基于CNN-S-GPR的宁夏枸杞高光谱影像估产方法

刘立波 王 涛 张 鹏

(宁夏大学信息工程学院， 银川 750021)

0 引言

枸杞作为宁夏特产，是国家地理标志产品，其产量是重要的经济信息，预测枸杞年产量对于种植管理和政府决策是一项难度较大但必不可少的研究。传统的农作物估产采用人工区域调查方法，从农学、气象学等不同角度建立作物估产模式。该方法速度慢、工作量大、成本高，且不利于时空动态监测。近年来，遥感技术因其覆盖范围广、重返周期短、获取成本低等优势，被广泛运用于农作物估产，成为遥感与农业交叉的研究重点[1]。

目前，遥感估产方法按模式不同主要分为2类：基于机器学习和基于深度学习的估产方法[2]。基于机器学习的方法将归一化植被指数(Normalized difference vegetation index, NDVI)、增强型植被指数(Enhanced vegetation index, EVI)等植被指数作为表征作物产量的重要信息，进而采用机器学习方法建立模型预测农作物产量。文献[3-5]利用NDVI、EVI等植被指数，基于Lasso、支持向量回归(Support vector regression, SVR)等回归算法构建了产量预测模型，达到了比较高的估产精度。但这类方法NDVI、EVI等指数仅利用近红外波段和红外波段2个波段计算得到，忽略了其余波段重要信息，不利于估产精度的提升，对作物长势反映不理想。所以，文献[6-10]提出了作物生长模型(Aquacrop、WheatSM、WOFOST等)，根据气象条件、土壤条件、叶面积指数(Leaf area index, LAI)以及作物栽培管理措施，定量描述作物生长、发育、籽粒形成等动态过程。这类模型虽然模拟作物生长机理性强、估产精度高，但过分依赖实测数据和人工经验，使工作量和成本大幅度上升，具有一定局限性。基于深度学习的估产方法，利用神经网络非线性映射能力拟合样本与观测值之间的关系，从而预测出产量。如高云[11]采用DNN(深度神经网络)准确预测了春小麦产量。针对DNN等前馈神经网络难以处理时序数据的问题，文献[12-13]采用长短期记忆网络(Long short-term memory, LSTM)对长期信息依赖进行处理，大大提升了估产精度，但仍需大量手工时间序列特征参与计算。文献[14-15]采用卷积神经网络(Convolutional neural networks, CNN)直接从时序影像中提取特征，摒弃了原有手工制作特征，进而利用全连接网络预测了中国北方冬小麦产量，并取得了很好的效果，但是仍存在以下问题：①基于CNN的遥感估产方法虽然简化了特征提取操作，却未对通道间依赖关系进行建模，即无法自动获取每个特征通道的重要程度，以强调有效信息并抑制无效信息。②CNN回归模型仅能提取遥感影像时间序列特征，却无法整合特定位置空间下对作物生长有影响的特征，因此具有一定局限性。

综上所述，本文在CNN模型的基础上，提出一种基于CNN-S-GPR的高光谱影像估产模型。通过CNN提取影像特征，融合通道注意力机制SENet(Squeeze-and-Excitation network)对CNN卷积层通道进行加权，以对通道间的依赖关系建模。接着在网络中引入高斯过程回归(Gaussian process regression, GPR)方法整合影像特征和空间位置特征，以宁夏回族自治区中宁县、贺兰县等16个县为估产单元，基于多级Modis高光谱影像数据，对各县域枸杞产量进行预测。

1 研究区概况

研究区域位于宁夏回族自治区(35°14′～39°23′N， 104°17′～107°39′E)，地处黄河水系，属北温带大陆性季风气候区，是全国主要枸杞种植地[16]，研究区概况如图1所示。宁夏枸杞种植主要分布在中南部地区，北部地区种植面积较少。此外，枸杞每年的休眠期从当年11月至次年4月，生长期为7个月，即萌芽期为4月上旬，开花初期为5月上旬，果熟期为6月中旬，连续开花，连续结果，直至10月下旬落叶。

