1 引言

目前,三维激光扫描已经得到了比较普遍的应用,特别是在构造物的模型重建领域。特征线是描述构造物的主要特征,已经在很多领域得到了应用,如道路提取^[1]、建筑物轮廓线的描述^[2]等。特征点都可以通过特征线得到。特征线在图像领域已经得到了明确的定义,如图像的边界线是两个不同灰度区域之间的界限^[3],通过构建霍夫变换来实现对图像边界线的提取^[4]。在点云领域,特征线分为边界线和折边^[5],边界线主要是非连接区域线,折边主要包含曲率变化区域大的特征点构成的线^[6]和折线^[7]。

点云特征线的提取主要包含以下两个方法^[8]:1) 将点云转为图像,依据图像特征线方法来实现点云特征线的提取;2) 通过分析点云几何特征,如曲率来判断点的归属,从而实现点云特征线的提取。第1种方法需要高质量的图像,第2种方法在计算曲率的过程中耗时比较长。针对此,文献[ 9]利用协方差分析点归属不同特征的概率,但无法实现对目标物突兀区域特征线的判断。文献[ 10-11]利用领域点法向量特征和张量投票方法来实现对点归属特征的判断。文献[ 12]采用泰森多边形协方差法来估计主曲率及法向量方向,进而在点云中判断点是不是属于特征线。文献[ 13]通过主曲率方向来实现对脊线的判断,但需要耗费大量的时间。为了缩短时间,文献[ 14-15]基于区域增长算法采用面状区域间交集来实现特征线的判断,但该方法无法准确地识别小区域和狭长区域面状。随机抽样一致性(RANSAC)算法被提出用于小区域平面的探测,根据相邻平面的交集便可实现特征线的判断^[16],但该方法却无法对边界点进行准确探测。为准确探测边界点,文献[ 5]通过结合RANSAC和邻近点角差方法来实现边界的判断。该方法比较适用于较大平面的交集线,但同时比较尖锐的特征可能会丢失。文献[ 17-19]给出了改进的特征面探测方法,从而实现了敏感特征线的探测,然而,该方法受噪声点影响较大,因为噪声点的梯度变化和相邻平面交集的变化相似。

目前,点云特征线的探测,主要是基于图像特征线提取方法、曲率方法、RANSAC方法及张量投票方法,这些方法没有顾及算法时间效率和点云边界及折边的区别。因此,不能高效地提取点云特征线。本文在顾及时间效率和点云特征线的基础上,分步对边界线和折边线进行提取。通过采集不同类型目标物点云数据,验证了所提方法的有效性。

2 算法模型

考虑到时间效率,基于点云特征线提出一种特征线提取方法,该方法包含边界线和折边的提取。主要步骤如下:

1) 原始点云数据的预处理,包括点云配准及较大噪声的剔除。

2) 利用最邻近点搜索算法,对邻近点进行搜索,利用整体最小二乘(TLS)构建基准面,将邻近点投影到基准面上,构建探测点和投影邻近点之间的向量,通过计算向量方位角G来实现对边界线的探测和提取。

3) 计算邻近点的法向量,并利用k-均值聚类算法^[20]来实现对邻近法向量的聚类,根据最佳的聚类数来实现折边点的判断。

其具体的流程如图1所示。

图 1. 特征线探测流程

Fig. 1. Flow chart of feature line detection

下载图片 查看所有图片

2.1 边界线提取

利用KNNS(k邻近点搜索算法)^[21]来实现对探测点P₀的k邻近点的搜索,根据加权主成分分析算法(WPCA)实现法向量n=(A,B,C)的计算,其中A、B、C为向量分量,如图2所示。

图 2. 基于邻近点的法向量

Fig. 2. Normal vector according to adjacent points

下载图片 查看所有图片

如图3所示,根据法向量和探测点P₀(X₀,Y₀,Z₀)可以直接确定投影基准面(x,y,z),即

A(x-X0)+B(y-Y0)+C(z-Z0)=0,(1)

将邻近点P_j=(X_j,Y_j,Z_j)投影到该基准面上,则投影点p_j=(x_j,y_j,z_j)可以表示为

A(xj-X0)+B(yj-Y0)+C(zj-Z0)=0xj-XjA=yj-YjB=zj-ZjC=λ,(2)

图 3. 投影平面上的法向量方位角。(a)边界上的点;(b)平面上的点

