在森林经营中，蓄积量预测是一项重要而复杂的工作，一直是测树学、森林资源清查等工作的重点之一，也是林业生产单位最为关注的指标(吴中伦，1984；孟宪宇，1996)。杉木(Cunninghamia lanceolata)作为我国南方的主要用材树种之一，在满足国民经济发展和人民群众生活对森林多种效益的需求上具有重要的地位和作用。因此，有必要对杉木人工林的蓄积收获进行研究，以便及时、准确地获得杉木资源的数据，指导林业生产和森林经营活动，并做出科学合理的决策(杜纪山等，2000)。

在构建林分生长和收获模型时，经常要用到联立方程组模型(Borders et al., 1986；Hasenauer et al., 1998)。目前联立方程组模型已经在林业上有所应用(Borders，1989; Fang et al., 2001a; 郎奎建，2004)，这有2方面的原因：一是基于相容性需要，即生长模型中的输出变量要与收获模型中的相一致；二是在估计模型参数时，一些变量在某一方程中是因变量而在另一个方程中就成了自变量，因此多个模型相互之间共有的一些变量存在着连续相关的问题(Fang et al., 2001a)。为了避免在估计中产生联合偏差，很多学者在进行参数估计时采用了很多方法，比如三阶段最小二乘方法、矩法和信息似然方法等(Borders，1989；LeMay，1990)。这些方法都没有完全很好地解决模型偏差问题，在模型估计精度上也没有大的提高。混合效应模型已经被证明在提高模型估计精度上有优势，因此可探索在联立方程估计中考虑参数的随机效应(Lynch et al., 2005；Uzoh et al., 2006；Sharma et al., 2007)。Hall等(2004)利用多元多层次的非线性混合效应模型方法(NLME)描述了美国乔治亚州和佛罗里达州湿地松(Pinus elliottii)的木材收获，作者首先利用一个三元的NLME方法来模拟优势木高、公顷断面积和公顷株数，然后利用这几个变量的预测值输入到蓄积模型中，最后对林分未来的生长和收获进行预测。Zhao等(2005)利用多层次的NLME方法模拟了经过4种不同经营措施处理的火炬松(Pinus taeda)林分的蓄积生长，认为控制采伐和施肥措施导致了最大的蓄积生长量，NLME方法能够很好地解决不均衡和不完整的重复测量数据问题。Fang等(2001a)基于混合效应模型方法构建了美国乔治亚州和佛罗里达州湿地松林分优势木平均高、林分断面积和林分蓄积作为内生变量的联立方程组，然后考虑了样地层次的随机效应及重复观测相关性，预测了经过不同经营措施处理的蓄积生长趋势，结果认为与传统方法相比，蓄积的误差从52.6%减小到5.8%。

上述学者只有Fang等(2001a)的蓄积估计是基于联立方程组方法，并且没有很好地解决内生变量间误差在多个模型间的传递问题。基于此，本文利用江西杉木样地数据，探讨基于样地层次的混合效应，采用模拟数据构建林分优势木平均高、林分断面积及林分蓄积的联立方程组模型。最后选择验证数据对传统模型和混合效应模型方法进行验证和预测。

1 数据来源

数据来源于江西省568块人工杉木固定样地数据，面积为0.08 hm²，从1991年开始，每隔5年调查1次，共进行了4次调查。由于部分样地在不同调查期间进行了各种方式和各种强度的采伐，对林分断面积和蓄积，尤其是优势木平均高的模拟有很大的影响，因此去掉了上层疏伐以及择伐措施的样地，只保留下层疏伐及采伐强度在20%以下的样地，最后选择了符合条件的365块样地(1 210个样本)。随机地把数据分成模拟数据和验证数据，其中模拟数据包括851个样本，验证数据包括359个样本。样地基本调查因子包括海拔、坡向、坡位、坡度、地貌特征、平均年龄、平均胸径、平均树高、样木总株数等。此外计算和统计了每个样地的公顷株数、公顷断面积、公顷蓄积和优势木平均高，样地因子统计结果见表 1。

