基于Logistic算法与遥感影像的棉花虫害监测研究

地力夏提·依马木，周建平,2，许 燕，樊湘鹏，亚里坤·沙吾提

(1 新疆大学 机械工程学院，新疆 乌鲁木齐 830047； 2 新疆大学 新疆维吾尔自治区农牧机器人及智能装备工程研究中心，新疆 乌鲁木齐 830047； 3 新疆维吾尔自治区农牧业机械化技术推广总站，新疆 乌鲁木齐 830054)

近年来，提升大田作物品质与保障农田作物产量是农业机械化领域的热点研究问题，给农业生产造成损失的主要原因是虫害的发生与蔓延。据现有统计信息，2020年中国棉花总产量为591万吨，因病虫害损失15%～20%[1]。作为中国最大的优质棉生产区域，新疆棉花播种面积占全国的76%左右，总产占全国棉花产能的84.9%[2-3]。棉花高产量与优产率得益于新疆独特的地域优势与气候环境优势[4-5]。棉蚜虫、棉红蜘蛛、棉铃虫是棉花的主要害虫，对棉花生产造成的危害较大。精准、快速地对棉花虫害发生情况进行监测，及时了解棉花虫害动态变化，对于棉花大面积稳健生长具有重要意义[6-7]。传统的虫害识别监测、定位方法需要耗费大量的人力、物力，仅适用于小片区棉田，难以满足新疆棉田规模化植保防治的要求[8]。赵亮等[9]基于便携式光谱仪的光谱信息分析螨害严重度与棉叶螨叶片光谱反射率的相关性，构建了基于敏感波段的虫害诊断模型，从而实现了对棉花螨害的识别提取。Prabhakar等[10]利用叶片波谱分析和受害叶片症状解读方法，构建叶蝉敏感波段，并利用比值指数组合构建叶蝉指数(Leafhopper index，LHI)有效识别叶蝉发生区域；竞霞等[11]利用卫星遥感数据提取归一化差值植被指数(Normalized difference vegetation index，NDVI)、重归一化差值植被指数 (Renormalized difference vegetation index，RDVI)、增强植被指数 (Enhanced vegetation index，EVI)、差值植被指数 (Difference vegetation index，DVI)、全球环境监测植被指数(Global environment monitoring index，GEMI)和修改型土壤植被指数 (Modified soil-adjusted vegetation index，MSAVI)共6种影响因子，使用偏最小二乘法 (Partial least squares,PLS) 和变量投影重要性(Variable importance projection，VIP)准则构建棉花黄萎症影响度估测模型，有效识别了病害发生区域并分析了其精度；胡根生等[12]利用植被中对蚜虫影响较大的植被指数生成8种特征因子，基于环保卫星远程遥感影像资料，提取相应光谱特征信息，利用最小二乘孪生支持向量机建立小麦蚜虫识别模型，且模型精度达到86.4%，较好地实现了对小麦蚜虫的监测。Stella等[13]基于农作物受病虫害侵害的时间构建作物病虫害的Logistic模型，以总农作物种群为常数，以病虫害媒介种群为变量，以基本病虫繁殖数为平衡点分析病虫害发生情况，并在棉花、小麦、玉米等农作物上进行模型验证；冯炼等[14]在冬小麦病虫害面积监测研究中选择二值Logistic回归法，分别建立植被指数与病虫害发生区域的Logistic模型，从模型识别结果来看，提取的冬小麦受病虫害胁迫面积与实测值结果相符，能满足农作物病虫害发生区域遥感监测要求，该模型在棉田病虫害识别方面也得到了较好的精度。学者们在预测虫害发生趋势研究方面利用Logistic模型确定性和随机性、虫害敏感波段、植被指数影响因子等方法，对比基于无人机遥感获取的高光谱、多光谱传感器影像数据，识别提取作物虫害发生区域并预测虫害发生趋势。结合无人机光谱遥感精准、高时效、高分辨率的优势，可以实现虫害区域预测研究[9,11,15-17]。而 Logistic模型可同时引用敏感波段、植被指数等数据[14,18]构建作物虫害识别模型，且从精度数据来看，对作物虫害识别的准确率均可满足虫害发生区域监测要求，故本研究采用Logistic模型实现棉花虫害监测与分类。现阶段无人机遥感影像数据在识别监测棉田虫害研究领域中仅停留在对单一虫害的识别监测，鲜见对多种虫害的识别监测，且鲜见对影响遥感影像精度的因素及其解决方案的系统性研究。借助遥感影像处理软件分析影响遥感影像精度的因素并对所获取棉田遥感影像进行几何校正、辐射校正，提升影像精度，其次利用Logistic模型、自相关性时间序列检测的互相关函数(Cross correlation function，CCF)等算法，实现棉田多种虫害快速精准识别监测，以期为棉花稳产高质提供技术服务，为虫害识别提供数据支撑。

