基于MODIS数据的山东省2014-2016年干旱监测分析

王正东, 郭 鹏, 万 红, 杨 纲

(山东农业大学 信息科学与工程学院, 山东 泰安 271018)

干旱是我国乃至世界上许多国家主要的自然灾害之一，不仅会对生态系统和环境造成极大破坏，还会严重影响社会经济活动以及居民生活。在全球变暖的大趋势下，干旱的发生变得越来越频繁，干旱监测也成为全球关注的一个科学热点问题。山东省地处黄淮海地区，是主要的产粮大省，受季风影响明显，属于水资源短缺省份，近几年，出现了持续的少雨及异常气候事件频率增多的现象，受旱灾影响严重。因此，如何及时、有效、准确地监测黄淮海地区干旱的发生、旱情的发展已成为有关部门关注的焦点。

常规的土壤水分监测手段通过人工测墒完成，耗时、耗力。利用遥感卫星监测干旱则具有周期短、范围广等优点，通过对研究区的遥感影像进行分析，可以快速获得干旱数据，确定干旱等级，从而为防灾减灾工作提供依据[1]。温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)是遥感监测干旱的常用方法之一，利用该方法可以及时、准确和有效地监测干旱情况，能较好地监测不同地表面，且具有计算较简单，容易实现等优点[2]，许多学者已进行过相关研究。如齐述华等[3]利用NOAA-AVHRR数据，采用温度植被干旱指数法反映表层土壤水分变化情况。王婷婷等[4]利用趋势线法对2002—2009年每年8月份的TVDI值进行回归分析，来研究2002—2009年松辽平原的干旱变化趋势。

本文以山东省为研究区，利用MODIS的NDVI，EVI，LST产品，在使用Savitzky-Golay(S-G)滤波方法进行重构，填补缺失数据的基础上，分别构建NDVI-LST和EVI-LST特征空间，建立基于时间序列的温度植被干旱指数对山东省2014—2016年干旱变化的时空演变特征进行分析，并探讨其与气象因子之间的关系。

1 数据与方法

1.1 研究区概况

山东省位于中国东部沿海，地处114°20′—122°43′E， 34°23′—38°23′N，地幅辽阔，陆地南北最长距离约420 km，东西最宽约700 km，全省总面积为15.8万km2，境内最高处位于泰山，海拔为1 524 m，最低处位于东北部的黄河三角洲，海拔为2 m左右。山东的气候类型属暖温带季风气候，降水集中，四季分明，冬夏较长，全年平均温度为10～14℃，年平均降水量为550～950 mm，但是降水季节分布不均衡，对农业生产影响较大，因此寻求一种大范围、实时的旱情监测手段对于山东省的社会经济发展具有现实推动意义。

1.2 数据来源与预处理

本研究所使用的MODIS遥感数据从MODIS官网(https:∥ladsweb.modaps.eosdis.nasa.gov/)下载得到，主要有8 d合成的地表温度产品(MOD11A2)和16 d合成的植被指数产品(MOD13A2)，分辨率为1 km，时间范围为2014年1月—2016年12月，MOD11A2和MOD13A2数据的原始格式都为EOS-HDF。在对数据预处理时，首先利用MRT(MODIS Reprojection Tools)对两种产品数据进行格式转换、拼接、重投影等批处理，然后根据山东省的行政区划边界裁剪得到研究区影像，最后得到研究区3年的归一化植被指数产品(NDVI)数据、增强型植被指数产品(16 d合成数据)和陆地表面温度数据(8 d合成数据)以及相应的质量控制文件。

温度和降水等气象数据来源于中国气象科学数据共享服务网(http:∥cdc.nmic.cn/)提供的中国地名气候资料月值数据集，本文选取惠民、定陶、济南、兖州、潍坊5个国家级气象站的每月月平均降水量和月平均温度数据，起止时间为2014—2016年，研究区气象站点的空间分布如图1所示。

图1 研究区范围与气象站点位置

1.3 数据重建

