基于CNN-RNN网络的中国冬小麦估产

赫晓慧，罗浩田，乔梦佳，田智慧※，周广胜

（1. 郑州大学地球科学与技术学院，郑州 450001； 2. 郑州大学信息工程学院，郑州 450001；3. 中国气象科学研究院郑州大学生态气象联合实验室，郑州 450001）

0 引 言

农业是人类社会发展和稳定的基础，而中国是人口大国，粮食安全尤为重要[1]。小麦是中国的主要粮作物之一，据统计，2020年小麦播种面积为23.38万hm2，总产量13 425万t，占全国粮食总产量的1/5，因此，准确、及时地预测中国小麦产量对全国粮食安全乃至世界农业发展至关重要[2-3]。

近几十年来，国内外学者围绕着作物产量预测、粮食安全等问题展开了深入而广泛的研究。许多学者基于经验统计模型[4]和作物生长模拟模型[5]等传统的估产模型预测作物产量，这些模型通常对于某一区域的特定农作物产量预测精度较高，但很难推广到拥有大规模种植区域以及需要进行多农作物产量预估的地区和国家[6-8]。

随着传感器的发展，遥感卫星可以获得大量高质量、高时空分辨率的图像，为大尺度农作物产量预测提供了可能。遥感卫星捕获的光谱信息能充分反映农作物的生长状态，以往的研究者大多倾向于从影像中提取作物相关指标来建立与作物产量之间的关系[9-11]，如归一化差异植被指数（Normalized Difference Vegetation Index，NDVI）、增强植被指数（Enhanced Vegetation Index，EVI）等[12-14]，由于农作物生长环境复杂，植株状态多变，仅依靠植被指数可能无法完全反映影像中农作物的生长特征。

为了充分利用遥感图像中的光谱信息，提升农作物估产的准确性，研究者逐渐将深度学习应用到了农作物估产中。如Qi等[15-16]通过构建两个独立分支的卷积神经网络（Convolutional Neural Networks，CNN）结构分别处理高分辨率RGB图像和多光谱影像，整合了影像中的空间- 光谱信息进行估产，并发现水稻成熟期有较高精度的产量预测；Nevavuori等[17]构建了基于无人机获取的NDVI和RGB图像的CNNs作物产量预测模型，结果表明基于RGB图像的CNN模型优于NDVI数据。尽管基于CNN的方法取得了较好的结果，但是这些研究仅仅基于农作物某个特定生长期的影像对产量进行预测，如Son等[18-20]发现在农作物的生育后期的预测精度最高。然而，农作物的生长是一个生物量长时间累积的过程，包括多个生育期。已有研究表明，在作物生长期内的不同时相遥感图像中，作物信息与其产量关系的显著性具有较大差别[21-25]。因此，综合利用农作物生长过程中的多时相遥感影像进行估产，有助于掌握不同时间阶段的农作物生长状况，提升产量预测精度[26-30]。虽然CNN被广泛用于从遥感影像中提取光谱特征[31]。但是由于结构的限制，CNN无法对时间序列进行建模[32-34]。因此，仅利用CNN无法从多时相影像提取中作物生长的特征。因此，如何在提取遥感影像空谱特征的基础上，为长时间序列的遥感影像进行建模，是提升农作物产量预测的关键问题。

为解决上述问题，本文将全国冬小麦主产区作为研究区域，并选取研究区2001－2018年县级冬小麦生育期的时间序列遥感影像作为数据源，同时设计了一个端到端的深度学习估产网络从多时相遥感影像中自动学习空间-光谱-时间联合特征。具体来说，首先将每个时间节点的冬小麦影像独立地嵌入卷积子网模块，在卷积子网充分提取影像中的空间-光谱特征之后，将其按时间序列传递到递归子网单元，整合冬小麦整个生育期的时间信息。本文所提出的农作物产量预测网络旨在能够从时间序列遥感影像中多尺度、全方位地挖掘作物生长特征与冬小麦产量之间的复杂关系，以实现大尺度的冬小麦产量预测。

1 研究区及数据

1.1 研究区

