李震 田邦森 張平 趙常軍 孫光德 劉暢

摘要積雪是影響輻射平衡、水資源等全球能量與水循環(huán)的重要參數(shù)，遙感是監(jiān)測積雪時空分布特征及其變化趨勢的重要手段.合成孔徑雷達（SAR）以其全天時、全天候能力，成為積雪遙感中重要的研究方向.本文從SAR與積雪的特性、積雪的散射模型、基于強度和相位信息的積雪參數(shù)反演等方面，針對SAR積雪參數(shù)反演方法的理論、技術進行總結(jié)和分析，指出了目前利用SAR數(shù)據(jù)反演積雪的研究現(xiàn)狀和存在的問題，并對未來的研究方向進行了討論.關鍵詞合成孔徑雷達;積雪;參數(shù);模型;反演;進展

中圖分類號P426.63+5;P415.2

文獻標志碼A

0引言

積雪影響全球輻射平衡、地下水、徑流、冰川、動植物以及人類活動.在所有的冰凍圈要素中，積雪覆蓋的范圍最廣，冬季北半球40%以上的地區(qū)被雪覆蓋.由于新雪的反照率通常在0.8到0.9之間，而大多數(shù)陸地表面的反照率在0.1到0.3之間，因此，積雪對太陽輻射的高反射能影響區(qū)域和全球的能量平衡.并且，積雪直接影響生態(tài)和社會經(jīng)濟系統(tǒng).融雪徑流不僅在水量方面，而且在分布、質(zhì)量和季節(jié)性分布方面，支配著下游地區(qū)的徑流狀況.此外，部分區(qū)域積雪和融雪也可能導致自然災害，如雪崩和春季洪水等.

由于積雪覆蓋范圍廣、時變性強、可及性差、分布位置偏遠等特點，使得遙感成為監(jiān)測大范圍積雪時空分布特征及其變化趨勢的唯一有效技術.利用光學/多光譜遙感數(shù)據(jù)監(jiān)測積雪的歷史比SAR（合成孔徑雷達）長，在積雪范圍提取方面有一定優(yōu)勢，但會受到云層覆蓋、極夜和穿透性的制約.星載SAR數(shù)據(jù)提供了一個有價值的替代方案，它具有全天時、全天候和穿透性的特點，尤其是近些年來新的SAR衛(wèi)星不斷升空，SAR積雪參數(shù)反演算法的發(fā)展，基于SAR數(shù)據(jù)來監(jiān)測積雪收到了越來越多的關注.

目前，基于SAR的積雪參數(shù)反演的理論、技術局限性、關鍵輔助數(shù)據(jù)、發(fā)展軌跡以及未來的發(fā)展前景，還沒有全面的論述.因此，本文總結(jié)了目前利用SAR數(shù)據(jù)反演積雪的技術，比較了它們的優(yōu)缺點，并對未來的研究方向進行了討論.

1SAR與積雪的特性

相比光學和被動微波遙感，SAR是一種主動式微波遙感系統(tǒng)，發(fā)射波長在0.01至1 m之間的電磁波，并接收地球表面發(fā)出的回波，其獲取的圖像距離向的空間分辨率大約等于實際天線長度的一半，并且與平臺高度無關，故能夠以較高的空間分辨率進行積雪遙感.由于微波可以穿透云層，不受天氣和光照條件的限制，因此具有全天時、全天候成像能力.就積雪遙感而言，這些能力尤為重要，因為積雪覆蓋的山區(qū)經(jīng)常被云層覆蓋，且高緯度地區(qū)在冬季易受極夜的影響.此外，SAR可穿透一定積雪厚度，潛在地探測有關積雪狀態(tài)的信息，如積雪粒徑和雪深等.SAR傳感器的另一個優(yōu)點來自它的相位信息，不同的極化（水平、垂直）允許檢測表面的物理特性，如目標的形狀、質(zhì)地或朝向，而基于SAR記錄的相位信息，可以進行相干和干涉測量以提取積雪的范圍和深度信息.

然而，目前SAR在積雪遙感應用中存在著一些問題.首先，由于軌道設計的原因，目前星載SAR具有相對較低的時間分辨率（即回訪周期，通常長于5 d）（表1）.根據(jù)全球氣候觀測系統(tǒng)（GCOS）給出的報告，水文和氣候研究對星載積雪產(chǎn)品的時間分辨率的最低要求是1至5 d.其次，與光學傳感器不同，SAR傳感器觀測的幾何結(jié)構更為復雜，因為它們的合成多波束回波在兩個方向（側(cè)視距離向和方位向）發(fā)射和接收，由此產(chǎn)生的變形影像因景觀和傳感器角度的不同而異，從而產(chǎn)生了透視收縮、疊掩和陰影等幾何變形，在山區(qū)尤其嚴重.此外，由于SAR圖像的每個像素接收到的信號是所有地面分辨單位內(nèi)反射信號相干的總和，因此而產(chǎn)生的斑點噪聲會顯著降低圖像質(zhì)量，使得SAR圖像的解釋和分析比傳統(tǒng)的光學觀測圖像更具挑戰(zhàn)性.

受SAR信號的穿透特性的影響，干濕雪在SAR圖像中的差別很大.另一方面，SAR信號穿透積雪的深度取決于信號的波長.例如，C波段SAR在觀測干雪時的潛在穿透深度約為20 m.由于雪的粒徑在0.1至0.3 mm之間，波長較長的合成孔徑雷達信號幾乎不受阻礙地穿過積雪.在高于1 MHz的頻率下，空氣、冰和水的介電常數(shù)分別為1.0、3.0和80.當雪中的空氣、冰和水的比例隨雪的融化發(fā)生變化時，介電常數(shù)也隨之變化，從最初的1.2到2.0之間的低值變得越來越高，SAR信號的穿透深度減小到3 cm左右.因此，積雪的液態(tài)水含量會影響到散射機理，在干燥積雪條件下，主要的散射過程是積雪的體積散射和積雪/地面界面的表面散射之和;當積雪變濕時，空氣/雪界面的表面散射是其主要散射機制.表2概述了各因素對干濕雪散射的影響以及后向散射特性.

