张玉存， 李亚彬， 付献斌

(1.燕山大学 电气工程学院，河北 秦皇岛 066004；2.河北环境工程学院 信息工程系，河北 秦皇岛 066102)

1 引 言

随着三维扫描技术的发展与广泛应用，三维扫描点云数据处理技术的研究也日益兴起[1]。由于扫描过程中测量设备的限制和环境干扰，导致测得的三维点云数据不可避免地存在噪声，点云消噪已成为三维几何数据处理领域的重要研究方向[2]。

目前，众多学者提出了多种点云模型去噪算法，主要有通过网格滤波算法、局部突出特征信息等去噪方法。其中根据网格滤波算法去噪的方法中以双边滤波光顺去噪最为普遍，Fleishman S等[3]介绍了一种各向异性的网格去噪算法，通过对图像去噪的双边过滤算法采用它来进行3D网格化，消除了网格中的噪声。Zheng等[4]提出一种双边标准滤波器，通过空间距离和信号距离，该特征保护方案在视觉上和数字上都比以前的方法更有效果，特别是在具有尖锐特征或不规则采样的区域；江亮亮等[5]提出一种用于网格去噪的自适应双边滤波器，根据特征强度自适应地调整双边滤波器的滤波参数进一步提高了去噪性能；Lu 等[6]引入有效光滑面法线的规范双边滤波器和一个与过滤面法线估计顶点法线的有效方案，在顶点更新中，提出了一种新的迭代顶点更新算法来有效地更新顶点位置；崔鑫等[7]利用主元分析法估算点云的曲率和法向量，根据曲率特征标识点云数据的特征区域，使用双边滤波器对点云模型进行平滑去噪；Zhang等[8]建立了基于扩散滤波一般原理的三维扩散滤波方程，推导了尺度参数和最大迭代次数的理论公式，实现了自适应降噪。但是由于双边滤波的平滑滤波因子和特征保持因子在不同特征区域的权重不能自适应调整，所以在滤波过程中可能会出现过光顺，因此研究者根据点云数据局部突出特征信息，如将模型的尖锐部分与其他部分分割，将点云数据分为平坦区与非平坦区等特征进行去噪。吴禄慎等[9]根据点云平均曲率的局部特征权值，针对不同特征区域分别采用邻域距离平均滤波算法和自适应双边滤波算法进行去噪滤波；Zheng等[10]根据特征类型为多个法线分配特征点提取尖锐的特征结构，根据过滤的法线更新点位置，可以稳健地滤除高水平的噪声，同时保持重要的几何特征完整；Mattei等[11]提出了一种移动鲁棒主成分分析(mrpca)，点的估计位置通过局部平均来计算，通过加权最小化来保留尖锐特征，利用权重测量本地邻域中的法向矢量之间的相似性进行消噪；Wang等[12]提出了一种有效的方法来处理基于新型滚动制导正常滤波器不同尺度的几何特征，平滑小规模的几何特征，同时保留大规模的功能；樊宇等[13]运用激光扫描光刀法扫描得到点云数据，提出了一种基于几何关系的三角形滤波法则，能够较好地进行去噪处理。

为避免对曲率较大的数据点的影响以保证边界特征得以保持，本文提出基于曲率约束的点云分割去噪方法。利用最小二乘拟合曲面与过邻域点作拟合曲面的法截面，结合相交曲线的旋转角估算邻域内点的曲率，以邻域点为根节点搜索下一邻域点，根据曲率的约束对点云数据分割；对分割后点集合利用一种映射关系对噪声点进行线性回归，移动噪声点到新位置，实现点模型去噪。

2 曲率约束分割算法

对于三维点云数据集P={Pi=(xi,yi,zi)T,i=1,2…N}，Pi为点云中一个采样点，N为采样点中心，建立坐标轴，以数据的第一列为x轴，第二列为y轴，第三列为z轴，进行曲面分割。

(1)

选取Pi点及其邻域内点，图1为同面内的邻域点，图2为不同曲面内的邻域点。根据Pi点与邻域点所形成三角面片结构，构建不同三角面片的平面拟合，即：

y=AKx+BKy+CKz+DK

图1 选取点P及其邻域点Fig.1 Point P and its neighborhood point

图2 邻域点在不同的曲面Fig.2 Neighborhood points are on different surfaces

设三角面片的法向量a=(Aa，Ba，Ca)，b=(Ab，Bb，Cb)，AK、BK、CK、DK为每一个三角面片所生成平面的系数，其中K为邻域点的个数。一个固定三角面片的法向量与其它法向量的夹角为θ，即：

(2)

当在θn的矩阵中所有的角度差值在一定的阈值内，则Pi的邻域内没有两曲面相交的点或两曲面夹角相对较小；当角度差值超过阈值时，则该邻域内存在曲面相交角度较大的交界点。