本文选取全国冬小麦主产区作为研究区域，包括14个省和自治区：河北、陕西、江苏、浙江、安徽、山东、河南、湖北、四川、云南、陕西、甘肃、宁夏回族自治区和新疆维吾尔族自治区，共1 713个区县。

1.2 数据源

1.2.1 MODIS影像

本文采用中等分辨率成像光谱仪（MODIS）卫星图像，该卫星影像有较高的时间分辨率，能获取足够多检测冬小麦产量的多时相影像。数据下载自Google Earth Engine（https://earthengine.google.com/）。本文主要使用3种MODIS产品数据，其时间覆盖为2001—2018年。MOD09A1 V6产品提供了500 m和8 d的地表反射率，包括7个波段，如红波段、蓝波段和近红外波段等，其采集的地表反射率影像可以反映作物的生长状态和环境[35-36]；MOD11A2 V6产品可在1 km×1 km的网格中提供每8 d的昼夜地表温度影像，昼夜地表温度与作物冠层温度联系紧密，有利于进行产量预测[37]。MCD12Q1 V6产品以年为单位提供全球土地覆盖类型，其中包括草原、农田、城市、水体等17类土地类型，其空间分辨率为500 m；数值比例表示把原始光谱值从整型转化为浮点型的比例系数；数值范围表示原始数据中光谱值的最大值和最小值范围，超出数值范围的则是无效数据。

1.2.2 产量数据

冬小麦产量数据是以区县为单位的统计数据，其中2001－2008年的冬小麦各区县年平均产量及播种面积数据来源于中国科学院资源学科创新平台[38]（http://www.data.ac.cn），2009—2018年数据收集自各省、市、区县统计年鉴。

1.3 数据预处理

为统一影像分辨率，将MOD11A2的空间分辨率重采样为500 m，再利用MCD12Q1中农作物掩膜数据提取各区县的农作物种植面积，由于农作物掩膜不能区分具体的农作物类型，因此需要结合作物产量数据，剔除产量较低和影像中农作物像素点较少的县，即将县域作物像素小于2×2的影像剔除，当影像中作物像素点太少，估产网络难以从中提取有用信息，反而会被当成噪声数据，进而影响网络的性能。对MOD09A1和MOD11A2各个波段进行提取、融合，即利用GDAL库分别读取MOD09A1和MOD11A2影像，将其波段按顺序排列，共计9个波段，包括MOD09A1的7个波段和MOD11A2的2个波段，最后输出融合后的影像。

2 研究方法

2.1 估产样本构建

本文选取2001—2018年的冬小麦影像数据作为样本。一般来说，冬小麦在10月中旬播种，次年7月初收获，由此确定影像数据的时间范围从每年的第289天到次年的第180天，共30个时间步。冬小麦生育期包括：播种期、出苗期、分蘖期、越冬期、返青期、起身期(生物学拔节)、拔节期、孕穗期、抽穗期、开花期、灌浆期、成熟期，每个生育期约21 d，把每两个生育期划分一个时间节点，即每两个生育期对应5个时间步：播种-出苗期（0～5）、分蘖-越冬期（5～10）、返青-起身期（10～15）、拔节-孕穗期（15～20）、抽穗-开花期（20～25）、灌浆-成熟期（25～30）。结合县边界数据，裁剪每个县

所有时间步长t∈(0,30)的影像其中，b代表参与农作物估产的波段数，id代表区县的国家行政区划代码，y代表预测年份。由于各县级单位冬小麦种植面积不同，为统一网络的输入，统计影像内农作物像素点个数，取32×32作为阈值。当县域作物像素小于阈值时，将除作物像素的其余点填充为0；当县域作物像素大于阈值时，剔除零星分布的农作物像素，这部分像素点可能会被当成噪声，取农作物比较集中的像素作为网络输入；以确保影像整体信息完整。为提升模型精度和收敛速度，对冬小麦数据集进行标准化。

式中μ表示均值，σ表示方差。计算出整个数据集的均值和方差，并应用上述公式对影像中每个像素点进行标准化。

2.2 CNN-GRU融合网络

