基于多时相无人机遥感植被指数的夏玉米产量估算

韩文霆 彭星硕 张立元,3 牛亚晓

(1.西北农林科技大学机械与电子工程学院， 陕西杨凌 712100；2.西北农林科技大学水土保持研究所， 陕西杨凌 712100；3.科罗拉多州立大学土木与环境工程系， 柯林斯堡 CO 80523)

0 引言

玉米是我国第一大粮食作物，种植面积最大可达3.811 93×107hm2，年产量2.25×108t，分别占粮食作物总面积和粮食总产量的33.6%和36.1%[1]。及时、准确预测玉米产量对农业经营管理、粮食政策制定、农业保险、农业用水效率评价具有重要意义，是发展精准农业的迫切需求[2-3]。

玉米产量的传统估算方法多采用地面采样估算，该方法费时、费力、具有破坏性，难以扩展到大尺度[4]。遥感技术在作物估产方面具有快速、无损、尺度大的优势，已成为精准农业的研究热点[5-6]。在作物遥感估产中，基于卫星遥感的作物产量估算具有传感器种类多、监测区域面积广等优势。文献[7]以美国玉米为研究对象，以各州为估产区，通过归一化植被指数(Normalized difference vegetation index，NDVI)与玉米单产间的最佳估产模型对2011年各州玉米单产进行了估算，并推算出全国玉米单产，结果表明，全国玉米单产的相对误差仅为2.12%；文献[8]利用重采样粒子滤波算法同化条件植被温度指数(Vegetation temperature condition index，VTCI)和叶面积指数(Leaf area index ,LAI)，并基于组合熵的方法构建加权VTCI 和LAI 与冬小麦单产的线性回归模型，同化VTCI和LAI的单产估测结果(R2=0.531)优于单独同化VTCI(R2=0.475)或LAI(R2=0.428)的估测结果。但是，对田块尺度的作物生长监测，在作物关键生育时期，要求遥感影像的采集频率为1～3 d/次[9]。卫星遥感存在时空分辨率低、易受气象条件限制等问题[10-11]，难以进行田块尺度的作物生长监测和产量估算。

无人机遥感技术具有时空分辨率高、成本小等优势[12]，提供了一种快速、准确获取田块尺度遥感影像的方法，可以对田块尺度的作物产量进行有效估算[13]。文献[14]建立了小麦单生育期植被指数和产量的线性模型，筛选小麦估产的最佳植被指数和最佳生育期，决定系数R2最高达到0.7。文献[15]结合产量与植被指数，建立了指数模型和Richards非线性回归模型，偏差信息准则(Deviance information criterion,DIC)最低为62.87。此外，还可将遥感数据与作物生长模型相结合[16]，或者结合多源的辅助信息，如作物氮含量、作物水分含量、辐射利用效率、作物高度、叶面积指数[17]、雷达数据、高重叠率可见光图像、倾斜成像构造的3D 点云等[18-20]。除结合辅助信息外，还可以综合作物不同时域的遥感影像。文献[21]建立了基于多时期水稻遥感数据的多元线性回归模型，其精度高于单一时期模型，决定系数R2最高达0.75。无人机遥感估产模型中加入辅助信息，融合多源和多时相的遥感数据，不仅提高了统计经验模型的可靠性与可拓展性[22-23]，也提升了作物长势诊断和估产的精度[24-25]。目前，尚未见针对玉米估产的最优植被指数和无人机遥感最佳监测时期方面的研究。

本文以内蒙古自治区规模化种植的夏玉米为研究对象，建立各生育期的6种典型植被指数和夏玉米产量的经验统计回归模型，确定适用于夏玉米估产的最优无人机遥感监测时期和最适植被指数。运用阈值滤波法，剔除影响遥感估产精度的土壤像元，评价土壤噪声对夏玉米无人机遥感估产的精度影响，以期探索简便、快速、灵活、准确的夏玉米遥感估产方法，建立适合研究区域以及周边地区的夏玉米遥感估产模型，为规模化农业经营管理提供决策辅助信息。

1 材料与方法

1.1 实验地概况

