基于NDWI和卷积神经网络的冬小麦产量估测方法

刘峻明 周 舟 和晓彤 王鹏新 黄健熙

(1.中国农业大学土地科学与技术学院， 北京 100083； 2.中国农业大学信息与电气工程学院， 北京 100083)

0 引言

冬小麦是世界三大粮食作物之一，准确预测冬小麦单产及其空间分布对保障国家粮食安全和挖掘区域可利用的农业资源具有重要意义。目前利用遥感数据进行冬小麦产量预测主要有基于作物敏感波段反射率的统计回归方法[1]、基于遥感信息和作物模型的数据同化方法[2]等。统计回归方法中应用最广泛的是植被指数法，植被指数是由多光谱数据经线性或非线性组合构成的对植被有一定指示意义的数值，被广泛用于植被长势监测和作物估产领域。然而采用植被指数的经验模型通常是基于某一作物在特定条件下建立，例如归一化差值植被指数(Normalized difference vegetation index, NDVI)在植被达到一定覆盖度后增长缓慢，在高植被覆盖地区存在过饱和现象，往往导致模型缺乏普适性。数据同化技术的思想最早由文献[3]提出并在作物估产领域得到应用。通过数据同化技术使遥感数据参与到作物模型模拟过程中，可以提升模型在区域尺度的应用精度，但作物模型的参数标定以及数据同化系统的运行效率仍是当前应用的难点，在实际应用中很难快速得到大范围产量预测结果且预测精度还有待提高[4]。

随着机器学习技术不断发展，机器学习和深度学习方法被成功应用到多个领域，如图像识别、机器翻译、信号处理等[5]。传统的机器学习方法如支持向量机、决策树及随机森林等能够较好地解决非线性问题，并且有较好的效果[6-9]。在作物估产领域，研究已表明[10-12]深度学习方法通常能够得到比传统机器学习方法更高的精度。其中，卷积神经网络是使用最广泛的深度学习方法。相较于其他神经网络结构，需要的参数相对较少，且可直接将多维图像数据作为网络输入特征，信息损失小。目前在作物产量预测方面，文献[13]利用卷积神经网络和遥感数据对玉米产量进行回归预测，相比于支持向量回归算法更有优势并取得了较好的预测结果。文献[14]提出了一种新的卷积神经网络估产框架，使用遥感影像的直方图信息作为模型输入，既达到了数据降维的效果，同时也取得较高的产量预测精度，此方法具有迁移学习的能力[15]，在产量数据较少的地区同样有较好的效果，为基于卷积神经网络的估产方法在区域扩展性上提供了参考。

归一化差值水指数(Normalized difference water index, NDWI)是基于MODIS地表反射率2个近红外波段计算得到，它能快速有效地提取植被冠层的水分含量，并及时地响应植被冠层受水分胁迫的影响[16]。本文基于MODIS数据，以遥感植被特征直方图分布信息作为输入变量，应用卷积神经网络对冬小麦产量进行回归预测，并通过分析不同遥感植被特征在冬小麦产量估测上的表现，探讨NDWI在区域冬小麦产量估测上的应用，同时结合田间实测霜冻害资料进一步探究模型在霜冻害发生年份的表现，以期为县域尺度农作物产量预测提供有效方法。

1 研究区与实验数据

1.1 研究区域

河南省地处31°23′0″～36°22′0″N，110°22′0″～116°38′0″E(图1)，属北亚热带湿润气候与暖温带半湿润季风气候间的过渡气候，年平均降水量为500～1 000 mm，研究区冬小麦多种植冬性、弱冬性品种，一般在9月中旬至10月上旬播种，翌年5月底至6月初成熟[17]。独特的地理位置和气候条件使其成为我国冬小麦霜冻害高发区之一，霜冻害年际发生率高达60%[18]，严重地区可减产60%～70%[19]，对冬小麦的高产稳产构成威胁。

1.2 数据来源与处理

采用MODIS 8d合成地表反射率产品(MOD09A1)作为遥感数据源，分别提取其红光波段、近红外波段、蓝光波段、绿光波段的地表反射率数据作为遥感输入特征，空间分辨率为500 m，时间分辨率为8 d(2004-10-08—2017-06-18)，应用MODIS重投影工具软件(Modis reprojection tool, MRT)对影像的研究区域进行裁剪，投影转换为WGS84坐标系统，并利用MODIS土地覆盖类型产品(MCD12Q1)植物功能型分类方案(Land cover type 5)识别冬小麦种植区。产量数据来自《河南省农村统计年鉴》[20]中县级冬小麦单产数据，在2005—2017年间连续种植冬小麦的县市共有103个，共获得1 339个有效单产数据。冬小麦受灾、成灾、绝收面积等霜冻害资料来自于农业农村部种植业管理司历年自然灾害数据库(http:∥sjcx.fldj.agri.cn/moazzys/zaiqing.aspx)。

