Develop Annual Stand Volume Growth Model of Chinese fir Including Different Stand Density Indices

样地设置在福建武夷山北部的邵武市，位于117°43′ E， 27°05′ N。地貌特征主要是高山和低丘陵地区，海拔为250～700 m，坡度为25°～35°，为亚热带季风气候，年平均温度17.7 ℃，1月平均温度6.8 ℃，7月平均温度28 ℃，最低极端气温为−7.9 ℃，年日照时间1 740.7 h，平均霜冻期为95 d。年降水量1 768 mm，年平均相对湿度为82%。气候条件适合杉木生长。

试验林使用1年生苗木于1982年造林，完全随机区组设计，分5种造林密度：A：2 m × 3 m（1 667 株·hm⁻²），B：2 m × 1.5 m（3 333 株·hm⁻²），C：2 m × 1 m（5 000 株·hm⁻²），D：1 m × 1.5 m（6 667 株·hm⁻²）和E：1 m × 1 m（10 000 株·hm⁻²）。每个样地大小为20 m × 30 m，每种造林密度均重复3次，总计15个样地。标记了一共4 800棵树，在冬季测量树高和胸径，从1984至1990年，每年进行1次测量；从1992年至2010年，每隔1或者2年进行1次测量。在每个样地中，选取最高的6株树，并计算其平均树高作为林分优势高。具体统计数据见表1。研究数据随机抽取60%用于模型构建，剩余40%用于模型验证。

利用常用的5种林分密度指标，分别是每公顷株数密度N、林分密度SDI指数、SD指数、优势高—营养面积比Z指数和相对植距RS指数构建杉木林分蓄积量生长模型。

　（1） SDI指数

SDI指数是Reineke基于完满立木度提出的一种林分密度指标^[1]。其优点是不仅能反映林地内的株数，而且能反映林木的大小，同时不受林龄和立地条件的影响，具有容易测算、使用简单的特点^[11]。其数学表达式为：

式中：Dq为林分平方平均胸径cm；D0为林分标准直径，杉木一般取值20 cm；β为Reinek自稀疏系数。Sun等^[12]认为杉木同龄林中林木的自然死亡率达到2%时才开始启动林分自稀疏，因此在剔除自稀疏死亡率小于2%的样本后，根据Reineke的自稀疏理论，对单位面积株数N和林分平均胸径Dq进行对数转化，即lnN和lnDq计算杉木自然稀疏线斜率β。据此分析，计算得到自然稀疏斜率β=−1.897 2。

（2） SD指数

根据Nilson^[2]对林木之间的平均距离L和平方平均胸径Dq的定义，林分相对密度SD指数的公式如下：

式中：a，b为参数，k=−a，D₀取值20 cm。

（3）优势高—营养面积比Z指数

张连金等^[3]采用单位面积内株数与优势木的平均高来反映林分中林木的相对密度的大小。其表达式如下：

式中： H_d为优势木平均高。该指数与林龄、立地指数等无关，计算相对较为简单。

（4） 相对植距RS指数

相对植距包含了林分优势高和每公顷株数，其表达式如下^[11]：

本研究选取了优势高-营养面积比Z指数、每公顷株数密度N、Reineke密度指数SDI、Nilson密度SD指数和相对植距RS构建杉木林分蓄积量模型；

式子中：M₁，M₂分别为第1期和第2期的林分蓄积量；

A₁为第1期的林分年龄；

Hd₁第1期的林分平均优势高；

K₁为第1期时的5种林分密度指标，包括：优势高-营养面积比Z指数，株数N，林分密度指数SDI，林分密度SD指数，相对植距RS指数；

x₁～x₄,α₁～α₃为模型参数；

