木材密度(wood density 或 specific gravity)是指单位体积的木材质量，是木材性质的一项重要指标。木材密度大小直接影响木材的物理和力学性质(MacDonald et al.，2002)。木材密度与木材的许多性质都存在着一定相关性，如木材密度与微纤丝角间存在一定的相互关系，并与木材强度和硬度密切相关(Mäkinen et al.，2002 ; Lundgren，2004 ;Jaakkola et al.，2005 ; Jordan et al.，2005)。过去人们对木材密度的研究主要是从定性分析的角度来研究木材密度的变异规律(骆秀琴等，1999 ; 马丽娜等，2003 ; 罗建勋等，2004 ; 刘青华等，2010 ; 王秀花等，2011)。木材密度在径向变异模式上大致划分为3种类型: 1)从髓心至树皮逐渐增加(直线增加或曲线增加); 2)从髓心向外逐渐降低，然后又升高(靠近树皮的木材密度比髓心的密度高或者低); 3)从髓心至树皮逐渐减小(直线或曲线减小)。木材密度在纵向变异模式上也大致划分为3种类型: 1)沿树干方向从下向上木材密度逐渐变小; 2)从下向上木材密度逐渐变大; 3)从下向上木材密度逐渐变小，而到树冠基部又逐渐变大。

过去木材科学工作者主要是利用线性和非线性回归来构建木材密度模型(姜景民等，1999; 徐有明等，2002; Pérez Cordero et al.，2002; Espinoza，2004)。由于木材密度有较大变异的特点，近年来混合模型方法在国外已用于构建木材密度模型(Molteberg et al.，2007; Mäkinen et al.，2007; Jordan et al.，2008; Schneider et al.，2008; Fujimoto et al.，2010)。本文以兴安落叶松(Larix gmelinii)人工林为例，考虑2层次效应，即单木效应和树干高度效应，构建落叶松木材密度2层次混合效应模型，并与单层混合效应模型及传统模型拟合效果进行对比分析。

落叶松试样采自七台河市林业局金沙林场3块立地条件相似的样地。在每块样地内伐取优势木、平均木和被压木各1株，共9株。标准木伐倒后从树干的根部、胸高处(1.3 m)、树高的20 % 、40 % 、60 % 、80 % 处截取5 cm 厚的圆盘各1个。在圆盘年轮比较清晰的部位自髓心向外取30 °楔形木块，在每个楔形木块角度平分处画一条直线。根据这条直线把楔形木块划分成相等的8份，并测量每部分的宽度和年轮数，然后将每部分切割下来。为了尽量保证试验精度，切割主要沿年轮方向切割。由于样品形状不规则，因此用排水法得到各样品的体积，即将样品浸泡在水中，根据测得的浮力来计算样品体积，然后用圆盘绝干质量[(103 ± 2)℃]与试样饱和水分时的体积之比得到木材基本密度。每株树共截6个圆盘，每个圆盘沿年轮方向共收集8个样本，因此每株树共测量48个样本的基本密度。把木材密度数据分成建模数据和检验数据，各样木生长轮年龄和木材密度统计见表 1 。

表 1 生长轮年龄和木材密度统计 Tab.1 Summary statistics for ring age and wood density

表 1 生长轮年龄和木材密度统计 Tab.1 Summary statistics for ring age and wood density 数据 Data 树号 Tree No． 样本数 Number 生长轮年龄 Ring age / a 木材密度 Wood density /(g·cm -3) 平均值 Mean 最小值 Min． 最大值 Max． 标准差 Std． 平均值 Mean 最小值 Min． 最大值 Max． 标准差 Std． 建模数据 2 48 15.62 3 41 9.54 0.52 0.31 0.71 0.09 Fitting data 4 48 14.87 3 39 9.48 0.54 0.21 0.70 0.09 5 48 14.60 3 41 9.15 0.49 0.36 0.83 0.08 7 48 13.87 1 43 9.76 0.52 0.38 0.77 0.08 8 48 14.77 3 39 8.81 0.50 0.39 0.67 0.07 9 48 14.71 2 43 9.93 0.57 0.37 0.99 0.11 检验数据 1 48 14.06 1 43 9.53 0.51 0.32 0.74 0.08 Validation data 3 48 14.75 3 39 9.74 0.50 0.35 0.71 0.07 6 48 14.32 3 41 9.23 0.50 0.34 0.77 0.07

线性混合模型包括单水平和多水平形式，根据 Pinheiro等(2000)，单水平线性混合模型的形式为:

