基于多光谱影像特征的小麦锈病反演模型构建

王志业, 刘涛,, 张寰, 张全国, 张海涛, 何玉红

(1.河南财经政法大学资源与环境学院, 河南 郑州 450002;2.河南农业大学机电工程学院,河南 郑州 450002;3.濮阳职业技术学院,河南 濮阳 457000)

小麦是中国主要的粮食作物之一,在国家粮食安全和国民经济发展中有着重要的地位。由于中国小麦种植区域范围广、传统农业管理粗放等因素,小麦锈病、白粉病、赤霉病等病虫害频繁发生,严重影响小麦的产量和品质。在众多作物病虫害类型中,锈病作为中国小麦常见病害中的一种,具有影响范围广、演变性强、传播速度快等特点,是小麦减产的重要原因之一[1]。传统的小麦锈病监测方法主要依靠人工田间调查,该方法费时费力,难以快速、精准获取受灾区域范围和受灾程度。此外,为了防治小麦锈病,往往过量使用农药,导致生产经营成本增加,甚至形成有毒残留、造成环境污染。因此,针对小麦锈病发生的区域范围、灾害等级进行及时有效的监测,是促进中国智慧农业、精准农业快速发展中的关键环节。

遥感技术通过解析地物反射或辐射的电磁波信息,可以获得地物类型[2]、生态环境[3]、理化参量[4]等信息的空间分布状况,以用于田间作物的生长情况监测,具有实时高效、大范围、无损等特点[5]。已有研究表明,卫星遥感能够获取大范围区域的影像数据,并实现作物生长参量反演和病虫害监测应用[6-8],但是其较长的访问周期、较低的分辨率及混合像元等问题,导致难以及时、有效获取作物病虫害监测结果。近年来,随着无人机技术的迅速发展,通过在无人机上搭载可见光、多光谱、高光谱相机可以实现对作物长势参量与病虫害等级的反演。目前已有学者开展了基于无人机遥感的小麦锈病反演研究,例如,HEIDARIAN等[9]基于无人机获取4个不同地点的患病小麦可见光影像,通过对不同波段的像元亮度值(digital number,DN)及其组合研究发现,绿、红波段的DN值适宜反演小麦条锈病的病叶面积,绿、红、蓝波段的DN值组合易于估算小麦叶锈病的病叶面积。SU等[10]使用多光谱无人机获取多个关键黄锈发育阶段的小麦田影像,基于互信息法筛选敏感指数,通过随机森林分类器对患病程度不同的小麦植株进行分类,结果表明,利用随机森林分类器开发的黄锈监测系统可以实现良好的分类性能,平均精度、召回率和准确率分别为89.2%、89.4% 和89.3%。GUO等[11]使用高光谱无人机获取小麦种植区影像,分别提取植被指数和纹理特征,通过筛选特征变量和偏最小二乘回归,实现了小麦锈病反演与制图。综合目前研究,无人机遥感技术可以及时获取高分辨率影像数据,为小麦锈病监测数据获取提供了方便、快捷的技术手段[12],但是针对如何从无人机影像中提取有效的特征,同时避免出现特征冗余等问题仍然缺乏理论性的探讨。

从无人机遥感影像上提取、筛选特征是构建作物病虫害监测模型的关键环节,选择合适的特征变量不仅能够减少模型复杂度、提高运算效率,还可以有效提升模型精度。目前已有学者利用无人机影像提取的光谱指数实现对患锈病小麦的分类识别[9-11],为了进一步挖掘多光谱影像蕴含的潜在特征、提高对患病小麦的识别精度,本研究以小麦锈病田块为研究对象,利用无人机获取实验区域多光谱影像,提取影像的光谱特征和纹理特征,同时考虑多个变量共同参与模型构建时可能出现的特征冗余等问题,通过相关系数-赤池信息准则(Pearson correlation coefficient-akaike information criterion,PCC-AIC)、灰色关联度-赤池信息准则(grey relational analysis-akaike information criterion,GRA-AIC)和变量投影重要性-赤池信息准则(variable importance projection-akaike information criterion,VIP-AIC)3种方法筛选特征变量,分别构建PLSR、BPNN、RF模型,以期为作物病虫害监测应用中的特征变量筛选与建模提供方法参考。

1 数据获取

1.1 无人机影像数据获取

本研究选取河南省长葛市某小麦锈病试验田区域(见图1)。长葛市位于河南省中部,属于暖温带季风区,年平均日照2 280 h,年平均气温15.0 ℃,年降水量579 mm,适宜于小麦生长。本试验田南北长115 m,东西长93.6 m,航向由东向西,区域内小麦患有不同程度的锈病。

图1 研究区示意图

