在森林生长模型中，由于径阶转移模型(矩阵模型)结构简单、易校检和检验、易与优化模型相结合(Lu et al., 1993)，自Usher(1966)第1次把预测动物种群未来年龄结构的矩阵模型应用于异龄林的生长动态分析以来，已得到广泛应用并在不断发展(Buongiorno et al., 1980; 1994; Solomon et al., 1986; Lin et al., 1996; Favrichon，1998; Kolbe et al., 1999; Sist et al., 2003; Liang et al., 2005; Hao et al., 2005a; Shao et al., 2006; Bollandsås et al., 2008)。矩阵模型的关键是转移概率矩阵的确定，包括进阶生长、进界生长、枯损和保留过程。根据转移概率将矩阵模型分为固定参数和可变参数2类(Lin et al., 1997) :前者将各径阶的转移概率视为常数; 后者则将转移概率表示为林分密度、期初胸径、立地因子以及树种和结构多样性等林分因子的函数，因为与林分状态有关，也称为与密度有关的非线性模型。如Buongiorno等(1980)将进界生长项由一个常数变成受林分因子约束的变量; Solomon等(1986)把生长表达为林分密度的函数。与密度有关的非线性模型能更准确地表达林分的真实生长情况，因此目前普遍被采用(Kolbe et al., 1999; Sist et al., 2003; Namaalwa et al., 2005)。在模型参数的估计方法上，因为3个转移概率间的相关性，也有采用联立方程组进行估计(Yang et al., 2008)。除了模拟不同经营方案对经济(如木材产量)和生态目标(物种和结构多样性)的影响外(Lu et al., 1993; Lin et al., 1996; Shao et al., 2006)，矩阵模型还和优化模型相结合，用于森林经营规划(Ingram，1996; Sist et al., 2003; Liang et al., 2005; Hao et al., 2005b; López et al., 2007)。我国学者也对矩阵模型进行了大量研究(郑耀军等，1987; 阳含熙等，1988; 宋铁英等，1989; 殷传杰，1989; 曾伟生等，1991; 谢哲根等，1993; 李荣伟，1994; 邵国凡等，1995; 王飞等，2004; 曲智林等，2006)。在经营模拟和优化方面，陆元昌等(2002)模拟了采伐对热带林的生长影响，并建立了采伐损伤模型; 王飞等(2005)、Shao等(2006)分析了不同采伐强度对红松(Pinus koraiensis)阔叶林年生长量、年收获量和恢复时间的影响; Hao等(2005a; 2005b)建立了长白山区混交异龄林的矩阵生长模型并确定了最优择伐方案; 郝清玉等(2006)建立了以矩阵生长模型为约束的林分结构及择伐周期优化模型。但研究存在以下问题:1)对于混交林，较少将树种分开进行建模; 2)模拟的经营目标仅限于木材生产; 3)转移概率模型求解没有考虑模型误差间的相关性。本文以落叶松云冷杉林为对象，建立分树种(组)的矩阵模型，并模拟不同采伐方案对木材生产、树种和林分结构多样性及碳贮量的综合影响，为多目标森林经营决策提供方法和依据。

1 研究地点概况

研究地区为吉林省汪清林业局金沟岭林场，位于吉林省汪清县境内东北部，所处的地理位置为130°05'—130°19'E，43°17'—43°25'N。属长白山系老爷岭山脉雪岭支脉，地貌属低山丘陵，海拔300 ~ 1 200 m，坡度一般在5° ~ 25°，个别陡坡在35°以上。林区属季风型气候，全年平均气温为3. 9 ℃左右，年平均降水量600 ~ 700 mm。其中5—9月的降水量为438 mm，占全年总降水量的80%。土壤主要是玄武岩中低山灰化土灰棕壤类型，平均厚度在40 cm左右。该区植被属长白山植物区系。

2 数据与方法 2.1 数据

用来建立模型的数据为金沟岭林场的20块固定样地。其起源为1964—1967年间营造的有部分针阔保留树种的人工落叶松林，经过多年的演变，大部分已成为落叶松云冷杉针阔混交林，具有天然林的部分特征，称为近天然林。以长白落叶松(Larix olgensis)、鱼鳞云杉(Picea jezoensis var. microsperma)和臭冷杉(Abies nephrolepis)为优势树种，其他树种有红松、色木(Acer mono)、水曲柳(Fraxinus mandshurica)、白桦(Betula platyphylla)、紫椴(Tilia amurensis)、枫桦(Betula costata)和春榆(Ulmus propinqua)等阔叶树种。样地面积在0. 077 5 ~ 0. 25 hm²之间。数据来自于1987—2007年的调查数据(1987, 1992,1997, 2002, 2007)，间隔期5年，其中2007年仅有4块调查样地。调查因子除每木检尺记录树种和胸径(≥5 cm)外，还包括立地因子如海拔、坡向、坡度，样地建立时基本概况见表 1。

