北麓河流域不同植被覆蓋條件下土壤凍融過程中地表土壤含水量變化

胡越中， 王廣軍， 杜海波， 梁四海， 羅銀飛， 董高峰， 彭紅明

（1.中國地質大學（北京）土地科學技術學院，北京 100083；2.內蒙古煤田地質局勘測隊，內蒙古 呼和浩特 010010；3.中國地質大學（北京）水資源與環境學院，北京 100083；4.青海省環境地質勘查局，青海 西寧 810007；5.青海省環境地質重點實驗室，青海 西寧 810007；6.青海省地質環境保護與災害防治工程技術研究中心，青海 西寧 810007）

0 引言

地表土壤水分是大氣降水、植物含水和地下水相互轉化的紐帶，是植物生長和生存的決定性要素之一［1］，已成為農業、水文、氣候、生態等研究領域的重要基礎研究內容［2-4］。世界上三分之一的國家和地區，包括我國青藏高原的絕大部分［5］，處于干旱或半干旱地帶，水分的匱乏直接制約了當地的生態環境和經濟發展。因此世界各國十分關注土壤水分的研究，并為其投入頗多［6］，如2009年歐空局（ESA）發射的SMOS與2015年美國宇航局（NASA）發射的SMAP均為專注于全球土壤水分監測的衛星。

自20 世紀80 年代，微波遙感技術開始廣泛應用于區域土壤含水量監測研究［7-8］。由于土壤介電常數和土壤粗糙度影響土壤的后向散射系數，研究者根據地表電磁波傳播理論，先后建立起描述后向散射系數與地表物理與幾何參數之間的理論模型，為土壤含水量的反演提供了理論基礎，如Kirchhoff Approximation（基爾霍夫近似），其又根據土壤的粗糙程度分為POM（物理光學模型），GOM（幾何光學模型）兩種［9］、SPM（小擾動）［10］、IEM（積分方程）模型［11］等。同時一些學者基于理論模型，結合陸基后向散射計的實測數據，發展出各種適用于主動微波遙感的經驗、半經驗模型，如Oh 模型［12］，Dubois 模型［13］等。這類模型大多針對裸露地表而言，除土壤紋理結構外，制約其反演精度的重要因素是地表粗糙度［14-16］，通過正確選擇雷達參數或多頻率、多極化、多入射角數據，可以在一定程度上降低或削弱地表粗糙度的影響；另外，當地表被植被覆蓋時，其反演的土壤含水量往往會偏低，使之無法直接應用于有植被覆蓋的復雜地表［17-18］。因此在植被覆蓋區反演土壤含水量往往需要借助光學遙感提供的參數去除植被層對后向散射系數的影響，如MIMICS（植被微波散射模型）模型［19］，水云模型等［20］。

但總體而言，相較于傳統的點尺度物理觀測、面尺度的光學或被動微波遙感數據土壤含水量反演方法，主動微波遙感技術，特別是合成孔徑雷達（SAR）技術，以其宏觀、高效、全天候、高分辨率和穿透性強等特點已經被廣泛地應用于農業、水文和生態等領域［21-25］。

北麓河流域位于青藏高原腹地多年凍土區，陸面生態系統以高寒植被為主。已有的研究結果表明，青藏高原土壤凍融對土壤水熱變化影響較大［26］，凍結過程有利于土壤維持其水分，而土壤含水量反過來又影響土壤凍融過程及熱量分布［27］，進而影響植被的發育程度［28］。在北麓河左冒西孔曲小流域內的定點觀測結果表明，土壤凍融過程、土壤含水量和植被生長呈現出復雜的關系［27］。多年凍土活動層開始凍結和消融時間隨著植被蓋度的減少不斷提前，且隨著植被蓋度減小，活動層含水量變化速率增大，植被起到抑制土壤含水量變化速率的作用［29］。李元壽等［30］的研究結果進一步表明，北麓河左冒西孔曲小流域內多年凍土活動層淺層土壤凍結和融化對植被覆蓋度的響應程度較強，接近深層土壤凍結和融化對植被覆蓋度的響應程度降低，且伴隨草地退化、植被覆蓋度降低，土壤含水量對植被覆蓋響應越加強烈，植被覆蓋度越高，土壤含水量響應越是滯后且土壤含水量變化越加平緩。Wang 等［31］在忽略植被作用的前提下，利用2014—2016 年TerraSAR-X 多角度時間序列數據去除地表粗糙度的影響，基于雙組分土壤含水量檢索模型完成了北麓河流域內小范圍（＜100 km2）土壤含水量反演，但文中沒有對年內土壤含水量變化情況進行進一步分析。北麓河多年凍土區活動層凍融過程中，植被覆蓋和土壤含水量的響應仍有待于進一步深入研究。

