基于地理栅格变量与机器学习的松材线虫病扩散风险分析

刘婷婷，杨晋帆，周汝良，刘 琳

（1.西南林业大学 地理与生态旅游学院，云南 昆明 650224；2.西南林业大学 理学院，云南 昆明 650224）

由松材线虫Bursaphelenchusxylophilus引起的松材线虫病是目前松树上最具破坏性的病害之一。自1905年日本首次报道发生松材线虫病以来，亚欧其他国家相继出现松材线虫病疫区[1-2]。中国的松材线虫病最早报道于1982年，出现在南京中山陵的黑松Pinusthunbergii上，之后该病害快速扩散，据国家林业和草原局统计，中国现共有松材线虫病县级疫区726 个、乡镇级疫点5 479 个，发生面积达180.9 万hm2，包括轻型疫区206 个(含2021年初新增疫情发生区)、重型疫区518 个[3-4]。目前松材线虫病已成为中国最为严重的林业病害之一[5-6]。为实现《全国松材线虫病疫情防控五年攻坚行动计划(2021—2025)》中新发疫情“早发现、早报告、早除治、早拔除”的目的，达到控制增量的总体目标，就需要把握疫情的总体趋势，精准定位和定量各疫区的疫情，提高防疫针对性[7]。

已有研究发现：全球变暖环境下，气温增高会影响松材线虫病的地理分布[8]；人类活动是影响松材线虫病远距离、跳跃式、梅花式扩散传播的重要因素[9-11]；年均气温、年均降水量等是影响松材线虫定殖的主要气候因素[12]；中国松材线虫病的发生受到寄主植物分布的强烈影响[13-14]；媒介昆虫对寄主的选择性也影响松材线虫病的发生[15-17]。当前对松材线虫病传播风险的研究多集中于采用传统的地统计学方法来分析松材线虫病在省、市范围的风险格局[18-20]，而将机器学习引入病害风险测报研究较少见。机器学习方法中随机森林、支持向量机具有需要样本量少，能够同时利用多变量进行非线性回归、分类预测并且结果精度高等优点，使风险预测精确化分析成为可能[21-23]。

本研究综合考虑松材线虫病扩散的影响因素，利用地理环境、气候环境、卫星遥感、道路交通网络等构成地理栅格数据集。对比随机森林和支持向量机2 种算法建立松材线虫病风险扩散模型，探究机器学习算法在松材线虫病扩散传播预测研究中的可行性和准确性，寻求最佳预测方案，将潜在疫区根据等级进行风险区划分，并进行栅格化分析，将全国松材线虫病测报落实到“山头地块”上，以期为林业生态系统中松科Pinaceae 植物的动态监测提供有效的手段。

1 数据与方法

1.1 疫区样本数据

以中国县区为行政单元，以国家林业和草原局2018—2021年关于松材线虫病疫区公告为基础[4]，收集整理了中国松材线虫病发生区731 个以及未发生区2 131 个，以松材线虫病发生(1)或未发生(0)作为模型响应变量，由点及面构建了2 862 个样点。其中随机选取30%样本点(954 个)作为模型独立检验样本。总样本点分布见表1。

表1 松材线虫病疫区Table 1 Pine wilt disease endemic area

1.2 测报变量空间模拟

1.2.1 气候环境数据 根据以往研究[8,12]，选取1970—2000年共30 a 逐月气象数据获得年均气温、年均降水量以及年均最高气温、年均最低气温、年均最低降水量与年均最高降水量数据。空间分辨率为 250 m的气候数据集来源于https://www.worldclim.org/。

1.2.2 传播媒介数据 选取影响昆虫媒介墨天牛属Monochamus的地理因子，包括高程数据，以及根据高程数据计算提取的坡度与坡向指数。选取海拔、坡度、植被覆盖度以及温度构建阻力面，表征其自然传播的能力。根据全国高速公路、国道、省道以及一级水系等交通网络，利用欧氏距离来模拟人类活动导致的有害生物空间传播度量[9]。道路、水系矢量数据等人为影响数据来源于Bigemap 平台，利用ArcGIS 处理形成综合交通便利性数据；土地利用类型数据来源于GLOBELAND 30 平台。