z₁～z₄为5种不同造林密度设置的哑变量(造林密度1 667 株·hm⁻²：z₁=1, z₂=z₃=z₄=0；造林密度3 333 株·hm⁻²：z₂=1, z₁=z₃=z₄=0；造林密度5 000 株·hm⁻²：z₃=1, z₁=z₂=z₄=0；造林密度6 667 株·hm⁻²：z₄=1, z₁=z₂=z₃=0；造林密度10 000 株·hm⁻²：z₁=z₂=z₃=z₄=0)。

由于调查的间隔并不都是1年1次，因此为了提高模型的预测精度，可变生长率法被引入到本研究中。可变生长率法考虑间隔期内，林分优势高的变化引起的林分蓄积量生长的变化。杉木林分蓄积量年生长模型利用递归的方式推导过程如下：

(t + 1)年时:

(t + 2)年时:

(t + q)年时:

式中，t为调查年，q为调查间隔期。

研究中构建的杉木优势高模型形式如下：

其中β₁、β₂为待估参数。

本研究中, 杉木林分蓄积量模型的参数估计均利用 SAS中非线性回归模块来完成。利用可变生长率法估计模型参数时, 应用了循环运算。

林分优势高模型和林分蓄积生长模型通过统计量平均绝对偏差(MAD)、均方根误差( RMSE)和决定系数(R²) 进行评价。

杉木林分优势高模型的参数估计结果以及模型的均方根误差、决定系数见表2。结果表明，模型的参数估计值均为有效值，并且模型的决定系数很高，达到了0.953 9，均方根误差RMSE较小，为1.343 1。因此模型估计得到的林分优势高是可信的，可以用于构建杉木林分蓄积量模型。

分别用5种不同林分密度指标，利用可变生长率法建立了6种杉木林分蓄积量年生长模型，其中包括一种不含密度指标的生长模型。参数估计和标准误差及模型评价见表3。从表3中可以看到，5类含有密度指标的模型决定系数R²均在0.979以上，模型精确度较高，且大于不含密度指标的对照组模型（模型决定系数为0.972 8）。林分密度对于杉木蓄积量具有不可忽视的影响，因此将林分密度指标引入杉木林分蓄积量模型可以提升模型精确度，优化模型表现。

从模型决定系数方面来看，所有模型的R²数值由高到低顺序为：包含每公顷株数N的林分蓄积量模型（0.979 9）、相对植距模型（0.979 9）、林分密度指数SDI模型（0.979 4）、优势高营养面积比Z模型（0.979 3）、Nilson密度指数模型（0.979 0）以及不含密度指标模型（0.972 8）。除去因为参数估计异常而被舍弃的2号和3号模型，杉木林分蓄积量模型中表现最好的是以林分密度指数SDI为密度指标的模型。同时，SDI模型也取得了最小的平均绝对偏差和均方根误差。因此，从本研究结果中可以认为，林分密度指数SDI模型是估计杉木林分蓄积量的最佳模型。而不含密度指数模型则是表现最差的林分蓄积量模型，这也说明了在拟合和估计杉木林分蓄积量时，缺少密度指数将会影响模型的准确度。SD指数模型也因为在所有密度指数模型中表现最差，在估计杉木林分蓄积方面成为最不适合的密度指标。我们还发现，根据哑变量参数的估计值可得，在低造林密度（1 667～3 333 株·hm⁻²），林分蓄积生长量要大于中高造林密度的林分（5 000～10 000 株·hm⁻²）。

基于研究得到的预估杉木林分蓄积量的最佳模型，SDI模型，得到了模型残差分布图（图1）。从图1可以看到，都倾向于低估杉木的林分蓄积量。在林分蓄积量处于0～100 m³·hm⁻²的时候，残差点基本分布于0～−20的范围内；当实际林分蓄积量达到100 m³·hm⁻²以上时，残差点分布范围则有所扩大，基本上均匀分布在−10～−50 m³·hm⁻²的区间中。当然在实际林分蓄积量处于1～10 m³·hm⁻²和较大蓄积的时候，模型也有少量高估蓄积量的情况发生。

