基于Landsat 时序影像和非线性边界监测土壤干旱研究——以1986—2020 年黄河内蒙古段生长季为例

张新平，权 全

(1.西安理工大学 艺术与设计学院,陕西 西安 710054；2.西安理工大学 省部共建西北旱区生态水利国家重点实验室,陕西 西安 710048)

土壤水分（soil moisture,SM）是干旱胁迫生态系统生产的主要驱动因素[1]。目前，基于遥感反演SM 的方法主要有：基于反射率的光学方法（ORM）、被动微波（MWP）、主动微波（MWA）和协同方法（光学+热红外，MWP+MWA，微波+光学）。实际工作中，通常依据SM 与土壤温度植被干旱指数（TVDI）间的关系，通过遥感反演表层土壤的干旱程度来间接地表征土壤的相对含水量[2-3]。目前，监测土壤水分与旱情的遥感指数主要有：水分亏缺指数（WDI）、温度植被干旱指数（TVDI）、土壤水分亏缺指数（SWDI）、土壤水分指数（SMI）、垂直干旱指数（PDI）、改进的垂直干旱指数（MPDI）、第二种改进的垂直干旱指数（MPDI1）和温度植被土壤水分干燥指数（TVMDI）。其中，TVDI 仅需通过多光谱卫星影像定量计算得出的地表温度（LST）和植被指数（NDVI），因没有其他限制条件，成为应用最普遍的干旱指数。关于长时间、广空间高精度遥感反演TVDI 的研究，主要集中在多源数据融合[4-5]、反演模型改进[6]和下垫面分类考虑[7]。三者同时兼顾的研究较少，仅见于基于Landsat 和Sentinel 在单景影像空间尺度上和LST-NDVI 和STR-NDVI 特征空间中开展的土壤水分与干燥度的遥感反演研究[4]。已有研究认为，土壤含水量等值线与模拟干边呈非线性相交，特别是高NDVI 地区的蒸发修正量随NDVI 增大而增长更迅速，而低NDVI 则增长缓慢，所以遥感反演TVDI 的过程中需要进行非线性干湿边界的拟合修订[8]。前期研究中发现在地形复杂、地表覆盖物异质性高的区域，梯形或三角模型的干湿边界通常表现为曲线，且非线性边界的预测精度高于线性边界[9]。基于Landsat 卫星的多景影像拼接后，在不同的土地利用/覆被类型下，通过非线性干湿边界遥感反演土壤水分与干燥度的研究，尚未见报道。鉴此，本研究提出以下2 个假设并进行证实：（1）基于非线性干湿边界反演的TVDI 数据是否可靠？（2）长时间、广空间的TVDI 遥感反演是否可以采用如下流程：同年度植被生长季内不同DOY（day of year）同源影像的LST、NDVI 各自像元水平上的均值合成→边缘非数值（NaN）区域裁剪、多景空间拼接→生成各个土地利用/覆被类型的LST-NDVI 特征空间的干湿边界→反演出各地类的TVDI→合并成整幅TVDI。本研究以黄河内蒙古段为例，以1986—2020 年的步长为4～5 a 的8 期Landsat-5/8 时序影像为数据源，开展了TVDI 估算、精度评价、空间分异等方面的研究，以期为干旱半干旱地区的土壤干旱监测与预报提供借鉴。

1 研究区概况与数据来源

黄河内蒙古河段地处黄河流域最北端，起自宁夏的石嘴山，止于内蒙古伊克昭盟准格尔旗的马栅乡，全长约820 km，为典型的水沙异源、径流量季节性波动河段[10]，其中巴彦高勒至头道拐约520 km 河道为强冲积性河道（图1）。行政区划上涵盖了呼和浩特、包头、鄂尔多斯、巴彦淖尔和乌海5 座城市，占地面积约198 671 km2，年均降雨量介于80～500 mm，为干旱半干旱区域。从区位经济体量上看，该区域是黄河“几”字湾都市圈的核心地段，区内有中国“三大一首”自流引水灌溉的河套灌区，盐碱化严重。据统计，该区域2019 年底总人口达1 038.12 万人，建成区面积为597.79 km2，市区平均人口密度5 310 人/km2，人均年收入为5 438.3 元，水资源总量为10 198 760 t，相对匮乏。王国庆等研究发现黄河流域未来30 年水资源量将减少，全流域经济社会的可持续发展将受到威胁[11]。因此，亟待开展高分辨率长时间序列的地表干旱监测研究，寻找气候适应对策，指导该区域水资源的精准管控和永续利用。

