基于修正NDVI时间序列的大区域冬小麦全生育期墒情监测

王金鑫，赵光成，张广周，于百顺，罗蔚然，张成才，李 颖

(1.郑州大学水利与环境学院，郑州 450001；2. 中国气象局·河南省农业气象保障与应用技术重点实验室，郑州 450003；3. 河南省气象科学研究所，郑州 450003)

0 引 言

遥感墒情监测是目前遥感技术应用的前沿领域[1,2]，也是公认的世界性研究难题之一[1-4]。究其原因，影响墒情的因素包括土壤类型、降水、太阳高度、大气、植物长势、地表状况以及田间管理等等，使土壤墒情遥感反演异常困难[1,5,6]。基于植被指数的墒情监测或反演是遥感墒情监测的常用方法。在一定区域、一定时间内，气温、土壤、地形等因素相对稳定，只降雨量变化会对作物生长造成较显著影响[7]，因而，植被指数的变化与土壤墒情之间有较显著相关性。但植被受缺水胁迫在短期内不会在生长状态上明显表现出来，具有一定的时滞性[8]，无法及时地反映出土壤水分状况(后面有详细讨论)，所以植被指数是一个灵敏度较差的水分胁迫指标[9]，但在农作物的整个生育期内，墒情的波动与植被指数的变化频率与趋势应该是一致的。除早期的基于归一化植被指数(NDVI)的距平植被指数(AVI)、条件植被指数(VCI)[10]外，考虑地表温度影响(其是十分复杂的[11])，人们又提出了植被供水指数(VSWI)[12,13]、温度植被旱情指数(TVDI)[14]、条件植被温度指数(VTCI)[15,16]以及特征空间法[14]等新的概念和方法用于土壤墒情监测，相应于研究案例都取得了较好的效果。经过多年的研究积累，墒情遥感监测也达成一些基本共识：①从土壤深度角度，以10～20 cm左右效果最好[1-3,17]，表层和深层土壤水分与遥感资料的相关关系往往不佳[3]；②从地表覆盖角度，裸土和低植被以热惯量模型和微波遥感效果最好[3,4,18]；完全植被覆盖以作物缺水指数(CWSI)、距平植被指数、条件植被指数、特征空间模型、能量指数法和植被供水指数法等效果较好[1-3,19]；部分植被覆盖采用农田蒸散双层模型(但模型复杂)、改进求解的热惯量模型和温度调节指数效果较好[3,18]。概而言之，裸土时，土壤的热特性起决定作用，基于温度的方法较好；植被茂密时，植被的生长状态是重要的指示因子，指数法占优；半植被覆盖时，要同时考虑土壤和植被热特性差异的双重效应，综合指数或方法最好。充分证明了地表温度和植被指数是农业墒情光学遥感监测中的主要指示因子[4]。

国内外已有研究表明，任何单一的指数或方法都有自己的特点和适用条件，地表生态的复杂性和遥感技术的特点也决定了不存在一种普适的墒情监测指数和方法。大范围农田墒情遥感监测更是存在遥感数据的适用性选择以及光热条件的区域分异问题。本文拟提出一种基于高时间分辨率、中空间分辨率的遥感时序数据，考虑地域和物候期差异，利用距平植被指数，实现大区域冬小麦全生育期土壤相对湿度的监测方法。以河南省中东部黄淮海平原冬小麦主产区为研究对象，以MODIS NDVI为基础墒情指数，针对不同地区、不同生长阶段引进相应因子进行自适应改进，实现对农田墒情及其时空分布规律的监测。

1 研究区域与数据源

河南省黄淮海平原位于北纬32°08′～36°21′和东经112°51′～116°35′之间，北起豫冀交界，南至淮河一线，西起海拔100 m等高线和豫西北丘陵边缘，东至豫鲁、豫皖分界线(为方便数据处理，本文以该区域的行政界线为准)，面积约8.7 万km2，占黄淮海平原总面积的25.6%，河南省总面积的52.4%，属于南暖温带半湿润大陆性季风气候区[20]。该地区是河南省的主要农业区，冬小麦主产区，年农业总产值占全省农业总产值的70%以上。

