基于多光谱卫星影像与机器学习算法的松材线虫病受害林分识别研究

邱万林,宗世祥

(北京林业大学林木有害生物防控北京市重点实验室,北京 100083)

松材线虫病(Pine Wilt Disease,PWD)被称为“松树癌症”,源自于北美,在日本、韩国和中国均有过疫情暴发,于1982年首次在中国被发现(叶建仁,2019)。目前中国是松材线虫病发生最严重的国家(叶建仁和吴小芹,2022)。松树从感染松材线虫病至植株死亡只需要2～3个月,最快40 d左右就会呈现出枯死状(杨宝君,2002)。2021年,松材线虫病在我国的发生面积为171.65万ha,全国共有19个省(自治区、直辖市)742个县级行政区发生疫情(李硕等,2022)。面对严峻的松材线虫病防控形势,及时清理枯死松树和有效管控疫木源头是遏制疫情蔓延的最好办法(国家林业和草原局,2022)。

我国有大量的松林分布于松材线虫的适生区范围内,面临着松材线虫病的潜在威胁,急需高效、实时、大面积监测疫情的技术手段。过去的松材线虫病疫木监测主要以人工地面调查为主,但是传统的监测手段面临着山区地形复杂、林分分布不均等问题,且人工调查成本高、效率低,难以覆盖所有受害林分,导致遗漏现象普遍,难以精确掌握病害发生情况(张红梅和陆亚刚,2017)。遥感监测作为一种重要的森林病虫害监测手段,具有灵活、高效、大面积监测病害发生区域等特点,可以作为传统地面人工监测的有效辅助手段。无人机遥感在森林病虫害监测已有广泛的应用,如基于Fast R-CNN和无人机遥感实现枯死木的快速定位(黄华毅等,2021)。但是无人机遥感无法获取历史影像,成本较卫星遥感更高而监测面积更小。卫星遥感作为一种传统的遥感技术,弥补了无人机遥感无法获取历史影像数据和难以形成长时序影像的缺点,并且在相同监测面积下较无人机平台成本更低、监测的范围更广(史舟等,2015)。

针叶光谱特征的变化是松材线虫病遥感监测的基础。松木在感染松材线虫病后其生理指标(含水量、含氮量和叶绿素含量等)发生变化,导致针叶颜色发生变化(理永霞和张星耀,2018),进而使针叶光谱反射率产生变化。所以通过遥感手段实现松材线虫病发生监测是可行的。目前松材线虫病卫星遥感监测是变色木遥感监测中的热点,国内外已有大量研究和应用结果表明应用卫星遥感能有效监测松材线虫病发生。有研究基于机器学习算法实现松材线虫病疫木信息提取,利用WorldView卫星全色增强融合图像和支持向量机提取受害松树(Takenakaetal.,2017)。有研究基于阈值分割法实现松材线虫病疫木信息提取,采用IKONOS和QuikBird等高分辨率卫星,基于局部最大值滤波器实现受害木信息提取(Leeetal.,2003,2007)。有研究利用长时序的MODIS卫星生成NDVI序列实现了松材线虫病受害林分提取和未来发生区域的预测(Haoetal.,2021)。也有许多研究基于深度学习实现松材线虫病受害木提取。LeNet(Zhouetal.,2022)、HRNet、DANet(Wangetal.,2022)和MA-UNet(Yeetal.,2022)等多种深度学习算法都被应用于松材线虫病疫木信息提取。

目前,有关松材线虫病卫星监测的研究集中于高分辨率卫星的单木尺度受害木提取,如用采用IKONOS(Leeetal.,2003)、Quikbird(Leeetal.,2007)和WorldView(Zhangetal.,2021)等高分辨率卫星实现松材线虫病疫木信息提取。但是高分辨率卫星存在成本高、历史数据少等问题。中分辨率卫星相较于高分辨卫星具有经济实惠、易获取多时相数据等优点,在松材线虫早期监测与受害林分识别等领域具有不错的应用前景。过去的研究中缺少对中分辨率卫星松材线虫病疫木监测能力的探究。而哨兵-2号和Landsat-8是免费卫星中分辨率最高、历史数据最丰富的多光谱卫星。目前利用哨兵-2号对松材线虫病进行监测的研究较少,有学者基于哨兵-2号卫星,利用辐射传播模型反演松树健康状态,进而识别受害松木(Lietal.,2022)。但是对哨兵-2号本身的波段信息及衍生信息的挖掘不够深入。基于哨兵-2号实现松材线虫病受害林分识别仍需进一步探究。同样作为中分辨率卫星的Landsat-8,具有免费、监测范围广和历史数据丰富等与哨兵-2号相同的特点,可以作为哨兵-2号的可靠对照。所以本试验的主要目的是:(1)比较Landsat-8与哨兵-2号卫星识别松材线虫病受害林分的能力;(2)比较不同机器学习模型与不同植被指数建立松材线虫病受害林分监测模型的准确率,筛选出最适合的植被指数和机器学习模型。

