近年来全球气候变化背景下大范围树木死亡逐渐引起人们的广泛关注^[1-2]。地区水平上大面积的树木死亡不仅影响碳平衡, 还影响森林生态系统服务功能^[3]。我国的三北地区气候也正经历着暖干化变化^[4]。该地区三北防护林经过40多年的建设, 为治理风沙、改善京津冀环境质量发挥了重要作用^[5]。但是近十年来, 杨树防护林出现了大面积的枯梢和死亡现象。枯梢林分面积接近80%, 枯死和濒死的约占总面积的1/3^[6], 防护功能明显下降。水分是该区域植物生长的主要限制因子^[7], 同时粗放经营也是导致防护林生长发育不良甚至衰退的重要原因^[4]。研究杨树林的死亡原因和影响因素, 对全球气候变化背景下三北防护林的林分经营管理是十分必要的。

国内外很多学者从造林方式^[8]、林分结构^[9-11]、立地条件^[12-13]、气候环境条件^[14-16]等多方面对树木死亡的原因开展了大量研究。Qiu等^[9]建立了林木相对断面积、胸径、密度、土壤理化性质和气候因子与树木死亡关系模型, 认为林分密度过大, 或者林木相对断面积增长量达不到一定值时, 树木死亡率增加。陈清等^[17]研究了林龄和林分断面积对树木死亡的影响。Fan等^[18]研究了胸径、密度、坡度、坡向、海拔等因子与树木死亡的关系, 发现胸径是影响树木死亡最重要的因子。De Toledo等^[11]总结了影响树木死亡的土壤和地质因子。以上研究表明影响树木死亡的原因有很多, 包括胸径、密度、林木相对断面积、土壤理化性质等。然而不同径级的林木, 导致其死亡的原因是否相同尚不清楚。

本文以张北坝上地区的小叶杨(Populus simonii Carr)为研究对象, 以林木胸径、林木相对断面积、林龄、林分密度、不同土层土壤含水量和土壤容重等11个因子作为输入变量, 运用随机森林建立不同径级树木死亡模型, 对影响树木死亡的变量重要性进行评估, 揭示不同径级杨树死亡的主要原因, 这对于今后坝上地区防护林的恢复重建和经营管护具有重要的意义。

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

试验地点位于河北省坝上地区的张北县, 属温带大陆性气候, 多年平均气温3.2℃, 极端最高气温30.8℃, 极端最低气温-30.3℃, 年降水量不足400mm。每年约有50—70d风力达到六级以上, 以偏西风和西北风为主, 土壤风蚀严重。地带性土壤为栗钙土和沙壤土。该地区杨树防护林自20世纪60年代开始营造, 小叶杨为主要树种。近年来, 杨树出现大面积的枯梢枯死现象。乔木树种还有榆树、落叶松等, 主要灌木有沙棘、柠条等。

1.2 数据来源与处理 1.2.1 数据来源

张北的防护林以杨树纯林为主。2016年7—8月和2017年7—8月, 在对张北全县18个乡镇进行踏查的基础上, 目测立地类型、胸径、树高、林龄、密度等主要调查因子后, 以此为轮廓设置了100块20m×20m的标准地进行调查, 基本上涵盖了张北防护林不同立地类型、林分密度、林龄、胸径分布范围等。当调查对象为带状防护林时, 为减少边缘效应, 设置矩形标准地, 并根据实际情况改变标准地的大小, 使得标准地内的林木特征和立地条件尽量一致。

记录标准地的地理位置、海拔、地形等因子。由于张北的防护林都在平坦区域, 因此没有记录坡度、坡向、坡位等信息。对标准地内的杨树进行每木检尺, 记录胸径、树高、林木枯梢情况、林木死亡情况等。杨树防护林大多为同龄林, 因此每块标准地选取3株标准木, 应用生长锥测定林木的年龄, 取其平均值代表该林分的年龄。在每个标准地内按照一定方向选取3个有代表性的位置挖掘土壤剖面, 分别在每一剖面按照0—10、10—20、20—60cm采集土壤样本, 并把同一层次土壤样本进行混合。利用烘干法测定土壤含水量, 环刀法测定土壤容重、孔隙度等土壤物理性质。由于土壤含水量是变化值, 因此对样地进行了多次调查并取其平均值。该地区降水频率和强度都较小, 为研究一般情况下土壤含水量对树木死亡的影响, 调查时间选择在长期未降水的时期, 避免降水对土壤含水量造成的影响。

