張鈞泳 丁建麗 譚 嬌

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

0 引言

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

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

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

1 研究區(qū)概況和數(shù)據(jù)資料

1.1 研究區(qū)概況

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

圖1 研究區(qū)位置圖Fig.1 Location of study area

1.2 實測數(shù)據(jù)

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

1.3 遙感數(shù)據(jù)

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

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

1.4 遙感影像預(yù)處理

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

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

2 研究方法

2.1 不同溫度-植被指數(shù)構(gòu)建

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

(1)

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

Tsmax=a1+b1VI

(2)

Tsmin=a2+b2VI

(3)

式中TVDI——溫度植被指數(shù)

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

VI——植被指數(shù)

Tsmin——相同植被指數(shù)值對應(yīng)的最小地表溫度

Tsmax——相同植被指數(shù)值對應(yīng)的最大地表溫度

a1、b1——干邊線性擬合方程的系數(shù)

a2、b2——濕邊線性擬合方程的系數(shù)

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

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

2.2 土壤后向散射系數(shù)計算

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

(4)

(5)

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

(6)

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

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

θ——雷達(dá)波入射角

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

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

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

VWC=2.15NDMI+0.32

(7)

其中

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

(8)

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

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

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

(9)

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

(10)

(11)

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

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

2.3 SVM算法構(gòu)建

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

(11)

而對應(yīng)的最優(yōu)分類函數(shù)形式可寫為

(12)

式中n——訓(xùn)練樣本個數(shù)

xi、yi、yj——訓(xùn)練樣本向量

b*——偏置項系數(shù)

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

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

3 結(jié)果與分析

3.1 地下水埋深和土壤含水率的相關(guān)性分析

野外實際考察發(fā)現(xiàn)，農(nóng)業(yè)生產(chǎn)中，當(dāng)?shù)赜?月初進(jìn)行冬灌和4月底進(jìn)行春灌，以改善土地的水分條件。2015年1— 4月期間，沒有灌溉用水、降雨量極少(圖3)，且綠洲蒸發(fā)量極大，此時土壤中的水分來自地下水埋深的補給。

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

通過分析不同深度的土壤含水率與地下水埋深之間的相互關(guān)系來探求兩者之間的線性關(guān)系，通過相關(guān)性分析得出(表1)，0～10 cm的土壤含水率與地下水埋深的決定系數(shù)達(dá)到了0.792 2，擬合效果相對較好，其次是40～60 cm的土壤含水率，其決定系數(shù)為0.722 3，10～20 cm和20～40 cm的土壤含水率與地下水埋深的決定系數(shù)分別為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水平上顯著。

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

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

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

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

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

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

3.3 土壤水分空間分布特征

圖5 基于TVDI和參數(shù)反演的土壤含水率空間分布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 反演結(jié)果驗證

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

圖7 預(yù)測值與實測值結(jié)果比較Fig.7 Comparison of predicted and measured results

4 結(jié)論

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