黄志鑫，邢 涛，邢艳秋，王 强

（东北林业大学 森林作业与环境研究中心，黑龙江 哈尔滨 150040）

骨架曲线以紧凑和直观的方式共同描述实体形状的几何、拓扑和对称等特性，不仅有助于分析空间结构信息，而且能够行之有效的降低处理过程中的复杂度[1-3]。由于森林调查涉及大量数据，在不影响林木基本信息的情况下，引入单木骨架曲线概念可望精简数据量[4]。它是树木三维几何模型构建的基础，在精准林业等研究领域有重要的意义[5]。

激光雷达（Light detection and ranging，LiDAR）技术是最近20年快速发展的新型遥感技术，能够探测森林植被空间结构特征，特别是森林垂直结构的测量能力，具有以往光学遥感设备难以比拟的优势[6]。近年来，背包激光扫描（Backpack laser scanning，BLS）凭借其灵活机动性，能够弥补现有激光雷达技术不能或不便深入复杂样地进行测量的局限性[7-8]，具有提高数据采集效率的潜力[9]，作为一种新型扫描概念被Kukko 等[10]和Hyyppä 等[7]引入，并于2014年首次应用于森林调查[6]。背包激光雷达采集数据通常以点云形式存在，针对点云骨架曲线提取，主要研究方法具体可分为：Runions 等[11]通过使用点云作为局部吸引子，在封套内生长骨架曲线结构。封套方法能够很好提取具有单一结构的点云形式，但是对于复杂结构进行骨架曲线提取不易实现。Bucksch等[12]利用八叉树对点云进行分割，根据八叉树各节点的相邻信息建立八叉图得到点云骨架曲线。此方法在中，八叉树深度对该骨架曲线提取质量起着关键作用。Cao 等[13]将拉普拉斯算子应用到点云模型，利用拉普拉斯收缩算法对点云收缩拓扑细化成骨架曲线。该算法对骨架曲线提取参数非常敏感，曲线质量依赖参数调节。Sheng 等[14]通过对自适应采样选择的骨架关键点连接形成语义器官组成骨架曲线。但关键点的选取建立在拉普拉斯收缩算法基础之上。

随着背包激光扫描的兴起，针对背包激光点云数据处理的研究不断深入[15]。上述的实验研究大多针对地面激光扫描(Terrestrial laser scanning，TLS)数据的骨架曲线提取，TLS 数据采用固定式扫描，数据质量相对较高，能够很好的用于骨架曲线提取。而BLS 设备通常在行进中对目标进行数据获取，由于目前采集过程中定位技术的局限性，扫描数据投影到二维平面表现为相当厚度的环状点云，并且存在大量离散点，使得数据质量不佳。鉴于此，本研究以背包激光雷达扫描获取的活立木数据为基础，着眼于林木单枝点云近似圆柱特性，重建与活立木形态结构相一致的骨架曲线，提出一种基于动态调整拟合圆柱半径的骨架曲线提取方法。

1 研究区概况及研究方法

1.1 研究区概况

研究区位于广西壮族自治区南宁市西北郊高峰林 场（108°8′～108°53′E，22°49′～23°15′N），如图1所示。地处南宁盆地边缘，地势东高西低。地貌特征主要以丘陵（占地面积55.5%）和山地（占地面积38.7%）为主，海拔高度150～485.8 m。属温润的亚热带季风气候，常年气候温和，全年雨量充足。林场主要树种为马尾松Pinus massoniana、速生桉Eucalyptusrobusta smith、杉木Cunninghamia lanceolate等。

图1 广西南宁高峰林场研究区域Fig.1 Study area in Gaofeng forest farm,Nanning,Guangxi

1.2 实测数据

本研究于2017年12月对样地进行外业测量，实测数据如表1所示。选取林下环境良好，地势相对平坦的杉木林区域，规划样地大小为30 m×30 m 矩形，样地内有杉木48 株，人工测量胸径均值17.22 cm，最大树高18.7 m，最小8.2 m，平均树高14.47 m。

表1 样地实测数据Table 1 Data of field measurement

