2016, Vol. 52
表(Table )
大兴安岭不同区域兴安落叶松可变指数削度方程
文章信息
- 姜立春, 马英莉, 李耀翔
- Jiang Lichun, Ma Yingli, Li Yaoxiang
- 大兴安岭不同区域兴安落叶松可变指数削度方程
- Variable-Exponent Taper Models for Dahurian Larch in Different Regions of Daxing'anling
- 林业科学, 2016, 52(2): 17-25
- Scientia Silvae Sinicae, 2016, 52(2): 17-25.
- DOI: 10.11707/j.1001-7488.20160203
-
文章历史
- 收稿日期：2015-03-19
- 修回日期：2015-11-20
-
作者相关文章
引用本文
姜立春, 马英莉, 李耀翔. 2016. 大兴安岭不同区域兴安落叶松可变指数削度方程[J]. 林业科学, 52(2): 17-25.
Jiang Lichun, Ma Yingli, Li Yaoxiang. 2016. Variable-Exponent Taper Models for Dahurian Larch in Different Regions of Daxing'anling. Scientia Silvae Sinicae, 52(2): 17-25. DOI: 10.11707/j.1001-7488.20160203
大兴安岭不同区域兴安落叶松可变指数削度方程
1. 东北林业大学林学院 哈尔滨 150040;
2. 东北林业大学工程技术学院 哈尔滨 150040
2. 东北林业大学工程技术学院 哈尔滨 150040
摘要：[目的] 基于林业上广泛应用的Kozak(1988),Kozak(1994),Kozak(2001)和Kozak(2002)可变指数削度方程,构建适合大兴安岭不同区域兴安落叶松树干削度方程,同时检验不同区域的显著性。[方法] 采用大兴安岭3个区域的兴安落叶松样木干形数据,利用SAS软件的非线性回归SUR法拟合4个削度方程,采用调整确定系数(Radj2)、平均误差(MAB)、均方根误差(RMSE)、相对误差(MPB)、预估精度(P%)和多重共线性指标(CN)以及相应的残差分布图等对削度方程进行综合比较分析。区域性检验采用非线性额外平方和方法,该方法需要拟合完整模型和简化模型。[结果] 1) Kozak(1988)和Kozak(1994)削度方程拟合和预估精度较高,但是存在较高的多重共线性问题;其他削度方程降低了模型的多重共线性,但可提高对树干上部的预测能力。2)不同区域模型的F检验分析发现,区域3与区域1和区域2的树干削度相差较大,区域1和区域2的树干削度相差较小,任何2个区域的对比检验都是显著的(P<0.0001),说明模型在不同区域不能共用一套参数,应有不同的区域参数估计。3)从不同区域各模型的干曲线模拟可以看出,同一模型在3个区域的干曲线模拟结果不同,尤其模型(3)在3个区域中的干曲线更明显地体现出3个区域的不同。模型(1)、模型(2)和模型(4)在3个区域的干曲线模拟体现在区域1和区域2的模拟结果比较接近,区域3则与其有明显不同,该结果与F检验结果一致。不同区域对树木干曲线有显著影响,削度方程参数估计值在不同区域的错误应用会导致较大误差。[结论] Kozak(2002)可变指数削度方程在拟合统计量、残差分布图和多重共线性等方面都表现出了一致性,预估精度达99%以上,可作为大兴安岭3个区域兴安落叶松的最优削度方程。
关键词：
可变指数削度方程 多重共线性 哑变量 非线性回归
Variable-Exponent Taper Models for Dahurian Larch in Different Regions of Daxing'anling
1. College of Forestry, Northeast Forestry University Harbin 150040;
2. College of Engineering and Technology, Northeast Forestry University Harbin 150040
2. College of Engineering and Technology, Northeast Forestry University Harbin 150040
Abstract: [Objective] The aim of this study was to develop taper functions of dahurian larch(Larix gmelinii)based on widely used variable-exponent taper models developed by Kozak(1988),Kozak(1994),Kozak(2001)and Kozak(2002),and the significances of different regions were also tested.[Method] Stem taper data of dahurian larch was collected from three regions of Daxing'anling. SUR method in SAS software was used to fit the four taper equations. Coefficient determination(Radj2), mean absolute bias(MAB), root mean square error(RMSE), mean percentage of bias(MPB)and prediction accuracy(P%)were employed to evaluate the precision of different models combining multicollinearity (CN)and graph of the residuals. The non-linear extra sum of squares method was used for regional comparison. The method required the fitting of full model and reduced model.[Result] 1)Variable-exponent taper models developed by Kozak(1988)and Kozak(1994)were better based on goodness of fit statistics and prediction precision, but they had high multicollinearity. The other models decreased multicollinearity and improved prediction precision of upper stem diameters. 2)F-test of regional comparison showed that stem taper in region 3 had large differences with region 1 and region 2. Stem taper in region 1 had small difference with region 2. Any two regional test showed significant difference(P<0.0001). It was indicated that any two regions could not use the same parameter estimates and separated parameter estimates were necessary. 3)Through the stem taper curve simulations in three regions, it could be seen that the same model showed the different taper curve simulations in three regions, especially the taper curves of model 3 showed significant difference in three regions. The taper curve simulations of other models showed small differences with region 1 and region 2, and significant difference was found for region 3. The results of simulations were consistent with the results from F-test. Different region had a significant effect on tree stem curve. Incorrectly application of taper models in different regions would result in larger prediction error.[Conclusion] Taper model developed by Kozak(2002)showed the consistent performance on fitting statistics, residual distribution, and multicollinearity. The prediction precision was more than 99%, and could be selected as the best taper equation for the three regions of Daxing'anling.
Key words:
variable-exponent taper models multicollinearity dummy variables nonlinear regression
干曲线的变化直接影响林木材积和材种出材量(孟宪宇，1982；胥辉等，1995；Li et al., 2010)，因此树干形状研究从干曲线方程到削度方程一直受到林学界的关注。削度是描述树干直径沿其树干向上随干径位置升高而逐渐减小变化程度的指标，削度方程是指树干上任意部位的直径为该位置距地面高、全树高及胸高直径的函数(Jiang et al., 2007；姜立春等，2011)。利用削度方程可以估计树干上任意部位的直径、任意既定直径部位距树基的长度、树干上任意分段的材积以及全树干的材积(Jiang et al., 2008； Yang et al., 2009)，而这些指标正是编制材种出材率表的重要依据。削度方程的形式有多种，按其发展顺序大致可分为简单削度方程、分段削度方程和可变指数削度方程三大类。林业专家最早用一个简单的数学函数来描述树干干形的变化(Kozak et al., 1969)，即简单削度方程，但是简单削度方程不能完整地表达树干形状。此后，研究者把树干分为了几部分，每部分都作为不同的几何体用数学函数来描述，一般把下部看作凹面体，中部看作抛物线，顶部则视为圆锥体，这就是分段函数的思想(Max et al., 1976； Leites et al., 2004； Jiang et al., 2005； Trincado et al., 2006)。近年来，一些学者提出了可变指数(variable exponent)模型，也称为变量形式(variable form)模型，通过变化的指数或变量描述干形从基部至树梢的连续变化，显示了更多的灵活性(Kozak，1988；Newnham，1992；Bi，2000； Sharma et al., 2004)。
本文以大兴安岭3个区域的兴安落叶松(Larix gmelinii)天然林为研究对象，选择Kozak(1988)，Kozak(1994)，Kozak(2001)和Kozak(2002)4个可变指数削度方程，采用SAS软件的SUR(seemingly unrelated regression)法对方程进行拟合，利用常用的统计检验指标调整确定系数(Radj2)、平均误差(MAB)、均方根误差(RMSE)、相对误差(MPB)和预估精度(P%)比较方程对不同区域树干任意处直径预测精度的高低，并结合削度方程的多重共线性指标和残差分布图等对方程进行综合比较分析。
1 数据与方法 1.1 数据本文所用数据为2004年12月—2005年2月在大兴安岭伊勒呼里山北坡西北部立地亚区(西林吉、图强、阿木尔、呼中林业局，漠河林场)、伊勒呼里山北坡东南部立地亚区(塔河、新林、十八站、韩家园、呼玛县林业局，二十二站林场)和大兴安岭北部东坡立地亚区(松岭、加格达奇林业局)采集的兴安落叶松样本数据。树木伐倒后，测量胸径、树高、冠幅、树冠高度、第一活枝高和第一死枝高，并以1m长为区分段进行干形测量。共收集了大兴安岭3个区域的样木分别为344，40和626株。将所收集的全部干形数据按75%和25%的比例分成建模数据样本和独立检验样本。兴安落叶松人工林各样木测树因子的统计量见表 1。
表1 大兴安岭3 个区域兴安落叶松各样木调查因子统计
Tab.1 Descriptive statistics for dahurian larch sample trees in three regions of Daxing’anling
|
选用Kozak(1988)，Kozak(1994)，Kozak(2001)和Kozak(2002)4个削度方程作为备选方程，具体形式如下：
Kozak(1988) d=b0Db1b2D
| $${\left( {\frac{{1 - \sqrt q }}{{1 - \sqrt p }}} \right)^{[{b_3}{q^2} + {b_4}\ln \left( {q + 0.001} \right) + {b_5}\sqrt q + {b_6}{e^q} + {b_7}D/H]}};$$ | (1) |
| $${\left( {\frac{{1 - \sqrt q }}{{1 - \sqrt p }}} \right)^{\left\{ {{b_3} + {b_4}{q^{1/4}} + {b_5}{q^{1/3}} + {b_6}{q^{1/2}} + {b_7}\arcsin (1 - {q^{1/2}}) + {b_8}\left( {1/\left[ {D/H + q} \right)} \right] + {b_9}H} \right\}}};$$ | (2) |
| $${\left( {\frac{{1 - {q^{1/4}}}}{{1 - {m^{1/4}}}}} \right)^{[{b_2} + {b_3}(1/{e^{D/H}}) + {b_4}{D^{(\frac{{1 - {q^{1/4}}}}{{1 - {m^{1/4}}}})}} + b5{{(\frac{{1 - {q^{1/4}}}}{{1 - {m^{1/4}}}})}^{D/H}}]}};$$ | (3) |
| $${\left( {\frac{{1 - {q^{1/3}}}}{{1 - {t^{1/3}}}}} \right)^{[{b_3}{q^4} + {b_4}(1/{e^{D/H}}) + {b_5}{{(\frac{{1 - {q^{1/3}}}}{{1 - {t^{1/3}}}})}^{0.1}} + {b_6}(1/D) + {b_7}{H^1} - {q^{1/3}} + {b_8}(\frac{{1 - {q^{1/3}}}}{{1 - {t^{1/3}}}})]}} \circ $$ | (4) |
模型(1)是Kozak(1988)提出的，可以对不同树种的材积进行无偏估计。模型(2)，(3)和(4)在模型(1)的基础上进行了改进，主要是为了降低模型(1)的多重共线性，提高树干上部的预测能力。
1.2.2 模型评价和检验指标
采用SAS软件对模型进行拟合得到参数估计值。拟合结果采用均方根误差(RMSE)和调整确定系数(Radj2)进行评价； 检验结果通过平均误差(MAB)、均方根误差(RMSE)、相对误差(MPB)和预估精度(P%)来对比不同模型。其相应的数学表达式为：
| $$R_{adj}^2 = 1 - \frac{{n - 1}}{{n - p - 1}}(1 - {R^2});$$ |
| $${R^2} = 1 - \left[ {\frac{{\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}} }}{{\sum\limits_{i = 1}^n {{{({y_i} - {{\bar y}_i})}^2}} }}} \right];$$ |
| $${\text{MAB}} = \frac{{\sum\limits_{i = 1}^n {\left| {{y_i} - {{\hat y}_i}} \right|} }}{n};$$ |
| $${\text{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}} }}{{n - 1}};} $$ |
| $${\text{MPE}} = 100 \times \sqrt {\frac{{\sum\limits_{i = 1}^n {\left| {{y_i} - {{\hat y}_i}} \right|} }}{{\sum\limits_{i = 1}^n {{y_i}} }};} $$ |
| $$P\% = (1 - {{{t_{0.05}}{{\rm{S}}_{{\rm{\bar y}}}}} \over {}}{\rm{)}} \times {\rm{100}}\% $$ |
区域性检验是为了比较一个区域的参数估计是否适合于另一个区域。如果没有显著不同，说明2个区域可以共用一套参数； 否则，就必须分别得到参数估计。区域性检验主要采用林业上常用的非线性额外平方和方法(Neter et al., 1996； Huang et al., 2000)，该方法需要拟合完整模型(full model)和简化模型(reduced model)，完整模型是利用哑变量方法引入区域参数，而简化模型则是所有区域共用一套参数。以模型(1)为例，用哑变量方法构造区域1和区域2的完整模型如下：
d=(b0+b00×r1)D(b1+b11×r1)(b2+b22×r1)D
| $${\left( {\frac{{1 - \sqrt q }}{{1 - \sqrt p }}} \right)^{\left[ {\matrix {({b_3} + {b_{33}} \times {r_1}){q^2} + ({b_4} + {b_{44}} \times {r_1}){\text{ln}}(q + 0.001) + } \\ {({b_5} + {b_{55}} \times {r_1})\sqrt q + ({b_6} + {b_{66}} \times r1){e^q} + ({b_{7 + }}{b_{77}} \times {r_1})D/H}} \right]}} \circ $$ | (5) |
零假设和备择假设如下：
H0： b00=b11=b22=b33=b44=b55=b66=b77=0；
Ha： 至少有一个参数不等于0。
F检验公式如下：
| $$F = \frac{{({\text{SS}}{{\text{E}}_{\text{R}}} - {\text{SS}}{{\text{E}}_{\text{F}}})/({\text{d}}{{\text{f}}_{\text{R}}} - {\text{d}}{{\text{f}}_{\text{F}}})}}{{{\text{SS}}{{\text{E}}_{\text{F}}}/{\text{d}}{{\text{f}}_{\text{F}}}}} \circ $$ | (6) |
如果F>F1－α;dfR－dfF，dfF，则拒绝零假设，即模型在不同区域间有显著不同； 如果F≤F1－α;dfR－dfF，dfF，则接受零假设，即模型在不同区域间没有显著不同。一般利用F值直接计算P值，P=FDIST(1－α;dfR－dfF，dfF)，若P<0.05，则区域间有显著不同。
2 结果与分析 2.1 模型的参数估计利用SAS统计软件，采用非线性回归方程分别拟合4个备选模型(1)～(4)，结果见表 2。在拟合过程中发现一些参数不显著，说明这些参数对模型没有显著影响，因此对模型进行修正，即假定模型中不显著的参数为0，再重新进行拟合，最终使各方程中的参数都显著(P<0.05)。
表2 大兴安岭3个区域兴安落叶松削度方程的参数估计值①
Tab.2 Parameter estimates of taper models for dahurian larch in three regions of Daxing’anling
|
|
表 3给出了不同模型拟合的统计量： 均方根误差(RMSE)、调整确定系数(Radj2)和多重共线性诊断指标条件数(condition number，CN)。Belsey(1991)指出，当CN＜30时，共线性不是主要问题； 当30＜CN＜100时，共线性存在，但模型可以接受； 当CN＞100时，模型存在严重的共线性。从表 3的评价指标可以看出，模型(2)的拟合精度最高，但是3个区域中模型的条件数分别为3 331，3 113和21 151，说明模型(2)存在严重的多重共线性；模型(1)的精度在4个模型中排在第2位，但是3个区域中模型的条件数分别为477.25，478.71和460.24，都大于100，说明模型(1)也存在严重的多重共线性；模型(4)的拟合精度优于模型(3)，虽然模型(4)的条件数大于模型(3)，但都在可接受的范围内。
表3 大兴安岭3个区域兴安落叶松削度方程的拟合统计量
Tab.3 Goodness-of-fit statistics of taper models for dahurian larch in three regions of Daxing’anling
|
|
模型的总体评价只能反映模型总体的拟合效果，不能反映模型的异方差性。为了全面评价各模型的拟合效果，分区域绘制各模型的残差分布图及趋势线(图 1～3)。可以看出，模型(1)，(2)和(4)在3个区域具有较高的等方差性和无偏性，模型(3)在3个区域都表现出较差的等方差性和无偏性。
|
图1 区域1各模型的残差分布 Fig.1 Graph of the residuals for taper models in region |
|
图2 区域2各模型的残差分布 Fig.2 Graph of the residuals for taper models in region |
|
图3 区域3各模型的残差分布 Fig.3 Graph of the residuals for taper models in region |
利用各区域的检验数据，基于表 2的参数估计值，采用SAS软件计算各模型的平均误差(MAB)、均方根误差(RMSE)、相对误差(MPB)和预估精度(P%)的统计量，结果见表 4。可以看出，各模型对树干不同高度处直径的预估精度都达到了99%以上，都可以很好地预估大兴安岭不同区域兴安落叶松树干不同高度处的直径。从偏差统计量来看，与拟合数据基本一致，模型(1)和(2)的偏差统计量较低，模型(4)的偏差统计量优于模型(3)。
表4 大兴安岭3个区域兴安落叶松削度方程的独立性检验
Tab.4 Validation for taper models for dahurian larch in three regions of Daxing’anling
|
|
将3个区域的数据两两进行合并，用合并后的数据分别拟合4个模型，利用式(6)计算F值及相应的P值，对比F检验结果见表 5。从模型(1)的检验结果可以看出，区域1和区域2有显著不同(F=8.94，P<0.000 1)，区域1和区域3有显著不同(F=51.19，P<0.000 1)，区域2和区域3有显著不同(F=56.75，P<0.000 1)。模型(2)，(3)和(4)的区域性检验都得到了相似的结果。针对每个模型的F检验分析发现，任何2个区域的对比都是显著的(P<0.000 1)。各模型对比检验F值的波动范围分别为5.71～8.94(区域1和区域2)、38.94～51.19(区域1和区域3)、45.42～56.75(区域2和区域3)，可以看出，区域3与区域1和区域2的树干削度相差较大，区域1和区域2的树干削度相差较小。
表5 不同区域模型对比F检验
Tab.5 F-test for models for different regions
|
|
为了更直观地分析不同区域兴安落叶松树干干形的变化，利用模型(1)～(4)进行模拟。参数估计值采用表 2中各区域的参数估计值，模拟树木的胸径假定为26 cm，树高为18 cm。图 4给出了4个模型在不同区域的干曲线模拟结果，可以看出，同一模型在3个区域的干曲线模拟结果是不同的，尤其是模型(3)在3个区域中的干曲线更明显地体现出3个区域的不同。模型(1)，(2)和(4)在3个区域的干曲线模拟体现在 区域1和区域2的模拟结果比较接近，区域3则与它们有明显不同。总体而言，不同区域对树木干曲线有一定影响，在选择削度方程时要区别对待。
|
图4 3个区域各模型干曲线模拟 Fig.4 Stem taper curve simulation for models in three regions |
本研究在大兴安岭3个区域分别选取样木测定干形数据，选取4个可变指数削度方程Kozak(1988)，Kozak(1994)，Kozak(2001)和Kozak(2002)，利用SAS统计软件中的非线性拟合方法拟合方程及其修正式，结果表明，虽然Kozak(1988)和Kozak(1994)削度模型拟合精度较高，但是这2个模型存在严重的多重共线性，与国外文献的研究结果一致(Kozak，2004； Rojo et al., 2005； Corral-Rivas et al., 2007)。残差分布图表明，模型(3)在3个区域都表现出了较差的等方差性和无偏性。模型(4)在拟合统计量和残差分布图上都表现了一致性，且预估精度达到99%以上，适合于大兴安岭兴安落叶松干形的预测。模型在不同区域之间的F检验表明，4个模型在3个区域检验中的P值都小于0.000 1。各模型对比检验F值的波动范围表明，区域2和区域3的树干削度相差较大，其次是区域1和区域3，区域1和区域2的树干削度相差较小。大兴安岭不同区域由于气候(温度和降水量)、土壤潮湿度、植被演替和类型等方面的不同导致了树木干形的不同，环境条件不仅影响树木的生长，也影响树木的干形和木材产量。由于数据没有林分变量，本研究没有探讨林分因子对干形的影响，随着数据的积累，这方面将进一步深入研究。
参考文献(References)
| [1] |
姜立春, 刘瑞龙. 2011. 基于非线性混合模型的落叶松树干削度模型. 林业科学, 47(4):101-106.. (Jiang L C,Liu R L.2011. A stem taper model with nonlinear mixed effects for dahurian larch. Scientia Silvae Sinicae,47(4):101-106.[in Chinese])(  1)
|
| [2] |
孟宪宇. 1982. 削度方程和出材率表的研究. 南京林产工业学院学报,(1):122-133.. (Meng X Y.1982. Studies of taper equations and the table of merchantable volumes. Journal of Nanjing Technological College of Forest Products,(1):122-133.[in Chinese])( 1)
|
| [3] |
胥辉, 孟宪宇. 1995. 天山云杉削度方程与材种出材率表的研究. 北京林业大学学报, 18(3):21-30.. (Xu H,Meng X Y.1995. Study on taper function and merchantable volume yielding rate tables of Picea schrenkiana var. tianshanica. Journal of Beijing Forestry University,18(3):21-30.[in Chinese])( 1)
|
| [4] |
Belsey D A. 1991. Conditioning diagnostics, collinearity and weak data in regression. John Wiley & Sons Inc., New York.(1)
|
| [5] |
Bi H. 2000. Taper equation for second-growth mixed conifers of northern California. For Sci, 46(3):397-409.(1)
|
| [6] |
Corral-Rivas J J, Diéguez-Aranda U, Rivas S C, et al. 2007. A merchantable volume system for major pine species in EI Salto, Durango(Mexico). For Ecol Manage, 238(1/3):118-129.(1)
|
| [7] |
Huang S, Price D, Titus S J. 2000. Development of ecological based height-diameter models for white spruce in boreal forests. For Ecol Manage, 129(1/3):125-141.(1)
|
| [8] |
Jiang L, Brooks J R. 2008. Taper, volume, and weight equations for red pine in West Virginia. North J Appl For, 25(3):151-153.(1)
|
| [9] |
Jiang L, Brooks J R, Hobbs G R. 2007. Using crown ratio in yellow-poplar compatible taper and volume equations. North J Appl For, 24(4):271-275.(1)
|
| [10] |
Jiang L, Brooks J R, Wang J. 2005. Compatible taper and volume equations for yellow-poplar in West Virginia. For Ecol Manage, 213(1/3):399-409.(1)
|
| [11] |
Kozak A, Munro D O, Smith J H G. 1969. Taper functions and their application in forest inventory. For Chron, 45(4):278-283.(1)
|
| [12] |
Kozak A.1988. A variable-exponent taper equation. Can J For Res, 18(11):1363-1368.(5)
|
| [13] |
Kozak A. 2004. My last word on taper functions. For Chron, 80(4):507-515.(1)
|
| [14] |
Leites L P, Robinson A P. 2004. Improving taper equations of loblolly pine with crown dimensions in mixed-effects modeling framework. For Sci, 50(2):204-212.(1)
|
| [15] |
Li R, Weiskittel A R. 2010. Comparison of model forms for estimating stem taper and volume in the primary conifer species of the North American Acadian region. Ann For Sci,67(3):1-16.(1)
|
| [16] |
Max T A, Burkhart H E. 1976. Segmented polynomial regression applied to taper equations. For Sci, 22(3):283-289.(1)
|
| [17] |
Neter J, Kutner M H, Nachtsheim C J, et al. 1996. Applied linear statistical models. New York:McGraw-Hill, 1048.(1)
|
| [18] |
Newnham R M. 1992. Variable-form taper functions for four Alberta tree species. Can J For Res, 22(2):210-223.(1)
|
| [19] |
Rojo A, Perales X, Sanchez-Rodriguez F, et al. 2005. Stem taper functions for maritime pine(Pinus pinaster Ait.)in Galicia(Northwestern Spain). Eur J For Res, 124(3):177-186.(1)
|
| [20] |
Sharma M, Zhang S Y. 2004. Variable-exponent taper equations for jack pine, black spruce, and balsam fir in eastern Canada. For Ecol Manage, 198(1/3):39-53.(1)
|
| [21] |
Trincado G, Burkhart H E. 2006. A generalized approach for modeling and localizing stem profiles curves. For Sci, 52(6):670-682.(1)
|
| [22] |
Yang Y, Huang S, Trincado G, et al. 2009.Nonlinear mixed-effects modeling of variable-exponent taper equations for lodgepole pine in Alberta, Canada. Eur J For Res, 128(4):415-429.(1)
|