全 文 ：第 51 卷 第 3 期
2 0 1 5 年 3 月
林 业 科 学
SCIENTIA SILVAE SINICAE
Vol. 51，No. 3
Mar．，2 0 1 5
doi: 10．11707 / j．1001-7488．20150304
收稿日期: 2014 － 05 － 20 ; 修回日期: 2014 － 09 － 19。
基金项目: 国防科工局重大专项项目(21 － Y30B05 － 9001 － 13 /15)。
* 倪成才为通讯作者。
基于混合效应模型及 EBLUP预测美国黄松
林分优势木树高生长过程*
祖笑锋1，2 倪成才2 Gorden Nigh3 覃先林1
(1．中国林业科学研究院资源信息研究所 北京 100091; 2．北华大学林学院 吉林 132013;
3． British Columbia Ministry of Forests，Lands and Natural Ｒesources Operations，Forest Analysis and Inventory Branch，
P． O． BOX9512，Stn． Prov． Govt． Victoria，B． C． V8W 9C2，Canada)
摘 要: 【目的】基于加拿大哥伦比亚省美国黄松 79 株解析木数据，研究如何用经验线性无偏最优预测法
(EBLUP)预测优势木树高生长过程，并分析预测精度与观测次数、观测间隔和预测时长的关系。【方法】随机抽取
49 株解析木数据拟合树高生长混合效应模型，30 株解析木数据用于 EBLUP 的预测分析。树高生长模型以
Ｒichards，Logistic，Korf 等为基础模型，选用 AIC，BIC 及 Loglik 3 个统计量评价模型的拟合效果。模型拟合用 Ｒ 软件
的 nlme 函数实现，预测分析以预测误差均方(MSPE)为评价标准。在分析观测间隔、观测次数和预测时长对 MSPE
的影响时，为分离出 1 个因素的影响效果，将 2 个因素保持不变，以分析第 3 个因素的影响作用。在 Ｒ 软件拟合结
果的基础上，用 SAS 的 IML 过程进行 EBLUP 预测分析。【结果】拟合结果表明，Logistic 方程的拟合精度最高，选为
EBLUP 预测分析的基本模型。预测分析结果表明，观测次数、观测间隔和预测时长对预测精度均有显著影响。随
着观测次数的增加，MSPE 一般表现出减少的趋势，但下降幅度与观测间隔有关:当间隔较大时，不同的观测值可以
提供更充分的生长过程信息，因而可以显著降低 MSPE 值;但当间隔较小时，观测值所提供的生长信息相互重叠，
对提高预测精度的增益有限。从预测时长角度看，在观测值附近一定区域内，EBLUP 预测结果非常精确，但随着预
测时长增加，预测误差呈逐渐增加的趋势。【结论】EBLUP 预测相当于两阶段拟合过程的第二阶段。第一阶段拟
合为估计混合参数模型确定参数的过程，而第二阶段则是在第一阶段拟合结果的基础上，依据一个特定林分的若
干树高观测值用 EBLUP 法预测此林分的随机效应值，并进一步预测树高生长过程。
关键词: 混合效应模型; 经验线性无偏最优预测法; 树高生长模型; 美国黄松
中图分类号:S758 文献标识码:A 文章编号:1001 － 7488(2015)03 － 0025 － 09
Based on Mixed-Effects Model and Empirical Best Linear
Unbiased Predictor to Predict Growth Profile of Dominant Height
Zu Xiaofeng1，2 Ni Chengcai2 Gorden Nigh3 Qin Xianlin1
(1 ． Institute of Forest Ｒesources Information Techniques，CAF Beijing 100091; 2 ． College of Forestry，Beihua University Jilin 132013;
3 ． British Columbia Ministry of Forests，Lands and Natural Ｒesources Operations，Forest Analysis and Inventory Branch，P． O． BOX9512，
Stn． Prov． Govt． Victoria，B． C． V8W 9C2，Canada)
Abstract: 【Objective】 This study analyzed prediction accuracy of empirical best linear unbiased predictor(EBLUP)，
and effects of previous observations，age interval of observations and prediction span on prediction accuracy，based upon
height data from 79 dominant trees of ponderosa pine in British Columbia，Canada． 【Method】We randomly selected 49
trees for fitting mixed-effects models and 30 trees for validating EBLUP． The base models were Ｒichards，Logistic，and
Korf． Fit statistics，AIC，BIC and Loglik，were used as evaluation criteria，and mean squared prediction error (MSPE)
for analyzing effects of previous observations，age interval of observations and prediction span on prediction accuracy． We
used the nlme function in Ｒ for model fitting，and the IML procedure in SAS for analyzing EBLUP prediction． To isolate
the effect of one factor，we kept two other factors fixed．【Ｒesult】Fitting results showed the Logistic model had the best
criteria among the three models of under investigation，indicating that it was the best-fitted model and was chosen for
林 业 科 学 51 卷
EBLUP prediction analysis． In the analysis of EBLUP prediction，we first introduced how to use EBLUP to predict random
effects associated with a stand through a detailed example． Data from six trees，which deviated significantly from
population-mean growth process，were used to present relationships among individual growth，population-mean growth，
and adjusted values given by EBLUP． The results indicated that EBLUP prediction could fully follow individual growth
process，given that there were multiple previous observations with long-enough age intervals． EBLUP analysis results also
presented the number of previous observations，age interval of observations and prediction span significantly affected
prediction accuracy． MSPE decreased as the number of previous observations increased，particularly when observations
separated long enough in age so that they could give more efficient growth information． With respect to prediction span，
prediction accuracy decreased as prediction span extended further away from the ages at which observations were obtained．
【Conclusion】We concluded that EBLUP could be taken as the second stage of a two-stage fitting process． The first stage
was used to estimate fixed model parameters，whereas the second stage to predict random effects associated with a stand on
the basis of parameter estimators obtained in the first stage，and then to predict the height growth process of the stand．
Key words: mixed-effects model; EBLUP; dominant height growth model; ponderosa pine
林分优势木树高生长模型是描述优势木树高生
长过程的统计模型，是林分生长与收获模型系统中
的重要组成部分。林分优势木树高的预测为适地适
树、森林经营活动和林分生长收获预测提供重要的
基础数据(孟宪宇，1994)。
由于林分立地质量不同，同一树种的不同林分
优势木树高生长过程也不同，表现为生长模型的参
数不同。基本有 2 种方法可以处理模型参数的变
化:方法一是将模型参数设定为其他林分测树因子
的函数，但这一方法模型结构复杂，参数过多，在预
测树高时首先要预测树高模型中的其他测树因子，
因此可能会带来更大的预测误差;另一种方法是将
一个或多个参数设定为随林分变化而变化的随机参
数，主要包括差分生长模型、广义差分生长模型和混
合效应模型。
差分生长模型 ( algebraic difference approach，
ADA)由 Clutter(1963)提出，并由 Bailey 等 (1974)
进行了系统化研究，成为立地指数模型的主流方法。
ADA 模型一般仅能将一个模型参数设定为随机参
数，用一个观测值即可估计此随机参数;但在有多个
观测值时，ADA 模型无法有效利用所有观测值中蕴
涵的生长过程信息。此外，ADA 模型的树高生长过
程，要么为单型曲线族，要么为多型曲线族，无法生
成具有可变水平渐进极值的多型曲线族。针对这一
不足，Ciezewski 等(2000a; 2000b)提出了广义差分
生长模型( generalized ADA，GADA)，GADA 模型可
以指定多个参数为随机参数，一定程度上克服了
ADA 模型的不足。但与 ADA 模型一样，GADA 模
型在模型拟合时，仍然无法考虑林分个体与总体间
的相关性和差异性;在模型预测时，无法有效利用所
有观测值中蕴涵的生长过程信息，而且也无法给出
总体平均生长趋势。
混合模型方法则在很大程度上弥补了 ADA，
GADA 模型的不足:在模型拟合时，可以指定任意多
个模型参数为随机参数，既可以反映总体的平均变
化趋势，又可以反映个体之间的差异;在模型预测
时，可以有效利用多个观测值估计特定林分的随机
效应参数值，并给出精确的预测。
Lappi 等(1988)开创性地建立了 Ｒichards 混合
效应优势木树高模型。Fang 等(2001)建立了美国
乔治亚州和佛罗里达州湿地松 ( Pinus elliottii) 的
Ｒichards 混合效应模型，考虑了经营措施(如采伐、
施肥和火烧)的影响和样地的随机效应，认为非线
性混合模型要比传统方法预测的精度高。Calegario
等(2005)考虑了样地的随机效应，为无性系杂交桉
树( eucalyptus hybrid)构造了优势树高非线性混合
模型，结果表明采用三参数 Logistic 非线性混合模型
具有相当大的灵活性，产生了依赖于无性系和环境
的多形地位指数曲线。目前国内对混合效应模型的
拟合方法进行了系统的研究( Tang et al．，2001; 李
永慈等，2004; 李春明等，2010)，但是混合模型在
林业上的预测应用却不多见。
混合模型主要用于分析定性、定量和随机因子
对应变量的影响，特别是各种随机因子对应变量的
影响程度;另一个重要用途就是在有多个观测值时，
用经验线性无偏最优预测法 ( empirical best linear
unbiased predictor，EBLUP)提高预测精度。EBLUP
最早应用于动物育种学中，与地统计学的 Kriging、
时间序列的卡尔曼滤波及小域估计法 ( small area
estimation)(吕萍等，2009)的数学原理一致。目前
国内林业上对混合效应模型的应用基本限于与传统
回归模型拟合效果的对比，对预测的研究成果较少。
62
第 3 期 祖笑锋等: 基于混合效应模型及 EBLUP 预测美国黄松林分优势木树高生长过程
本文基于加拿大哥伦比亚省美国黄松 ( Pinus
ponderosa)79 株解析木数据，研究如何用 EBLUP 对
优势木树高生长过程进行预测，使用 49 株解析木的
树高连年生长数据拟合 Ｒichards-Chapman，Logistic，
Korf 等混合效应模型，利用 30 株解析木数据进行
EBLUP 预测分析。模型拟合用 Ｒ 软件的 nlme 函数
完成，并用 SAS 的 IML 过程对 EBLUP 预测进行
分析。
1 数据来源与方法
1. 1 数据来源
数据 来 源 于 加 拿 大 哥 伦 比 亚 省 ( British
Columbia)美国黄松林分的临时样地数据。哥伦比
亚省林业局于 2000 和 2001 年夏天建立了 100 块临
时样地，样地设置时充分考虑了现有林分的立地质
量、年龄、密度以及物种丰富度等因素，从中选取具
有代表性的林分作为样地。样地形状为圆形，面积
为 0. 01 hm2(半径 5. 64 m)。每个样地挑选 1 株未
受损或未受压的健康优势木或亚优势木作为样木，
最终得到 79 株优势木的解析木。样木伐倒后，对于
树桩和树干部分采取纵切的方法，由髓心直接查出
树高 －年龄值;对于不适合做纵切的树梢部分采取
轮生枝法确定每年的树高生长量。随机选取 49 株
优势木进行拟合，30 株优势木进行预测。表 1 分别
列出拟合数据(49 株)和预测数据(30 株)的描述性
统计量，数据的详细介绍见文献 Nigh(2004)。
1. 2 树高生长模型的选择与建立
林分优势树高生长模型种类很多，如 Logistic，
Mitscherlich，Gompertz，Ｒichards 和 Korf 方程等，其中
Ｒichards-Chapman，Logistic 和 Korf 方程应用较为广
泛，其模型结构见表 2。表中的 yi，j 表示 i 个林分在
年龄 xi，j时的树高; β1 ，β2 ，β3 为方程的确定效应参
数; bi，1 ，bi，2 ，bi，3 为与 β1 ，β2 ，β3 对应的随机效应
参数，可假设服从三维正态分布; ei，j 为方程的误差
项，且与随机效应参数 bi，1 ，bi，2 ，bi，3 任何一个均相
互独立。
表 1 拟合和预测优势木树高生长模型的数据统计
Tab． 1 The summary statistics of the fitting and calibration data for the dominant tree
株数
Number
均值
Mean
最小值
Minimum
最大值
Maximum
标准差
SD
拟合数据 Fitting data
优势木树高 Dominant height /m 49 24. 2 16. 9 35. 8 5. 6
年龄 Age / a 49 135 83 243 38
验证数据 Calibraton data
优势木树高 Dominant height /m 30 23. 3 14. 7 39. 5 6. 0
年龄 Age / a 30 117 84 206 25
表 2 优势木树高生长方程和混合效应模型
Tab． 2 The equations of dominant tree and mixed-effects model
模型
Model
基础模型
Base model
混合效应模型
Mixed-effects model
Ｒichards-Chapman yi，j = β1［1 － Exp( － β2 xi，j)］β3 + ei，j y i，j = ( β1 + bi，1 ) {1 － Exp［－ ( β2 + bi，2 ) xi，j］}
( β3 + b i，3) + ei，j
Logistic yi，j =
β1
1 + Exp［β2 + β3 lg( xi，j)］
+ ei，j y i，j =
( β1 + bi，1 )
1 + Exp( β2 + bi，2 ) + ( β3 + bi，3 ) lg( xi，j)
+ ei，j
Korf yi，j = β1 Exp
－ β2
xβ3i，
( )
j
+ ei，j y i，j = ( β1 + bi，1 )Exp
－ ( β2 + bi，2 )
x( β3 + b i，3)i，
[ ]
j
+ ei，j
1. 3 EBLUP 的基本原理
EBLUP 的推导方法有多种，下面介绍易于理
解、基于多元正态分布理论的推导方法。
以 b i = 0 为基点，用一阶泰勒公式可以将非线
性混合模型
y i = f(A iβ + B ib i，X i) + e i
近似地表达为:
y i ≈ f(A iβ，X i) + Z ib i + e i。 (1)
式中: X i 为非线性混合模型的预测变量矩阵; y i 为
与 X i 的元素 xi1 ，xi2 ，…，xin i 相对应的 ni 个树高观
测值构成的向量，y i的第 j个元素为第 i个林分的第
j个观测值; e i 为与 y i 相关联的随机误差项向量，服
从 ni维正态分布，其数学期望为 0 ; β 为 p维确定参
数列向量; b i 为 q 维随机参数列向量，一般假设 b i
服从 q 维的多元正态分布，其数学期望为 0 ，而且与
e i 相互独立; A i 和 B i 为具有适当维数的关联矩阵
(即由 0 或 1 构成的矩阵)。以表 2 中的混合效应模
72
林 业 科 学 51 卷
型为例，对应的 A i 和 B i 均为三阶单位阵，模型参数
均可表示为 β k + bi，k( k = 1，2，3 )。以上所述矩阵可
分别表示如下:
β =
β1
β2

