基于MaxEnt模型預測黃脊竹蝗在中國的適生區(qū)

溫玄燁，王 越，姜 璠，唐 健，林 曉

(國家林業(yè)和草原局森林和草原病蟲害防治總站，林業(yè)有害生物監(jiān)測預警國家林業(yè)和草原局重點實驗室，沈陽 110034)

黃脊竹蝗CeracriskiangsuTsai，又稱竹蝗，屬直翅目Orthoptera網(wǎng)翅蝗科Arcypteridae，以跳蝻和成蟲取食葉及嫩梢，主要危害剛竹屬Phyllostachys、箣竹屬Bambusa竹種，也可危害玉米Zeamays、水稻Oryzasativa、棕櫚Trachycarpusfortunei等農(nóng)作物及近百種雜草(鐘洪武等, 2010)。黃脊竹蝗食性雜、取食量大、危害時間長、防治困難，對竹子產(chǎn)業(yè)危害嚴重，在中國南部、越南、老撾等熱帶和亞熱帶地區(qū)均有分布(伍遠平等, 2005; Yuetal., 2011; 孔志堅和寸佳蒞, 2018)。

黃脊竹蝗是我國記錄最早的林業(yè)害蟲之一，早在明朝嘉慶年間就曾大規(guī)模發(fā)生(程佳等, 2010)。上世紀50年代，黃脊竹蝗曾在我國南方地區(qū)暴發(fā)，造成了嚴重的經(jīng)濟損失。近年來，我國經(jīng)濟林面積不斷擴大，黃脊竹蝗危害又有抬頭之勢。2020年6月以來，云南省中老邊境發(fā)生境外黃脊竹蝗遷飛入境事件，對相關(guān)地區(qū)的農(nóng)林業(yè)生產(chǎn)造成危害，引起了農(nóng)業(yè)農(nóng)村部和國家林業(yè)和草原局的高度重視。我國學者對黃脊竹蝗開展了大量研究，在生物學和防治技術(shù)等領(lǐng)域取得了豐富的成果和進展(程佳等, 2009; 張威等, 2016; 高尚等, 2020)，但關(guān)于其在我國的潛在分布區(qū)研究仍屬空白。因此，開展黃脊竹蝗在我國的潛在分布區(qū)預測及適生性分析研究，對做好災害預防和風險評估有重要意義。

生態(tài)位模型是通過物種分布及環(huán)境變量，利用數(shù)學模型描述物種生態(tài)需求并進行空間投影，分析物種適生性的方法，目前已應用于珍稀植物分布、外來物種入侵、野生動物保護等領(lǐng)域，均取得較好效果(遲翔文等, 2009; 李望軍等, 2019; 馮春慧等, 2020)。最大熵(MaxEnt)模型通過已知空間分布，尋找關(guān)鍵環(huán)境變量，構(gòu)建約束集合，模擬二者關(guān)系的生態(tài)位模型，其優(yōu)點是預測結(jié)果準確、需求樣本數(shù)量少，可較好地預測物種的適生區(qū)，國內(nèi)現(xiàn)已基于MaxEnt模型對紅脂大小蠹Dendroctonusvalens(王濤等, 2018)、美國白蛾Hyphantriacuneain(李濤等, 2018)、松針紅斑病Dothistromapini(王曉偉等, 2019)等林業(yè)有害生物適生區(qū)進行預測，均得到了較好的效果，準確地反映了其潛在適生區(qū)，為林業(yè)有害生物的防控提供了依據(jù)。

本研究以黃脊竹蝗為研究對象，基于本次云南黃脊竹蝗遷飛入境點位和我國已報道的物種分布信息，在當前氣候條件下，通過MaxEnt模型結(jié)合ArcGIS對黃脊竹蝗在我國的適生區(qū)進行預測，分析影響其生存的關(guān)鍵環(huán)境變量，以期為黃脊竹蝗的監(jiān)測預警和綜合防治提供理論依據(jù)。

1 材料與方法

1.1 材料

