1 引言

近年来,地面三维激光扫描作为一种新的测量技术已被广泛应用于城市建模和森林调查等领域^[1-5]。地面三维激光扫描仪采用非接触的方式获取物体表面采样点的距离信息,通过计算距离与角度信息获得点云数据,并进行数据处理与分析^[6]。地面三维激光扫描仪因数据精度高和分辨率高等特点也常被用于林业测量,可以辅助人工测量,提高树木信息(如胸径)的采集效率^[7]。树木骨架是对树木三维结构的精简描述,其可表示树木的形状和拓扑结构,根据生成的树木骨架,可进一步计算出树木的几何结构、分枝模式和分枝角度等信息,从而体现树木的结构和生长特点。

目前,已有不少学者在基于三维点云生成骨架方面进行了一些研究工作^[8-16],其中包括一些生成树木骨架的研究成果。树木骨架的生成方法大致可以分为两种。一种是通过分割空间,得到体素、元胞等新的计量单位,进行骨架节点的定位与计算,最后连接节点生成树木骨架。Bucksch等^[17-18]提出了一种基于八叉树的骨架生成方法,生成八叉树之后,从八叉树元细胞中提取一个由点组成的连接图,然后消除图中的闭环,形成一条直线段,最后得到树木骨架。Ramamurthy等^[19]提出了一种基于体素的植物骨架生成算法。首先设定体素的体积大小,并计算点云数据所属的体素,通过细化算法删除冗余的体素后,利用近邻搜索连接得到植物骨架。另一种是基于点云局部分割,计算骨架节点,然后根据图论或法向量等方法连接骨架节点,生成树木骨架。Yan等^[20]将提取的树木骨架分成分割与重建两个部分,该方法基于K-均值聚类算法实现区域识别和定位,然后建立邻域图提取骨架节点和树木分支,最后通过B样条模型计算生成树木骨架。Xu等^[21]提出了一种利用点云生成数骨架的基于最短路径和聚类的算法,该方法首先根据每个点的邻接关系构造邻域图,然后用最短路径算法计算各点与根点之间的距离,并基于距离聚类生成聚类块,最后利用聚类块内的点云拟合骨架节点并连接,从而得到树骨架。Delagrange等^[22]扫描获得没有叶子的树木点云数据,并沿着Z轴提取210层宽度为1 cm的切片,从每个切片中计算出质心,再根据最小生成树算法连接质心生成树木骨架。唐丽玉等^[23]基于水平采样和最小二乘法获取枝干骨架点,根据骨架节点与特征点的拓扑结构重构树木三维模型。赵春江等^[24]基于立体视觉原理,通过两幅特定角度的图像先提取二维骨架,然后通过极线约束算法和合理的匹配准则,最终生成玉米植株的三维骨架信息。李巍岳等^[25]提出了一种通过稀疏点云生成树木骨架的方法。该方法主要通过模糊C-均值分割点云,然后使用B样条拟合聚类团质心,根据法向量方向变化率确定节点位置,从而生成树木骨架。Li等^[26]基于图论的概念对点云数据进行建图,采用近邻搜索算法对每个点进行标号,再进行非线性量化处理,使数据分层形成聚类块,然后根据聚类块中的点云拟合骨架节点,最后根据深度搜索算法连接骨架节点,生成树木骨架。

通过对上述算法的研究分析,结合树木点云的特点,本文基于体素空间细化处理的思想,提出了一种体素细化改进的树木骨架生成方法:1)通过树木点云的范围和体素的数目构建体素空间;2)根据滤除阈值对噪点体素进行滤除;3)对目标体素进行细化处理,并利用细化后的体素拟合骨架节点;4)基于树木的空间结构和深度优先搜索算法连接骨架节点,生成骨架。该方法提出了一种体积形状可变的体素,用于构建体素空间,同时增加了对体素孤立分布程度的分析,然后利用体素细化算法提取树木的骨架,大大提高了树木骨架生成的运行效率。

2 材料与方法

本文方法处理流程如图1所示,主要包括体素空间构建、噪点体素滤除、细化算法处理和树木骨架生成4部分。

1) 体素空间构建。计算树木点云所占空间大小,记为体素空间。按照输入参数将体素空间(X、Y、Z三个维度上的体素个数)分成若干个体素。算法中约定包含树木点云的体素为目标体素,包含非树木点云的体素为噪点体素。

2) 噪点体素滤除。根据体素中点云的密集程度和体素的分布情况,通过设定滤除点数阈值和孤立体素滤除阈值,对噪点体素进行判断与滤除。

3) 细化算法处理^[27-32]。噪点体素滤除完成之后,对目标体素进行细化处理,通过迭代删除边缘的目标体素,最后只保留中间呈线状连接的目标体素。