Fig. 3. Azimuthal angle of normal vectors on projection plane. (a) Point at boundary; (b) point in plane

下载图片 查看所有图片

由(1)式可以得到基准面参数

λ=AX0+BY0+CZ0-AXj-BYj-CZjA2+B2+C2,(3)

将(3)式代入(2)式,投影点p_j表示为

xj=Aλ+Xjyj=Bλ+Yjzj=Cλ+Zj,(4)

探测点和投影点之间的向量表示为

pj=xj-X0,yj-Y0,zj-Z0。(5)

在投影面上构建二维坐标系统o'-x'y',其坐标为(x'_j,y'_j),则二维坐标系统中探测点和邻近点方位角为

τj=arctany'jx'j。(6)

对方位角进行排序,得到新的方位角τ'_j,则排序后的相邻方位角差值为

Δτ'j=τ'j-τ'j-1。(7)

根据邻近方位角差值比较的最大值max(Δτ'_j)和阈值Δτ'_thr进行比较便可确定边界线。

1) 如果max(Δτ'_j)>Δτ'_thr,则探测点P₀属于边界点。

2) 如果max(Δτ'_j)<Δτ'_thr,则探测点P₀属于非边界点。

为了确定阈值Δτ'_thr,分别对不同距离和不同邻近点的阈值进行判断,结果如图4所示。由图4(a)可知,该阈值不受距离分辨率的影响,其值在85°左右;而由图4(b)可知,阈值受邻近点的影响较大,并随着邻近点的增加而减少,而根据扫描特点,其每个点的最少的邻近点是8个,而在邻近点为8时,其对应的阈值为90°左右,因此,邻近方位角差值的阈值近似为90°,最终将阈值确定为90°。

图 4. 邻近法向量方位角的最大差值。(a)不同距离分辨率下;(b)不同邻近点下

Fig. 4. Maximum difference of azimuthal angles of adjacent normal vectors. (a) Under different range resolution; (b) under different neighborhood points

下载图片 查看所有图片

2.2 折边提取算法

首先对点云法向量进行计算,再利用k-均值聚类算法对法向量进行聚类,根据最佳的聚类数来判断点的归属情况,其具体步骤如下:

1) 利用最邻近点搜索算法^[21]确定点云中每个探测点的最邻近点。将探测点和任意两个邻近点组成三角形,如图5所示。

图 5. 探测点和邻近点组合成的三角形

Fig. 5. Triangles composed by detection point and nearest neighbor points

下载图片 查看所有图片

2) 利用探测点P₀(x₀,y₀,z₀)和邻近点P_i(x_i,y_i,z_i)确定三角形法向量n_j:

ni-1=xi-1-x0,yi-1-y0,zi-1-z0ni=xi-x0,yi-y0,zi-z0nj=ni-1×ni。(8)

3) 对法向量n_j进行单位化,得到领域单位法向量

ej=njnj。(9)

4) 利用k-均值聚类算法^[20]对领域单位法向量e_j进行聚类,根据最终的最佳聚类数来判断探测点的归属,由于单位法向量存在两个方向,则在探测点是平面的情况下,最佳的聚类数是2。

5) 目前,已有多种指标来判断最佳聚类数,如文献[ 22]提出的CH(Calin-ski-Harabasz)指标,文献^[23]提出的Wint(Weighted inter-intra)指标和文献[ 24]提出的Sil(Silhouette)指标(I_sil)。本文将利用目前已经广泛应用的Sil指标来判断最佳聚类数。

最佳聚类数的基本判断情况如下:

1) 如果I_sil=2,则探测点P_i,属于平面点,不属于折边点,如图6(a)所示。

2) 如果I_sil=4,则探测点P_i,属于折边点,如图6(b)所示。

3) 如果I_sil>4,则探测点P_i,属于曲面点,不属于折边点,如图6(c)所示。

图 6. 探测折边点的聚类结果。(a)平面上点法向量的聚类;(b)折边上点法向量的聚类;(c)曲面上点法向量的聚类

Fig. 6. Cluster results of fold point detection. (a) Cluster of point normal vectors in plane; (b) cluster of point normal vectors on fold; (c) cluster of point normal vectors on curve

下载图片 查看所有图片

3 实验分析

3.1 边界提取