1.2.3 寄主植物数据 全国1∶100 万植被分布数据来源于资源环境数据云平台http://www.Resdc.cn/，从中提取出易受松材线虫病感染的针叶林植被。根据前人研究结果[24]将不同树种的寄主易感指数赋分为0～100，数值越高表示寄主易感性越高，划分标准见表2。

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

综上，选择16 个环境变量因子作为主要环境数据(表3)。环境数据集的空间分辨率通过重采样统一为250 m。地理数据为1∶320 万的中国国界、省界以及县界行政区划图，来自于自然资源部标准地图服务网http://bzdt.ch.mnr.gov.cn/。

表3 松材线虫病预测模型自变量表Table 3 Independent variables of the prediction model of pine wilt disease

2 扩散风险模型构建

2.1 模型

随机森林(Random Forest Classifier，RF)是一种监督学习的机器学习算法，即将多个决策树分类器集合在一起进行数据的分类和回归[25-27]。模型从原始样本抽取的随机性，以及利用小于样本集特征的数据建立训练集，增加了不同决策树结果的差异性，中和了过拟合现象，使得模型的泛化能力提高使其具有较高的分类准确性[28-29]。支持向量机(Support Vector Machine，SVM)算法能够解决机器学习中样本量小、非线性和高维模式识别中“维数灾难”和“过学习”等方面的问题[30-31]。该方法选择适当的核函数通过非线性变换将输入空间变换到一个高维空间，然后在这个高维空间求取最优分类面，找到输入变量和输出变量之间的一种非线性关系。

2.2 精度评价

ROC-AUC 特征曲线分析方法用来评价二值分类器的优劣[32]。ROC 曲线反映了输出概率的准确性，曲线上每个点反映对同一信号刺激的感受性。横轴为负正类率(false postive rate，FPR)特异度，纵轴为真正类率(true postive rate，TPR)，灵敏度均为[0,1]。ROC 曲线在斜对角线以下，表示该分类器效果差于随机分类器，反之，效果好于随机分类器[33]。ROC 曲线下的面积(area under curve，AUC)作为数值可以直观地评价分类器的好坏。

2.3 扩散风险模拟

RF 和SVM 模型均以松材线虫病在全国县域尺度中存在(1)或不存在(0)作为响应变量，对所有环境变量在输入前进行归一化处理，以消除由于纲量导致的差异。数据集以7∶3 的比例随机分成两部分进行模型训练和测试。在建模过程中16 个环境变量根据模型重要性及其影响，筛除部分低贡献度变量，并将针叶林分布变量提出作为重要变量与预测风险区叠加分析。模型预测精度和ROC-AUC 曲线用来评价每个模型的性能，利用预测结果与独立检验样本通过混淆矩阵进行精确性检验。通过Python 转为0～1 之间的发生概率，输入ArcGIS 中进行风险可视化处理。除此之外，全国栅格尺度更能够体现空间连续变化，在进行运算操作、叠加分析方面也比县级矢量数据更具优势[34]。利用从全国1∶100 万植被分布类型中提取出来的寄主松科植物，根据松科植物不同的易感性重分类，栅格化后与模型预测的松材线虫病潜在发生区域叠加分析。

3 结果与分析

3.1 模型精度验证

通过模型精度和ROC-AUC 检验结果(表4)可知：2 个模型均有较好的预测能力，其预测精度和模型稳定性均在75%以上。模型精度是通过RF 和SVM 模型对独立检验样本数据进行松材线虫病的发生风险概率值预测，设定临界值为0.5，将大于此值县域视为松材线虫病“存在”(y=1)，小于此值归为“不存在”(y=0)，得出分类结果与独立检验样本的实际值来验证模型预测的精度。结果表明：构建模型的测试集预测精度在RF 和SVM 模型中总精度分别为83.95%和77.97%。

