基于无人机影像的SEGT 棉花估产模型构建

2020-10-21祁佳峰汪传建张国顺

农业工程学报 2020年16期

徐权，郭鹏※，祁佳峰，汪传建，张国顺

（1. 石河子大学理学院，石河子 832003；2. 兵团空间信息工程技术研究中心，兵团空间信息工程实验室，石河子 832003）

0 引言

作物估产对保障国家粮食安全、经济政策制定和农业可持续发展具有重要意义[1]。棉花是中国重要的经济作物，及时了解棉花的长势和产量，不仅对加强其生产管理和生产计划安排有利，也对棉花外贸和进出口计划制定有直接帮助，便于农业部门提前实施相应的管理与决策，获得更好的经济及环境效益[2-3]。

随着信息技术的不断发展，对作物产量的估算研究已经从传统的地面测量发展到了多维时空的遥感估算[4]。遥感技术由于具有宏观性强、获取信息快、感测范围广等特点，已被广泛应用于农业研究当中[5-8]。针对棉花产量的估算，国内外众多学者也展开了卓有成效的研究[9-13]。Haghverdi 等[14]利用 Landsat 8 遥感数据，通过人工神经网络和作物物候指标预测了美国田纳西州西部的棉花皮棉产量。 Dalezios 等 [15] 利用NOAA/AVHRR 卫星数据，建立了棉花生长季内时间序列的NDVI 模型以估算产量。庄丽等[16]通过分析棉花产量与冠层光谱指数之间的相关性，建立了基于高光谱成像的棉花估产模型。刘焕军等[17]利用时间序列的Landsat遥感影像，采用相关分析和线性回归方法对植被指数与实测产量数据的关系进行了分析，建立了棉花产量预测模型，对美国加州圣华金河谷种植区的棉花产量进行了预测。由此可见，利用卫星影像和时间序列数据对棉花产量估算是一种行之有效的方法，但在实际估产应用中，高空遥感影像受下垫面、大气效应、作物物候等影响较大，估产精度有待进一步提高[18]。而无人机遥感空间分辨率高、获取数据快、操作简单、成本低，能够快速针对某一区域进行影像采集，获取更精确的作物分布信息，成为航空遥感和卫星遥感的重要补充，对作物监测技术的发展和应用具有重大意义，众多学者也利用无人机开展了大量研究工作[19-21]。

但遥感估产方法易受大气效应和作物物候等影响，估产精度有待提高，结合棉花出苗和生长状态监测的估产研究相对较少，本研究利用无人机高分辨率遥感影像，提出一种基于出苗株数并结合生长过程状态的棉花估产理念和方法，构建了基于苗铃生长趋势的SEGT（Seedling Emergence and Growth Trend）估产模型，以期为棉花估产提供新的研究思路。

1 研究区概况

研究区位于新疆生产建设兵团第八师石河子垦区（86°01′00″～86°01′50″E，44°29′30″～44°30′00″N），如图 1 所示。第八师石河子垦区地处天山北麓中段，古尔班通古特大沙漠南缘，该地区地势平坦，平均海拔450.8 m，自东南向西北倾斜，属于典型的温带大陆性气候，冬季长而严寒，夏季短而炎热，年平均气温在6.5～7.2 ℃之间，气温北部低，南部高，年降水量在125.0～207.7 mm 之间，无霜期为168～171 d，日照充沛，年日照时数为2 721～2 818 h。该地区生态气候条件独特，耕地平整连片，条田建设规范，棉花种植的机械化和规模化水平较高，适宜实施遥感估产和精准农业。

图1 研究区示意图Fig.1 Sketch diagram of study area

2 研究方案

首先借助无人机可见光影像数据，通过计算植被指数与大津法、形态学滤波相结合的方法，获取研究区内每株棉花的空间位置；然后利用无人机多光谱时间序列影像数据，分析各时期归一化差异植被指数 NDVI（Normalized Difference Vegetation Index）与实际产量的相关特征，构建棉花生长状态图，并对生长状态进行等级划分，反演每株棉花的预测成铃数；最后结合棉花单铃质量构建 SEGT 模型进行棉花产量估算，并根据实测产量进行模型精度验证。研究方案如图2 所示。

图2 技术路线图Fig.2 Technology roadmap

2.1 无人机影像数据获取与处理