如图1所示，所提出的CNN-GRU估产模型由卷积子网、递归子网和全连接层组成。为提取多时相影像的联合光谱-空间-时间特征，卷积子网从每个时间步的影像中提取空间-光谱特征映射，最终得到30个时间步的特征向量。在卷积子网的基础上，递归子网以卷积层的特征作为输入，利用门控循环单元（Gated Recurrent Unit，GRU）挖掘不同层次的特征间的依赖关系，整合CNN框架部分输出的特征向量。最后利用全连接层进行产量预测。

2.2.1 基于卷积子网的空间光谱特征提取

农作物的生长主要伴随光谱变化，因此从遥感影像中准确地挖掘光谱信息是农作物产量预测地关键。同时，由于农作物分布通常呈现聚集性，农作物像素之间的空间相关性也十分重要，而CNN作为深度学习的一个重要分支，其独特的卷积结构能提取遥感影像中丰富的光谱-空间特征[39-40]。

卷积子网包括4个卷积层、4个池化层和一个全连接层，如图2所示。输入数据为每个时间步的冬小麦影像。在每一张影像上利用3×3的卷积核的点积之和来提取固定大小的光谱特征，并利用窗口滑动以覆盖整个空间维度，以提取冬小麦的空间-光谱信息。本文选择线性整流函数（Rectified Linear Unit，ReLu）作为激活函数，不仅能加速网络收敛，还能避免梯度爆炸或梯度消失。最后输出特征图表示为

式中f(·)表示ReLu激活函数，xi表示上一个卷积层的输出，wi表示第i-1层和第i层之间的权重，bi表示贝叶斯偏置项。

通过卷积层提取特征后，为减少计算量，本文选择最大池化函数来处理从卷积运算中获得的特征映射结果，其中池化窗口为2×2，池化步长为2。

式中xpool表示池化层的输出。卷积子网最终将每个时间步的冬小麦影像映射为高维的空间-光谱特征{f1,f2,…,f30}，并将其作为递归子网的输入。

2.2.2 基于递归子网的时间相关性建模

在农作物生长过程中，其植株形体特征会随生育期变化而改变，卷积神经网络（CNN）仅能提取单一生育期的农作物生长特征，而通过农作物不同生育期植株状态的变化来挖掘深层次的时序特征尤为重要。循环神经网络（Recurrent Neural Network，RNN）可以很好地关联上下文信息，尤其擅长针对时间序列数据建模[41-42]。但传统RNN面临着梯度消失的问题，难以对长时间序列数据建模。因此，本文使用GRU代替传统RNN网络，GRU是循环神经网络的变体，适用于挖掘长时间序列数据的内在联系，并且有效解决了梯度消失的问题[43]。

因此，在CNN的基础上，本文进一步使用GRU网络构建农作物生长过程中的时间关系。GRU由重置门rt、更新门zt和隐藏状态ht组成，将从不同时间步获取的高维空谱特征xt作为输入，其中隐藏状态ht用来处理连续的冬小麦时间特征，在考虑先前冬小麦生长状态的情况下学习当前的特征；重置门用于控制前一时刻的隐藏单元ht1-对当前输入xt的影响；更新门确定在当前隐藏状态ht之前，保留多少先前时刻的状态信息。第t∈(0,30)时间步各个门的迭代参数可表示为

式中xt为第t时间步的输入分别表示更新门、重置门、候选状态、隐藏状态，σ是sigmoid激活函数，表示权重矩阵和偏执向量，☉表示矩阵运算。这种结构使GRU能够根据冬小麦先前的生长信息，结合新的输入进行自我更新。因此，GRU可以充分挖掘各个时间步之间的内在联系，从而有效反映农作物生长过程中的相互关系。

2.3 loss函数和评价指标

在网络最后，将所有的时间输出展平成一个向量，利用一个神经元的全连接层输出冬小麦预测产量，并在每个池化层和全连接层添加比例为0.5的dropout结构，以避免模型出现过拟合现象。本文使用自适应矩估计（adaptive moment estimation，Adam）优化器训练网络，并通过L2损失函数[44]来计算预测产量（preid）与真实产量（realid）之间的训练损失值，如公式所示：

本文使用均方根误差（RMSE）、平均绝对误差（Mean Absolute Error，MAE）和决定系数（R2）来评估模型预测精度[45-47]。

