全 文 ：第 49 卷 第 7 期
2 0 1 3 年 7 月
林 业 科 学
SCIENTIA SILVAE SINICAE
Vol. 49，No. 7
Jul．，2 0 1 3
doi:10．11707 / j．1001-7488．20130713
收稿日期: 2012 － 08 － 15; 修回日期: 2012 － 11 － 16。
基金项目: 中央高校基本科研业务费专项资助(DL12DA01) ;国家“十二五”科技支撑课题(2012BAD22B0202) ;国家自然科学基金项目
(31170591)。
* 姜立春为通讯作者。
基于 2 层次线性混合模型的落叶松木材密度模拟*
李耀翔1 姜立春2
(1． 东北林业大学工程技术学院 哈尔滨 150040; 2． 东北林业大学林学院 哈尔滨 150040)
摘 要: 以黑龙江省七台河市林业局金沙林场 9 株人工落叶松 432 个样品密度数据为例，利用逐步回归技术构
建落叶松木材密度模型:WD = β1 + β2 RN + β3 RN
2 + β4 h。利用 S-PLUS 软件中的 LME 过程，分别考虑单水平和多水
平效应，拟合线性木材密度混合效应模型。结果表明:基于单水平和多水平效应的混合模型拟合精度高于传统的
基本模型，并且考虑单水平树高效应和 2 层次效应时的混合模型精度高于考虑单水平样木效应影响的混合模型。
模型检验结果表明:混合效应模型不但能反映总体平均木材密度变化趋势，还能反映分组之间的差异。
关键词: 木材密度; 两水平线性混合模型; 落叶松
中图分类号: S757 文献标识码:A 文章编号:1001 － 7488(2013)07 － 0091 － 08
Modeling Wood Density with Two-Level Linear Mixed Effects
Models for Dahurian Larch
Li Yaoxiang1 Jiang Lichun2
(1 ． College of Engineering and Technology，Northeast Forestry University Harbin 150040;
2 ． College of Forestry，Northeast Forestry University Harbin 150040)
Abstract: In this study，the sample data was based on 432 samples of 9 trees from dahurian larch( Larix gmelinii)
plantations located in Qitaihe Forest Bureau in Heilongjiang Province． The stepwise regression techniques were used to
develop wood density model: WD = β1 + β2 RN + β3 RN
2 + β4 h． Then，the developed model was fitted using single level
and multilevel linear mixed-effects modeling approach based on LME procedure of S-PLUS software． The mixed effects
models showed better model fitting results than basic model whatever considering single level and multilevel linear mixed
effects． Moreover，the mixed effects model considering height effects and both effects showed more precision than that
considering individual tree effects． Model test indicated that mixed effects models not only showed the mean trends of wood
density，but also showed the variations among groups．
Key words: wood density; two-level linear mixed model; dahurian larch (Larix gmelinii)
木材密度(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; 徐有明
林 业 科 学 49 卷
等，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 层次混合效应模型，并与单层混
合效应模型及传统模型拟合效果进行对比分析。
1 数据与方法
1. 1 数据
落叶松试样采自七台河市林业局金沙林场 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
数据 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
1. 2 方法
1. 2. 1 线性混合模型 线性混合模型包括单水平
和多水平形式，根据 Pinheiro 等(2000)，单水平线性
混合模型的形式为:
yi = Xiβ + Zibi + ε i，i = 1，…，m，
bi ～ N(0，D)，
ε i ～ N(0，Ri
{
)。
(1)
式中: yi 是第 i 株树或树干不同高度中木材密度的
测量值; m 是分组数; Xi 是已知设计矩阵; β是固定
参数向量; Zi 是随机效应的设计矩阵; bi 是随机参
数向量; D 是随机参数的方差协方差矩阵; Ri 是树
木内方差协方差结构; ε i 是模型的误差项，随机效
应参数 bi 和误差项 ε i 假定相互独立。
多水平线性混合模型 ( multilevel linear mixed
effects model，MLMEM)是基于单水平线性混合模型
基础上增加随机效应因子构成的。由多个相互嵌套
随 机 效 应 因 子 产 生 的 MLMEM 称 为 多 水 平
MLMEM。对于嵌套 2 水平的 MLMEM，模型表达
式为:
yij = Xijβ + Zi，j b i + Zij b ij + ε ij，
i = 1，…，m，j = 1，…，mi，
bi ～ N(0，D1)，bij ～ N(0，D2)，
ε ij ～ N(0，Ri





)。
(2)
式中: yij是第 i株树中树干高度 j处的木材密度测量
值; m 为第 1 水平的分组数; mi 为对应于第 1 水平
的第 2 水平的分组数; Xij 是已知设计矩阵; β 是固
定参数向量; Zi，j 和 Zij 分别是 1 水平和 2 水平的随
机效应设计矩阵; bi 和 bij 分别是 1 水平和 2 水平的
随机参数向量且相互独立; D1 和 D2 分别是 1 水平
和 2 水平随机参数的方差协方差矩阵; Ri 是方差协
方差结构; ε ij 是模型的误差项。
根据 Pinheiro 等(1998)，构建混合模型需要确
定 3 个步骤:1) 确定参数效应。本文主要考虑不同
随机参数的组合来比较(Pinheiro，2000)。2) 确定
29
第 7 期 李耀翔等: 基于 2 层次线性混合模型的落叶松木材密度模拟
树木内方差协方差结构( Ri )。采用式(3)来描述
(Davidian et al．，1995; Calama et al．，2004; Trincado
et al．，2006):
Ri = σ
2G0 ． 5i Γ i G
0 ． 5
i (3)
式中: σ2 指模型的误差方差值; Γ i 为树木内误差相
关性结构; Gi 为描述方差异质性的对角矩阵。
3) 确定随机效应的方差协方差结构(D)。比
较林业上常用的 3 种方差协方差结构:复合对称
(CS)、对角矩阵(diagonal matrix)、正定矩阵( general
positive-definite matrix) ( Fang et al．，2001; Jordan et
al．，2005)，以包括 3 个随机参数( b1，b2，b3 )的方差
协方差结构为例，
复合对称:
D =
σ2b1 + σ
2 σ2b1 σ
2
b1
σ2b1 σ
2
b1
+ σ2 σ2b1
σ2b1 σ
2
b1 σ
2
b1
+ σ

