CN103017696A - 一种基于曲线模拟测定单木材积量的方法 - Google Patents

Abstract

本发明公布了一种基于曲线模拟测定单木材积量的方法。包括如下步骤：首先，在树干的0.5米、1.3米、2.5米处，沿3个方向(南北、东北至西南、西北至东南)钻取3根树芯；然后，测定树芯型值点至髓心的宽度值，确定型值点的二维坐标；其次，将型值点坐标代入三次参数样条曲线方程，模拟树干横断面曲线；再次，根据3个高度的横断面曲线以及树高信息，利用抛物样条曲线方程模拟树干纵断面曲线，形成树干立体形状；最后，将树干平分为9000段，分别计算每一段的体积并累加得出单木材积量。

Description

一种基于曲线模拟测定单木材积量的方法

技术领域

本发明涉及一种测量单株树木体积的新方法，该方法通过计算机模拟出的树干曲线进行树木体积的计算。

背景技术

树木材积的测量在林业***有着广泛的应用，例如林地生产量的计算、原木使用率和出材率的计算、立木的生长监测等。传统的材积测量方法主要有每木检尺、角规测树、全站仪测树等，这些方法假设树干横断面是规则圆形或者椭圆形，然后通过测量树高、胸径来计算树木的材积量。然而在实际测量中，树干横断面并不是规则的，因此传统测量法得到的材积量与实际材积量之间会存在一定的误差，使得测量结果不够准确。

发明内容

本发明的设计分为3个部分，具体技术方案如下：

第一部分是树干横断面曲线的模拟：

首先，在树干0.5米、1.3米、2.5米3个高度上各钻取3根树芯，获取到3组6个不同方位的型值点，型值点在树干曲线上的分布如图1所示；

然后，选择三次参数样条曲线方程进行树干横断面曲线的模拟，曲线方程表达式为 $P_i\left(t\right)=P_i+P_i^{\′}t+[\frac{3(P_{i+1}-P_i)}{{t_{i+1}}^2}-\frac{2P_i^{\′}}{t_{i+1}}-\frac{{P_{i+1}}^{\′}}{t_{i+1}}]t^2+[\frac{2(P_i-P_{i+1})}{{t_{i+1}}^3}-\frac{P_i^{\′}}{{t_{i+1}}^2}-\frac{{P_{i+1}}^{\′}}{{t_{i+1}}^2}]t^3,$ 式中P_i为第i个型值点；P_i′是在第i个型值点处的切向矢量；t是自变量表示分段曲线的弦长，0≤t≤t_i+1；i＝1，2，…，n-1；t_i+1为第i+1个型值点处的弦长；