实验区域位于内蒙古自治区鄂尔多斯市达拉特旗昭君镇(40°26′0.29″N，109°36′25.99″E，海拔1 010 m)。属于温带大陆性气候，种植作物为夏玉米，属于春播中晚熟型，一年一熟。播种时间为2018 年5 月中旬，收获时间为9月。行距58 cm，株距23～30 cm，采用中心轴式喷灌机进行灌溉，实验小区布置如图1所示。实验地划分为5个扇形区域，每个扇形区域划分有3个样方(长×宽为6 m×6 m)[26]。

图1 无人机作业田块及采样监测点分布Fig.1 Distribution of UAV field and sampling points

1.2 实验设计及地面数据采集

本研究对5个扇形区域采用不同的灌溉条件，以田间持水率的不同百分比在营养生长期(Vegetative stages，V期)和生殖生长期(Reproductive stages，R期)进行不同梯度的水分处理，作物长势和产量存在较大的空间差异性，具有良好梯度差异。收获时，按图1a中标记的位置进行采样测产，采样面积为3 m2，将收获的玉米籽粒干燥8 h，得到玉米的干质量，测算玉米单产。灌溉和降雨量以及产量如表1所示。

表1 实验区降水和灌溉量Tab.1 Precipitation and irrigation in experimental area

1.3 无人机多光谱系统及遥感图像采集

多光谱图像采集采用六旋翼无人机搭载Red Edge五波段多光谱相机，如图2所示，具体波段如表2所示，相机焦距为5.5 mm，视场角为47.2°，图像分辨率为1 280像素×960像素。相机配备了光强传感器和2个3 m×3 m的灰板。光强传感器用于校正无人机航拍过程中外界光线变化对光谱影像造成的影响，灰板具有固定的反射率，可对航拍影像进行反射率的校正，从而生成反射率影像，进行植被指数的提取。

图2 无人机及相机Fig.2 UAV and camera

表2 波段中心波长和带宽
Tab.2 Wavelength and bandwidth

nm

实验时晴朗无云，平均气温28.6℃，相对湿度61.96%，平均风速1.12 m/s，微风。实验时间在11：30—14：30。多光谱无人机飞行高度为70 m，飞行方向为南北方向，航向、旁向重叠度分别为80%和70%。将每次实验按照固定航线拍摄的多幅图像，以日期为索引导入到瑞士Pix4D 公司的Pix4D mapper软件中，以实时动态测量(Real time kinematic，RTK)的方法获取地面相控点，导入小图对应的POS 数据，在软件中进行初始化处理，几何校正，构建三维模型，提取纹理以及构造地物特征，最终生成高清正射多光谱影像。

数据于2018年5月29日开始采集，2018年9月7日结束，共16次实验，期间覆盖夏玉米的V3(出苗期)、V6(拔节期)、VT(抽雄期)、R1(吐丝期)、R2(籽粒建成期)、R3(乳熟期)、R5(蜡熟期)6个关键生育期。由于拼接预处理后的原始图像包含除研究区域以外很大的区域，为了更加突显遥感影像的作物特征，需在ENVI 软件中裁剪处理。根据可见光影像中的样方对应裁剪出多光谱影像的15块实测区域，每一块实测区域对应实验样地的每一小样方，裁剪后取每一小样方的植被指数均值。

1.4 夏玉米估产方法

研究方法如图3所示，遥感影像拼接校正后，进行波段运算得到多种植被指数，再筛选最优植被指数；基于夏玉米的生长规律将夏玉米分为不同生育期，比较单生育期和多生育期的模型精度，筛选最佳监测生育期；最终采用阈值滤波法剔除土壤像元的干扰，得到无人机夏玉米最优遥感估产模型。

图3 研究方法流程图Fig.3 Flow chart of research method

1.4.1植被指数

结合无人机遥感平台所搭载的传感器，本文选取了6种常用的无人机遥感植被指数和作物的产量构建估产模型，即无蓝色波段增强型植被指数(EVI2)[27]、次生修正土壤调整植被指数(MSAVI2)[28]、优化土壤调节植被指数[29]、归一化植被指数[30]、土壤调整植被指数[31]、绿度归一化植被指数[32]。

1.4.2单生育期的玉米估产模型

