区域尺度油松毛虫灾害发生相关气象因子筛选1)

宋雄刚 王鸿斌 李国宏

(国家林业局森林保护学重点实验室(中国林业科学研究院森林生态环境与保护研究所))，北京，100091)

于志军 陈国发 孔祥波 张真

(国家林业局森林病虫害防治总站) (国家林业局森林保护学重点实验室(中国林业科学研究院森林生态环境与保护研究所))

责任编辑:程 红。

农林病虫害是影响农林生产及其生态系统稳定的 重 大 自 然 灾 害 之 一［1－2］。松 毛 虫(Dendrolimus spp．)是我国松属树种主要的食叶害虫，大发生时会取食全部针叶，并引起人畜中毒;油松毛虫(D．tabulaeformis)是松毛虫中危害严重者之一，主要在北方发生［3］5，［4］。昆虫生活史的各个阶段对温度都较敏感［5］，种群数量的变化及暴发规律除受到自身遗传特性的影响外，还受到环境因素的重大影响［3］19，［6］166－192，［7－8］;气候因素是影响松毛虫种群发生规律的主要环境因素，不同的气候因子在松毛虫发育的不同时期可单独或交互的起到重要作用［3］22－23，［8］。在气象因素和松毛虫的发生具有显著相关性的基础上，根据多种统计学或统计力学原理可建立各种监测预报模型。张玉书等运用气象期距及有效积温法对阜新县油松毛虫的始见期预测，其预测结果比实际上树期晚一天［9］;周广学等运用相关分析和主成分分析法建立朝阳市油松毛虫越冬死亡率和发生量的预测模型，结果显示对死亡率影响最大的是春季温度、对发生量影响最大的是前一年10月份降水和当年1月中旬平均温度［10］;陈绘画等运用逐步回归法和人工神经网络模型建立了仙居县灾害预测模型，其准确性能达到100%［11］;Liang et al．用一般线性模型(GLM)研究了山松大小蠹(Dendroctonus ponderosae)在不同气候情形下死亡率分布［12］。上述均以气象因子的月或旬值为基础展开研究，而对气候因子动态变化及气候因子在昆虫不同发育期的累积效应考虑较少;费海泽等考虑了上述影响后用人工神经网络法建立广西马尾松毛虫(Dendrolimus punctatus)不同发生程度的预测模型且准确性较高［8］。气候因子变化对昆虫的影响还体现在灾害发生期的提前或分布区的扩大上［13－16］。上述研究大多以县级或具有相同立地条件的区域为研究对象，在小区域或相同的某一栖息环境中，可得到较准确的结果;而在包括多种立地、林型的更大尺度上或更宽区域，则其预测结果代表性下降。故本研究以辽宁、北京、河北和山西等主要发生区为代表，以油松毛虫2002—2011年的灾害监测数据(包括监测面积和不同发生程度)和相应的气候数据为基础，结合气候因子的动态变化和油松毛虫的不同生活阶段，运用典型相关分析法和主成分分析法，分析北方区域尺度上影响油松毛虫发生的主要物候因子。这些因子在气候变化背景下分析预测油松毛虫的灾情变化有重要意义。

1 研究区域概况

研究区由辽宁、北京、河北及山西4个省行政区组成。该区域为油松毛虫在我国北方的主要发生区。我国近12年来，油松毛虫的发生面积共158．69 hm2，其中2002年最为严重。从2009年开始油松毛虫发生面积围绕10 hm2波动(图1)(数据来源于森防总站的历年监测数据和总站的内部资料《林业有害资料监测预警年报》)，到2014年辽宁凌源大面积发生，松林镇成灾，至2015年凌源持续大发生，单株虫口密度达到2 000头(数据由辽宁凌源森防站提供)，所以油松毛虫的防治依然面临严峻形势。油松毛虫的灾害性暴发主要出现于我国北方的9个省市，且以辽宁省发生面积最多，达到33%，天津最少，接近0(表1)。研究区的气候因子、地形地势、寄主情况能反应油松毛虫在我国北方分布与发生的一般特征(图2)。

图1 2002—2013年全国油松毛虫发生面积变化趋势

图2 研究区地理地势、松属分布示意图

表1 2002—2013年各省市油松毛虫平均灾害面积及发生比例

2 数据资料来源