β


p
，b i =
bi，1
bi，2

bi，


q
，X i =
xi，1
xi，2

xi，n


i
，
y i =
yi，1
yi，2

yi，n


i
，e i =
ei，1
ei，2

ei，n


i
。
随机向量 b i ～ N(0，D)，e i ～ N(0，Ｒ i);e i 和 b i
协方差矩阵分别为 Ｒ i 和 D ;矩阵 Ｒ i 是一个维数为
ni 的对称矩阵，而 D 为 q 维对称矩阵。式(1)中的
矩阵 Z i 定义如下:
Zi =
zi，1，1 zi，1，2 … zi，1，q
z i，2，1 zi，2，2 … zi，2，q
  … 
zi，n i，1 zi，n i，2 … zi，n i，


q
，
zi，j，k
f(A iβ + B ib i，xi，j
bi，k ( b i = 0)
。 (2)
此处 k = 1，2，…，q 。对于元素 zi，j，k ，对式(1)求关
于随机效应参数 bi，k (此处 k = 1，2，…，q )的偏导
数，求出偏导数后，令所有随机效应参数取值为 0 即
可求出。
由于 f(A iβ，X i) 为确定的量，且由于 e i 和 b i 相
互独立，因此对式(1)求 y i 的协方差矩阵为:
Cov( y i) = Z iDZ
T
i + Ｒ i。
此处上标 T 表示矩阵转置。同样根据 e i 和 b i 相互
独立性，b i 和 y i 间的协方差矩阵为:
Cov(b i，y i) = Cov［b i，f(A iβ，X i) + Z ib i + e i］ =
Cov(b i，Z ib i) = DZ
T
i。
随机向量 b i 和 y i 的联合分布为:
b i
y( )i ～ N 0f(A iβ，X i( )) ， D DZ
T
i
Z iD
T Z iDZ
T
i + Ｒ( )[ ]i 。
假设 y i 的观测值向量 珓y i 含有 ni 个观测值，分别
为 珓yi，1，珓yi，2，…，珓yi，n i ，依次对应于 xi，1 ，xi，2 ，…，xi，n i，
根据多元正态分布理论，当 y i = 珓y i 时，b i 的条件分
布数学期望为:
E(b i | y i = 珓y i) = E(b i) + DZ
T
i (Z iDZ
T
i + Ｒ i)
－1
［珓y i － f(A iβ，X i)］。
对随机变量取值的最优预测值为其数学期望
值，因此随机变量 b i在 y i = 珓y i时的取值可以用 b i的
条件分布数学期望表示为:
b^ i = DZ
T
i (Z iDZ
T
i + Ｒi)
－1［珓yi － f(Aiβ，Xi)］。(3)
由式(3)可见，矩阵 Z i 与 (A iβ，X i) 一样，均为
参数向量 β 的函数。如果获得了 β 的估计值，即可
求出 Z i 与 f(A iβ，X i) 的估计值。D 在模型拟合时，
可由统计软件给出估计值。Ｒ i 的值随林分变化而变
化，在模型拟合时可以指定描述 ei，j 变化规律的方差
函数和自相关函数，在统计软件获取此类函数的参
数估计值后，可计算求出 Ｒ i 估计值。假设矩阵 β ，
D 及 Ｒ i 均已获得估计值，那么将式(3)中的所有未
知量如 β，D，Ｒ，Z i 及 f(A iβ，X i) 用相对应的估计量
取代后，则有:
b^ i = D^Z^
T
i ( Z^ i D^Z^
T
i + Ｒ^ i)
－1［珓y i － f(Aβ^ i，X i)］。(4)
式 ( 4 ) 称 为 经 验 线 性 无 偏 最 优 预 测 法
(EBLUP)，详细推导介绍请参考 Ni 等 (2011)。对
于一阶泰勒级数的基点选择，有 2 种基本方法:一是
本文介绍的 b i = 0 为基点;二是以 b i 的估计值为基
点。对 于 2 种 方 法 的 讨 论，可 参 考 Meng 等
(2009a)。混合效应模型在林业中的应用多以第 1
种方法为主 ( Fang et al．，2001; Hall et al．，2001;
Huang et al．，2009b; Jiang et al．，2010)，对于第 2 种
方法的应用可参考 Hall 等 ( 2004 ) 和 Huang 等
(2009a)。一般认为第二种方法迭代过程繁琐费
时，比第一种方法仅存在微小的优势。
2 结果与分析
2. 1 树高生长模型的比较
采用 Ｒichards-Chapman，Logistic，Korf 3 个生长
模型作为基本模型，利用 Ｒ 软件的 nlme 过程对 49
株解析木进行拟合，其拟合统计量见表 3。
使用 AIC，BIC 和 Loglik 3 个模型评价指标对
拟合结果进行评价。由表 3 可以看出，Logistic 方
程的 AIC，BIC 最小且 loglik 值最大，表明 Logistic
方程拟合效果最好，因此选为混合效应模型的基
本模型。
表 3 拟合统计量
Tab． 3 Fit statistics
模型
Model
AIC BIC Loglik
Ｒichards-Chapman 606. 348 638. 109 － 293. 174
Logistic 602. 390 634. 151 － 291. 195
Korf 639. 391 671. 153 － 309. 695
2. 2 拟合非线性混合效应模型
Logistic 方程不仅拟合精度高，而且具有良好的
82
第 3 期 祖笑锋等: 基于混合效应模型及 EBLUP 预测美国黄松林分优势木树高生长过程
解析性、预测性以及广泛的适应性。按照 Pinheiro 等
(2000)、Davidian 等(1995)以及 Littell 等(1996)，表 2
中 Logistic 混合效应方程的 β1 + bi，1 =  i，1 表示当年
龄趋向无穷时应变量的水平渐近线; β2 + bi，2 =  i，2
是年龄在  i，2 /2 时的应变量的取值; β3 + bi，3 =  i，3
是规模参数，表示拐点与 0. 73  i，1 点在 x 轴上的
距离。
利用 Ｒ 软件的 nlme 函数对 Logistic 模型参数进
行估计，用于 EBLUP 运算的统计量列于表 4。由于
假设随机误差项 ei，j 为独立等方差结构，因而 ei，j ～
N(0，σ2) ，仅需估计 σ2 的值即可。nlme 函数给出
σ2 估计值为 0. 558 2，但未列于表 4。
表 4 模型拟合结果
Tab． 4 The statistical results of the model fitting
固定效应参数
Fixed-effects parameter
随机效应参数协方差矩阵
Variance -covariance matrix D
^
参数
Parameter
估计值
Estimate
bi，1 bi，2 bi，3
β1 35. 023 6 bi，1 87. 455 3 － 4. 1842 1. 013 0
β2 9. 105 3 bi，2 － 4. 184 2 2. 124 0 － 0. 364 2
β3 － 2. 070 9 bi，3 1. 013 0 － 0. 364 2 0. 070 8
已知一个林分，假设在 X i 年时已经观测了树高
值 珓y i ，可以利用式(4)预测此林分的 b i 值。y i 的第 j
个元素可以表示为:
yi，j =
(β1 + bi，1)
1 + Exp［(β2 + bi，2) + (β3 + bi，3)lg(xi，j)］
+ ei，j。
(5)
式中: β1 ，β2 ，β3 可用于描述总体平均树高生长过
程; bi，1 ，bi，2 ，bi，3 为与 β1 ，β2 ，β3 相关联的随机效
应参数，表示第 i 个林分树高生长过程参数与 β1 ，
β2 ，β3 间的随机偏离值。bi，1 ，bi，2 ，bi，3 的方差及协
方差估计值，可参考表 4。bi，1 的标准误较 bi，2 ，bi，3
的大很多，说明 49 株优势木的水平渐进线值变化幅
度较大，仅用确定效应，难以充分代表各株树的特定
生长过程。
2. 3 林分生长预测分析
2. 3. 1 单株示例分析 式(2)和(4)中 zi，j，k 的估计
值可以分别为式 (6 ) ～ (8)，此处 k = 1，2，3 。式
(4)中，列向量 f(A i β^，X i) 的第 j 个元素可由式(9)
表示:
z^ i，j，1 =
1
1 + Exp( β^2 + β^3 lgxij)
; (6)
z^ i，j，2 =
－ β^1
［1 + Exp( β^2 + β^3 lgxi，j)］
2
Exp( β^2 + β^3 lgxi，j); (7)
z^ i，j，3 =
－ β^1 lgxi，j
［1 + Exp( β^2 + β^3 lgxi，j)］
2
Exp( β^2 + β^3 lgxi，j); (8)
β^1
1 + Exp( β^2 + β^3 lgxij)
; (9)
在获得 bi，1 ，bi，2 ，bi，3 的预测值 b
^
i，1 ，b
^
i，2 ，b
^
i，3
后，即可用式(10)预测在年龄 xi，j 时的树高值。式
(6) ～ (10)中的 β^1 ，^β2 ，^β3 为参数 β1 ，β2 ，β3 的估
计值，同样由 Ｒ 软件给出，参考表 4。z^ i，j，k 可由式
(6) ～ (8)求出。本文用式 (10)进行预测，是由采
用的 EBLUP 算法所决定的［参考式(4)］，具体的原
因分析可参考 Ni 等(2011)。
y^i，j =
β^1
1 + Exp(β^2 + β^3 lgxij)
+Σ
3
k = 1
z^i，j，k b
^
i，k。 (10)
选取一样木的 3 对观测值用于预测 b i ，各项取
值及相应的计算结果列于表 5。表 5 中，珓y i 为在 X i
时的优势木树高观测值; f(A i β^，X i) 用式(9)求出，
所需参数估计值 β^ k( k = 1，2，3) 列于表 4，Z^ i 各列
元素依次按式 (6) ～ (8)计算。将 X i 列的值 (50，
80，120)代入式(6) ～ (8)可获得 Z^ i 各列的取值; D^
估计值由表 4 给出;在模型拟合时，仅指定 ei，j 为独
立的同质方差结构，故 Ｒ^ i 由 σ^
2 (值为 0. 558 2)与三
维单位阵相乘而得。根据表 5 数据和式(4)所计算
的 b i 值，同列于表 5。
表 5 EBLUP 预测的矩阵和向量
Tab． 5 The matrices and vectors used for EBLUP prediction
j 珓y i /m X i / a f(A i β^，X i) /m Z^ i Ｒ^ i b
^
i
1 6. 49 50 9. 392 3 0. 268 1 － 25. 631 3 － 100. 270 3 0. 558 2 0 0 － 11. 104 4
2 12. 18 80 17. 243 9 0. 492 3 － 17. 779 6 － 77. 910 9 0 0. 558 2 0 0. 153 9
3 16. 86 120 24. 233 3 0. 691 9 － 10. 790 2 － 51. 658 31 0 0 0. 558 2 － 0. 039 9
根据表 5 的 b
^
i ，预测了此样木年龄 80 ～ 83 间
的 4 个树高值，并与实测值同列于表 6。总体平均
值由式(9)计算而得，对应于 f(A i β^，X i) ，而 EBLUP
预测值由式(10)计算。校正值为 EBLUP 预测值与
总体平均值的差，是根据表 5 中的 3 个观测值所提
供个体生长过程信息，对总体平均值进行校正，从而
给出了此样木的树高预测值。以 82 年时的优势木
树高为例，其校正值约为 － 5. 301 1。
92
林 业 科 学 51 卷
表 6 单株优势木树高观测值、总体平均值、EBLUP 预测值及校正值
Tab． 6 Observations of dominant tree，population mean，and EBLUP predicted，and adjusted values for a selected single trees
年龄
Age / a
优势木树高
Dominant tree /m
总体平均值
Population means /m
EBLUP 预测值
EBLUP prediction
校正值
Adjusted values
80 12. 18 17. 243 0 12. 076 7 － 5. 167 2
81 12. 25 17. 469 2 12. 234 6 － 5. 234 5
82 12. 43 17. 691 7 12. 390 5 － 5. 301 1
83 12. 54 17. 911 4 12. 544 6 － 5. 366 8
选取验证数据中与总体平均生长过程有较大差
异的 36，39，43，45，46，53 号解析木，利用 SAS 9. 3
中的 IML 过程进行 EBLUP 预测，并与总体平均值进
行对比(图 1)。由图 1 可见，单株生长过程与总体
平均生长过程存在较大差异，并随着年龄的增加而
增加。除了 45 号解析木 (图 1 中的 plot45 ) 外，
Logistic 混合效应模型的 EBLUP 预测值可以充分预
测单株生长过程。45 号样木除了在幼龄阶段表现
出曲线特征外，已观测到的生长过程似乎用线性模
型描述更为精确，故 EBLUP 预测结果不佳。
2. 3. 2 EBLUP 预测的经验性分析 以验证数据的
30 株解析木为基础，经验性地分析 EBLUP 的特点。
用于 评 价 预 测 质 量 的 指 标 为 预 测 误 差 均 方
(MSPE)，其定义可参考式(11):
MSPE = 1
Σ
m
i = 1
ni
Σ
m
i = 1
Σ
n i
j = 1
( yij － y^ ij)
2。 (11)
观测间隔、观测次数和预测时长 3 个因素对
MSPE 的影响见表 7。观测间隔指每隔多少年进行
1 次树高生长观测，而观测次数则指进行了几次树
高观测，MSPE 列的数字代表起始观测年龄，亦即从
何年龄开始对树高进行观测。基本分析方法为，将
3 个因素中的 2 个保持不变，以比较出第 3 个因素
的影响作用。根据验证数据的年龄分布特征，用
30 ～ 100 年间的数据进行了预测验证。以第 4 行的
MSPE 值 1. 568 7 为例，以林龄为 50，55，60 时树高
观测值为基础，利用式(4)预测 b i 值。在获得了 b i
的预测值后，对 30 ～ 100 年间的树高值逐年进行预
测并求总 MSPE 值，即为 1. 568 7。
随着观测次数的增加，MSPE 表现出减少的趋
势(表 7)，这一特征同样与一般回归分析相类似。
从第 3 和 5 行、第 4 和 7 行以及第 5 和 8 行的逐对
比较可见，虽然观测间隔相同，但随着观测次数从 3
次增加到 5 次，MSPE 值均有 50% ～ 80% 以上的降
幅。在仅有一个观测值时，MSPE 值非常大，与
ADA，GADA 模型相比(赵磊等，2012)，并无明显的
优越之处。表 7 中的第 1，2，3 和 6 行的预测间隔均
为 2 年，但随着观测次数从 1 次增加到 5 次，MSPE
值显著降低。以起始观测年龄 30 年为例，MSPE 从
16. 694 2 降低至 4. 082 0。
在观测次数相同时，观测间隔增加同样可以降
低 MSPE 值。以第 6，7 和 8 行为例，均有 5 个观测，
但观测间隔不同，依次分别为 2，5 和 10 年，其 MSPE
值随着观测间隔的延长显著下降，此结论同样可从
表 7 的第 3 至第 5 行得出。在观测次数相同时，增
加观测值间年龄间隔，可以提供更为充分的树高生
长过程信息，因而可以显著降低 MSPE 值。
在观测间隔和观测次数相同的前提下，表 7 中
不同起始观测年龄间的 MSPE 显著不同，这主要由
预测时长的变化引起的，不足以归结为起始观测年
龄对 MPSE 的影响。
图 2 给出了表 7 中第 4 行的 MSPE 值(以林龄
为 50，55，60 时观测值求解 b i 值，并预测 30 ～ 100 年
树高值)，在观测值临近拟合数值时的预测非常精
确。在 35 ～ 70 年林龄的范围内，MSPE 值不超过
1. 0。随着预测偏离用以预测 b i 值的观测值，MSPE
值逐渐上升。
混合效应模型的 EBLUP 预测可以理解两阶段
拟合的第二阶段。第一阶段为混合效应模型的拟
合，用于估计总体分布特征，如随机效应的分布特征
和确定效应参数值;而第二阶段则在第一阶段拟合
的基础上，根据一个特定林分的已知观测值，对此林
分单独进行拟合，亦即以此林分的观测值为基础，预
测此林分特定的 b i 值，继而获得该林分生长过程模
型的参数值并进行预测。将预测 b i 值视为回归模
型的拟合，有助于理解 EBLUP 预测的特点。在一般
的回归分析中，MSPE 将随着预测点远离拟合数据
而表现出逐渐增加的趋势。与此相对应，EBLUP 预
测同样表现出在观测值附近的一定区域内，预测结
果非常精确，但随着预测时长增加，MSPE 值呈逐渐
增加的趋势。
03
第 3 期 祖笑锋等: 基于混合效应模型及 EBLUP 预测美国黄松林分优势木树高生长过程
图 1 预测单株解析木生长曲线
Fig． 1 The chart of predicting individual sample growth
表 7 验证数据统计分析
Tab． 7 The statistical analysis of data verification
序号
No．
观测间隔
Interval
观测次数
Observations
不同起始观测年龄时 MSPE
MSPE of observing different starting ages /m2
30 40 50 60
1 2 1 16. 694 2 10. 727 4 7. 255 4 6. 477 1
2 2 2 13. 082 3 8. 913 6 5. 882 2 5. 759 3
3 2 3 8. 981 8 6. 106 7 3. 878 2 4. 218 3
4 5 3 5. 583 2 3. 122 8 1. 568 7 1. 782 1
5 10 3 4. 092 3 1. 970 6 1. 008 9 0. 754 9
6 2 5 4. 082 0 3. 195 0 1. 263 6 1. 963 7
7 5 5 2. 894 8 1. 396 1 0. 876 7 0. 867 0
8 10 5 1. 185 4 0. 420 7 0. 221 8 0. 493 7
13
林 业 科 学 51 卷
图 2 30 ～ 100 逐年预测的 MSPE 值
Fig． 2 The value of 30 － 100 yearly forecast MSPE
3 结论与讨论
本文基于 3 种树高生长模型，在均衡考虑了统
计学和生物学特征的基础上，建立三参数 Logistic 混
合效应模型。在有多个以往观测数据的前提下，用
EBLUP 可以大大提高预测精度。对 EBLUP 进行深
入分析，得出以下几个结论:
1) 在有多个观测值且观测值时间间隔适宜的
情况下，EBLUP 可以充分预测树高生长过程。
2) 观测次数对预测精度的变化有很大影响。
随着观测次数的增加，预测误差表现出减少的趋势。
本研究中随着观测次数从 3 次增加到 5 次，预测误
差(MSPE)值均有 50% ～ 80% 以上的降幅，所以观
测次数越多，预测越准确。
3) 观测间隔不同，预测精度也有很大变化。观
测间隔的增加同样可以降低预测误差。在观测次数
相同的前提下，加大观测间隔可以大幅度降低预测
误差，提高预测精度。
4) 减小预测时长可提高 EBLUP 预测精度。
EBLUP 预测表现出在观测值附近的一定区域内，预
测结果非常精确，但随着预测时长的增加，预测误差
呈逐渐增加的趋势。
混合效应模型可以指定任意多个参数为随机效
应参数，此特征优于 ADA，GADA 模型。一般而言，
GADA 模型参数间必须假定存在特定的关系，以便
于用韦达定理求解。从已发表的研究成果看，一般
至多可以指定 2 个随机效应参数(Cieszewski et al．，
2000a; 赵磊等，2012)。
另一方面，从预测误差 MSPE 角度看，在仅存在
一个观测值的前提下，EBLUP 与 ADA /GADA 模型
相比，虽无定论但一般认为并无明显的优越之处
( Jiang et al．，2010; Meng et al．，2009b; Yang et al．，
2011);但在有多个观测值时，EBLUP 可以充分地利
用观测所涵盖的生长过程信息，极大地提高预测准
确性。与此相反，ADA，GADA 模型却仅能利用一个
观测值，即使有多个观测值也仅能选用一个而已。
参 考 文 献
李春明，张会儒 ． 2010． 利用非线性混合模型模拟杉木林优势木平
均高 ．林业科学，46(3) : 89 － 95．
(Li C M，Zhang H Ｒ． 2010． Modeling dominant height for Chinese fir
plantation using a nonlinear mixed-effects modeling approach．
Scientia Silvae Sinicae，46(3) : 89 － 95． ［in Chinese］)
李永慈，唐守正 ． 2004． 用 Mixed 和 Nlmixed 过程建立混合生长模型 ．
林业科学研究，17(3) : 279 － 283．
(Li Y C，Tang S Z． Establishment of tree height growth model based on
mixed and nlmixed of SAS． Forest Ｒesearch，17 (3) : 279 － 283．
［in Chinese］)
吕 萍，朱 钰 ． 2009． 基于最佳线性无偏估计的模型权数的小域
估计 ．统计与决策，(1) : 9 － 11．
(Lü P，Zhu Y． 2009． Estimate based on the number of best linear
unbiased estimate of the model right small domains． Statistics ＆
Decision，(1) : 9 － 11． ［in Chinese］)
孟宪宇 ． 1994． 测树学 ． 2 版 ．北京:中国林业出版社 ．
( Meng X Y． 1944． Forest mensuration． 2 nd edition． Beijing: China
Forestry Publishing House． ［in Chinese］)
赵 磊，倪成才，Gordon Nigh． 2012． 加拿大哥伦比亚省美国黄松广
义代数分型地位指数模型 ．林业科学，48(3) : 74 － 81．
(Zhao L，Ni C C，Nigh G． 2012． Generalized algebraic difference site
index model for ponderosa pine in British Columbia，Canada．
Scientia Silvae Sinicae，48(3) : 74 － 81． ［in Chinese］)
Bailey Ｒ L， Clutter J L． 1974． Base-age invariant polymorphic
sitecurves． Forest Science，20(2) : 155 － 159．
Calegario N，Daniels Ｒ F，Maestri Ｒ，et al． 2005． Modeling dominant
height growth based on nonlinear mixed-effects model: a clonal
eucalyptus plantation case study． Forest Ecology and Management，
204(1) : 11 － 21．
Cieszewski C J，Bailey Ｒ L． 2000a． Generalized algebraic difference
approach: theory based derivation of dynamic site equations with
polymorphism and variable asymptotes． Forest Science，46 ( 1 ) :
116 － 126．
Cieszewski C J， Harrison M，Martin S T． 2000b． Examples of
practicalmethods for unbiased parameter estimation in self-
referencingfunction ∥ Cieszewski C J． Proceedings of the first
international conferenceon measurements and quantitative methods
management． Jekyll Island，Georgia，November 17 － 18．
Clutter J L． 1963． Compatible growth and yield models for loblolly pine．
Forest Science，9(3) : 354 － 371．
Davidian M， Giltinan D M． 1995． Nonlinear models for repeated
measurement data． Chapman and Hall，London，62: 359．
Fang Z，Bailey Ｒ L． 2001． Nonlinear mixed effects modeling for slash
pine Dominant height growth following intensive silvi cultural
treatments． Forest Science，47(3) : 287 － 300．
Hall D B，Bailey Ｒ L． 2001． Modeling and prediction of forest growth
23
第 3 期 祖笑锋等: 基于混合效应模型及 EBLUP 预测美国黄松林分优势木树高生长过程
variables based on multilevel nonlinear mixed models． For Sci，
47(3) :311 － 321．
Hall D B，Clutter M． 2004． Multivariate multilevel nonlinear mixed
effects models for timber yield predictions． Biometrics，60 ( 1 ) :
16 － 24．
Huang S，Meng S X，Yang Y． 2009a． Estimating a non-linear mixed
volume-age model with and without taking into account serially-
correlated errors: differences and implications． Modern Appl Sci，
3(5) :3 － 20．
Huang S，Wiens D P，Yang Y，et al． 2009b． Assessing the impacts of
species composition，top height and density on individual tree height
prediction of quaking aspen in boreal mixedwoods． For Ecol
Manage，258(7) :1235 － 1247．
Jiang L，Li Y． 2010． Application of nonlinear mixed-effects modeling
approach in tree height prediction． J Computers， 5(10) :
1575 － 1580．
Lappi J，Bailey Ｒ L． 1988． A height prediction model with random stand
and tree parameter: an alternative to traditional site methods． Forest
Science，34(4) : 907 － 927．
Littell Ｒ C，Milliken G A，Stroup W W，et al． 1966． SAS system for
mixed models． SAS Institute Inc． Cary，NC，633．
Meng S X，Huang S． 2009． Improved calibration of nonlinear mixed-
effects models demonstrated on a height growth function． For Sci，
55(3) :238 － 248．
Meng S X，Huang S，Yang Y，et al． 2009． Evaluation of population-
averaged and subject-specific approaches for modeling the dominant
or codominant height of lodgepole pine trees． Can J For Ｒes，
39(6) : 1148 － 1158．
Ni C C，Nigh G． 2011． An analysis and comparison of predictors of
random parameters demonstrated on planted loblolly pine diameter
growth prediction． Forestry，85(2) : 271 － 280．
Nigh G． 2004． A comparison of fitting techniques for ponderosa pine
height-age models in British Columbia． Ann For Sci，61(7) : 609 －
615．
Pinheiro J C，Bales D M． 2000． Mixed effects models in S and S-plus
Springer Verlag． New York．
Tang S Z，Meng F Ｒ． 2001． Analyzing parameters of growth and yield
models for Chinese Fir —145． provenances with a inearmixed
approach． Silvae Genetica，50(3 /4) : 140 － 145．
Yang Y，Huang S． 2011． Comparison of different methods for fitting
nonlinear mixed forest models and for making predictions． Can J For
Ｒes，41(8) : 1671 － 1686．
(责任编辑 石红青)
33