1 引言

激光三维成像是指系统主动发射激光照射目标,结合激光飞行时间或相位得到目标的三维点云。其具有距离远、精度高、不受光照影响等特点,可用于目标识别、逆向工程及相对导航等,具有很高的**和民用应用价值。点云分割依据几何结构特征将目标点云分割为若干个有意义的部分,是点云处理的重要步骤,单个点云分割后可用于区域特征识别,两个部分重叠的点云分割后对应的重叠区域可用于点云配准。然而,受测距精度和后向散射等因素影响,激光三维成像得到的目标点云存在严重的高斯噪声和离群点,且目标复杂多样,无法预定义分割规则,这给激光三维成像的点云分割造成了一定的困难。

国内外学者对点云分割进行了大量研究,点云分割方式主要有:在空域依据点的局部几何特征连续性或相似性聚类,以及在频域依据点的局部几何特征相似性进行谱聚类。

Koschan^[1]采用区域生长方法分割点云,该方法对生长种子点的预设选取参数敏感;Rabbani等^[2]依据点的法向量和曲率等局部几何特征连续性分割点云,其对含噪声点云表现不好。依据点的局部几何特征连续性的分割方法对噪声和预设参数较敏感,对激光点云的适用性并不好。

Shlafman等^[3]较早将k-means聚类应用于点云分割;任同群等^[4]采用聚类方法将点云中非平坦区域作为匹配特征,拟合后实现点云配准;袁小翠等^[5]利用k-means聚类将点云分割为平坦区域和非平坦区域,实现了点云精简的密度自适应;卢维欣等^[6]将非地面点按照距离准则聚类分割,并对每个点云簇进行特征分析,实现了特征识别;王肖等^[7]在相对导航应用中采用聚类方法对象化分割目标点云,提高了效率和准确率;周文振等^[8]利用k-means聚类对迭代最近点配准算法(ICP)进行改进,完成了简单室内激光点云的配准。采用中心聚类的点云分割方法具有一定的抗噪性,在处理良分离目标应用中表现较好,但是依然存在对目标的几何结构解析能力较差的问题。

Golovinskiy等^[9]借助谱聚类中的Min-cut准则实现了点云前景和背景的自适应分割;马腾等^[10]指出点在谱空间内能够表现出更好的类簇性,并使用归一化的几何矩和角度距离构造邻接矩阵,实现了谱聚类点云分割;韩丽等^[11]提取点的测地线等多特征构造邻接矩阵,并通过拉普拉斯矩阵特征值本征间隙自适应确定聚类数目,完成了点云的自适应分割。可见,谱聚类分割方法将点云特征转换到谱空间,保持了聚类算法一贯较好的抗噪优势,且从全局角度考虑分割问题,能更好地解析目标点云几何结构。

但目前的谱聚类分割方法大多依靠邻域点提取点的局部几何特征,在特征提取及分割数目自适应方面仍有改进空间。本文首先采用混合概率主成分分析(MPPCA)法提取点的局部几何特征构造邻接矩阵,然后利用N-cut实现点云特征降维嵌入,最后利用类间类内划分(BWP)算法进行自适应分割。

2 本文算法

点云的谱聚类分割是将点云Xx1,x2,…,xN(N为点云中点的数量)按照局部几何特征相似性分割为k类,同一类点的特征具有较大的相似性,而不同类点的特征具有较大的相异性。谱聚类分割步骤可分为:1)定义点的局部几何特征描述符,构造点云的邻接矩阵;2)利用拉普拉斯变化将点云的邻接矩阵映射到谱空间,降维得到描述每个点分割特征的低维向量;3)依据每个点的分割特征进行聚类,得到分割结果。

2.1 点的局部几何特征提取

点云谱聚类分割需要首先定义点的局部几何特征描述符,然后计算点云中所有点的局部几何特征之间的相似关系,构造邻接矩阵W(i,j)(i,j=1,2,…,N)表征点x_i、x_j之间的相似性,为了避免过拟合,一般认为点x_i仅与其最近邻点k_NN有关,将W(i,j)定义为

W(i,j)=pij, xi,xjbeingnearestneighbortoeachother0, otherwise,(1)

