基于綜合干旱指數的淮河流域土壤含水量反演

張 文， 任 燕， 馬曉琳， 胡藝杰

(1.武漢大學遙感信息工程學院，武漢 430079； 2.國家知識產權局專利局專利審查協作河南中心，鄭州450046； 3.河南省白沙水庫管理局，禹州 461670； 4.華北水利水電大學資源與環境學院，鄭州 450045)

0 引言

陸表土壤水不僅是旱情監測的重要指標，還是氣候、水文、生態、農業等領域的重要參數，也是全球氣候變化的重要組成部分[1]。傳統測量方法只能采集到點和小范圍的墑情信息[2]，而遙感技術的發展為獲取大范圍地表土壤水分提供了有效手段。遙感的主要特點是探測范圍廣、速度快、周期短、受地面條件限制少，能夠頻繁持久地提供地表的信息[3]。因此，利用遙感技術來進行大區域的土壤水分信息提取具有極其廣闊的應用前景。

從20世紀60年代開始，遙感技術已經開始應用到土壤水分監測中，歷年來，國內外科研工作者提出過許多模型和方法。1969年，Jordan[4]提出了最早的一種植被指數，比值植被指數(ratio vegetation index，RVI)； 1977年，Kahle[5]得到第一個遙感熱慣量影像; Price給出了用遙感數據計算熱慣量的一般方法[6]，并于1985年提出了表觀熱慣量(apparent thermal inertia，ATI)模型[7]； Jackson等[8]綜合分析研究土壤水分、葉片溫度和植被指數相互之間關系,于1981年提出了作物缺水指數法(crop water stress index，CWSI)； Carlson等[9]于1994年提出了綜合考慮植被指數和植被冠層溫度的植被供水指數(vegetation supply water index，VSWI)； 2002年，Sandholt等[10]運用陸地表面溫度(land surface temperature，LST)與歸一化植被指數(normalized differential vegetation index，NDVI)，建立了NDVI-Ts特征空間，并提出了溫度植被干旱指數(temperature vegetation dryness index，TVDI)。各種干旱指數都有其優缺點與不同的適用范圍，綜合多種干旱指數對土壤水分進行反演會得到更可靠的效果。因ATI模型與VSWI模型計算簡便，效果良好，本文選取這2種模型進行實驗。已有研究表明ATI模型只適用于低植被覆蓋區域[11]，而VSWI模型在作物覆蓋度較高時比較有效[9]，本文通過NDVI來區分地表植被覆蓋度，結合ATI與VSWI模型，建立綜合干旱指數(comprehensive drought index，CDI)模型，進行土壤水分的初步反演。

但是僅利用遙感數據來監測土壤水分也有不足，因為利用上述模型基于遙感影像提取到的土壤水分信息都只能反映土壤的相對干濕情況，并不是真實的土壤含水量值，對土壤含水量的定量分析和使用有較大的限制。因此，本研究引入實測的土壤含水量數據，探究CDI結果與實測土壤含水量之間的相關關系，選出一種最佳的相關模型，則可利用該模型把CDI結果轉化為真實的土壤含水量數據，以期提高大區域的陸表土壤含水量監測效率，便于土壤水分產品的業務化生產。

1 研究區概況及數據源

1.1 研究區概況

淮河流域地處我國東部，介于長江和黃河2流域之間，位于E111°55′～121°25′，N30°55′～36°36′，面積約27萬km2(圖1)。淮河流域西部、西南部及東北部為山區、丘陵區，其余為廣闊的平原[12]。平原面積約占總面積的2/3，在我國農業生產中占有重要地位。

圖1 研究區及實測站點示意圖

淮河流域地處我國南北氣候過渡帶，冬春干旱少雨，夏秋悶熱多雨，冷暖和旱澇轉變急劇，多年平均降水量約為920 mm，由南向北遞減，山區多于平原，沿海大于內陸，旱澇時有發生。因此，對淮河流域進行業務化的土壤含水量監測具有極大的價值。因其流域面積大，地貌特征復雜，利用MODIS影像選取ATI和VSWI模型進行綜合研究，是一種有效可行的方法。