为了能够迅速、准确地对玉米产量进行估测，本研究采用线性模型和最小二乘法构建玉米在不同生育期内的植被指数和作物产量之间的经验模型，公式为

Y=aVt+b

(1)

式中Y——产量V——植被指数

t——植被指数标号

a、b——待定系数

1.4.3多生育期的玉米估产模型

作物产量是一个有机物不断积累的结果，为了能够更好地体现相关机理，减少偶然误差，本文将玉米植被指数和玉米生育期进行牛顿-梯形积分，然后基于其积分值和作物产量构建模型

Y=cTg+d

(2)

其中

(3)

式中i——采集次序n——采集总次数

g——植被指数类型

ti——采集植被指数时距播种时间，d

Vi——ti时植被指数

Tg——植被指数基于全生育期的梯形积分

c、d——待定系数

1.4.4土壤像元背景干扰剔除

无人机遥感影像的精度较高，土壤像元在研究区域的比重较大。为了提取植被像元即分离土壤像元与植被像元，本研究采用最大类间方差法OTSU(Maximization of interclass variance)对拔节期的NDVI植被指数影像进行阈值分割。OTSU算法是通过统计整个图像的直方图特性来实现全局阈值的自动选取，计算得到分割阈值为0.4[33-34]。根据阈值通过掩膜处理得到对应的掩膜，并将掩膜应用于遥感背景中的土壤像元剔除，技术流程如图4所示。

图4 阈值滤波技术流程图Fig.4 Threshold filtering process

1.5 精度评价

选取15个采样点的数据对模型精度进行验证，采用两个指标对估产模型的精度进行评价，即决定系数R2和均方根误差RMSE。较高的决定系数代表较好的拟合程度，但不能完全体现精度；均方根误差用来衡量观测值同真值之间的偏差，越小的均方根误差代表越高的精度。

2 结果分析

2.1 基于夏玉米不同生育期和多生育期的估产模型评价

选取了玉米4个生育期(拔节期、抽雄期、吐丝期和蜡熟期)的植被指数进行单生育期作物产量建模，选取了多生育期的植被指数梯形积分进行多生育期作物产量建模。

图5为4个单生育期和多生育期的决定系数R2和均方根误差RMSE评价结果。图中S1为拔节期，S2为抽雄期，S3为吐丝期，S4为蜡熟期，S5为多生育期。从图5a可以看出，拔节期的R2都在0.63以上，最高达到0.72。单生育期的估产模型精度由高到低依次为抽雄期、吐丝期、蜡熟期、拔节期。多生育期的模型R2普遍高于单生育期的估产模型，最高能够达到0.80。图6为基于无人机遥感植被指数的夏玉米4个单生育期和多生育期的估算产量与实测产量的比较。图6显示无论是单生育期还是多生育期的估算产量，大部分都集中在6 000～8 000 kg/hm2，当实测产量低于7 200 kg/hm2时，模型均会出现高估的现象。当实测产量高于7 200 kg/hm2时，模型均会出现低估现象。综上所述，在4个单生育期内，抽雄期的植被指数作物估产模型决定系数最高，误差最小；多生育期的植被指数估产模型精度高于单生育期的模型。

图5 基于无人机遥感的4个单生育期和多生育期的植被指数估产模型评价Fig.5 Evaluation of yield estimation models based on vegetation index derived from UAV measurements for four single and whole growth periods

2.2 基于夏玉米不同植被指数的估产模型评价

图6 基于无人机遥感4个单生育期和多生育期的植被指数估算产量和实测产量比较Fig.6 Comparisons of maize yields as actually measured and estimated by models based on vegetation indices of four single and full growing periods from UAV remote sensing

图5和图6显示,基于不同植被指数的玉米估产模型精度不同。在拔节期，各植被指数的R2在0.1～0.25之间，精度较差，不具备任何实用意义；在抽雄期只有植被指数GNDVI和EVI2的R2高于0.7，其余植被指数的R2均高于0.63；吐丝期各个植被指数的模型精度仅次于抽雄期，模型精度由高到低依次为GNDVI、EVI2、MSAVI2、NDVI、OSAVI、SAVI；蜡熟期植被指数的模型精度较低，R2均在0.42以下。基于多生育期模型精度最高的植被指数为EVI2和GNDVI，R2大于0.78，其他植被指数的R2也均高于0.63。综上所述，单生育期精度最高的植被指数为GNDVI，R2为0.72，RMSE为487.20 kg/hm2；多生育期精度最高的植被指数为EVI2，R2为0.80，RMSE为413.29 kg/hm2，估产模型公式为

