浅析理查德方程拟合树木高生长模型及参数求定

向南海

（旌德县蔡家桥林场，安徽 宣城 242600）

浅析理查德方程拟合树木高生长模型及参数求定

向南海

（旌德县蔡家桥林场，安徽 宣城 242600）

在林学方面，描述树木及林分生长过程时，理查德方程是近代应用最为广泛、适应性较强的一类生长曲线方程。笔者就速生树种杉木进入中龄林后用理查德方程拟合树木高生长模型并求定参数，在此与大家做一交流。

理查德方程；参数

1 树木生产量测定意义和树木生长方程概念

1.1 树木生长量测定意义

找出树木生长规律测定树木生长量，在森林经营管理上有着很重要的意义，它既反映立地条件的好坏，又可以作为判断营林效果和森林生产能力以及确定年伐量和主伐年龄的重要依据。

在现代林业特别是集体林权制度深化改革中，森林资源作为资产纳入资本运作，给林业发展增添了新的活力，且森林资源资产评估工作已经作为必不可少的工作环节纳入多个经济领域。在一个森林经理期内，测定林木生长量的预期收获值直接影响森林资产价值评估结果，给林业资本运作也带来一定的牵引作用。

1.2 树木生长方程的基本概念

树木生长方程是描述某树种（组）各调查因子总生长量y（t）随年龄（t）生长变化规律的数学模型。它是该树种某调查因子的平均生长过程，也就是在均值意义上的生长方程。主要特点如下：

当t＝0时y（t）＝0。此条件称之为树木生长方程应满足的初始条件。

y（t）存在一条渐进线y（t）＝A，A是该树木生长极大值。

树木的生长是不可逆的，使得y（t）是关于年龄（t）的单调非减函数。

y（t）是关于t的连续且光滑的函数曲线。

2 理查德方程产生背景与应用范围

2.1 理查德方程产生背景

理查德（Richards）方程是基于著名的Bertalanffy生长理论发展而来。方程的基本形式为：y=A（1±Bexp-kt）1/（1-m）（当m＞1时取“+”号，0＜m＜1时取“-”号）。该方程的参数或参数的组合有明显的生物学意义，而且随着方程中参数m的取值不同，它构成了不同的生物生长或种群动态方程。因此，理查德方程有着广泛的适用性，逻辑性强，参数可由独立的试验加以验证，即参数可作出生物学解释，从理论上对未来的趋势进行预测。其产生背景简述如下：

Bertalanffy通过分析动物的生长，发现在动物生长期间，动物的体重增长速率为同化速率与消耗速率之差，而后两者分别和同化器官的大小以及动物体重成比例，即：

式中：F—同化器官重，W—体重，Ra—同化速率，Rt—消耗速率。

由相对生长关系，有F=rWm，因此：

式中：η=ar

理查德（Richards）认为，m的取值为2/3，从而得到著名的理查德生长方程：

2.2 理查德方程在林业中的应用

（1）在胸径和断面积生长中的应用：林木胸径和断面积生长曲线满足t=t0（生长至1.3 m 所需的年龄），y＝0的初始条件。

（2）在树木高生长中的应用满足当t=0时、y=0的初始条件。

因此，通过引入参数使对树木生长具有广泛的适应能力，主要应用在以下几个方面：1）与林分密度有关的林分断面积生长模型；2）树高生长及地位指数方面研究；3）林分直径结构研究方面。

3 理查德方程应用实例

3.1 拟合对象

拟合对象为杉木中龄林高生长过程。林业调查中，按龄组、龄级划分条件定义杉木中龄林生长阶段为林龄10至20年，5年1个龄级跨Ⅲ和Ⅳ二个龄级。

3.2 拟合模型

将方程两边各取对数得：

式中：t为树木年龄观测值，H为树高，a、k、c为参数，其中k为模拟生长规律时制定的修正值，一般在1.0以内。

3.3 参数求定

3.3.1 参数k值求定

利用计算机应用程序求定参数k值，应用程序环境为：Microsoft Excel 2010。

第一步，列出杉木中龄林年龄与树高观测值于Excel表中。在实际调查中，可根据不同坡向坡位等不同立地条件调查多组数据，这里仅列出三组调查数据如下：