1.2 數據源

1.2.1 MODIS數據

MODIS數據是美國國家航空航天局(National Aeronautics and Space Administration,NASA)對地觀測系統系列遙感衛星平臺上的主要傳感器，它具有36個光譜通道，每1～2 d可以獲取一次全球地表數據[13]。因其波段范圍廣，數據更新頻率快，在全球范圍免費接收等優點，本實驗選用MODIS數據為主要數據源[14]。

選取空間分辨率為1 000 m的MODIS 1B產品數據，即MOD021KM, 以及其對應的地理定位文件，即MOD03。使用ENVI+IDL進行編程，完成數據預處理和相關指數的計算。

1.2.2 實測數據

地面實測土壤含水量數據來自中國水利部數據中心，其墑情數據主要由水利部門的土壤墑情監測站提供。墑情監測站點以縣為單元，根據氣候類型、地形地貌、作物布局、灌排條件和土壤類型等因素，設置在區域范圍內代表性較強的地塊。墑情采集方法主要為固定監測法，即埋設固定式自動監測設備，傳感器分別埋入土層深度10 cm，20 cm和40 cm處，按0～10 cm，10～20 cm和20～40 cm共3個層次監測土壤含水量，并于每月的1日、11日和21日各觀測一次。

淮河流域耕地面積為1 333 hm2，農作物以冬小麥、水稻、棉花和油菜等為主。每年的5月上中旬是冬小麥的抽穗揚花期，5月中下旬為灌漿乳熟期，5月底到6月上旬為收割期，因此5月份是冬小麥生長需水的最關鍵時期，一旦出現旱情，則會抑制子粒灌漿及干物質向子粒的運輸與積累，從而直接影響到小麥單產水平。故本文選取了2014年5月1日、11日、21日和6月1日4天作為實驗日。其中5月1日、11日屬于小麥生長的中后期，葉面基本覆蓋農田； 5月21日和6月1日屬于小麥成熟期，部分地區小麥開始收割，地表覆蓋情況不一。

2 研究方法

2.1 干旱指數模型選擇

目前利用衛星影像反演土壤水分的方法己經有很多，例如ATI、距平植被指數、CWSI、條件溫度植被指數(vegetation-temperature condition index,VTCI)、VSWI等。現將使用比較廣泛的幾種方法進行介紹與分析。

2.1.1 ATI模型

土壤熱慣量隨土壤含水量的變化而變化，因此可建立土壤熱慣量模型來監測土壤含水量。Price[7]提出的ATI是在熱慣量定義的基礎上，不考慮太陽高度角和緯度等因素，其形式為

ATI=(1-A)/(Tmax-Tmin)，

(1)

式中：ATI為表觀熱慣量；A為全波段反照率，可由 MODIS數據 1和2 通道的反射率得到；Tmax和Tmin分別為一天中最高和最低溫度，可分別由MODIS數據 31 通道的地表溫度得到。ATI值越高，表示土壤含水量越大，反之亦然。ATI模型參數都可從遙感影像中獲取，且計算簡便，但僅適用于低植被覆蓋區域[10]。

2.1.2 CWSI模型

CWSI以能量平衡為基礎，是最常用的植被蒸散法之一，定義為[15]

CWSI=1-ET/ETp，

(2)

式中：CWSI為作物缺水指數；ET為水分的日蒸散量；ETp為在水分供應充分條件下的日潛在蒸散量。ET值越小，CWSI值越大，土壤的含水量也越少，反之亦然。但是ET與ETp值不能夠從遙感影像上直接獲取，需要大量的地面實測資料，應用起來也比較困難。

2.1.3 VTCI模型

在NDVI-Ts構成三角形空間的基礎上，王鵬新等[16]提出了VTCI模型 ，其計算公式為

(3)

式中：VTCI為條件植被溫度指數；LSTNDVIi,max和LSTNDVIi,min分別為研究區域內具有相同NDVI值的像元的最高溫度和最低溫度;LSTNDVImax為NDVI最大值相應像元的溫度；LSTNDVIi為NDVI值為NDVIi的相應像元的溫度。VTCI的值越小代表土壤含水量越低，反之亦然。該模型對大區域的旱情監測效果較好，但計算比較復雜。