多水平线性混合模型(multilevel linear mixedeffects model，MLMEM)是基于单水平线性混合模型基础上增加随机效应因子构成的。由多个相互嵌套随机效应因子产生的 MLMEM 称为多水平MLMEM。对于嵌套2水平的 MLMEM，模型表达式为:

根据 Pinheiro等(1998)，构建混合模型需要确树木内方差协方差结构(R_i)。采用式(3)来描述(Davidian et al.，1995 ; Calama et al.，2004 ; Trincado et al.，2006):

3)确定随机效应的方差协方差结构(D)。比较林业上常用的3种方差协方差结构:复合对称(CS)、对角矩阵(diagonal matrix)、正定矩阵(generalpositive-definite matrix)(Fang et al.，2001; Jordan et al.，2005)，以包括3个随机参数(b₁，b₂，b₃)的方差协方差结构为例，

复合对称:

对角矩阵:

广义正定矩阵:

拟合和检验结果通过以下指标评价:绝对误差(Bias)、均方根误差(RMSE)、确定系数(R²)。

$ {\rm{Bias = }}\frac{{\left| {\sum\nolimits_{i = 1}^m {\sum\nolimits_{j = 1}^{{m_i}} {\sum\nolimits_{k = 1}^{{m_{ij}}} {\left({{y_{ijk}} -{{\hat y}_{ijk}}} \right)} } } } \right|}}{n}; $

$ {\rm{RMSE = }}\sqrt {\frac{{{{\sum\nolimits_{i = 1}^m {\sum\nolimits_{j = 1}^{{m_i}} {\sum\nolimits_{k = 1}^{{m_{ij}}} {\left({{y_{ijk}} -{{\hat y}_{ijk}}} \right)} } } }^2}}}{{n -1}}} ; $

$ {R^2}{\rm{ = 1 -[}}\frac{{{{\sum\nolimits_{i = 1}^m {\sum\nolimits_{j = 1}^{{m_i}} {\sum\nolimits_{k = 1}^{{m_{ij}}} {\left({{y_{ijk}} -{{\hat y}_{ijk}}} \right)} } } }^2}}}{{{{\sum\nolimits_{i = 1}^m {\sum\nolimits_{j = 1}^{{m_i}} {\sum\nolimits_{k = 1}^{{m_{ij}}} {\left({{y_{ijk}} -{{\bar y}_{ijk}}} \right)} } } }^2}}}{\rm{]}}。 $

模型检验使用 Vonesh等(1997)的方法来计算随机参数b_k:

初步分析木材密度与年龄、不同树干高度、胸径、树高、冠幅等林木变量的关系，研究发现木材密度与年龄和不同树干高度有密切关系，并随着年龄的增加而增加，达到最大值之后有所下降。采用逐步回归技术建立如下木材密度模型:

利用最小二乘法对模型(8)进行拟合，利用t检验和F检验法分别对模型参数和模型进行检验。参数t检验证明模型(8)所有参数估计值都是显著的(P < 0.000 1)，说明这些参数显著影响模型的变化。模型F检验(F=87.54，P < 0.000 1)都达到了极显著水平。

将木材密度模型(8)写成2水平混合效应模型形式为:

单水平是在2水平混合模型基础上去掉随机参数b_i，j。混合模型在拟合时，需要确定分类变量。根据木材密度数据的特征，分类变量包括单木、树干不同高度及二者同时作用。因此本研究分别考虑单水平和2水平，单水平主要包括单木或树干不同高度效应，2水平指将二者同时考虑。当考虑单木效应影响时，单木水平中每组对应的观察次数都是48 ;当考虑树干不同高度效应时，每个树干高度水平对应的观察次数都是8 ;当同时考虑2水平时，每个分组对应的观察次数也是8。

拟合结果见表 2。可以看出，能够收敛的混合模型有6种。当模型没有随机参数时，AIC 值等于-859.781 2，BIC值等于-840.695 6。当考虑1个随机参数时，模拟1的 AIC和BIC 值最小; 当考虑2个随机参数时，模拟5的 AIC和BIC 值最小;当考虑3个随机参数时，模型不能收敛。LRT 检验表明:基本模型和模拟1有显著不同(P < 0.000 1)，模拟1和模拟5没有显著不同(P=0.389 7)，这说明增加随机参数个数不能提高模型的拟合精度。综合以上分析，把含有随机参数 b1的模型选为基于单木效应的最优混合模型。

表 2 基于单木效应混合模型拟合结果^① Tab.2 Fitting results of mixed model based on individual tree effects