2 基于非线性混合效应的蓄积联立系统 2.1 联立方程组模型的确定

在构建蓄积联立方程系统中，选择了3个内生变量，即优势木平均高、林分断面积和林分蓄积。优势木平均高形式见Fang等(2001b)，林分断面积形式见李春明等(2004)，林分蓄积形式见Fang等(2001a)。为了确定这3个模型在模拟时是否有意义，利用SAS软件采用最小二乘回归参数估计方法分别进行拟合。由于立地条件对杉木生长有很大的影响，因此在模拟时考虑了立地因子，以哑变量的形式加入到模型中去。经过多次模拟，选择模拟精度最高并且参数的回归系数显著，即P < 0.05的模型形式，结果发现只有地貌特征对林分断面积生长有显著影响(方差分析差异显著)，其他立地因子影响不显著，模拟结果见表 2。

利用F检验方法对模型进行检验。结果表明:林分优势木平均高模型的F=26 395.5>F_0.95临界值(P < 0.000 1)，林分断面积模型的F=1 063.68>F_0.95临界值(P < 0.000 1)，林分蓄积的F=10 606.6>F_0.95临界值(P < 0.000 1)，这说明3个模型对于各自描述的生长趋势有显著意义。通过t检验来判断模型参数是否显著，结果所有参数P值均小于0.05，说明这些参数描述模型差异显著。因此最后确定的联立方程组形式如下：

(1)

式中：HD为优势木平均高；BA为林分断面积；V为林分蓄积；N为单位面积株数；dm为地貌特征的量化指标；t为林龄；ε为误差；其余为待估参数。

2.2 混合模型模拟方法 2.2.1 基于混合模型的联立方程求解方法

2.1部分已经确定了联立方程组形式，因此如何考虑混合参数是本文首先要解决的问题。联立方程的基本原理是：

式中：；i为样地数；t_k为测量次数。

在利用混合模型方法模拟时，这3个向量(林分断面积、林分优势木平均高和林分蓄积)可合并成一个长向量，在长向量中，这3个变量可通过哑变量或不同的函数指数来区别。本文中，Level 1表示林分优势木平均高、Level 2表示林分断面积、Level 3表示林分蓄积，不同的Level有不同的数学期望及方差。合并成一个长向量后，就可以按普通混合模型方法求解固定和混合参数，计算主要通过S-PLUS软件中的NLME模块进行。

2.2.2 混合参数的确定

基于样地水平，分别对林分优势木平均高、林分断面积和林分蓄积单独进行混合效应模型的拟合，基于模型的效果评价指标(AIC, BIC及-2lg Likelihood)，最后确定β₁和β₂作为林分优势木平均高公式中的混合参数、b₃和b₅作为林分断面积模型中的混合参数、c₁和c₂作为林分蓄积模型中的混合参数时模拟精度最高。在(1)式中，混合参数表示为：

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

对式(1)的联立方程组进行拟合，结果β₁和β₂同时为混合参数时，模型不能收敛，c₁和c₂同时为混合参数时，模型也不能收敛。林分优势木平均高模型中考虑β₁为混合参数要比β₂为混合参数的模拟效果好，林分蓄积模型中c₁作为混合参数要比c₂为混合参数的模拟效果好。因此最后确定β₁作为林分优势木平均高模型中的混合参数、b₃和b₅作为林分断面积模型中的混合参数、c₁作为林分蓄积模型中的混合参数同时进行模拟。在模拟过程中，这4个参数同时进行模拟时，模型不能收敛，因此随机去掉1个参数，选择精度最高的作为最后结果，模拟结果显示β₁，b₃和b₅同时作为混合参数时模拟精度最高, 然后再随机去掉1个参数，发现模拟效果明显降低，且差异显著(P < 0.000 1)，因此最后确定的混合参数为β₁，b₃和b₅，通过S-PLUS软件的NLME模块进行模拟，并与传统回归方法进行比较，其结果见表 3。

