松毛虫属(Dendrolimus)世界上已知30余种，主要分布于欧亚大陆，是针叶树的主要食叶害虫，该属害虫经常暴发成灾。如西伯利亚松毛虫(D. sibiricus)从19世纪到20世纪末毁坏俄罗斯林分400万hm²、西伯利亚西部契他州林分700万hm²(OEPP/EPPO，2005)；落叶松毛虫(D. superans)在20世纪中期危害库页岛林地达到15 356 hm²，1986—1989年俄罗斯和芬兰等地由落叶松毛虫造成的死亡林分达1 800 hm²(Pechenina，1991)；类似的情况在日本、韩国也有发生(Maeto，1991; OEPP/EPPO，2005)。而松毛虫作为我国主要林业害虫类群，具有分布广、种类多等特点，有29种松毛虫，其中6种经常暴发成灾(陈昌洁，1990)，每年平均有超过200万hm²的松林遭到各种松毛虫的危害，损失森林生物量超过360万m³。松毛虫在造成经济损失的同时也会导致原有森林生态系统的稳定性发生改变甚至对其造成破坏，使林分受到次期性害虫尤其是干部害虫攻击的风险性提高(Pechenina，1991; OEPP/EPPO，2005; 曾菊平等，2010)。研究环境因素和松毛虫生物学之间的关系，对种群数量和灾害发生作出有效的监测预报，是制定该类害虫科学有效管理措施的依据。

经典数学统计模型经常用于揭示灾害发生与各种环境生物因子之间关系。如应用回归模型分析英国冬季温度和蚜虫迁飞期的关系(Zhou et al.，1995)，建立西班牙多种蝴蝶出现高峰期和温度的关系模型(Stenfanescu et al.，2003)，应用逐步回归模型分析影响马尾松毛虫(D. punctatus)发生的主要气象因子(薛贤清等，1984；夏乃斌等，1990)；灰色系统和人工神经网络法也得到应用(陈绘画等，2003；2004)。而将数学模型和气候情景模拟相结合，可用来预测未来气候情景下的害虫发生情况。Yamamura等(2002)用1945—2001年的诱捕数据，建立了日本水稻二化螟(Chilo suppressalis)、叶蝉(Nephotettix cincticeps)和灰飞虱(Laodelphax striatellus)与温度的回归模型，并结合气候情景数据预测了2031—2050年每年害虫的数量变化趋势；Harrington等(2007)预测2050年蚜虫的发生期会随着温度的升高提前8天；结合有效积温模型也可对未来气候情景下的种群暴发和分布区域作出趋势预测(Morimoto et al.，1998; Yamamura et al.，2002; Bryant et al.，2002; 陈瑜等，2010)。

结合多种数学关系和物候关系的物种潜在适生区预测模型MaxEnt在物种潜在适生区预测上得到广泛应用(Pearson et al.，2004；赵文娟等，2009；Peterson et al.，2008)。MaxEnt是基于熵最大即误差最小的原理开发的一种计算机通用学习方法，起源于统计物理力学，通过最大熵法和贝叶斯原理运用不完整的有限信息对未知态作出最优预测或估计(Jaynes，1957)，在计算机科学和统计学领域，尤其是自然语言处理和识别方面得到广泛应用(Phillips et al.，2004；Berger et al.，1996)。Phillips等(2006)基于生态位理论并考虑多种环境因子，利用最大熵原理构建了物种地理潜在分布的生态位模型，在动物(齐增湘等，2011)、植物(Pearson et al.，2004；郭水良等，2011)、病菌(赵文娟等，2009)、入侵物种(雷军成等，2010)和昆虫(Peterson et al.，2008; Kroschel et al.，2013；李勤文等，2013)等方面得到应用，而在林业害虫发生预测方面仅见对美国西部的山松大小蠹(Dendroctonus ponderosae)、西部松大小蠹(Dendroctonus brevicomis)和松齿小蠹(Ips pini)在未来气候情景下灾害暴发的模拟预测(Paul et al.，2011)。MaxEnt模型在物种预测上的优势：其一，可以同时运用连续和离散数据并体现不同变量的组合关系，连续的输出结果能满足不同区域建模的需要，可满足有条件和无条件建模，当训练数据有限时其内在的通用性会有明显优势；其二，能利用独立变量和目前数据来界定环境范围并预测分布；其三，输出结果为连续结果，能很好地适应不同区域建模，模型有效性的检验函数可直观地表述模拟的准确度，即接受者操作特征曲线(receiver operating characteristic curve，ROC)曲线和Jackknife检验(Phillips et al.，2004; Phillips et al.，2006)。

