基于无人机高清影像的棉花单产预测

祁佳峰，刘笑，杜文玲，郭鹏*

（1. 石河子大学理学院，新疆 石河子 832003；2. 绿洲城镇与山盆系统生态兵团重点实验室，新疆 石河子 832003）

新疆棉花种植面积占新疆耕地面积的50%左右[1]，占全国棉花种植总面积的78%左右，总产量占全国产量的80%以上。 棉花产量的准确预测对保证区域粮食安全，以及当地政府和农业部门制定相关政策等都具有指导意义[2]。

近年来遥感技术不断发展和进步，由于其获取信息速度快、观测范围广而广泛应用于农业生产的各个领域，而作物遥感估产也已逐步取代了传统的人工地面测量和统计的估产方式[3-4]。 目前遥感估产主要借助高空卫星和无人机两种工具进行，主要包括两种估产方式，一是根据地面统计数据与遥感指数之间的相互关系构建估产模型实现估产。 刘姣娣等[5]以Landsat7 和Spot4 为数据源，分析棉花不同生育时期的光谱特征后构建多个植被指数，通过与实际产量建立回归分析模型实现了对新疆生产建设兵团棉花产量的预测。 而程乙峰等[6]通过结合棉花花铃期的LAI 和归一化植被指数（normalized difference vegetation index,NDVI）与实际产量之间的关系，构建多元复合估产模型，实现了对新疆北疆地区棉花产量的预测。 卫星遥感棉花估产虽可以实现大尺度范围内的产量预测， 但估产精度有待进一步提高。为此，借助无人机图像空间分辨率高、获取速度快、影像成本低等优势，以无人机影像和光谱指数为基础进行了棉花估产研究[7]。白丽等[8]通过野外光谱仪获取了棉花冠层不同生育阶段的光谱曲线， 根据棉花冠层光谱特征构建光谱指数，建立了棉花产量的高光谱估算模型，实现了对棉花产量的预测。 二是将遥感反演的多种参数通过数据同化技术[2]融合到作物生长模型中实现产量模拟[9-10]。 戴建国等[11]利用高清无人机影像数据，通过采集棉花不同生育时期的时间序列数据，将不同时期的影像输入到卷积神经网络和双向长短期记忆的混合串行结构中实现了棉花产量估算。柏军华等[12]利用多时相遥感数据，通过叶面积指数（leaf area index, LAI）动态与棉花实际产量之间的关系， 构建叶面积指数-产量模型实现了对棉花产量的定量遥感监测。而在其他作物的估产研究中，何亚娟等[13]、徐春萌等[14]、王鹏新等[15]通过拟合地面统计数据与遥感数据， 构建估产模型，分别实现了对甘蔗、大豆和玉米的产量预测。Perros 等[16]利用无人机获取了2 个水稻品种4 个时期的无人机影像，基于植被指数建立了水稻估产模型。 Hassanzadeh 等[17]通过低空无人机系统，利用连续2 年的夏季遥感数据加地面点实测数据，使用2 种不同的植被指数检测算法，以偏最小二乘法建立估产回归模型，实现了对蚕豆产量的估算。

3～4 叶期之前，棉花植株较小，不同棉花之间的形态差异较小，无法进行有效分级；不同棉株的生长状态有较大差异；在此之后，棉花生长形成对地表的全覆盖，无法进行单株棉花的提取和分级。 已有研究根据作物幼苗的直径、株高和根系数量等特征指标实现对作物等级的评定和划分，从而实现对作物后期成活率、长势和产量等的预测[18-21]。 经试验发现，基于生长前期长势归类的同等级棉花在同等生长条件（土质、土壤肥力、水分、土壤盐分和土壤酸碱度）下，其生长状态在后期生长过程中不会发生较大改变，并且结铃期不同生长状态棉花的结铃数具有较大差异。以往的作物遥感估产研究多以拟合地面实测产量与遥感影像的方式构建估产模型，而且多以NDVI（normalized difference vegetation index, NDVI）和LAI 为关键指标实现估产，也有少量研究通过时间序列影像结合深度学习技术实现估产，但以棉花出苗和生长状态进行估产的研究相对较少。 鉴于此，本文在前人研究的基础上，借助无人机高空间分辨率的优势，根据棉花生长特性，根据3～4 叶期棉苗直径对棉花进行等级划分， 并结合后期统计的单株结铃数和铃重实现对田块尺度上棉花产量的预测。

