利用FY-3卫星MWRI数据探测海冰分布

吴展开，王星东, ，王 成

(1. 河南工业大学信息科学与工程学院，河南 郑州 450001；2. 中国科学院遥感与数字地球研究所，北京 100094)

北极是全球气候系统的重要组成部分，是全球气候变化的指示器。作为全球气候的冷源，其大气、海洋、海冰的运动直接或间接影响着全球水汽环流的强度、全球热平衡和气候变化，其中海冰是影响全球气候系统的关键因素[1-3]。近几十年由于全球变暖的影响，北极海冰呈逐年退缩、海冰厚度变薄、海冰范围的最小记录不断被刷新、夏季出现北极航道通航等现象，使得连接亚洲、欧洲和北美洲的北极航线有望全线开通。北极航线作为北美洲、欧洲和东亚地区之间潜在的最快捷的海上通道，具有重要的航运价值和战略价值。相比传统航线，北极航线可规避海峡拥堵及海盗等多种问题，且大大缩短我国至欧洲地区的航程，因此对我国具有重要的经济、军事和战略价值，也将有利于国家“一路一带”倡议的推进[4]。

北极航道的开通受多种因素的制约及影响，其中受海冰的影响极为严重[4-5]。另外，海冰动态与全球气候变化间的相互影响机制十分复杂，海-冰-气相互作用过程的反馈机理还不十分清楚。但是，通过对历史海冰观测资料和同期海冰参数的比较，进行逐旬、月、季和年间距分析可以获取海冰的动态变化信息；结合气象观测数据，可对海冰变化趋势做出预测。

由于微波辐射计能全天候工作，且受天气影响较弱，因此常用于海冰的研究。Andersen等[6]总结了针对 SSM/I 数据的7种海冰密集度的反演算法。其中，NASA-Team 算法[7]和Bootstrap算法[8]主要是基于低频数据对海冰密集度进行反演，这类方法受天气影响较弱，但空间分辨率也相对较低；NASA-Team2 算法[9-10]、SEA LION 算法[11]和 ARTSIST Sea Ice (ASI)[12]等算法均使用高频数据，可得到空间分辨率相对较高的结果。ASI算法具有不需要额外数据源的优点，而且与利用其他通道的海冰密集度算法有相似的结果[13]。同样，近年来国内学者，对海冰分布及海冰密集度的探测也做出了诸多研究成果，目前大多是基于AMSR-E数据及SSM/I进行[5,14]，而利用FY-3 MWRI数据对海冰的研究较少。

利用AMSR-E 数据基于ASI[15-16]、NASA-Team2 和Boot-strap Algorithm (ABA)[17]算法3种算法分别对北极海冰密集度进行反演，并将所得结果与2003年3—4月和2004年7—8月“北极星”号考察船船舶观测数据进行对比分析，结果表明3种算法得到的结果与观测数据的相关系数分别为0.80、0.79和 0.81。过去30年来，北极海冰范围减少明显，约以每年4.51×104km2的速度消退，且季节平均及月平均海冰范围均表现出下降趋势，4个季节中，夏季减少最为明显[18]。1979—2006年3月平均海冰范围以每10年2.8%的速度减少，而9月以每10年8.6%的速度减少[19]。通过对1978—2000年的微波数据进行分析研究，结果表明：夏季出现北极海冰范围最小时的海冰主要由多年冰构成，且前后两个10年间的平均海冰范围及平均海冰面积分别减少了6.3%和11.0%[20]。1979—2008年北极平均无冰季节长度表现出了较为明显的增加趋势，其中1979—2006年间，以无冰季节长度每年1.1 d的速度增长，2007和2008年无冰季节长度分别达到了168和158 d，且2007年的无冰季节长度为1979—2008年中的最高纪录[21]。

在分析海冰(一年冰、多年冰)与海水的辐射率变化特征的基础上，提出基于微波成像仪(MWRI) 19 GHz水平极化和垂直极化的极化差(极化差最大)结合Otsu算法(确定海冰和海水的分类阈值)来探测北极海冰分布的新方法，利用2016年1月共31 d的数据得到北极海冰月平均分布图。

1 研究区域与数据源

北极地区一般是指北纬66°34′以北的地区，由北冰洋、欧亚大陆的边缘陆地及众多的岛屿岛礁组成，面积约2.1×107km2。北冰洋面积约1.4×107km2，约占世界大洋总面积的3%，约占北极地区面积的67%，其大部分区域常年被海冰覆盖。白令海峡是北冰洋与太平洋连通的路径，在格陵兰岛西侧的加拿大北极群岛有许多与巴芬湾相连接的通道，北极西北航道就在此片岛屿中穿行。

