2016, Vol. 52
文章信息
- 姜立春, 马英莉, 李耀翔
- 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
-
作者相关文章
2. 东北林业大学工程技术学院 哈尔滨 150040
2. College of Engineering and Technology, Northeast Forestry University Harbin 150040
干曲线的变化直接影响林木材积和材种出材量(孟宪宇，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。
|
选用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)。
|
|
表 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)，但都在可接受的范围内。
|
|
模型的总体评价只能反映模型总体的拟合效果，不能反映模型的异方差性。为了全面评价各模型的拟合效果，分区域绘制各模型的残差分布图及趋势线(图 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)。
|
|
将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的树干削度相差较小。
|
|
为了更直观地分析不同区域兴安落叶松树干干形的变化，利用模型(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的树干削度相差较小。大兴安岭不同区域由于气候(温度和降水量)、土壤潮湿度、植被演替和类型等方面的不同导致了树木干形的不同，环境条件不仅影响树木的生长，也影响树木的干形和木材产量。由于数据没有林分变量，本研究没有探讨林分因子对干形的影响，随着数据的积累，这方面将进一步深入研究。
| [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)
|