一种基于激光点云数据的树干材积获取方法与流程

本发明涉及一种树干参数提取技术，尤其是涉及一种基于激光点云数据的树干材积获取方法。

背景技术：

地面三维激光扫描技术(terrestriallaserscanningtechnology，简称tls)是一种新型的测绘技术，产生于20世纪90年代。tls可快速准确的获得目标物体表面的点云数据信息，从而可以快速构建目标物体的三维模型。经过近20多年来的发展，三维激光扫描仪已经可以连续快速地对被观测物体进行非接触式测量，其通过获取物体表面至扫描仪的距离和发射强度获取大量物体表面的三维点云数据。使用tls提取林业相关参数是近10年来的研究热点。

树干材积是测树学中的基本测树因子，是树木经济价值的直接体现，也是计算树木生物量与碳储量的基础。传统上将树干简化为一个绕中心轴旋转的旋转体，以此为基础估算树干材积。然而树干形状是不规则的，旋转体的简化结构过于简化了树干的几何形状，在此基础上计算得到的树干材积必定与树干的真实材积存在着差异，也正因如此，如何准确计算树干材积是林业工作者一直关注的问题。

技术实现要素：

有鉴于此，本发明的目的是针对现有技术的不足，提供一种基于激光点云数据的树干材积获取方法，用以根据树干点云获取精确的树干材积。

为达到上述目的，本发明采用以下技术方案：

一种基于激光点云数据的树干材积获取方法，包括如下步骤：

步骤1，对树干点云垂直分段；

步骤2，逐一计算各个垂直分段的体积；

步骤3，计算所有垂直分段体积的累加和。

优选的，所述步骤1，对树干点云垂直分段，包括：

以高度为标准，计算树干点云中的最低点plowest，其中zlowest为最低点plowest的高度值，并以所述最低点plowest为基础计算树干的高度；

构建平面π0与π1，其方程分别为z＝zlowest与z＝zlowest+thickness，其中thickness是每一个垂直分段的高度；则位于平面π0与π1间的点云形成树干点云的一个垂直分段q0；继续向上构建平面π2，其方程为z＝zlowest+thickness*2，则位于平面π1与π2间的树干点云形成另一个垂直分段q1；则，平面πi的方程可表示为z＝zlowest+thickness*i，位于平面πi与πi+1间的树干点云形成点云垂直分段qi；其中，平面πi与平面πi+1分别称为垂直分段qi的下平面与上平面；按照上述垂直分段的构建方法直至树梢部位得到若干个树干点云的垂直分段。

优选的，所述步骤2，逐一计算各个垂直分段的体积，包括：

计算垂直分段在下平面的投影点集；

计算垂直分段的断面轮廓点；

构建断面轮廓曲线并计算树干断面的断面积；

计算垂直分段的体积。

优选的，所述计算垂直分段在下平面的投影点集，包括：

计算垂直分段qi在其下平面的投影点集qi′；将垂直分段qi中的树干点云投影至垂直分段qi的下平面πi得到一个平面点集qi′；计算平面点集qi′的凸包多边形的质心点作为平面点集qi′的几何中心点ci。

优选的，所述计算垂直分段的断面轮廓点，包括：

设点pm为垂直分段qi中的一个点，pm′表示点pm在垂直分段qi下平面的投影，则几何中心点ci与投影点pm′构成的向量与x轴的夹角可表示为：根据角度分区的角度值θ，将平面点集qi′中的点分为个角度分区，每个角度分区qij中的点可表示为：计算每一个角度分区qij中点的重心点，所述重心点为所述角度分区中所有点的几何平均值；并将所述重心点作为所述角度分区的轮廓点，由此得到垂直分段qi的个轮廓点。

优选的，所述构建断面轮廓曲线并计算树干断面的断面积，包括：设垂直分段qi的断面轮廓点集合用rprofile表示，首先采用非有理三次b样条曲线构建一条插值于点集rprofile的闭合曲线，所述闭合曲线即为断面轮廓曲线；断面轮廓曲线所围区域的面积即为树干断面的断面积。

优选的，所述构建断面轮廓曲线并计算树干断面的断面积，包括：

将rprofile首部的3个点依次加入到rprofile的尾部同时也将尾部的3个点依次加入到rprofile的首部，构成插值点r＝{rk},k＝0,1,...,n；

使用弦长参数法计算插值点r的节点参数值令d为总弦长，节点参数值定义为：

根据节点参数值使用平均值法计算节点向量u＝{u0,u0,...,um}，有

其中，p＝3表示样条曲线的次数，n＝|r|表示插值点的个数，m＝n+p-1；

采用全局曲线插值的方法求解非有理三次b样条曲线的控制点，得到系数矩阵为(n+1)×(n+1)的线性方程组：

其中，ni,p(u)表示为