本文采用的遥感数据包括：研究区域2011-2015小麦生产年度的b1、b2和b7波段的Terra MODIS数据、MODIS的每日陆面温度产品等，从共享网站(https:∥ladsweb.nascom.nasa.gov/)下载。 b1、b2和b7波段数据的数据等级是L2级的陆地表面反射率数据，已经过大气纠正和几何纠正。b1和b2波段的空间分辨率为250 m，共下载1367幅；b7波段的分辨率为500 m，根据需要下载了75幅；MODIS每日陆面温度数据(简称遥感地表温度数据，包括每天白天和夜晚两幅数据，这里取其平均值)是从共享网站下载的成品数据，空间分辨率是1 000 m，数量和使用情况同b1、b2波段。主要问题是有一些日期的数据“漏洞”很多(没有温度反演数据)。MODIS源数据为Sinusoidal projection，利用MRT软件统一转换为WGS84基准的Albers等积投影，并将地温数据和b7波段数据重采样为250 m，以与b1、b2波段匹配，进行像元级计算。

本研究所涉及的墒情及其他数据：研究区域2014- 2015年度10和20 cm土壤水分含量实测数据、2011-2014年2月上旬和11月下旬的10和20 cm土壤水分含量实测数据、豫南、中、北地区的多个监测站1981-2010年(30年)逐月平均气温数据、2014-2015年河南省冬小麦种植区域数据以及河南省行政界线数据等，从河南省气象局调研获取。墒情实测数据是研究区域均匀分布的100多个站点的10和20 cm每日土壤体积含水量的平均值。以各观测站30年的逐月平均气温进行内插(本文的内插方法均采用基于球面函数的普通克里金插值)，分别得到各区域每个像元的常年逐月平均气温。将2014- 2015年冬小麦种植区域数据、河南省界线数据也转换为与上述遥感数据统一的基准与投影，然后对影像数据进行剪裁和植被指数计算。由于冬小麦是河南省的主要粮食作物，种植范围年际变化不大，所以本文假设2014-2015年的小麦种植范围同前三年的种植范围相同，并用之对MODIS数据做掩膜，得到最终的研究区域。

2 研究方法与过程

2.1 研究思想阐述

中原地区的冬小麦生长周期较长，从头年十月到来年六月，跨深秋、冬季、春季和初夏几个季节，长达八个月，温度起伏大，生态因素多变。要进行大尺度麦田墒情监测，既要考虑地域差异，又要考虑冬小麦的物候特征。根据相关研究成果，黄河和驻(马店)漯(河)交界是河南省重要的综合地理分区和小麦生态分区界线[21,22]。本文以33°N纬线(驻马店市与漯河市交界处)和黄河为界，把河南省黄淮海平原冬小麦主产区划分为豫南、豫中、豫北三个地区(图1)；大致以返青和乳熟期为基点(咨询河南农业科学院和国家小麦工程中心专家)把小麦全生育期划分为：前期(从播种到返青，为裸土到半植被覆盖期)、中期(从返青到乳熟，全植被覆盖期)和后期(从乳熟到收割，为麦子变黄时期)，见表1。

表1 研究区冬小麦生长期划分Tab.1 Division of growing period of winter wheat in research

图1 河南省冬小麦主产区区域划分Fig.1 Regional division of main wheat producing areas in Henan Province

在小麦生长前期，麦苗尚不能完全覆盖裸土，NDVI不能完全反映墒情状况，必须考虑土壤温度因子。根据土壤的热特性，在同一(季节)时间(光热条件下)，若土壤含水量高，其温度要降低；若土壤含水量低，其温度必升高。考虑到NDVI为一个无量纲的量，这里对原始NDVI进行以下修正：

(1)

在小麦生长中期，麦子生长茂盛，基本全部覆盖地表，NDVI可以指示根部墒情。但正如前面所述，NDVI相对于墒情波动具有滞后性。图2是豫中龙城站中期NDVI与10和20cm土壤含水量的曲线图。NDVI与10、20cm土壤含水量的相关系数分别为-0.218，-0.213；将植被指数曲线前移5d，NDVI与10、20cm土壤含水量的相关系数分别为0.664，0.467。可见，NDVI滞后于墒情波动约4～6d。此外，众所周知，在高植被覆盖期，NDVI具有饱和性。需注意的是，这里的NDVI是利用遥感波段计算的原始值，由于受到太阳光照角度、观测角度、云、气溶胶等因素的影响，其往往包含着很多噪声，并不都是植被指数的真实值(一般比实际值要小)，尤其是低值部分。