本研究以我国北方油松(Pinus tabulaeformis)林内危害严重的油松毛虫(D. tabulaeformis)为例，基于历史灾情和气象因子运用MaxEnt模型对油松毛虫的灾害发生进行模拟预测。在分析判断MaxEnt训练模型对松毛虫历史灾害模拟拟合效果的基础上，结合不同气候变化背景进行未来变化分析，对油松毛虫未来灾害发生情况作出趋势预测。本研究不仅可探索预测预报松毛虫大尺度灾害发生的有效手段，也可扩展物种潜在适生区预测模型MaxEnt在森林保护学上的适用范围。

1 材料与方法 1.1 研究区概况

以山西省油松分布区为研究区域。山西省位于太行山以西的黄土高原上，地理坐标为34°34′—40°44′N，110°14′—114°33′E，总面积15.67万km²，山区面积占80%以上。海拔180~3 061.1 m，有山地、平原、丘陵等多种复杂地貌。属温带大陆性季风气候，年降水量400 mm，多季节性河流，河水流量随时间的变化较明显。复杂的地理环境孕育了丰富的植物资源和植被类型，南部以落叶阔叶林和针阔混交林为主；北部植被较少，主要以草原为主；中部以针叶林为主阔叶林为辅。针叶林多以松纯林为主，易受病虫害侵染危害，油松毛虫、天牛和红脂大小蠹(Dendroctonus valens)为主要害虫(宋玉双等，2000)。油松毛虫在山西油松林区均有分布，危害严重区主要在太原、大同、阳泉、长治、晋城、晋中、忻州、临汾和五台、太岳、中条和吕梁山林区等地。1978—2014年松毛虫仅在山西太原就出现了3次大发生高峰，1978年发生面积4 000 hm²，1989年85 733 hm²，2001年8 000 hm²(郎东升，2004； 国家林业局森林病虫害防治总站等，2013)，在山西的暴发期距为10年左右，对林业安全和生产造成严重损失。

1.2 松毛虫灾害数据

数据均来源于国家林业局森林病虫害防治总站历年的监测调查统计，对其中2002—2011年松毛虫全国发生不同灾害程度(轻、中、重)的县市级数据进行筛选，取得10年间山西有重度发生记录的数据，并以县级为最小数据单位定义灾害发生点，结合Googleearth 7.1.2.2042确定其几何中心的经纬度坐标，10年中油松毛虫重度发生点合计为69个。将69个点分为2组，即训练数据层(A)和预测检验数据层(B，C，D)(表 1)。

将A，B，C和D 4个数据层以*.csv格式文档保存。MaxEnt运用的最小凸多边形(MCPs)方法通过从发生点分析得到的背景候选点来提高模型的准确性(Vanderwal et al.，2009)。

1.3 物候数据 1.3.1 物候因子衍生

油松毛虫的生长、发育、繁殖、迁徙、种群数量的变化除受到种群自身暴发规律影响外，还受环境因子的影响，其中的气候因子对生物种群和暴发周期的影响更为显著(陈昌洁，1990；韩瑞东等，2005；何善勇等，2012)。温度、湿度、降水、风速和光照等气象因子通过对松毛虫生活史各个阶段的影响来调节种群数量(陈昌洁，1990)。结合油松毛虫在山西的生活史和已有的研究成果，共衍生出和油松毛虫幼虫期、蛹期、成虫期、卵期、越冬和上下树期相关的物候因子80个，包括温度、湿度、降雨、风速和光照等(表 2)。