2.4 模型训练

本文使用python对影像数据进行预处理，在基于TensorFlow的深度学习框架下进行模型的搭建、训练和测试，并安装CUDA等运算平台以搭建GPU加速环境。实验的硬件仿真环境为Intel(R) Core(TM) i7-6900K CPU，使用NVIDIA GeForce GTX 1080 GPU对模型进行加速，显存大小为8G。

为实现全国冬小麦主产区的产量预测，本文对CNN-GRU估产模型进行训练。在训练阶段，以全国冬小麦主产区整个生育期的遥感影像作为训练数据，年平均产量为标签，年份为2001－2015年，时间步t∈(0,30)。为确保无偏性，每个批次从训练集中随机抽取32个样本进行训练。初始学习率为0.000 1，分阶段逐步递减为初始值的10倍，达到迭代次数时，训练停止。在测试阶段以2016－2018年全国冬小麦主产区的年平均产量和影像作为测试数据。

3 结果与分析

3.1 数据贡献度分析

为评估MODIS影像中地表反射率（7个波段）和地表昼夜温度数据对农作物估产的影响程度。本文对温度数据和反射率数据分别进行训练，并将反射率7个波段依次去除一个，使用其余波段进行训练。如表1所示（以RMSE为例，单位kg/hm2），其中删除波段1表示删除第一个波段，仅使用第2～7个波段，以此类推，最后一列波段1～7表示使用地表反射率的所有波段。结果表明：1）只利用反射率训练的CNN_GRU估产模型RMSE年平均值分别为841 kg/hm2，优于仅使用温度数据的预测结果，说明反射率数据相较于温度数据包含了更多的农作物信息，如红光波段对绿色植被具有强吸收性，对农作物的覆盖度、生长状态比较敏感。2）使用反射率所有波段训练的模型效果最好，说明反射率7个波段均对均对农作物产量预测起决定性作用。在分别删除第1、2、7个波段后，RMSE下降了8%左右，而删除第3个波段时，RMSE仅下降了3%，说明第1（红光波段）、2（近红外波段）、7（短波红外波段）相较于其他波段对农作物估产更为重要，而第3个波段（蓝光波段）主要反映植被的下垫面土壤与植被之间的差异，可能相较于其他波段不能更直接的反映作物的生长状态。

表1 地表反射率（7个波段）和地表昼夜温度数据对农作物估产的影响Table 1 The effect of surface reflectance (7 bands) and surface day and night temperature data on crop yield estimation kg·hm-2

3.2 估产模型性能检测

为验证CNN-GRU估产模型的优越性，将其与使用最广泛的多个农作物产量预测模型进行比较，如深度学习框架下的CNN、GRU和传统机器学习算法中的SVR、DT和RF。其中所有估产模型均使用同一个数据集。但CNN难以处理长时间序列数据，则将时间序列和波段相叠加作为CNN的输入，即（32,32,30×9）。而GRU网络需要将输入数据展平为时间序列向量（30,32×32×9）。此外，由于传统机器学习模型本身的结构限制，在输入之前需要把输入数据展平为一维向量，但30×32×32×9=276480维度太高，容易导致“维度灾难”，因此，使用影像的NDVI作为传统机器学习模型的输入。各个模型的关键参数设置如下：

1）卷积神经网络：包含4个卷积层、4个池化层、1个全连接层，批处理大小为32，初始学习率为0.000 1。

2）门控循环单元模型：由128个神经元的GRU和全连接层组成，批处理大小为32，初始学习率为0.001。

3）支持向量回归：使用径向基（RBF）作为SVR的核函数，惩罚因子C=2,5,10,100，松弛变量γ=10-1,10-2,10-3。

4）决策树：包括树的最大深度max_depth=3,4,5,6,7，叶子节点最少样本数min_samples_leaf=10,15,20,30，内部节点再划分所需最小样本数min_samples_split=5,8,10。

5）随机森林：包括子树的数量num_estimators=50,80,100，每棵树的最大深度max_depth=3,6,9。