在获取遥感影像时，因为存在传感器设计、大气传输、太阳光照角度、观测视角等随机干扰因素，影像中会存在很多噪声，当有云覆盖时，会造成数据的缺失，所以获取的影像数据有时无法准确表达地表情况，因此有必要进行数据重建。由Savitzky等[5]提出的移动窗口的加权平均算法S-G滤波重建能够较好地完成数据平滑，内插，反映真实的植被生长、气象变化等地表情况，适用于NDVI，EVI等数据的重建[6]。S-G滤波是一种特殊的低通滤波器，用于数据的噪声平滑，具有实现简单、先验知识少等优点，可以较真实地体现植被的生长变化情况[7]。

S-G滤波是以平滑时间序列数据为移动窗口的加权平均算法[8]，该方法中的加权系数取决于一个滤波窗口内给定高阶多项式的最小二乘拟合次数。若某个数据集内包含n个点，半窗口的宽度为m，滤波器的长度等于滑动数组的宽度等于N。Ci为平滑窗口内对应每个点的卷积系数，则对于第j个点来说，平滑后的值即为：

(1)

式中：Yj*为拟合后的结果。

根据像元可信度值(pixel reliability)和权重选择质量控制(quality control)文件分别对影像的每个像元赋予权重，通过编程分别对NDVI，EVI和LST进行S-G时序滤波处理。除因时间而异的观测条件会使时间序列呈现类似于锯齿状的不规则波动，在其中的某一次观测时由于不同设备间的电磁干扰，或者观测设备的位置偏移等原因会产生高频噪声与干扰，所以还需要进行空间滤波，目的是提高单次观测获取的影像的质量，因此本研究使用中值滤波消除空间上的异常值。

1.4 数据重建效果分析

图2 S-G滤波重构值、站点实测值与原始值的对比

1.5 植被指数-地表温度特征空间

Lambin等[9]指出了地表温度与归一化植被指数两者存在的关系，建立TVDI模型对于消除土壤背景的影响效果较好。如果研究区是从裸土到植被完全覆盖，土壤湿度从极湿润变化到极干旱时，以LST和VI为纵横坐标构建的散点图呈三角形。植被覆盖区的TVDI理论取值应处于(0，1)之间[10]。TVDI值越低则表明干旱程度越轻，反之TVDI值越高，则干旱程度越严重。TVDI由植被指数(NDVI/EVI)和地表温度(LST)计算得到，其定义为公式:

(2)

式中：TS为在一定的分辨率条件下任意像元的地表温度；TSmin表示某一NDVI对应的最低地表温度值，即湿边；TSmax为某一NDVI对应的最高地表温度值，即干边。

TSmax=a1+b1×NDVI

TSmin=a2+b2×NDVI

(3)

式中：a1，b1是干边拟合方程的系数；a2，b2是湿边拟合方程的系数[10]。由地表温度最大值拟合而成的称为干边，最小值拟合而成的则称为湿边。由TVDI的原理可知，地表温度最大值与植被指数呈负线性关系。若植被覆盖度小于20%，低植被覆盖度对应的NDVI不能准确反映出区域内植被生长情况；若植被覆盖度大于80%时，其对应的NDVI增长速度会逐渐减慢并趋于平缓，出现一种饱和状态，并且对于植被监测的灵敏度也下降。由此得出，NDVI更适用于中等植被覆盖率的情况。因此，在拟合特征空间中每一步长NDVI对应LST最大值与最小值时，选择处于0.2～0.8范围内的VI值[11]。

2 结果与分析

2.1 TVDI计算

