留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于混合效应模型及EBLUP预测杉木树高生长过程

王明初 孙玉军

王明初, 孙玉军. 基于混合效应模型及EBLUP预测杉木树高生长过程[J]. 浙江农林大学学报, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003
引用本文: 王明初, 孙玉军. 基于混合效应模型及EBLUP预测杉木树高生长过程[J]. 浙江农林大学学报, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003
WANG Mingchu, SUN Yujun. Based on mixed-effects model and empirical best linear unbiased predictor predicting growth profile of height for Chinese fir[J]. Journal of Zhejiang A&F University, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003
Citation: WANG Mingchu, SUN Yujun. Based on mixed-effects model and empirical best linear unbiased predictor predicting growth profile of height for Chinese fir[J]. Journal of Zhejiang A&F University, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003

基于混合效应模型及EBLUP预测杉木树高生长过程

doi: 10.11833/j.issn.2095-0756.2017.05.003
基金项目: 

基于FORPLAN的森林多功能经营技术引进项目 2015-04-31

国家林业局重点资助项目 2012-07

林业科技成果国家级推广项目 [2014]26

详细信息
    作者简介: 王明初, 从事林业资源调查与监测研究.E-mail:wangmc@bjfu.edu.cn
    通信作者: 孙玉军, 教授, 博士, 博士生导师, 从事森林资源监测与模型等研究.E-mail:sunyj@bjfu.edu.cn
  • 中图分类号: S797.27

Based on mixed-effects model and empirical best linear unbiased predictor predicting growth profile of height for Chinese fir

图(2) / 表(9)
计量
  • 文章访问数:  1131
  • HTML全文浏览量:  123
  • PDF下载量:  358
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-12-16
  • 修回日期:  2017-02-18
  • 刊出日期:  2017-10-20

基于混合效应模型及EBLUP预测杉木树高生长过程

doi: 10.11833/j.issn.2095-0756.2017.05.003
    基金项目:

    基于FORPLAN的森林多功能经营技术引进项目 2015-04-31

    国家林业局重点资助项目 2012-07

    林业科技成果国家级推广项目 [2014]26

    作者简介:

    王明初, 从事林业资源调查与监测研究.E-mail:wangmc@bjfu.edu.cn

    通信作者: 孙玉军, 教授, 博士, 博士生导师, 从事森林资源监测与模型等研究.E-mail:sunyj@bjfu.edu.cn
  • 中图分类号: S797.27

摘要: 基于福建省将乐县国有林场15块标准地的30株杉木Cunninghamia lanceolata标准木的解析数据,首先对5个生长方程运用非线性最小二乘法进行拟合,选出拟合效果最好的模型作为基础模型,利用解析木数据构建非线性混合效应树高生长模型。以单株树木作为随机效应,通过变换混合效应参数个数,利用R软件选择赤池信息准则(AIC),贝叶斯信息准则(BIC)最小,对数似然函数(Loglik)值最大的混合效应模型作为最优模型,基于混合效应模型研究经验线性无偏最优预测法(EBLUP)预测树高生长过程的特点。结果表明:Weibull方程中,β1β2β3等3个参数都作为混合效应参数的模型模拟精度最高。观测次数相同时,延长观测间隔能够降低预测误差,提高预测精度;观测间隔相同时,增加观测次数,预测精度会提高。

English Abstract

