基于地理柵格變量與機器學(xué)習(xí)的松材線蟲病擴散風(fēng)險分析

劉婷婷，楊晉帆，周汝良，劉 琳

（1.西南林業(yè)大學(xué) 地理與生態(tài)旅游學(xué)院，云南 昆明 650224；2.西南林業(yè)大學(xué) 理學(xué)院，云南 昆明 650224）

由松材線蟲Bursaphelenchusxylophilus引起的松材線蟲病是目前松樹上最具破壞性的病害之一。自1905年日本首次報道發(fā)生松材線蟲病以來，亞歐其他國家相繼出現(xiàn)松材線蟲病疫區(qū)[1-2]。中國的松材線蟲病最早報道于1982年，出現(xiàn)在南京中山陵的黑松Pinusthunbergii上，之后該病害快速擴散，據(jù)國家林業(yè)和草原局統(tǒng)計，中國現(xiàn)共有松材線蟲病縣級疫區(qū)726 個、鄉(xiāng)鎮(zhèn)級疫點5 479 個，發(fā)生面積達180.9 萬hm2，包括輕型疫區(qū)206 個(含2021年初新增疫情發(fā)生區(qū))、重型疫區(qū)518 個[3-4]。目前松材線蟲病已成為中國最為嚴(yán)重的林業(yè)病害之一[5-6]。為實現(xiàn)《全國松材線蟲病疫情防控五年攻堅行動計劃(2021—2025)》中新發(fā)疫情“早發(fā)現(xiàn)、早報告、早除治、早拔除”的目的，達到控制增量的總體目標(biāo)，就需要把握疫情的總體趨勢，精準(zhǔn)定位和定量各疫區(qū)的疫情，提高防疫針對性[7]。

已有研究發(fā)現(xiàn)：全球變暖環(huán)境下，氣溫增高會影響松材線蟲病的地理分布[8]；人類活動是影響松材線蟲病遠(yuǎn)距離、跳躍式、梅花式擴散傳播的重要因素[9-11]；年均氣溫、年均降水量等是影響松材線蟲定殖的主要氣候因素[12]；中國松材線蟲病的發(fā)生受到寄主植物分布的強烈影響[13-14]；媒介昆蟲對寄主的選擇性也影響松材線蟲病的發(fā)生[15-17]。當(dāng)前對松材線蟲病傳播風(fēng)險的研究多集中于采用傳統(tǒng)的地統(tǒng)計學(xué)方法來分析松材線蟲病在省、市范圍的風(fēng)險格局[18-20]，而將機器學(xué)習(xí)引入病害風(fēng)險測報研究較少見。機器學(xué)習(xí)方法中隨機森林、支持向量機具有需要樣本量少，能夠同時利用多變量進行非線性回歸、分類預(yù)測并且結(jié)果精度高等優(yōu)點，使風(fēng)險預(yù)測精確化分析成為可能[21-23]。

本研究綜合考慮松材線蟲病擴散的影響因素，利用地理環(huán)境、氣候環(huán)境、衛(wèi)星遙感、道路交通網(wǎng)絡(luò)等構(gòu)成地理柵格數(shù)據(jù)集。對比隨機森林和支持向量機2 種算法建立松材線蟲病風(fēng)險擴散模型，探究機器學(xué)習(xí)算法在松材線蟲病擴散傳播預(yù)測研究中的可行性和準(zhǔn)確性，尋求最佳預(yù)測方案，將潛在疫區(qū)根據(jù)等級進行風(fēng)險區(qū)劃分，并進行柵格化分析，將全國松材線蟲病測報落實到“山頭地塊”上，以期為林業(yè)生態(tài)系統(tǒng)中松科Pinaceae 植物的動態(tài)監(jiān)測提供有效的手段。

1 數(shù)據(jù)與方法

1.1 疫區(qū)樣本數(shù)據(jù)

以中國縣區(qū)為行政單元，以國家林業(yè)和草原局2018—2021年關(guān)于松材線蟲病疫區(qū)公告為基礎(chǔ)[4]，收集整理了中國松材線蟲病發(fā)生區(qū)731 個以及未發(fā)生區(qū)2 131 個，以松材線蟲病發(fā)生(1)或未發(fā)生(0)作為模型響應(yīng)變量，由點及面構(gòu)建了2 862 個樣點。其中隨機選取30%樣本點(954 個)作為模型獨立檢驗樣本?？倶颖军c分布見表1。

