多种植被指数与土壤墒情响应关系对比分析

2020-06-11王金鑫于百顺李聪玲

节水灌溉 2020年3期

王金鑫，于百顺，李聪玲，姚静

(郑州大学水利科学与工程学院，河南郑州 450001)

0 引言

土壤墒情一般指地表土壤的湿度状况。是水文学、气象学以及农业科学研究领域中的一个重要指标参数[1]，是土壤-植物-大气连续体的一个重要因子[2]，对气候、农业、旱情监测都具有极为重要的意义[3]。传统台站式离散监测方法难以表达大区域连续的土壤墒情，已逐渐被能够快速进行周期性、大面积同步观测的遥感方法以及基于传感器的物联网自动监测方法所取代。早在1965年，西方地理学家Bowers与Hanks[4]就指出，将土壤光谱反射特性和热红外数据信息相结合的方法在探究土壤信息方面是可行的；1971年，Waston与 Rowen[5]使用地表温度日较差推导了热惯量公式，创立了热惯量墒情监测方法，成为日后该领域研究的理论基础；Kogan[6]在归一化植被指数(NDVI)的基础上提出了植被状态指数(VCI)，减弱了由于地理位置和土壤条件等因素导致的差异；Sandholt[7]等人运用陆地表面温度(LST)与归一化植被指数(NDVI)，建立了NDVI-Ts 特征空间，并在此基础上提出了温度植被干旱指数(TVDI)，用来估算土壤的表层的水分分布情况；Ulaby等人[8]最早开始研究主动微波遥感，通过实验得出，微波后向散射系数相对于土壤的介电常数有很大的相关性。我国开展土壤水分遥感监测的研究起步较晚，张仁华[9]提出一个改进的热惯量模型，该模型综合概括了地表显热通量和潜热能量等影响因素，并利用改进模型来反演土壤含水量；隋洪智[10]等人在直接利用卫星资料基础上，依据能量平衡方程的原理，进行简化推算出表面热惯量，通过经验模型拟合热惯量与土壤水分之间的相互关系对旱灾进行探测分析；姚春生[11]通过对各种遥感定量反演土壤含水量的模型进行了综合比较分析，提出对于裸地或作物的生长初期应该运用表观热惯量模型，对于稍高植被覆盖度选择TVDI或Ts/NDVI模型，以此来反演全国土壤含水量的分布情况；胡蝶等人[12]以甘肃定西为研究区，用Radarsat-2数据和土壤水分分析仪，较精确地解析了共极化反射系数与不同深度土壤水分之间的关系，结果表明共极化散射系数和 10～20 cm 土壤的湿度之间具有最高的相关性。墒情遥感监测从原理上可分为土壤热惯量法、土壤波谱特征法、能量平衡模型法和植被光谱特征法等四类方法[13]，它们各有其适用条件。基于植被指数(即植被光谱特征)的土壤墒情监测方法利用光学和红外遥感实现，其在时空间分辨率方面有较大优势，数据类型多、数量足，且方便实现大区域土壤墒情的遥感监测，一般适用于较高植被覆盖期[14]。植被指数直接反映植被的生长状况，影响植被生长状况的因子十分复杂，温度、湿度和肥力无疑是最主要的影响因子。在一定区域、一定时间内，气温、土壤、地形等因素相对稳定，只降雨量变化会对作物生长造成较显著影响[15]，因而，植被指数的变化与土壤墒情之间有较显著相关性[16]。虽然植被的缺水胁迫相对于其生长状态有一定的时滞性[17](研究表明，NDVI滞后于墒情波动约4～6 d[16])，但在农作物的整个生育期内，墒情的波动与植被指数的变化频率与趋势应该是一致的[16]。由于每种植被指数对水分的响应存在差异，而且植物在不同的季节、不同的生长阶段，具有不同的水分利用特征[18]，所以植被指数作为水分响应指标，本身存在适用性及其适宜的时空尺度[19]。因而，在利用遥感数据进行长时间、大范围农业干旱遥感监测过程中，如何根据不同区域、不同作物生长时期，综合考虑作物与其生长环境之间的关系，选择最合适的监测指标[19]就成为问题的关键。本文以河南省中东部黄淮海平原冬小麦主产区的农田墒情为研究对象，对归一化植被指数(NDVI) 、增强型植被指数(EVI) 、植被供水指数(VSWI) 、温度植被干旱指数(TVDI)这几种典型的植被指数与土壤含水量进行相关性分析，对比不同深度、不同植被指数对土壤含水量的响应，并且分别选用这四种植被指数与对应的土壤墒情数据进行回归分析，得到最优的响应深度及响应指数，为土壤墒情的监测提供指导。

