基于TVDI 的旱情时空动态变化监测——以神东矿区为例

邵天意，包斯琴*，王 楠，韩阿茹汗

（1.内蒙古农业大学 沙漠治理学院，呼和浩特 010020；2.中国地质大学（北京）土地科学技术学院，北京 100089）

0 引 言

【研究意义】干旱是自然生态系统产量的重要输入因子，同时也是植被发挥蒸腾作用及光合作用的限制因子，其在一定程度上能够决定一个地区的植被类型和植被生长结构[1]，更是衡量荒漠化程度以及指导生态恢复措施实施的重要影响因子，对准确且快速的掌握旱情状况对土地复垦、生态修复和农业生产具有指导意义[2]。【研究进展】传统长时间序列的旱情监测方法需要消耗巨大的人力、物力和财力。相比于传统的监测方法，遥感监测具有监测范围广、时效性长等特点[3]。诸多学者[4-6]针对遥感监测领域进行的大量研究表明，表观热惯量模型进行干旱遥感监测需要地表反照率和太阳辐射作为参数，其更适合于裸土的研究区，同时还需考虑土壤类型变化对监测结果的影响。基于冠层温度的遥感监测方法在实际监测中易受自然气候条件，卫星传感器等影响会产生不确定的误差[7-8]。综合作物长势进行干旱监测以植被长势状态表征干旱程度存在相对大的滞后性，预警时效较差。微波遥感监测旱情对卫星影像数据的精度要求极高[9-10]。综上所述，利用研究区归一化植被指数（Normalized Difference Vegetation Index，NDVI）与地表温度（Land Surface Temperature，LST）建立温度干旱植被指数（Temperature Drought Vegetation Index，TVDI）干旱监测模型能够消除了单一因素对监测结果的影响，更具时效性，且该模型监测适用于中高植被覆盖地区[11]。TVDI模型是一种采用光学和热红外遥感通道数据进行地表水分的遥感反演方法[12]。对于植被覆盖区而言，地表水分一定程度上决定着植被的冠层温度，在一定条件下，植被的冠层温度能间接反映植被的供水状况[13-15]。据此，Sandholt等[16]基于NDVI与LST建立了Ts-NDVI特征空间，计算得到的TVDI能够很好地反演出地表水分，其反演精度也得到了验证[17]。诸多学者[18-20]也基于不同数据源利用TVDI模型进行地表水分反演并进行精度验证，发现TVDI在0～10 cm 土层深度具有良好的反演效果。【切入点】国内基于TVDI模型进行地表水分动态监测的研究很多，但针对地势复杂沟壑纵横的矿区进行的研究尚不多见。神东矿区受地理位置的影响干旱程度高，加之矿区内常年进行煤炭开采作业，导致地表结构被破坏，土壤水分流失严重，干旱程度加剧。【拟解决的关键问题】本文基于TVDI模型，利用神东矿区1991—2018 年植被生长季遥感数据进行干旱程度动态监测，以期为研究区的生态恢复和土地复垦提供指导建议。

1 材料与方法

1.1 研究区概况

神东矿区位于陕西省榆林市北部和内蒙古自治区鄂尔多斯市南部（109°51′—110°46′E，38°52′—39°41′N），区内大部分为典型的风成沙丘及沙滩地地貌。海拔800～1 385 m，沟壑纵横，地表支离破碎，水土流失十分严重。神东矿区地处毛乌素沙地与黄土高原的过渡地带，属半干旱大陆性季风气候，温度年较差较大，年平均气温6.6 ℃，1 月平均气温-10.1 ℃，7 月平均气温20.5 ℃。区内干旱少雨，降水集中，汛期为6—9 月，占全年降水量的76%，降水年际变化较大，最大年达919.10 mm，最小量年仅为108.60 mm。区内以沙生植物为主，植被每年4 月返青，10月叶落，郁闭较差。区内风沙土（Aeolian sandy soil）占矿区总面积的50%，黄土性土（Loess soil）占总面积的30%，红土性土（Red native soil）占总面积的10%，土壤具有质地较粗，结构不良，肥力低，抗侵蚀性差的特点。研究区位置如图1 所示