1 材料与方法

1.1 试验设计

1.1.1试验区。 试验区位于新疆生产建设兵团第八师（后文简称“八师”）石河子垦区的石河子总场三分场二连（85°58′30″E，44°22′41″N），面积为26 779 m2。 该地区地势平坦，平均海拔450 m，地势东南高、西北低，冬季长而严寒，夏季短而炎热， 是典型的温带大陆性气候区， 年平均气温6.5～7.2 ℃，年降水量125.0～207.7 mm，年日照时间为2 721～2 818 h[22]。 该地区热量充足，光照丰富，非常适合棉花生长，是我国棉花的主要产区之一[23]。

1.1.2试验流程。 整个试验主要分5 步进行，依次为：第一，计算各类植被指数，选取具有最佳图像分割效果的指数， 并确定相应的分割阈值；第二， 根据棉花最大和最小直径实现对3～4 叶期棉花图像的分割，棉花幼苗提取和计数；第三，根据直径对棉花等级进行评价和划分；第四，根据实地调查数据， 对每株棉花的单株结铃数进行统计分析和等级划分，并对铃重进行统计；第五，构建估产模型，实现估产。

1.2 数据获取

1.2.1无人机影像。 试验使用的是深圳市大疆创新科技有限公司生产的大疆御2 四旋翼无人机，其整体重量约1 kg， 综合续航时间约20 min，本次试验获取的数据包含红、绿、蓝和近红外四个通道。在本次试验中，无人机飞行高度为30 m，航向重叠率为60%，旁向重叠率为80%。 影像拍摄时镜头垂直向下、悬停拍摄[1]。

影像拍摄时间为2021 年5 月19 日13：00，拍摄时天气晴朗、光线充足、无风，具备良好的飞行和拍摄条件。 最终获取的影像空间分辨率为0.007 5 m，拍摄的单幅照片均以JPEG 格式储存，拍摄后的图片由Pix4Dmapper 软件进行拼接和正射校正， 最终以WGS_1984_UTM_zone_45N为投影方式，并以TIFF 格式储存。

1.2.2地面点数据。 试验区的棉花品种为新陆早64 号， 是八师种植面积最广的棉花品种之一，生育期一般为123 d，属典型的早熟棉，该品种茎秆较硬，抗倒伏，适宜机采。试验区棉花于2021 年4月18 日播种，5 月初开始出苗， 至8 月中旬逐渐开始吐絮。

为建模和验证模型精度，对试验区实际出苗数、单株结铃数、铃重和实际产量数据进行采集。出苗数据采集时间为2021 年5 月19 日10：00－14：00，单株结铃数及铃重数据于2021 年10 月8日采集，实际产量数据在棉花采摘之后，咨询农户获取。

试验区棉花以宽窄行行距（66＋10）cm，1 膜6 行的方式种植，地膜宽1.82 m，地膜间距30 cm。利用GPS 进行定位， 共在试验区设置20 个1.82 m×1.82 m 的地面数据采样区统计单株结铃数， 同时又随机设置了18 个1.82 m×1.82 m的采样区进行估产模型精度的检验。 为保证试验的准确性并验证不同尺度上的精度，根据地膜总宽和膜间距， 分别设置3.94 m×3.94 m，6.06 m×6.06 m，8.18 m ×8.18 m，10.3 m ×10.3 m，12.42 m×12.42 m 的地面采样点各3 个，进行棉花提取精度和出苗率的验证。 同理，为确定铃重，在整个试验区随机采集200 个成熟棉铃进行称量统计。 试验区位置及采样区的分布如图1 所示。

图1 试验区示意图Fig. 1 Schematic diagram of the test area