4) 树木骨架生成。根据细化处理后的体素内的点云数据拟合骨架节点,并根据树木结构与近邻搜索算法连接骨架节点,得到最后的树木骨架。

图 1. 方法流程图

Fig. 1. Flow chart of proposed method

下载图片 查看所有图片

2.1 数据获取

本次实验使用的地面三维激光扫描仪器为RIEGL VZ-400型地面三维激光扫描仪,扫描距离最远可以达到400 m,每秒最多可以获取4.2×10⁴个点,垂直扫描角度范围为-40°~+60°,水平扫描角度范围为360°。数据采集于北京林业大学校内的一棵银杏树(Ginkgo)和一棵重瓣榆叶梅树(Amygdalus triloba f. multiplex tree),采集时间为2018年1月23日。实验过程中,对银杏树和重瓣榆叶梅树分别进行了多次扫描,并通过配套的软件RISCAN PRO手动截取样本树木的点云数据,删除周围噪点。数据对应的具体角分辨率和点数如表1所示。

2.2 构建体素空间

通过遍历点云,获取数据在X、Y、Z三个维度上的最大值(X_max,Y_max,Z_max)和最小值(X_min,Y_min,Z_min),从而构建体素空间。输入参数M、N、K,M、N 、K分别代表体素划分的行数、列数和层数,即将体素空间划分为M×N×K个体素。体素大小(V_{x_size}、V_{y_size}、V_{z_size})计算公式为

Vx_size=(Xmax-Xmin)MVy_size=(Ymax-Ymin)/NVz_size=(Zmax-Zmin)K。(1)

显然,树木点云的位置在空间三个维度上不同,故对于可调整的输入参数M、N、K,V_{x_size}、V_{y_size}、V_{z_size}可以设置为相同,也可设置为不同,从而保证体素划分的灵活性。

确定体素大小后,设定体素所在的行数、列数和层数作为体素的坐标。则点云所属体素计算公式为

Vx=int(x/Vx_size)Vy=int(y/Vy_size)Vz=int(z/Vz_size),(2)

式中:(V_x,V_y,V_z)为点云所在体素的位置坐标;(x,y,z)为点云的空间坐标;int()函数为取整函数。

2.3 噪点体素滤除

扫描设备、外部环境、树木特征和枝叶遮挡等会导致噪点的产生,这可能造成骨架的冗余与缺失,影响树木骨架的生成。本研究从两个方面对包含噪点的体素进行判断与滤除。

1) 基于体素内点云数量滤除体素。体素内点云数量过少,说明这些点可能是噪点。在算法中设定滤除阈值ε₁,当体素中点数S_num小于ε₁,则认为该体素是噪点体素并删除。根据角分辨率和点到原点的距离可以计算出体素内充满点时点的个数,选取一个小于1的比例系数,即得滤除阈值ε₁,可表示为

ε1=a×Vx_sizeRx×Vz_sizeRz,(3)Rx=2L×sinθ2Rz=2L×sinϕ2,(4)

式中:R_x和R_z分别为点云中X,Z方向上两点间距;L为点到原点的距离;θ和ϕ分别为水平方位角和俯仰角的扫描角分辨率,a(a<1)为滤除比例。

比较阈值和体素内的点数,将体素区分为目标体素和噪点体素,并进行标记。l_abel为1,表示该体素为目标体素;l_abel为0,表示该体素为噪点体素:

label=0,Snum<ε11,Snum>ε1。(5)

2) 根据目标体素分布的孤立程度滤除体素。若一个目标体素周围26邻域中的目标体素过少,则该目标体素有较大可能是噪点体素。计算目标体素周围26邻域中目标体素数量v_num。设定孤立程度判断阈值为ε₂,如果v_num小于ε₂,视为噪点体素并删除。l_abel为1,表示该体素为目标体素;l_abel为0,表示该体素为噪点体素:

label=0,vnum<ε21,vnum>ε2。(6)

2.4 细化算法处理

本文算法采用文献[ 32]中提到的细化策略,对三维空间中的目标体素进行细化处理。设定一组体素空间中的三维模板,用于判定边缘体素,删除模板需要考虑上下左右前后6个方向,“上”方向的6个删除模板如图2所示。

图 2. “上”方向的6个删除模板

Fig. 2. Six delete templates in upper direction

下载图片 查看所有图片

删除模板是以当前目标体素P为中心的3×3×3的体素集合,噪点滤除处理之后,每一个体素都有一个l_abel值,通过判断体素的l_abel值,可以将模板中的体素分为A类和B类两种,A类体素指该体素的l_abel为固定值,标号为0或1;B类体素指l_abel不固定,标号也为0或1。由于体素是三维立体的,每个模板可以顺时针旋转90°、180°和270°,得到的模板仍属于删除模板。为保证各向同性,需考虑上下左右前后6个方向(U、D、W、E、N、S),不同方向的模板可以通过相互翻转得到。若目标体素符合6个方向中的任何一个删除模板,即可进行删除,l_abel值改为0。图2为“上”方向的6个删除模板(不包含旋转之后的模板),其中中心点为需要判别的体素。