2.2 研究方法 2.2.1 模型结构

模型结构采用Buongiorno等(1980)的形式，在考虑采伐和多树种的情况下，模型为:

式中: ，s = 1…m，m为树种组个数，i = 1…n，n为径阶数，y_s，i，t+1为s树种组t+1时刻i径阶的单位面积株数; Y_t = [y_s，i，t]，y_s，i，t为s树种组t时刻i径阶的单位面积株数; H_t = [H_s，i，t]，H_s，i，t为s树种组t时刻i径阶的单位面积采伐株数。当H_t为0时，为自然生长。

G_t为t时刻到t+1时刻的径阶转移概率矩阵，即

式中:a_s，n，t为第s树种组n径阶在t+1时刻保留在原径阶内的概率即保留率，b_s，n-1，t为第s树种组n-1径阶在t+1时刻向上生长到第i径阶的概率即进阶率。本研究中采用期初为第i径阶到期末上升到第j(j＞i)径阶的株数除以期初为第i径阶的株数，表示第i径阶到第j径阶的转移概率。

向量I_t = [I_st ]，I_st为s树种的进界向量，i_st为s树种从t时刻到t+1时刻进入起测径阶的株数:

本研究中划分为4个树种组即m = 4，分别为:落叶松、红松云冷杉、慢阔(色木、水曲柳、紫椴和枫桦)和中阔(白桦、春榆和杂木)。调查间隔期为5年，样地树木直径在5 ~ 60 cm之间，径阶宽度为5 cm，因此共划分为11个径阶即n = 11，保证树木只向上生长1个径阶，进界生长只生长到第1个径阶。最大径阶为60 cm，即60 cm以上的树木枯损概率为1。样地中的各自变量通过面积的线性关系化为每公顷的数值。

2.2.2 转移概率模型自变量的选择及参数估计

根据文献(Solomon et al., 1986; Liang et al., 2007; Lexerød，2005; Namaalwa et al., 2005; Woodall et al., 2005)，确定的备选自变量主要包括与径阶有关的因子[径阶中值、林分最小径阶株数、各树种(组)各径阶株数及其占林分总株数的比]、林分因子[林分断面积、林分株数、树种(组)断面积及其在林分中所占的比例]、多样性(树种和大小多样性)和立地因子，其中立地因子采用海拔、坡向、坡度及其交互作用来表示(Stage et al., 2007)。采用多元逐步回归来建立进阶生长、进界、枯损概率3个子模型，并对备选自变量各种变形(平方、导数、对数等)进行试验。自变量及其定义如表 2所示。

由于一些自变量间存在共线性，估计量的方差79林业科学47卷很大，相应地会产生较大的参数标准误，因此采用方差膨胀因子(VIF)来判断自变量间的多重共线性。一般认为，当VIF＞10时，有严重的共线性，此时，标记共线性较大的自变量，保留共线性弱而对因变量贡献大的自变量。此外，由于3个子模型的自变量部分相同，模型之间误差项可能相关，在多元逐步回归把自变量选出后，采用似不相关线性回归(seemingly unrelated linear regression) (唐守正等，2002)来进行参数的联合估计，并和单个模型普通最小二乘(OLS)回归估计进行比较。概率显著性水平取0. 05。

2.2.3 模型拟合优度检验

由于受建模样本数量限制，除第4期(2007年)调查的数据外，所有数据参加建模。对各径阶株数的预测值和实测值做卡方检验，检验分树种(组)和全林分进行。卡方值为:

式中:n为样本数，y_i为观测值，为预测值，i = 1… n，n为径阶数。概率显著性水平取0. 05。

2.2.4 采伐方案模拟

利用建立的矩阵生长模型模拟不同采伐周期和采伐强度对经营目标(见2.2.5)的影响。模拟林分各树种的起始径阶状态见图 1，即2002年所有样地径阶分布的平均值。模拟周期为50年，模拟分期设为5年，采伐强度(蓄积)分别为5%，10%，15%，20%，采伐周期分别为5，10和15年，最小采伐径阶为30 cm，采伐向量H即各径阶的采伐量按各树种的株数在各径阶所占的比例分配，采伐可分为只采针叶树、只采阔叶树和针阔叶都采3种情况，但因为本研究中样地阔叶树种(组)的比例很少，只采针叶树种(组)和其余2种采伐方法差别不大。因此只模拟采针叶树(按比例分配)情况。共得到13种采伐方案(包括自然生长)，用采伐周期+采伐强度来表示，例如0Y0P，5Y5P …15Y20P。

