Effects of Topography and Forest Vegetation Factors on Occurrence and Population Size of Overwintering Dendrolimus superans

研究地点位于内蒙古大兴安岭绰源林业局 (120°17′52″~121°40′19″ E，47°35′21″~48°37′50″ N)，总面积约为3.13 × 10⁶ hm²，森林覆盖率83.7%，海拔750～1 450 m。该地区属大兴安岭山脉南麓，寒温带大陆性气候，夏季气候凉爽湿润，冬季冰雪覆盖期长，年平均气温-2.1℃，年均降水量461.9 mm，无霜期87 d。该区域常见植物60科214属390种，具有寒温带针叶林生态系统的典型特征，主要建群种为兴安落叶松 (Larix gmelinii (Rupr.) Kuzen)和白桦 (Betula platyphylla Suk.)，伴生乔木主要包括蒙古栎 (Quercus mongolica Fisch.ex Ledeb.)、黑桦 (Betula dahurica Pall.)、山杨 (Populus davidiana Dode)等，下层灌木主要包括杜香 (Ledum palustre L.)、越桔 (Vaccinium vitis-idaea L.)等。

2019年在绰源林业局所辖林场 (翠岭、梨子山、青岭、苏格河和育林林场)进行样地的选取，依据绰源林业局森林资源二类调查数据提供的各小班林组信息，筛选林组为近熟林的小班，随机设立142个20 m × 20 m不同林型的临时调查样地。调查每个样地的坡向、坡度、坡位、海拔。坡向分为平地、阴坡 (北坡和东北坡，337.5°～22.5°和22.5°～67.5°) 、半阴坡 (东坡和西北坡，67.5°～112.5°和292.5°～337.5°) 、半阳坡 (西坡和东南坡，112.5°～157.5°和247.5°～292.5°) 和阳坡 (南坡和西南坡，157.5°～247.5°)。坡度分为< 6° (平坡)、6°～15° (缓坡)、16°～25° (斜坡)和 >25° (陡坡)。坡位分为谷地、坡下、坡上和山脊。

2020年4月下旬，在落叶松毛虫幼虫上树为害前，采用五点取样法，在每个样地的4角和中心附近各选取2株兴安落叶松，共10株，清除目标树木胸径以下的侧枝和与其树干相连的灌木枝条，在每株树的胸径处环刮去树干的粗皮，用塑料胶带在此处粘闭合的胶带环，再往胶带上均匀涂抹粘虫胶，形成一个宽度15 cm的闭合隔离环。幼虫开始上树后，每隔3天记录闭合环上的幼虫数量，并清除杀死，直到连续3天未见有新虫上树，观察结束，统计每株的幼虫上树数量。测量每株树的胸径、树高，记录各样地林分密度、郁闭度、草本盖度等林地环境因子。其中，林分密度通过计数样地内胸径≥5 cm的林木株数换算得到。

林分密度指单位面积上的林木株数，可利用下式得到：

式 (1)中，S为林分密度 /(株·hm⁻²)，N为样地内林木株数，A为样地面积，本研究中A为0.04 hm²。

调查142个样地的松毛虫数量，统计发生落叶松毛虫样地中10株树木的落叶松毛虫幼虫总数，计算平均值作为该样地落叶松毛虫的发生量，依据《松毛虫监测预报技术规程》^[22]，利用落叶松毛虫的发生量将各个样地的松毛虫危害程度进行分级，分别为极轻度 (0～20 条·株⁻¹)、轻度 (20～40 条·株⁻¹)、中度 (40～70 条·株⁻¹)、重度 (>70 条·株⁻¹)危害。

利用二元Logistic回归模型，将样地中是否存在落叶松毛虫作为因变量，以地形和森林植被因子为自变量，其中坡向、坡度等级、坡位和林分类型为定性变量称为项目，需将定性变量进行划分称为类目(表1)，并以指数函数(0, 1)进行定量处理，只要定性变量中有一个类目的影响效果显著(P<0.05)，则把该变量代入到模型中。通过判别预测结果，得到混淆矩阵，然后计算敏感性和特异性，做ROC曲线图，并计算AUC作为模型精度的预测指标。二元Logistic回归模型的基本公式为：

对式 (2)进行转化得出发生概率的函数与自变量之间的线性回归模型公式为：

式 (2)和 (3)中，P_i为落叶松毛虫发生的概率；y为因变量，表示该样地中落叶松毛虫是否发生，发生=1，未发生=0；x_i为自变量，表示第i种影响因素；β_i为自变量的回归系数。

在落叶松毛虫发生的126个样地中，对各定量变量进行分级处理，利用单因素方差分析（one-way ANOVA）探讨地形和森林植被对落叶松毛虫发生量的影响，并利用 LSD 检验进行多重比较。采用对数转化以保证其发生量的正态分布和方差齐性。

