2016, Vol. 36
文章信息
- 王少杰, 邓华锋, 黄国胜, 王雪军
- WANG Shaojie, DENG Huafeng, HUANG Guosheng, WANG Xuejun
- 基于哑变量的油松人工林和天然林生长模型
- Dummy variables models in Pinus tabulaeformis artificial forest and natural forest growth
- 森林与环境学报, 2016, 36(3): 325-331
- Journal of Forest and Environment, 2016, 36(3): 325-331.
- http://dx.doi.org/10.13324/j.cnki.jfcf.2016.03.012
-
文章历史
- 收稿日期: 2015-11-27
- 修订日期: 2015-12-26
2. 国家林业局规划调查设计院, 北京 100714
2. Academy of Forest Inventory and Planning, State Forestry Administration, Beijing 100714, China
油松(Pinus tabulaformis Carr.)，阳性树种，喜生长在排水良好的酸性土壤中，是北京市主要的造林绿化树种之一。依据林分起源，油松有人工林和天然林2种，人工林相对较多，油松在保护生态环境安全方面发挥着至关重要的作用[1]。一般来说，针叶人工林的生产力比天然林的高，但天然林具有较强抵抗病虫害和维持土壤肥力的能力，在满足物种资源的保存及环境保护的要求等方面，天然林都优于工业人工林[2]。为了优化林分的质量，提高收益，管理好林分资源，建立林分生长模型来全面掌握林分动态变化规律是一个非常有效的途径。国内外学者已经构建了大量的林分断面积和蓄积量生长模型[3, 4]，李春明等[5]借助基础经验模型方程，通过对模型再次参数化并选择合适的自变量拟合长白落叶松(Larix olgensis Henry.)和杉木[Cunninghamia lanceolata (Lamb.)Hook.]林分断面积，取得了较好的预测精度。李忠国等[6]通过引入区域特征哑变量，建立了北亚热带高山区和暖温带中山区日本落叶松的材积生长模型。高东启等[7]在研究油松断面积生长模型时，在经验模型方程基础上，引入2个分别代表间伐与否的哑变量和代表间伐指标变量建立模型，取得了较好的预估精度。
在生长模型的研究中，所建立的模型一般是基于同一林分起源单独建立的，因而无法用同一个模型来表达不同林分起源林分的生长状况。理论上，相同的树种不同林分起源的林分是需要分开单独建模的，但人工林和天然林存在一定相关关系，分开建模增加了工作量，且对于不同林分起源合并建立断面积和蓄积量生长模型很少有报道。目前研究主要是用在生物量建模上，如符利勇等[8]通过把不同起源作为哑变量，针对不同起源的马尾松(Pinus massoniana Lamb.)构建生物量相容性生长模型，随后引入2个哑变量分别代表不同地域和不同林分起源，构建了东北落叶松地上生物量通用方程[9]。利用北京市油松人工林和天然林的定期清查数据，在相关的树木生长方程中引入代表林分起源的哑变量，将人工林和天然林合并起来建立油松林分断面积和蓄积量生长模型，并从断面积和蓄积量模型中分别选择出拟合较好的生长模型，最后用检验样本对其预估结果进行适应性分析，从而建立精度高、且能适用于人工林和天然林的油松林分生长模型。
1 研究区概况北京市地处华北平原北端，北纬39°28′-41°05′，东经115°25′-117°30′，与天津市、河北省毗邻，属于暖温带半湿润大陆性季风型气候，气候干燥，降水分布不均匀。北京市下辖有16个区，在怀柔、密云、平谷、海淀、延庆和平谷区都分布有油松。
2 材料与方法 2.1 数据来源试验采用北京市2001、2006、2011年组成的3期一类连续清查的油松固定样地数据，样地面积0.066 7 hm2。剔除数据缺失、记录不详和异常的样地，从3期固定样地中选出133块样地，其中，人工林样地103块，天然林样地30块。人工林样地林龄最大57 a，最小林龄7 a；天然林样地林龄最大63 a，最小林龄23 a。样地基本情况见表 1，建模和检验样地的变量因子基本情况见表 2。对133块样地，随机抽样选出88个样本用于建模，其中人工林73块、天然林15块；剩下的45块样地用于检验数据，其中人工林30块、天然林15块。
| 样地 Plot | 人工林Artificial forest | 天然林 Natural forest | |||
| 样地数量 Number of plots | 所占比例 Percentage/% | 样地数量 Number of plots | 所占比例 Percentage/% | ||
| 幼龄林 Sapling forest | 23 | 23 | 3 | 10 | |
| 中龄林 Middle aged forest | 25 | 24 | 21 | 70 | |
| 近龄林 Nearmature forest | 25 | 24 | 5 | 17 | |
| 成龄林 Mature forest | 30 | 29 | 1 | 3 | |
| 合计 Total | 103 | 100 | 30 | 100 | |
| 数据类别Type ofdata | 样本数Number ofsamples | 平均林龄Averageage/a | 平均胸径AverageDBH/cm | 平均树高Averageheight/m | 株数Number oftrees | 林分蓄积量Stand volume/(m3·hm-2) | 林分断面积Stand basal area/(m2·hm-2) |
| 建模数据 Fit data | 88 | 33 | 12.3 | 6.9 | 783 | 38.800 | 9.772 |
| 检验数据 Validation | 45 | 36 | 12.3 | 7.0 | 621 | 42.143 | 10.276 |
哑变量又称为虚拟变量，是自变量经过量化后所得的，一般取值0或者1。引入哑变量[10]，可以将人工林和天然林2种类型的林分用定性代码0或1表示。将其中一种起源的林分样地定义为Zi，将定性数据Zi转化为(0，1)形式。
| $ {Z_{\rm{i}}}{\rm{ = }}\left\{ {\begin{array}{*{20}{c}} {1,当{Z_i}为第i个类型林分时,}\\ {0,否则。} \end{array}} \right. $ |
式中：i=1，2；Z1,Z2分别表示人工林分和天然林分的定性代码。
2.3 地位级指数的计算由于采用的一类连续清查数据中没有直接提供林分优势木的平均树高数据，故以林分平均树高及林龄为自变量建立方程模型，把经过拟合得到的地位级指数作为立地质量指标[11]。
利用Richards模型拟合林分平均树高生长过程，其具体形式如下
| $ H = a{[1 - exp\left( {bT} \right)]^c} $ | (1) |
式中：H为林分平均树高(m)；T为平均林龄(a)；a 、b 、c为模型的待估测参数。
用林分平均树高及平均林龄计算地位级指数
| $ I = H\frac{{1 - \exp {{\left( { - b{T_0}} \right)}^c}}}{{1 - exp{{\left( { - bT} \right)}^c}}} $ | (2) |
式中：I为地位级指数；H为平均树高(m)；b 、c为方程的参数；T为平均林龄(a)；T0为油松林分的标准林龄，取25 a[12, 13]。
2.4 林分密度指数的计算林分密度指数(stand density index，SDI)是一个适用性较广、评价较高的密度指标。
| $ SDI = M{\left( {{D_g}/{D_0}} \right)^\alpha } $ | (3) |
式中： SDI为林分密度指数；N为单位面积株数(株·hm-2)；Dg为林分的平均胸径(cm)；D0为标准平均直径(cm)，取10 cm；α为自然稀疏率，取1.601 2[13, 14]。
2.5 林分断面积生长模型林分断面积是预估林分收获量的重要变量，具有较高的稳定性和可预测性，因此国内外很多学者在这方面作了较多研究，Richards和Schumacher模型是目前主要的2种树木生长方程形式[15, 16, 17]。这2类模型参数中都有与断面积最大值相联系的参数，模型中参数的意义较明确，林分断面积的生长快慢也可以在模型参数中体现，拟合时比较容易收敛，故适合用来拟合人工林和天然林分断面积的生长。以Richards模型为基础模型，拟合油松断面积，在模型方程上继续参数化，林分立地质量通过地位级指数间接反映，密度指标用林分密度指数或者林分每公顷株数来表示，同时在模型中引入一个代表林分起源的哑变量，最终模型的形式如下。
| $ G = \left( {{a_0}{Z_1} + {a_1}{Z_2}} \right){I^{a2}}{\left\{ {1 - \exp [- {a_3}{{\left( {SDI/1000} \right)}^{a4}}\left( {T - {T_{1.3}}} \right)]} \right\}^{a5}} $ | (4) |
| $ G = \left( {{a_0}{Z_1} + {a_1}{Z_2}} \right){I^{a2}}{\left\{ {1 - \exp [- {a_3}{{\left( {N/1000} \right)}^{a4}}\left( {T - {T_{1.3}}} \right)]} \right\}^{a5}} $ | (5) |
式中：G为林龄T时的林分断面积(m2·hm-2)；Z1 、Z2分别为天然林和人工林的定性代码；I为地位级指数；N为单位面积株数(株·hm-2)；T为林分平均林龄(a)；T1.3为林木生长到1.3 m时的林龄；a0-a5为待估参数。
2.6 林分蓄积量生长模型通常情况下，林分立地质量(地位级指数)的大小与模型中的渐近值参数有着密切的相关关系，且林分密度主要控制拟合模型的曲线形状[18]。同样拟合林分蓄积量模型以Richards模型为基础模型，选用林分断面积和林分密度指数作为密度指标，并引入林分起源作为哑变量，模型最终的形式如下
| $ V = \left( {{a_0}{Z_1} + {a_1}{Z_2}} \right){I^{a2}}{\left\{ {1 - \exp \left( { - {a_3}{G^{a4}}T} \right)} \right\}^{a5}} $ | (6) |
| $ V = \left( {{a_0}{Z_1} + {a_1}{Z_2}} \right){I^{a2}}{\left\{ {1 - \exp [- {a_3}{{\left( {SDI/1000} \right)}^{a4}}T]} \right\}^{a5}} $ | (7) |
式中：V为林分蓄积量(m3·hm-2)。
2.7 参数估计与模型检验应用Forstat 2.2软件对油松林分断面积和蓄积量模型采用最小二乘回归参数估计方法进行拟合，并用Excel 2010进行简单的数据处理，采用平均偏差(mean deviation，MD)、平均绝对偏差(mean absolute deviation，MAD)、决定系数(R2)及预估精度(P)4个评价指标对模型的结果进行适应性检验分析。
| $ {\rm{MD = }}\sum {\left( {{y_i} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over y} }_i}} \right)/n} $ | (8) |
| $ {\rm{MAD = }}\sum {\left| {{y_i} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over y} }_i}} \right|/n} $ | (9) |
| $ {R^2} = 1 - \sum {{{\left( {{y_i} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over y} }_i}} \right)}^2}/\sum {{{\left( {{y_i} - {{\bar y}_i}} \right)}^2}} } $ | (10) |
| $ P=\left[ 1-{{t}_{0.05}}\sqrt{\sum{{{\left( {{y}_{i}}-\overset\frown{{{y}_{i}}} \right)}^{2}}/n\left( n-p \right)}}/\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\bar{y}} \right]\times 100%$ | (11) |
式中：yi为模型断面积实际值(m2·hm-2)或蓄积量的实际值(m3·hm-2)； $ {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over y} }_i} $为模型断面积预测值(m2·hm-2)和蓄积量的预估值(m3·hm-2)；$ {\bar y} $为实测值的平均值(m2·hm-2或m3·hm-2)； $ \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\bar{y}} $为模型预估值的平均值(m2·hm-2或m3·hm-2)；n为样本数；t0.05为置信水平0.05时的t分布值；p为模型中参数的个数。
3 结果与分析选用Richards方程拟合林分平均高生长过程结果H=10.321[1-exp(-0.037T)]1.125，决定系数R2为0.677 0，待估参数a，b，c分别为10.321，0.037，1.125。将结果代入公式(2)中，求出每个样地的地位级指数。当油松的树高为1.3 m时，由拟合的公式得出，T1.3=4.0，基本符合油松的生长习性，将T1.3代入公式(4)和公式(5)，对油松林分断面积生长模型继续进一步拟合，林分断面积和蓄积量生长模型参数如表 3所示。拟合油松断面积模型时，当密度指标采用林分密度指数时，油松断面积模型R2为0.938 0；而密度指标采用每公顷株数时，油松断面积模型的R2只有0.878 3，说明所拟合的林分断面积模型效果都较好，且以林分密度指数为密度指标更能反映林分断面积生长变化规律，即模型(4)拟合林分断面积更好。拟合油松林分蓄积量模型，密度指标采用林分断面积时，拟合的模型R2为0.991 8；而密度指标采用林分密度指数时，拟合的油松蓄积量模型R2为0.919 3，选用2种密度指标的模型拟合效果都较好，且以林分断面积为密度指标的模型(6)比模型(7)拟合的林分蓄积量效果更好。因此，选用模型(4)、模型(6)继续进一步检验分析。
| 参数Parameter | 林分断面积生长模型Stand basal area growth model | 林分蓄积生长模型Stand volume growth model | |||
| 模型(4)Model (4) | 模型(5)Model (5) | 模型(6)Model (6) | 模型(7)Model (7) | ||
| a0 | 1.600 2 | 0.411 4 | 15.440 9 | 4.528 4 | |
| a1 | 1.796 2 | 0.484 5 | 15.914 7 | 5.302 0 | |
| a2 | 0.206 6 | 0.793 3 | 0.094 4 | 0.351 0 | |
| a3 | 0.017 1 | 0.035 9 | 0.000 0 | 0.025 3 | |
| a4 | 2.360 4 | 0.670 9 | 4.577 6 | 1.314 0 | |
| a5 | 0.389 2 | 1.726 0 | 0.218 7 | 0.777 4 | |
| R2 | 0.938 0 | 0.878 3 | 0.991 8 | 0.919 3 | |
对拟合的模型进行适应性检验分析，综合统计建模数据和检验数据的各项评价指标，结果如表 4所示。在建模数据中，林分断面积和蓄积量的R2均在0.9以上，P均超过96%，MD 、MAD等误差都较小，说明所建立的模型对林分断面积和蓄积量拟合精度高。在检验数据中，4个评价指标的数值都较优，进一步证明所拟合的油松林分断面积、蓄积量生长模型预测效果较好。
| 数据类别 Types of data | 模型 Model | MD | MAD | R2 | P/% |
| 建模数据 Fit data | 林分断面积生长模型Stand basal area growth model(4) | -0.019 5 | 0.973 0 | 0.938 0 | 96.46 |
| 林分蓄积量生长模型Stand volume growth model(6) | -0.042 0 | 1.755 6 | 0.991 8 | 98.63 | |
| 检验数据Validation data | 林分断面积生长模型Stand basal area growth model(4) | -0.239 9 | 1.484 3 | 0.922 0 | 92.99 |
| 林分蓄积量生长模型Stand volume growth model(6) | -0.779 6 | 2.143 9 | 0.993 5 | 97.53 |
表 5列出了模型对人工林和天然林的预测效果，对于人工林，检验指标中P均在96%以上，MD和MAD均较小，绝对值均小于2，说明模型拟合人工林分断面积和蓄积量精度高。对比断面积和蓄积量检验指标知，蓄积量的拟合精度高于断面积；对于天然林，蓄积量的拟合效果较好，R2和P均较高，断面积的拟合效果稍差，但到达了理想精度。综合分析人工林和天然林的拟合效果，所建模型能够较好地预测人工林和天然林的生长过程，且模型对人工林的拟合准确性高于天然林。
| 林分起源 Stand origin | 模型 Model | MD | MAD | R2 | P/% |
| 人工林Artificial forest | 林分断面积生长模型 Stand basal area growth model(4) | 0.027 0 | 0.836 6 | 0.929 0 | 96.93 |
| 林分蓄积量生长模型 Stand volume growth model(6) | -0.004 5 | 1.550 1 | 0.994 1 | 98.82 | |
| 天然林 Natural forest | 林分断面积生长模型 Stand basal area growth model(4) | -0.506 0 | 2.205 7 | 0.843 8 | 91.88 |
| 林分蓄积量生长模型 Stand volume growth model(6) | -1.284 3 | 3.043 6 | 0.986 3 | 97.11 |
分别建立检验数据、人工林与天然林林分断面积和蓄积量实测值、预测值线性相关方程。检验数据的林分断面积和蓄积量R2分别为0.922 0、0.993 5(图 1)；图 2中的散点基本分布在回归线上，说明模型对于人工林断面积和蓄积量的拟合效果较好；图 3中，天然林的断面积和蓄积量R2分别为0.843 8、0.986 3，相对于人工林来说偏低。综合分析3组回归图可知，林分蓄积量的R2均高于林分断面积的R2，实际值的分布更加靠近回归线，预估效果更好。3组回归图中，方程的决定系数R2都接近1，常数项接近0，表明所拟合的生长模型的预测效果均较理想，模型具有很好的参考价值。
 |
