赵朝方， 徐 锐， 赵 可

(1. 中国海洋大学信息科学与工程学院，山东 青岛 266100； 2. 青岛海洋科学与技术国家试点实验室，区域海洋动力学与数值模拟功能实验室，山东 青岛 266011)

海冰是一种敏感的气候指标，限制海洋和大气之间的能量、动量和物质的交换。同时，海冰的存在影响海表反照率，减少地球表面吸收的太阳辐射量。海冰在形成过程中析出盐可以导致冰下海水密度增加，有助于垂直方向的海洋的温盐环流，而海冰融化会增加上层海洋的稳定性。此外，海冰直接影响高纬度地区船舶航行和近海作业。因此，在高纬或极地海域，准确确定海冰的范围，对航行安全十分重要[1]。同时，研究极地海冰的时空特性对于监测全球气候变化、天气预报、渔业捕捞和海洋运输具有重要意义[2]。

早期的海冰观测方式包括测量船、飞机探测等，这两种观测方式不仅效率低、成本高，而且对于位置偏远且环境恶劣的极地海冰，更是难以直接观测。随着遥感技术的发展，卫星遥感已经成为研究极地海冰的有效手段。用于海冰检测的遥感传感器主要包括可见、红外传感器和微波传感器。可见光与近红外波段传感器包括NOAA/AVHRR传感器、EOS/MODIS传感器、以及近年来发射的高分辨率可见和红外传感器等，它们所获取的各波段遥感资料被广泛应用于海冰范围及海冰密集度的监测研究[3]。但是可见和红外传感器在多云和极夜等条件下容易受到云或雾等天气的影响，造成遥感图像模糊影响观测结果。微波可以穿透云层，具备全天时、全天候的工作能力，因此微波传感器成为观测海冰的重要工具。目前有两种微波传感器用于海冰观测，被动微波传感器和主动微波传感器。被动微波传感器主要是指多波段微波辐射计如雨云系列卫星Nimbus-5搭载的电子扫描微波辐射计(Electronically Scanning Microwave Radiometer,ESMR)，主要通过测量亮温获取海冰边缘等参数。搭载在Nimbus-7卫星的扫描多通道微波辐射仪(Scanning Multichannal Microwave Radiometer,SMMR)，将SMMR海冰观测与飞机观测结果进行比较，SMMR对总冰密集度的估计精度为±3%，多年冰密集度的估计精度为+10%[4]。1987年开始的美国国防气象卫星计划(Defence Meteorological Satellite Program,DMSP)系列卫星先后搭载了专用微波成像仪和专用微波成像探测仪(Special Sensor Microwave/Imager and Special Sensor Microwave Imager/Sounder,SSM/I和SSMIS)，成为海冰数据的主要来源[5]。2002年EOS-Aqua卫星搭载改进型多通道扫描辐射计(Advanced Microwave Scanning Radiometer-Earth Observing System,AMSR-E)，比SSM/I具有更高的空间分辨率和更大的光谱范围，能够更好地评估不同类型海冰的空间分布[6]。

观测海冰的另一种微波传感器是主动微波传感器，主要包括合成孔径雷达(Synthetic Aperture Radar,SAR)和散射计等。1978年，海洋卫星Seasat-A成为第一颗提供高分辨率海冰SAR影像的卫星，SAR传感器以其高分辨率成像能力、不受多云天气影响，可以获取海冰边缘线和冰块缝隙等信息，在极地和高纬度海冰检测中发挥着重要作用。但SAR的覆盖范围小，资料成本高，其时空覆盖率难以满足极地海冰和中大尺度海洋现象的研究[7]。1990年代以来，微波散射计开始用于海冰探测，包括C波段的ERS-1、ERS-2、ASCAT以及Ku波段的NSCAT、QUIKSCAT、HY-2A/SCAT等。微波散射计是非成像传感器，其设计初衷是通过测量后向散射系数(σ0)反演风场，但是近年来在海冰方面的探测也逐渐成熟起来。