1 研究区概况及数据源

1.1 研究区概况

河南省黄淮海平原(图1)位于北纬32°08′～36°21′和东经112°51′～116°35′之间，地势西高东低。西起海拔100 m等高线和豫西北丘陵边缘，东至豫鲁、豫皖分界线，北起豫冀交界，南至淮河一线，面积约8.7 万km2，占河南省总面积的52.4%。自北向南横跨海河、黄河、淮河三大流域。地处暖温带和北亚热带，属于湿润至半湿润季风气候。多年平均气温在12～16 ℃之间，气温年较差、日较差均较大，全年无霜期从北向南约为180～240 d。年平均降水量在500～900 mm，受到季风影响，年降水量时空分布不均，由北向南逐渐递增，全年降水量主要集中在夏季，大约占全年降水量的一半左右，冬季降水较少。日照时数2 100～2 600 h，光热资源较充足，基本可以满足一年两熟的农作物需求。年均蒸发量1 300～1 600 mm，由北向南逐渐递减[20]。该区自然条件优越，土壤类型众多，农业发达，是我国主要的冬小麦生产区。

1.2 数据源及数据预处理

本研究所涉及的遥感数据包括：研究区域2015年2-6月的MODIS红光波段、近红光波段和蓝光波段的地表反射率数据，MODIS的每日地表温度数据(数据来源： https:∥ladsweb.nascom.nasa.gov/)。其中，从MOD09GQK/全球250 m地表反射率产品中得到红光波段、近红光波段的反射率，分辨率为250 m；从MOD09GHK/全球500 m地表反射率产品中得到蓝光波段的反射率，分辨率为500 m；从MOD11_L2/1KM地表温度产品中得到每日地表温度，分辨率为1 000 m，该地表温度产品分为白天和晚上两幅影像，这里取其平均值作为当天的地表温度。以上数据都已经经过几何校正和大气校正。利用MRT软件统一转换为WGS84基准的Albers等积投影，并将地表温度数据和蓝色波段地表反射率数据重采样为250 m的分辨率。

本研究所涉及的其他数据包括：研究区域2015年2-6月的0～10 cm和10～20 cm土壤水分含量实测数据(进而可得到0～20 cm实测数据)、2014-2015年河南省冬小麦种植区域数据以及河南省行政界线数据，从河南省气象局调研获取。土壤含水量实测数据是研究区102个站点的0～10 cm和10～20 cm每日土壤体积含水量的平均值，所有站点均位于小麦大田内。其中豫北区域29个站点，豫中区域57个站点，豫南区域15个站点，分布如图1所示。每个站点通过经纬度坐标与MODIS图像像元进行匹配。将2014-2015年冬小麦种植区域数据、河南省界线数据也转换为与上述遥感数据统一的基准与投影，然后对影像数据进行剪裁。

图1 研究区站点分布图Fig.1 Site map of the study area

2 研究方法与过程

2.1 冬小麦生长区域及生长期的划分

本文以33°N纬线(驻马店市与漯河市交界处)和黄河为界，将河南省主要种植冬小麦的区域分为豫北，豫中，豫南3个地区[16]。

中原地区冬小麦生长的茂盛期为拔节到灌浆之间，不同地区略有差异，如表1[16]所示。

表1 研究区域冬小麦生长茂盛期的划分Tab.1 Division of the flourishing period of winter wheat in the study area

2.2 土壤墒情监测指标

归一化植被指数(Normalized Difference Vegetation Index，NDVI)：

NDVI=(NIR-RED)/(NIR+RED)

(1)

式中：NIR表示近红外波段的反射值；RED表示红光波段的反射值。NDVI越大，植被长势越好；但其在较高植被覆盖度和较高产量水平下，易出现饱和现象[21-24]。

地表温度(Land Surface Temperature，LST)是反映土壤墒情的重要参数。裸地区域，上层土壤温度代表了地表温度。其含水量越低，土壤热容量越小，在接收太阳辐射后表面温度上升得越高；在植被覆盖区域，遥感监测到的地表温度实际上是植被的冠层温度[25]。当土壤含水量下降时，作物为了减少由于蒸腾作用所导致的水分损失而闭合叶片的部分气孔，引起地表潜热通量的降低和感热通量的增加，从而导致冠层温度的升高[26]。

增强型植被指数(Enhanced Vegetation Index，EVI)是对NDVI饱和性的一种改进。它引入了蓝光波段和土壤调节因子，有效地消除了大气和植被的影响，有研究表明EVI和NDVI存在很大的相关性[27]，其计算公式如下：

