基于多时相无人机遥感生育时期优选的冬小麦估产

王晶晶 李长硕 卓 越 檀海斌 侯永胜 严海军,4

(1.中国农业大学水利与土木工程学院， 北京 100083； 2.国家半干旱农业工程技术研究中心， 石家庄 050051； 3.中国农业机械化科学研究院集团有限公司土壤植物机器系统技术国家重点实验室， 北京 100083； 4.农业节水与水资源教育部工程研究中心， 北京 100083)

0 引言

小麦是中国三大粮食作物之一，种植面积达到2.338×107hm2，占全国粮食作物总面积的22.0%[1]。在小麦成熟收获前进行精准产量估算对保持粮食市场稳定、确保国家粮食安全具有重要意义。传统的粮食测产采用抽样方法，依赖大量的实地采样，成本高且费时费力，无法保证测产结果的时效性及准确性。卫星遥感可进行大面积的作物生长监测[2-4]，但是其估测精度易受云层、天气等影响，且遥感图像的时空分辨率均较低。

无人机遥感具有覆盖范围广、测量周期短、成本低和机动性强等特点，能够快速、准确获取作物表型关键信息，提供区域尺度的高时空分辨率图像。基于无人机遥感平台，利用解析图像获取植被指数，并结合机器学习算法，可对作物进行高效快捷的大面积产量估算[5-7]，以克服传统方法的不足[8]。文献[9-11]利用无人机搭载各种传感器，采用不同算法构建产量估算模型，相比于传统回归方法，机器学习算法提高了产量估算精度，但在采用神经网络等算法建模时需要大量的训练样本，且极易出现过拟合问题。偏最小二乘法(PLS)将预测变量和观测变量投影到一个新空间，建立线性回归模型，可以避免数据非正态分布、因子不确定性和模型不能识别等潜在问题[12]。随机森林(RF)是一种基于多决策树分类的新型机器学习算法，具有很好的抗噪性能和极强的拟合能力，且不会产生过拟合现象等[13]。套索算法(LASSO)是一种同时进行特征选择和正则化的回归分析方法，可增强统计模型的预测准确性和可解释性[14]。然而作物的光谱特征受生物因子和非生物因子的共同影响[15]，在不同生长阶段具有不同的冠层特点，其相应的敏感植被指数有所不同，产量估算模型精度亦有所不同。

为提高作物产量估算精度，相关学者开展了多时相无人机遥感的研究。ZHOU等[16]与韩文霆等[17]采用梯形公式对多时相植被指数积分，发现相比于单生育时期，采用多时相无人机遥感植被指数可提高水稻与夏玉米的产量估算精度。上述研究在建模时选择单一植被指数进行线性回归，然而不同生育时期的最优植被指数不同且存在饱和问题，同时输入参数与产量之间存在非线性关系，因此采用单个植被指数线性建模会降低模型的普适性和预测精度。程千等[18]采用冬小麦连续4个生育时期遥感数据进行产量估算，结果显示多生育时期的冬小麦产量估算模型精度均高于单生育时期，但其建模方法为多个生育时期植被指数的简单合并，并不能解释小麦全生育阶段的日植被指数的动态变化趋势。综上可见，多生育时期遥感可以提高农作物产量估算精度，然而随着时间的推移，多次收集信息会增加工作量和成本，因此评估最优生育时期及采集次数至关重要。

针对多生育时期遥感信息，为解决光谱在时间序列上不均匀问题，杨雪[19]与ZHANG等[20]使用Logistic回归构建了水稻冠层植被指数NDVI随时间的动态变化趋势，但有学者认为相比于Logistic和Gompertz，三次B样条曲线集成了三次样条和曲率最小化的优点，能够解释作物营养生长与生殖生长时期的变化趋势[21-22]，更适用于描述作物动态生长变化规律。然而不同作物的植被结构存在差异，会影响冠层光谱信息，采用样条曲线描述冬小麦逐日植被指数的规律以及对小麦产量估算精度的影响还有待研究。