在VI为0.2～0.8的范围内，以0.01为步长，提取每一步长NDVI，EVI像元对应的所有LST像元中的最大值与最小值，并采用最小二乘拟合方法对每期特征空间的干湿边进行线性拟合，得到干湿边的拟合方程及其相关系数，结果见表1—2。从干湿边的拟合结果来看，NDVI与LST的最大值呈负相关关系，与最小值呈正相关关系。EVI-LST构成的特征空间干湿边斜率与NDVI-LST一致，但EVI与LST构建的特征空间中干湿边拟合的相关系数更高，主要因为EVI加入了蓝色波段来增强植被信号，消除土壤背景和气溶胶散射的影响，因此能够克服NDVI在高植被覆盖地区增长缓慢趋于平缓，容易饱和，在低植被覆盖地区容易受到土壤植被影响等缺点。进入7月份后，相关系数均降低，原因可能是进入夏季后，山东大部分地区地表温度普遍较高，因此随着NDVI，EVI的变化，地表温度变化较小。

图3为2015年1月济南市NDVI，EVI与LST构建的散点图，由图3可知，相比NDVI，EVI总体偏低，适用于高植被覆盖区域，在构建特征空间时，拟合的干湿边更稳定，趋势更明显，干湿边的相关系数分别达到0.84，0.86，而NDVI的干湿边相关性仅为0.79，0.74，因此，本文选择EVI和LST构建TVDI模型进行下一步的研究。

表1 NDVI-LST干湿边拟合方程及其相关系数

表2 EVI-LST干湿边拟合方程及其相关系数

2.2 研究区干旱状况时空分析

根据EVI-LST计算得到山东省2014—2016年每个月TVDI图像，再以3—5月份为春季，6—8月份为夏季，9—11月份为秋季，12月—翌年2月份为冬季，制作每个季度平均TVDI图像。然后按照整个研究区域内像元TVDI值的分布直方图和中国土壤湿度界定干旱的标准对研究区干旱情况进行分级[12]，将干旱划分为5个等级，分别是潮湿(0≤TVDI≤0.2)、湿润(0.2≤TVDI≤0.4)、正常(0.4≤TVDI≤0.6)、干旱(0.6≤TVDI≤0.8)、重旱(0.8≤TVDI≤1.0)，分级结果如图4所示。

为了分析山东省2014—2016年土壤水分时空分布变化趋势，统计山东省每年各干旱等级土壤面积以及百分比(表3)。从分析结果看，2014年上半年旱情较严重，进入秋季之后，重旱与干旱区域的面积显著下降，表明干旱得到一定程度的缓解；2015年旱情为3年中最严重的一年，重旱与干旱区域的面积比例分别高达7.46%，34.24%，表明在一年中山东省有近一半面积处于干旱状态，全年各干旱等级土壤面积比例与2014年基本相一致，但是2015年是在进入冬季之后，干旱才得到减弱，东北部地区表现为相对湿润，其他地区基本属于正常状态，只有局部地区存在干旱情况；2016年正常和湿润区域面积占比分别达到44.52%，42.54%，相对于前两年，干旱面积显著下降。

图3NDVI，EVI-LST特征空间干湿边方程

图4研究区3年旱情等级时空变化

表3 不同干旱等级面积及比例

2.3 TVDI与气象数据的相关性分析

湿润指数中涉及多种气象因子，为了进一步探究气象因子与TVDI变化的关系，本文研究了主要气象因子气温和降水的年际变化及其与TVDI的关系。利用山东省的济南、兖州、潍坊、定陶和惠民5个国家级气象站收集的每月平均温度和月均降水量与TVDI进行相关性分析，根据温度、降水量与TVDI的相关系数(表4)可以得出，TVDI与温度均通过0.01水平的显著性检验，表明相比于降水，TVDI对温度更为敏感，在温度升高时，作物生长需水量大，地表水分蒸发速度加快，蒸发量大于降雨量，土壤含水量降低，使得地表较为干旱，导致TVDI值升高；TVDI与降水也均通过0.05水平的显著性检验，且TVDI与降水量存在显著负相关关系，当存在连续降水时，TVDI的值会呈降低趋势，但两者并不同步，TVDI的变化存在一定的滞后性，当降水量变化稳定时，TVDI也随之趋于稳定。以潍坊市为例，2014—2016年3年中的降水量明显少于其他区域，从气象站统计3年年均温度与年均降水量中可以看出3年年平均降水量仅为378，435，478 mm，相比于同时期其他地区，明显偏少，所以潍坊在3年中一直处于TVDI偏高的情况。

