CN102682475A - 一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法 - Google Patents

Abstract

本发明涉及一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法。基于Dijkstra算法和最小生成树，提出了距离最小生成树方法(DMst)重构树木点云间的几何拓扑关系，获得近似树三维骨架，利用该骨架计算其上每点的主方向，采用物体表面点受力平衡的假设建立优化函数，以此修补缺失的树木点云信息，用DMst构建完整的树骨架结构，基于树异向生长特性模型计算树干的半径，实现树骨架的三维重建。本发明能很好的处理缺失点云数据，在大量点云缺失的情况下，仍能精确重建树的骨架结构。本发明具有很好的鲁棒性和数据自适应性。

Description

一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法

一、技术领域

本发明涉及一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法，属于空间信息技术领域。 

二、背景技术

植物生长的研究涉及到众多的科学研究领域，如植被遥感、生态等。地面激光雷达能够快速、精确的获取植被三维点云信息，从点云数据可以提取树木上的点集，利用该点集能够重建物体的表面结构。但是树木的几何结构较为复杂，尤其树木较多时，树木间相互遮挡，造成激光点云数据部分缺失。如何从地面激光点云甚至数据缺失的激光点云中重建树木的几何结构信息是空间信息技术领域亟待需要解决的问题之一。 

本发明的目的是在单一方向地面激光雷达扫描的情况下，采用物体表面点受力平衡的假设建立优化方程，让点云能够自动的填充被遮挡的部分，并且远离一些可能错误相交的部分，并且用距离最小生成树(DMst)构建的树型结构，既反映局部的点密度又兼顾点到根部的距离。 

目前基于点云对树三维骨架建模方面，常用的有基于启发式方法，该方法用单方向的激光点云重建树的主要树枝，在构建初始的树以后，用辐射传输模型对树被遮蔽的地方进行模拟，结合L系统来填充这些地方，但该方法建模效率较低，并且难以处理点云缺失的情况。也有学者采用八叉树的数据结构，通过设置假设，去除那些不符合假设的边，最后获得树形结构。但树结构复杂和相互遮挡，并且激光点云往往存有噪声，不能刚性的设置条件来获取树的结构，而是尽可能用可变形的方法来对树进行建模。 

本发明采用单方向的激光点云来构建树的三维结构，主要特色体现在：针对树的相互遮挡等原因造成点云数据缺失，甚至缺失比较严重的情况，基于物体表面连续且粒子在物体表面是受力平衡的思想，提出了优化方法，构建优化函数，可以较好的表达点云在树表面的特征，实现树木三维骨架的建模，而以往树建模方法几乎没有涉及过数据有较大缺失的情况。 

三、发明内容

1、目的：本发明的目的是从地面激光点云数据中自动重建三维树木骨架。由于树木几何形态比较复杂，待处理的点云数据往往不完整或存有噪声，本发明在单一方向的地面激光雷达扫描的情况下，采用物体表面点受力平衡的假设建立优化方程，自动填充缺失的区域，用距离最小生成树构建三维树骨架结构。本发明在植被遥感、生态等领域具有很好的应用价值。 

2、技术方案：本发明一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法，具体步骤如下(如图1)： 

步骤一：初步建立树木三维几何结构 

对于包含多棵树的地面激光点云，需要先识别树的株数及其对应的底部。获取一定高度层的点，然后聚类提取单棵树的点云集合。计算每棵树点云集合的中点，作为该树的底点。对于复杂的区域需要交互式指定底部来确定树的棵数和底点。 

在树木点云中，利用k-Nearest Neighbor(KNN)算法计算每个点最邻近的n个点，一般的，n＝50；设其中两点分别为i点和j点，KNN算法计算出的最邻近点i的最相邻点包含j，而j的最邻近不包括i，所以让i的最相邻点包含j，且点的相邻点包含i，点与点连线形成边的长度作为该边的权重，以此构建一个无向图，通过该无向图构建最小生成树(MST)。用Dijkstra算法计算每个底点到MST各点的最短距离，把每个点归属到离所对应底点最近的树上去。这样就获得了一株树的点云以及其连接方式。 

步骤二：构建距离最小生成树 