2
; (4)
对角矩阵:
D =
σ2b1 0 0
0 σ2b2 0
0 0 σ2b

3
; (5)
广义正定矩阵:
D =
σ2b1 σ b1b2 σ b1b3
σ b1b2 σ
2
b2 σ b2b3
σ b1b3 σ b2b3 σ
2
b

3
。 (6)
1. 2. 2 模型评价和检验指标 拟合和检验结果通
过以下指标评价:绝对误差 ( Bias)、均方根误差
(RMSE)、确定系数( R2 )。
Bias =
∑
m
i = 1∑
mi
j = 1∑
mij
k = 1
yijk － y^( )ijk
n
;
RMSE =
∑
m
i = 1∑
mi
j = 1∑
mij
k = 1
yijk － y^( )ijk
2
n －槡 1 ;
R2 = 1 － ∑
m
i = 1∑
mi
j = 1∑
mij
k = 1
yijk － y^( )ijk
2
∑
m
i = 1∑
mi
j = 1∑
mij
k = 1
yijk － 珋( )y

2
。
其中: yijk 为观测值; y^ ijk 为预测值; 珋y 为观测值的平
均值;m 为 1 水平的样木株数; mi 为 1 水平下 2 水
平的分组数; mij 为 1 水平下 2 水平对应研究对象的
观测次数;n 为样本数。
1. 2. 3 模型检验 模型检验使用 Vonesh 等(1997)
的方法来计算随机参数 bk:
b^ k ≈ D
^ Z^ Tk Z
^
kD
^ Z^ Tk + R
^( )
k
－1 e^ k。 (7)
式中:变量定义如前所述。
2 结果与分析
2. 1 基本模型的构建
初步分析木材密度与年龄、不同树干高度、胸
径、树高、冠幅等林木变量的关系，研究发现木材密
度与年龄和不同树干高度有密切关系，并随着年龄
的增加而增加，达到最大值之后有所下降。采用逐
步回归技术建立如下木材密度模型:
WD = β1 + β2 RN + β3 RN
2 + β4h。 (8)
式中: WD 为木材密度( g·cm － 3 );RN 为年轮数; h
为树干不同高度 ( m ); β1，β2，β3，β4 为模型预估
参数。
利用最小二乘法对模型(8)进行拟合，利用 t 检
验和 F 检验法分别对模型参数和模型进行检验。
参数 t 检验证明模型(8)所有参数估计值都是显著
的(P ＜ 0. 000 1)，说明这些参数显著影响模型的变
化。模型 F 检验(F = 87. 54，P ＜ 0. 000 1)都达到了
极显著水平。
2. 2 混合效应模型拟合
将木材密度模型(8)写成 2 水平混合效应模型
形式为:
WD ijk = (β1 + b1 i + b1 i，j) + (β2 + b2 i + b2 i，j)RN ijk +
(β3 + b3 i + b3 i，j)RN ijk
2 + β4hijk + ε ijk，
bi =
b1 i
b2 i
b3