EVI=2.5(NIR-RED)/(NIR+6.0RED-7.5BLUE+1)

(2)

式中：NIR、RED、BLUE分别代表近红外波段、红光波段、蓝光波段的反射值。

由于植被指数对作物水分的反映具有一定的延迟，而冠层温度又容易受大气环境状况的影响，单一考虑植被指数或冠层温度来反映作物水分不够敏感和有效[28]。植被供水指数(Vegetation Supply Water Index，VSWI)是综合归一化植被指数和地表植被冠层温度的一种用于监测土壤缺水状况的综合指数。其定义如下：

VSWI=NDVI/Ts

(3)

式中：VSWI为植被供水指数；NDVI为归一化植被指数；Ts为植被冠层温度。

当植被供水正常时，归一化植被指数在一定的生长期内稳定在一个范围，卫星遥感的植被冠层温度也保持在一个稳定的范围；若土壤含水量下降，植被供水不足，一方面植被的生长受到影响，归一化植被指数降低，另一方面植被为避免过多失水而关闭部分的叶片的气孔，导致植被冠层温度升高。因其利用的是作物缺水时在不同反射波段上的反应，所以在植被覆盖度比较高的地区尤其是在作物生长茂盛期更为适用[29]。由于VSWI数值很小，为方便计算比较在实验中将VSWI放大100倍。

Price[30]、Carlson等[31]研究发现利用遥感反演地表温度与植被指数所构成的散点分布图呈三角形，而Nemani等[32]、Moran等[33]研究认为LST和NDVI构成的散点分布图为梯形，即LST-NDVI特征空间。Sandholt等[7]在基于LST-NDVI特征空间的基础上提出了温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)。计算公式如下：

TVDI=(Ts-Tsmin)/(Tsmax-Tsmin)

(4)

式中：Ts表示待测像元的地表温度；Tsmax和Tsmin表示待测像元NDVI在特征空间中对应的最大地表温度(干边温度)和最小地表温度(湿边温度)，其计算公式如下：

Tsmin=a1+b1NDVI

(5)

Tsmax=a2+b2NDVI

(6)

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

实际应用中，应根据研究区植被覆盖状况和遥感数据源来选择LST-NDVI特征空间的形状[34,35]。在植被覆盖度较高或者使用卫星遥感数据时，采用梯形特征空间较好[34]。在冬小麦主要生长期内，TVDI与土壤水分具有显著的相关性[36]。

2.3 研究过程

考虑到地域和物候期的差异，在选择数据样本时，本文将豫北、豫中、豫南3个地区按照植被生长状态相近的原则进行匹配，并考虑数据质量因素，每隔10 d左右的时间选择一组数据，共选择8组数据进行研究。表2为各组数据对应的日期(均为2015年)。

确立好八组样本数据后，首先对每个站点各相关植被指数进行计算提取；接着进行各植被指数与各层土壤墒情进行相关分析，得出最佳响应指数；最后建立各植被指数与各层土壤墒情的回归方程，实现墒情的定量反演，进一步验证实验结果。图2为研究区LST-NDVI特征空间实例。

3 结果与分析

3.1 植被指数与土壤墒情的相关性分析

表3所示的是土壤墒情与植被指数的相关情况。可以看出：

表2 各组数据对应日期Tab.2 Date corresponding to each group of data

图2 豫中地区第60天LST-NDVI特征空间Fig.2 LST-NDVI feature space on the 60th day of center Henan area

(1)各指数对不同深度土壤墒情在置信度(双测)为 0.01 时显著相关；

(2)NDVI、EVI、VSWI与土壤墒情呈现正相关，TVDI与土壤墒情呈现负相关；

(3)横向上，相关性呈现出NDVI

(4)纵向上，相关性呈现出10～20 cm>0～20 cm>0～10 cm的规律。NDVI与10～20 cm土壤墒情的相关性最佳的组数为6组，NDVI与0～20 cm土壤墒情相关性大于0～10 cm土壤墒情的组数为7组，总体呈现出10～20 cm>0～20 cm>0～10 cm的趋势；EVI与10～20 cm土壤墒情的相关性最佳的组数为5组，EVI与0～20 cm土壤墒情相关性大于0～10 cm土壤墒情的组数为7组，总体呈现出10～20 cm>0～20 cm>0～10 cm的趋势；VSWI与10～20 cm土壤墒情的相关性最佳的组数为7组，VSWI与0～20 cm土壤墒情相关性大于0～10 cm土壤墒情的组数为6组，总体呈现出10～20 cm>0～20 cm>0～10 cm的趋势；TVDI与10～20 cm土壤墒情的相关性最佳组数为4组，TVDI与0～20 cm土壤墒情相关性大于0～10 cm土壤墒情的组数为4组，不同深度土壤墒情与TVDI的相关性较均匀。可以看出，4种植被指数均对10～20 cm深度土壤墒情的反应最敏感。