2.1.4 VSWI模型

VSWI是通過計算NDVI和植物冠層溫度的比值得到的，其公式為[17]

VSWI=NDVI/T，

(4)

式中：VSWI為植被供水指數；T為植被的冠層溫度。VSWI值越大，表明土壤含水量越大； VSWI值越小則代表土壤含水量越小。但VSWI值是根據植被覆蓋狀況的變化來進行反演的，所以不適用于低植被覆蓋地區[18]。

2.2 建立綜合干旱指數模型

淮河流域面積大，地貌特征復雜，僅用一種干旱指數無法滿足監測復雜地表土壤含水量的需求。通過2.1節對各干旱指數模型的優缺點與適用性特征的分析，可知ATI與VSWI模型效果良好且計算簡便，但ATI模型比較適用于低植被覆蓋地區和裸土區域，而VSWI模型適用于植被覆蓋度較高的地區，所以本文結合這2種模型來建立一種新的干旱指數模型——CDI模型。

建立CDI模型需要對植被覆蓋情況進行分類，而NDVI反映了植被覆蓋情況，故本文利用NDVI設置閾值來區分植被覆蓋類型，以便合理使用適合不同覆蓋類型的指數來反演土壤水分，以提高整個區域土壤水分反演精度。已有實驗證明，在NDVI>0.3的情況下VSWI非常有效[19]，而在NDVI<0.35的情況下ATI非常有效[20]，因此本研究閾值選為0.33。

由于ATI和VSWI的取值范圍不在同一個尺度上，因此在綜合使用這2種模型之前，需要把它們歸一化到同一個區間范圍內，以便最終建立CDI，其公式為

(5)

式中：CDIi為在任意像元點i的綜合干旱指數；ATIi為任意像元點i的表觀熱慣量；VSWIi為任意像元點i的植被供水指數；ATImax和ATImin分別為表觀熱慣量的最大值和最小值；VSWImax和VSWImin分別為植被供水指數的最大值和最小值。

2.3 基于實測數據建立相關模型

根據2.2節所述的原理與公式，利用ENVI+IDL編程可以得到每個實驗日的CDI結果。但CDI的結果是個無量綱常數，只能反映土壤相對的干濕情況，并不能代表土壤含水量的定量結果。為了得到CDI與實際土壤含水量之間的關系，需要在這2種數據上分別取樣點進行分析，建立兩者之間的函數關系，分別用經典的線性、指數和對數函數進行擬合，根據精度選擇一個最佳的函數對應關系。

實測站點分布如圖1所示，選擇2014年5月1日的數據進行建模。因為CDI值是遙感影像反演得到的，而理論上遙感只能穿透土壤表層，所以選用0～10 cm深度的實測數據來與CDI值進行建模。根據2.2節所述，首先用ENVI編程計算得到5月1日的CDI結果； 然后用ArcGIS把實測點對應的CDI值提取出來，去掉其中CDI值為空的數據，再去掉實測點為0或極端大的異常數據； 最后留下了81組樣點。這些樣點均勻分布在研究區域內，基本上能代表整個研究區域的情況。

常見的用于統計建模的數學模型有3種： 線性模型、對數模型和指數模型。為了建立選擇的81組數據之間的函數關系，分別用這3種模型進行擬合，然后通過比較其精度來選擇最優模型。相同的置信度下，用于判定擬合優度的指標相關系數R2越大越好。通過實驗，3種模型結果如圖2所示。由圖可知選取線性模型函數進行擬合時，R2值最高。于是擬建立線性回歸關系，進一步做置信度分析，結果如表1所示。

(a) 線性關系 (b) 對數關系 (c) 指數關系

圖2 擬合關系對比Fig.2 Comparison of fitting relationships表1 淮河流域土壤含水量與CDI值之間的線性統計關系Tab.1 Linear statistical relationship between measured soilmoisture and CDI value in Huaihe River Basin