使用大疆精灵4多光谱版无人机(Phantom4-M,P4M)采集试验区域多光谱影像数据,影像获取时间为2021-04-15。如表1所示,P4M集成了一个可见光相机和5个多光谱相机,相机有效像素200万像素,焦距5.74 mm。试验通过DJITerra软件进行航线规划,计划飞行高度30 m,航向和旁向重叠率分别设置80%和75%,单张影像的空间分辨率为每像素1.6 cm。试验通过Pix4Dmapper进行影像拼接,获得该区域的正射影像,使用ENVI5.3软件对正射影像进行辐射校正。

表1 多光谱相机参数

1.2 地面实测数据

采集无人机影像的同时,随机选取87个地面采样点(如图1所示),使用10 cm×10 cm大小的单位框标注采样点区域,分别记录采样区域内的小麦叶片总数和患病叶片个数,通过千寻位置(www.qxwz.com)提供的厘米级差分定位服务(real time kinematic,RTK),标注采样点中心位置坐标并记录编号。依据采样点坐标使用ENVI软件在多光谱影像中提取采样点光谱值,并且每个点重复3次以上,选取多次重复提取光谱值的平均值作为采样点光谱值。本研究使用患病叶片数与叶片总数的比值RD(病叶比率)作为采样点的实测患病程度指标,计算公式如式(1)所示:

RD=d/S

(1)

式中:d为患病叶片个数;S为小麦叶片总数。

2 研究方法

2.1 指数构建

本研究在获取5个波段的遥感多光谱影像基础上,分别提取采样点位置的经验植被指数、计算光谱指数构建光谱指数候选集,用于后续特征筛选。如表1所示,使用前人研究中常用于病虫害监测的植被指数作为小麦锈病识别反演的经验植被指数(empirical vegetation indexes,EmVI),具体包括:花青素反射指数(anthocyanin reflectance index,ARI),优化土壤调整植被指数(optimization of soiladjusted vegetation index,OSAVI)三角植被指数(triangular vegetation index,TVI),宽动态范围植被指数(wide dynamic range vegetation index,WDRVI),土壤调节植被指数(soil adjusted vegetation index,SAVI),红边叶绿素指数(red edge position index,CIre),修正型三角植被指数(modified triangular vegetation index,MTVI2),改良非线性植被指数(improved non-linear vegetation index,MNLI),绿叶指数(green leaf index,GLI)绿色叶绿素植被指数(green position index,CIG)[13-22]共计10种光谱指数。如式(2)～(4)所示,计算任意2个不同波段间的差值光谱指数(difference spectral index,DSI)、比值光谱指数(ratio spectral index,RSI)和归一化光谱指数(normalized spectral index,NDSI),作为本研究的计算光谱指数(calculate spectral index,CaVI)。

DSIi,j=Ri-Rj

(2)

(3)

(4)

式中:Ri、Rj分别表示蓝、绿、红、红边、近红外波段中的任意一个波段;DSIi,j表示i波段与j波段的差值光谱指数;RSIi,j表示i波段与j波段的比值光谱指数;NDSIi,j表示i波段与j波段的归一化光谱指数。

无人机影像空间分辨率高,除了蕴含光谱信息外,还包含着丰富的空间信息,而纹理指数能够反映影像灰度特征值的空间变化,并且能够反映出作物遭受病害后的参量变化[23]。所以,本研究在提取单波段影像和光谱指数影像基础上,利用灰度共生矩阵提取出目标影像的8种纹理特征,分别为均值(mean),协同性(homogeneity,Hom),相异性(dissimilarity,Dis),对比度(contrast,Contr)信息熵(entropy,En),二阶矩(secondmoment,Sm),相关性(correlation,Cor),对比度(contrast,Contr),方差(variance,Var)。计算窗口为7×7,计算方向为90°,为了方便表示,用“_”连接原始影像和纹理特征,如B_Hom和B_Cor分别表示从Blue波段影像提取的协同性特征和相关性特征。

2.2 特征排序方法

分别使用皮尔逊相关系数(Pearson correlation coefficient,PCC)、灰色关联度分析(grey relation analysis,GRA)和变量投影重要性(variable importancein projection,VIP)3种方法对经验植被指数、计算光谱指数和纹理指数进行排序。其中,皮尔逊相关系数是一种统计学方法,用来衡量2个变量之间的相关性,其值介于-1到1之间,当相关系数的绝对值越靠近1,说明2个变量的相关性越大,具体计算公式如式(5)所示:

(5)

灰色关联度分析是灰色数学中的一种方法,用于度量自变量和因变量之间的关联性大小,在本实验计算出采样点小麦病叶比率RD与光谱指数之间的GRA值,GRA值越大,说明小麦病叶比率和该光谱指数的内在联系更紧密[24]。灰色关联度的计算公式如式(6)～(7)所示:

