全球气候变化对森林生态系统动态有着重要影响(Henttonen et al., 2014; Guillemot et al., 2015)。森林生长对气候变化如何响应，而人类对此应采取怎样的应对措施受到广泛关注。树木生长受树种、年龄等内在因素(Gerendiain et al., 2012)及外部环境条件如林分密度、立地因子和气候因子等影响。雷相东等(2009)发现林分密度、立地因子和单木竞争对树木生长均有显著影响。气候因子对树木生长的影响较为复杂，不仅生长季内的温度和降水等气候因子会影响树木生长(Wimmer et al., 1997; 李宗善等, 2010; Chen et al., 2012)，生长季以前的气候因子也会对树木生长产生影响(Lebourgeois, 2000; Linderholm, 2001; Mäkinen et al., 2003)。这种现象与在生长季之前的气候因子对土壤水分、树木的生理活动等方面的影响有关(Rolland, 1993; Jyske et al., 2012)。树木生长与气候因子的响应关系会受树种(Frank et al., 2005; 王丽丽等, 2005)、年龄(Szeicz et al., 1994; Carrer, 2004)、地理位置(Savva et al., 2006; Henttonen et al., 2014)等因素的影响。为保证树木生长对气候因子变化的敏感性，研究地多选在寒冷、干旱地区的林线处(kirdyanov et al., 2003; 于健等, 2016)。在林分内部树木生长的研究中，通常认为在一个区域内气候是较为稳定的，因而不考虑气候因子对树木生长的影响。然而在全球气候变化的背景下，气候的稳定性正在发生改变，有必要考虑到气候因子对森林内部树木生长的影响。

目前，人类活动干扰也已成为影响树木生长的因素。间伐是常用的营林措施，可以降低林内竞争(Millar et al., 2007; 贾忠奎等, 2012a; Magruder et al., 2013)，改善林内土壤、光照等条件(Barbour et al., 1994; 李国雷等, 2008; 贾芳等, 2009; 贾忠奎等, 2012b)，提高保留木对气候因子的适应能力(Misson et al., 2003; Kerhoulas et al., 2013; Lebourgeois et al., 2014; Oliver et al., 2014)，促进保留木生长(Koga et al., 1997; 雷相东等, 2005; Guller, 2007)。不同强度的间伐对保留木的促进作用也有所差异，通常，中、高强度间伐对树木生长的促进作用较强，而低强度间伐对树轮宽度的影响较小或不显著(马履一等, 2007; 莫日根等, 2013; 龚固堂等, 2015)。间伐对树木生长的影响会在伐后数年间逐渐增强后减弱(Misson et al., 2002; Koga et al., 2002; Mehtatalo et al., 2014)。研究不同间伐强度对林木生长的影响，有助于为特定林分找出合适的经营方案。

树木年轮记录了树木逐年的径向生长，具有分辨率高和易于获取等特点(吴祥定等, 1996)。利用年轮宽度数据研究树木生长时，由于同一样地或同一个体间存在相关性，因此不能满足独立同分布的假设。混合模型可通过规定不同的协方差结构来表示相关的误差，允许数据间存在相关性和方差异质性，可提高预测精度和解释随机误差来源，因此被广泛应用于树木生长的研究中(雷相东等, 2009; 李春明等, 2010)。长白落叶松(Larix olgensis)以其生长快、适应性强的特点，作为先锋树种，大量用于我国东北地区人工林的营造，其生长与收获受到广泛关注。本研究以长白落叶松为对象，建立树轮宽度指标的非线性混合模型，分析间伐与气候因子对树木生长的影响，为全球气候变化条件下的林业生产等营林活动、森林生态系统管理提供参考依据。

1 研究区概况

研究地位于吉林省汪清林业局金沟岭林场(130°05′—130°20′E，43°17′—43°25′N)，属长白山系老爷岭山脉雪岭支脉。低山丘陵，海拔300~ 1 200 m，平均坡度10°~25°。温带大陆性季风气候，年均气温3.9 ℃，1月份气温最低，平均最低气温－32 ℃，7月份气温最高，平均最高气温32 ℃；全年无霜期138天；年均降水量600~700 mm，多集中在7月。土壤种类以暗棕壤为主，平均厚度40 cm左右。该地区植被属长白山植物区系。

