一种稀疏点云环境下的单株树骨架提取算法

李巍岳，刘春，吴杭彬，孙伟伟

(1.同济大学 测绘与地理信息学院，上海 200092；2.现代工程测量国家测绘地理信息局 重点实验室，上海 200092；3.矿山空间信息技术国家测绘地理信息局重点实验室，焦作 454000)

1 引 言

树木三维建模一直是计算机图形图像学及仿真学研究的热点之一，其模型重构技术已在虚拟现实、林业调查、生态环境建设、古树名木保护以及农林业等方面有着广泛的应用。多数学者主要用生物形态法及几何构造法两种方法来构建模型。前者主要是从植物学的角度模拟树木形成、生长及衰败的过程[1]，但该方法没有考虑树木的形态及结构信息，所以很难构建出真实树木的三维模型；后者通过计算机图形学的手段，参照树木的拓扑结构信息来生成模型[2]，弥补了前者在建模上的不足，但该方法没有严格遵循植物学的规律。树木重构首先要提取树木的骨架信息[3]，树木的骨架信息对描述形状具有重要的意义，不仅能够保持树木拓扑结构、方向及局部特征，而且可以以直观、紧凑的形式有效地表示3D模型，是树木建模的关键问题之一[4-6]。骨架提取目前没有严格的概念定义，无法确定提取算法的优劣，这是因为骨架提取针对不同的研究领域，质量的评价目前还没有一个统一的标准[7]。在树木骨架提取的研究中，应满足骨架的拓扑一致性、居中性及细性等特点。已有的研究[8-10]多数基于图像及手工绘制草图的交互式建模，提取树木的骨架信息。这部分研究中，除了基于图像的建模外，大部分的研究都是建立虚拟三维树木，没有实现真实树木的数字化。由于许多树木形体高大，枝条众多，树叶与枝干间存在遮挡，单纯的利用图像很难实现树木准确的骨架提取。

激光扫描测距技术(Light Detection and Ranging，LiDAR)是一种快速、直接地获取地形表面模型的技术；与传统的光学及微波遥感不同，LiDAR能够快速、精确地获取地面特征在水平和垂直方向上的位置，克服了传统观测手段中树木部分结构枝叶间遮挡的问题，是目前能测定森林覆盖地区地面高程的可行技术之一[11-12]。近些年来，利用点云信息提取树木骨架的研究逐渐增多，Pyysalo等采用矢量多边形的方法，在点云中识别出云杉、松树、桦树三类树种，分别计算点云的算数中心及各点与中心点的连线、夹角，按照逐渐增大的角度顺次连接各点形成单木树的骨架[13]；Gorte等基于点云与栅格影像，分别选择滤波、数学形态学、最短路径等算法，在构建的体素模型中提取3D树木骨架[14]；Delagrange等在垂直方向上将点云按照1cm间隔进行水平切片，找出每一层数据点构成的中心，然后从下向上通过最小生成树将这些中心点连接得到树木骨架[15]；Bucksch等利用八叉树数据结构来组织点云数据，通过体素模型标注物体的延伸方向，再结合图论方法及约束规则来提取骨架[16]；Raumonen等采用圆柱体模型，在收集到的地面LiDAR数据中，快速自动地构建出单株树的树枝及树干[17]。目前针对激光扫描数据提取树木骨架的方法中，多数建立在三维网络模型或体素化模型的基础上。由于点云模型获取困难，噪声多，缺少模型的连接信息，在提取过程中，需要进行点云去噪，体素剖分等处理，过程费时，结果较难预料[18]；特别是在稀疏点云的条件下，影响的因素更多，严重地影响到树木的三维重建，因此有必要寻找一种方法能够解决上面提到的问题。本文根据树木的生长规律，分别阐述了单株树的主干、节点及枝条信息的获取方法，总结出在稀疏点云(地面及机载)条件下提取树木骨架的最优化方法。

2 单株树骨架提取算法

树木信息一般较为复杂，采用激光扫描技术所获得的点云受到漫反射和其他信息的干扰，点云比较稀疏且精度有限。虽然所获取的点云数据分布稀疏离散，但点云之间存在密度关联与较好的层次性。鉴于此，本文主要的思想是：大多数的激光点反映的是树木不同位置的叶子信息，叶子点云内具有团聚的特点，将每一个聚类中心取代周围的叶子点云，作为树木骨架的末端；还有部分点云信息反映了树木的主干及枝条，采用自由型曲线进行模拟，可以最佳的还原树木自身的结构，同时，控制自由型曲线走向的法向量得到了树木枝条的节点；最终得到单株树的骨架(图1)。

图1 树木骨架提取流程

