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

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

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

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

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

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

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

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

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

1.2 測(cè)報(bào)變量空間模擬

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

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

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

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

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

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

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

2.1 模型

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

2.2 精度評(píng)價(jià)

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

2.3 擴(kuò)散風(fēng)險(xiǎn)模擬

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

3 結(jié)果與分析

3.1 模型精度驗(yàn)證

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

表4 模型精度驗(yàn)證數(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 個(gè)模型在預(yù)測(cè)松材線蟲病縣級(jí)風(fēng)險(xiǎn)分布過(guò)程中具有邏輯可行性，預(yù)測(cè)結(jié)果具有可信性，但RF 的效果、模型穩(wěn)定性、模型分類結(jié)果泛化能力均高于SVM 模型。

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

3.2 變量重要性評(píng)價(jià)

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

圖2 RF 模型變量貢獻(xiàn)度Figure 2 Contribution of RF model variables

3.3 松材線蟲病縣域風(fēng)險(xiǎn)測(cè)報(bào)

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

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

3.4 全國(guó)柵格尺度松材線蟲病風(fēng)險(xiǎn)測(cè)報(bào)

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

表7 RF 與寄主易感性的松材線蟲病地理柵格潛在風(fēng)險(xiǎn)預(yù)測(cè)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)計(jì)方法，利用地理柵格信息描述驅(qū)動(dòng)變量，以RF 模型建模預(yù)報(bào)風(fēng)險(xiǎn)的方法，實(shí)現(xiàn)了測(cè)報(bào)變量的數(shù)字矩陣描述，解決了測(cè)報(bào)結(jié)果的空間分異性與連續(xù)化描述的難題。全國(guó)柵格尺度的松材線蟲病風(fēng)險(xiǎn)測(cè)報(bào)能夠更好地反映人類活動(dòng)、地理環(huán)境變化、寄主易感性對(duì)松材線蟲病傳播軌跡的影響，對(duì)基層單位開(kāi)展檢疫、預(yù)防和治理等工作有指示和支持作用。

降水、氣溫和海拔是影響松材線蟲病生存的主要因素，統(tǒng)計(jì)分析得出松材線蟲病最低生存年均降水量范圍為110～121 mm。影響松材線蟲病發(fā)生的年均低溫范圍為10～13 ℃，這與葉建仁[35]以年均氣溫10～14 ℃區(qū)域?yàn)樗刹木€蟲病可以發(fā)生區(qū)，大于14 ℃為流行區(qū)劃分的預(yù)測(cè)區(qū)域基本一致。影響松材線蟲病發(fā)生的海拔范圍為4～247 m，這與LEE 等[21]預(yù)測(cè)的松材線蟲病最適生存海拔小于200 m 基本一致。在傳播影響中，除了媒介昆蟲通過(guò)自然傳播，人類活動(dòng)對(duì)松材線蟲病的擴(kuò)散也具有重要的影響。松材線蟲病早期發(fā)生地點(diǎn)多位于城市化活躍、人類活動(dòng)強(qiáng)烈的區(qū)域內(nèi)，其原因?yàn)槌鞘杏玫亍⒌缆酚玫卦黾右约叭斯び痔娲贾脖唤档土松鷳B(tài)穩(wěn)定性[36]，這為松材線蟲病爆發(fā)與傳播埋下巨大隱患。

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

5 結(jié)論

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

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

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