本文使用2016年1月的FY-3 卫星的MWRI数据对北极地区的海冰分布进行分析研究。FY-3气象卫星能够全天候连续获取全球范围内的大气、地表及海洋表面多种特征参数，并广泛应用于三维大气探测、地表遥感特征提取等各方面研究。其中，MWRI是装载于FY-3卫星上的重要遥感仪器。MWRI可用于地表热辐射的监测，通过使用天线绕轴旋转形成圆锥形跨轨的扫描方式，在10.6～89 GHz间的多个固定频段上对来自大气及地表的双极化电磁辐射进行接收，并通过对数据的加工及分析，用于降雨反演、水汽监测、地表及海表上各类目标物的研究，且可为天气预报及环境监测提供参考数据。MWRI传感器具体参数见表1。

表1 FY-3 MWRI传感器若干参数

2 海冰分布探测方法

2.1 极化差

地物波谱特征是遥感反演的基础，不同地物的波谱特征曲线反映出了不同地物对电磁波反射的差异。其中，一年冰、多年冰和开阔水的光谱特征曲线如图1所示[22-23]。由亮温的定义可知，当物体的物理温度一定时，一定频率下亮温大小只与物体的辐射率有关。因为垂直极化和水平极化的电磁波在海冰表面或海水表面同时发出时，对同一物体的物理温度相同，因此极化差只受辐射率大小的影响。通过图1中的仿真结果可知，在微波成像仪(MWRI)的这几个频段中，19 GHz的一年冰和多年冰中的极化差值相差不多且与开阔水面相比差别最大，即可利用19 GHz的极化差P=19V-19H识别冰和水进而得到海冰的分布。即如果P>T，那么此像素为开阔水，否则为冰。

图1 海冰和海水发射率与频率的关系

2.2 Otsu算法

Otsu算法，又被称为大津法。该方法是根据图像数据值之间的方差选取阈值。其原理如下：

首先假设阈值具有L级数据的图像划分为两类：C0∈[0,T]及C1∈[T+1,L+1]图像各数据级对应的概率为

(1)

(2)

(3)

C0和C1类的均值分别为

(4)

(5)

(6)

海冰和海水的最佳分类阈值T*应使类间方差最大，即

(7)

2.3 极化差结合Otsu算法的海冰分布探测

利用Otsu算法对19 GHz极化差P进行处理，得到海冰和海水的分类阈值，进而得到北极海冰的分布信息。图2为海冰分布探测方法的流程图，基本步骤如下：

(1) 数据预处理：辐射定标、掩膜、异常数据点处理(当某个数据点与以其为中心的窗口内的数据点的中值相差太大时，该数据点被定义为异常数据点，用窗口内的数据点的中值代替)。

(2) 计算极化差：利用19 GHz的垂直极化和水平极化计算极化差P=19V-19H。

(3) 确定海冰和海水分类的最优阈值：选择样本点，基于Otsu算法对极化差P进行处理得到海冰和海水分类的最优阈值T。

(4) 得到海冰分布图：根据阈值T对极化差P图像进行分类，得到北极海冰分布图。

(5) 验证：所得结果与美国冰雪数据中心结果进行对比验证。

3 结果与验证

以2016年1月共31 d的FY-3 MWRI 19 GHz水平极化数据和垂直极化数据进行海冰分布探测，对19 GHz水平极化数据和垂直极化数据进行上述处理，得到19 GHz下的极化差异P，再选取样本区域，应用Otsu算法得到阈值为48.79 K。根据此阈值对图像进行分类得到2016年1月平均北极海冰分布图(如图3(a)所示)。获取同期1月NSIDC发布的北极海冰分布图(如图3(b)所示)。对比图3 可以看出，2个结果海冰外缘线分布状况基本一致；其中，MWRI数据反演得到的结果在高纬度地区的陆地边缘的海冰数量要多于NSIDC结果，由于陆地边缘在一定程度上会对海冰的判别造成误差，初步判断可能是美国冰雪数据中心为了消除陆地边缘造成的误差，对所有的陆地边缘进行了处理，但本文认为处于冬季的高纬度陆地边缘存在一些数量的海冰也是合理的，因此没有对其做进一步的处理；在鄂霍次克海区域(图幅的左上偏右区域)NSIDC结果海冰数量多于MWRI数据得到的结果，造成这种差异的原因可能如下：①由于微波数据是单轨成像，进行数据拼接，2种卫星过境时间不同，处于鄂霍茨克海区域轨道成像时间可能差异较大，因此造成这种误差；②该区域位于西北太平洋，与开阔水域链接畅通，可能受季风及洋流的影响较大，因此可能引起海冰移动，进而造成两者的误差；③2种微波数据的本身特性存在某种差异。当然，为了进一步地分析验证FY-3 MWRI数据及本文提出的算法应用到海冰识别方面的精确度，对于以上几个问题的深入研究，将会是笔者接下来努力的方向。