2 研究方法 2.1 样地设置和调查

本研究样地数据采自金沟岭林场的20块固定样地。其起源为1964—1967年营造的长白落叶松人工林，经多年演变，大部分成为长白落叶松云冷杉针阔混交林。林内除长白落叶松外，还有鱼鳞云杉(Picea jezoensis)、臭冷杉(Abies nephrolepis)、红松(Pinus koraiensis)、色木槭(Acer mono)、白桦(Betula platyphylla)、水曲柳(Fraxinus mandschurica)、椴树(Tilia tuan)和榆树(Ulmus pumila)等。20块样地按照区组试验设计，分为5个区，每区包含的4块样地在1987年分别按胸高断面积的百分比进行轻度(20%)、中度(30%)、重度(40%)下层疏伐及对照处理，分别标记为A1, A2, A3, A4, B1, B2, B3, B4, C1, C2, C3, C4, D1, D2, D3, D4, D5, E1, E2, E3, E4；其中D1和D2样地于1993年再次分别进行了强度为30%和20%的间伐。各样地概况与间伐情况见表 1。样地在1987—2012年间进行过13次调查，间隔期为2~3年，内容包括树种、胸径和立地因子等。对于其间未调查年份的数据通过线性插值补全(Filipescu et al., 2014)。

2.2 数据来源

年轮宽度数据来源为2016年夏季采集的长白落叶松生长锥年轮条。年轮条采样前，依据样地调查数据中的胸径转化为径阶将林木划分为3个等级：1级为样地内径阶不低于四分位点的林木，2级为径阶不足1级木但高于林分平均胸径所在径阶的林木，3级为低于林分平均胸径所在径阶的林木。通过简单随机抽样的方式从各样地内1，2，3级林木中分别抽取4株长白落叶松钻取年轮条。共采集243株长白落叶松林木的年轮条，其中，4条样品(1条采自2级木，3条采自3级木)断裂难以恢复。对剩余的年轮条通过Lintab 6树轮宽度测量仪测量年轮宽度，并将年轮宽度数据通过COFECHA程序进行交叉定年。交叉定年结果显示8条样品年轮宽度数据与主序列相关性不显著，确认测量无误后舍弃。最终231条样品年轮宽度数据通过交叉定年。保留的各样地年轮条样品信息见表 1。年轮宽度数据的80%用于建模，20%用于检验。

气候数据为样地附近气象站提供数据的线性插值结果(Shen et al., 2015)，包括1964—2010年间每月最高气温、最低气温、平均气温、降水量及生长季大于5 ℃积温等。

本研究使用的1987—2010年的样地、气候与树轮宽度数据。

2.3 模型建立

本研究采用非线性混合模型进行建模。以树轮宽度值(RW)作为因变量时，模型残差分布表现出明显趋势，对树轮宽度值的多种变型形式试验后，发现$ \sqrt {{\rm{RW}}} $结果最好，因此选择树轮宽度指标$ \sqrt {{\rm{RW}}} $作为因变量。模型形式如下：

$ \sqrt {{\rm{R}}{{\rm{W}}_{ijk}}} = f\left({\alpha, A} \right) + g\left({\beta, T} \right) + \gamma X + {b_i} + {b_{ij}} + {\varepsilon _{ijk}}。$

(1)

式中：$ \sqrt {{\rm{R}}{{\rm{W}}_{ijk}}} $为i号样地j号树k年的树轮宽度指标；f(α, A)为包括固定参数向量α和形成层年龄A的函数，表示因变量随形成层年龄变化的趋势；g(β, T)为包括固定参数向量β和表示间伐因素的自变量矩阵T的函数；γ为固定效应参数向量；X为不包含形成层年龄与间伐因素的自变量矩阵；b_i为样地水平随机效应向量；b_ij为树木个体水平随机效应向量；ε_ijk为随机误差项。

对于树轮宽度随形成层年龄变化的趋势，选用三参数Logistic方程模拟，形式如下：

$ f\left({\alpha, A} \right) = {\alpha _0} + {\alpha _1}/[1 + {\text{e}^{({\alpha _2} \times A)}}]。$