分布數(shù)據(jù)：黃脊竹蝗的分布數(shù)據(jù)來源于實地調(diào)查、檢索中國林業(yè)有害生物防治信息管理系統(tǒng)(http://os.forestpest.cn/)和相關(guān)文獻資料(黃煥華等, 2003; 陳良昌等, 2013; 方蓉等, 2015; 楊華等, 2018; 宋玉雙, 2019)。

環(huán)境變量：從世界氣象WorldClim Version 2.0數(shù)據(jù)庫(http://www.worldclim.org)中下載1970-2000年19個柵格大小為25 km2生物氣候環(huán)境變量數(shù)據(jù)Bio1～Bio19，分別為年均氣溫、晝夜溫差月均值、等溫性、溫度變化方差、最熱月份最高溫、最冷月份最低溫、年溫變化范圍、最濕季度平均溫度、最干季度平均溫度、最暖季度平均溫度、最冷季度平均溫度、年降水量、最濕月份降水量、最干月份降水量、降水量變化方差、最濕季度降水量、最干季度降水量、最暖季度降水量、最冷季度降水量。

地圖地形：來源于國家自然資源部標準地圖服務系統(tǒng)(http://bzdt.ch.mnr.gov.cn/)下載的1 ∶1 600萬標準中國地圖，審圖號是GS(2016)2923號，通過ArcGIS空間分析工具獲取3個地形因子(海拔、坡向和坡度)。

1.2 方法

1.2.1數(shù)據(jù)處理

通過數(shù)據(jù)庫或文獻檢索到的物種分布記錄，部分數(shù)據(jù)缺失經(jīng)緯度坐標，需要通過Google Earth v7.1.8輸入地名查詢。為避免地理空間采樣過于密集導致生態(tài)空間分布點重復(朱耿平等, 2014)，將所得數(shù)據(jù)利用分辨率為2.5 arc-min的柵格去除重復，降低采樣點的空間偏差。最終獲得物種分布點合計267個，其中，侵入點25個，已發(fā)生數(shù)據(jù)242個(圖1)，根據(jù)MaxEnt軟件要求，將所得數(shù)據(jù)整理成“.CSV”格式備用。

1.2.2主要環(huán)境變量的篩選

由于黃脊竹蝗種群的分布與存活與環(huán)境變量存在相關(guān)性，為減少環(huán)境因子的多重貢獻性對模型的影響，避免空間相關(guān)性對建模產(chǎn)生誤差，利用SPSS 14.0軟件對19個環(huán)境因子進行Pearson共線性檢驗。當2個環(huán)境變量相對系數(shù)≥0.8時，表明2個變量間存在共性關(guān)系，選取貢獻較高的環(huán)境變量，當2個環(huán)境變量相對系數(shù)<0.8時，表明變量間不存在共性關(guān)系，保留該變量用于模型分析。

篩選得到等溫性(bio3)、最熱月最高溫(bio5)、最冷月最低溫(bio6)、溫度年較差(bio7)、最干月降水量(bio14)、降水量季節(jié)性變異系數(shù)(bio15)、最濕季降水量(bio16)等7個環(huán)境變量(表1)和3個地形因子海拔(Altitude)、坡向(Aspect)和坡度(Slope)用于進一步分析。

圖1 黃脊竹蝗分布點位Fig.1 Distribution point of Ceracris kiangsu注：該圖基于自然資源部標準地圖服務系統(tǒng)下載的審圖號為GS(2016)2923號的標準地圖制作，底圖無修改。Note: The drawing was based on the standard map with drawing No. GS (2016) 2923 downloaded by the standard map service system of the Ministry of Natural Resources, and the bottom map was not modified.

表1 環(huán)境變量的相關(guān)性矩陣

1.2.3參數(shù)設置

將黃脊竹蝗地理分布點數(shù)據(jù)和主要環(huán)境因子數(shù)據(jù)導入MaxEnt軟件，選擇隨機25%數(shù)據(jù)坐標點作為測試數(shù)據(jù)，75%坐標點為訓練數(shù)據(jù)，在默認收斂域限10-5,最大迭代次數(shù)500狀態(tài)下，進行10次重復，選取平均值作為黃脊竹蝗分布預測結(jié)果，輸出結(jié)果以ASC文件形式顯示(Peterson and Eaton, 2010)。

1.2.4黃脊竹蝗適生區(qū)預測與評價

將MaxEnt模型得到的黃脊竹蝗潛在分布結(jié)果輸入ArcGIS，轉(zhuǎn)化為柵格形式(趙健等, 2019)。按照模型生成的適生概率(P)，利用自然間斷點分級法(Jenks)，結(jié)合黃脊竹蝗在發(fā)生區(qū)所造成的危害嚴重程度，將黃脊竹蝗的適生區(qū)分為4級：非適生區(qū)(0≤P<0.05)，低適生區(qū)(0.05≤P<0.25)，中適生區(qū)(0.25≤P<0.5)，高適生區(qū)(0.5≤P<1)，利用ArcGIS中重分類功能獲取黃脊竹蝗在我國的適生區(qū)等級分布圖(謝岷等, 2019)。

1.2.5環(huán)境變量對黃脊竹蝗生存率的影響

采用Jackknife(刀切法)探討環(huán)境變量和地形因子對黃脊竹蝗分布區(qū)域和生存概率的影響，通過1.2.2獲得10個環(huán)境變量(包括地形因子的3個變量)，依次去除1個環(huán)境變量，利用剩余變量建模，分析去除該變量后的相關(guān)性，從而確定10個環(huán)境變量與黃脊竹蝗生存的相關(guān)程度，即生存率。

1.2.6模型評價

通過ROC(受試者工作)曲線對模型精度進行評價。以特異度(假陽性率)為X軸，靈敏度(遺漏率)為Y軸，計算ROC曲線下積值即AUC值對模型的精度進行評價，當AUC值<0.7時，預測結(jié)果準確性較差，當0.7≤AUC值<0.9時，預測結(jié)果準確性一般，當AUC>0.9時，預測結(jié)果準確性高(Merowetal., 2013)。

2 結(jié)果與分析

2.1 模型精度評價

對MaxEnt模型進行準確性檢驗，得出訓練集和測試集10次重復結(jié)果(圖2)。AUC平均值為0.958，平均標準誤差為0.005，表明此模型預測準確度高，利用現(xiàn)有分布點數(shù)據(jù)和環(huán)境變量可有效預測黃脊竹蝗在我國的潛在分布范圍。

圖2 黃脊竹蝗適生區(qū)的ROC驗證曲線Fig.2 ROC curve of potential distribution of Ceracris kiangsu

2.2 黃脊竹蝗在我國的適生區(qū)

MaxEnt模型預測結(jié)果顯示，黃脊竹蝗在我國江淮流域、長江中下游地區(qū)、華南及西南等地有不同程度的適生區(qū)域，其適生區(qū)面積占我國陸地總面積的18.9%，其中，高適生區(qū)占3.0%，中適生區(qū)占5.6%，低適生區(qū)占10.3%。高適生區(qū)主要分布在長江中下游的江西、湖南、浙江、安徽，華南沿海的廣東、廣西，以及西南的云南等省區(qū)。以湘贛地區(qū)的鄱陽湖平原、兩湖平原為中心向四周發(fā)散式分布，南至兩廣丘陵北部，東至浙西丘陵地帶，北至鄂皖地區(qū)的大別山南麓，向西跳躍式分布于云南中老邊境及重慶中部地區(qū)。中適生區(qū)主要分布于安徽中部、湖北西北部、湖南西部、廣西大部、重慶西北部及云南西南部地區(qū)，以及上述高適生區(qū)的周邊區(qū)域。低適生區(qū)廣泛分布于我國淮河以南各主要省區(qū)，向南可至海南省北部及臺灣省南部地區(qū)(圖3)。

2.3 環(huán)境變量對黃脊竹蝗分布的影響

應用刀切法對環(huán)境變量進行分析，結(jié)果顯示，最干月降水量(bio14)和最冷月最低溫(bio6)貢獻率分別51.5%和22.3%，正規(guī)化訓練增益值>1.4，對模型提供了最大的效益，說明使用最干月降水量和最冷月最低溫時比其他變量包含更多有用信息，對黃脊竹蝗分布的影響最大。最熱月最高溫(bio5)、最濕季降水量(bio16)、降水量季節(jié)性變異系數(shù)(bio15)、溫度年較差(bio7)正規(guī)化訓練增益值>1.0，四個環(huán)節(jié)變量累計貢獻率為22.1%，說明對黃脊竹蝗的分布產(chǎn)生較大影響。而等溫性(bio3)、海拔(altitude)、坡向(aspect)和坡度(slope)等環(huán)境變量對黃脊竹蝗的分布影響相對較小(圖4)。

圖3 黃脊竹蝗在我國的適生區(qū)預測分布圖Fig.3 Predicted potential suitable habitats of Ceracris kiangsu注：該圖基于自然資源部標準地圖服務系統(tǒng)下載的審圖號為GS(2016)2923號的標準地圖制作，底圖無修改。Note: The drawing was based on the standard map with drawing No. GS (2016) 2923 downloaded by the standard map service system of the Ministry of natural resources, and the bottom map was not modified.

圖4 刀切法分析環(huán)境變量的重要程度Fig.4 Jackknife test for evaluating the relative importance of environmental variables

影響黃脊竹蝗分布的主要環(huán)境變量是降水量和溫度，從降水量來看，當最干月降水量不足50 mm時，黃脊竹蝗的生存率不足0.4，當降水量超過200 mm時，生存率上升至0.66(圖5-A)，按其適生區(qū)等級劃分方法，以P=0.05為閾值，得到降水量變化范圍為7～2 619 mm(圖5-B)，最低適宜值為72.16 mm。從溫度來看，最冷月最低溫極值為1℃(圖5-C)，低于或高于該值時，黃脊竹蝗生存率降低，最熱月最高溫度不足25℃時，黃脊竹蝗生存率不足0.1，隨著溫度上升其生存率也逐漸上升，33.0℃以上時，黃脊竹蝗的生存率大于0.5。

圖5 不同環(huán)境變量與黃脊竹蝗生存概率間的關(guān)系Fig.5 Relationship between different environmental variables and survival probability of Ceracris kiangsu注：圖層中的溫度數(shù)值(℃)為實際數(shù)值×10。Note: Temperature value (℃) in the layer was the actual value×10.

3 結(jié)論與討論

黃脊竹蝗屬寄主主導型食葉害蟲，其生長發(fā)育與寄主種類密切相關(guān)(張守科等, 2017)，在景觀尺度上，昆蟲的適生區(qū)范圍與其寄主分布直接相關(guān)(王夢琳等, 2017)。我國第九次森林資源清查顯示，湖南、江西、福建、浙江、安徽、四川、浙江、廣東等地為我國竹資源主產(chǎn)區(qū)(李玉敏和馮鵬飛, 2019)，與本研究所得適生區(qū)分布結(jié)果基本一致，結(jié)合已有的發(fā)生區(qū)報道，模型輸出結(jié)果較為準確，可為黃脊竹蝗在我國的風險分析提供參考。研究發(fā)現(xiàn)廣西西部、貴州東部為黃脊竹蝗潛在適生區(qū)，當?shù)刂脖缓蜌夂驐l件較為適宜，建議當?shù)亓謽I(yè)主管部門加強竹林監(jiān)測，以防傳入危害，造成損失。

環(huán)境主成分分析表明，影響黃脊竹蝗分布的最主要因子為最干月降水量和最冷月最低溫。結(jié)合黃脊竹蝗的生物學特性發(fā)現(xiàn)，6-8月是跳蝻羽化至成蟲的關(guān)鍵時期，降水量不足會加劇黃脊竹蝗失水風險(張威等, 2017)。黃脊竹蝗以卵越冬，室內(nèi)試驗發(fā)現(xiàn)8℃處理90 d以上卵解除滯育，冬季氣溫過高或過低會對胚胎發(fā)育造成損害，降低其孵化率(趙琴, 2009)。黃脊竹蝗產(chǎn)卵喜選擇背北向陽的山腰或竹坡產(chǎn)卵(孫立峰等, 2016)，但坡度、坡向因子在本模型中貢獻并不突出，推測可能與其成蟲在接近交配時期常進行長距離的群飛遷移,另尋產(chǎn)卵地有關(guān)。

就此次中老邊境黃脊竹蝗遷飛過境事件而言，筆者認為黃脊竹蝗在侵入普洱市、紅河州、西雙版納州等地后傳播擴散風險較小，主要原因如下：1)云南省竹資源主要分布于瀾滄江流域滇南地區(qū)，后續(xù)向北擴散寄主數(shù)量和范圍明顯減少；2)結(jié)合模型預測結(jié)果，個舊、玉溪、楚雄一代氣溫及降水量不適宜黃脊竹蝗生存，環(huán)境條件不利于蟲情擴散；3)黃脊竹蝗防治技術(shù)較為成熟，通過積極除治可有效降低蟲口密度。

MaxEnt模型顯示結(jié)果為物種分布可能性，不代表物種實際分布。本研究所得黃脊竹蝗適生區(qū)分布，主要基于19個環(huán)境變量和3個地形地貌因子，但我國竹產(chǎn)區(qū)主要受東南季風影響，云南等地局部也受西南季風影響，加之黃脊竹蝗遷飛能力較強，因此氣流可能是影響其遷飛傳播的關(guān)鍵因素。此外，本研究采取了默認參數(shù)，未選擇最低復雜度的模型參數(shù)來建模(朱耿平和喬慧捷, 2016)，可能會造成模型的轉(zhuǎn)移能力偏低，是否會對模型結(jié)果的準確性造成影響有待進一步研究。