研究区位于新疆生产建设兵团第八师145团六分场六连, 中心位置为86° 3'30.85″E, 44° 32'30.23″N, 属中高纬度的典型大陆性气候, 光照充足, 雨量稀少, 气候干燥, 冬季严寒, 夏季严热, 研究区气候条件适宜棉叶螨正常生长^[9]。 研究区棉田面积34.7万m², 约520亩, 种植品种为新陆早166号和新陆早163号, 采用1膜3管6行种植模式, 蕾期一般处于当地的6月上中旬至7月上旬。

1.2.1 棉花冠层高光谱数据采集与处理

2018年于北京时间7月12日12:00— 14:00对蕾期棉花进行不同螨害程度的冠层光谱采集, 采集设备为ASD FieldSpec Pro FR2500型背挂式野外高光谱辐射仪, 光谱波长范围为350~2 500 nm, 采样间隔为1 nm, 光谱仪视场角为25° 。 采集时保持探头垂直向下, 光谱仪距离棉花冠层顶垂直高度约为0.5 m, 每次测量前后均进行标准白板校正, 每个采样点采集5组光谱数据, 取平均值作为该采样点的光谱值。

1.2.2 无人机多光谱数据采集与处理

2019年棉花蕾期内, 选择晴朗无云、 无风或风力较小的天气, 分别于6月22日、 27日、 29日和7月3日、 9日、 13日在研究区内进行6次多光谱数据采集, 采集设备为瑞士eBeeSQ固定翼无人机搭载ParrotSequoia五通道多光谱相机。 该相机镜头焦距为3.98 mm, 分辨率为1 280× 960, 可同时获取4个多光谱波段和可见光数据。 4个多光谱波段分别为绿光(中心波长550 nm, 波段宽度40 nm)、 红光(660 nm, 40 nm)、 红边光(735 nm, 10 nm)和近红外光(790 nm, 40 nm)。

无人机起飞前用Sequoia多光谱相机配备的参考板进行辐射定标处理, 根据预设航线垂直拍摄, 影像范围覆盖整个研究地块。 设置无人机飞行航线为S型, 飞行高度为120 m, 飞行速度为13 m· s^-1, 旁向重叠度为60%, 航向重叠度为75%, 影像空间分辨率为0.12 m。 数据采集后, 将影像数据导入Pix4Dmapper软件进行辐射定标和拼接处理。

1.2.3 螨害实地调查数据采集

在每期无人机多光谱影像采集的时间点附近由人工进行实地叶片螨害抽样调查。 在研究区棉田内使用VIVO X20手机拍摄健康和螨害叶片图像, 采样比例1: 1, 共135个采样点。 采样时记录叶片受害程度并根据手机自带的GPS精确定位功能记录采样点经纬度信息。 螨害等级参照国家标准《棉花叶螨测报技术规范》进行划分为4级: 0级为完害; 1级为叶面有零星黄斑; 2级为红斑占叶面积1/3以下; 3级为红斑占叶面积1/3以上。 4种不同螨害程度棉叶示例图像如图1。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 不同螨害程度棉叶图像 (a): 0级; (b): 1级; (c): 2级; (d): 3级Fig.1 Images of cotton leaves with different degree of mite damage (a): Level 0; (b): Level 1; (c): Level 2; (d): Level 3

1.2.4 环境数据采集与处理

无人机多光谱数据采集期内, 通过棉田研究区搭建的气象监测系统和土壤墒情设备, 获取2019年6月— 2019年7月期间逐日每小时的气象数据和土壤温湿度数据共13种田间环境数据。 再以4 d作为一个周期进行统计, 计算出9种气象数据, 即每周期最高温度、 平均温度、 最低温度、 平均湿度、 累计降雨量、 平均风力、 温湿系数、 温雨系数、 积温(逐日平均温度≥ 25 ℃); 4种土壤数据, 即10 cm土壤最高温度、 10 cm土壤最低温度、 10 cm土壤平均温度和10 cm土壤平均湿度。 同时考虑到农户施药对螨害的影响, 在施药日期(6月24日和7月8日)节点附近将周期长度细化为2 d。