1 研究平台与方法

1.1 研究区概况

研究区位于新疆维吾尔自治区巴音郭楞蒙古自治州尉犁县境内，地处北纬 41°29′27″～41°29′58″，东经 86°15′30″～86°16′3″，占地 127.3 hm2。属暖温带大陆性荒漠气候，冷热差异悬殊，全年平均气温10.1 ℃，年平均降水量 43 mm，年平均蒸发量2 700 mm。本次研究选择‘新陆中68号’棉花作为研究对象，生育期134～143 d，研究区棉花采取2膜12行种植模式，行距配置为(10+66+10+66+10) cm，种植行间距为宽行 66 cm，窄行 10 cm，株距10 cm。种植区域采用滴灌模式，其他管理措施按照当地大田棉花高产栽培模式进行。

1.2 试验设计

棉田遥感影像数据获取试验于2020年7月15日早间进行，天气晴，气温22 ℃，风向西北，风速3.1 m/s，适宜无人机航拍作业。无人机遥感平台选用DJI大疆经纬 Matrice 600多旋翼无人机，搭载MicaSense RedEdge-M 5通道光谱相机构成多旋翼无人机遥感系统。RedEdge-M 5通道光谱相机分辨率为 1 280×960，光谱波段包括蓝光 475 nm、绿光 560 nm、红光 668 nm、近红外 840 nm 和红边 717 nm，其对应的带宽分别为 20、20、10、40 和 10 nm。

数据采集时，根据预定航线采取垂直拍摄模式，覆盖整片研究区地块，飞行高度50 m，影像航向重叠率80%，旁向重叠率70%，满足影像处理需求。从研究区棉田划分4.67 hm2棉田作为验证区，取20个采样节点，每个采样节点为3 m×3 m空间范围，采样点地理位置信息采取空地设备结合方式，地面选取轻测Lite RTK设备获取地理位置信息，该设备满足5 Hz刷新率，RTK精度达到±(10+1×10-6D)mm，其中D为被测点间距；选用本研究执飞无人机GPS模块获取空间地理位置信息。获取的地理位置数据经由几何校正后，生成试验区正射影像图，从而保证采样点定位精度。选取此空间范围内所有棉花植株作为采样对象，采样对象包含棉花植株叶片、花蕾、花铃所发生的虫害，采样节点位置示意图如图1所示。采样过程样叶分成左基部、右基部两端，每端各测2次，每次测2条光谱，生成4组光谱波长数据，取其平均值作为样叶光谱值，棉花叶片采样位置如图2所示。在20个采样节点中选取300株棉花验证模型准确度，从中选取总样本数2/3的棉花植株(200株)用于模型训练，其中，100个样本选取受棉蚜虫、棉红蜘蛛、棉铃虫虫害的棉花植株，100个样本为健康棉花植株；剩余1/3的棉花植株样本数据用来测试精准度，其中，50个样本为受棉蚜虫、棉红蜘蛛、棉铃虫虫害的棉花植株，50个样本为健康棉花植株。本研究中健康棉叶与受虫害棉叶区分标准为：受虫害棉叶指叶片受虫害侵蚀区域面积大于单个叶片面积的5%，健康棉叶指叶片受虫害侵蚀区域面积少于单个叶片面积的5%。

图1 采样节点位置示意图Fig.1 Schematic diagram of sampling node location

图2 棉花叶片采样位置示意图Fig.2 Schematic diagram of sampling position of cotton leaves

1.3 研究方法

针对现阶段无人机遥感在规模化农田虫害识别研究中存在的遥感影像质量严重依赖天气、图片质量影响虫害识别精度及当前研究停留于单一虫害识别等问题，本研究利用几何校正、辐射校正方法，提升图像地理位置精度与像素覆盖密度，有效提升遥感影像质量；利用逻辑回归(Logistic regression，LR)方法识别虫害发生区域，对虫害识别精度及虫害发生区域特征进行分析。