1.3 BLS 数据采集

本研究BLS 数据于2017年12月使用背包激光雷达扫描设备采集，设备参数见表2，采集路线如图2所示。采集过程中选取样地一角点（五角星处）作为扫描起始点，测量人员背负设备以“S”形在样地中行进扫描，回到起点使扫描路线形成闭合，完成采集获取样地背包激光雷达点云数据。

图2 数据采集路线Fig.2 Data acquisition route

表2 背包激光雷达扫描技术参数Table 2 Technical parameters of BLS

1.4 研究方法

本研究单木骨架曲线提取技术路线如图3所示，研究方法具体分为3 个阶段。阶段1：获取待提取的树状点云数据。从背包点云数据集中分离出单株杉木点云，对获取的单木点云进行树叶点与地面点滤波，并对单木枝干点进行离散点去除。阶段2：从枝干点云中提取单枝骨架点。采用动态圆柱拟合林木单枝，对拟合出的单枝点云沿主方向进行逐层聚类，获取相应的单枝骨架点。阶段3：获取单木骨架曲线。采用3 次B 样条曲线对单枝骨架点进行平滑处理，再利用同向向量近邻搜索，对单枝骨架曲线实现拼接，最终得到完整单木骨架曲线。

图3 研究技术路线Fig.3 The workflow of research

1.4.1 数据预处理

选用研究区单株杉木点云数据进行处理，由于地面点和树叶点与骨架曲线提取无关，对其进行相应滤波。局部地面点云表现为面状，采用Sahin 等[16]和Dong 等[17]提出的自动3D 平面分割方法实现地面点的去除。由于树木枝干与树叶对激光反射强度不同，获取的点云密度不同，则利用Louhichi 等[18]提出的基于密度的聚类算法实现树叶点云的滤除。

1.4.2 动态圆柱拟合

根据单枝点云近似圆柱特征，基于RANSAC(Random sample consensus)进行圆柱拟合[19]，该方法能够从一组观测数据中，通过多次迭代方式滤除“局外点”，进而估计出圆柱模型的参数。本研究根据树干近地面处的最大半径确定圆柱模型的半径起始值，限制每个局内点到圆柱模型距离阈值，设置法线影响权重值，实验采用的圆柱拟合模型如式（1）所示。

式(1)中：(x0，y0，z0)为圆柱轴线L上一点，(l，m，n)为圆柱轴线L方向向量，r为圆柱半径，此7 参数可以确定唯一圆柱方程，（x，y，z）为点云的三维坐标。

对于整株单木点云，如图4a所示，由于单枝粗细大小不同，圆柱拟合所需半径亦不同，无法采用固定半径进行单木全局单枝拟合。此外，由于背包激光雷达目前技术的局限性，如图4b所示，背包激光扫描单枝点云数据横截面呈近圆环状分布，故需要对圆柱拟合半径进行动态调整，从而实现圆柱拟合之后获得粗细大小不同的单枝。

图4 动态参数调整条件Fig.4 Conditions of dynamic parameter adjustment

本研究提出动态圆柱拟合单枝方法，通过设置半径递减参数δ，在拟合过程中动态调整圆柱半径大小实现单枝点云拟合。为达到充分利用待处理点云的目的，多次设置δ值进行单枝圆柱拟合，通过对拟合单枝圆柱点云结果进行判定，从而确定最优参数值。

BLS 点云数据呈圆环状分布，横截面如图5a所示。若一次将圆柱完全拟合，即拟合出半径为r的圆柱R 可将枝干点云完全拟合。但由于存在不同粗细枝干，需要不同拟合半径，故无法一次性实现圆柱拟合，需要对点云进行迭代拟合。

单枝点云迭代拟合过程中，首先根据最大拟合圆柱半径r1，如图5b所示，通过设置内点到圆柱模型的距离最大值d/2，获得圆环状柱体R1，即圆C1与圆C2构成的圆环。

图5 单枝截面图Fig.5 Cross section of sticks

