云南松(Pinus yunnanensis Franch)又称飞松，分布的地理位置约在北纬23°~ 32°，东经98°30′~ 106°之间，构成一个不规则的多角形分布区，其中云南省分布最为密集^[1]。云南松不仅是我国西南地区的乡土树种，也是该地区更新造林的先锋树种，自然条件下常形成大面积纯林，在经济、社会和生态效益方面具有一定价值，在林业生产中也占有重要地位。因此，为实现对云南松林分的科学经营管理，准确预估林分的生长、收获量显得尤为重要，必须对林分的生长变化规律有充分的了解。

近年来，许多学者对云南松进行了大量的研究，主要集中在遗传多样性、群落结构及演替、生长生理响应特征、生物量等方面^[2-3]，如蔡年辉等^[4]、王健敏等^[5]分别就云南松天然林与人工林的群落结构和物种多样性进行了对比研究；孙宝刚等^[6]、杨阳等^[7]、冉啟香等^[8]对云南松生物量模型进行了相关研究。但是，对云南松蓄积生长模型的研究较少，尤其是间伐措施对拟合云南松生长模型的影响少有报道。众所周知，现实林分的生长经常会受到抚育间伐的影响，而在林分研究中忽略间伐措施对林分生长的效应可能会在一定程度上影响研究结果的准确性，因此，有必要通过对间伐及未间伐林分的研究来深入了解间伐对林分生长的影响。近年来，国内外已有诸多关于抚育间伐的报道^[9-10]，如胥富强等^[11]分析了抚育间伐对华山松林木胸径、树高和材积的影响；鲍斌等^[12]以马尾松为研究对象，探讨了间伐在林分生长方面的影响。但目前间伐林分的研究中较少以云南松作为研究对象。从理论上来说，为保证模型的适用性及较高的预估精度，不同类型的林分需要分开单独建模，但这会导致工作量及成本增加等问题，而哑变量模型为合并建模提供了可能途径。为此，笔者在前人研究的基础上，利用云南省云南松一类清查数据，以Richards和Schumacher模型为基础，选取常用的两个林分密度指标，即林分密度指数和林分断面积，并利用地位级指数反映立地质量，建立云南松蓄积生长模型，在最优的传统生长模型上，通过引入哑变量建立同时适用于间伐林分、未间伐林分的林分蓄积生长模型。

1 材料与方法 1.1 数据来源

本研究所采用的数据为云南省第7、8次森林资源一类清查数据，来源于国家林业局调查规划设计院，每个样地面积为0.066 7 hm²。主要对每个样地的林分平均胸径、林分平均树高、林分蓄积、林分平均年龄、郁闭度、坡向、坡位及海拔高度等因子进行了记录。从2期数据中剔除数据缺失、记录不详及明显有误的样地，最终选出135块样地，其中，间伐样地94块，未间伐样地41块。对135个样地进行随机抽样，选出85个样本作为建模数据，其中间伐样地63块，未间伐样地22块，剩余50个样本用于检验，包括31块间伐样地和19块未间伐样地。样地分布和样地基本情况分别见表 1、表 2。

1.2 林分密度指数的确定

林分密度指数(stand density index，SDI)是指现实林分的株数换算到比较直径时所具有的单位面积林木株数^[13]，计算方法如下：

${I_{{\rm{SD}}}} = N{(D/{D_0})^b}$

(1)

式中：I_SD为林分密度指数；N为现实林分每公顷株数；D为现实林分平均直径；D₀为标准平均直径；b为自然稀疏率。参考郎荣等^[13-14]的研究，云南松D₀取15 cm，b为1.373 27。

1.3 地位级指数的计算

由于一类调查数据中的树高是林分平均高而不是优势木平均高，故通过拟合林分平均高和林分平均年龄的关系，将建模所得的地位级指数作为评价立地质量的指标^[15]。通过对几个常用模型的拟合对比，最终选用Mistcherlich式拟合林分平均高生长曲线，其基本形式如下：

$H = a[1 - {\rm{exp}}\left( { - bT} \right)]$

(2)

式中：H为林分平均高(m)；T为林分平均年龄(a)；a、b为参数。

地位级指数计算公式如下：

${I_{{\rm{SC}}}} = H[1 - {\rm{exp}}( - b{T_0})\left] / \right[1 - {\rm{exp}}\left( { - bT} \right)]$

(3)

式中：I_SC为地位级指数；T₀为基准年龄，云南松取20 a^[14]。

1.4 哑变量的定义

哑变量(dummy variable)又称虚拟变量，是量化了的自变量，用以反映质的属性，通常取值为0或1。哑变量的定义为：对于定性数据x，用变量δ(x, i)表示为：