2.1 模糊c-均值聚类

叶子点云数据较为离散，不同的样本类属存在不确定性，模糊c-均值分类方法使用目标函数将数据分成互相分离的类别，其主要算法为：

设数据集X={x1,x2,…,xn}⊂Rr，每个数据样本xi由m个特征定义，m选择根据诸克军等提出的遗传-迭代自组织分析技术(GA-ISODATA)来进行优选[19]。

对数据集X进行模糊分区，主要算法如下：

(1)

(2)

vi由m个特征来表述，对于点云中每个类的中心坐标按照下面公式计算：

(3)

j是特征空间中的一个变量，即j=1，2，…，m。通过经典模糊聚类方法，可以识别出三维离散点云中密度相对集中的区域(图2)，将离散的点组织为若干个组别。在树木识别中，通过距离的选择，可以初步判断骨架点云与叶子点云的分布。叶子点云中每个类的中心能够反映枝条的末端。

图2 基于模糊c-均值点云聚类

2.2 B样条曲线拟合

树木主干及枝干结构较为复杂，单纯采用直线段来模拟，不美观也不符合自然规律，常采用带有弧段变化的曲线来表现。同时，选取的曲线需要满足控制树木形态及易于实现的特点。B样条方法能够表示与设计物体自由型的曲线与曲面(图3)，是逆向工程中广泛应用的数学描述方法[20]，其曲线方程可写为：

(4)

其中，pi(i=0，1，…，n)为控制顶点，顺序连接成的折线成为B样条控制多边形；Ni，k(u)(i=0，1，…，n)称为k次规范B样条基函数[21]。

B样条基是多项式样条空间具有最小支承的一组基，可通过改变控制点个数设计一条曲线，不需要改变多项式的次数。鉴于上述B样条的优势，选用其作为主干及枝干重建曲线。

2.3 法矢变化

点云根据B样条拟合得到的任一点p(u，w)处，构造一个垂直与样条曲面的矢量，该矢量称为法矢。单位法矢可以根据该点处切矢pu和pw的矢量积得到。

单位法矢在几何造型中是不可缺少的，在本研究中主要用其方向衡量B样条曲线的弯曲程度，点云法矢方向变化率(Δn)大，曲线波动大。主干及枝干在分支或者节点处，Δn较大；根据给定的初始阈值，可识别出法矢方向变化率较大的位置(图4)。此外，规定法线方向在样条曲线走向的顺时针方向上。

图3 离散点B样条曲线拟合

图4 法矢变化及节点识别

3 试验及其结果分析

3.1 试验1：地面稀疏点云树木骨架提取

3.1.1 数据获取及处理

试验采用FARO公司的Focus 3D高精度三维激光扫描仪，该仪器以每秒最大976000点的速率可扫描最长为153.49m，具有三维虚拟重现、速度控制、高精确度及视野范围大等特点[22]。首先，设计各扫描站及标靶的位置，满足每测站之间至少有3个标靶重合，使点云数据能够统一到相同的仪器坐标系下；然后，对各标靶进行精确扫描，统一每个标靶的标识，满足相同的标靶在不同测站中的标识必须相同。通过这几个原则，对实地树木进行采集，得到点云图(图5)。数据中包括大量冗余，采用系统抽稀的方法，间隔50个点作为抽样间隔，其原理是在50个样本点中随机选择一个样本数据点，图6为压缩后的点云数据。通过这种方法来随机模拟稀疏点云的环境。

图5 地面LiDAR树木点云图

图6 抽稀后的树木点云图

图7 聚类后的树木点云图

3.1.2 骨架提取

利用模糊c-均值聚类方法对抽稀后的点云数据进行分割，根据式(3)，当距离取为2cm，GA-ISODATA算法叶子点云被分为50类时，得到较为理想的分割密度(图7，不同的颜色表示不同的聚类团)。其中，骨架点云与叶子点云通过不同的颜色初步得以区分；另外，不同叶子点云表现出一定的团聚特点。不同的聚类团由相应的枝条相串联，根据树木生理特点，其枝条的末端应分布在叶子点云聚类团的质心处(图8，质心用红色点表示)，提取每个聚类团的质心点为枝条节点，最外围的质心数据作为枝条末端。

图8 聚类团质心