数据资料包括气象数据和虫情数据。其中气象数据由中国气象科学数据共享服务网(http://cdc．nmic．cn/home．do)和中国农业科学院农业环境与可持续发展研究所(http://www．isi－mip．org)提供，包括日、月、年的温度、相对湿度、降雨、日照和风速等;2002—2011年的油松毛虫县市级虫情数据(包括监测面积和各种程度的发生面积)由中国森防总站提供。

3 研究方法

3．1 衍生物候因子

气候因子影响油松毛虫生活史的各个期。油松毛虫在研究区1年1代，且跨年度完成发育，即油松毛虫的灾害情况和前一年、当年的环境因素相关，为了更科学的分析气候条件对油松毛虫世代生长发育和其暴发的影响，结合研究区油松毛虫的各个生长发育期对自然年界进行如下划分:首先将油松毛虫的生活史划分为4个生长发育期，即羽化—产卵—孵化—幼龄幼虫期、越冬期、上树期和老熟幼虫食叶—化蛹期;然后将这4个期松毛虫活动特征和气候结合并分别对应为其生长发育的Ⅰ、Ⅱ、Ⅲ和Ⅳ期(表2)。结合油松毛虫生活史不同时期对不同气象因子不同的敏感性和相关的研究文献［3］5－25、［5］、［7］，衍生物候因子。

表2 油松毛虫生长期对应时间

3．2 虫情灾害数据资料

在2002—2011年的县级数据库中，选择研究区中具有完整数据记录的县市级灾害发生点77个，并分别计算轻、中和重度发生程度的面积比例(p)。p为各个发生程度的发生面积和监测面积的比值。计算公式如下:

式中:Pij为i点j灾害程度比例;i为发生点;j为不同的发生程度，其中1代表轻度发生、2代表中度发生、3代表重度发生;Sij为i点j发生程度灾害面积;Si为i点的监测面积。

3．3 物候因子筛选

典型相关分析法是利用综合变量对之间的相关关系来反应2组变量间的整体相关性的多元降维统计方法［17］，可用来筛选与油松毛虫灾害发生情况相关的衍生物候因子。考虑到单一统计方法在相关性分析中存在的缺陷，同时运用主成分分析法对与油松毛虫灾害发生的相关气候因子作分析，经过2种方法的结合，得到主要物候因子。其一般过程为:首先对所有的衍生因子作正态性检测;其次用典型相关分析筛选衍生因子，并对显著性进行检验;最后用主成分分析法得出相关结果，并与典型相关分析结果做比较分析。用SPSS 19．0和Correl函数进行数据分析。用典型相关分析筛选物候衍生因子的原理过程详见文献［17］。

4 结果与分析

4．1 衍生的物候因子

根据相关研究和本研究设定的条件，共衍生出可能和油松毛虫灾害发生相关的物候因子73个。这些因子有温度、湿度、降水、日照时间和风速等(表3)，包括了油松毛虫在北方区域中完成生活史不同阶段的基础物候因子和极端物候条件，如Ⅰ期平均温度(X6)指前一年的8、9、10月份的平均温度，此期间是油松毛虫幼龄幼虫进食期，适宜的温度是其能否完成各个期的基础条件，而期间的湿度、风速和最大日降雨(X43、X59、X62等)等条件在很大程度上影响了成活率、扩散范围［3］19－25，可直接造成其种群分布和数量的波动，从而影响其灾害的暴发。

表3 区域尺度油松毛虫预报衍生因子

4．2 衍生因子典型相关性

用SPSS 19．0对衍生物候数据和油松毛虫灾害发生比率数据矩阵块分析计算，得到物候因子和油松毛虫灾害情况的均值和标准差，同时分析得出典型相关系数及其检验结果(表4)。由表4可见，第1、2典型相关系数分别为1．000和0．988，其统计检验达到极显著水平(P1＜0．001、P2＜0．001)。第3典型相关系数虽然也较大但其统计检验无法达到显著水平，故剔除。因此，选前2对典型变量来分析物候因子与油松毛虫灾害发生之间的相关关系。

表4 预报因子对油松毛虫灾害发生影响的典型相关分析

由于物候因子和虫情数据的量纲不同，为了降低分析的误差，故在相关分析中对变量进行标准化处理。标准化处理后的第1组典型相关变量为:

在第1组典型相关变量中，U1主要受到温度因素(X6、X16、X26、X28、X30)、湿度因素(X2、X17、X22、X29、X33、X35、X37)和风速因子(X9)的影响。由偏相关系数分析可得Ⅲ期平均温度(X16)、Ⅰ期平均温度(X6)及构成Ⅲ期各月份平均温度(X26、X28、X30)为其主要影响因子。其中，X16、X26、X28、X30呈正相关，X6呈负相关。即在一定范围内油松毛虫上树活动期温度越高，灾害发生的概率就越大;反之，油松毛虫灾害发生的概率就越小。V1主要受到轻度发生比率y1的影响。在上树期温度因子的作用下，油松毛虫灾害的发生主要为轻度发生。

标准化处理后的第2组典型相关变量为:

在第2组典型相关变量中，U2主要受到Ⅳ期平均湿度X22及Ⅳ期各月份平均湿度(X33、X35、X37)的影响，且呈负相关性。即Ⅳ期的湿度越高，油松毛虫灾害发生的概率就越低，Ⅳ期的湿度越低则其灾害发生的概率就会越高。因此，将Ⅳ期平均湿度称为影响油松毛虫灾害发生的第2因子。V2则主要决定于轻度发生比例和中度发生比例。此外，典型相关表达式中也反映出Ⅰ期风速对油松毛虫灾害发生的影响较大且呈现正相关。

综合第1、第2典型相关变量，Ⅲ期平均温度、2月份平均温度、3月份平均温度、4月份平均温度、Ⅰ期平均温度、Ⅰ期平均风速构成影响油松毛虫灾害发生的主要物候因子。

4．3 衍生因子的主成分分析

典型相关分析可得出和油松毛虫灾害发生最相关的物候因子，但同时也有可能漏掉一些相关因子而使得筛选的相关因子不全面，故对衍生因子运用主成分分析法对其进行补充。主成分分析法在保证信息可靠性的前提下，可通过少数因子反映全部因子的大部分信息(＞85%)。而当分析过程中设定抽取主成分和衍生因子相等时，可保证筛选物候因子的全面性，故而可补充典型相关分析的不足。对上述73个因子运用SPSS 19．0做主成分分析，当选取前11个主成分时，累计贡献率达到86．46%(表5)，故可用前11个成分来代替分析原来的73个因子。

表5 前11个主成分方差

将前11个主成分的载荷矩阵利用SPSS 19．0的Transfrom计算得到11个主成分的特征向量。结果表明，第1主成分中变量系数绝对值在0．30～0．45的变量共33个，是对温度和风速的综合。温度与第1主成分呈正相关，相关顺序由大到小依次为X1、X6、X16、X73、X60、X21、X11;风速和第1主成分呈负相关，相关顺序由大到小依次为X4、X9。第2主成分集中反应了相对湿度的影响，其变量系数为0．30～0．37，由大到小依次为X2、X44、X12、X22、X17。第3主成分到第11主成分中未集中反应某种因子的相关性，但反应了多种因子的综合影响，其偏相关系数都较小(＜0．1)。故在分析影响油松毛虫灾害发生的相关因子时，以第1和第2主成分为主。

在衍生物候因子的基础上，综合典型相关分析和主成分分析结果，选择温度因子中的年平均温度、Ⅰ期平均温度、Ⅳ期平均温度、Ⅱ期平均温度、Ⅱ期平均最高温、9月份最高温，风速因子中的年平均风速、Ⅰ期平均风速，湿度因子中的年平均湿度、9月份平均湿度、Ⅱ期平均湿度、Ⅳ期平均湿度、Ⅲ期平均湿度为我国北方油松毛虫灾害发生预测预报上的主要因子。

4．4 筛选物候因子的Correl函数

由典型相关分析和主成分分析筛选出和油松毛虫灾害发生关系最大的13个物候因子，并用Correl函数对这些因子与不同程度的灾害发生作相关性分析。

表6 筛选的物候因子与油松毛虫不同发生比例间的相关系数(r)

筛选物候因子和不同灾害程度的发生相关性变化范围较大(－0．74≤r≤0．93)。表6显示，相关性最强的是Ⅱ期平均最高温度X73和重度发生(r=0．93)，并呈正相关;相关性最弱的是Ⅱ期平均湿度X12和轻度发生(r=－0．01)，并呈负相关;同时，年平均温度(X1)和轻度发生比例、Ⅱ期平均温度(X11)和重度发生比例、Ⅲ期平均湿度(X17)和重度发生比例、9月份最高温(X60)和重度发生比例都有较强的相关，其相关系数都在0．75附近。

综上分析:与油松毛虫轻度灾害发生呈正相关的主要物候因子由强到弱依次为X1、X17、X60、X21、X4、X9。与中度灾害发生呈正相关的主要物候因子由强到弱依次为X11、X6、X9。与重度灾害发生呈正相关的主要物候因子由强到弱依次为X73、X60、X11、X6、X1、X44，呈负相关的物候因子由强到弱依次为X17、X2，即越冬期高温和上树期低温对油松毛虫重度灾害的发生有重要影响。

5 结论与讨论

研究表明，在众多的物候因子中，温度、风速、湿度因子是影响我国北方油松毛虫灾害发生的主要物候因子。这些因子包括基础物候因子和极端物候因子，基础物候因子有年平均温度、Ⅰ期平均温度、Ⅱ期平均温度、Ⅳ期平均温度、年均风速、Ⅰ期平均风速、年平均湿度、9月份平均湿度、Ⅱ期平均湿度、Ⅲ期平均湿度、Ⅳ期平均湿度，极端物候因子Ⅱ期平均最高温、9月份最高温。这些物候因子对油松毛虫不同程度灾害发生响应也不同。轻度发生与年平均温度和上树期的平均相对湿度相关性较强。中度发生与羽化—产卵—孵化—幼龄幼虫期的平均温度和越冬期平均温度相关性较强。重度发生与越冬期最高温、9月份最高温和上树期平均相对湿度相关性较强，即暖冬和早春干旱均会提高油松毛虫重度灾害发生的可能性。本研究结果与前人众多的物候因子方面的研究［10，13－14，18－20］各有异同。就松毛虫而言，薛贤清［18］通过对25个县、市马尾松毛虫发生情况及气象因子分析，得出和马尾松毛虫灾害发生的主要气象因子由强到弱为平均气温、降水量、相对湿度;而谢伟忠［19］则认为3、8月份的平均气温，4、5、6月份的平均湿度及其温湿度系数是影响马尾松毛虫灾害发生的主要因子，5—7月份连续的暴雨可减少松毛虫的危害。刘志明等［20］研究结果表明，相对湿度是影响落叶松毛虫灾害发生的首要因子。对其他虫种类而言，Flower et al．的研究证明前一年的低湿对云杉卷叶蛾的发生有促进作用，而高湿则很少导致灾害的发生［13］;Bentz et al．则发现温度的升高对红脂大小蠹和山松大小蠹行为造成显著影响［14］;Liang et al．报道冬季高温、夏季少雨可提高山松大小蠹的死亡率［12］。可见，气候条件尤其是温湿度的变化会对不同林业昆虫造成不同的影响。综上所述，虫种不同，气候条件对其的作用亦不尽相同，但对于松毛虫灾害的发生都受到温湿度的显著影响。

油松毛虫灾害发生与温度、相对湿度和风速相关。油松毛虫在北方以2～4龄幼虫越冬，越冬期和早春复苏期其幼虫的自我调节能力及抗环境胁迫能力都较差，越冬期高温可提高越冬代的成活率。高湿环境利于昆虫病害的流行，从而在很大程度上对松毛虫种群密度的增长起到抑制作用。在高湿环境下，赤松毛虫(D．spectabilis)CPV在林间使用良好［21］，思茅松毛虫(D．kikuchii)CPV病毒感染率可达到97%以上［22］。所以，经过越冬期复苏的幼虫生命力和抵抗力不强极易感染病害，而上树期的低湿可抑制昆虫病害的流行提高虫口基数。不同的风速对油松毛虫的灾害发生的影响主要体现在幼龄幼虫期，过大的风速如日风速＞5 m·s－1或月风速＞10 m·s－1则会造成幼虫的大量死亡而降低幼虫成活率［3］24，而适当的风速会对松毛虫的传播起到促进作用［23］。

松毛虫灾害的暴发不仅和种群自身遗传有关，还取决于其栖息地诸多环境因素。而众多环境因素中，气象因子的可量化性使其在模拟预测中地位突出。森林昆虫完成生活史的各个时期，都要以适宜的气候环境作为依托。故本研究在考虑油松毛虫和气象因子的相关性时，为了保证其生长发育的各个时期和对应的气象因子的最大同步性，最大程度上结合了我国北方区域实际气候特征和油松毛虫在该区不同时间的发育情况，将前一年8—10月份定义为Ⅰ期，11月份至次年1月份为Ⅱ期，2—4月份为Ⅲ期，5—7月份为Ⅳ期。同时，将各月的气象因子以月为单位分别处理，并结合油松毛虫生长发育不同时期的不同敏感性［5］，加入其极端气候变量。从而可更好地筛选影响其灾害发生的主要因子，除上述筛选的主要物候因子外，降雨和光照也有一定作用。降雨的影响取决于降雨的强度和降雨的时间，月总降雨量大于500 mm或日降雨量超过100 mm虫口密度会明显降低［3］22，［6］169，［23］;松毛虫对光周期的敏感性会引起松毛虫滞育［24］，有利于其越冬，从而保证其虫口密度并对灾害发生起积极作用;反之，则起到消极作用。但由于缺乏油松毛虫长期的监测数据，只能用2002—2011年相对完整的77组历史数据进行分析，其结果可能存在一定的偏差，为了更好地了解研究气象因子对森林病虫害的影响，应加强监测样点布置的代表性和监测数据量的积累。建立气候因子与有害生物灾情之间的关系模型将是未来林业有害生物监测预报的研究重点内容之一。同时，以温度和降水变化为主要特征的气候变化背景下，结合气象关键因子和其他的环境因子对林业有害生物灾害暴发［12］的评价和预报也是研究的重要内容。

［1］段永春，袁洪刚，张永亮，等．山东茶树虫害的种类及主要品种的发生规律研究［J］．中国农学通报，2010，26(11):284－289．

［2］王志伟，张东霞，马雅丽，等．山西山西省冬小麦主要病虫害气象等级预报模型［J］．中国农学通报，2010，26(11):267－271．

［3］陈昌洁．松毛虫综合管理［M］．北京:中国林业出版社，1990．

［4］曾菊平，戈峰，苏建伟，等．我国林业重大害虫松毛虫的灾害研究进展［J］．昆虫知识，2010，47(3):451－459．

［5］Yamamura K，Yokozawa M，Nishimori M，et al．How to analyze long-term insect population dynamics under climate change:50-year data of three insect pests in paddy fields［J］．Population Ecology，2006，48(1):31－48．

［6］侯陶谦．中国松毛虫［M］．北京:科学出版社，1989．

［7］韩瑞东，何忠，戈峰，等．马尾松毛虫黄河流影响松毛虫种群动态的因素［J］．昆虫知识，2004，41(6):504－511．

［8］费海泽，王鸿斌，孔祥波，等．马尾松毛虫发生相关气象因子筛选及预测［J］．东北林业大学学报，2014，42(1):136－140．

［9］张玉书，冯锐，陈鹏狮，等．松毛虫发生期与气象条件关系［J］．中国农业气象，2004，25(3):26－28．

［10］周广学，张国林，梁群，等．气象条件对油松毛虫的影响及其预测模型的构建［J］．东北林业大学学报，2012，40(11):131－134．

［11］陈绘画，朱寿燕，崔相富，等．基于人工神经网络的马尾松毛虫发生量预测模型的研究［J］．林业科学研究，2003，16(2):159－165．

［12］Liang L，Hawbaker T J，Chen Y L，et al．Characterizing recent and projecting future potential patterns of mountain pine beetle outbreaks in the southern rocky mountains［J］．Applied Geography，2014，55:165－175．

［13］Flower A，Gavin D G，Heyerdahl E K，et al．Drought-triggered western spruce budworm outbreaks in the Interior Pacific Northwest:A multi-century dendrochronological record［J］．Forest Ecology and Management，2014，24(3):16－27．

［14］Bentz B J，Régnière J，Fettig C J，et al．Climate change and bark beetles of the western United States and Canada:direct and indirect effects［J］．BioScience，2010，60(8):602－613．

［15］Musolin D L．Insects in a warmer world:ecological，physiological and life-history responses of true bugs(Heteroptera)to climate change［J］．Global Change Biology，2007，13(8):1565－1585．

［16］Kroschel J，Sporleder M，Tonnang H E Z，et al．Predicting climate-change-caused changes in global temperature on potato tuber moth(Phthorimaea operculella)distribution and abundance using phenology modeling and GISmapping［J］．Agricultural and Forest Meteorology，2013，17(6):228－241．

［17］田兵．典型相关分析及其应用［J］．阴山学刊，2013，27(4):11－22．

［18］薛贤清．在马尾松毛虫测报中应用逐步回归电算方法的研究［J］．林业科学，1984，20(1):42－49．

［19］谢伟忠．温湿度对广东阳春市马尾松毛虫种群的影响［J］．中国森林病虫，2005，24(4):17－19．

［20］刘志明，晏明，张旭东，等．气象因子对落叶松毛虫的影响分析［J］．灾害学，2002，17(4):34－37．

［21］陈昌洁，王志贤，刘革，等．赤松毛虫质型多角体病毒的引进和利用研究［J］．林业科学研究，1988，1(1):14－27．

［22］曾述圣，杨桂，杨中学，等．思茅松毛虫质型多角体病毒的应用［J］．森林病虫通讯，2000，2(1):28－29．

［23］孙思清，邵守英，刘鸿岩．台风暴雨对赤松毛虫虫口密度的影响［J］．山东林业科技，1997(6):21－23．

［24］李兆麟，贾凤友．油松毛虫的光周期反应［J］．昆虫学报，1989，32(4):410－417．