图2 豫中龙城监测站NDVI相对于墒情变化的滞后性Fig.2 The lag of NDVI relative to soil moisture change in Longcheng monitoring station

基于上述考虑，对于中期NDVI修改如下：

(2)

式中：KMNDVI指中期修正NDVI(MediumNDVI)；β为修正加系数；Tn为某像元某天的实际地表温度；Tn-1为相应像元该天的前一天的实际地表温度；其他符号同前。

这里修改的前半部分(分子第一项)主要作用是将NDVI的滞后性提前。其理由是如果当天下雨或大面积灌溉，地表温度肯定低于前一天的温度，温差较大时，修改效果明显；如果不下雨或没有灌溉，那么二者的温度接近，修改效果可以忽略。后一部分(分子第二项)的作用与前期类似，这里主要用于改善NDVI的饱和性。

在小麦生长后期，麦子将逐渐变黄，NDVI的值要逐渐变小，但并不一定意味着土壤湿度低。因此也要对后期NDVI进行修正。杜晓等[23]依据对水的吸收率曲线和植被与土壤反射率曲线的分析，指出MODIS的b6、b7波段反映了含水土壤和植被的波谱信息，并构造了地表含水量指数(SWCI)模型。在此基础上，张红卫等[24]等考虑绿色植被叶绿素对光谱的吸收作用，引入b1波段信息，以b1/b7为修正系数，构造了增强型土壤表层水分含水量指数(ESWCI)模型，并应用于植被茂盛期的墒情监测。小麦生长成熟期是一个叶绿素含量逐渐减少，作物由茂盛过渡到黄枯的阶段，所以，与上述修正系数相反，本文对该阶段的植被指数修改如下：

(3)

式中：KLNDVI指后期修改NDVI(LateNDVI)；γ为修正乘系数；b1、b2、b7为MODIS相应波段光谱反射率。

将修正NDVI与实测墒情数据进行相关分析，验证其有效性。在通过验证的基础上，构建修改NDVI时间序列，并与前三年的平均值进行距平分析，得到监测年与前三年平均相比的墒情分布时空规律。

2.2 技术路线与过程

综上所述，本文的技术路线如图3所示。

图3 研究技术路线Fig.3 Research technical route

本文通过旬最大值合成(MVC)和改进S-G滤波法[25]，去除噪声数据，构建反映冬小麦生长实际规律的、较光滑连续的修正NDVI生长曲线。

3 结果与分析

3.1 修正前后像元生长曲线对比

图4表示豫中淮阳站2014-2015年植被指数修改前后的情况。

图4 豫中淮阳监测站2014-2015年NDVI修正前后的效果Fig.4 The effect comparisons of NDVI before and after modification of Huaiyang observation station at central Henan province in 2014-2015

可以发现，修正前后，植被指数曲线的走势基本是一样的，但修正后植被指数的变化幅度一般较小。

3.2 修正NDVI与实测墒情的相关分析

为了验证修正NDVI的有效性，分别在3个地区均匀选择8个监测站，将修正前后的植被指数与10和20cm实测墒情进行了相关分析。其结果如表2所示。

由表2可知，原始NDVI与实测墒情的相关情况，后期比前期和中期要好，中期比前期要好。修正后NDVI与墒情相关的规律与原始NDVI基本类似，但修正NDVI在各时期的相关系数都有不同程度提高，说明了修正系数的有效性，修正后的NDVI更适宜作为农田的墒情指数。

表2 NDVI修正前后与实测墒情的相关分析对比Tab.2 Comparison of correlation between NDVI before and after modified and measured soil moisture

续表2 NDVI修正前后与实测墒情的相关分析对比

注：*表示在0.05水平(双侧)上显著相关，**表示在0.01水平(双侧)上显著相关。

3.3 经MVC和S-G滤波后像元生长曲线情况