2014年4月3日，ESA 哥白尼計劃（Global Monitoring for Environment and Security）中的地球觀測衛星Sentinel-1A 衛星發射升空，并于2017 年4 月對外共享數據，為主動微波反演土壤含水量提供了新的數據支持。針對Sentinel-1A 數據源，El Hajj 等［32］結合IEM 與水云模型建立帶有噪聲干擾的模擬C波段SAR 數據庫，并利用其進行神經網絡的訓練與驗證。文中嘗試多種神經網絡的輸入數據組合，充分研究了土壤含水量與C 波段SAR 數據和植被指數的相關性。結果表明單獨使用VV 極化或協同使用VV 與VH 極化反演土壤含水量均能取得較好的效果；Bao 等［33］基于水云模型發展出一種土壤含水量與雷達后向散射系數和植被指數關系的半經驗模型，并通過在英國與西班牙地區的實測數據進行擬合來獲得模型所需的參數。結果表明該模型可應用于復雜的植被覆蓋地表。

本文結合主動微波遙感與光學遙感兩種數據源，在改進的“水云模型”去除植被覆蓋影響的前提下，協同利用Dubois 模型與Oh 模型以消除難以獲取的土壤參數，利用模型正向模擬建立數據庫并通過代價函數最小化的方式反演土壤含水量。通過分析時序反演結果，從面尺度研究北麓河流域2018年年內2 037.94 km2土壤凍融過程、土壤含水量和植被覆蓋的關系。

1 研究區域及數據

1.1 研究區域概況

北麓河流域位于青藏高原腹地長江源區，地理位置如圖1 所示，坐標范圍為：92°54′～92°56′ E，34°48′～35°00′ N。研究區海拔為4 510～4 680 m，多年凍土極為發育，活動層厚度一般為2～3 m，地表土層以細粒土（黏土、亞砂土）為主，局部分布少量粗顆粒土（角礫土、碎石土）［34］，屬青藏高原可可西里區沖、洪積高平原地貌［35］。

圖1 北麓河流域位置示意Fig.1 Overview of Beiluhe basin

研究區植被生長發育情況較好，主要的植被類型為高寒草甸，主要植物種為矮嵩草、短穗兔兒草、藏嵩草等。在高寒荒漠區間或分布著火絨草、墊狀點地梅等。該地區屬于亞寒帶半干旱氣候，年均降水量在290～300 mm 之間，集中在6—8 月，占全年80% 以上；年最高氣溫為19.2 ℃，最低氣溫為-27.9 ℃，年均氣溫為-5.0～-3.8 ℃［35］，土壤凍結期為9月—次年4月［34］。

1.2 數據源和預處理

Sentinel-1A衛星處于高度為693 km的太陽同步軌道，軌道傾角為98.18°，軌道周期約為96分鐘，重訪周期為12天。Sentinel-1A衛星載有C波段合成孔徑雷達（SAR），使用TOPSAR的掃描方式，可較好的解決ScanSAR 幅度不勻調制的問題，保證了掃描區內的影像質量均勻［36］。其中心頻率為5.405 GHz，極化方式為VV+VH 與HH+HV，入射角度為20°～45°，輻射分辨率為1 dB·（3σ）-1，最大噪聲等效散射系數（Noise Equivalent Sigma Zero）為-22 dB。包括4 種工作模式：條帶模式（Strip Map Mode，SM），超寬幅模式（Extra Wide Swath，EW），寬幅干涉模式（Interferometric Wide Swath，IW）和波模式（WaveMode，WV）［37］。