根据图2给出的删除模板,可以计算得到其余5个方向上,以及顺时针旋转90°、180°和270°后的所有模板。每个目标体素均与这些三维模板进行空间匹配,若与任意一个模板匹配,则删除该目标体素。当该处理多次循环迭代直至目标体素数量不再减少,即认为体素细化完毕。细化算法流程图如图3所示。

图 3. 细化算法流程图

Fig. 3. Flow chart of thinning algorithm

下载图片 查看所有图片

2.5 生成骨架

细化处理后的目标体素保留了树木的骨架结构,但相互连接关系尚未建立。骨架生成包括三步:

1) 通过计算体素内的所有点的坐标平均值,拟合出每一个体素内所有点云的几何中心点S_P,将其作为树木骨架中的骨架节点,可表示为

Sp(x,y,z)=(∑0nxin,∑0nyin,∑0nzin),(7)

式中:S_P(x,y,z)为体素内部拟合的点的坐标;(x_i,y_i,z_i)为体素内第i个点的坐标;n为体素内所有点云的总个数。

2) 对体素进行26邻域搜索,判断哪些体素相互连接。

3) 确定体素之间的连接关系后,检索连接相邻体素内的骨架节点,构建树木骨架。图4为骨架节点连接示意图。

图 4. 骨架节点连接

Fig. 4. Connection of skeleton nodes

下载图片 查看所有图片

3 分析与讨论

本研究使用获取的银杏树和重瓣榆叶梅树的点云数据对算法进行验证。图5为银杏树与重瓣榆叶梅树的照片及骨架生成结果。图5(a)左侧为银杏树照片,右侧为重瓣榆叶梅照片;图5(b)左侧为由地面三维激光扫描仪获取的银杏树的原始点云,右侧为由地面三维激光扫描仪获取的重瓣榆叶梅树的原始点云;图5(c)左侧为本文方法生成的银杏树骨架,右侧为本文方法生成的重瓣榆叶梅树骨架。树木骨架与点云数据在形态结构上相对一致,说明所提算法有一定的可行性和稳健性。

图 5. 银杏树和重瓣榆叶梅树的骨架生成结果。(a)左侧为银杏树实物图,右侧为重瓣榆叶梅树实物图;(b)左侧为银杏树的点云数据,右侧为重瓣榆叶梅树的点云数据;(c)左侧为银杏树生成骨架,右侧为 重瓣榆叶梅树生成骨架

Fig. 5. Generated skeletons of Ginkgo tree and Amygdalus triloba f. multiplex tree. (a) Picture of Ginkgo on left and picture of Amygdalus triloba f. multiplex on right; (b) point cloud data of Ginkgo on left and point cloud data of Amygdalus triloba f. multiplex on right; (c) generated skeleton of Ginkgo on left and generated skeleton of Amygdalus triloba f. multiplex on right

下载图片 查看所有图片

3.1 数据扫描精度对骨架生成的影响

本研究分别采用扫描角分辨率为0.02°、0.05°和0.10°时得到的银杏树木点云进行树木骨架生成。骨架生成结果如图6所示,角分辨率设置如表2所示。随着扫描精度的降低,得到的点云数据量减少,相应的运算时间也随之减少。采用本文方法处理3个角分辨率下采集到的数据,所得到主枝干的骨架无较大差别,但细小枝干处存在差异。显然,扫描角分辨率会影响点云数据量,从而影响骨架生成中的节点数量和骨架的精细程度,但对树干和主要枝条的骨架生成影响不大。

3.2 体素划分规格对骨架生成的影响

以重瓣榆叶梅树的数据进行参数分析。不同的(M,N,K)参数对骨架生成有一定影响。图7为不同划分规格下生成的重瓣榆叶梅树骨架,参数设置如表3所示。可以看出,划分的体素数量越多,得到的骨架节点越多,在骨架结构上表现有更多的枝干细节。同时,因为原始点云数据量一定,划分规格会出现饱和现象,即随着划分规格的增加,树木形态的保留程度不会出现明显的提高,图7(d)的划分程度比图7(c)的高,但生成骨架的形态无较大差别。

3.3 滤除参数对骨架生成的影响