1 材料与方法

1.1 试验材料

1.1.1试验样地情况

试验地位于辽宁省东部的抚顺市章党镇大伙房林场(E124°13′06″,N41°56′15″)。研究区属于典型的温带大陆性季风气候,夏季高温多雨,冬季干旱少雨。年平均气温为6.8℃,年平均降水量为789.5 mm。林场内森林资源丰富,主要以人工种植的50年生红松PinuskoraiensisSiebold &Zuccarini、油松PinustabuliformisCarriere和华北落叶松Larixgmeliniivar.principis-rupprechtii(Mayr)Pilger为主。经过为期一年的实地调查(2020年7月-2021年10月),共选取了5块感染松材线虫病的纯林作为试验样地(图1),其中3块为红松纯林样地,2块为油松纯林样地。除红松2号样地受害较轻外(样地内松木感染率不超过20%),其余4块样地受害较为严重,林地内松木感染率均超过了50%。

图1 研究区样地图(辽宁抚顺)Fig.1 Study area(Fushun, Liaoning)

1.1.2无人机遥感数据

在2020年7-11月使用无人机采集受害红松纯林和油松纯林的可见光数据(图2)。无人机的型号为大疆精灵-4,飞行高度240 m,采集样地面积为50.7 ha。每个月进行2～3次无人机可见光数据采集,作为基于卫星数据识别受害松林的基础数据。

图2 无人机图像Fig.2 UAV images

1.1.3卫星遥感数据

