基于HY-2A/SCAT數據極地海冰檢測方法研究*

趙朝方， 徐 銳， 趙 可

(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在海冰分類方面的研究，了解不同算法的物理機制，以及尋找更多更適合的海冰檢測特征，為極地海冰監測與全球變化研究服務。