图1 研究区位置示意图Fig.1 Schematic diagram of study area location

1.2 数据来源

遥感影像数据来源于地理空间数据云（http://www.gscloud.cn/）中的美国国家航空航天局（National Aeronautics and Space Administration,NASA）发射的陆地系列卫星（Landsat）所搭载的30 m 空间分辨率数据，包括 Landsat5-TM 和Landsat8-OLI，时间分辨率为16 d 的2 种传感器。数据时相选取神东矿区1991、2002、2007、2010、2014、2018 年植被生长季（5—10 月）共36 景影像，在1991—2018 年数据质量上均选择云量小于10%反演效果好的数据单元，借助ENVI 5.3 进行辐射定标、大气校正、镶嵌等预处理，然后利用数据估算得到空间分辨率30 m 的NDVI及LST数据集。

1.3 研究方法

1.3.1TVDI计算方法

以提取得到研究区的归一化植被指数（NDVI）为横坐标，地表温度最大值和最小值为纵坐标进行线性拟合，形成关于Ts-NDVI的三角形特征空间，将该特征空间的斜率作为TVDI的等值线，TVDI自下而上升高，TVDI的绝对值越大，干旱程度越大。因此，在研究区拟合出其特征空间的干湿边方程，即可得到每个像元的干旱指数，其计算式为：

式中：TVDI为温度干旱植被指数；Ts为地表真实温度；Ts_max为地表最高温度；Ts_min为地表最低温度。TVDI值越大，代表地表水分越低，干旱程度越大。TVDI的干湿边拟合方程为：

式中：a1、b1、a2、b2分别为干湿边线性拟合方程的系数。

在干湿边方程中，其斜率在生态学意义上代表着土壤含水率在饱和与不足时的地表温度值，将各像元的NDVI值带入所对应的干湿边方程，计算出最高地表温度（Ts_max）与最低地表温度（Ts_min），利用式（1）计算各像元对应的TVDI值。本文采取适用干旱半干旱区的TVDI分级标准进行分级，如表1 所示[21]。

表1 TVDI 分等定级Table 1 TVDI classification and classification

1.3.2 偏差分析法

本文采用偏差分析法[22]分析研究区1991—2018年植被生长季的TVDI时间分布特征。偏差分析法可以在时间尺度上通过某年TVDI与多年平均的TVDI距离，表示某一时期TVDI偏离多年平均TVDI的程度。TVDI偏离值为正值，表示年均TVDI高于多月平均水平，TVDI在时间序列上呈上升趋势，TVDI偏离值为负值，表示年均TVDI低于多年平均水平，TVDI在时间序列上呈下降趋势。

1.3.3 Mann-Kendall 检验法[23]

Mann-Kendall 检验法（M-K 检验法）是世界气象组织推荐并广泛使用的非参数检验方法。其不需要服从特定的分布，亦不受样本值的影响，因此广泛应用于分析一般数据在时间序列上的趋势检验和突变点检验。对于时间变量（x1,x2, …,xn），n为时间序列的长度，M-K 法定义了其统计量S：

式中：Sgn（）为函数符号，规则如下：

S为正态分布，其均值为 0 ，方差Var(s)=(n-1)(2n+5)/18，当n＞10 时，正态分布统计量计算为：

若Z＞0，表明TVDI在时间序列上呈上升趋势，若Z＜0，反之，绝对值越大，趋势越明显；且Z≥1.28、1.96、2.32 时，分别通过了90%、95%、99%水平的信度检验，在该水平上显著。UF曲线与UB曲线的交点，为TVDI在时间序列上的突变点。

1.3.4 趋势分析和F检验

线性倾向率（Bslope）能够反映TVDI在空间上随时间变化的的上升或下降趋势。因此，本文逐像元计算1991、2002、2007、2010、2014、2018 年的TVDI均值，得到相应年份的TVDI空间分布，并利用线性倾向率计算研究区1991—2018 年植被生长季的TVDI在空间上的线性变化趋势。其计算式为：