Y=14 054TGNDVI+4 632

(4)

Y=136TEVI2+2 527

(5)

2.3 过滤土壤像元后的估产模型评价

为了更好地提高夏玉米的无人机遥感估产模型，选取了4个单生育期和多生育期的遥感数据，进行阈值滤波法过滤土壤像元处理。表3、4显示，经过土壤像元滤波处理之后，拔节期和多生育期的R2提升效果显著，其中MSAVI2由0.17提升到了0.31；多生育期的模型精度最好，其中 GNDVI的R2最高为0.89，有3个植被指数的估产模型高于0.87，依次为GNDVI、MASVI2、EVI2；其他的单生育期的R2没有明显的差别。RMSE的变化和R2的变化类似，拔节期和多生育期的RMSE明显下降，其他单生育期的RMSE变化不明显。综上所述，经过阈值滤波法过滤土壤像元之后的估产模型，在拔节期和多生育期的精度提升最高，尤其是多生育期，精度最高的植被指数为GNDVI，决定系数R2为0.89，RMSE为299.35 kg/hm2，估产模型为

Y=192.5TGNDVI+5 127

(6)

基于该估产模型的产量示意图如图7所示。

3 讨论

3.1 不同生育期对夏玉米估产模型的影响

作物对不同光谱的反射率随着作物生育期的不同而变化，通过波段计算得到的植被指数可以很好地反映不同生育期内作物的特征，这就为通过植被指数进行作物估产奠定了理论基础。但是由于各个生育期内作物的生长状态不一样，基于各个生育期构建的估产模型的精度也就存在显著差异。玉米在拔节期主要为营养生长，该时期的生长特征并不能完全反映后期产量器官的有机物质积累，因此该生育期的模型精度较低。抽雄期是夏玉米由营养生长转为生殖生长的关键时期，玉米从抽雄前10 d到抽雄后25～30 d是玉米干物质积累最快、吸肥最多的阶段，从抽雄期到吐丝期期间，玉米的生长特征直接反映作物产量器官的干物质积累，在这两个生育期之间的夏玉米单生育期的估产模型精度最高。在夏玉米成熟期，玉米的营养物质向籽粒转移，叶片中的叶绿素含量下降，基于不同波段的植被指数和作物产量的相关性下降，所以在玉米的成熟后期估产模型精度下降。

表3 不同生育期和全生育期土壤像元过滤前后R2(n=15)对比Tab.3 Comparison of R2 (n=15) before and after soil pixel filtration in different growth stages and whole growth stages

表4 不同生育期和全生育期土壤像元过滤前后RMSE(n=15)对比Tab.4 Comparison of RMSE (n=15) before and after soil pixel filtration in different growth stages and whole growth stages kg/hm2

图7 产量示意图Fig.7 Output diagram

此外，本研究也针对多种植被指数构建了多元线性回归模型，但模型的误差很大，这可能是由于模型的过拟合以及参数的多元共性造成的。如果采用人工神经网络以及深度学习等算法，由于样本数过少，导致过拟合现象严重。同时，图6表明产量大部分都集中在6 000 kg/hm2到8 000 kg/hm2之间，当实测产量低于7 200 kg/hm2时，模型均会出现高估的现象。当实测产量高于7 200 kg/hm2时，模型均会出现低估现象。综上所述，本研究为了突出无人机快速、便捷的特点，折中选用了线性模型进行建模，在今后的研究中可以根据植被指数、纹理特征等参数结合实测产量构建非线性模型。

3.2 多生育期夏玉米估产模型精度