1.3.2 物候因子数据及其筛选

气象数据均来自国家气象数据共享网(http://cdc.cma.gov.cn/ho-me.do)，包括2002—2011年山西19个气象站点的年、月和日值地面气象数据资料。由于油松毛虫在山西全境内均有分布，即在山西境内有油松分布的地方就有松毛虫分布，故气象站点的分布可代表研究区的气候特征及变化趋势(图 1)(王鸿斌等，2007)。从气象数据中统计提取出10年间80个物候因子数据值。结合油松毛虫重度发生情况，运用主成分分析和逐步回归法筛选出与油松毛虫重度发生相关的前8个物候因子，即 X₂₉，X₄₃，X₅₄，X₅₅，X₅₆，X₆₂，X₆₃，X₆₇。其中X₂₉，X₄₃，X₅₄，X₅₆与油松毛虫重度灾害呈正相关，而X₅₅，X₆₂，X₆₃，X₆₇则呈负相关。

1.3.3 MaxEnt模拟预测数据层

将上述8个物候因子数据按照1.2中的分组方法分为模型训练数据层和预测检验数据层，每个数据层也由8个物候因子组成。训练数据层定义为a，检验预测数据层定义为b，c，d(表 3)。对每个数据层的所有因子用下式去量纲化：

$x=\frac{{{x}_{i}}-{{x}_{\min }}}{{{x}_{\operatorname{m}\text{ax}}}-{{x}_{\min }}}\circ $

式中：x_i为因子原始值，x_min为因子最小值，x_max为因子最大值。

将去量纲化后各数据层运用ArcGIS 9.3空间分析模块(spacial analysis)中的样条插值函数进行空间插值，最后得出空间分辨率为5.0 km×5.0 km的连续表面栅格数据层，并将其转化为*.asc格式保存。

1.4 MaxEnt模拟训练

本研究运用的MaxEnt software v.3.3.3k(http:www.cs.princeton.edu/～schapire/MaxEnt)为目前发布的最新版本。MaxEnt基于熵最大原理和生态学理论，可在生态学研究中运用(Phillips et al.，2004； 2006)。就本文而言，将研究区定义为由地理点元素构成的集合X，对集合中每个点灾害发生的概率为π。π构成和为1的非负分布π(x)，则π的最优估计为灾害发生的概率，记为${\hat{\pi }}$，则其原理函数为：

$H(\hat{\pi })=-\sum\limits_{x\in X}^{{}}{\hat{\pi }(x)\ln \hat{\pi }}(x)\circ $

式中，$H(\hat{\pi })$为熵值，满足最大熵原则的概率分布为${{{\hat{\pi }}}^{*}}\text{=}\arg \max H(\hat{\pi })$。

1.4.1 模型训练

由A代表的灾害数据层和a代表的物候训练层模拟其灾害发生情况。模型输出结果为ASCII码栅格图层(*.asc)。数据中75%用于训练，25%用于训练中检验。重复10次迭代运算并输出平均模拟结果。选择ROC曲线分析模型的精度，选择刀切法(Jackknife)分析物候因子对灾害影响的大小，其他参数保持默认设置。

1.4.2 模型预测检验

以上述训练模型为基础，分别以2003，2007，2011年物候数据层来模拟对应各年份的灾害发生情况，而后将各年份灾害数据点层和对应的模拟结果层叠加并结合油松毛虫寄主油松的分布来判断模拟的准确性。

1.4.3 未来灾害趋势模拟

全球气候变化尤其是温度和降雨的变化将改变物种栖息地种群数量的变化规律。温、湿度的改变和油松毛虫的生活史关系密切，气候变化也会对其灾害情况造成影响，使其灾害分布在现有基础上发生改变。以上述模型为基础，在世界气象数据共享网(http://www.worldclim.org)上获得IPCC5最新发布的2种外排情景(RCP4.5和RCP6.0)、3种气候模式(BCC-CSM1-1、HadGEM2-ES和MIROC-ESM-CHEM)下的与昆虫生活史有关的衍生物候因子数据作为其未来环境变量，以1950—2000的数据作为其当前环境变量，用2002—2011年我国油松毛虫重度发生点作为当前灾害数据，对我国未来2种外排情景下的油松毛虫灾害的发生趋势进行预测。

2 结果与分析 2.1 模型训练

将油松毛虫在山西的灾害模拟结果利用ArcGIS中Natural Breaks分级方法重新进行分类，将油松毛虫灾害在山西发生的风险指数分为5个等级，分别为无风险区(0.0~0.17)、低风险区(0.17~0.29)、中风险区(0.29~0.40)、高风险区(0.40~0.53)和超高风险区(0.53~0.82)(图 2a)。超高风险区主要集中在广灵、灵丘、五台山和晋东南等地。对油松毛虫灾害严重发生的地点以县市为单位做频度分析(图 2b)，其结果与训练结果相同。

ROC检验可反映模型的模拟预测精度，AUC(area under the ROC curve)值越大，表示随机分布相距越远，即物候变量与灾害发生模拟相关性越大，表现为模拟预测效果越好。ROC曲线的评估标准根据AUC值可分为5级，失败(0.50~0.60)、较差(0.60~0.70)、一般(0.70~0.80)、较好(0.80~0.90)、非常好(0.90~1.0)(王运生等，2007)。本文训练模型的ROC值如图 3a所示，AUC值为0.820，SD为0.019，表明模拟效果较好。

刀切法检验分析结果能反映物候因子对油松毛虫灾害发生影响的大小。如图 4所示模型模拟中，对模拟结果贡献最大的是 X₄₃，X₅₄和X₅₅，而贡献最小的是X₂₉和X₆₂。X₄₃为7月日均相对湿度大于75%的天数，7月正是山西松毛虫羽化产卵期，老熟幼虫化蛹羽化时，35%以下的湿度会使其逐渐死亡而不能化蛹，随着环境湿度的升高，化蛹成功率也会逐渐升高，当相对湿度达到85%以上时，均可正常化蛹羽化；成虫交配产卵后，湿度75%以上的场所中的卵，孵化率在90%以上，而70%以下的湿度则不利于孵化，且湿度愈小孵化率愈低(陈昌洁，1990)。X₅₄为3月平均风速、X₅₅为4，5，6月平均风速，前者会通过影响越冬幼虫的上树期来影响幼虫数量，后者适当的风速可使幼龄幼虫迅速传播而对灾害的扩大起到促进作用(陈昌洁，1990)。

2.2 检验预测

在训练模型的基础上，分别对3个检验预测层进行预测模拟后，用与训练结果相同的方法分级处理，并将预测结果和代表当年的灾害发生县市的数据点叠加，由叠加结果可判断预测的准确性(图 5)。预测结果均显示油松毛虫灾害发生都以灵丘、广灵、五台山和晋东南为代表的区域最为严重。分别对不同风险区的灾害发生点统计，并对2003，2007和2011年灾害发生风险程度变化做卡方检验(χ²=6.000，df=4，P=0.199>0.05)，结果表明预测结果具有相似性。

模型预测的ROC曲线检验值AUC在2003年为0.817，2007年为0.797，2011年为0.804(图 3b，c，d)，即表明模型预测2003年和2011年的预测结果较好，而2007年模型的预测结果一般。刀切法检验结果与训练结果相一致，表明筛选的物候因子可用于山西油松毛虫灾害发生预测。

2.3 未来气候情景趋势预测

图 6a为根据历史气候条件得到的当前气候情形下灾害模拟预测情况，叠加油松分布情况后，可见油松毛虫主要发生在北京、河北、山西，辽宁、内蒙古与河北交界处，而在2002—2011年我国油松毛虫以省为单位的灾害严重发生频度统计中，以上区域每年都有严重灾害发生，说明模型有相对的准确性。图 6b，c分别是未来2050年时2种外排模式(RCP4.5，RCP6.0)下3种模拟气候情景(BCC-CSM1-1、HadGEM2-ES和MIROC-ESM-CHEM)的平均预测结果。图 6b说明2050年在RCP4.5下，油松毛虫灾害发生会整体向北移动，灾害发生将以有油松分布的北京、河北及河北与内蒙古交界处为中心，而山西南部及陕西等地灾害发生概率会降低；图 6c则显示2050年在RCP6.0下灾害发生趋势为在当前基础上河北北部和山西中南部有所加重。气候变化尤其是温度升高使昆虫和寄主的分布界限北移，不规律的降水和极端气候造成灾害发生的适宜环境可能随时变化(Harringtion et al.，2007)。

未来预测的ROC检验中，AUC为0.984，SD为0.009，模型拟合性非常好。MaxEnt模拟预测以训练模型为基础，训练模型精度的稳定性与训练数据样本的大小直接相关，当样本较大时，模型精度会升高(陈新美等，2012)。

3 结论与讨论 3.1 结论

物种地理潜在适生区预测模型MaxEnt可利用 气象物候因子对油松毛虫灾害发生进行模拟预测。本模拟预测中，7月平均湿度>75%的天数(X₄₃)为最主要因素，其次是3月和4，5，6月平均风速(X₅₄，X₅₅)，而10月的低温天数和风速在模型预测中作用最小(X₂₉，X₆₂)。

在未来气候情形下，油松毛虫重度灾害发生情况与油松分布区叠加结果显示，不同的外排模式下油松毛虫的灾害发生趋势不同，具体表现为：在油松分布区RCP4.5下灾害发生中心会向北京、河北、河北与内蒙古交界处集中，山西等地灾害发生概率则会降低；在RCP6.0下向山西中南部扩散。即环境条件的变化会引起灾害发生的相应变化(Musolin，2007; Kroschel et al.，2013)。

3.2 讨论

油松毛虫种群数量呈现出周期性暴发的特点，表现为周期不等的暴发和崩溃规律，一般是7~10年暴发1次。1个周期可分为4个阶段，即潜伏期、上升期、猖獗期和衰退期。影响油松毛虫种群变化的因素有环境因素和自身因素，各种环境因素在影响油松毛虫的各个阶段和暴发规律中发挥重要作用(Levins，1968)。因此，基于生态位理论构建的物种地理潜在适生区模型适用于松毛虫灾害暴发的模拟预测，但本研究中，由于模型模拟训练的样本量较少(仅为山西省)，故模型精度稳定性不够，增加模型训练的灾害暴发点样本量，或可提高模型的准确性(陈新美等，2012)。

灾害预测时本研究用到一个假设，即在未来条件下，对油松毛虫灾害发生起作用的是目前使用的关键气候因子而不是其他气候因子，在这个假设成立的前提下对未来的灾害趋势预测才是有意义的。同时模型在时间和空间尺度上存在无法检验的局限(Heikkinen et al.，2006)。

本研究只分析了气象因子在油松毛虫灾害发生模拟预测中的作用，而灾害的发生同时也受到其他环境因子的影响，包括种间关系、适应性、地理环境、经营措施、林火、天敌、社会经济发展水平以及其他影响因素(Heikkinen et al.，2006; Sinclair et al.，2010)。今后可在更小尺度上对油松毛虫的灾害发生情况做模拟预测，将寄主林分的树种、树龄组成及林分的抗病性等都纳入模型中，得到更加周密的预测结果，这对经营森林资源、监测防治害虫发生更具参考价值。