基于SVM的綠洲荒漠交錯帶土壤水分與地下水埋深反演

張鈞泳 丁建麗 譚 嬌

(1.新疆大學資源與環境科學學院， 烏魯木齊 830046；2.新疆大學智慧城市與環境建模自治區普通高校重點實驗室， 烏魯木齊 830046；3.新疆財經大學計算機科學與工程學院， 烏魯木齊 830032)

0 引言

地下水分布范圍廣泛，其可持續利用性強，在水資源中占有重要地位。淺層地下水作為地下水組分之一，是干旱植被生長環境的限制變量，地下水資源循環周期長，一旦賦存環境遭到破壞，再生能力減弱。隨著人類社會發展，長期不合理、無規則開采使地下水水位持續下降，地下水系統遭受破壞[1]。大多數干旱內陸河流域面臨嚴峻的大型灌溉農業用水[2-3]，地下水是水源補給的主要方式，故干旱區綠洲生態環境的調控作用顯得尤為重要[4]。合理的淺層地下水位能滿足干旱區綠洲植被的正常生長需求，這使得淺層地下水、植被與土壤之間的相互作用和植被群落分布特征與淺層地下水的相互響應研究尤為重要[5]。遙感技術具有宏觀動態快速監測的特點，避免了傳統檢測方法中時間長和精力消耗多的缺點，在監測地下水位應用中可以快速獲得地下水位，實現大面積監測[6]。隨著衛星遙感技術不斷進步，定量反演地表特征參數，結合相關實測數據，建立地下水相關模型，可快速獲得地下水相關信息[7-11]。

目前，遙感手段在大尺度區域反演地下水埋深研究中已經成為主流趨勢，但同時地下水遙感反演也是一個技術難題，地下水遙感監測主要是通過地下水與地表水[12]、地表溫度[13]、土壤水分[14]和植被[15]等相關性進行研究。目前，對地下水進行遙感探測有以下4種類型：①熱紅外遙感監測法，早期該方法利用熱紅外波段判斷地下水的存在，現今姜紅等[12]分析了新疆開都河流域周圍的綠洲，監測和分析該區域17年間的地下水位變化；付俊娥等[13]利用以MODIS遙感影像計算植被指數(NDVI)以及地表溫度(LST)數據，通過SVM定量估算了河西走廊疏勒河流域的地下水位；尹濤等[14]通過1 d中不同時間段DNVI、LST進行數學變換，反演了在植被生長期的地下水埋深情況。②環境因素遙感信息法，根據與地下水相關的環境因子與地下水間的關系，來監測和反演區域地下水狀況。田凱等[15]利用MODIS數據并計算NDVI以及條件植被覆蓋率(VCR)，與野外實測的地下水埋深相結合，發現黃河中上游地區的植被生長所適合的地下水埋深。③遙感信息定量反演分析法，該方法案例較多，如毛德發等[16]在北京市大興區消耗水用量以及本著地下水可持續利用的條件下，利用MODIS數據以及氣象降水徑流等數據，得出該區域綜合節水的有效方法。④微波遙感監測分析法，RODELL等[17]利用GRACE衛星數據，對密西西比河主要的子流域地下水儲量變化進行了估算。KOMAROV等[18]在水文因素、土壤介電常數和土壤對雷達波段的反射特性分析的基礎上,建立土壤水分和地下水之間的實驗方程。尹楠等[19]在不同極化模式下，研究了后向散射系數對壟行方向變化的敏感性及去除壟和地表粗糙度對土壤含水率反演的影響。綜上所述，利用微波遙感、土壤水分、植被指數和熱紅外等基于遙感信息定量監測地下水位已經取得了較大的進展，其次，干旱植被指數在土壤水分和地下水的監測中有著至關重要的作用[20-21]。

本文以新疆渭干河-庫車河綠洲(以下簡稱：渭-庫綠洲)為研究區，選擇Sentinel-1A影像和Landsat-OLI影像為數據源，通過微波和光學遙感數據以及土壤含水率、地下水埋深數據，結合植被以及土壤條件，利用支持向量機算法進行定量反演研究區土壤含水率以及地下水埋深，為監測綠洲的地下水埋深提供技術支持。

1 研究區概況和數據資料

1.1 研究區概況