表 2 基于单木效应混合模型拟合结果① Tab.2 Fitting results of mixed model based on individual tree effects 模拟 Simulation 随机参数Random parameter 参数个数 Number of parameters AIC BIC Log-likelihood LRT P ①参数个数包括固定参数个数、随机参数方差协方差参数个数及模型偏差参数。下同。 The number of parameters includes a fixed number ofparameters，random pammeer variance and covariance parameters and model error parameters． The same below． 基本模型 Basic model None 5 -859.781 2 -840.695 6 434.890 6 1 b1 6 -888.953 2 -866.050 6 450.476 6 31.172 0 < 0.000 1 2 b2 6 -883.660 6 -860.757 9 447.830 3 3 b3 6 -874.111 6 -851.208 9 443.055 8 4 b1，b2 6 -886.334 4 -855.797 5 451.167 2 5 b1，b3 8 -886.838 1 -856.301 2 451.419 1 1.884 8 0.389 7 6 b2，b3 8 -881.846 7 -851.309 8 448.923 4 7 b1，b2，b3 11 不收敛 No convergence

拟合结果见表 3 。可以看出，能够收敛的混合模型有7种。当考虑1个随机参数时，模拟1的 AIC和BIC 值最小; 当考虑2个随机参数时，模拟4的 AIC和BIC 值最小; 当考虑3个随机参数时，模拟7优于模拟1和4 。LRT 检验表明:基本模型和模拟1有显著不同(P < 0.000 1)，模拟1和模拟4有显著不同(P < 0.000 1)，模拟4和模拟7有显著不同(P < 0.000 1)。综合考虑不同随机参数的组合，把模拟7，即把基本模型含有随机参数 b₁，b₂，b₄的模型选为基于树干不同高度效应的最优混合模型。然后又测试了林业上常用的3种方差协方差结构。通过表 4的 AIC，BIC，Log-likelihood和似然比检验等统计量指标表明:正定矩阵结构优于对角矩阵和复合对称结构。

表 3 基于高度效应混合模型拟合结果 Tab.3 Fitting results of mixed model based on height effects

表 3 基于高度效应混合模型拟合结果 Tab.3 Fitting results of mixed model based on height effects 模拟 Simulation 随机参数Random parameter 参数个数 Number of parameters AIC BIC Log-likelihood LRT P 基本模型 Basic model None 5 -859.781 2 -840.695 6 434.890 6 1 b1 6 -888.545 9 -865.643 2 450.272 9 30.7647 < 0.000 1 2 b2 6 -884.167 6 -861.264 9 448.083 8 3 b3 6 -869.446 1 -846.543 5 440.723 1 4 b1，b2 6 -899.044 2 -868.507 3 457.522 1 14.498 3 0.000 7 5 b1，b3 8 -885.840 9 -855.304 0 450.920 5 6 b2，b3 8 -895.268 7 -864.731 8 455.634 4 7 b1，b2，b3 11 -915.689 2 -873.700 9 468.844 6 22.644 9 < 0.000 1

表 4 基于不同方差协方差结构混合模型拟合结果 Tab.4 Fitting results of mixed model based on different variance-covariance structures

表 4 基于不同方差协方差结构混合模型拟合结果 Tab.4 Fitting results of mixed model based on different variance-covariance structures 方差协方差结构 Variance-covariance structure 参数个数Number of parameters AIC BIC Log-likelihood LRT P 正定矩阵 General positive-definite matrix 11 -915.689 2 -873.700 9 468.844 6 对角矩阵 Diagonal matrix 8 -888.857 6 -858.320 8 452.428 8 32.831 5 < 0.000 1 复合对称 CS 7 -867.864 9 -841.145 1 440.932 4 55.824 3 < 0.000 1

综合考虑单木效应和树干高度效应，这属于多层次混合模型问题，与单层混合模型的研究方法相同，拟合结果见表 5 。可以看出，能够收敛的混合模型共有6种。当同时考虑6个随机参数时，模型不收敛;当考虑5个随机参数组合时，模拟2的 AIC和BIC 值最小;当考虑4个随机参数组合时，模拟5的 AIC和BIC 值最小。综合考虑不同随机参数的组合，模拟5的效果最好。 LRT 检验表明:基本模型和模拟2有显著不同(P < 0.000 1)，模拟2和模拟5没有显著不同(P=0.61)。综合以上分析，把模拟5，即把样木效应含有随机参数b_1i和树干不同高度效应含有随机参数b_1i，j，b_2i，j ，b_3i，j 的模型选为基于2层次混合模拟的最优混合模型。