LR又称为广义线性回归，属于判别分析，用于处理因变量为分类变量的数据。在实际应用过程中响应变量为二值变量，为解决二分类问题，延伸出二值逻辑回归法[19]。在解决二分类问题过程中，因变量为二分类变量，利用0与1代替表达不同状态。本研究中棉蚜虫、棉红蜘蛛、棉铃虫3种虫害发生与否为非线性变化求解问题，因此，以上3种虫害监测问题属于非线性二分类问题，基于此建立Logistic回归判定模型：

式中，β0为常数项或截距，β1，β2，…，βp为模型回归系数，x1，x2，…，xp代表p个自变量，P表示在p个自变量作用下的概率。也可表示为在影响Y取值范围的p个自变量条件作用下棉蚜虫、棉红蜘蛛、棉铃虫发生的概率为P=P(Y=1|x1，x2，…，xp)。

对式(1)进行Logistic变换，其线性表达式如下：

从式(2)可以看出，回归系数β在控制自变量条件下反映变量影响力，响应变量Y求解二分类问题时分别取值0和1，表示发生与不发生，P为发生概率，以P=0.5作为界限，P≥0.5表示逻辑为真，即棉花植株发生虫害，当P＜0.5表示棉花植株健康(未发生虫害)。

2 结果与分析

2.1 遥感影像预处理

影像拼接及几何校正过程使用瑞士Pix4D Mapper软件进行，可实现空三加密，结合控制点信息与所采集图像地理位置信息数据进行几何校正，生成研究区正射影像[20-22]，几何校正过程利用归一化函数提高配准点精度，从而有效提高图像精度。

经几何校正后，有效提升图像最优模型精度，模型精度达到（0.059±0.012） m。所生成正射影像如图3所示，几何校正后绝对位置和方向不确定性参数如表1所示。

图3 研究区正射影像图Fig.3 Orthophoto map of the study area

表1 绝对位置和方向不确定性参数Table 1 Absolute geographic location and directional uncertainty parameters

由于受光照强度及其他因素影响，遥感影像光谱响应误差表现为像元值像素覆盖密度数值波动，造成图像解析不准确，无法精准有效反映光谱反射率。为消除数值波动造成的像素密度误差，本研究利用Back propagation(BP)神经网络模型对棉田遥感影像数据进行监督分类。辐射校正过程使用ENVI软件Band Math工具定义辐射定标公式：

式中，Lλ为单元像元 (Digital number，DN)内像素覆盖密度，LMAXλ为单元像元最大值光谱反射率，LMINλ为单元像元最小值光谱反射率，QCAL为单元像元值，QCALMAX是最大量化定标像素值。红光波段、绿光波段、蓝光波段、近红外波段辐射校正前单位像元内覆盖密度分析结果如图4a所示。

图4 辐射校正前（a）、后（b）单位像元内像素覆盖密度Fig.4 Pixel coverage density in unit pixel before (a) and after radiation collection (b)

对所获取无人机遥感棉田影像提取的红光波段 (Red)、绿光波段 (Green)、蓝光波段 (Blue)、近红外波段(Near-infrared，Nir)组合作为输入层构建Red-BP、Green-BP、Blue-BP、Nir-BP 多层感知器模型优化波段反射率，光谱波段反射率作为输出层。以辐射校正前影像像素反射率数值作为预测值，该预测值数据由Micro-Hyperspec小型光谱检测仪实测叶片波长数据计算得出，辐射校正后BP模型反射率数值作为实测值，选用Sigmoid函数作为激活函数，实数值自变量区间为[0，1]，采用梯度下降法训练数据集。经辐射校正后红光波段、绿光波段、蓝光波段、近红外波段单位像元内覆盖密度结果如图4b所示。由辐射校正模型预测值-实测值散点图(图5)可知，红光波段、蓝光波段、近红外波段和绿光波段基于BP法的Red-BP、Blue-BP、Nir-BP和Green-BP模型光谱反射率预测值与实测值分布较为均匀，其预测值与实测值有较好的拟合效果，且决定系数(R2)分别为0.945、0.906、0.913和0.875，模型精度较高，满足图像解析要求。

图5 BP神经网络模型对辐射校正模型光谱反射率估测结果Fig.5 Estimation results of spectral reflectance of image pixel after radiation correction by BP neural network model

根据图4提取辐射校正前、后单元DN内像素覆盖密度数据，再针对像素覆盖密度数据绘制ROC曲线(图6)，对比数据预测价值比较如表2所示。

图6 4种波段像素覆盖密度Fig.6 ROC curves of pixel coverage densities of four bands