为了更为直观地对比2种海冰面积结果在1月份的变化趋势及差异，统计得到海冰范围统计表(见表2)及变化趋势图(如图4、图5所示)。

图2 19 GHz数据的海冰分布探测流程

图3 2016年1月北极海冰平均分布

图4 基于FY-3 2016年1月海冰范围变化趋势

图5 NSIDC 2016年1月海冰范围变化趋势

日期NSIDC海冰范围(×106)/km2FY-3海冰范围(×106)/km2FY-3与NSIDC海冰范围同比差值/(%)2016-01-0112.656212.18713.70612016-01-0212.746212.33093.25832016-01-0312.724312.37322.75922016-01-0412.820612.47002.73482016-01-0512.836812.39403.44952016-01-0613.011212.41044.61742016-01-0713.143112.43905.35572016-01-0813.015612.59453.23522016-01-0913.135012.57814.23962016-01-1013.348012.67205.06952016-01-1113.505012.65936.26152016-01-1213.392512.76714.66912016-01-1313.513712.95954.10112016-01-1413.455613.05732.95992016-01-1513.460613.07072.89612016-01-1613.521813.01603.74042016-01-1713.577513.02844.04392016-01-1813.663713.04074.55922016-01-1913.650613.08954.11032016-01-2013.625613.12073.70512016-01-2113.715613.19983.76052016-01-2213.704313.17463.86512016-01-2313.646213.20463.23572016-01-2413.748713.27573.44002016-01-2513.742513.27343.41322016-01-2613.866213.22854.59862016-01-2713.863113.21004.71122016-01-2813.861213.22144.61602016-01-2913.958713.25265.05842016-01-3013.928713.43962.79322016-01-3113.933713.33604.2892平均值13.493312.90564.3550极大值13.958713.43963.7185极小值12.656212.18713.7061极差1.30251.25253.8387

由表2及图4、图5可知，基于FY-3 MWRI数据利用Otsu算法结合19 GHz频段极化差得到的2016年1月海冰范围，与NSIDC提供的结果在日变化趋势上基本一致，其中MWRI得到的结果，以平均每天0.038 4×106km2的速度增加，NSIDC结果以平均每天0.041 9×106km2的速度增加，二者的增长速度仅相差0.003 5×106km2；NSIDC 1月海冰范围平均为13.493 3×106km2，MWRI数据结果1月海冰范围平均为12.905 6×106km2，二者同比相差4.355%；NSIDC 1月海冰范围极差为1.302 5×106km2，MWRI数据结果1月海冰范围极差1.252 5×106km2，二者同比相差3.838%；NSIDC 1月海冰范围单日极大值为13.958 7×106km2，出现日期为1月29日，MWRI 1月海冰范围单日极大值为13.439 6×106km2，出现日期为1月30日，造成二者的差异可能是由于2个卫星过境时间不同，成像不同步引起的，但综合来看二者的极大值比较接近，同比相差仅为3.718%；NSIDC 1月海冰范围单日极小值为12.656 2×106km2，出现日期为1月1日，MWRI 1月海冰范围单日极小值为12.187 1×106km2，出现日期为1月1日，二者较为同步，且同比相差仅为3.706%。因此，结合以上几个指标来看，基于本文算法得到的北极海冰分布范围与美国冰雪中心的海冰分布范围及日变化趋势基本一致，这在一定程度上说明本文提出的方法可行。

4 结 语

以2016年1月份共31 d FY-3 MWRI数据为例，基于海水与海冰在不同频段的极化差异的特性不同，选取19 GHz频段的数据，通过计算极化差，然后利用Otsu算法得到了海冰和海水的分类阈值，进而获取了北极的海冰分布信息。并与美国冰雪数据中心海冰分布的业务化产品进行对比分析，结果表明二者海冰外缘线分布基本一致，海冰面积也相差不大，海冰范围日变化趋势也较为同步。相对于其他光学及高频数据海冰分布的反演算法，本文提出的基于19 GHz频段极化差异特性的算法，更能用于大范围长时间序列的海冰监测与研究；同样在算法的实现过程中仅需确定一个分类阈值即可，具有较高的运算效率，对于即时性的大量数据的海冰信息发布与监测研究，具有明显的优势。因此表明本文提出的这种算法可行，尤其是在更高效率的运用国产卫星数据的层面具有一定价值与意义。