通过设置递减参数δ，能够实现对圆柱半径调整。进行第一次圆柱拟合之后，圆柱半径r1已经无法满足利用剩余点云进行拟合，对r1与δ进行做差，所得结果即为第二次圆柱拟合半径r2，最终获取第二个圆柱体模型R2。通过此种方法对剩余点云迭代拟合，充分利用点云数据，直至将单枝点云全部拟合。

1.4.3 单枝点云优化

在圆柱拟合过程中，δ是实现拟合圆柱半径动态调整的重要参数，其取值将直接影响单枝拟合效果。为充分实现单木单枝圆柱提取，采用较小的递减参数值，以期达到更精细拟合单枝的目的，但随之而来产生两大问题：同一单枝多次拟合与不同单枝误分，如图6所示，图7中相同颜色的点云为同一圆柱拟合结果。这些现象会严重影响后续骨架点的提取。为解决这一问题，需要对拟合完成的单枝点云进行以下优化。

图6 单枝拟合结果（红色与蓝色点云分别代表不同单枝）Fig.6 Results of sticks’ fitting(red and blue point clouds represent different branches respectively)

1）圆柱合并。对于执行多次拟合的同一单枝，其拟合后的点云呈嵌套分布并且依旧保持圆柱的特性，这些圆柱半径大小存在差异，但具有相同或相近的轴线方向。依此属性，对嵌套圆柱利用参数信息计算圆柱轴线夹角实现多个拟合圆柱合并，夹角方程如式(2)所示。

2）误分单枝剔除。为解决非同一单枝错误分为同一单枝点云问题，引入Zhao 等[20]提出的基于密度的聚类算法(DBSCAN)，通过该算法能够找到任意形状的聚类。通常情况下，误分点云所占数量比例很小，并且其原有单枝已能够实现相应单枝的骨架点的提取，可直接对误分重合点云部分进行剔除。

1.4.4 骨架曲线提取

骨架点的提取采用沿单枝点云主方向基于体元逐层聚类[21]。通过对单枝整体点云的XYZ 坐标进行主成分分析可得到3 个特征值[22]（λ1≥λ2≥λ3≥0），最大特征值λ1所对应的特征向量即为该单枝的主方向。利用模糊C均值算法对底层体元内的点云数据进行分类以确定体元质心，以此作为起始聚类中心，沿单枝圆柱主方向逐层聚类，每层聚类中心即为该层骨架点。

采用3 次B 样条曲线拟合方法[23]对单枝骨架点进行自然衔接，能够使得拟合后的骨架曲线具有C²连续性[24]。从而满足骨架曲线对平滑的要求，有利于更准确的表现单木骨架曲线的流畅性。

由于单枝拟合得到的骨架曲线并没有进行相互连接，不具有整体效果，未能达到反映出各级单枝之间的关系，故需要对各个骨架曲线进行连接。采用基于同向向量近邻查找的方法将所有单枝骨架曲线连接成完整单木骨架曲线，从而实现整体骨架线的提取。

2 结果与分析

2.1 骨架曲线提取

本研究对研究区单株杉木进行骨架曲线提取，如图7所示。实验结果表明，利用该方法能够提取出完整的骨架曲线模型。对比原始点云，所得骨架曲线具有很好的一致性与还原性。图7a 为单木原始点云数据，存在大量的地面点云与树叶点云。图7b 为单枝圆柱拟后的点云，图7c 为单枝骨架曲线构成的完整单木骨架曲线，提取的骨架曲线基本还原图7b所示单木的结构特征，并且能够直观的获得枝干之间的拓扑特征。

图7b 中不同颜色的点云代表不同的枝干类型，对于存在较大弯折部位以及明显未连接的单枝，由于拟合圆柱轴线不同，圆柱拟合算法将点云拟合成两个不同的单枝点云。单枝点云分布明显，没有错分为同一枝干的点云。图7c 中不同颜色的曲线代表不同的单枝骨架曲线，单枝骨架曲线提取良好，未有错误拼接状况，单枝骨架曲线之间连接准确自然，整体曲线能够表征完整单木的骨架曲线。

