1 引言

三维激光扫描技术是测绘领域继全球定位系统(GPS)技术之后又一项测绘新技术,其具有扫描精度高、速度快、易携带和非接触式采集操作等优点,已成为空间数据获取的重要技术手段,可以快速地从模型表面获取数十万甚至数百万的测量数据,并可准确的反映空间位置,实现了信息的可视化,原始点云在建筑的三维模型化设计^[1]以及变形监测^[2-3]中已经发挥了不可忽视的作用。但是因为数据量大和处理效率低等问题,而催生出三维模型的特征线提取技术。建筑物的特征线是建筑物形状描述、三维重建的重要参数,要建立一个三维建筑物模型首先要做特征线的提取。目前,特征线的提取已广泛用于逆向工程、曲面重建、数据分割、模型简化、曲线匹配和拼接等几何处理领域^[4-5],通过提取树木特征信息确定其位置分布的方法甚至已经被应用于森林生态、管理林区等工作^[6]。

对于三维点云模型,其特征曲线提取过程就是对模型上的数据点进行分析计算,找出其中的特征点,并连接构成光滑特征曲线。针对点云特征线的提取,许多学者做了大量的研究并且提出了相关的算法:1) 基于曲率值、法向量以及投影等几何特征的方法。钱锦峰等^[7]提出了一种利用点云数据小邻域内点的相邻关系检测边界特征点的算法,该算法通过最小生成树算法对点云数据进行处理,构造目标点邻域的最小二乘拟合平面,计算出投影点的角度标准差阈值,从而确定特征点并排序;Daniels等^[8]对每个点的邻域内使用稳健移动最小二乘法(RMLS)进行曲面拟合,并计算当前点在曲面上的投影残差值,将投影残差值较大的点定义为潜在特征点;陈朋等^[9]提出了一种建筑物点云特征线提取方法。2) 基于主成分分析的方法。Wang等^[10]提出了一种体素维度分析法,该方法以主成分分析(PCA)降维为基础,计算体素点云协方差矩阵,对协方差矩阵进行分解得到特征值和特征向量,根据特征值大小判断出显著特征向量,从而实现特征点的保留;Demarsin等^[11]首先用PCA法求出点云数据的法向量,再对其变化大的点采用区域生长算法,并聚类分割数据,从而得到特征线。史皓良等^[12]提出一种双阈值点云特征信息提取算法,采用PCA法和局部二次曲面拟合法对点云模型进行微分几何信息估算,得到k 邻域内采样点平均法矢夹角和平均曲率的特征权值,利用双阈值检测方法对散乱点云的特征信息进行提取。3) 基于映像和模型构建的特征线提取方法。王大莹等^[13]运用数学形态学中的灰度形态学,将光探测和测距(LIDAR)数据通过C++平台处理得到深度影像,并在Matlab平台下进行特征提取,从而得到建筑物的边缘信息;邹万红等^[14]提出了一种新的点云数据骨架提取方法,通过对点云数据的空间层次剖分后建立其简化模型,根据离散Morse理论提取主要的特征点,得到模型的初步骨架,再将初步骨架内推至模型内部,对内推后的骨架光顺及聚类后形成最终骨架。经过研究学习以及对已有算法的实现得出已有算法存在算法复杂、容易受到实验中误差的影响、无法有效地剔除平面上的特征点和提取折边的特征点从而导致特征线不完整以及无法准确提取出不规则点云特征线等问题。

针对以往算法的不足,从聚类的角度出发,本文提出了一种基于高斯映射和K均值(K-means)聚类判断特征点的特征提取方法,该方法首先对大量原始点云数据进行简化,可以弥补之前算法的特征点提取数量大、噪声多、折边丢失、特征线不完整、仅对规则图形提取效果好等缺陷,较之已有算法有以下优点:1) 简化数据之后大大降低了数据处理量,减少了噪声点,提高了处理效率;2) 将邻域三角形集合法向量投影在高斯球上,根据投影点在曲面上的分布规律,进行K-means聚类,根据聚类有效性指标筛选出最佳聚类个数,从而得到特征线,实验证明,该算法不仅可以完整地提取规则点云数据的特征线,还可以高效精准的提取不规则点云数据的特征线。

2 算法概述

所提算法对于特征线的提取首先用k近邻搜索寻找点云数据每个目标点的k个临近点,然后将目标点与其任意两个邻域点排列组合构造成三角形,并求出其单位法向量,再将其进行高斯映射,最后利用聚类的方法达到特征点判别的目的。算法流程如图1所示。

图 1. 算法流程

Fig. 1. Flow chart of algorithm

下载图片 查看所有图片

2.1 k近邻算法