本文使用的SAR 數據為IW 模式的Level-1 級地距多視GRD 產品。極化方式為VV+VH，空間分辨率為5 m×20 m，幅寬為240 km，重訪周期12 天，成像時間為2018年1月3日—2018年12月29日，共29 景。為抑制SAR 影像斑點噪聲，利用GammaMAP［38］濾波對影像進行去噪處理，濾波窗口大小為3×3；再使用SARspace 軟件進行輻射定標，地理編碼，重采樣與裁剪，最終獲得研究區空間分辨率為10 m的后向散射系數影像。

此外，采用Landsat-8衛星上搭載的陸地成像儀（OLI）數據反演植被覆蓋度，以消除植被冠層散射對土壤含水量反演結果的影響。Landsat-8 OLI 數據包括9 個波段，波長范圍為0.43～12.51 μm，重訪周期為16天（成像時間與Sentinel-1A最近影像最多相差8 天），空間分辨率為30 m，包括一個15 m 的全色波段。所使用的L1TP 產品已經經過系統輻射校正、幾何校正和地形校正，因此只需要進行輻射定標，將原始DN 值轉換為大氣外層表面反射率，之后利用ENVI軟件提供的FLAASH 大氣校正模塊消除大氣影響，并重采樣為10 m 分辨率的影像與Sentinel-1A 影像對應。使用Landsat-8 紅光（band4）和近紅外（band5）波段計算得到歸一化植被指數（NDVI）并應用像元二分模型反演植被覆蓋度［39］。此外，考慮到沒有實測葉面積指數（LAI）數據，無法由NDVI通過線性回歸求解LAI，以及采用其它經驗模型由于應用區域不同會引入新的誤差等問題，本文在前述植被覆蓋度反演的基礎上，借鑒陳麗等［40］的研究方法，采用多次散射過程冠層模型求解迭代的方法來逐步逼近精確的LAI 數值，以作為后續水云模型的輸入參數之一。

為獲取研究區表層土壤含水量數據，于2019年8月27日對北麓河流域進行了星地近似同步實地勘測（圖1），當日天氣陰，氣溫約為4～19 ℃，表層土壤未凍結。Sentinel-1A 衛星過境時間為北京時間2019 年8 月27 日6 時52 分，實地采樣時間為當日08：00—16：00?？紤]到地形與交通條件，沿109 國道選擇地勢平坦，植被長勢均勻的地方布設土壤含水量采樣點（圖1）。采用土壤含水量記錄儀（Uni1000土壤含水量記錄儀）獲取了地表0～8 cm深度的土壤含水量實測數據。采集樣方4個角點和中心點共5 個土壤含水量數據，通過距離加權平均計算得到與臨近影像像元地理位置相同的土壤含水量。依照上述方法，共采集16個樣方的土壤含水量均值數據。

2 研究方法

2.1 改進的水云模型

水云模型是一種適用于低矮植被的半經驗后向散射模型，已被證明在C、X 波段和20°～40°入射角的范圍內模擬精度較好［41］。本文在考慮單次散射的前提下，將給定像元（像元后向散射系數為植被和裸土混合，即混合像元）的植被覆蓋度引入水云模型以進行改進，使混合像元內總后向散射系數分解為植被覆蓋區地表貢獻和無植被覆蓋地表裸土地表貢獻，則雷達后向散射分量可以寫成：

式中：為給定混合像元總后向散射系數，上標pp為極化方式（本文為同向極化VV 或交叉極化VH）；f為植被覆蓋度。為更好描述冠層的散射特性，Prevot等［42］建議使用葉面積指數LAI代替植被含水量，則為植被冠層（水云層）后向散射系數，τ2為雙向植被透射率或衰減因子，分別由下式［42-43］給出：

式中：a，b為經驗系數，因研究區內土壤種類較多，難以獲取空間分辨率足夠高的土壤類型數據，因此采用Bindlish 等［44］的研究成果，分別取a，b為allland類型0.0012和0.091；θ為雷達入射角。

