文章信息
- 姜立春, 张锐, 李凤日
- Jiang Lichun, Zhang Rui, Li Fengri
- 基于线性混合模型的落叶松枝条长度和角度模型
- Modeling branch length and branch angle with linear mixed effects for Dahurian Larch
- 林业科学, 2012, 48(5): 53-60.
- Scientia Silvae Sinicae, 2012, 48(5): 53-60.
-
文章历史
- 收稿日期:2011-04-11
- 修回日期:2011-08-30
-
作者相关文章
树冠是树木进行光合作用和呼吸作用等生理活动的主要场所,其形状决定了树木的生长状况和林分动态(Chmura et al., 2007; 涂洁, 2007)。树冠形状和结构能用来估计树冠盖度、体积、表面积、树叶分布和透光率(Deleuze et al., 1996;Waguchi, 2004;廖彩霞等, 2007),树冠宽度决定了树木之间的竞争能力及林分郁闭程度(Baldwin et al., 1997;Marshall et al., 2003;Jordan et al., 2007),因此树冠形状的准确表述是许多生物过程模型的重要组成部分。
树冠模型构建过程大致可划分为直接和间接方法。直接方法是把树冠设想成一种或几种不同的几何形状,如圆锥体(cone)、圆柱体(cylinder)、凹曲线体(neiloid)和抛物线体(paraboloid),以林木变量(如胸径、树高、冠长等)、着枝深度(DINC)和相对着枝深度等为自变量来直接构建树冠形状预测模型(Hann, 1999)。直接方法虽然能较好地预测树冠形状,但是不能预测树冠内枝条大小和分布。间接方法是先预测枝条大小及其空间变化特征,如枝长(BL)、着枝角度(BA)等变量,然后利用枝长和着枝角度的三角关系来计算树冠的半径,进而实现树冠形状的预测(Cluzeau et al., 1994; Roeh et al., 1997;Mäkinen et al., 2002;Hein et al., 2008)。间接方法不但能预测树冠形状、了解树冠动态变化,还能作为解释树木生长机制的基本依据之一。
在我国,李凤日(2004)以长白落叶松(Larix olgensis)为对象, 用间接方法建立枝条基径、枝长、着枝角度和弦长等预估模型来预测树冠形状,在技术方法上做出了有益贡献;刘兆刚等(2008)构建了樟子松(Pinus sylvestris var. mongolica)人工林一级枝条基径和枝长模型,为樟子松人工林树冠形状预测奠定了基础。
由于枝条特征数据取自不同的林木和不同的样地,近年来混合模型技术在国外已应用于枝条特征模型中(Hein et al., 2007;Weiskittel et al., 2007)。与传统的回归分析相比,混合模型能得到渐进无偏的参数估计,通过引入随机参数能提高模型的拟合精度,并用方差协方差结构来反映数据间的相关性及异质性。目前国内应用混合模型研究枝条长度和角度模型还未见报道。本文以兴安落叶松(Larix gmelinii)人工林为例,采用线性混合模型技术建立落叶松枝条长度和角度混合效应模型,包括确定参数效应、树木内方差协方差结构、随机效应参数的方差协方差结构及随机参数的校正过程,然后对混合效应模型与传统模型拟合效果进行检验及比较分析。
1 数据与方法 1.1 数据研究地位于黑龙江省伊春市五营林业局境内。五营林业局位于黑龙江省伊春市北部、小兴安岭南坡腹部,129°06′—129°30′E,47°54′—48°19′N。五营区属中温带大陆性湿润季风气候,除受纬度、地理条件和大气环流控制外,还受森林和局部地形影响。五营林区四季气候变化很大, 冬季在极地大陆气团控制下气候严寒、干燥并漫长;夏季受副热带海洋气团的影响,降水集中,雨水充沛,气候湿热,日照时间长,适宜作物生长。年平均气温为0.2 ℃左右,无霜期为111天,全年平均降雨量为626.9 mm,年降雨量最大值为832.7 mm。土壤种类以山地暗棕色森林土为主,少量草甸土、沼泽土、石质土。
用来建立模型的数据来自该局丽林实验林场的10块落叶松人工林样地。样地面积为0.04 hm2, 实测样地内每株林木的胸径和树高。在每块标准地内随机选取优势木、平均木、被压木各1株,树木被伐倒后,测量以下指标: 1)树高; 2)第一活枝高和第一死枝高;3)树冠高度和冠幅。将解析木的树干按1 m区分段进行区分,并在每个区分段的中央位置锯取树干解析圆盘。树冠部分也按1 m区分段进行区分。在每个区分段内对枝条进行编号,并测定每个枝条的总着枝深度(DINC)、枝条的方位角(A)、着枝角度(BA)、基径(BD)、枝长(BL)、弦长(BCL)及弓高(BAH)等枝条特征变量。将所收集全部枝条数据,按80%和20%的比例分成建模数据样本和独立检验样本,分别用于建立和检验枝条长度和角度模型。落叶松人工林样木和枝条特征调查因子的统计量见表 1。
固定效应线性模型形式为:
(1) |
单水平线性混合模型的形式为(Pinheiro et al., 2000):
(2) |
式中:yi是第i株树中(ni×1) 维枝条长度或角度观测值;m是样木数量;ni是第i株数的观察值数量;Xi是(ni×p)维已知设计矩阵;β是(p×1) 维固定参数向量;Zi是(ni×q)维随机效应的设计矩阵;bi是(q×1) 维随机参数向量;D是随机参数的方差协方差矩阵;Ri是树木内方差协方差结构;εi是(ni×1) 维模型的误差项。
根据文献(Pinheiro et al., 1998), 构建以上混合模型需要确定以下3个步骤。
1)确定参数效应。Pinheiro等(2000)建议如果模型能收敛,将不同随机参数组合的模型进行拟合,比较模型拟合的统计量,即比较AIC、BIC和对数似然值指标,AIC和BIC越小越好,对数似然值越大越好。为了避免过多参数化问题,具有不同参数个数的模型还要进行似然比检验(LRT),即利用LRT和P值进行检验,如P值小于0.05即认为差异显著。本文按照此方法进行随机效应参数确定。
2)确定树木内方差协方差结构(Ri)。为了确定树木内方差协方差结构,必须解决树木内误差相关性和异方差问题。目前,在统计和林业上基本都采用式(3)来描述(Davidian et al., 1995;Calama et al., 2004; Trincado et al., 2006):
(3) |
式中:σ2指模型的误差方差值;Γi为树木内误差相关性结构;Gi为描述方差异质性的对角矩阵。
3)确定随机效应的方差协方差结构(D)。树木间的方差协方差结构反映了树木间的变化。本研究检验了林业上3种常用的方差协方差结构,包括复合对称(CS)、对角矩阵(diagonal matrix)、广义正定矩阵(general positive-definite matrix)(Fang et al., 2001; Jordan et al., 2005)。
1.2.2 模型评价和检验指标拟合和检验结果通过以下统计量和指标评价:AIC、BIC和Log Likelihood 3个统计量指标,绝对误差(Bias)、均方根误差(RMSE)、确定系数(R2)。
式中:yij为观测值;ŷij为预测值;ȳ为观测值的平均值;ni为树木内观察值数量;m为树木株数;n为样本数。
1.2.3 模型检验模型的独立性检验采用建模时未使用的独立样本(检验样木)数据,对所确定模型的预测性能进行综合评价。混合效应模型中固定效应部分的检验与传统的检验方法相同;但是,随机效应部分的检验需要二次抽样来计算随机参数值。本研究使用Vonesh等(1997)的方法来计算随机参数bk:
(4) |
式中:D̂为随机效应参数的方差协方差矩阵;R̂k为树木内方差协方差结构;Ẑk为设计矩阵;êk为实际值减去用固定效应参数计算的预测值。
2 结果与分析 2.1 固定模型建立通过绘制散点图分析枝条长度和角度与胸径、树高、冠幅、冠长、总着枝深度等林木变量的关系。研究发现,枝条长度与枝条的总着枝深度有密切关系,并随着总着枝深度的增加而增加,到达最大值之后有所下降;枝条角度与枝条的总着枝深度也有一定相关性;枝条长度和角度与林木胸径、树高、冠幅、冠长等变量组也有一定的相关关系。基于这些变量,采用逐步回归技术建立如下枝条长度和角度模型:
(5) |
(6) |
式中: BL为枝条长度(cm);BA为枝条角度(°);DINC为枝条的总着枝深度(m);DBH为胸径(cm);b1,b2,b3,b4为模型预估参数。
2.2 随机效应的方差协方差结构(D)混合模型构建首先要确定随机效应的方差协方差结构。本研究主要考虑了林业上常用的3种方差协方差结构:复合对称、对角矩阵、广义正定矩阵。以包括3个随机参数(b1,b2,b3)的方差协方差结构为例,复合对称:
(7) |
对角矩阵:
(8) |
广义正定矩阵:
(9) |
方差协方差结构比较过程中假设误差效应满足独立同分布的假设,随机效应默认包括所有的随机参数,结果见表 2。从表 2可以看出:无论考虑枝条长度还是枝条角度模型,广义正定矩阵结构始终显示了最高的拟合效果。因此,把随机效应具有广义正定矩阵结构选为枝条长度和角度最优随机效应的方差协方差结构。
考虑树木效应和不同随机参数的组合,误差效应满足独立同分布的假设,利用S-PLUS软件的LME模块对模型(5)进行拟合。利用AIC、BIC、Log Likelihood及似然比检验等统计量指标对模型的拟合优度进行比较,拟合结果见表 3。从表 3可以看出:混合模型收敛的情况共有7种,7次混合模型模拟结果的AIC和BIC值都比基于最小二乘法的模型(5)值小。当考虑随1个随机参数时,模型(5.2)优于模型(5.1)和(5.3);当考虑随2个随机参数组合时,模型(5.6)优于模型(5.4)和(5.5);当考虑随3个随机参数时,模型(5.7)优于模型(5.6)。
为了避免过多参数化问题的产生,还要对不同混合模型进行显著性检验,即利用LRT和P值进行方差分析,如P < 0.05即认为差异显著。本文主要对具有代表性的基本模型(5)、模型(5.2)、模型(5.6)和模型(5.7)进行似然比检验,检验结果表明基本模型(5)和模型(5.2)有显著不同(P < 0.000 1),模型(5.2)和模型(5.6)有显著不同(P < 0.000 1),模型(5.6)和模型(5.7)有显著不同(P < 0.000 1)。综合以上分析,模型(5.7)拟合效果最好。
2.4 枝条角度混合模型拟合结果考虑树木效应和不同随机参数的组合,利用S-PLUS软件的NLME模块对模型(6)进行拟合, 拟合结果见表 4。从表 4可以看出:混合模型收敛的情况共有7种,7次混合模型模拟结果的AIC和BIC值都比基于最小二乘法的模型(6)值小。当考虑随1个随机参数时,模型(6.1)优于模型(6.2)和(6.3);当考虑随2个随机参数组合时,模型(6.4)优于模型(6.5)和(6.6);当考虑随3个随机参数时,模型(6.7)优于模型(6.4)。对基本模型(6)、模型(6.1)、模型(6.4)和模型(6.7)进行似然比检验,检验结果表明基本模型(6)和模型(6.1)有显著不同(P < 0.000 1),模型(6.1)和模型(6.4)有显著不同(P < 0.000 1),模型(6.4)和模型(6.7)有显著不同(P < 0.000 1)。综合以上分析,模型(6.7)拟合效果最好。
混合模型在参数效应和随机效应的方差协方差结构确定后,还必须确定误差的自相关性和异质性。由于枝条数据不是连续观测数据,因此本研究误差自相关性不考虑。通常判断误差的异质性主要通过残差分布图进行比较(图 1)。从图 1可以看出:枝条长度混合模型的残差分布显示了极为明显的喇叭状,即异方差性。本研究选用混合模型中常用的幂函数和指数函数来描述混合模型产生的异方差现象,结果见表 5。从表 5可以看出:幂函数和指数函数都能提高枝条长度混合模型的拟合效果,3个评价指标都说明指数函数优于幂函数。从残差分布图 1对比可以看出:加入指数函数的枝条长度混合模型显示了均匀的方差分布。用相同的方法对枝条角度混合模型进行了测试,结果见表 5和图 2。从表 5和图 2可以看出:幂函数和指数函数都能提高枝条角度混合模型的拟合效果,3个评价指标都说明幂函数优于指数函数,但残差分布图变化并不明显。指数函数和幂函数的方程如下:
(10) |
(11) |
式中:μij为基于固定效应参数的预测值;α, β为指数函数和幂函数的参数。
表 6给出了枝条长度和角度最优混合模型的固定参数估计值、随机参数的方差协方差组成、指数函数和幂函数的参数估计值,并给出了最优混合模型与传统基本模型拟合统计量的比较。结果表明:无论考虑枝条长度和角度模型拟合时,混合模型的确定系数大于基本模型的确定系数,而混合模型的绝对误差和均方根误差均小于基本模型的绝对误差和均方根误差,这些评价指标说明引入随机参数提高了模型的拟合精度。
对剩余的落叶松枝条长度和角度数据进行验证。用式(4)计算随机参数值,具体计算采用SAS软件的PROC IML模块计算。最后利用绝对误差和均方根误差2个指标与传统的最小二乘法进行预测精度比较,检验结果如表 7所示。从表 7可以看出:枝条长度和角度混合模型绝对误差和均方根误差均小于基本模型的绝对误差和均方根误差,这些评价指标说明引入随机参数提高了模型的预测精度。
混合模型在应用预测时包括平均预测和个体预测。平均预测相当于传统的非线性回归方法,即用固定效应参数来预测;而个体预测在固定效应的基础上还要计算随机参数值。本研究以枝条长度的一个实际例子来详细说明。为了减小矩阵的维数,从检验数据中随机抽取着枝深度和枝条长度的6次测量值(0.9, 98),(1.9, 106),(3.5, 255),(4.1, 251),(4.6, 271),(5.7, 245)来计算,首先计算êk,R̂i和Ẑk,然后利用式(4)计算随机参数b̂k:
因此,基于随机参数值和固定参数值计算混合模型枝条长度预测值为(74.1, 152.5, 233.8, 250.3, 258.2, 257.1),绝对误差和均方根误差为(19.516 0,24.114 0);而只利用固定效应参数方法计算枝条长度预测值为(68.5, 122.4, 187.7, 205.6, 217.7, 235.5),绝对误差和均方根误差为(36.907 8,42.157 9)。由此可以看出:考虑随机效应参数混合模型的预测精度明显高于只考虑固定效应参数的方法,通过计算随机参数值,枝条长度的个体预测优于平均预测。
3 结论与讨论本研究采用逐步回归技术建立了落叶松人工林枝条长度和角度模型:BL=b1+b2DINC+b3DINC2+b4DBH·DINC2,BA=b1+b2DINC+b3DINC2+b4DBH·DINC。
基于所建基本模型,利用线性混合模型方法构建了枝条长度和角度混合模型。利用AIC、BIC、对数似然值及似然比检验评价不同线性混合模型的效果表明,混合模型的拟合精度都比基本模型的拟合精度高。无论拟合枝条长度模型还是枝条角度模型,b1,b2,b3同时作为混合参数时模型拟合最好。为了解决树木内方差异质性问题,把指数函数和幂函数加入到枝条长度和角度混合模型中,指数函数较好地描述了枝条长度混合模型的异方差现象,幂函数较好地描述了枝条角度混合模型的异方差现象。虽然从统计上指数函数和幂函数都能显著提高枝条长度和角度混合模型的拟合效果,但实际的残差图对比并不明显,尤其是枝条角度残差图的变化,因此枝条长度和角度混合模型在应用时异方差现象可以不考虑。为了进行混合模型检验,本研究还探讨了随机参数值的校正过程,通过校正随机参数值混合模型能提高模型的预测精度。因此混合模型在应用上不仅能反映总体平均枝条长度和角度预测,而且能通过方差协方差结构校正随机参数来反映树木之间的差异。
树冠枝条特征变化受年龄、林分、立地条件、经营措施及树木生长水平等多种因子的影响与制约,因此本文所构建的落叶松枝条长度和角度预测模型仅适用于与本研究立地条件相似的地区,大尺度模型的建立还需考虑其他因子的影响,随着数据的积累,这方面将进一步深入研究。
[] | 李凤日. 2004. 长白落叶松人工林树冠形状的模拟. 林业科学, 40(5): 16–24. DOI:10.11707/j.1001-7488.20040503 |
[] | 廖彩霞, 李凤日. 2007. 樟子松人工林树冠表面积及体积预估模型的研究. 植物研究, 27(4): 478–483. |
[] | 刘兆刚, 舒扬, 李凤日. 2008. 樟子松人工林一级枝条基径和枝长模型的研究. 植物研究, 28(2): 244–248. |
[] | 涂洁, 刘琪. 2007. 中亚热带红壤丘陵区湿地松枝条生长规律. 浙江林学院学报, 24(2): 162–167. |
[] | Baldwin V C, Peterson K D. 1997. Predicting the crown shape of loblolly pine trees. Can J For Res, 27(1): 102–107. DOI:10.1139/x96-100 |
[] | Calama R, Montero G. 2004. Interregional nonlinear height-diameter model with random coefficients for stone pine in Spain. Can J For Res, 34(1): 150–163. DOI:10.1139/x03-199 |
[] | Chmura D J, Rahman M S, Tioelker M G. 2007. Crown structure and biomass allocation patterns modulate aboveground productivity in young loblolly pine and slash pine. For Ecol Manage, 243(2/3): 219–230. |
[] | Cluzeau D, LeGoff N, Ottorini J M. 1994. Development of primary branches and crown profile for Fraxinus excelsior. Can J For Res, 24(12): 2315–2323. DOI:10.1139/x94-299 |
[] | Davidian M, Giltinan D M. 1995. Nonlinear models for repeated measurement data. Chapman and Hall, London, UK. |
[] | Deleuze C, Herve J C, Colin F, et al. 1996. Modeling crown shape of Picea abies: spacing effects. Can J For Res, 26(11): 1957–1966. DOI:10.1139/x26-221 |
[] | Fang Z, Bailey R L. 2001. Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments. For Sci, 47(3): 287–300. |
[] | Hann D W. 1999. An adjustable predictor of crown profile for stand-grown Douglas-fir trees. For Sci, 45(2): 217–225. |
[] | Hein S, Mäkinen H, Yue C, et al. 2007. Modeling branch characteristics of Norway spruce from wide spacings in Germany. For Ecol Manage, 242(2/3): 155–164. |
[] | Hein S, Weiskittel A R, Kohnle U. 2008. Branch characteristics of widely spaced Douglas-fir in south-western Germany: comparisons of modelling approaches and geographic regions. For Ecol Manage, 256(5): 1064–1079. DOI:10.1016/j.foreco.2008.06.009 |
[] | Jordan G J, Ducey M J. 2007. Predicting crown radius in eastern white pine (Pinus strobus L.) stands in New Hampshire. North J Appl For, 24(1): 61–64. |
[] | Jordan L, Daniels R F, Clark Ⅲ, et al. 2005. Multilevel nonlinear mixed effects models for the modeling of earlywood and latewood microfibril angle. For Sci, 51(4): 357–371. |
[] | Marshall D D, Johnson G P, Hann D W. 2003. Crown profile equations for stand-grown western hemlock trees in northwestern Oregon. Can J For Res, 33(11): 2059–2066. DOI:10.1139/x03-126 |
[] | Mäkinen H, Song T. 2002. Evaluation of models for branch characteristics of Scots pine in Finland. For Ecol Manage, 158(1/3): 25–39. |
[] | 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. |
[] | Roeh R L, Maguire D A. 1997. Crown profile models based on branch attributes in coastal Douglas-fir. For Ecol Manage, 96(1/2): 77–100. |
[] | Trincado G, Burkhart H E. 2006. A generalized approach for modeling and localizing stem profile curves. For Sci, 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. |
[] | Waguchi Y. 2004. Accuracy and precision of crown profile, volume, and surface area measurements of 29-year-old Japanese cypress trees using a Spiegel relascope. J For Res, 9(2): 173–176. DOI:10.1007/s10310-003-0060-0 |
[] | Weiskittel A R, Maguire D A, Monserud R A. 2007. Modeling crown structural responses to competing vegetation control, thinning, fertilization, and Swiss needle cast in coastal Douglas-fir of the Pacific Northwest, USA. For Ecol Manage, 245(1/3): 96–109. |