年龄( t ) 1 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0树高( h ) 6 . 5 7 . 6 8 . 2 9 . 0 9 . 9 1 1 . 0 1 1 . 7 1 2 . 3 1 3 . 2 1 4 . 0 1 4 . 6年龄( t ) 1 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0树高( h ) 6 . 4 7 . 5 8 . 0 9 . 1 1 0 . 0 1 0 . 6 1 1 . 3 1 2 . 4 1 3 . 0 1 3 . 8 1 4 . 7年龄( t ) 1 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0树高( h ) 6 . 8 7 . 3 8 . 0 9 . 1 9 . 5 9 . 9 1 0 . 6 1 1 . 4 1 1 . 5 1 1 . 9 1 3 . 0

第二步，绘制树高曲线图。选择年龄与树高观测值数据区域，点击插入“散点图”，然后添加趋势线。趋势线格式选择“多项式”，并勾选“显示公式”和“显示R平方值”，见下图1。图2为树高曲线图绘制结果。

图2 中，公式y=-0.008 2x2+0.987 4x-2.467为Excel自动拟合的树高与年龄的回归方程，R2值是回归方程拟合优度，取值范围是[0，1]，越接近1拟合程度越好，反之越差。

第三步，分别取三个年龄观测值t1、t2、t3，将它们代入回归方程中，求出三个相对应的理论高H1、H2、H3。在这里，t1、t2、t3分别取值为10、15、20，即杉木中龄林的下限、中限、上限三个年龄值。理论高求算结果为6.6、10.5、14。

第四步，根据算出的H1、H2、H3，求观测比。

第五步，由（5）式可得：

因此，观测比λ可改写为：

已知t1＝10、t2＝15、t3＝20，代入（7）式如下：

已知λ＝0.617，令e-10k=x，这里x＞0，且x＜1，则：

第六步，在Excel中规划求解（15）式中x值。方法如下：

首先，装载规划求解加载宏。依次点击“文件－选项－加载宏－管理Excel加载项转到……”，将分析工具库和规划求解加载项一并勾选，见图3，点击“确定”；然后，点击“数据”菜单“规划求解”按钮，弹出对话框进行设置，见图4。

设置目标单元格是要显示（15）式中左边公式内容所在的单元格，目标值选定（15）式中等号右边值，即1。在这里我们已经于规划求解前在$T$8单元格中预先用函数写下了公式“=POWER（1+$T $9，0.617）*（1-$T$9）+POWER（$T$9，3/2）”，式中$T $9是（12）式中变量x引用的单元格绝对位置，也是要规划求解的值，并添加约束规则2项：$T$9＜=1，$T$9＞=0。勾选“使无约束变量为非负数”，选择求解方法为“非线性GRG”。接着点击“选项”按钮，见图5、图6。在所有方法选项卡中，如果想要看到迭代计算结果就勾选“显示迭代结果”，要直接看到计算结果就取消勾选，其他可不用填写；在非线性GRG选项卡中派生选择“中心”，并勾选“使用多初始点”，确定退出。

回到图4界面，点击“求解”按钮，得到全局最优解，见图7对话框。

关闭对话框，回到Excel工作表中，见图8。在$T $9单元格中已经看到我们要的x变量结果：x= 0.62224955。如果最优解得到的是1或0，重复求解过程便可。

第七步，将x值代入定义的e-10k=x式中，求得：k=0.047 441 406。

3.3.2 参数a、c值求定

将年龄观测值t1、t2、t3和对应的用回归方程求解的树高理论值H1、H2、H3以及参数k值一并代入式（7）、（8）、（9）中，在Excel中用函数轻松求解，得到：a=29.918 477 22，c=1.552 514 895。

3.4 杉木中龄林高生长方程模型拟合结果

将上述参数求定值代入理查德方程模型，拟合结果为：

将年龄观测值代入公式验证：当t1=10时，解得H1＝6.6。

当t2=15时，解得H2＝10.5。

当t3=20时，解得H3＝14.0。验证通过。

4 理查德方程应用注意事项

一般的直线方程或曲线方程有a、b两个待定参数，而理查德方程 H=a(1-e-kt)c中有三个待定参数a、k、c，问题显然复杂。我们在上述实例中，设定k值为模拟生长规律时制定的修正值，单独求算出，然后再代入方程求算出其他两个参数，这样便于求解。在这里，需要特别注意的是：k值模拟树木生长规律求算时，要尽量满足树木生长函数曲线要求，缩短拟合区间范围，以各龄组分段拟合为最优。

S711

A

2095-0152（2015）03-0071-04

2015-03-08

2015-05-07

向南海(1967- )，男，林业工程师，主要从事林业调查规划设计和森林资源资产评估咨询工作。E-mail：376141298@qq.com

杨婷婷）