基于MTVDI與DDI二元回歸模型對毛烏素沙地腹部土壤表層水分的研究

王思楠，李瑞平，韓 剛，田 鑫，王耀強，胡勇平

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

土壤水分是監測土地退化的重要指標，同時也是衡量土壤干旱程度的重要指標[1]。及時知曉土壤水分狀況，可以了解旱情程度[2]。遙感技術能夠有效、大面積、實時動態地獲取干旱地區旱情資料[3-5]，為各級政府與農業生產部門提供決策依據。因此，利用遙感進行旱情監測是一個研究和應用的熱點[6-8]。作為同時與歸一化植被指數(Normalized vegetation index，NDVI)和陸地表面溫度(Land surface temperature，LST)相關的溫度植被干旱指數(Temperature vegetation drought index，TVDI)可用于干旱監測[9-12]，尤其是監測特定年內某一時期整個區域的相對干旱程度。

目前利用地表溫度-植被指數(Ts-NDVI)特征空間進行旱情監測已取得一定的進展[13]，但是使用歸一化植被指數容易造成紅光飽和，背景的土壤噪聲也在一定程度上損害了NDVI的空間一致性。鑒于此，利用改進的TVDI模型降低實際計算過程中的偏差，同時提高土壤水分反演計算的準確性就顯得尤為重要。伍漫春等[14]在此基礎上采用土壤調節植被指數(Modified soil adjusted vegetation index，MSAVI)對其進行改進，證明地表溫度-土壤調節植被指數(Ts-MSAVI)能夠更好地反映區域土壤水分狀況，是一種更有效的土壤水分監測方法[15]。

雖然TVDI 以及改進的MTVDI都可以較好地監測土壤水分，但是精度還不是很高。針對這一問題，本研究根據研究區的特點采用修正土壤調整植被指數MSAVI構成的MTVDI指數，并利用野外同步實測的土壤水分數據分析和比較不同土壤深度的監測性能，探討荒漠化差值指數(Desertification difference index，DDI)對毛烏素沙地腹部土壤水分狀況的監測精度的影響，能夠對毛烏素沙地農牧業活動提供一定的指導，同時也為改善大尺度區域的氣候及改善毛烏素沙地生態環境打下基礎。

1 研究方法

1.1 研究區概況

研究區為毛烏素沙地腹地的烏審旗,它位于鄂爾多斯高原西南部，處于蒙、陜、寧經濟發展的“金三角”地帶。東經108°17′-109°40′,北緯37°38 ′-39°23′，面積11 645 km2。其所處位置為年平均氣溫6.8℃，多年平均降水量350～400 mm，多年平均蒸發量2 443 mm，海拔一般在1 300～1 400 m的干旱區，生態環境脆弱，沙生植被是本地區植被主體。

圖1 采樣點分布Fig.1 Sampling site distribution

圖2 子樣點采樣設計Fig.2 Sampling design at each site

1.2 數據來源

1.2.1 遙感數據 研究用的2016年4月21日和9月28日兩期軌道號為128/33、128/34,覆蓋烏審旗的Landsat8 OLI遙感數據由地理空間數據云提供。對Landsat 影像數據進行輻射定標,利用Modtran 4模型進行大氣校正以及裁剪鑲嵌等一系列的處理。

1.2.2 土壤含水量數據 參照土地利用類型圖、荒漠化程度分布圖等資料，在2016年4月20日、21日與9月27日、28日借助手持GPS布設了23個樣區(圖1)，每個樣區中包含5個子樣點A1～A5共5組(1 km×1 km像元，如圖2)共計115個子樣點，重復0～10 cm、10～20 cm、20～30 cm的土層深度，逐點用土鉆采集樣品，并刮去樣本點土壤表層的浮土后迅速封裝。在實驗室利用烘干法計算土壤含水量數據。

1.3 數據處理

1.3.1 構建Albedo-NDVI特征空間 不同荒漠化信息相對應的地表反照率Albedo與歸一化植被指數NDVI的特征空間中具有顯著的線性關系，計算公式如下：

(1)

A=0.356α2+0.130α4+0.373α5+0.085α6

+0.072α7-0.0018

(2)

式中,A為地表反照率，α2、α4、α5、α6、α7為經過大氣校正后的第2、4、5、6、7波段的反射率。