表 5 基于单木效应和高度效应混合模型拟合结果^① Tab.5 Fitting results of mixed model based on individual tree and height effects

表 5 基于单木效应和高度效应混合模型拟合结果① Tab.5 Fitting results of mixed model based on individual tree and height effects 模拟 Simulation 单木效应 Individual tree effects 高度效应 Height effect AIC BIC Log-likelihood b 1 i b2 i b 3i b 1 i，j b 2i，j b 3i，j ①带△变量为包括这个随机参数。 With △ variables including the random parameters． 1 △ △ △ △ △ △ 不收敛 No convergence 2 △ △ △ △ △ -865.593 6 -812.321 7 446.796 8 3 △ △ △ △ △ -865.384 8 -812.112 9 446.692 4 4 △ △ △ △ △ -865.515 3 -812.243 4 446.757 6 5 △ △ △ △ -868.604 9 -822.943 2 446.302 4 6 △ △ △ △ -867.631 0 -821.969 3 445.815 5 7 △ △ △ △ -862.443 9 -816.782 3 443.222 0

误差的异质性主要通过残差分布图进行比较。可以看出，与基本模型的残差分布图相比(图 1A)，基于单木效应混合模型残差分布图变化不明显(图 1B)，这也说明，单木效应对木材密度影响不大。而基于树干高度效应和基于2层次混合效应模型的残差分布图(图 1C和D)显示了更均匀的分布。误差相关性测试了林业上常用的4种描述时间序列的自相关结构[CS，AR(1)，MA(1)和ARMA(1，1)]，结果显示没有提高混合模型的精度，因此时间序列自相关结构在本研究中没有进一步考虑。

	图 1 基于传统模型(A)、单木效应混合模型(B)、高度效应混合模型(C)和2层次效应混合模型(D)的残差分布 Fig.1 Residual plots of traditional model(A)，mixed model based on individual tree effect(B)，mixed model based on height effect(C)and two-level effects mixed model(D)

表 6给出了落叶松不同木材密度模型的固定参数估计值、随机参数方差协方差的参数估计值。

表 6 不同模型参数和方差组成估计值^① Tab.6 Parameter and variance components estimates of different models

表 6 不同模型参数和方差组成估计值① Tab.6 Parameter and variance components estimates of different models

本研究以1.3m处木材密度预测的一个例子来说明:从检验数据中树高1.3m处随机抽取年龄和密度的6次测量值(5，0.35)，(9，0.50)，(11，0.53)，(15，0.59)，(19，0.63)，(27，0.64)来计算。由于树干高度效应和2层次效应对混合模型拟合结果影响相似，因此参数和随机效应方差协方差估计值采用基于高度效应的参数估计值。首先计算${_{{{\hat e}_k}}}$，${{{\hat R}_i}}$${{{\hat Z}_k}}$，然后利用式(7)计算随机参数${_{{{\hat b}_k}}}$:

$ {{\hat e}_k} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {0.35}\\ {0.50}\\ {0.53} \end{array}}\\ {0.59}\\ {0.63}\\ {0.64} \end{array}} \right] -\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {0.48725}\\ {0.52765}\\ {0.54485} \end{array}}\\ {0.57325}\\ {0.59365}\\ {0.61045} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} { -0.13725}\\ { -0.02765}\\ { -0.01485} \end{array}}\\ {0.01675}\\ {0.03635}\\ {0.02955} \end{array}} \right], $

$ {{\hat R}_i} = {\sigma ^2}I(6)= 0.0025\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {1\;\;0\;\;0\;\;0\;\;0\;\;0}\\ {0\;\;1\;\;0\;\;0\;\;0\;\;0} \end{array}}\\ {0\;\;0\;\;1\;\;0\;\;0\;\;0}\\ {0\;\;0\;\;0\;\;1\;\;0\;\;0}\\ {\begin{array}{*{20}{c}} {0\;\;0\;\;0\;\;0\;\;1\;\;0}\\ {0\;\;0\;\;0\;\;0\;\;0\;\;1} \end{array}} \end{array}} \right], $

