1 引　言

极地是地球热机冷源地，极地海冰的变化通过海洋和大气影响全球能量和物质的循环。随着人们对海冰认识的不断深入，极地海冰研究已成为气候变化和可持续发展研究的热点之一（IPCC，2013；秦大河等，2014）。北极和南极海冰的变化通过气候系统成员间的相互作用，可以影响欧亚大陆特别是中国的天气、气候（卞林根等，2008；武炳义等，2011）。海冰密集度（Sea Ice Concentration，SIC）作为描述极地海洋表层热状况最为重要的参数之一，被广泛地应用于上层海洋过程、海气热量交换、海洋大气数值模拟与预报等研究和业务工作中（国巧真等，2006；蒋兴伟等，2008）。因此，研制高质量的海冰密集度数据，对于研究极地环境和大气环流等具有非常重要的意义（效存德，2008；马丽娟等，2004；张雷等，2017；张琦等，2018；Hui，et al，2017）。

国际上针对科研、业务需要，研发了多种高质量的海冰密集度产品，其中美国国家海洋和大气管理局（NOAA，National Oceanic and Atmospheric Administration）研制的OISST （Banzon，et al，2016；Reynolds，et al，2007）海冰密集度数据具有时间序列长、时效性高等特点，被广泛使用。在海冰密集度数据的处理过程中，早期（1981—2004年）数据源使用了美国国家航空航天局戈达德航天中心（NASA's Goddard Space Flight Center）的海冰密集度产品（Cavalieri，et al，1999），近期（2005年至今）数据源使用NCEP Gridded Ses Ice Concentra-tion产品（Grumbine，1996）。采用中值滤波方法生成了全球0.25°逐日海冰密集度产品。英国气象局（Met Office，National Meteorological Office）的OSTIA数据（Donlon，et al，2012），使用了OSI SAF被动微波海冰密集度卫星产品（Tonboe，et al，2017），通过双线性插值等方法，生成了全球0.05°逐日海冰密集度产品。英国Met Office Hadley Centre为气候变化研究研制了HadISST资料（Titchner，et al，2014），其中海冰资料源包括数字化图表和被动微波卫星数据，为保证融合资料的均一性，利用同期的数字化图表订正微波遥感数据，使卫星资料和图表资料统一，在无观测数据的网格上使用气候态海冰资料进行插补。美国国家冰中心（NIC）为北半球海冰监测研制的IMS（Helfrich，et al，2007）海冰覆盖数据，主要通过人机交互的方式融合了多种卫星资料（可见光、被动微波、主动微波）和冰雪数字化图表等数据源资料。IMS的海冰数据代表网格的地表类型而非该网格内的海冰密集度数据。

从数据源看，在OISST和OSTIA资料中使用的海冰密集度数据均为被动微波卫星遥感数据产品，在HadISST和IMS资料中使用了更多的数据源。数据处理方法上，对于单一数据源资料的处理相对简单。目前，大尺度观测海冰主要利用卫星遥感手段，包括可见光、红外、被动微波和主动微波遥感等。可见光遥感分辨率高，但在夜晚、极夜不能对海冰进行监测，另外云层对可见光和红外遥感的影响非常大。微波遥感数据具有不受昼夜限制，受云雾影响较小以及较好的时、空连续性等特点，但被动微波遥感目前分辨率较低，主动微波分辨率高却存在回访周期长和成本较高等不足。因此，每种观测方式都有其长处与不足，开展海冰密集度多源数据融合，可以吸收各种数据源的优势提高数据质量，也可以增强数据的稳定性（师春香等，2019；潘旸等，2018；Shen，et al，2014）。

为满足海冰监测、数值天气预报等实时业务与科研的需要，开展了逐日全球多源海冰密集度资料的融合试验。本研究发展了一套全球多源海冰密集度融合方案，采用包括风云气象卫星产品在内的多种数据源，通过综合考虑各数据源的优势和误差特征，订正各数据源资料的系统误差并生成超级观测场，利用多尺度变分技术将超级观测场与模式背景场融合生成2016年逐日0.25°×0.25°分辨率的融合试验产品，并与国际上广泛认可的海冰密集度数据产品对比评估。融合方案考虑了卫星观测数据的实时可获取性，以满足实时业务需要。

2 资料和预处理方法

根据海冰密集度融合试验的目的，融合数据源需要考虑资料的实时性和可获取性。全球海冰密集度融合数据源采用了欧洲气象卫星应用组织海洋海冰应用中心OSI SAF的海冰密集度卫星产品、美国国家冰雪数据中心（NSIDC）的NISE海冰密集度资料、中国国家卫星气象中心（NSMC）的MWRI和VIRR海冰密集度卫星资料、美国国家冰中心（NIC）的IMS海冰覆盖资料以及ERA-Interim的海冰密集度数据（表1）。