1.2.2 数据处理

通过实地调查发现, 张北的杨树林以成熟林和过熟林居多, 即大径级林木较多。为了使不同径级的林木样本量相近, 去除掉部分大径级的林木和部分生长健康无枯死木的林分, 只保留部分生长健康的林分, 最终选取26块样地, 基本涵盖了张北地区杨树林的径级范围, 并考虑到不同的立地类型和林分密度等因子。

研究显示与树木死亡有显著相关性的因子有林木径级大小、林龄、密度、土壤理化性质、立地等^{[9, 18]}。张北的杨树防护林分布于平坦区域, 海拔的分布范围较窄。研究表明张北防护林的退化死亡情况与土壤养分没有明显的相关关系^[6]。因此选择表征林龄、林木径级大小、林木疏密程度、土壤物理性质的因子。表征林木径级大小的因子有胸径和林木相对断面积, 表征林木疏密程度的因子有林分密度和林分断面积。土壤物理性质包括0—20cm土壤含水量、20—60cm土壤含水量、土壤平均含水量、0—20cm土壤容重、20—60cm土壤容重、土壤平均容重。林木相对断面积即林木断面积与林分平均断面积的比。

应用主成分分析(PCA)的方法对11个解释变量降维, 经过降维后得到4个主成分, 累计方差贡献率为87%。PC1主要由土壤含水量和土壤容重决定, 说明土壤物理性质；PC2主要由林木相对断面积决定；PC3主要由密度决定；PC4说明林龄。将林木分为5个径级(< 10cm、10—15cm、15—20cm、20—25cm、>25cm), 利用逻辑斯蒂广义线性混合效应模型分别不同径级建立树木死亡与4个主成分因子的回归模型。因变量为二元分类变量, 枯死记为1, 无枯死记为0。由于空间自相关性的影响, 距离相近的林木其生长与死亡可能互相影响, 因此将林木所属的样地号作为模型的随机效应。固定效应系数b表示模型中固定效应对因变量影响大小的参数, 将固定效应系数和模型显著性相近的划分为一类, 以此作为划分径级大小的标准。

应用随机森林研究解释变量对树木死亡的相对重要性, 建立不同径级树木死亡预测模型。

所有的统计分析应用R软件实现。

1.3 研究方法 1.3.1 逻辑斯蒂广义线性混合模型

设树木死亡的概率为P(二项分类因变量Y=1), 树木无枯死的概率为1-P(二项分类因变量Y=0)。对P进行Logit变换,

(1)

概率P与解释变量的关系为

(2)

1.3.2 随机森林算法

随机森林算法是由Breiman于2001提出的一种基于分类回归树的机器学习算法。如果因变量有n个观测值, 有k个解释变量与之相关, 它利用Bootstrap重抽样的方法从n个原数据中随机抽取n个观测值, 选择部分变量作为分类树节点, 从而构建几百甚至几千个分类树组合成随机森林, 将重复程度最高的树作为最终结果。利用随机森林算法可以对解释变量的相对重要性进行评价, 对解释变量进行预测。平均准确率降低度是衡量将一个变量的取值变为随机数后随机森林预测准确性降低程度的指标, 该值越大, 说明该变量的重要性越大。

2 结果与分析 2.1 不同径级树木死亡与环境因子关系

根据环境因子对各径级树木死亡影响的相似性, 可将树木径级划分为 < 10cm、10—25cm、>25cm(表 1)。各径级林木死亡率与土壤含水量、土壤容重(PC1)和林龄(PC4)的关系较弱。径级为10—25cm的林木, 其死亡率与林木相对断面积(PC2)具有极显著的正相关关系(P < 0.001), 与密度(PC3)具有显著的负相关关系(P < 0.05), 固定效应系数分别为1.011和-0.499。径级 < 10cm和>25cm的林木, 死亡率与解释变量之间无显著的相关关系。径级为10—25cm时, 固定效应解释的变异最大, R²为0.297。不同径级的回归模型中, 固定效应解释的变异小于随机效应解释的变异。