其次，设定曲线的连接条件使得每条分段曲线能够首尾相连，且保证连接后的整条曲线是光滑的，本发明设定的条件是每段曲线在端点连接处二阶连续，得到以下矩阵方程： $\left(t_{i+2}2(t_{i+2}+t_{i+1})t_{i+1}\right)\left(\begin{array}{c}{P}_i^{\′}\\ {P_{i+1}}^{\′}\\ {P_{i+2}}^{\′}\end{array}\right)=3[\frac{t_{i+1}}{t_{i+2}}(P_{i+2}-P_{i-1})+\frac{t_{i+2}}{t_{i+1}}(P_{i+1}-P_i)],$ 对于由n个型值点确定的曲线，可以得到以下方程组，共n-2个方程，有n个未知数(型值点切矢量) $\left(\begin{array}{ccccccc}{t}_3& 2(t_3+t_2)& t_2& 0& 0& \·\·\·& 0\\ 0& t_4& 2(t_4+t_3)& t_3& 0& \·\·\·& 0\\ 0& 0& t_5& 2(t_5+t_4)& t_4& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& 0& t_n& 2(t_n+t_{n-1})& t_{n-1}\end{array}\right)\left(\begin{array}{c}{p}_1^{\′}\\ p_2^{\′}\\ p_3^{\′}\\ \·\\ \·\\ \·\\ p_n^{\′}\end{array}\right)=\left(\begin{array}{c}3[\frac{t_2}{t_3}(p_3-p_2)+\frac{t_3}{t_2}(p_2-p_1)]\\ 3[\frac{t_3}{t_4}(p_4-p_3)+\frac{t_4}{t_3}(p_3-p_2)]\\ 3[\frac{t_4}{t_5}(p_5-p_4)+\frac{t_5}{t_4}(p_4-p_3)\\ \·\\ \·\\ \·\\ 3[\frac{t_{n-1}}{t_n}(p_n-p_{n-1})+\frac{t_n}{t_{n-1}}(p_{n-1}-p_{n-2})]\end{array}\right)$ 通过解算以上方程可以求得P₂′，P₃′，...，P_n+1′的取值；

再次，设定端点条件，即闭合曲线的首末端点二阶导数相同，本发明假设该二阶导数为0，即得到以下运算方程

和

可以求得P₁′，P_n′，由此n个型值点处的切向矢量都被计算出来，曲线方程系数即被确定；

最后，按照以上方法，利用0.5米、1.3米、2.5米3个高度的型值点数据模拟出树干3个高度的横断面曲线，如图2所示。

第二部分是树干纵断面曲线的模拟：

首先，选择抛物样条曲线进行树干纵断面曲线的模拟，曲线方程表达式为：P(t)＝A₁+A₂t+A₃t²，式中，A₁、A₂、A₃分别为树干立体曲线段拟合的3个待定系数，0≤t≤1；

接着，建立抛物样条曲线的3个独立条件：(1)设抛物线过P_i、P_i+1、P_i+2三个点，P_i为始点，当参变量t＝0时，曲线过P_i点；(2)P_i+2为终点，当参变量t＝1时，曲线过P_i+2点；(3)当参变量t＝0.5时，曲线过P_i+1点，且切矢量平行于P_i+2-P_i。则有：

$S_i\left(t_i\right)=(2t_i^2-3t_i+1)P_i+(4t_i-4t_i^2)P_{i+1}+(2t_i^2-t_i)P_{i+2};$

然后，通过上式得出的曲线段集合在每两段曲线之间的搭接区间，两条抛物线是不会自然重叠成一条曲线的，为了把两条曲线用一条光滑的曲线连接，需要加权合成，即选择两个合适的权函数，分别为f(T)＝1-T；g(T)＝T(0≤T≤1)，加权合成后的曲线为：P_i+1(t)＝f(T)S_i(t_i)+g(T)S_i+1(t_i+1)。模拟中设t的取值范围为：0≤t≤0.5，则T＝2t；t_i＝0.5+t；t_i+1＝t；则合成曲线(两条抛物线的重叠段)可以写成：P_i+1(t)＝(-4t³+4t²-t)P_i+(12t³-10t²+1)P_i+1+(-12t³+8t²+t)P_i+2+(4t³-2t²)P_i+3；

其次，建立端点条件。n个型值点只能得到n-3段曲线，但n个型值点之间应有n-1个区段，设计中添加两个端点P₀和P_n+1，且满足P₀＝P₁；P_n+1＝P_n。联合前面的n-3段曲线就可得到n-1个曲线段。这样利用上述的抛物样条方程即可模拟树干曲线。

再次，利用已模拟的树干横断面曲线，量测树的胸径，根据树高-胸径关系式计算树高坐标，并从3个高度的树干横断面曲线中选取3个对应的型值点，保证3个型值点在同一树干纵断面曲线上，从而形成抛物样条曲线的4个型值点，代入抛物样条曲线方程模拟树干纵断面曲线；

最后，将每个横断面曲线按角度平均分成1000份，按照上一步方法，形成1000条树干纵断面曲线，从而模拟出整个树干的立体形状，如图3所示。

第三部分是树干材积量的计算：

首先，沿树干方向将树干平均分成9000段，如图4所示；

然后，计算每个分割段的下底面面积。由于树干立体曲线是由1000条树干纵断面曲线构成，因此每个分割段实际也被分割成1000份，计算出所有1000份的面积并累加，即可得到下底面面积，如图5所示；

最后，由于完整树干被平均分割成9000段，因此在计算体积时可以假设每一小段的上、下底面面积是近似相等的，体积计算就是其中一个底面积与高的乘积。通过步骤2计算出的底面积与分割段的高度可以得到每个分割段的体积，累加9000段体积即可得到整个树干的体积，即树干材积量。

与现有技术相比，本发明具有如下优点：

(1)本发明将模拟树干真实形状作为主要的研究对象，使得测量结果相对于将树干横断面近似为圆形或者椭圆形的传统方法更准确；

(2)本发明方法以计算机处理为主，人工测量为辅，因此更加智能，大大降低了人的劳动强度，提高了人的工作效率；

(3)本发明方法为了获取树干型值点数据，只用在特定位置钻取树木树芯即可，不会破坏活立木的生长；

(4)本发明适合于杉木、湿地松、马尾松等多林种材积量的测定，具有非常广泛的应用价值。

附图说明

图1为型值点方位角度示意图；

图2为在0.5米、1.3米和2.5米3个高度的横断面曲线示意图；

图3为通过树干纵断面曲线模拟树干立体形状的示意图；

图4为将树干平分为9000份的示意图；

图5为将树干底面平分为10000份计算底面积的示意图。

具体实施方式

实施例1：树干曲线型值点的提取

首先，选取森林样方，并用测绳将其围起，进行每木检尺，记录树种名称和胸径值，并同时记录样方名称、采样日期、坡度、坡向、经纬度、海拔等信息。

然后，在样方中选取典型的树干，进行树芯数据的采集。分别在距地面0.5米、1.3米、2.5米处各钻取3根生长树芯，每木采用相同取样策略进行取样，即在南北、东北至西南、西北至东南方向各取一芯。树芯采集完成后，用树枝条将钻孔填塞严实。

最后，将样芯立即装入纸吸管内并在吸管上用油性笔标注代码，待样芯干燥后，用乳胶将其粘在样本木槽上。为避免在乳胶干燥过程中发生样芯变形，用细线绳将其按十字形式捆绑固定，待乳胶干后方可拆线。然后用120、320、600目的砂纸依次对样芯进行打磨，直到在显微镜下能够清晰地分辨细胞的大小。以靠近树皮的部分为采样年份，从***向髓心进行定年。用树轮宽度量测仪(LINTAB)对树木样芯进行宽度量测，量测精度为0.01mm。

实施例2：树干横断面曲线的模拟

首先，根据每个型值点P_i至髓心的宽度值及型值点所在的方向角，计算出以髓心为坐标原点的型值点二维坐标值；

序号	树芯高度(m)	型值点方向角(°)	至髓心宽度(cm)	X坐标	Y坐标
						1	0.5	-90	3.584	0	-3.584
2	0.5	-30	3.8905	3.369272	-1.94525
						3	0.5	30	3.4165	2.958776	1.70825
4	0.5	90	3.774	0	3.774
						5	0.5	150	3.462	-2.99818	1.731
6	0.5	-150	3.654	-3.16446	-1.827
						7	1.3	-90	3.0445	0	-3.0445
8	1.3	-30	2.838	2.45778	-1.419
						9	1.3	30	2.858	2.475101	1.429
10	1.3	90	2.913	0	2.913
						11	1.3	150	2.954	-2.55824	1.477
12	1.3	-150	3.024	-2.61886	-1.512
						13	2.5	-90	1.6635	0	-1.6635
14	2.5	-30	1.765	1.528535	-0.8825
						15	2.5	30	1.799	1.55798	0.8995
16	2.5	90	2.1515	0	2.1515
						17	2.5	150	1.6635	-1.44063	0.83175
18	2.5	-150	1.837	-1.59089	-0.9185

然后，建立模拟树干横断面曲线的三次参数样条曲线方程，方程表达式为 $P_i\left(t\right)=P_i+P_i^{\′}t+[\frac{3(P_{i+1}-P_i)}{{t_{i+1}}^2}-\frac{2P_i^{\′}}{t_{i+1}}-\frac{{P_{i+1}}^{\′}}{t_{i+1}}]t^2+[\frac{2(P_i-P_{i+1})}{{t_{i+1}}^3}-\frac{P_i^{\′}}{{t_{i+1}}^2}-\frac{{P_{i+1}}^{\′}}{{t_{i+1}}^2}]t^3,$ 式中0≤t≤t_i+1；i＝1，2，…，n-1；取t为分段曲线的弦长。

其次，为了使模拟的曲线能够平滑过渡，设定相应的端点条件和连接条件，得到 $\left(t_{i+2}2(t_{i+2}+t_{i+1})t_{i+1}\right)\left(\begin{array}{c}{P}_i^{\′}\\ {P_{i+1}}^{\′}\\ {P_{i+2}}^{\′}\end{array}\right)=3[\frac{t_{i+1}}{t_{i+2}}(P_{i+2}-P_{i-1})+\frac{t_{i+2}}{t_{i+1}}(P_{i+1}-P_i)],$ $2P_1^{\′}+P_2^{\′}=\frac{3(P_2-P_1)}{t_2}$ 和

的方程组，将型值点坐标信息代入以上方程组，求出P₁′，P₂′，P₃′，...，P_n+1′，P_n′等变量，从而求解出三次参数样条曲线的系数；

再次，对于每个三次样条曲线段，在程序中设定把每个曲线段的弦长分割成9000份，即在每段曲线中***8999个插值点，根据求解出的三次参数样条曲线方程计算各插值点的坐标值，然后连线绘图，绘制出本段样条曲线；

最后，绘出全部分段曲线，即完成了一条树木横断面曲线的模拟。根据在树的0.5米、1.3米和2.5米处的型值点数据，模拟出3条树干横断面曲线。

实例3：树干纵断面曲线的模拟

首先，在3个不同高度上的横断面曲线上选择3个对应的型值点，保证这3个点在树干的同一纵断面上；

然后，根据型值点数据测量树的胸径值，根据树高-胸径关系式计算树干的高度；

其次，建立模拟树干纵断面曲线的抛物样条曲线方程，其表达式为 $S_i\left(t_i\right)=(2t_i^2-3t_i+1)P_i+(4t_i-4t_i^2)P_{i+1}+(2t_i^2-t_i)P_{i+2}.$ 为了使n个型值点用一条光滑曲线连接，本发明选择一种加权合成方法，对曲线方程进行改造，使曲线方程变为P_i+1(t)＝(-4t³+4t²-t)P_i+(12t³-10t²+1)P_i+1+(-12t³+8t²+t)P_i+2+(4t³-2t²)P_i+3

再次，将第一步获取的3个型值点和第二步获取的型值点组成抛物样条曲线的4个型值点，将4个点的坐标代入第三步方程，得到抛物样条曲线方程

最后，将3条横断面曲线分别按照角度平均分成1000份，得到每条横断面曲线上的1000个型值点，从而在纵向上形成1000组抛物样条型值点集合，按照以上操作步骤模拟出1000条树干纵断面曲线，形成树干的立体形状。

Claims (5)

1.一种基于曲线模拟测定单木材积量的方法，包括以下步骤：

第一步，树干横断面曲线的模拟：

首先，在树干0.5米、1.3米、2.5米3个高度上各钻取3根树芯，获取到3组6个不同方位的型值点；

然后，选择三次参数样条曲线方程进行树干横断面曲线的模拟，曲线方程表达式为 $P_i\left(t\right)=P_i+P_i^{\′}t+[\frac{3(P_{i+1}-P_i)}{{t_{i+1}}^2}-\frac{2P_i^{\′}}{t_{i+1}}-\frac{{P_{i+1}}^{\′}}{t_{i+1}}]t^2+[\frac{2(P_i-P_{i+1})}{{t_{i+1}}^3}-\frac{P_i^{\′}}{{t_{i+1}}^2}-\frac{{P_{i+1}}^{\′}}{{t_{i+1}}^2}]t^3,$ 式中P_i为第i个型值点；P_i′是在第i个型值点处的切向矢量；t是自变量表示分段曲线的弦长，0≤t≤t_i+1；i＝1，2，…，n-1；t_i+1为第i+1个型值点处的弦长；

其次，设定曲线的连接条件使得每条分段曲线能够首尾相连，且保证连接后的整条曲线是光滑的，本发明设定的条件是每段曲线在端点连接处二阶连续，得到以下矩阵方程： $\left(t_{i+2}2(t_{i+2}+t_{i+1})t_{i+1}\right)\left(\begin{array}{c}{P}_i^{\′}\\ {P_{i+1}}^{\′}\\ {P_{i+2}}^{\′}\end{array}\right)=3[\frac{t_{i+1}}{t_{i+2}}(P_{i+2}-P_{i-1})+\frac{t_{i+2}}{t_{i+1}}(P_{i+1}-P_i)],$ 通过解算方程可以求得P₂′，P₃′，…，P_n+1′的取值；

再次，设定端点条件，即闭合曲线的首末端点二阶导数相同，本发明假设该二阶导数为0，即得到以下运算方程和可以求得P₁′，P_n′，由此每个型值点处的切向矢量都被计算出来，曲线方程系数即被确定；

最后，按照以上方法，利用0.5米、1.3米、2.5米3个高度的型值点数据代入曲线方程模拟出树干3个高度的横断面曲线；

第二步，树干纵断面曲线的模拟：

首先，选择抛物样条曲线进行树干纵断面曲线的模拟，曲线方程表达式为：P(t)＝A₁+A₂t+A₃t²，式中，A₁、A₂、A₃分别为树干立体曲线段拟合的3个待定系数，0≤t≤1；

接着，建立抛物样条曲线的3个独立条件：(1)设抛物线过P_i、P_i+1、P_i+2三个点，P_i为始点，当参变量t＝0时，曲线过P_i点；(2)P_i+2为终点，当参变量t＝1时，曲线过P_i+2点；(3)当参变量t＝0.5时，曲线过P_i+1点，且切矢量平行于P_i+2-P_i，则有：

$S_i\left(t_i\right)=(2t_i^2-3t_i+1)P_i+(4t_i-4t_i^2)P_{i+1}+(2t_i^2-t_i)P_{i+2};$

然后，通过上式得出的曲线段集合在每两段曲线之间的搭接区间，两条抛物线是不会自然重叠成一条曲线的，为了把两条曲线用一条光滑的曲线连接，需要加权合成，即选择两个合适的权函数，分别为f(T)＝1-T；g(T)＝T(0≤T≤1)，加权合成后的曲线为：P_i+1(t)＝f(T)S_i(t_i)+g(T)S_i+1(t_i+1)，模拟中设t的取值范围为：0≤t≤0.5，则T＝2t；t_i＝0.5+t；t_i+1＝t；则合成曲线(两条抛物线的重叠段)可以写成：P_i+1(t)＝(-4t³+4t²-t)P_i+(12t³-10t²+1)P_i+1+(-12t³+8t²+t)P_i+2+(4t³-2t²)P_i+3；

其次，建立端点条件，n个型值点只能得到n-3段曲线，但n个型值点之间应有n-1个区段，研究中添加两个端点P₀和P_n+1，且满足P₀＝P₁；P_n+1＝P_n，联合前面的n-3段曲线就可得到n-1个曲线段，这样利用上述的抛物样条方程即可模拟树干曲线；

再次，利用已模拟的树干横断面曲线，量测树的胸径，根据树高-胸径关系式计算树高坐标，并从3个高度的树干横断面曲线中选取3个对应的型值点，总共形成抛物样条曲线的4个型值点，代入抛物样条曲线方程模拟树干纵断面曲线；

最后，将每个横断面曲线按角度平均分成1000份，按照上一步方法，形成1000条树干纵断面曲线，从而模拟出整个树干的立体形状；

第三步，树干材积量的计算：

首先，沿树干方向将树干平均分成9000段；

然后，计算每个分割段的下底面面积，由于树干立体曲线是由1000条树干纵断面曲线构成，因此每个分割段实际也被分割成1000份，计算出所有1000份的面积并累加，即可得到下底面面积；

最后，由于完整树干被平均分割成9000段，因此在计算体积时可以假设每一小段的上、下底面面积是近似相等的，体积计算就是其中一个底面积与高的乘积，通过步骤2计算出的底面积与分割段的高度可以得到每个分割段的体积，累加9000段体积即可得到整个树干的体积，即树干材积量。

2.如权利要求1所述的单木材积测量方法，其特征在于：所述的单木材积测量包括树干横断面曲线模拟、树干纵断面曲线模拟和树干体积测量。

3.如权利2所述的树干横断面曲线模拟方法，其特征在于：所述树干横断面曲线模拟，首先，从树干上钻取3根树芯，根据树芯数据计算型值点二维坐标，然后将二维坐标代入三次参数样条曲线方程模拟相邻型值点之间的曲线段，其次通过设定端点条件和连接条件将各独立曲线段进行光滑连接处理，形成树干横断面曲线。

4.如权利2所述的树干纵断面曲线模拟方法，其特征在于：所述树干纵断面曲线模拟，首先在树干0.5米、1.3米和2.5米3个高度模拟树干横断面曲线，然后从3条曲线上分别选取3个对应的型值点，并通过查找胸径-树高一元表得到树高信息，将型值点和树高信息代入抛物样条曲线方程模拟相邻型值点之间的曲线段，最后通过加权合成将独立的曲线段进行光滑连接处理，形成树干纵断面曲线。

5.如权利2所述的树干体积测量方法，其特征在于：所述树干体积测量，首先将树干平均分成9000份，然后假设每一份树干的上下底面积是近似相等，通过底面积与高度乘积得到每一份树干的体积，最后通过累加9000份树干体积得到整个树干体积。

Images