式中p_ij为点x_i、x_j之间的局部几何特征相似性参数,一般情况下可由其邻域点局部采样协方差矩阵进行奇异值分解得到,或者用测地线及法向量等描述。但在点云分割的应用中,空间距离较近的点x_i、x_j对应邻域重合率很大,处于两个分割边界的点与周围相邻区域存在很大的相似性,极易出现过分割或欠分割;并且该方法容易受到邻域中噪声和离群点干扰,而对点云进行滤波处理后再分割则会对目标的细节造成一定影响^[12],因此该特征提取方法并不完全适用于激光点云分割。

考虑到成像目标表面由平面和曲面构成,目标点云为分布于三维欧氏空间的点簇,从流形学习角度可以视点云为分布于三维空间的线性与非线性混合流形,不同的几何特征对应不同的流形结构,处于相同流形上的点几何特征相似,处于不同流形上的点几何特征不同。由此可以不依赖邻域点而完成点的局部几何特征描述,避免了基于邻域点的几何特征提取方法对噪声和参数的敏感性。

本文采用常用的流形学习方法——MPPCA^[13]分析点的局部几何结构特征,主要基于以下几个原因:1)非线性流形在局部可以被多个线性流形组合很好地近似表示^[14],因此其具备分析描述复杂几何结构点云的能力;2)点所在流形主子空间的主角度相似性能够描述点的局部几何特征相似性;3)MPPCA由主成分分析(PCA)法改进得到,具备PCA克服噪声影响提取主成分的能力,通过全局迭代最优解避免了PCA对邻域点选取的依赖,因此,MPPCA具备良好的抗噪性和对预设参数的稳定性。借助MPPCA将处于不同流形的点分离,利用点所在流形的几何特征描述点的几何特征,处于相似流形的点之间的p_ij较大,而处于不同流形的点之间的p_ij较小,这样可以得到符合目标几何结构特征的邻接矩阵。

MPPCA的主要思路是训练M个线性主成分分析器对每个点的局部几何结构进行估计,用最大期望(EM)算法得到最优描述点云分布的M个线性主成分分析器和每个点x_i对应的概率密度分布。假设在第m个线性主成分分析器中,每个点x_i对应于一个降维后的点y_i,则有

xi=Vmyi+μm+εm,(2)

式中V_m为第m个分析器的主子空间,μ_m为分布于第m个分析器的所有点均值,ε_m为噪声。利用高斯混合模型将y_i和ε_m估计为y_i~N(0,I)和ε_m~N(0, σm2I), σm2为噪声标量。结合(2)式和贝叶斯公式,可得x_i关于m和y_i的条件分布为

p(xi|yi,m)=(2πσm2)-3/2exp-12σm2‖xi-Vm-μm‖,(3)

y_i关于m的条件分布为

p(yi|m)=(2π)-1exp(-12yiTyi)。(4)

由(3)式和(4)式可得

p(xi|m)=(2π)-3/2σm2I+VmVTm-1/2exp-12xi-μmTσm2I+VmVTm-1xi-μm。(5)

构造目标函数ζ,即第m个主成分分析器上所有点x_i的对数最大似然估计,得到

ζ=∑i=1Nln∑m=1Mπmp(xim,(6)

式中π_m为各个主成分分析器的混合比例,π_m≥0且 ∑m=1Mπ_m=1。利用EM算法最大化目标函数ζ,即可得到各个主成分分析器的参数。EM算法主要为E步和M步两个部分的交替迭代,令第t次迭代的参数为π_m、 Vm,μm,σm2,第t+1次的估计值为 π~m、 V~m,μ~m,σ~m2,则有

1) E步

π~m=∑i=1Np(m|xi)/N,(7)μ~m=∑i=1Nxip(m|xi)/∑i=1Np(m|xi),(8)

式中m对于x_i的后验概率分布为

p(m|xi)=[πmp(xi|m)]/p(xi);(9)

2) M步

V~m=SmVm(σm2I+Tm-1VTmSmVm)-1,(10)σ~m2=trSm-SmVmTm-1(V~m)T/2,(11)