本研究分别采用无人机可见光和多光谱影像开展试验，其中可见光影像由大鹏CW-10 垂直起降固定翼无人机平台搭载佳能EF-M18-55 相机获得，佳能相机可获得RGB 三通道灰度图像，空间分辨率为2.5 cm；多光谱影像由eBee SQ 无人机搭载Parrot Sequoia 传感器获得，空间分辨率为 10 cm。棉花苗的提取对影像分辨率要求极高，因此试验中没有使用Parrot Sequoia 自带的可见光影像，而是选用空间分辨率更高的佳能相机。无人机可见光影像于2018 年5 月23 日上午11 点拍摄，多光谱影像分别于2018 年6 月至8 月拍摄，共10 期数据，数据采集时天气情况良好，无风，无人机飞行高度150 m，航向重叠度为80%，旁向重叠度为 60%。无人机影像数据通过 Pix4Dmapper 软件进行影像拼接和正射校正，处理后影像以TIFF 格式存储各通道的灰度信息，每个通道包含8 bit 信息，投影方式为WGS_1984_UTM_zone_45N。

2.2 地面实测数据获取

地面实测数据主要包括棉花的真实出苗信息和实际产量，用于模型的构建和结果的验证。试验在研究区内均匀划定了60 个3 m × 3 m 的样地，其中40 个试验样区，20 个验证样区（详细分布见图1）。为准确记录各个样区的位置，分别在每个样区中心插入一根花杆，并在花杆顶端放置一个红色圆盘，以便可以更准确地在影像中找到样区位置。试验于2018 年5 月23 日进行，利用手持GPS 采集并记录60 个样区内棉花苗的数量和空间位置，共计采集5 291 株棉花苗。将实测数据与提取结果进行对比，对棉花苗提取结果进行精度评价。试验选用的棉花品种为新陆早 64 号，是新疆兵团第八师种植面积最广的品种之一，该品种属于早熟陆地棉常规产品，生育期 123 d，株型为塔型，茎秆坚硬抗倒伏，生长稳健，叶片中等大小，叶上举，通透性好。铃卵圆形，较大，吐絮畅，宜机采，衣分 43.6%左右，纤维长度约30 mm，马克隆值4.2，整齐度85.5%。棉花产量在吐絮期已经基本形成，测产一般都在该时期内进行，试验于2018 年9 月20 日开展测产工作，利用天平称量并记录 60 个样区收获的籽棉质量，产量调查结果用于构建模型和验证估产精度。

2.3 棉花出苗位置与株数计算

鉴于棉花植株与裸土在影像中存在明显差异，本试验利用植被指数和大津法相结合的阈值检测方法确定棉花苗的空间位置，并利用数学形态学滤波进行降噪处理。

2.3.1 植被指数的筛选

植被指数可以根据绿色植物在各个波段不同的反射和吸收特性，对所获取的波段进行运算，以增强植被信息，其本质是在综合考虑各光谱信号的基础上，把多波段反射率进行数学变换，从而形成一些表征植被参数的数值[22]。相对于卫星遥感影像，搭载佳能数码相机获取的无人机影像由于缺少近红外波段，只能借助可见光波段构建植被指数来提取棉花信息。常用于表征植被信息的指数有过绿减过红指数（Excess Green-Excess Red）[23]，可见光波段差异植被指数（Visible-Band Difference Vegetation Index，VDVI）[24]，归一化绿红差异指数（Normalized Green-Red Difference Index，NGRDI）[25]，归一化绿蓝差异指数（Normalized Green-Blue Difference Index，NGBDI）[26]等，具体计算如表1 所示。

表1 常用可见光植被指数Table 1 Usual visible vegetation indices

2.3.2 灰度阈值法与降噪

大津法能使类间方差达到最大时的阈值自动确定为最佳阈值，本研究利用ENVI 软件中的Class Activation to Pixel ROI 工具，选择OTSU 算法，对植被指数灰度处理后的图像进行二值化处理，获取棉花的出苗信息。

图像滤波可以简化图像数据，保持图像基本形状特征，并除去不相干的结构，可进行噪声抑制、边缘检测、图像分割和形状识别等处理。通过ExG-ExR 指数和大津法相结合可以得到较好的棉花苗提取结果，但图像仍存在少量噪点，本研究采用数学形态学滤波来解决这个问题，利用ENVI 软件中的Convolutions and Morphology 工具进行卷积运算，卷积核大小为3（Kernel size=3×3，像素），滤波的重复次数为1（Cycles=1）。

2.3.3 精度评价

采用精确率（Precision）、召回率（Recall）、F1值3 个指标对植被指数的提取结果进行精度评价。

式中 TP 为正确提取的棉花苗数量，FP 为错误提取的棉花苗数量，FN 为未被提取的棉花苗数量。

2.4 棉花生长状态与结铃数等级划分