表1 松材線蟲病疫區(qū)Table 1 Pine wilt disease endemic area

1.2 測報變量空間模擬

1.2.1 氣候環(huán)境數(shù)據(jù) 根據(jù)以往研究[8,12]，選取1970—2000年共30 a 逐月氣象數(shù)據(jù)獲得年均氣溫、年均降水量以及年均最高氣溫、年均最低氣溫、年均最低降水量與年均最高降水量數(shù)據(jù)?？臻g分辨率為 250 m的氣候數(shù)據(jù)集來源于https://www.worldclim.org/。

1.2.2 傳播媒介數(shù)據(jù) 選取影響昆蟲媒介墨天牛屬Monochamus的地理因子，包括高程數(shù)據(jù)，以及根據(jù)高程數(shù)據(jù)計算提取的坡度與坡向指數(shù)。選取海拔、坡度、植被覆蓋度以及溫度構(gòu)建阻力面，表征其自然傳播的能力。根據(jù)全國高速公路、國道、省道以及一級水系等交通網(wǎng)絡(luò)，利用歐氏距離來模擬人類活動導(dǎo)致的有害生物空間傳播度量[9]。道路、水系矢量數(shù)據(jù)等人為影響數(shù)據(jù)來源于Bigemap 平臺，利用ArcGIS 處理形成綜合交通便利性數(shù)據(jù)；土地利用類型數(shù)據(jù)來源于GLOBELAND 30 平臺。

1.2.3 寄主植物數(shù)據(jù) 全國1∶100 萬植被分布數(shù)據(jù)來源于資源環(huán)境數(shù)據(jù)云平臺http://www.Resdc.cn/，從中提取出易受松材線蟲病感染的針葉林植被。根據(jù)前人研究結(jié)果[24]將不同樹種的寄主易感指數(shù)賦分為0～100，數(shù)值越高表示寄主易感性越高，劃分標(biāo)準(zhǔn)見表2。

表2 松材線蟲病寄主量化Table 2 Quantification of pine wilt disease hosts

綜上，選擇16 個環(huán)境變量因子作為主要環(huán)境數(shù)據(jù)(表3)。環(huán)境數(shù)據(jù)集的空間分辨率通過重采樣統(tǒng)一為250 m。地理數(shù)據(jù)為1∶320 萬的中國國界、省界以及縣界行政區(qū)劃圖，來自于自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)http://bzdt.ch.mnr.gov.cn/。

表3 松材線蟲病預(yù)測模型自變量表Table 3 Independent variables of the prediction model of pine wilt disease

2 擴散風(fēng)險模型構(gòu)建

2.1 模型

隨機森林(Random Forest Classifier，RF)是一種監(jiān)督學(xué)習(xí)的機器學(xué)習(xí)算法，即將多個決策樹分類器集合在一起進行數(shù)據(jù)的分類和回歸[25-27]。模型從原始樣本抽取的隨機性，以及利用小于樣本集特征的數(shù)據(jù)建立訓(xùn)練集，增加了不同決策樹結(jié)果的差異性，中和了過擬合現(xiàn)象，使得模型的泛化能力提高使其具有較高的分類準(zhǔn)確性[28-29]。支持向量機(Support Vector Machine，SVM)算法能夠解決機器學(xué)習(xí)中樣本量小、非線性和高維模式識別中“維數(shù)災(zāi)難”和“過學(xué)習(xí)”等方面的問題[30-31]。該方法選擇適當(dāng)?shù)暮撕瘮?shù)通過非線性變換將輸入空間變換到一個高維空間，然后在這個高維空間求取最優(yōu)分類面，找到輸入變量和輸出變量之間的一種非線性關(guān)系。

2.2 精度評價