求解方程组得到n+1个控制点{pi}，根据非有理三次b样条曲线的表达式p＝3；可得到一条连续曲线，设置rprofile的第一个节点r0为首节点，从首节点处开始构建，当曲线环绕一周到达首节点时停止曲线的构建，即可得到完整的闭合曲线，所述的闭合曲线即为断面轮廓曲线；

断面轮廓曲线是一条闭合曲线，根据格林公式将该闭合曲线所围区域的面积转换为第一类曲线积分的计算，从而得到树干断面的断面积si。

优选的，所述计算垂直分段的体积，包括：

计算垂直分段qi的体积vi，将垂直分段近似为一个柱体，断面形状为断面轮廓曲线所围区域的形状，则根据柱体体积公式，垂直分段qi的体积vi可表示为vi＝si*thickness。

优选的，所述步骤3，计算所有垂直分段体积的累加和，包括：所有垂直分段体积的累加和即为树干材积，树干材积的计算公式可表示为v＝thickness*∑si，v为树干材积。

优选的，所述垂直分段的高度值thickness为0.5厘米。

本发明的有益效果是：

传统上使用区分求积法计算伐倒木的树干材积，其建立在将树干简化为树干横断面是圆，树干中心轴为直线的旋转体的基础上，该方法将树干区分成若干等长或不等长的区分段，使各区段的干形更接近于正几何体，再根据近似求积式计算各分段材积，合计各分段材积得到树干材积。计算立木的树干材积的传统方法是采用经验公式，将立木高度与直径的测量值代入材积经验公式得到立木的树干材积。由于材积经验公式建立在伐倒木材积计算的基础上，因此伐倒木与立木的材积计算都建立在树干横断面为圆，树干中心轴为直线的规则旋转体的几何理论基础之上。然而树干形状是不规则的，树干横断面也不是规则的圆形，因此在此简化基础上计算的材积与树干的真实材积存在着差异。

本发明建立在精确计算树干断面积的基础上，通过精确获取树干断面的轮廓点，进而构建一条闭合的断面轮廓曲线，闭合曲线所围面积即为断面面积，该段断面积计算方式建立在获取树干断面轮廓的基础上，贴近于树干断面的真实轮廓；而且使用的thickness值为0.5厘米，对于一段长度为0.5厘米的树干，其在各处的断面几乎相同，因此本发明计算的垂直分段的体积非常接近于垂直分段的真实体积。相对比传统上伐倒木和立木的不同材积计算方法，本发明在获取树干点云的基础上，无需区分伐倒木还是立木，均可使用本发明的方法计算树干材积。

本发明在精确计算树干断面断面积的基础上精确计算树干材积，相比于传统的将树干断面简化为圆来计算树干断面的断面积的方法，本发明在获取树干断面轮廓的基础上计算树干断面的断面积，具有较高的精确性，从而为准确提取树干生物量、碳储量与树干量材优化等与树干材积的相关研究与应用提供基础数据。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

图1是本发明一种基于激光点云数据的树干材积获取方法的流程图。

图2是本发明一种基于激光点云数据的树干材积获取方法步骤2的流程图。

图3是树干点云垂直分段的示例图。

图4是断面轮廓曲线的示例图。

具体实施方式

下面结合附图和实施例对本发明作进一步描述。

如图1所示，本发明提供一种基于激光点云数据的树干材积获取方法，包括如下步骤：

步骤1，对树干点云垂直分段；

步骤2，逐一计算各个垂直分段的体积；

步骤3，计算所有垂直分段体积的累加和。

步骤1，对树干点云垂直分段，包括：

根据z轴坐标，也即高度值，计算树干点云中的最低点plowest，设最低点plowest的z轴坐标为zlowest，并根据最低点为基础计算树干的高度；

构建平面π0与π1，其方程分别为z＝zlowest与z＝zlowest+thickness，其中thickness是每一个垂直分段的高度值。则位于π0与π1间的点云形成树干点云的一个垂直分段q0；继续向上构建平面π2，其方程为z＝zlowest+thickness*2，则位于平面π1与π2间的树干点云形成另一个垂直分段q1；则，平面πi的方程可表示为z＝zlowest+thickness*i，位于平面πi与πi+1间的树干点云形成点云垂直分段qi；其中，平面πi与平面πi+1分别称为垂直分段qi的下平面与上平面；按照上述垂直分段的构建方法直至树梢部位得到若干个树干点云的垂直分段；如图3所示，图中上平面1与上平面2之间是垂直分段qi。

thickness的值可根据实际需求预先设定，其值的大小可调，其值越小越好，但考虑到点云密度的影响，在一个实施例中，thickness的值为0.5厘米。thickness的值越小，计算得到的体积越精确，但受点云密度的影响，若thickness的值太小，则垂直分段中包括的点云太少，容易导致计算的断面轮廓点无法反映断面的真实轮廓，若thickness的值太大，树干在此段范围内干形有较大变化，则导致体积计算精度降低。