(2)

式中：α₀表示生长稳定时的树轮宽度指标；α₁为靠近髓心的树轮宽度指标；α₂为尺度参数。

间伐对树轮宽度的影响通常表现为在伐后数年间先增大后逐渐减弱的趋势(Misson et al., 2002; Mehtatalo et al., 2014)，因此在比较不同方程形式后选择公式(3)拟合树轮宽度指标受间伐影响并随伐后时间变化的趋势。不同间伐强度对树轮宽度指标的影响以哑变量的形式表现。

$ \begin{array}{l} g\left({\beta, T} \right) = ({\beta _{111}}{T_{11}} + {\beta _{112}}{T_{12}} + {\beta _{113}}{T_{13}}){t_1}/\\ \;\;\;{\text{e}^{({\beta _{21}} \times t_1^2)}} + ({\beta _{121}}{T_{21}} + {\beta _{122}}{T_{22}}){t_2}/{\text{e}^{({\beta _{22}} \times t_2^2)}}。\end{array} $

(3)

式中：T₁₁为第1次轻度间伐的哑变量，T₁₂为第1次中度间伐的哑变量，T₁₃为第1次重度间伐的哑变量，T₂₁为第2次轻度间伐的哑变量，T₂₂为第2次中度间伐的哑变量，所有哑变量均是为1，否为0。β₁₁₁为第1次轻度间伐对树轮宽度指标影响程度的固定效应参数；β₁₁₂为第1次中度间伐对树轮宽度指标影响程度的固定效应参数；β₁₁₃为第1次重度间伐对树轮宽度指标影响程度的固定效应参数；β₁₂₁为第2次轻度间伐对树轮宽度指标影响程度的固定效应参数；β₁₂₂为第2次中度间伐对树轮宽度指标影响程度的固定效应参数；β₂₁为第1次间伐对树轮宽度指标影响程度与持续时间的固定效应参数；β₂₂为第2次间伐对树轮宽度指标影响程度与持续时间的固定效应参数；t₁为第1次间伐后年数；t₂为第2次间伐后年数。

2.4 模型自变量选择与参数估计

影响树木生长的因子较多，为研究方便，将除形成层年龄与间伐因素以外的自变量分为4组：立地因子组、林分因子组、单木竞争组和气候因子组。各变量组包含变量及说明见表 2。考虑到应用，单木竞争组变量采用较为方便获取的与距离无关的指数。林分密度指数(SDI)通过公式(4)计算。

$ {\rm{SDI}} = {N_{\rm{s}}}{({D_{\rm{g}}}/{D_0})^a}, $

(4)

式中：N_s为林分每公顷株树；D_g为林分平均胸径；D₀为标准胸径；a为自然稀疏率；本研究中长白落叶松标准直径与自然稀疏率的取值与李春明等(2010)的研究一致：D₀=20 cm，a=1.58。

本研究自变量较多，且气候因子间存在相互作用，各自变量间可能存在共线性，影响参数估计的精确性。研究中利用方差膨胀因子(VIF)判断自变量间的多重共线性。仅保留VIF < 5且回归系数显著的自变量(雷相东等, 2009)。本研究中各自变量量纲不一致，为通过模型参数比较各自变量对因变量影响的大小，将表 2中自变量标准化，此时模型参数代表在自变量变化1个标准差时因变量的变化情况。

对于可能存在的误差的方差异质性和序列相关性，将优选出合适的方差函数(幂函数、常数幂函数)及协方差结构(自回归、滑动平均、自回归滑动平均)。模型的选择基于Akaike’s information criterion (AIC值)(Akaike, 1974)。对于嵌套模型，AIC能够辨识出较好的模型，其值小者为优。模型参数估计采用极大似然法。所有计算通过R软件的nlme包完成。

2.5 模型评价及检验

通过以下指标对模型拟合结果进行评价：决定系数R²、绝对偏差(Bias)、相对偏差(Bias%)、均方根误差(RMSE)和相对均方根误差(RMSE%)，见公式(5)-(9)(雷相东等, 2009)。

$ {R^2} = r_{{y_i}{{\hat y}_i}}^2; $