图7b 与图7c 对比，单枝点云与单枝骨架曲线一一对应，每一个单枝点云都能够正确地提取出骨架曲线。通过对单枝骨架曲线进行拼接，构成单木完整骨架曲线，连接部位准确自然，结果表明，该方法具有一定的鲁棒性，可以看出单木骨架的拓扑结构与树状原始点云结果基本保持一致。

图7 骨架提取过程Fig.7 Process of skeleton extraction

2.2 动态调节参数

动态圆柱拟合所得单枝点云是后续骨架提取关键，故在本实验中对拟合圆柱半径递减参数δ进行研究。图8分别显示不同递减参数的拟合结果，其中蓝色点表示圆柱拟合前的枝干点云数据，绿色点表示圆柱拟合之后的结果。图中绿色点所占比例越小，即点云利用率越低，表示圆柱拟合对点云剔除越明显，单枝骨架点提取效果越差；绿色点所占比例越大，待处理点云利用率越高，则单枝拟合效果越好，骨架点提取准确度也越高。

图8a 为参数0.3 cm 时拟合出的结果，绿色点与蓝色点数量百分比为95.04%，点云利用率最高，图中绿色点基本将蓝色点云完全遮挡，说明拟合结果能够最大程度利用待处理点云，点剔除现象不明显。对于单枝点云，设置递减参数为0.3 cm，能够非常完整地将各个单枝进行圆柱拟合。图8b为参数0.5 cm 拟合结果，数量百分比为92.70%，点云利用率较高，绿色点基本覆盖蓝色点，单枝点云存在部分点被剔除现象，但是每个单枝绿色点所占比例较大，被剔除的部分点并未对单枝骨架点提取造成影响。图8c 为参数0.8 cm 拟合结果，绿色点数量与蓝色点百分比为79.70%，点云利用率低，图中蓝色点明显，并且存在单枝点完全为蓝色现象，大量点因被剔除而不能参与骨架点提取，并且存在一些单枝无法提取骨架点的情况。

图8 不同参数圆柱拟合点云（绿色）与未处理点云（蓝色）数量百分比Fig.8 Percentage of cylindrical fitting point clouds with different parameters(green) and unprocessed point clouds (blue)

由图8可知，参数δ值为0.5 cm 时，圆柱拟合单枝点云能够基本保留原始单枝，在后续骨架曲线提取时，能够相对完整的提取出较多细节的单木骨架曲线；若充分获得圆柱拟合单枝点云，可以依靠过拟合实现，即设置递减参数δ值为0.3 cm；若在处理过程中过分考虑运算效率，可以通过欠拟合实现，即设置递减参数δ值为0.8 cm，但是这种情况存在部分单枝点云缺失严重现象。

2.3 参数对比结果

本研究设置过拟合参数0.3 cm和参数0.5 cm、欠拟合参数0.8 cm，对3 种情况进行详细对比。对同一株杉木点云数据进行3 组递减参数圆柱拟合，待处理点数量为62 967，具体数据见表3。如图8b所示，选用参数值为0.3 cm 时，即实现单枝点云过拟合，圆柱拟合后点云数量为59 843，运算时间为515.05 s，能够最大限度的保留最初单木点云数据。但由表3可知，良好的拟合效果建立在增加运算时间的基础上。图8b 已基本能够保留待处理单木点云单枝信息，当参数为0.5 cm 时，拟合后点云数量为58 372，运算时间为369.34 s，相比于参数为0.3 cm 的结果，拟合后的单枝点云能够满足骨架曲线提取，且节约了28.29%运算时间。与图8c 相比，取参数值为0.8 cm 时，对单枝点云进行欠拟合，拟合后点云数量为50 184，运算时间为349.35 s，相比于参数为0.3 cm，运算时间节约了31.17%。但是存在较多单枝点云缺失的情况，采用δ值为0.8 cm 无法满足提取骨架曲线的要求。

表3 不同半径递减参数拟合结果Table 3 Fitting results of different radius parameters

3 讨 论