ROC-AUC 特征曲線分析方法用來評價二值分類器的優(yōu)劣[32]。ROC 曲線反映了輸出概率的準(zhǔn)確性，曲線上每個點反映對同一信號刺激的感受性。橫軸為負(fù)正類率(false postive rate，F(xiàn)PR)特異度，縱軸為真正類率(true postive rate，TPR)，靈敏度均為[0,1]。ROC 曲線在斜對角線以下，表示該分類器效果差于隨機分類器，反之，效果好于隨機分類器[33]。ROC 曲線下的面積(area under curve，AUC)作為數(shù)值可以直觀地評價分類器的好壞。

2.3 擴散風(fēng)險模擬

RF 和SVM 模型均以松材線蟲病在全國縣域尺度中存在(1)或不存在(0)作為響應(yīng)變量，對所有環(huán)境變量在輸入前進行歸一化處理，以消除由于綱量導(dǎo)致的差異。數(shù)據(jù)集以7∶3 的比例隨機分成兩部分進行模型訓(xùn)練和測試。在建模過程中16 個環(huán)境變量根據(jù)模型重要性及其影響，篩除部分低貢獻度變量，并將針葉林分布變量提出作為重要變量與預(yù)測風(fēng)險區(qū)疊加分析。模型預(yù)測精度和ROC-AUC 曲線用來評價每個模型的性能，利用預(yù)測結(jié)果與獨立檢驗樣本通過混淆矩陣進行精確性檢驗。通過Python 轉(zhuǎn)為0～1 之間的發(fā)生概率，輸入ArcGIS 中進行風(fēng)險可視化處理。除此之外，全國柵格尺度更能夠體現(xiàn)空間連續(xù)變化，在進行運算操作、疊加分析方面也比縣級矢量數(shù)據(jù)更具優(yōu)勢[34]。利用從全國1∶100 萬植被分布類型中提取出來的寄主松科植物，根據(jù)松科植物不同的易感性重分類，柵格化后與模型預(yù)測的松材線蟲病潛在發(fā)生區(qū)域疊加分析。

3 結(jié)果與分析

3.1 模型精度驗證

通過模型精度和ROC-AUC 檢驗結(jié)果(表4)可知：2 個模型均有較好的預(yù)測能力，其預(yù)測精度和模型穩(wěn)定性均在75%以上。模型精度是通過RF 和SVM 模型對獨立檢驗樣本數(shù)據(jù)進行松材線蟲病的發(fā)生風(fēng)險概率值預(yù)測，設(shè)定臨界值為0.5，將大于此值縣域視為松材線蟲病“存在”(y=1)，小于此值歸為“不存在”(y=0)，得出分類結(jié)果與獨立檢驗樣本的實際值來驗證模型預(yù)測的精度。結(jié)果表明：構(gòu)建模型的測試集預(yù)測精度在RF 和SVM 模型中總精度分別為83.95%和77.97%。

表4 模型精度驗證數(shù)據(jù)表Table 4 Model accuracy verification data sheet

RF 和SVM 模型的ROC-AUC 曲線值(圖1)分別為0.89 和0.77 (平均標(biāo)準(zhǔn)誤差為0.01)。這證明2 個模型在預(yù)測松材線蟲病縣級風(fēng)險分布過程中具有邏輯可行性，預(yù)測結(jié)果具有可信性，但RF 的效果、模型穩(wěn)定性、模型分類結(jié)果泛化能力均高于SVM 模型。

圖1 ROC-AUC 曲線Figure 1 ROC-AUC curve

3.2 變量重要性評價

根據(jù)RF 模型對數(shù)據(jù)貢獻度(圖2)結(jié)果得出，重要性排序中位于前列的為氣候變量，包括年均最低降水量(0.302 59)、年均降水量(0.170 33)、年均低溫(0.102 98)；其次是海拔變量(0.150 84)，這2 類是影響松材線蟲病的重要環(huán)境變量。松材線蟲病發(fā)生概率在低海拔地區(qū)較高，隨著年均低溫的升高風(fēng)險概率增加，而高溫、干旱、低濕氣候直接影響松材線蟲病的爆發(fā)。在傳播媒介影響中，人為因素中各級交通網(wǎng)絡(luò)影響力重要性更為突出(0.067 58)，說明密集的人流和物流是松材線蟲病遠(yuǎn)距離傳播的重要媒介。在重要性排序中松科植物類型數(shù)據(jù)貢獻度與實際影響松材線蟲病的情況不符，寄主是影響松材線蟲病存在和發(fā)展的必然條件，故將針葉林分布單獨提出柵格化，與預(yù)測風(fēng)險區(qū)疊加分析其影響。