在本节中，使用不同的估产模型分别预测2016—2018年冬小麦产量，为确保比较的公平性，所有估产模型均在相同的软硬件环境下训练和测试。分析表2可知：1）所提出的CNN-GRU估产模型年平均RMSE为818.3 kg/hm2，相较于CNN、GRU、SVR、RF和DT分别降低了20.13%、18.81%、29.51%、34.84%和36.57%，年平均MAE为560 kg/hm2，相较于其他模型分别降低了25.46%、22.39%、37.09%、39.95%和42.26%。表明CNN-GRU模型能准确地预测冬小麦产量。但对于部分低产量的区县（如新疆维吾尔自治区的伊宁县和和硕县）和部分高产量的区县（陕西省西安市的灞桥区和河南省南阳市的新野县）CNN-GRU模型预测误差偏大（如图3所示）。产生这种误差的原因可能是冬小麦数据集中低产量和高产量的样本数占偏低，导致模型对于部分低产量和高产量的地区泛化能力较差，预测的误差较高；对于低产量的新疆维吾尔自治区海拔高、环境复杂等不确定性因素导致预测精度降低；对于高产量的灞桥区和新野县来说，地处平原，均为冬小麦气候适宜区，种植面积较大，而模型输入中没有表达区域异质性的参数，导致低估了部分高产量地区的冬小麦产量。2）CNN和GRU的性能相当，年平均RMSE分别为920、996 kg/hm2，均高于所提出的CNN-GRU估产模型，表明单一的神经网络在作物估产领域效果不佳。CNN能有效地提取影像中空间-光谱特征，但作物生长随时间变化，农作物在不同的生育期有不同的形态特征，时间信息的缺失严重影响了CNN的整体性能；而GRU能整合长时间序列的作物生长信息，但不能捕捉复杂的空间关系。因此所提出的CNN-GRU估产模型可以有效地提取冬小麦生长期的空间-光谱-时间信息，并将其应用于冬小麦产量预测中。3）深度学习网络在性能上均优于传统的机器学习模型，在冬小麦产量预测这一复杂问题的研究中，深度学习网络更能挖掘冬小麦影像中的时空特征。

表2 不同估产模型性能对比Table 2 Performance comparison of different yield estimation models

3.3 不同生育期冬小麦预测

本节分别对作物累计生育期的时间步（0～5，0～10，0～15，0～20，0～25，0～30）和单一生育期的时间步（0～5，＞5～10，＞10～15，＞15～20，＞20～25，＞25～30）进行单产预测，以评估CNN-GRU网络模型预测冬小麦产量的最佳时间窗。

如图4所示，CNN-GRU估产模型在生育期初期性能较差，2016－2018年平均RMSE和MAE分别为1 109、840.6 kg/hm2，这是由于估产模型在生育期初期缺乏足够的作物生长信息和环境信息，但随着时间的推移，越来越多的影像信息被整合到估产网络中，模型的性能逐步提高。1）单一生育期的预测精度很低，相较于全生育期影像（0～30）低了25%至28%，这是由于单一生育期仅使用了部分时间步的影像数据，缺乏前后生育期冬小麦的生长状态，不能充分发挥CNN-GRU估产网络整合长时间序列影像时空特征的优势。2）返青-起身期和拔节-孕穗期两个时间步的RMSE相差最大，约10%左右，拔节-孕穗期相较于返青-起身期预测精度明显提升。在返青-起身期时，麦苗开始生长春生第一叶，茎部第一节开始伸长但为完全露出地面。而在拔节期，冬小麦植株茎部露出地面1.5～2 cm，生长特征明显。表明CNN-GRU模型能有效地捕捉冬小麦的各个生育期的特点及其生长状态的变化。3）在灌浆-成熟期时，CNN-GRU估产模型精度最高，2016－2018年的平均RMSE和MAE分别为817、556 kg/hm2，R2为0.76，均优于其他各个生育期。在冬小麦生育期后期，所有的作物生长信息均被整合到CNN-GRU估产模型中，该估产模型能精确地预测全国冬小麦主产区产量。4）2016－2018年冬小麦拔节-孕穗期和抽穗-开花期估产模型预测精度相当，年平均RMSE分别为849、823 kg/hm2，仅比成熟期低3.9%和0.7%，说明CNN-GRU估产模型在冬小麦生育期的中后期已经能准确地预测其产量，即所提出的估产模型有能力提前2个月预测全国冬小麦主产区产量。