植被指数能够综合不同遥感波段下的植被反射特性，与作物生物量、叶面积指数等存在极强的相关关系，因此许多研究经常使用植被指数评价作物的生长状态[27]，其中NDVI（Normalized Difference Vegetation Index）是监测植被覆盖情况和生长状况的最佳遥感指数[28]。因此本文利用无人机多光谱时间序列影像，通过分析各期棉花 NDVI 与实际产量之间的相关特征，确定各期 NDVI在生长状态评级中的权重，构建能够整体反映棉花生长状态的综合植被指数CNDVI（Comprehensive Normalized Difference Vegetation Index）CNDVI 值越大，棉花生长状态越好，CNDVI 定义如下：

式中ai代表不同时期NDVI 的权重；j代表多光谱影像期数，j=10。根据每期NDVI 与实际产量的相关系数大小占相关系数总和的比例获得。

根据得到的CNDVI 指数，参考李敏等[29]研究成果，利用自然间断点分级法对棉花生长状态进行等级划分，并结合当地实际情况，将棉花分为 3 个等级，分别为一等苗、二等苗、三等苗。一等苗指发育正常且健壮的苗；二等苗指比正常苗偏小偏弱的苗；三等苗指病苗、小苗和弱苗。不同的生长状态等级，对应不同的棉铃数，也对应了不同的预估产量区间。针对分级后的CNDVI数据，于2019 年8 月29 日棉花吐絮期通过GPS 定位在研究区60 个样区内实地测算不同生长状态等级下的棉铃数，根据实际调查情况，CNDVI 值越大，棉花结铃数越多，结铃数大于12 个和小于6 个的植株数量相对较少，约占总数的 5%，不具有代表性，因此最终确定一等苗对应 12个棉铃，二等苗对应9 个棉铃，三等苗对应6 个棉铃。

3 估产模型构建与精度验证

3.1 估产模型构建

棉花产量估算模型定义如下：

式中SEGT 为区域内棉花单位面积估测产量，kg/hm2；n代表棉花生长状态等级，n=3；Ci代表不同生长状态等级下棉花苗数量；Qi代表对应等级下单株棉铃数量；W代表单个棉铃质量，kg；S代表区域面积，hm2。

试验用棉花品种的单铃实测质量约6.3 g，因此W取值为0.006 3 kg，棉花产量估算公式为

式中C1、C2、C3分别代表一等苗、二等苗、三等苗的数量。

3.2 估产模型验证

用实测产量对估产结果进行验证，选取统计学中常用的决定系数（R2）、均方根误差（RMSE）和相对误差（ME）3 个指标评价模型估产效果。

4 结果与分析

4.1 棉花出苗提取结果分析

棉花苗提取是构建SEGT 模型的一个关键环节，在很大程度上影响最终的估产结果，因此棉花苗的提取结果准确性至关重要。本试验利用无人机获取的高分辨率可见光影像数据，通过计算植被指数提取每株棉花苗的空间位置，统计区域内总株数，局部提取结果如图3 所示。最终，在研究区42.47 hm2面积内，共计提取出4 364 255 株棉花苗，其种植密度为102 761 株/hm2。将提取结果与实测数据（2.2 节地面实测棉花苗数据）对比可知，该方法提取的精度较高，精确率为93%，召回率为92.33%，F1值为92.66%。评价结果表明，该方法能够有效提取棉花出苗信息，可为后续棉花估产提供良好的基础。

ExG-ExR、VDVI、NGRDI 和 NGBDI 指数对棉花苗识别精度如表2 所示。由表2 可知，ExG-ExR 指数针对棉花苗的提取在精度和稳定性方面最好，其中 3 块样地的平均精确率、平均召回率、平均F1值分别达到了93%、92.33%和 92.66%，明显优于其他植被指数。VDVI 指数提取效果相对较好，但是在稳定性方面不如ExG-ExR 指数，NGRDI 与NGBDI 指数在精度和稳定性方面表现都相对较差。这些指数对植被识别都具有较好的效果，但对于单株棉花苗的提取识别，研究认为ExG-ExR 指数效果最佳。

图3 基于无人机影像的棉花苗提取结果Fig.3 Extraction results of cotton seedling based on UAV images

表2 基于不同植被指数的棉花苗提取精度评价Table 2 Evaluation of extraction accuracy of cotton seedling based on different vegetation indices %

4.2 棉花生长状态评级结果分析