(5)

$ {\rm{Bias}} = \sum\limits_{i = 1}^n {({y_i} - {{\hat y}_i})/n}; $

(6)

$ {\rm{Bias}}\% = \sum\limits_{i = 1}^n {({y_i} - {{\hat y}_i})/\sum\limits_{i = 1}^n {{{\hat y}_i}} ;} $

(7)

$ {\rm{RMSE}} = \sqrt {\sum\limits_{i = 1}^n {{{({y_1} - {{\hat y}_i})}^2}/\left({n - 1} \right)} } ; $

(8)

$ {\rm{RMSE}}\% = \sqrt {[\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}/\left({n - 1} \right)]/(\sum\limits_{i = 1}^n {{{\hat y}_i}/n)} } } 。$

(9)

式中：y_i为观测值；$ {\hat y_i} $为预测值；$ r_{{y_i}{{\hat y}_i}}^2 $为观测值y_i与预测值$ {\hat y_i} $的相关系数；n为样本数。

3 结果与分析 3.1 模型结果

除形成层年龄和间伐因素外，将9个自变量经标准化后引入模型，9个自变量分别为林分密度指数(SDI)、对象木胸径与林分平均胸径的比值(D_r)、生长季大于5 ℃积温(D_gs)、夏季热湿指数(H_s)、冬季降水量(P_w)、7月最高气温(T_7max)、9月最低气温(T_{9 min})、5月降水量(P₅)和4月平均温(T_4mean)。最终模型见式(10)。模型的参数估计结果见表 3。除表示第1次强度为轻度的间伐造成影响的参数β₁₁₁不显著外，其他参数均在95%置信水平上显著(P < 0.05)。模型中加入样地与单木水平随机效应后，AIC值降低，因此将样地与单木水平的随机效应引入模型。协方差结构选择一阶自回归结构AR(1)，此时模型AIC值最小。

$ \sqrt {{\rm{R}}{{\rm{W}}_{ijk}}} = f\left({\alpha, A} \right) + g\left({\beta, T} \right) + \gamma X + {b_i} + {b_{ij}} + {\varepsilon _{ijk}}; $

其中：

$ \begin{array}{*{20}{c}} {f\left({\alpha, A} \right) = {\alpha _0} + {\alpha _1}/[1 + {{\rm{e}}^{({\alpha _2} \times A)}}], {\rm{ }}}\\ {g\left({\beta, T} \right) = ({\beta _{111}}{T_{11}} + {\beta _{112}}{T_{12}} + {\beta _{113}}{T_{13}}){t_1}/}\\ {{{\rm{e}}^{({\beta _{21}} \times t_1^2)}} + ({\beta _{121}}{T_{21}} + {\beta _{122}}{T_{22}}){t_2}/{{\rm{e}}^{({\beta _{22}} \times t_2^2)}}, }\\ {\gamma X = {\gamma _1}{\rm{SDI}} + {\gamma _2}{D_{\rm{r}}} + {\gamma _3}{D_{{\rm{gs}}}} + {\gamma _4}{H_{\rm{s}}} + }\\ {{\gamma _5}{P_{\rm{w}}} + {\gamma _6}{T_{{\rm{7max}}}} + {\gamma _7}{T_{{\rm{9min}}}} + {\gamma _8}{P_5} + {\gamma _9}{T_{{\rm{4mean}}}}{\rm{, }}}\\ {{\varepsilon _{ijk}} \sim (0, {\sigma ^2}{R_{ij}}), }\\ {{\sigma ^2}{R_{ij}} = {\sigma ^2}{\mathit{\Gamma} _{ij}}\left(\varphi \right), }\\ {{\mathit{\Gamma} _{ij}}\left(\varphi \right) = {\rm{AR}}\left(1 \right)}。\end{array} $

(10)

式中：D_r为对象木胸径与林分平均胸径的比值；D_gs为生长季大于5 ℃积温；H_s为夏季热湿指数；P_w为冬季降水量；T_7max为7月最高气温；T_{9 min}为9月最低气温；P₅为5月降水量；T_4mean为4月平均气温；γ₁, γ₂, γ₃, γ₄, γ₅, γ₆, γ₇, γ₈和γ₉均为待估参数；Γ_ij(φ)为序列相关函数；AR(1)为一阶自回归结构。