圖2 RF 模型變量貢獻度Figure 2 Contribution of RF model variables

3.3 松材線蟲病縣域風(fēng)險測報

利用Python 將預(yù)測結(jié)果提取為0～1 之間的概率數(shù)據(jù)，依據(jù)離散程度劃分為5 個等級，由高到低依次為極高風(fēng)險、高風(fēng)險、中等風(fēng)險、低風(fēng)險、極低風(fēng)險。因樣本原因，不涉及中國香港、臺灣、澳門行政區(qū)。對比RF 與SVM 模型的分析結(jié)果(表5和表6)，SVM 對重要變量的模糊表達導(dǎo)致其預(yù)測結(jié)果等級差距小、重點疫區(qū)不突出、預(yù)測精確性低、全國大部分地區(qū)呈同一等級，這使得對松材線蟲病疫區(qū)的監(jiān)測與防治難度增大。相比而言，RF 預(yù)測結(jié)果具有更高的準(zhǔn)確性，其潛在疫區(qū)分布與原有疫區(qū)重合度高，危險等級一致，并能夠明顯表達城市、道路以及地形對潛在疫區(qū)分布的影響。

表6 SVM 松材線蟲病縣域潛在風(fēng)險等級Table 6 SVM potential risk levels of pine wilt disease in counties

3.4 全國柵格尺度松材線蟲病風(fēng)險測報

寄主植物根據(jù)其易感性劃分等級，松材線蟲病易感性松科植物多分布于浙江、福建、湖南、廣西等省，另外，云南、貴州、四川以及東北林區(qū)也有大量寄主植物分布。利用RF 模型疊加寄主易感性分布數(shù)據(jù)對全國31 個省進行地理柵格尺度的疫區(qū)劃分與分析，其中概率大于0.8 以上的柵格點為極易感染區(qū)域。結(jié)果(表7)表明：松材線蟲病的極高風(fēng)險區(qū)為華東地區(qū)的浙江、江西、福建，華南地區(qū)的廣西、廣東以及華中地區(qū)的湖南，高易感性寄主分布與極高風(fēng)險區(qū)預(yù)測結(jié)果高度重合，疫區(qū)內(nèi)多為赤松、黑松、馬尾松；高風(fēng)險區(qū)包括華東地區(qū)的安徽，華中地區(qū)的湖北，以及西南地區(qū)的重慶、貴州，高風(fēng)險區(qū)寄主也呈高易感性，但其分布分散且多于其他植被混生，森林具有類型多樣化；中等風(fēng)險區(qū)包括華中的河南、江蘇，東北地區(qū)遼寧，華東地區(qū)山東，寄主植物易感性降低，多呈中等易感性且面積小，多為散布狀態(tài)，其中很多省份為糧食大省，以種植農(nóng)作物為主，林業(yè)分布較少；無風(fēng)險區(qū)包括西北地區(qū)甘肅、寧夏、青海、新疆，華北地區(qū)的內(nèi)蒙古、山西和西南地區(qū)的西藏，缺少寄主植物的分布，人流、物流相較于其他地區(qū)也較少，難以形成松材線蟲病潛在風(fēng)險區(qū)；其余省份劃分為低風(fēng)險區(qū)，其中需注意的是云南、四川大部分地區(qū)存在高易感性的思茅松、云南松等優(yōu)勢種，風(fēng)險等級容易上升。

表7 RF 與寄主易感性的松材線蟲病地理柵格潛在風(fēng)險預(yù)測Table 7 Prediction of the potential risk of geographic grids for RF and host susceptibility to pine wilt disease

4 討論

相比于縣域行政單元尺度的傳統(tǒng)地學(xué)統(tǒng)計方法，利用地理柵格信息描述驅(qū)動變量，以RF 模型建模預(yù)報風(fēng)險的方法，實現(xiàn)了測報變量的數(shù)字矩陣描述，解決了測報結(jié)果的空間分異性與連續(xù)化描述的難題。全國柵格尺度的松材線蟲病風(fēng)險測報能夠更好地反映人類活動、地理環(huán)境變化、寄主易感性對松材線蟲病傳播軌跡的影響，對基層單位開展檢疫、預(yù)防和治理等工作有指示和支持作用。