选取哨兵-2号和Landsat-8作为数据源。其中哨兵-2号数据从欧洲航天局官方网站下载(https://www.esa.int/),拍摄日期为2020年 9月 30日,影像获取时天气状况良好,影像云层覆盖率小于10%,研究区内无云覆盖;Landsat-8影像数据从美国地质勘探局官方网站下载(https://www.usgs.gov/),拍摄日期为2020年10月3日,影像获取时天气状况良好,影像云层覆盖率小于10%,研究区无云覆盖。

1.2 试验方法

1.2.1受害木等级划分

在试验地中,容易遭受松材线虫病侵染的树种是油松和红松。通过贝尔曼漏斗法抽检变色的油松和红松样本均发现有松材线虫存在。而落叶松变色木数量低于5%且抽检结果没有松材线虫存在,故本研究选择3块红松纯林和2块油松纯林作为试验样地(图1)。

本研究中将树木的健康状况划分为了两个等级:(1)健康木,当枝梢变色率低于10%,判定为健康木;(2)感病木,当枝梢变色率大于等于10%,针叶褪绿变为黄色,即判定为感病木。

1.2.2无人机数据处理

在采集完成无人机数据之后,通过目视解译的方法提取出感染松材线虫病疫木的位置并制作感病木分布图。哨兵-2号和Landsat-8的多光谱传感器分辨率分别为20 m和30 m,而试验地中油松和红松的树冠直径平均值为5 m,所以同一像元内存在多棵松木。因此本研究将试验地的卫星影像像元划分为两种类型:(1)健康像元,感病木所占像元面积低于10%;(2)感病像元,感病木所占像元面积大于等于10%。以感病木分布图为基础,将试验地的卫星影像像元划分为健康像元和感病像元。

1.2.3卫星数据的预处理

卫星数据预处理的目的是为了将大气表层反射率转化为地表反射率并提取出研究区域内的森林植被。哨兵-2号和Landsat-8的预处理过程主要包括:(1)辐射定标;(2)大气校正;(3)研究区的裁剪;(4)林分掩膜。哨兵-2号的预处理过程在Snap中完成,Landsat-8的预处理过程在ENVI中完成,具体处理方法见(陈文静,2021)。

基于Snap中的Sen2Res插件可实现影像的超分辨率合成,能将哨兵-2号卫星影像的分辨率由20 m提升至10 m。在完成预处理后,使用Sen2Res插件生成分辨率为10 m的哨兵-2号影像。最终得到10 m、20 m的哨兵-2影像数据和30 m的Landsat-8影像数据。

1.2.4松材线虫病监测模型构建

松树在感染松材线虫病后,外部形态上表现为针叶褪色和枯萎,内部生理状态上表现为含水量和叶绿素含量降低(徐华潮,2013)。由于植株的含水量和叶绿素含量发生变化,会导致其光谱反射率产生变化,这一现象为松材线虫病疫木遥感监测奠定了理论基础。

本研究通过提取试验地像元的地理信息和特征参数(光谱反射率、植被指数),构建松材线虫病疫木监测模型。基于不同多光谱传感器的哨兵-2号和Landsat-8分别选取了50个和27个特征参数(表1)。在10 m影像中共提取了3 250个健康像元和1 287个感病像元的特征信息;在20 m影像中共提取了614个健康像元和520个感病像元的特征信息;在30 m影像中共提取了220个健康像元和300个感病像元的特征信息。

机器学习是一种常用的构建病虫害监测模型的方法(Durgabai and Bhargavi,2018)。本研究中分别对3种分辨率的影像数据(10 m、20 m的哨兵-2号和30 m的Landsat-8)采用4种机器学习算法(随机森林、决策树、支持向量机和极端梯度提升)进行拟合,最终构建出多种松材线虫病监测模型。具体模型拟合方法见Scikit-Learn网站(https://scikit-learn.org/)。

表1 特征参数名称

续表1 Continued table 1

续表1 Continued table 1

1.2.5松材线虫病监测敏感特征选择

特征选择是机器学习中常见的提升模型效能的方法。特征选择的意义有:(1)能有效减少模型运行时间;(2)避免特征的冗余、共线性和模型的过拟合等问题;(3)从不同角度筛选出对松材线虫病监测最为敏感的特征。

在特征参数筛选之前,为比较不同特征参数的贡献值,会对特征参数进行重要性排序。本研究选取了两种不同的重要性排序方法,分别是随机森林重要性排序与排列重要性(Permutation Importances)。随机森林的重要性排序结果会在模型拟合之后由模型自动生成。而排列重要性方法是输入不同的特征参数,通过多次迭代模型,得出重要性排序结果。具体重要性排序计算方法见Scikit-Learn网站(https://scikit-learn.org/)。

本研究只对10、20和30 m遥感数据中表现最优的数据集进行特征参数的筛选,因为试验的主要目的是通过筛选特征参数,进一步提升最优的模型性能。本研究中共采用了3种特征选择方法:卡方检验法、互信息法和递归消除特征法。3种特征选择方法都在Python中运行实现。对特征筛选后的数据集进行模型再拟合,最终得到不同筛选方法的模型优化能力。

1.2.6模型精度验证

为比较不同的监测模型性能,需要计算模型的性能度量指标。本研究中选取的性能度量指标为准确率(Accuracy)、Kappa系数(Kappa coefficient)、K折交叉验证(k-fold cross validation)和受试者特征曲线(Receiver Operation Characteristic,ROC)。计算方法:

(1)

TP是真阳性数量,TN是真阴性数量,FP是假阳性数量,FN是假阴性数量。

(2)

P0是总体分类精度,Pe是每一类的真实样本个数乘以预测样本个数的加和除以样总体样本数的平方。

本文中K折交叉验证的结果和ROC曲线图像通过Python计算得到,具体计算方法见Scikit-Learn网站(https://scikit-learn.org/)。

2 结果与分析

2.1 松材线虫病监测模型参数重要性

在不同的数据集中,相同特征参数重要性得分不同(图3、图4)。在随机森林重要性排序中,不同数据集的重要性得分最大值之间的差值均超过0.06,其中不同数据集的得分最高值由大到小排序为30、20、10 m,重要性得分越高表明其对模型的贡献值越大,说明得分最高的特征在30、20、10 m数据集中对模型拟合结果的影响是逐渐降低的。植被指数中只有NGRDI在3个数据集中重要性得分排名均为前2,说明其对松材线虫病疫木监测的能力强于研究中其余植被指数。光谱反射率中只有红波段在3个数据集中重要性排名均为前2,说明其对松材线虫病疫木监测的能力强于研究中其余光谱反射率。其余植被指数中NBRI在10 m、20 m影像中重要性得分高于0.06,但在30 m影像中低于0.04,且在10 m、20 m影像中NBRI重要性得分总和仅次于NGRDI,说明NBRI在10 m、20 m影像中表现最佳并对10 m、20 m影像拟合模型的重要性仅次于NGRDI。GLI和GNDVI与NBRI相反,在30 m影像中得分高于0.08而在10 m、20 m影像中得分低于0.04,说明GLI和GNDVI在30 m影像中表现最佳。而TVI、NDVI、RVI和PSSR等特征仅在10 m影像中得分高于0.06,说明TVI、NDVI、RVI和PSSR在10 m影像中表现最佳。排列重要性结果与随机森林重要性排序相似。也表现出不同数据集的相同特征参数重要性得分不同。不同数据集的得分最高值的差值均超过0.02,重要性的得分最大值由大到小排序也为30、20、10 m。NGRDI、NBRI、TVI、NDVI、RVI和PSSR等特征参数的排名在排列重要性排序中与在随机森林重要性排序中的结果相同(图4)。

2.2 松材线虫病监测模型性能

两种不同卫星的建模结果见图5,发现基于哨兵-2号影像数据建立的松材线虫病监测模型除20 m的决策树模型外准确率均高于基于30 m数据集拟合的模型。而基于哨兵-2号的10 m和20 m影像数据集,表现出在相同算法的条件下10 m监测模型准确率高于20 m监测模型。总的来说,在相同算法条件下,遥感影像的分辨率越高拟合出的模型效能也会越强。

图3 随机森林重要性得分结果Fig.3 Random Forest Importance scores result

图4 排列重要性得分结果Fig.4 Permutation Importance scores result

图5 建模结果Fig.5 Modeling results注:图中准确率均为模型在相同随机数条件下计算得出,没有进行迭代,为单次建模计算结果。Note: Accuracy rates in the figures were all calculated by the model under the same random number condition without iterations and were the results of a single modeling calculation.

4种不同的机器学习算法中,决策树算法在相同影像条件下拟合出的模型准确率均低于其余3种算法。在不同影像中,随机森林、支持向量机和极端梯度提升3种机器学习算法中的最优算法不同。在10 m影像中拟合模型准确率最高的算法是随机森林算法,而在20 m和30 m影像中拟合模型准确率最高的算法是支持向量机算法。

本研究通过计算ROC曲线中的AUC值来比较不同模型的拟合效果。曲线下面积(Area Under Curve,AUC)是ROC曲线下的面积。AUC值越大说明模型的拟合效果越好。图6中展示了10、20和30 m 三种影像构建模型的拟合效果,不同影像中随机森林、支持向量机和极端梯度提升等算法拟合的模型AUC值均大于0.80且差值均小于0.03,说明模型拟合效果接近且拟合结果优异,分类器预测结果可靠。决策树算法在3种不同分辨率影像中拟合的模型AUC值均低于0.80,说明决策树模型拟合效果较差。

松材线虫病监测模型是基于卫星影像数据与机器学习算法共同构成的,所以选择表现最佳的影像数据和机器学习算法能有效提升模型准确率。研究中发现在相同算法条件下基于哨兵-2号遥感影像拟合的模型准确率高于基于Landsat-8影像拟合的模型,而在10 m和20 m哨兵-2影像中,基于10 m哨兵-2号影像拟合的模型准确率最高,所以下一步特征筛选是基于10 m哨兵-2号影像实现的。从模型的准确率和ROC曲线综合评估4种机器学习算法构建模型的效能。除决策树外,在相同影像条件下随机森林、支持向量机与极端梯度提升等机器学习算法拟合的模型准确率差异不超过1%,因此在下一步再拟合过程中可采用任意机器学习算法。本研究中选取的再拟合机器学习算法为极端梯度提升算法。

2.3 松材线虫病监测敏感参数筛选

在特征筛选后利用机器学习算法再拟合模型(表2)。3种特征筛选方法减少的特征数量均少于16种且相互之间的差值小于3种,说明3种筛选方法在减少特征的能力上是相近的。在特征筛选后基于新的数据集进行再建模。根据2.2得到的结果,选取极端梯度提升作为再建模中采用的机器学习模型。比较基于不同筛选方法建立的模型准确率、Kappa系数、AUC值以及特征参数减少数量,发现只有递归消除特征法在降低特征数量的同时提升了模型准确率和Kappa系数。所以递归消除特征法是最优参数筛选方法。

图6 ROC曲线Fig.6 ROC curve

表2 不同特征筛选方法再建模结果

3 结论与讨论

在单时相卫星监测中,卫星空间分辨率越高,得到监测模型的准确率也越高。在过去的研究中,有许多学者比较不同卫星在同一病虫害监测中的准确率差异。发现在落叶松松毛虫受害区识别研究中(陈文静,2021)、山松甲虫受害木识别(Abdullahetal.,2018)、松材线虫病受害木识别中(Lee,2007)都表现出基于更高空间分辨率影像构建的模型准确率越高。卫星的分辨率还包括了时间分辨率与光谱分辨率。已有研究证明,更高的时间分辨率和光谱分辨率也能提升监测模型的准确率并能实现病害的早期监测(Abdullahetal.,2018)。本研究中仅利用单时相卫星监测受害松林,缺乏对哨兵-2号长时序特点的探究,也缺乏对松材线虫病早期感染阶段(受害林分外部形态尚未表现出感病特征)的探究(Chengetal.,2010)。通过对早期感染阶段的监测,进一步实现疫木早发现、早伐除的目标,最大程度的遏制松材线虫病的蔓延。有许多学者基于哨兵-2号长时序影像数据的特点,对山松甲虫的早期感染阶段展开研究(Bártaetal.,2021;Huoetal.,2021)。但是在松材线虫病疫木识别方向缺少类似研究。在未来的研究中,应该把基于多时相卫星监测松材线虫病早期感染阶段作为研究的重点。

红波段、绿波段、近红外波段和短波近红外波段等光谱反射率及其衍生植被指数在松材线虫病受害松木识别中具有重要意义。在过去的研究中发现,基于地面高光谱的松材线虫病疫木监测中绿波段、红波段、近红外波段和短波近红外波段的峰值变化与感病针叶光谱变化表现出强相关性(徐华潮,2013; Juetal.,2014)。植被指数NGRDI、NBRI等在过去的研究中展现出对松材线虫病疫木监测的敏感性。其中归一化绿红差异指数(NGRDI),对松木针叶颜色变化十分敏感,具有反应松叶变色的能力,对描述松木树冠颜色变化(Zhangetal.,2021)、调查松材线虫病发生区(国家林业和草原局,2022)起到了重要作用。而归一化燃烧率指数(NBRI)对含水量的减少十分敏感,由于受松材线虫感染后松树针叶表现为水分含量减少,该指数对受害松木监测模型构建具有高贡献值。在未来的松材线虫病受害松木研究中,应该加强对红波段、绿波段、近红外和短波近红外等光谱特征及其衍生光谱植被指数的研究,研发出更多对松材线虫病受害松木监测敏感的植被指数,从而提升遥感监测松材线虫病发生的准确率。

采用不同的信息提取方法、不同的遥感特征,以及将多种信息提取方法结合应用等都能实现模型准确率提升。现已有多种松材线虫病受害木信息提取的方法,包括机器学习(Takenakaetal.,2017)、深度学习(Zhouetal.,2022)、阈值分割(Lee,2007)、模型反演(Lietal.,2022)和长时序NDVI提取(Haoetal.,2021)等方法。而遥感特征则包括光谱特征、地形特征、物候特征和图像纹理特征等。为实现受害木监测模型准确率的提升通常会将多种不同的特征与不同的信息提取方法结合。所以在未来的研究中,为实现模型准确率的提升,应该着力于算法改进(图像增强算法、识别算法),加入更多的遥感特征数据和辅助数据(地形要素、树种分布),进一步提高松材线虫病受害林分识别准确率。