式中：Bslope为线性倾向趋势；i为36 个月变量（i=1,2, 3, …, 36），n=36，TVDIi为第i月的TVDI值。若Bslope＞0，表示干旱程度在空间上随时间变化呈上升趋势；Bslope＜0，表示研究区干旱程度在空间上随时间变化呈下降趋势。用F检验对TVDI的空间变化趋势进行显著性分析，根据F检验显著性临界值（α=0.05），将F值划分为显著下降、显著上升、不显著下降、不显著上升4 种趋势。

1.3.5 空间转移矩阵

空间转移矩阵来源于系统分析中对系统状态与状态转移的定量描述，可以定量识别TVDI不同等级在某一时间间隔的空间格局变化，不仅可以反映TVDI不同等级的面积变化，还可以直观反映TVDI各等级面积转入转出情况[25]。

2 结果与分析

2.1 神东矿区植被生长季TVDI 时间变化特征

2.1.1TVDI月际变化特征

利用ArcGis10.7 像元统计工具，统计TVDI值，得到神东矿区1991—2018年植被生长季TVDI月际变化如图2 所示。由图2 可知，神东矿区1991—2018年TVDI月际变化非常明显，月均TVDI在0.35～0.85之间波动。2007 年6 月TVDI均值最高为0.85，说明地表干旱程度最大。2018年9月TVDI均值最低为0.39，说明地表干旱程度最小。月均TVDI在每年的6—9 月出现最低值，这是由于神东矿区降水集中导致的，与神东矿区6—9 月为汛期的自然条件相吻合。

图2 神东矿区1991—2018 年植被生长季TVDI 月际变化Fig.2 Intermonthly TVDI variation of vegetation growing season in Shendong Mining Area from 1991 to 2018

2.1.2TVDI突变分析

采用M-K 突变检验法对神东矿区1991—2018 年植被生长季的月均TVDI进行检验，如图3 所示。图3 检验结果表明，UF～UB曲线几乎处于显著性水平α=0.05 置信区间内，TVDI值呈下降趋势（UF＜0），与神东矿区月际变化规律相吻合。1991 年7—10 月存在4 个突变点，2014 年9 月存在1 个上升趋势突变点，2018 年5 月、7 月左右存在一个上升趋势、一个下降趋势突变点。

图3 1991—2018 年神东矿区植被生长季TVDI M-K 检验Fig.3 TVDI M-K test of vegetation growing season in Shendong Mining Area from 1991 to 2018

2.1.3TVDI年际变化特征

计算得到神东矿区各年植被生长季TVDI均值如图4 所示。从年际变化来开，年均值波动于0.56～0.69之间，按照地表干旱等级类型划分，神东矿区干旱程度基本处于正常、较干旱等级，这与TVDI的空间分布相吻合。在研究时段内，2002 年TVDI达到峰值为0.686，自2002 年后TVDI逐年下降，2018 年达到最低的0.561，下降速率为0.02/30 a。该结果在实地调研中得到了验证，这是由于神东矿区在煤炭开采后进行了大量的土地复垦工作，水土流失强度减弱，种植大量的防风固沙植被，神东矿区地表涵养水源的能力得到提升，因此地表水分上升，地表干旱程度呈下降趋势。

图4 神东矿区1991—2018 年植被生长季TVDI 年际变化Fig.4 Interannual TVDI variation during the vegetation growing season in Shendong Mining Area from 1991 to 2018

2.1.4TVDI偏差分析

利用偏差分析法分析神东矿区1991—2018 年年际TVDI变化趋势，如图5 所示。神东矿区1991、2002、2007、2010 年TVDI偏离值＞0，说明1991—2010 年TVDI均值呈上升趋势，干旱程度具有增大趋势。其中1991—2002 年TVDI偏离程度最大，上升趋势最明显，高于TVDI均值0.05。2014、2018 年TVDI偏离值＜0，说明2010—2018 年TVDI呈下降趋势，干旱程度减小。2014—2018 年偏离程度最大为0.56，低于TVDI均值（0.08）。该结果与神东矿区TVDI年际变化趋势一致。