2 研究方法

2.1 遥感植被指数计算

遥感植被指数的选取主要围绕植物光合作用和冠层水分条件，因此对于这2个参数敏感的MODIS可见光、近红外波段都被纳入到候选波段中，其中包括计算NDVI所需的红光波段、短波近红外波段，以及计算NDWI所需的短波近红外波段、长波近红外波段。通过将提取的可见光、近红外波段的反射率数据作为输入特征，分别计算得到NDVI[21]、NDWI[16]、绿红植被指数(Green red vegetation index, GRVI)[22]、绿色归一化植被指数(Green normalized difference vegetation index, GNDVI)[23]、调整土壤亮度植被指数(Optimal soil adjusted vegetation index, OSAVI)[24]、土壤调节植被指数(Soil adjusted vegetation index, SAVI)[25]、改进型土壤调节植被指数(Modified soil-adjusted vegetation index, MSAVI)[26]8 d间隔的时间序列，时间序列长度为32，覆盖了冬小麦的整个生育周期(10月至次年6月中旬)。

2.2 直方图信息提取

因各县域遥感影像形状、像元数量差异较大，为实现样本结构的标准化，统计各县域不同遥感植被指数在冬小麦生育期时间序列上的像素分布直方图作为样本特征。NDVI、GRVI、GNDVI、OSAVI、MSAVI的像素值范围为(0，1)，NDWI的范围为(-0.25，0.5)，SAVI的范围为(0，1.5)，使用最大最小值法将NDWI和SAVI归一化至0到1之间，计算式为

式中Inorm——归一化后NDWI或SAVI值

I——NDWI或SAVI值

Imin——NDWI或SAVI的最小值

Imax——NDWI或SAVI的最大值

将所有植被指数值等间隔划分至32个区间内，统计各区间像素百分比，得到频率直方图。

2.3 单产去趋势处理

2003—2017年河南省小麦产量总体呈增长趋势(图2)，主要得益于小麦品种的改进、管理技术的提高和农业政策的改革[27-29]，为去除由技术进步、经济社会发展等因素引起的冬小麦单产变化，本文采用5 a滑动平均法计算得到其趋势单产，用实际单产减去趋势单产对研究区单产数据进行去趋势处理。

2.4 卷积神经网络

卷积神经网络是一种前馈型的神经网络，其主要组成结构包括输入层、卷积层、池化层、全连接层、激活函数等。以研究区逐年各县的不同遥感植被指数的像素直方图及其对应的县域冬小麦单产作为模型的样本数据集，采用卷积神经网络模型进行训练和验证。本文的卷积神经网络输入层为32×32×7的矩阵，各卷积层的卷积核个数依次是128、128、256、256、512、512、512，卷积核尺寸都是3×3，滑动步长分别为1、2、1、2、1、1、2。同时，在每一个卷积层上进行批归一化和线性整流函数(ReLU)激活操作，并在全连接层加入随机失活(Dropout)操作。具体参数设置为：使用方差缩放方法初始化网络权重，偏差初始化为0，初始学习率为0.001，Drouput设置为0.5，使用自适应矩估计(Adam)优化器，运行一次输入的样本数为32。

2.5 模型精度评价与验证

采用决定系数(Coefficient of determination，R2)、均方根误差(Root mean square error，RMSE)和平均绝对误差(Mean absolute error，MAE)3个指标对模型拟合程度优劣进行评价[30]。

3 结果与分析

3.1 不同植被特征估产精度比较

为评估NDWI在冬小麦产量估测上的表现，利用经遥感植被指数计算直方图信息提取得到的NDVI、NDWI、GRVI、GNDVI、OSAVI、SAVI、MSAVI 7个遥感特征，以2005—2014年1 030组数据作为训练样本，分别将7个遥感特征的直方图信息作为特征集，实际单产作为目标变量构建基于CNN的回归预测模型，以2015—2017年309组数据作为验证样本输入模型，对比分析不同遥感植被指数对预测精度的影响并对去趋势前后模型预测精度进行比较，结果见表1。从表1中可以看出，相对于植被指数NDVI、SAVI、OSAVI、GNDVI、MSAVI、GRVI，NDWI表现出更好的预测效果，单产去趋势前后的NDWI对产量的预测精度均高于NDVI、SAVI等植被指数。

