极地海冰对大气和海洋环流有着重要的影响, 在全球气候变化中起着重要作用^[1]。海冰的存在改变了大气和海洋之间的热交换, 影响了海洋的热平衡, 阻碍了大气和海洋之间的动量传递。与开阔水域相比, 海冰的反照率更大, 海冰将更多的太阳辐射反射回空间, 吸收较少的短波辐射。特别是海冰范围的季节性变化很大, 对全球气候系统具有重要影响^[2-3]。随着全球变暖, 卫星观测显示, 自1979年以来, 北极海冰密集度(Sea Ice Concentration, SIC)一直在下降^[4]。北极海冰在2022年9月18日达到其年度最小面积, 海冰面积缩小到4.67×10⁶ km², 比1981—2010年的平均最小面积6.22×10⁶ km²少了约1.55× 10⁶ km²。极地海冰的减少是国际社会面临的一个紧迫问题, 需要进一步优化海冰数据集的精度、时间序列的连续性和空间分辨率的一致性。SIC是海冰变化的主要参数之一, 海冰的变化与温度、风场、海水温度和盐度密切相关^[5]。SIC被定义为给定区域海冰的面积分数, 能反映海冰的空间分布特征。通过这个参数, 我们可以对大气和海洋预测等业务进行研究^[6]。因此, 研制高质量的SIC数据, 对于极地研究、人文活动等方面都有着重要的意义^[7-9]。

目前, 国际上发展出以微波遥感监测为主, 光学遥感监测为辅的海冰监测体系。在海冰研究中, 微波遥感和光学遥感各有优势。光学遥感海冰数据具有丰富的多光谱信息且具有较高的空间分辨率, 然而其易受到大气衰减和天气条件的影响^[10-11]。相比之下, 微波遥感可以穿透云和雨, 具有全天候的海冰监测能力且可以获得丰富的纹理特征, 但其分辨率较低。因此, 充分利用不同数据源数据的特征信息, 实现优势互补, 克服单数据源的局限性, 本研究提出了一种基于微波遥感和光学遥感的多源SIC融合方法。

针对多源SIC数据的融合, 2001年, 吴奎桥等^[12]基于NOAA/ARHRR和Radarsat/ScanSAR数据, 利用小波多分辨率分析的方法对辽东湾海域的海冰遥感数据进行了融合分析, 证明了此融合方法可以有效提高海冰遥感图像的解析能力。2003年, Shokr^[13]提出了一种将粗分辨率传感器的每个网格点与细分辨率传感器的相应网格点的视场共同定位的方法, 有助于解决Radarsat海冰图像视觉分析中检索海冰参数的一些模糊问题。2007年, Helfrich等^[14]采用人机交互的方法, 将冰雪数字化图表等数据与多种卫星数据(可见光, 被动微波, 主动微波)进行融合, 为监测北半球海冰研制了IMS(Ice Mapping System)海冰覆盖数据。2009年, Yu等^[15]基于快速散射仪和先进微波扫描辐射计海冰图像数据, 利用主成分分析的方法进行融合, 使得北极西部海冰类型数据精度得到提高。2015年, Liu等^[16]利用小波变换方法、HIS (Hue-intensity-saturation)变换方法分别对海冰SAR (Synthetic Aperture radar)图像和光学图像进行了处理, 从而保留了原始图像的空间细节和光谱信息, 然后利用主成分分析方法得到融合图像, 有效地融合了海冰SAR图像和光学图像的优点。2021年, Yan等^[17]基于卷积神经网络在提取深层特征方面的优势, 针对SAR图像和光学图像的卷积神经网络结构, 提取异质多尺度特征和多层次特征, 利用特征级融合实现海冰图像分类, 设计了一种深度学习和异构数据融合的海冰图像分类方法, 明显提高了海冰图像分类精度。同年, 张雷等^[19]采用了多尺度三维变分的融合方法^[18], 研制出了全球逐日0.25°分辨率SIC融合试验数据, 并研发了全球SIC融合业务软件(v1.0)。2022年, Tsung等^[20]提出了一种极地自主船舶近场海冰感知的雷达图像融合处理方法, 将海洋雷达和冰雪雷达的图像进行融合, 得到了一套适用于低海冰浓度和高海冰密集度情况的近场海冰风险评估和预警流程。