渭-庫綠洲位于亞歐大陸腹地。本文所界定的野外實地調查采樣區間為81°26′～83°17′E，41°06′～41°56′N。從地貌條件看，北部地區海拔最高，其次為內部平原區，南部沙漠地區海拔最低，是一個典型的沖積扇平原。氣候干燥少雨，氣溫呈北低南高的特點。綠洲內部的光照和熱量資源十分充足，土壤類型豐富，主要農業經濟作物為小麥、玉米、棉花、石榴和無花果等。植被類型多為鹽生植被，呈片狀分布，多分布于綠洲外部、綠洲東北部荒漠帶及綠洲荒漠交錯帶[22]。土壤水分反演區域為渭-庫綠洲全區，如圖1a所示，地下水埋深反演區域為庫車綠洲的東北角，如圖1b所示。

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

1.2 實測數據

2015年3月21日—4月30日對研究區開展了一次野外調查，調查的內容主要是了解春灌之前稀疏植被情況以及采集土壤樣本，獲取土壤含水率、溫度、質地、容重、介電常數以及地下水埋深數據。根據實驗室多年數據的累積，選擇了38個具有代表性的采樣點，均勻分布在綠洲、綠洲荒漠交錯帶以及荒漠帶區域。其采用五點法進行土壤樣本的采集，使用W.E.T型傳感器、Hydra型土壤測試儀測量土壤含水率、溫度以及介電常數等數據，土壤電導率(EC)和總溶解固體(TDS)含量由inolab cond 7310型臺式電導率測試儀(德國WTW公司)測定，pH值由pH7310型臺式酸度計(德國WTW公司)進行測定，土壤容重使用環刀法進行測定。使用美國Onset公司生產的HOBO型自動記錄水位計(以下簡稱水位計)記錄地下水埋深數據，該儀器可根據實驗的要求設定水位計的記錄周期，水位計產生的誤差最大為1 cm，降低了實驗誤差。

1.3 遙感數據

Sentinel-1A是合成孔徑雷達(Synthetic aperture radar，SAR)的一種，其主要的特點是空間分辨率高，工作狀態不受云雨天氣影響和波動，重訪周期短，可全天時、全天候監測，具有單雙極化方式，較強的干涉能力等優點[23]，在全球范圍內可用于森林監測、農業監測、地震監測、陸地冰雪提取、洪水監測以及土壤水分、土壤濕度監測等[24-29]。本文Sentinel-1A數據來源于歐空局的哨兵科學數據中心(Sentinels Scientific Data Hub)，影像成像時間為2015年4月21日，為地距影像(Ground range detected，GRD)Level 1級產品，極化方式為VV，工作方式為干涉寬幅模式(Interferometric wide，IW)，分辨率為5 m×20 m，與太陽同步軌道，載波波段為C波段，工作頻率為5.044 GHz，軌道高度為693 m，重訪周期為12 d，入射角為20°～40°，幅寬240 km。

使用遙感數據為對地觀測衛星Landsat OLI的影像，獲取時間為2015年4月26日，空間分辨率為30 m，數據來源于http:∥glovis.usgs.gov/。

1.4 遙感影像預處理

利用歐空局(ESA)提供的Sentinel-1 Toolbox(S1TBX)軟件對微波圖像進行預處理：①輻射校正。利用S1TBX軟件提供的輻射校正工具完成輻射校正，得到影像后向散射系數。②噪聲處理。首先對圖像進行熱噪聲移除；其次，為了降低影像中的斑點噪聲，采用Refined Lee濾波進行處理，濾波器窗口設置為7像素×7像素。③幾何校正。先將處理好的雷達圖像轉換成后向散射系數，然后進行幾何校正。將處理之后的Sentinel-1A SAR數據輸出為ENVI數據格式，并將其重采樣到與Landsat 8數據具有相同的空間分辨率。

Landsat OLI需要進行輻射校正、大氣校正、幾何校正以及圖像裁剪等預處理操作，具體操作參照文獻[30]。

2 研究方法

2.1 不同溫度-植被指數構建

SANDHOLT等[31]發現不同植被指數與溫度指數均呈不規則簡化三角形分布，在此基礎上提出了溫度植被干旱指數的概念。Ts-VI由植被指數和地表溫度計算得到，其定義為

(1)

其中，對最大、最小地表溫度散點進行線性擬合，得到TVDI特征空間的干濕邊方程

Tsmax=a1+b1VI

(2)

Tsmin=a2+b2VI

(3)

式中TVDI——溫度植被指數

Ts——地表任意像元的地表溫度

VI——植被指數

Tsmin——相同植被指數值對應的最小地表溫度

Tsmax——相同植被指數值對應的最大地表溫度

a1、b1——干邊線性擬合方程的系數

a2、b2——濕邊線性擬合方程的系數