从表 3可看出，基于混合效应模型的联立方程组比传统回归方法3个指标值要小得多，说明考虑参数的随机效应后，模拟效果要好得多，利用LRT(likelihood ratio test, 似然比检验)进行方差分析，LRT=6 284.162(P < 0.000 1)，说明差异显著。

2.2.3 考虑变量的相关性及异方差问题

由于优势木平均高与断面积和蓄积、断面积与蓄积之间存在着一定的相关性，在模拟时会产生联合偏差，因此要考虑三者误差之间的序列相关性，这样就可以有效地降低产生的联合偏差。为了表达3个变量之间的相关性，选择了一阶自回归矩阵模型[AR(1)]、一阶自回归与滑动平均模型相结合的矩阵模型[ARMA(1, 1)]及复合对称矩阵模型(CS)作为相关性矩阵结构来描述(Wang et al., 2007)。最后拟合的比较结果见表 4。

从表 4可看出，同不考虑误差结构矩阵相比，3个误差结构矩阵的模拟效果都有很大的提高，3个指标值均有下降，LRT分析表明差异显著。CS结构在3个指标值上比其他2个更加合理，因此可作为优势木平均高、断面积和蓄积三者之间的误差结构矩阵。

模型是否存在异方差问题可通过残差分布图来分析。图 1是传统回归方法和基于混合效应模型方法的残差分布图(图 1A1, B1, C1为传统回归模型方法，图 1A2, B2, C2为基于混合效应模型方法)。

图 1表明，与传统回归方法相比，基于混合效应模型方法残差不仅分布范围大幅减小，而且分布也更均匀。传统回归方法存在着明显的异方差，而基于混合效应模型方法异方差并不明显。利用幂函数和指数函数来描述基于混合效应模型方法异方差，结果模型不收敛。因此，最终综合考虑样地效应和误差结构矩阵的模型基本形式如式(2)。

(2)

式中：D为样地间随机效应方差协方差矩阵；R_i为样地内误差效应方差协方差矩阵；u_i为随机效应参数；Ψ_i^0.5为异方差矩阵；Γ_i(θ)为自相关矩阵；i为样地数。其余变量和参数意义同式(1)。

对式(2)进行模拟，具体结果见表 3的c部分。在实际预测过程中，如果已知前几期林分优势木平均高、林分断面积和林分蓄积的测量结果，则可采用表 3的模拟结果对林分未来优势木平均高、林分断面积和林分蓄积进行估计和预测。

3 模型验证

对基于传统回归估计方法和混合效应模型方法的联立方程组模拟完成后，利用验证数据对2种方法的模拟精度进行验证，采用的评价指标包括确定系数(R²)、平均绝对残差(|E|)及均方根误差(RMSE)3个模型精度评价指标对模拟结果进行效果评价，|E|及RMSE越小而R²值越大说明模型的精度越高(Calama et al., 2004)。具体的验证结果见表 5。

从表 5可看出，无论是林分优势木平均高、林分断面积还是林分蓄积，考虑样地的混合效应及误差方差矩阵的3个评价指标值都要优于传统回归估计方法，说明混合效应模型的精度较传统的回归估计方法好。这一结论与前面模拟效果结论是一致的。

4 模型预测

建立联立方程组的主要目的就是对未来林分的生长进行预测，在预测过程中至少要考虑2种形式，即不考虑随机效应和考虑随机效应。