2.2.5 经营目标

1) 经济目标  经济目标为木材产量，即模拟期采伐木的蓄积量，由一元材积式计算获得。2)生态目标  ①树种和大小多样性用Shannon多样性指数表示(Buongiorno et al., 1995)。②地上碳贮量采用陈传国等(1989)的生物量模型计算地上部分生物量，碳贮量通过生物量乘以0. 5得到。需要说明的是，由于研究将树种划分为树种组，但各个树种有各自的生物量方程，因此树种组的生物量方程由某一树种的生物量方程替代，具体为:红云冷树种组由冷杉替代、慢阔树种组由色木替代、中阔树种组由白桦替代。

2.2.6 采伐方案评价

以上3个方面的目标既相互依赖又可能相互排斥，要求各子目标同时达到最优是困难的。为了对各种经营方案对各种经营目标的影响进行综合分析，需要构造一个总目标函数。由于各目标的量纲不一致，首先对其进行标准化处理。本文采用线性变换法即目标值与最大值之比。不失一般性，各目标权重都为1，目标函数为:

式中:T为标准化后的木材产量，ΔBsp，ΔBsz，ΔC分别为标准化后的树种多样性指数、大小多样性指数和地上碳贮量的变化量。

3 结果与分析 3.1 基本因子统计量

从表 3可以看出:落叶松的枯损概率最高，平均(每5年)为0. 127 6;慢阔的枯损概率最低，平均为0. 060 6;红云冷和中阔的枯损率分别为0. 072 4和0. 095 0，但中阔的变动较大。红松云冷杉有最大的进阶率，为0. 336 4;中阔的进阶率最低，为0. 295 0，但中阔变动较大; 落叶松和慢阔的进阶率分别为0. 319 7和0. 187 7。虽然落叶松的进界率最高，为0. 003 9，但标准差也最大，为0. 019 0，中阔的进界率最低，为0. 001 0。

3.2 子模型的参数估计

枯损、进阶和进界3个子模型最小二乘法的参数估计值如表 4，5，6所示。可以看出:影响枯损概率的变量包括径阶大小、林分断面积、树种多样性、最小径阶株数和海拔。径阶大小只对落叶松和阔叶树种组的枯损有显著影响，且随着胸径的增大树木枯损的比例会迅速降低，当胸径达到一定程度时树木的枯损下降程度开始减缓。针叶树种枯损率和胸高断面积呈正相关，尤其是落叶松，主要因为落叶松是阳性树种。对落叶松和慢阔树种，树种多样性与树木枯损呈现负相关; 同时随着海拔的增高，落叶松和红云冷的枯损率逐渐变大，而最小径阶株数只对落叶松枯损有显著影响，因为落叶松不能自然更新。

进阶模型中，落叶松径阶中值的倒数与进阶率呈现负相关，即随着其值的增加进阶率会增大。胸高断面积在针叶树种及中阔中系数为负，表明随着断面积的增大进阶率会有一定的下降。落叶松树种多样性与进阶率呈现正相关，说明混交对针叶树种有一定的促进作用。红松、云冷杉中大小多样性越大，该树种的进阶率越高。在针叶树种中，各树种(组)最小径阶的株数则与进阶呈正相关，其中慢阔的进阶没有自变量选入模型，因此采用该树种组建模数据中各径阶的进阶率平均值作为转移概率进行计算。

进界模型中，各树种(组)的胸高断面积系数为负，表明随着林分胸高断面积增加，进界率会下降。各树种(组)最小径阶株数的系数为正，其值较大时有较大的进界生长率。落叶松树种多样性越大其进界生长越小。海拔对落叶松和阔叶树种组的进界有显著影响，但系数很小。

不分树种模型表明:枯损模型中，在林分断面积较大或海拔较高时枯损率较大，树种多样性越大枯损率越小，较小径阶的树木有较大的枯损率; 在进阶模型中，林木较小或林分断面积越大时树木进阶越困难，而当树种多样性和最小径阶的株数变大时进阶率亦会越大; 林分的胸高断面积平方和树种多样性与进界量呈负相关，最小径阶的株数越大时进界生长越大。