然而，由于積雪是一個復雜的多層結(jié)構，雪的粒度、密度、深度、分層、雜質(zhì)含量和表面粗糙度都會影響其后向散射.除積雪特性外，許多外部因素影響散射機制，如雷達波長、極化、入射角、表面粗糙度和介電性能.

2積雪的散射模型

目前關于積雪特性、積雪參數(shù)由底層到表層的變化規(guī)律等研究，多集中在積雪演化、積雪熱力學相關研究中，如Lehning團隊[1-2]基于有限元方法的物理演化模型SNOWPACK.

積雪這種占空比在10%～40%之間的致密性介質(zhì)，粒子之間的散射不是相互獨立的，聚集性散射效應和近場散射波的相干特性必須考慮.基于聚集性散射效應和近場散射波的相干特性，Tsang等[3]最早提出了密集介質(zhì)輻射傳輸理論（Dense Media Radiative Transfer，DMRT），QCA/DMRT理論很好地考慮了雪粒子間的相干散射項.基于該理論的遙感探測積雪方法得到廣泛研究，目前，這些工作主要集中在利用星載微波輻射計數(shù)據(jù)（SMM/I、AMSR、AMSR-E）獲取雪水當量和雪深等參數(shù)[4-6].

但是，由于自然界中的積雪受重力、風以及融凍等因素的影響，在不同積雪深處，雪粒子的大小、密度、含水量等會明顯不同，具有典型的分層特性，這會使衛(wèi)星觀測數(shù)據(jù)與QCA（Quasi Crystalline Approximation）/DMRT理論的模擬結(jié)果產(chǎn)生較大誤差[7].因此，Tedesco等提出一個包含濕雪層和干雪層的雙層DMRT微波輻射模型[5]，但是模型中濕雪層輻射被認為是黑體輻射，所以該模型實際上并不能反映多層積雪的散射和輻射特性.最近，Liang等[8]發(fā)展了一個基于QCA/DMRT理論的多層積雪輻射模型，每層中的輻射特性遵循QCA/DMRT，不同層間入射角-反射角-透射角的關系遵循Snell法則，雪層-地表之間的輻射量用經(jīng)典的Q-H模型[9]表示.該模型的模擬數(shù)據(jù)與CLPX實驗中的實測數(shù)據(jù)比較，能得到比傳統(tǒng)單層密集介質(zhì)模型更好的結(jié)果[8].Brandt等[10]利用時域有限差分（FDTD）模型，研究了探地雷達對由于季節(jié)性侵滲產(chǎn)生的粒雪響應，他們發(fā)現(xiàn)冰層厚度、各層之間的距離以及是否出現(xiàn)橢圓形冰棱都與雷達響應密切相關.Gherboudj等[11]利用致密介質(zhì)相和幅度校正理論（Dense Medium Phase and Amplitude Correction Theory，DM-PACT）和MD（Matrix Doubling）方法，針對淡水面上形成的多層河冰開發(fā)了多層后向散射模型，并分析不同冰類型對雷達后向散射信號的影響.Du等[12]利用準晶近似相干位（QCA-CP）和橢圓形瑞利散射相矩陣，通過MD和AIEM模型發(fā)展了考慮多次散射的多層積雪散射模型，研究顯示出多層多次散射對HH極化和交叉極化表達具有明顯的改進.上述結(jié)果表明了密集介質(zhì)散射、輻射理論模型的研究已經(jīng)發(fā)展到了多層模型階段.

InSAR正演模型依賴于介質(zhì)與電磁波的相互作用.早期發(fā)展的積雪電磁散射模型是隨機離散介質(zhì)模型，該理論假設各散射粒子的散射是相互獨立的，散射體之間的相互作用發(fā)生在單個散射體散射的遠場且相干作用很小，只考慮輻射強度的疊加，而不考慮輻射場的相關性.該理論下有Rayleigh（瑞利）散射和Mie散射，前者發(fā)生在散射粒子的尺度大小遠小于波長時，后者發(fā)生在粒子大小與波長尺度相當時[4，13].2016年基于QCA理論的致密介質(zhì)InSAR正演模型[14]可以較為準確地刻畫單層積雪，該模型假定積雪粒子為球形粒子，用QCA理論求解積雪介質(zhì)內(nèi)的相干場，并將雷達觀測到的總電場近似為相干場，然后結(jié)合兩次觀測的觀測條件計算復相干系數(shù).然而，由于積雪的分層特性，電磁波在多層積雪中的傳播路徑較單層積雪復雜且與積雪子層參數(shù)序列有關，而傳播路徑的長度、子層折射率則直接影響InSAR相干系數(shù)，因此單層積雪InSAR正演模型無法模擬多層積雪的相干系數(shù).隨后，李震等[15]建立了一種新型多層積雪模型，該模型基于QCA近似，考慮單多層積雪的差異對積雪子層分別建模，利用菲涅爾定律描述多層積雪子層間的電磁波透射，結(jié)合InSAR相干過程依次求出雷達收到的單個雪粒子的相干場、子層積雪的相干散射場共軛積，構建了多層積雪InSAR正演模型，并開展地基SAR實驗驗證其有效性.

3SAR積雪參數(shù)反演

SAR積雪反演主要包括兩類算法：一是依賴于多頻多極化后向散射強度觀測與眾多積雪、土壤和植被覆蓋等參數(shù)關系的經(jīng)驗算法[16];另一種是利用干涉SAR測量數(shù)據(jù)在觀測積雪時產(chǎn)生的干涉相移，直接估算雪水當量[17].

3.1基于強度的積雪參數(shù)反演

3.1.1采用L、C、X波段進行積雪參數(shù)反演

