基于无人机RGB图像的棉花产量估算

张静，郭思梦，韩迎春，雷亚平，邢芳芳，杜文丽，李亚兵，*，冯璐，*

（1.郑州大学农学院，棉花生物学国家重点实验室郑州大学研究基地，郑州 450000；2.中国农业科学院棉花研究所，棉花生物学国家重点实验室，河南 安阳 455000）

棉花是我国最重要的经济作物之一，棉花产量的快速准确估算对其精准管理和生产决策十分重要。传统的作物产量估算主要依靠田间测产调查，不仅费时费工，而且准确性较低[1]。近年来，遥感技术在农业生产监测中应用越来越广泛。卫星遥感技术由于易受天气影响不能满足田块尺度监测的需要[2-3]，而低空无人机遥感技术具有便于携带、低成本、高分辨等优点[1,4]，弥补了这一缺点，已成为精准农业研究的热门方向。

目前，国内外学者基于无人机搭载不同图像采集器（可见光相机、热红外相机、多光谱相机）对作物表型信息开展了大量研究[5-6]，实现了对作物生长参数，如株高[7]、生物量[8]、叶面积指数[9-10]等的有效评估。近年来，作物产量估算成为无人机生产监测的重要研究内容。王来刚等[11]基于多光谱图像提取植被指数建立玉米估产模型，发现与单生育期相比，多生育期估算效果更好。申洋洋等[12]基于无人机的多光谱图像利用结合机器学习算法提取植被指数，建立了小麦产量估算模型，在抽穗期估算效果最好。此外，基于可见光相机采集的RGB图像提取颜色、纹理等[13]特征，在小麦产量估测[14-15]中也得到了较好的验证，且多特征融合比单一特征估算精度高；多类型图像特征融合（如多光谱和热红外）较单一类型图像特征估算精度高[16]。因此，无人机遥感技术在作物产量估算中具有良好的应用前景。然而，现有的图像获取大多是在同一飞行高度进行，研究发现3 m高度下的水稻估算产量模型优于6和9 m[17]，但是不同高度的比较研究仍然较少，此外，基于无人机估算作物产量的研究多集中于小麦、玉米等粮食作物，对棉花的研究相对较少，且已有关于棉花产量估算方面的研究多是基于多年气象和棉花物候信息对全省范围的产量建立模型[18-19]，或者是通过光谱图像对中型田块进行棉花产量估算[1,20]，少见基于无人机RGB图像对田块尺度的棉花产量估算。

因此，本研究以棉花为研究对象，通过无人机搭载可见光相机采集不同生育时期和不同高度的RGB图像，为解决单一因素的预测精度不高的问题[4,6]，提取颜色指数特征和纹理特征并将其结合，通过逐步回归分析和因子分析2种方法处理数据，建立多元回归模型估算棉花产量，对比分析不同生育时期和不同高度建立的棉花产量估算模型，筛选出建模的最佳生育时期和最优高度。旨在为利用无人机快速估算棉花产量提供参考。

1 材料与方法

1.1 研究区域与试验设计

本试验于2020年7月至2020年9月在河南省安阳市白璧镇中国农业科学院棉花研究所东场试验基地（36°06′N,114°21′E）开展。供试棉花品种为鲁棉研28(SCRC28)。试验田共设24个小区，各小区面积为64 m2（8 m×8 m），种植行间距0.8 m，共10行。每组6个小区，均设置6个密度处理（1.5、3.3、5.1、6.9、8.7、10.5万株·hm-2），共4次重复（重复Ⅰ~Ⅳ）。棉花田间管理采用常规高产管理方式。

1.2 棉花产量测定

在收获期按小区实收测产，2020年分别在9月28日和11月2日收取2次棉絮，10月29日收获青铃，将其相加计算实收籽棉产量，公式如下。

式中，Y1为收取棉絮后计算的产量，kg·hm-2；Y2为收取的青铃产量，Y为总籽棉产量，W1为籽棉毛重，W2为口袋重，W3为小区青铃重，单位均为kg；S指小区面积，单位为hm2。

1.3 图像获取与处理