图5 神东矿区1991—2018 年TVDI 变化趋势Fig.5 TVDI trend of Shendong Mining Area from 1991 to 2018

2.2 神东矿区TVDI 空间变化特征

2.2.1 神东矿区植被生长季TVDI干旱等级空间分布

神东矿区1991—2018 年植被生长季干旱等级的空间分布如图6 所示。1991—2002 年研究区西部干旱等级下降，东部区域的干旱等级上升。2002—2007年研究区中部、南部干旱等级上升。2007—2010 年研究区中部、西部干旱等级基本达到正常级别，东部地区干旱等级上升。2010—2014 年东部区域干旱等级下降，西部区域干旱等级略有提升。2014—2018年研究区西南部干旱等级增大，中部及东北部的干旱等级有小范围的下降。从干旱等级划分来看，神东矿区干旱等级始终表现为西南高于东北，其主要原因为神东矿区的降水量从西南到东北呈递减趋势，空间上基本处于正常、较干旱和干旱3 种等级。

图6 神东矿区1991—2018 年植被生长季干旱等级的空间分布Fig.6 Spatial distribution of drought grades in the vegetation growing season from 1991 to 2018 in Shendong Mining Area

2.2.2 神东矿区线性倾向趋势分析

利用ArcGIS10.7，基于神东矿区1991—2018 年植被生长季TVDI的空间分布，运用式（7）计算其线性倾向率（Bslope），分析其线性倾向趋势，并对其线性倾向趋势空间分布在α=0.05 显著性水平上进行F检验，如图7 所示。由图7 可知，神东矿区1991—2018年TVDI线性变化趋势下降趋势面积大于上升趋势面积，Bslope＜0 的区域集中在研究区的大部分地区，该区域TVDI呈下降趋势，地表干旱程度变小。Bslope＞0的区域零散分布于研究区的北部、东部、南部，这些区域TVDI呈上升趋势。地表干旱程度变大。根据F检验的空间分布得知，在TVDI呈下降趋势的区域中，大部分地区下降趋势显著，下降趋势不显著的地带则零星分布于其中。在TVDI呈上升趋势的区域中，上升显著的地区大部分集中在神东矿区的东北部，该区域内分布着神东矿区的矿井群，这是由于采煤后塌陷，造成植被破坏严重，地表水分下渗速度加快，地表水分减小，所以干旱程度增大。神东矿区TVDI空间变化结果与时间变化的趋势一致

图7 1991—2018 年神东矿区植被生长季TVDI 线性变化趋势空间分布及F 检验Fig.7 Spatial distribution and F-test of TVDI linear trend in vegetation growing season in Shendong Mining Area from 1991 to 2018

通过对神东矿区TVDI线性倾向趋势面进行计算发现：神东矿区1991—2018 年植被生长季TVDI下降趋势面积为3 788.43 km2，上升趋势面积为1 215.33 km2，下降趋势面积大于上升趋势面积。从F显著性检验结果来看，神东矿区显著下降面积最大为2 394.74 km2，显著上升面积最小为425.91 km2，说明研究区的地表干旱程度得到缓解，该地区生态恢复效果显著。

2.2.3 神东矿区TVDI转移矩阵

利用ArcGIS10.7 空间分析工具，选择1991、2002、2010、2018 年的TVDI干旱等级分布进行空间转移矩阵的计算，得到神东矿区1991—2018 年植被生长季干旱等级转移情况，如图8 所示。由图8 可知，1991—2002 年，整个研究区的西北、西南、东南大范围区域由干旱变化为其他4 种干旱等级。较湿润、正常、较干旱变为干旱零星分布于中东部，该部分区域旱情有所加重，其他区域基本保持原有的干旱等级。2002—2010 年，研究区大部分区域干旱等级由较干旱向正常、较湿润等级转化，只有西南部小范围由干旱变为其他等级，东部由较湿润、正常、较干旱转变为干旱。2010—2018 年，研究区中部、东部干旱等级由干旱转为其他等级，旱情程度下降，中部及北方大部分干旱等级则由较干旱变为正常，西南小部分区域由较干旱转换为干旱，这一时段研究区旱情转移基本表现为由干旱变为较干旱、较干旱变为正常，研究区大部分旱情得到一定程度的缓解。