綜上，得到了把CDI轉換為真實土壤含水量值的模型,即

SM=0.265 6CDI+0.048 6，

(6)

式中SM為土壤含水量。

對該模型的可靠性進行檢驗，顯著水平取α=0.01，R=0.87>0.283=Rn-2，α，檢驗通過。F=434.67>7.08=F1-α(1，n-2)，檢驗通過。SignificanceF是在顯著性水平下的Fα臨界值，等于P值，故本例中，P=7.43e-34<0.001，故置信度達到99.9%以上。因此可以使用該線性模型把淮河流域的CDI結果轉化為實際土壤含水量。

3 結果與分析

3.1 土壤含水量反演結果

如上述方法，利用MODIS影像計算出CDI指數后，根據線性模型式(6)，反演出2014年5月11日、21日和6月1日3 d的淮河流域土壤含水量結果，如圖3所示。

(a) 5月11日 (b) 5月21日 (c)6月1日

圖3淮河流域土壤含水量分布

Fig.3DistributionofsoilmoistureinHuaiheRiverBasin

由圖3可見，土壤含水量分布大致趨勢為南部地區略大于北部地區，東部地區略大于西部地區。江蘇省土壤含水量普遍較其他區域略高，可能與靠近海洋空氣濕潤有關。而河南和安徽等地區土壤含水量則相對略低，這主要是因為華北地區蒸發較強，夏季風弱，地下水位低。初步判斷土壤含水量反演結果比較合理。整體看淮河流域土壤含水量大多在0.1～0.3 m3/m3之間，5月11日整體相對較高，5月21日整體相對較低。且5月11日土壤含水量分布不均勻，浮動較大，土壤含水量高與低的區域過度較不平緩。

3.2 精度驗證及分析

為了驗證本文提出的反演模型的可靠性，結合地面實測土壤含水量數據，對反演得到的土壤含水量結果進行精度驗證。實測數據點去除數據為 0 和異常大值的記錄，最終每個實驗日選擇了90個左右質量比較好的數據點。

圖4為實測0～10 cm、10～20 cm深度土壤含水量與反演土壤含水量的對比。可看出反演的土壤含水量與實測土壤含水量無論是從絕對值還是變化趨勢上都大體一致，特別是0～10 cm深度土壤含水量與反演結果相似度很高，只是反演值比實測值略偏低。0～20 cm深度的土壤含水量也與反演結果比較相關，但是相較而言相關性不如0～10 cm深度數據。這主要是由于遙感監測穿透土壤深度有限，所以反演結果與淺層土壤含水量的相關性更好。可以看出3個實驗日中2014年5月21日土壤含水量較低，大致在0.1～0.2 m3/m3之間，在圖3(b)上也有表現，其他實驗日土壤含水量均在0.2 m3/m3左右浮動。

(a) 5月11日 (b) 5月21日 (c) 6月1日

圖4淮河流域實測土壤含水量與反演土壤含水量對比

Fig.4ComparisonofmeasuredsoilmoistureandestimatedsoilmoistureinHuaiheRiverBasin

具體精度需要進行誤差統計分析，定義地面觀測土壤含水量為X，反演的土壤含水量為X0，N為參與比較的總樣本數。通過對比分析，統計出最大誤差(MaxE)、絕對誤差(ABVR)和均方根誤差(RMSE)。

3個實驗日的實測0～10 cm、10～20 cm土壤含水量數據與反演得到的土壤含水量數據的誤差分析如表2所示。

表2 實測土壤含水量與反演土壤含水量誤差分析Tab.2 Error analysis of measured soil moisturewith estimated soil moisture