如图2所示，步骤2，逐一计算垂直分段qi的体积vi，包括如下步骤：

步骤201，计算垂直分段在下平面的投影点集；

步骤202，计算垂直分段的断面轮廓点；

步骤203，构建断面轮廓曲线并计算树干断面的断面积；

步骤204，计算垂直分段的体积。

对于一个特定的垂直分段qi，将qi中的树干点云投影至垂直分段的下平面得到平面投影点集qi′，计算平面点集qi′构成的凸包多边形的质心点ci，根据点ci和qi′中点构成的向量与角度分区的角度值θ度对投影点集qi′进行角度分区，再计算每一个角度分区的重心点，个角度分区的重心点构成断面的轮廓点rprofile；再采用非有理三次b样条曲线构建一条过点集rprofile的闭合曲线，该曲线为断面轮廓曲线；然后再使用格林公式计算闭合曲线所围区域的面积的积分转换为第一类曲线积分，计算得到树干断面的断面积si；垂直分段qi的体积vi可表示为vi＝thickness*si平方厘米。

在一个实施例中，θ值为1度(θ的值越小，断面轮廓点越多，断面轮廓会越准确，但受点云密度的影响，当θ太小，容易导致角度分区中无点云，当θ值为1度时，可确保断面在360个方位上都有轮廓点，基本上能真实反映树干断面的轮廓特征)。该实施例中，采用平面点集凸包多边形的质心点作为平面点集的几何中心点；根据角度分区的角度值与平面点集的几何中心点对平面点集进行角度分区的方法。

步骤2具体包括如下步骤：

计算垂直分段qi在其下平面的投影点集qi′，将垂直分段qi中的树干点云投影至垂直分段qi的下平面πi得到一个平面点集qi′；计算平面点集qi′的凸包多边形的质心点作为平面点集qi′的几何中心点ci；

设点pm为垂直分段qi中的一个点，pm′表示点pm在垂直分段qi下平面的投影，则几何中心点ci与投影点pm′构成的向量与x轴的夹角可表示为：

根据角度分区的角度值θ，将平面点集qi′中的点分为个角度分区，每个角度分区qij中的点可表示为：pm′∈qi′}，计算每一个角度分区qij中点的重心点(该角度分区中所有点的几何平均值)，并将该重心点作为该角度分区的轮廓点，由此得到垂直分段qi的个轮廓点；

构建断面轮廓曲线并计算树干断面的断面积，设垂直分段qi的断面轮廓点集合用rprofile表示，首先采用非有理三次b样条曲线构建一条插值于点集rprofile的闭合曲线，该闭合曲线即为断面轮廓曲线；断面轮廓曲线所围区域的面积即为树干断面的断面积。具体方法为：

将rprofile首部的3个点依次加入到rprofile的尾部同时也将尾部的3个点依次加入到rprofile的首部，构成插值点r＝{rk},k＝0,1,...,n；

使用弦长参数法计算插值点r的节点参数值令d为总弦长，参数值定义为：

根据节点参数值使用平均值法计算节点向量u＝{u0,u0,...,um}，有

其中，p＝3表示样条曲线的次数，n＝|r|表示插值点的个数，m＝n+p-1；

采用全局曲线插值的方法求解非有理三次b样条曲线的控制点，得到系数矩阵为(n+1)×(n+1)的线性方程组：

其中，ni,p(u)表示为

求解方程组得到n+1个控制点{pi}，根据非有理三次b样条曲线的表达式p＝3；可得到一条连续曲线，设置rprofile的第一个节点r0为首节点，从首节点处开始构建，当曲线环绕一周到达首节点时停止曲线的构建，即可得到完整的闭合曲线，所述的闭合曲线即为断面轮廓曲线；图4断面轮廓曲线的示例图；

断面轮廓曲线是一条闭合曲线，根据格林公式可将该闭合曲线所围区域的面积转换为第一类曲线积分的计算，从而得到树干断面的断面积si；

计算垂直分段qi的体积vi，将垂直分段近似为一个柱体，断面形状为断面轮廓曲线所围区域的形状，则根据柱体体积公式，垂直分段qi的体积vi可表示为vi＝si*thickness。

步骤3，计算所有垂直分段体积的累积和，包括：

合计步骤2计算的每一个垂直分段的体积即得到树干的材积，其公式为v＝thickness*∑si。在一个实施例中，thickness的值为0.5厘米，则公式可表示为v＝0.5*∑si。

需要说明的是，以上实施例仅用以说明本发明的技术方案而非限制，本领域普通技术人员对本发明的技术方案所做的其他修改或者等同替换，只要不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。