干雪深度反演的同極化相位差模型

宋依娜，肖鵬峰,2,3，張學良，卓 越，馬 威

1.南京大學地理與海洋科學學院，江蘇 南京 210023；2.自然資源部國土衛星遙感應用重點實驗室，江蘇 南京 210023；3.江蘇省地理信息技術重點實驗室，江蘇 南京 210023

積雪是冰凍圈和水圈重要的組成部分，在全球變化研究中具有重要地位[1]。一方面，積雪作為一種特殊的地表覆蓋物，通過反射大部分太陽輻射以及在消融過程中吸收大量能量，對全球氣候系統產生了顯著的降溫作用[2-3]。同時，積雪融水作為全球水資源的重要組成部分，為全球1/6人口提供了水源[4]，并進一步產生了水力發電、農業灌溉等其他用途[5]。另一方面，積雪具有時空不穩定性，過量或非適時的積雪也會引發自然災害，給人類生命和財產安全帶來巨大的威脅和損失[6]。而積雪深度作為積雪的重要因子，表征局部氣候環境特征與水資源條件，是水文預測、氣候模擬、雪災監測與評價等模型的重要輸入參數[7]。獲取高精度雪深空間分布信息可為水資源管理、氣候變化研究和防災減災工程等提供科學支撐。

合成孔徑雷達(synthetic aperture radar,SAR)因其有較高的空間分辨率[8]且對積雪特性有著較高的敏感度等優勢[9-10]，在流域尺度的雪深反演中具有較強的適用性[11-12]。極化SAR通過回波獲取關于積雪的散射特征[13]，被廣泛應用于積雪相關研究，如干濕雪的分類、雪濕度和雪密度等參數的反演[14-16]，但是目前的研究對極化SAR的極化相位信息的分析較少[17]。極化SAR計算同一天線、不同極化方式的極化相干，可獲取極化相干系數和極化相位信息，包括同極化相位差(co-polarized phase difference,CPD)。積雪為各向異性介質，導致HH和VV極化信號在雪層中具有不同的傳播速度，因此產生CPD。文獻[18]根據此原理建立積雪CPD模型，定量分析了CPD與積雪各向異性結構、雪深的關系。但是，根據CPD模型反演干雪深度需要冰粒的各向異性結構(以冰粒軸間比表示)和積雪密度兩個參數，其中冰粒軸間比參數比較抽象，難以通過野外觀測獲得。而且，目前國際上基于CPD模型進行雪深反演的研究主要使用X波段，暫未見使用C波段極化SAR數據的研究報道。文獻[19]假設研究區積雪各向異性介電常數不變，使用一個實測站點處的積雪密度0.07 g/cm3和假定軸間比1.5作為雪深反演算法的輸入，利用全極化TerraSAR-X數據反演整個影像的雪深，并以該實測站點的雪深進行驗證，獲得較高的反演精度。但是，由于積雪的各向異性、密度空間分布不均，僅使用一個實測點的積雪密度和假定軸間比反演整個影像的雪深，將給反演結果帶來較大誤差。文獻[20]根據實測點處CPD的正負，結合該點處實測雪深、雪密度等參數，代入CPD正演模型，求得冰粒軸間比，再借助積雪密度反演算法[21]獲取全局雪深分布，但此反演方法較為復雜，且積雪密度反演算法可能會帶來新的誤差。

本文基于CPD模型提出了新的雪深反演方法，將反演所需要的參數減少為遙感獲取的CPD數據以及實測雪深數據，根據全部實測點的雪深與CPD的關系，擬合得到最適合研究區積雪特性的反演模型系數，用它們反演研究區的雪深。為了驗證反演方法在C波段、不同積雪條件下的應用能力，使用GF-3全極化數據計算CPD，以新疆阿爾泰山南坡克蘭河上游為研究區，根據實測雪情差異劃分深雪區與淺雪區進行雪深反演，分析反演結果的差異，以期為山區雪深反演提供新思路。

1 研究區與數據

1.1 研究區概況

克蘭河流域位于新疆北部，阿爾泰山南麓，在北溫帶大陸性氣候條件影響下，所形成的積雪為典型的“干寒型”積雪[22]，年平均雪深在40 cm以上。克蘭河屬融雪徑流補給為主的河流，年平均徑流量6.0×108m3[23]。積雪融水對當地農業灌溉及畜牧業發展具有重要影響，研究該流域的積雪深度具有重要意義。