(6)

(7)

式中:x0为参考序列;xi为特征序列;ρ为分辨系数,一般在0～1中选取,取0.5。

变量投影重要性是指自变量在解释因变量时的作用和重要性[25],其计算过程是通过自变量的主成分传递来解释因变量;一般情况下,VIP值越高,该变量对模型的贡献越大,具体计算公式如式(8)所示:

(8)

式中:m为变量的个数;F为偏最小二乘法提取的主成分个数;wji为变量j中第i个主成分的权重值;SSYi为主成分解释方差的平方和,代表了主成分对因变量的解释能力;SSYtotal为因变量平方总和。

2.3 特征筛选方法

对特征指数进行排序后,分别使用判定系数(coefficientof determination,R2)和赤池信息准则(akaike information criterion,AIC)进行特征筛选。其中,R2是用来测试在因变量中回归模型能够解释的比例,当R2越靠近1时,可以说明变量间拥有较好的线性相关关系。具体计算公式如式(9)所示:

(9)

AIC=aln(L)+2(k+1)

(10)

式中:a为样本数;L为模型的残差平方和;k为模型中的变量数。

2.4 建模方法与精度检验

从样本中随机选取70%数据作为训练集、30%数据作为测试集,并基于筛选得到的特征,分别使用PLSR、BPNN和RF等回归方法对小麦病叶比率RD进行建模,并依据R2和均方根误差(root mean square regression,RMSE)检验模型精度,公式如(11)。其中,PLSR是一种集成主成分分析、线性相关分析和多元线性回归的多变量回归分析方法,可以同时实现数据降维、回归建模和分析数据间的相关性。BPNN是一种多层前馈神经网络,具有信号向前传播、误差反向传播的特点,主要是以误差逆向传播算法训练模型[27];本研究中构建BP神经网络时包含了1个输入层、10个隐藏层和1个输出层,学习率设置为0.01。RF是一种多因子机器学习算法,中心思想是利用Bootstrap重抽样算法从训练集中随机又放回的抽取样本,对每个随机抽取的样本集进行决策树建模,将这些决策树组成随机森林,并将多棵决策树投票决定的结果作为最终的预测结果[28],本研究建模时,设置决策树的个数为500。具体公式如式(11)所示:

(11)

式中:N为样本数;Yi为样本真实值;f(xi)为样本预测值。

3 结果与分析

3.1 光谱指数筛选结果

根据2.2节中介绍的方法,对光谱指数候选集中所有的指数进行排序,得到的结果如表2所示。可以看出,PCC值最大的计算光谱指数和经验植被指数分别是DSI1,3和ARI,值分别为-0.28和-0.26;GRA值最大的计算光谱指数是RSI3,1,值为0.69;GRA值最大的经验植被指数是MNLI,值为0.68;VIP值最大的经验植被指数是MTVI2,其值为2.58;VIP值最大的计算光谱指数是RSI4,5,其值为2.26。

表2 不同特征排序方法对光谱指数排序结果

依据第2.3节中的方法对排序结果进行光谱指数筛选,逐步增加光谱指数进行PLSR回归分析,理论上当AIC值越低、R2越高时,模型效果较好,实验结果如图2所示(图2-a为R2变化结果,图2-b为AIC值变化结果)。可以看出,利用PCC-AIC方法筛选光谱指数时,由排序前10位指数构建的模型具有最高精度,AIC值为198.79,R2值为0.664,RMSE值为0.22;GRA-AIC方法筛选光谱指数时,利用前5个指数构建的模型精度最高,AIC值为194.46,R2值为0.68,RMSE值为0.21;利用VIP-AIC方法筛选光谱指数时,由排序前8位指数构建的PLSR模型具有最高精度,此时AIC值为194.6,R2值为0.69,RMSE值为0.21。对比3种方法,利用VIP-AIC方法筛选变量并构建PLSR模型效果最好。

图2 光谱指数筛选结果

3.2 光谱+纹理指数筛选结果