利用数量化理论Ⅰ将样地内落叶松毛虫的发生量作为因变量，以地形和森林植被因子为自变量纳入数量化模型中，其中，定性变量需以指数函数 (0, 1)进行定量处理，例如，将谷地定量为 (1, 0, 0)、上坡定量为 (0, 1, 0)、下坡定量为 (0, 0, 1)、那么山脊则定量为 (0, 0, 0)，并带入到线性模型中，因此，每个定性变量的最后一个类目 (山脊)是多余的参数，得分为0，并筛选出显著性影响因子 (P<0.05)，筛选方法与二元Logistic回归模型相同。数量化模型公式为：

式 (4)中：y_i为因变量，是第i块样地的落叶松毛虫发生量；$ {\widehat{b}}_{0} $为常数；$ {\delta }_{i\left(i,k\right)} $为类目的反应值；$ {b}_{jk} $为j项目k类目的得分值，(j=1，2，…，m；k=1，2，…，r_j)；x_j是定量变量的值。

方程模型由随机抽取的96个样本数据拟合而成，剩余30个未参加模型建立的数据用于该模型精度的检验和评价。模型精度检验采用决定系数 (R²)、均方根误差 (RMSE)对模型进行评价，以确定模型的可靠性和精度。本研究中二元Logistic回归模型、ROC曲线、单因素方差分析和多重比较在SPSS Statistics 26.0中完成，数量化理论Ⅰ是在自变量中包含定性变量的多元线性回归，通过Forstat 2.1完成^[23]，随机抽样和绘图采用Excel 2010软件。

坡向和草本盖度对落叶松毛虫的发生影响显著 (P<0.05)。将表2中常量和回归系数代入公式 (3)得到模型：

模型 (6)表明，落叶松毛虫发生在阴坡的概率最低，发生在阳坡的概率最高，且随着草本盖度的增加，发生概率增加。AUC为0.826，范围为0.717～0.939，说明模型的预测结果较好，精度较高(图1)。

Figure 1.  ROC curve of the occurrence probability of Dendrolimus superans

在落叶松毛虫发生的126个样地中，落叶松毛虫的发生数量受到地形和森林植被因子的显著影响 (P<0.05)。图2A表明，落叶松毛虫在阳坡的幼虫数量显著最多，在阴坡的幼虫数量显著最少 (P<0.05)。图2B表明，落叶松毛虫数量随坡度等级的增大而显著减少 (P<0.05)。图2C和图2D显示，落叶松毛虫在山脊的数量显著多于其他相应坡位，且在海拔<900 m的样地发生数量显著多于海拔>1 100 m的样地 (P<0.05)。图2E表明，草本盖度在60%～70%的样地，虫数显著少于另外3个草本盖度等级的样地。图2F表明，在林分密度大于1200株·hm⁻²的样地虫数显著少于另外3个林分密度等级的样地。

Figure 2.  Influences of topography and forest vegetation on the occurrence quantity of Dendrolimus superans

将表3中各变量的得分代入到公式 (4)中，得到落叶松毛虫发生数量与地形和森林植被因子的方程，模型为：

坡向、坡度等级、坡位、郁闭度和林分密度的影响效果显著 (P<0.05)，为影响落叶松毛虫发生数量的关键因子 (表3)。

利用这5项关键因子重新建立其与落叶松毛虫发生数量的方程，模型为：

落叶松毛虫的发生数量在阴坡最少，阳坡最多。随着坡度等级的增加，其数量随之减少。在山脊的落叶松毛虫发生量最多。森林郁闭度增加，或林分密度降低，落叶松毛虫的发生数量随之增多 (表4)。

模型 (7)和 (8)的R²分别为0.685和0.635，RMSE分别为0.936和0.964，且回归方程极显著 (P<0.001)，表明2个模型均较可靠，用5个关键因子可以达到全部地形和森林植被因子的拟合效果 (表5)。

利用模型 (8)计算得到未参加方程建立的剩余30个样地的落叶松毛虫发生数量预测值，与其实测值作比较，有6个样地的落叶松毛虫预测与实测危害等级不同，且多为危害等级为中度的预测。预测危害等级的准确率为80.0% (图3)。

Figure 3.  Predicted and measured occurrence quantity of Dendrolimus superans

坡向和草本盖度是决定大兴安岭南部落叶松毛虫是否发生的关键因子 (表2)，随着从阴坡到阳坡的光照增强和水分减少，落叶松毛虫的发生概率逐渐增加。这主要是因为温度的升高有利于虫卵的发育和幼虫摄食强度提高；而湿度的增加会打破落叶松毛虫体内水分平衡，同时使白僵菌等微生物的繁殖加快，提升其死亡率^{[20, 24]}。对西伯利亚松毛虫 (Dendrolimus sibiricus)和欧洲松毛虫 (Dendrolimus pini)的研究也发现，夏季高温和少雨造成的干旱生境会缩短这两类松毛虫的生命周期，是造成它们大规模暴发的主要诱因^{[15, 25-26]}。由此可以推测，随着气候变暖，本研究区域的落叶松毛虫暴发风险将持续提升。另外，草本盖度对该虫害发生概率的作用途径很可能是通过增加土壤呼吸速率，土壤温度和肥力的增加更有助于落叶松毛虫产卵，并且幼虫在枯枝落叶层生存以林内灌草为掩护，不易被天敌捕食和人工防治。本研究中发生概率模型的AUC值为0.826 (图1)，与于跃等基于大兴安岭南部落叶松毛虫建立模型的AUC值相近^[8]。然而，相对于Kumar等^[27]建立的苹果蠹蛾 (Laspeyresia pomonella Linnaeus)全球分布预测模型 (AUC=0.94)以及吴思俊等^[28]对云南松毛虫 (Dendrolimus houi)适生区预测模型 (AUC=0.93)略低，这可能由于本研究选择小班为单位，空间尺度较小、为避免多重共线性采用的变量也较少，因而预测精度会随之降低，如能综合考虑气候、天敌等更多维度的因素，模型精度可能会进一步提高。