其中

Sm=∑i=1Np(m|xi)(xi-μ~m)(xi-μ~m)TNπ~m,(12)Tm=σm2I+VTmVm。(13)

σm2的初始化参数为高斯分布中的随机数,π_m的初始化参数为均匀分布,μ_m的初始化参数由数次k-means迭代得到。当两次迭代之间的参数变化小于预设值时,认为算法收敛,此时得到的参数即为M个主成分分析器的参数以及x_i对应的概率密度分布。将概率最大的主成分分析器作为点的所在流形,则点x_i,x_j之间的相似关系为^[15]

pij=∏svd(VmxiTVmxj)。(14)

至此,由(14)式结合(1)式即可得到点云中各点基于MPPCA的邻接矩阵W。

2.2 特征的谱空间低维嵌入

得到邻接矩阵W后,采用谱聚类方法进行切割,以子类内部相似且子类之间距离较远为切割准则。作为一种基于图论的聚类方法,谱聚类可以构造一个全局评价函数,并得到分割最优解,具体方法有Min-cut、Ratio-cut及N-cut等,其中,N-cut最常用且表现最好,其评价函数为^[16]

Ncut(A1,A2,…,Ak)=∑l=1k∑xi∈Al,xj∈A-lW(i,j)/∑xi∈Al,xj∈XW(i,j),(15)

式中A_l(l=1,2,…,k)为第l个子类, A-l为A_l的补集。该评价函数中的惩罚因子包含了对子类数目的约束,避免孤立点形成子类,所以能使分割结果中子类数目较平均,且能够结合非规范化拉普拉斯矩阵L,将NP-complete问题松弛化后转化为P问题。L定义为L=D-W,其中,对角矩阵D为待切割图的度矩阵,对角元素定义为d_i=∑j=1NW(i,j)。根据拉普拉斯矩阵的性质可知,最小化N-cut值的问题可以转化为最小化L矩阵的特征值λ的问题,即取L矩阵的前几个最小特征值。从第二小特征值开始,每个特征值对应的特征向量均可表示点的分类关系。按照对应特征值由小到大从第二个特征向量开始,取q个特征向量,q值由矩阵最大本征间隙决定,将q个N维特征向量构成N×q矩阵,矩阵中第i行的q维向量即对应点x_i的分割特征描述符。

谱聚类把点云转换到谱空间,并在谱空间对点云进行分析,然后将分割特征描述符嵌入低维向量,相较于直接对矩阵W进行聚类分割,谱空间能够更好地表达点云的类别特征,得到更好的聚类结果。

2.3 自适应分割

点云分割过程期望实现无监督分割,即根据点云的几何结构自适应分为k类。为此,结合BWP算法实现点云分割中分类数k的自适应选择^[17]。BWP算法根据子类距离构建对聚类结果的评价函数,并对一定范围内的k值聚类结果进行评价,将评价结果最优的k值作为自适应分割的最优k值。具体表达式为

BWP(l,i)=[b(l,i)-w(l,i)]/[b(l,i)+w(l,i)],(16)

式中b(l,i)和w(l,i)分别为第l个子类中点x_i的最小类间距离和所有点的平均距离,两者定义分别为

b(l,i)=min1≤g≤k,g≠l∑j=1Ngxj-xi2/Ng,(17)w(l,i)=∑p=1,p≠iNlxi-xp2/(Nl-1),(18)

式中点x_p为第l个子类中不同于x_i的点,N_l为第l个子类中点的数量,点x_j为第g个子类中的点,N_g为第g个子类中点的数量。BWP算法能够按照子类内部相似且子类之间距离较远的标准评价聚类结果,BWP值越大说明聚类效果越好。在合适的评价指标之下,将分割的子类数k限定在[k_min,k_max],对每一个k值,计算当前点云中所有点的BWP均值

kBWP_avg=∑l=1k∑i=1NlBWP(l,i)/N,(19)

在k的取值范围内,满足k_{BWP_avg}取最大值的k值即为最佳聚类参数,对应的聚类结果即为点云的最佳分割结果。

3 实验与分析

3.1 仿真点集测试与分析