图1 研究区概况图Fig.1 Overview of study area

2 研究方法

2.1 数据准备及预处理

2.1.1数据准备

本文采用的实验数据由宁夏回族自治区遥感影像、枸杞种植区域矢量图和年际枸杞产量3类数据组成。其中，遥感影像采用Modis高光谱影像数据，来源于EARTHDATA网站，行列号为h26v04，其波段信息丰富，光谱分辨率可达纳米级，因而可提取农作物的反射峰、吸收谷等特征，更好地刻画农作物长势及信息差异；同时，Modis遥感数据具有高时间分辨率特性，可获得充足的多时相估产数据。根据枸杞每年4月上旬萌芽，10月下旬落叶，选择每年第97天至第297天，生长季总计201 d，26幅时相影像(时间分辨率为8 d)来构建CNN-S-GPR网络模型训练样本。每个时相影像包括MOD09A1、MOD13A1、MYD11A2、MCD15A2H 4类Modis产品数据，共13个波段，影像及矢量数据如表1所示。MOD09A1为地表反射率数据，其1～7波段可反映作物生长环境及土壤含水率；MOD13A1为植被指数产品数据，其中NDVI和EVI 2个波段可用于作物产量预测[17-18]；MYD11A2为地表温度数据，包含白天和夜间地表温度波段，与作物冠层温度密切相关[19]；MCD15A2H为叶面积指数和光合有效辐射产品数据，LAI和FPAR(光合有效辐射)是大量作物生长模型的基础，与作物产量之间的关系更为直接[20]。此外，枸杞种植区域矢量数据来源于宁夏农林科学院，为shp文件，由工作人员于2017年实地考察记录枸杞种植区并绘制而成；年际枸杞产量数据由宁夏回族自治区统计局提供，包括2010—2019年宁夏22个县(县级市区)枸杞种植面积和实际产量，因统计数据的不连续性，筛除西吉、隆德等6个县后，保留了同心、中宁等16个县作为估产区域(图1)，2019年产量数据如表2所示，研究区总面积为27 960 hm2，总产量为94 843 t，平均产量为3.39 t/hm2。

表1 影像及矢量数据Tab.1 Image and vector data

表2 2019年产量数据Tab.2 Yield data of 2019

2.1.2数据预处理

数据预处理过程主要分为4步：①利用MRT工具将Modis遥感影像和枸杞种植区域掩膜数据重投影为基于WGS-84椭球体的UTM投影，并将MYD11A2影像和枸杞矢量数据重采样的空间分辨率为500 m，保证数据空间位置的一致性。②由于MOD13A1时间分辨率为16 d，其余Modis数据产品为8 d，为了保证时序的完整性，采用上下影像求算术平均值的方法对影像进行补充，如将MOD13A1第97天和第113天影像的均值对第105天进行补充。③根据宁夏行政区划图分别将枸杞种植区域掩膜和Modis遥感影像数据裁剪为16个县枸杞种植区域掩膜及216 320幅遥感影像，并提取各县中心经纬度。④利用不同县域枸杞种植矢量图提取Modis高光谱遥感数据的感兴趣区域，以中宁县枸杞种植感兴趣区域提取为例，图2为中宁县枸杞种植区域矢量图，绿色区域为枸杞种植地块。图3为经步骤③处理所得的中宁县Modis高光谱影像MOD09A1产品数据，调用GDAL(Geospatial data abstraction library)库中的warp函数以实现图2对图3的裁剪，从而得到中宁县枸杞种植感兴趣区域图(图4)，并对同一天不同产品进行波段提取和融合处理后得到16 640幅包含13个波段的影像。

图2 中宁县枸杞种植区域矢量图Fig.2 Wolfberry planting area in Zhongning

图3 中宁县Modis高光谱影像MOD09A1图像Fig.3 Modis HSI MOD09A1 in Zhongning

图4 中宁县枸杞种植感兴趣区域图Fig.4 ROI for wolfberry planting in Zhongning

2.2 数据集构建

2.2.1直方图降维与归一化