k近邻算法是一种机器学习算法,对于点云数据集p_i=(x_i,y_i,z_i)∈R³,i=1,2,…,n,其中R为实数集;n为点云集合中点云的个数;i为正整数集合,表示点云顺序;点p_i的k邻域是指集中在点p_i周围,并与p_i距离最短的k个数据点,按与当前采样点的距离d从小到大排列的点云集合 pjk=(x_j,y_j,z_j),j=1,2,…,k,k为正整数集,表示邻域点个数;j为正整数集合,表示k个邻域点的顺序。利用距离度量来进行临近点的搜索,如欧氏距离、曼哈顿距离、闵科夫斯基距离,本文选用最常用的欧氏距离作为距离度量,表示为

d=(xi-xj)2+(yi-yj)2+(zi-zj)2,(xi,yi,zi)∈p,(xj,yj,zj)∈pk,(1)

式中:p为点云集合;p^k为某一点的k邻域点集合。

k近邻搜索实现的步骤中最重要的是k值的选择与搜索过程。

1) k值过小容易受到噪点的影响,导致特征点判断不准确,k值过大,虽然会减小噪点导致的误差,但是会增加计算量,降低运算速度,通过对多个实验点云模型计算,确定邻域大小k的取值一般取为8~32。

2) k近邻算法搜索的过程是通过构造K维度(Kd)树实现的,Kd-树中每个节点即为一个K维的点,变量K表示空间维度数。每个非叶子节点可以想象为一个分割超平面,利用垂直于坐标轴的超平面将空间分为两部分,这样递归的从根节点不停地划分,直到没有实例为止,如图2所示。

图 2. Kd-树的递归过程

Fig. 2. Recursive process of Kd-tree

下载图片 查看所有图片

经典的构造Kd-树的规则如下:

1) 随着树的深度增加,对于数据点集及其所在的空间,计算每个维度上的方差,数据方差最大时表明沿该坐标轴方向上的数据点分散的比较开,在此进行的数据分割可以获得最好的分辨率;循环的选取坐标轴,作为分割超平面的法向量。对于三维树来说,根节点选取x轴,根节点的孩子选取y轴,根节点的孙子选取z轴,根节点的曾孙子选取x轴,这样循环下去;

2) 每次均为所有对应实例的中位数作为切分点,切分点作为父节点,左右两侧为划分的左右两子树。如此循环往复直至空间中只包含一个数据点。通过上述搜索过程,对于三维点云模型可表示为p_i=(x_i,y_i,z_i)∈R³,i=1,2,…,n,寻找到每个采样点p_i的k邻域为 pjk=(x_j,y_j,z_j)∈R³,j=1,2,…,k,p_i的k个近邻点按距离值d由小到大依次排序。

2.2 高斯映射理论

高斯映射已广泛应用于曲面识别、点云拟合、点云分割等领域。高斯映射的定义为:对一个分段光滑的空间三维曲面∑,计算其上所有点的单位法矢量n,并将法矢的起点平移到单位球S²的球心上,则法矢量的终点落在单位球面上,形成投影点的集合称为该曲面∑的高斯图,曲面上的点向球面上映射的过程,称为高斯映射G_p:∑→S²,该单位球称为高斯球,如图3所示,其中O为原点。

图 3. 曲面∑的高斯映射

Fig. 3. Gauss map of surface ∑

下载图片 查看所有图片

使用上述原理,以点云数据中取当前采样点p_i为顶点,与其邻域内的k个近邻点构成k(k-1)/2个三角形,设这些三角形的集合为T={Δ_mn=(p_i,p_m,p_n)∈R³,1≤m,n≤k},m、n表示邻域任意两个临近点编号;Δ_mn表示以p_i、p_m、p_n3点为顶点所构成三角形;T为三角形的集合;每个三角形Δ_mn的法向量为两条边的叉乘即N_mn=p_ip_m×p_ip_n,然后将这些法向量单位化,p_i是单位高斯球的圆心,法向量的起点即单位球原点,则点p_i的k邻域内的高斯映射点的定义为G_p…T→S²,某一目标点的邻域三角形法向量集合在单位高斯球的投影如图4所示。

Δmn→xmn:pi+Nmn‖Nmn‖。(2)

图 4. 某一目标点的邻域法向量投影

Fig. 4. Neighborhood normal vector projection of certain target point

下载图片 查看所有图片