i
～ N 0，D( )1 ，bi，j =
b1 i，j
b2 i，j
b3 i，

j
～ N(0，D2)，
ε ijk ～ N(0，Ri





)。
(9)
式中: WD ijk 是第 i 株树中树干高度 j 处的木材密度
的 k 次测量值;RN ijk是第 i 株树中树干高度 j 处的生
长轮的 k 次观察值; β1，β2，β3 是固定参数; bi 和 bi，j
分别是 1 水平和 2 水平的随机参数向量;其他变量
定义如前所述。
单水平是在 2 水平混合模型基础上去掉随机参
数 bi，j 。混合模型在拟合时，需要确定分类变量。
根据木材密度数据的特征，分类变量包括单木、树干
不同高度及二者同时作用。因此本研究分别考虑单
水平和 2 水平，单水平主要包括单木或树干不同高
度效应，2 水平指将二者同时考虑。当考虑单木效
应影响时，单木水平中每组对应的观察次数都是
48;当考虑树干不同高度效应时，每个树干高度水平
对应的观察次数都是 8;当同时考虑 2 水平时，每个
分组对应的观察次数也是 8。
2. 2. 1 基于单木效应混合模型拟合结果 拟合结果
39
林 业 科 学 49 卷
见表 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
模拟
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. 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 8 － 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
①参数个数包括固定参数个数、随机参数方差协方差参数个数及模型偏差参数。下同。The number of parameters includes a fixed number of
parameters，random pammeer variance and covariance parameters and model error parameters． The same below．
2. 2. 2 基于树干不同高度效应混合模型拟合结果
拟合结果见表 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，即把基本模型含有随机参数 b1，b2，b3 的模
型选为基于树干不同高度效应的最优混合模型。然
后又测试了林业上常用的 3 种方差协方差结构。通
过表 4 的 AIC，BIC，Log-likelihood 和似然比检验等
统计量指标表明:正定矩阵结构优于对角矩阵和复
合对称结构。
表 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 8 － 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
方差协方差结构
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
49
第 7 期 李耀翔等: 基于 2 层次线性混合模型的落叶松木材密度模拟
2. 2. 3 基于 2 层次混合模型拟合结果 综合考虑
单木效应和树干高度效应，这属于多层次混合模型
问题，与单层混合模型的研究方法相同，拟合结果见
表 5。可以看出，能够收敛的混合模型共有 6 种。
当同时考虑 6 个随机参数时，模型不收敛;当考虑 5
个随机参数组合时，模拟 2 的 AIC 和 BIC 值最小;当
考虑 4 个随机参数组合时，模拟 5 的 AIC 和 BIC 值
最小。综合考虑不同随机参数的组合，模拟 5 的效
果最好。LRT 检验表明:基本模型和模拟 2 有显著
不同(P ＜ 0. 000 1)，模拟 2 和模拟 5 没有显著不同
(P = 0. 61)。综合以上分析，把模拟 5，即把样木效
应含有随机参数 b1 i 和树干不同高度效应含有随机
参数 b1 i，j，b2 i，j，b3 i，j 的模型选为基于 2 层次混合模拟
的最优混合模型。
表 5 基于单木效应和高度效应混合模型拟合结果①
Tab． 5 Fitting results of mixed model based on individual tree and height effects
模拟
Simulation
单木效应 Individual tree effects 高度效应 Height effect
b1 i b2 i b3 i b1 i，j b2 i，j b3 i，j
AIC BIC Log-likelihood
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
①带△变量为包括这个随机参数。With △ variables including the random parameters．
2. 2. 4 误差的异方差性和自相关性 误差的异质性
主要通过残差分布图进行比较。可以看出，与基本模
型的残差分布图相比(图 1A)，基于单木效应混合模型
残差分布图变化不明显(图 1B)，这也说明，单木效应
对木材密度影响不大。而基于树干高度效应和基于 2
层次混合效应模型的残差分布图(图 1C 和 D)显示了
更均匀的分布。误差相关性测试了林业上常用的 4 种
描述时间序列的自相关结构［CS，AR(1)，MA(1)和
ARMA(1，1)］，结果显示没有提高混合模型的精度，因
此时间序列自相关结构在本研究中没有进一步考虑。
表 6 给出了落叶松不同木材密度模型的固定参
数估计值、随机参数方差协方差的参数估计值。
图 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)
59
林 业 科 学 49 卷69
表
6
不
同
模
型
参
数
和
方
差
组
成
估
计
值
①
Ta
b．
6
Pa
ra
m
et
er
an
d
va
ri
an
ce
co
m
po
ne
nt
s
es
tim
at
es
of
di
ff
er
en
t
m
od
el
s
参
数
Pa
ra
m
et
er
基
本
模
型
Ba
sic
m
od
el
混
合
模
型
( 单
木
)
M
ix
ed
m
od
el
(t
re
e)
混
合
模
型
( 高
度
)
M
ix
ed
m
od
el
(h
ei
gh
t)
混
合
模
型
(2
层
次
)
M
ix
ed
m
od
el
(t
wo
le
ve
l)
固
定
参
数
Fi
xe
d
pa
ra
m
et
er
β 1
0.
44
3
42
(0
.0
12
4)
*
0.
44
0
92
(0
.0
14
7)
*
0.
43
1
70
(0
.0
18
7)
*
0.
43
9
14
(0
.0
20
7)
*
β 2
0.
01
1
46
(0
.0
01
4)
*
0.
01
1
87
(0
.0
01
3)
*
0.
01
3
60
(0
.0
02
1)
*
0.
01
3
68
(0
.0
02
0)
*
β 3
－
0.
00
0
19
(0
.0
00
04
)*
－
0.
00
0
20
(0
.0
00
04
)*
－
0.
00
0
25
(0
.0
00
05
)*
－
0.
00
0
26
(0
.0
00
05
)*
β 4
－
0.
00
4
51
(0
.0
00
6)
*
－
0.
00
4
54
(0
.0
00
6)
*
－
0.
00
5
00
(0
.0
00
9)
*
－
0.
00
6
00
(0
.0
00
7)
*
随
机
效
应
方
差
协
方
差
结
构
Ra
nd
om
ef
fe
ct
va
ria
nc
e-
co
va
ria
nc
e
D
=
［
0．
00
0
6］
D
=
0．
00
8
7
－
0．
00
0
9
2．
1E
－
05
－
0．
00
0
9
0．
00
0
1
－
2．
5E
－
06
2．
1E
－
05
－
2．
5E
－
06
5．
[
]
6E
－
08
D
1
=
0．
[
]
00
0
5
D
2
=
0．
01
0
0
－
0．
00
1
0
2．
3E
－
05
－
0．
00
1
0
0．
00
0
1
－
2．
5E
－
06
2．
3E
－
05
－
2．
5E
－
06
5．
[
]
5E
－
08
模
型
方
差
M
od
el
va
ria
nc
e
σ2
0.
00
4
4
0.
00
3
8
0.
00
2
5
0.
00
2
5
①
括
号
内
为
标
准
差
In
pa
re
nt
he
se
s
in
th
e
sta
nd
ar
d
de
vi
at
io
n;
*
:
P
＜
0.
00
0
1．
第 7 期 李耀翔等: 基于 2 层次线性混合模型的落叶松木材密度模拟
2. 3 模型应用
本研究以 1. 3 m 处木材密度预测的一个例子来
说明:从检验数据中树高 1. 3 m 处随机抽取年龄和
密度的 6 次测量值 ( 5，0. 35 )，( 9，0. 50 )，( 11，
0. 53)，(15，0. 59)，(19，0. 63 )，(27，0. 64 )来计
算。由于树干高度效应和 2 层次效应对混合模型拟
合结果影响相似，因此参数和随机效应方差协方差
估计值采用基于高度效应的参数估计值。首先计算
e^k，R
^
i 和 Z
^
k，然后利用式(7)计算随机参数 b
^
k:
e^k =
0 ． 35
0 ． 50
0 ． 53
0 ． 59
0 ． 63
0 ．

