半干旱草原型流域表層土壤飽和導水率傳遞函數及遙感反演研究*

黎明揚 劉廷璽，2? 羅艷云，2 段利民，2 張俊怡 周亞軍 Scharaw Buren

（1 內蒙古農業大學水利與土木建筑工程學院，呼和浩特 010018）

（2 內蒙古自治區水資源保護與利用重點實驗室，呼和浩特 010018）

（3 Application Center for System Technologies，Fraunhofer IOSB，llmenau 98693，Germany）

土壤是植被生活的基質，不僅為植物生長提供必需的養分，同時也在整個生態水文循環中起到傳輸和貯存水分、物質和能量的重要作用。土壤傳遞函數（Pedo-Transfer Functions，簡稱PTFs）是土壤科學中被廣泛應用的一種方法，可利用有機質含量、容重、粒徑組成等土壤基本特性指標，間接計算出土壤飽和導水率（KS）、水分特征曲線等常用的參數和曲線［1］。國內外學者通過大量室內和野外實驗對KS和PTFs進行了廣泛的研究［2-5］，Aimrun和Amin［6］在Sempadan水稻種植區域的研究表明通過以自然對數形式輸入土壤干容重、顆粒百分比、有機質和幾何平均粒徑等因變量擬合出來的KS是最佳的PTFs形式；Xu等［7］將土壤光譜信息作為附加輸入，使用人工神經網絡模型（ANN）重新評估了PTFs在估算水力傳導率方面的優越性。通過各種方法建立的PTFs雖然能夠較為準確地預測試驗點上的KS，但要估算一個流域或區域不同空間點的KS，如果去實地測試每個點上的土壤基本特性指標，再通過PTFs估算各點的KS，這無疑是費時、費力、不科學的。合成孔徑雷達（SAR）問世于20世紀60年代，是目前反演和解譯流域或區域土壤特征參數方面最為有效的技術［8-9］。Radarsat-2是一款搭載C波段傳感器的全天時、全天候的空間微波遙感成像雷達，其作用距離遠，抗干擾性能好，對云層和地面植被有一定穿透能力，該雷達被大量用于表層土壤體積含水率的反演中［10-12］，但目前利用其來反演土壤KS的研究還不多，具有重要的學術意義。

干旱半干旱地區草原生態系統十分脆弱，無論是單純采用大面積土壤采樣的PTFs還是直接利用全極化SAR數據反演的方法來預測，均無法在面尺度上準確地體現區域空間變化，因此研究一種通過PTFs與SAR相結合且對生態干擾破壞最小的新方法變得至關重要。本文采用野外原位采樣與同時期Radarsat-2衛星數據相結合的方法，首先擬合了研究區表層土壤KS的PTFs，之后結合全極化SAR數據與表層土壤參數均值建立函數關系，最終預測出雷達數據研究區范圍內表層土壤的KS。

1 材料與方法

1.1 研究區概況

選擇內蒙古高原半干旱的錫林河流域為研究區，其位于內蒙古自治區中部錫林郭勒盟，發源于赤峰市克什克騰旗寶爾圖山，流經錫林郭勒盟阿巴嘎旗，在貝力克牧場轉向西北流經錫林浩特市，最終注入查干諾爾沼澤地自然消失［13］。開展野外采樣的區域位于錫林河支流浩勒圖郭勒與干流錫林高勒河交匯處以上流域，43°24″～44°4″N，116°17″～117°15″E，面積約為1 852 km2，區域地勢三面環山，植被多為天然牧草，以茅針和羊草最為常見。研究區屬中溫帶半干旱大陸性季風氣候，多年平均降水量為266.8 mm，其中，6—8月降水量占年總降水量的50%以上［14-15］。

圖1 研究區位置、采樣點分布及雷達影像Fig. 1 Location and sampling points in and radar image of the study area

1.2 試驗設計

首先，利用全國1∶100萬土壤類型圖（http://westdc.westgis.ac.cn）和全國DEM高程圖（http://www.gscloud.cn）細致刻畫出研究區五種土壤類型分區圖和等高線分布圖，在此基礎上，按照每種土壤類型所占研究區的面積比例，選擇多個具有代表性的樣地，在雷達影像范圍內的采樣點共布設42個，其中分為32個建模采樣點，10個驗證采樣點，土壤采樣點垂直于等高線布設，包含研究區90%以上的海拔高度、坡度、坡向及不同地上生物生長狀況。使用手持GPS記錄采樣點坐標，在保證選取土壤質地均一且完整的前提下，去除地表浮土及植被后，每個樣點處挖50 cm深的剖面，從表層向下取樣，每隔10 cm一層，每層采用100 cm3的環刀重復采集土樣2個和自封袋重復樣3個，分別用于測定水力參數、干容重、土壤粒徑分布及其他物理化學特性。反演土壤參數和驗證模型時均使用0～30 cm表層土壤數據。本次試驗于2017年7月23日至30日共采集原狀環刀土樣420個，自封袋土樣630個。