目前来看，针对HY-2A/SCAT数据进行海冰检测算法还较少，缺少更加全面的比较。贝叶斯海冰检测算法虽然已经被很多散射计所应用，却还未用HY-2A/SCAT数据进行过实验；此外，李明明等[2]已经用FLDA算法对HY-2A/SCAT数据进行了海冰检测，模式识别中广泛应用的支持向量机和神经网络算法在很多方面都比FLDA表现出独特的优越性，因此将这两种算法加入HY-2A/SCAT数据的海冰检测研究具有重要意义。此外，参考文献[14]在参数选择方面的启发，本文在输入参数方面基于主成分分析(Principal Component Analysis,PCA)对特征参量进行选择，在此基础上利用BP神经网络检测算法来进行海冰检测。本文使用贝叶斯海冰检测算法、线性判别算法、支持向量机算法和基于PCA的BP神经网络算法对HY-2A/SCAT进行极地海冰检测，分析比较不同算法的结果，并在此基础上分析极地海冰的季节变化。

1 数据介绍

1.1 主要数据

1.2 辅助数据

本文所用辅助数据是搭载于DMSP的24通道无源微波辐射计SSMIS海冰密集度数据，用来与HY-2A/SCAT海冰检测结果进行对比分析。SSMIS频率范围为19～183 GHz，可获得每日海冰密集度，是目前公认的海冰数据集，该海冰密集度数据可以在美国数据冰雪中心获取[17]。

图1 散射计观测几何[16]

SSMIS极地网格海冰密集度数据的空间分辨率为25 km，数据来源：ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051gsfc_NASAteam_seaice/。

1.3 数据时空范围

本文使用HY-2A/SCAT全球数据对南北极地区的海冰进行探测。作为标准数据的SSMIS极地网格数据采用NASA 和 NSIDC提出的极地椭球投影方法(正轴等角割方位投影)[18]，为了便于观察和比较，对HY-2A/SCAT南北极地区数据采用相同的方法进行投影，投影前后分辨率均为25 km。投影后北极和南极平面栅格大小分别为448×304和332×316。投影后北半球经纬度范围：30.98°N～90°N，180°W～180°E。投影后南半球经纬度范围：39.23°S～90°S，180°W～180°E。

本文所研究的时间范围：2013年9月—2014年8月。图2为极地平面地图网格，其中lon表示经度(Longitude),lat表示纬度(Latitude)。

2 研究内容

2.1 贝叶斯海冰检测算法

2.1.1 贝叶斯海冰检测算法原理 贝叶斯海冰检测算法在海水和海冰地球物理模式函数(GMF)的基础上，通过统计到海水和海冰两个GMF对应的后向散射系数(σ0)分布的概率距离,利用贝叶斯公式来实现海冰检测。后向散射系数的大小主要取决于海面粗糙度，在无冰的海面上海面粗糙度主要由风应力决定[10]，因此海水GMF即海面风场的GMF。

根据贝叶斯公式，海冰的条件概率可以写成：

(1)

式中：p(ice|σ0)代表待识别单元分类为海冰的后验概率，p0(ice)与p0(wind)的初始值设定为0.5,条件概率p(σ0|ice)与p(σ0|wind)可用下式表示：

p(σ0|wind)=p(MLEwind)；

(2)

p(σ0|ice)=p(MLEice)。

(3)

MLEwind和MLEice分别为特征值到海面风场GMF(相当于没有结冰的海面)和冰GMF对应的σ0分布函数的最大似然估计(Maximum Likelihood Estimators, MLE)，表示特征值相对于两个模型σ0分布的最小平方距离。参考Seawinds海面风场反演时MLE的定义[9]，MLEwind和MLEice分别可以表示为：

(方框内为研究区域。Rectangle region is the research area.) 图2 极地平面地图网格

(4)

(5)

MLEwind和MLEice函数的期望分布分别可以近似为具有N-2和N-1个独立自由度的卡方分布[20],概率分布表示为：

(6)

(7)

为了充分利用历史先验信息，每处理一轨数据都会更新p0(ice)与p0(wind)值，根据Seawinds进行贝叶斯海冰检测时的p0(ice)与p0(wind)的更新原则[9]，对HY-2A/SCAT进行处理时也采用相同的阈值，即：

(8)

2.1.2 贝叶斯海冰检测算法参数获取

2.1.2.1 海面风场GMF与海冰GMFσ0分布函数参数 本文使用的海面风场GMF 为NSCAT-4 GMF，该GMF是在NSCAT 双极化后向散射测量与ECMWF(European Centre For Medium Range Weather Forecasts)模式基础上建立的[20]，将后向散射系数(σ0)表达为风速、风向、入射角以及极化方式的函数[21]。