測量積雪的密度和雪深（SWE）是基于后向散射雷達信號與雪量的相互作用建立的.Shi等[6，18]利用SIR-C/X-SAR獲得的基于多頻率（L、C和X波段）雙極化（VV和HH）SAR數(shù)據(jù)，估算了SWE.該方法是使用L波段來估計雪密度和下層地表介電和粗糙度特性，利用C和X波段，在最小化地表后向散射信號的影響下，估算雪深和粒徑.

該方法通過使用IEM模型[19]基于多種入射角、介電特性、粗糙度以及入射波數(shù)（對應100～550 kg/m3的雪密度范圍）建立后向散射系數(shù)數(shù)據(jù)庫來完成.然后，通過回歸分析，確定了每個入射角和波數(shù)下，HH和VV后向散射特征與表面介電常數(shù)和粗糙度之間的關系，最終導出僅使用SAR測量值估計雪密度的算法.其中雪密度可根據(jù)Looyenga的半經(jīng)驗介電公式[20]估算：

εs=1.0+1.599 5ρs+1.861ρ3s .（1）

用半經(jīng)驗模型來描述積雪-地面的相互作用，這種模型比獨立散射假設下的公式更能真實地描述雪地相互作用的成分.高山地區(qū)的自然表面相當粗糙，非相干分量對積雪-地面的相互作用具有重要貢獻.通常根據(jù)積雪體散射反照率、光學厚度、地面反射率和地面均方根高度來建立半經(jīng)驗模型.

建立半經(jīng)驗模型來描述C波段和X波段地面后向散射分量之間的關系，需要參數(shù)化C波段和X波段與積雪消光特性之間的關系.因為消光性質(zhì)是高度相關的，發(fā)展出C波段和X波段消光關系可以將積雪后向散射分量中的未知量可以減少到兩個：體散射反照率ω和光學厚度τ.最終未知數(shù)僅包括τ，ω和表面分量σsvv（X）.因此，這3個未知數(shù)可以用3個SAR測量值在每個像素上進行數(shù)值求解：σtvv（C）、σthh（C）和σtvv（X）.

為了估計雪深d，需要在估計的τ中分離消光系數(shù)ke和雪深d，需要確定消光系數(shù)的絕對值.這可以通過估計吸收系數(shù)來實現(xiàn)：

kaX=1.334+1.218 2lgVi-

3.421 7lg（τX[1-ω（X）]τC[1-ω（C）]） ，（2）

其中Vi是冰的體積分數(shù)，可通過使用L波段從每個像素處估計的雪密度得出測量.

τ1-ω=ka·d=τa是除去散射后每個頻率光學厚度的吸收部分.最終，雪的深度可以通過公式估計得到：

d=τX[1-ω（X）]ka（X） .（3）

隨后Thakur等[21-22]又對該方法進行了實驗和改進.

3.1.2采用X、Ku波段進行積雪參數(shù)反演

在X波段，高達60%的信號來自于較厚的積雪.因此，SAR測量結(jié)果對積雪更加敏感，使用X波段或更高頻率的SAR數(shù)據(jù)反演積雪參數(shù)更加可靠.選擇合適的雷達頻率，通過遙感觀測積雪特性，一方面要考慮到穿透較厚積雪的能力，同時要考慮到雙向傳播路徑，另一方面，利用散射源的期望強度從積雪中獲得適合反演的回波信號.

“寒區(qū)水文高分辨率觀測計劃”（CoreH2O），作為歐空局地球探測器計劃[16，23]的一個候選衛(wèi)星計劃，提出了一種在Ku波段（17.2 GHz）和X波段（9.6 GHz）頻率下工作的雙頻合成孔徑雷達.較短的波長（Ka波段）不能提供足夠的穿透力，信號主要由雪層頂部的微觀結(jié)構（粒度和形狀）控制[24].在較長波段（C波段、L波段）的雷達回波主要是干雪覆蓋下的地面的后向散射.

來自發(fā)射和接收極化組合pq的積雪地面的總后向散射σt，由以下貢獻組成：

σtpq=σaspq+σvpq+σgvpq+σg′pq ，（4）

其中，σas表示空氣-積雪界面的散射，σv是積雪直接體散射項，σgv是地面-雪和雪-地面的相互作用，σg′是通過積雪后的地表后向散射（圖1）[23].單個貢獻的大小取決于積雪的散射和吸收特性以及背景介質(zhì)的后向散射信號.

干雪后向散射的主要貢獻主要來自積雪和積雪下面的地面，空氣/積雪界面的散射最多達到百分之幾.圖2顯示了使用致密介質(zhì)輻射傳輸理論（DMRT）模擬的地面干雪的后向散射，該模型使用DMRT來描述雪的微觀結(jié)構并導出散射特性[25].模擬結(jié)果表明：在9.6 GHz（X波段），對80 cm以下的積雪，土壤的后向散射貢獻占主導地位（對應于自由模型參數(shù)（SWE，雪水當量）=144 mm）;而在17.2 GHz（Ku波段），積雪的后向散射貢獻占主導地位.另外，理論和實驗數(shù)據(jù)都表明了積雪的微觀結(jié)構對后向散射敏感度的重要性，在反演雪深和SWE時需要考慮后向散射敏感度.

Lemmetyinen等在芬蘭北部的Sodankyl利用散射計（SnowScat）觀測了X波段和Ku波段不同粒徑條件下， SWE與散射觀測之間的關系[26].結(jié)果表明在40°的入射角下，對于SWE≤100 mm，Ku波段的靈敏度約為30 mm（SWE）/dB;而對于較高的SWE值，Ku波段的靈敏度下降到40 mm（SWE）/dB.另一方面，X波段的后向散射對SWE的敏感性相對較小.在后向散射反演積雪SWE時，要考慮積雪粒徑的影響.粗顆粒的積雪的后向散射對SWE的敏感性較高，而細顆粒的后向散射對SWE的敏感性較低（圖3）.King等[27]在加拿大哈德遜灣西南海岸一個積雪厚度較淺（<25 cm）、粒徑較大的實驗區(qū)測得Ku波段后向散射測量的靈敏度為15 mm（SWE）/dB.