点云骨架曲线提取是林业可视化过程中的基础研究内容之一。本研究提出了一种基于动态调整圆柱半径进行圆柱拟合的骨架曲线提取方法。首先对单枝点云进行骨架点提取，构建单枝骨架曲线；再将单枝骨架曲线进行向量搜索点拼接，构成整体单木骨架曲线，所提取的骨架曲线基本恢复了单木树状点云的拓扑形态特征。

对单株杉木点云采用动态圆柱拟合方法进行骨架曲线提取，所得效果如图7所示。结果表明，提取出的单木整体骨架曲线不仅形状上与单木点云具有很好的结构吻合度，而且骨架曲线具有良好的整体性。依据实验结果，本研究所提出的骨架曲线方法切实可行。

本研究提出动态圆柱拟合单枝为前提的骨架曲线提取方法对圆柱半径递减参数进行研究。采用不同参数值进行圆柱拟合，对结果进行分析，在充满考虑算法运行效率、圆柱拟合结果的前提下，参数为0.3 cm 时，拟合后点云数量为59 843，相比待处理点云数量，点云利用率为95.04%，点云利用率非常高，拟合效果最好，但是运算时间为515.05 s，处理时长表现不佳。参数为0.8 cm 时，运算时间为349.35 s，相比参数0.3 cm 减少了31.17%，运算效率表现最好，但点云利用率为79.70%，且因部分单枝缺失而无法提取骨架曲线。参数为0.5 cm 时，点云利用率为92.70%，运算时间相比参数0.3 cm 降低了28.29%。相比于0.8 cm，点云利用率能够保证所有单枝的骨架曲线提取，满足研究所需条件。

与现有骨架提取算法对比，基于动态圆柱拟合单枝进行骨架曲线的提取方法的算法相对简单。相比于Alexander 等方法，本研究方法不需要预先考虑八叉树深度，充分利用PCL 库即可实现该算法；本实验采用背包激光雷达林地点云数据，数据质量相对较好，但也会存在单枝密集处点云稀疏的状况，Sheng 等[14]需要足够稠密的点云，通过自适应采样确定骨架关键点形成语义器官，本算法充分利用单枝点云结构，不需要稠密的单木点云数据就能够基本保证单木单枝拓扑结构的合理性。

基于圆柱拟合单枝点云进行骨架曲线提取直接作用于背包激光雷达点云，现有算法大多对点云进行网格化处理。Cao 等[13]将点云网格化，利用拉普拉斯收缩并拓扑细化形成骨架曲线。但是，在网格化过程中，网格插值方法导致的不确定性将影响单枝细化的准确度；其次，网格化骨架曲线提取没有利用激光雷达点云完整的三维结构。本研究将林木单枝点云包含于其对应的拟合圆柱体内，不必对点云进行网格化处理，并且通过对阈值参数的动态调整，降低离散点对提取准确度的影响，从而达到完全利用目标点云提取骨架曲线三维结构的目的。

在验证算法可行性的前提下，本研究使用单株杉木作为研究对象进行骨架曲线的提取，研究样本过于单一。因此，在接下来的研究中，我们将继续使用多种林木背包激光雷达点云提取骨架曲线，以期能够更好的验证本算法的鲁棒性。

4 结 论

本研究使用单株杉木背包激光雷达点云进行骨架曲线提取，利用林木单枝近似圆柱生长特性，基于动态圆柱拟合提取骨架曲线，得出如下结论：

1）利用枝干自然生长近似圆柱的特性，本研究基于动态圆柱拟合方法能够实现骨架曲线的提取，提取结果具有很好的一致性与还原性，能够恢复单木树状点云的拓扑形态特征。

2）相比于地面激光雷达，背包激光雷达点云具有更高的离散性，数据质量弱于地面激光雷达点云。为克服数据本身的问题，通过设置动态调节参数，从单个枝干点云中拟合多个圆柱体点云以获取完整的单枝点云，能够实现背包激光雷达点云数据的单木骨架曲线。

3）本研究方法直接作用于单木枝干点云，不必对点云进行八叉树及网格处理，半径递减参数为0.5 cm，点云利用率为92.70%时，所得骨架曲线具有良好的整体性。