為了進一步確定公式(4)中的K值，在研究區選擇分布于不同沙漠化類型的300個點，并進行歸一化處理，然后利用Albedo和NDVI兩組數據進行回歸擬合(如圖3)，得到相應的方程：

A=0.8438-0.5531×NDVI

(3)

1.3.2 荒漠化差值指數提取 Verstrate等[16]研究發現在代表荒漠化變化趨勢的垂直方向上劃分Albedo-NDVI特征空間，荒漠化差值指數DDI能夠有效地區分出不同程度的荒漠化土地，計算公式如下：

DDI=K×NDVI-Albedo

a×K=-1

(4)

式中,a為式(3)中垂線的斜率，由此，確定荒漠化差值指數DDI的最終表達式為：

DDI=1.808×NDVI-Albedo

(5)

1.3.3 植被指數-地表溫度特征空間 Sandholt等[17]在研究NDVI和Ts的散點圖時發現二者呈現三角形分布，在此基礎上，前人還發現NDVI和Ts構成的特征空間可以間接表示土壤含水狀況，因此提出了溫度植被干旱指數的概念。本文使用TVDI由改進型土壤調整植被指數和地表溫度[18]計算得到，其定義為公式：

MTVDI=(Ts-Tsmin)/(Tsmax-Tsmin)

(6)

圖3 Albedo與NDVI的回歸分析Fig.3 Albedo and NDVI regression analysis

Tsmax=MSAVIa+b
Tsmin=MSAVIc+d

(7)

式中,Ts為任意像元地表溫度；Tsmax為干邊上的地表溫度；Tsmin為濕邊上的地表溫度。a、b、c、d是干、濕邊通過線性擬合的模型參數。

2 結果分析

2.1 基于Landsat8數據的Ts和MSAVI構建特征空間

通過Landsat8數據計算的植被指數對應地表溫度的最大值和最小值，構建毛烏素沙地腹部2期不同時相的Ts-MSAVI特征空間。根據公式(6)并擬合干、濕邊方程(表1)。結果表明，MSAVI對應Ts的最大值和最小值呈近似線性關系，隨著植被指數的增大，地表溫度的最大值呈減小趨勢，地表溫度的最小值呈增大趨勢，二者差值呈減小趨勢。從圖4可以看出2期不同時相的散點圖的形狀相似[19]，其中9月的干、濕邊的R2都比4月的大。

2.2 Landsat8-MTVDI與土壤水分相關性分析

利用野外實測采樣點0～10 cm、10～20 cm、20～30 cm土層的土壤含水量值與MTVDI值進行回歸分析用于反演精度驗證。從圖5和圖6中可以看出不同遙感數據獲取的不同月份MTVDI和土壤各層含水量具有一定的負相關。即MTVDI越高，土壤含水量越低，滿足MTVDI指數值越大土壤水分越低的原理。4月份MTVDI指數與0～10 cm、10～20 cm、20～30 cm土壤含水量的R2值分別為0.656、0.646、0.637，整體高于9月份R2值0.457、0.436、0.431。不管是低植被覆蓋度的4月還是高植被覆蓋度的9月，在0～10 cm深度的R2值均大于10～20 cm、20～30 cm深度的R2值。說明MTVDI能夠較好地反映土壤表層的含水量狀況。

圖4 2016年4月(a)、9月(b)Landast8的Ts-MSAVI特征空間Fig.4 Ts-MSAVI feature space of Landast8 in April(a)and September(b)2016

表1 不同時相干濕邊擬合結果

圖5 2016年4月不同深度下土壤水分與MTVDI回歸分析Fig.5 The regression analysis of soil moisture and MTVDI at different depths in April 2016

圖6 2016年9月不同深度下土壤水分與MTVDI回歸分析Fig.6 The regression analysis of soil moisture and MTVDI at different depths in September 2016

2.3 基于MTVDI與DDI的二元線性回歸干旱監測模型的分析

毛烏素沙地荒漠化程度的加劇，在土壤水分補給和散失的過程當中又受到土地利用、地貌、地形等多種因素的影響。本研究考慮荒漠化因素的條件下，利用Landsat8遙感數據提取2016年4月、9月不同時相的MTVDI指數和DDI指數并與0～10 cm土層深度含水率建立二元線性監測模型(SMC=a×MTVDI+b×DDI+c)，利用野外實測土壤0～10 cm土層含水量進行驗證。最后和單一使用MTVDI指數與0～10 cm土層深度含水量建立的一元線性監測模型(SMC=a×MTVDI+c)進行精度對比分析，結果如表2、3所示。