为了单独对该方法边界提取的可靠性进行验证,对图7(a)进行提取,结果如图7(c)所示,并将该提取方法与改进的RANSAC方法进行比对,结果如图7(b)所示。

图 7. 边界的提取。 (a)原始点云;(b)改进的RANSAC方法;(c)所提方法

Fig. 7. Extraction of boundary. (a) Original point cloud; (b) improved RANSAC method; (c) proposed method

下载图片 查看所有图片

由图7(b)可知,改进的RANSAC方法提取的边界含有较多的噪声,造成该现象的原因是邻近点平面拟合容易受点云噪声的影响,从而探测到错误的边界点信息。所提方法没有利用邻近点平面拟合方法,而是利用了邻近点法向量夹角信息,从而避免点云噪声的影响,由图7(c)可知,所提出方法可以准确地提取边界线,且该方法明显优于改进的RANSAC方法。

3.2 折边提取

利用所提算法对某建筑物进行提取,并将提取结果与基于曲率方法^[12]和基于平面拟合方法^[15]进行比对分析,结果如图8所示。

由图8可知,基于平面拟合方法提取效果最好,所提方法和基于平面拟合方法相似,基于曲率方法最差,主要原因是点云存在误差,平面区域曲率变化较大,误将平面认为折边。

为了比较全面地判断所提方法的普适性,对非平面构造的建筑物进行折边的提取,结果如图9所示。

由图9可知,所提方法主要是根据邻近点法向量聚类来判断折边情况,其提取的折边能够比较清晰地概括建筑物的轮廓,而基于平面拟合方法和曲率方法丢失了一部分折边点,因此,在非平面情况下,所提方法在折边提取效果方面要优于其他两种方法,说明所提方法不仅能够适用于平面构造物,也适用于非平面构造物的折边提取。

图 8. 不同方法对平面构造物折边的提取效果。(a)原始点云;(b)基于曲率方法;(c)基于平面方法;(d)所提方法

Fig. 8. Fold lines extracted by different methods under condition of plane construction. (a) Original point cloud; (b) based on curvature method; (c) based on plane method; (d) based on proposed method

下载图片 查看所有图片

图 9. 不同方法对非平面构造物折边的提取效果。(a)原始点云;(b)基于曲率方法;(c)基于平面方法;(d)所提方法

Fig. 9. Fold lines extracted by different methods under condition of non-plane construction. (a) Original point cloud; (b) based on curvature method; (c) based on plane method; (d) based on proposed method

下载图片 查看所有图片

3.3 特征线提取

利用RIEGL-VZ400对某桥梁进行两期扫描,扫描时间分别为2017年8月及2017年11月,两期扫描位置及参数设置相同,距离分辨率设置为1 cm,点云平均密度为4000 m^-2。扫描距离约为80 m,采集的原始点云数据如图10所示。利用所提方法得到两期桥梁点云的特征线如图11所示。由图可知,所提方法提取的特征线保留了桥梁的所有特征,特别是对于桥梁缆线的提取,没有造成缆线的截断和损失,同时桥梁的主体边缘特征也得到了保留,从而满足了特征线提取的要求,同时,由图11(a)和(b)对比可知,两期桥梁点云特征提取效果基本相同,即所提方法适用于不同时间复杂度扫描的点云数据。

图 10. 桥梁点云

Fig. 10. Point cloud of bridge

下载图片 查看所有图片

图 11. 桥梁特征线的提取。(a)第一期点云提取的特征线; (b)第二期点云提取的特征线

Fig. 11. Extracted feature lines of bridge. (a) Feature lines extracted from first phase point cloud; (b) feature lines extracted from second phase point cloud

下载图片 查看所有图片

4 结论

通过分析基准面上探测点和领域点法向量方位角的分布情况来判断目标物点云的边界归属情况。通过计算邻近点法向量聚类情况来判断目标物点云的折边归属情况。利用三维激光扫描仪采集的平面和非平面构造物来对其进行验证,所提出的边界提取算法要优于改进的RANSAC方法。对于平面构造物,所提出的折边提取算法和基于平面拟合算法对于折边提取效果相似,且都优于基于曲率方法的折边提取;对于非平面构造物,所提方法提取的折边效果明显优于其他两种方法。最终利用所提方法对不同时期扫描的桥梁点云特征线进行了提取,用最少的点比较清晰地再现了桥梁轮廓,证实所提方法适用于不同时期扫描的点云。