图5是豫北某像元修正NDVI曲线情况。

图5 豫北某像元经MVG合成和改进S-G滤波后的修正NDVI生长曲线Fig.5 The growth curve of a Pixel modified NDVI after MVC synthesis and improved S-G filtering in the north of Henan province

3.4 基于修正NDVI距平分析的冬小麦墒情监测

基于生长曲线得到距平值系列(图6)，并进行距平值等级划分，最终得到监测年(2014-2015年)与前三年均值相比，研究区域冬小麦田相对湿度的时空分布规律(图7)。

图6 研究区像元的距平值时间序列Fig. 6 The time series of the anomalies of the study area

图7 2014-2015年河南省冬小麦主产区土壤湿度相对于前三年均值偏离情况的时空分布Fig.7 Temporal and spatial distribution of soil moisture relative to the mean of the previous three years in main winter wheat production areas in Henan province from 2014 to 2015

这里的关键是距平值(距平植被指数)的等级划分。实验数据(图6)表明，距平值的分布区间约为[-1,1]，集中在[-0.2,0.2]之内。如果等级划分过粗，那么就看不出墒情的空间分布结构；如果划分过细，则会引起数据失真，不能反映真实情况。综合相关研究文献[18,26]和实验结果，本文把距平值划分为7个等级：<-0.5，[-0.5,-0.15), [-0.15,-0.07)，[-0.07,0.07)，[0.07,0.15)，[0.15,0.5)，>0.5，分别表示相对湿度从干到湿的分布序列。其中，[-0.07,0.07)为与往年持平。需要说明的是，图7所示的情况是监测年相对于前三年平均的湿度偏离情况，并不是绝对的墒情等级。

从图7中可以看出，2014-2015年的10月下旬、3月中旬、4月中旬至5月上旬，研究区冬小麦田的墒情基本与前三年平均持平；3月下旬、4月上旬和5月中旬至下旬稍偏干；11月中旬至12月上旬偏干；12月下旬至3月上旬偏湿。

为了验证上述结果的正确性，我们选取了明显偏干的2014年11月下旬和明显偏湿的2015年2月上旬作为实证。将上述两个时期前三年的10和20 cm实测墒情(取最大值)的平均值与2015年相同时期实测墒情(取最大值)的值做比较(均值为被减数)，其结果如图8所示(等级划分的相对比例与图7一致)。从图8(a)、8(c)与8(b)，8(d)、8(f)与8(e)的比较可以看出，监测值与实测值分布规律基本是一致的，误差基本都在一个量级内；比较而言，监测值与10 cm实测土壤水分符合度更好。需要说明的是，图8的色标适用于除8(b)、8(e)以外的其他图。

图8 相对湿度验证实例Fig.8 Comparison of relative moisture verification

4 结 论

土壤墒情是作物旱情诊断和灌溉管理的基础参数，是农作物产量的限制因子，监测土壤墒情对提高水分利用效率和效益具有重要意义[27]。本文充分利用MODIS遥感的高时间分辨率特性，构建基于区域和作物生育期自适应修改的植被指数生长曲线，通过距平植被指数，实现大区域全生育期的冬小麦墒情监测。研究表明：

(1)与NDVI值相比，基于地域和物候期分异的修正NDVI值与土壤水分的相关性更好，更适合作为农田墒情的标识指数；

(2)基于自适应修正NDVI作物生长曲线的遥感墒情监测是一种成本低、精度较好、适宜性较广的有效方法；

(3)农作物生长曲线是一条隐含诸多农情信息的特征曲线，对其进行农情信息挖掘是一个很有价值的研究方向。

本文提出对中期NDVI修正的初衷之一是将其滞后性提前，但实验中我们发现，恰是在降雨的时候，往往也是云层遮挡较为严重的时候，MODIS的每日陆面温度就存在较为严重的漏洞。为保证结果的科学性，这时我们采取的措施是不修改。从图4可以看出，滞后性修改的效果不明显，这在一定程度上影响了实验结果。此外，距平植被指数法虽然不能准确地划定墒情的等级，但如果时间序列积累足够长，则可直接得到监测年相对于常年的偏离情况，具有一定的实用参考价值。