为进一步验证最终的估产精度，设置由不同农户种植的与试验区相邻的地块为验证地块，总面积约2 500 m2， 棉花品种和种植密度均与试验地块一致。 验证地块的实际产量同样在采摘之后咨询农户获取；棉花铃重、不同等级棉株对应的单株结铃数等参数均采用与试验地块相同的设置。

1.3 数据处理与分析方法

1.3.1指数筛选。 本次试验中使用的四通道数据可以利用Envi 5.6 软件中的Spectral Indices 工具构建34 种植被指数。 为确定最佳的分割指数，分别利用34 种指数进行整幅图像的分割试验，其中，15 种指数无法对目标和背景进行分割，19 种指数可以对目标和背景进行分割，但部分指数在分割后产生大量噪点，试验所选取的指数及计算公式如表1 所示。

表1 植被指数列表Table 1 List of vegetation indices

1.3.2杂草去除方法。 经实地调查，杂草大多数生长在膜之间的未铺膜区，由于地膜的覆盖和遮挡，生长在地膜覆盖区的杂草极少，可忽略不计。本试验根据杂草与棉花的相对位置关系，以距离法对杂草进行去除，即在第1 行左侧和第6 行右侧10 cm 处划线，将10 cm 线之外的提取物确定为杂草，并将其去除。 杂草去除原理如图2 所示。

图2 去除杂草示意图Fig. 2 Schematic diagram of weed removal principle

1.3.3棉花提取与分级标准。3～4 叶期的棉花在形态方面已有较大差异，利用Envi 5.6 精准农业工具包 （Envi Crop Science） 中的作物计数工具（count crops） 对图像中棉株进行提取和计数，用测量工具（mensuration）测定棉花的直径。 首先均匀地在整个试验区内选择100 株较小的棉株测量其直径， 取平均值作为该棉花群体的最小直径，以同样的方式测量该群体的最大直径，最终确定该试验区棉株直径范围为4～14 cm。 以0.5 cm 为增量，从4～14 cm 共划分出21 个等级，最后以作物计数工具实现不同等级所有目标的提取。

然后，以试验区域内幼苗的最大、最小直径为限， 利用ArcGis 工具将棉花幼苗按直径大小以自然断点法分为4 个等级，对应为甲等苗直径范围为11.5～14.0 cm， 乙等苗直径范围为9.0～11.5 cm，丙等苗直径范围为6.5～9.0 cm，丁等苗直径范围为4.0～6.5 cm。

采用SPSS11.5软件包用于统计分析,等级资料的比较采用Ridit检验,对分级数据进行比较,P<0.05为差异有统计学意义。

1.3.4精度检验。采用统计学中常用的决定系数（coefficient of determination,R2）、 均 方 根 误 差（root mean squared error,RMSE） 来衡量模型精度。 其计算公式为：

式中：n代表样本数，yi代表预测产量，xi代表实际产量，y軈i代表预测产量平均值，x軈i代表实际产量平均值。

对于棉花幼苗提取结果精度的检验，分别在不同的尺度上验证其准确率（accuracy,A），其计算公式为：

式中：A代表棉苗提取的精度，T表示提取的棉苗数量，N表示遗漏未提取的棉苗数量， 本试验中未出现棉苗错提的情况。

2 结果与分析

2.1 棉苗提取

2.1.1影像分割。 棉花幼苗、地膜及裸土在形态和颜色方面有巨大差异，在遥感影像中可以根据绿色植物典型的光谱反射特性，利用植被指数实现对棉花植株的提取。 因此，试验以高清无人机影像为基础，利用各指数对影像进行分割，从而实现对棉花幼苗的提取，局部分割效果如图3 所示。 大部分指数在进行图像分割时会产生大量的噪点， 并且易对棉花幼苗产生的阴影进行误提，从而在进行图像分割时对分割效果产生巨大影响， 进而对之后棉花幼苗等级的划分造成影响，因此不适宜作为分割指数。 在选取的众多植被指数中，仅有绿叶指数（Green leaf index,GLI）和可见大气阻抗菌被指数（Visible atmospherically resistant index, VARI）可以将目标和背景进行有效区分，并且产生的噪点较少，分割效果较好。 但VARI 的分割阈值难以确定， 当其阈值过小时不能对棉花叶片的整个外在轮廓进行有效提取，而当其阈值调整至可以提取整个叶片轮廓时，会产生一定的噪点。 GLI 对整幅影像分割的总体效果优于VARI，在细节的提取和分割上效果更佳，分割时产生的噪点极少，图像分割也更加完整。