图 1为模型参数估计结果代入公式(1)所绘制曲线，显示形成层年龄在15~45年的时间段内，树轮宽度指标逐渐减小并趋于稳定。

在间伐效应方面，模型中第1次轻度间伐的参数没达到95%置信水平显著，表明该次间伐对树轮宽度的影响并不显著，其余间伐措施均对长白落叶松树轮宽度指标表现出显著影响。图 2为将模型参数估计结果代入公式(2)绘制的曲线，图中显示树轮宽度指标在不同间伐时间与间伐强度间有所差异：在第1次间伐中，重度间伐引起树轮宽度指标增量最多，中度间伐次之，轻度间伐最少，在伐后1~2年达到峰值，分别为0.038，0.029和0.014个单位，后逐渐减弱并持续至约伐后4年；第2次间伐中，轻度间伐引起的树轮宽度指标增量比中度间伐大，且在伐后4~5年达到峰值，分别为0.069和0.057个单位，后逐渐减弱约持续至伐后15年。

描述立地因子的备选变量均未对长白落叶松树轮宽度指标表现出显著影响；林分因子组中SDI与单木竞争组中D_r对长白落叶松树轮宽度指标表现出显著影响。气候因子组中有7个变量对长白落叶松树轮宽度指标表现出显著影响。其中生长季大于5 ℃积温和7月最高气温每升高1个标准差所引起树轮宽度指标的变化量最多，分别为-0.023 6和0.022 0个单位；夏季热湿指数的影响次之，每升高1个标准差引起树轮宽度指标降低0.010 0个单位；4月平均温、5月降水量、9月最低气温和冬季降水量对树轮宽度指标的影响相对较小。表 4给出了模型中的自变量标准化计算时所用到的均值和标准差。

3.2 模型评价与检验

将检验数据代入模型，计算出评价统计量，与建模数据拟合结果的评价统计量进行对比(表 5)，显示模型检验与拟合的各评价统计量基本一致。引入随机效应参数后，模型的决定系数(R²)有大幅提高，对于建模和检验数据而言，分别从0.37和0.33提高至0.73和0.73，表明样地及树木个体间的差异是随机误差的重要来源。绝对偏差(Bias)、相对偏差(Bias%)、均方根误差(RMSE)和相对均方根误差(RMSE%)在建模数据和检验数据的拟合结果中相差不大，认为模型稳定可靠。

图 3显示模型残差分布均匀，未见明显趋势，因此不考虑方差异质性。

4 讨论

本研究利用非线性混合模型气候因子及间伐对长白落叶松树轮宽度的影响。模型的因变量$ \sqrt {{\rm{R}}{{\rm{W}}_{ijk}}} $是树木年轮宽度的变型，可反映出树木逐年生长量的相对多少。影响长白落叶松生长的因素包括形成层年龄、间伐、林分密度、单木竞争和气候因子。描述立地因子的变量并没有对树轮宽度指标表现出显著影响，原因可能是各样地立地因子较为均一，相差不大。模型结果显示高密度林分内长白落叶松生长情况相对较差，而在林分中处于优势竞争地位的林木生长的更好，这与雷相东等(2009)和Filipescu等(2014)的研究结果一致。