2015年降水量相比于2014年有所增加，但是受到“历史最强”厄尔尼诺的影响，干旱情况更加严重，2015年TVDI与温度的相关系数较高，除定陶外均达到0.76以上，而同年降水与TVDI的相关系数仅有0.5左右，最高仅达到0.69，因此温度成为2015年TVDI值升高、干旱加重的主要影响因素，所以即使在降水量增加的情况下，干旱程度并没有减弱反而增强。2016年年平均降水量明显升高，潍坊、定陶、兖州等地降水量均达到3年最高水平，TVDI反演的季度干旱情况也表明全省干旱情况有所减弱。

表4 温度、降水量与TVDI的相关系数

注：*表示达到0.01的显著性水平；**表示达到0.05的显著性水平；剩余表示未通过显著性检验。

3 结论与讨论

(1) 利用S-G滤波与空间滤波对植被指数与地表温度进行重建，可以有效消除卫星数据在时间和空间上存在的缺失值和噪声的影响，对滤波后的温度值与实测值做相关性分析，相关系数均在0.91以上，说明经过重构后的数据与站点实测数据在趋势上基本一致，因此可以利用重构后的数据进行TVDI模型的构建。

(2) NDVI/EVI-LST特征空间的干边随着植被指数的增大而呈递减趋势，而湿边在气温较低月份与植被指数关系明显，另外NDVI与EVI在构建特征空间上也有所不同，NDVI与LST构成的特征空间更平缓；而EVI与LST构成的特征空间范围更集中，在干湿边拟合时，趋势更稳定，相关性更高，因此运用EVI-LST构建TVDI更为合理。

(3) 从时间上看，山东省2014年干旱主要集中在春季和夏季，全年平均干旱面积占到全省面积的37.62%，在进入春季之后，温度迅速回升，植被生长迅速，蒸发旺盛，干旱面积显著提高，进入8月之后，除局部发生重旱外，大部分地区的旱情得到缓解。2015年干旱分布情况基本与2014年一致，但持续时间更长，全年平均干旱面积比例增加，占全省面积的41.7%。只有当进入12月份之后，干旱才减弱。2016年由于降水量显著增多，只有在2,3月份有一次范围较大的旱情，其余月份除局部地区外全省大部分地区相对湿润。

(4) 从空间上看，2014年和2015年年均重旱面积比例分别达到4.93%，7.46%，其中鲁中与半岛交界处、鲁西南地区是干旱较为严重的地区，根据3年的变化情况分析得出，干旱分布体现为由东南向西北转移的趋势。湿润以及正常的区域主要分布在中部的山区、最北部的平原地区以及胶东地区。

(5) 对TVDI与降水量、温度数据做相关性分析，均达到显著相关，并发现TVDI与温度的相关性更强，且通过0.01水平的显著性检验，表明改进的TVDI指数能够较好地反映山东省2014—2016年的旱情变化情况，对研究区旱情的快速准确检测和干旱演变过程的研究具有一定的参考价值，有助于制定相应的防灾和减灾决策。

受限于MODIS数据产品本身的问题，本文仍存在一些不足之处：温度是TVDI的影响因子之一，而高程的变化是影响地表温度的重要条件。山东省地处丘陵地带，境内高程变化较大，因此还应该考虑高程对温度变化的影响，利用高程对地表温度进行校正从而计算TVDI。另外TVDI与温度、降水数据与有些气象站点的相关关系较弱，其原因可能是：① 由于范围不同，整个行政区观测数据与3 km×3 km缓冲区数据的相关性不可能达到很高的水平；② 虽然气象资料是准确的，但参与TVDI建模的地表温度和植被指数都是按固定期数合成的数据，尽管假设两者在时相上是相一致的，但实际上可能并不完全一致，存在部分偏差；③ TVDI与气象因子对干旱的反映侧重点不同，TVDI是地面植被、土壤等受干旱影响的响应程度，气象因子则比较直观地表现地面湿润情况。