降水、氣溫和海拔是影響松材線蟲病生存的主要因素，統(tǒng)計分析得出松材線蟲病最低生存年均降水量范圍為110～121 mm。影響松材線蟲病發(fā)生的年均低溫范圍為10～13 ℃，這與葉建仁[35]以年均氣溫10～14 ℃區(qū)域為松材線蟲病可以發(fā)生區(qū)，大于14 ℃為流行區(qū)劃分的預(yù)測區(qū)域基本一致。影響松材線蟲病發(fā)生的海拔范圍為4～247 m，這與LEE 等[21]預(yù)測的松材線蟲病最適生存海拔小于200 m 基本一致。在傳播影響中，除了媒介昆蟲通過自然傳播，人類活動對松材線蟲病的擴散也具有重要的影響。松材線蟲病早期發(fā)生地點多位于城市化活躍、人類活動強烈的區(qū)域內(nèi)，其原因為城市用地、道路用地增加以及人工育林替代原始植被降低了生態(tài)穩(wěn)定性[36]，這為松材線蟲病爆發(fā)與傳播埋下巨大隱患。

全國地理柵格擴散風(fēng)險格局分析表明：松材線蟲病在華北平原和長江中下游平原持續(xù)擴散，并且存在3 條傳播廊道：一是沿著長江水系向西南地區(qū)及四川盆地擴散，現(xiàn)以水熱條件優(yōu)越、高城市化的重慶潛在風(fēng)險等級為最高；二是沿著黃渤海臨海城市向東北平原的擴散廊道，這個廊道中北京和天津的潛在風(fēng)險等級突出，可能是受到高城市化發(fā)展和密集交通網(wǎng)絡(luò)的影響；三是沿金沙江、珠江流域形成了一條由廣東、海南、廣西向云貴高原擴散的通道，雖西南地區(qū)海拔較高，但其豐富的松科植物仍為松材線蟲病的生存和傳播提供基礎(chǔ)。

5 結(jié)論

基于多類型數(shù)據(jù)建立的RF 與SVM 模型均有較好的預(yù)測性能(ROC-AUC＞75%)。RF 模型具有更高的精度為83.95%，并且預(yù)測的極高風(fēng)險區(qū)與寄主高易感性分布高度重合，證明RF 模型分類結(jié)果更符合實際。

影響松材線蟲病發(fā)生的敏感變量是年均降水量、年均低溫和海拔，而坡度小、人類活動強度高是導(dǎo)致松材線蟲病快速傳播的重要因素。潛在擴散區(qū)位于人類活動密集的低海拔地區(qū)、道路通達的林區(qū)、城市城鎮(zhèn)分布區(qū)和人工林分布區(qū)。其中，極高風(fēng)險分布地區(qū)主要位于華東地區(qū)的浙江、江西及福建，華南地區(qū)的廣西、廣東以及華中地區(qū)的湖南。松材線蟲擴散過程中存在明顯傳播廊道與地理阻隔。松材線蟲在華北平原的擴散傳播受到秦嶺和太行山脈的阻隔，并未侵入西北地區(qū)的黃土高原，而高風(fēng)險疫區(qū)多集中在黃河水系以南地區(qū)。

本研究僅使用了RF 及SVM 共2 種機器學(xué)習(xí)算法，樣點數(shù)據(jù)僅為國家林業(yè)和草原局2018—2021年松材線蟲病疫區(qū)數(shù)據(jù)，也還未考慮環(huán)境動態(tài)變化導(dǎo)致的松材線蟲病疫區(qū)發(fā)生的風(fēng)險變化。未來的研究應(yīng)以松材線蟲病實際病害發(fā)生位置為基礎(chǔ)，通過更多的機器學(xué)習(xí)算法或深度學(xué)習(xí)算法進行對比研究，充分考慮與研究區(qū)相適應(yīng)的自然環(huán)境、人類活動等時空數(shù)據(jù)，以期為松材線蟲病疫區(qū)的預(yù)測提供更加準(zhǔn)確、精細(xì)的風(fēng)險測報。