OSI SAF海冰密集度产品（Tonboe，et al，2017）是基于DMSP（Defense Meteorological Satellite Program）系列卫星上SSMIS（Special Sensor Microwave Imager Sounder）传感器，采用OSHD（OSI SAF Hybrid Dynamic）算法，并使用ECMWF数值模式结果进行大气校正，得到的空间分辨率为10 km的海冰密集度数据产品，覆盖范围为南、北半球高纬度区域。NISE数据（Brodzik，et al，2016）同样是由搭载于美国DMSP系列卫星上SSMIS传感器的数据，采用NASA Team算法生成的分辨率为25 km的海冰密集度产品。MWRI微波成像仪和VIRR可见光红外扫描辐射计两个传感器均搭载在FY-3系列极轨气象卫星上（杨军等，2011）。相应海冰密集度数据产品由中国气象局国家卫星气象中心提供，MWRI极区海冰密集度产品为12.5 km分辨率通用横球面投影极区日产品，VIRR海冰密集度产品为0.1°分辨率等经纬度日数据。IMS高分辨率北半球海冰密集度融合资料前文已经进行了详细介绍，在该次试验中其主要用于极点附近缺少卫星观测区域的数据补充。ERA-Interim资料（Dee，et al，2011）中的海冰密集度数据作为融合背景场。

数据在融合前，首先要对各网格点资料进行质量控制，包括（1）空间检查：核对数据的经纬度信息，并利用融合背景场海陆分布数据检查观测点是否在海洋上；（2）时间检查：保证其时间在24 h范围内；（3）阈值检查：海冰密集度记录应该在0%至100%范围内。多种来源的资料存在资料类型、分辨率和投影方式等差异，需要进行统一。对不同来源不同分辨率的数据进行插值（钱辉等，2009），将等面积数据转换为等经纬度数据。在此过程中，利用网格面积加权平均的方法将判识值统一为范围值，生成0.25°分辨率的单数据源观测场。

3 融合试验 3.1 海冰密集度卫星资料的系统误差订正

海冰密集度卫星资料的系统误差表现为明显的非独立系统误差特征。这种系统误差随观测值的大小、空间区域和时间范围而变化，通常的线性方法难以订正，而概率密度函数（PDF）匹配法能够较好地订正此类系统误差（宇婧婧等，2013；Li，et al，2015）。中外研究机构对OSI SAF海冰密集度卫星资料的评估表明，该资料的精度和空间分辨率均较高（Tonboe，et al，2017；Lavelle，et al，2017；Lam，et al，2018；Blunden，et al，2020；郝光华等，2018；赵杰臣等，2017），且获取途径相对稳定，OSI SAF海冰密集度数据可作为参考，订正其他卫星资料的系统误差。

利用概率密度函数匹配方法，以OSI SAF资料为基准，通过调整非基准资料，使两种资料的概率密度分布一致，从而实现对NISE、MWRI和VIRR卫星资料的系统误差订正。订正方案中匹配样本的时间滑动窗口采用的是当天并向前取29 d（共30 d）；匹配样本的空间滑动窗口是以分析点为中心向外40个网格内所有有效格点；要求时、空窗口内匹配样本超过300对才进行订正。

从系统误差订正前、后观测−参考的概率分布（图1）可以看出，PDF订正后的曲线更倾向于正态分布，说明系统误差订正具有较好的效果。从订正前、后定量统计结果（表2、表3）可以看出，系统误差订正后各数据源的偏差（Bias）和均方根误差（RMSE）较订正前明显减小，说明PDF订正方法能够有效地减小海冰密集度卫星产品的系统性误差。

3.2 超级观测场的生成

超级观测场“Super Observation”（Reynolds，et al，2007）是将不同来源的多种观测资料进行有效结合，生成统一的观测场。由于海冰密集度观测种类较多，且分布非常不均匀，融合时易产生求解的“不适定问题”，即紧邻的观测点的观测值相差较大，求解通常为病态解，而超级观测场方案可以有效解决这一问题。

根据系统误差订正后各数据源误差的大小，生成超级观测场。确定超级观测每个网格海冰密集度值的策略如下：对于每个分析网格，如果观测误差相同，则超级观测值是其平均观测值；如果观测误差不同，则取误差小的观测值作为超级观测值，误差大的观测值则丢弃。计算公式如下

