基于Freeman全極化分解的干雪識別指數模型構建

黃 林，李 暉，康 璇

（廈門理工學院計算機與信息工程學院，福建 廈門 361024）

積雪是冰凍圈內分布最為廣泛且最為活躍的要素，在水文循環、雪災預測等方面具有重要作用［1］。目前，使用光學遙感數據進行積雪研究易受多云多雨的影響［2］，而合成孔徑雷達（synthetic aperture radar，SAR）能夠穿透云層，具有全天候、全天時成像的特點，在山區積雪識別及分類應用中具有極大的潛力［3］。全極化SAR數據可以反映積雪的極化散射特性，區分地物的散射機制，為積雪識別提供了豐富的極化信息，有助于提高積雪識別的精度［4］。從SAR圖像中提取有效的極化特征對積雪特性研究和分析具有重要意義。同時，積雪極化特征的提取也是SAR 圖像識別的關鍵技術，對積雪的準確識別起著重要作用［5］。

國內外學者利用遙感數據對積雪識別方法進行了大量的研究。龐海洋等［6］基于Landsat8 OLI數據源，綜合考慮藍紫光、藍光、綠光和短波紅外波段，提出了增強型歸一化差值雪指數（enhanced normalized difference snow index，ENDSI）。高揚等［7］對青藏高原地區不同土地覆蓋類型下判識積雪的最優歸一化積雪指數（normalized difference snow index，NDSI）閾值進行了研究，表明不同土地覆蓋類型下的NDSI 閾值優化可以有效提高青藏高原積雪判別精度。馬騰耀等［8］利用GF-3 全極化SAR 數據，提出了一種基于特征優選的積雪識別方法。Varade等［9］提出了一種使用全極化SAR數據評估熵（H）和各向異性（A）繪制積雪的新指標，并將其稱為雷達雪分數（radar snow fraction，RSF）。王光遠等［10］基于高分六號（GF-6）多光譜相機（PMS）影像的目視解譯結果，結合光譜特征規律，探討了深度學習、決策樹等算法的積雪信息識別?？佃龋?1］對SAR 數據極化特征和散射特征進行組合，采用隨機森林算法對積雪進行了識別提取。目前，基于積雪指數的研究主要使用光學遙感數據，而使用SAR 數據構建指數的研究較少。其他領域的SAR 遙感指數研究主要是通過不同極化后向散射系數的數學組合［12-14］，或是基于SAR 數據的后向散射系數結合光學遙感指數的反演模型研究［15-19］，其中使用的方法主要還是機器學習方法［20-21］，但利用極化分解技術獲得不同散射機理的后向散射分量，進而構建物理意義更加明確的SAR 遙感指數的研究還較少［22］。因此，對于單一極化分解獲取的極化特征的信息提取還有很大的提升空間，如何利用較少的散射信息識別積雪范圍仍需進一步探索。

基于模型的分解直接與地物真實散射機制相關聯，具有一定的物理意義。Freeman 分解是第一個被提出的基于模型的分解［23］，相對于其他分解方法，Freeman分解更符合地物的散射機制，對地物散射特征描述得更加充分［24］?；诖?，本文對全極化SAR數據進行Freeman分解，通過極化特征值在積雪與非積雪面的特征，構建一種新的指數進行干雪識別。

1 研究區概況及數據源

1.1 研究區概況

以新疆天山中段北麓瑪納斯河流域（43.79°N～44.09°N，85.74°E～86.13°E）作為研究區，地表覆蓋物主要為草地，氣候屬于中溫帶大陸性干旱氣候，總體降水量較少，受水汽來源、地形和高程影響較大。根據自然高度帶，研究區可劃分為3 帶：50%以上區域海拔高度在0～1 300 m，為草原至半灌木過渡帶；南部約25%地區海拔在1 300～2 000 m，為地草甸草原帶；最南部少部分地區海拔在2 000～2 700 m，為云杉林帶［25］。2013 年12 月12 日至2013 年12 月16 日期間，在研究區內開展了為期5 d 與Radarast-2 衛星SAR 數據同步的地面積雪實地觀測，獲取研究區觀測點的雪深、海拔、坡度、氣溫、下墊面等信息，根據實地觀測，研究區積雪深度范圍介于0～20 cm。