Figure 1.  Residual distribution of stand volume growth model including density index SDI

林分密度是林分生长模型预估的一个重要指标。不同的密度指标其应用的范围不一样，因此也受到了不同学者的争议。冉啟香等^[13]在对北京地区油松（Pinus tabuliformis Carr.）林分蓄积生长模型研究的过程中指出，选用不同的密度指标构建模型会直接影响模型的预估效果，同时对每公顷胸高断面积（BA）、每公顷株数（N）、林分密度指数（SDI）等不同林分密度指标进行比较，最终得出用每公顷胸高断面积（BA）作为密度指标预估林分蓄积量生长模型时效果更好。张雄清等^[14]利用林分株数密度N作为密度指标构建杉木林分蓄积量模型，并且模型的精度也较高。此外，许多学者在研究林分生长时也经常将林分密度指数SDI当作较好的密度指标纳入模型。吴宏炜等^[15]基于福建湿地松（Pinus elliottii Engelm.）人工林调查数据，选取林分密度指数SDI构建湿地松生长模型，取得了较好的模型精确度。由此可以看出，选取合适的林分密度指标构建林分蓄积量模型很有必要。

但是对于SDI指数，争议较多的还是自然稀疏线斜率是否为−1.605。多数学者认为自然稀疏线斜率不是固定不变−1.605，而是随着树种、初始密度、立地以及气候的变化而变化^[16-19]。Pretzsch和Biber^[20]在构建挪威云杉（Picea abies (L.) Karst.）林分生长收获表时发现，用固定的−1.605斜率会导致SDI指数偏差，进而影响收获表的准确性。Zhang 等^[21]基于气候敏感的林分密度指数SDI构建了杉木林分生长模型，发现该SDI指数很好地反映了密度的变化情况，使得林分生长模型精度得到提高。总之，基于我们的研究，对于杉木林分蓄积量建模预测，结果显示林分密度指数SDI为最佳的密度指标。同时，密度效应在杉木林分蓄积量生长中起到了很重要的作用，不含有密度指标的模型将在精度和预测能力上有所降低。

此外，对于定期（不定期）的调查数据，如何分析并构建年生长量模型对于精准预测很关键。近年来，基于已有的连续调查数据，通常用固定生长率法，其前提是假设在整个生长期，林分的年生长量不变，常用定期平均生长量来代替年生长量^[11]。很明显该方法并不符合林分的实际生长规律。因为随着森林的演替过程，一些林分的特征（立地，年龄，林分竞争等）都会随着林龄变化而变化，因此也必然会导致林分的年生长量也发生变化。之后Ochi和Cao^[22]提出利用可变生长率法建立林分年生长量相容性模型，得到了很好的拟合结果。之后张雄清等^[23]基于北京山区油松定期调查数据，运用可变生长率法构建了北京油松全林分年生长模型。模型结果显示，利用可变生长率法构建的油松生长模型更加符合其生长规律，同时解决了模型在预测油松生长阶段中的无偏性，提升了模型预测性能。

基于可变生长率法，建立了含5种林分密度指标的杉木林分蓄积年生长模型，含有密度指标的模型决定系数R²均在0.979以上，精度高于不含密度指标的对照组模型。从模型的参数估计和模型精度上来看，精度最高的是包含每公顷株数密度N和相对植距RS的林分蓄积量模型。但由于这两个模型中的密度指标的参数估计不显著，模型表现不稳定，因此不适宜预测。而包含林分密度指数SDI的蓄积量模型相较于优势高营养面积比Z、Nilson密度SD指数和不含有密度指标的蓄积量模型，取得了最大的R²和最小的均方根误差和平均绝对误差，因此SDI指数在本研究中是建立杉木林分蓄积量模型中最好的密度指标。其次，本研究还发现在低造林密度（1 667～3 333 株·hm⁻²），林分蓄积生长量要大于中高造林密度的林分（5 000～10 000 株·hm⁻²）。