表3 土壤墒情与植被指数Pearson相关性Tab.3 Pearson correlation between soil moisture and vegetation index

注：** 表示在置信度(双测)为 0.01 时，相关性是显著的。

图3为不同深度土壤墒情与各植被指数相关系数的八组平均趋势(为方便比较，对TVDI的相关系数取绝对值)。

图3 不同深度土壤墒情与各植被指数相关性比较Fig.3 Correlation between soil moisture and vegetation index at different depths

3.2 土壤墒情对植被指数的回归分析

在相关分析的基础上，对不同深度的土壤墒情与各植被指数进行回归分析。经过选择的站点按照4∶1的比例进行关系建立与验证。选择两组数据进行结果展示。

表4和表5分别为是第4组和第6组样本不同深度的土壤墒情和各植被指数的回归方程情况。

表5 第6组样本土壤墒情与植被指数的回归方程Tab.5 Regression equations for soil moisture and vegetation index of group 6

图4为第六组样本10～20 cm深度土壤墒情与各植被指数的回归情况，可代表本次实验不同深度土壤墒情与各植被指数的回归趋势。

图4 第6组样本10～20 cm土壤墒情与各植被指数的回归曲线Fig.4 Regression curve of 10～20 cm soil moisture and vegetation index of group 6

图5展示了这两组数据不同深度土壤墒情与各植被指数回归R2的变化趋势。

图5 第4组和第6组土壤墒情与植被指数回归R2变化趋势Fig.5 Regression trend of soil moisture and vegetation index in Group 4 and Group 6

从上面图表中，可以得出：各层土壤墒情与各植被指数的回归结果与前面相关性分析结果一致。横向上，回归方程的R2都基本呈现出NDVI0～20 cm>0～10 cm的规律。

3.3 预测精度分析

表6为利用植被指数预测土壤墒情的精度情况统计。

表6 第4组与第6组样本各指数预测墒情与实际墒情差的标准偏差Tab.6 Standard deviations of the indices of the 4th and 6th samples for predicting and actual soil moisture

从表6中可以看出：各组土壤墒情定量预测的标准偏差横向上符合NDVI>EVI>VSWI>TVDI，纵向上符合10～20 cm<0～20 cm<0～10 cm的规律。与前面的实验结一致。

4 结论与讨论

4.1 结论

本文以中原地区大区域冬小麦土壤墒情的遥感监测为研究对象，利用共享遥感数据和实测墒情数据，选择4种典型的植被指数，考虑地域和农作物物候期差异，探究不同深度土壤墒情与各植被指数的响应关系。主要结论如下。

(1)从相关系数的角度，各植被指数与不同深度土壤含水量的相关系数呈现出NDVI0～20 cm>0～10 cm的规律性。其中，NDVI、EVI、VSWI为正相关，TVDI为负相关(这里TVDI的相关系数取其绝对值)；

(2)从回归显著性检验的角度，各植被指数与不同深度土壤含水量的回归R2呈现出NDVI0～20 cm>0～10 cm的规律性；

(3)从土壤墒情定量反演的角度，各植被指数与不同深度土壤含水量定量预测标准偏差呈现出NDVI>EVI>VSWI>TVDI、10～20 cm<0～20 cm<0～10 cm的规律性。

NDVI、EVI与VSWI直接反映植被的生长状况，与植被根部的墒情密切相关。墒情越好，植被生长就越茂盛。冬小麦的根部大致位于10～20 cm深处，所以其上述指数与10～20 cm的土壤墒情成较显著的正相关。相反，一定条件下，土壤含水量越高，其温度就越低，所以，主要反映土壤温度因子的TVDI与土壤墒情成负相关。

综上，针对中原地区的冬小麦而言，4种植被指数与不同深度土壤墒情均具有较显著的相关性。其中，TVDI是最佳的墒情响应指数，10～20 cm是植被指数最佳的墒情响应深度。利用TVDI进行遥感农田墒情反演是一种成本低、适应性广、精度较高的有效方法。本文的研究结果对于冬小麦田间管理、产量估测以及粮仓安全保障具有重要参考价值。