表2 4种波段像素覆盖密度预测价值比较Table 2 Comparison of the prediction values for the pixel coverage densities of four bands

辐射校正后单元DN内红光波段、蓝光波段、近红外波段和绿光波段像素覆盖密度的ROC曲线下面积 (The area under the ROC curve，AUC)分别为0.875、0.750、0.688和0.875；在最佳点位置，像素覆盖密度灵敏度分别为88.3%、79.3%、80.9%和87.1%。从灵敏度、ROC曲线下面积及校正情况来看，辐射校正可有效降低光谱响应误差，增强影像像素密度。

2.2 指数模型构建

为更好地反映出大田棉花种植区域、棉花长势情况及棉叶光谱反射率基本情况，为Logistic模型提供数据，本研究构建NDVI、EVI、SAVI模型。NDVI是体现农作物植被覆盖度的参数，是描述植被生长状态的一个重要指标[23-25]，对于本研究识别棉花虫害区域有着重要的意义。为增强植被反射率信息，EVI在归一化指数基础上加入蓝光波段数据，矫正土壤参数和气溶胶散射数据的影响，EVI常用于植被覆盖密度较大区域[24,26]。SAVI是基于NDVI和大量遥感数据观测修正后提出用以减小土壤背景影响的指数模型。利用ENVI软件Spectral Math工具构建指数模型，包括NDVI模型、EVI模型和SAVI模型，3种植被指数模型公式如下：

式中，L是随着植被覆盖度变化的参数，取值阈值为0～1，当植被密度很高时为0，反之为1。很明显，如果L=0，SAVI=NDVI。对于棉花植株，L取0.5时SAVI消除土壤背景影响的效果较好[27]。

NDVI模型可直观地表述棉花与土壤分布情况，EVI模型可直观地表述棉花生长趋势情况，SAVI模型可直观地表述土壤分布及棉花生长趋势；3种模型均可精准提取棉花长势信息及分布面积，在棉田虫害模型构建过程中提供棉花分布情况数据及土壤边界。

2.3 棉花虫害敏感波段提取

使用Micro-Hyperspec小型光谱检测仪对采样节点受虫害区域及健康棉花区域进行光谱值测定，利用Excel软件对样叶光谱波长数据进行统计分析，利用反射率公式生成虫害反射率及健康棉花叶片波长数据，为Logistic模型提供敏感波段数据公式如下：

式中，LMAX、LMIN分别为像元辐射亮度最大值和最小值；QCAL为某一像元DN值；QCALMAX、QCALMIN分别为像元所取最大值和最小值。

使用Origin软件构建棉田虫害光谱敏感波段模型及健康棉叶光谱波长信息模型，由图7可以看出，棉蚜虫的严重程度与713～831 nm波段的棉花蚜虫叶片光谱反射率显著相关；棉红蜘蛛的严重程度与605～724 nm波段的棉花红蜘蛛叶片光谱反射率显著相关；棉铃虫的严重程度与709～794 nm波段的棉花铃虫叶片光谱反射率显著相关，此范围内遥感因子可以作为辨认棉花棉蚜虫、棉红蜘蛛、棉铃虫3种虫害的敏感波段。

图7 健康棉叶及棉花虫害光谱波长反射率曲线Fig.7 Spectral wavelength reflectance curves of healthy cotton leaves and cotton pests

2.4 Logistic算法虫害识别研究

2.4.1 基于 Logistic 算法虫害识别模型 模型构建前，为使构建模型特征相同，对特征变量进行归一化分析。基于此，利用SPSS软件分析300株棉花数据，将本研究中应用的3种植被指数作为回归模型输入变量因子，构建虫害识别Logistic回归模型，同时计算均方根误差 (Root mean square error,RMSE)。当引入1个输入变量因子时，回归模型为y=10.25+9.4SAVI，RMSE为0.193；输入2个变量因子时，回归模型为y=-13.16+155.1SAVI-24.6NDVI，RMSE为0.187；输入3个变量因子时，回归模型为y=-9.85+138.38SAVI-75.7NDVI-44.8EVI，RMSE 为 0.376。NDVI、EVI、SAVI与棉花受害情况的相关系数分别为-0.914，-0.813和0.937，分别达到0.01、0.05和0.01水平。综上，构建虫害识别模型选取2个变量因子的识别效果最佳，由SAVI模型和NDVI模型建立的Logistic回归模型是识别棉蚜虫、棉红蜘蛛、棉铃虫的最优模型。