則根據式1，裸土地表或植被下墊面后向散射系數可由下式得到：

依照式（2）、（3）、（4），利用Sentinel-1A 影像處理得到的VV，VH 極化后向散射系數值以及與其時間最接近的Landsat-8 影像處理得到的植被覆蓋度和葉面積指數，計算得到土壤表層的后向散射系數值。由此，土壤含水量的反演轉化為去除植被覆蓋影響后土壤表層后向散射系數刻畫的問題，可以使用針對裸露地表反演土壤含水量的模型代替水云模型中的土壤表層后向散射系數函數。目前常用的經驗、半經驗模型有Dubois 模型，Oh 模型，Zribi模型等［45］。本文在前述考慮單次散射和混合像元求得的裸土地表或植被下墊面后向散射系數的基礎上，綜合Dubois 模型和Oh 模型，進行土壤含水量的反演。

2.2 模型的綜合使用

Dubois 模型是基于地面實測數據建立的經驗模型。在1.5～11 GHz 頻率范圍，30°～65°入射角范圍內，Dubois模型可以較好地模擬VV與HH后向散射。但是在地表較為粗糙［ks＞2.5，k為自由空間波束，s為均方根高度（cm）］，土壤水分較高時（ms＞35%，ms為土壤體積含水量）Dubois 模型將不再適用。由Dubois經驗模型［13］：

式中：ε為土壤介電常數；λ為波長（cm）。則可推導出：

Oh 提出了同向極化和交叉極化的后向散射系數經驗模型［12］。模型的適用的條件為：4%＜ms＜29.1%，0.13＜ks＜6.98，10°≤θ≤70°。

將式（6）代入Oh模型：

則可去除地表粗糙度的影響，式（7）中，ms為土壤體積含水量。進一步結合TOPP公式［46］：

可以建立起土壤體積含水量ms和土壤同向極化VV和交叉極化VH的關系模型，因公式過于復雜，為簡化起見，本文表示為：

本文利用查找表法［47］解決式（9）ms反演問題。對于某一像元，利用式（4）去除植被影響后，根據，ms與的關系模型對進行正向模擬（為真實值，ms范圍設定為0.01～0.99 cm3·cm-3），從而建立起查找表。并利用與的平方誤差代價函數S查找結果，即使S最小化時相對應的ms最優解作為最終的反演結果。S表達如下：

式中：為觀測VH 極化后向散射系數；為使用關系模型模擬的VH極化后向散射系數。

2.3 反演流程

本文反演土壤含水量的流程如圖2所示：

圖2 反演流程Fig.2 The inversion process

（1）首先選擇Sentinel-1A 影像及與其時間最為接近的Landsat-8 影像，并對其進行預處理。利用Landsat-8 影像計算得到研究區的植被覆蓋度及葉面積指數。

（2）逐像元進行判斷，若其植被覆蓋度大于0，則利用改進水云模型進行處理，最終獲得VV 與VH極化土壤表層的后向散射系數。

（3）綜合Dubois、Oh、Topp模型，消除中間量ks，建立ms、和三者間的關系式。并利用查找表法對每一像元進行反演，最終獲得研究區土壤含水量影像。

3 結果與分析

3.1 精度驗證

圖3為實測土壤含水量與反演土壤含水量的散點圖：實線為回歸方程，虛線為1∶1 直線，反演土壤水含量誤差棒為1.96 倍標準差（95%置信區間），實測土壤含水量誤差棒為儀器誤差；回歸方程調整決定系數（AdjustedR-square，Adjusted-R2）為0.6848，均方根誤差（Root Mean Square Error，RMSE）為0.039 cm3·cm-3。當土壤含水量較高時，會出現明顯低估的現象，從而造成擬合曲線斜率小于1。Dubois模型與Oh模型的最佳適用范圍分別為：土壤含水量在0～0.35 cm3·cm-3與0.04～0.291 cm3·cm-3以內。因此可歸因于土壤含水量超出最佳適用范圍，從而造成的低估現象。此外，在本文中多種模型運用過程中的累積誤差也是造成實測土壤含水量與反演土壤含水量存在差異的原因。