$ {S_{{{ij}}}} = \left\{ {\begin{array}{*{20}{l}} {v^{{k}}_{{{ij}}} }&{\;\varepsilon ^{{k}}_ {{{ij}}} = {\rm{min}}\left({\varepsilon ^ 1_{{{ij}}}\,,\,\varepsilon^2_{{{ij}}} \,, \cdots \cdots,\,\varepsilon^{{{K}} - 1}_{{{ij}}} \,,\,\varepsilon ^{{K}}_ {{{ij}}} } \right)}\\ {\dfrac{\,1\,}{\,P\,}{{\displaystyle\mathop \sum \limits_{{{p}} = 1}^{{P}} v^{{p}}_ {{{ij}}} }}}&{\underbrace {\varepsilon^{{m}}_{{{ij}}} = \cdots \cdots = }_{P {\text{种}}}\varepsilon^{{r}}_ {{{ij}}} = {\rm{min}}\left({\varepsilon^1_ {{{ij}}} \,,\,\varepsilon^2_{{{ij}}} \,, \cdots \cdots,\,\varepsilon^{{{K}} - 1}_ {{{ij}}} \,,\,\varepsilon^{{K}}_ {{{ij}}} } \right)} \end{array}} \right. $

(1)

式中，S_ij表示第i行第j列网格的超级观测场的值，K为观测资料种类数，k、m、r为资料的种类，当第k类资料在该网格的误差 $\varepsilon _{{{ij}}}^{{k}} $ 最小时，则此类资料的观测值 $ v _{{{ij}}}^{{k}} $ 作为该网格的超级观测值。当在同一网格存在P种资料的误差相同时，则超级观测值是P种资料的算术平均值 $ \dfrac{\,1\,}{\,P\,}{{\displaystyle\mathop \sum \limits_{p = 1}^{{P}} v^{{p}}_{{{ij}}}}} $ 。另外，由于卫星轨道与地轴之间存在一定的夹角，常用的极轨卫星观测不到极点周围的海冰密集度，利用IMS数据填补极点附近的超级观测场。

3.3 多尺度三维变分方法融合

利用多尺度三维变分（STMAS）方法（Xie，et al，2011）将超级观测场和ERA-Interim海冰密集度背景场进行融合，生成全球海冰密集度融合试验产品。STMAS方法通过依次对不同尺度的观测信息进行分析，消除了长波信息与短波信息的混淆，能够克服由背景场误差协方差矩阵不准确引入的误差。这相当于通过改变网格的粗、细来调整观测信息的影响范围，从而起到了调整背景场误差协方差矩阵的作用，该方法将背景场误差协方差矩阵设为简单的单位矩阵，从而降低对运算资源的要求。STMAS方法在粗网格的目标函数对误差长波进行分析，细网格的目标函数对误差短波进行分析，最终分析结果为各重网格分析结果的叠加。该方法在大尺度的计算上节约了时间，具有较高的计算效率，并能够详细提取小尺度的信息，具有较低的均方根误差。该融合技术具有计算高效、运行稳定等特点（Xie，et al，2011；Li，et al，2008；徐宾等，2018；韩帅等，2018）。在每一重网格上的目标泛函采用如下形式

$\begin{split} {J^{({{n}})}} =& \frac{\,1\,}{\,2\,}{{\boldsymbol X}^{{{({{n}})}^{\rm{T}}}}}{{\boldsymbol X}^{({{n}})}} + \frac{\,1\,}{\,2\,}{\left[ {{{\boldsymbol H}^{({{n}})}}{{\boldsymbol X}^{({{n}})}} - {{\boldsymbol Y}^{({{n}})}}} \right]^{\rm{T}}}{{\boldsymbol O} ^{{{({{n}})}^{ - 1}}}}{\text •}\\&\left[ {{{\boldsymbol H}^{({{n}})}}{{{\boldsymbol X}}^{({{n}})}} - {{\boldsymbol Y}^{({{n}})}}} \right]\quad{{n}} = 1,2,3 \cdots ,{{N}} \end{split}$

(2)

式中，J为目标函数，X=X_a−X_b，Y=Y_obs−HX_b，O是观测误差协方差矩阵，H是从模式网格到观测点的双线性插值算子，n表示第n重网格，N为网格重数，T表示矩阵的转置。其中，X_a是分析场，X_b为背景场，Y_obs为观测场，X为相对背景场的修正，它由变分同化系统计算而来，Y是观测场与背景场的差值，不同尺度网格Y值的表达式如下：

$ \left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol Y}^{(1)}} = {{\boldsymbol Y}_{{\rm{obs}}}} - {{\boldsymbol H}^{(0)}}{{\boldsymbol X}_{\rm{b}}}}&{{{n}} = 1}\\ {{{\boldsymbol Y}^{({{n}})}} = {{\boldsymbol Y}^{({{n}} - 1)}} - {{\boldsymbol H}^{({{n}} - 1)}}{{\boldsymbol X}^{({{n}} - 1)}}}&{{n}} = 2,3 \cdots,{{N}} \end{array}} \right. $

