基于MaxEnt模型预测黄脊竹蝗在中国的适生区

温玄烨，王 越，姜 璠，唐 健，林 晓

(国家林业和草原局森林和草原病虫害防治总站，林业有害生物监测预警国家林业和草原局重点实验室，沈阳 110034)

黄脊竹蝗CeracriskiangsuTsai，又称竹蝗，属直翅目Orthoptera网翅蝗科Arcypteridae，以跳蝻和成虫取食叶及嫩梢，主要危害刚竹属Phyllostachys、箣竹属Bambusa竹种，也可危害玉米Zeamays、水稻Oryzasativa、棕榈Trachycarpusfortunei等农作物及近百种杂草(钟洪武等, 2010)。黄脊竹蝗食性杂、取食量大、危害时间长、防治困难，对竹子产业危害严重，在中国南部、越南、老挝等热带和亚热带地区均有分布(伍远平等, 2005; Yuetal., 2011; 孔志坚和寸佳莅, 2018)。

黄脊竹蝗是我国记录最早的林业害虫之一，早在明朝嘉庆年间就曾大规模发生(程佳等, 2010)。上世纪50年代，黄脊竹蝗曾在我国南方地区暴发，造成了严重的经济损失。近年来，我国经济林面积不断扩大，黄脊竹蝗危害又有抬头之势。2020年6月以来，云南省中老边境发生境外黄脊竹蝗迁飞入境事件，对相关地区的农林业生产造成危害，引起了农业农村部和国家林业和草原局的高度重视。我国学者对黄脊竹蝗开展了大量研究，在生物学和防治技术等领域取得了丰富的成果和进展(程佳等, 2009; 张威等, 2016; 高尚等, 2020)，但关于其在我国的潜在分布区研究仍属空白。因此，开展黄脊竹蝗在我国的潜在分布区预测及适生性分析研究，对做好灾害预防和风险评估有重要意义。

生态位模型是通过物种分布及环境变量，利用数学模型描述物种生态需求并进行空间投影，分析物种适生性的方法，目前已应用于珍稀植物分布、外来物种入侵、野生动物保护等领域，均取得较好效果(迟翔文等, 2009; 李望军等, 2019; 冯春慧等, 2020)。最大熵(MaxEnt)模型通过已知空间分布，寻找关键环境变量，构建约束集合，模拟二者关系的生态位模型，其优点是预测结果准确、需求样本数量少，可较好地预测物种的适生区，国内现已基于MaxEnt模型对红脂大小蠹Dendroctonusvalens(王涛等, 2018)、美国白蛾Hyphantriacuneain(李涛等, 2018)、松针红斑病Dothistromapini(王晓伟等, 2019)等林业有害生物适生区进行预测，均得到了较好的效果，准确地反映了其潜在适生区，为林业有害生物的防控提供了依据。

本研究以黄脊竹蝗为研究对象，基于本次云南黄脊竹蝗迁飞入境点位和我国已报道的物种分布信息，在当前气候条件下，通过MaxEnt模型结合ArcGIS对黄脊竹蝗在我国的适生区进行预测，分析影响其生存的关键环境变量，以期为黄脊竹蝗的监测预警和综合防治提供理论依据。

1 材料与方法

1.1 材料