由于積雪SWE的反演需要依據(jù)積雪的體散射和積雪下介質(zhì)的后向散射貢獻，并受積雪的粒徑等參數(shù)的影響，需要的觀測量較多，通常需要一些先驗知識參與反演.如從Ku和X波段后向散射值中反演積雪SWE的基本算法，需要關于雪特性的統(tǒng)計先驗信息進行正則化.采用約束最小化方法，使用成本函數(shù)迭代兩個自由變量[16，23]：

J（x）=∑4i=112vari[Fx1，x2;c1，…，cr-σ0i]2+∑2j=11λ2j（xj-〈xj〉）2 ，（5）

其中F（·）是計算積雪-地面后向散射的前向模型;i是雷達信息參數(shù)（i=1，2，3，4→X-vv，X-vh，Ku-vv，Ku-vh）;x1，x2是自由模型參數(shù)（SWE和有效晶粒半徑（RE））;c1，…，cr是正向模型的配置參數(shù);σ0i是雷達觀測i中測量的后向散射系數(shù);vari是由于測量不確定度引起的后向散射系數(shù)的方差;j是自由模型參數(shù)的指數(shù);〈xj〉是正則化參數(shù)的平均值;λj是正則化參數(shù)的標準差.

該算法對單個雪層使用輻射傳輸模型，模型將后向散射特征與積雪的物理特性相關聯(lián).對兩個積雪參數(shù)SWE和RE執(zhí)行迭代.迭代法將測量和正演模擬后向散射數(shù)據(jù)之間的差異最小化.所需的先驗信息包括結(jié)構參數(shù)（雪密度和溫度）和從降雪開始前的后向散射測量獲得的無雪地面背景信號.對于自由變量SWE和RE的平均值和標準差，需要初始值.結(jié)構參數(shù)和自由變量的初始值可從氣候積雪統(tǒng)計中獲得，由數(shù)值氣象數(shù)據(jù)驅(qū)動的積雪過程模型提供的時間序列數(shù)據(jù)提供了較好的估計[28-29].

3.2基于InSAR積雪參數(shù)反演

3.2.1利用干涉相位進行積雪參數(shù)反演

星載合成孔徑雷達差分干涉（DInSAR）技術利用電磁波穿透雪層的傳播延遲來反演雪水當量（SWE）.采用降雪前后的SAR圖像做干涉處理，獲取的相位、相干系數(shù)等參數(shù)與積雪的密度、粒徑、深度信息相關，可反算雪水當量和雪深.Guneriussen等最早利用InSAR反演雪水當量[17]，該方法對雪水當量變化十分敏感，如Envisat C波段波長5.62 cm，入射角23°，雪密度為0.3 g/cm3的情況下，雪深10.2 cm就會發(fā)生相位纏繞，對于L波段，其對應的雪深為58 cm，長波段、小入射角可以反演更大的雪深和雪水當量.

圖4描述了雷達波的傳播路徑，T1時刻為無雪地面，T2時刻覆蓋有厚度為ds并且介電常數(shù)為ε的雪層.θi為入射角，θr為折射角.兩次入射波的光程長度差異如下：

ΔRop=εΔRs-ΔRa.（6）

光程長度差異是因為兩次的傳播路徑不同，這是因為積雪和空氣的介電常數(shù)不同導致折射率和傳播速度有差異.假設雪層對雷達波的吸收和損失可以忽略，對于一塊均勻覆蓋厚度為ds且反射率為n=ε的雪層，其傳播路徑延遲如下：

ΔRop=dscosθi-ε-sin2θi，（7）

其中，θi為雷達波在雪層上表面的入射角.

由于雷達波具有周期性，相差整周波長時相位相同，這將產(chǎn)生2π的相位模糊[30].因此，為了避免相位纏繞，兩次雷達影像期間的雪水當量的變化量需要小于2π相位差對應的量.當入射角θi=30°，L波段的2π相位差對應于SWE=128 mm，而C波段對應于SWE=31 mm.隨著入射角的增加和波長的減短，整周模糊所對應的雪水當量會減少，但InSAR相位對雪水當量的敏感性增加了，表明在選取最優(yōu)波段時需要在雪水當量的變化量和探測的敏感性上綜合考慮.

Engen等[31]在Guneriussen等[17]工作的基礎上，采用了InSAR子孔徑分解技術（delta-K）發(fā)展了一種新的InSAR雪水當量反演方法.該方法將傳統(tǒng)的通道分解為兩個通道和，克服了傳統(tǒng)InSAR雪水當量反演算法中存在的相位纏繞問題，其采用有雪和無雪的兩景ERS-1影像對雪水當量估算提高了估算精度.此后，該方法又進行了不斷實驗和改進[32-37]，結(jié)果表明該方法的精度主要受限于失相干和相位纏繞，失相干將導致干涉相位精度變差.

為了克服InSAR雪水當量估算算法中的相位纏繞問題，多位學者又進行了通道分割、多波段組合法等方法研究解決相位纏繞[38-39]，但積雪整體垂直結(jié)構還不能利用干涉測量獲得.為解決此問題，目前正在發(fā)展InSAR正演模型，表達了積雪在一定觀測條件（波長、極化方式、基線等）下，積雪參數(shù)、地面參數(shù)和InSAR復相干系數(shù)之間的定量關系，解決微波穿透積雪內(nèi)部、冰粒子之間發(fā)生多次散射、積雪-地面之間的散射、多層積雪雪層之間的透射和散射問題[15].對積雪子層分別建模，聯(lián)合InSAR相干過程，構建出多層積雪InSAR正演模型，分析積雪分層對相干場共軛積的影響，發(fā)展多層積雪InSAR反演算法是重要的研究方向[14，40].

3.2.2利用相干系數(shù)進行積雪參數(shù)反演