圖3 實測與反演土壤含水量散點圖Fig.3 Scatter plot of soil moisture content from measurements and inversions

3.2 土壤含水量年度變化分析

整體上看，北麓河流域1 月上旬（1 月3 日）—3月下旬（3 月28 日）土壤含水量基本無變化；從4 月初（4 月9 日）開始緩慢增加，5 月初（5 月3 日）—6 月中旬（6 月20 日），土壤含水量開始迅速增加并達到最高；6 月中旬后（6 月20 日），除個別地區略有波動外，北麓河流域土壤含水量在在整個夏季保持較穩定的狀態；從9 月上旬（9 月12 日）開始，土壤含水量開始下降，經過兩個月的下降期至11 月中上旬（11月11日）達到穩定，土壤含水量不再有明顯的變化，并一直持續到年末（12月29日）（圖4）。

圖4 2018時序土壤含水量分布Fig.4 Time series of soil moisture content in 2018

北麓河地區5 月中旬開始進入夏季融化期［26］，且6 月中旬以后，北麓河地區的降水類型主要以降雨為主［30］，因此，降雨和升溫導致的活動層融化是6月土壤含水量激增的重要原因；夏季土壤含水量的波動變化則主要是由于降雨量變化造成的；而9 月—11 月的土壤含水量下降則是由于降水減少且主要為固態降水造成。

但需要注意的是，北麓河流域部分地區的土壤含水量時序變化不盡相同，這可能與植被覆蓋或土壤質地有關，因此有必要根據植被覆蓋和土壤質地等選取典型區進行對比分析。

為弱化地形對雷達后向散射系數的影響，且考慮到研究區植被覆蓋普遍較低（夏季最高植被蓋度普遍在30%以下，個別水熱及土壤養分較好區域夏季最高植被蓋度在90%左右）選取地形平坦、植被覆蓋度能進行顯著區分及土壤質地差異明顯的4個區域作為典型研究區（圖1），分析研究區土壤凍融過程、土壤含水量和植被覆蓋的關系。為減弱土壤含水量空間變異性的影響，每個區域大小為3×3像元。

圖5為各典型區的年內土壤含水量變化。對比各典型區土壤含水量變化折線，可以發現：

圖5 典型區土壤含水量變化折線圖Fig.5 Temporal variations in soil moisture content in the four typical areas

（1）相較于其它典型區，90%植被蓋度區多年凍土活動層在11 月—4 月期間含水量基本無變化，保持穩定，但隨著植被覆蓋度的降低，在此期間，其它典型區（沙化地表除外）的土壤含水量則波動較大。

分析原因，這可能與地表均質性或土壤質地對雷達后向散射系數的影響有關。90%植被蓋度區（主要是高寒草甸）土壤地表土質為亞砂土、亞黏土，且地表腐殖質較厚［5］，地表均質性較好，進而減少了對雷達后向散射系數的干擾，且在此期間，降水主要以固態水形式，活動層凍結后未凍水分很小。

活動層與大氣的水汽交換主要在土壤表層以升華等形式進行，水汽交換量很?。?6］，而植被層又具有保溫作用［25］，從而使高植被覆蓋區土壤含水量保持較好的穩定性。但低植被覆蓋區礫石較多且分布不均，或裸土地表景觀破碎化較大（圖6），而對雷達后向散射系數影響較大，進而使反演的多年凍土活動層土壤含水量波動性較大。沙化地表由于主要是細砂、粉砂組成，均質性較好，故沙化地表土壤含水量在此期間也保持了較好的穩定性。

圖6 低植被覆蓋區地表實地照片Fig.6 Photos of the low vegetation cover areas

90%植被蓋度的典型區土壤含水量從4 月9 日開始穩步緩慢增加至5 月3 日，推測該時間段內土壤升溫較慢［26］?？紤]到衛星重訪周期，該區多年凍土活動層消融時間應在4 月9 日之后，這與沱沱河埋深4 cm 活動層消融時間為4 月15 日的野外觀測結果［27］吻合（沱沱河距北麓河直線距離＜90 km），進一步驗證了本文觀測結果的準確性。同理，根據雷達反演的土壤含水量年內曲線，推測該區凍結期應在10月18日后。