植被指数(vegetation index, VI)通过不同波段反射率的线性或非线性组合变化, 在一定程度上可以突显植被的不同特征, 削弱背景等因素对植被光谱特征的干扰。 工作中根据地面尺度的棉花冠层高光谱数据提取出多个螨害敏感波段, 从中选取与低空尺度的无人机多光谱波段位置相对应的4个敏感波段, 初步构建23种螨害监测植被指数。 植被指数及其计算公式如表1所示。

表1

Table 1

表1(Table 1)

表1 植被指数及其计算公式 Table 1 Calculation formulas for vegetation index 植被指数 公式 NDGI (G-R)/(G+R) TVI 0.5× [120× (NIR-G)-200× (R-G)][10] MSR (NIR/R-1)/[(NIR/R)0.5+1][10] MSAVI 0.5× {(2× NIR+1)-[(2× NIR+1)2- 8× (NIR-R)]0.5}[11] RDVI (NIR-R)/[(NIR+R)0.5][10] RVI NIR/R[11] REWDRVI (0.15× NIR-RE)/(0.15× NIR+RE)[11] DVI NIR-R[10] GNDVI (NIR-G)/(NIR+G)[11] SAVI 1.5× (NIR-R)/(NIR+R+0.5)[12] ARI (1/G)-(1/RE) OSAVI 1.16× (NIR-R)/(NIR+R+0.16)[12] NDVI (NIR-R)/(NIR+R)[12] EVI 2.5× [(NIR-R)/(NIR+6× R-7.5× G+1)][13] GRVI NIR/G[13] GDVI NIR-G[13] GOSAVI (1+0.16)× (NIR-G)/(NIR+G+0.16)[13] REDVI NIR-RE[13] RERVI NIR/RE[13] RESAVI 1.5× [(NIR-RE)/(NIR+RE+0.5)][13] RERDVI (NIR-RE)/[(NIR+RE)0.5][13] RENDVI (NIR-RE)/(NIR+RE)[13] REOSAVI (1+0.16)× (NIR-RE)/(NIR+RE+0.16)[13] 注: G, R, RE, NIR分别为绿, 红, 红边和近红外波段光谱反射率 Note: G, R, RE, NIR represent reflectance of green, red, red-edge and near-infrared bands respectively

表1 植被指数及其计算公式 Table 1 Calculation formulas for vegetation index

通过相关性分析筛选反映螨害最敏感的植被指数和田间环境数据作为建模特征值, 分别用0(无螨害)和1(有螨害)作为目标值, 将135个样本数据集输入SVM分类器进行训练和分类, 采用单一环境数据、 单一植被指数、 环境数据与植被指数相结合的3种方法建立基于SVM棉田螨害发生监测模型。 模型训练和测试使用Matlab R2019b和LIBSVM 3.24工具箱, 并选取准确率(Accuracy)、 精确率(Precision)、 召回率(Recall)和对精确率和准确率整体评价的F₁值作为模型精度的评价指标, 确定出最优螨害发生监测模型。

使用田间环境数据分别与螨害面积占比值作相关性分析, 筛选出显著特征值, 再基于多元逐步回归分析法选出与螨害面积值关系最密切的环境特征建立棉田螨害面积预测模型, 选取决定系数R²作为模型评价指标。

棉花叶片在遭受不同程度的棉叶螨危害后, 棉花冠层光谱反射率曲线呈现规律性变化, 呈现明显的波峰和波谷等特征。 图3为不同螨害程度的棉花冠层光谱反射率曲线。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同螨害程度地面冠层光谱反射率曲线Fig.3 Spectral reflectance curves of ground canopies with different degrees of mite damage