为此，本文基于无人机遥感平台采集多时相冬小麦遥感影像，通过PLS、RF与LASSO构建单生育时期冬小麦产量估算模型，获得最优生育时期及算法；采用三次B样条曲线与牛顿复合梯形公式，利用最优算法构建5种特定生育阶段的产量估算方案，探究遥感采集次数对冬小麦产量估算精度的影响，以期优选出无人机遥感冬小麦产量估算的生育时期。

1 材料与方法

1.1 研究区域与试验设计

试验地点位于河北省涿州市中国农业大学教学实验场(北纬39°27′，东经115°21′)，海拔42 m，属于暖温带半湿润季风区，季节温差变化大且四季分明，多年平均降雨量563.3 mm，年平均气温11.6℃。试验地土壤类型(国际制)以砂土为主，0～80 cm土层的砂粒、粉粒和黏粒平均质量分数为87.92%、8.40%和3.68%，有机质、有效磷、速效钾、硝态氮、铵态氮含量(质量比)分别为11.72 g/kg、32.45 mg/kg、54.98 mg/kg、14.01 mg/kg、4.56 mg/kg。

研究作物冬小麦(品种为农大212)于2020年10月12日播种，2021年6月8日收获。播种量为270 kg/hm2，行距15 cm，播种时统一施入底肥，合计纯N为54 kg/hm2、P2O5为138 kg/hm2、K2O为81 kg/hm2。在冬小麦起身期至成熟期设置施氮处理和灌水处理两个试验。施氮处理采用人工撒施尿素(总氮含量大于等于46%)，设置5种追施氮肥水平，即N0(0)、N1(66.7 kg/hm2，1/3 N3)、N2(133.4 kg/hm2，2/3 N3)、N3(200 kg/hm2)、N4(266.7 kg/hm2，4/3 N3)。每个水平处理重复3次，共15个样区，各样区的总灌水量均为270 mm(含越冬灌水量30 mm)。灌水处理试验地块呈约60°扇形，划分成12个试验区，对应设置4个灌溉水平，即W1(150 mm，极度亏缺)、W2(190 mm，中度亏缺)、W3(230 mm，轻度亏缺)、W4(270 mm，充分灌溉)，均含越冬灌水量30 mm。灌水试验区采用水肥一体化追施尿素(总氮含量大于等于46%)溶液，追施氮量均为200 kg/hm2，每个试验区内划分3个样区，共36个样区。施氮和灌水处理各样区随机布置，面积均为6 m×6 m，统一采用一台3跨、半径为140 m的圆形喷灌机进行灌溉，其他田间管理相同，试验布置图如图1所示。

图1 研究区域和试验设计Fig.1 Experimental area and design

1.2 数据获取与处理

1.2.1多光谱图像

使用无人机(精灵Phantom 3 Advanced，深圳大疆创新科技有限公司)搭载多光谱相机(RedEdge，Micasense，美国)获取冬小麦多光谱图像。该相机由5个波段组成，波长范围分别为蓝(465～485 nm)、绿(550～570 nm)、红(663～673 nm)、红边(712～722 nm)和近红外(820～860 nm)。选择在11:00—13:00晴朗低风条件下飞行无人机，飞行高度为50 m，航向及旁向的重叠度均为75%。分别在起身期(4月6日)、拔节前期(4月14日)、拔节后期(4月22日)、孕穗期(4月30日)、抽穗期(5月7日)、开花期(5月11日)、灌浆前期(5月18日)、灌浆后期(5月27日)进行了8次影像采集。使用Pix4D Mapper 4.4软件对多光谱图像进行几何校正和辐射校正，利用ArcGIS 10.2软件对拼接预处理的遥感影像进行裁剪处理，得到每个试验样区光谱图像，提取各样区光谱反射率并求得植被指数。本文选取常见的15种无人机遥感植被指数构建冬小麦产量估算模型，相关植被指数的计算公式如表1所示。

表1 植被指数计算公式Tab.1 Equations for vegetation indices

1.2.2产量