通常，相干系數(shù)的幅度γ可以采用最大似然評估法求取窗口中N個像素得到[30]：

γ=∑Nn=1Vn1V*（n）2∑Nn=1Vn12∑Nn=1Vn22，（8）

其中，V1和V2是兩景SAR影像中相同位置像素的復數(shù)值， *表示復共軛.γ的范圍在0到1之間.在InSAR中，去相干是指誤差源導致的噪聲，誤差源的相關長度小于相干評估窗口.總相干系數(shù)可以表示為以下幾個相關因子的乘積：

γtot=γthermal×γspatial×γvolume×γprocessing×γtemporal.（9）

相關因子γthermal是由于系統(tǒng)噪聲導致的去相干，γspatial和γvolume是由于兩次雷達觀測不同導致的空間和體散射去相干，γprocessing是由于SAR影像處理和配準過程帶來的去相干.前4個因子的影響可以通過選取合適的傳感器參數(shù)以及衛(wèi)星較小的基線來控制到最小[30].InSAR雪水當量反演的一個重要影響因素是時間去相干，它是由于觀測目標的物理特性發(fā)生變化所導致的.影響積雪相干性的主要因素是表面融雪、降雪以及風吹和沉積導致的雪粒移動.融雪時雷達波的穿透性降低，因此InSAR雪水當量的反演不適用于濕雪.

Zebker等[41]分析了格林蘭島冰蓋積雪干涉的失相干源，推導出體失相干與積雪穿透深度（圖5）的反演表達式，具體如下：

對SAR圖像上每一個復像素點對應的分辨單元，天線測量的信號可表達為

s1=f（y，z）exp-j4πλr0+ysin θ1-zcos θr1·

Wr（y，z）dydz， （10）

其中，f（y，z）是散射體內(nèi)每個散射單元的復后向散射函數(shù)，λ是波長，r0是SAR天線至分辨單元中心的斜距，θ1和θr1分別為入射角和折射角，Wr（y，z）為距離相關的脈沖響應函數(shù).因此，干涉或交叉相關兩次觀測信號可寫為

s1s*2=f（y，z）f*y′，z′e-jβ（y，z）·Wr（y，z）W*ry′，z′dydzdy′dz′， （11）

其中，β（y，z）=4πλ[y（sin θ1-sin θ2）-z（cos θr1-cos θr2）].假設積雪由均勻分布的非相干散射點構成，則，〈f（y，z）f*（y′，z′）〉=σ0（y，z）δ（y-y′，z-z′）.其中，σ0（y，z）為歸一化的雷達散射截面積.

假設分辨單元內(nèi)散射體的歸一化雷達散射截面僅與積雪深度有關，即σ0（y，z）=σ0z.

假設散射介質(zhì)是均勻有耗的，并且在雷達所能探測到的深度范圍內(nèi)擴展，即：

σ0z=σ0ae2zdcos θr，（12）

其中z在介質(zhì)內(nèi)部為負值，d為穿透深度，cos θr考慮了無線電波在介質(zhì)中的非垂直傳播.將上式代入傅里葉積分式并歸一化，可得：

γvolume=11+2πdB⊥λr0tanθ2.（13）

因此，如果式（9）中其他失相干系數(shù)已知，通過式（13）即可從無偏相干系數(shù)γ反推出雷達波的穿透深度d.

該方法的有效性在地基SAR系統(tǒng)中得到了很好的驗證.Leinss等[39]在芬蘭北部采用工作頻率在10～16 GHz的地基散射計得到了高時間分辨率的干涉數(shù)據(jù)，獲得了4個冬季積雪變化的時間序列.由于采用了4 h的時間分辨率，時間去相干得到了很好的抑制.反演得到的小于200 mm的雪水當量變化與實測數(shù)據(jù)結(jié)果一致，RMSE位于±5到±15 mm之間.同樣的方法用5.9 GHz（C波段）和2.2 GHz（S波段）的地基SAR系統(tǒng)開展的實驗，也得到了成功的驗證[42].同時，實驗也發(fā)現(xiàn)這2 d頻率之間的時間去相干差異顯著，從相干像素的統(tǒng)計結(jié)果來看，C波段相干性2 d內(nèi)從平均值γ=0.83降到γ=0.28，然而S波段相干性從γ=0.91降到γ=0.67，這與波長對積雪地表相干性影響的理論預期是一致的[43].

3.3層析SAR積雪參數(shù)反演

層析SAR技術于20世紀90年代末期發(fā)展起來，1999年德國宇航局（DLR）實現(xiàn)了機載L波段實測數(shù)據(jù)的三維成像，用德國Oberpfaffcnhofen地區(qū)的14軌L波段全極化數(shù)據(jù)開展實驗，證明了層析SAR三維成像的可行性和應用價值[44].近年來隨著機載、星載SAR系統(tǒng)的不斷成熟和SAR層析三維成像技術的發(fā)展，高質(zhì)量多基線SAR影像的成功獲取，更適合SAR層析三維成像實際應用的成像算法被發(fā)展出來.

層析成像通常在目標不同的高度上開展多次觀測，圖6顯示的是層析SAR成像幾何模型，共有Nel條基線的數(shù)據(jù)，方位向合成孔徑長度為Laz[45].采用時域反投影算法得到距離向壓縮的信號為

Srr0，p=∫S（p，k）exp（-ikr0-p）dk，（14）

其中，r0和p表示目標和傳感器的位置，S（p，k）為雷達記錄的回波信號，k=4πλ為波數(shù).接下來再對方位向和高度向合成孔徑就可以得到三維成像結(jié)果：

Ir0=∑Nelnz=1∑Naznx=1Srr0，p（nz，nx）2.（15）

此時域反投影三維成像方法在穿透積雪介質(zhì)時，由于積雪的分層不均一性，會造成電磁波在不同層的傳播方向和速度的變化.由于分層折射率的差異造成的傳播時間延遲如下：