其中，Tsmin對應Ts-VI特征空間的濕邊；Tsmax對應Ts-VI特征空間的干邊。Ts-VI的閾值在0～1之間，當TVDI值逐漸向1靠近時，土壤干旱情況越來越嚴重；反之，當TVDI值越靠近0時，土壤濕度越高。因此Ts-VI與土壤濕度的相關性，在兩種極端的情況下可以反映干、濕情況。

為消除大氣散射的影響以及地表相鄰點反射光折射的影響，利用ENVI圖像處理軟件通過數字高程(DEM)數據提取坡度、坡向數據。其操作步驟參照文獻[32]，進而對LST數據進行地形校正。

2.2 土壤后向散射系數計算

為了得到比較準確的土壤表面后向散射系數，對有植被影響的后向散射系數進行處理，以消除植被所引入噪聲帶來的負面影響，從而獲取區域表層土壤真實后向散射系數。根據研究區當時實際情況，選擇適合本區域的Water-Cloud算法

(4)

(5)

γ2(θ)=exp(-2Bmvegsecθ)

(6)

γ2(θ)——農作物層的雙層衰減因子

mveg——植被含水量，kg/m3

θ——雷達波入射角

其中，A和B分別為依賴于植被類型的參數，可以通過多次迭代衰減或最小二乘法獲得。

估算土壤水分時應規避植被覆蓋影響，則需要應用Water-Cloud模型消除植被層在土壤水分后向散射中的干擾。研究學者對于A、B參數的合理取值及地域性分析進行了研究，本文采用周鵬等[33]所標定的符合渭-庫綠洲的最優解，即參數A=0.001 9，B=0.137。

在利用Water-Cloud模型進行有植被區的土壤水分定量估算時，其中一個至關重要的輸入參數為植被含水量(VWC)。研究區植被多為棉花及低矮的植被，為了獲取研究區的植被含水量，利用改進型歸一化植被水分指數NDMI[34]，并根據經驗方程[35]，則VWC計算式為

VWC=2.15NDMI+0.32

(7)

其中

NDMI=(ρNIR-ρMIR)/(ρNIR+ρMIR)

(8)

式中ρNIR——遙感影像近紅外波段反射率

ρMIR——遙感影像中紅外波段反射率

將式(7)代入式(5)與式(6)，運用Water-Cloud模型消除植被覆蓋影響從而得到裸土后向散射系數，表示為

(9)

γ2(θ)=exp(-2×0.137mvegsecθ)

(10)

(11)

通過式(5)和式(11)可以計算出研究區有無規避植被影響的土壤后向散射系數，結果如圖2所示。

圖2 植被效應與去除植被效應后向散射系數比較Fig.2 Comparison of backscatter coefficients before and after removal of surface vegetation cover in April

2.3 SVM算法構建

已有研究證明，與其他核函數相比，徑向基核函數(RBF)在土壤含水率預測模型中效果更優[36]。本文選擇RBF核函數建立土壤含水率預測模型，其功能選擇為回歸功能。目前研究發現微波遙感數據對土壤水分變化的敏感度也非常強，而且微波數據具有實時觀測能力，但僅考慮后向散射系數和土壤介電特性等因素的微波遙感反演土壤水分是遠遠不夠的[37]，因其未考慮地表上的植被覆被等因素。因此，本文選擇國內外水鹽遙感研究最常用的多光譜數據(Landsat OLI)和合成孔徑雷達數據(Sentinel-1A)，選擇支持向量機(Support vector machine，SVM) 算法，此方法在土壤水分監測方面反演精度較高[38]。將干旱指數、土壤含水率、后向散射系數進行耦合，區域地下水埋深與影響因子之間的相互關系可以看作復雜的非線性函數問題。支持向量機非線性基本思想是：選擇適合的核函數，通過空間變換(將輸入向量z映射到一個高維的特征空間上)，實現其線性可分性，使得到的結果風險小。根據泛函數的相關理論，只要核函數K(xi，yi)滿足Mercer條件，它就對應于某一變換空間中的內積，公式為

(11)

而對應的最優分類函數形式可寫為

(12)

式中n——訓練樣本個數

xi、yi、yj——訓練樣本向量

b*——偏置項系數

2.4 支持向量機回歸算法操作流程