采用似不相关回归估计3个子模型参数，发现3个子模型间误差相关很小(相关系数r＜0. 12)，可认为它们是独立的。因此，采用最小二乘法得到参数估计值。

3.3 矩阵模型的检验

基于1987年数据预测20块样地15年后及4块样地20年后(未参加建模)的径阶分布，各树种的预测值和观测值的卡方值及显著度为0. 05情况下的临界值如表 7所示。

从表 7可以看出:除慢阔树种外，其他树种(组)都通过卡方检验。主要原因是慢阔的进界模型中预测的进界值较高，且其第1径阶的进阶率偏低，使得第1径阶的树木在以后的分期内大多保留在第1径阶。但从总体来说，检验说明建立的矩阵生长模型可以用来预测落叶松云冷杉林的生长。

3.4 采伐方案模拟

利用建立的矩阵生长模型模拟不同采伐方案对木材生产、树种和大小多样性及地上碳贮量的影响。图 2~5给出了不同采伐方案下50年间木材生产、树种和大小多样性及地上碳贮量的变化。

3.4.1 木材生产

从不同采伐方案下的蓄积量图(图 2)可以看出:高采伐强度导致林分蓄积的持续下降且采伐周期越短下降越快(方案3，4，5)，50年时都恢复不到原来的水平(方案3，4，5，8，9);长周期和低采伐强度的采伐方案(方案6，7，10，11，12，13)则能在采伐后得到生长和恢复，林分蓄积到50年时达到或超过原来的水平。

不同采伐方案50年间的生长量无明显差异(表 8)，最小为方案2，即5年5%采伐强度的方案，生长量为363. 55 m³·hm^-2; 最大为方案5，即5年20%采伐强度的方案，生长量为394. 75 m³ ·hm^-2; 而方案4为次优，其生长量为392. 43 m³·hm^-2，年生长量在7. 27 ~ 7. 90 m³ · hm^-2 a^-1之间。但枯损量有着明显的区别，都随着采伐强度的增加而减小，以自然生长的枯损量为最大，以方案5为最小。方案5的总收获量最大，因此，单从收获木材的经营目标来看，方案5即短周期大强度的方案为最优方案。

3.4.2 树种和大小多样性

树种多样性表现出先增加后减少的趋势(图 3)，主要是因为采伐方案设计的是只采针叶树。整体来看(表 9)，高强度采伐方案3，8，9，13导致树种多样性的减少，而自然生长和低强度采伐(方案6，10)情况下树种多样性略有增加; 对于大小多样性，采伐方案3，4，5，8，9造成大小多样性的降低，而自然生长和长周期方案10使大小多样性增加较大。这些与采伐树种和采伐径阶的选择有关，因为本研究采用的采伐方式是采伐针叶树种组并且从最大径阶起伐，采伐周期变短、采伐强度(方案3，4，5，8，9)增大会导致大径阶的林木缺失。自然生长和采伐方案10在50年间多样性的增加最大。因此，从维持和增加多样性的角度，以自然生长和长周期小强度为最佳方案。

此外，该林分的起源是在伐后有部分保留树种下造的落叶松人工林，但已具有天然林的部分特征。从50年时不同采伐方案的各树种组成来看(表 9) :自然生长(方案1)落叶松比期初状态下降，而红松云冷杉、慢阔和中阔的比例略有增加; 采伐总体上导致落叶松所占比例的增加，红松云冷杉的比例既有增加也有减少，慢阔和中阔的比例略有增加。因此，该林分的自然演替是一个长期的过程。

3.4.3 地上碳贮量

从不同采伐方案碳贮量表中地上碳贮量的变化来看(图 5) :对短周期的采伐方案(方案2 ~ 5)，同一采伐周期内，碳贮量增加量基本随采伐强度的增加而减少，且不同采伐强度间差别较大，这是因为采伐除去了部分生物量，且在短期内难以通过生长进行恢复; 10年的采伐周期，弱度采伐的碳贮量增加量(方案6和7)高于强度采伐(方案9);对于15年的采伐周期，不同采伐强度间差别不大，且没有表现出一致性的规律。综合来看(表 8)，除短周期高采伐强度(方案5)外，其余方案都能保持地上碳贮量的增长。

3.4.4 综合目标

从以上可以看出: 3个目标间存在相互冲突，要满足多目标经营的要求，需要对各个目标进行折衷。将各目标标准化后加权(权重都相同)得到不同采伐方案的总目标，如表 10所示。可以看出: 13个方案中，在3个目标同等重要的前提下，以方案10(采伐周期15年采伐强度5%)为最优方案，方案1(自然生长)为次优，但显然方案1忽略了木材生产目标，属于保护型方案。因此，长周期低强度采伐可以同时满足人们对木材生产、保护多样性和增加碳贮量多目标的需要，使森林经营的总目标值达到最大。这说明合理的森林经营可以实现森林的多个目标。