从树木点云中用Dijkstra算法获得树的三维结构缺乏相邻点的信息，而MST也不能很好描述树的伸展形状。树枝通常是向上生长的，小范围内树干的分布应该是相对光滑的，树上的点被认为是以离根节点最近的方式连接，并且越相临的点在同一树枝上的可能性越大。本发明提出距离最小生成树方法(DMst)来构建子树。DMst的方法流程如下： 

(1)对单株树的点云构建无向图，其中顶点集合为V。设每个顶点到根节点的距离D_dij，每个顶点到当前DMst(此时，DMst为空)的最小距离为D_mst，每个顶点指向为P，每个顶点的MST指向为P_mst，每个顶点Dijkstra树的指向为P_dij； 

(2)把底点做为DMst的根节点，初始化V，除底点以外的所有顶点放入V中，计算所有点的D_dij和D_mst，把根节点相邻点的P、P_dij和P_mst指向根节点。 

(3)在V中，选取D_mst+D_dij之和最小的点，把该点从V中移出。 

(4)更新P_mst和D_mst，设点i、j之间的距离为d_ij，如果D_mst[i]＞d_ij，那么P_mst指向j，D_mst[i]＝d_ij，更新P_dij和D_dij；如果D_dij[i]＞D_dij[j]+d_ij，那么P_mst指向j，D_dij[i]＝D_dij[j]+d_ij；当D_mst[i]+(D_mst[i]+D_dij[j])＜＝D_dij[i]-D_dij[k]+D_dij[i]时，P[i]指向距树最小的位置j，否则，P[i]指向距根节点距离最小的位置。 

(5)重复步骤(3)和(4)，直到V中没有点为止，通过P获取所有点的组织结构，完成DMst构建。 

构建DMst时，用D_mst[i]+w×D_dij[i]代替D_mst+D_dij控制入选的点。更新时，因为在DMst中每次获取i点和j点的公共前向，计算复杂度大，而用D_mst[i]+(D_mst[i]+D_dij[j])-D_dij[k]＜＝w×(D_dij[i]-D_dij[k]+D_dij[i]-D_dij[k])作为判定条件，可以获得较好效果。 

DMst获得的树形结构中，大量的树节点点只有少量的子节点，去除那些只有少量子节点的树节点，剩下的点构成一个近似的树骨架。设置每个点的权值c_i，对没有子节点的点来说，c_i＝0.1，对有子节点的点， 

$c_i=\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}(c_j+d_{\mathrm{ij}})---\left(1\right)$

式(1)中，Ω是i点相邻子节点的集合，j点是i点相邻的子节点，d_ij是i点和j点之间的距离。 

设阈值为δ，剔除c_i＜δ的点而保留c_i＞δ的点。保留的点构成了树的近似骨架，并且骨架点是树的表面点，骨架节点之间的方向很好反映了树枝的伸展方向。DMst构成树的连接方式具有方向性，因为考虑了每个点到根部的距离，形成树枝的方向主要是向上的，符合树的生长特征，这也为缺失点云区域的修补提供了帮助。 

步骤三：优化树骨架结构 

一般的，树的表面是连续的，树枝的变化是平滑的，不会出现粗枝瞬间变细枝或者细枝瞬间变成粗枝的情况，所以用如下假设来描述树的表面：地面激光雷达扫描仪扫描树木表面形成点云的过程可看作连续表面离散化的过程，点云类似于粒子，它们之间有相互的斥力，所在位置对其本身有一定的束缚作用，这些粒子所构成的表面处于受力平衡的状态，因此有点云缺失的树木区域，本应该存在的表面点没有了，整个表面的受力平衡状态就破坏了，粒子会产生运动，最后会获得一个新的受力平衡状态，而这个平衡状态就是修复了缺失的状态。(1)提取枝干主方向 

因为点密集分布在树的表面，相邻点之间的方向并不能代表树枝的伸展方向，所以点与点之间的方向不能代表表面上每一个点受力方向，用主方向来反映树枝的伸展方向。 

为了获得主方向，对去除c_i＜δ点的点云进行重采样。将c_i＞δ且所有后向点的c_i＜δ的点作为末端点，如果一个点有两个或两个以上后向点c_i＞δ，这个点是树的分叉点。在所有采样点构成的树骨架上，每隔距离d再采一个点，d是点与点之间树干上的距离而不是直线距离，把这些点按DMst的结构连接起来，d＝0.2。设相邻两点为i，j，离根节点近的点设为i，另一点设为j，由i指向j的方向作为这两点之间的方向，该方向作为i点的所有子节点的主方向；取末枝点的前向节点到末枝点的方向作为末枝点的所有子节点的主方向。 