对于这2种形式，在联立方程组模型预测过程中，每一种形式又可再分为2种情况：一种是方程右侧的所有内生变量和外生变量值在当前预测期被观测了；另一种情况是方程右侧仅部分变量在当前预测期被观测了。对于形式2来说，还有2个与形式1不同的特殊情况：一种是过去所有观测时期的观测记录是完整的，即所有变量包括外生变量在过去被完整地测量和记录了；另一种情况是对于过去的所有观测，某些时期的观测记录完整，另一些时期的观测数据不完整，只是部分数据被记录了。这是林业调查过程中最普遍的事情。例如，在林业调查或经营中，对某一林分未来蓄积和断面积进行预测时，一种情况是，过去某一时期所调查的记录很全面，包括蓄积、断面积以及优势木平均高、林分平均胸径以及林分年龄等因子；另一情况是，过去某一时期仅仅调查了优势木平均高、林分平均直径和林龄，或者是仅仅调查了林龄。这些情况就会造成预测结果的不同。一般来讲，对于某一样地，在联立方程组模型中，假如过去的观测值已知，则这个样地的随机参数可估，这样就可被用来提高估计精度。

正常情况下，林分的生长过程主要包括3种情况，即自然稀疏生长过程、先等株数生长后自然稀疏生长过程以及等株数生长过程。由于大部分杉木样地都经过了1次或多次采伐，很难确定未来林分具体株数，也很难建立林分株数与其他林分因子的关系模型，因此，本文针对数据不同信息情况按等株数生长过程进行蓄积预测。

1) 假如知道林分当前的年龄是13年，并且在年龄为13年的林分密度是2 463株·hm^-2，而没有其他测量信息。这种情况就无法计算此样地的随机效应参数，只能够考虑参数的固定效应，按以下方法预测：

估计出优势木平均高的值后就可以把估计值代入断面积公式中进行断面积估计，估计值为：

把断面积和优势木平均高估计值代入蓄积公式，最后求出的蓄积值为：

最后计算出的蓄积为49 m³。同样，根据这些方法可以预测出18，23，28年及以后各林龄的估计值。

2) 假如只知道林龄13年时的优势木平均高是8.1 m，在估计13年时的林分断面积则可把优势木平均高值直接代入到断面积公式中，则估计的断面积值为11.33 m²。把优势木平均高值和估计的断面积值代入蓄积公式中，最后估计的蓄积值为38.84 m³。

3) 假如只知道林龄13年时的林分断面积(其值为10.46 m²)，其他因子未知，则需要预测林分的优势木平均高和蓄积。

根据公式计算出的优势木平均高为10.8 m，把估计的优势木平均高值和实测的林分断面积代入到蓄积公式中，最后计算出的蓄积为39.48 m³。

4) 假如测量了林龄13年时的优势木平均高为8.1 m，林分单位断面积为10.46 m²，其他因子同上，则可以直接把优势木平均高值和林分单位断面积值代入到蓄积公式中，最后计算出的蓄积为35.58 m³。

以上的计算都是基于没有考虑随机效应的情况，因此只能预测当前林龄的蓄积，本样地实测蓄积为31.33 m³。从上述结果可看出只知道林分年龄和株数情况的误差最大，达到了17.67 m³；其次是知道林分年龄、株数和断面积，误差达到了8.15 m³；再次是知道林分年龄、株数和优势木平均高，误差为7.51 m³；误差最小的为知道林分年龄、株数、优势木平均高和断面积，只有4.25 m³。从而可以说明，获取样地的信息越多，预测的精度就会越高。下面是考虑样地随机效应的预测。

5) 如果测量了13，18和23年时的优势木平均高，则可根据连续重复观测通过随机效应参数公式计算出样地在估计优势木平均高公式中的随机效应，进而调整计算公式，然后计算出其他林龄时的优势木平均高，最后再计算蓄积。例如要估计林龄为28年时的蓄积则首先要确定优势木平均高公式的样地效应参数值，然后再估计28年时的优势木平均高，把优势木平均高值代入到断面积公式中，最后把优势木平均高值和林分断面积值代入到蓄积公式中估计出蓄积。具体计算如下。