1.3 土壤理化指標測定

土壤粒徑利用HELOS & RODOS激光粒度分析儀進行干法測定，分級標準采用美國制：黏粒（＜0.002 mm），粉粒（0.05～0.002 mm），砂粒（2～0.05 mm）；土壤干容重采用環刀法測定；土壤有機質采用濃硫酸-重鉻酸鉀外加熱法測定；飽和含水率采用環刀底面纏紗布，由底部浸水測定，每隔6 h稱重一次，待兩次稱重誤差小于2%時，取平均值為其試驗值；土壤飽和導水率利用定水頭下的馬氏瓶滲透儀進行測定。

1.4 雷達影像的處理

為保證野外采樣時間與雷達成像時間對應一致，提前購置2017年7月26日RADARSAT-2精細四極化SLC（single look complex）格式雷達影像一景，相幅為25 km×25 km，分辨率為8 m。使用ENVI軟件的SARscape模塊對原始SLC影像進行多視（Multilooking）、Refined-Lee濾波（Refined-Lee Filtering）、地理編碼和輻射定標（Geocoding and Radiometric Calibration）等處理，在Google影像上選取GCP（Ground Control Point）對影像進行幾何配準后，得到標準四極化后向散射系數影像圖，利用Output ROIs to ASCII功能將所需的采樣點后向散射系數輸出成文本文檔以供使用。

由于試驗前期極度干旱，連續4個月總有效降水不超過20 mm，導致試驗期間研究區平均植被蓋度僅為46%，平均植株高度低于10 cm，且植株長勢差，各采樣點平均地表生物量為58 g·m-2，故研究認為植被對雷達返回信息的影響可以忽略，未進行植被的去極化處理。具體研究區位置及高程信息、土壤類型分區、采樣布設及處理后雷達影像如圖1所示。研究采用同極化差SAR數據計算組合地表粗糙度參數ZS［8］，同極化差Δσo=σovv-σohh與ZS的關系可表示為

式中，A(θ)和B(θ)的表達式如式(2)及式(3)所示，θ為本地入射角，為了簡化計算，統一取雷達入射角40.917°為本地入射角。

1.5 研究方法

將研究區32個建模采樣點0～30 cm每隔10 cm一層的表層土壤數據分別代入Cosby［16］、Saxton［17］、Wosten［18］及非線性多元經驗回歸模型［19］，建立錫林河流域表層土壤KS的PTFs，通過32個建模采樣點表層土壤參數均值對PTFs進行驗證，并與雷達后向散射系數建立多元線性模型，利用10個未參與建模的驗證采樣點雷達后向散射系數反演表層土壤參數，代入PTFs中進行驗證，最終選擇模擬效果最佳的PTF進行區域表層土壤KS的預測。研究采用判定系數R2、均方根誤差RMSE及 t 檢驗對模擬值進行精度分析。

2 結果與討論

2.1 土壤參數的統計特征

研究區主要被草甸地為主的大粒徑砂質土壤覆蓋，地貌多變，根據各采樣點所屬的土壤類型，將42個參與模型計算的采樣點實驗數據劃分為厚栗黃土、草甸沼澤土、荒漠風沙土、石灰性草甸砂土、淡黑土5種類型，并按照土層的不同深度進行分類統計，其結果如圖2所示。

研究區土壤砂粒含量極高，幾乎不存在黏粒，草甸沼澤土、厚栗黃土及石灰性草甸砂土砂粒含量最高，各層砂粒含量皆高于80%，荒漠風沙土最低，其平均含量仍達80.19%。在隨深度變化方面，5種土壤的表層粒徑變化均不明顯，這說明長期的降水會讓向下運移的水分不斷攜帶體積小、粒徑細的顆粒一同運動，最終導致土壤在運移較劇烈的表層達到粒徑分配方面的平衡。土壤較高的砂粒含量導致其保水能力差，水分下滲速度快，這使得研究區植被稀松，有機質含量較低，這與Zhang等［20］和Yu等［21］對半干旱區植被與土壤理化性質關系的研究結果近似。在地形平坦的區域，土壤有機質含量隨深度增加呈下降趨勢，而主要分布在地勢陡峭的厚栗黃土，其土壤有機質含量在20 cm土層有明顯增高，這是因為植物為了生長在畜水困難的坡地上，深層表土根系更加發達，擁有大量毛細根的20 cm土層經過常年的積累有機質含量有明顯提高。