为了体现本文算法在全局角度分析目标几何结构的能力,分别利用二维空间中两条螺旋线合成的二维点集和三维空间中相交的两条直线及一个平面合成的三维点集对算法进行测试,结果如图1所示。可以看出,本文算法没有在局部曲率变化最大的地方切割点集,而是考虑了几何结构的连续性,从全局结构进行了切割,这一特征使得该算法能够更好地体现目标的几何结构,而不是单纯地依靠曲率变化等局部特征进行分割,该切割方式更符合对目标的结构认知。

图 1. 仿真点集的分割结果。(a)两条曲线构成的螺旋线;(b)相交的直线与平面

Fig. 1. Segmentation results of simulated point sets. (a) Spirals made by two curves; (b) crossed lines and plane

下载图片 查看所有图片

3.2 算法主要参数分析

需要预设的主要参数为主成分分析器数量M和近邻点数k_NN。预设参数的调节是点云分割的一大难题,根据目标点云调节参数会带来巨大的工作量,自适应点云分割期望在较宽预设参数范围内得出理想的分割结果,本文参数分析部分主要考察算法对预设参数的稳定性。实验采用了三组激光点云Rabbit、Horse和Armadillo^[18-20]对算法进行测试。图2所示的分割结果为实验中出现最多的结果,即本文算法收敛于全局最优解时得到的分割结果,可以看出,该结果较好地符合目标几何结构。但是不同的参数会影响得到此结果的概率。当参数不合适时,自适应分割结果会收敛于存在严重分割错误的局部最优解。以图2所示的结果为基准,由分割结果一致性评价算法对参数的稳定性,其中分割结果一致性为算法在不同参数条件下收敛于全局最优解的概率。

图 2. 符合目标点云几何结构的分割结果。(a) Rabbit; (b) Horse; (c) Armadillo

Fig. 2. Segmentations of point clouds based on target geometry. (a) Rabbit; (b) Horse; (c) Armadillo

下载图片 查看所有图片

实验过程如下:首先选出一组分割结果一致性较好的M和k_NN,固定k_NN,将M以5为步长在5~100之间变化,若分割结果与图2中分割数目相同且各部分点与点对应比例均在95%以上,则认为该分割结果与图2结果一致,测试结果记为1,反之则记为0。由于MPPCA中的EM学习是迭代最优化的,因此完全相同参数条件下分割结果也有可能不同,故对每个M值进行100次重复实验,统计测试结果为1的百分比,将其作为当前M值的一致性评价标准。k_NN的测试按相似的方式进行。

1) 分割结果一致性与M的关系

M值与目标的复杂程度有关,需要对目标进行具体的调试才能得到最佳值。但是从图3可以看出,在M值偏小或偏大不严重的情况下,本文算法均能得到较理想的结果。其中Rabbit由于结构明显,在合适的M取值情况下,分割结果的一致性达到了100%,Horse和Armadillo的分割结果一致性也均达到了80%以上。但M的取值和分割结果一致性并不成正相关,这是因为适当数量的主成分分析器已经可以较好地描述点云结构,M过大会造成点云中的点在众多主成分分析器上分布过于离散,点与点之间的几何结构关系提取不准确。

同时,M值与EM算法中的混合概率模型复杂程度有关,其对程序运行时间会造成影响,测试结果如图4所示。其中Rabbit、Horse、Armadillo图中三维点数量分别为3921、5091、5370。综上,推荐M值大于30即可。

图 3. 点云分割结果一致性与M之间的关系

Fig. 3. Relationship between consistency of segmentation and M

下载图片 查看所有图片

图 4. 程序运行耗时与M之间的关系

Fig. 4. Relationship between time consuming and M

下载图片 查看所有图片

2) 分割结果一致性与k_NN的关系

点的局部几何结构由M个主成分分析器的混合概率描述,与k_NN的取值没有直接关系。但是为了避免欠分割或将局部几何结构相似而不相邻的点分割到一个部分,将Wi,j在x_i、x_j不互为k_NN点时赋值为0。因此,k_NN值与分割结果一致性也存在一定的相关性,测试结果如图5所示。可见该算法对k_NN的取值也存在较大的宽容度。