(3)

多尺度变分STMAS方法从最粗网格到细网格依次进行分析，当式（3）中n=1时， $ {{\boldsymbol H}}^{\left({0}\right)} $ 表示从背景场到观测点的双线性插值算子， $ {{\boldsymbol Y}}^{\left({1}\right)} $ 为观测和背景场到观测位置映射的差值，代入式（2）中求得 $ {{J}}^{\left({1}\right)} $ 的极小化解 $ {{\boldsymbol X}}^{\left({1}\right)} $ ；当n≥2时，上一层 $ {{J}}^{\left({n}{-1}\right)} $ 的极小化解 $ {{\boldsymbol X}}^{\left({n}{-1}\right)} $ 被插值到更细的第n层网格上， $ {{\boldsymbol Y}}^{\left({n}{-1}\right)} $ 可从上一层获得，从而计算得到第n层网格的 $ {{\boldsymbol Y}}^{\left({n}\right)} $ ，代入式（2）求得 $ {{J}}^{\left({n}\right)} $ 的极小化解（ $ {{\boldsymbol X}}^{\left({n}\right)} $ ）。通过迭代计算，求得每重网格分析增量（ $ {{\boldsymbol X}}^{\left({n}\right)} $ ）的值。由式（4）得到分析场（X_a），即为背景场（X_b）与各重网格分析增量叠加 $\displaystyle\mathop \sum \limits_{{{n}} = 1}^N {{\boldsymbol{X}}^{{{(n)}}}}$ 的和。

${{\boldsymbol{X}}_{\rm{a}}}{\rm{ = }}{{\boldsymbol{X}}_{\rm{b}}}{\rm{ + }}\mathop \sum \limits_{n = 1}^{{N}} {{\boldsymbol X}^{(n)}} $

(4)

4 融合结果 4.1 全球范围评估

将融合试验产品与OISST和OSTIA海冰密集度资料进行对比，评估融合试验产品的质量。图2给出了2016年2月1日融合试验产品与OISST和OSTIA海冰密集度资料空间分布的对比。从图中可以看出，融合试验产品与OISST和OSTIA的海冰数据在空间分布上高度相似，在北半球和南半球高海冰密集度区域几乎一致。但3种海冰密集度资料在部分海冰边缘的低海冰密集度海域存在一定的差异，该日在南极地区威德尔海海域OISST海冰资料的低海冰密集度范围明显大于OSTIA和融合试验产品。在此海域，这种大范围的低海冰密集度覆盖在同日OSI SAF、NISE、MWRI和VIRR海冰密集度数据中均未发现。对于低海冰密集度区域各种传感器观测数据的反演能力均较差，造成了海冰边缘海域低海冰密集度区域卫星资料的不确定性较大。因此，在该区域不同产品之间存在较大差异。

北半球3月前后海冰范围最大，9月前后海冰范围最小，因此选取了3月15日、9月15日海冰密集度融合试验产品与同期OISST和OSTIA资料进行对比，分析融合试验产品与两种常用资料空间分布的差异。图3和4分别为2016年3月15日和9月15日北极、南极地区融合试验产品与OISST、OSTIA海冰密集度差值的空间分布。北极地区融合产品与OISST和OSTIA在3月和9月主要差别出现在海冰边缘区域，在空间分布上融合试验产品与OSTIA资料更接近。南极（图4）9月融合产品与OISST相比低海冰密集度区域有明显的负偏差，但与同期的OSTIA相比在边缘区域还出现了正偏差，3月融合产品与OISST和OSTIA资料差异的空间分布存在偏差相反的情况。整体看，融合产品与两种资料的空间差异主要表现在海冰边缘区域，空间差异也随时间发生变化，融合产品与OSTIA资料一致性更高，这与融合试验产品使用的数据源及偏差订正选取的标准有关。

图5为融合产品与OISST、OSTIA海冰密集度资料2016年逐日的比较结果。可以看到，融合试验产品与OSTIA全年的偏差均在0值附近，与OISST相比6—8月为负偏差，11月—12月中旬存在正偏差。在均方根误差和相关系数方面，融合产品与OSTIA海冰密集度资料相关要高于与OISST的相关，其主要差别出现在7—11月。

