-
削度方程能预测树干总材积，还能预测任意给定小头直径或一定高度的商品材积[1]。在森林经营规划、木材收获与利用、估算出材率和树干生物量等方面有着广泛的应用[2–5]。因此，削度方程是林业生产活动中重要的工具。由于树干的削度变化受树种[6]、遗传[7]、营林措施[8–10]、气候[11]以及立地条件[12]等多种因素的影响，目前没有任何一个模型可以适用于所有树种或不同区域的同一树种的干形的描述[13]。当前大多数削度模型的构建都是基于参数法[14–16]，如树高、胸径及树干不同高度构建的削度方程。一些学者在削度模型构建中尝试引入其他林分变量或气候变量，但有研究表明加入的额外变量并没有使削度模型的精度显著提高[17]，另外，参数模型的构建经常要考虑方差的同质性、多重共线性以及空间自相关性等假设前提，这增加了模型构建和拟合的难度。
随着广义加性模型理论和计算机技术的发展，广义加性模型在林业上得到了一定的应用[18–24]。广义加性模型对模型的前提假设要求不严，对数据的分布没有要求，主要通过数据来驱动模型，具有较好的稳健性[18]。因此，广义加性模型是对参数模型的一种扩展，具有较强的适用性和适应性。目前国内利用广义加性模型研究削度方程以及与参数削度方程的比较还未见报道，因此本研究以大兴安岭樟子松为研究对象，采用广义加性模型构建削度方程，基于R软件mgcv包中gamm函数的6个常见光滑样条函数来拟合削度方程，并与传统的3个参数削度方程曾伟生等（1997）[25]，Bi（2000）[26]和Kozak（2004）[1]进行比较。该研究对探索森林经营管理预测的新方法具有重要的理论和现实意义。
-
数据来源于黑龙江省大兴安岭地区呼玛县韩家园林业局（125° 09′ ~ 127°10′ E，50°45′ ~ 52°19′ N）。林业局地处大兴安岭东南坡、小兴安岭过渡带的西北部，境内平均海拔350 m。属寒温带大陆性气候，冬季严寒、漫长，夏季温暖，春秋两季短、持续大约6周。1月平均气温−26.2 ℃，7月平均气温20.6 ℃，年平均气温−1.01 ℃。全年降水量为471 mm，全年无霜期100 d左右。
本研究选择干形较好的樟子松天然林样木，在样木伐倒前对其胸径进行测量，伐倒后使用皮尺测量树干实际高度。在每株样木的相对树高为0、0.02、0.04、0.06、0.08、0.10、0.15、0.20、0.30、0.40、0.50、0.60、0.70、0.80、0.90处测量带皮直径，共收集180株樟子松样木的数据。所收集样木测树因子统计量如表1所示。
表 1 样木调查因子统计量
Table 1. Descriptive statistics for sample trees
变量 Variable 样木数 Sample tree number 最小值 Min. value 最大值 Max. value 平均值 Mean 标准差 SD 变异系数 CV 胸径 DBH/cm 180 5.2 54 28.48 10.42 36.60 树高 Tree height/m 180 6.2 25.4 18.93 4.09 21.61 在数据收集过程中，树干上有些部位有大的结节或者凸起，通过画相对树高和相对直径散点图能直观发现异常点，可以逐个剔除，但是这种方法费时费力。Bi（2000）采用非参数LOESS方法拟合相对树高和相对直径并系统剔除异常值的方法，该方法能快速有效地剔除异常数据，提高工作效率[26]。本研究使用Bi（2000）中的方法剔除异常点72个，占总数据的2.67%，剔除的数据不再计入后续的计算中。图1表示数据处理前（a）和处理后（b）的相对树高和相对直径散点图。
图 1 非参数回归LOESS曲线剔除异常值前后相对树高和相对直径散点图
Figure 1. Scatter plots of relative tree height and relative diameter before and after eliminating outliers with nonparametric regression LOESS curve
-
广义加性模型（generalized additive model）由不同的加性项组成，每个加性项使用单个光滑样条函数来估计，每一加性项中可以解释因变量如何随自变量变化而变化。其结构形式如下[27]：
$${Y_i} = \alpha + {S\!_1}\left( {{x_1}} \right) + {S\!_2}\left( {{x_2}} \right) + \cdots + {S\!_n}\left( {{x_n}} \right) + \varepsilon $$ (1) 式中：Yi表示因变量，α是截距，x1，···，xn表示自变量，n为变量个数，S1，···，Sn表示非参数光滑样条函数，ε表示随机误差项。
广义加性模型的拟合主要考虑使用不同的光滑样条函数，本研究利用R语言mgcv软件包gamm函数[27]，选择6个常见的光滑样条函数对广义加性削度模型进行拟合。所选样条函数有B-样条函数（B-spline function，BS），三次回归样条函数（cubic regression spline function，CR），Duchon样条函数（Duchon spline function，DS），高斯过程平滑样条函数（Gaussian process smooth spline function，GP），P-样条函数（P-spline function，PS），薄板回归样条函数（thin plate regression spline function，TP）[27-28]。
-
为了与广义加性削度方程比较，本研究选用国内外林业上拟合精度较高的3个变指数削度方程[1, 25-26]进行比较，模型具体形式如下：
曾伟生等（1997）提出的削度方程[25]：
$$\frac{d}{D} = {\left( {\frac{{H - h}}{{H - 1.3}}} \right)^{{b_1} + {b_2}{T^{1/4}} + {b_3}{T^{1/2}} + {b_4}\frac{\scriptstyle{D}}{\scriptstyle{H}}}}$$ (2) Bi（2000）提出的削度方程[26]：
$$ d = D{\left( {\frac{{\ln \sin \left(\dfrac{{\text{π}} }{2}T\right)}}{{\ln \sin \left( {\dfrac{{\text{π}} }{2} \cdot \dfrac{{1.3}}{H}} \right)}}} \right)^{{b_1} + {b_2}\sin \left( {\frac{\scriptstyle{\text{π}} }{\scriptstyle{2}}T} \right) + {b_3}\cos \left( {\frac{{3{\scriptstyle{\text{π}}} }}{\scriptstyle{2}}T} \right) + {b_4}\sin \left( {\frac{\scriptstyle{\text{π}} }{\scriptstyle{2}}T} \right)/T + {b_5}D + {b_6}T\sqrt D + {b_7}T\sqrt H }}$$ (3) Kozak（2004）提出的削度方程[1]：
$$d = {b_1}{D^{{b_2}}}{H^{{b_3}}}{x_2}^{{b_4}{T^4} + {b_5}\left( {\frac{\scriptstyle{1}}{\scriptstyle{{{\rm{e}}^{D/H}}}}} \right) + {b_6}{x_2}^{0.1} + {b_7}\left( {\frac{\scriptstyle{1}}{\scriptstyle{D}}} \right) + {b_8}{H^Q} + {b_9}{x_2}}$$ (4) 式中：
$T = \dfrac{h}{H}$ ，${x_2} = \dfrac{{1 - {{\left( {h/H} \right)}^{1/3}}}}{{1 - {{\left( {{{1.3} / H}} \right)}^{1/3}}}}$ ，$Q = 1 - \sqrt[3]{T}$ ，d表示不同部位的直径，cm；h表示树干某一位置处高度，m；D表示胸径，cm；H表示树高，m；bi为参数。 -
评价参数模型和广义加性模型的统计量有：平均绝对误差（mean absolute bias，MAB），相对百分误差（relative percentage bias，MPB），均方根误差（root mean square error，RMSE），确定系数（coefficient of determination，R2）。模型检验使用留一交叉检验法（leave-one-out cross-validation）[29]，检验结果通过MAB、MPB和RMSE来比较。各统计量数学表达式如下：
$${\rm{MAB}} = \frac{1}{n}\mathop \sum \limits_{i = 1}^n \left| {{y_i} - {{\hat y}_i}} \right|$$ (5) $${\rm{MPB}} = \frac{{\displaystyle\sum\limits_{i = 1}^n {\left| {{y_i} - \hat y} \right|} }}{{\displaystyle\sum\limits_{i = 1}^n {{y_i}} }} \times 100$$ (6) $${\rm{RMSE}} = \sqrt {\frac{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}{{n - p}}} $$ (7) $${R^{\rm{2}}} = 1 - \left[ {\frac{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{y_i} - \bar y} \right)}^2}} }}} \right]$$ (8) 式中：
${y_i}$ 为实测值，${\hat y_i}$ 为预测值，n为样本数，$\bar y$ 为实测值的均值，p表示模型参数个数。 -
广义加性削度方程的变量选取与参数模型一样，都是选用不同部位直径d或相对直径d/D作为因变量，自变量选用树高H、胸径D和相对树高h/H等变量。同时考虑变量否需要进行变量转换，如对数、平方、开根号等。经过初步筛选和分析本文构建如下广义加性削度模型：
$$\frac{d}{D} = \alpha + s\left( {{D^2}} \right) + s\left( H \right) + s\left( {\sqrt {\frac{h}{H}} } \right)$$ (9) 式中：α表示截距，s(·)表示光滑样条函数，其他变量定义如上。
-
利用R软件mgcv包的gamm函数对模型（9）进行拟合，在模型拟合时分别考虑6种光滑样条函数（BS，CR，DS，GP，PS，TP），广义加性模型由含有截距的参数项和含有平滑项的非参数项组成。通过拟合，结果如表2所示。可以看出，影响拟合结果的变量有胸径的平方，树高，相对高度的算数平方根。通过对各模型的非参数项进行显著性检验发现P值均小于0.001，这说明这些变量对于树干干形变化有显著的意义。
表 2 广义加性模型拟合结果
Table 2. Fitting results of generalized additive models
模型 Model 参数项 Parametric term 非参数项 Nonparametric term 光滑样条
Smooth spline截距
Intercept估计值
Estimated value标准误
SEt值
t valueP值
P value平滑项
Smooth term自由度
dfF值
F valueP值
P valueB-样条函数
B-spline function (BS)${\alpha _1}$ 0.796 4 0.000 8 933.7 < 0.001 ${s_1}({D^2})$ 4.449 28.46 < 0.001 ${s_1}(H)$ 8.646 12.16 < 0.001 ${s_1}(\sqrt {h/H} )$ 8.040 12 645.02 < 0.001 三次回归样条函数
Cubic regression spline function (CR)${\alpha _2}$ 0.796 4 0.000 9 816.8 < 0.001 ${s_2}({D^2})$ 5.297 17.497 < 0.001 ${s_2}(H)$ 6.384 7.763 < 0.001 ${s_2}(\sqrt {h/H} )$ 3.991 19 245.985 < 0.001 Duchon样条函数
Duchon spline function (DS)${\alpha _3}$ 0.796 4 0.000 8 935.7 < 0.001 ${s_3}({D^2})$ 5.378 22.83 < 0.001 ${s_3}(H)$ 11.078 11.35 < 0.001 ${s_3}(\sqrt {h/H} )$ 10.702 6 285.95 < 0.001 高斯过程平滑样条函数
Gaussian process smooth spline function (GP)${\alpha _4}$ 0.796 4 0.000 8 927.7 < 0.001 ${s_4}({D^2})$ 5.243 23.06 < 0.001 ${s_4}(H)$ 7.178 11.19 < 0.001 ${s_4}(\sqrt {h/H} )$ 8.788 11 422.21 < 0.001 P-样条函数
P-spline function (PS)${\alpha _5}$ 0.796 4 0.000 8 924.6 < 0.001 ${s_5}({D^2})$ 4.053 27.25 < 0.001 ${s_5}(H)$ 5.220 13.32 < 0.001 ${s_5}(\sqrt {h/H} )$ 7.762 12 843.74 < 0.001 薄板回归样条函数
Thin plate regression spline function (TP)${\alpha _6}$ 0.796 4 0.000 8 928.3 < 0.001 ${s_6}({D^2})$ 5.183 23.51 < 0.001 ${s_6}(H)$ 7.388 11.34 < 0.001 ${s_6}(\sqrt {h/H} )$ 8.751 11 485.39 < 0.001 注：α1 ~ α6分别表示BS、CR、DS、GP、PS、TP对应的截距，s1 ~ s6分别表示BS、CR、DS、GP、PS、TP对应的光滑样条函数。
Notes: α1−α6 represent the intercept of BS, CR, DS, GP, PS and TP, respectively. s1−s6 represent smooth splines of BS, CR, DS, GP, PS and TP, respectively.此外，广义加性模型的总体拟合效果也可以从偏残差（partial residuals）分布图看出。图2给出了广义加性模型各组成的偏残差图，第i个加性项的偏残差指当前模型的因变量减去不包含该加性项的其他加性项的估计值所得的残差，如果每个加性项的偏残差均匀分布在图中曲线的周围表明模型拟合效果较好[27]。从图2可以看出，所有光滑样条函数不同加性项的偏残差分布均匀，且具有相似的趋势。广义加性模型的相对直径拟合值随着相对树高算数平方根的增大而减小，且表现出较强的非线性关系，这符合树木干形的变化规律，即相对树高值越大，相对直径值越小，这与表2中相对高度的平方根的F检验相一致，即与其他2个变量的F值相比，相对高度平方根的F值都较大。将相对高度的算数平方根换算为相对高度，在相对高度为0.04 ~ 0.4的范围内，偏残差分布的范围明显变窄，这说明在这个范围内模型的预测效果较好。虽然胸径的平方和树高与干形表现出的非线性关系不强，但对这两个变量的显著性检验发现P值均小于0.001（表2），所以这两个变量在建模中作为自变量也具有显著意义。
图 2 广义加性模型各组成的偏残差
Figure 2. Partial residuals of components for generalized additive models
-
通过调用R软件nlme软件包对3个参数削度方程进行拟合，参数估计结果见表3。曾伟生等（1997）和Bi（2000）模型都有不显著的参数（P > 0.05），如果模型删除不显著的参数重新拟合会轻微增大RMSE值，因此本研究保留不显著参数。
表 3 模型参数估计
Table 3. Parameter estimate of models
模型 Model b1 b2 b3 b4 b5 b6 b7 b8 b9 曾伟生等(1997)
Zeng et al.(1997)3.068 2
(0.049 3)−5.225 5
(0.125 6)2.837 7
(0.088 1)−0.006 5
(0.007 7)nsBi(2000) −0.673 7
(0.051 2)0.532 6
(0.029 1)0.151 8
(0.006 6)0.413 1
(0.031 1)−0.000 7
(0.000 1)0.012 7
(0.003 6)0.008 6
(0.007 1)nsKozak(2004) 0.847 3
(0.024 0)0.997 1
(0.004 3)0.057 4
(0.010 7)0.480 0
(0.014 3)0.237 1
(0.037 9)0.364 7
(0.008 9)−2.088 8
(0.248 3)−0.012 3
(0.001 4)0.164 9
(0.020 3)注：括号内为标准误差，ns 表示参数估计在0.05水平不显著。Notes: standard error is showed in bracket, ns represents that parameter is not significant at 0.05 level. -
为了比较模型预测效果，本研究使用留一交叉检验法对各模型进行检验。结果如表4所示。可以看出，参数模型中依然是Kozak（2004）模型预测效果较好，广义加性模型中除CR外，其他样条函数预测效果相似，且均优于参数模型，其中BS光滑样条的预测效果最好。
表 4 参数模型和广义加性模型交叉检验统计量
Table 4. Cross-validation for parametric models and generalized additive models
模型
Model平均绝对误差
Mean absolute
bias (MAB)/cm相对百分误差
Relative
percentage
bias (MPB)均方根误差
Root mean
square error
(RMSE)/cm曾伟生等(1997)
Zeng et al. (1997)1.048 8 4.956 9 1.474 1 Bi (2000) 1.018 6 4.814 2 1.450 8 Kozak (2004) 0.957 4 4.525 1 1.361 0 BS 0.911 9 4.309 8 1.310 9 CR 1.040 0 4.915 4 1.468 5 DS 0.923 3 4.364 1 1.331 8 GP 0.915 7 4.327 9 1.318 3 PS 0.923 8 4.366 2 1.323 9 TP 0.916 0 4.329 6 1.317 7 为了进一步分析3个参数模型和BS样条在预测树干不同高度处直径的精度，图3给出了不同模型MPB和RMSE随相对高度变化。通过MPB的变化情况发现，BS模型在预测树干不同位置的精度高于3个参数模型，特别是随着相对高度的增加，BS预测精度更高。通过比较RMSE随着相对高度的变化也发现BS在预测不同高度时的精度高于其他模型。通过图3说明BS模型在预测树干不同部位的精度优于参数模型。
图 3 模型相对百分误差和均方根误差随相对树高变化
Figure 3. Changes of MPB and RMSE with relative tree height for models
为了更直观地比较广义加性模型和变指数模型，本研究选择广义加性模型中预测效果最好的BS和参数模型中预测效果最好的Kozak（2004）模型对不同大小的树干进行模拟。从已知数据中分别选择胸径为9.2 cm、树高10.1 m的一株小树和胸径为43.6 cm、树高23.6 m的一株大树进行模拟如图4所示。图4a表示对小树的模拟效果，图4b表示对大树的模拟效果。可以看出对于小树，BS模型的模拟效果更好，而Kozak（2004）模型随着高度的增加偏离了实际值。通过对大树进行模拟发现Kozak（2004）和BS的模拟效果相似，都能有效模拟大树的干形变化。
图 4 削度模型对不同大小树干的模拟
Figure 4. Estimated stem taper equation for trees with different sizes
-
本研究以大兴安岭地区樟子松天然林为研究对象，用构建削度方程常用的变量胸径D、树高H和不同部位高度h构建了估计h对应位置直径d的广义加性模型。该削度方程的构建考虑了不同变量的变换形式，选取不同的变量形式对模型的预测精度影响较大，特别是使用直径d和相对直径d/D时模型的拟合效果差异较大。当以d为因变量时，使用不同变量形式和不同光滑样条函数进行组合，发现无论哪种组合方式，模型的预测精度都没有提高，且有的组合方式会使d的预测值出现负数，这严重不符合干形的变化规律，所以将d作为因变量不予以考虑。设置D取平方，相对树高开根号为自变量时，模型精度最高，最终确定式（9）为基础广义加性模型。然后使用R软件mgcv软件包gamm函数的6个光滑样条函数（BS，CR，DS，GP，PS，TP）对模型进行拟合。当考虑使用不同的样条函数对每一个加性项进行组合，发现所有加性项使用同一光滑样条函数比使用不同样条函数组合的精度高，其中广义加性模型中CR模型拟合精度没有其他光滑样条函数高。本研究构建的削度形式和Robinson等[30]构建的加拿大6个树种的广义加性削度方程形式一致，但是他们只考虑了CR光滑样条函数。本研究比较了6种光滑样条函数，发现CR并不适合预测该地区樟子松的削度方程，同时将广义加性削度方程和3个拟合精度较高的变指数削度方程进行了比较，发现除CR外，其他广义加性模型的精度均高于参数削度方程，其中BS模型的预测精度最高。因此本研究最终确定的广义加性削度模型如式（10）：
$$\frac{d}{D} = 0.796 \; 4 + s\left( {{D^2}} \right) + s\left( H \right) + s\left( {\sqrt {\frac{h}{H}} } \right)$$ (10) 式中：各符号代表含义同式（9）。
为了更加准确估计树干的干形，林业工作者不断提出许多不同形式且越来越复杂的参数或分段削度方程，有时使得这些方程在使用中不太方便。而广义加性削度方程只需进行简单的编程就能进行树干干形的预测，没有复杂的模型构建过程，且精度能够达到要求。此外，广义加性削度方程通过简单编程也能实现树干材积的预测。在实际应用中，为了确保预测精度，所预测数据的范围尽量在拟合数据的分布范围内。
-
本研究构建了以相对直径为因变量，胸径的平方、树高、相对树高的算术平方根为自变量的樟子松树干削度方程，并使用gamm函数常见的6个光滑样条函数进行了拟合，通过对6个光滑样条函数（BS，CR，DS，GP，PS，TP）进行比较，发现除了CR外，其他光滑样条函数均能较好地描述樟子松干形的变化规律，同时和林业上精度较高的变指数削度方程曾伟生等（1997）、Bi（2000）和Kozak（2004）进行了比较，通过拟合和交叉检验发现本文所构建的广义加性模型均优于所选的3个参数削度方程，其中BS的精度最高。因此本研究所构建的广义加性削度方程BS可替代传统参数削度方程来预测该地区樟子松的干形和材积。
Research on stem taper equation of Scots pine based on generalized additive model
-
摘要:
目的 基于广义加性模型理论，构建樟子松的广义加性树干削度方程，并和林业上精度较高的变指数削度方程曾伟生等（1997）、Bi（2000）以及Kozak（2004）进行预测精度比较。 方法 以大兴安岭樟子松为研究对象，使用胸径、树高和不同部位高度及该部位树干直径及其变形构建广义加性削度方程，利用R软件mgcv软件包gamm函数对广义加性模型进行拟合，拟合过程中采用6种样条函数：B样条函数（BS）、三次回归样条函数（CR）、Duchon样条函数（DS）、高斯过程平滑样条函数（GP）、P样条函数（PS）和薄板回归样条函数（TP）。使用留一交叉检验法对模型进行检验。 结果 （1）将相对直径作为因变量，将胸径的平方、相对树高的算术平方根和树高作为自变量构建了最优的广义加性削度方程结构。（2）拟合结果表明，除CR外，其他光滑样条函数表现了相似的拟合结果，且均优于变指数削度方程的统计指标。（3）交叉检验结果表明，除CR光滑样条函数外，广义加性模型（BS，DS，GP，PS，TP）总体与拟合结果基本一致，即预测精度都优于曾伟生等（1997）、Bi（2000）和Kozak（2004）模型，其中广义加性模型中BS模型的预测精度最高，变指数削度方程中Kozak（2004）预测精度最高。（4）通过对比BS和Kozak（2004）模型的干曲线模拟发现，Kozak（2004）在预测小树树干上部时误差较大，而BS在模拟小树和大树上都具有较高的精度。 结论 广义加性模型是构建削度方程的一种非参数方法，基于BS样条函数的广义加性削度方程预测精度最高，适合大兴安岭地区樟子松的干形预测。 Abstract:Objective Based on the theory of generalized additive model, stem taper equation was constructed for Scots pine (Pinus sylvestris). Accurate variable exponent taper equations such as Zeng et al. (1997), Bi (2000) and Kozak (2004) in forestry were used for comparison. Method The generalized additive taper equation was constructed using DBH, tree height, the height at different stem parts, the diameter at different tree heights and their transformation based on taper data of Scots pine. The model was fitted using the gamm function in mgcv library of the R software. Six smooth splines were chosen for fitting process, i.e. B-spline function (BS), cubic regression spline function (CR), Duchon spline function (DS), Gaussian process smooth spline function (GP), P-spline function (PS) and thin plate regression spline function (TP). The models were validated using leave-one-out cross-validation method. Result (1) The optimal generalized additive model form of taper equation was constructed by response variable relative diameter and explanation variable square of DBH, total height and the square root of relative height. (2) The fitting results showed that the generalized additive taper equations were similar and better than parametric taper equation except for CR function. (3) The cross validation results showed that the generalized additive models (BS, DS, GP, PS and TP) were basically consistent with the fitting results except for CR, i.e. they were superior to parametric taper equation of Zeng et al. (1997), Bi (2000) and Kozak (2004). The BS model had the highest prediction accuracy in the generalized additive models. Kozak (2004) had the highest prediction accuracy in the variable exponential taper equations. (4) Through the simulation of stem curves of BS and Kozak (2004) models, it was found that Kozak (2004) had a large error in predicting the upper part of the stem of a small tree. However, BS had higher accuracy in simulating small tree and large tree. Conclusion The generalized additive model is a nonparametric method for constructing taper equation. The generalized additive taper equation based on BS spline has the highest prediction accuracy. It’s suitable for the prediction of the shape of Scots pine. -
Key words:
- generalized additive model /
- spline function /
- Scots pine /
- taper equation
-
图 1 非参数回归LOESS曲线剔除异常值前后相对树高和相对直径散点图
光滑系数 = 0.3 Smoothing coefficient = 0.3
Figure 1. Scatter plots of relative tree height and relative diameter before and after eliminating outliers with nonparametric regression LOESS curve
图 2 广义加性模型各组成的偏残差
在同一个光滑样条函数对应的偏残差图中，横坐标中sq.D表示胸径的平方，H表示树高，sqrt.rh表示相对高度的平方根，垂直于x轴的小线段表示变量的数值。y轴表示偏残差。纵坐标括号内的数值表示自由度。图中灰色区域表示95%的置信区间。 In the partial residual graph corresponding to the same smooth spline function, sq.D represents square of DBH, H represents tree height, sqrt.rh represents square root of relative tree height. The ‘rug plots’, along the x-axis bottom of each plot show the values of variable. y-axis represents partial residuals. The number in each y-axis caption is the degree of freedom (df) of the term being plotted. The grey area in the graph represents 95% confidence interval.
Figure 2. Partial residuals of components for generalized additive models
图 3 模型相对百分误差和均方根误差随相对树高变化
Figure 3. Changes of MPB and RMSE with relative tree height for models
图 4 削度模型对不同大小树干的模拟
小树、大树的胸径分别为9.2、43.6 cm，树高分别为10.1、23.6 m。The DBH for small tree and large tree are 9.2 and 43.6 cm, respectively. The tree height for them are 10.1 and 23.6 m, respectively.
Figure 4. Estimated stem taper equation for trees with different sizes
表 1 样木调查因子统计量
Table 1. Descriptive statistics for sample trees
变量 Variable 样木数 Sample tree number 最小值 Min. value 最大值 Max. value 平均值 Mean 标准差 SD 变异系数 CV 胸径 DBH/cm 180 5.2 54 28.48 10.42 36.60 树高 Tree height/m 180 6.2 25.4 18.93 4.09 21.61 
下载: 导出CSV
表 2 广义加性模型拟合结果
Table 2. Fitting results of generalized additive models
模型 Model 参数项 Parametric term 非参数项 Nonparametric term 光滑样条
Smooth spline截距
Intercept估计值
Estimated value标准误
SEt值
t valueP值
P value平滑项
Smooth term自由度
dfF值
F valueP值
P valueB-样条函数
B-spline function (BS)${\alpha _1}$ 0.796 4 0.000 8 933.7 < 0.001 ${s_1}({D^2})$ 4.449 28.46 < 0.001 ${s_1}(H)$ 8.646 12.16 < 0.001 ${s_1}(\sqrt {h/H} )$ 8.040 12 645.02 < 0.001 三次回归样条函数
Cubic regression spline function (CR)${\alpha _2}$ 0.796 4 0.000 9 816.8 < 0.001 ${s_2}({D^2})$ 5.297 17.497 < 0.001 ${s_2}(H)$ 6.384 7.763 < 0.001 ${s_2}(\sqrt {h/H} )$ 3.991 19 245.985 < 0.001 Duchon样条函数
Duchon spline function (DS)${\alpha _3}$ 0.796 4 0.000 8 935.7 < 0.001 ${s_3}({D^2})$ 5.378 22.83 < 0.001 ${s_3}(H)$ 11.078 11.35 < 0.001 ${s_3}(\sqrt {h/H} )$ 10.702 6 285.95 < 0.001 高斯过程平滑样条函数
Gaussian process smooth spline function (GP)${\alpha _4}$ 0.796 4 0.000 8 927.7 < 0.001 ${s_4}({D^2})$ 5.243 23.06 < 0.001 ${s_4}(H)$ 7.178 11.19 < 0.001 ${s_4}(\sqrt {h/H} )$ 8.788 11 422.21 < 0.001 P-样条函数
P-spline function (PS)${\alpha _5}$ 0.796 4 0.000 8 924.6 < 0.001 ${s_5}({D^2})$ 4.053 27.25 < 0.001 ${s_5}(H)$ 5.220 13.32 < 0.001 ${s_5}(\sqrt {h/H} )$ 7.762 12 843.74 < 0.001 薄板回归样条函数
Thin plate regression spline function (TP)${\alpha _6}$ 0.796 4 0.000 8 928.3 < 0.001 ${s_6}({D^2})$ 5.183 23.51 < 0.001 ${s_6}(H)$ 7.388 11.34 < 0.001 ${s_6}(\sqrt {h/H} )$ 8.751 11 485.39 < 0.001 注：α1 ~ α6分别表示BS、CR、DS、GP、PS、TP对应的截距，s1 ~ s6分别表示BS、CR、DS、GP、PS、TP对应的光滑样条函数。
Notes: α1−α6 represent the intercept of BS, CR, DS, GP, PS and TP, respectively. s1−s6 represent smooth splines of BS, CR, DS, GP, PS and TP, respectively.
下载: 导出CSV
表 3 模型参数估计
Table 3. Parameter estimate of models
模型 Model b1 b2 b3 b4 b5 b6 b7 b8 b9 曾伟生等(1997)
Zeng et al.(1997)3.068 2
(0.049 3)−5.225 5
(0.125 6)2.837 7
(0.088 1)−0.006 5
(0.007 7)nsBi(2000) −0.673 7
(0.051 2)0.532 6
(0.029 1)0.151 8
(0.006 6)0.413 1
(0.031 1)−0.000 7
(0.000 1)0.012 7
(0.003 6)0.008 6
(0.007 1)nsKozak(2004) 0.847 3
(0.024 0)0.997 1
(0.004 3)0.057 4
(0.010 7)0.480 0
(0.014 3)0.237 1
(0.037 9)0.364 7
(0.008 9)−2.088 8
(0.248 3)−0.012 3
(0.001 4)0.164 9
(0.020 3)注：括号内为标准误差，ns 表示参数估计在0.05水平不显著。Notes: standard error is showed in bracket, ns represents that parameter is not significant at 0.05 level.
下载: 导出CSV
表 4 参数模型和广义加性模型交叉检验统计量
Table 4. Cross-validation for parametric models and generalized additive models
模型
Model平均绝对误差
Mean absolute
bias (MAB)/cm相对百分误差
Relative
percentage
bias (MPB)均方根误差
Root mean
square error
(RMSE)/cm曾伟生等(1997)
Zeng et al. (1997)1.048 8 4.956 9 1.474 1 Bi (2000) 1.018 6 4.814 2 1.450 8 Kozak (2004) 0.957 4 4.525 1 1.361 0 BS 0.911 9 4.309 8 1.310 9 CR 1.040 0 4.915 4 1.468 5 DS 0.923 3 4.364 1 1.331 8 GP 0.915 7 4.327 9 1.318 3 PS 0.923 8 4.366 2 1.323 9 TP 0.916 0 4.329 6 1.317 7
下载: 导出CSV
-
[1] Kozak A. My last words on taper equations[J]. The Forestry Chronicle, 2004, 80(4): 507−515. doi: 10.5558/tfc80507-4. [2] Özçelik R, Crecente-Campo F. Stem taper equations for estimating merchantable volume of Lebanon cedar trees in the Taurus Mountains, southern Turkey[J]. Forest Science, 2016, 62(1): 78−91. doi: 10.5849/forsci.14-212 [3] 姜立春, 刘瑞龙. 基于非线性混合模型的落叶松树干削度模型[J]. 林业科学, 2011, 47(4):101−106. Jiang L C, Liu R L. A stem taper model with nonlinear mixed effects for Dahurian larch[J]. Scientia Silvae Sinicae, 2011, 47(4): 101−106. [4] Schröder T, Costa E A, Valério A F, et al. Taper equations for Pinus elliottii Engelm. in southern Paraná, Brazil[J]. Forest Science, 2015, 61(2): 311−319. doi: 10.5849/forsci.14-054 [5] 姜立春, 李凤日, 刘瑞龙. 兴安落叶松树干削度和材积相容模型[J]. 北京林业大学学报, 2011, 33(5):1−7. Jiang L C, Li F R, Liu R L. Compatible stem taper and volume models for Dahurian larch[J]. Journal of Beijing Forestry University, 2011, 33(5): 1−7. [6] Özçelik R, Karatepe Y, Gürlevik N, et al. Development of ecoregion-based merchantable volume systems for Pinus brutia Ten. and Pinus nigra Arnold. in southern Turkey[J]. Journal of Forestry Research, 2016, 27(1): 101−117. doi: 10.1007/s11676-015-0147-4. [7] Gomat H Y, Deleporte P, Moukini R, et al. What factors influence the stem taper of Eucalyptus: growth, environmental conditions, or genetics?[J]. Annals of Forest Science, 2011, 68(1): 109−120. doi: 10.1007/s13595-011-0012-3. [8] Jiang L, Liu R. Segmented taper equations with crown ratio and stand density for Dahurian larch (Larix gmelinii) in northeastern China[J]. Journal of Forestry Research, 2011, 22(3): 347−352. doi: 10.1007/s11676-011-0178-4. [9] Wiklund K, Konöpka B, Nilsson L O. Stem form and growth in picea abies (L.) karst, in response to water and mineral nutrient availability[J]. Scandinavian Journal of Forest Research, 1995, 10(1−4): 326−332. doi: 10.1080/02827589509382899 [10] Sharma M, Oderwald R G. Dimensionally compatible volume and taper equations[J]. Canadian Journal of Forest Research, 2001, 31(5): 797−803. doi: 10.1139/x01-005. [11] Nigh G, Smith W. Effect of climate on lodgepole pine stem taper in British Columbia, Canada[J]. Forestry, 2012, 85(5): 579−587. doi: 10.1093/forestry/cps063. [12] Lee W K, Biging G S, Son Y, et al. Geostatistical analysis of regional differences in stem taper form of Pinus densiflora in central Korea[J]. Ecological Research, 2006, 21(4): 513−525. doi: 10.1007/s11284-006-0152-3 [13] Newnham R M. A variable-form taper function[M]. Chalk River: Petawawa National Forestry Institute, Forestry Canada, 1988. [14] Rojo A, Perales X, Sánchez-Rodríguez F, et al. Stem taper functions for maritime pine (Pinus pinaster Ait.) in Galicia (Northwestern Spain)[J]. European Journal of Forest Research, 2005, 124(3): 177−186. doi: 10.1007/s10342-005-0066-6. [15] Corral-Rivas J J, Diéguez-Aranda U, Corral R S, et al. A merchantable volume system for major pine species in El Salto, Durango (Mexico)[J]. Forest Ecology and Management, 2007, 238(1−3): 118−129. doi: 10.1016/j.foreco.2006.09.074. [16] Li R, Weiskittel A R. Comparison of model forms for estimating stem taper and volume in the primary conifer species of the North American Acadian Region[J]. Annals of Forest Science, 2010, 67(3): 302. doi: 10.1051/forest/2009109. [17] Gómez-García E, Crecente-Campo F, Diéguez-Aranda U. Selection of mixed-effects parameters in a variable-exponent taper equation for birch trees in northwestern Spain[J]. Annals of Forest Science, 2013, 70(7): 707−715. doi: 10.1007/s13595-013-0313-9. [18] Adamec Z, Drápela K. Generalized additive models as an alternative approach to the modelling of the tree height-diameter relationship[J]. Journal of Forest Science, 2015, 61(6): 235−243. [19] Kuželka K, Marušák R. Input point distribution for regular stem form spline modeling[J/OL]. Forest Systems, 2015, 24 (1): e008 [2019−11−01]. https://revistas.inia.es/index.php/fs/article/view/5815. DOI: 10.5424/fs/2015241-05815. [20] Schmidt M, Kiviste A, Von Gadow K. A spatially explicit height-diameter model for Scots pine in Estonia[J]. European Journal of Forest Research, 2011, 130(2): 303−315. doi: 10.1007/s10342-010-0434-8. [21] Wang Y, Raulier F, Ung C H. Evaluation of spatial predictions of site index obtained by parametric and nonparametric methods, a case study of lodgepole pine productivity[J]. Forest Ecology and Management, 2005, 214(1−3): 201−211. doi: 10.1016/j.foreco.2005.04.025. [22] Guisan A, Edwards T C, Hastie T. Generalized linear and generalized additive models in studies of species distributions: setting the scene[J]. Ecological Modelling, 2002, 157(2−3): 89−100. doi: 10.1016/S0304-3800(02)00204-1. [23] M’Hirit O, Postaire J G. A nonparametric technique for taper function estimation[J]. Canadian Journal of Forest Research, 1985, 15(5): 862−871. doi: 10.1139/x85-139. [24] 余黎, 雷相东, 王雅志, 等. 基于广义可加模型的气候对单木胸径生长的影响研究[J]. 北京林业大学学报, 2014, 36(5):22−32. Yu L, Lei X D, Wang Z Y, et al. Impact of climate on individual tree radial growth based on generalized additive model[J]. Journal of Beijing Forestry University, 2014, 36(5): 22−32. [25] 曾伟生, 廖志云. 削度方程的研究[J]. 林业科学, 1997, 33(2):32−37. Zeng W S, Liao Z Y. Research on taper equation[J]. Scientia Silvae Sinicae, 1997, 33(2): 32−37. [26] Bi H. Trigonometric variable-form taper equations for Australian eucalypts[J]. Forest Science, 2000, 46(3): 397−409. [27] Wood S N. Generalized additive models: an introduction with R[M]. 2 ed. Boca Raton: CRC Press, 2017. [28] Wood S N. P-splines with derivative based penalties and tensor product smoothing of unevenly distributed data[J]. Statistics and Computing, 2017, 27(4): 985−989. doi: 10.1007/s11222-016-9666-x. [29] Rodríguez F, Lizarralde I, Bravo F. Comparison of stem taper equations for eight major tree species in the Spanish Plateau [J/OL]. Forest Systems, 2015, 24(3): e034 [2019−12−02]. http://dx.doi.org/10.5424/fs/2015243-06229. [30] Lane S E, Robinson A P, Baker T G. The functional regression tree method for diameter distribution modelling[J]. Canadian Journal of Forest Research, 2010, 40(9): 1870−1877. doi: 10.1139/X10-119. -
点击查看大图
图(4) / 表 (4)
计量
- 文章访问数: 177
- HTML全文浏览量: 44
- PDF下载量: 29
- 被引次数: 0