σV,dB=A×σH,dB+B。

(9)

计算海冰GMF σ0分布函数的根本是确定线性方程的斜率A和截距B。海冰GMF σ0分布线性方程的斜率受混合冰水以及冰架等很多因素的影响，为了保证海冰GMF σ0分布线性方程的一致性，应选取最具代表性时间段的纯海冰拟合这条直线的斜率[9]。

图与之差的绝对值分布(单位：dB)

图4是不同海冰密集度下北极每日线性模式方程斜率。参考Seawinds海冰GMFσ0分布函数的参数计算过程，结合斜率曲线，选择1—3月份，60°N以上范围内的海冰的平均斜率为最终的海冰GMFσ0分布函数斜率。

最终确定的北极海冰GMFσ0分布函数的斜率为1.095 8，截距为-0.192 2。

得到海冰GMFσ0分布拟合函数的斜率和截距后，可以最终确定NSCAT-4 GMF 与海冰GMF在原始后向散射测量空间中的分布。图5为两个模型在HH和VV两个极化后向散射测量空间中的σ0分布，虚线框选范围为海面风场GMFσ0分布，直线为HY-2A/SCAT冰GMFσ0分布的线性拟合。

(北极参考斜率用红色虚线标出。The reference Arctic winter slope is highlighted in red dashed line.)

图4 不同海冰密集度下北极每日海冰 GMFσ0分布线性函数拟合斜率

Fig.4 Arctic daily fitting slopes to sea ice GMFσ0distribution linear model of three different sea ice concentration

图5 NSCAT-4GMF与HY-2A/SCAT海冰 GMF在垂直和水平极化测量空间中的σ0 分布

图6 后向散射系数观测值到海冰GMF σ0 分布线性函数的距离(标准差)示意图

2.2 线性判别(FLDA)海冰检测算法

假设有n组样本数据，每一组样本数据可以表示为四维参数矢量xi，由于海冰检测只有冰和水两个类别，因此只需要把样本数据投影到一条直线上。对任意一组样本矢量,它在直线w上的投影可以表示为wTxi,设两个类别的中心点μ0,μ1， 在直线w的投影分别为wTμ0和wTμ1。根据线性判别原理，需要让冰水两个类别中心的距离尽可能大，亦即最大化‖wTμ0-wTμ1‖2,同时要求同一类别数据的投影点尽可能的接近，相当于要求协方差wT∑0w和wT∑1w尽可能小，即最小化wT∑0w+wT∑1w。综上所述，优化目标函数为：

(10)

类间离散矩阵表示为

SB=(μ0-μ1)(μ0-μ1)T。

(11)

类内离散矩阵表示为

图7 海冰后向散射到海冰GMFσ0 分布线性函数的距离频率直方图分布

SW=Σ0+Σ1=

(12)

因此优化目标函数可以重写为

(13)

采用拉格朗日乘子法解出上式的最佳投影方向

(14)

则判别函数可以写为

y=w*Tx。

(15)

然后确定分类阈值w0,本文采用下式来获得w0:

(16)

根据这个原则将所有样本数据投影至最优一维矢量空间后，比较投影后待识别单元的y值与阈值w0的大小,判断该待识别单元属于海冰还是开阔海水。

2.3 支持向量机(SVM)海冰检测算法

由于线性判别分析(FLDA)只能处理线性可分问题，但是数据样本通常并不一定线性可分。支持向量机(Support Vector Machine,SVM)算法通过引入核函数解决非线性可分问题。

wT×φ(xi)+b=0 (i=1,…,n)。

(17)

(18)

s.t.((wT*φ(xi))+b)≥1 (i=1,…,n)。

(19)

式中：y是类别标签；yi是y的第i个分量。本文仅有海冰和海水两个类别，即y设置为+1和-1。这个公式可以使用拉格朗日乘数法来求解。即先对每条约束添加拉格朗日乘子βi≥0，则该问题的拉格朗日函数可写为：

(20)

用序列最小优化算法(SMO)可得w的最优解w*：

(21)

得到w*后可代入式(17)得到的b最优解b*。则对于任意测试样本点m，其决策函数可以写为：

f(m)=

(22)

式中：m表示投影后测试样本点参数；φT(xi)φ(m)即为投影后测试样本点与训练样本的内积，由于直接计算这个内积比较困难，所以用核函数代替内积：