在原始数据单波段影像和光谱指数影像(EmⅥ、CaⅥ)基础上,分别提取8个纹理特征(Mean,Hom,Dis,En,Sm,Cor,Contr and Var),将纹理特征与光谱指数合并作为特征候选集,使用3种方法进行排序和筛选,排序结果如表3所示。可以看出,在3种排序方法中,纹理特征均位于前列,这说明纹理特征能够更好地反映小麦叶片病变导致的生理变化。根据排序结果,逐步添加特征变量进行PLSR回归,结果如图3所示。可以看出,利用PCC-AIC方法筛选光谱指数时,由排序前19位指数构建的模型具有最高精度,AIC值为186.958,R2值为0.799,RMSE值为0.22;GRA-AIC方法筛选光谱指数时,利用前7个指数构建的模型精度最高,AIC值为179.548,R2值为0.708,RMSE值为0.21;利用VIP-AIC方法筛选光谱指数时,由排序前15位指数构建的PLSR模型具有最高精度,此时AIC值为178.187,R2值为0.802,RMSE值为0.21。试验结果表明,在优选植被指数的基础上,添加纹理特征作为候选集,使用3种筛选方法得到的特征参量构建PLSR模型精度都要高于仅使用光谱指数。

图3 光谱指数+纹理特征组合筛选结果

表3 不同特征排序方法对混合多元特征排序结果

3.3 模型评估与精度验证

为了进一步验证本研究筛选特征方法的有效性,使用筛选出的光谱+纹理特征组合,构建BPNN和RF等机器学习回归模型,训练集的结果如表4所示。在BPNN建模结果中,利用VIP-AIC方法筛选出的特征指数组合具有最高精度,其R2为0.918,RMSE为0.128;在RF模型建模结果中,利用VIP-AIC方法筛选出光谱指数精度最高,R2值为0.879,RMSE为0.114。筛选光谱指数在测试集中性能如图4所示,图4中红色区域代表模型的95%置信带,样本点越多位于置信带内说明模型拟合效果越好;可以看出,利用PLSR、BPNN和RF等3种模型预测样本点大都位于95%置信带中,利用VIP-AIC方法筛选的指数构建的BP模型精度最高(R2为0.93,RMSE为0.075),利用该方法反演制图结果如图5所示。

图4 测试集验证结果

图5 小麦影像反演结果

表4 机器学习模型结果

4 结论与讨论

农作物病虫害监测是保证农业生产安全的重要一环,实现及时准确的病虫害分类识别能够有效的推动农业生产者对农田的精准化生产管理。现有针对小麦锈病反演的研究,主要利用影像光谱指数进行患病小麦的分类识别,使用相关系数作为特征筛选指标。本研究通过多种特征筛选方法对光谱指数和纹理指数构成的混合多元特征集进行优选,构建出精度较高的小麦锈病反演模型,并通过AIC与R2相结合的方法来减少特征冗余问题。

前人研究中,纹理指数通常从原始影像中提取。在本研究中,进一步从优选的光谱指数影像中提取纹理特征,试验结果发现,在3种筛选方法得到的优选混合特征集中,纹理特征相比光谱特征更具数量优势,进一步说明纹理特征对于小麦锈病导致的作物表型变化更加敏感,同时也与无人机多光谱影像较高的空间分辨率有关。此外,在筛选出的纹理特征中有部分来自原始影像,更多的是通过EmVI、CaVI等优选光谱指数影像中提取得到,表明利用光谱指数筛选后的特征影像提取纹理特征,能够更好反映出作物表型、理化参量之间的差异,有助于提升后续反演模型精度。

为了进一步提升小麦锈病无人机遥感反演的模型精度,本研究尝试将线性模型筛选出的特征拓展应用于神经网络、随机森林等机器学习模型。结果表明,相比于偏最小二乘回归模型,机器学习模型的精度均有所上升,其中RF模型相比于PLSR模型,精度提升最大的是使用GRA-AIC方法筛选出的特征集,R2提高了0.16;BPNN模型相比于PLSR模型,精度提升最大的是使用VIP-AIC方法筛选的特征集,R2提高了0.11。结果证明,BP神经网络和随机森林等机器学习模型更适合于小麦锈病反演建模,原因可能是机器学习模型在原理上可以模拟非线性拟合结果,进而更适合反映出作物遭受病虫害后表型参量发生的复杂变化,使非线性模型在小麦病害反演建模中更具优势。

本研究针对小麦锈病田块开展了基于无人机多光谱的病害识别与反演制图试验,选取10个经验植被指数、40个计算光谱指数和192个纹理特征,通过PCC、GRA和VIP等方法排序特征指数,使用AIC和R2等指标进行特征筛选,分别构建PLSR、BPNN、RF等回归模型对小麦病叶比率反演建模。结果表明,VIP-AIC方法可以筛选出较好的特征集;相比其他2种回归模型,BP模型具有更好的拟合精度。相比仅使用光谱指数构建的反演模型,使用光谱+纹理组合特征构建的模型精度更高。本研究不足之处在于,筛选特征变量时仅用了PLSR线性模型,没有进一步拓展到非线性模型中;下一步需要尝试使用更多非线性模型组合,用于筛选特征变量和反演建模。