对于存在曲面相交点的邻域点进行曲面拟合，如图3、图4所示。

图3 采样点邻域数据点Fig.3 Sampling point neighborhood data point

图4 邻域数据点拟合曲面Fig.4 Neighborhood data point fitting surface

设拟合曲面方程为：

f(x,y)=a1x2+a2y2+b1xy+b2yz+c1x+c2y+d

将各点投影至曲面上，通过微分几何计算曲面上各点的高斯曲率K、平均曲率H：

(3)

(4)

(5)

曲率的差值的数据为：

(6)

取两组数据μj=(X1，X2)做为初始聚类中心，使得差值数据分为两类，则目标函数为：

(7)

μj为分类的簇中心，更新迭代使得簇成员到簇中心差值和为最小。

提取两类中最大值与最小值：

(8)

(9)

则定义采样点的特性为：

(10)

当一曲面中不存在相交点但在某段的曲率偏高时，不利于后期对数据的拟合处理，在曲率较高的情况下运用区域生长的方法对曲率较高的曲面进行分割，可减小数据后期处理时的拟合偏差。

对于除去相交点后的点云数据，对Pi点邻域内的点进行曲面拟合得曲面S0，将Pi点的邻域点投影至拟合曲面，拟合曲面方程为：

f(x,y,z)=a1x2+a2y2+a3z2+b1xy+b2yz+
b3xz+c1x+c2y+c3z+d

图5 曲面的法截面Fig.5 The normal section of a surface

(11)

记法截面的方程为g(x,y,z)=0，法截面S1与曲面S0的相交曲线为l0，则l0的方程表示为：

(12)

(13)

(14)

(15)

(16)

(17)

(18)

则法截面与曲面相交的弧线的旋转角为：

(19)

旋转θ值越大则曲线l0的平均曲率越小，旋转θ值越小则曲线l0的平均曲率越大。

图6 同一曲率范围内的点Fig.6 Points in the range of the same curvature

图7 不同曲率范围内的点Fig.7 Points in the range of different curvature

3 曲率约束消噪方法

以点云数据点为球心，半径为r的球体范围内的点云数据点数为该点的密度值[14]。在以点云数据点为球心，半径为r的球体范围内寻找点云个数。当该数据点为大尺度噪声点或离群点时，以该点为球心半径为r的球体内的点云个数较少，即该点的密度值偏小；当该数据点为正常值时，以该点为球心半径为r的球体内的点云个数较多，即密度值较大。通过设定球体内的点云数的阈值范围，可将大尺度噪声及离群点去除。

对于集合MN{N=1,2,3,…}中含有噪声的点云数据，其中不乏含有噪声点距离真实点较远，偏离主体点云。对集合MN中的点Pi，i为集合中点云的个数，计算其到所有临近点的距离值，设定阈值r，得到在阈值r的范围的点云数据个数成为该点基于距离r的密度，记为density(Pi,r)。根据点云密度设定在距离为r的范围内含有的数据点个数的阈值ξ。当Pi点的密度值density(Pi,r)<ξ，则判定点Pi为噪声点，应予剔除；当Pi点的密度值density(Pi,r)≥ξ，则判定点Pi为合格点。

在集合MN中含有噪声的点云数据中，建立一个三次曲面来逼近数据点，三次曲面方程为：

f(x,y)=a9x3+a8y3+a7x2y+a6xy2+a5x2+
a4y2+a3xy+a2x+a1y+a0

对比图8、图9可知邻域点在经过三次曲面拟合后，曲面在噪声点处的曲率变化优于二次曲面拟合，根据拟合的曲面，求取各点在曲面的投影点，即每一点距离曲面最近的点为该点的投影点，(x0，y0，z0)为数据点，(x，y，f(x，y))为曲面上的点，则距离d满足下式：

d2=(x-x0)2+(y-y0)2+(f(x,y)-z0)2

图8 二次曲面集合内的点Fig.8 Points in a set of two surfaces

图9 三次曲面集合内的点Fig.9 Points in a set of three surfaces

求d的最小值，则：

(20)

解得投影点坐标为：

(21)

取曲率较大的数据点构成协方差矩阵：

(22)

(23)

式中f()为一种映射。对噪声点光顺平移，设噪声点移动距离为(Δxi，Δyi，Δzi)，则：

(24)

(25)

式中：ω为各分量方向的权重系数。

移动值与真实值之间距离插值最小时，移动后的点最接近真实点，设t*=(txi，tyi，tzi)，u*=(Δxi，Δyi，Δzi)，ω*=(ωxi，ωyi，ωzi)表示不同系数权值，则其目标函数为：

(26)

求得系数：