根據表2中的數據可知，反演精度比較理想，特別是與0～10 cm深度土壤含水量相關性較好，平均RMSE為0.022 3。而反演結果與10～20 cm深度土壤含水量的相關性略差，因為遙感只能穿透表層土壤，隨著土壤深度增加，遙感反演的精度會降低。反演結果與0～10 cm和10～20 cm深度的RMSE最小值都出現在5月21日，從圖3(b)和圖4(b)中可以看出，當日土壤含水量在0.1～0.2m3/m3左右，相對較低。所以推斷在土壤水分較低的區域，該反演模型的結果相對更加精確。而5月11日RMSE較大，從圖3(a)和圖4(a)中可以看出該日土壤含水量相對較高且浮動范圍比較大，特別是從圖3(a)可以看出土壤含水量高與低的區域過度不太平滑，可能是由于局部灌溉造成，在此種狀況下模型反演精度略有降低。

(a) 5月11日 (b) 5月21日 (c) 6月1日

圖50～10cm實測土壤含水量與反演土壤含水量相關性分析

Fig.5Correlationanalysisofmeasured0～10cmsoilmoistureandestimatedsoilmoisture

圖5為0～10 cm深度的實測站點土壤含水量與反演土壤含水量之間的線性回歸分析結果。因為上文結果已經證實，遙感反演的土壤含水量與0～10 cm深度的實際土壤含水量關系比較密切，所以在此著重分析反演結果與0～10 cm深度土壤含水量的關系。圖中紅線是擬合的直線，黑線是斜率為1的輔助線，理論上，紅線應該與黑線重合。紅線在黑線下的部分說明反演的土壤含水量較實測值低估，紅線在黑線上的部分說明反演的土壤含水量較實測值高估。由圖5中的相關系數可以看出，3個實驗日中，相關系數R2均在0.7左右，說明反演的土壤含水量結果是比較可靠的，可以用反演的土壤含水量來反應實際的土壤含水量情況。

4 結論與展望

本文基于MODIS數據與實測數據，反演了淮河流域的土壤含水量。該反演方法具有以下優點：

1)提出的CDI模型考慮了低植被覆蓋區與高植被覆蓋區的不同適用性特點，綜合ATI和VSWI這2種模型來進行研究。CDI模型可適用于復雜的植被覆蓋區域，相較使用單一的反演模型，反演精度將大大提高。

2)建立的CDI與實測土壤含水量轉換模型，可對遙感反演的無量綱結果進行轉換，從而得到真實的土壤含水量值。這對于土壤含水量的定量監測與分析有較大的意義。

經過精度驗證，反演的土壤含水量與0～10 cm和10～20 cm深度土壤含水量的平均均方根誤差分別為0.022 3和0.036 0，特別是與0～10 cm深度數據的相關性較高，R2在0.7左右，精度較好，可滿足日常業務化監測需要，且簡便快捷，具有較大的應用潛力。

本次實驗中也存在一些誤差，在2.3節建立相關模型時，選取的是2014年5月1日的樣本數據，樣本所選取的日期不同，可能會對數據的相關性結果產生影響，這也會對反演的精度造成影響。將來可考慮進行更多后續實驗，以求獲得具有更高精度、更高效率的土壤含水量成熟產品。

參考文獻(References)：

[1] Chen C F,Son N T,Chang L Y,et al.Monitoring of soil moisture variability in relation to rice cropping systems in the Vietnamese Mekong Delta using MODIS data[J].Applied Geography, 2011,31(2):463-475.

[2] 劉 敏,王亮亮,蔡秋鵬.FDR和TDR測定幾種典型土壤含水量的對比分析[J].水利信息化,2016(6):32-36.

Liu M,Wang L L,Cai Q P.Compaveral typical soil’s water content measuring with FDR and TDR[J].Water Resources Informatization,2016(6):32-36.

[3] Wagner W,Hahn S,Kidd R,et al.The ASCAT soil moisture product:A review of its specifications,validation results,and emerging applications[J].Meteorologische Zeitschrift,2013,22(1):5-33.

[4] Jordan C F.Derivation of leaf-area index from quality of light on the forest floor[J].Ecology,1969,50(4):663-666.

[5] Kahle A B.A simple thermal model of the Earth’s surface for geologic mapping by remote sensing[J].Journal of Geophysical Research,1977,82(11):1673-1680.

[6] Price J C.Thermal inertia mapping:A new view of the earth[J].Journal of Geophysical Research,1977,82(18):2582-2590.