图 1 检验数据林分断面积、蓄积量实测值与预测值 Fig.1 Correlations of observed and predicted values on validation data for stand basal area and stand volume |
 |
图 2 人工林分断面积、蓄积量实测值与预测值 Fig.2 Observed and predicted values of stand basal area and stand volume on artificial forest |
 |
图 3 天然林分断面积、蓄积量实测值与预测值 Fig.3 Observed and predicted values of stand basal area and stand volume on natural forest |
利用一类连续清查数据，基于Richards模型方程引入一个代表林分起源的哑变量，将人工林和天然林合并，使得2种不同起源的林分用一个统一的林分断面积、蓄积量生长模型来表示，这样有利于减少工作量，提高工作效率，同时也可以解决模型不相容的问题，使油松人工林和天然林合并一起建立生长模型。 采用MD、MAD 、R2 、P四个数学指标统计量及通过分析比较这些指标的大小，作为建模精度和模型适应性检验的标准。以往所建立林分生长模型能很好地说明单一起源的林分生长情况，而引入哑变量建立一个包含不同林分起源的生长模型，使所建立的模型同时适用于人工林和天然林，即减少了工作量，又能很好地解决不同起源林分合并建模不相容问题。
Richards模型是应用最为广泛，适用性较强的一种生长曲线方程[19, 20]。因此研究选用Richards方程为基础模型拟合林分断面积和林分蓄积量[21]。研究结果表明，模型拟合精度都较高，且拟合模型经检验对人工林和天然林断面积和蓄积量预测都很好。另外，无论是林分断面积还是林分蓄积量，人工林的拟合精度均高于天然林。林分断面积生长模型中，模型的密度指标采用林分密度指数比采用每公顷株数的预估精度更高，能更准确地反映油松林分断面积生长变化规律，这与高东启等[22]拟合蒙古栎断面积生长模型时的结果一致。在拟合林分蓄积量时，林分断面积作为密度指标，模型的拟合效果更佳。由于拟合的林分平均高生长方程的决定系数偏低，可能会降低油松生长模型预测结果精度。为了方便模型能够更好地在今后的林业发展中应用，须在增加样地数据的同时考虑时间效应对模型的影响，建立混合模型进行深入地研究，为下一步科学经营油松人工林和天然林提供重要的理论依据。
| [1] | 高东启,邓华锋,王海宾,等. 基于哑变量的蒙古栎林分生长模型[J]. 东北林业大学学报,2014,42(1):61-64. |
| [2] | 徐化成. 人工林和天然林的比较评价[J]. 世界林业研究,1991(3):50-56. |
| [3] | BREIDENBACH J,KUBLIN E,MCGAUGHEY R,et al. Mixed-effects models for estimating stand volume by means of small footprint airborne laser scanner data[J]. The Photogrammetric Journal of Finland,2008,21(1):4-15. |
| [4] | ANTÓN-FERNÁNDEZ C,BURKHART H E,AMATEIS R L. Modeling the effects of initial spacing on stand basal area development of loblolly pine[J]. Forest Science,2012,58(2):95-105. |
| [5] | 李春明,杜纪山,张会儒. 间伐林分的断面积生长模型研究[J]. 林业资源管理,2004(3):52-55. |
| [6] | 李忠国,孙晓梅,陈东升,等. 基于哑变量的日本落叶松生长模型研究[J]. 西北农林科技大学学报(自然科学版),2011,39(8):69-74. |
| [7] | 高东启,邓华锋,蒋益,等. 油松林分断面积生长预估模型研究[J]. 西南林业大学学报,2015,35(1):42-46. |
| [8] | 符利勇,雷渊才,孙伟,等. 不同林分起源的相容性生物量模型构建[J]. 生态学报,2014,34(6):1461-1470. |
| [9] | 符利勇,唐守正,张会儒,等. 东北地区两个主要树种地上生物量通用方程构建[J]. 生态学报,2015,35(1):150-157. |
| [10] | 唐守正,郎奎建,李海奎. 统计与生物数学模型计算:ForStat教程[M]. 北京:科学出版社,2009. |
| [11] | 江传阳. 福建柏地位级指数曲线模型的研制[J]. 林业勘察设计,2014(2):5-8,13. |
| [12] | 李佩萍,武建林. 晋中东部山区油松人工林地位指数表的编制[J]. 山西林业科技,1999(4):19-23,27. |
| [13] | 孟宪宇. 测树学[M]. 北京:中国林业出版社,2006. |
| [14] | 刘永霞. 北京山地油松林分生长过程数量化模拟研究[D]. 北京:北京林业大学,2007. |
| [15] | 杜纪山. 用二类调查样地建立落叶松单木直径生长模型[J]. 林业科学研究,1999,12(2):160-164. |
| [16] | 李春明. 利用非线性混合模型进行杉木林分断面积生长模拟研究[J]. 北京林业大学学报,2009,31(1):44-49. |
| [17] | 王霓虹,杨英奎,戴巍,等. 基于Richards方程的落叶松人工林断面积生长模型[J]. 森林工程,2015,31(1):22-25,53. |
| [18] | 郎荣,许建初,TENNIGKEIT T,等. 基于样方数据的云南松林分生长模型研究:以云南省保山市杨柳白族彝族乡为例[J]. 植物分类与资源学报,2011,33(3):357-363. |
| [19] | 杜纪山,唐守正. 杉木林分断面积生长预估模型及其应用[J]. 北京林业大学学报,1998,20(4):1-5. |
| [20] | 池上评,陈金章,江传阳,等. 基于间隔期的福建柏人工林动态生长模型及应用[J]. 福建林学院学报,2014,34(4):304-308. |
| [21] | 胡晓龙. 长白落叶松林分断面积生长模型的研究[J]. 林业科学研究,2003,16(4):449-452. |
| [22] | 高东启,邓华锋,程志楚,等. 蒙古栎间伐林分和未间伐林分生长模型研究[J]. 中南林业科技大学学报,2014,34(2):50-54. |