利用ENVI软件基于规则面向对象工作流程工具，波谱类型选取SAVI模型、NDVI模型与棉蚜虫、棉红蜘蛛、棉铃虫3种虫害光谱反射率曲线，基于Logistic回归判定的虫害识别模型，生成虫害识别区域如图8所示，由图8可明显看出受虫害影响区域。

图8 3种害虫为害区域识别结果Fig.8 Identification results of three pests infested arer

2.4.2 模型分类结果 基于RMSE选择具有代表性的Logistic回归虫害识别模型进行评价比较，具体分类结果分为选择1、2、3个输入变量因子，模型分类比较结果如表3所示。从比较结果中可以看出，输入变量为2的虫害识别模型识别效果最佳，而输入变量为3的虫害识别模型在准确率、召回率及F1值上均低于其他2个虫害识别模型。结合模型构建过程，由SAVI模型和NDVI模型建立的棉蚜虫、棉红蜘蛛、棉铃虫识别模型为最优模型，其训练样本准确率和测试样本准确率分别达到93.7%和90.5%，召回率、F1值分别为96.6%和93.5%。

表3 模型分类比较结果1）Table 3 Comparison results of model classification

分别对3种虫害发生区域像元数据参数估算，分析基于Logistic模型识别的棉田虫害区域预测值-实测值数据，验证Logistic模型在棉田虫害区域监测的精准度。3种虫害预测值-实测值相关性分析如图9所示。由图9可知，Logistic模型在棉田虫害区域监测的精准度较高，对棉蚜虫、棉红蜘蛛、棉铃虫发生区域识别模型的决定系数分别为0.942、0.851和 0.663。

图9 3种虫害发生区域识别精准度Fig.9 Accuracy of the occurrence area of three types of pests

2.4.3 虫害发生区域特征分析 对棉田受棉蚜虫、棉红蜘蛛、棉铃虫灾害区域进行自相关性时间序列检测(CCF)，分析随时间变化的棉田受虫害面积的相关性，结果见图10。由图10可知，棉蚜虫、棉红蜘蛛与棉铃虫的发生与蔓延均具有相似的条件。随时间推移，棉蚜虫与棉红蜘蛛的发生蔓延在1周内相关系数达到0.863，随后逐渐降低，棉蚜虫与棉铃虫的发生蔓延相关系数从开始的0.394逐渐降低，而棉红蜘蛛与棉铃虫的发生蔓延相关系数在2周内保持为0.35～0.40；从CCF相关性系数变化范围与对比数据来看，3种虫害发生初期其发生条件受彼此影响。

图10 3种棉田虫害受灾面积相关性分析Fig.10 Correlation a Analysis of the areas affected by three types of pests in cotton field

本研究主要开展棉田虫害识别监测与分类，对于虫害识别，健康棉叶与受虫害棉叶光谱反射率有明显偏差，受虫害侵蚀越严重，其反射率偏差越明显。虫害发生早期症状不明显，其反射率数值偏差较低，但在开展虫害识别监测与分类研究时，经量化细分，反射率数值偏差仍能较好反映，因此对样叶的采样标准会对虫害识别模型结果造成影响，但仅表现于虫害发生早期，不会造成显著影响。

3 结论

利用Logistic回归判定方法提取棉蚜虫、棉红蜘蛛、棉铃虫发生区域，由SAVI模型和NDVI模型建立的棉蚜虫、棉红蜘蛛、棉铃虫识别模型为最优模型，其训练样本准确率和测试样本准确率分别达到93.7%和90.5%，召回率、F1值分别为96.6%和93.5%。本研究在识别3种虫害发生区域过程中使用了Logistic模型，并利用自相关性时间序列检测法对虫害发生区域进行动态识别监测，从试验数据来看，Logistic模型在棉田虫害区域监测的精准度达到较高水平，该方法可应用于多种虫害的监测识别，并可基本满足棉田精准植保作业相关要求。借助识别区域Shapefile矢量文件，可为精准植保作业提供经纬度等数据，便于实现智慧农田管理的相关作业。

本研究仍存在一定的局限性，遥感影像数据优化过程中使用方法仍存在改进空间，此外对于会造成虫害发生与蔓延的人为土地管理干预、虫害天敌、极端气候气象条件等影响因子未做进一步研究，后续开展研究时，应多方面考虑虫害发生与蔓延的相关因子，以提高虫害动态监测模型精度。