研究區范圍如圖1所示，位于克蘭河流域上游東支小東溝源頭區域。研究區西南部為戈壁，地形起伏較為平緩；東北部是阿爾泰山南麓侵蝕中山山地，海拔較高，地形起伏較大，是典型的山區環境[24]。

1.2 研究數據

(1)SAR數據。CPD為HH和VV極化通道進行極化相干得到復相干系數的相位角，SAR衛星必須具有HH和VV極化通道極化相干成像的能力才能計算CPD，目前常見的可以計算CPD的SAR衛星有Radarsat-2、TerraSAR-X/TanDEM-X、ALOS-PALSAR等。研究表明C波段是CPD對雪深變化仍有較高敏感性的頻率下限[25]，本文研究選取我國的GF-3衛星數據，探討基于CPD模型使用C波段極化SAR數據反演雪深的能力。GF-3衛星于2016年8月10日發射，最高分辨率為1 m，設計有12種成像模式[26]。綜合考慮GF-3衛星數據的成像時間和圖像覆蓋區域、成像質量、同步觀測試驗可行性等，選取2018年1月17日獲取的阿勒泰地區GF-3衛星全極化條帶數據作為實驗數據，影像范圍為47°49′—48°07′N,88°03′—88°24′E。研究所用的SAR影像相關參數見表1。

(2)野外實測數據。CPD反演雪深的可行性分析、模型參數的確定與反演結果的精度驗證均需地面實測數據的支持，因此于2018年1月17日對研究區進行星地同步積雪觀測，設計了積雪樣帶觀測和積雪屬性剖面觀測相結合的觀測方案。此時為研究區的積雪穩定期，觀測期間最高氣溫在0℃以下，未發生降水。積雪樣帶觀測共計44個樣點，主要測量雪深、積雪密度、雪層溫度等參數。積雪剖面觀測共計8個剖面點，主要測量雪深、雪層密度、雪層溫度、雪層介電常數、液態水含量等參數。雪深直接由鋼制直尺測量，取觀測點周圍3 m的范圍內3次測量的平均值，測量精度為0.1 cm；積雪密度由鋼制直尺、托板和自制的量雪筒測量，測量精度為1 g/cm3；自積雪表層向下每10 cm劃分雪層，測得積雪參數代表該雪層積雪參數，其中雪層溫度利用經過標定的針式溫度計進行測量，測量精度為0.1℃，雪層介電常數、密度和液態水含量通過雪特性分析儀SnowFork測量3次取平均值得到[27]。

表1 所用的GF-3影像的成像參數Tab.1 Parameters of the used GF-3 image

雪深測量結果如圖2所示，可見研究區雪深空間分布不均，西南部戈壁區域積雪較淺，東北部山地區域積雪較深。這是因為來自西南方向的暖濕氣團在東移過程與阿爾泰山成近正交角，地形對氣流的抬升作用顯著，有利于山區降水，使得山麓的降水多于河谷平原地區[28]。為了探究模型在不同積雪條件下的適用性，將研究區分為深雪區與淺雪區兩個區域，分別進行反演與分析。淺雪區的雪深觀測點16個，平均雪深27.5 cm；深雪區的雪深觀測點36個，平均雪深72.6 cm。

深雪區觀測剖面5個，平均積雪密度0.23 g/cm3，平均體積含水量0.64%；淺雪區觀測剖面3個，平均積雪密度0.17 g/cm3，平均體積含水量0.27%。雪層含水量都接近于0，可認為積雪為干雪，滿足CPD模型條件。分層積雪密度與液態水含量的觀測結果如圖3所示。可以看出淺雪區密度與含水量分層垂直差異小，深雪區分層積雪密度存在明顯的垂直差異，雪層密度呈中間大、底部和頂部小的特征。深雪區底層積雪的含水量高于表層積雪，垂直分布表現出較大的差異。圖4為研究區雪粒圖像，可以看出深雪區較淺雪區雪粒垂直分異更明顯，積雪重結晶變質作用更強。

圖1 研究區位置Fig.1 Location of the study area

圖2 觀測雪深點的空間分布Fig.2 Position of the measured snow depth points

圖3 積雪剖面觀測結果Fig.3 Distribution of snow depth,snow density,and snow wetness of the measured snow samples

圖4 各雪層的雪粒觀測結果Fig.4 Observed results of snow particle of each layer