图1 研究区地理信息和所用Landsat 影像行列号Fig.1 Geographical information of study area and the path and row numbers of Landsat images

研究数据来源于美国NASA Landsat 5/8多光谱与热红外影像（https://glovis.usgs.gov）。影像分带信息见图1，时间跨度为35 a（1986—2020 年），步长为4～5 a，8 期Landsat 系列卫星数据，共816 景无云影像（TM 612 景、OLI&TIRS 204 景），在影像云量覆盖度较高的情况下，用相邻年份的年积累日（DOY）接近的影像替代，NDVI 与LST 详细的反演过程见“2.1 TVDI 计算方法”部分。

与研究期一致的8 期（1980s,1990,1995,2000,2005,2010,2015 和2020 年）30 m 土地利用栅格数据集来自中国资源环境科学与数据中心（https://www.resdc.cn），1980—2020 年每月0.500°×0.625°土壤湿度数据MERRA-2(M2T1NXLND 5.12.4)，获取自：https://disc.gsfc.nasa.gov/datasets/M2TMNXLND_5.12.4/summary?keywords=https:%2F%2Fdisc.gsfc.nasa.gov；该表层土壤湿度数据为netCDF 格式，其中soil moisture L4，GWETTOP 变量用于本研究TVDI 估算结果的验证。采用ArcGIS10.2，ENVI5.3，OriginPro2015 处理数据。

2 研究方法与思路

2.1 TVDI 计算方法

图2 呈现了依据Sandholt 等[12]提出的植被指数-地表温度特征空间反演TVDI 的光学原理，数学计算过程见式(1)，在三角形UVW中，U点为干燥裸土，V点为湿润裸土，W点为湿润密闭冠层；m1表示（t-tmin），代表某一像元与相同植被覆盖情况下最湿像元的温度差；m2表示（tmax-tmin），代表在一定植被覆盖条件下最大温度差。特征空间内的斜线可看作θ相同的等值线，θ自下而上升高，斜率绝对值大的相对于斜率绝对值小的较为干旱，因此，在一定面积区域拟合出其特征空间的干边（LSTmax）、湿边（LSTmin），即可得到每个像元的干旱指数。

图2 被指数与地表温度特征空间和非线性干湿边界遥感反演TVDI 原理示意Fig.2 Vegetation index and surface temperature space and schematic diagram of TVDI remote sensing version based on the nonlinear dry and wet edges

式中：θ为温度植被干旱指数（TVDI）；t为地表温度；tmin、tmax分别为给定NDVI 值下的最小地表温度、最大地表温度。本研究中温度变量的单位均为℃。研究区海拔较高，平均1 309 m，地表起伏明显，气温差异大，地表温度受高程影响显著，需用式（2）对LST 数据进行高程订正[3,13]。

式中：tdem为修正后地表温度；t0为遥感反演的地表温度；订正系数α=0.006 ℃/m；h为30 m DEM。

本研究参照覃艺等[3]提出的自变量等间隔区间的LSTmin-NDVI 与LSTmax-NDVI 的矩阵散点数据获取方法，在ArcGIS 10.2 获取点对数据，然后在Origin Pro 2015 中绘制成2D 散点图，并根据拟合模型的决定系数（R2）显著性和实际符合情况，在3～9 次幂多项式函数中选择最佳的拟合方程，分别获取LSTmin与LSTmax的关于NDVI 的非线性方程。

地表温度通过Landsat 卫星的热红外波段反演得到，其中TM-6 采用单通道算法[14]；TIRS-10/11 采用广义单通道[14]与劈裂窗协方差-方差比率（SWCVR）[15]相结合的算法[16]；依据式（3）和（4）获得地表温度。

式中：γ与δ为中间反演参数[12,17]；ε为地表比辐射率；ψ1、ψ2和ψ3为大气参数；d为Landsat 卫星的热红外波段的数字量化值；G和O为增益与偏置值（头文件中获得）。

NDVI 是植被覆盖的一种表征，在生态和环境领域的研究中得到了广泛应用，其计算式为：