表1 不同植被指数去趋势前后验证结果Tab.1 Verification results of different vegetation indexes

在冬小麦主要生育期，水分是限制冬小麦叶片生长的重要因素，通过影响冬小麦的光合作用来限制其子房发育灌浆后期，茎叶中的营养物质通过植株体内的水运输到籽粒中，水分会直接影响营养物质的输送，进而影响冬小麦产量[31]。相比于去趋势前模型预测结果，NDWI在单产去趋势后R2提高了0.05，实测单产和预测单产的散点图如图3所示，趋势线与1∶1线交于0点附近，大部分样本聚集在1∶1线周围，R2最高达到0.79，MAE和RMSE分别为482、637 kg/hm2，主要是由于去除了产量年际间社会经济因素的影响。

3.2 不同生育阶段NDWI对产量的影响

为进一步分析不同生育阶段NDWI对产量的影响，将冬小麦全生育期划分为6个时间段，分别为10月8日—11月25日、12月3日—2月26日、3月6日—3月30日、4月7日—4月30日、5月1日—5月17日、5月25日—6月18日，大致对应冬小麦出苗—越冬、越冬—返青、返青—拔节、抽穗—灌浆、乳熟—成熟阶段，分别以各时间段NDWI作为样本特征输入，去趋势单产为目标变量，验证结果见图4。从图4中可以看出，抽穗—灌浆阶段模型预测精度最优，MAE和RMSE分别为552 kg/hm2和759 kg/hm2，R2最高达到0.74，说明该阶段影像反映的植被状况对产量的影响最大，乳熟—成熟以及返青—拔节阶段次之，出苗—越冬阶段预测效果相对较差。冬小麦在返青—拔节阶段主要进行营养生长，该阶段是决定穗数和粒数的关键时期，但其生长特征并不能完全反映产量形成器官的干物质积累过程[32]，因此该生育阶段的模型精度较低。抽穗—灌浆阶段有机物从营养器官转移到籽粒，该阶段NDWI与冬小麦千粒质量密切相关[33]，故此阶段估产精度最高。在冬小麦乳熟—成熟阶段，冠层和茎秆的营养物质向籽粒转移，叶片中的叶绿素含量下降，与产量的相关性变弱[34]，故在成熟后期估产模型精度下降。

为确定研究区冬小麦产量估测的最佳时间，根据NDWI在不同生育阶段的模型验证结果，对抽穗—灌浆阶段进一步划分为3个时间段，分别为4月7—14日、4月15—23日、4月23—30日，分别对各时间段的NDWI进行训练并预测2015—2017年对应时段的冬小麦产量，结果见表2。从表2中可以看出，4月23—30日的NDWI对产量的决定系数可达到0.72，MAE和RMSE分别为566、763 kg/hm2。这主要是由于籽粒最终产量主要来源于抽穗—成熟阶段叶片的光合产物，而灌浆后期是籽粒干物质积累最旺盛的时期，地上干生物量中籽粒比重较大[35]。

表2 抽穗—灌浆生育阶段验证结果Tab.2 Verification results of NDWI in heading—filling growth stage

3.3 模型预测误差空间分布

对NDWI模型精度进行逐年验证，2015—2017年的精度验证结果如图5所示，每一年份分为4幅子图，包括单产估测值与实测值散点图、县域实测单产分布、模型估测单产分布以及误差分布。从空间分布上来看，模型估测单产与实测单产图中高产区和低产区分布基本一致，东部单产最高，中部次之，西部最低，整体呈东高西低。从误差分布图中可以看出，大部分区域误差在±300 kg/hm2内，估测误差大于900 kg/hm2的县主要分布在西部和北部山区与东部黄淮海平原交界处，这些县的单产较低，对应于散点图中的低产部分，低产县的模型估测单产普遍高于实测单产。造成此误差的主要原因有：①低产区的单产样本数量较少，从散点图中可以看出，中产和高产数据的密度较大，而低产数据则相对较少，而样本分布不平衡是造成机器学习和深度学习预测偏差的主要原因之一。②低产区域多为山区-平原过渡地带，地势西高东低，地形较为复杂，农田小气候与地形因素会对冬小麦产量产生一定的影响，从而影响到模型估测效果。