首先利用固定效应参数估计出优势木平均高的预测值，进而用测量值减去预测值计算出ê_i，Z_i^T是设计矩阵(具体是优势木平均高公式的偏导数)，和在模拟中给出，最后算出b_i。如果不考虑优势木平均高、断面积和蓄积之间误差存在的相关性，则可考虑为单位矩阵。

则蓄积为83.5 m³。

6) 如果只测量了13，18和23年时的林分断面积，则可根据连续重复观测通过随机效应参数公式计算出样地在估计林分断面积公式中的随机效应，进而调整计算公式，然后计算出其他林龄时的林分断面积，最后再计算蓄积。例如要估计林龄为28年时的蓄积则首先要确定林分断面积公式中的样地效应参数值，然后再估计28年时的林分断面积，最后把估计的优势木平均高值和林分断面积值代入蓄积公式中估计出蓄积。具体计算如下：

利用固定效应参数估计出林分断面积的预测值，进而用测量值减去预测值计算出ê_i，Z_i^T是设计矩阵(具体是林分断面积公式的偏导数)，和在模拟中给出，最后算出b_i。

则BA=exp(b₀+b₀₁·dm+b₁/t)N^b₂+b₄/tHD^{(b₃+0.242 7)+b₃₁·dm+(b₅-3.806 9)/t}=26.87 m²。

则最后蓄积为：

则蓄积为121.1 m³。

7) 如果测量了13，18和23年时的优势木平均高和林分断面积，则可估计出优势木平均高的样地效应参数(即β₁的随机效应值)和林分断面积的样地效应参数(即b₃和b₅的随机效应值)。最后再估计林龄28年的蓄积。计算如下：

最后的蓄积为：

则蓄积为110.1 m³。

8) 如果测量了13，18和23年时的优势木平均高、林分断面积和蓄积，则可估计出优势木平均高的样地效应参数(即β₁的随机效应值)、林分断面积的样地效应参数(b₃和b₅的随机效应值)。与前7个预测情况不同的是这种情况可获得3个变量之间的误差自相关矩阵。最后再估计林龄28年的蓄积。根据表 5的模拟结果计算出的蓄积为107.3 m³。

5 结论与讨论

1) 与传统回归估计方法相比，混合效应模型方法在模拟和估计精度上有很大的提高，并且能够很好地解释和减小误差传递问题。

2) 优势木平均高是联立方程组中最基本的组成部分，因为优势木平均高是断面积模型和蓄积模型中的自变量。另外，优势高的预测误差是断面积和蓄积模型中的主要误差来源，因此优势木平均高的准确预测是联立方程组的关键之一。断面积的准确预测也是联立方程组模型中十分重要的内容，断面积的准确与否与预测蓄积的准确性有直接的关系。

3) 综合模拟结果表明，在优势木平均高模型中β₁作为混合参数，b₃和b₅作为断面积模型的混合参数的联立方程组模拟效果最好。确定好合理的混合参数及联立误差自相关矩阵后，解释蓄积的随机效应被证明是没有必要的。

4) 优势木平均高、林分断面积既作为自变量又作为因变量，因此考虑三者之间的误差结构能够提高估计精度，并且CS结构描述三者的误差自相关性效果最好。

5) 在联立方程组中，考虑异方差时模型没有收敛，但这并不意味着异方差不重要。实际上，参数的随机效应和3个变量(优势木平均高、林分断面积和蓄积)的相关性能够解决或减少异方差问题。

6) 本研究在预测未来林分生长时，限于数据本身问题，只假定林分株数为等株数生长，没有考虑其他情况，如以后有更充分的数据则可相应地构建林分株数与其他林分因子关系模型来进行预测。

6) 本研究在预测未来林分生长时，限于数据本身问题，只假定林分株数为等株数生长，没有考虑其他情况，如以后有更充分的数据则可相应地构建林分株数与其他林分因子关系模型来进行预测。