表4 模型精度验证数据表Table 4 Model accuracy verification data sheet

RF 和SVM 模型的ROC-AUC 曲线值(图1)分别为0.89 和0.77 (平均标准误差为0.01)。这证明2 个模型在预测松材线虫病县级风险分布过程中具有逻辑可行性，预测结果具有可信性，但RF 的效果、模型稳定性、模型分类结果泛化能力均高于SVM 模型。

图1 ROC-AUC 曲线Figure 1 ROC-AUC curve

3.2 变量重要性评价

根据RF 模型对数据贡献度(图2)结果得出，重要性排序中位于前列的为气候变量，包括年均最低降水量(0.302 59)、年均降水量(0.170 33)、年均低温(0.102 98)；其次是海拔变量(0.150 84)，这2 类是影响松材线虫病的重要环境变量。松材线虫病发生概率在低海拔地区较高，随着年均低温的升高风险概率增加，而高温、干旱、低湿气候直接影响松材线虫病的爆发。在传播媒介影响中，人为因素中各级交通网络影响力重要性更为突出(0.067 58)，说明密集的人流和物流是松材线虫病远距离传播的重要媒介。在重要性排序中松科植物类型数据贡献度与实际影响松材线虫病的情况不符，寄主是影响松材线虫病存在和发展的必然条件，故将针叶林分布单独提出栅格化，与预测风险区叠加分析其影响。

图2 RF 模型变量贡献度Figure 2 Contribution of RF model variables

3.3 松材线虫病县域风险测报

利用Python 将预测结果提取为0～1 之间的概率数据，依据离散程度划分为5 个等级，由高到低依次为极高风险、高风险、中等风险、低风险、极低风险。因样本原因，不涉及中国香港、台湾、澳门行政区。对比RF 与SVM 模型的分析结果(表5和表6)，SVM 对重要变量的模糊表达导致其预测结果等级差距小、重点疫区不突出、预测精确性低、全国大部分地区呈同一等级，这使得对松材线虫病疫区的监测与防治难度增大。相比而言，RF 预测结果具有更高的准确性，其潜在疫区分布与原有疫区重合度高，危险等级一致，并能够明显表达城市、道路以及地形对潜在疫区分布的影响。

表6 SVM 松材线虫病县域潜在风险等级Table 6 SVM potential risk levels of pine wilt disease in counties

3.4 全国栅格尺度松材线虫病风险测报

寄主植物根据其易感性划分等级，松材线虫病易感性松科植物多分布于浙江、福建、湖南、广西等省，另外，云南、贵州、四川以及东北林区也有大量寄主植物分布。利用RF 模型叠加寄主易感性分布数据对全国31 个省进行地理栅格尺度的疫区划分与分析，其中概率大于0.8 以上的栅格点为极易感染区域。结果(表7)表明：松材线虫病的极高风险区为华东地区的浙江、江西、福建，华南地区的广西、广东以及华中地区的湖南，高易感性寄主分布与极高风险区预测结果高度重合，疫区内多为赤松、黑松、马尾松；高风险区包括华东地区的安徽，华中地区的湖北，以及西南地区的重庆、贵州，高风险区寄主也呈高易感性，但其分布分散且多于其他植被混生，森林具有类型多样化；中等风险区包括华中的河南、江苏，东北地区辽宁，华东地区山东，寄主植物易感性降低，多呈中等易感性且面积小，多为散布状态，其中很多省份为粮食大省，以种植农作物为主，林业分布较少；无风险区包括西北地区甘肃、宁夏、青海、新疆，华北地区的内蒙古、山西和西南地区的西藏，缺少寄主植物的分布，人流、物流相较于其他地区也较少，难以形成松材线虫病潜在风险区；其余省份划分为低风险区，其中需注意的是云南、四川大部分地区存在高易感性的思茅松、云南松等优势种，风险等级容易上升。