上述过程虽然得到了法向量,但是法向量会出现反向簇,需要将法向量的方向调整一致,才不会影响后续的处理,通常点云都有主扫描面,以该面上z坐标最大的点为基准点,计算每一点与该点法向量的夹角,如果夹角>90°,则将该点的法向量朝向改为反方向。对于夹角<90°的点不用调整,直到所有的点都调整完毕。x_mn即为p_i的高斯映射点,不同的曲面在高斯球上的映射点有不同的分布规律:1) 零维分布,对于平面上的点,其法向量n为常数,投影点重合于一点,如图5(a)所示;2) 一维分布,对于圆锥、圆柱面和一般拉伸面,任何一个点p_i的法向量n和轴线l成一个固定夹角α,投影点形成一个二次圆弧曲线,即高斯映射像呈一维分布,如图5(b)所示;3) 二维分布,对于球面、圆环面、一般直纹面上的投影点在球面上均匀分布,形成面形分布点,如图5(c)所示。

图 5. 高斯映射点的分布规律。(a)点型聚类;(b)线型聚类;(c)面型聚类

Fig. 5. Distribution law of Guassian map points. (a) Point clustering; (b) linear clustering; (c) surface clustering

下载图片 查看所有图片

2.3 高斯映射点的K-means聚类

对点云数据中每个采样点p_i邻域三角形的单位法向量进行单位高斯球映射之后,运用K-means对高斯映射点进行聚类。K-means聚类是应用十分广泛的聚类算法之一,其计算过程十分直观:在每一个采样点p_i高斯映射点的集合中,选取k'个元素作为k'个簇各自的中心(此处为了与上文变量(邻域点个数k和维度数K)做出区分,选用变量k'作为聚类中心个数)。再分别计算剩下的元素到k'个簇中心的相异度为

dii'=xi-xi'+yi-yi'+zi-zi',(xi,yi,zi)∈p,(xi',yi',zi')∈pk',(3)

式中:d_ii'为两点间的曼哈顿距离;i'为整数集合,表示k'个聚类中心的顺序编号;p^k'为k'个聚类中心组成的集合。将这些元素分别划归到相异度最低,即曼哈顿距离d_ii'最小的簇;根据聚类结果,重新计算簇中所有元素各自维度的算术平均数,作为k'个簇新的中心;然后将集合中全部元素按照新的中心重新聚类;重复上一步,直到聚类结果不再变化,将结果输出。

其中,确定最佳聚类数的聚类有效性指标是衡量聚类算法产生的聚类结果是否达到最优的标准,该指标将最优聚类结果所对应的聚类数作为最佳聚类数。目前已有的聚类有效性指标主要有CH (Calin-ski-Harabasz)^[15]指标、Wint (Weighted inter-intra)^[16]指标和Sil (Silhouette)^[17]指标等。在对点云数据的多次对比以及重复性实验的过程中得出,与其他几个指标相比,Sil指标简单易用且具有良好的评价能力^[18],轮廓系数的定义为:1) 计算样本ii同簇到其他样本的平均距离a_ii,其值越小,说明样本ii越应该聚类到该簇,因此将a_ii称为样本ii的簇内相异度。某簇中所有样本的a_ii均值称为该簇的簇相异度。2) 计算样本ii到其他某簇的所有样本的平均距离b_iij',j'为其他簇顺序编号,j'=1,2,…,k';称为样本ii与其他簇的相异度。定义为样本ii的簇间相异度b_ii=min{b_ii1,b_ii2,…,b_iik'},b_ii越大,说明样本i越不属于其他簇。3) 根据样本ii的簇内相异度a_i和簇间相异度b_i,定义样本ii的轮廓系数为

sii=(bii-aii)maxaii,bii=1-aiibii,aii<bii0,aii=biibiiaii-1,aii>bii。(4)

当s_ii接近1时,说明样本ii聚类合理;当s_ii接近-1时,说明样本ii更应该分类到另外的簇;当s_ii近似为0,则说明样本ii在两个簇的边界上。求出所有样本的轮廓系数后再求平均值就得到了平均轮廓系数。其取值范围为[-1,1],且簇内样本的距离越近,簇间样本距离越远,平均轮廓系数越大,聚类效果越好。Sil指标反映了聚类结构的簇内紧密性和簇间分离性,那么,平均轮廓系数最大的k'便是最佳聚类数,接下来根据最佳聚类数进行特征点的判断。

2.4 特征点的判断

通过以上聚类算法,对当前采样点p_i邻域的高斯映射点进行聚类,得出每个目标点的最佳聚类数,根据上述不同曲面高斯映射点的分布规律可知:若以p_i为中心的单位高斯球上的仅保留一个簇,其潜在区域则是一个平面,高斯映射点为零维分布,则认为不是特征点;若保留了2~4个簇,其潜在区域则为圆锥、圆柱面或一般拉伸面,高斯映射点呈一维分布,最多可以聚类成4个簇,则该点是特征点;若是折边(两个面或者多个面相交),则交线上样本点的高斯映点,将在高斯球上聚类成两个或者多个簇,则也为特征点;若高斯球上没有保留簇或者大于4个簇,则认为该点不是特征点。

3 实验结果及其分析

为了验证所提算法的有效性,采用三维激光扫描仪(RIEGLVZ-400,Riegl公司,奥地利)进行数据采集,对武汉长江二桥进行扫描,截取了索塔的立柱部分,采用Cloud Compare的“edit-octreeresample”命令对其进行均匀抽稀(Cloud Compare 是一个开源点云处理软件,主要用于多组点云的对比显示、拼接和分割等),将原始点云由400万简化到约22.5万,既可以保持原始点云的特征又可以提高数据处理效率,得到如图6(a)所示的效果图。

为了证明所提算法提取特征线的优越性,引入了对比实验,图6为对长江二桥立柱使用不同方法进行实验的结果,其中图6(b)为所提算法近邻点个数k=14,最佳聚类数2~4时,得到的特征线;图6(c)为文献[ 9]中提出的点到面的距离均值法k=14, D>0.04时,D表示目标点到平面的距离,得到的特征线;图6(d)为文献[ 10]中提出的特征值法,k=14时筛选出来的特征点。对比图6(b)~(d)可以看出,文献[ 9]算法平面上残留的噪声点较多,文献[ 10]算法虽然噪声少,但同时也过滤掉了大量的特征点,而所提算法可以清晰完整地提取出大桥立柱的特征轮廓线,且噪声少。

图 6. 不同算法对规则点云特征线提取效果对比。(a)简化后点云;(b)所提算法;(c)文献[ 9]中的算法;(d)文献[ 10]中的算法

Fig. 6. Extraction effect comparison of feature lines of regular point clouds by different algorithms. (a) Simplified point cloud; (b) proposed algorithm; (c) method in Ref. [9]; (d) method in Ref. [10]

下载图片 查看所有图片

由图6的对比实验可以看出,所提算法可以精准地提取出规则点云的特征线,高效地剔除非特征点,保留了完整清晰的特征轮廓线。为了验证所提算法对不规则点云的提取效果,选取某广场雕塑进行对比实验。

图7(a)为用Cloud Compare简化后的原始点云;图7(b)为所提算法实验结果;图7(c)为文献[ 9]算法的结果;图7(d)为文献[ 10]算法的结果。文献[ 9]算法对规则部分提取效果较好,不规则部分只保留少部分散乱点云, 文献[ 10]算法的效果优于文献[ 9]算法,但是不规则部分仍过滤掉许多特征点。通过比较可见,3种算法在提取雕塑上部不规则部分的效果具有明显区别,而对于不规则部分, 所提算法的效果明显优于另外两种算法,既可以完整清晰地提取出不规则曲面的特征线,又可以保留细部的特征点。

表1为3种算法对点云模型特征点和运行时间的统计对比,从运行速率和特征点保留两方面,对不同算法进行了定量地分析。针对特征点保留,结合表1以及图7(b)~(d)可以得出以下结论:所提算法对特征点的保留率为14.21%,完整清晰地保留了原始的特征点;文献[ 9]算法的特征点保留率为12.64%,但是由图7(c)可以看出,这些保留的点中包含有大量的非特征点,提取效果并不是很好;文献[ 10]算法与所提算法保留的特征点数目相差无几,保留率为13.11%,但是结合图7(d)可以看出,特征值法虽然具有良好的特征点提取以及非特征点去除性,但是会造成特征线的部分缺失。从时间角度分析,所提算法运行效率高。综上所述,所提算法可以清晰完整地保留点云数据的原始特征点,准确地去除非特征点,而且运行效率相对较高,对于海量点云数据的运算能够节省大量的运行时间。

图 7. 不同算法对不规则点云特征线提取效果对比。(a)简化后点云;(b)所提算法; (c)文献[ 9]中的算法;(d)文献[ 10]中的算法

Fig. 7. Extraction effect comparison of feature lines of irregular point clouds by different algorithms. (a) Simplified point cloud; (b) proposed algorithm; (c) method in Ref. [9]; (d) method in Ref. [10]

下载图片 查看所有图片

4 结论

提出了基于高斯映射的三维激光特征线提取方法,经实验证明,与已有算法相比,该方法可以清晰完整地提取出规则以及不规则点云模型的特征线,噪声少,算法运行时间短且精确、高效,还能够快速高效地获取物体或者建筑的特征模型,对于逆向工程、建筑物三维重建工程等,可以大大提高工作效率和精确度。该算法不仅可以很好地提取不规则点云的特征模型,还对安全领域变形监测(滑坡、建筑变形监测等)具有相当的意义。