分布数据：黄脊竹蝗的分布数据来源于实地调查、检索中国林业有害生物防治信息管理系统(http://os.forestpest.cn/)和相关文献资料(黄焕华等, 2003; 陈良昌等, 2013; 方蓉等, 2015; 杨华等, 2018; 宋玉双, 2019)。

环境变量：从世界气象WorldClim Version 2.0数据库(http://www.worldclim.org)中下载1970-2000年19个栅格大小为25 km2生物气候环境变量数据Bio1～Bio19，分别为年均气温、昼夜温差月均值、等温性、温度变化方差、最热月份最高温、最冷月份最低温、年温变化范围、最湿季度平均温度、最干季度平均温度、最暖季度平均温度、最冷季度平均温度、年降水量、最湿月份降水量、最干月份降水量、降水量变化方差、最湿季度降水量、最干季度降水量、最暖季度降水量、最冷季度降水量。

地图地形：来源于国家自然资源部标准地图服务系统(http://bzdt.ch.mnr.gov.cn/)下载的1 ∶1 600万标准中国地图，审图号是GS(2016)2923号，通过ArcGIS空间分析工具获取3个地形因子(海拔、坡向和坡度)。

1.2 方法

1.2.1数据处理

通过数据库或文献检索到的物种分布记录，部分数据缺失经纬度坐标，需要通过Google Earth v7.1.8输入地名查询。为避免地理空间采样过于密集导致生态空间分布点重复(朱耿平等, 2014)，将所得数据利用分辨率为2.5 arc-min的栅格去除重复，降低采样点的空间偏差。最终获得物种分布点合计267个，其中，侵入点25个，已发生数据242个(图1)，根据MaxEnt软件要求，将所得数据整理成“.CSV”格式备用。

1.2.2主要环境变量的筛选

由于黄脊竹蝗种群的分布与存活与环境变量存在相关性，为减少环境因子的多重贡献性对模型的影响，避免空间相关性对建模产生误差，利用SPSS 14.0软件对19个环境因子进行Pearson共线性检验。当2个环境变量相对系数≥0.8时，表明2个变量间存在共性关系，选取贡献较高的环境变量，当2个环境变量相对系数<0.8时，表明变量间不存在共性关系，保留该变量用于模型分析。

筛选得到等温性(bio3)、最热月最高温(bio5)、最冷月最低温(bio6)、温度年较差(bio7)、最干月降水量(bio14)、降水量季节性变异系数(bio15)、最湿季降水量(bio16)等7个环境变量(表1)和3个地形因子海拔(Altitude)、坡向(Aspect)和坡度(Slope)用于进一步分析。

图1 黄脊竹蝗分布点位Fig.1 Distribution point of Ceracris kiangsu注：该图基于自然资源部标准地图服务系统下载的审图号为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 环境变量的相关性矩阵

1.2.3参数设置

将黄脊竹蝗地理分布点数据和主要环境因子数据导入MaxEnt软件，选择随机25%数据坐标点作为测试数据，75%坐标点为训练数据，在默认收敛域限10-5,最大迭代次数500状态下，进行10次重复，选取平均值作为黄脊竹蝗分布预测结果，输出结果以ASC文件形式显示(Peterson and Eaton, 2010)。

1.2.4黄脊竹蝗适生区预测与评价

将MaxEnt模型得到的黄脊竹蝗潜在分布结果输入ArcGIS，转化为栅格形式(赵健等, 2019)。按照模型生成的适生概率(P)，利用自然间断点分级法(Jenks)，结合黄脊竹蝗在发生区所造成的危害严重程度，将黄脊竹蝗的适生区分为4级：非适生区(0≤P<0.05)，低适生区(0.05≤P<0.25)，中适生区(0.25≤P<0.5)，高适生区(0.5≤P<1)，利用ArcGIS中重分类功能获取黄脊竹蝗在我国的适生区等级分布图(谢岷等, 2019)。

1.2.5环境变量对黄脊竹蝗生存率的影响

采用Jackknife(刀切法)探讨环境变量和地形因子对黄脊竹蝗分布区域和生存概率的影响，通过1.2.2获得10个环境变量(包括地形因子的3个变量)，依次去除1个环境变量，利用剩余变量建模，分析去除该变量后的相关性，从而确定10个环境变量与黄脊竹蝗生存的相关程度，即生存率。

1.2.6模型评价

通过ROC(受试者工作)曲线对模型精度进行评价。以特异度(假阳性率)为X轴，灵敏度(遗漏率)为Y轴，计算ROC曲线下积值即AUC值对模型的精度进行评价，当AUC值<0.7时，预测结果准确性较差，当0.7≤AUC值<0.9时，预测结果准确性一般，当AUC>0.9时，预测结果准确性高(Merowetal., 2013)。

2 结果与分析

2.1 模型精度评价

对MaxEnt模型进行准确性检验，得出训练集和测试集10次重复结果(图2)。AUC平均值为0.958，平均标准误差为0.005，表明此模型预测准确度高，利用现有分布点数据和环境变量可有效预测黄脊竹蝗在我国的潜在分布范围。

图2 黄脊竹蝗适生区的ROC验证曲线Fig.2 ROC curve of potential distribution of Ceracris kiangsu

2.2 黄脊竹蝗在我国的适生区

MaxEnt模型预测结果显示，黄脊竹蝗在我国江淮流域、长江中下游地区、华南及西南等地有不同程度的适生区域，其适生区面积占我国陆地总面积的18.9%，其中，高适生区占3.0%，中适生区占5.6%，低适生区占10.3%。高适生区主要分布在长江中下游的江西、湖南、浙江、安徽，华南沿海的广东、广西，以及西南的云南等省区。以湘赣地区的鄱阳湖平原、两湖平原为中心向四周发散式分布，南至两广丘陵北部，东至浙西丘陵地带，北至鄂皖地区的大别山南麓，向西跳跃式分布于云南中老边境及重庆中部地区。中适生区主要分布于安徽中部、湖北西北部、湖南西部、广西大部、重庆西北部及云南西南部地区，以及上述高适生区的周边区域。低适生区广泛分布于我国淮河以南各主要省区，向南可至海南省北部及台湾省南部地区(图3)。

2.3 环境变量对黄脊竹蝗分布的影响

应用刀切法对环境变量进行分析，结果显示，最干月降水量(bio14)和最冷月最低温(bio6)贡献率分别51.5%和22.3%，正规化训练增益值>1.4，对模型提供了最大的效益，说明使用最干月降水量和最冷月最低温时比其他变量包含更多有用信息，对黄脊竹蝗分布的影响最大。最热月最高温(bio5)、最湿季降水量(bio16)、降水量季节性变异系数(bio15)、温度年较差(bio7)正规化训练增益值>1.0，四个环节变量累计贡献率为22.1%，说明对黄脊竹蝗的分布产生较大影响。而等温性(bio3)、海拔(altitude)、坡向(aspect)和坡度(slope)等环境变量对黄脊竹蝗的分布影响相对较小(图4)。

图3 黄脊竹蝗在我国的适生区预测分布图Fig.3 Predicted potential suitable habitats of Ceracris kiangsu注：该图基于自然资源部标准地图服务系统下载的审图号为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 刀切法分析环境变量的重要程度Fig.4 Jackknife test for evaluating the relative importance of environmental variables

影响黄脊竹蝗分布的主要环境变量是降水量和温度，从降水量来看，当最干月降水量不足50 mm时，黄脊竹蝗的生存率不足0.4，当降水量超过200 mm时，生存率上升至0.66(图5-A)，按其适生区等级划分方法，以P=0.05为阈值，得到降水量变化范围为7～2 619 mm(图5-B)，最低适宜值为72.16 mm。从温度来看，最冷月最低温极值为1℃(图5-C)，低于或高于该值时，黄脊竹蝗生存率降低，最热月最高温度不足25℃时，黄脊竹蝗生存率不足0.1，随着温度上升其生存率也逐渐上升，33.0℃以上时，黄脊竹蝗的生存率大于0.5。

图5 不同环境变量与黄脊竹蝗生存概率间的关系Fig.5 Relationship between different environmental variables and survival probability of Ceracris kiangsu注：图层中的温度数值(℃)为实际数值×10。Note: Temperature value (℃) in the layer was the actual value×10.

3 结论与讨论

黄脊竹蝗属寄主主导型食叶害虫，其生长发育与寄主种类密切相关(张守科等, 2017)，在景观尺度上，昆虫的适生区范围与其寄主分布直接相关(王梦琳等, 2017)。我国第九次森林资源清查显示，湖南、江西、福建、浙江、安徽、四川、浙江、广东等地为我国竹资源主产区(李玉敏和冯鹏飞, 2019)，与本研究所得适生区分布结果基本一致，结合已有的发生区报道，模型输出结果较为准确，可为黄脊竹蝗在我国的风险分析提供参考。研究发现广西西部、贵州东部为黄脊竹蝗潜在适生区，当地植被和气候条件较为适宜，建议当地林业主管部门加强竹林监测，以防传入危害，造成损失。

环境主成分分析表明，影响黄脊竹蝗分布的最主要因子为最干月降水量和最冷月最低温。结合黄脊竹蝗的生物学特性发现，6-8月是跳蝻羽化至成虫的关键时期，降水量不足会加剧黄脊竹蝗失水风险(张威等, 2017)。黄脊竹蝗以卵越冬，室内试验发现8℃处理90 d以上卵解除滞育，冬季气温过高或过低会对胚胎发育造成损害，降低其孵化率(赵琴, 2009)。黄脊竹蝗产卵喜选择背北向阳的山腰或竹坡产卵(孙立峰等, 2016)，但坡度、坡向因子在本模型中贡献并不突出，推测可能与其成虫在接近交配时期常进行长距离的群飞迁移,另寻产卵地有关。

就此次中老边境黄脊竹蝗迁飞过境事件而言，笔者认为黄脊竹蝗在侵入普洱市、红河州、西双版纳州等地后传播扩散风险较小，主要原因如下：1)云南省竹资源主要分布于澜沧江流域滇南地区，后续向北扩散寄主数量和范围明显减少；2)结合模型预测结果，个旧、玉溪、楚雄一代气温及降水量不适宜黄脊竹蝗生存，环境条件不利于虫情扩散；3)黄脊竹蝗防治技术较为成熟，通过积极除治可有效降低虫口密度。

MaxEnt模型显示结果为物种分布可能性，不代表物种实际分布。本研究所得黄脊竹蝗适生区分布，主要基于19个环境变量和3个地形地貌因子，但我国竹产区主要受东南季风影响，云南等地局部也受西南季风影响，加之黄脊竹蝗迁飞能力较强，因此气流可能是影响其迁飞传播的关键因素。此外，本研究采取了默认参数，未选择最低复杂度的模型参数来建模(朱耿平和乔慧捷, 2016)，可能会造成模型的转移能力偏低，是否会对模型结果的准确性造成影响有待进一步研究。