王明初, 孙玉军. 基于混合效应模型及EBLUP预测杉木树高生长过程[J]. 浙江农林大学学报, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003
引用本文: 王明初, 孙玉军. 基于混合效应模型及EBLUP预测杉木树高生长过程[J]. 浙江农林大学学报, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003
WANG Mingchu, SUN Yujun. Based on mixed-effects model and empirical best linear unbiased predictor predicting growth profile of height for Chinese fir[J]. Journal of Zhejiang A&F University, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003
Citation: WANG Mingchu, SUN Yujun. Based on mixed-effects model and empirical best linear unbiased predictor predicting growth profile of height for Chinese fir[J]. Journal of Zhejiang A&F University, 2017, 34(5): 782-790. doi: 10.11833/j.issn.2095-0756.2017.05.003
  • 林木树高是基本测树因子之一,能够反映林木生长状况,是森林资源经营管理和林木生长收获研究所必需。树高生长模型是描述树高生长过程的统计模型,也是林木生长收获模型系统中一个重要的组成部分,同时林分优势木树高的预测为适地适树、森林经营活动和林分生长收获预测提供重要的基础数据[1]。因此,树高生长模型的研究对于建立林分生长模型系统及评价立地质量具有重要意义。模型的建立方法很多,混合效应模型是近代发展起来的新统计方法,混合效应模型由固定效应和随机效应两部分组成,既可以反映总体变化趋势,又可以提供方差、协方差等多种信息来反映个体之间的差异[2]。国外开展了大量树高生长的混合效应模型研究[3-12],国内对混合效应模型的拟合方法也进行了系统的研究[13-15],但是目前国内林业上对混合效应模型的应用基本限于与传统模型拟合效果上的比较,对预测的研究还很少。祖笑锋等[16]基于混合效应模型及经验线性无偏最优预测法(EBLUP)预测美国黄松Pinus ponderosa林分优势木树高生长过程,很好地解释了EBLUP,并且深入地分析了EBLUP的特点。EBLUP最早用于动物育种学中,与地统计学的Kriging,时间序列的卡尔曼滤波及小域估计法的数学原理一致[16]。本研究以福建将乐林场杉木Cunninghamia lanceolata人工林为研究对象,基于30株解析木数据,建立拟合最优的非线性混合效应树高生长模型,研究如何利用EBLUP预测树高生长过程,并通过设置不同的观测次数、观测间隔研究与预测精度之间的关系,并探究提高预测精度的方法。

    • 研究区位于福建省将乐县国有林场,26°26′~27°04′N,117°05′~117°40′E。将乐县位于福建省西北部,地处武夷山脉东南部,以中、低山为主,最高峰海拔为1 640.2 m。属亚热带季风气候,具有海洋性和大陆性气候特点,年平均气温为19.8 ℃,年平均降水量为2 027.0 mm。境内气温较高,夏季时间长,冬季较温暖,霜冻较少,生长期长。土壤以红壤为主,并分布有黄红壤,土层深厚,土质较好的一般为沙壤土或轻壤土,水分充足,土壤肥沃。植被以亚热带植物区系为主,植被种类非常丰富。其中人工次生林的乔木主要有杉木、马尾松Pinus massoniana等;灌木主要有粗叶榕Ficus hirta,黄毛楤木Aralia decaisneana,檵木Loropetalum chinensis等;草本主要有乌毛蕨Blechnum orientale,乌蕨Stenoloma chusanum,铁线蕨Adiantum capillus-veneris等。

    • 2010-2012年,根据林分不同年龄和密度,以典型抽样原则设置了15块20 m × 30 m标准地。对标准地内的林木进行每木检尺,选取标准木2株·标准地-1,共30株,伐倒并进行树干解析,解析木编号为F01~F30。根据解析木内业表内插法获得各年龄对应的树高值。随机选择20株进行拟合,10株进行预测(表 1)。

      表 1  建模数据和检验数据概况

      Table 1.  Summary statistics of the fitting and calibration data

      项目拟合数据(20)验证数据(10)
      树高/m年龄/a树高/m年龄/a
      平均值14.32315.124
      最大值25.03322.338
      最小值6.664.17
      标准差5.47.77.512.5
          说明:表头括号中为随机选择的株数。
    • 选择Richards方程、Weibull方程、Korf方程以及Logistic方程的2种变化形式等5个方程作为基础模型,对树高生长进行模拟(表 2)。运用非线性最小二乘法对5个树高生长方程进行拟合,选出拟合效果最好的模型作为基础模型,构建非线性混合效应模型[17]

      表 2  树高生长方程

      Table 2.  The equations of height growth

      编号模型方程
      1Richards $H = {\beta _1}{\left[ {1 - {\text{exp}}\left( {{\beta _2}t} \right)} \right]^{{\beta _3}}}$
      2Weibull $H = {\beta _1}\left[ {1 - {\text{exp}}\left( {{\beta _2}{t^{{\beta _3}}}} \right)} \right]$
      3Korf $H = {\beta _1}{\text{exp}}\left( {\frac{{ - {\beta ^2}}}{{{t^{{\beta _3}}}}}} \right)$
      4Logistic $H = \frac{{{\beta _1}}}{{1 + {\text{exp}}\left( {{\beta ^2} - {\beta _3}t} \right)}}$
      5Logistic2 $H = \frac{{{\beta _1}}}{{{\text{exp}}\left( {{\beta ^2} + {\beta _3}{\text{lg}}\left( t \right)} \right)}}$
    • 非线性混合效应模型主要特点是参数分为固定效应和随机效应,固定效应可以反映研究对象的总体变化规律,随机效应随个体的不同而变化,反映总体中不同个体的变化规律,从而获得较好的拟合效果。笔者主要研究利用EBLUP预测树高生长过程,并探究提高预测精度的方法。因此将样木作为随机效应,建立混合效应模型后,通过选取生长过程与总体平均生长过程差异较大的样木进行预测分析。非线性混合模型的一般表达式为:

      ${y_i} = f\left( {{A_i}\beta + {B_i}{b_i},{X_i}} \right) + {e_i}。$

      (1)

      式(1)中:βp×1维的固定效应向量;p为固定参数个数,biq×1维的随机效应向量;q为随机参数个数;AiBi为具有相应维度的设计矩阵,其元素通常为0,1或与固定效应和随机效应相关的协方差值;ei为随机误差项向量;f为非线性函数;Xi为年龄;yi为树高。其中随机向量biN(0, D),eiN(1, Ri),并且相互独立;eibi协方差矩阵分别为RiD,且为ni维和q维的对称矩阵。

      构建非线性混合效应模型还要确定以下3种结构:① 混合效应参数。依据PINHEIRO等[18]的研究,将基础模型中参数进行组合作为混合参数进行模拟,选择收敛的模型,通过比较其赤池信息量准则(AIC),贝叶斯信息准则(BIC)以及对数似然函数(Loglik)模型评价指标。AIC,BIC值越小,Loglik值越大,模型拟合效果越好。② 随机效应内的方差协方差矩阵。为了确定随机效应内的方差协方差矩阵,需要解决异方差和自相关性结构。以指数函数(Exp),幂函数(Power)和常数加幂函数(ConstPower)3种结构形式,消除数据间的异方差。表达式为:

      ${R_i}\left( {\beta ,{b_i},p} \right) = {\sigma ^2}G_i^{0.5}{I_{ni}}G_i^{0.5};$

      (2)

      ${\text{VarExp}}\left( {{e_i}} \right) = {\sigma ^2}{\text{exp}}\left( {\beta t} \right);$

      (3)

      ${\text{VarPower}}\left( {{e_i}} \right) = {\sigma ^2}{t^\beta };$

      (4)

      ${\text{VarConstPower}}\left( {{e_i}} \right) = {\sigma ^2}\left( {{\beta _1} + {t^{{\beta _2}}}} \right)。$

      (5)

      式(2)~式(5)中:σ2为误差方差值,由模型残差方差值给出;Gini×ni维对角矩阵来解释方差异质性,对角元素为相应误差项的标准差;由于本研究数据间不存在自相关性,因此,Inini×ni维的单位矩阵。ββ1β2为参数,t为年龄[19-22]。从自变量为年龄的3种函数中,根据AIC,BIC,Loglik值以及似然比检验,确定残差方差模型。③ 随机效应间的协方差矩阵。随机效应间的协方差矩阵反映了随机效应之间的变化性。本研究包含3个随机参数,因此,矩阵为3 × 3维的方差协方差矩阵,结构例如下:

      $D = \left[ {\begin{array}{*{20}{l}} u \\ v \\ w \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {\sigma _u^2}&{{\sigma _{uv}}}&{{\sigma _{uw}}} \\ {{\sigma _{vu}}}&{\sigma _v^2}&{{\sigma _{vw}}} \\ {{\sigma _{wu}}}&{{\sigma _{wv}}}&{\sigma _w^2} \end{array}} \right]。$

      (6)

      式(6)中:σu2σv2σw2分别为随机参数uvw的方差,σuv=σvu为随机参数uv的协方差,σuw=σwu为随机参数uw的协方差,σvw=σwv为随机参数wv的协方差。

    • 推导方法基于多元正态分布理论,以bi=0为基点,用一阶泰勒公式将非线性混合效应模型式(1)近似地表达为:

      ${y_i} \approx f\left( {{A_i}\beta ,{X_i}} \right) + {Z_i}{b_i} + {e_i}。$

      (7)

      式(7)中:矩阵Ziq × ni维矩阵,其元素可以通过对混合效应模型式(1)分别求关于随机效应bi的偏导数,随后令所有bi取0即可求出。前文提到biN(0, D),eiN(1, Ri),且相互独立,因此对式(7)求yi的方差矩阵为:

      ${\text{Var}}\left( {{y_i}} \right) = {Z_i}DZ_i^T + {R_i}。$

      biyi间协方差矩阵为:

      $\begin{gathered} {\text{Cov}}\left( {{b_i},{y_i}} \right) = {\text{Cov}}\left[ {{b_i},f\left( {{A_i}\beta ,{X_i}} \right) + {Z_i}{b_i} + {e_i}} \right] \hfill \\ = {\text{Cov}}\left( {{b_i},{Z_i}{b_i}} \right) = DZ_i^T。 \hfill \\ \end{gathered} $

      biyi的联合分布为:

      $\left( {\begin{array}{*{20}{l}} {{b_i}} \\ {{y_i}} \end{array}} \right) \sim {\text{N}}\left[ {\left( {\begin{array}{*{20}{c}} 0 \\ {f\left( {{A_i}\beta ,{X_i}} \right)} \end{array}} \right),\left( {\begin{array}{*{20}{c}} D&{DZ_i^T} \\ {{Z_i}{D^T}}&{{Z_i}DZ_i^T + {R_i}} \end{array}} \right)} \right]。$

      根据多元正态分布理论,当 ${y_i} = {{\tilde y}_i}$ ( ${{\tilde y}_i}$ 为观测值)时,bi的条件分布数学期望为:

      $E\left( {{b_i}|{y_i} = {{\tilde y}_i}} \right) = E\left( {{b_i}} \right) + DZ_i^T{\left( {{Z_i}DZ_i^T + {R_i}} \right)^{ - 1}}\left[ {{{\tilde y}_i} - f\left( {{A_i}\beta ,{X_i}} \right)} \right]。$

      对于随机变量取值最优预测值是其数学期望值,所以随机变量biyi取观测值时可以用bi的条件分布数学期望表示:

      ${{\hat b}_i} = DZ_i^T{\left( {{Z_i}DZ_i^T + {R_i}} \right)^{ - 1}}\left[ {{{\tilde y}_i} - f\left( {{A_i}\beta ,{X_i}} \right)} \right]。$

      将式中βDRZif(Ai β, Xi)用相应的估计值取代后则:

      ${{\hat b}_i} = \hat D\hat Z_i^T{\left( {{{\hat Z}_i}\hat D\hat Z_i^T + {{\hat R}_i}} \right)^{ - 1}}\left[ {{{\tilde y}_i} - f\left( {{A_i}\hat \beta ,{X_i}} \right)} \right]。$

      (8)

      式(8)称为经验线性无偏最优预测法(EBLUP)[16]

      设定 ${y_i} = \left( {{array}{*{20}{l}} {{y_{i,1}}} \\ {{y_{i,2}}} {array}} \right)$ ,yi, 1为观测值,yi, 2为预测值。则yi, 1yi, 2的联合分布为:

      ${y_i} = \left( {\begin{array}{*{20}{l}} {{y_{i,1}}} \\ {{y_{i,2}}} \end{array}} \right) \sim N\left[ {\left( {\begin{array}{*{20}{c}} {f\left( {{A_i}\beta ,{X_{i,2}}} \right)} \\ {f\left( {{A_i}\beta ,{X_{i,1}}} \right)} \end{array}} \right),\left( {\begin{array}{*{20}{c}} {{W_{22}}}&{{W_{21}}} \\ {{W_{12}}}&{{W_{11}}} \end{array}} \right)} \right]。$

      其中:W21=Cov(yi, 2, yi, 1)=Zi, 2DZi, 1T+Cov(ei, 1T, ei, 2),W11=Var(yi)=Zi, 1DZTi, 1+Ri, 1,则yi, 2的条件分布数学期望:

      ${{\tilde y}_{i,2}} = f\left( {{A_i}\tilde \beta ,{X_{i,2}}} \right) + {{\tilde W}_{21}}\tilde W_{11}^{ - 1}\left[ {{y_{i,1}} - f\left( {{A_i}\tilde \beta ,{X_{i,1}}} \right)} \right]。$

      将各项参数用相应估计值带入的:

      ${{\tilde y}_{i,2}} = f\left( {{A_i}\tilde \beta ,{X_{i,2}}} \right) + {{\hat Z}_{i,2}}\underbrace {\hat D\hat Z_{i,1}^T{{\left( {{{\hat Z}_{i,1}}\hat D\hat Z_{i,1}^T} \right)}^{ - 1}}\left[ {{y_{i,1}} - f\left( {{A_i}\hat \beta ,{X_{i,1}}} \right)} \right]}_{{{\hat b}_i}}。$

      简化为:

      ${{\tilde y}_{i,2}} = f\left( {{A_i}\hat \beta ,{X_{i,2}}} \right) + {{\hat Z}_{i,2}}{{\hat b}_i}。$

      (9)

      结合公式(8)和公式(9)即可利用混合效应模型进行EBLUP预测。详细推导过程参考文献[23]。

      根据EBLUP预测原理,式(9)中Zi, k可通过求关于随机效应bi的偏导数,求出偏导数后令所有bi取0即可。本模型有3个随机效应,故k取1,2,3。Zi, k的估计值为:

      ${{\hat Z}_{i,1}} = 1 - {{\text{e}}^{ - {{\tilde \beta }_2}{t_i}^{{{\tilde \beta }_3}}}};$

      (10)

      ${{\hat Z}_{i,2}} = {{\hat \beta }_1}{{\text{e}}^{ - {{\tilde \beta }_2}{t_i}^{{{\tilde \beta }_3}}}}{t_i}^{{{\tilde \beta }_3}};$

      (11)

      ${{\hat Z}_{i,3}} = {{\hat \beta }_1}{{\text{e}}^{ - {{\tilde \beta }_2}{t_i}^{{{\tilde \beta }_3}}}}{{\tilde \beta }_2}{t_i}^{{{\tilde \beta }_3}}{\text{lg}}{t_i}。$

      (12)

      列向量 $f\left( {{A_i}\hat \beta ,{X_i}} \right)$ 表达式为:

      $f\left( {{A_i}\hat \beta ,{X_i}} \right) = {{\hat \beta }_1}\left( {1 - {{\text{e}}^{ - {{\tilde \beta }_2}{t_i}^{{{\tilde \beta }_3}}}}} \right)。$

      (13)

      式(11)~式(14)中: ${{\hat \beta }_{i,1}},{{\hat \beta }_{i,2}},{{\hat \beta }_{i,3}}$ 为βi, 1βi, 2βi, 3的估计值,通过R语言软件得出。随机效应参数的预测值 ${{\hat b}_{i,1}},{{\hat b}_{i,2}},{{\hat b}_{i,3}}$ 根据式(8),利用SAS 9.3的IML过程计算获得。将各参数估计值代入式(9)进行预测。

    • 利用R语言软件基于表 1中拟合数据,对表 2中6个树高生长方程进行拟合,拟合结果见表 3。根据决定系数R2越接近于1,均方根差Rmse以及平均绝对残差|E|越接近于0,其拟合效果越好的原则,选取方程1,方程2,方程3作为基础模型。

      表 3  各树高生长方程拟合结果

      Table 3.  Simulant result of height growth

      方程参数确定系数R2均方根差(Rmse平均绝对残差|E|
      β1β2β3
      145.517 70.024 11.192 10.819 92.698 80.042 7
      241.994 90.013 21.157 90.819 92.699 10.076 2
      32 647.267 08.858 00.177 10.820 12.697 20.002 7
      420.734 92.434 60.169 20.808 32.784 42.976 8
      565.986 04.821 0-2.730 00.819 92.986 00.042 1
    • 基于基础模型,将单木考虑作为随机效应,利用R语言软件中的nlme函数对模型进行估计,改变随机参数个数比较拟合效果。基于不同随机参数组合的混合模型的拟合效果见表 4。从表中拟合情况可以看出模型2.5的AIC,BIC最小且Loglik值最大,即当模型2以β1β2β3同时作为混合效应参数时,模型拟合效果最佳。因此,选择模型2.5(Weibull方程)进一步构建混合效应模型。

      表 4  基于不同随机效应参数组合的树高生长模型拟合精度比较

      Table 4.  Comparison of models' fitting precisions with different combinations of random effects parameters

      模型编号随机效应参数AICBICLoglik
      1.1β11 516.717 01 537.274-753.358 4
      1.2β21 461.503 01 482.060-725.751 3
      1.3β31 600.736 01 621.293-795.368 1
      1.4β1, β31 207.652 01 236.432-596.826 0
      1.5β2, β31 147.023 01 175.803-566.511 4
      2.1β11 516.378 51 536.936-753.189 3
      2.2β21 455.967 71 476.525-722.983 8
      2.3β31 382.081 01 402.638-686.040 5
      2.4β2, β31 158.258 61 187.039-572.129 3
      2.5β1, β2, β3973.412 21 014.527-476.706 1
      3.1β21 513.753 01 534.311-751.876 7
      3.2β31 484.579 01 505.136-737.289 6
      3.3β1, β21 253.315 01 282.096-619.657 7
      3.4β2, β31 241.380 01 270.160-613.690 0
      其他情况不收敛
    • 加权残差方差模型的评价指标见表 5。考虑残差加权的3种模型中幂函数的AIC和BIC最小且Loglik值最大,故选它作为模型的残差方差模型。

      表 5  各残差方差模型的模拟结果

      Table 5.  Simulation results for each residual variance model

      残差方差模型AICBICLoglik
      指数函数975.447 11 020.673-476.723 5
      幂函数975.070 01 020.296-476.535 0
      常数加幂函数977.447 01 026.785-476.723 5
    • 利用R语言软件拟合模型的各参数列于表 6,基于非线性混合效应模型的杉木树高生长方程表达式为:

      表 6  模型拟合结果

      Table 6.  Statistical results of the model fitting

      固定效应参数随机效应参数协方差矩阵
      参数估计值参数bi, 1bi, 2bi, 3
      β11.954 00bi, 15.936 40.016 60.473 7
      β20.018 96bi, 20.016 6-0.302 0-0.251 0
      β31.631 90bi, 30.473 7-0.251 0-0.655 0

      $\begin{gathered} H = 1.954\;0\left( {1 - {{\text{e}}^{ - 0.018\;96{t^{1.6319}}}}} \right); \hfill \\ {b_i} = \left[ {\begin{array}{*{20}{c}} {{u_i}} \\ {{v_i}} \\ {{w_i}} \end{array}} \right] \sim N\left( {0,D} \right){e_i} \sim N\left[ {0,{R_i}} \right]; \hfill \\ R = 0.535\;2G_i^{0.05}{\text{Var}}\left( {{e_i}} \right) = 0.535\;2{t^{ - 0.004\;446}}。 \hfill \\ \end{gathered} $

      (14)
    • 由于EBLUP算法的原因,本研究利用式(9)预测,具体原因参考NI等[23]。从验证数据中选取F04解析木作为预测数据,选3对观测值用于预测bi,结果统计于表 7。表中 ${{\tilde y}_i}$ 为样木在Xi时树高观测值, $f\left( {{A_i}\hat \beta ,{X_i}} \right)$ 通过式(13)计算得出, ${{\tilde Z}_i}$ 通过式(10)~式(12)计算得出, ${{\tilde b}_i}$ 是通过SAS 9.3的IML过程获得。利用表 7的 ${{\tilde b}_i}$ 预测F18样木19~22 a间的树高值,并与观测值比较列于表 8。总体估计值通过式(13)计算得到,预测值通过式(9)计算获得,校正值为预测值减总体估计值的差。可见根据3对观测值获得的对样木进行EBLUP预测树高,能较准确地反映样木生长过程。在验证数据中选择样木生长过程与总体平均生长过程有差异的F04,F09,F11,F18解析木,进行EBLUP预测,将预测生长曲线、解析木树高曲线和总体平均生长过程进行比较(图 1)。通过图 1可看出:当单木生长过程和总体估计值差距较大时,基于Weibull混合效应模型的EBLUP预测能充分地预测单木生长过程。本研究依据EBLUP预测的算法,利用样木的3组观测值,基于混合效应模型使每个样木获得对应的随机效应参数,进而充分地预测样木的生长过程。

      表 7  预测的向量和矩阵

      Table 7.  The matrices and vectors used for prediction

      ${{\tilde y}_i}/{\text{m}}$ ${X_i}/{\text{a}}$ $f\left( {{A_i}\hat \beta ,{X_i}} \right)/{\text{m}}$ ${{\hat Z}_{i,1}}$ ${{\hat Z}_{i,2}}$ ${{\hat Z}_{i,3}}$ ${{\hat b}_i}$
      16.602318.713 90.957 7137.837 18.196 0-2.412 3
      15.001817.195 60.880 0262.154 114.369 5-0.016 4
      12.001313.923 50.712 6369.254 317.961 20.312 6

      表 8  F18样木树高观测值、总体估计值、EBLUP预测值比较

      Table 8.  Comparison of observations of height, population mean, and EBLUP predicted for F18

      年龄/a树高观测值/m总体估计值/mEBLUP预测值校正值
      2216.2818.509 216.533 9-1.975 3
      2115.9618.261 816.312 6-1.949 2
      2015.6417.965 116.0342-1.930 9
      1915.3217.612 115.689 3-1.922 8

      图  1  树高生长过程比较

      Figure 1.  Comparison of the height growth

    • 在进行单木树高生长预测时发现,观测值的选择方式对预测精度高低有较大的影响,因此本研究通过设置不同的观测值间隔、观测值选取次数以及预测时长,比较预测结果的精度高低。评价预测精度的指标为均方误差(EMS):

      ${E_{{\text{MS}}}} = \frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} 。$

      以解析木F04,F18,F15和F11为例(由于这几株样木生长过程与总体平均生长过程差距较大,通过不同处理方式能较明显反应预测精度的提高),通过不同的选择观测值方式来研究对EMS的影响,其中观测间隔指每隔多少年对树高进行一次观测。本研究分别选取间隔1,3,5 a;观测次数是指观测多少个树高值,由于本研究基于具有3个随机效应参数的混合模型,预测随机效应参数至少需要3对观测值,故分别选取3,6,9次观测次数;相应地,因为每3对观测值会求得1组随机效应参数,故观测次数为6,9次时,需要分别取2,3组数据,分别按照1/2,1/3将解析木数据分段,年龄向上取整,每年龄段从最大值开始向下取观测值。以样木F15(35年生)为例,观测间隔5 a观测3次时,使用的是35,30,25年生的树高观测值;观测间隔为3 a观测6次时,分别使用35,32,29年生以及18,15,12年生这2组共6个树高观测值;观测间隔为1 a观测9次时,分别使用35,34,33年生,24,23,22年生以及12,11,10年生这3组共9个树高观测值。由于F04,F18年龄较小,故没有进行间隔5 a 9次的观测。以树高观测值为基础,根据式(8)预测bi值,再对树高逐年预测最后求总均方误差。结果列于表 9

      表 9  验证数据统计分析

      Table 9.  The statistical analysis of data verification

      观测间隔/a观测次数/次均方误差
      F04(23年生)F18(28年生)F15(35年生)F11(38年生)
      135.546 04.178 712.396 022.533 3
      332.866 11.092 43.388 92.310 8
      532.746 00.391 72.704 30.941 7
      165.218 63.982 95.238 33.607 2
      360.332 10.21452.212 00E925 9
      560.271 60.195 60.284 00E707 1
      194.256 12.401 72.631 71E020 8
      390.16700.189 40.200 30.379 7
      590.094 60E247 0

      观测次数相同时,随着观测间隔的增大,均方误差降低。观测次数均为3,观测间隔取1,3,5 a,各解析木的均方误差随着观测间隔的增大而明显降低。观测间隔相同时,随着观测次数的增加,均方误差降低。观测间隔均为3 a时,随着观测次数从3次增加到9次,各组的均方误差均相应降低,观测间隔为1或5 a按同样方式比较,也可以得出相同的结论。以F15为例,观测间隔均为1 a时,随着观测次数从3次增加到9次,均方误差从12.396 0降低到2.631 7。

      图 2中可以看出:“1 a 3次”在其观测值28,27,26年生附近均方误差较低,远离观测值时,均方误差开始变大。同样,“3 a 3次”在远离20年生时均方误差开始变大,而“5 a 3次”的观测值为28,23,18年生,故在图 2中均方误差较低。所以在EBLUP预测时,随着预测远离用来预测bi的观测值,均方误差逐渐变大。

      图  2  解析木F18逐年预测的均方误差

      Figure 2.  The value of F18 yearly forecast EMS

    • 本研究选用5个生长方程,选取3个拟合效果较好的方程作为基础模型,进一步改变不同的随机效应组合,最终以式(2)为基础的Weibull方程中,3个参数都作为混合模型的模型精度最高,拟合效果最好。通过EBLUP预测表明:本研究建立的非线性混合效应模型在考虑样木随机效应时,预测效果明显优于固定模型,能更好地反映树高生长过程。将随机效应考虑到样木水平,就相当于为每株样木建立树高生长模型,可以充分预测样木的生长过程,能够为将乐林场的杉木林经营决策提供一定的依据。进行预测时,观测次数相同,增加观测间隔能够降低预测误差,提高预测精度。由于增加观测值间年龄间隔,可以更充分提供树高生长过程信息,因此可以显著降低预测误差。观测次数不同,预测精度变化较大。观测间隔相同时,增加观测次数,预测精度会提高。因此在预测杉木树高生长过程时,要根据其年龄选择适当的观测间隔和观测次数。EBLUP可以充分利用观测值所涵盖的生长过程信息降低预测误差,随着预测远离用来预测bi的观测值,均方误差逐渐变大。

      通过本次研究得到的数据不能明确观测间隔和观测次数对预测准确度影响,有关研究还需要进一步进行。

参考文献 (23)

目录

    /

    返回文章
    返回