冬小麦成熟时，在每个样区内选取3个代表性的1 m×1 m样方，对小麦脱粒后进行称量。通过谷物水分测定仪(LDS-1G型)测量冬小麦籽粒含水率，并将其换算为含水率13%的产量，3个样方产量平均值即为该样区产量实测值。施氮和灌水试验的产量数据如表2所示，51个样区冬小麦产量均值为5 675 kg/hm2，各处理间产量差异显著(P<0.01)。

表2 不同施氮和灌水处理的冬小麦产量Tab.2 Winter wheat yield of different fertilization and irrigation treatments

1.3 建模方法

首先将本文计算的15种植被指数与冬小麦产量采用SPSS 20.0进行相关性分析，然后选取显著相关(P<0.01)的植被指数，作为产量估算模型的输入参数，使用机器学习算法构建单生育时期以及多生育时期产量估算模型。从51组样本数据中随机选取35组作为模型训练集，剩余16组作为验证集，通过Python软件编写回归程序。

1.3.1单生育时期产量估算模型

为探究产量估算模型的最优生育时期和算法，分别采用PLS、RF、LASSO算法构建无人机遥感采集的冬小麦8个单生育时期产量估算模型。

1.3.2多生育时期产量估算模型

为分析多生育时期光谱采集次数对产量估算精度的影响，采用三次B样条拟合多时相植被指数，将采集时间点的植被指数数据作为节点，应用广义

交叉验证防止样条参数过拟合，寻找最优平滑水平，得到日植被指数，并通过复合梯形公式进行时间积分，复合梯形公式为

式中VI——植被指数积分值

f(t1)、f(ti)、f(tn)——采集时间为t1、ti、tn的植被指数

采用单生育时期产量估算模型的最优算法构建多生育时期产量估算模型，输入参数为各植被指数积分值，输出参数为产量。根据计算的日植被指数，定义5种特定生育阶段的日植被指数积分方案，其中方案1为采集的8个生育时期遥感数据，方案2至方案4依次向前推进一个生育时期，方案5为从8个生育时期中间隔挑选5个。其中4月6日(第1天)开始采集数据，5月27日(第52天)结束采集，具体方案如表3所示。

表3 5种多时相植被指数积分方案Tab.3 Five multi-temporal vegetation indices integration schemes

1.4 精度评价

植被指数拟合值与实测值的吻合程度采用相对误差(RE)进行评价，RE越接近0，表明拟合效果越好。采用决定系数(Coefficient of determination，R2)、平均绝对误差(Mean absolute error，MAE)、均方根误差(Root mean square error，RMSE)、标准均方根误差(Normalized root mean square error，NRMSE)来评价产量估算结果的准确性。R2越接近于1，MAE、RMSE、NRMSE越接近0，即表示模型预测精度越高；根据NRMSE标准分为4个等级，即极好(≤10%)、较好((10%，20%])、一般((20%，30%])、较差(>30%)[35]。

2 结果与分析

2.1 植被指数与产量相关性分析

冬小麦各生育时期的植被指数与产量相关性分析结果如表4所示。可以看出，大部分植被指数与产量的相关系数随着冬小麦的生长呈增大趋势，在灌浆后期达到最大值，其中NDVI、NDRE、OSAVI、GNDVI、PVI、GOSAVI的相关系数均达到0.90以上。在8个生育时期内，EVI、DVI、MCARI及EVI2与产量未全部达到显著性水平，因此选取RVI、NDVI、NDRE、OSAVI、GNDVI、CIRE、Vplot、PVI、GOSAVI、MSR和MSAVI共11种与产量显著相关(P<0.01)的植被指数作为冬小麦产量估算模型的输入参数。

表4 植被指数与产量相关系数Tab.4 Correlation coefficients between vegetation indices and yield

2.2 基于单生育时期的产量估算精度