将提取出的部分骨架点云，利用1.2节及1.3节的办法，进行法矢方向变化率判断及B样条曲线拟合。由于地面LiDAR主要以近水平方式扫描树木，树木枝干点云存在大量的冗余，骨架起始节点选择树木点云最下端形成的流形曲面中心(图9中A点)。调整曲线法向量变化的阈值，在搜索中遵循点与点之间的最短距离，选择进行B样条曲线拟合的点，完成该部分点云树木骨架的提取(图9)，在节点处(图9中B、C、D、E、F、G、H、I点)，法矢方向变化率较大。按照同样的方法，将部分点云树木骨架的末端与聚类团的质心进行法矢方向变化率判断及B样条曲线拟合，形成完整的树木骨架(图10)。该方法能够解决点云数据冗余及遮挡的问题，提取出的骨架能够真实地反映树木的拓扑结构。

图9 部分点云树木骨架提取

图10 完整的树木骨架

3.2 试验2：机载点云树木骨架提取

3.2.1 数据获取及处理

机载LiDAR数据采用LiteMapper 5600系统获取，其中激光扫描仪为Riegl LMS-Q560，波长为1550nm，激光脉冲的长度为3.5ns，激光脉冲发散角小于等于0.5mrad。获取的离散点云数据采用WGS84坐标系，UTM投影北半球6度带第50分带。其中，点云之间的平均点间隔为0.15m，平均点密度为6.5个/m2。机载点云数据反映了物体垂直范围的情况，在树木信息提取中，主要可以获取单木树高、冠幅、林区的平均数高、体积、密度等。

通过点云数据的滤波处理，分离出地面点与非地面点；根据植被与人工地物的反射特性及空间几何关系，利用数学形态学方法将植被及人工地物区分开；同时利用分水岭算法，分割了树木冠层高度模型，提取出表述单株树的点云数据。

3.2.2 骨架提取

机载LiDAR通过较远的距离及垂直扫描方式获取树木信息，其点云数据普遍稀少，且由于叶子的遮挡，缺少对枝干特征的点云数据描述。试验中利用相关文献[23-24]里提到的树木位置、树干的提取方法：通过低通滤波处理的树木冠层高度模型中观察局部最大点，从而得到树木的位置；每棵树的树冠由平行的水平多边形包围，树干由地面到树顶的垂直矢量表示(图11)。利用模糊c-均值对提取的树木点云进行聚类，当距离取50cm叶子点云被分为30类得到较为理想的分割结果(图12)。

图11 树木位置、树干识别

图12 聚类后的树木点云分布图

选择离树干最近的4个点A、B、C、D作为B样条曲线模拟的起始点，参照1.3节的方法追踪法矢方向变化率，在搜索中遵循点与点之间的最短距离，与聚类团的质心拟合得到树木一级枝条；然后从一级枝条的下一个节点再进行法矢方向变化率及最短距离搜索，与聚类团的质心拟合次一级枝条，最终遍历整个点云数据，得到了机载条件下的树木骨架(图13)。

图13 机载条件下提取的树木骨架

3.3 定性分析

通过地面及机载方式获取的树木信息方式不同，前者以近水平的方式扫描树木，激光波长较小，树木的骨架信息比较完整，但存在数据冗余及部分点云信息遮挡的问题；后者以近垂直的方式扫描树木，激光波长较大，主要可以反映树木的高度、冠幅信息，骨架信息严重缺失，影响到后期的树木建模。

树木骨架提取目前较为成熟的方法是参数化L-系统[25]，传统的参数化L-系统是根据树木的类型选定合理的公理和产生式，设计不同的初始状态与树木描述规则，通过不断地迭代产生一系列字符串，然后逐一读取字符，按照执行规则的不同逐一对读取的字符串进行解释。针对机载LiDAR树木建模，Pyysalo[13]等在2002年提出了一种构建树木矢量模型的方法，依次按照树冠之间的1m间隔将整个点云分成若干层，在每个高度层内的点按递增的方位角排列，从而得到树木的水平多边形，该方法已经投入到了大范围林木林地模型的生产中。将本文提出的方法与这两种方法进行比较，结果如表1所示。本文的方法在算法上简单、易于实现；支持地面及机载LiDAR数据的树木骨架提取，骨架具有一定的弯曲度，接近真实树木的情况；适用于压缩后的数据，可以为单木树及大范围树木建模提供一定的理论基础。

表1 不同树木骨架提取方法比较

4 结束语

骨架提取在图形可视化及建模的应用中是一个比较基本的问题，它能够反映物体的结构特征及拓扑信息。稀疏点云数据的精度较低，存在树木数据结构缺失的情况，选择一种方法最大限度地提取出树木骨架信息至关重要。

本文首先基于模糊c-均值聚类的方法，对点云数据进行分割，根据树木的枝叶分布特性，选取聚类团的质心拟合枝干形状的自由型结构，并通过曲线的法矢变化率判断树木的节点，最终实现了树木骨架的提取，并分别通过地面及机载LiDAR数据的试验，验证了该方法适用于稀疏点云条件。