图3 不同植被指数分割效果图Fig. 3 Images of segmentation effect of different vegetation indices

综合比较各植被指数提取分割效果， 选取GLI 作为本试验中提取棉花幼苗的最终分割指数，经反复尝试和调整，当分割阈值为0.041 时，影像分割效果较好， 可以有效区分背景和幼苗。整个试验区域中共提取389 126 个目标。

2.1.2幼苗提取和出苗率精度检验。 为检验在不同尺度上对棉花幼苗的提取精度，在试验地块中随机设置不同面积的区域进行精度验证（图1），每个尺度上共进行3 次重复试验，不同尺度上的棉苗提取精度变化如表2 所示。 经人工实地计数与影像提取对比发现，当棉花幼苗叶片之间重叠部分过大时，存在错将2 株幼苗当作1 株进行提取的情况， 导致提取的棉花幼苗数量比实际数量偏少。但随着验证区域面积的增大，相对于整个验证样本的快速增加，漏提棉苗的数量增加相对缓慢，因此，随着面积的增大，棉苗提取精度逐渐提高。

表2 不同田块尺度棉花幼苗提取精度变化Table 2 Variation of extraction accuracy of cotton seedlings at different scales

2.2 棉苗提取与分级结果

图4 去除杂草前（A）后（B）效果图Fig. 4 Images before (A) and after (B) weed removal

不同生长状态的棉花幼苗分级结果如图5所示。 甲等苗为生长状态极佳、极为健硕的棉苗；乙等苗为生长状态略逊色于甲等苗、比甲等苗略小的棉苗； 丙等苗为生长状态略优于丁等苗，且比丁等苗略大的苗；丁等苗为发病、特别弱小的棉苗。

图5 棉花幼苗实景图（A）和示意图（B）Fig. 5 Actual view (A) and schematic diagram (B) of cotton seedlings

按照棉苗提取和分级结果，在成功提取并有效去除杂草后剩余的380 715 株棉花中， 甲等苗2 657 株，乙等苗103 753 株，丙等苗214 691 株，丁等苗59 614 株。

2.3 数据统计分析与模型构建

2.3.1数据统计与棉苗分级精度检验。 为确定模型构建所需参数，随机对20 个采样区内的1 483株棉花单株结铃数进行调查，经统计，在1 483株棉花中， 单株结铃数为4～13 的棉株数最多，占比超过96%； 而单株结铃数超过15 的共11株，单株结铃数大于13 小于15 的共21 株；而单株结铃数低于4 的共17 株。 按照采样区内棉花实际结铃数的统计结果， 以13 和4 分别作为单株结铃数的最大值和最小值， 按相等间隔划分，确定4 个等级棉花对应的结铃数，最终将甲等苗确定为13 个铃，乙等苗确定为10 个铃，丙等苗确定为7 个铃，丁等苗确定为4 个铃。 单株结铃数大于13 或小于4 的棉花不足5%， 没有代表性，忽略不计。