文献[35]指出，在对作物产量进行遥感估测时，遥感数据如果不是来自于决定作物产量的生育期，最终的估产模型是不准确的。作物产量是一个干物质累积的结果，单生育期的作物长势并不能很好地代表作物的产量。因此，本文考虑到作物的营养生长和生殖生长的过程，利用无人机遥感监测作物全生育期长势，基于夏玉米的生长期，对夏玉米的植被指数进行积分处理。为了凸显模型的实用性和无人机方便、快捷的特点，采用了牛顿-梯形积分和线性模型。通过对比发现，基于多生育期的作物估产模型精度高于任何一个单生育期的估产模型，基于多生育期的夏玉米估产模型不仅提高了模型的实用性，还提升了模型的机理性。在今后的工作中，可以将作物的营养生长和生殖生长分别处理，结合作物生长机理，构建精度更高、机理性更强的估产模型。

3.3 不同植被指数的夏玉米估产模型精度

基于植被指数对作物的相关特征进行反演，是一种常用的作物监测手段。植被指数多采用差值和比值的形式，可以抵消大部分单一波段的误差，获得更佳的植被信息。其中NDVI使用最广，植物叶片组织对蓝光和红光辐射有强烈的吸收，但是对于近红波段具有强烈的反射作用，NDVI采用非线性拉伸的方式增强了近红波段和红光波段反射率的对比度。在植被覆盖度比较低时能够准确指示作物的长势和生物量，但是当监测地区的植被量为高饱和状态时，NDVI的指示能力会下降。

农田尺度的夏玉米为植被覆盖度较高的区域。文献[32]指出，当NDVI对植被组成和光合太阳辐射的吸收率很敏感时，应添加一个绿色植被指数，以感知叶绿素的浓度，测量光合作用的速率以及监测植物的胁迫。文献[36]研究发现，GNDVI的估产效果要好于NDVI、TSAVI，尤其是抽雄期到灌浆期时GNDVI的估产效果最佳。文献[37]研究表明，在非干旱地区，EVI2 的估产效果优于NDVI，且对于产量差异幅度较大区域的反演效果更佳。SAVI通过添加土壤调整系数，减弱了土壤背景对NDVI的影响，OSAVI和MSAVI2是对SAVI的进一步优化，其中MSAVI广泛应用在作物的生长和产量的估测。文献[38]研究指出，MSAVI2能够更好地消除土壤背景对冠层的反射率影响。以上的研究成果多基于卫星遥感数据，本研究也发现GNDVI、EVI2、MSAVI2在玉米生育期的模型精度较高，决定系数R2在0.87～0.89之间，均方根误差在300～330 kg/hm2之间。针对不同的作物和不同的地区，还需进一步对相关的植被指数和遥感参数进行评选和优化。

3.4 土壤像元对夏玉米估产模型的精度影响

通过阈值滤波的方式可以有效地减少土壤像元噪声对估产模型精度的影响。表3和表4显示土壤像元噪声对夏玉米拔节期和全生育期的估产精度影响显著，对于抽雄期、吐丝期和蜡熟期的影响较小，这是由于玉米在抽雄灌浆之后叶面积达到峰值，冠层的郁闭度很高，土壤像元噪声很少。在今后的工作中，可以发挥无人机的高重叠率的特点，构建作物冠层的三维立体曲面，从而将玉米的冠层结构信息结合到玉米作物的估产模型中去。

4 结束语

通过实验和数据分析，利用无人机多光谱系统建立了夏玉米估产模型，并分析了土壤像元对不同生育期估产精度的影响，同时针对玉米单生育期和多生育期分别建立估产模型，并分析其精度。研究发现：不同生育期的玉米估产模型精度存在显著差异。单生育期中，精度由高到低依次为：抽雄期、吐丝期、蜡熟期、拔节期。多生育期的估产模型精度高于单生育期的估产模型精度。对遥感数据进行土壤像元滤波处理之后，拔节期和多生育期的模型精度得到了显著提升，而抽雄期、吐丝期和蜡熟期的模型精度并没有显著变化。基于土壤像元滤波之后的多生育期，运用植被指数GNDVI时，估产模型的结果最优，决定系数R2能够达到0.89。故采用无人机遥感植被指数进行作物估产，可以迅速、有效地对作物的长势和产量进行估测，为规模化农业经营提供了便捷，同时也为评价用水效率提供了一种估算指标。