τ=1c0（d0+∑Ni=1nidi），（16）

其中，c0為空氣中電磁波的傳播速度，d0是到達積雪界面之前空氣中的傳播距離，di和ni分別對應每層的厚度和折射率.如果對此不做校正，會造成層析成像結(jié)果與實際位置出現(xiàn)偏差，影響后續(xù)的三維積雪參數(shù)反演.

Rekioua等[46]基于地面平行的假設，采用了一種迭代方法得到每一層的厚度和折射率.其基本思路是，當采用正確的折射率做三維成像時得到與地面平行的成像結(jié)果.據(jù)此，從最頂層開始，依次迭代校正所有層，完成此校正后就可以得到正確位置的三維成像結(jié)果，并且得到了積雪的子層厚度以及反射率.依據(jù)反射率可以計算分層介電常數(shù)和分層雪密度.該實驗通過與地面雪坑實測數(shù)據(jù)做比對，表明用層析SAR反演積雪分層厚度和分層介電常數(shù)等參數(shù)是可行的.

SAR層析技術用于積雪方面可以得到積雪內(nèi)部分層結(jié)構信息，進而反演獲得多層積雪參數(shù).Tebaldini等[47]于2010年12月開展了地基合成孔徑雷達（Ground-Based Synthetic Aperture Radar，GBSAR）積雪實驗并成功獲得三維成像結(jié)果.Famil等[48]于2013年2月采用X和Ku波段GBSAR對積雪做三維成像，層析SAR成像結(jié)果能夠得到清晰的積雪分層.

此外，SAR層析成像技術應用于冰川領域，可以得到冰川表面下一定深度的內(nèi)部結(jié)構信息.Banda等[49]在2012年5月采用P波段全極化數(shù)據(jù)在格陵蘭西南地區(qū)對冰川做三維成像，得到P波段對冰川三維成像的穿透深度大約為20～60 m.Tebaldini等[50]在2014年3月采用L波段全極化數(shù)據(jù)對奧地利Mittelbergferner冰川做三維成像，結(jié)果與地面探冰雷達的結(jié)果對比，有很好的一致性.

4結(jié)論

積雪是控制全球輻射平衡、水文現(xiàn)象、植被覆蓋和影響人類活動的重要組成部分之一.星載合成孔徑雷達的極化和相位數(shù)據(jù)提供了關于積雪特性的獨特信息，使定量反演積雪參數(shù)成為可能.但是，由于SAR探測積雪所依賴的積雪后向散射系數(shù)與積雪參數(shù)（雪密度、介電常數(shù)、粒徑、濕度、深度、粗糙度等）以及雪下地表之間復雜的相互關系，目前并無成熟穩(wěn)定的運行化方法.另外，植被和山區(qū)地形的影響，也制約了算法的應用.

近30年來，基于SAR的積雪參數(shù)反演算法得到了迅速而深入的發(fā)展，包括設計新的星載SAR傳感器、新的積雪反演算法、提高衍生產(chǎn)品的空間和時間分辨率、提高精度以及加深對傳輸機制的理解.當前，基于C波段SAR的算法占據(jù)了研究的主導地位，但X波段SAR的發(fā)展為其提供了一個很有前途的選擇.由于C波段SAR衛(wèi)星歷史長，與L波段SAR相比具有更好的積雪探測能力，故比其他傳感器更常被用于雪地覆蓋探測.然而，最近的許多研究已經(jīng)證明，X和Ku波段更適合于探測干雪.考慮到新的和計劃的X波段任務的數(shù)量，基于X波段的積雪探測算法將會有更大的發(fā)展.

此外，隨著最近發(fā)射的SAR衛(wèi)星的技術進步，如更寬的探測范圍、更短的回訪時間和全極化，使得基于干涉合成孔徑雷達的積雪參數(shù)反演更加有效.因為擴大的覆蓋范圍可以提高反演的效率，縮短的重訪時間可以支持基于InSAR的方法保持更有用的相干性.全極化衛(wèi)星數(shù)據(jù)豐富了基于極化信息分解的PolSAR技術的應用，相關的應用有待于進一步的探索，其大范圍、長時間序列的適用性有待驗證，新的普遍適用的模型需要進一步發(fā)展.層析SAR技術用于積雪探測可以獲得多層積雪參數(shù)，但數(shù)據(jù)獲取、姿態(tài)參數(shù)處理和影響要素去除方法等方面還需要進一步解決.

另一方面， SAR與其他傳感器的協(xié)同使用會有助于算法的提高和應用程度，但SAR與其他傳感器（如光學和被動微波）協(xié)同提高積雪參數(shù)反演質(zhì)量的研究還不成熟，需要進一步研究.

目前，一些針對積雪特性與模型的研究的序列觀測計劃正在實施，如ESA的“寒區(qū)水文高分辨率觀測計劃”（CoreH2O），NASA的SnowEx計劃、GCPEX計劃等，目的是發(fā)展針對積雪空間觀測的多種設備與技術.其中，SnowEx是一個延續(xù)5年的項目，該項目采用機載和地面同步測量實驗，已經(jīng)在2016—2017年、2018—2019年冬季利用5種型號的飛行器搭載了10種不同的傳感器獲取美國典型區(qū)域的積雪數(shù)據(jù)，后續(xù)還將在2019—2020年冬季繼續(xù)獲取數(shù)據(jù)，以研究全球積雪SWE反演的最佳方法.

參考文獻

References

[1]Bartelt P，Lehning M.A physical SNOWPACK model for the Swiss avalanche warning[J].Cold Regions Science and Technology，2002，35（3）：123-145

[2]Lehning M，Bartelt P，Brown B，et al.A physical SNOWPACK model for the Swiss avalanche warning[J].Cold Regions Science and Technology，2002，35（3）：147-167

[3]Tsang L，Kong J A，Shin R T.Theory of microwave remote sensing[M].Wiley Interscience，1985