2.2 随机森林模型预测精度

运用随机森林分别不同径级建立树木死亡预测模型。表 2的混淆矩阵显示了随机森林模型的预测结果。不同径级模型预测的总体误判率为10.7%—31.25%, 解释了树木死亡70%—90%的变异。其中, 模型对于杨树防护林无枯死的误判率为4.8%—13.1%, 对于枯死的误判率为30.2%—44.8%。预测精度较为满意。

2.3 不同径级树木死亡影响因子重要性评估

分别不同径级, 运用随机森林评估11个影响因子对树木死亡的相对重要性, 并比较哪些变量的影响是相似的。由图 1可知, 胸径 < 10cm时, 林木相对断面积对死亡的影响远远大于其他因子, 其次为胸径和0—20cm土壤水分, 密度和林分断面积的影响相对较小。当林木相对断面积 < 0.22或>0.45时, 树木死亡率较大(图 2)。

胸径为10—25cm时, 密度对树木死亡的影响最大, 其次20—60cm土壤水分、林木相对断面积、林分断面积对树木死亡的影响相似。当密度>600株/hm²时, 树木死亡率随着密度的增加明显上升(图 2)。

胸径>25cm时, 20—60cm土壤水分和林分断面积对树木死亡的影响最大且相近, 其次为林分密度、林木相对断面积、0—20cm土壤容重的影响相近, 胸径对树木死亡的影响最小。当20—60cm土壤质量含水量>5%时, 树木死亡率明显下降。

3 讨论 3.1 不同径级树木死亡对环境因子的响应

本文的研究表明林木相对断面积、密度和土壤水分对该区域树木死亡的影响较大, 可能与该区域林木可获得的资源及受到的大风、冻害等恶劣的自然环境有关。林分中每株树木获得的资源一方面由其获得资源的能力即竞争力决定, 一方面由环境可提供的资源多少决定。胸径和林木相对断面积较大时, 林木竞争资源的能力较大^[19]。林分密度和林分断面积较大时, 环境可提供给每株树木的资源较少^[20]。本文的研究显示不同径级的林木, 林木相对断面积对死亡的影响大于胸径对死亡的影响, 说明林木相对断面积可以更好的描述林木之间的竞争力, 相对越大的树木竞争力越强, 这与Luo等^[21]的研究结果相近。

不同径级的林木, 各影响因子对树木死亡的影响权重不同。胸径 < 10cm时, 林木相对断面积对树木死亡的影响最大, 其次为胸径和0—20cm土壤水分, 密度和林分断面积的影响最小。说明该阶段树木大小是影响树木死亡最主要的因子。Zhu等^[12]的研究也发现树木大小对小径级树木的死亡有显著影响。本研究中当0.23 < 林木相对断面积 < 0.43时, 树木死亡率最低(图 2)。当林木相对断面积 < 0.23时, 树木死亡率较高。坝上地区属半干旱地区, 水资源严重缺乏, 当林木相对较小时, 一方面其竞争资源的能力较弱, 包括水资源、土壤资源、光照资源等。杨树属于速生物种, 木材密度较低^[22], 若得不到充足的水分、营养物质和光照, 将难以维持碳平衡, 根系和主干的碳储存减少, 树干密度降低, 抵御病虫害的能力随之降低, 死亡率增加。另一方面, 周围高大的林木会造成阴影, 不利于矮小的林木获得光照。同时, 有研究表明当树木年均断面积增长量达不到一定值时, 死亡率会增加^[9]。本研究中当林木相对断面积>0.43时, 林木死亡率增大。原因可能是坝上地区大风天气频繁, 相对较高的林木更易受大风冻害的影响, 大风造成落叶增加, 叶损失会减少树木光合作用产物, 增加非结构性碳消耗, 减少树木体内的碳储备, 当非结构性碳消耗到一定程度时会导致树木因碳饥饿而死亡^[23]。这与Coomes^[24]和De Toledo^[11]的研究结果相似, 大径级林木更易受大风、冻害等外界环境干扰, 死亡率相对较高。