(23)

本文使用高斯核函数，可以表示为：

(24)

式中,σ是达到率，即函数值跌落到0的速度参数。

2.4 基于主成分分析的BP神经网络海冰检测算法

主成分分析是通过正交变换将一组可能存在相关性的变量转换为线性不相关的变量。主成分分析可以减少特征数量，有助于降低数据变量的复杂性。BP(Back Propagation)神经网络是采用误差逆向传播算法训练的多层前馈神经网络[22]，利用BP神经网络对目标进行分类，若输入太多的样本特征，则会降低网络的训练速度与效率，严重时会导致网络不收敛。所以结合PCA对原始数据进行主成分提取可以提高神经网络训练速度与效率。

3 结果分析

为了对比不同算法检测海冰的能力，本文从冰水误判、海冰边界、海冰面积变化趋势三个方面对四种算法的海冰检测结果进行分析。

3.1 海冰检测结果误判率对比

图8为利用四种算法得到的2013年9月16日北极海冰范围检测结果，以及与同时期SSMIS 15%海冰密集度的海冰范围(见图8(e))进行比较。从图中可以看出，FLDA和SVM算法得到的海冰检测结果在区域C处都存在大量海冰误判，贝叶斯算法与基于PCA的BP神经网络算法在这个位置的误判相对较少。四种方法的误判率(噪声点与全部测量点的比)进行统计见表1。

从表1可知，贝叶斯算法的总误判率最小，其次是基于PCA的BP神经网络算法，SVM算法和FLDA算法的误判率较大。Rivas[20]利用贝叶斯算法对Seawinds数据进行冰水检测时，发现高风速会影响检测结果，因此利用数值预报NWP(Numerical Weather Prediction)风速进行校正；此外，应用RL-N算法海冰检测也发现，高风速影响检测效果，尤其是在格陵兰岛和威德尔海域[9]。其原因在于高风速的海面，海冰和开阔海水的后向散射系数非常相近，导致海冰和开阔海水的特征区分不明显。

(a 贝叶斯算法 Bayesian； b 基于PCA的BP神经网络算法 BP neural network based on PCA； c 线性判别算法 FLDA； d 支持向量机算法 SVM； e SSMIS海冰范围 SSMIS sea ice extent。区域C处表示四种算法的误判差异。Area C indicates misjudgment difference among four algorithms.)

图8 不同算法得到的2013年9月16日 北极海冰范围和SSMIS结果

Note： ① Water to ice; ② Ice to water; ③ All misjudgment points; ④ Bayesian; ⑤ FLDA; ⑥ SVM; ⑦ BP neural network based on PCA

为了验证区域C处产生误判的原因，对2013年9月16日区域C内后向散射系数分布(见图9)和ECMWF风速(见图10)进行分析。

从图中可以看出，误判位置的后向散射系数比周围开阔海水区高的多，可以达到-10 dB，风速也都在15 m/s以上。因此，该位置产生的过多误判与风速有着很大联系。

图9 北极地区分布(2013年9月16日)

图10 北极地区ECMWF风速(m/s)(2013年9月16日)

贝叶斯方法由于结合了先验概率信息，其冰水检测不完全依赖当前观测的后向散射系数，结合海面风场和海冰GMF并加入了NWP风速，可以校正风速过高引起的误差，因此误判相对较低，优于FLDA和SVM算法；基于PCA的BP神经网络算法，在主成分分析时仅保留最主要的特征信息，并且BP神经网络通过反复的后向反馈，理论上可以实现任何非线性映射的功能，所以误判率也相对较低，但仔细分析冰水的边界，发现这种算法得到的海冰边界比较粗糙。

3.2 海冰边界对比

SSMIS海冰密集度数据与HY-2A/SCAT海冰范围采用相同的极地投影方式，并且具有相同的分辨率，可以直接进行海冰边界的比较。将四种算法的海冰边界与SSMIS海冰边界进行对比后得到，四种算法的海冰边界均介于SSMIS 0%和30%海冰密集度海冰边界之间。

为了对HY-2A/SCAT南北极海冰边界进行具体分析，本文只列出了贝叶斯算法和SVM算法的海冰边界与SSMIS 0%、30%海冰密集度海冰范围边界的对比图(见图11)。为了清楚地展示比较结果，海冰轮廓叠加在SSMIS 0%密集度图像上。