已有的SIC数据融合研究主要针对小尺度海域, 且数据时间连续性较差, 可应用范围较窄, 针对大尺度海域SIC数据融合的研究, 其研究成果较少, 资料分辨率较低, 误差较大, 不能满足北极海冰监测, 预报等实时业务的需要。为解决此类问题, 本文提出了一种北极多源SIC融合方法, 采用OSI SAF等多种数据源, 通过概率密度匹配法(Probability Density Function, PDF)订正各数据源的系统误差, 利用小波变换结合卡尔曼(Kalman)滤波的融合算法得到2016年空间分辨率为0.05°的逐日融合资料, 并与OISST和OSTIA SIC数据进行对比, 验证了产品的可靠性。本文结构如下, 第一节介绍了数据和预处理方法, 第二节介绍本文提出的SIC融合方法, 第三节对融合结果进行了对比验证, 第四节对研究结果进行讨论与总结。

1 资料介绍

考虑到SIC数据源资料需具备实时性和可获取性。因此本论文采用了欧洲气象卫星应用组织(European Organization for the Exploitation of Meteorological Satellites, EUMETSAT)、美国国家冰雪数据中心(National Snow and Ice Data Center, NSIDC)、中国国家卫星气象中心(National Satellite Meteorological Centre, NSMC)等的数据资料。EUMETSAT研制的OSI SAF海冰密集度资料^[21]是基于美国DMSP(Defense Meteorological Satellite Program)系列卫星上的SSMIS(Special Sensor Microwave Imager Sounder)传感器, 使用ECMWF(European centre for Medium-Range Weather Forecasts)数值模式结果进行大气校正, 并采用OSHD (OSI SAF Hybrid Dynamic)算法, 生成空间分辨率为10 km的SIC资料, 覆盖了南北半球高纬度地区。NSIDC研制的NISE (Near-Real-Time Ⅰ-SSMIS EASE-Grid) SIC资料^[22]也是以美国DMSP系列卫星上的SSMIS传感器为基础, 通过NASA(National Aeronautics and Space Administration) Team方法, 得到空间分辨率为25 km的SIC资料。NSMC研制的微波成像仪(Microwave Radiation Imager, MWRI)和可见光红外扫描辐射计(Visible and Infrared Radiometer, VIRR) SIC资料, 分别以FY-3号系列极地轨道气象卫星上的MWRI和VIRR传感器为基础^[23], 得到空间分辨率为12.5 km的MWRI极区SIC日产品和空间分辨率为0.1°的VIRR SIC日产品。NIC研制的IMS海冰覆盖资料在本研究中, 主要为北极极点区域缺失数据进行补充, 在前文已作详细介绍。各资料如表 1所示。

2 方法

本文对各资料进行预处理之后, 通过PDF匹配法订正各资料的系统误差, 然后利用小波变换将各资料分解为低频信息和高频信息, 低频信息计算各资料的融合权重进行融合, 高频信息利用Kalman滤波进行处理筛选, 然后利用小波反变换得到融合资料。流程图如图 1所示。

2.1 数据预处理

由于多源SIC数据存在资料类型、分辨率、阈值范围和投影方式等方面的不同, 甚至有的数据源存在缺失经纬度信息的情况, 因此在进行实验之前需要对各数据源进行预处理。首先, 对不同类型的数据源进行读取; 其次, 对没有经纬度信息的数据进行处理, 确保所有数据都有经纬度信息; 然后, 将所有的SIC数据源阈值范围统一调整至0~100%; 最后, 采用最近邻点插值法, 将所有数据源统一至0.05°分辨率。

2.2 卫星资料的误差订正

由于SIC数据存在着明显的非独立性系统误差, 即误差会随着资料值大小、空间范围等的不同而变化。针对此类系统误差PDF匹配法可以较好地订正^[24]。因此, 本研究以OSI SAF SIC资料为基准, 对其他资料进行订正。

从表 2的均方根误差(RMSE), 偏差(Bias)和相关系数(CC)可以看出, 经过PDF匹配法订正后, NISE、MWRI和VIRR SIC数据的偏差和均方根误差都得到了显著的降低, 说明该方法能有效地降低各数据源的非独立性系统误差。

2.3 融合权重的确定

由于各数据源的精度不同, 所以每个数据源在融合时的权重也不同, 因此就需要计算每个数据源的融合权重, 原理为: 数据源的误差越小, 一致性测度数值就越大, 权重值也就越大。所以通过计算数据源的一致性测度值, 就可以得到低频信息的融合权重值。