本研究采用大疆精灵四旋翼无人机平台，搭载2 048万像素CMOS传感器，FOV84°镜头，自动对焦。分别在棉花花铃期和吐絮期采集冠层图像，采集时间分别为7月10日至8月21日和8月25日至9月25日，各采集10次，图像采集当天时间为9∶00—10∶00，选择晴朗无风天气，阴雨天气顺延。飞行高度分别设置为20、30和40 m，每隔4 d采集1次，各小区在每个高度下各采集1幅图像，总计72幅图像。镜头垂直向下，获取图像像素大小为5 472×3 078，同步存储于Micro SD卡，以JEPG格式导入计算机中。

借助Photoshop CC2019和PyCharm 2020.2.3软件对RGB棉花冠层图像进行预处理和特征提取。预处理包括图像裁剪、分割、灰度化处理等。裁剪是按照小区划分裁剪成单一密度图像，并将裁剪后的图像保存为JPG格式；通过HSV阈值分割方法将棉株和背景（土壤、传感器等）分离；采用视觉法对分割后的图像进行灰度化进而提取纹理特征。经预处理后，棉花花铃期和吐絮期各飞行高度分别获得240幅冠层图像。

1.4 图像特征参数计算

1.4.1 颜色指数特征提取图像有多种颜色空间，无人机获取的原始图像是基于RGB彩色空间的[21]。本研究基于分割后的棉花冠层图像提取RGB颜色特征值，为了减少阴影的影响，首先对其进行归一化处理并计算出r、g、b值，然后根据不同组合计算12个可见光颜色指数（color index,CI）[22-26]，分别是可见光大气阻抗植被指数（visible atmospherically resistant index,VARI）、超绿植被指数（excess green vegetation index,EXG）、超红植被指数（excess red vegetation index,EXR）、超蓝植被指数（excess blue vegetation index,EXB）、绿红差值指数（excess green minus excess red vegetation index,EXGR）、归一化差分指数（normalized greenred difference index,NDI）、改良绿红植被指数（modified green red vegetation index,MGRVI）、Woebbecke指数（woebbecke index,WI）、川岛指数（kawashima index,IKAW）、绿叶指数（green leaf algorithm,GLI）、红绿蓝植被指数（red green blue vegetation index,RGBVI）、植物指数（vegetativen,VEG）。

式中，r、g、b分别为色度坐标，R、G、B分别为各分量的最大值。

1.4.2 纹理特征提取纹理反映了物体本身的属性，整体上表现出某种规律性，其灰度分布具有周期性[21]。数字图像中的纹理是相邻像素的空间位置变化的视觉表现[21]。本研究基于灰度共生矩阵（gray-level co-occurrence matrix,GLCM）[14,27]和Tamura纹理特征[28-29]提取灰度图像的纹理特征（textural feature,TF）。其中，灰度共生矩阵选择了4个常用参数，分别是能量（energy）、对比度（contrast）、熵（entropy）、逆差分矩（inverse different moment,IDM）；而Tamura纹理特征选择的3种常见指标是粗糙度（coarseness）、对比度（contrast）、方向度（directionality）。

1.5 产量估算

本文主要使用逐步回归分析（stepwise regression,SWR）和因子分析（factor analysis,FN）方法对RGB图像特征数据进行处理并建立多元回归模型，主要步骤为：①各生育时期采集的10次图像提取特征值计算平均值；②结合图像特征参数（12个颜色指数特征、7个纹理特征），利用逐步回归分析方法对所有特征值进行筛选，同时利用因子分析方法处理特征值，得到方差贡献率高于90%的因子；③分别探讨图像特征参数与产量建立的模型，确定估算棉花产量的最佳生育时期和最佳高度；④选择重复Ⅰ~Ⅲ为建模数据集，对建立的模型进行分析，用未参与建模的重复Ⅳ样本验证产量估算模型的精度。

逐步回归分析原理是依次不断反复引入自变量直到将所有显著变量引入，采用F检验和t检验进行显著性检验，最终将不显著自变量删除后得到最优解释变量集[30]；因子回归是利用降维方法将多变量归结为少数几个新变量，即因子的统计[31-32]。

1.6 评价指标

本研究采用决定系数（R2）、平均绝对百分比误差（mean absolute percentage error，MAPE）、均方根误差（root mean square error，RMSE）和标准均方根误差（normal root mean square error，NRMSE）4个指标来评价预测值与实测值拟合性。R2越大模型拟合效果越好；MAPE、RMSE越小模型精度越高；NRMSE越小，预测值和实测值一致性越好。