由图3可见, 在绿光波段范围内553 nm波长附近出现第一个反射波峰后反射率逐渐下降, 至红光波段范围内680 nm波长附近出现反射波谷-红谷, 此时具有较低的反射率。 反射率于705 nm波长处开始呈现加速上升趋势, 且在近红外波段760~1 350 nm范围内反射率较高, 同时在此范围内光谱反射率随着螨害程度的逐渐加重表现出较为明显的下降趋势, 在中红外波段1 550~1 800 nm范围内的光谱反射率也表现出了同样的规律。

红边是绿色植物最明显的光谱特征之一。 为分析不同螨害程度棉叶冠层光谱的红边特征, 对红边(680~750 nm)范围内的光谱反射率求一阶微分, 如图4所示。 可以看出680~750 nm范围内的一阶微分光谱曲线图像呈现“ 双峰” 现象, 且随着螨害加重, 反射率一阶微分值逐渐减小, 表明红边波段对不同程度的螨害较为敏感。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同螨害程度地面冠层一阶微分光谱曲线Fig.4 First derivative spectral curve of ground canopy with the different degrees of spider mites damage

分析结果表明, 棉花叶片在遭受不同程度的棉叶螨为害后, 可通过棉花冠层光谱反射率及反射率一阶微分的变化进行表征, 绿光波段553 nm附近、 红光波段680 nm附近、 红边波段680~750 nm范围、 近红外波段760~1 350 nm范围都是与螨害发生相关的敏感波段, 与文献[8]研究结果一致。

3.1节中高光谱的螨害敏感波段分析结果显示, 本实验中无人机携带的多光谱相机波段范围属于螨害敏感波段, 为进一步确定无人机携带多光谱相机进行虫害监测的可行性, 对无人机多光谱遥感影像进行了不同螨害程度采样点的反射率提取, 图5为反射率曲线。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 不同螨害程度的多光谱影像反射率曲线Fig.5 Reflectance curve of multispectral image with different degrees of spidermite damage

由图5可见, 绿光波段550 nm附近、 550~660 nm波段范围、 660~735 nm波段范围处不同等级螨害发生点的光谱反射率基本重叠在一起, 可能的原因是无人机的飞行高度影响影像分辨率, 造成不同螨害程度的采样点在光谱响应上的变化差异较小, 光谱差异很难区分, 不能进行受害等级区分, 但有害和无害的光谱反射率差异比较明显。 故利用无人机多光谱遥感进行螨害监测模型构建时, 只构建基于SVM的二分类模型, 将螨害1级及以上均视为螨害发生, 即未发生螨害记为0, 发生螨害记为1。

将选取的23种植被指数和预处理后的13种田间环境数据分别与地面采样点的螨害发生数据做相关性分析, 用显著性(sig)作为评价指标, 当sig≤ 0.05时, 植被指数与螨害发生达到显著相关水平, 当sig≤ 0.01时, 植被指数与螨害发生达到极显著相关水平。 从表2看出, 所选的植被指数中, SAVI, OSAVI, TVI和NDGI与螨害发生呈极显著相关, RDVI, RVI和MSR与螨害发生呈显著相关。 从表3可以看出, 平均湿度、 温湿系数和10 cm土壤平均温度与螨害发生呈极显著相关, 最高温度、 平均温度、 积温、 10 cm土壤最高温度和10 cm土壤平均湿度与螨害发生呈显著相关。 研究中选取sig≤ 0.05的植被指数及环境特征作为螨害监测建模的特征值。

表2

Table 2

表2(Table 2)

表2 植被指数与螨害发生之间相关性 Table 2 Correlation between the occurrence of spider mite and vegetation indices 植被指数 RDVI SAVI OSAVI TVI NDGI RVI MSR 相关系数 0.407* 0.222* * 0.304* * 0.492* * -0.304* * 0.210* 0.197* 注: * * 和* 分别表示0.01, 0.05的显著水平 Note: * * and* represent significant at the 0.01 and 0.05 levels of probability respectively

表2 植被指数与螨害发生之间相关性 Table 2 Correlation between the occurrence of spider mite and vegetation indices

表3

Table 3

表3(Table 3)