除坡向和草本盖度外，坡度等级、坡位、海拔和林分密度也是影响落叶松毛虫发生数量的关键因子。其中，坡度越大落叶松毛虫的发生数量越少，这与对马尾松毛虫 (Dendrolimus punctatus)的研究发现类似，坡度较大时，林内小气候和微环境复杂，虫害在环境和天敌的作用下不易暴发^[29]。从山谷到山脊土壤含水量不断下降，因此落叶松毛虫更适合在干燥的山脊处生存繁衍 (图2C)。落叶松毛虫的发生数量在海拔小于900 m的样地显著多于海拔大于1 100 m的样地 (图2D)，这间接反映了落叶松毛虫对温度的适应性趋势^[3]。需要注意的是，虫害的数量多不完全等同于危害大。对更高纬度的西伯利亚松毛虫的研究发现，其对森林危害 （主要指虫害导致的树木死亡）最严重的地形条件主要为海拔400 m左右，坡度处于11°～13°的缓坡地带^{[6, 30]}。当林分密度大于1 200株·hm⁻²时 (图2F)，落叶松毛虫发生数量显著减少，而其数量在落叶松纯林和混交林之间无显著差异，此时资源稀释效应更加突显^[12]，也可能是由于降低的林内风速会阻碍该虫的远距离扩散和定向转移^[31]。落叶松是落叶松毛虫的优先寄主，其长势与抗虫能力一致^[10]，然而，本研究中落叶松胸径和树高并没有显著影响该昆虫的发生数量。这可能受选取样地时落叶松龄组均为近熟林所致，落叶松代谢出对害虫生长起抑制作用的植物单宁含量在不同林龄阶段存在显著差异^[13]。

在数量化模型 (公式7)中，筛选出坡向、坡度等级、坡位、郁闭度和林分密度这5个关键因子。其中，落叶松毛虫的发生数量随郁闭度的增加而显著增多，高的郁闭度使林内通风不畅，该虫尤其是幼虫多生存在背风区域，这可能由于成年雌蛾飞行能力受腹中卵粒数量影响，导致背风区域其飞行距离有限而繁殖在此环境下^[31]。徐海舟等在湘北地区的研究发现，思茅松毛虫 (Dendrolimus kikuchii)发生量随郁闭度的增加而减少^[32]，本研究与其不一致的原因可能归因于气候条件的差异性。

在预测和防治方面，应用马尔科夫链模型对虫害进行预测较为广泛，然而，这种方法只能预报虫情的大致趋势。有关松毛虫的研究表明，该模型精度为76.9%，且划分的虫害发生量等级越多，精度越低，预测条件至少需10年的虫情资料^[33]。而本研究运用数量化理论Ⅰ建立的模型，精度较高且减少了数据收集的工作难度。但实测与预测危害等级的样地出现了偏差，尤其是对虫害较高样地的预测。胡瑞瑞等^[5]运用数量化模型建立了立地因子与虫基指数的关系，其对预测昆嵛山腮扁叶蜂危害等级的准确率达到87.5%，本研究为80.0%。在模型中加入土壤因子，即考虑土壤质地对落叶松毛虫发生量的影响，有可能提高预测危害等级的准确率。虽然气候变暖在时间和空间上对松毛虫生活史和主要寄主的影响仍存在很大的不确定性^{[15, 25]}，但目前的防治重点应放在干旱时期海拔较低的阳坡和半阳坡区域，通过增加乔灌混交度降低林下草本盖度，从而最大程度的降低该类虫害暴发的风险。此外，本文仅着重分析了地形和森林植被对落叶松毛虫发生和数量的影响，而有关于落叶松毛虫幼虫不同生活史阶段对地形和森林植被的响应有待于长期的监测和研究。

本研究分析了地形和森林植被因子与落叶松毛虫发生概率和数量的关系，证实了立地条件和森林植被特征均对其发生概率和数量具有显著影响，并发现坡向和草本盖度为影响其发生概率的关键因子。根据坡向和草本盖度，做好发生概率的预测可以准确的为落叶松毛虫的防治节省大量资源。不同坡向、坡度等级、坡位、海拔、草本盖度和林分密度之间的落叶松毛虫发生数量存在显著差异。利用坡向、坡度等级、坡位、郁闭度和林分密度可以有效预测该区域的松毛虫发生数量。为防止落叶松毛虫暴发成灾，应对处于阳坡、平坡和山脊的森林区域重点防治，通过合理的补植、疏伐和定期除草等森林经营措施将虫害控制在较低水平。