圖2 5種土壤粒徑、容重及有機質含量隨深度變化規律Fig. 2 Variation of soil particle size, bulk density and organic matter content with depth relative to soil type

2.2 飽和導水率土壤傳遞函數的適用性

土壤的結構、物理化學性質直接或間接影響土壤KS，將建模采樣點表土各層數據輸入到Cosby、Saxton、Wosten三種傳統模型及多元經驗回歸模型（Empirical model）中建立研究區表層土壤KS的PTFs，代入建模采樣點0～30 cm表層土壤參數的均值進行模型的檢驗。各PTFs形式、建模及檢驗的模擬實測值對比如表1和圖3所示，圖3中層值與均值的線性結果分別位于各分圖的左上及右下。

表1 土壤飽和導水率的傳遞函數Table 1 Transfer functions of soil saturated hydraulic conductivity

注: 方程中，Ks為飽和導水率(m·d-1)，c1、c2、c3分別為土壤的黏粒、粉粒和砂粒含量(%)，ρ和c4分別為土壤容重(g·cm-3)和土壤有機質含量(g·kg-1)，θs為土壤飽和含水量(%)，dg為平均粒徑(mm)，σg為粒徑標準偏差(%) Note: In the equation, Ksstands for saturated hydraulic conductivity (m·d-1); c1for clay content (%); c2for silt content (%); c3for sand content (%); ρ for soil bulk density(g·cm-3); c4for soil organic matter content (g·kg-1); θsfor saturated water content (%); dgfor average particle size (mm); σgfor standard deviation in particle size (%)

從圖3中4種PTFs的模擬值與實測值擬合關系可以發現，Cosby模型效果最差，這與Mermoud和Xu［22］的研究結果相同；其余3種模型的擬合及檢驗結果均高于0.958，Saxton模型擬合效果最佳，建模精度可達0.985，RMSE為0.262，Wosten與多元經驗回歸模型次之；各模型置信度均大于95%，說明PTFs模型在多種土壤類型模擬中不存在顯著性差異并具有普遍適用性，這與前人對多種類型PTFs的研究結果相近［23-25］。

PTFs的不同形式側面說明了模擬結果的差異性，由于Cosby模型只涉及到土壤的顆粒組成，所以擬合結果偏差較大，其主要問題表現為當實測土壤砂粒含量大于80%且土壤容重超過1.6 g·cm-3時，模型KS的模擬值會出現較大幅度的偏差，這說明Cosby模型在擬合砂粒含量和土壤容重同時較高的土壤時，具有一定的局限性。

2.3 合成孔徑雷達對土壤參數的反演

圖3 飽和導水率土壤傳遞函數模擬值與實測值對比Fig. 3 Comparison between measured and PTFs-estimated values of soil transfer function of saturated hydraulic conductivity

輻射校準后的SAR圖像是目標特征的綜合反映，SAR能同時實現對地面目標的距離向和方位向的高分辨率成像，全極化SAR系統交替發射兩路相互正交的極化脈沖，并同時用兩路相互正交的極化通道接收信號，其獲得的類似光學照片的目標圖像與傳統光學遙感相比，具有一定的優越性［8］。使用位于雷達影像范圍內建模采樣點0～30 cm表層土壤參數的平均值與微波雷達RADARSAT-2四極化后向散射系數HH、HV、VH、VV及其兩種組合HH/VV、HV/VH以及組合地表粗糙度ZS建立多元線性回歸方程，以刻畫雷達影像研究區范圍內各土壤參數的分布情況。各土壤參數的結果及反演形式如表2及式(4)至式(9)所示。

表2 0～30 cm表層土壤參數反演結果Table 2 Inversion of the parameters of the 0～30 cm surface soil layer