分析表2、3可知，2016年2個不同時相在0～10 cm土層深度的含水量的反演過程中，單獨利用MTVDI監測含水量的相關系數分別為0.810、0.676，相對誤差分別為12.93%、14.35%，平均相對誤差為13.64%。當引進了DDI指數影響因子之后，二元線性監測模型中相關系數分別為0.842、0.734，相對誤差分別為10.26%、11.64%，平均相對誤差為10.95%。二元線性監測模型提高了反演土壤含水量的精度。說明在2個不同時相土壤表層0～10 cm深度的含水量受荒漠化信息的影響較大。二元線性回歸模型監測0～10 cm土層深度的相對含水量空間分布，如圖7所示。

由圖7可知：上述2期土壤含水量反演圖均可反映烏審旗當時的土壤含水量分布情況。分析表4可知，2016年2期數據中可得到該區域0～10 cm表層土壤含水量在4月份5%～10%之間的區域，占總面積的53.72 %以上，達到了6 256 km2，含水量偏低。0～10 cm表層土壤含水量在9月份10%～15%之間的區域，均占總面積的51.55 %以上，達到了6 003 km2，含水量偏高。與此同時根據研究區氣象局2016年氣象資料可知該年9月份降水偏高，與本文研究的結果基本一致。

表2 MTVDI模型對0～10 cm土層深度含水量的回歸分析

注：a,c為一元線性監測模型的系數。

Note: “a” and “c” are the coefficients of a linear monitored model.

表3 MTVDI與DDI二元監測模型對0～10 cm深度含水量的回歸分析

注：a,b,c為二元線性監測模型的系數。

Note: “a”,“b” and “c” are the coefficients of binary linear monitored model.

圖7 2016年MTVDI與DDI二元監測模型反演0～10 cm深度含水量空間分布Fig.7 Inversion of spatial distribution of water content in 0～10 cm soil depth in 2016 modeled by MTVDI and DDI

表4 土壤表層含水量面積分布

3 結論與討論

本研究建立了毛烏素沙地腹部MTVDI與DDI二元線性監測模型與土壤表層0～10 cm含水率的關系模型，對該地區土壤表層0～10 cm含水率進行反演與精度分析并統計表層土壤含水率的面積分布。

1)不同植被覆蓋度下的MTVDI均能反映不同深度的土壤含水量，且呈現不同程度的負相關，R2值4月份整體高于9月份，在0～10 cm深度的R2值均大于10～20 cm、20～30 cm深度的R2值，并且都高于0.4。說明MTVDI可以作為有效指示地表土壤水分狀況的指標。

2)荒漠化程度加劇，地表覆蓋程度下降，地表能量與水分平衡發生變化，均可導致土壤水分發生改變。荒漠化指數DDI與MTVDI結合建立二元線性回歸模型監測區域0～10 cm土層深度含水率，平均相對誤差值要比單獨TVDI小2.69 %。

3)通過二元線性回歸監測模型反演的土壤表層含水量分布圖發現：該區域0～10 cm表層土壤含水量在4月份5%～10%之間的區域，占總面積的53.72%以上，達到了6256km2，含水量偏低；0～10 cm表層土壤含水量在9月份10%～15%之間的區域，達到了6 003 km2，均占總面積的51.55%以上，含水量偏高。

本研究在MTVDI指數的基礎上，考慮了毛烏素沙地腹部荒漠化信息，建立MTVDI指數和DDI指數的二元線性回歸監測模型，在一定程度上提高了MTVDI指數指示土壤水分的合理性。從另一個角度看，MTVDI指數指示土壤水分的能力與DDI指數有很大關系，使用單一的MTVDI指數無法精確達到反映干旱特征的目的，這也促使干旱監測由單因素向多因素綜合發展，綜合多指數的干旱監測模型是研究復雜的干旱監測問題的新途徑，在解決干旱監測的復雜問題中有著較大的應用潛力。本研究尚屬可行性研究，但單期的影像不能說明研究結果的普遍性，在后續的研究中將會進一步考慮采用多期遙感數據做動態分析。