(3)其他數據。根據阿勒泰國家基準氣候站提供的氣溫和雪深數據，在星地同步觀測前，1月13日前后發生了一次較大的降雪，自1月13日至2月底無較大降雪天氣過程，平均雪深約為14.5 cm。此外，采用國家基礎地理信息中心發布的全球30 m地表覆蓋數據GlobeLand 30掩膜森林等復雜散射機理的區域，使其不參與雪深反演；采用30 m分辨率的SRTM DEM數據作為數字高程數據。

2 本文研究方法

2.1 積雪各向異性結構

積雪的物理特性與其微觀幾何結構有關。在理想狀態下，干雪可以看作是冰粒和空氣的二相混合物[29]。為了模擬雪的各向異性結構，假定干雪冰粒為規則的橢球體，以橢球體中心為原點，建立三維笛卡兒坐標系(x,y,z)，z軸為冰粒晶體對稱軸，平行于重力方向，x和y平面平行于平坦地表。ax、ay和az是橢球體在x、y、z方向上的半軸長，假設ax=ay，冰粒軸間比r=ax/az，表示冰粒的形狀。r=1時，為各向同性隨機結構的新降雪；r>1時，積雪在自身重力壓實作用下迅速密實化[30]，逐漸轉變為各向異性、圓盤狀水平排列結構的雪粒[31]，r越大，各向異性越強；r<1時，積雪受溫度梯度作用發生重結晶和變質作用，形成各向異性、針狀垂直排列的變質雪粒[32-33]，r越趨近于0，各向異性越強。由新降雪到發生重結晶與變質作用變成陳雪的過程，伴隨冰粒晶體形狀及軸間比r的演變如圖5所示。

圖5 積雪冰粒演化過程Fig.5 Evolution of snow particle shape

2.2 同極化相位差模型

當電磁波穿透干雪層時，積雪粒子的大小遠小于電磁波波長，干雪可認為是均質的、體散射可忽略不計的介質[34-35]，因此后向散射信號主要來源于下墊面，因雪層折射率而延遲。文獻[36]基于Maxwell-Garnett理論將冰川冰的各向異性結構與觀測到的負相位差聯系起來，建立了CPD模型。文獻[37]進一步基于干雪微結構及其微波極化特性，將CPD模型與雪的微觀結構聯系起來。

如圖6所示，當電磁波以θ入射角穿透一個深度為SD的干雪層，冰粒在電磁波極化電場作用下極化，依據Laplace方程及其邊界條件，計算出3個方向的退極化因子Ni，i∈{x,y,z}[38]

(1)

Nx=Ny=0.5(1-Nz)

(2)

圖6 極化信號穿透雪層產生相位差示意圖(改編自文獻[37])Fig.6 Electromagnetic wave penetrates the snow layer and results in the phase difference (adapted from reference[37])

根據Maxwell-Garnett等效介質理論計算冰粒各向異性介電常數εeff,i，i∈{x,y,z}

(3)

式中，Ni為退極化因子；εice為冰的等效介電常數；εair為空氣的等效介電常數；fvol為積雪與冰的密度比，故各向異性介電常數可由積雪密度ρsnow、冰粒軸間比r計算得到。

因為折射率n與介電常數εe存在εe=n2的關系，故HH極化和VV極化的折射率nH和nV與3個軸向上的等效介電常數εeff,i、入射角θ滿足以下關系

(4)

(5)

HH和VV極化信號在雪層中的折射率不同，導致后向散射信號所傳播的路徑和距離不同，從而產生了相位差異。利用波長與相位的關系，以及極化信號傳播路徑和雪層的幾何關系，得到不同雪深引起的同極化相位差

(6)

2.3 雪深反演方法

根據式(6)可計算厚度為SD的積雪導致的相位差，由此可以推導雪深反演模型為

(7)

根據nV、nH與εeff,i的關系、εeff,i與積雪密度ρsnow、Ni的關系以及Ni與冰粒軸間比r的關系可知，在確定SAR數據計算CPD之后，輸入積雪密度和積雪冰粒軸間比參數即可得到雪深。

本文研究在文獻[19]的基礎上進一步深入，通過模擬CPD與雪深的關系可知，在假定積雪各向異性介電常數不變的情況下，CPD與雪深呈線性關系，如圖7所示。當研究區有多個實測點時，可根據實測點雪深數據與從SAR數據計算的CPD擬合得到最適合研究區積雪特性的線性模型經驗系數a、b，從而得到該區域全部像元的雪深。即

(8)

(9)