表3 环境数据与螨害发生之间相关性 Table 3 Correlation between the occurrence of spider mite and environmental data 环境数据 最高温度 平均温度 平均湿度 温湿系数 积温 10 cm土壤最高温度 10 cm土壤平均温度 10 cm土壤平均湿度 相关系数 -0.218* -0.178* 0.221* * 0.243* * -0.213* -0.201* -0.258* * 0.213* 注: * * 和* 分别表示0.01, 0.05的显著水平 Note: * * and* represent significant at the 0.01 and 0.05 levels of probability respectively

表3 环境数据与螨害发生之间相关性 Table 3 Correlation between the occurrence of spider mite and environmental data

选取表2及表3中的7种植被指数和8种田间环境数据作为建模特征值, 在135个地面采样点中, 随机选取出总样本数的2/3作为训练集, 其余的1/3作为测试集进行模型构建, 分别基于单一环境数据、 单一植被指数、 环境数据与植被指数相结合的3种方法建立SVM棉田螨害发生监测模型M1, M2及M3。 SVM模型的精度主要与核函数和SVM内的参数选取有关, 研究中选择径向基函数(radial basis function, RBF)作为SVM训练模型的核函数, 使用交叉验证的方法调整SVM内的惩罚参数c和gamma参数, 得到各类评价指标效果最优的监测模型。 模型结果如表4所示。

表4

Table 4

表4(Table 4)

表4 不同螨害监测模型分类结果对比 Table 4 Comparison of classification results of the different degrees of spider mite monitoring models 建模方式 样本 训练样本 测试样本 健康 螨害 准确率/% 健康 螨害 准确率/% 精确率/% 召回率/% F1值/% 单一环境数据M1 健康 螨害 26 9 23 32 64.44 13 10 8 14 60 56.52 61.9 59.09 单一植被指数M2 健康 螨害 27 12 9 42 76.67 15 8 6 16 68.89 65.22 71.43 68.18 环境数据与植被 指数结合M3 健康 螨害 26 3 13 48 82.22 12 4 5 24 80 75 70.59 72.73

表4 不同螨害监测模型分类结果对比 Table 4 Comparison of classification results of the different degrees of spider mite monitoring models

表4中, M3效果最好, 测试集分类准确率为80%, 精确率为75%, 召回率为70.59%, F₁值为72.73%; M2次之, M1分类效果最差。 表明, 同时使用环境数据和植被指数综合考虑了影响螨害发生的主要因素, 从而使模型具有更高的精度和监测效果。

基于棉田螨害监测模型M3, 分别绘制6月22日、 27日、 29日和7月3日、 9日、 13日共计6景无人机遥感螨害监测空间分布图[图6(a— f)]。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 不同时期(6月22日, 27日, 29日, 7月3日, 9日, 13日)棉田螨害监测空间分布图(a, b, c, d, e, f)Fig.6 Spatial distribution map (a, b, c, d, e, f) of cotton field spider mite monitoring in (June 22, 27, 28, Jule 3, 9, 13) different periods

图6中可见, 6月22日开始, 研究区螨害成片发生且基本分布在棉田东部以及中路以南区域[图6(a)], 棉田实际调查发现棉田东部及南部区域靠近路边, 地边杂草较多且农户活动多聚集于两地, 同时发现监测期间研究区中部以南地区因棉田灌溉设施漏水, 棉花长势较弱, 这些因素综合导致了棉田边缘易发生的螨害扩散到了棉田中部及以南。 根据6月22日棉田螨害监测情况, 农户于6月24日使用25%阿维乙螨唑悬浮剂对叶螨进行施药防治, 从图6(b)可以看出, 6月27日分布图中螨害面积逐渐减小, 叶螨数量较之前有所下降。 但由于棉田内施药时间不同及叶螨抗药性强, 观察6月29日棉田螨害监测图[图6(c)]发现棉田北部螨害扩散较其他地区略微严重。

