基于线性混合模型的红松人工林一级枝条大小预测模拟

基于线性混合模型的红松人工林一级枝条大小预测模拟*
董灵波　刘兆刚**
　李凤日　姜立春
(东北林业大学林学院,哈尔滨150040)
摘　要　基于黑龙江省孟家岗林场60株人工红松955个标准枝数据,采用线性混合效应模型理论和方法,考虑树木效应,利用SAS 软件中的MIXED 模块拟合红松人工林一级枝条各因子(基径㊁枝长㊁着枝角度)的预测模型.结果表明:通过选择合适的随机参数和方差协方差结构能够提高模型的拟合精度;把相关性结构包括复合对称结构CS ㊁一阶自回归结构AR (1)及一阶自回归与滑动平均结构ARMA (1,1)加入到一级枝条大小最优混合模型中,AR (1)可显著提高枝条基径和角度混合模型的拟合精度,但3种结构均不能提高枝条角度混合模型的精度.为了描述混合模型构建过程中产生的异方差现象,把CF 1和CF 2函数加入到枝条混合模型中,CF 1函数显著提高了枝条角度混合模型的拟合效果,CF 2函数显著提高了枝条基径和长度混合模型拟合效果.模型检验结果表明:对于红松人工林一级枝条大小预测模型,混合效应模型的估计精度比传统回归模型估计精度明显提高.
关键词　红松人工林　枝条基径　枝条长度　枝条角度　线性混合模型
文章编号　1001-9332(2013)09-2447-10　中图分类号　S758.1　文献标识码　A
Primary branch size of Pinus koraiensis plantation :A prediction based on linear mixed effect model.DONG Ling⁃bo,LIU Zhao⁃gang,LI Feng⁃ri,JIANG Li⁃chun (College of Forestry ,North⁃east Forestry University ,Harbin 150040,China ).⁃Chin.J.Appl.Ecol .,2013,24(9):2447-2456.
Abstract :By using the branch analysis data of 955standard branches from 60sampled trees in 12sampling plots of Pinus koraiensis plantation in Mengjiagang Forest Farm in Heilongjiang Province of Northeast China,and based on the linear mixed⁃effect model theory and methods,the models for predicting branch variables,including primary branch diameter,length,and angle,were devel⁃oped.Considering tree effect,the MIXED module of SAS software was used to fit the prediction models.The results indicated that the fitting precision of the models could be improved by choosing appropriate random⁃effect parameters and variance⁃covariance structure.Then,the correlation structures including complex symmetry structure (CS ),first⁃order autoregressive structure [AR(1)],and first⁃order autoregressive and moving average structure [ARMA(1,1)]were added to the optimal branch size mixed⁃effect model.The AR(1)improved the fitting precision of branch diameter and length mixed⁃effect model significantly,but all the three structures didn’t improve the precision of branch angle mixed⁃effect model.In order to describe the heteroscedasticity during building mixed⁃effect model,the CF 1and CF 2functions were added to the branch mixed⁃effect model.CF 1function improved the fitting effect of branch angle mixed model significantly,whereas CF 2function improved the fitting effect of branch diameter and length mixed model significantly.Model validation confirmed that the mixed⁃effect model could improve the precision of prediction,as compare to the traditional regression model for the branch size prediction of Pinus koraiensis planta⁃tion.
Key words :Pinus koraiensis plantation;branch diameter;branch length;branch angle;linear mixed model.
*林业公益性行业科研专项(201004002007)㊁ 十二五”国家科技计划项目(2012BAD22B0202)㊁长江学者和创新团队发展计划项目(IRT1054)资助.
**通讯作者.E⁃mail:lzg19700602@ 2012⁃11⁃12收稿,2013⁃06⁃03接受.
应用生态学报　2013年9月　第24卷　第9期 Chinese Journal of Applied Ecology,Sep.2013,24(9):2447-2456
树冠是林木进行光合作用㊁呼吸作用等一系列生理活动的主要场所,其大小㊁结构及其分布形式直接决定树木的生长活力和生产力[1].一级侧枝作为树冠重要组成部分,是构成树的骨架,它的长短及其空间排列对树冠的形状起着支配作用,定量模拟其数量㊁大小㊁分布和生长可以很好地预测树冠动态㊁树木的生长和木材质量,有利于评价和选择合理的经营措施.
目前,一些学者对不同树种树冠内枝条大小特征进行了研究,但建立的枝条预测模型集中在采用传统线性和非线性回归模型上,以林分变量(立地㊁密度等)㊁林木变量(胸径㊁树高㊁冠长等)和着枝深度(DINC)为自变量来直接构建枝条各变量(如基径㊁枝长㊁着枝角度㊁数量等)的预测模型[1-6];采用枝解析方法,使同一林分不同样地㊁同一样地不同树木间具有相关性,同时也使同一树木各枝条之间具有相关性.因此,这种 样地效应”或 树木效应”使误差不满足独立同分布的假设.由于样地和树木是随机选择的,所以随机误差至少应包括样地间㊁树木间和树木个体(即枝条间)内3个随机效应部分[7].而混合模型通过规定不同的协方差结构来表示相关的误差情况,消除了层次结构数据产生的误差相关性和异质性,从而提高模型的预估精度并解释随机误差来源[8].此方法既能反映总体的平均变化趋势,又可以提供数据方差㊁协方差等多种信息来反映个体之间的差异.近年来,混合模型技术已广泛应用于单木胸径[9-10]㊁树高[11-12]㊁材积[13]㊁干形[14],以及林分的断面积[8]㊁树高曲线[15]㊁枯损[16]等模型的研究中,取得了很好的拟合效果.
国外已有学者利用混合模型方法对树木枝条特征进行研究.如Meredieu等[17]应用线性混合模型建立欧洲黑松(Pinus nigra)枝条长度和角度模型; Beaulieu等[18]采用线性混合效应模型方法,考虑了区域㊁样地和单木的随机效应,建立了北美短叶松(Pinus banksiana)枝条特征模型;Hein等[19]采用线性混合模型方法模拟了不同密度(1600㊁700和350株㊃hm-2)下挪威云杉(Picea abies)枝条预测模型.以上研究都表明:在模拟枝大小条特征时,混合模型方法的精度比传统回归方法高.迄今,国内应用混合模型研究枝条特征模型还鲜有报道,仅姜立春等[20]采用线性混合模型方法建立了落叶松(Larix gmeli⁃nii)枝条长度和角度模型,而且未包含枝条基径模型.由于不同树种枝条特征的变化既受自身遗传因素的制约,也受外界环境条件的影响,所以有必要开展不同树种枝条属性变化规律的研究.因此,本文应
用黑龙江省东北部地区孟家岗林场60株红松枝解
析数据,采用线性混合模型建立红松(Pinus koraien⁃sis)人工林一级枝条大小(基径㊁长度和着生角度)混合效应预测模型,利用验证数据进行验证,选择确
定系数(R2)㊁平均绝对误差(MAE)和均方根误差(RMSE)3个评价指标,对考虑单木效应㊁时间序列相关性及方差异质性的混合模型模拟结果与传统线
性回归模拟方法进行精度比较,以检验混合模型在
枝条特征预测模型中的应用.
1　研究地区与研究方法
1.1　研究区概况
研究区域位于黑龙江省佳木斯市孟家岗林场(46°20′30″ 46°30′50″N,130°32′42″ 130°52′36″E).该地区地处完达山西麓余脉,以低山丘陵为主,坡度较平缓,平均海拔250m;属东亚大陆性季风气候,年均气温2.7℃,年均降水量550mm,年全日照时数1955h,年无霜期120d;土壤种类以暗棕壤为主,另有少量白浆土㊁草甸土㊁沼泽土和泥炭土.该林场总经营面积16274hm2.其中,天然林面积3597 hm2,主要是以柞树(Quercus mongolica)㊁椴树(Tilia tuan)㊁白桦(Betula platyphylla)㊁山杨(Populus da⁃vidiana)为主的次生落叶阔叶混交林;人工林面积9482hm2,主要为红松㊁樟子松(Pinus sylvestris var. mongolica)㊁落叶松(Larix gmelinii)等,其中,红松人工林面积1835hm2,蓄积量1.83×105m3,分别占整个林场总面积和蓄积的11.3%和13.0%.
1.2　研究数据
为揭示红松人工林果实产量与林分因子㊁林木
因子和树冠结构的关系,结合研究区域红松人工林
资源现状,于2010年在研究区域具有代表性的红松
人工林中设置12块标准地,样地面积一般为0.06 hm2,最大0.09hm2,实测样地内每株林木的胸径㊁树高㊁冠幅㊁枝下高和林木坐标等(表1).每块标准地按照等断面积径级标准木法,将林木大小划分5个等级,将各等级林木的平均胸径㊁树高和冠长作为选取解析木的标准,并在每块标准地附近选伐样木.将解析木的树干按1m区分段进行区分,并在每个区分段的中央位置锯取树干解析圆盘.树冠部分也按1m区分段进行区分,树梢在一般情况下是不足1m的梢头.在每个区分段内对一级枝条进行逐枝编号,并测定每个枝条的总着枝深度(DINC)㊁枝条的基径(BD)㊁枝长(BL)㊁弦长(BCL)㊁弓高(BAH)㊁
8442应　用　生　态　学　报 24卷
表1　标准地林分调查因子概况
Table1　Summary of stand description attributes of sample plot
林分变量
Stand variable标准地数
Number of plot平均值
Mean最小值
Minimun最大值
Maximum标准差SD变异系数
CV(%)年龄Age(a)1237.583246 3.669.7
平均胸径Mean DBH(cm)1218.7715.1721.46 2.3912.7
平均树高Mean tree height(m)1212.069.6514.64 1.3511.2林分密度Stand density(plant㊃hm-2)121085.676501650349.5832.2林分蓄积量Stand volume(m3㊃hm-2)12137.17113.23180.6220.0014.6
表2　红松人工林样木和一级枝条因子概况
Table2　Summary of sample tree and primary branch attributes for Pinus koraiensis plantation
变量Variable
拟合数据Modeling data
平均值
Mean最小值
Minimum最大值
Maximum标准差SD
检验数据Validation data
平均值
Mean最小值
Minimum最大值
Maximum标准差SD
树木变量Tree variable(n=60)n=48n=12
胸径DBH(cm)19.8812.5027.30 3.7719.6512.3023.40 3.51树高Tree height(m)12.189.5016.10 1.5011.929.7014.60 1.22冠幅Crown width(m) 2.28 1.33 3.700.66 2.14 1.30 3.130.50枝条变量Total branch(n=955)n=772n=183
总着枝深度Depth into crown(m) 2.610.078.92 1.72 2.420.12 6.95 1.58方位角Azimuth angle(°)162.320350.00107.28147.600350.00102.95着枝角度Branch angle(°)64.6920.0090.0012.3566.8930.0090.0011.68基径Branch diameter(mm)23.95 1.4074.9111.1222.37 5.4253.0810.43枝长Branch length(cm)180.5111.00506.00109.40160.8211.00493.00103.95弦长Branch chord length(cm)166.1811.00491.0099.52144.469.00469.0095.10弓高Branch arc height(cm)22.460106.0017.6720.66090.0016.42
着枝角度(BA)和方位角(A)等枝条特征变量.按系统抽样选取48棵解析木的772个枝条数据用于建模,其他12棵解析木的183个枝条数据用于模型检验.红松人工林样木和枝条特征调查因子的统计量见表2.
1.3　研究方法
1.3.1基础模型　通过绘制散点图分析枝条基径㊁长度和着枝角度与树木因子和林分因子的关系.研究发现,枝条的基径㊁长度和着枝角度与枝条的总着枝深度密切相关,并随着总着枝深度的增加而持续增加;同时,其与林木胸径㊁树高㊁冠幅㊁冠长等变量组也存在一定的相关关系,但与林分因子没有明显的相关性.基于以上关系,本研究中模型自变量选择过程用方差膨胀因子(variance inflation factor,VIF)来判断自变量间的多重共线性.一般认为,当VIF>
10时,有严重的共线性,此时剔除共线性较严重的自变量,保留共线性弱而对因变量贡献大的自变量[9-10,21].只有回归系数显著(P<0.05)㊁且方差膨胀因子<10的变量才能进入模型.本文最终建立的枝条基径㊁长度和着枝角度预测模型如下:
ln(BD)=a0+a1ln(DINC)+a2ln(DBH)(1)
ln(BL)=a0+a1ln(DINC)+a2DINC2+
a3ln(DBH)(2) BA=a0+a1ln(DINC)+a2DBH+a3BD(3)式中:BD为枝条基径(mm);BL为枝条长度(cm); BA为枝条着枝角度(°);DBH为胸径(cm);DINC 为枝条的总着枝深度(m);a0㊁a1㊁a2㊁a3均为模型预估参数.
1.3.2线性混合模型　线性混合模型是在线性模型的基础上,通过引入随机效应而构建.在本研究中,单水平线性混合效应模型可表示为[22]:
y ijk=X xjkβ+Z ijkμ+εijk
εijk~N(0,R ij),μij~N(0,D)
μij=^D^Z ij T(^R ij+^Z ij^D^Z ij T)-1^εijk
i=1, ,m;j=1, ,n i;k=1, ,f
ì
î
í
ï
ïï
ï
ï
ij
(4)
式中:y ijk为一级枝条大小因子的观测值;m为样地数量;n i为样地j中解析木的数量;f ij为第i样地第j 株样木中标准枝的数量;X ijk为(f ijk×p)维已知设计矩阵;β为(p×1)维固定参数向量;Z ijk为(f ijk×q)维随机效应设计矩阵;μ为(q×1)维随机参数向量;D为(q×q)随机效应方差协方差矩阵;R ij为第i样地第j 株树木(k×k)维误差效应方差协方差结构;^D为预
9442
9期 董灵波等:基于线性混合模型的红松人工林一级枝条大小预测模拟
估的随机效应参数方差协方差结构矩阵;^R ij为预估的树木内方差协方差结构;^Z ij为设计矩阵;^εijk为(k×1)维固定效应模型的残差向量;εijk为(n ij×1)维混合效应模型的误差项.
在基本模型的基础上构建混合模型,包括固定参数㊁随机参数和方差协方差结构.混合模型的构建需要确定以下3个结构[11,23]:
1)确定参数效应.混合模型中参数效应的确定一般依赖于所研究的数据.Pinheiro和Bates[22]建议,如果模型能够收敛,首先应把模型中所有的参数都看成是混合参数,然后将不同随机参数组合的模型进行拟合,比较模型拟合的统计量,即比较AIC㊁BIC和⁃2log Likelihood指标,3个指标越小越好.为了避免过多参数化问题,具有不同混合参数个数的模型还要进行似然比检验(likelihood ratio test, LRT),即利用LRT和P值进行检验,如果P<0.05即认为差异显著.
2)树木内方差协方差结构(R ij).为了确定树木内方差协方差结构,必须解决异方差和自相关结构两方面的问题.由于枝解析数据通常都存在自相关和异方差问题,为了解决树木内误差的异质性和自相关问题,在统计和林业上基本上都采用下式进行描述[11,24-26]:
R ij=σ2G ij0.5Γij G ij0.5(5)式中:R ij为模型的方差协方差矩阵;σ2为模型的残差方差值;Γi j为描述误差效应自相关结构矩阵,即描述同一树木内不同测量值相关性;G i j为描述方差异质性的对角矩阵.
3)确定随机效应的方差协方差结构(D).树木间的方差协方差结构反映了树木之间的变化和差异.本研究检验了林业上常用的3种方差协方差结构:无结构矩阵模型(UN)㊁复合对称矩阵模型(CS)和对角矩阵模型(UN(1)).以包含3个随机参数(a1㊁a2㊁a3)的方差协方差结构为例,无结构矩阵模型为:
D=
σa
1
2σ
a1a2
σa
1a3
σa
1a2
σa
2
2σ
a2a3
σa
1a3
σa
2a3
σa
3
æ
è
ç
ç
çç
ö
ø
÷
÷
÷÷
2
(6)
复合对称矩阵模型为:
D=
σ12+σ2σ12σ12
σ12σ12+σ2σ12
σ12σ12σ12+σ
æ
è
ç
ç
ç
ö
ø
÷
÷
÷2
(7)
对角矩阵模型为:
D=
σa
1
200
0σa220
00σa3
æ
è
ç
ç
çç
ö
ø
÷
÷
÷÷
2
(8)
式中:σa
1
2为随机参数a
1的方差;σa2
2为随机参数a
2
的方差;σa32为随机参数a3的方差;σa1a2为随机参数
a1和a2的协方差;σa
1a2
为随机参数a1和a3的协方
差;σa
2a3
为随机参数a2和a3的协方差;σ2为模型的
方差.
1.3.3模型评价和检验　本研究采用的3个模型模
拟效果评价指标为AIC㊁BIC㊁⁃2Log likelihood.利用
确定系数(R2)㊁平均绝对误差(MAE)和均方根误差
(RMSE)3个模型精度评价指标对验证数据进行模
拟精度评价.其中,RMSE和MAE值越小㊁R2值越
大,说明模型的模拟精度越高.
R2=1-
∑m i=1∑n i j=1∑f ij k=1(y ij-^y ij)2
∑m i=1∑n i j=1∑f ij k=1(y ij-⎺y)
é
ë
ê
ê
ê
ù
û
ú
ú
ú2
(9)
MAE=∑m i=1∑n i j=1∑f ij k=1|y ij-^y ij|
é
ë
ê
ê
ù
û
ú
ú
n
(10)
RMSE=∑m i=1∑n i j=1∑f ij k=1(y ij-^y ij)2
n-
é
ë
ê
ê
ù
û
ú
ú
1
(11)
式中:y ij为观测值;^y ij为预测值;⎺y为观测值的平均
值.
2　结果与分析
2.1　基础模型的拟合结果
利用SAS9.2软件对式(1)~(3)3个模型进行
拟合,选择确定系数㊁平均绝对误差和均方根误差对
模拟结果进行比较分析.从表3可以看出,枝条基径
和长度模型的拟合效果较好,其模型的R2分别达到
0.64和0.89,并且模型MAE和RMSE值也较小;模
型(1)和(2)的参数t检验结果表明,所有参数估计
均极显著(P<0.01);模型检验(F=462.65㊁P<
0.0001,F=1969.25㊁P<0.0001)表明,模型(1)和
(2)对描述红松人工林一级枝条基径和长度的变化
具有显著意义.模型3的模拟效果相对较差,R2仅为
0.16,但模型(3)的参数t检验也极显著(P<0.01),
并且模型检验(F=48.41,P<0.0001)结果表明,该
模型对描述红松人工林一级枝条着枝角度的变化规
0542应　用　生　态　学　报 24卷
表3　红松人工林枝条大小基础模型的拟合结果
Table3　Simulation results of the primary branch size models for Pinus koraiensis plantation(n=772)
模型
Model参数估计值Parameter estimate
a0a1a2a3
拟合统计量Good⁃of⁃fit statistics
确定系数
R2
平均绝对误差
MAE均方根误差
RMSE
1 2.0650.5130.2140.6960.2210.295
2 3.780 1.009-0.0090.1880.8850.1960.284 361.9309.0410.429-0.5190.1598.77511.349
律具有一定意义.因此,本文选择上述模型作为构建混合模型的基础模型.
2.2　基于单木效应混合模型的随机参数确定
通常,确定参数效应和随机效应的方差协方差结构须同时进行,首先把随机效应的方差协方差结构设定为无结构矩阵,考虑树木效应,利用SAS软件的MIXED模块对上述模型的不同随机参数组合进行拟合.利用AIC㊁BIC㊁-2Log Likelihood及似然比检验等统计量指标对模型的拟合优度进行比较,3个指标越小,说明拟合的精度越高.为了避免过多参数化问题,具有不同参数个数的模型还采用似然比检验来比较,如P<0.05即认为差异显著.拟合结果见表4.限于篇幅,仅给出了具有相同参数个数中的最优模型,即AIC㊁BIC㊁-2Log Likelihood最小.
从表4可以看出:1)对于枝条基径模型,混合模型收敛的情况共有5种(未给出),其中具有混合效应参数a0㊁a1的模型1.2(模型1.2是在模型1的基础上,通过在参数a0㊁a1上添加随机项组成,是混合模型的基本方法,由于本文涉及到的这类模型较多,以下不一一列举)的3个指标最小.2)对于枝条长度模型,混合模型收敛的情况共有16种(未给出),其中,模型2.1~2.4模拟结果的3个指标值都小于传统最小二乘法,与模型2.3相比,模型2.4的拟合效果没有显著提高(P=0.192).因此,把具有混合效应参数a0㊁a1㊁a2的模型2.3作为枝条长度的最优混合模型.3)对于枝条着枝角度模型,混合模型收敛的情况共有11种(未给出),当考虑1个随机参数时,模型3.1的拟合效果最好;当考虑2个随机参数组合时,模型3.2拟合效果最好;当考虑3个随机参数时,虽然模型3.3的拟合效果比模型3.2有所提高,但并不显著(P=0.284).因此,把具有混合效应参数a0㊁a1的模型3.2作为枝条角度基于树木效应的最优混合模型.
2.3　方差协方差结构
当考虑单木效应时,随机参数反映了树木间的变化,通过添加随机参数能提高模型的拟合精度,似然比的显著性检验说明,树木间的随机效应是显著的.然后又测试了其他2种方差协方差结构:复合对称CS和对角矩阵UN(1).从表5的拟合结果可以看出:对于枝条基径模型,3种方差协方差结构均能提高模型的拟合精度,但采用UN和CS结构的混合模型似然比检验不显著(P=0.646),说明采用CS 结构更适合于描述枝条基径在不同树木间的变化;对于枝条长度模型,采用UN和CS结构的拟合效果好于UN(1),而且采用UN结构矩阵的混合模型显著优于CS结构(P=0.022);对于枝条角度模型,3
表4　基于不同随机效应参数混合的枝条大小模型拟合结果比较
Table4　Comparisons of mixed BD,BL and BA models on different random effects parameters
模型代码
Model code随机参数
Random parameter参数个数
Parameters number赤池信息准则
AIC贝叶斯信息准则
BIC对数似然值
⁃2Log likelihood似然比检验
LRT P 1无None4320.7325.3318.7
1.1a15304.4308.1300.418.3<0.001 1.2a0,a17288.5296.0280.519.9<0.001 2无None5273.0277.6271.0
2.1a16228.7232.4224.746.3<0.001 2.2a0,a1814
3.4150.9135.489.3<0.001 2.3a0,a1,a211140.4153.5126.49.00.029 2.4a0,a1,a2,a315142.3162.9120.3 6.10.192 3无None55945.25949.95943.2
3.1a065896.55900.25892.550.7<0.001 3.2a0,a185879.75887.25871.720.80.000 3.3a0,a1,a3115881.95895.05867.9 3.80.2841542
9期 董灵波等:基于线性混合模型的红松人工林一级枝条大小预测模拟
个评价指标值及似然比检验结果都说明采用UN结构的混合模型拟合效果最好.综上,本文把随机效应具有复合对称矩阵的混合模型选为枝条基径最优混合模型,同时把随机效应具有无结构矩阵的混合模型选为枝条长度和角度最优混合模型.
2.4　基于混合模型的误差结构矩阵
由于本研究所用的枝条数据是采用枝解析技术获取,数据间存在着时间序列相关性.本研究采用复合对称结构CS㊁一阶自回归结构AR(1)以及一阶自回归与滑动平均结构ARMA(1,1)来描述样木内树木内枝条的时间序列相关性,把这3个函数分别加入以上所得枝条基径㊁长度和角度最优混合模型中.
表6的拟合结果表明:当尝试把CS模型加入到枝条基径混合模型中时,模型的拟合效果没有显著提高(P=0.646);当把AR(1)和ARMA(1,1)结构加入到基径混合模型中时,3个评价指标及似然比检验说明,AR(1)和ARMA(1,1)都能显著提高混合模型的拟合效果,但ARMA(1,1)和AR(1)的似然比检验不显著(LRT=0.0,P=1.0),因此AR(1)更适合于描述基径模型的误差相关性.当尝试把CS 结构加入到枝条长度混合模型中时,模型拟合效果没有提高;当把AR(1)和ARMA(1,1)结构加入到长度混合模型中时,3个评价指标及似然比检验说明,模型的拟合效果显著提高,但ARMA(1,1)和
AR(1)的似然比检验不显著(LRT=0.4,P= 0.527),说明AR(1)结构同样适合于描述长度模型的误差相关性.当把3个自相关结构加入角度混合模型中时,并不能提高模型的拟合精度.这可能是不同树种的分枝角度主要受自身遗传因素的控制,对外界不同环境条件的响应较小,因此本研究不考虑枝条着枝角度混合模型的时间序列结构,即Γij=I ij.
考虑单木效应和时间序列相关性后仍不能完全解决模型的异方差问题.通常判断误差的异质性主要通过残差分布图进行比较(图1a~c).从图1b可以看出,枝条长度混合模型的残差分布呈较明显的
表5　基于随机参数不同方差协方差结构的混合模型模拟结果比较
Table5　Comparisons of mixed model based on different variance⁃covariance structure of random parameters
模型
Model方差协方差结构
Variance⁃covariance参数个数
Parameters number赤池信息准则
AIC贝叶斯信息准则
BIC对数似然值
⁃2Log likelihood似然比检验
LRT P值
P value 枝条基径对角矩阵UN(1)6292.6298.2286.6
Branch diameter复合对称矩阵CS6286.7292.3280.7 5.9-无结构矩阵UN7288.5296.0280.50.20.655枝条长度对角矩阵UN(1)8164.8170.5158.8
Branch length复合对称矩阵CS7141.9147.5137.920.9<0.001无结构矩阵UN11140.4153.5126.411.50.022枝条角度对角矩阵UN(1)75894.95900.55888.9
Branch angle复合对称矩阵CS75887.95893.55881.97.0-无结构矩阵UN85879.75887.25871.710.20.001 -不进行似然比检验No likelihood ratio test.下同The same below.
表6　考虑自相关矩阵模型的红松人工林枝条大小模拟结果
Table6　Comparisons of mixed models with different error autocorrelation matrix for Pinus koraiensis plantation
模型
Model时间序列相关结构
Time series
correlation structure 参数个数
Parameters
number
赤池信息
准则
AIC
贝叶斯信息
准则
BIC
对数似然值
-2Log
likelihood
似然比检验
LRT P值
P value
枝条基径无None6286.7292.3280.7
Branch diameter复合对称矩阵CS7288.5296.0280.50.20.655一阶自回归结构AR(1)7274.3281.8266.314.40.000
一阶自回归与滑动平均结构ARMA(1,1)8276.3285.6266.314.40.001
枝条长度无None11140.4153.5126.4
Branch length复合对称矩阵CS12142.4157.3126.40.0-一阶自回归结构AR(1)12132.7148.7117.78.70.003
一阶自回归与滑动平均结构ARMA(1,1)13135.3152.2117.39.10.011
枝条角度无None85879.75887.25871.7
Branch angle复合对称矩阵CS95881.75891.15871.70.0-一阶自回归结构AR(1)95881.75891.15871.70.0-
一阶自回归与滑动平均结构ARMA(1,1)105883.75894.95871.70.0-2542应　用　生　态　学　报 24卷
倒喇叭状,存在着异方差问题.本研究选用林业中常用的2种函数来描述混合模型产生的异方差现象[27-28],方程如下:
CF (1,ijk )=1+s ijk 2/2
(12)CF (2,ijk )=exp(s ijk 2/2)
(13)
式中:CF 1和CF 2均为方差函数;s ijk 为模型误差.
2个异方差函数都能提高枝条基径混合模型的
拟合效果,3个评价指标都说明CF 2函数优于CF 1函数(表7);加入CF 2函数的枝条基径混合模型显示了更均匀的残差分布(图1a).2个异方差函数同样提高了枝条长度混合模型的拟合效果,而且CF 2函数显著优于CF 1函数(LRT =32.5);采用混合模型方法,枝条长度模型的残差并不随预测值的增大而减小,分布虽然并不完全均匀,但异方差现象明显减小(图1b).当把2个异方差函数加入枝条角度混合模型中时,3个评价指标的结果都说明模型的拟合效果得到明显提高,但残差分布图的变化并不明显(图1c).
综上,本文确定的红松人工林一级枝条大小因子混合预测模型如公式(14)~(16)所示:
枝条基径:
Y =^Y +εijk ,^Y =(a 0+u 0)+(a 1+u 1)ln(DINC )+a 2ln(DBH )εijk ~N (0,R ij ),u 0ij u 1æèçö
ø÷ij ~(0,D ),D =CS ,Γij (θ)=AR (1)
R ij =σ2G ij 0.5Γij G ij 0.5=σ2CF 1 ρf ij -1┇⋱┇ρf ij -1 æèççöø÷
÷1CF ,CF =1/CF (2,ijk )i =1, ,m ;j =1, ,n i ;k =1, ,f ìîíïïïïï
ïïïïï
ij
(14)
枝条长度:
Y =^Y +εijk ,^Y =(a 0+u 0)+(a 1+u 1)ln(DINC )+(a 2+u 2)DINC 2+a 3ln(DBH )εijk ~N (0,R ij ),u 0ij u 1ij u 2æèççöø÷÷ij ~(0,D ),D =UN
Γij (θ)=AR (1)
R ij
=σ2G ij 0.5Γij G ij
0.5=σ2CF 1 ρf ij -1
┇⋱┇ρf ij -1 æèççöø÷÷1CF
CF =1/CF
(2,ijk )
i =1, ,m ;j =1, ,n i ;k =1, ,f ìîíï
ïï
ï
ïï
ï
ï
ï
ïïïï
ïïï
ij
(15)
枝条角度:
Y =^Y+εijk ,^Y =(a 0+u 0)+(a 1+u 1)ln(DINC )+a 2DBH+a 3BD εijk ~N (0,R ij ),u 0ij u 1æèçöø÷ij ~(0,D ),D =UN R ij =σ2G ij 0.5Γij G ij 0.5=σ2CF ×I f ij ×CF ,CF =1/CF
(1,ijk )
i =1, ,m ;j =1, ,n i ;k =1, ,f ìîíï
ïï
ïï
ïïï
ij
(16)
式中:Y 为枝条大小因子实测值;^Y
为枝条大小因子估计值;CF 为描述方差异质性的对角阵;ρ为待估参数.
2.5　模型评价
利用SAS 软件中的MIXED 模块对上述混合模型进行参数估计.表8给出了枝条基径㊁长度和角度最优混合模型的固定参数估计值㊁随机参数的方差协方差组成㊁时间序列函数的参数估计值,以及各最优混合模型的拟合统计量指标.与表3对比可以看表7　考虑不同异方差函数的红松人工林枝条大小模拟结果*
Table 7　Comparison of mixed models with different variance functions for Pinus koraiensis plantation
模型Model
方差函数Variance function
参数个数Parameters number
赤池信息准则
AIC
贝叶斯信息准则
BIC
对数似然值-2Log likelihood
似然比检验
LRT 枝条基径无None 7274.3281.8266.3Branch diameter 函数CF 17172.5180.0164.5101.8函数CF 27
152.6160.1144.6121.7枝条长度无None 12132.7148.7117.7Branch length 函数CF 112 2.2
17.2
-13.8131.5函数CF 212-30.3-15.3-46.3
164.0枝条角度无None 85879.75887.25871.7Branch angle
函数CF 185072.75080.15064.7807.0函数CF 2
85541.2
5548.7
5533.2338.5
*参数个数相同,不进行似然比检验The same parameters number,not carrying on the likelihood ratio test.
3
5429期 董灵波等:基于线性混合模型的红松人工林一级枝条大小预测模拟
图1　基于传统回归(A)和混合模型方法(B)的枝条基径(a)㊁长度(b)和角度(c)模型残差分布
Fig.1　Residual distribution of branch diameter (a),lenglth (b)and angle (c)based on ordinary regression (A)and mixed⁃effect model (B).
表8　红松人工林枝条大小预测模型参数估计结果及拟合统计量
Table 8　Parameter estimates and goodness⁃of⁃fit statistics of branch size prediction models for Pinus koraiensis planta⁃tion
变量Variable 参数Parameter 枝条基径Branch diameter 枝条长度Branch length 枝条角度Branch
angle
固定参数
a 0 2.091 3.95661.097Fixed parameter a 10.518 1.07010.762a 20.208-0.0130.540a 30.129-0.598误差方差Error variance σ
2
0.0650.047 1.740随机效应方差
σ2a 00.0050.01121.422Random effect variance
σ2a 10.0050.021 5.500σ2a 20.000σ2a 0a 1-0.003-0.013-9.064σ2a 0a 20.000σ2a 1a 2
-0.001时间序列相关结构
Time series correlation structure ρ0.1790.219拟合统计量
R 2
0.7460.9290.345Goodness⁃of⁃fit statistics
MAE 0.1930.1567.615RMSE
0.265
0.221
10.038
出,对于枝条基径㊁长度和角度模型,混合模型的确定系数高于固定效应模型,而平均绝对误差和均方根误差均小于固定效应模型,一方面说明混合模型的模拟精度高于固定效应模型,另一方面计算出的
随机参数也能够体现各样木之间枝条大小因子的差异性.
2.6　模型检验
应用表2中的模型独立样本检验数据,根据式
(1)中的μij =^D ^Z ij T (^R ij +^Z ij ^D ^Z ij T )-1^εijk 和表8中各式的模拟结果,求出12棵验证样木数据的随机效应参
数μij ,进而计算出12棵样木枝条大小因子预测值,具体计算利用Matlab2010a 软件实现.最后利用确定系数㊁平均绝对误差和均方根误差3个指标与传统的最小二乘法进行精度比较.从表9可以看出,枝条基径㊁长度和角度混合模型的确定系数均大于基本模型的确定系数,混合模型的平均绝对误差和均方根误差均小于基本模型,尤其是枝条角度模型的
4542应　用　生　态　学　报 24卷。