式中，xi为实测的棉花产量，为实测的棉花产量平均值；yi为预测的棉花产量；n为样本数。

2 结果与分析

2.1 特征值筛选和因子选择结果

将颜色指数特征和纹理特征结合起来，分别进行逐步回归分析和因子分析。逐步回归结果如表1所示，可以看出，WI、方向度和IKAW、方向度分别是高度20 m花铃期和吐絮期筛选后的特征值；WI、ENG和WI分别是高度30 m花铃期和吐絮期的影响因素；而高度40 m花铃期和吐絮期筛选的特征值分别是IDM和粗糙度。表2中FA1、FA2、FA3分别表示不同高度、不同生育时期经因子分析后贡献率超90%的3个因子，经计算得出各因子的系数值如表2所示。

表1 逐步回归分析结果Table 1 Results of stepwise regression analysis

表2 因子分析结果Table 2 Results of factor analysis

2.2 不同产量模型拟合效果分析

分别将逐步回归筛选出的特征值以及因子分析得出的因子与产量建立多元回归模型，对同一高度下花铃期和吐絮期的产量估算模型进行对比分析。由表3可知，在20 m高度下，花铃期的R2比吐絮期高，RMSE、NRMSE比吐絮期低；在30 m高度下，花铃期的R2均在0.7以上（=0.784 9，=0.757 4），而吐絮期的R2仅为0.415 9和0.505 2，RMSE、NRMSE均比花铃期低；在40 m高度下，花铃期和吐絮期的R2都在0.5左右，RMSE、NRMSE值较高，模型拟合效果一般。综合来看，20和30 m的高度下，花铃期产量估算模型拟合度均优于吐絮期且拟合效果较好，而40 m高度花铃期和吐絮期拟合效果一般。

表3 不同产量模型对比分析Table 3 Comparative analysis of different yield models

2.3 不同生育时期产量模型验证分析

2.3.1 高度20 m模型分析对建立的产量估算模型进行验证分析，结果表明，对于高度20 m（图1）建立的模型，在花铃期通过逐步回归方法建立的模型验证结果显著（P<0.05）而通过因子分析方法建立的估产模型不显著（P>0.05）；在吐絮期两种方法建立的估产模型验证结果都显著（P<0.05），其中，通过因子分析方法的R2更高、拟合效果更好，但通过逐步回归方法验证的R2与之仅相差约0.05，而MAPE、RMSE和NRMSE更小，其模型具有更高的精度，与实测产量的一致性更好。综合对比逐步回归方法在花铃期和吐絮期的拟合模型和验证结果，花铃期的产量估算模型拟合效果更好，模型精度更高。

图1 高度20 m花铃期和吐絮期产量估算模型验证结果Fig.1 Verification results of flowering and boll development stage and boll-opening stage yield estimation model at the height of 20 m

2.3.2 高度30 m模型分析 在30 m高度（图2）下观察到了类似的结果，无论是花铃期还是吐絮期，验证结果均显著（P<0.05）。所有验证的R2均达到了0.9以上，通过逐步回归验证的4个评价指标（R2、MAPE、RMSE和NRMSE）均优于因子分析方法的验证结果；对比通过逐步回归方法在花铃期和吐絮期的模型验证结果，花铃期的R2更高、拟合效果更好，RMSE和NRMSE更小，说明其模型具有更高的精度、与实测产量的一致性更好。因此，综合对比产量估算模型和验证结果，通过逐步回归方法在花铃期建立的产量估算模型拟合效果更好，模型精度更高。

图2 高度30 m花铃期和吐絮期产量估算模型验证结果Fig.2 Verification results of flowering and boll development stage and boll-opening stage yield estimation model at the height of 30 m

2.3.3 高度40 m模型分析40 m建立的产量估算模型验证结果如图3所示，可以看出，吐絮期验证的R2在0.8以上（=0.801 2，=0.946 9），但与拟合模型相差较大（=0.506 0，=0.448 1）（表3），尽管通过逐步回归方法验证的R2比因子分析方法低，但MAPE、RMSE和NRMSE均比因子分析小，因此，与通过因子分析方法预测的结果相比，逐步回归方法预测值更接近实测值。