算法中的滤除参数对树木骨架的生成也有一定影响。本研究根据体素内的点数与体素孤立程度进行判定,设定不同的阈值时,相应保留的目标体素数目也不同,进而生成的树木骨架效果也不尽相同。本研究中ε₁分别取10和20,分析其对骨架生成的影响;树木骨架中应至少存在两个相连体素,所以体素邻域内目标体素小于2的一定是噪点体素,本研究中ε₂分别取2和4,分析其对骨架生成的影响。图8为不同(ε₁, ε₂)参数设置生成的骨架,参数设置如表4所示。可以看出,滤除参数越大,得到的骨架节点越少,运算时间越短。图8(a)和图8(b)为骨架节点孤立程度分析,可以看出,ε₂阈值越大,滤除的点越多,但会在一定程度上破坏树骨架本身的结构;阈值较小时,在滤除孤立点的同时能够较好地保留骨架完整性。图8(c)和图8(d)为骨架节点点数阈值

图 6. 不同角分辨率下点云生成的银杏树骨架。(a) 0.02°;(b) 0.05°;(c) 0.10°

Fig. 6. Skeletons of Ginkgo tree generated by point clouds at various angular resolutions. (a) 0.02°; (b) 0.05°; (c) 0.10°

下载图片 查看所有图片

图 7. 不同划分规格的重瓣榆叶梅树骨架生成结果。(a) (40,40,40);(b) (60,60,60);(c) (80,80,80);(d) (100,100,100)

Fig. 7. Generated skeletons of Amygdalus triloba f. multiplex tree with different partition specifications. (a) (40, 40, 40); (b) (60, 60, 60); (c) (80, 80, 80); (d) (100, 100, 100)

下载图片 查看所有图片

图 8. 不同滤除参数的重瓣榆叶梅树骨架生成结果。(a) (0,4);(b) (0,2);(c) (10,0);(d) (20,0)

Fig. 8. Generated skeletons of Amygdalus triloba f. multiplex tree with different filtering parameters. (a) (0,4); (b) (0,2); (c) (10,0); (d) (20,0)

下载图片 查看所有图片

判定,可以看出,ε₁阈值越大,保留的骨架节点越少,滤除效果越明显。

3.4 算法效率比较

对本文方法(简称IVTA方法)与文献[ 26]中的方法(简称为GSA方法)进行比较,从生成的骨架节点数量与方法运行时间两个方面进行分析。图9所示为两种方法分别处理角分辨率为0.02°的银杏树和重瓣榆叶梅树的点云数据时生成的骨架结果。图9(a)为GSA方法生成的银杏树骨架,图9(b)为GSA方法生成的重瓣榆叶梅树骨架,图9(c)为IVTA方法生成的银杏树骨架,图9(d)为IVTA方

法生成的重瓣榆叶梅树骨架。表5中列出了IVTA方法和GSA方法处理点云数据生成树木骨架获得的骨架节点数和运行时间。显然,IVTA方法相比GSA方法,在完整描述骨架的同时,能够生成更多的骨架节点,保留更多的枝干细节。而且,从表5可知,IVTA方法的运行时间远低于GSA方法。对于实验中的银杏树,GSA方法生成骨架的运行时间约为IVTA方法的31倍;而对于数据量更大的重瓣榆叶梅树,GSA方法的运行时间约为IVTA方法的67倍。显然,相比GSA方法,IVTA方法在运行效率上有很大提高。

图 9. GSA和IVTA方法下银杏树与重瓣榆叶梅树骨架生成结果。(a) GSA方法生成的银杏树骨架;(b) GSA方法生成的 重瓣榆叶梅树骨架;(c) IVTA方法生成的银杏树骨架;(d) IVTA方法生成的重瓣榆叶梅树骨架

Fig. 9. Tree skeletons of Ginkgo tree and Amygdalus triloba f. multiplex tree generated by GSA and IVTA methods. (a) Skeleton of Ginkgo tree generated by GSA method; (b) skeleton of Amygdalus triloba f. multiplex tree generated by GSA method; (c) skeleton of Ginkgo tree generated by IVTA method; (d) skeleton of Amygdalus triloba f. multiplex tree generated by IVTA method

下载图片 查看所有图片

4 结论

提出了一种通过点云体素细化生成树木骨架的方法。通过构建体素空间,判断包含树木点云的目标体素和包含非树木点云的噪点体素;根据体素内的点数与体素邻域内目标体素的数量,对噪点体素进行滤除;利用细化算法对体素进行处理,得到树木骨架的单体素连接结构;最后通过现存的目标体素拟合骨架节点并连接,生成树木骨架。

基于获取的银杏树和重瓣榆叶梅树的点云数据对该方法进行验证。实验设置不同参数验证其对骨架生成的影响。实验结果表明,划分规格越大,单个体素的数量越多,保留的树木细节越多。体素内的点数与体素的孤立分布程度对噪点的滤除都有一定影响。点数滤除阈值比孤立分布阈值的滤除效果略好,对树木的枝干骨架生成效果的帮助更大。实验成功生成了树木骨架,且生成时间都在4 s以内。与GSA方法相比,本文方法能够大大减少生成树木骨架的运算时间,说明本文方法具有一定的有效性和可行性。