$\delta \left( {x,{\rm{ }}i} \right) = \left\{ {\begin{array}{*{20}{l}} {1,}&{当x为第i等级时,}\\ {0,}&{否则。} \end{array}} \right.$

(4)

这种方法叫做定性因子(0，1) 化展开，变量δ(x, i)称为哑变量^[16]。通过引入间伐措施作为哑变量，将间伐林分和未间伐林分用定性代码0或1表示，第i种林分类型表示为S_i，定性数据S_i转化为(0，1) 形式:

${s_i} = \left\{ {\begin{array}{*{20}{l}} {1,}&{当x为第i种林分类型时,}\\ {0,}&{否则。} \end{array}} \right.$

(5)

式中：i=1，2；S₁，S₂分别表示未间伐林分和间伐林分的定性代码。

1.5 生长模型的选择

目前对生长模型的研究很多，模型种类、形式多样，但应用较多的主要有Logistic模型、单分子式、Richards方程和Schumacher模型等。其中，Richards方程由于其广泛的适应能力及较好的解释性，得到了大范围的应用，而Schumacher模型因形式简单、计算简便等优点，在实际应用中得到了众多研究者的偏爱。因此，基于这两种模型，通过引入地位级指数、选取林分密度指数和林分断面积作为林分密度指标来拟合林分蓄积生长模型，模型形式如下：

$V = {a_1}I_{SC}^{{a_2}}{({I_{{\rm{SD}}}}/1\,000)^{{a_3}}}{\rm{exp}}( - {a_4}/T)$

(6)

$V = {a_1}I_{{\rm{SC}}}^{{a_2}}{G^{{a_3}}}{\rm{exp}}( - {a_4}/T)$

(7)

$V = {a_1}I_{{\rm{SC}}}^{{a_2}}{\{ 1 - {\rm{exp}}[ - {a_3}{({I_{{\rm{SD}}}}/1\,000)^{{a_4}}}T]\} ^{{a_5}}}$

(8)

$V = {a_1}I_{{\rm{SC}}}^{{a_2}}{[1 - {\rm{exp}}( - {a_3}{G^{{a_4}}}T)]^{{a_5}}}$

(9)

1.6 参数估计及模型评价

通过ForStat、Excel实现数据处理及模型参数估计。采用平均绝对偏差(mean absolute deviation, MAD)、均方根误差(root-mean-square error, RMSE)、模型准确性(v)^[17]、决定系数(R²)和预估精度(P)等指标对模型进行评价和对比。几项指标中，R²和P越接近1越好，其余3项值越小说明模型模拟精度越高。各指标计算表达式如下：

${D_{{\rm{MA}}}} = \sum \left| {\frac{{{y_i} - {{\hat y}_i}}}{n}} \right|$

(10)

${E_{{\rm{RMS}}}} = \sqrt {\frac{{\sum {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}} }}{n}} $

(11)

$v = \frac{{\sum {{{({y_i} - {{\hat y}_i})}^2}} - \sum {{{[({y_i} - {{\hat y}_i})]}^2}/n} }}{{n\left( {n - 1} \right)}}$

(12)

${R^2} = 1 - \frac{{\sum {{{({y_i} - {{\hat y}_i})}^2}} }}{{\sum {{{({y_i} - y)}^2}} }}$

(13)

$P=[1-{{t}_{0.05}}\frac{\sum{{{({{y}_{i}}-{{{\hat{y}}}_{i}})}^{2}}}}{n\left( n-p \right)}/\hat{\bar{y}}]\times 100%$

(14)

式中：D_MA为平均绝对偏差；E_RMS均方根误差；y_i为林分蓄积量实测值(m³·hm^-2)；${{{\hat y}_i}}$为林分蓄积量预估值(m³·hm^-2)；y为实测值的平均值(m³·hm^-2)；${\hat{\bar{y}}}$为模型预估值的平均值(m³·hm^-2)；n为样本数；t_0.05为置信水平为0.05时的t分布值；p为模型中参数的个数。

2 结果与分析 2.1 传统模型的拟合结果

采用Mistcherlich式拟合林分平均高生长结果为H=18.913 319[1-exp(-0.020 188T)]，决定系数R²为0.858 4。将结果代入公式(3)，计算每个样地的地位级指数。表 3列出了模型6~9的参数估计值及各项评价指标。所有模型的决定系数R²都大于0.95，且模型9的决定系数最高，达到0.971 6，模型6的决定系数最小；模型决定系数从大到小排序为:模型9＞模型7＞模型8＞模型6，其中模型9和模型7选用的林分密度指标为林分断面积，而模型8和模型6选用的林分密度指标为林分密度指数，这表明相对于林分密度指数来说，林分断面积更能反映林分蓄积量的生长规律。当选用相同密度指标时，Richards模型的决定系数均大于Schumacher模型，说明Richards模型可能更适合模拟云南省云南松天然林的蓄积生长。从模型预估精度及各项评价指标来看，4个模型的拟合效果都较好，预估精度P值均超过95%，最高达到96.79%。比较模型评价指标可知，模型9无论R²还是P值都比其他3个模型高，且其余3项误差指标均较其他3个模型更好。因此，选择模型9作为最优传统模型来构造哑变量模型。

2.2 哑变量模型及其拟合结果

根据传统模型拟合结果，在最优传统模型中加入哑变量。经过尝试发现，在不同参数中加入哑变量时，模型的决定系数有差异，其中，在参数a₄中加入哑变量后所得模型的决定系数最高，因此，选择该模型作为哑变量模型，其形式如下：

$V = {a_1}I_{{\rm{SC}}}^{{a_2}}{[1 - {\rm{exp}}( - {a_3}{G^{{a_4}{S_1} + {a_5}{S_2}}}T)]^{{a_6}}}$

(15)

式中：S₁、S₂分别代表未间伐林分和间伐林分。利用建模数据拟合模型15，参数a₁~a₆分别为79.077 9、0.660 4、0.000 7、1.070 1、1.014 5、0.939 0，R²为0.973 0，P为96.84%。

2.3 模型拟合结果

综合利用建模数据和检验数据对最优传统模型(模型9) 和哑变量模型(模型15) 进行模型精度的验证，统计各项模型评价指标，结果如表 4所示。无论建模还是检验数据，相比传统模型，引入哑变量后的模型的决定系数和预估精度值都有所提高，平均绝对偏差、均方根误差和模型准确性均有所下降，模型的拟合优度得到提高。

区分间伐林分和未间伐林分来进一步检验哑变量模型对两种类型林分的预测效果(表 5)。结果表明：所建哑变量模型对两种类型的林分蓄积量预测的决定系数均大于0.980 0；对于间伐林分，平均绝对偏差、均方根误差和模型准确性分别为1.506 0、1.757 9、0.099 8，预测精度超过98%；对于未间伐林分，平均绝对偏差、均方根误差和模型准确性分别为3.094 2、4.342 2、0.982 8，预测精度相对间伐林分略低，但也超过92%，达到了理想精度，说明所建哑变量蓄积生长模型对间伐及未间伐林分的预测效果均较好，且模型对间伐林分的拟合准确性更高。

图 1为检验数据两种类型林分的蓄积量实测值与预测值散点图。对于检验数据总体、间伐林分和未间伐林分其蓄积实测值和估计值的线性方程的R²分别为0.995 8、0.998 8、0.987 3。检验数据总体和间伐林分的常数项及回归系数均分别比较接近于0和1；未间伐林分其回归系数接近于1，而常数项偏离0较远，但散点基本分布在回归线上，拟合效果也较好。综合分析三者，进一步说明所建的含哑变量的蓄积生长模型的预估效果较好，模型具有一定的参考价值，可以用于云南省云南松间伐及未间伐林分的蓄积预估。

3 结论与讨论

利用云南省一类清查数据，以Schumacher和Richards模型为原型，考虑林分立地质量差异，选取林分密度指数、林分断面积作为林分密度指标，建立了云南省云南松天然林蓄积生长模型。结果表明：4个传统模型的拟合效果都较好，模型的决定系数均在0.950 0以上，不同密度指标对模型的决定系数有直接影响，且林分断面积更适合作为密度指标来反映林分蓄积生长规律，这与高东启等^[18]、冉啟香等^[19]的研究结果一致。在2个选用林分断面积作为林分密度指标的模型中，以Richards模型为原型的模型9的拟合效果更佳，因此，选择该模型作为最优传统模型。

在模型9中引入代表间伐的哑变量，使所建模型同时适用于两种类型的林分，这样一方面能减少外业调查工作量，另一方面不同类型林分合并建模不相容的问题也能得到有效解决。研究表明：加入哑变量能提高模型的决定系数及预估精度，各项误差也相应减少，哑变量模型的拟合效果更好。区分林分类型对哑变量模型进行适应性检验，结果表明：哑变量模型对两种类型林分预测效果都较好，且对间伐林分的预估效果更好。

由于所用数据中未间伐样地数量有限，可能对模型的预估精度产生影响，且在研究中未考虑时间效应的影响，为方便模型在林业生产实践中更好的应用，须获取更多数据做进一步研究。在今后云南松林分的研究中，可以考虑时间效应、林分起源等因素对林分生长的影响，采用混合模型的方法，以期为今后云南松林分的经营管理提供参考。