归一化水体指数用于河南省干旱监测适用性分析

2020-02-03谷佳贺薛华柱董国涛程结海

干旱地区农业研究 2020年6期

谷佳贺，薛华柱，董国涛，程结海

(1.河南理工大学测绘与国土信息学院，河南焦作 454003；2.黑河水资源与生态保护研究中心，甘肃兰州 730030)

干旱是一种发展缓慢的自然灾害，是我国乃至世界上许多国家主要的自然灾害之一，对生态系统和社会经济以及居民生活造成严重的影响[1]。在全球变暖的背景下，干旱有加剧态势，因此干旱监测也成为全球关注的重点问题[2-3]。卫星遥感具有覆盖范围广、持续时间长等特点[4]，是近些年来较有效的技术手段之一[5]，可用于监测大范围的土壤湿度和植被生长状况[6]。MODIS(moderate-resolution imaging spectroradiometer)数据具有高时间分辨率、高光谱分辨率、适中空间分辨率的特点[7]，近年来在旱情监测中得到了广泛的应用[8]。

国内外学者利用MODIS数据构建了不同的指数对干旱进行监测，彭擎等[9]采用MODIS归一化植被指数(NDVI)和地表温度(LST)数据分析了新疆2000—2015年生长季3个阶段NDVI与LST的时空变化特征及相关关系；Seo-Yeon Park等[10]利用MODIS地表温度(LST)、植被健康指数(VHI)、蒸散发以及降水数据应用于干旱监测并评估了其适用性，确定了不同时间尺度标准化降水指数(SPI)与干旱指数之间的相关性，可间接应用于农业或水文干旱监测；刘英等[11]以MODIS NDVI和LST数据构建了温度植被干旱指数(TVDI)模型并分析了2000—2016年间陕西省旱情时空分布特征和规律。Zhang等[12]根据我国2000—2014年的MOD13A3 NDVI和MOD16A2 ET/PET数据集，计算了干旱严重程度指数(DSI)并将其用于农业干旱监测，为DSI在中国乃至世界其他地区农业监测中的应用提供了理论基础；刘一哲等[13]采用模糊数学法建立了基于MODIS TVDI的干旱等级划分标准，实现对藏北地区春夏旱情的动态连续监测并分析旱情的时空变化特征。Khan等[14]利用MODIS NDVI与LST数据计算得到植被温度条件指数(VTCI)并采用地理空间近实时耦合(NRTC)方法对巴基斯坦Punjab平原干旱进行了研究。杨波等[15]利用MODIS增强植被状态指数(EVCI)和温度状态指数(TCI)建立了干旱状态指数(DCI)，将该指数应用于湖南省农作物的旱情监测并得到了旱情等级的空间分布图。

NDVI在进行干旱监测时得到了较为广泛的应用，但NDVI在监测干旱时不能及时反映土壤水分含量，只有当水分胁迫十分严重进而阻碍了作物生长时才会引起NDVI值的显著变化，这表明NDVI对重旱有较好的反映[16]，因此NDVI对土壤水分的反映以及应用于干旱监测具有一定的滞后性[17-18]。有研究表明，在高纬度、高海拔地区，尤其是作物生长的前期和后期，基于NDVI-LST特征空间的TVDI模型可能并不适用[19]。归一化水体指数(NDWI)一般用于识别水体目标，NDWI对植被冠层含水率比NDVI更为敏感，在短期干旱监测中，NDWI能及时地反映旱情的时空变化[20]。

河南省是我国的农业大省，粮食年产量占全国总产量的10%左右，是我国粮食重要生产地区[21]。近年来，在全球变暖的背景下，该地区的干旱灾害日益严重，对粮食生产构成巨大威胁，干旱问题研究亟待加强[22]。准确监测干旱的发生时间、发展程度和影响范围,对保障社会经济发展、促进生态环境恢复、维持区域和谐稳定具有重要意义。2014年河南省大部分地区遭受百年不遇的大旱，因此本文选取2014年MODIS观测数据从时间和空间尺度上分别分析NDWI与实测土壤相对湿度之间的相关性，根据相关系数来评价NDWI用于河南省干旱监测的适用性。

1 研究区概况

河南省位于我国中东部、黄河下游，地处北纬31°23′～36°22′，东经110°21′～116°39′之间，总面积达16.7万km2。河南省地势西高东低，由平原、盆地、山地、丘陵等构成，具有独特的地理位置和复杂的地貌特征。大部分地区气候处于暖温带，南部跨亚热带，属于大陆性季风气候。在北亚热带向暖温带气候过渡、自东向西由平原向丘陵山地气候过渡过程中，受季风型气候的影响，降雨分布不均匀[23]，具有四季分明、雨热同期、复杂多样和旱涝灾害频繁的特点。全省由北向南平均气温为12.1～15.7℃，年均降水量532.5～1 380.6 mm，降雨以6—8月份为主[24]。河南省主要土地利用类型以耕地为主。