间伐对长白落叶松的生长表现出先增强后减弱的促进作用(图 2)，这样的趋势也见于Misson等(2002)、Koga等(2002)和Mehtatalo等(2014)的研究中。这种促进效果受间伐强度影响，并且在2次间伐中表现出差异。第1次间伐时林龄为20年，间伐后树木生长加快，重度间伐下林木径向生长最快，中度间伐次之，轻度间伐与对照组差异不显著，与其他学者的研究结果一致(马履一等, 2007; 莫日根等, 2013; 龚固堂等, 2015)。第2次间伐中轻度间伐比中度间伐对树木生长表现出了更强的促进作用，而其他学者关于间伐的研究中通常为轻度间伐对树木生长的促进效果弱于中度间伐(马履一等, 2007; 莫日根等, 2013; 龚固堂等, 2015)，原因可能是本研究中2次间伐相隔时间较短，仅有6年，树木生长还没有受到林内资源的强烈限制。在这样的情况下，30%强度的间伐与20%强度相比，对林内环境的破坏更多，从而抵消了部分间伐对林木生长带来的好处，因此表现出中强度间伐比低强度间伐对树木生长的促进作用弱。第2次间伐对长白落叶松生长的影响比第1次间伐更强且持续时间更长。原因可能是第1次间伐时林龄较小，树木生长旺盛，间伐对林内竞争的缓解由于林木的快速生长而迅速减弱。随着林龄增长，第2次间伐时林木的生长速度减缓(图 1)，因此该次间伐对林木竞争的缓解可以持续更长时间。然而，本研究中仅2块样地进行了第2次间伐，样本量小，有偶然因素造成这些现象的可能。

本研究中对长白落叶松生长表现出显著影响的气候因子与其他学者的相关研究结果有所差异(于大炮等, 2005; 于健等, 2016; 陈力等, 2014)。原因可能是不同研究地生长季的起止时间不同，各时间段的气候因子对树木生长的生物学意义也不一致(Savva et al., 2006)。同时不同研究地的气候条件存在差异，树木生长的限制因子不尽相同(Mäkinen et al., 2003; Wilmking et al., 2004)。

模型结果显示温度对长白落叶松树木生长影响较大，其中影响程度最大的是7月最高气温和生长季大于5 ℃积温。7月是该地区一年中温度最高的月份，较高的最高气温会促进林内的蒸发散，加剧干旱胁迫，高温又会促进呼吸作用，加速营养物质消耗，不利于树木生长(Lebourgeois, 2000; Wilmking et al., 2004)。Wang等(2002)在研究中发现树木早材宽度对生长季早期积温敏感，且早材比晚材在树轮宽度的变化中贡献更大，因此，生长季积温的提高对树木生长有促进作用。此外，夏季热湿指数对长白落叶松生长也表现出较大影响，较高的夏季热湿指数意味着在树木生长最旺盛的夏季，温度相对较高而降水不充足，干旱条件会限制树木生长(Wilmking et al., 2004; Henttonen et al., 2014)。4月平均温、5月降水量、9月最低气温和冬季降水量的变化所引起的树轮宽度指标的变化量相对较小，但也表现出显著影响。温暖的4月有助于土壤解冻及激活形成层，影响木质部生长的开始和持续时间，延长生长季，促进树木生长(Savva et al., 2006; Jyske et al., 2012)。5月降水对长白落叶松生长表现出抑制作用，推测原因是此时正处于春季，已消融的积雪补充了土壤水分，因而5月的降水不是树木生长的限制因子。而降水的增多意味着阴天与太阳辐射减少，不利于树木新枝条的生长与光合作用，并且寒冷潮湿的环境会抑制落叶松的生长(Chen et al., 2012)。9月进入了生长季末期，温度较低，树木蒸腾作用减缓，水分对树木生长的限制减少，温度成为树木生长的限制因子，因此表现出9月最低气温升高会促进树木生长(史江峰等, 2006)。冬季降水会在林内形成积雪，可以隔离树根并降低冰冻深度，且积雪融化可以补充土壤水分，有助于树木抵御干旱，有利于其生长(Linderholm, 2001)。

5 结论

本研究表明对长白落叶松云冷杉混交林内的长白落叶松生长有显著影响的因素有：形成层年龄、林分密度、单木竞争、间伐和气候因子。在林龄20年左右，中高强度的间伐对长白落叶松树木生长有较强的促进作用，促进效果可持续4年。低海拔长白落叶松生长受大量气候因子综合作用，其生长的主要限制因子是7月最高气温和生长季大于5 ℃积温。生长于林分内部的树木对气候因子变化敏感。在全球气候变化的背景下，在对树木生长进行研究时，为使研究结果更精确，应当同时考虑气候因子的影响。

关于气候因子及间伐的交互作用，包括不同间伐措施对树轮宽度与气候因子的响应关系的影响和在不同气候条件下间伐对树轮宽度影响的变化在本研究中并没有很好地体现，将在以后的研究中加以完善。