$ {{\hat Z}_k} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\frac{{\partial f\left({R{N_{ij1}},{\beta _1}} \right)}}{{\partial {\beta _1}}}\;\;\;\frac{{\partial f\left({R{N_{ij1}},{\beta _2}} \right)}}{{\partial {\beta _2}}}\;\;\;\frac{{\partial f\left({R{N_{ij1}},{\beta _3}} \right)}}{{\partial {\beta _3}}}}\\ { \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots } \end{array}}\\ { \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots } \end{array}}\\ { \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots }\\ { \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \cdots }\\ {\frac{{\partial f\left({R{N_{ij6}},{\beta _1}} \right)}}{{\partial {\beta _1}}}\;\;\;\frac{{\partial f\left({R{N_{ij6}},{\beta _2}} \right)}}{{\partial {\beta _2}}}\;\;\;\frac{{\partial f\left({R{N_{ij6}},{\beta _3}} \right)}}{{\partial {\beta _3}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {1\;\;\;\;5\;\;\;25}\\ {1\;\;\;\;9\;\;\;81} \end{array}}\\ {1\;\;11\;\;121} \end{array}}\\ {1\;\;15\;\;225}\\ {1\;\;19\;\;361}\\ {1\;\;27\;\;729} \end{array}} \right], $

$ {{\hat b}_k} = \left({\begin{array}{*{20}{c}} { -0.127265}\\ {0.012936}\\ { -0.000236} \end{array}} \right)。 $

因此，基于混合模型密度预测值为(0.419 415，0.498 343，0.531 975，0.587 575，0.627 623，0.661063)，绝对误差和均方根误差为(0.016 5，0.029 6)。而基于固定效应密度预测值为(0.489 55，0.523 95，0.538 75，0.563 55，0.581 95，0.599 55)，绝对误差和均方根误差为(0.047 8，0.064 2)。可以看出，考虑混合效应模型的预测精度明显高于只考虑固定效应的方法。

基于以上计算过程，利用独立检验数据进行验证。传统最小二乘法的3个指标值计算结果为平均误差0.048 0 、均方根误差0.066 3 、确定系数0.441 6 。当只考虑单木效应时的混合模型平均误差为0.043 2 、均方根误差为0.061 4 、确定系数为0.5219 ;当只考虑树高效应时的混合模型平均误差为 0.029 6 、均方根误差为0.044 7 、确定系数为0.746 3 ;当2个层次都考虑时的混合模型平均误差为0.000 4 、均方根误差为0.044 1 、确定系数为0.7503 。可以看出，引入随机参数提高了模型的预测精度。

1)本研究建立了落叶松木材密度模型，该模型能实现落叶松树干木材密度的径向和纵向预测;同时利用不同层次线性混合模型方法构建了木材密度混合模型。当只考虑单木效应时，b1作为随机参数时模型拟合最好;当只考虑高度效应时，b₁，b₂，b₃同时作为随机参数时模型拟合最好;当2个层次效应都考虑时，把样木效应含有随机参数b1i和树干不同高度效应含有随机参数b_1i，j，b_2i，j，b_3i，j的模型选为基于2层次混合模拟的最优混合模型。

2)对不同层次的混合效应比较表明:2层次效应和树干高度效应模拟与检验精度都高于单木效应。基于2层次效应和树干高度效应模拟结果相似，但基于2层次效应模拟精度略高于基于树干高度效应模拟精度。这也说明在拟合混合效应木材密度模型时，树木之间木材密度差异较小，误差的主要来源是树干不同高度之间木材密度的差异。这也符合数据的特征，因为本研究木材密度的样木来自立地条件相似的样地。

3)混合模型在实际应用中需要计算随机参数值，因此可应用基于树干不同高度效应模拟的参数估计值来代替2层次效应的参数估计值，即用单水平代替2水平，这样可以简化计算步骤。

4)Espinoza(2004)用相同数据收集的方法获取了云南石梓(Gmelina arborea)木材密度数据，研究了云南石梓木材密度的径向和纵向变异规律，但这只是从定性分析角度和简单线性回归分析来研究，虽然利用构建的简单线性模型能预测树干5个不同高度的木材密度值，但对其他树干高度应用时显示了模型的局限性，更没有利用混合模型技术进行木材密度模拟。

5)本研究将混合模型技术应用到落叶松木材密度模型构建中。这与以往的研究基本一致，如 Jordan等(2008)利用混合模型技术构建了美国南部6个区域火炬松(Pinus taeda)木材密度混合模型。Fujimoto等(2010)利用相同的方法研究了初植密度对日本落叶松(Larix kaempferi)木材密度的影响。但他们都只利用混合模型技术研究了树干胸高处木材密度变化规律，并没有考虑树干不同高度的随机效应，更没有对混合模型应用过程中随机参数计算过程进行详细阐述。