图8 神东矿区1991—2018 年植被生长季TVDI 转移矩阵Fig.8 TVDI transfer matrix of vegetation growing season from 1991 to 2018 in Shendong Mining Area

综上所述，神东矿区在研究区间内，干旱等级变化均有不同程度的波动，但在空间上基本表现为由干旱等级向较干旱、正常等级转变，较干旱等级向正常、较湿润等级进行转变，说明神东矿区进行的一系列生态恢复及水土保持措施达到了一定成效。神东矿区地处我国西北内陆，空间跨度较大，因此水热状况差异较大，且随着地表与地下的采煤扰动，地表覆盖变化复杂，加之自然因素及人类活动的影响导致了地表水分在时间上和空间上变化明显[26]。

3 讨 论

本文基于TVDI模型对神东矿区1991—2018年植被生长季的旱情进行监测。从TVDI的月际变化来看，神东矿区每年植被生长季内6—9 月TVDI出现最低值，该结果与神东矿区降水集中的特点相吻合[27]。TVDI的年际变化显示，神东矿区自1991—2018 年，TVDI整体呈下降趋势，下降速率为0.02/10 a，地表水分逐年上升，旱情逐步得到缓解，该结果与刘英等[28]基于梯度结构相似度的神东矿区土壤湿度空间分析得到的神东矿区土壤湿度在空间上呈增加趋势，旱情缓解的结果一致。对TVDI时间变化进行趋势分析发现，神东矿区2018 年TVDI下降趋势最为显著，这是由于神东矿区在2018 年前后实施了大量的水土保持措施，在不考虑自然因素对各年TVDI变化影响的情况下，研究结果与实地调查结果一致。神东矿区1991—2018 年干旱等级的空间分布表明，研究区干旱程度，始终为西南部高于东北部，这是由于该地区西南部地处毛乌素沙地边缘[29]。干旱等级面积减少，正常及较湿润面积增加，该结果与TVDI的时间变化特征相吻合。从线性倾向趋势的空间分布来看，研究区整体TVDI下降趋势面积远远大于TVDI上升趋势面积，且TVDI呈上升趋势的区域集中在神东矿区的矿井群附近，这是由于采煤后地表塌陷，地表发生变形，地表水下渗速度加快，加快该区域的水土流失[30]。利用转移矩阵分析发现，整个研究区在1991—2018年间，干旱等级向较干旱、正常转化，而较干旱等级大多向正常、较湿润转化，旱情逐步得到缓解。

本文是在保证TVDI模型对神东矿区旱情监测具有较高的精度下进行研究，研究发现的TVDI在时间序列中存在多个突变点问题，还需要在后续研究中对突变点的变化情况进行相关验证。下一步还需要对造成其发生时空演变的自然及人为因素进行更深一步的研究，以期为神东矿区的水土保持措施及生态恢复提供参考。

4 结 论

1）1991—2018 年神东矿区植被生长季TVDI在0.35～0.85 之间波动，每年6—9 月出现低值。TVDI月际变化整体呈下降趋势。

2）年际变化表明，1991、2002、2010 年TVDI偏离值＜0，TVDI呈下降趋势，2014、2018 年TVDI偏离值＞0，TVDI呈上升趋势。TVDI在时间序列上整体呈下降趋势，下降速率为0.02/10 a。

3）神东矿区1991—2018 年植被生长季TVDI干旱等级空间分布，始终表现为西南高于东北。研究区干旱等级面积逐年减少，正常、较湿润等级面积逐年增加。

4）研究区TVDI下降趋势面积远大于TVDI上升趋势面积。显著下降趋势面积最大为2 394.74 km2，显著上升面积最小为425.91 km2。

5）1991—2018 年，TVDI的干旱等级逐渐向较干旱、正常等级转变。TVDI的较干旱等级逐渐向正常和较湿润等级转变。向干旱及较干旱等级转变的范围较小。

（作者声明本文无实际或潜在的利益冲突）