采用PLS、RF和LASSO构建的冬小麦单生育时期产量估算模型训练集结果如表5所示。可以看出，在起身期时3种模型的R2较低，均未超过0.40，MAE、RMSE也均在1 000 kg/hm2以上，NRMSE均大于20%，产量估算精度较差，不具备实用意义。3种模型的产量估算精度在拔节前期与拔节后期均有所提升。在孕穗期和抽穗期，PLS和LASSO模型的精度接近，而RF模型的精度优于PLS和LASSO，但优势不明显。开花期之后，3种模型的精度进一步提升，其中PLS和RF模型均在灌浆前期取得最佳值，相应的R2分别为0.81、0.88，MAE分别为509、454 kg/hm2，RMSE分别为658、527 kg/hm2，但此时PLS模型的NRMSE为11.87%，而RF的NRMSE为9.51%，达到了极好水平。LASSO模型的精度在灌浆后期达到了最佳值，R2为0.81，MAE、RMSE分别为570、666 kg/hm2，NRMSE为12.01%。综上所述，单生育时期产量估算模型精度随着冬小麦的生长总体呈递增趋势；除拔节前期外，3种算法中RF的产量估算效果最好。

表5 单生育时期产量估算模型训练集结果Tab.5 Results of training set of yield prediction model for single growth stage

图2为单生育时期冬小麦产量估算模型的验证集结果，可以看出验证集的产量估算效果和训练集保持一致，即从起身期到灌浆后期，3种模型的预测精度逐渐提高，其中PLS和RF在灌浆前期预测效果最好，验证集R2分别为0.90、0.91，NRMSE 分别为10.26%、9.91%。而LASSO模型在灌浆后期预测效果最好，验证集R2为0.87，NRMSE为12.00%。比较各生育时期3种算法的产量估算精度，最终选择RF作为冬小麦单生育时期产量反演的最优算法，并用于构建多生育时期冬小麦产量估算模型。

图2 冬小麦单生育时期产量估算模型验证集精度Fig.2 Verification accuracy of yield of winter wheat for single growth stage

2.3 基于多生育时期的产量估算精度

2.3.1植被指数样条曲线拟合结果

选出与冬小麦产量显著相关的11种植被指数，采用三次B样条曲线进行拟合，可以得到植被指数在冬小麦起身期至灌浆后期的变化趋势。由于方案2、3、4植被指数变化规律与方案1一致，受篇幅所限，本节比较了方案1与方案5的拟合结果。图3仅列出有代表性的施氮试验N3处理和灌水试验W4处理拟合结果。可以看出，冬小麦起身期植被覆盖度较小，冠层光谱影像中包含大量的土壤背景信息，植被指数较小；拔节期至孕穗期冬小麦快速生长，植被指数呈快速增长趋势；孕穗期至开花期植被指数变化较稳定，并在抽穗期前后达到峰值；开花期至灌浆后期为小麦生殖生长阶段，植被指数呈现略微下降趋势。

图3 方案1与方案5中施氮处理N3、灌水处理W4植被指数样条曲线拟合结果Fig.3 Vegetation indices fitted by spline curve of nitrogen treatment N3 and irrigation treatment W4 in scheme 1 and scheme 5

图4为方案1与方案5中施氮处理N3、灌水处理W4的11种植被指数拟合值与实测值的相对误差。可以看出，方案1中N3处理的RE在拔节前期达到最大值9.36%，W4处理的RE在灌浆前期达到最大值11.52%；方案5中两种处理的RE均在拔节前期达到最大值，分别为10.41%、12.23%。两种方案不同生育时期的RE均小于15%，表明采用三次B样条曲线拟合计算的日植被指数，可用于构建多时相植被指数积分的冬小麦产量估算模型。

图4 方案1与方案5中施氮处理N3、灌水处理W4植被指数拟合值与实测值的相对误差Fig.4 Relative errors of vegetation indices of nitrogen treatment N3 and irrigation treatment W4 in scheme 1 and scheme 5

2.3.2积分预测模型精度