随着7月新疆地区温度逐渐升高, 高温条件对棉叶螨扩散极其有利, 同时药剂对叶螨的防治效果下降, 导致7月3日螨害逐渐加重, 呈现局部棉田叶螨暴发状态且多集中于东、 西部地块边缘区域[图6(d)]。 因此, 农户于7月8日再次对棉叶螨进行施药防治, 但随着植保人员和农户在田中不断移动进行实地调查等操作, 导致叶螨随之迁移扩散, 通过图6(e)可以看出, 中部以北区域出现局部田块爆发成灾现象, 防治难度增加, 易由逐渐扩散转变成迅速扩散。 7月13日再次对棉田进行监测, 发现螨害中部以北区域已达到大面积爆发状态, 而中部以南区域螨害防治程度较好, 无局部爆发状态[图6(f)]。

监测图整体趋势与实际调查情况较为相符, 通过基于环境数据与植被指数相结合建立的SVM棉田螨害监测模型可以反映田间螨害变化, 基于无人机的低空遥感可以辅助用来实现棉田大范围的螨害发生监测。

首先基于模型分类结果计算螨害面积占比值, 考虑到实际应用环节, 环境数据更易获取, 故只将初选的13种田间环境数据与棉田螨害面积值作相关性分析, 选取sig≤ 0.05的田间环境数据, 结果如表5所示。

表5

Table 5

表5(Table 5)

表5 环境数据与棉田螨害面积之间相关性 Table 5 Correlation coefficient between environmental data and cotton field spider mite area 环境数据 最高温度 平均湿度 温湿系数 积温 10 cm土壤最高温度 10 cm土壤平均温度 相关系数 0.707* 0.844* 0.931* * 0.837* 0.856* 0.974* * 注: * * 和* 分别表示0.01, 0.05的显著水平 Note: * * and * represent significant at the 0.01 and 0.05 levels of probability yrespectively

表5 环境数据与棉田螨害面积之间相关性 Table 5 Correlation coefficient between environmental data and cotton field spider mite area

分别使用SVM回归和多元逐步回归方法构建螨害面积模型, 以决定系数R²作为指标来衡量模型的性能。 在Matlab下调用LIBSVM工具包建立SVM螨害面积预测模型M₁, 并通过交叉验证方式寻优参数, 得到最终模型R²=0.796。 采用SPSS进行多元逐步回归分析, 从表5中逐步选出与螨害面积值关系最密切的田间环境数据作为建立预测模型M₂的特征值, 最终得到多元逐步回归方程为

y=110.826-0.193x1-0.089x2-0.5x3+0.423x4-4.022x5(1)

式(1)中, x₁是最高温度, x₂是平均湿度, x₃是温湿系数, x₄是10 cm土壤最高温度, x₅是10 cm土壤平均温度, y是预测的螨害面积占比值。 模型评价结果见表6。

表6

Table 6

表6(Table 6)

表6 模型评价结果 Table 6 Model evaluation results 模型 R R2 调整后R2 标准估算的误差 M2 0.874 0.848 0.835 0.914 9

表6 模型评价结果 Table 6 Model evaluation results

由表6可知, 多元逐步回归预测模型的决定系数R²=0.848, 说明该模型拟合度优于SVM回归预测模型, 并对多元逐步回归预测模型进行显著性检验, 计算出P值为0.001, 小于0.05。

本研究最终使用多元逐步回归模型对测试样本进行预测, 将螨害面积占比值单位转换成亩数, 结果如表7所示。

结果表明模型对棉田螨害发生面积的预测效果较好。

表7

Table 7

表7(Table 7)

表7 棉田螨害面积预测模型预测结果对比 Table 7 Comparison of prediction results of cotton field spider mite area prediction models 日期 实际值/亩 预测值/亩 6.22 64.168 58.374 6.27 18.361 18.591 6.29 20.348 18.101 7.3 50.055 41.666 7.9 43.087 42.39 7.13 72.831 69.83

表7 棉田螨害面积预测模型预测结果对比 Table 7 Comparison of prediction results of cotton field spider mite area prediction models