3.4 CNN-GRU模型鲁棒性检验

为对提出的CNN-GRU估产模型进行鲁棒性检验，将2008－2018年估产样本逐年作为验证集，其余年份作为训练集进行训练，如表3所示。结果表明，2008－2018年CNN-GRU估产模型的RMSE和MAE平均值为727、534 kg/hm2，R2为0.805。在不同训练集和验证集下，预测产量与统计产量值离散程度基本一致，表明CNN-GRU估产模型有较高的鲁棒性。

表3 CNN-GRU模型逐年精度评估Table 3 Annual accuracy evaluation of CNN-GRU model

综合分析发现，2010年产量预测误差最大，RMSE和MAE分别为982、752 kg/hm2，R2仅0.617；2012年模型精度最高，RMSE和MAE分别为604、446 kg/hm2，R2为0.867。2009、2011、2013和2014年模型预测精度较高，RMSE和MAE在650、500 kg/hm2左右；2016－2018年预测效果一般，RMSE在820 kg/hm2左右。不同年份模型性能不同的主要原因有：1）在2010年，中国西南地区遭遇特大干旱，北方冬小麦主产区出现持续低温天气，在极端天气情况下，统计数据极值差异较大[41]。如图5a所示，2010年全国冬小麦主产区县级尺度产量预测出较多离群值，而多数离群值的预测产量都大于统计产量，导致模型整体精度降低。2）2012年的RMSE和MAE分别为604、446 kg/hm2，R2为0.867。如图5b所示，预测产量与统计产量接近，仅有零星几个离群值。3）2016－2018年的收集的产量数据较少，可用于训练和测试的区县数量仅约650个，而其余每年的区县数据约1 500个。并且随着基因技术和种植技术的发展，近年来作物单产飞速提升，2016－2018年的低产量地区数量较少，导致数据分布不均衡。而2009、2011、2013和2014年的产量分布比较均衡，导致模型在不同年份预测的精度不同。

4 结 论

作物产量预测在农业管理中具有重要意义。近年来的研究表明，遥感是一种有效的估产方法。本文提出了一种基于多时相遥感影像的CNN-GRU估产模型，利用CNN从多时相遥感影像中提取丰富的空间-光谱特征，在此基础上，GRU自适应学习冬小麦生育期各阶段之间的时间依赖。试验证明该模型能从多时相遥感影像中提取联合空间-光谱-时间特征，挖掘作物生长特征与产量之间的复杂关系，并能够对收获前期的冬小麦产量进行准确预测。结论如下：

1）深度学习技术在产量预测上有巨大潜力，对比结果表明深度学习方法优于传统机器学习方法。同时，与单一的深度学习网络CNN和GRU相比，本文所提出的CNN-GRU估产模型精准预测了2016－2018年中国冬小麦主产区的产量，其精度均优于其他估产模型。说明该模型能有效地提取反映作物生长信息的空间-光谱-时间特征，解决了单一神经网络提取作物生长特征不充分问题，进一步提高了大范围冬小麦单产估算效率。

2）基于混合神经网络的遥感估产模型鲁棒性较高。本文对2008－2018年数据集逐年样本训练，所提出的CNN-GRU估产模型RMSE和MAE平均值分别为727、534 kg/hm2，R2平均值为0.805。该模型在不同数据集上预测产量与统计产量值离散程度基本一致，CNN-GRU估产模型性能稳定，预测精度高，可以作为一种高效的方法为全国提供常规可行的冬小麦产量监测和预报。

3）在冬小麦成熟期CNN-GRU估产模型性能最高，但拔节-孕穗期和抽穗-开花期同样取得较高精度，RMSE分别为849、823 kg/hm2，R2均高于0.7，与成熟期预测结果几乎没有差异。说明本文提出的CNN-GRU估产模型有能力提前2个月预测县级尺度的全国冬小麦主产区产量，为冬小麦生育早期产量预测提供了新的思路和方法。