(2)优化点云缺失的树木 

对于完整点云集合来说，其中每个点处于受力平衡状态，所以一个点所连接的所有点对该点的合力在主方向上的投影F_r与其所在位置对该点的束缚力F_s之和F＝F_r+F_s，作为该点所受的外力，显然F＝0时，该点处于平衡状态。类似于同性电荷之间的力，点之间的斥力与其距离成反比，通过公式(2)计算每个点的电荷量q_i， 

$q_i=\frac{1}{r_i}---\left(2\right)$

式(2)中，q_i是对应点i的电荷量，r_i是对应点i到n个最临近点的平均距离，相连i与点j之间的斥力是， 

$F_r(i,j)=\frac{q_i{q}_j}{{\left|\right|P_i-P_j\left|\right|}^2}(P_i-P_j)---\left(3\right)$

式(3)中P_i，P_j是优化后对应点i和j的坐标，q_i，q_j是对应点的电荷量。 

i点所在的原始位置对其的束缚力应该和i点偏离原始位置的距离成正比。通过公式(3)求当前点i的阻尼系数K_i。 

$K_i=\frac{\left(\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}{q}_j\right)q_i}{{r_m}^2}---\left(4\right)$

式(4)中，r_m是树上两点之间的平均距离，Ω是与i点相连接的所有点集合，j是i相连的点，i所在位置对i的束缚力为， 

F_s(i)＝K_i(Q_i-P_i)(5) 

式(5)中，P_i是优化后对应点i的坐标，Q_i是为优化前点i的坐标。最后的优化函数： 

$F=\min \underset{i\∈\mathrm{tree}}{\mathrm{\Σ}}{\left|\right|\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left(O_i^T{F}_r(i,j)O_i\right)+w{F}_s\left(i\right)\left|\right|}^2---\left(6\right)$

式(6)，Ω是i点所有相连接点的集合，j是i相连的点，O_i是i点的主方向，w是权重，w的改变可以控制点的偏离程度。 

对DMst进行重采样，并且对重采样的每一个末枝点进行延伸，寻找末枝点的后向子节点中c_i最大的点，作为待处理点，再搜素待处理点的后向子节点中c_i最大的点，如果点到末枝点在树干上的距离超过d，那么把这个点作为新的末枝点加到骨架中，直到找到一个没有后向点的点为止，并把该点加入到骨架。这样就得到最后的树骨架。 

步骤四：重建树木三维骨架 

用拉普拉斯方程(7)对上面获取的树骨架进行平滑。 