假设有n个传感器, 则传感器i在时刻Q处所获得的值为X_i(q) (q=1, 2, …, s; q代表像素的编号; s是观察区域像素的数量), 并且在Q时刻传感器j获得的值为X_j(q), 此时, 在q位置上的传感器i和j的支撑力度矩阵表示为:

$ {\delta _{ij}}(q) = \frac{{\min \left\{ {{x_i}(q),{x_j}(q)} \right\}}}{{\max \left\{ {{x_i}(q),{x_j}(q)} \right\}}} . $

(1)

在时间是Q的情况下, 每个传感器在像素位置q处的总支撑值表示为:

$ \sum {(q) = \left( {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1 \\ {{\delta _{21}}(q)} \\ \vdots \\ {{\delta _{n1}}(q)} \end{array}}&{\begin{array}{*{20}{c}} \cdots \\ \cdots \\ \vdots \\ \cdots \end{array}}&{\begin{array}{*{20}{c}} {{\delta _{1n}}(q)} \\ {{\delta _{2n}}(q)} \\ \vdots \\ 1 \end{array}} \end{array}} \right)} . $

(2)

在整个支撑矩阵中, 各行支撑值的总和表示相应传感器的最终支撑值。所以, 在时刻Q的情况下, 每个传感器在像素位置q的一次性测度表示为:

$ {r_i}(q) = \frac{{\sum\limits_{j = 1}^n {{\delta _{ij}}(q)} }}{n} , $

(3)

式中, 一致性测度r_i(q)是指在q像素位置, 传感器i得到的数据与其他传感器得到的数据的接近程度。单独一个像素的r_i(q)值, 并不能完全反映观察区域传感器的可靠度, 所以应将观察区域内所有像素的一致性测度计算在内。因此, 在时间为Q时, 第i个传感器的支撑值的均差和方差分别表示为:

$ {R_i}(Q) = \frac{{\sum\limits_{q - 1}^s {{r_i}(q)} }}{s} , $

(4)

$ \delta _i^2(Q) = \frac{{\sum\limits_{q = 1}^s {{{\left[ {{r_i}(q) - {R_i}(Q)} \right]}^2}} }}{s} . $

(5)

如果一个传感器观察区域内的像素平均值越大, 并且其灰度值越低, 那么这个传感器就越稳定, 其支撑值就越高, 权重就越大。所以, 这个传感器的支撑值可以用信号噪声比来形容, 用平均值除以方差就是信号噪声比的值。当时间是Q时, 传感器i的权重值ω_i(Q) 表示为:

$ {\omega _i}(Q) = \frac{{{R_i}(Q)}}{{\delta _i^2(Q)}} . $

(6)

归一化之后的值表示为:

$ {W_i}(Q) = \frac{{{\omega _i}(Q)}}{{\sum\limits_{j = 1}^n {{\omega _j}(Q)} }} . $

(7)

2.4 小波变换与Kalman滤波技术

本文使用的多源SIC数据融合技术是采用小波变换结合Kalman滤波来实现。小波变换能够将原始图像分解为具有不同频率特性、空间分辨率和方向特性的子图像, 从而能够更好地保留原始图像的光谱信息^[25]。该方法的基本思想类似于一对共轭镜象正交滤波器组, 将影像分解成1条低频带和3条高频带(水平, 垂直, 对角线3条)^[26], 其中, 低频段代表影像的总体特性, 而高频段则代表影像的详细特性。同时, 小波变换算法中的双正交或正交多分辨率融合算法已经广泛应用于多传感器图像融合领域, 其融合效果得到了广泛认可。在本文提出的融合技术中, 我们将Kalman滤波应用于小波变换后的高频部分。Kalman滤波是一种高效的自回归滤波器, 能够寻求最优的预测结果, 在效率、通用性等方面都较好。

