基于类主干点的三维树木几何模型重建方法与流程
本发明涉及的是一种基于类主干点的三维树木几何模型重建方法，属于城市地面植被测量及建模技术领域。
背景技术：
在数字化城市场景的过程中，由于树木、灌木等植被在城市场景中随处可见，快速、高效地对树木三维模型进行重建也是城市建模不可或缺的重要组成部分，高精度且逼真的树木模型能够使得数字化的城市场景更富有真实感与沉浸感。此外，随着当前大数据与信息化进程的发展，数字林业建设过程也进一步加快，构建更加真实的树木模型也是推动林业数字化进程的基础和关键。
激光雷达技术（lightdetectionandranging，lidar）以作业场景大、空间数据获取速度快、自动化程度高、精度高、时效性好、通用性强等特点，逐渐成为城市场景和林业进行树木三维空间数据采集的重要方式，同时也为进行高精度、大范围场景林木模型的重建工作研究提供了数据保障。目前通过地面、机载、移动式激光扫描获取的激光点云数据，均可用于树木模型的重建工作。其中地面激光雷达技术（terrestriallaserscanning，tls）获取的点云数据精度已能达到毫米级，且得到的点云数据密度大（单株树木经单站扫描已能获取近百万数量的点）、精度高，单棵树木的完整性较好，因此在不结合树木照片、精确点位等异源数据信息的情况下，单独利用地面激光扫描技术，已能实现较高精度的树木三维模型重建；机载激光雷达（airbornelaserscanning，als）技术是利用航空飞机或无人机平台，实施大范围林木场景、快速高效模型重建，广泛应用于树冠、林木顶端结构的模型构建。当前的机载激光点云精度与密度都有了极大提高，2017年纽约大学发布的dublinals点云数据集已经达到300点/平方米，为实现大范围林木场景的模型重建提供了良好的数据源；移动式激光扫描（mobilelaserscanning，mls）技术是一种新兴的、综合的、髙效灵活的空间信息获取手段，在实施高效获取树木点云数据的同时，也存储了自身运动轨迹信息与相关全景影像最具有代表性特点。此外，移动式激光扫描通常搭载了由全球卫星导航系统(gnss-globalnavigationsatellitesystem)和惯性导航系统(ins-inertialnavigationsystem)组合而成的定位定姿系统(pos)，为获取的点云数据提供了精确的位置信息。然而，无论是地面、机载还是移动式激光扫描技术，其获取的点云数据，往往存在着点密度不均匀、噪声点与离值点以及树冠内部遮挡严重、数据缺失等问题；此外树木本体具有形态各异、无规则可循、几何结构复杂，这种点云自身缺陷与树木几何形态固有的特点为树木模型的重建，带来了巨大的挑战。近二十年来，模式识别、人工智能与虚拟现实技术的全面发展，在促进激光雷达与点云数据建模研究发展的同时，也对二者提出了更高的要求。如何自动化、高效化、高保真度、轻量化地实现树木模型重建，是数字城市与数字林业建设过程中必须重视与解决的问题。
当前主流的植被地面点云建模方法主要包括：基于聚类思想的建模、基于图论方法的建模、基于先验假设的约束建模、基于拉普拉斯算子的建模、以及基于轻量化的树木建模；上述建模方法均存在不同程度的缺陷：
（1）基于聚类思想的建模，理论上可以提取任意树木的骨架点，但聚类算法的时间复杂度较高，并随着点数增加而呈现指数变化，不适合大数据量点云的建模研究；此外，聚类算法对点邻域的界定存在很大的依赖，点的邻域图构建直接影响后续聚类过程的进行。
（2）基于图论方法的建模，往往对体素的分辨率有很大的依赖性，体素的规格大小会对点云的计算精度与邻域范围分析产生影响：体素过大，骨架节点提取精度降低，体素过小又会加大计算成本、削弱计算效率。因此，合理的体素大小选择也是利用该方法进行建模时所必须解决的问题；最小生成树极其改进方法，在抗噪声与数据缺失上比较鲁棒，能够很好地对缺失点云数据进行弥补，但优化迭代过程需要进行进一步研究。
（3）基于先验假设的约束建模方法依赖去噪的预处理过程，且由于现实场景中，由于自然风、光照等环境因素的诱导使得树木结构并非完全规则的圆柱体或截面并非规则椭圆体，严格的先验假设约束在一定程度上损失了树木模型的精度。
（4）基于拉普拉斯算子的建模策略中算子的平滑过程时间复杂度高，建模效率受到影响，利用拉普拉斯算子建模时，往往未兼顾算子的时间复杂度与骨架线提取精度。
（5）基于轻量化的树木建模策略中点云三维空间信息的完整性与全面性被削弱，不考虑骨架特别是树冠内部树枝骨架，必然会造成枝条拓扑连接的不准确与细枝模型的丧失，在充分表达树木特别是植冠内部一些细节处的信息上，也存在部分拓扑链接的缺陷。
技术实现要素：
本发明的目的在于克服上述现有树木三维建模方法存在的种种缺陷，提出一种基于类主干点的三维树木几何模型重建方法，该方法对点云数据缺失和点云的密度不敏感，对试验区不同的树种建模具有较高的鲁棒性，通过对不同密度的点云和不同抽稀层次下的点云进行效率测试，结合类主干点表达和点云抽稀共同确保对大场景树木点云数据的重建。
本发明的技术解决方案：基于类主干点的三维树木几何模型重建方法，具体包括如下步骤：
（1）单株树木分割；
（2）构建树木初始骨架；
（3）优化初始骨架；
（4）确定枝干半径；
（5）添加树叶.
所述的步骤（1）单株树木分割，具体包括如下步骤：
对实验场景进行多站点扫描，扫描时合理布设标靶，利用标靶信息配准多站扫描数据，避免多株树木之间的遮挡和单株树木枝干之间自遮挡造成的数据缺失，从而保证扫描树木点云的完整性;在此基础之上，进一步利用cloudcompare软件，从扫描场景点云数据中手动分割单株树木，以此作为树木几何建模算法的输入。
所述的步骤（2）构建树木初始骨架依据激光反射强度值的差异，从预处理原始树木点云中，提取出树木的“类主干点”，然后采用图论中的“最小生成树”算法组织“类主干点”和剩余的“非主干点”，形成树木的初始骨架，最后通过“点密度调整”、“树枝平滑”对初始树木骨架进一步优化，得到比较精确且分布合理的树木骨架；具体包括如下步骤：
①提取“类主干点”
从原始树木点云中粗略区分出“类主干点”和“非主干点”，其中，“类主干点”主要指位于树主干和树冠中优势枝干上的点，而“非主干点”主要指位于树冠细枝和树叶上的点。首先从预处理后的原始树木点云中，依据激光点云反射强度差异，采用“直方图统计分析”方法，通过实验得到经验阈值
②构建树木的粗略骨架
采用无向图g(v,e)组织点云，并采用数组表示法存储点云，其中数组vertex存储数据元素信息（即顶点），数组edge存储数据元素之间的关系（即边）；提取出的“类主干点”中的每一个点作为图中的一个顶点，搜索每一个顶点周围最邻近的k个点，连接该点和这k个点形成边，并将每条边两顶点间的欧氏距离度量作为这条边的权重；采用prim算法，从上述无向图中构建mst，从而提取树木的初始骨架。
所述的步骤（3）优化初始骨架，具体包括如下步骤：
①点密度调整
基于ε邻近图进行计算，ε邻近图的定义为：两点i~j之间具有相互连接的关系，如果j∈b(i;ε)，其中b(i;ε)表示以i为中心以ε为欧式度量距离的原始点云的集合。我们基于ε邻近图组织原始树木点云，并根据公式（1）计算点云的密度：
di=n(b(i;ε))/ε（1）
其中n(b(i;ε))表示ε邻近图中与i点直接连接的点的数目，ε取值为0.1m；生成mst时，采用欧式距离作为边的权重；
树点云的分布特征：（i）在不存在数据缺失的情况下，mst上两点间的点云密度是相似的，在较小的范围内，点密度具有连续性且过渡平滑；（ii）树枝区域点密度大，说明该位置存在粗壮的树干的可能性较大，反之，该位置存在一些细枝或数据缺失的可能性较大；基于上述点云分布特征，修改无向图g(v,e)中每条边的权重，以适应点密度的变化：对应性质（i），相邻节点往往拥有相似的点密度，所以具有相似点密度节点间的权重小于非相似节点间的权重；对应性质（ii），优先连接点密度大的区域点云，继而判断点密度区域较小的点云是否与之连接，也就是说要优先在具有较大点密度或相似点密度的区域构建树枝骨架；此时，图g(v,e)中边的权重由公式（2）决定：
eij=dij/(di×dj)×|di-dj|（2）
其中参数di和dj分别是与点i和j相连的所有点组成的点集密度，dij是点i与j的欧式距离。
②树枝平滑
采用式（3）laplacian方程对点云密度调整后的骨架进一步优化。
参数wl和wh是对角矩阵，wl控制平滑的程度，wh控制平滑后的点云与原始点云相似的程度。wl中第i个对角元素的值定义成wl,i，wh中第i个对角元素的的值定义成wh,i。如果wl中的元素较大，那么平滑的程度较大；如果wh中的元素较大，平滑后的点云和原始点云相似度越高，平滑的程度越小。v是原始的点集，v'是优化过后的点集，l是n×n的拉普拉斯算子，通过公式（4）定义：
式中，e表示骨架结构上的边的集合，i、j、k表示骨架结构上的点。将对角矩阵wl的对角线上的元素值设为1，使每个点具有相同的收缩能力，,通过设置wh控制最后平滑的效果，矩阵wh对角线上的值是通过下述定义给出：
式中，k是i的父节点，j是i的子节点,θkij是∠kij。
所述的步骤（4）确定枝干半径，具体包括如下步骤：
分别采用“圆柱体拟合”和“异速生长模型”计算构架结构上某一特点pi的半径，结果分别为ric和ria，然后检查ric是否满足如下公式：
如果满足，则令ric为当前点pi的真值半径，该过程计算完成后，本发明优化方程（7）计算骨架上其它未知点的半径，该式用已知点的真值半径作为正则项，采用“异速生长模型”构建约束项。
所述的步骤（5）添加树叶，具体包括如下步骤：
获得每一节点的半径后，骨架模型膨胀成三维树模型，构建成一棵没有叶子的树木三维几何模型；针对每一个节点，根据公式（7）判断该节点是否存在叶子。
gi=di/ri（7）
参数ri代表每个节点的半径，di是与i点相连的所有点的点云密度；所有节点获得了对应的gi值后，对这些值从大到小排序，设定一个百分比的阈值控制整棵树存在叶子节点的比例，获得这些点的位置以后，在这些点上随机加上一些叶片模型，则完成完整的三维树木几何模型重建。
本发明的优点：（1）基于反射强度直方图分析的初始树木骨架提取算法：先提取“类主干点”再进行骨架优化的方法能提高骨架生成的效率，对较小的数据缺失和噪声不敏感，既能灵活地适应树枝的伸展方向，又能良好地描述树枝的细节，在散乱的树冠点云中，避免了树的主干部分的连接方式受细小的树枝影响。基于初始骨架的提取算法，相当于使得输入的植被点云被抽稀为一定比例的关键点，保证算法可以应用于处理大场景植被点云数据。
（2）结合“点密度调整”和“树枝平滑”的树木骨架的优化算法：由于面激光扫描数据一般存在大量的缺失和点密度极其不均质，利用此类树木点云数据进行建模是一项具有挑战性的工作，通过初始骨架的优化，恰恰可以处理点云缺失数据，在一定程度上可以缓解由于单站扫描对树木建模带来的挑战。
（3）本发明充分耦合“圆柱拟合”和“异速生长模型”这两类算法，计算树木枝干的半径。在树木枝干的精确性和不确定性两方面取得平衡，从而实现对枝干半径的合理估算。
附图说明
附图1是本发明基于类主干点的三维树木几何模型重建方法整体流程图。
具体实施方式
下面结合附图对本发明的技术方案作进一步说明。
如图1所示的基于类主干点的三维树木几何模型重建方法，首先从预处理完毕的树木点云中依据激光反射强度值的差异提取“类主干点”，采用最小生成树算法组织“类主干点”和剩余点，形成树的初步骨架并对骨架作进一步优化；然后在优化过的骨架上计算树枝半径并添加树叶；最后输出绘制好的三维树木几何模型。
（1）单株树木分割
本发明首先对实验场景进行多站点扫描，扫描时合理布设标靶，利用标靶信息配准多站扫描数据，有效避免多株树木之间的遮挡和单株树木枝干之间自遮挡造成的数据缺失，从而保证扫描树木点云的完整性。在此基础之上，进一步利用cloudcompare商业软件，从扫描场景点云数据中手动分割单株树木，以此作为树木几何建模算法的输入。
（2）构建树木初始骨架
树木骨架是构建三维树木模型的基础，本发明着重阐述如何从地面激光点云中提取树木的初始骨架，进而优化初始骨架，最终实现树木骨架的提取。为此，首先依据激光反射强度值的差异，从预处理原始树木点云中，提取出树木的“类主干点”，然后采用图论中的“最小生成树”算法组织“类主干点”和剩余的“非主干点”，形成树木的初始骨架，最后通过“点密度调整”、“树枝平滑”对初始树木骨架进一步优化，得到比较精确且分布合理的树木骨架。
由于树干、树冠细枝和树叶等部分的几何、材质显著不同，因此激光脉冲从树木不同部位得到的后向散射强度值（光谱）存在显著差异。据此，本发明设计一种提取树木“类主干点”的方法，从原始树木点云中粗略区分出“类主干点”和“非主干点”。其中，“类主干点”主要指位于树主干和树冠中优势枝干上的点，而“非主干点”主要指位于树冠细枝和树叶上的点。“非主干点”往往对树木骨架提取构建起到干扰作用，因此为充分区分这两类点集，本发明首先从预处理后的原始树木点云中，依据激光点云反射强度差异，采用“直方图统计分析”方法，通过实验得到经验阈值
获得树木的“类主干点”以后，本发明利用图论中的最小生成树算法mst（minimumspanningtree，mst）组织“类主干点”和部分“非主干点”，形成树木的初始骨架。本发明采用无向图g(v,e)组织点云，图的存储结构分为数组表示法、邻接表法、十字链表法、邻接多重表法。本发明采用数组表示法存储点云，其中数组vertex存储数据元素信息（即顶点），数组edge存储数据元素之间的关系（即边）。本发明提取出的“类主干点”中的每一个点作为图中的一个顶点，搜索每一个顶点周围最邻近的k个点，连接该点和这k个点形成边，并将每条边两顶点间的欧氏距离度量作为这条边的权重。考虑到单株树木点云数量较多，所建立的无向图为稠密图，本发明采用prim算法，从上述无向图中构建mst，从而提取树木的初始骨架。mst是一种所包含边的权重之和是所有生成树中边的权重之和最小的树组织形式。因此要使得权重之和最小，mst算法会优先连接短边，可以有效地描述区域内部主干点的连接关系。
（3）优化初始骨架
由于点云存在噪声、离值点和数据缺失，初始树木骨架中存在不少细枝与主干的连接错误，需要对得到的初步骨架进一步优化，得到反映真实树木几何外观且具有较高保真度的树木骨架，本发明主要采用“点密度调整”和“树枝平滑”实现上述目标。
点云的密度我们基于ε邻近图进行计算，ε邻近图的定义为：两点i~j之间具有相互连接的关系，如果j∈b(i;ε)，其中b(i;ε)表示以i为中心以ε为欧式度量距离的原始点云的集合。我们基于ε邻近图组织原始树木点云，并根据公式（1）计算点云的密度：
di=n(b(i;ε))/ε（1）
其中n(b(i;ε))表示ε邻近图中与i点直接连接的点的数目，在本发明中ε取值为0.1m即可较好的平衡点云的噪声和骨架的精度。如果ε取值过小，则树木的骨架会受到部分误分的类主干点的影响，使得所生成的枝干中会产生过多的毛刺，为后续树枝平滑带来一定的压力；相反ε取值过大，则会对骨架产生过渡抽象表达，从而损耗几何表达的精度。生成mst时，采用的是欧式距离作为边的权重，点云密度的变化会影响mst的连接方式，因此需要修改图g(v,e)中边的权重。一般来说，一棵树的表面形态以及树枝伸展变化是连续且平滑的，例如在细小枝条的基础上很难生长出比其更加粗壮的枝条，相反，一根粗壮枝条具有较小的可能长出比它细小很多的枝条。根据这些特点，得到树点云的分布特征：（i）在不存在数据缺失的情况下，mst上两点间的点云密度是相似的，在较小的范围内，点密度具有连续性且过渡平滑。（ii）树枝区域点密度大，说明该位置存在粗壮的树干的可能性较大，反之，该位置存在一些细枝或数据缺失的可能性较大。
基于上述点云分布性质，本发明重新调整图g(v,e)中边的权重。首先修改无向图g(v,e)中每条边的权重，以适应点密度的变化。对应性质（i），相邻节点往往拥有相似的点密度，所以具有相似点密度节点间的权重小于非相似节点间的权重。对应性质（ii），优先连接点密度大的区域点云，继而判断点密度区域较小的点云是否与之连接，也就是说要优先在具有较大点密度或相似点密度的区域构建树枝骨架。基于上述思想将会极大地降低由于树枝遮挡或自遮挡造成的数据缺失对生成树木骨架的敏感性，同时也会使得错综复杂的小枝对生成树木主干结构的影响显著降低，此时，图g(v,e)中边的权重由公式（2）决定。
eij=dij/(di×dj)×|di-dj|（2）
参数di和dj分别是与点i和j相连的所有点组成的点集密度，dij是点i与j的欧式距离。一个点所在区域的点密度较大时，周围点到该点的距离就越小，反之就越大。i和j周围的点密度差距越大，|di-dj|差距就越大，边的权重eij越大；i和j周围的点密度差距越小，1/(di×dj)越大，边的权重eij就越大。公式（2）很好地满足树点云分布的两条性质。
点密度调整以后树的骨架结构虽然得到了一定程度的调整，但是树枝的伸展方式并不自然，需要对这些树枝进行平滑，从而符合真实树枝的伸展方式。本发明采用式（3）laplacian方程对点云密度调整后的骨架进一步优化。
参数wl和wh是对角矩阵，wl控制平滑的程度，wh控制平滑后的点云与原始点云相似的程度。wl中第i个对角元素的值定义成wl,i，wh中第i个对角元素的的值定义成wh,i。如果wl中的元素较大，那么平滑的程度较大；如果wh中的元素较大，平滑后的点云和原始点云相似度越高，平滑的程度越小。v是原始的点集，v'是优化过后的点集，l是n×n的拉普拉斯算子，通过公式（4）定义：
式中，e表示骨架结构上的边的集合，i、j、k表示骨架结构上的点。将对角矩阵wl的对角线上的元素值设为一致，这样可以使每个点具有相同的收缩能力，或者说，有相同的平滑能力，本发明中设为1。通过设置wh控制最后平滑的效果，矩阵wh对角线上的值是通过下述定义给出：
式中，k是i的父节点，j是i的子节点,θkij是∠kij。特别地，当点i表示的是末端点、根节点或交叉点时，将i单独设置成较大的值，使得这几类点和原始点相比变化最小，因为这些点的位置基本没有变化，保持了树骨架原来的形状，并很好平滑了这些点中间的点。
（4）确定枝干半径
提取得到最优树木骨架后，需要计算枝干的半径，从而利用opengl绘制具有真实感的三维树木几何模型。当前主流的计算枝干半径的方法主要有“圆柱拟合”和“异速生长模型”两类，本发明充分耦合这两类方法，在树木枝干的精确性和不确定性两方面取得平衡，从而实现对枝干半径的合理估算。
“圆柱体拟合法”的关键是确定圆柱体轴线的位置和方向以及圆柱体的半径。由于树木点云数据存在噪声、离值点，往往难以保证树枝局部范围内的所有点分布在一个圆柱体的表面。另外由于树枝之间存在一定程度的遮挡和自遮挡，这使得面对扫描仪方向的树木枝干获得的点云多，背对着扫描仪的方向枝干获得的点云较少，导致局部范围内树干点云难以拟合圆柱体。上述两方面的缺陷共同导致基于数据驱动思想的“圆柱体拟合法”对树木点云质量具有极强的敏感性。“异速生长模型”主要通过异速生长理论，计算出骨架节点的半径，从而完成树干的重构。该方法需要提供较准确的树种和胸径等系列参数才有望构建具有较高逼真度和几何外观的三维树木模型。由于“异速生长模型”是基于模型驱动思想，同时耦合了特定树种的异速生长理论，因此得到的树木枝干半径是个近似值，其精度往往没有“圆柱体拟合法”高。
在现实中，如果点云扫描的质量较高，树干的半径可以用圆柱体拟合的方法来进行较为精确的估计，通过骨架结构上相邻两个点之间mst上所有点进行拟合。但是如果这些点分布不均匀，且存在大量噪声、离值点和数据缺失，采用圆柱体拟合难以获得良好的结果，此时“异速生长模型”则为估算树枝半径提供了一种替代方案。基于上述分析，本发明结合“圆柱体拟合”和“异速生长模型”共同确定骨架结构上每一点对应的半径。具体而言，分别采用“圆柱体拟合”和“异速生长模型”计算构架结构上某一特点pi的半径，结果分别为ric和ria，然后检查ric是否满足如下公式:
如果满足，则令ric为当前点pi的真值半径，该过程计算完成后，本发明优化方程（7）计算骨架上其它未知点的半径，该式用已知点的真值半径作为正则项，采用“异速生长模型”构建约束项。
（5）添加树叶
获得每一节点的半径后，骨架模型就膨胀成三维树模型，就构建完成了一棵没有叶子的树木三维几何模型。通常情况下，如果一棵树存在叶子，叶子位于点密度相对较高，但半径却相对较小的位置的可能性较高。为了识别这些位置，针对每一个节点，根据公式（7）判断该节点是否存在叶子。
gi=di/ri（7）
参数ri代表每个节点的半径，di是与i点相连的所有点的点云密度。从上式发现，半径小的时候gi大，密度大的时候gi大，所以gi用于判断节点i是否存在叶子。所有节点获得了对应的gi值后，对这些值从大到小排序。设定一个百分比的阈值控制整棵树存在叶子的节点的比例，比如前30%的点存在叶子节点，在筛选后点上随机加上一些叶子，就建立了一棵完整的三维树模型。
- 还没有人留言评论。精彩留言会获得点赞！