64
－
0 ． 487 25
0 ． 527 65
0 ． 544 85
0 ． 573 25
0 ． 593 65
0 ．

610 45
=
－ 0 ． 137 25
－ 0 ． 027 65
－ 0 ． 014 85
0 ． 016 75
0 ． 036 35
0 ．

029 55
，
R^ i = σ
2 I(6) = 0 ． 002 5
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0

0 0 0 0 0 1
，
Z^ k =
f(RNij1，β1)
β1
f(RNij1，β2)
β2
f(RNij1，β3)
β3
… … …
… … …
… … …
… … …
f(RNij6，β1)
β1
f(RNij6，β2)
β2
f(RNij6，β3)
β


3
=
1 5 25
1 9 81
1 11 121
1 15 225
1 19 361

1 27 729
，
b^ k =
－ 0 ． 127 265
0 ． 012 936
－ 0 ．

000 236
。
因此，基于混合模型密度预测值为(0. 419 415，
0. 498 343，0. 531 975，0. 587 575，0. 627 623，0. 661
063)，绝对误差和均方根误差为(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. 521
9;当只考虑树高效应时的混合模型平均误差为
0. 029 6、均方根误差为 0. 044 7、确定系数为 0. 746
3;当 2 个层次都考虑时的混合模型平均误差为
0. 000 4、均方根误差为 0. 044 1、确定系数为 0. 750
3。可以看出，引入随机参数提高了模型的预测
精度。
3 结论与讨论
1) 本研究建立了落叶松木材密度模型，该模型
能实现落叶松树干木材密度的径向和纵向预测;同
时利用不同层次线性混合模型方法构建了木材密度
混合模型。当只考虑单木效应时，b1 作为随机参数
时模型拟合最好;当只考虑高度效应时，b1，b2，b3 同
时作为随机参数时模型拟合最好;当 2 个层次效应
都考虑时，把样木效应含有随机参数 b1 i 和树干不同
高度效应含有随机参数 b1 i，j，b2 i，j，b3 i，j 的模型选为基
于 2 层次混合模拟的最优混合模型。
2) 对不同层次的混合效应比较表明:2 层次效
应和树干高度效应模拟与检验精度都高于单木效
应。基于 2 层次效应和树干高度效应模拟结果相
似，但基于 2 层次效应模拟精度略高于基于树干高
度效应模拟精度。这也说明在拟合混合效应木材密
度模型时，树木之间木材密度差异较小，误差的主要
来源是树干不同高度之间木材密度的差异。这也符
合数据的特征，因为本研究木材密度的样木来自立
地条件相似的样地。
3) 混合模型在实际应用中需要计算随机参数
值，因此可应用基于树干不同高度效应模拟的参数
估计值来代替 2 层次效应的参数估计值，即用单水
平代替 2 水平，这样可以简化计算步骤。
4) Espinoza(2004)用相同数据收集的方法获取
了云南石梓 (Gmelina arborea)木材密度数据，研究
了云南石梓木材密度的径向和纵向变异规律，但这
只是从定性分析角度和简单线性回归分析来研究，
虽然利用构建的简单线性模型能预测树干 5 个不同
高度的木材密度值，但对其他树干高度应用时显示
了模型的局限性，更没有利用混合模型技术进行木
材密度模拟。
79
林 业 科 学 49 卷
5) 本研究将混合模型技术应用到落叶松木材
密度模型构建中。这与以往的研究基本一致，如
Jordan 等(2008)利用混合模型技术构建了美国南部
6 个区域火炬松(Pinus taeda)木材密度混合模型。
Fujimoto 等(2010)利用相同的方法研究了初植密度
对日本落叶松( Larix kaempferi)木材密度的影响。
但他们都只利用混合模型技术研究了树干胸高处木
材密度变化规律，并没有考虑树干不同高度的随机
效应，更没有对混合模型应用过程中随机参数计算
过程进行详细阐述。
参 考 文 献
姜景民，孙海菁，吕本树 ． 1999． 火炬松木材基本密度的株内变异 ．
林业科学研究，13(1) :97 － 102．
刘青华，张 蕊，金国庆，等 ． 2010． 马尾松年轮宽度和木材基本密
度的种源变异及早期选择 ． 林业科学，46(5) : 49 － 54．
骆秀琴，管 宁，文小明，等 ． 1999． 木材材性株内径向变异模式初
探 VI． 19 个杉木种源木材密度径向变异模式的研究马尾松年轮
宽度和木材基本密度的种源变异及早期选择 ． 林业科学，35
(6) : 86 － 92．
罗建勋，李晓清，孙 鹏，等 ． 2004． 云杉天然群体管胞和木材基本
密度性状变异的研究 ． 北京林业大学学报，26(6) : 80 － 85．
马丽娜，付孝德，张 明，等 ． 2003． 人工林杨树木材密度变异规律
的研究 ． 安徽农业大学学报，30(4) : 410 － 413．
徐有明，林 汉，江泽慧，等 ． 2002． 橡胶树生长轮宽度、木材密度变
异及其预测模型的研究，林业科学，38(1) : 95 － 102．
王秀花，马丽珍，马雪红，等 ． 2011． 木荷人工林生长和木材基本密
度 ． 林业科学，47(7) : 138 － 144．
Calama R，Montero G． 2004． Interregional nonlinear height-diameter
model with random coefficients for stone pine in Spain． Canadian
Journal Forest Research，34(1) :150 － 163．
Davidian M， Giltinan D M． 1995． Nonlinear models for repeated
measurement data． Chapman and Hall，London，UK．
Espinoza J A． 2004． Within-tree density gradients in Gmelina arborea in
Venezuela． New Forest，28(2 /3) : 309 － 317．
Fang Z，Bailey R L． 2001． Nonlinear mixed effects modeling for slash
pine dominant height growth following intensive silvicultural
treatments． Forest Science，47(3) :287 － 300．
Fujimoto T，Koga S． 2010． An application of mixed-effects model to
evaluate the effects of initial spacing on radial variation in wood
density in Japanese larch ( Larix kaempferi ) ． Journal of Wood
Science，56(1) : 7 － 14．
Jaakkola T，Mkinen H，Saranp P． 2005． Wood density in Norway
spruce: changes with thinning intensity and tree age． Canadian
Journal Forest Research，35(7) :1767 － 1778．
Jordan L，Daniels R F，Clark A，et al． 2005． Multilevel nonlinear
mixed effects models for the modeling of earlywood and latewood
microfibril angle． Forest Science，51(4) :357 － 371．
Jordan L，Clark A，Schimleck L R，et al． 2008． Regional variation in
wood specific gravity of planted loblolly pine in the United States．
Canadian Journal Forest Research，38(4) :698 － 710．
Lundgren C． 2004． Microfibril angle and density patterns of fertilized and
irrigated Norway spruce． Silva Fennica，38(1) : 107 － 117．
MacDonald E，Hubert J． 2002． A review of the effects of silviculture on
the timber quality of Sitka spruce． Forestry，75(2) : 107 － 138．
Mkinen H，Saranp P，Linder S． 2002． Wood density variation of
Norway spruce in relation to nutrient optimization and fibre
dimensions． Canadian Journal Forest Research，32(2) :185 － 194．
Mkinen H，Jaakkola T，Piispanen R，et al． 2007． Predicting wood and
tracheid properties of Norway spruce． Forest Ecology and
Management，241(1 /3) : 175 － 188．
Molteberg D，Hib O． 2007． Modelling of wood density and fibre
dimensions in mature Norway spruce． Canadian Journal Forest
Research，37(8) :1373 － 1389．
Pérez Cordero L D，Kanninen M． 2002． Wood specific gravity and
aboveground biomass of Bombacopsis quinata plantations in Costa
Rica． Forest Ecology and Management，165(1 /3) : 1 － 9．
Pinheiro J C，Bates D M． 1998． Model building for nonlinear mixed
effects models． Tech Rep Dep of Stat，Univ of Wisconsin．
Pinheiro J C，Bates D M． 2000． Mixed-effects models in S and S-PLUS．
Springer，New York，528．
Schneider R，Zhang S Y，Swift D E，et al． 2008． Predicting selected
wood properties of jack pine following commercial thinning．
Canadian Journal Forest Research，38(7) :2030 － 2043．
Trincado G，Burkhart H E． 2006． A generalized approach for modeling
and localizing stem profile curves． Forest Science， 52 ( 6 ) :
670 － 682．
Vonesh E F，Chinchilli V M． 1997． Linear and nonlinear models for the
analysis of repeated measurements． Marcel Dekker，New York．
(责任编辑 石红青)
89