胸径为10—25cm时, 密度对树木死亡的影响最大, 其次为20—60cm土壤水分、林木相对断面积、林分断面积对树木死亡的影响相近。说明林分所能提供的资源多少成为该阶段树木生长的主要限制因子。由于根系的生长, 这一阶段林木主要利用20—60cm土壤水分。苗博等人^[25]的研究中, 在生长季节胸径在该范围的杨树主要利用30—80cm土壤水分, 与本文的研究结果相似。随机森林的结果显示, 在这一阶段密度越大, 林木的死亡率越高。因为密度越大, 来自周围林木的资源竞争越强, 林分所能提供给单株林木的资源越少, 同时林木也会缺乏生长所需的空间。说明这一阶段的杨树都已成林, 需要进行疏伐。然而逻辑斯蒂广义线性模型的拟合结果显示, 这一阶段密度与树木死亡呈现负相关关系。可能是由于本研究是对树木死亡和4个主成分之间进行了拟合回归, PC3不能完全代表密度。同时也可能是由于该区域在冬春季节大风、冻害频繁, 当林分的密度较大时, 有利于林分抵御不利的自然条件。因此该径级范围内, 密度与树木死亡的关系还需进一步的探究, 从而得到树木生长的最适密度。

胸径>25cm时, 20—60cm土壤水分和林分断面积对树木死亡的影响最大且相近。说明该阶段林木主要利用20—60cm土壤水分, 且林分所能提供的水资源多少成为导致树木死亡的主要因子。在该阶段, 密度是导致树木死亡的第三个主要因子, 其影响小于林分断面积, 因为林分断面积是一个随着林分生长而变化的密度指标, 不仅能表示林木的拥挤程度, 还能反应林分对资源的需求。

在本研究中林龄对不同径级树木死亡的影响都较小, 影响树木生长的主要原因是水资源缺少, 大风、冻害等恶劣的自然环境, 以及人类不合理的经营。这与部分学者认为的坝上地区树木死亡是树龄老化的观点不一致^[26], 在以后需要对杨树是否达到老龄化这个问题进行进一步的研究。

3.2 模型预测精度

本研究运用随机森林算法建立的树木死亡预测模型可以解释树木死亡70%—90%的变异, Botter等^[27]、De Toledo等^[11]、Luo等^[21]建立的简单线性模型或回归树模型可以解释树木死亡20%—80%的变异, 说明本研究中建立的树木死亡模型拟合效果较好。逻辑斯蒂回归对自变量的多元共线性非常敏感, 要求自变量之间相互独立, 因此本研究在运用模型之前将11个影响因子做了主成分分析, 而随机森林算法可以评估所有变量的重要性, 不需要顾虑一般回归问题面临的多元共线性问题。

本研究只选取了表示林分结构和土壤性质的几个因子作为解释变量, 在以后的研究中可以加入气候因子。另外本文只研究了0—60cm土壤水分与树木死亡的关系, 以后可以增加更深层次土壤水分对树木死亡的影响。

4 结论

本研究基于随机森林, 分析了坝上张北地区导致不同径级树木死亡的主要影响因子, 并得到了关键因子的影响阈值。结果表明, 林木相对断面积、密度、林分断面积、土壤水分是影响该区域树木死亡的主要因子。林木胸径 < 10cm时, 当林木相对断面积在0.23—0.43之间时, 树木死亡率较低；胸径为10—25 cm时, 密度>600株/hm², 树木死亡率明显上升；胸径>25cm, 当20—60cm土壤质量含水量>5%时, 树木死亡率明显下降。不同径级树木死亡的主要影响因子不同, 但都表明水资源是影响该区域树木生长的关键限制因子。不同阶段的林木, 有其不同的适宜密度。基于水资源承载力的林分密度动态调控是该区域防护林经营的关键。要加强森林的经营管护, 且要遵循自然规律, 不同生长发育阶段的林木经营措施不同, 才能获得森林的可持续发展。