基于Matlab平臺，以土壤含水率預測為因變量，使用支持向量機回歸算法進行模型的構建，其基本思想步驟如下：①建模集樣本和驗證集樣本的劃分，本文選擇26個樣點用于模型的訓練，12個樣本用于模型的驗證。②SVM回歸模型的創建，利用svmtrain函數進行樣本訓練，在訓練時考慮數據的實際情況對數據進行歸一化處理。選擇高斯徑向基核函數(RBF)，并通過交叉驗證網格搜索法(Grid-search)，利用訓練樣本反復地調整懲罰因子c和參數g，獲得最佳的參數組合。③利用svmpredict函數實現回歸模型的預測。④精度評價通過svmpredict函數進行預測時產生的均方根誤差RMSE和決定系數R2等評價標準，即可對建立的SVM回歸模型的預測能力進行綜合評價。

3 結果與分析

3.1 地下水埋深和土壤含水率的相關性分析

野外實際考察發現，農業生產中，當地于1月初進行冬灌和4月底進行春灌，以改善土地的水分條件。2015年1— 4月期間，沒有灌溉用水、降雨量極少(圖3)，且綠洲蒸發量極大，此時土壤中的水分來自地下水埋深的補給。

圖3 月平均降雨量Fig.3 Monthly average precipitation

通過分析不同深度的土壤含水率與地下水埋深之間的相互關系來探求兩者之間的線性關系，通過相關性分析得出(表1)，0～10 cm的土壤含水率與地下水埋深的決定系數達到了0.792 2，擬合效果相對較好，其次是40～60 cm的土壤含水率，其決定系數為0.722 3，10～20 cm和20～40 cm的土壤含水率與地下水埋深的決定系數分別為0.311 8和0.103 1，擬合效果相對較差。地下水埋深和不同深度的土壤含水率的擬合效果由高到低依次為0～10 cm、40～60 cm、10～20 cm、20～40 cm。

表1 地下水埋深與不同深度的土壤含水率的回歸分析Tab.1 Regression analysis of groundwater depth and soil moisture content at different depths

注：** 表示在0.01水平上極顯著，*表示在0.05水平上顯著。

結果表明，0～10 cm的土壤含水率與地下水埋深的擬合效果最優，土壤表層的土壤含水率與地下水埋深之間的相互關系最為密切，利用土壤表層的含水率來反演地下水埋深是可行的。

3.2 土壤水分反演模型精度比較

從圖4a～4d可看出，當只考慮TVDI作為SVM模型參數時，TVDIMSAVI、TVDINDVI、TVDIEVI、TVDIDVI均可以反映土壤表層的土壤濕度狀況，可用于土壤表層的水分反演。建模集R2分別為0.74、0.66、0.71和0.55，RMSE分別為4.66%、4.96%、4.71%和4.83%，TVDIMSAVI在SVM建模的結果R2為0.74，RMSE為4.66%，驗證集中的R2分別是0.70、0.69、0.74和0.59；RMSE分別為4.65%、4.91%、4.83%和4.56%。驗證結果與建模結果基本相同，在4種不同植被指數構建特征空間得到的干旱指數中，TVDIMSAVI精度最優，其次為TVDIEVI，最差的為TVDIDVI。建模集和驗證集不同植被指數的TVDI預測出來的土壤含水率相關性由大到小排序為：TVDIMSAVI、TVDIEVI、TVDINDVI、TVDIDVI。

表2 不同參數的SVM模型預測方案Tab.2 SVM model prediction scheme with different parameters

圖4 基于不同參數的SVM建模精度和驗證精度Fig.4 SVM modeling accuracies and verification accuracies with different parameters

表3 SVM回歸模型的最佳參數組合Tab.3 SVM regression model of the best combination of parameters

3.3 土壤水分空間分布特征

圖5 基于TVDI和參數反演的土壤含水率空間分布Fig.5 Spatial distribution of soil moisture content based on TVDI and parameter inversion

3.4 地下水埋深反演

圖6 土壤含水率反演的地下水埋深Fig.6 Retrieved result of groundwater depth by soil moisture content

3.5 反演結果驗證

在得到研究區地下水埋深分布圖后，選取HOBO水位計的數據與反演結果進行驗證，對水位計實測數據和反演結果進行比較，如圖7所示。SVM模型反演的地下水埋深與實測值的平均誤差為8.23%，優于李相等[40]得出的18.06%的結果。

圖7 預測值與實測值結果比較Fig.7 Comparison of predicted and measured results

4 結論

(1)通過對不同深度土層的土壤含水率與地下水埋深進行相關性分析，得出與地下水埋深相關性最優的為0～10 cm的土壤含水率，40～60 cm土層的土壤含水率次之，10～20 cm和20～40 cm的土層相關性較低。