为验证棉苗分级的准确性，用20 个采样区内不同等级棉花的结铃数对分级准确性进行检验。 在20 个采样区的1 483 株棉花中，3～4 叶期的无人机图像中共划分出甲等苗75 株， 乙等苗415 株，丙等苗756 株，丁等苗237 株。 通过统计棉花结铃时的实际生长状况发现，27 株甲等苗生长状态欠佳，单株结铃数较少，与划分的等级有略微差异，占采样区内甲等苗的36%、占采样区内全部棉花的1.82%。 在27 株生长欠佳的棉花中，3 株的单株结铃数低于4，17 株的单株结铃数在4～7 之间， 其余7 株的单株结铃数在7～10之间。 此外，有23 株丁等苗生长情况有所改善，占采样区内丁等苗的9.70%， 占采样区内全部棉苗的1.55%。 在丁等棉花中，单株结铃数为4～10的共21 株， 剩余2 株的单株结铃数在10 以上。经检验，在1 483 株棉花中，共有50 株棉花的结铃数与苗期预测的结铃数有差异，生长状态发生改变的棉株仅占采样区中棉株总数的3.37%。 因此， 根据3～4 叶期棉花幼苗直径大小对棉花的分级取得了较好效果。

2.3.2估产模型构建。 经称量，200 个成熟棉铃中，铃重最大为7.8 g，最小为3.9 g，众数为6.4 g，平均为6.2 g，即0.006 2 kg。本试验以铃重平均值作为参数进行产量预测。 根据不同等级棉花幼苗数量、 单株结铃数以及铃重构建的NDCS（num ber and diameter of cotton seedlings）棉花产量预测模型定义如下：

式中：Y表示总产量，li表示划分的不同等级棉株的数量，ni表示不同等级单株结铃数，m表示铃重（kg）。

2.3.3模型检验和产量估计。 为验证模型的精度， 利用18 个采样区的实际产量对模型的验证结果如图6 所示。 决定系数R2越趋近于1，表示模型拟合程度越好， 而均方根误差RMSE越小，则说明模型越稳定，可信度也越高。 通过实际产量拟合，NDCS 模型的R2为0.919 2，采样区内产量的RMSE为0.168 7，表明模型精度较好，可以应用于棉花的估产研究。

图6 模型产量预测验证图Fig. 6 Validation graph of model yield prediction

在试验年间，试验地未遭遇极端天气和病虫害，试验材料生长状况良好，根据农户提供的产量数据可知，样地实际总产量为18 419.07 kg，结合样地区域面积，棉花单产为6 872.79 kg·hm－2。该试验通过对棉苗总数、不同等级棉花单株结铃数和铃重进行准确提取后，根据估产公式（5）对样地产量进行估计。 根据估产公式预测的样地区域总产量为17 442.86 kg，结合样地面积，预计试验区的棉花单产为6 508.53 kg·hm－2， 略低于当地实际产量。 该试验的预测产量是实际产量的94.70%，取得了较好的试验结果。

为了进一步验证模型在不同地块的估产效果，在试验区相邻地块（约2 500 m2）进行估产模型验证。 在验证地块中共提取甲等苗361 株，乙等苗10 919 株，丙等苗15 757 株，丁等苗2 985 株，经预测得到的总产量为1 463.96 kg，经咨询农户获取的该验证地块的实际总产量为1 581.29 kg，预测产量是实际产量的92.58%，同样取得了较好的估产效果。

3 讨论

3.1 基于无人机影像数据的棉苗提取与杂草去除方法

作物估产是现代精准农业中不可或缺的一部分，在众多估产方法中，多时相遥感估产受作物物候和大气效应的影响较大，从而易对估产精度造成影响[22,42-43]。 诸多作物估产研究中，利用棉花生长早期表现进行估产的研究还较少， 因此，本试验主要通过对棉花幼苗的提取和生长状态的分级并结合后期地面统计数据实现对田块尺度上的棉花产量预测。 本估产模型的关键在于棉花幼苗的准确提取和生长状态的合理划分。 对于作物幼苗的提取，闫春雨等[44]以面向对象的分类思想方法，充分利用影像的颜色和形状等纹理特征，以5 种不同的分类器实现了棉苗的提取并进行了精度对比分析，确定支持向量机为最佳的支持学习分类器。 魏文丽等[45]以地块的归一化植被指数均值和均方差反映了地块内玉米出苗的基本情况，最终精度达到80%以上。 而本试验直接利用最新的Envi 5.6 农业工具包实现棉苗的提取，经验证，利用GLI 指数提取棉苗，当验证区域面积为（12.42×12.42）m2时提取精度为96.2%。