鉴于训练样本的稀疏性，无法采用端到端方式训练估产模型。因此，本文将各县枸杞种植感兴趣区域的高维遥感影像映射到直方图中，使得产量与不同灰度的像素数量相关，而与影像像素位置无关。然而遥感影像中的像素值是离散的，假设每个波段最多可以取b个不同值，那么对于d个波段的影像将得到具有bd个像素区间的直方图。如在可见光图像中，每个波段可以取b=256个不同值，波段数d=3，则直方图像素区间个数约为1.67×107，显然将b个不同的取值离散至b个区间不可取，会导致维数灾难问题。因此，本文划定x(x

(1)

式中Hi——直方图归一化结果

hi、hj——像素直方图中每个区间像素数

2.2.2时间序列融合与维度转换

经上述直方图降维与归一化处理后，使得每个影像变成尺寸为32×13的矩阵。为了进一步融合高光谱影像不同时相上的枸杞生长信息，对上述直方图统计结果在时间序列上进行融合，从而得到时序样本(x轴为直方图区间维度，y轴为图像波段维度，z轴为时间维度)。然而，不同于3D CNN提取3维时序数据样本特征的方式，本文采用2D CNN获取3维时序数据样本的植被长势信息。为了提取时序样本中枸杞生长时间特征，需对图像波段维度和时间维度进行转换，最终形成32×26×13的矩阵作为网络的输入。将高、中、低产量不同波段的数据融合结果(x轴为像素区间，y轴为时间维度，z轴为影像波段维度)进行可视化对比，如图5所示。从图5可以看出，在高产、中产和低产中波段1、波段7、NDVI、EVI等波段明显存在视觉差异，表明特征抽取网络能够提取到有效特征。

图5 数据融合可视化结果Fig.5 Visualization results of data fusion

2.3 枸杞估产模型构建

2.3.1网络结构设计

高光谱影像波段宽而广，包含丰富的光谱信息和图像信息，能较好地反映作物长势。本研究首先通过构建CNN对影像特征进行提取，接着为了将注意力集中在对枸杞估产具有重要作用的通道上，在CNN后融合了通道注意力(SENet)模块，以建模通道间的依赖关系；最后引入高斯过程回归方法(GPR)整合影像特征和空间位置特征，对各县枸杞年际产量进行预测。综上所述，采用CNN作为骨干网络，通过引入SENet模块和GPR方法，提出一种基于CNN-S-GPR网络的高光谱影像估产模型用于枸杞估产，网络结构如图6所示。其中，第1层为输入层，为32×26×13的矩阵。第2层为特征提取层，由6个卷积层构成，其中淡青色为步长为1的卷积层，蓝色为步长为2的卷积层，边缘采用1个像素零填充，卷积核数量分别为128、256、256、512、512、512，卷积核大小全为3×3。在每个卷积层上进行批归一化处理后，采用ReLU函数激活，并加入Dropout层和L2正则化，以避免“梯度弥散”和模型过拟合问题[21]。第3层为SENet模块，对骨干网络最后一层进行Squeeze和Excitation 2个操作，以学习到不同通道的权重。最后为全连接层，将提取特征全连接为2 048维向量，并融合位置特征后作为高斯过程回归模型的输入。

图6 CNN-S-GPR网络结构图Fig.6 CNN-S-GPR network structure diagram

2.3.2高斯过程回归

基于预处理后的遥感影像数据集，首先采用CNN网络对2019年宁夏16个县的枸杞产量进行预测；然后利用半变异函数将各县距离与产量绝对误差相关联，将计算结果进行可视化，如图7所示。由图7可知，当空间上距离越接近时，半变异函数值越小。由此可见，枸杞种植空间位置对枸杞估产精度影响较大。若将遥感影像特征与枸杞种植空间位置结合起来，将对提升枸杞估产精度具有重要意义。

图7 半变异函数计算结果可视化Fig.7 Visualization of semivariogram results

因此，采用高斯过程回归预测枸杞产量，不仅整合了影像特征和空间位置特征，而且还解决了小样本模型训练问题[22]，其线性模型为

y(x)=h(x)Tβ+
(β～N(w,σwI),～N(0,σ2))

(2)

式中y(x)——观测值(产量预测值)

h(x)——SENet加权后的CNN特征向量

w——h(x)的权重向量

σw——超参数

I——主对角线为1的方阵

由于观测值y的绝对误差与x分布有关，即2个x距离越近，y的绝对误差越小。所以引入高斯过程f(x)替换了式(2)中的残差函数。不同于深度学习回归模型，该方法隐式地将输入特征映射到不同维度的空间，从而求解到模型最优权重，GPR方法直接对函数建模，可显式整合影像特征和空间位置特征，公式为

y(x)=h(x)Tβ+f(x)
(f(x)～N(0,k(x,x′)))

(3)

式中f(x)——核函数为k(x,x′)、零均值的高斯过程函数

核函数k(x,x′)采用径向基函数核，加入高斯噪声后，其公式为

(4)

式中gloc、g′loc——各县归一化后的中心经纬度

rloc、σ、σn——超参数

经上述推导后，整合影像特征和位置特征的GPR公式为

(5)

为了使高斯过程回归拟合效果达到最优，基于交叉验证，采用网格搜索进行最优超参数值的搜索，最终将σ设为1，σw为0.01，σn为0.32，rloc为0.5。

2.4 损失函数

在回归问题中，损失函数常被用于描述模型真实值与预测值的接近程度，常见的损失函数有L1损失和L2损失，具体公式为

(6)

(7)

式中n——样本总数

yi——第i个样本的真实值

L1loss——L1损失L2loss——L2损失

L1损失表示真实值和预测值之间绝对误差的平均值。然而，由于L1损失梯度不变，即使损失很小，梯度也非常大，不利于模型收敛。所以本文采用梯度变化的L2损失作为损失函数，以表示真实值和预测值之间的距离平方和。

2.5 评价指标

采用平均相对误差(Mean relative error, MRE)、均方根误差(Root mean square error, RMSE)和决定系数(Coefficient of determination,R2)作为评价指标，以验证基于CNN-S-GPR模型准确性。

3 实验结果与分析

3.1 模型训练

通过对宁夏16个县2010—2019年Modis遥感影像、枸杞种植区域矢量图和枸杞产量统计数据进行预处理、直方图降维、数据融合和维度转换后，共获得160个数据样本。实验将2010—2016年共112个样本作为训练集，2017—2018年32个样本作为验证集，2019年16个样本作为测试集，训练估产模型，训练集、验证集和测试集划分比例为7∶2∶1。

本文实验环境为Ubuntu 16.04.12操作系统，深度学习框架选用 Pytorch，软件环境为Python 3.6.10和Cuda 9.0，GPU 采 用 NVIDIA GeForce RTX 2080 Ti。模型采用分批训练，每个批次随机选择26个样本训练模型，共训练30 000次，并利用Adam优化器优化。学习率初始化为0.001，为防止模型因学习率较大导致跳过全局最优点的问题，采用分段常数衰减策略降低学习率，即当训练次数分别达到4 000和15 000时，学习率降为原来的10%。同时设置patient值，当验证集损失累计10次未发生改变时，提前终止训练。

训练效果如图8所示，模型逐渐收敛，直至迭代次数为20 000时，训练集和测试集损失达到最小并趋于稳定。保存模型并在测试集上进行测试，测试结果如图9所示，R2达到0.91。同时采用直线回归拟合真实值与预测值间的关系，直线斜率为0.81，接近于1。由此可见，经直方图降维和归一化处理后的高光谱影像信息丢失较少，且通过高斯过程回归将影像信息与位置信息整合后，能够准确拟合枸杞高产量和低产量样本。

图8 训练效果Fig.8 Training effect diagram

图9 测试结果Fig.9 Test result graph

3.2 实验对比

3.2.1消融实验

为了验证SENet模块和GPR方法对本文所提模型性能的提升效果，在自建数据集上，通过以下3个不同的模型进行消融实验：①用CNN(baseline)表示原卷积神经网络模型。②用CNN-SENet表示在模型①中融合通道注意力后的网络。③用CNN-S- GPR表示在模型②中添加GPR后的网络，实验结果如表3所示。

表3 消融实验结果Tab.3 Results of ablation experiments

从表3可以看出，各个改进后的模型在性能上均有一定程度的提升。其中，仅融合通道注意力后的网络(CNN-SENet)与基线网络模型相比，在MRE和RMSE上分别下降了0.26个百分点和15.39 t，且R2达到0.85，说明产量预测值中85%变异可由加权图像特征解释，通过通道注意力机制对CNN特征进行加权后，强调了有效信息，能更好地表征枸杞生长机理。融合了通道注意力和高斯过程回归后的网络(CNN-S- GPR)，针对影像特征及枸杞种植位置特征进行整合，可避免因缺少空间位置信息产生的估产误差，MRE、RMSE和R2分别达到13.57%、776.58 t和0.91，比其他2种方法效果更好。与CNN-SENet相比， MRE和RMSE分别下降了0.69个百分点和67.26 t，R2提升了0.06；与基线模型相比， MRE和RMSE分别下降了0.95个百分点和82.65 t，R2提升了0.08。由此验证了本文方法的有效性。

3.2.2对比实验

为了使实验结论更具有准确性与说服力，将本文模型与Lasso、岭回归(Ridge regression，RR)、SVR、DNN和 LSTM 5种近年来用于估产的主流回归模型进行对比。其中，Lasso、RR和SVR为机器学习回归方法；DNN和LSTM为深度学习回归方法，均有3个隐藏层和256个神经元。与高斯过程回归不同，对比的回归模型输入均采用手工特征，由高光谱影像感兴趣区域的平均NDVI、EVI、地表温度、LAI和FPAR值序列组成，共计16 640条数据。对比实验结果如表4所示。

由表4可知，本文方法在自建数据集上，MRE、RMSE和R2均优于近年来用于估产的主流回归方法。与表中对比机器学习回归方法不同，本文方法一方面运用时间序列融合方法融合了不同时相上的时间信息，另一方面采用CNN提取了NDVI、EVI等手工特征不具有的其他波段信息，使特征丰富度进一步增强，所以本文R2高于其他方法。虽然LSTM也能很好地处理时间序列，但是忽略了空间位置对枸杞生长的影响，所以GPR方法在整合影像特征和空间位置特征后，MRE和RMSE相对于LSTM下降了0.44个百分点和52.48 t。总之，与上述所提方法相比，本文方法提升了枸杞估产的准确性。

表4 对比实验结果Tab.4 Comparative experimental results

3.3 各县估产准确性分析

为了更进一步评估CNN-S-GPR模型在各县区枸杞估产的准确性，对2019年宁夏16个县(县级市区)枸杞统计产量与预测产量进行比较，各县估产准确性分析如图10所示。误差较大的区域主要为灵武市、盐池县、大武口区和贺兰县，误差产生的原因可能是没有足量的一年一度枸杞种植区域矢量图，枸杞种植面积的年际变化性使得模型预测产量数据偏差较大；此外，灵武市、盐池县和大武口区种植面积过小且存在混合像元现象，CNN所提取的枸杞生长信息过少，给模型准确预测带来了一定的困难。与之相反，模型在中宁县、同心县、沙坡头区等县表现优秀，统计产量曲线与预测产量曲线几乎重合。总之，CNN-S-GPR模型在研究区的估产精度较高，可实现宁夏枸杞年际估产。

图10 各县估产准确性分析Fig.10 Analysis chart of yield estimation accuracy by counties

4 结束语

提出的CNN-S-GPR高光谱影像年际作物估产模型，在CNN特征提取网络后融合通道注意力机制，表征不同通道的重要性，增强了枸杞生长机理性模拟，在此基础上，引入高斯过程回归方法，解决了小样本估产问题，同时整合影像特征和空间位置特征，使估产特征更加丰富。与其他估产模型相比，该模型平均相对误差和均方根误差下降了0.44～0.95个百分点和52.48～82.65 t，且R2达到0.91，验证了本文方法的有效性。