4 结论与讨论 4.1 结论

1) 影响枯损、进阶和进界概率的主要变量包括径阶中值、公顷断面积、树种多样性、海拔和最小径阶株数，这和大多数研究一致(Lin et al., 1996; 1997; Ralston et al., 2003; Sist et al., 2003; Namaalwa et al., 2005; Liang et al., 2005; 2007; Lexerød，2005)。建立的矩阵生长模型可以用来预测落叶松云冷杉林的生长动态。

2) 建立的矩阵生长模型能灵敏地反映不同经营方案的差异。在综合考虑木材生产、树种和大小多样性和地上碳贮量多个目标的情况下，长周期低强度(15年5%采伐强度)为最优方案，表明合理的森林经营可以实现森林经营的多个目标。对于木材生产、树种和多样性和碳贮量，延长采伐周期、降低采伐强度可以同时增加木材产量、树种和大小多样性及地上碳贮量，这与Buongiorno等(1995)、Boscolo等(1997)、Favrichon (1998)的研究结论类似。该研究为模拟不同经营方案的多目标经营效果提供了一种方法。

3) 本研究中的落叶松云冷杉林是一种人工林与天然林混交林，落叶松作为阳性树种，不能天然更新，其最终会消失，但从模拟树种组成看，50年后落叶松仍占有一定的比例，因此该林分的发展演替是一个长期的过程，其发展演替有待进一步验证。

4.2 讨论

1) 关于转移概率子模型  转移概率子模型的模型形式，除了采用线性外，也有采用Logistic模型，它可以非常灵敏地反映因变量在0 ~ 1区间内的变化(Orois et al., 2002; Bollandsås et al., 2008)。对于影响转移概率的因子，本研究将树种和大小多样性作为自变量，发现其对进界、生长和枯损都有显著的影响，这对异龄混交林的生长模拟有一定的启示。在参数估计方法上，虽然研究中采用似不相关回归和普通最小二乘结果无显著差异，但前者仍是理想的方法。如Yang等(2008)的研究表明:采用似不相关回归有转移概率方程进行联合估计较单独估计单个方程有显著改善。本研究中的落叶松云冷杉林是一种介于人工林和天然林之间的近天然林，子模型的决定系数从0. 08到0. 66，但总体模型基本通过了卡方检验。模型的决定系数低并不意味着模型不能应用，只能说明选出的自变量只解释了部分转移概率的变异，但所有的自变量统计上均显著，且模型的常数项的标准差很小。大多数研究的转移概率子模型的决定系数从0. 004到0. 83，最小值从0. 004到0. 06，尤其是天然林(Orois et al., 2002; Ralston et al., 2003; Sist et al., 2003; Namaalwa et al., 2005; Liang et al., 2005; 2007)。若样本数量充足，可分径阶建立转移概率模型(Hao et al., 2005)。

2) 树种分组  将林分组成树种划分为树种组进行模拟，是矩阵生长模型的通用作法; 但这导致林分蓄积量、碳贮量等的计算存在偏差，也是矩阵模型的一个局限。因为以径阶为单元建模，要求每个组成树种都要有足够的样本量。

3) 采伐方案的设计  采伐设计是一个复杂的过程，需要考虑采伐周期、采伐强度及采伐木，对于混交林更是如此。本文仅考虑了采针按比例方案，这使采伐方案受到局限。若资料充足，以后的研究可以根据不同的需要继续丰富采伐方案。像针阔(组)都采、只采阔叶树种(组)、保留珍贵树种(组)采伐、确定不同的最小采伐直径等多种方式进行模拟比较，为经营者实施更灵活的经营策略提供理论依据。

4) 多目标经营  本研究对多目标经营模拟进行了尝试，除木材生产外，还考虑了生物多样性和碳贮量目标。但碳贮量仅考虑了地上碳贮量，未包括根、土壤及地表枯落物中的碳。最新的研究已经可以跟踪产品碳(Baskent et al., 2008)，因为林木被采伐后碳并未完全释放，而是以不同的林产品形式固定下来，对不同的林木产品有不同碳吸存值，这就需要对林木的出材量及木材产品对木材的利用率进行分析。此外，模拟的经营方案有限，下一步可和优化模型相结合，得到最优解指导森林经营规划。