定量评估显示（表4），2016年全年融合试验产品与OISST匹配的样本量达到44321040个，偏差为−1.170%，均方根误差为9.835%，相关系数为0.985。与OSTIA的偏差仅为0.276%，均方根误差为7.427%，相关系数达到0.999。融合试验产品较OSTIA为正偏差，与OISST为负偏差，试验产品整体偏差介于二者之间。

4.2 关键区海冰密集度的评估

巴伦支海—喀拉海—拉普捷夫海海冰通过海冰、海水和大气之间的相互作用，对中国的气候产生明显的影响，成为中国海冰监测的关键区。关键区域覆盖的范围定义为（72°—83°N，30°—150°E）（Zuo，et al，2016），紧靠亚欧大陆北端。

3和9月前后分别是北极海冰覆盖范围最大和最小值时期。图6为2016年巴伦支海—喀拉海—拉普捷夫海海冰关键区融合试验产品、OISST和OSTIA海冰密集度数据3、4、9和10月各月15日的空间分布。从图中可以看到，融合产品与OISST、OSTIA资料都能反映出巴伦支海—喀拉海—拉普捷夫海区域海冰密集度的变化，3种资料均能很好地展现北极海冰的气候特征。在巴伦支海—喀拉海—拉普捷夫海区域，3种海冰密集度资料的空间分布具有较高的一致性，但在部分陆地边缘海域的低海冰密集度区域也存在差异。

在海冰关键区的年内变化中（图7），融合试验产品与OSTIA资料的偏差一直较小，与OISST资料相比在5月中旬至9月中旬存在一定的负偏差，9月下旬至11月存在正偏差。试验产品与OSTIA资料的RMSE较小，与OISST资料相比，6月至12月的RMSE值较大。海冰关键区，融合试验产品与OSTIA海冰密集度资料的相关整体也高于与OISST资料的相关。

表5定量给出了关键区2016年逐日融合试验产品与OISST、OSTIA海冰密集度资料的对比分析。与OISST和OSTIA的偏差分别为−2.155%和−0.475%，相应的均方根误差为11.560%和9.770%。试验产品与OSTIA海冰密集度资料的相关系数为0.924。

针对气候关键区开展的评估分析表明，融合试验产品在该区域与OISST和OSTIA的空间分布具有较高的一致性，与两种资料的偏差较小，在年内变化上试验产品与OSTIA的差异较小且具有很好的稳定性。

5 结　论

本研究设计了一套全球多源海冰密集度数据融合方案，即首先对各数据源资料进行质量控制，利用概率密度函数法订正数据源的系统误差，根据订正后各数据源误差的大小生成超级观测场，最后利用STMAS法将超级观测场和背景场融合，生成融合试验产品。基于此方案研制生成了2016年全球逐日0.25°分辨率的海冰密集度融合试验产品，将融合产品与国际广泛使用的海冰密集度资料进行对比评估。得到以下结论：

融合方案中概率密度函数订正法能够有效地减小海冰密集度卫星资料的系统误差，是误差订正的一种有效方法。STMAS融合方法，在大尺度网格上抓住海冰密集度空间分布的整体特征，在小尺度网格上提取细节信息，不但实现了超级观测场与背景场的融合，还提升了计算效率。本方案考虑了各卫星数据源的时效性、可获取性，能够保障实时产品稳定生成。因此，上述方案的提出对于全球海冰密集度实时产品的研制具有较好的适用性。

融合试验产品与OISST、OSTIA海冰密集度资料的偏差分别为−1.170%和0.276%，试验产品整体偏差介于两种对比资料之间，产品相关系数均超过0.980。融合试验产品整体偏差介于OISST与OSTIA之间，反映了试验产品系统误差较小的良好特性。试验产品、OISST和OSTIA三种资料在南、北极海冰密集度空间分布上高度一致，表明融合试验产品具有较高的精度。

文中探讨了一种适用于多源海冰密集度数据的融合方案，基于该方案研制的全球逐日海冰密集度融合试验产品与OISST和OSTIA资料高度一致，但在海冰边缘区3种资料之间均存在一定差异，对此问题还需要更长期的检验与改进工作。

根据文中所述融合方案，已研发了全球海冰密集度融合业务软件（v1.0），由该软件实时生成全球0.25°×0.25°逐日海冰密集度融合实况分析产品，上述产品可通过国家气象业务内网http://10.1.64.154/idata/获取。

致　谢：感谢国家卫星气象中心武胜利正研级高级工程师、郑照军高级工程师和中山大学杨清华教授在海冰密集度资料处理方面给予的指导。