根据样条曲线拟合结果，将日植被指数按照设定的5种方案采用复合梯形公式进行积分，采用单生育时期优选的RF算法构建模型，模型的训练集与验证集的结果如图5所示。5种方案的预测值均与实测值呈良好线性相关，均匀分布在1∶1线两侧。相比于单生育时期，方案1、2、3、5的R2、MAE、RMSE、NRMSE有明显改善，均优于单生育时期产量估算模型精度。与起身期至灌浆后期8次遥感(方案1)相比，起身期至灌浆后期5次遥感(方案5)的训练集模型精度R2从0.96降为0.94，MAE、RMSE分别从212、299 kg/hm2增至318、378 kg/hm2，NRMSE由5.39%增至6.82%，可见方案1与方案5的模型预测效果均为极好，验证集各参数结果与训练集结果一致，表明适当减少采集次数并没有明显降低模型预测精度。

图5 5种积分方案产量预测结果Fig.5 Yield prediction results of five integral schemes

随着多光谱遥感采集生育时期前移，方案1、2、3、4的训练集模型精度呈逐渐下降趋势。方案1、2、3的NRMSE均低于10%，而方案4的NRMSE高达11.05%。表明采用起身期至开花期6次遥感(方案3)的产量估算模型精度达到极好的水平，可用于冬小麦产量的早期估算。

3 讨论

本研究优选的RVI、NDVI、NDRE、OSAVI、GNDVI、CIRE、Vplot、PVI、GOSAVI、MSR和MSAVI 11种植被指数与冬小麦产量呈现显著相关(P<0.01)，随着冬小麦的生长，大部分植被指数与产量的相关性呈增大趋势。采用PLS、RF和LASSO构建单生育时期植被指数的产量估算模型精度在冬小麦营养生长阶段呈快速递增趋势，在生殖生长阶段呈现缓慢递增趋势，这是由于随着冬小麦生长表型差异逐渐显著。3种模型产量估算的最优生育时期不同，PLS和RF模型为灌浆前期，LASSO模型为灌浆后期。这与文献[36-37]提出在小麦灌浆期进行产量估算有较高的预测精度相一致，因此灌浆期可能是冬小麦产量估算的最佳时期。本研究发现除拔节前期外，RF模型精度均高于PLS和LASSO，与文献[18]发现RF构建的冬小麦产量估算精度优于PLS和支持向量机(SVR)的结果相似，这是由于RF采用多个决策树进行样本训练，在处理复杂输入参数中优势明显[13]。

本研究采用RF构建的多生育时期冬小麦产量估算模型精度高于单生育时期，从起身期至灌浆后期8次遥感产量估算精度最高(R2=0.96，MAE为212 kg/hm2、RMSE为299 kg/hm2、NRMSE为5.39%)。与 8次遥感数据相比，从起身期至灌浆后期5次遥感数据的模型精度NRMSE为6.82%，表明减少遥感采集次数并未明显降低产量预测精度；而采用起身期至开花期6次遥感的模型精度NRMSE为9.16%，亦达到极好效果，这可以提前估测冬小麦产量。由于冬小麦不同生育时期的光谱反射特性会不同程度地受到外界环境影响，需要进一步有效融合多传感器数据或集成多种机器学习算法，以消除外界因素干扰，提高产量估算精度。

4 结论

(1)单生育时期产量估算模型精度随着冬小麦生育时期的发展总体呈现递增趋势。PLS、RF、LASSO模型的最优生育时期分别为灌浆前期、灌浆前期、灌浆后期。除拔节前期外，RF模型优于PLS和LASSO。

(2)采用三次B样条曲线拟合获得冬小麦起身期至灌浆后期的8次遥感(方案1)与5次遥感(方案5)的植被指数拟合值，与实测值的RE均小于15%，可用于构建植被指数积分的冬小麦产量估算模型。

(3)采用RF构建的多生育时期产量估算模型精度高于单生育时期，最优产量估算方案为起身期至灌浆后期8次遥感数据(方案1，NRMSE为5.39%)；与之相比，起身期至灌浆后期的5次遥感(方案5)减少了采集次数，可获得极好的产量预测精度(NRMSE为6.82%)；起身期至开花期的6次遥感(方案3)亦达到极好的预测效果(NRMSE为9.16%)。