通過表2可以看出，由于微波雷達對地物的高穿透性，各表層土壤參數與全極化SAR數據的擬合結果較為理想。土壤的粒徑分布也可以較好地表示，C1作為粒徑處于0.002～0.05 mm的粉粒含量，在反演中效果最差，本文通過反演效果較好的C1及C3推求獲得，這樣做既滿足粒徑分布的平衡要求，其結果也有一定的提升；θs與ρ有較好的線性關系，這與前人的研究結果一致［3,26-29］，可通過式(10)計算獲得。如式(9)所示，利用SAR數據也可以直接推求表層土壤的KS，相較前文中的四種模型，擬合效果僅優于Cosby模型，這說明在結合SAR數據對KS進行面尺度推廣時，使用PTFs可以獲得更精準的結果。

2.4 土壤傳遞函數精度檢驗及飽和導水率的預測

利用未參與模型運算的10個驗證采樣點表層土壤數據的平均值進行模型的驗證，結果如圖4所示。研究表明，Saxton模型的雷達預測檢驗效果最好，模擬值最為收斂且誤差最小，這說明Saxton模型更傾向于黏粒含量低的砂質土壤，這與Buccigrossi等［29］的結論相同，相較其他模型而言，Saxton模型中的參數并不是最多的，這說明參數多的模型不一定較簡單的模型效果好。直接利用全極化SAR數據進行KS預測的方法精度有待提高，其模擬值收斂性不高，容重高于1.7 g·cm-1的模擬值偏低，甚至出現負值；Cosby和Wosten模型的檢驗結果較差，當利用SAR估計的表層土壤容重偏高時，原本對容重變化不敏感的Cosby模型模擬效果會更差，Wosten模型包含了土壤有機質含量一次方及二次方項，當表層土壤有機質含量的預測出現問題時，KS模擬值便會受到一定程度的影響，這導致在建模階段表現尚好的Wosten模型在驗證階段失準；多元經驗回歸模型由于其計算了平均粒徑及粒徑的標準偏差，這可以在一定程度上修正SAR對表層土壤粒徑分布方面的預測偏差。

選擇效果最佳的Saxton模型，使用ENVI軟件的波段運算工具（band math）對雷達影像研究區范圍的表層土壤KS進行預測，并利用彩虹條從高到低表示。圖5為雷達影像研究區范圍內表層土壤KS的預測分布圖。

圖5反映了Saxton模型模擬的雷達影像研究區范圍內表層土壤KS的分布情況，可以看出，人類活動對表層土壤的KS影響較大。圓形的滴灌區顏色深而純，說明砂質土壤在經過農業改造后表層土壤的KS可以獲得大幅度降低，最低不足1 m·d-1，表層土壤較高的持水保水能力將更有利于植物的生長；灌區附近分布著少量村莊，磚瓦和混凝土影響了水分的下滲，導致該區域土壤KS雖大于滴灌區但仍下滲較慢。河流兩側山地地勢起伏大，為表面裸露、植物稀疏的巖漿巖山區，土壤中包含大量孔隙發達的火成巖，導致土壤下滲速度較快，最高可達25 m·d-1；靠近河流的河谷地區水分充足，水分下滲速度較緩，使得該區域生長了大量的水生植物，含水量較多的地表植物其葉片形狀和植株倒伏特性是隨機的，這會在一定程度上影響雷達對土壤參數反演的精確性，這與Millard等［30］和Liao等［31］在濕地中植被對RADARSAT-2反演土壤參數時產生較大影響的觀點一致。通過預測可以發現研究區表層土壤KS的分布規律，即：山區裸地＞沙丘沙地＞河間濕地＞村鎮建筑用地＞滴灌區。

圖4 雷達反演飽和導水率精度分析Fig. 4 Accuracy analysis of the inversion of saturated hydraulic conductivity by radar

圖5 雷達影像研究區范圍內Saxton模型表層土壤飽和導水率的遙感預測Fig. 5 Remote-sensing based prediction of saturated hydraulic conductivity of the soil surface layer with Saxton model within the radar image of the study area

3 結 論

在利用SAR大面積預測土壤KS時，使用PTFs預測的結果明顯優于直接利用SAR數據反演的結果；PTFs模型的擬合效果與其中包含參數的個數無關，Saxton模型更適用于砂粒含量較高、有機質含量較低的半干旱草原型流域。環境退化的草原表層土壤的KS普遍為4～8 m·d-1，自然條件下缺少植被、砂粒含量高是導水較快的主要原因，在灌區、村鎮等經過開墾或建設開發的地區，表層土壤的KS大幅降低，說明人類活動是其變化的重要影響因素。而利用SAR進行區域表層土壤KS的大面積預測技術尚不成熟，有待進一步驗證。