$L_{\mathrm{ij}}=\left\{\begin{array}{cc}{w}_{\mathrm{ij}}=1& \mathrm{if}(i,j)\∈E\\ {\mathrm{\Σ}}_{(i,k)\∈E}^k-w_{\mathrm{ik}}& \mathrm{if\; i}=j\\ 0& \mathrm{otherwise}\end{array}\right.---\left(7\right)$

L是n×n的拉普拉斯算子矩阵，E是i相邻点的集合，j和k是i的相邻点。 

$\left[\begin{array}{c}{W}_{L}L\\ W_H\end{array}\right]V^{\′}=\left[\begin{array}{c}0\\ W_{H}V\end{array}\right]---\left(8\right)$

公式(8)中，W_L和W_H是对角矩阵，用来平衡收缩力和吸引力；V’是优化后的点，V是初始点；W_L对角线上的值设为1，根据需要设置平滑的大小。 

因为点云数据包含噪声或存在缺失，难以确定树枝干的宽度。一般的，树底部附近的点比较密集，并且周围干扰点较少，可以估算出树底部的半径，根据树的异速生长模型计算所有树枝的半径。如果父节点只有一个子节点，那么， 

$r_c=r_p{\left(\frac{l_c}{l_p}\right)}^{\frac{3}{2}}---\left(9\right)$

r_c是子节点的半径，r_p是父节点的半径，l_c是子节点的c，l_p是父节点的c。 

如果一个父节点i有多个子节点，那么， 

$r_{\mathrm{ci}}=r_p{\left(\frac{l_{\mathrm{ci}}}{\underset{j}{\mathrm{\Σ}}{l}_{\mathrm{cj}}}\right)}^{\frac{3}{2}}---\left(10\right)$

j是子节点的个数。 

至此计算了每个节点的半径。对于有叶子的树，认为树干半径小但周围点云多的位置长叶子的可能较大，所以，在优化后点云中，树骨架上的每点的q_i除以该点处树干半径r_i，对这些值进行降序排序，根据树叶的疏密，添加叶子。 

3、优点及功效：本发明基于优化的思想构建目标函数，把实体表面点云看做粒子，认为粒子处于受力平衡状态，能够很好的修复由于树相互遮挡等原因造成的树木部分激光点云数据缺失，利用DMst重建树的三维骨架。本发明方法具有很好的鲁棒性、数据自适应性和较高的数据建模效率，能很好的处理缺失点云数据，在大量点云缺失的情况下，仍能精确重建树的骨架结构。 

四、附图说明

图1本发明方法的流程示意图 

图2(a)无叶子树点云(4,503个点) 

图2(b)树的三维模型 

图3(a)激光雷达扫描的方向的照片 

图3(b)长有密集树梢的树的点云 

图3(c)该树的三维模型 

图4(a)树枝的部分点云缺失 

图4(b)该树的三维模型 

图5(a)激光雷达扫描方向的照片 

图5(b)树相应的点云 

图5(c)构建的3颗交叠树的三维模型 

五、具体实施方式

本发明涉及一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法，该方法具体步骤如下(如图1)： 

步骤一：初步建立树木三维几何结构 

对于包含多棵树的地面激光点云，需要先识别树的株数及其对应的底部。获取一定高度层的点，然后聚类提取单棵树的点云集合。计算每棵树点云集合的中点，作为该树的底点。对于复杂的区域需要交互式指定底部来确定树的棵数和底点。 

在树木点云中，利用k-Nearest Neighbor(KNN)算法计算每个点最邻近的n个点，一般的， n＝50；设其中两点分别为i点和j点，KNN算法计算出的最邻近点i的最相邻点包含j，而j的最邻近不包括i，所以让i的最相邻点包含j，且点的相邻点包含i，点与点连线形成边的长度作为该边的权重，以此构建一个无向图，通过该无向图构建最小生成树(MST)。用Dijkstra算法计算每个底点到MST各点的最短距离，把每个点归属到离所对应底点最近的树上去。这样就获得了一株树的点云以及其连接方式。 

步骤二：构建距离最小生成树 

从树木点云中用Dijkstra算法获得树的三维结构缺乏相邻点的信息，而MST也不能很好描述树的伸展形状。树枝通常是向上生长的，小范围内树干的分布应该是相对光滑的，树上的点被认为是以离根节点最近的方式连接，并且越相临的点在同一树枝上的可能性越大。本发明提出距离最小生成树方法(DMst)来构建子树。DMst的方法流程如下： 

(6)对单株树的点云构建无向图，其中顶点集合为V。设每个顶点到根节点的距离D_dij，每个顶点到当前DMst(此时，DMst为空)的最小距离为D_mst，每个顶点指向为P，每个顶点的MST指向为P_mst，每个顶点Dijkstra树的指向为P_dij； 

(7)把底点做为DMst的根节点，初始化V，除底点以外的所有顶点放入V中，计算所有点的D_dij和D_mst，把根节点相邻点的P、P_dij和P_mst指向根节点。 

(8)在V中，选取D_mst+D_dij之和最小的点，把该点从V中移出。 

(9)更新P_mst和D_mst，设点i、j之间的距离为d_ij，如果D_mst[i]＞d_ij，那么P_mst指向j，D_mst[i]＝d_ij，更新P_dij和D_dij；如果D_dij[i]＞D_dij[j]+d_ij，那么P_mst指向j，D_dij[i]＝D_dij[j]+d_ij；当D_mst[i]+(D_mst[i]+D_dij[j])＜＝D_dij[i]-D_dij[k]+D_dij[i]时，P[i]指向距树最小的位置j，否则，P[i]指向距根节点距离最小的位置。 

(10)重复步骤(3)和(4)，直到V中没有点为止，通过P获取所有点的组织结构，完成DMst构建。 

构建DMst时，用D_mst[i]+w×D_dij[i]代替D_mst+D_dij控制入选的点。更新时，因为在DMst中每次获取i点和j点的公共前向，计算复杂度大，而用D_mst[i]+(D_mst[i]+D_dij[j])-D_dij[k]＜＝w×(D_dij[i]-D_dij[k]+D_dij[i]-D_dij[k])作为判定条件，可以获得较好效果。 

DMst获得的树形结构中，大量的树节点点只有少量的子节点，去除那些只有少量子节点的树节点，剩下的点构成一个近似的树骨架。设置每个点的权值c_i，对没有子节点的点来说，c_i＝0.1，对有子节点的点， 

$c_i=\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}(c_j+d_{\mathrm{ij}})---\left(1\right)$