(a 贝叶斯算法(南极) Bayesian(Antarctic); b 贝叶斯算法(北极) Bayesian(Arctic); c 支持向量机(南极) SVM(Antarctic); d 支持向量机(北极) SVM(Arctic)。绿色线条表示SSMIS 0%密集度边界；黄色线条表示SSMIS 30%密集度边界；红色线条代表贝叶斯和SVM得到海冰边界。Green line indicates SSMIS 0% concentration edge; yellow line indicates SSMIS 30% concentration edge; red line indicates Bayesian and SVM edge.)

图11 基于贝叶斯和SVM算法得到的 2013年9月20日极地海冰边界

Fig.11 Polar sea ice edge based on Bayesian and SVM on Sept.20,2013

从图11可以看出，两种算法获得的海冰边界都在海冰密集度0%～30%之间，由此说明这两种算法都可以用于海冰边界的检测。

3.3 南北极海冰范围季节趋势分析

考虑到PCA的BP神经网络算法运行效率，以及FLDA误判较多，本文选取贝叶斯算法和SVM算法海冰检测结果，对HY-2A/SCAT 2013年9月到2014年8月份极地海冰变化趋势进行分析(见图12)。其中，对SVM算法结果进行了图像腐蚀膨胀技术处理，贝叶斯海冰检测算法中α取3，将SSMIS 0%、15%、30%密集度以上的海冰面积作为比较标准。

从整体来看，贝叶斯算法与SVM算法海冰面积变化趋势图都能很好的展现海冰面积的季节性变化。在海冰结冰期(冬季)，与SSMIS 15%、30%密集度海冰面积相比，贝叶斯和SVM算法得到的海冰面积变化都与其比较一致，且都介于SSMIS 0%、30%密集度海冰面积曲线之间；而在海冰融冰期(夏季)，贝叶斯算法和SVM算法探测的海冰面积都要比SSMIS 30%海冰密集度以上的海冰面积大，且贝叶斯算法探测的海冰面积超过SSMIS 30%海冰密集度以上的海冰面积的程度更大，而SVM算法结果与SSMIS 15%海冰密集度以上的海冰面积比较一致,该结果与之前研究[9，23]基本相符。

图12 2013年9月1日～2014年8月31日 基于HY-2A/SCAT的极地日海冰范围

4 结语

本文针对HY-2A/SCAT可以进行海冰检测的有用参数，详细研究了贝叶斯、FLDA、SVM以及基于PCA的BP神经网络的海冰检测算法，分析HY-2A/SCAT在极地海冰检测中的应用效果，以及在冰水误判、海冰面积变化趋势等方面的表现。

从海冰检测结果来看，四种方法都能基于HY-2A/SCAT进行海冰检测。初步结果表明，高风速引起的误判对贝叶斯算法和基于PCA的BP神经网络算法影响较小,FLDA以及SVM算法受高风速影响较大。其中针对贝叶斯算法，本文基于HY-2A/SCAT数据获取了适合HY-2A/SCAT的海冰GMF σ0分布拟合参数及其方差，实现了HY-2A/SCAT贝叶斯海冰检测算法。

将贝叶斯和SVM算法得到的海冰边界与SSMIS 0%、30%海冰密集度的海冰边界相比，结果显示这两种算法的海冰边界都介于SSMIS 0%、30%海冰密集度的海冰边界之间。

对贝叶斯和SVM算法得到的海冰面积季节变化趋势与SSMIS 0%、15%、30%海冰密集度的海冰面积变化相比较，它们在结冰期(冬季)吻合度较高；而在融冰期(夏季)，基于散射计的两种算法得到的海冰面积比SSMIS 30%海冰密集度以上的海冰面积要大，但SVM算法探测的海冰面积与SSMIS 15%海冰密集度以上的海冰面积比较一致。

综上所述，四种算法均适用于HY-2A/SCAT的海冰检测研究，证明了HY-2A/SCAT具有良好的海冰检测能力，不失为一种海冰探测的有效手段。在未来的工作中，将进一步探索HY-2A/SCAT在海冰分类方面的研究，了解不同算法的物理机制，以及寻找更多更适合的海冰检测特征，为极地海冰监测与全球变化研究服务。