文章信息
- 李耀翔, 姜立春, 李凤日
- Li Yaoxiang, Jiang Lichun, Li Fengri
- 基于混合模型的落叶松微纤丝角模型
- Modeling Microibril Angle with Mixed Models for Dahurian Larch
- 林业科学, 2012, 48(4): 81-86.
- Scientia Silvae Sinicae, 2012, 48(4): 81-86.
-
文章历史
- 收稿日期:2010-11-30
- 修回日期:2011-02-06
-
作者相关文章
2. 东北林业大学林学院 哈尔滨 150040
2. College of Forestry, Northeast Forestry University Harbin 150040
微纤丝角(microfibril angle, MFA)为细胞壁中纤维素纤维与细胞纵向轴线之间的夹角(Lichtenegger et al., 1999),其角度越小则细胞的抗张强度越大。微纤丝角是木材机械性能的主要决定因子之一,特别是弹性模量、断裂模量及纵向和切向收缩性,微纤丝角与木材密度存在一定的相互关系,并与木材强度和硬度密切相关(Cave et al., 1994; Walker et al., 1995; MacDonald et al., 2002; Fang et al., 2004; Lundgren, 2004; Jordan et al., 2005;2007)。目前微纤丝角已成为林产品加工业木材质量评价的一个重要指标,以往主要是从定性分析的角度研究木材微纤丝角随生长轮年龄的变异规律(Addis et al., 1995; 洑香香等,2002; 汪贵斌等,2007; 曹琳等,2009),研究结果表明:微纤丝角在髓心处最大,从髓心到树皮微纤丝角逐渐变小。
长时间以来,木材科学工作者一直在探索从定量分析的角度来研究木材微纤丝角随生长轮年龄变异的规律,通过建立数学模型来模拟微纤丝角的变化,主要是用线性和非线性回归过程来分析微纤丝角模型(鲍甫成等,1999; 杨文忠等,2004; 马顺兴等,2006)。由于微纤丝角在不同树木之间及同一树木不同生长轮之间有较大变异,因此近年来混合模型技术在国外已应用于微纤丝角模型中(Jordan et al., 2005)。与传统的回归分析相比,混合模型能得到渐进无偏参数估计,通过引入随机参数能提高模型的拟合精度,并用方差协方差结构来反映数据间的相关性及异质性。混合模型包括总体和个体预测,总体预测只能反映平均水平而不能体现树木之间的差异,而个体预测能通过二次抽样校正随机参数来反映树木之间的差异。目前国内应用混合模型研究微纤丝角还未见报道。本文以兴安落叶松(Larix gmelinii)人工林为对象,采用非线性混合模型方法建立落叶松木材微纤丝角径向混合效应模型,包括确定固定效应参数、随机效应参数、方差协方差结构及相关性结构,并对混合效应模型中随机参数值的校正进行探讨,然后对混合效应模型与传统模型拟合效果进行检验及比较分析。
1 数据与方法 1.1 数据研究地点位于黑龙江省七台河境内。七台河市位于黑龙江省东部的张广才岭与完达山脉两大山系衔接地带,地处低山丘陵带。地理坐标为130°06′—131°58′E,45°16′—46°37′N。七台河市处于中温带湿润气候区,四季分明,降水各季分布不均。冬季长而干燥寒冷,夏季短而湿热多雨,春、秋两季气候多变,春季回暖快、风大而少雨干旱,秋季降温快、来霜早。年平均温度4.0 ℃,最冷月(1月)平均气温-18.3 ℃, 极端最低气温-39.0 ℃,最热月(7月)平均气温21.9 ℃, 极端最高气温37.4 ℃。年平均无霜期128天,年平均降水量549 mm。土壤以暗棕壤为主。
落叶松试样采自七台河市金沙林场3块立地条件相似的样地。在每块样地内伐取优势木、平均木和被压木各1株,共9株。标准木伐倒后,测定胸径、树高、冠幅等因子,并从树高的1.3 m处截取5 cm厚的圆盘各1个。在每个圆盘的中心位置过髓心取2 cm宽的长条,从髓心到树皮,沿南向逐年轮取宽度为1.5 cm的木块试样。由于标准木年龄较大,一些年轮处晚材不好区分,为了避免试验误差,本研究只对早材切取弦切面切片,切片厚度6~10 μm,用10%铬酸和10%硝酸等量混合液浸泡离析后,置于偏光显微镜下测定微纤丝角。每个年轮内测取10个微纤丝角数据。把微纤丝角数据分成建模数据和检验数据,各样木生长状况及相关指标见表 1。
综合国内外参考文献,研究木材微纤丝角的模型主要有直线方程、对数方程、多项式方程、乘幂式方程、指数式方程和修正的Logistic方程(杨文忠等,2004; Jordan et al., 2005)。本研究采用这6种常用的模型形式:
(1) |
(2) |
(3) |
(4) |
(5) |
(6) |
式中: y为微纤丝角; x为年龄; β1, β2, β3为模型待定参数; ε为模型的误差项。
1.2.2 非线性混合模型非线性混合模型是在非线性模型基础上通过引入随机效应而构建的。在本研究中,单水平非线性混合模型的形式为(Pinheiro et al., 2000):
(7) |
式中: yij是第i株树中第j次微纤丝角观测值; m是样木数量; ni是第i株数的观察值数量; xij是相应于微纤丝角观测值的生长轮年龄; β是(p×1)维固定参数向量;bi是(q×1)维随机参数向量; Aij和Bij是包含有0和1的设计矩阵; D是随机参数的方差协方差矩阵; Ri是树木内方差协方差结构; εij是模型的误差项。
根据文献(Pinheiro et al., 1998), 构建以上混合模型需要确定3个步骤。
1) 确定参数效应。混合模型构建过程中最重要的一步是确定参数效应,即哪个参数为固定参数,哪个参数为混合参数。Pinheiro等(2000)建议:如果模型能收敛,首先应把模型中所有参数看成是混合参数,然后将不同随机参数组合的模型进行拟合,比较模型拟合的统计量,即比较AIC、BIC和对数似然值指标,AIC和BIC越小越好,对数似然值越大越好。为了避免过多参数化问题,具有不同参数个数的模型还要进行似然比检验(LRT),即利用LRT和P值进行检验,如P<0.05即认为差异显著。本文按照此方法进行随机效应参数确定。
2) 确定树木内方差协方差结构(Ri)。为了确定树木内方差协方差结构,必须解决树木内误差相关性和异方差问题。目前,在统计和林业上基本都采用式(8)来描述(Davidian et al., 1995; Calama et al., 2004; Trincado et al., 2006):
(8) |
式中: σ2为模型的残差方差值; Γi为树木内误差相关性结构; GI为描述方差异质性的对角矩阵。
3) 确定随机效应的方差协方差结构(D)。树木间的方差协方差结构反映了树木间的变化。本研究检验了林业上4种常用的方差协方差结构,包括无结构(UN)、复合对称(CS)、对角矩阵(diagonal matrix)、广义正定矩阵(general positive-definite matrix)(Fang et al., 2001; Jordan et al., 2005)。
1.2.3 模型评价和检验指标拟合和检验结果通过以下统计量和指标评价: AIC、BIC和Log Likelihood 3个统计量指标, 确定系数(R2), 绝对误差(Bias), 均方根误差(RMSE)。
式中: yij为观测值;
模型的独立性检验采用建模时未使用的独立样本(检验样木)数据,对所确定模型的预测性能进行综合评价。混合效应模型中固定效应部分的检验与传统的检验方法相同; 然而,随机效应部分的检验需要二次抽样来计算随机参数值。本研究使用Vonesh等(1997)的方法来计算随机参数bk:
(9) |
式中:
利用S-Plus软件对式(1)~(6)6个模型采用最小二乘回归参数估计法进行拟合,采用绝对误差、均方根误差及确定系数对模拟结果进行比较分析,结果见表 2。从表 2的拟合结果可以看出:模型(6)的R2值比其他5个模型都高,并且RMSE的值比其他5个模型都低。虽然模型(1)和(2)的Bias值都小于模型(6),但残差分布图显示模型(6)的分布要优于模型(1)和(2)的分布(未给)。模型(6)参数t检验证明模型(6)所有参数估计都是显著的(P<0.000 1)。模型检验(F =9 593.11,P<0.000 1)表明模型(6)对于描述落叶松微纤丝角的径向变化有显著意义。
通常确定参数效应和随机效应的方差协方差结构须同时进行,首先把随机效应的方差协方差结构设为广义正定矩阵,利用S-Plus软件的NLME模块对模型(6)的不同随机参数组合进行拟合。利用AIC、BIC、Log Likelihood及似然比检验等统计量指标对模型的拟合优度进行比较,具有相同参数个数的模型,AIC、BIC指标越小和Log Likelihood指标越大说明拟合的精度越高。为了避免过多参数化问题,具有不同参数个数的模型还采用似然比检验来比较,如P<0.05即认为差异显著。从表 3的拟合结果可以看出:混合模型收敛的情况共有7种,当考虑1个随机参数时,模型(6.3)优于模型(6.1)和(6.2), 模型(6.3)和模型(6)的似然比检验的显著性(P<0.000 1)说明增加1个随机参数的拟合效果要比传统最小二乘法好; 当考虑随2个随机参数组合时,模型(6.5)优于模型(6.4)和(6.6), 模型(6.3)和模型(6.5)的似然比检验的显著性(P<0.000 1)说明增加2个随机参数的拟合效果要比增加1个随机参数的拟合效果好;当考虑随3个随机参数组合时,模型(6.7)和模型(6.5)的似然比检验的显著性(P<0.000 1)说明增加3个随机参数的拟合效果要比增加2个随机参数的拟合效果好。当考虑树木效应时,随机参数反映树木间的变化,通过添加随机参数能提高模型的拟合精度,似然比检验的显著性说明树木间的随机效应是显著的。然后又测试了其他3种方差协方差结构:无结构、复合对称、对角矩阵。广义正定矩阵结构始终显示了最高的拟合精度,因此,把随机效应具有广义正定矩阵结构的混合模型选为最优混合模型。随机效应(b1,b2,b3)的方差协方差结构如下:
(10) |
式中: σb12为随机参数b1的方差; σb22为随机参数b2的方差; σb32为随机参数b3的方差; σb1b2为随机参数b1和b2的协方差; σb1b3为随机参数b1和b3的协方差; σb2b3为随机参数b2和b3的协方差。
2.3 方差协方差结构混合模型在参数效应确定后,还必须确定误差的异质性和树木内误差相关性。通常判断误差的异质性主要通过残差分布图进行比较(图 1)。从图 1可以看出:非线性混合模型的残差分布没有显示极不规则的形状,如喇叭状、抛物线状等,而分布范围要比传统的非线性回归模型小得多,因此异方差的影响在本研究中不考虑,即Ri=σ2ΓiIi。
为了表达树木内的误差相关性,复合对称模型(CS)、一阶自回归模型AR(1)、一阶移动平均模型MA(1)及一阶自回归与移动平均模型[ARMA(1, 1)]经常被用在混合模型中(Fang et al., 2001; Calama et al., 2004; Jordan et al., 2005)。本研究采用这4个模型来描述树木内误差相关性,结果见表 4。从表 4可以看出:CS结构不能提高模型的拟合精度。当把AR(1),MA(1)和[ARMA(1, 1)]模型加入到混合模型中时,3个评价指标及似然比检验都说明AR(1),MA(1)和[ARMA(1, 1)]模型能显著提高原有混合模型的精度。似然比检验证明[ARMA(1, 1)]模型优于AR(1)模型(P<0.000 1),因此本研究采用[ARMA(1, 1)]模型来描述微纤丝角的树内径向误差相关性。
本研究所建立的最终混合模型参数估计如下。
固定参数: b1=24.245 81,b2=0.064 15,b3=7.939 17。
方差协方差: σ2=0.405 31,σ2b1=0.000 07,
σ2b2=0.000 03,σ2b3=0.154 87,
σb1b2=-0.000 03,σb1b3=0.002 49,
σb2b3= -0.002 06。
[ARMA(1, 1)]模型参数: ρ=0.993 14,θ=-0.800 96。
2.4 模型评价混合模型的检验首先需要计算随机参数值,这是混合模型在实际预测过程中最重要的一步。本研究以一个实例来详细说明。为了减小矩阵的维数,从检验数据中随机抽取年龄和微纤丝角的6次测量值(4, 17),(10, 12.5),(15, 10.5),(23, 10),(27, 8.4),(31, 7.7)来计算,首先计算
因此,混合模型微纤丝角的预测值为(17.03, 14.35, 12.45, 10.17, 9.35, 8.71),绝对误差和均方根误差为(0.994 4,1.239 2);而传统非线性回归方法估计值为(18.61, 16.27, 14.56, 12.38, 11.55, 10.88),绝对误差和均方根误差为(3.025 7,3.136 0)。由此可以看出:混合模型的预测精度明显高于传统非线性回归方法。
基于以上计算步骤,对所有检验数据进行验证,具体采用SAS软件的PROC IML模块来计算。最后利用确定系数、绝对误差和均方根误差3个指标与传统的最小二乘法进行预测精度比较。计算及比较结果如表 5。从表 5可以看出:混合模型的确定系数大于基本模型的确定系数,而混合模型的误差和均方根误差均小于基本模型的误差和均方根误差,这些评价指标说明引入随机参数和误差相关性结构提高了模型的预测精度。
本研究采用非线性混合模型的方法建立落叶松人工林微纤丝角模型。将不同随机参数组合的模型进行拟合,利用AIC、BIC、对数似然值及似然比检验评价不同非线性混合模型的效果表明:混合模型的拟合精度都比基本模型的拟合精度高,b1, b2, b3同时作为混合参数时模型拟合效果最好。为了解决树木内误差相关性问题,把相关性结构包括复合对称结构(CS)、一阶自回归结构AR(1)、一阶移动平均结构MA(1)及一阶自回归与移动平均结构[ARMA(1, 1)]加入到微纤丝角最优混合模型中,一阶自回归与移动平均结构[ARMA(1, 1)]显著提高了微纤丝角混合模型的拟合精度。这与以往的研究结果基本一致,如Jordan等(2005)利用混合模型提高了美国南部火炬松(Pinus taeda)微纤丝角的预测,但在表达树内的误差相关性时,Jordan发现一阶自回归结构AR(1)更好地解释了长叶松树木内误差相关性。而本研究发现[ARMA(1, 1)]模型更适合描述落叶松树内误差相关性结构。为了进行混合模型检验,本研究还探讨了随机参数值的校正过程,通过校正随机参数值混合模型能提高模型的预测精度。因此混合模型在应用上不但能反映总体微纤丝角预测,还能通过方差协方差结构和相关性结构校正随机参数来反映个体微纤丝角差异,而传统的非线性回归分析只能反映总体微纤丝角的平均变化规律。
准确预测树木的微纤丝角可以提高木材质量的预测和使用,且在实际应用中,通过对活立木在胸高处(1.3 m)用生长锥钻树心的方法容易得到这些数据,因此本研究建立的落叶松微纤丝角混合模型具有一定的实用性。由于样本数量有限,本研究只考虑了单水平的混合效应模型,即只考虑了树木间的随机效应,随着数据的积累,在多层混合微纤丝角模型方面将进一步深入研究。
[] | 鲍甫成, 江泽慧, 刘盛全. 1999. 人工林杨树材性与生长轮年龄和生长速度关系的模型. 林业科学, 35(1): 77–82. |
[] | 曹琳, 赵广杰. 2009. 毛白杨微纤丝角在株内的变异. 北京林业大学学报, 31(1): 67–70. |
[] | 洑香香, 杨文忠, 方升佐. 2002. 木材微纤丝角研究的现状和发展趋势. 南京林业大学学报:自然科学版, 26(6): 83–87. |
[] | 马顺兴, 王军辉, 张守攻, 等. 2006. 日本落叶松无性系微纤丝角遗传变异的研究. 林业科学研究, 19(2): 188–191. |
[] | 汪贵斌, 曹福亮, 柳学军, 等. 2007. 落羽杉种源木材微纤丝角和纤维形态的变异. 林业科学, 43(6): 117–122. |
[] | 杨文忠, 方升佐. 2004. 杨树无性系微纤丝角的时空变异模式. 东北林业大学学报, 33(1): 25–28. |
[] | Addis T, Buchanan A H, Walker J C F. 1995. A comparison of density and stiffness for predicting wood quality or density: the lazy man's guide to wood quality. Journal of the Institute of Wood Science, 13(6): 539–543. |
[] | 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. DOI:10.1139/x03-199 |
[] | Cave I D, Walker J C F. 1994. Stiffness of wood in fast-grown plantation softwoods: the influence of microfibril angle. Forest Products Journal, 44(5): 43–48. |
[] | Davidian M, Giltinan, D M. 1995. Nonlinear models for repeated measurement data. Chapman and Hall, London, UK. |
[] | Fang S Z, Yang W Z, Fu X X. 2004. Variation of microfibril angle and its correlation to wood properties in poplars. Journal of Forestry Research, 15(4): 261–267. DOI:10.1007/BF02844949 |
[] | Fang Z, Bailey RL. 2001. Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments. Forest Science, 47(3): 287–300. |
[] | 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, He R, Hall D B, et al. 2007. Variation in loblolly pine microfibril angle in the southeastern United States. Wood and Fiber Science, 39(2): 352–363. |
[] | Lichtenegger H, Reiterer A, Stanzl-Tschegg S E, et al. 1999. Variation of cellulose microfibril angles in softwoods and hardwoods-A possible strategy of mechanical optimization. Journal of Structural Biology, 128(3): 257–269. DOI:10.1006/jsbi.1999.4194 |
[] | 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. DOI:10.1093/forestry/75.2.107 |
[] | 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. New York, Springer: 528. |
[] | 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. New York, Marcel Dekker. |
[] | Walker J C F, Butterfield B G. 1995. The importance of microfibril angle for the processing industries. New Zealand Forestry, 40(4): 34–40. |