图3 高度40 m吐絮期产量估算模型验证结果Fig.3 Verification results of boll-opening stage yield estimation model at the height of 40 m

综合来看，3个高度的2种方法对比结果均显示，通过逐步回归方法建立的产量估算模型拟合度和精度优于因子分析方法。在高度20 m花铃期建立的产量估算模型预测效果更好、精度更高；高度30 m在花铃期建立的产量估算模型预测效果和精度相当；而40 m高度下，与花铃期相比，在吐絮期建立的产量估算模型效果更好，但拟合效果一般，稳定性较差。

2.4 不同高度模型对比

进一步比较通过SWR方法在20（花铃期）、30（花铃期）和40 m（吐絮期）3个高度下建立的3个模型并筛选出产量估算模型的最佳高度和最佳生育时期。由表3可知，高度30 m在花铃期建立的模型R2（0.784 9）最高，20 m高度下在花铃期建立的模型的R2（0.733 3）次之，高度40 m在吐絮期的R2（0.506 0）最低，而RMSE和NRMSE的值从小到大依次是高度30、20、40 m。综合来看，在花铃期、高度为30 m建立的产量估算模型拟合度优于20和40 m。

分析3个拟合模型的验证结果可知（图1～3），以花铃期、高度30 m图像建立的模型验证结果R2（0.977 6）最高，MAPE（2.31%）、RMSE（93.7 kg·hm-2）和NRMSE（2.81%）最小，与拟合的指标相差最小。因此，无论是从模型拟合度还是验证集的评估指标来看，高度30 m下的模型拟合度最好、稳定性最强、精度最高。

综合不同生育时期、不同高度的分析结果，产量估算模型的最佳生育时期为花铃期，最优高度为30 m。

3 讨论

多项研究表明，作物产量和图像颜色特征、纹理特征存在相关性并能够拟合模型进而估测产量，如龚红菊[32]验证了基于纹理特征估算水稻产量的可行性。但是通过无人机获取图像时，光照条件不同可能导致基于单一特征指数的预测效果随机性较大，将颜色特征和纹理特征结合，其结果比单一颜色特征建立的小麦产量预测模型精度高[14]。本研究选择结合棉花冠层图像的颜色指数特征和纹理特征对产量进行估算并得到较稳定的模型，与万亮等[6]研究结果一致。

本研究探讨了不同棉花生育时期建立的产量估算模型，但在不同的生育时期预测效果不同，高度20和30 m下，花铃期建模效果均优于吐絮期，模型拟合程度更高、验证结果更显著。刘焕军等[33]基于时间序列Landsat影像建立棉花最优产量估算模型，发现其随着生长阶段表现出差异，得出花铃期更适合棉花产量估测的结果，与本研究结果一致。而在40 m高度下，吐絮期产量估算模型拟合度优于花铃期，但拟合模型的R2与验证数据集相差较大，尽管棉絮占产量的大部分，但吐絮期叶片大量脱落、棉株倒伏[34]等原因可能导致提取特征不精确，也可能与样本数量较少有关。本研究用一年数据进行建模，而田间试验受气温、降水等自然环境条件、管理措施等影响较大[35]，拟进一步增加样本数量以及用不同品种、不同年限的棉花数据进行分析，得到更普遍的产量估算模型。

本文探讨了不同高度的棉花产量模型对比分析，得到30 m的模型预测结果最好。一般来说，不同高度对试验地的图像分辨率不同，飞行高度越低精度越高[36]，但本文结果不同，可能与高度设置较小有关，有研究者在10[14]、30[37]、50[1]、80 m[11]都得到了较好的预测结果；也有可能与采集图像时的风力、光照变化[38]及预处理时使用的图像分割方法相关。因此，今后可以设置更多高度梯度以找出更合适的飞行高度进一步优化棉花产量估算模型，也可以基于聚类和深度学习[39]等方法对图像进行分割处理找出更合适的分割算法。本研究结果表明，通过提取图像特征建立产量估算模型是可行的，但仍需要进行更多的研究，例如气候、农艺性状等信息与图像结合进一步建立更复杂的模型。