式(1)中，Ω是i点相邻子节点的集合，j点是i点相邻的子节点，d_ij是i点和j点之间的距离。 

设阈值为δ，剔除c_i＜δ的点而保留c_i＞δ的点。保留的点构成了树的近似骨架，并且骨架点是树的表面点，骨架节点之间的方向很好反映了树枝的伸展方向。DMst构成树的连接方式 具有方向性，因为考虑了每个点到根部的距离，形成树枝的方向主要是向上的，符合树的生长特征，这也为缺失点云区域的修补提供了帮助。 

步骤三：优化树骨架结构 

一般的，树的表面是连续的，树枝的变化是平滑的，不会出现粗枝瞬间变细枝或者细枝瞬间变成粗枝的情况，所以用如下假设来描述树的表面：地面激光雷达扫描仪扫描树木表面形成点云的过程可看作连续表面离散化的过程，点云类似于粒子，它们之间有相互的斥力，所在位置对其本身有一定的束缚作用，这些粒子所构成的表面处于受力平衡的状态，因此有点云缺失的树木区域，本应该存在的表面点没有了，整个表面的受力平衡状态就破坏了，粒子会产生运动，最后会获得一个新的受力平衡状态，而这个平衡状态就是修复了缺失的状态。 

(1)提取枝干主方向 

因为点密集分布在树的表面，相邻点之间的方向并不能代表树枝的伸展方向，所以点与点之间的方向不能代表表面上每一个点受力方向，用主方向来反映树枝的伸展方向。 

为了获得主方向，对去除c_i＜δ点的点云进行重采样。将c_i＞δ且所有后向点的c_i＜δ的点作为末端点，如果一个点有两个或两个以上后向点c_i＞δ，这个点是树的分叉点。在所有采样点构成的树骨架上，每隔距离d再采一个点，d是点与点之间树干上的距离而不是直线距离，把这些点按DMst的结构连接起来，d＝0.2。设相邻两点为i，j，离根节点近的点设为i，另一点设为j，由i指向j的方向作为这两点之间的方向，该方向作为i点的所有子节点的主方向；取末枝点的前向节点到末枝点的方向作为末枝点的所有子节点的主方向。 

(2)优化点云缺失的树木 

对于完整点云集合来说，其中每个点处于受力平衡状态，所以一个点所连接的所有点对该点的合力在主方向上的投影F_r与其所在位置对该点的束缚力F_s之和F＝F_r+F_s，作为该点所受的外力，显然F＝0时，该点处于平衡状态。类似于同性电荷之间的力，点之间的斥力与其距离成反比，通过公式(2)计算每个点的电荷量q_i， 

$q_i=\frac{1}{r_i}---\left(2\right)$

式(2)中，q_i是对应点i的电荷量，r_i是对应点i到n个最临近点的平均距离，相连i与点j之间的斥力是， 

$F_r(i,j)=\frac{q_i{q}_j}{{\left|\right|P_i-P_j\left|\right|}^2}(P_i-P_j)---\left(3\right)$

式(3)中P_i，P_j是优化后对应点i和j的坐标，q_i，q_j是对应点的电荷量。 

i点所在的原始位置对其的束缚力应该和i点偏离原始位置的距离成正比。通过公式(3)求当前点i的阻尼系数K_i。 

$K_i=\frac{\left(\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}{q}_j\right)q_i}{{r_m}^2}---\left(4\right)$

式(4)中，r_m是树上两点之间的平均距离，Ω是与i点相连接的所有点集合，j是i相连的点，i所在位置对i的束缚力为， 

F_s(i)＝K_i(Q_i-P_i)(5) 