1.2 數據源及預處理

獲取2013年12月13日的Radarsat-2全極化數據，包括HH、HV、VV、VH 4種極化方式，為單視復數數據（single look complex，SLC）C 波段產品，距離向和方位向分辨率分別為4.733 m 和4.799 m，中心入射角為43.45 °。結合實地觀測和天氣情況，對應積雪狀態為干雪。同時，選擇2013 年10 月6日的Landsat-8 OLI 數據（無雪期）、2013 年12 月14 日的GF-1 WFV 數據（干雪期）輔助參考。無雪期數據用于與積雪期數據對照判斷積雪區域及積雪區下墊面，干雪期數據用于同一時期SAR 數據的積雪判斷與精度評價參考。查詢歷史天氣記錄（https：//lishi.tianqi.com）可知，2013年12月13—14日溫度均在0 ℃以下，未有降水（降雪）記錄，說明干雪期SAR 數據和干雪期GF-1 WFV 數據在地表覆蓋類型上具有一致性。

為保證與光學遙感數據空間分辨率（20 m）相近，對Radarsat-2數據的距離向和方位向分別采用3×3 多視處理。對于全極化Radarsat-2 數據，相干斑濾波不僅需要考慮4 個通道（HH、HV、VH、VV）的濾波，還要考慮各個通道之間的相關性。因此，選用Refined Lee 濾波算法進行斑抑制，濾波窗口設置為5×5。最后，選擇影像對應的DEM數對影像進行地理編碼，本文選用NASA提供的90 m分辨率數據STRM4 DEM（https：//srtm.csi.cgiar.org/）對影像區域進行地理編碼操作，將影像從斜距幾何轉換為地圖投影坐標。

2 研究方法

2.1 Freeman極化分解、極化特征獲取與指數構建

在散射過程中，對SAR 目標極化效應進行定量描述的常用方法是運用極化散射矩陣，它包含目標的全部極化信息，其定義為

式（1）中：S表示單個像素散射特性；Shh、Svv表示同極化項；Shv、Svh表示交叉極化項。

極化協方差矩陣同極化散射矩陣一樣，包含雷達測量得到的全部目標極化信息。對于互易介質Shv= Svh，極化散射矩陣矢量可表示為以3×3 的視矩陣為例，協方差矩陣C 定義為

式（2）中：*表示共軛矩陣；T表示矩陣的轉置；|·|表示矩陣的模；〈·〉表示隨機散射介質在各向同性下的空間統計平均。

極化分解可以基于散射矩陣、相干矩陣或協方差矩陣實現［26］。全極化SAR 圖像通常使用9 個獨立參量表示接收到的數據，從而構成協方差矩陣或者相干矩陣，二者是等價的，都是非負定Hermitian矩陣［27］。Freeman分解以極化協方差矩陣為分解對象，將地物分為體散射、面散射和二面角散射3 類散射機制［28］。相對于其他分解方法，Freeman 分解更符合地物的散射機制，對地物散射特征描述得更加充分。3類散射方式的協方差矩陣分量為

式（3）～式（5）中：C1為面散射對應的協方差矩陣分量；C2為二面角散射對應的協方差矩陣分量；C3為體散射對應協方差矩陣分量；fs、fd、fv為各分量貢獻值；α、β為各種散射方式下極化散射矩陣分量中水平極化波與垂直極化波的比值?？偟臉O化協方差矩陣可表示為

由式（1）～（6）可解出未知量的值，并得到3類散射分量的功率為

式（7）～式（9）中：Ps、Pd、Pv分別為散射、二面角散射和體散射各分量的散射功率。

極化分解后獲得的極化參數，都可作為積雪識別過程中的一個特征［29］?；跇O化參數在覆蓋與非覆蓋區的不同散射機理，選取可以較大程度增加積雪與非積雪間差異的極化特征以構建歸一化指數。

2.2 Ostu閾值分割

Ostu 算法［30］是由日本學者大津提出的，也稱最大類間方差法，是圖像分割領域的經典算法。Ostu 算法的思想是：根據直方圖確定閾值t，將灰度圖像分成背景和前景2 部分，前景與背景的類間方差越大，說明兩者的差別越大［31］。若大小為m×n的圖像，設圖像中像素的灰度值小于閾值t的像素個數為n0，大于閾值t的像素個數為n1，則有