式中：v为NDVI；ρ3、ρ4分别为Landsat-5/8 的红、近红外波段的地表反射率。

依据均值合成原理，取各研究期的6—9 月份3 景LST 与NDVI 进行像元水平上的均值计算，将各自的均值合成结果作为各研究期的TVDI 估算的数据源。

2.2 水体与建设用地的剔除

水体与建设用地是地表土壤水分含量接近1 与0 的两种特殊地类，由于反演TVDI 指数的梯形或三角模型对湿度过度饱和与完全干燥的像元均比较敏感，所以待评估区域中不应含有滞留的水域和连片的不透水面[18]。借鉴Xu[19]提出的修正的归一化水体指数（MNDWI，简记为β）与Feyisa 等[20]提出的自动的水体提取指数（AWEI，简记为η），同时依据建筑指数（IBI，简记为λ）和土壤指数（SI，简记为ξ）的均值合成的归一化差异建筑与土壤指数（NDBSI，简记为μ）[21]：

式中：β为MNDWI；η为AWEI；λ为IBI；ξ为SI；μ为NDBSI；ρ1、ρ2、ρ3、ρ4、ρ5和ρ6分别为Landsat-5/8 的蓝、绿、红、近红外、短波红外1 波段和短波红外2 波段的地表反射率。计算AWEI 和NDBSI，以同期高分辨率Google Earth 和土地利用数据设定为阈值，剥除各个时期的水域与建设用地。

2.3 TVDI 估算结果验证

本研究借助估算模型精度评价指标[9,22]，即平均误差（ME，简记为Λ）、平均相对误差（MRE，简记为Δ）和均方根误差（RMSE，简记为Ω）来评价TVDI 的遥感反演结果的准确性。Λ=0，表示无偏差；Λ＞0，表示高估；Λ＜0，表示低估。Δ与Ω是对模型估算精度的度量，理论上，其值越低越好。

3 结果与分析

3.1 非线性干湿边界拟合结果

图3 中，遥感反演的8 期4 类土地利用/覆被类型下的研究区LST-NDVI 特征空间的非线性干、湿边方程结果如图所示。可见，干湿边界不是完美的线性，而是高次多项式曲线，包络呈不规则的近似梯形形状。黄河内蒙古段植物生长季（6—9 月）的LST-NDVI 特征空间的干燥边界（tmax）与湿润边界（tmin），除个别呈9 次多项式特例外，都呈7 次多项式曲线。干燥边界的R2显著高于湿润边界，且两者均达到0.05 水平上的显著，这表明拟合模型是合理的。

图3 1986—2020 年黄河内蒙古段4 类土地利用下LST-NDVI 特征空间干燥边界与湿润边界非线性拟合结果Fig.3 Non-linear fitting results of dry edge and wet edge of LST-NDVI space in Inner Mongolia section of the Yellow River under four land use types from 1986 to 2020

3.2 TVDI 时空变化

图4 呈现了8 期6—9 月平均TVDI 的空间分布格局。从空间分布规律看，依据TVDI 土壤湿度分级标准：0＜θ≤0.2，极湿润；0.2＜θ≤0.4，湿润；0.4＜θ≤0.6，正常；0.6＜θ≤0.8，干旱；0.8＜θ≤1.0，极干旱。自1986 年以来，干旱等级经历了如下变化：干旱主导→正常、干旱共存→正常兼有干旱→干旱兼有正常→干旱、极干旱平分→干旱主导、正常与湿润镶嵌。2015 年土壤干旱程度达到峰值，极干旱范围由研究区西北部的连片的未利用地（乌拉特后旗）向整个研究区扩散。2000 年出现了湿润半湿润带，主要分布在研究区东北部武川县大青山林区和黄河北岸河套灌区一带。由图4（i）可知，黄河内蒙古段土地利用/覆被以草地和未利用地为主体，近35 年来，建设农用地、林地面积呈稳步增加变化，草地面积呈先增后减变化，水体、农田和未利用地面积呈轻微波动变化。可见，未利用地（沙地、盐碱地、沼泽地、裸土地、裸岩等）是黄河内蒙古段土地干旱化的主要地类，需要加强监测和修复治理。另外，因受地形及资源环境的限制，农田的干旱化倾向，加重了对农业经济的影响，农业用水需求应考虑统筹调配。