式(5)中，P_i是优化后对应点i的坐标，Q_i是为优化前点i的坐标。最后的优化函数： 

$F=\min \underset{i\∈\mathrm{tree}}{\mathrm{\Σ}}{\left|\right|\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left(O_i^T{F}_r(i,j)O_i\right)+w{F}_s\left(i\right)\left|\right|}^2---\left(6\right)$

式(6)，Ω是i点所有相连接点的集合，j是i相连的点，O_i是i点的主方向，w是权重，w的改变可以控制点的偏离程度。 

对DMst进行重采样，并且对重采样的每一个末枝点进行延伸，寻找末枝点的后向子节点中c_i最大的点，作为待处理点，再搜素待处理点的后向子节点中c_i最大的点，如果点到末枝点在树干上的距离超过d，那么把这个点作为新的末枝点加到骨架中，直到找到一个没有后向点的点为止，并把该点加入到骨架。这样就得到最后的树骨架。 

步骤四：重建树木三维骨架 

用拉普拉斯方程(7)对上面获取的树骨架进行平滑。 

$L_{\mathrm{ij}}=\left\{\begin{array}{cc}{w}_{\mathrm{ij}}=1& \mathrm{if}(i,j)\∈E\\ {\mathrm{\Σ}}_{(i,k)\∈E}^k-w_{\mathrm{ik}}& \mathrm{if\; i}=j\\ 0& \mathrm{otherwise}\end{array}\right.---\left(7\right)$

L是n×n的拉普拉斯算子矩阵，E是i相邻点的集合，j和k是i的相邻点。 

$\left[\begin{array}{c}{W}_{L}L\\ W_H\end{array}\right]V^{\′}=\left[\begin{array}{c}0\\ W_{H}V\end{array}\right]---\left(8\right)$

公式(8)中，W_L和W_H是对角矩阵，用来平衡收缩力和吸引力；V’是优化后的点，V是初始点W_L对角线上的值设为1，根据需要设置平滑的大小。 

因为点云数据包含噪声或存在缺失，难以确定树枝干的宽度。一般的，树底部附近的点比较密集，并且周围干扰点较少，可以估算出树底部的半径，根据树的异速生长模型计算所有树枝的半径。如果父节点只有一个子节点，那么， 

$r_c=r_p{\left(\frac{l_c}{l_p}\right)}^{\frac{3}{2}}---\left(9\right)$

r_c是子节点的半径，r_p是父节点的半径，l_c是子节点的c，l_p是父节点的c。 

如果一个父节点i有多个子节点，那么， 

$r_{\mathrm{ci}}=r_p{\left(\frac{l_{\mathrm{ci}}}{\underset{j}{\mathrm{\Σ}}{l}_{\mathrm{cj}}}\right)}^{\frac{3}{2}}---\left(10\right)$

j是子节点的个数。 

至此计算了每个节点的半径。对于有叶子的树，认为树干半径小但周围点云多的位置长叶子的可能较大，所以，在优化后点云中，树骨架上的每点的q_i除以该点处树干半径r_i，对这些值进行降序排序，根据树叶的疏密，添加叶子。 

实施例1： 

在一台配置有Intel(R)Core^TM 2.4GHz 2处理器，2G内存，ATI Radeon HD图形显卡的计算机上进行了实施。采用lmd360型和z62型获取的地面激光雷达点云数据，并且z62型地面激光雷达获取点时，同时拍摄了对应的照片。 

对包含一棵基本完整树的z62型地面激光点云数据(图2a)进行建模(图2b)，可以看出树主干中端略有断开，采用上面方法构建的树形状得到较好保留，并且很好描述了树之间的交叉。图2和图3所示的点云都只有少量点云被其他树或者自身遮挡，可以较好处理。即使由于遮挡等原因造成有大范围缺失的树点云数据，本发明仍然可以从保留的部分枝干信息中有效的对缺失点云进行修复，从而得到精确的树的三维骨架。图3是z62型地面激光雷达扫描的结果，包含6,811个点。图4(a)中由于扫描过程中前面树的遮挡造成点云有两处明显的数据缺失，图4(b)是本发明方法对该点云进行建模的结果。树的建模结果较好贴近了树真实的三维结构，缺失的点云有较好的修复。图5(b)是lmd360型地面激光雷达扫描的三棵交叠长有叶子树的点云，包含220,915个点。图5(c)是本发明建模的结果。 

表1列出了树点云的个数、树的数目以及建模时间。本发明方法对树的骨架结构有较好的描述，并且具有很高的建模效率。 

表1.相关信息统计 

点的个数	树的棵数	时间
			10663	1	9秒
4,503	1	15秒
			11,856	1	21秒
111,833	1	2分2秒
			6,811	1	8秒
220,915	3	5分22秒
			1,173,621	22	12分

Claims (1)

1.一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法，步骤包括：

步骤一：初步建立树木三维几何结构

对于包含多棵树的地面激光点云，需要先识别树的株数及其对应的底部，获取一定高度层的点，然后聚类提取单棵树的点云集合，计算每棵树点云集合的中点，作为该树的底点，对于复杂的区域需要交互式指定底部来确定树的棵数和底点。

在树木点云中，利用k-Nearest Neighbor(KNN)算法计算每个点最邻近的n个点，以此构建一个无向图，通过该无向图构建最小生成树(MST)，用Dijkstra算法计算每个底点到MST各点的最短距离，把每个点归属到离所对应底点最近的树上去。

步骤二：构建距离最小生成树

从树木点云中用Dijkstra算法获得树的三维结构缺乏相邻点的信息，而MST也不能很好描述树的伸展形状；树枝通常是向上生长的，小范围内树干的分布应该是相对光滑的，树上的点被认为是以离根节点最近的方式连接，并且越相临的点在同一树枝上的可能性越大，本发明提出距离最小生成树方法(DMst)来构建子树，DMst的方法流程如下：

(1)对单株树的点云构建无向图，其中顶点集合为V，设每个顶点到根节点的距离D_dij，每个顶点到当前DMst的最小距离为D_mst，每个顶点指向为P，每个顶点的MST指向为P_mst，每个顶点Dijkstra树的指向为P_dij；

(2)把底点做为DMst的根节点，初始化V，除底点以外的所有顶点放入V中，计算所有点的D_dij和D_mst，把根节点相邻点的P、P_dij和P_mst指向根节点。

(3)在V中，选取D_mst+D_dij之和最小的点，把这个点从V中移出。

(4)更新P_mst和D_mst，设点i、j之间的距离为d_ij，如果D_mst[i]＞d_ij，那么P_mst指向j，D_mst[i]＝d_ij，更新P_dij和D_dij；如果D_dij[i]＞D_dij[j]+d_ij，那么P_mst指向j，D_dij[i]＝D_dij[j]+d_ij；当D_mst[i]+(D_mst[i]+D_dij[j])＜＝D_dij[i]-D_dij[k]+D_dij[i]时，P[i]指向距树最小的位置j，否则，P[i]指向距根节点距离最小的位置。

(5)重复步骤(3)和(4)，直到V中没有点为止，通过P获取所有点的组织结构，完成DMst构建。

构建DMst时，用D_mst[i]+w×D_dij[i]代替D_mst+D_dij控制入选的点，更新时，用D_mst[i]+(D_mst[i]+D_dij[j])-D_dij[k]＜＝w×(D_dij[i]-D_dij[k]+D_dij[i]-D_dij[k])作为判定条件。

DMst获得的树形结构中，去除那些只有少量子节点的树节点，剩下的点构成一个近似的树骨架，设置每个点的权值c_i，对没有子节点的点来说，c_i＝0.1，对有子节点的点，

$c_i=\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}(c_j+d_{\mathrm{ij}})---\left(1\right)$