由式（10）～式（13），類間方差g（t）可表示為

式（10）～式（14）中：ω0為背景像素點所占比例；ω1為前景像素點所占比例；μ 為圖像灰度值；μ0和μ1分別為背景和前景的均值。通過不斷對閾值t進行遍歷，當g（t）取得最大時，類間方差最大，即錯分概率最小，此時t即為最佳閾值。

本文中，通過Ostu 自動閾值分割算法確定閾值t，將小于閾值的像素點記為非積雪點，將大于閾值的像素點記為積雪點，對得到的積雪指數模型進行分割，以達到積雪識別的目的。

2.3 監督分類與精度評價

在遙感圖像監督分類中，距離法和最大似然法是較常用的2種方法。本文選用最小距離、馬式距離和最大似然分類方法對研究區進行分類，獲得積雪與非積雪。

在閾值分割和監督分類之后，本文以同一時期的GF-1 WFV 數據作為參考，隨機選擇積雪區域與非積雪區域各10 個進行精度評價，共130 541 個像元，其中積雪像元70 952 個，非積雪像元59 589個，對識別結果建立混淆矩陣。以總體精度（overall accuracy，OA）、Kappa系數（Kappa coefficient）、用戶精度（user’s accuracy，UA）和制圖精度（producer’s accuracy，PA）對識別結果進行精度評價?？傮w精度可以分析整體精度效果，Kappa系數可以確定積雪識別結果與真值的一致性，用戶精度和制圖精度可以體現誤分率和漏分率。

3 實驗結果與分析

3.1 Freeman極化分解特征分析

Freeman 極化分解后的極化特征為Freeman_Dbl、Freeman_Vol和Freeman_Odd，分別表示二面角散射、體散射和面散射，圖1 為研究區Radarsat-2 全極化數據Freeman極化分解后結果。

由于研究區積雪深度范圍介于0～20 cm，相對較淺，C波段的電磁波可以在雪表面發生折射，并穿透干雪。當有積雪覆蓋時，在干雪條件下，主要散射為雪-地界面的面散射，雪層體散射比總體面散射少［32］，結合光學遙感數據可知，積雪覆蓋區表現為藍紫色。無雪覆蓋區主要出現在高覆蓋度草地和林地，主要散射為植被體散射，圖像呈綠色。

3.2 基于Freeman極化分解的FSCI_Vol_Dbl構建

在Freeman 極化分解的圖像中，選擇干雪期積雪覆蓋區域與非積雪覆蓋區域，分別統計各散射特征的概率分布密度，結果如圖2 所示。由圖2（a）可知，在積雪覆蓋區，3 個散射分量的分布由小到大分別是二面角散射、體散射和面散射；二面角散射介于-29～-21 dB，主要集中在-25 dB；體散射介于-22～-14 dB，主要集中在-19 dB；面散射介于-20～-5 dB，分布特征呈高原形，數值較不集中。因此可知，二面角散射分布較為集中，且與體散射、面散射均存在較高的可分離性；體散射普遍小于面散射，且由于積雪覆蓋較淺，體散射和面散射無法直接分開，有較多區域混合在一起。由圖2（b）可知，2個散射分量的分布由小到大的分布是二面角散射、面散射和體散射；二面角散射介于-29～-15 dB，主要集中在-24 dB；體散射主要集中在-8 dB，且普遍比面散射大；面散射分布進一步擴大，分布范圍介于-20～-1 dB之間。

圖2 干雪期積雪區與非積雪區散射特征概率分布密度Fig.2 Sampling histogram of snow covered area and non snow covered area in dry snow period

對比圖2（a）和圖2（b）發現，積雪區體散射小于非積雪區的體散射，面散射在積雪區與非積雪區均呈高原型；由數值大小坐落范圍分析可知，二面角散射與體散射、面散射均存在較好的分離度，但二面角散射與體散射均分布較為集中。因此，基于二面角散射與體散射的規律，將二面角散射與體散射構建差值指數為