式中，γc為絕對相干系數；SVV和SHH分別為VV和HH極化單視復數影像；〈〉為濾波器。本文借助干涉雷達數據處理的軟件Gamma計算CPD，包括數據預處理、極化相干、空間濾波、地理編碼等步驟。考慮實測樣點的空間分布，設置地理編碼CPD的空間分辨率為10 m。

圖7 積雪各向異性介電常恒定下CPD與雪深的關系Fig.7 Relationship between CPD and snow depth based on constant snow anisotropic relative permittivity

3 反演結果分析

3.1 反演精度分析

深雪區與淺雪區積雪屬性差異大，區域內部積雪屬性差異小。為獲得更精確的反演結果，分區域進行雪深反演。對于同一實測點，實測雪深可以看作一個常數值，而不同濾波窗口大小計算得到的CPD不同，則反演得到的擬合關系、反演精度不同。本文使用留P法進行交叉驗證，以有效避免過擬合和欠擬合狀態，評價不同濾波窗口大小下深雪區、淺雪區的反演模型優劣。為了充分利用所有實測點數據，P取1，即每次假定一個點的雪深未知，使用其余點處實測雪深值與CPD擬合得到模型系數，反演得到該點雪深，然后使用相關系數R和均方根誤差RMSE來分析所有點上雪深實測值與反演結果的誤差，如圖8所示。

結果表明，淺雪區的最優濾波窗口為59×59像元，相關系數R為0.83，均方根誤差RMSE為2.72 cm；深雪區的最優濾波窗口為33×33像元，相關系數R為0.54，均方根誤差RMSE為11.69 cm。隨著濾波窗口大小的增加，淺雪區和深雪區的R均呈現出先上升再下降的趨勢，RMSE呈現先大幅上漲波動、然后減小、最后穩定的趨勢。在濾波窗口25×25像元以下的結果，R與RMSE出現較大波動，這是因為濾波窗口較小則不能有效抑制相干斑噪聲，相位的波動使相干性降低；當濾波窗口過大，則會平滑CPD的變化，削弱圖像的細節，也導致相干性降低。

圖8 反演結果R和RMSE與濾波窗口大小的關系Fig.8 Relationship between R,RMSE and the window size of filter

不同區域最優濾波窗口大小下雪深的空間分布如圖9所示。圖9(a)顯示深雪區反演得到的雪深自西向東逐漸減少，南部和北部的雪深小于中部區域。圖9(c)為淺雪區雪深反演的結果，可見除西北部反演雪深明顯大于其他區域以外，其他區域反演雪深呈相間分布。掩膜森林、水體等具有復雜的散射機制及受人類活動影響嚴重的人造地表區域的雪深結果如圖9(b)、(d)所示。

3.2 反演不確定性分析

地形通過影響水熱條件從而影響積雪空間分布，地形起伏帶來雷達構像幾何關系的變化從而影響雷達散射機制，對雪深的反演結果產生影響[39]。為分析地形對反演精度的影響，選取坡度作為變量與反演誤差進行相關性分析。如圖10所示，雪深反演誤差與坡度顯著相關，隨著坡度的增加，雪深的反演誤差呈現增加趨勢，且淺雪區的相關性強于深雪區。由淺雪區實測點坡度分布范圍可知，淺雪區的地形起伏較深雪區更大，由于構建CPD模型時假定冰粒的z軸平行于重力方向，復雜的地形會導致極化信號穿透雪層的幾何關系發生改變，影響模型的適用性，同時電磁波回波由于受到地形起伏干擾作用會引入過多的噪聲，從而增大反演誤差。

圖9 最優濾波窗口下的雪深反演結果Fig.9 Spatial distribution ofretrieved snow depth with optimal window size of filter

CPD模型以雪層均質、忽略體散射為基礎，然而嚴格來說，自然界中的媒介質都是不均勻的，不可避免會引入誤差。由實測結果可知，深雪區積雪密度較淺雪區大，由圖4實測雪粒圖像也可以看出深雪區經過重力壓實與變質作用導致粒子粗糙度更大，雪層均質性更差，極化信號穿透雪層時體散射增強，此時，后向散射信號不僅僅取決于下墊面，也由雪層中不同的散射體決定，有效散射中心向雪層上表面移動，因此同極化相干性受到體去相關的影響而降低[18]；根據電磁波在介質中的穿透深度公式[29]，針對本文所用的GF-3數據，對不同濕度下的積雪穿透深度進行了模擬，結果如圖11所示。可以看出C波段穿透深度隨著積雪濕度增加而顯著下降。由圖3的剖面雪深與濕度的關系可以看出，由于深雪區較強的土壤熱傳導作用導致積雪底層含水量較大，部分底層含水量超過3%，水的介電常數在很大程度上改變了雪的介電特性，大幅降低了極化信號對雪層的穿透能力，導致極化信號不能到達下墊面，因此深雪區雪深反演精度低于淺雪區。

