基于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指数无法精确达到反映干旱特征的目的，这也促使干旱监测由单因素向多因素综合发展，综合多指数的干旱监测模型是研究复杂的干旱监测问题的新途径，在解决干旱监测的复杂问题中有着较大的应用潜力。本研究尚属可行性研究，但单期的影像不能说明研究结果的普遍性，在后续的研究中将会进一步考虑采用多期遥感数据做动态分析。