（2）雖然30%植被蓋度的低植被覆蓋區11月—3 月期間波動較大，但從圖4 中仍然可以看出，3 月28 日—次年4 月21 日，土壤含水量有明顯增加趨勢。同因衛星重訪周期，推測30%植被蓋度區多年凍土活動層消融時間應在3 月28 日后。同時根據土壤含水量變化曲線，該區的凍結期應早于90%植被蓋度的典型區。

（3）裸土地表從3 月4 日開始，土壤含水量開始呈現波動性上升趨勢，因此推測裸土地表活動層的開始融化時間應在3 月4 日后。對比后續土壤含水量推測在3 月16 日后研究區有一次降溫過程，致使3 月28 日土壤含水量再次降低。而同期90%植被蓋度典型區含水量無變化，可能是由于植被保溫效果，導致氣溫波動對土壤含水量的影響較少所致。

（4）沙化地表在6 月8 日后土壤含水量才開始出現明顯攀升，10 月18 日后下降到低值點，推測埋深在8 cm 以上的沙化地表受凍融作用影響很小，其含水量主要受降水影響。

（5）總體上看，植被覆蓋對活動層凍融過程具有明顯的遲滯作用，即植被覆蓋度越高，活動層凍結和消融時間越滯后，這與劉光生等［29］和李元壽等［30］的研究結論一致。

4 結論與展望

本文基于2018 年Sentinel-1A 時序數據，輔以Landsat-8 光學遙感數據，利用改進的“水-云”模型，綜合Dubois 模型、Oh 模型和TOPP 公式，完成了北麓河流域2018年土壤含水量的反演，并基于實測數據驗證精度。通過對年內土壤含水量時序反演結果的分析，得到了土壤凍融過程、土壤含水量和植被覆蓋的關系。主要結論如下：

（1）本文提出的在改進的“水云”模型基礎上，綜合Dubois 模型、Oh 模型和TOPP 公式進行土壤含水量反演，有效去除了植被覆蓋和地表粗糙度的影響。16 個點的同步實測結果表明：回歸方程為Adjusted-R2為0.6848，RMSE為0.039 cm3·cm-3。

（2）地表均質性或土壤質地對雷達后向散射系數的影響較大。11 月—次年4 月期間，因高植被覆蓋度典型區地表均質性較好，加之降水主要以固態水形式，活動層與大氣的水汽交換量小，以及植被的保溫作用，含水量基本無變化，保持穩定，但隨著植被覆蓋度的降低，在此期間，其它典型區（沙化地表除外）的土壤含水量則波動較大。

（3）植被覆蓋對活動層凍融過程具有明顯的遲滯作用，即植被覆蓋度越高，活動層凍結和消融時間越滯后。推測裸地、30%植被蓋度和90%植被覆度的典型區2018 年消融開始時間分別在3 月4 日后、3月28日后和4月9日后。90%植被蓋度典型區凍結期應在10 月18 日后，而裸地和30%植被蓋度典型區則早于90%植被蓋度典型區。

（4）本文研究結果中北麓河地區年內土壤含水量變化趨勢與前人定點觀測所得到的趨勢基本吻合，表明基于Sentinel-1A 年度時序數據，從流域面尺度研究土壤凍融過程、土壤含水量和植被覆蓋的關系是可行的。

多年凍土活動層凍融期土壤水分以未凍水與冰兩種形式存在，散射機理較為復雜；同時因研究區凍融氣候環境惡劣，未能采集該時期的研究區土壤含水量數據，因此凍融土壤含水量反演的方法與精度有待進一步研究和驗證。另外，10 月30 日Sentinel-1A數據缺失，導致本文無法給出上述典型區的準確凍結日期。

土壤含水量與降水等密切相關，而北麓河流域又地處青藏高原腹地，氣候多變，多數情況下研究區內各地夏季降雨不同步，要準確地探究北麓河流域多年凍土區，特別是夏季土壤含水量與植被覆蓋的關系，仍是后續要持續探究解決的難題之一。