参考文献：

[1] 程章林，陈宝权.大规模树木三维建模研究[J].先进技术研究通报，2009，3(6)：40-44.

[2] SU Z X，ZHAO Y D，ZHAO C J，et al.Skeleton extraction for tree models[J].Mathematical and Computer Modelling，2011，54(3-4)，1115-1120.

[3] LIVNY Y，YAN F，OLSON M，et al.Automatic reconstruction of tree skeletal structures from point clouds[J].ACM Transactions on Graphics.2010，29(6)：1-8.

[4] AU O K C，TAI C L，CHU H K，et al.Skeleton extraction by mesh contraction[J].ACM Transactions on Graphics，2008，27(3)：1-10.

[5] HENNING J G，RADTKE P J.Detailed stem measurements of standing trees from ground-based scanning LiDAR[J].Forest Science，2006，52(1)：67-80.

[6] 张义宽，张晓鹏，查红彬，等.3维点云的拓扑结构表征与计算技术[J].中国图象图形学报，2008，13(8)：1576-1587.

[7] 黄文伟.基于拉普拉斯算子的点云骨架提取[D].大连：大连理工大学，2009：2-4.

[8] TENG C H，CHEN Y S，HSU W H.Constructing a 3D trunk model from two images[J].Graphical Models，2007，69(1)：33-56.

[9] CHENG Z L，ZHANG X P，THIERRY F.Tree skeleton extraction from a single range image[C].The 2ndInternational Symposium on Plant Growth Modeling，Simulation，Visualization and Application，2006，274-281.

[10] OKABE M，OWADA S，IGARASHI T.Interactive design of botanical trees using freehand sketches and example-based editing[J].Computer Graphics Forum，Proceedings of Eurographics，2005，24(3)：487-496.

[11] 刘春，陈华云，吴杭彬.激光三维遥感的数据处理与特征提取[M].北京：科学出版社，2010：8-14.

[12] 徐祖舰，王滋政，阳锋.机载激光雷达测量技术及工程应用实践[M].武汉：武汉大学出版社，2009：2-6.

[13] PYYSALO U，HYYPPA，H.Reconstruction tree crowns from laser scanner data for feature extraction[R].International Society for Photogrammetry and Remote Sensing-ISPRS Commission III Symposium，Graz，Austria，2002.

[14] GORTE B，PFEIFER N.3D image processing to reconstruct trees from laser scans[R].Proceedings of the 10thAnnual Conference of the Advanced School for Computing and Imaging(ASCI)，Ouddorp，the Netherlands，2004.

[15] DELAGRANGE S，ROCHON P.Reconstruction and analysis of a deciduous sapling using digital photographs or terrestrial-LiDAR technology[J].Annual Botany，2011，108(6)：991-1000.

[16] BUCKSCH A，LINDENBERGH R，MENENTI M.Skeltre：Robust skeleton extraction from imperfect point clouds[J].The Visual Computer，2010，26(10)：1283-1300.

[17] RAUMONEN P，KAASALAINEN M，AKERBLOM M，et al.Fast automatic precision tree models from terrestrial laser scanner data[J].Remote Sensing，2013，(5)：491-520.

[18] 刘雪梅，庄晋林，张树生，等.利用自适应模糊椭球聚类实现点云分区[J].计算机工程与应用，2007，43(15)：33-35.

[19] 诸克军，苏顺华，黎金玲.模糊C均值中的最优聚类与最佳聚类数[J].系统工程理论与实践，2005，(3)：52-61.

[20] 成思源，谢韶旺.Geomagic Studio 逆向工程技术及应用[M].北京：清华大学出版社，2010：19-21.

[21] WU Z K，ZHOU M Q，WANG X C.Interactive modeling of 3D tree with ball B-spline curves[J].The International Journal of Virtual Reality，2009，8(2)：101-107.

[22] AZMY S N，SAH S A，SHAFIE N J，et al.Counting in the dark：Non-intrusive laser scanning for population counting and identifying roosting bats[J].Scientific Reports，2010，(2)：524.

[23] YU X W，HYYPA J，HARRI K，et al.Automatic detection of harvested trees and determination of forest growth using airborne laser scanning[J].Remote Sensing of Environment，2004，90(4)：451-462.

[24] HYYPPA J，INKINEN M.Detecting and estimating attributes for single trees using laser scanner[J].The Photogrammetric Journal of Finland，1999，16：27-42.

[25] 石银涛，程效军，张鸿飞.基于参数L-系统的三维树木仿真[J].同济大学学报(自然科学版)，2011，39(12)：1871-1876.