圖10 反演誤差與坡度的關系Fig.10 Correlation between slope and inversion errors

此外，為了探討不同入射角觀測幾何條件下CPD模型反演能力的差異，使用入射角發生改變時單位雪深引起的CPD變化來表征，單位雪深引起的CPD變化越大則代表CPD對該入射角下的電磁波越敏感，反演能力越強。模擬CPD與不同入射角之間的關系如圖12所示。CPD隨SAR電磁波入射角的增加而增加，單位雪深引起的CPD變化也隨著SAR電磁波入射的增加而增加，但二者是非線性關系。當入射角增大時，HH極化與VV極化的折射率差增大，所引起的CPD變化也增大，即CPD對入射角較大的SAR電磁波更加敏感。故不同入射角觀測幾何條件下，本文方法的反演能力存在差異，它更適合具有大入射角的SAR衛星。

圖11 C波段極化信號積雪穿透深度模擬Fig.11 Simulation the penetration depth of C band polarimetric signal

圖12 模擬CPD與不同入射角之間的關系Fig.12 Simulation the relationship between CPD and incidence angle

3.3 反演方法比較

為驗證本文所提出的反演方法的優勢，與已有方法進行對比。根據文獻[19]提出的雪深反演方法，使用研究區平均積雪密度和假定軸間比1.5代入反演模型，求得雪深；根據文獻[20]提出的雪深反演方法，結合不同濾波窗口下的實測點處CPD、實測雪深、雪密度等參數求得軸間比，再結合該點處的實測積雪密度，輸入模型反演該點處的雪深。使用留P法(P=1)進行交叉驗證來評價3種反演方法優劣，比較結果見表2。根據文獻[19]提出的雪深反演方法得到掩膜后的深雪區、淺雪區的雪深分布結果分別如圖13(a)、13(c)所示，根據文獻[20]提出的雪深反演方法得到掩膜后的深雪區、淺雪區的雪深分布結果分別如圖13(b)、13(d)所示。

表2 本文方法與已有方法的雪深反演精度比較Tab.2 Comparison of the proposed method with the existing methods

圖13 已有雪深反演方法的反演結果Fig.13 Spatial distribution of retrieved snow depth based on the Patil’s and Majumdar’s methods

可以看出，對于本文研究所用GF-3數據，本文方法比另外兩種方法精度更高，且將反演所需要的參數減少為遙感獲取的CPD數據和進行模型訓練的實測雪深數據。比較結果證明了本文方法具有較好的反演積雪深度的能力。

4 結 論

本文以新疆阿勒泰克蘭河上游地區為研究區，在獲取C波段全極化GF-3數據及地面同步觀測數據的基礎上，基于CPD模型構建線性的半經驗雪深反演模型，根據實測積雪屬性差異，分深雪區、淺雪區反演雪深，探討了反演方法在不同積雪條件下的適用性，對反演不確定性進行了分析并與已有方法進行比較，主要結論包括：

(1)在積雪各向異性介電常數視為恒定的理想情況下，CPD僅是雪深的函數，可根據全部實測點的雪深數據與從SAR數據獲取的CPD數據，擬合得到最適合研究區積雪特性的模型經驗系數，反演研究區雪深。反演精度的高低與計算CPD過程中使用的濾波器的窗口大小有關，淺雪區的最優濾波窗口為59×59像元，此時R為0.83,RMSE為2.72 cm，深雪區的最優濾波窗口為33×33像元，此時R為0.54,RMSE為11.69 cm。

(2)雪深反演誤差與坡度顯著相關，隨著坡度的增加，雪深的反演誤差呈現增加的趨勢，雪深反演不確定性受雪層變質程度、含水量及衛星入射角觀測幾何條件影響，反演方法對于干燥、雪層變質結晶程度低、均質的積雪，以及具有大入射角的SAR衛星有更好的適用性。

(3)對比已有基于CPD模型的雪深反演方法，本文所提出的反演方法精度更高，并且將反演所需要的參數減少為遙感獲取的CPD數據和進行模型訓練的實測雪深數據，能夠可靠地反演積雪深度。