本文选取Haar小波基, 经过小波变换, 将各数据源分解为1个低频信息和3个高频信息, 这样做可以将原始数据分解成不同层数的小波系数, 从而实现对信号的频域分析, 同时也可以降低数据噪声和冗余, 提高数据的稳定性和可靠性, 分解后资料的分辨率为原来的一半。然后根据式(7)计算权重值, 按权重值加权处理低频信息, 留作小波反变换的低频信息进行融合; 接着对高频信息进行Kalman滤波处理, 假设有n个数据源, 进行Kalman滤波处理后总共生成3×n个高频信息, 在水平、垂直和对角线3个方向上, 每个方向都有n个高频信息, 按照最大值原则对3组高频信息进行筛选, 然后得到3个新的高频信息; 最后用小波反变换将1个新的低频信息和3个新的高频信息进行融合, 得到一个融合资料。算法对每天的数据都作上述处理; 最后得到相应的逐日融合资料。

3 融合结果

通过与国际上认可的OISST、OSTIA SIC数据相比较, 对融合资料的质量进行了检验。如图 2所示, 为2016年1月1日该融合资料与OISST、OSTIA SIC数据的空间分布图。由于北半球海冰范围在3月最大, 9月最小, 因此选取3月1日、9月1日海冰密集度融合数据与同期OISST、OSTIA资料进行对比, 图 3所示是该日融合资料与OISST和OSTIA资料差值分布, 结果表明, 该方法所得结果在空间分布上与OISST、OSTIA SIC数据有较好的一致性, 且与北半球SIC较高的地区基本吻合。但是, 3种类型的SIC数据在一些低SIC的边缘地区, 其密集度值有较大区别。这是因为各传感器对于低SIC区域的反演能力较差, 所以同日在NISE、MWRI、OSISAF和VIRR SIC资料中都未在相应海域出现海冰覆盖情况, 因此, 在该区域不同资料之间存在较大差异。

图 4为2016年度融合实验数据和OISST、OSTIA SIC数据的逐日对比, 结果表明, 该数据与OSTIA及OISST SIC数据的年度偏差在–1%~1%, 与OSTIA相比在6—9月存在正偏差, 其他时间均为负偏差, 从均方根误差和相关系数的角度看, 相对于OISST SIC数据, 该融合数据与OSTIA SIC数据相关性更高。

从表 3可以定量看出, 在2016年度, 融合试验数据与OISST数据的匹配样本为6 280 523 594个, 偏差、均方根误差和相关系数分别为–0.108%、8.010%和0.967。与OSTIA匹配的样本量达到6 293 207 725个, 与OSTIA的偏差为–0.156%, 均方根误差为5.140%, 相关系数达到0.989。且相对于前人的研究^[17], 本融合资料与OISST的偏差由–1.170%减少到–0.108%, 与OSTIA的偏差由0.276%减少到–0.156%; 与OISST和OSTIA的均方根误差分别由9.835%减少为8.010%以及由7.427%减少为5.140%。本融合资料的偏差以及均方根误差都得到了显著的提升, 具有较高的质量。为了验证数据的可靠性, 另做了2019年的融合资料, 如表 4所示, 可以看出融合资料在偏差与均方根误差方面都有着不错的结果。

4 结论

本文提出了一种北极多源SIC资料融合方案, 即首先, 每一个数据源都要经过质量控制; 利用PDF匹配方法对各个数据源进行系统误差校正; 然后, 利用小波变换结合Kalman滤波的融合方法将各订正后的数据源进行融合; 最后, 生成2016年北极逐日0.05°分辨率的海冰密集度融合资料。将融合资料与OISST和OSTIA数据进行对比, 得出以下结论: 小波变换结合Kalman滤波的融合方法, 在低频分量上抓住海冰密集度空间分布的整体特征, 在高频分量上抓住海冰数据的细节特征, 能够有效的提高融合资料的质量。该方法充分考虑了各资料的时效性和可获取性, 保证了融合资料的实时稳定产出。

融合资料与OISST和OSTIA SIC资料的偏差分别为–0.108%和–0.156%, 均方根误差分别为8.010%和5.140%, 相关系数均超过0.967。且与传统的SIC融合产品相比, 本融合资料与OISST的偏差由–1.170%减少到–0.108%, 与OSTIA的偏差由0.276%减少到–0.156%; 与OISST和OSTIA的均方根误差分别由9.835%减少到8.010%以及由7.427%减少为5.140%, 偏差和均方根误差都得到了明显的提升。

该资料在空间分布上与OISST、OSTIA资料有很好的相似性, 且与北半球SIC高的地区基本吻合, 验证了该融合结果具有较高的精度, 但是在北极低SIC区域, 3组数据之间有较大差别, 有待今后更深入地探讨。