4 讨论

4.1 霜冻害影响下NDWI模型预测精度

考虑到河南省霜冻害发生次数较频繁，造成小麦减产，有效地预测霜冻害影响下的冬小麦产量，对于冬小麦灾害预警、稳产高产具有重要意义。为探究NDWI模型在霜冻害影响下的预测效果，使用留一年法对模型精度进行逐年验证，结合2005—2017年农业农村部种植业管理司历年自然灾害数据库河南省霜冻害资料，验证结果如图6所示。由图6可以看出，模型预测精度在2005—2017年整体上呈现波动变化的趋势，有无霜冻发生年份均维持在一较高水平，在预测精度达到最高的2006年、2008年和2013年均有霜冻发生。霜年平均R2约为0.78，平均RMSE、MAE分别为682、527 kg/hm2，这表明模型在霜冻害影响下仍然能保持较好的预测效果。

2013年4月15—23日，河南省商丘地区发生春季霜冻事件，气温骤降15℃以上，此时小麦正处于籽粒形成的关键时期，受冻后会导致明显的缺粒现象，对当地冬小麦产量造成极大影响。根据田间霜冻害调查结果可知，该年商丘地区平均穗粒数减少率达到40%左右，平均减产率达到39.6%左右[36]。本次春霜事件为探索NDWI在冬小麦低温胁迫下的变化特征并进一步验证霜冻害影响提供了一个理想的案例。冬小麦发生冻害时，冠层含水量上升，植株体内发生了失水情形。当极端低温超过了小麦的耐寒能力时，植株细胞原生质体内以及细胞间隙间的水分发生放热凝固现象，植株体内的水由于固结失去流动性，无法将营养物质运送到各个器官，且细胞结构发生不可逆的损坏，对水分的控制能力下降，当白天气温上升时，植株体内的结冰开始融化，水分开始外渗故造成失水[37]。图7为2013年3月6日至5月17日商丘地区冬小麦种植区域的NDWI和NDVI的时间序列变化曲线，NDWI与NDVI整体变化趋势相同，均呈先上升后下降的变化态势。在返青初期，NDWI值较低，4月15日至4月23日和5月1日至5月9日间，NDWI有2次明显的上升过程，结合气象站提供的日最低温度和日降水量数据，5月6—9日间有明显降水现象，可以解释在此期间NDWI数值的突然升高，而4月15—23日间商丘地区出现了大幅度降温现象，4月21日最低气温达到0℃，低温持续近3 d，同时期内降水量较低，因此降水对于植被水分条件的影响可以忽略。在没有降水输入的情况下，NDWI数值的异常升高可认为是由冻害引起的。

4.2 模型局限性

从估产验证结果来看，NDWI能够很好地反映植被最终的生长状态，其R2与NDVI相比，提高了0.06，且在霜冻害影响下仍能保持较好的预测效果，但本研究仍存在一些不足以及值得进一步探索的地方：

(1)使用植被指数的直方图信息作为模型输入，模型在产量数据较少的地区同样有较好的效果，为基于CNN的估产方法在区域上的扩展性提供了参考。但直方图信息的提取需要样本区域内有足够多的有效像元个数，因此可能并不适用于像元尺度上的单产估测。

(2)采用的去趋势方法对各年份研究区全省范围的趋势产量进行计算，而非单独计算各年份各县的趋势，而不同县的单产增长趋势存在差异，故单产在空间上的变异性仍然存在。

(3)不同生育阶段划分的时间长度不一致，例如越冬—返青阶段时间约为90 d，而返青—拔节阶段、抽穗—灌浆阶段时间跨度相对较小，不足30 d，所包含的影像序列长度不一，模型存在不确定性。

5 结论

(1)NDWI能够很好地反映植被最终的生长状态，在冬小麦生育早期的产量预测上表现出更好的预测效果，R2最高可达0.79，MAE和RMSE分别为482、637 kg/hm2，因此适合作为冬小麦估产指标。

(2)NDWI在抽穗—灌浆阶段对冬小麦最终产量影响最大，NDWI在4月23—30日时间段内对产量的决定系数可达到0.72，综合对比可知抽穗—灌浆阶段NDWI对冬小麦最终产量影响较大。

(3)模型预测精度在霜冻发生年份，R2最高可达0.83，这表明模型在霜冻害影响下仍能保持较好的预测效果。