棉花生长状态评级是构建 SEGT 模型的另一关键环节。植被生长是一个极其复杂的过程，土壤、水分、叶绿素等成分的差异，都可能会对最终的产量造成影响，因此对棉花进行基于时间序列的监测是十分必要的。试验首先获取棉花6 至8 月的无人机多光谱影像，根据天气情况，每隔7～10 d 采集1 次数据，共计得到研究区内 10 个时期的 NDVI 数据，日期分别为2018-06-10、2018-06-19、2018-06-27、2018-07-08、2018-07-19、2018-07-23、2018-08-04、2018-08-14、2018-08-21、2018-08-28。然后将 NDVI 分别与试验样区的最终实际测产结果进行相关性分析以确定各时期NDVI 的权重。最后根据NDVI 权重得到综合反映棉花生长状况的CNDVI 指数，对棉花植株进行等级划分，以获取棉铃数。试验结果表明，根据相关性计算，10个时期 NDVI 与样区最终实际总产量的相关系数分别为 0.69，0.72，0.75，0.81，0.88，0.87，0.82，0.83，0.75 和0.72。可以看出，棉花铃期的NDVI 与样区最终实际产量之间的相关性较大，而蕾期和吐絮期相关性则相对较小。根据相关系数的大小，利用公式（4）对棉花 10 个时期的 NDVI 分别赋予 0.09，0.09，0.10，0.10，0.11，0.11，0.10，0.11，0.10 和 0.09 的权重，构建 CNDVI并进行棉花估产。

利用时间序列NDVI 构建的CNDVI 如图4 所示。根据CNDVI 计算结果得到，最终CNDVI 取值在0～0.70 之间。采用自然间断点分级法，选取0.29 和0.54作为分界点，将棉花划分为3 个等级，分别对应一等苗（0.54＜CNDVI≤0.70）、二等苗（0.29＜CNDVI≤0.54）和三等苗（0≤CNDVI≤0.29）。从图中可以看出，研究区内棉花整体生长状态较好，以一等苗和二等苗为主（图中绿色和黄色区域），三等苗相对较少，多分布在地块边缘处（图中红色区域），与实际调研情况相符。

图4 棉花综合植被指数分级结果Fig.4 Classification results of cotton Comprehensive Normalized Difference Vegetation Index(CNDVI)

4.3 棉花估产结果分析

获取棉花出苗信息和棉铃数后，根据公式（5）对研究区进行估产。研究区3 个样地面积分别为17.72、15.35和9.4 hm2，总面积约为42.47 hm2，共计提取出4 364 255株棉花苗。棉花预测单产分别为 6 078.12、6 386.47 和5 900.44 kg/hm2，预估总产量分别为107 704.32、98 032.27和55 464.16 kg，总计约为261 200.75 kg。棉花实际总产量分别为101 542.08、90 389.69 和61 626.85 kg，总计约为 253 558.62 kg。研究区 3 块样地的估产精度分别为94.28%、92.2%和90%。

为了进一步验证估产结果的精度，利用20 个验证样区（分布见图1）的实际产量对预测结果进行验证，选用决定系数R2、均方根误差 RMSE、相对误差 ME作为评价指标检验估产模型的精度和可信度。R2越接近于1，说明模型精度越高；RMSE 值越小表明模型预测能力越强，精度高且稳定性好、可信度高。通过模型验证结果可知（图 5），估产模型的决定系数R2达到了 0.92，均方根误差 RMSE 为 0.1，相对误差 ME 为3.47%，由此可见，利用SEGT 模型进行棉花估产具有较高的精度和可信度。

图5 棉花预测产量验证图Fig.5 Verification diagram of cotton predicted yield

5 讨论

棉花产量预测是一项复杂的工作，既要考虑技术的实用性和可行性，同时又要考虑预测结果的可信度和精度。本文充分发挥无人机空间分辨率高、获取数据快、操作简单、成本低的优势，构建了一种 SEGT 棉花产量估算模型，为无人机遥感在作物估产研究中的应用提供了一种新的思路，但试验中一些关键点仍值得进一步探讨：本研究进行了 1 个年度的棉花产量估算，但不同年份的棉花出苗情况以及生长状态存在较大差异，模型能否在不同年度、不同区域都取得较好的效果，需要在接下来的研究中，在更多区域、更多年度开展试验工作，验证并完善模型算法，进一步扩大估产区域，分析多年间产量的时空变化及其原因，提升模型与方法的适用性；棉花的生长状态较难准确评估，为避免人为主观因素的影响，本研究采用自然间断点分级法进行划分，但该分级方法与影像数据本身关联性较大，难以针对不同区域划分统一标准，后续研究中还应加强生长状态分级和棉铃数关系的研究，进一步提高量化标准，实现客观、准确的棉花生长状态分级和棉铃数获取，以提高模型的准确性；本研究通过获取苗数与每株棉花单产累加的方法得到整个区域的棉花产量，棉铃数和单铃质量是模型的2个重要参数，其与棉花品种存在较大关联性，因此在后续的研究中，需对更多棉花品种进行系统测量和参数修正，以提高模型的稳定性和适用性。