(27)

将ω*代入式(25)得到移动距离(Δxi，Δyi，Δzi)，光顺去除噪声点。当移动完成，对移动后的点进行二次曲面拟合，将各点投影至曲面，做各点的切线，以某一点为基准点根据式(18)判断角度θ。根据阈值η，当所有角度满足θ≤η，对下一数据点集合进行光顺处理；当不满足θ>η时，则对移动后的点进行三次曲面拟合，重复式(19)～(27)步骤，至所有点满足条件，完成点云数据的光顺去噪。

4 实验与分析

本文利用曲率的约束对点云数据进行分割处理，对每个区域进行去噪，在去噪的同时可以保持曲率变化较大区域的细节特征；设计了一种映射关系，使得噪声点利用曲率的变化回归至真实点区域。将用双边滤波的处理结果与本文处理结果进行了比较。

图10中展示点云模型去除大尺度噪声结果，其中，图10(a)为含有大尺度噪声点云模型，图10(b)为去除大尺度噪声点云模型的结果。通过对比图10(a)实线矩形框和虚线矩形框与图10(b)实线矩形框和虚线矩形框中的点云数据，可知通过设定点云密度的阈值可以将大尺度点云数据噪声进行去除，并且点云模型保持原来特征。

图10 点云模型去除大尺度噪声结果Fig.10 Point cloud model to remove large scale noise results

表1显示点云模型经过点云密度阈值去除大尺度噪声的结果，设定密度球的半径r为0.09 mm，原始点云数据为12 950个，经点云密度阈值去除大尺度噪声后的点云数为10 132个。去除大量大尺度噪声与离群点，并使得模型的顶点与边缘特征得到良好的保持。

表1 大尺度噪声去噪结果Tab.1 The result of large scale noise de-noising

Fandisk噪声模型在经过双边滤波处理与本文算法处理后的比较，如图11所示。在图11(d)与图11(e)中即图11(b)与图11(c)中的区域1，通过对比图11(d)与图11(e)可以看出在曲率变化较大的区域经过双边滤波处理后仍有明显噪声影响，表面凸起，有明显噪声点残留；本文算法对曲率变化较大的区域的特征保持效果较好，噪声去除效果较好。在图11(f)与图11(g)中即图11(b)与图11(c)中的区域2，通过对比图11(f)与图11(g)可以看出在两个曲面相交的邻域范围内，双边滤波处理后的曲面有明显噪声存在，模型存在较多的噪声棱角；本文处理后的曲面光顺效果较好，且模型无明显凸起棱角。在图11(h)与图11(i)中即图11(b)与图11(c)中的区域3，通过对比图11(h)与图11(i)可以看出，本文算法在边缘特征保持与曲面光顺方面效果较好。

图11 Fandisk噪声模型去噪结果比较Fig.11 Comparison of de-noising results of Fandisk noise model

图12为Fandisk噪声模型边缘及顶点去噪结果。通过对比图12(a)、图12(b)中的区域1和区域2中的边缘信息，可知双边滤波进行消噪时对模型的边缘信息保留不够完整，存在较多噪声棱角，边缘信息处理粗糙。对比图12(a)、图12(b)中的的区域3，可知在尖锐的顶点处，双边滤波处理将顶点数据平顺，顶点位置及边缘发生变化；本文算法的去噪结果对模型的边缘及顶点特征保存完整。

由图11、12及表2可以看出：在运用双边滤波进行消噪时，选取最佳参数σc=10，σs=20；运用本文的算法时选取密度球的半径r为0.09 mm。

表2 双边滤波与本文算法的去噪结果对比Tab.2 Comparison of denoising results between bilateral filtering and our algorithm

经过双边滤波的算法处理，能够去除部分噪声数据，但在模型的边缘特征及曲率较大区域的特征保持方面小规模稍差；本文算法处理后，能够有效地去除数据的冗余，同时能够保持模型的边缘特征及曲率较大区域的特征。

图12 Fandisk噪声模型边缘及顶点去噪结果Fig.12 Edge and vertex denoising results of Fandisk noise models

5 结 论

本文提出了一种基于曲率约束的点云分割去噪方法。根据采样点与其邻域点在局部拟合曲面上的投影点之间的关系构建曲率权重约束条件，下一采样点的与其邻域点在局部拟合曲面上投影点之间同样采用曲率权重约束条件，判断两邻域是否添加至同一点云数据集合中，对点云数据进行曲率分割；运用协方差分析法对每一点云数据集合得出噪声点的特征向量，构建一种映射关系，实现点云数据去噪。实验结果表明，本文算法在确保消除点云噪声的同时，能有效地保存点云曲面变化较大时的信息特征，提高去噪效果，防止消噪过程中点云模型的扭曲变形。