试验区共存在4 种不同的地物， 分别为棉苗、杂草、地膜和土壤，而试验主要是将棉苗作为提取目标从4 种地物中进行分离。 地膜、土壤与棉苗之间具有明显特征差异，利用GLI 即可对背景和目标进行分割。 但3～4 叶期的棉花与杂草的外在形态和颜色极为相似，利用植被指数无法对二者进行有效区分，是影响产量预测准确度的重要因素，因此，有效去除杂草至关重要。 在地膜覆盖的棉田，杂草和棉苗最大的不同在于相对位置的差异，棉花幼苗以1 膜6 行的播种方式整齐排列，全部生长在地膜覆盖区；而杂草呈无规则分布，并且由于地膜的覆盖和遮挡，主要生长在膜与膜的间隙。 因此，可以利用距离法将提取到的地膜覆盖区以外的目标去除，从而达到去除杂草的目的。

3.2 植株的分级与应用

根据植株地径[21]、株高[46]、直径、株重和不定根数[47]等指标对苗木进行分级的研究已取得了一定研究成果，并已经广泛应用于多个物种中[48]。李亚麒等[49]根据云南松幼苗的地径和株高为分级指标，将幼苗划分为6 个不同的等级，最终实现了生物量的估算。关参政等[50]以枸杞种苗的根数、根径和地径为分级依据，对其后期的成活率和长势进行了预测，经验证取得了较好的试验结果。 结合前人的研究结果，本试验以3～4 叶期棉苗直径对棉株进行等级划分具有一定的合理性，并且根据后期棉株实际结铃数，在1 483 株棉花中，仅有3.37%的棉株与苗期划分的等级有差异， 说明棉苗的分级取得了较好的试验效果。

3.3 研究局限与展望

本试验充分利用和发挥了无人机影像高分辨率的优势，主要基于棉花出苗情况和苗期状态构建估产模型实现产量预测。 试验虽取得了较好效果，但存在一些局限和不足需要在后续的研究工作中进一步完善：（1）不同品种棉花之间的单株结铃数和铃重有一定差异，而本试验只针对一个棉花品种进行估产。 （2）不同的种植密度对棉花铃重和单株结铃数均会造成一定影响，但本研究仅在一个播种密度的棉田进行了试验。 （3）暴风、 暴雨和病虫害等自然灾害会对棉花长势、铃重和结铃数产生影响，而试验区在试验期间未遭受到极端天气和病虫害，因此逆境胁迫生长条件下对利用3～4 叶期棉苗长势的估产模型精度的影响有待进一步的评估。 （4）大面积尺度上，外界条件（土壤质地、土壤肥力、水分、土壤盐分和土壤酸碱度）有明显差异，这些差异会对棉花生长产生重要影响；而在田块尺度上，这些外界条件基本一致，因此，本试验更适宜田块尺度上的估产研究。

针对模型存在的局限与不足，综合考虑田块尺度上出苗整齐度和后期干扰因素，对不同品种类型、不同种植密度的棉田进行多年多点连续估产试验是以后的研究目标和重点。

4 结论

本文以苗期棉花高清影像结合地面统计数据构建棉花估产模型NDCS，实现了对棉花产量的预测。试验结果主要包括：（1）试验共提取3～4叶期棉花幼苗380 715 株，经不同尺度重复验证，棉花幼苗的平均提取精度达96.2%， 共划分出甲等苗2 657 株，乙等苗103 753 株，丙等苗214 691株，丁等苗59 614 株。（2）以不同等级棉苗数量及相应的单株结铃数和铃重构建的模型取得了较好的估产效果，经验证，估产模型的决定系数R2为0.919 2， 采样区内产量的均方根误差RMSE为0.168 7，模型拟合效果较好。 利用估产模型预测的试验区产量为6 508.53 kg·hm－2， 实际产量为6 872.79 kg·hm－2，预测产量的精度为94.70%；模型验证地块的预测产量精度达到92.58%，也取得了较好的估产效果。