式（15）中：FSCI_Vol_Dbl表示由體散射、二面角散射計算的積雪覆蓋度結果；Freeman_Vol為Freeman 的體散射；Freeman_Dbl為Freeman 的二面角散射。FSCI_Vol_Dbl通過積雪區與非積雪區體散射與面散射的特點，計算出的范圍為-1.0～0，一般越接近于0，亮度值越高，積雪覆蓋度越高。

同時，為了進一步驗證選擇二面角散射與體散射優于其他分量構建指數，在干雪期，本文從其他分量呈現出的規律，構建了另外2種積雪范圍識別指數。

1）基于面散射、體散射在積雪與非積雪區的規律，選擇面散射、體散射分量做差值，構建指數為

式（16）中：FSCI_Vol_Odd表示由面散射、體散射計算的積雪覆蓋度結果；Freeman_Vol為Freeman 的體散射；Freeman_Odd為Freeman的面散射。

2）基于二面角散射與面散射的規律，選擇二面角散射、面散射做差值，構建指數為

式（17）中：FSCI_Odd_Dbl表示由面散射、二面角散射計算的積雪覆蓋度結果；Freeman_Odd為Freeman 的面散射；Freeman_Dbl為Freeman的二面角散射。

3.3 利用FSCI_Vol_DBl進行積雪識別的分析與比較

3.3.1 3種指數積雪識別對比

將3 種指數得到的結果與GF-1 WFV 影像、NDSI 進行比較，基于Freeman 分解特征分量的指數構建結果與對比如圖3 所示。NDSI 利用積雪在光譜的可見光部分比短紅外波段更高反射率的特點，使用綠色和短紅外波段進行計算，范圍為-1.0～1.0，值越接近于1，亮度值越高，積雪覆蓋度越高。由于同時期獲取的GF-1 WFV 數據只有4 個波段，缺少計算NDSI 所需的短紅外（SWIR）波段，且該時期GF-1 WFV 數據北部被大量云層覆蓋，因此，本文同時利用同時期的GF-1 WFV 影像（見圖3（a））、2013 年12 月9 日的Landsat-8 OLI 數據計算的NDSI（見圖3（b））與FSCI_Vol_Dbl結果進行比較。由于Landsat-8 OLI得到的NDSI比GF-1 WFV數據提前5 d，因此會存在一定誤差。

圖3 基于Freemman分解特征分量的指數構建結果與對比Fig.3 Index construction results and comparison based on Freeman decomposition feature component

圖3（c）為FSCI_Vol_Dbl的結果，可見，利用二面角散射與體散射構建的指數在積雪覆蓋區的灰度值遠遠高于非積雪覆蓋區，可以通過亮暗差異反映積雪覆蓋；圖3（d）為FSCI_Vol_Odd的結果，可見，構建的指數在積雪覆蓋區與無雪覆蓋區亮度差異并不明顯，在積雪區與非積雪區無法通過灰度值大小直觀反映出積雪覆蓋；圖3（e）為FSCI_Odd_Dbl的結果，可見，展現出的積雪覆蓋效果較不理想，主要原因是面散射的概率分布范圍較寬，密度較不集中，無法較統一地將積雪與非積雪的差異增大。

由圖3的結果對比可見，利用二面角散射與體散射構建的指數是積雪覆蓋度的最佳指示因子，增大了積雪與非積雪的差異，突出了積雪范圍。

FSCI_Vol_Dbl與Landsat-8 OLI數據計算的NDSI相比，兩者在研究區內積雪與非積雪區輪廓較為明顯，二者顏色分布較為相似。與GF-1 WFV數據相比，FSCI_Vol_Dbl中積雪區/非積雪區體現的差異比NDSI大，NDSI數據中整個研究區內亮度值均較高，存在較大的積雪高估范圍。FSCI_Vol_Dbl利用Freeman極化分解后的特征分量，通過非線性變換，增強了FSCI_Vol_Dbl低值部分，并抑制了高值部分，使得FSCI_Vol_Dbl數值較一地達到飽和狀態，從而對高覆蓋度積雪區敏感度降低。因二面角散射與體散射均為負值，計算出FSCI_Vol_Dbl的值為負數，在-1～0，避免了數據過大或過小給數據分析造成影響。此外，FSCI_Vol_Dbl還能夠消除部分山區積雪識別易受云層覆蓋限制的影響。由于Freeman 極化分解獲取的特征分量為體與面的散射，極易受地形因素影響，計算得出FSCI_Vol_Dbl的值接近-1，一般表示為地勢起伏較大的山地，地形復雜，使得出現小于-1 的異常值。FSCI_Vol_Dbl的值接近0，一般表示為該地區為積雪覆蓋。