式(1)中，Ω是i点相邻子节点的集合，j点是i点相邻的子节点，d_ij是i点和j点之间的距离。

设阈值为δ，剔除c_i＜δ的点而保留c_i＞δ的点，保留的点构成了树的近似骨架。

步骤三：优化树骨架结构

(1)提取枝干主方向

为了获得主方向，对去除c_i＜δ点的点云进行重采样，将c_i＞δ且所有后向点的c_i＜δ的点作为末端点，如果一个点有两个或两个以上后向点c_i＞δ，这个点是树的分叉点，在所有采样点构成的树骨架上，每隔距离d再采一个点，d是点与点之间树干上的距离而不是直线距离，把这些点按DMst的结构连接起来，设相邻两点为i、j，离根节点近的点设为i，另一点设为j，由i指向j的方向作为这两点之间的方向，该方向作为i点的所有子节点的主方向；取末枝点的前向节点到末枝点的方向作为末枝点的所有子节点的主方向。

(2)优化点云缺失的树木

对于完整点云集合来说，一个点所连接的所有点对该点的合力在主方向上的投影F_r与其所在位置对该点的束缚力F_s之和F＝F_r+F_s，作为该点所受的外力，显然F＝0时，该点处于平衡状态。

类似于同性电荷之间的力，点之间的斥力与其距离成反比，通过公式(2)计算每个点的电荷量q_i，