2 数据与方法

2.1 数据来源与预处理

2.1.1 遥感数据来源及预处理本文使用的MODIS数据来源于美国国家航天局NASA网站LAADS DAAC数据中心(https://ladsweb.modaps.eosdis.nasa.gov/)，时间为2014年，编号为H27V05。MOD09A1、MYD09A1反射率产品提供了波段1-7的500 m分辨率8 d合成的反射率数据，MOD13A1、MYD13A1数据为16 d合成的500 m空间分辨率植被指数产品。其中，MOD为Terra卫星产品，数据获取时间约为地方时上午10∶30，MYD为Aqua卫星产品，数据获取时间约为地方时13∶30。由于同一天Terra和Aqua两颗卫星过境时间不一样，可在MOD产品出现缺失值时利用MYD产品对其进行填补生成较为连续的数据集。以上MODIS数据利用MRT处理工具对其进行波段提取、投影及格式转换。

Landsat8 OLI数据来源于美国地质调查局(https://earthexplorer.usgs.gov/)，其可见光波段数据分辨率为30 m。本研究选取2013—2015年植被生长季(6—9月)，条代号为(122-126，35-38)共12景数据，生成河南省30 m分辨率的NDVI分布图，用于分析观测站点周围地表的异质性。

2.1.2 测站土壤水分数据土壤相对湿度数据来源于河南省气象局，监测深度为10、20、30、40 cm和50 cm，共有162个土壤墒情监测站数据。选取2014年每天的实测土壤相对湿度数据，剔除缺测土壤层的数据点和有明显异常的数据点剩余148个站点(图1)，根据MODIS产品时间尺度将其以8 d进行平均，用于验证NDWI的干旱监测精度。

图1 河南省气象站分布图

2.2 植被指数计算及预处理

2.2.1 归一化水体指数归一化水体指数(normalized difference water index,NDWI)是基于近红外波段与绿波段建立的归一化比值指数[25]，在利用遥感影像提取水体方面应用十分广泛。计算公式如下：

(1)

式中，ρNIR、ρSWIR分别为近红外波段、短波红外波段的反射率。本研究利用MODIS第2波段和第5波段数据计算NDWI。

2.2.2 增强型植被指数增强型植被指数(enhanced vegetation index,EVI)是目前应用比较广泛的植被指数，它利用背景调节参数L和大气修正参数C1、C2同时减少背景和大气的作用，对气溶胶等残留做了进一步改正，并具有消除土壤背景和大气影响的优势。其计算公式如下：

(2)

式中，ρBLUE和ρRED分别为蓝光和红光波段反射率，L为土壤调节参数，C1和C2为拟合系数。MODIS EVI产品数据反演过程中，L=1，C1=6，C2=7.5。EVI时间序列季节性明显，能够更好地反映高植被覆盖区的季节性变化特征，并且很少有突降现象，时间序列曲线较平滑[26]。

2.3 研究方法

2.3.1 图像信息熵信息熵常被用来作为一个系统信息含量的量化指标，从而可以进一步用来作为系统方程优化的目标或者参数选择的判据。其计算公式如下：

(3)

式中，H(X)是求得的信息熵的值；P(x)表示x在计算范围内出现的概率。由信息熵的定义可知，一个系统越有序，信息熵就越低；反之，一个系统越混乱，信息熵就越高[27]。对于图像而言，系统有序表示灰度值相同或相近，反之表示灰度值相差较大。对NDVI图像灰度图而言，若灰度值相同或相近，可表明地物相同或相似，灰度值不同则表明地物类型不同。

利用土壤水分实测数据对基于MODIS数据计算的NDWI用于农业干旱监测适用性进行验证时，应选择像元范围内下垫面较为单一的站点数据。由于建筑物与植被的反射率不同，通过30m分辨率NDVI影像很容易区分。本文利用NDVI图像信息熵方法对气象站点进行筛选，以气象站所在像元为中心计算33*33窗口内图像的信息熵值，根据熵值大小判断站点周围地表异质性。

2.3.2 时间序列谐波分析法时间序列谐波分析法(harmonic analysis of time series，HANTS)用于对时间序列数据进行平滑，它能够充分利用遥感影像的时间性和空间性，将空间上的分布规律和时间上的变化规律联系起来。本文利用ENVI提供的HANTS工具对EVI进行了时间序列平滑，该方法通过最小二乘法的迭代拟合去除时序EVI值中受云污染影响较大的点，借助于傅立叶在时间域和频率域的正反变换实现曲线的分解和重构可生成较为平滑的EVI数据集。HANTS方法充分考虑了植被生长周期性和数据本身的双重特点，可真实反映植被的周期性变化规律，平滑前后如图2所示。

图2 EVI平滑前后对比图

2.3.3 相关性分析本文通过计算气象站点的土壤湿度数据和NDWI的相关系数，评价NDWI在河南省的干旱监测能力。相关系数的取值范围在[-1,1]之间，其绝对值越大表明相关程度越高，正(负)号代表两个变量之间呈正(负)相关性。干旱指数与土壤湿度之间的相关程度，用于表征某种干旱指数是否能够准确反映土壤的干湿状况[28]。

3 结果与分析

3.1 观测站点异质性评价

本文采用NDVI图像信息熵方法分析了观测站点1 km2范围内的地表异质性，利用Google Earth查看气象站周围环境并与信息熵计算结果进行对比，结果发现NDVI图像信息熵的大小与站点周围的均匀程度具有较好的一致性，表1为辉县等6个观测站点的NDVI图像信息熵，各站点对应周边如图3所示。

图3 不同站点周围环境分布

表1 不同站点的信息熵值

3.2 NDWI与土壤水、植被覆盖关系分析

在干旱发生时，大气相对湿度、气温与植被受水分胁迫相关。气温越高，大气相对湿度越小，则植被受水分胁迫也越严重。图4为2014年开封市大后庄站点土壤水分、NDWI与EVI的变化趋势图。此站点位于开封市南部，地处中原腹地，地势平坦、土壤肥沃，适宜各类农作物种植，属温带季风气候，四季分明。由图4可知，NDWI与土壤相对湿度的变化呈现明显的负相关性。EVI随时间呈现波动性变化趋势，首次波峰出现在第73天，波谷出现在第129天，土壤水分在此区间变化状态总体上呈增长状态，波动情况比较明显，NDWI的增长趋势则较为平缓。EVI在第193、257天达到最大值和最小值，分别为0.70、0.09。由于受降雨的影响，土壤水分的波动变化较大，但NDWI对其负相关性依旧不变。总体来说，NDWI与土壤水分呈负相关性,在EVI值较大的情况下负相关性更加明显。

图4 大后庄站点NDWI、SM10、EVI年变化情况

3.3 时间序列相关性分析

为验证NDWI的干旱监测能力，选取地表异质性较小的站点，将EVI的阈值设为0.4并利用NDWI与土壤水分数据进行相关性分析(图5)。大后庄、后李、留固、正阳、董天龙气象站点均匀分布在河南省平原区域，受地势起伏影响较小且植被长势较好，NDWI与土壤相对湿度有较好的相关性，决定系数R2均能达到0.4以上。在EVI>0.4时，大后庄、后李和留固的NDWI与土壤相对湿度的相关性均大于EVI<0.4的情况，正阳和董天龙站在两种状态下的相关性基本相同。社旗站点位于河南省西南部南阳盆地东部，地形平缓，植被覆盖度相对较高，NDWI指数对土壤水分的反映较好。栾川、走马岭、梁庄等站点位于河南省西部，地形以山脉、丘陵为主，土壤相对湿度与NDWI的决定系数R2分别为0.253、0.175、0.084，说明地形的复杂多样以及小范围气候带的差异性往往会对遥感监测带来一定的影响，海拔较高地区对光谱反射率的准确性造成误差，从而影响NDWI指数对土壤水分的监测。整体来看，NDWI与土壤水分呈负相关性且在平原区域相关性较好，在植被覆盖高(EVI>0.4)的区域NDWI对土壤水分的反演更为敏感。

图5 NDWI与土壤相对湿度时间序列相关性

3.4 空间尺度相关性分析

河南省各地地形、气候、土壤类型均有所区别，豫南山地以及整个西部区域地形复杂、海拔较高、气候多样，再加上山地丘陵区与平原地区土壤持水能力的不同，对遥感数据所计算的NDWI指数有一定的影响，NDWI指数对土壤水分的相关性可能有一定的差异。为了减小气候、土壤类型等因素对实验结果的影响，根据地形特征将河南省分为北部、中部、南部、西部，安阳、濮阳、鹤壁、新乡、焦作为北部地区，郑州、开封、商丘、许昌、漯河、周口为中部地区，驻马店和信阳为南部区域，西部地区有三门峡、洛阳、南阳和平顶山。随机选取第121、201、313天数据计算不同区域NDWI与实测土壤相对湿度之间的相关性以分析NDWI在河南省干旱监测中的适用性。

图6展示了第121、201、313天空间上4个区域的NDWI指数与土壤水分的相关性。由图6可知，第121天时NDWI与土壤水分的相关性由高到低依次是北部、中部、南部、西部，决定系数R2依次为0.2998、0.2839、0.1372、0.0811。植被高覆盖区域主要集中在北部地区鹤壁、新乡和中部地区商丘、周口、漯河以及南部地区驻马店，这3个区域地势均相对较为平坦，NDWI与土壤水分的相关性较高。南部地区水系较其他区域发达，植被覆盖度高于其他区域，因此NDWI与土壤相对湿度的相关性也较高。而西部地区地形较为复杂，NDWI与土壤水分相关性低于其他区域，说明NDWI更适用于平原地区的干旱监测。河南省植被的生长季一般位于3—9月份，第201天正是植被快速生长阶段，北部和中部地区NDWI与土壤相对湿度之间仍然呈负相关，决定系数分别为0.206和0.255。然而在南部和西部地区NDWI与土壤水分呈现正相关性，与理论相矛盾，可见由于地形和其他未知因素影响了两者之间的相关性，NDWI在山地丘陵地进行干旱监测适用性较差。北部、中部、南部和西部地区第313天土壤水分与NDWI的决定系数依次是0.165、0.163、0.482、0.139，该时间4个区域NDWI与土壤水分的相关关系仍然为负相关，北部和中部地区相关性仍然高于西部地区，但南部地区相关性最高。由于受天气影响，该时间内南部地区无云覆盖的观测站点较少，该相关性不能完全代表整个区域。再一个可能的原因是南部地区第313天内温度和植被覆盖要高于其他地区，因而该区域内NDWI与土壤水分的相关性高于其他区域。

图6 NDWI与土壤相对湿度空间相关性

3.5 基于NDWI干旱情况分析

2014年为河南省特大干旱年，各月份NDWI的空间分布及变化趋势如图7所示。总体来说，在时间和空间上均体现出了全年干旱变化情况，4—8月份全省较为干旱，在8月26—30日期间，河南省大部分地区有较大降雨后，对比图7中8、9月份两个月的NDWI分布图可以看出，9月份的NDWI值明显降低，说明地表土壤水分含量增加，全省旱情得到缓解。

图7 河南省2014年不同时期NDWI空间分布图

在豫北山地以东、豫东平原、豫南山地以北区域NDWI值均呈现较低趋势，表明土壤含水量较大，地区较为湿润。11月、12月及翌年1月为冬季，土壤较为湿润的地区主要集中在豫北山地，即太行山脉附近、伏牛山、熊耳山、小秦岭、豫南山地等山地地区。此时其他地区庄稼均已收割，冬季作物刚播种下去，地表植被覆盖度较低，而山区植被类型大多为林地，此时植被覆盖度要高于其他平原地区，在相同的降雨条件下，山地林区比平原地区具有更强的保水能力；同时，低海拔地区地表蒸发和植物蒸腾作用强烈，对土壤水的消耗作用较大，而高海拔地区温度较低，对地表蒸散发具有一定的限制作用，土壤水分的消耗相对较少[29]。因此在11月、12月及翌年1月山地的NDWI明显低于其他平原地区。

4 讨论

利用地面站点观测数据验证卫星遥感数据时，两者尺度不同，直接验证会存在误差，地表异质性越强，误差越大。本文提出利用NDVI图像信息熵值对观测站点进行筛选，通过和Google Earth图像对比，该方法效果较好，可用于其他卫星遥感产品检验时站点异质性分析，其操作比半变异函数更为简单。

NDWI可有效指示植被冠层的含水量信息，在植被受水分胁迫时，指示效果比NDVI敏感，可用于进行大范围的干旱监测，选用归一化水体指数NDWI法，还可最大限度地消除植被和土壤等信息从而提高精度[30]。总体上，NDWI与冠层含水量呈负相关，且植被覆盖度越高相关性越好。本文利用地面观测土壤相对湿度验证NDWI用于干旱监测的精度，土壤相对湿度较低时，植被冠层开始受水分胁迫，但胁迫过程存在一定的延时性，而NDWI指示的是植被冠层的含水量信息，利用观测的土壤相对湿度验证时会产生一定的影响。在地形较为复杂的地区，卫星遥感观测精度会有所降低，NDWI用于干旱监测的精度比平原地区有所降低，应对遥感数据进行地形校正后再使用。河南省大部分地区为平原，西部三门峡及洛阳部分地区为山地，其余信阳、新乡、焦作等地有少量山地，利用NDWI进行干旱监测具有可行性。

2014年为河南省特大干旱年份，NDWI较好的反映出了该年份干旱的时空动态变化规律。本文研究结果可为不同地表大范围遥感干旱监测提供了参考，且更易于实施。