[7] Price J C.On the analysis of thermal infrared imagery:The limited utility of apparent thermal inertia[J].Remote Sensing of Environment,1985,18(1):59-73.

[8] Jackson R D,Pinter Jr P J.Detection of water stress in wheat by measurement of reflected solar and emitted thermal IR radiation[C]//Spectral Signatures of Objects in Remote Sensing.Versailles:Institut National de la Recherche Agronomique,1981:399-406.

[9] Carlson T N,Gillies R R,Perry E M.A method to make use of thermal infrared temperature and NDVI measurements to infer surface soil water content and fractional vegetation cover[J].Remote Sensing Reviews,1994,9(1/2):161-173.

[10] Sandholt I,Rasmussen K,Andersen J.A simple interpretation of the surface temperature/vegetation index space for assessment of surface moisture status[J].Remote Sensing of Environment,2002,79(2/3):213-224.

[11] 趙英時.遙感應用分析原理與方法[M].北京:科學出版社,2003.

Zhao Y S.Analysis Principle and Method of Remote Sensing Applications[M].Beijing:Science Press,2003.

[12] 郭月婷,徐建剛.基于模糊物元的淮河流域城市化與生態環境系統的耦合協調測度[J].應用生態學報,2013,24(5):1244-1252.

Guo Y T,Xu J G.Coupling coordination measurement of urbanization and eco-environment system in Huaihe River Basin of China based on fuzzy matter element theory[J].Chinese Journal of Applied Ecology,2013,24(5):1244-1252.

[13] 張 怡,何政偉,薛東劍.MODIS數據預處理方法[J].地理空間信息,2013,11(3):49-51.

Zhang Y,He Z W,Xue D J.Pre-processing method of MODIS data[J].Geospatial Information,2013,11(3):49-51.

[14] 刁 偉,孟令奎,張東映,等.MODIS數據在湖北旱情監測中的應用[J].水利信息化,2013(1):12-15.

Diao W,Meng L K,Zhang D Y,et al.Application of MODIS data on the drought monitoring of Hubei Province[J].Water Resources Informatization,2013(1):12-15.

[15] 申廣榮,田國良.基于GIS的黃淮海平原旱災遙感監測研究——作物缺水指數模型的實現[J].生態學報,2000,20(2):224-228.

Shen G R,Tian G L.Remote sensing monitoring of drought in Huanghe, Huaihe and Haihe Plain based on GIS：The calculation of crop water stress index model[J].Acta Ecologica Sinica,2000,20(2):224-228.

[16] 王鵬新,龔健雅,李小文.條件植被溫度指數及其在干旱監測中的應用[J].武漢大學學報,2001,26(5):412-418.

Wang P X,Gong J Y,Li X W.Vegetation temperature condition index and its application for drought monitoring[J].Geomatics and Information Science of Wuhan University,2001,26(5):412-418.

[17] Nichol J E,Abbas S.Integration of remote sensing datasets for local scale assessment and prediction of drought[J].Science of the Total Environment,2015,505:503-507.

[18] Cai G Y,Du M Y,Liu Y.Regional drought monitoring and analyzing using MODIS data：A case study in Yunnan Province[C]//Proceedings of the 4th IFIP TC 12 Conference on Computer and Computing Technologies in Agriculture IV.Nanchang,China:Springer,2010:243-251.

[19] 宋榮杰.基于遙感的旱區土壤濕度反演方法研究[D].楊凌:西北農林科技大學,2013.

Song R J.Study on Soil Moisture Inversion Method Based on Remote Sensing[D].Yangling:Northwest A and F University,2013.

[20] 吳 黎,張有智,解文歡,等.改進的表觀熱慣量法反演土壤含水量[J].國土資源遙感,2013,25(1):44-49.doi:10.6046/gtzyyg.2013.01.08.

Wu L,Zhang Y Z,Xie W H,et al.The inversion of soil water content by the improved apparent thermal inertia[J].Remote Sensing for Land and Resources,2013,25(1):44-49.doi:10.6046/gtzyyg.2013.01.08.