確定積雪與非積雪間的閾值，通過數值的比較，即可判斷積雪與非積雪，快速提取積雪空間分布信息。本文利用Ostu 閾值分割算法對FSCI_Vol_Dbl得到的閾值結果為-0.246。以-0.246 為閾值確定的積雪與非積雪區，具體如圖4 所示。由圖4 可見，FSCI_Vol_Dbl識別結果與GF-1 WFV影像中積雪范圍相似，可以識別積雪。

圖4 FSCI_Vol_Dbl閾值分割結果Fig.4 Threshold division result of FSCI_Vol_Dbl

3.3.2 與監督分類結果比較

圖5 為利用FSCI_Vol_Dbl進行監督分類的干雪識別結果。其中，圖5（a）、圖5（b）、圖5（c）分別為FSCI_Vol_Dbl以最小距離、馬氏距離和最大似然法得到的積雪識別結果。

圖5 監督分類積雪識別結果Fig.5 Snow cover recognize results of supervised classification

由圖5可見：整體上，3種監督分類方法都完成了積雪范圍的識別，且范圍與光學遙感影像相近。其中，在光學數據北部被云覆蓋的地區（圖3a），SAR 數據可以減輕云層的影響，達到了識別效果。從分類方法上分析，FSCI_Vol_Dbl的最小距離分類對積雪識別的范圍斑塊化較少，結果較平滑，但存在部分非積雪覆蓋地區誤分為積雪地區，易誤判的地區集中在山區的背陰面及地形起伏較大的地區；馬氏距離分類的結果較零碎，積雪識別的范圍在所有方法中最小，積雪地區存在較多非積雪斑塊；最大似然分類的結果與馬氏距離相近，但識別的積雪零碎斑塊比馬氏距離少。研究區南部地區地勢較高，地形復雜，該地區的識別結果在最小距離分類的結果中與馬氏距離、最大似然分類結果差異較大。

表1 為混淆矩陣精度評價結果，精度評價中將FSCI_Vol_Dbl閾值分割的結果（圖4）與所有監督分類結果（圖5）比較。從OA 和Kappa 系數上看，利用本文提出的FSCI_Vol_Dbl進行積雪識別的精度介于80%～90%，Kappa 系數介于0.60～0.80，呈現出高度的一致性。通過FSCI_Vol_Dbl閾值分割的精度比最小距離監督分類精度高3.73%，比馬氏距離監督分類精度高0.9%，與最大似然分類結果相當。從UA和制圖精度（PA）分析，FSCI_Vol_Dbl閾值分割的結果UA 高于PA，可見，通過FSCI_Vol_Dbl指數閾值分割識別積雪對積雪像元的判斷高于非積雪像元。

表1 混淆矩陣精度評價結果Table 1 Precision evaluation of confusion matrix

4 結論

對新疆瑪納斯河流域全極化SAR 數據通過極化分解技術，獲得不同散射機理的后向散射分量，根據概率分布密度在干雪期積雪覆蓋面與非積雪覆蓋面極化特征的規律，構建了基于極化分解的干雪識別指數模型。結果表明，基于極化分解的極化特征分量受到積雪覆蓋影響，通過增大不同散射機理的數值差構建指數模型，利用閾值可劃分積雪/非積雪范圍，總體精度為85.83%，Kappa 系數為0.716，具有一定的準確性，在進行山地積雪識別時具有一定的參考意義。

同時，由于受山地地形環境復雜、植被類型多樣、其他大氣狀況、積雪本身的晶體狀態及下墊面等因素的影響，SAR指數的適用范圍往往具有一定的地域性和時效性，本文指數是否適合其他區域的研究還需要進一步驗證。后續將綜合地形等因素，充分利用獲取的散射特征，構建從機理層次與積雪緊密聯系的指數進行研究。