三元立木材积模型的初步研究

苏振海

(辽宁省林业调查规划监测院,辽宁 沈阳 110022)

当今森林资源各类调查监测项目对调查精度要求越来越高,而立木材积则是其中最重要的指标之一。现在各地主要使用的是一元、二元材积表或材积式,虽然可以满足日常生产、监测需求,但精度由于受参与计算的变量的数量限制,无论采取何种材积式,其精度都不能有本质上的提高,因此采用三元材积模型是较为可行的。原本意义上的三元立木材积表也是林业调查中最常用的数表,是以胸径、树高以及干形某一指标作为分级因子而编制的立木材积表[1]。干形指标多采用形数、形率、直径差表示。原本三元材积表多用于测算林分蓄积量,简单说就是先测定干形分级因子,采用的还是二元立木材积表程序进行测算。也就是说,其根本还是二元材积表,只不过是多了一个分级因子,而且多采用平均干形因子确定,对单株立木估算仍存在系统误差。而本文所提到的三元立木材积模型,是指通过测量立木上部一个直径,通过材积模型计算,给出更准确的单株立木材积。

1 模型建立目标

建立模型的目的是提高估测精度；模型要有实际的意义,即模型要与立木材积有实际的联系；模型参数尽量不要太多,避免出现病态模型,即模型稳定性要好。

2 样本组织

辽宁柞树共采集1150株,其中集体392株,国有758株；柞树样本分布在抚顺市、本溪市、丹东市、铁岭市和鞍山市。样本采集后,采用散点图法和三倍标准差法剔除异常值。

柞树编表样木数据834个,剔除异常数据8个,剩余826个；柞树检验样木数据316个,剔除异常数据15个,剩余301个。

3 建模基本方程选择

由于笔者经验有限,因此只能在原有二元立木材积式的基础上改编得到三元立木材积式。二元立木材积经常使用山本式、奥盖尔、动态模型、动态改编模型、孟宪宇、斯泊尔、纳斯伦德、迈耶式[2]。三元材积模型要增加一个变量,因而会造成参数的增加,原有二元材积式中,动态模型、动态改编模型、孟宪宇、纳斯伦德、迈耶式原有参数就是5～9个,改编后可能增加到8～14个,不利于模型稳定性,因此只选择了山本式、奥盖尔、斯泊尔3个参数较少的二元材积式进行改编。

山本式：

V=C0DC1HC2

奥盖尔：

V=D2(C0+C1H)

斯泊尔：

V=C0+C1D2H

4 新增变量的选择

树干上部直径可分为固定高度和相对位置高度,固定位置如树干4m、7m位置[3],相对位置如主干的1/2、1/10、枝下高。笔者以为,固定位置比较受限制,如树高不足7m,或树高过高,4m位置代表性又不够强,而枝下高,就单株树木而言变化范围较大,代表性不稳定,因此采用相对位置作为新增变量。而相对位置有多个,从1/10～9/10,具体采用哪个可采用主因子分析确定,方法如下。

使用ForSatat统计之林软件的因子分析中的主成分分析工具,选入树高H、胸径及树高的1/10～9/10作为变量,使用相关矩阵进行标准化,避免量纲不同造成的影响。由于本文是建立三元材积模型,因此除胸径和树高外,只需要选择一个因子参与模型编制,所以主分量选择大于3个即可。通过因子负荷量来确定相对位置高度,分析结果如表1。

表1 因子负荷量表

从表1可知，主成分1负荷量中,除树高H、胸径外,树干直径D0.1～D0.9对描述树干都有很好的效果,但D0.5处更好。

除数据分析结果外,还要考虑实际测量难度,树冠高度直接影响主干上部直径的测量,越是向上,测量难度越大。综合以上因素,选择中央直径作为新增变量最为适合。

5 三元立木材积模型编制

依据前文3中选择的基本模型改编如下模型:

V=C1D1.3C2HC3+C4D0.5C5HC6

(1)

V=C1(D1.32+D0.52)H+C2D0.5C3HC4

(2)

V=(C1D1.32+C2D0.52)H+C3

(3)

V=C1(D0.5/D1.3)C2D1.3C3HC4

(4)

V=C1(D1.3C2+D0.5C2)HC3+C4D0.52H

(5)

V=C1(C2D1.3C3+D0.5C3)HC4

(6)

V=C1D1.3C2H(D0.5/D1.3)C3

(7)

V=D1.32(C1+C2H)(D0.5/D1.3)C3+C4

(8)

V=D1.32(C1+C2H)+D0.52(C3+C4H)+C5

(9)

式中,D1.3为胸径；D0.5为树干1/2直径；H为树高；C*为模型参数。

式(1)参考山本式,相当于将树干分为2段,分别用山本式进行拟合求材积；式(2)参考山本式,相当于树干分为2段,第1段(0m至树干中央高度,下同)采用2个位置断面积平均值作为断面积,第2段(树干中央高度至梢头,下同)用山本式求算；式(3)参考斯泊尔,相当于将树干分为2段,分别用斯泊尔式进行拟合求材积；式(4)参考山本式,采用D0.5与D1.3的比值作为调整系数；式(5)参考斯泊尔,相当于树干分为2段,类似于式(2),调整了参数与常数的位置；式(6)参考山本式,在式(1)基础上,统一了2个直径指数位置参数和树高的参数,增加1个胸径调整参数；式(7)在式(4)基础上减少H的参数；式(8)参考奥盖尔,采用D0.5与D1.3的比值作为调整系数；式(9)参考奥盖尔,相当于将树干分为2段,分别用奥盖尔式进行拟合求材积。

6 模型参数计算

将柞树编表样木(去掉异常值)826个数据代入上文模型公式,用麦夸特迭代方法计算各模型参数值。

表2 各模型参数参考值表

7 建模方程评价

评价回归模型的准则在统计文献中常见的有残差平方和Q、复相关系数R[4]、参数变动系数(稳定性)等,其中参数稳定性是评价通用性立木材积预测模型极为重要的准则。

7.1 各模型参数变动系数

一般情况下,参数变动系数超过50%模型就不稳定。因此式(1)、式(5)、式(9)不够稳定。

7.2 回归方程适应性检验

用样本实测材积V与样本回归估计值V’=f(D1.3,D0.5,H)比较。如果二者之间呈直线回归方程y=a+bx的系数(a,b)与(0,1)没有显著差异,则可以判断原回归方程f的选型是正确的[5]。

通过表4看出,当模型F

表4 模拟方程选优对照表

7.3 总相对误差

为检验选定模型的适用精度,除编表样本外另选取了检验样本301个；检验方法用总相对误差(RS),RS应介于±3%。二元山本式RS=1.12%,模型(6)RS=0.65%,总相对误差减少41.9%。

表3 各模型参数变动系数表

8 小结

本文通过分析验证,给出4个三元(胸径、树高、树干中央直径)立木材积式,分别是(3)、(4)、(6)、(7)。其中，式(3)模型参数少,精度和稳定性不及同样为3参数的式(7)；式(4)精度较高,稳定性好；式(6)精度最高,稳定性较好。因此,选择3参数,式(7)最优；选则4参数,式(6)精度更好。

本文样本测量值均为围径,而实际测量过程中,往往树干上部直径的测量值均为轮径,应建立围径与轮径的回归方程带入模型计算。树干上部直径测量方法及工具对测量精度各有不同,本文未进行分析,但其必然会对模型使用精度产生影响。