图4 1986—2020 年黄河内蒙古段TVDI 空间分布和土地利用/覆被面积变化Fig.4 Spatial distribution of TVDI and changes of land use/cove ratio in Inner Mongolia section of the Yellow River

3.3 TVDI 估算结果验证

图5（a）呈现了8 期54 对遥感反演TVDI 与数据集（M2TUNXLND）的干燥度（1-GWETTOP）2D 散点图对比结果。本研究反演的TVD I 与（1-GWETTOP）二维散点图相对均衡地分布在1∶1 参考线（蓝线）两侧，整体表明TVDI 反演方法可行、结果准确。从图5（b）中模型估算精度评价指数折线图可知，4 类土地利用/覆被下的TVDI 在近35 年间呈先增后减变化，变化的转折时间出现在2015 年，TVDI 值大小顺序为：未利用地＞农田＞草地＞林地，林地和草地在正常与干旱之间变化，未利用地与农田在干旱和极干旱之间变化。TVDI 标准偏差呈增减波动变化，表层土壤旱情的空间异质性也呈现出高低交错变化。8 个研究期平均误差分别为：0.030 6、0.004 3、0.060 9、-0.002 6、0.078 1、-0.003 7、-0.002 1 和-0.049 4，平均相对误差（MRE）依次为：0.075 1、0.021 2、0.192 0、0.482 2、0.232 3、0.018 6、0.020 7 和0.094 2，均方根误差分别为：0.064 1、0.021 5、0.127 3、0.234 4、0.144 3、0.018 2、0.018 3 和0.104 6，总体误差为2%～8%，相对准确。这证实了本研究引言部分提出的两个科学假设，即考虑LULC、非线性干湿边界的多景Landsat 影像反演TVDI 方法合理、结果准确。

图5 TVDI 遥感反演精度评价结果Fig.5 Evaluation results of TVDI remote sensing inversion accuracy

4 结语

本研究首先提出了基于Landsat 时序数据、土地利用/覆被数据和非线性干湿边界生成长时间、广空间的TVDI 的理论框架与技术流程，然后以黄河内蒙古段为案例进行了实证研究，并将遥感反演TVDI 与NASA 发布的0.500°×0.625°月尺度的地表诊断数据集（M2TUNXLND）衍生的表层土壤干燥度（1-SM）进行比对评估，证明了本研究所提出的两个科学假设，即（1）“季节内均值合成→分用地类型拟合高次多项式干湿边界→合并生成全区域TVDI”，该方法流程是可行的；（2）基于上述方法遥感反演的TVDI 数值是可靠的。

传统的TVDI 反演模型的假设前提是土壤干燥度（1-SM）与LST 呈线性关系，且干燥、湿润边界均为线性的。本研究则是通过更高次幂多项式（7 次或9 次）研究了三者之间的关系，干湿边界的64 个拟合曲线决定系数为0.497 6～0.920 9，均达到了p＜0.001 水平上的显著性要求，并获得了可靠的预测结果。这表明，考虑LST-NDVI 特征空间的干湿边界的非线性关系，在一定程度上可以改善区域性TVDI 的预测精度。

在利用TVDI 或改进TVDI 进行干旱分析时，大多数研究人员将研究区域作为一个整体考虑，较少考虑到不同的土地利用/覆被的影响。由于林地、草地、农田和未利用地之间的植被高度和盖度存在较大差异，相应的地表最高和最低温度也存在较大差异，导致TVDI 计算结果存在较大差异。同种类型土地利用/覆被提供了相对均质的地表环境，有利于准确表达TVDI 与LST 和NDVI 之间的关系，有效降低了TVDI 遥感反演的环境差异误差。本研究所得结论与SHI 等[7]的结果一致。

在后续的土壤旱情预测预报研究中，需要深入地开展干燥、湿润边界的线性与非线拟合效果评估研究，同时兼顾地貌、气候和土壤质地等因素。此外，应考虑卫星遥感数据的分辨率与下垫面的异质性相匹配。为了贯彻黄河流域生态保护与高质量发展重大国家战略，需要以水定人、以水定产，精准核算黄河“几”字湾区域的水资源承载力，助力该区域都市群高质量发展，实现黄河大保护生态目标。