基于多光譜衛星影像與機器學習算法的松材線蟲病受害林分識別研究

邱萬林,宗世祥

(北京林業大學林木有害生物防控北京市重點實驗室,北京 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)等方法。而遙感特征則包括光譜特征、地形特征、物候特征和圖像紋理特征等。為實現受害木監測模型準確率的提升通常會將多種不同的特征與不同的信息提取方法結合。所以在未來的研究中,為實現模型準確率的提升,應該著力于算法改進(圖像增強算法、識別算法),加入更多的遙感特征數據和輔助數據(地形要素、樹種分布),進一步提高松材線蟲病受害林分識別準確率。