图 5. 点云分割结果一致性与kNN之间的关系

Fig. 5. Relationship between consistency of segmentation and value of kNN

下载图片 查看所有图片

综上可知,算法在较宽的参数范围内能得到符合目标几何特征的全局最优解。虽然不适当的参数会影响得到正确结果的概率,但仍然可以通过多次运行得到最大概率下点云的正确分割。

3.3 点云测试与分析

受测距精度和后向散射等因素影响,激光点云一般存在高斯噪声严重且离群点较多的情况。为了测试算法在激光点云分割中的应用效果,分别对添加噪声的仿真数据和真实数据进行测试。

1) 添加高斯噪声和离群点的仿真点云分割

根据激光点云的噪声特性,对三种点云添加均值为0,标准差为点云最小外包盒对角线的0.01倍的噪声,并且额外添加数量为点数0.25倍的离群点,得到待测试点云如图6(a)所示,分别利用基于k-means的聚类分割算法^[3]、基于局部角距离的谱聚类分割算法^[21]以及本文算法,重复实验后选取出现次数最多的结果作为最终结果,如图6(b)~(d)所示。

图 6. 不同分割算法对含噪声的点云图像的分割结果。(a)添加噪声的点云;(b) k-means算法;(c)局部角距离谱聚类算法;(d)本文算法

Fig. 6. Segmentation results of point clouds with noises by different segmentation algorithms. (a) Point cloud with noise; (b) k-means algorithm; (c) local angular distance spectral clustering algorithm; (d) proposed algorithm

下载图片 查看所有图片

从图6(b)可以看出,基于k-means的分割算法没有得到正确的分割结果,这是因为基于k-means的分割算法受到点的距离影响较严重,对点云几何结构解析能力较差;从图6(c)可以看出,基于局部角距离的谱聚类分割算法在Rabbit点云中得到了错误的分割结果,且在Armadillo点云中对分割边界等细节的处理仍有一定欠缺,这是因为基于局部角距离提取特征会使得分割结果在本属于一个连续部分但曲率变化较大时容易断开,且在受到噪声干扰时不能准确解析几何特征;而本文算法在三种点云中均表现良好。

2) 切片式激光三维成像的真实点云分割

结合基于距离选通切片式激光三维成像得到的卫星模型点云进行实验,目标的灰度图像如图7(a)所示;切片式激光三维成像的点云具有激光三维成像的主要特征,即点云存在严重的噪声和大量离群点,图7(b)为目标点云侧视图;本文没有对点云进行去噪处理,直接对图7(b)中点云进行分割,结果如图7(c)、(d)所示,可见在点云受到严重噪声干扰的情况下,本文算法依然准确提取出了目标的结构特征,将卫星模型的目标点云分为了实验舱、左右太阳能电池板以及航天员模型,符合模型实际特征,且分割边界良好,充分说明本文算法对激光点云的自适应性。

图 7. 卫星模型点云的分割结果。(a)灰度图像;(b)切片式三维成像法得到的侧视图;(c)本文算法得到的正视图; (d)本文算法得到的侧视图

Fig. 7. Segmentation results of point cloud of satellite model. (a) Gray image; (b) side view of slice three-dimensional imaging; (c) front view of segmentation by the proposed algorithm; (d) side view of segmentation by the proposed algorithm

下载图片 查看所有图片

4 结论

结合激光点云的特点,提出了一种基于混合流形谱聚类的激光点云自适应分割方法,该算法具有良好的结构分析能力和抗噪性。实验结果表明,本文算法对预设参数稳定性好,在较宽参数范围内对三种受测点云均能得到以80%以上概率收敛于合理分割的结果。且在高斯噪声和离群点干扰的情况下表现稳定,适合于激光点云的自适应分割。对切片式三维成像点云进行分割,也取得了较好的结果。

本文算法仍然存在一些不足,例如对较复杂的非线性流形构成的混合点云分割准确率会有所下降,因此还待进一步改进。本文算法得到的分割结果可用于激光点云的配准与目标特征识别,对分割结果的应用也有进一步研究的价值。