[4]Chen C T，Tsang L，Guo J，et al.Frequency dependence of scattering and extinction of dense media based on three-dimensional simulations of Maxwell's equations with applications to snow[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41（8）：1844-1852

[5]Tedesco M，Kim E J.Retrieval of dry-snow parameters from microwave radiometric data using a dense-medium model and genetic algorithms[J].IEEE Transactions on Geoscience and Remote Sensing，2006，44（8）：2143-2151

[6]Shi J，Dozier J.Estimation of snow water equivalence using SIR-C/X-SAR.I.inferring snow density and subsurface properties[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38（6）：2465-2474

[7]Rosenfeld S，Grody N.Anomalous microwave spectra of snow cover observed from Special Sensor Microwave/Imager measurements[J].Journal of Geophysical Research：Atmospheres，2000，105（D11）：14913-14925

[8]Liang D，Xu X L，Tsang L，et al.The effects of layers in dry snow on its passive microwave emissions using dense media radiative transfer theory based on the quasi crystalline approximation （QCA/DMRT）[J].IEEE Transactions on Geoscience and Remote Sensing，2008，46（11）：3663-3671

[9]Wang J R，O'Neill P E，Jackson T J，et al.Multifrequency measurements of the effects of soil moisture，soil texture，and surface roughness[J].IEEE Transactions on Geoscience and Remote Sensing，1983，GE-21（1）：44-51

[10]Brandt O，Langley K，Giannopoulos A，et al.Radar response of firn exposed to seasonal percolation，validation using cores and FDTD modeling[J].IEEE Transactions on Geoscience and Remote Sensing，2009，47（8）：2773-2786

[11]Gherboudj I，Bernier M，Leconte R.A backscatter modeling for river ice：analysis and numerical results[J].IEEE Transactions on Geoscience and Remote Sensing，2010，48（4）：1788-1798

[12]Du J Y，Shi J C，Rott H.Comparison between a multi-scattering and multi-layer snow scattering model and its parameterized snow backscattering model[J].Remote Sensing of Environment，2010，114（5）：1089-1098

[13]Chen C T，Tsang L，Chang A T C，et al.Scattering by densely distributed moderate size particles with size distribution and applications to snow[C]∥IEEE International Geoscience & Remote Sensing Symposium，1999，DOI：10.1109/IGARSS.1999.771533

[14]Lei Y，Siqueira P，Treuhaft R.A dense medium electromagnetic scattering model for the InSAR correlation of snow[J].Radio Science，2016，51（5）：461-480

[15]李震，李治顯，田幫森，等.基于準晶體近似的多層積雪InSAR散射模型[J].中國科學（地球科學），2018，48（8）：1070-1084

LI Zhen，LI Zhixian，TIAN Bangsen，et al.An InSAR scattering model for multi-layer snow based on Quasi-Crystalline Approximation （QCA） theory[J].Scientia Sinica（Terrae），2018，48（8）：1070-1084

[16]Rott H，Yueh S H，Cline D W，et al.Cold regions hydrology high-resolution observatory for snow and cold land processes[J].Proceedings of the IEEE，2010，98（5）：752-765

[17]Guneriussen T，Hogda K A，Johnsen H，et al.InSAR for estimation of changes in snow water equivalent of dry snow[J].IEEE Transactions on Geoscience and Remote Sensing，2001，39（10）：2101-2108

[18]Shi J，Dozier J.Estimation of snow water equivalence using SIR-C/X-SAR.II.Inferring snow depth and particle size[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38（6）：2475-2488

[19]Fung A K.Microwave scattering and emission models and their applications [M].Boston：Artech House，1994

[20]Looyenga H.Dielectric constants of heterogeneous mixtures[J].Physica，1965，31（3）：401-406

[21]Thakur P K，Garg R D，Aggarwal S P，et al.Snow density retrieval using SAR data：algorithm validation and applications in part of North Western Himalaya[J].The Cryosphere Discussions，2013，7（3）：1927-1960

[22]Thakur P K，Aggarwal S P，Arun G，et al.Estimation of snow cover area，snow physical properties and glacier classification in parts of Western Himalayas using C-band SAR data[J].Journal of the Indian Society of Remote Sensing，2017，45（3）：525-539

[23]ESA.Report for mission selection：CoReH2O，ESA SP-1324/2[R].Noordwijk：European Space Agency，2012

[24]Strozzi T，Wiesmann A，Mtzler C.Active microwave signatures of snow covers at 5.3 and 35 GHz[J].Radio Science，1997，32（2）：479-495

[25]Xu X L，Tsang L，Yueh S.Electromagnetic models of Co/cross polarization of bicontinuous/DMRT in radar remote sensing of terrestrial snow at X-?and Ku-band for CoReH2O and SCLP applications[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2012，5（3）：1024-1032

[26]Lemmetyinen J，Kontu A，Pulliainen J，et al.Nordicsnow radar experiment[J].Geoscientific Instrumentation，Methods and Data Systems，2016，5（2）：403-415

[27]King J，Kelly R，Kasurak A，et al.Spatio-temporal influence of tundra snow properties on Ku-band （17.2 GHz） backscatter[J].Journal of Glaciology，2015，61（226）：267-279

[28]Brun E，Vionnet V，Boone A，et al.Simulation of northern Eurasian local snow depth，mass，and density using a detailed snowpack model and meteorological reanalyses[J].Journal of Hydrometeorology，2013，14（1）：203-219

[29]Carmagnola C M，Morin S，Lafaysse M，et al.Implementation and evaluation of prognostic representations of the optical diameter of snow in the detailed snowpack model SURFEX/ISBA-Crocus[J].The Cryosphere Discussions，2013，7（5）：4443-4500

[30]Hanssen R F.Radar interferometry：data interpretation and error analysis[M].Netherlands：Springer，2001

[31]Engen G，Guneriussen T，Overrein O.Delta-K interferometric SAR technique for snow water equivalent （SWE） retrieval[J].IEEE Geoscience and Remote Sensing Letters，2004，1（2）：57-61

[32]Deeb E J，F(xiàn)orster R R，Kane D L.Monitoring snowpack evolution using interferometric synthetic aperture radar on the North Slope of Alaska，USA[J].International Journal of Remote Sensing，2011，32（14）：3985-4003

[33]李暉，肖鵬峰，馮學智，等.基于重軌InSAR的積雪深度反演方法[J].冰川凍土，2014，36（3）：517-526

LI Hui，XIAO Pengfeng，F(xiàn)ENG Xuezhi，et al.Snow depth derived from repeat-pass InSAR sounding[J].Journal of Glaciology and Geocryology，2014，36（3）：517-526

[34]Li H，Wang Z，He G J，et al.Estimating snow depth and snow water equivalence using repeat-pass interferometric SAR in the northern piedmont region of the Tianshan mountains[J].Journal of Sensors，2017，2017：1-17

[35]Liu Y，Li L H，Yang J M，et al.Estimating snow depth using multi-source data fusion based on the D-InSAR method and 3DVAR fusion algorithm[J].Remote Sensing，2017，9（11）：1195

[36]劉洋，李蘭海，楊金明，等.D-InSAR技術的積雪深度反演[J].遙感學報，2018，22（5）：802-809

LIU Yang，LI Lanhai，YANG Jinming，et al.Snow depth inversion based on D-InSAR method[J].Journal of Remote Sensing，2018，22（5）：802-809

[37]Conde V，Nico G，Mateus P，et al.On the estimation of temporal changes of snow water equivalent by spaceborne SAR interferometry：a new application for the Sentinel-1 mission[J].Journal of Hydrology and Hydromechanics，2019，67（1）：93-100

[38]Evans J R，Kruse F A.Determination of snow depth using elevation differences determined by interferometric SAR （InSAR）[C]∥IEEE Geoscience and Remote Sensing Symposium，2014：962-965

[39]Leinss S，Wiesmann A，Lemmetyinen J，et al.Snow water equivalent of dry snow measured by differential interferometry[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2015，8（8）：3773-3790

[40]Li Z，Li Z X，Tian B S，et al.An InSAR scattering model for multi-layer snow based on Quasi-Crystalline Approximation （QCA） theory[J].Science China Earth Sciences，2018，61（8）：1112-1126

[41]Zebker H A，Weber Hoen E.Penetration depths inferred from interferometric volume decorrelation observed over the Greenland Ice Sheet[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38（6）：2571-2583

[42]Luzi G，Noferini L，Mecatti D，et al.Using a ground-based SAR interferometer and a terrestrial laser scanner to monitor a snow-covered slope：results from an experimental data collection in Tyrol （Austria）[J].IEEE Transactions on Geoscience and Remote Sensing，2009，47（2）：382-393

[43]Rott H，Nagler T，Scheiber R.Snow mass retrieval by means of SAR interferometry[C]∥3rd FRINGE Workshop，European Space Agency，Earth Observation.2003：1-6

[44]Reigber A，Moreira A.First demonstration of airborne SAR tomography using multibaseline L-band data[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38（5）：2142-2152

[45]Rekioua B，Davy M，F(xiàn)erro-Famil L.Snowpack characterization using SAR tomography：experimental results of the AlpSAR campaign[C]∥European Radar Conference （EuRAD），2015，DOI：10.1109/EuRAD.2015.7346230

[46]Rekioua B，Davy M，F(xiàn)erro-Famil L，et al.Snowpack permittivity profile retrieval from tomographic SAR data[J].ComptesRendus Physique，2017，18：57-65

[47]Tebaldini S，F(xiàn)erro-Famil L.High resolution three-dimensional imaging of a snowpack from ground-based SAR data acquired at X and Ku band[C]∥IEEE International Geoscience and Remote Sensing Symposium（IGARSS），2013：77-80

[48]Ferro-Famil L，Tebaldini S.Comparison of parametric and non-parametric approaches for the full-rank polarimetric SAR tomography of volumetric environments[C]∥IEEE International Geoscience and Remote Sensing Symposium（IGARSS），2013，DOI：10.1109/IGARSS.2013.6723295

[49]Banda F，Dall J，Tebaldini S.Single and multipolarimetric P-band SAR tomography of subsurface ice structure[J].IEEE Transactions on Geoscience and Remote Sensing，2016，54（5）：2832-2845

[50]Tebaldini S，Nagler T，Rott H，et al.Imaging the internal structure of an alpine glacier via L-band airborne SAR tomography[J].IEEE Transactions on Geoscience and Remote Sensing，2016，54（12）：7197-7209

Overview of the snow parameters inversion from synthetic aperture radar

LI Zhen1TIAN Bangsen1ZHANG Ping1ZHAO Changjun1SUN Guangde1LIU Chang1

1Aerospace Information Research Institute，Chinese Academy of Sciences，Beijing100094

AbstractSnowpack is an important factor affecting the global energy water cycle，such as radiation balance and water resources，while remote sensing is a key method to monitor the temporal and spatial distribution characteristics of snowpack and its changing trend.Synthetic aperture radar （SAR） has become an important research means in snow remote sensing because of its all-day and all-weather observation capability.In this paper，the theory and technology of snow parameter inversion methods are summarized and analyzed from aspects of the characteristics of SAR and snowpack，snow scattering models，snow parameter inversion based on SAR intensity and phase information.And the current research states as well as the problematic issues in snow inversion using SAR data are pointed out，and the future research direction is discussed.

Key wordssynthetic aperture radar;snowpack;parameters;models;inversions;progress

收稿日期2019-10-22

資助項目科技基礎資源調(diào)查專項資助（2017FY100502）;中國科學院空天信息創(chuàng)新研究院重點部署項目（Y950930Z2F）

作者簡介李震，男，博士，研究員，主要從事冰雪遙感研究.lizhen@radi.ac.cn

1中國科學院空天信息創(chuàng)新研究院，北京，100094