$q_i=\frac{1}{r_i}---\left(2\right)$

式(2)中，q_j是对应点i的电荷量，r_i是对应点i到n个最临近点的平均距离，相连点i与点j之间的斥力是，

$F_r(i,j)=\frac{q_i{q}_j}{{\left|\right|P_i-P_j\left|\right|}^2}(P_i-P_j)---\left(3\right)$

式(3)中P_i，P_j是优化后对应点i和j的坐标，q_i，q_j是对应点的电荷量。

i点所在的原始位置对其的束缚力应该和i点偏离原始位置的距离成正比，通过公式(3)求当前点i的阻尼系数K_i，

$K_i=\frac{\left(\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}{q}_j\right)q_i}{{r_m}^2}---\left(4\right)$

式(4)中，r_m是树上两点之间的平均距离，Ω是与i点相连接的所有点集合，j是i相连的点，i所在位置对i的束缚力为，

F_s(i)＝K_i(Q_i-P_i)(5)

式(5)中，P_i是优化后对应点i的坐标，Q_i是为优化前点i的坐标，最后的优化函数：

$F=\min \underset{i\∈\mathrm{tree}}{\mathrm{\Σ}}{\left|\right|\underset{j\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left(O_i^T{F}_r(i,j)O_i\right)+w{F}_s\left(i\right)\left|\right|}^2---\left(6\right)$

式(6)中，Ω是i点所有相连接点的集合，j是i相连的点，Oi是i点的主方向，w是权重。

对DMst进行重采样，对重采样的每一个末枝点进行延伸，寻找末枝点的后向子节点中c_i最大的点，作为待处理点，再搜素待处理点的后向子节点中c_i最大的点，如果点到末枝点在树干上的距离超过d，那么这个点作为新的末枝点加到骨架中，直到找到一个没有后向点的点为止，并把该点加入到骨架，这样就得到最后的树骨架。

步骤四：重建树木三维骨架

用拉普拉斯方程(7)对上面获取的树骨架进行平滑。

$L_{\mathrm{ij}}=\left\{\begin{array}{cc}{w}_{\mathrm{ij}}=1& \mathrm{if}(i,j)\∈E\\ {\mathrm{\Σ}}_{(i,k)\∈E}^k-w_{\mathrm{ik}}& \mathrm{if\; i}=j\\ 0& \mathrm{otherwise}\end{array}\right.---\left(7\right)$

L是n×n的拉普拉斯算子矩阵，E是i相邻点的集合，j和k是i的相邻点。

$\left[\begin{array}{c}{W}_{L}L\\ W_H\end{array}\right]V^{\′}=\left[\begin{array}{c}0\\ W_{H}V\end{array}\right]---\left(8\right)$

公式(8)中，W_L和W_H是对角矩阵，用来平衡收缩力和吸引力；V’是优化后的点，V是初始点；W_L对角线上的值设为1，根据需要设置平滑的大小。

根据树的异速生长模型计算所有树枝的半径，如果父节点只有一个子节点，那么，

$r_c=r_p{\left(\frac{l_c}{l_p}\right)}^{\frac{3}{2}}---\left(9\right)$

r_c是子节点的半径，r_p是父节点的半径，l_c是子节点的c，l_p是父节点的c。

如果一个父节点i有多个子节点，那么，

$r_{\mathrm{ci}}=r_p{\left(\frac{l_{\mathrm{ci}}}{\underset{j}{\mathrm{\Σ}}{l}_{\mathrm{cj}}}\right)}^{\frac{3}{2}}---\left(10\right)$

j是子节点的个数。

至此计算了每个节点的半径，对于有叶子的树，在优化后点云中，树骨架上的每点的q_i除以该点处树干半径r_i，对这些值进行降序排序，根据树叶的疏密，添加叶子。

Images