表7 RF 与寄主易感性的松材线虫病地理栅格潜在风险预测Table 7 Prediction of the potential risk of geographic grids for RF and host susceptibility to pine wilt disease

4 讨论

相比于县域行政单元尺度的传统地学统计方法，利用地理栅格信息描述驱动变量，以RF 模型建模预报风险的方法，实现了测报变量的数字矩阵描述，解决了测报结果的空间分异性与连续化描述的难题。全国栅格尺度的松材线虫病风险测报能够更好地反映人类活动、地理环境变化、寄主易感性对松材线虫病传播轨迹的影响，对基层单位开展检疫、预防和治理等工作有指示和支持作用。

降水、气温和海拔是影响松材线虫病生存的主要因素，统计分析得出松材线虫病最低生存年均降水量范围为110～121 mm。影响松材线虫病发生的年均低温范围为10～13 ℃，这与叶建仁[35]以年均气温10～14 ℃区域为松材线虫病可以发生区，大于14 ℃为流行区划分的预测区域基本一致。影响松材线虫病发生的海拔范围为4～247 m，这与LEE 等[21]预测的松材线虫病最适生存海拔小于200 m 基本一致。在传播影响中，除了媒介昆虫通过自然传播，人类活动对松材线虫病的扩散也具有重要的影响。松材线虫病早期发生地点多位于城市化活跃、人类活动强烈的区域内，其原因为城市用地、道路用地增加以及人工育林替代原始植被降低了生态稳定性[36]，这为松材线虫病爆发与传播埋下巨大隐患。

全国地理栅格扩散风险格局分析表明：松材线虫病在华北平原和长江中下游平原持续扩散，并且存在3 条传播廊道：一是沿着长江水系向西南地区及四川盆地扩散，现以水热条件优越、高城市化的重庆潜在风险等级为最高；二是沿着黄渤海临海城市向东北平原的扩散廊道，这个廊道中北京和天津的潜在风险等级突出，可能是受到高城市化发展和密集交通网络的影响；三是沿金沙江、珠江流域形成了一条由广东、海南、广西向云贵高原扩散的通道，虽西南地区海拔较高，但其丰富的松科植物仍为松材线虫病的生存和传播提供基础。

5 结论

基于多类型数据建立的RF 与SVM 模型均有较好的预测性能(ROC-AUC＞75%)。RF 模型具有更高的精度为83.95%，并且预测的极高风险区与寄主高易感性分布高度重合，证明RF 模型分类结果更符合实际。

影响松材线虫病发生的敏感变量是年均降水量、年均低温和海拔，而坡度小、人类活动强度高是导致松材线虫病快速传播的重要因素。潜在扩散区位于人类活动密集的低海拔地区、道路通达的林区、城市城镇分布区和人工林分布区。其中，极高风险分布地区主要位于华东地区的浙江、江西及福建，华南地区的广西、广东以及华中地区的湖南。松材线虫扩散过程中存在明显传播廊道与地理阻隔。松材线虫在华北平原的扩散传播受到秦岭和太行山脉的阻隔，并未侵入西北地区的黄土高原，而高风险疫区多集中在黄河水系以南地区。

本研究仅使用了RF 及SVM 共2 种机器学习算法，样点数据仅为国家林业和草原局2018—2021年松材线虫病疫区数据，也还未考虑环境动态变化导致的松材线虫病疫区发生的风险变化。未来的研究应以松材线虫病实际病害发生位置为基础，通过更多的机器学习算法或深度学习算法进行对比研究，充分考虑与研究区相适应的自然环境、人类活动等时空数据，以期为松材线虫病疫区的预测提供更加准确、精细的风险测报。