基于DBSCAN的改进RANSAC点云平面拟合算法

Improved RANSAC Point Cloud Plane Fitting Algorithm Based on DBSCAN

摘要

针对随机采样一致性（random sample consensus，RANSAC）算法对含有噪声的点云数据进行平面拟合时效果不佳和容易产生误识别的问题，对算法进行改进. 通过基于密度的噪声应用空间聚类（density-based spatial clustering of applications with noise，DBSCAN）算法改变RANSAC算法初始点集合的选择策略，并使用主成分分析法（principal component analysis，PCA）计算点云各点法向量，以点到平面距离以及点的法向量与平面法向量夹角两个约束条件同时作为RANSAC算法平面拟合模型内点判定的准则. 采用无噪声与分别含有300个噪声点和500个噪声点的点云仿真数据进行测试，本文算法拟合结果均接近理论值且内点距离标准差分别为1.007×10^-8、0.003、0.007，优于RANSAC算法. 采用实际工件点云数据在两种工况场景下进行测试，本文算法拟合平面内点比率相对于传统RANSAC算法分别提高24.7%和24.6%，平面提取完整度及准确率同样优于RANSAC算法. 仿真模拟及实例分析证明了本文算法的有效性.

Abstract

To address the problem of poor point cloud plane fitting effect and easy misidentification of the random sample consensus （RANSAC） algorithm for noisy point cloud data， improvements to the algorithm are necessary. The proposed algorithm employs density-based spatial clustering of applications with noise （DBSCAN） to modify the selection strategy of the initial point set in the RANSAC algorithm and uses principal component analysis （PCA） to compute the normal vectors of each point in the point cloud. Two constraint conditions， the distance from the point to the plane and the angle between the normal vectors of the point and the plane， are simultaneously used as the criteria for determining the points within the RANSAC algorithm model. The point cloud simulation data with no noise and 300 and 500 noisy points are used for testing. The fitting results of the proposed algorithm are all approximate to the theoretical values， and the standard deviations of the inner point distance are 1.007×10^-8， 0.003 and 0.007， respectively， better than those of RANSAC algorithm. Using actual workpiece point cloud data for testing in two operating scenarios， the proposed algorithm improves the fitting ratio of in-plane points by 24.7% and 24.6% compared to the traditional RANSAC algorithm， respectively. The completeness and accuracy of plane extraction are also superior to those of the RANSAC algorithm. Simulation and case analysis validate the effectiveness of the proposed algorithm.

Graphical abstract

关键词

点云平面拟合 / 随机采样一致性 / 噪声 / 密度聚类 / 主成分分析

Key words

point cloud plane fitting / random sample consensus / noise / density clustering / principal component analysis

随着人工智能的不断发展，三维空间信息的获取与处理显得越来越重要^［1］. 点云数据即空间中点的集合，是物体几何形状的一种具体数字化体现. 点云数据的获取方法主要有激光测距法、结构光法、RGB-D相机法、三维扫描法等，由于其相较二维图像，具有可量化复杂环境^［2］，提供三维空间信息的特点，被广泛应用于三维重建、增强和虚拟现实、智能导航等领域^［3］. 点云的丰富信息为研究者提供了深入理解和模拟真实场景的机会，近年来，对点云数据的处理成为研究热点.

点云场景中常包含大量平面特征，故点云平面拟合作为一项基础的任务，在实际应用中意义重大. 例如，在机器人导航与定位中，通过拟合地面点云平面，可以帮助机器人实现精准的定位和路径规划；在无人驾驶和智能车辆领域，点云平面拟合也被广泛应用于道路和交通场景的建模和理解，为自动驾驶系统提供重要的环境感知能力；在工业制造中，通过对工件表面进行点云数据采集和平面拟合，可以实现对工件表面的几何特征和质量状态进行检测和分析，为工业制造提供了高效、精准的质量控制解决方案. 综上，由于点云平面拟合对于场景理解、物体识别和机器人感知等应用起着至关重要的作用，故实现精准的点云平面拟合非常必要. 点云平面拟合的办法包括特征值法、最小二乘法、霍夫变换法、随机采样一致性（random sample consensus，RANSAC）等，其中，由于RANSAC算法较其他算法具有一定鲁棒性，故应用较广. 但RANSAC算法的鲁棒性并不是绝对的，在点云数据含有较多噪声点时，仍然容易导致平面拟合精度下降. Zheng等^［4］提出了一种结合RANSAC算法和改进特征值算法的点云平面拟合方法，提高了平面拟合的精度. Kang等^［5］基于贝叶斯定理，使用多个初始基元并行检测的策略优化RANSAC算法. Chen等^［6］采用局部采样方法改进RANSAC算法，在针对建筑物屋顶点云平面拟合提取时，有良好效果. Oehler等^［7］通过霍夫变换，生成多分辨率聚类，然后再使用RANSAC算法识别拟合平面，以此提升算法准确性，但霍夫变换需要耗费较大内存，有着较高的计算成本. Li等^［8］使用正态分布变换，通过在每次循环中选取一个平面无损检测单元作为最小样本的方式，提高了RANSAC算法在采样时的准确度，但单元大小的选取较为复杂烦琐. Yang等^［9］采用启发式搜索策略来构建初始点集合，能较为有效地提取室内空间环境的平面，但其内点判定仍然只以距离阈值为标准，存在误识别的可能性. 通过上述分析可知，现有改进算法大多面临着复杂的计算问题，且内点判定准则均较为单一，而错误的内点判定将对最终的拟合结果造成较大影响.

本文基于DBSCAN，在经过预处理的点云数据中搜索RANSAC平面拟合算法所需要的初始点集合，并使用主成分分析法计算点云数据中各点的法向量，增加法向量夹角这一约束条件，将其与点到平面的距离一同作为内点判定准则. 与其他改进算法相比，本文算法避免了高成本的计算，方便算法后续的硬件实现，且完善了内点判定准则，进一步降低误识别概率. 仿真试验与实物实验分析表明，本文算法较常规RANSAC算法，具有更好的拟合效果.

1 RANSAC算法

RANSAC算法是一种鲁棒性的估计分析方法，通过不断迭代的方式，估算拟合数据集的模型参数. 使用RANSAC算法可以将数据集划分为内点（inliers）和外点（outliers），内点即分布在所求模型上的点，外点即模型之外的点.

RANSAC算法原理本质上可概括如下：首先，从数据集中随机选取最小样本集，这个样本集的大小通常取决于所要拟合模型的参数个数，使用所选的样本集对模型进行拟合，求出模型的初始参数. 接着，给定一个合适的阈值，计算所有数据点到此前求出的拟合模型的误差，确定数据集中所有符合模型要求的数据点，拟合误差小于或等于阈值的点，被认为是内点，而其余点，则被认为是外点. 多次重复上述两个步骤，选择拥有最多内点或最小内点误差的模型作为最终估计.

当RANSAC算法具体应用到拟合空间点云平面时，其平面拟合示意图如图1所示.

基本步骤如下：

1）由于至少3个点才能确定一个平面，所以拟合平面时，需要从数据集中随机选取3个点作为初始的最小样本集，针对这3个点计算平面模型参数.

2）设置合适的距离阈值d_t，计算所有剩下的点到该平面的距离d_i，将符合d_i\led_{\mathrm{t}}的点视作内点，将其余点视作外点. 当内点数量满足要求时，保留此平面模型.

3）继续重复上述两步，如果当前模型的内点数量大于此前保存的最大内点数量，则更新模型参数. 不断迭代，直到达到迭代阈值，将内点个数最多的模型作为最终的平面模型.

2 改进的RANSAC算法

使用RANSAC算法进行点云平面拟合时，由于对原始点云数据集中3个初始样本点的选取完全随机，极大增加了取到数据集中平面外的点的可能性，由此得到的模型参数往往无法满足要求，当采样次数一定时，初始点选取平面外的点的次数越多，符合要求的模型集就越少，最后得到最优模型的概率就越小，平面拟合的准确度也将越低. 同时，RANSAC算法对内点的判定仅仅以数据集中的点到平面模型的距离为依据，只要符合距离要求的点就被视为内点，这也有可能导致误识别. 针对以上问题，本文提出一种基于DBSCAN的改进RANSAC算法，通过改变初始点的选择策略以及内点的判定方式，提高RANSAC算法的准确度. 改进RANSAC点云平面拟合算法流程图如图2所示，首先对原始点云进行预处理，接着使用主成分分析法（principal component analysis，PCA）计算各点法向量^［10］，通过DBSCAN寻找RANSAC算法的初始点集合，在初始点集合中随机选取3个点进行平面拟合，将满足点到平面距离阈值和点与拟合平面法向量角度阈值的点判定为内点，保留内点个数满足要求的模型，迭代后选择内点个数最多的模型作为平面拟合结果.

2.1 点云滤波与下采样

一般来说，3D点云原始数据往往包含较多的噪声以及离群点，且数据量庞大，为了进一步提高平面拟合的准确度和算法效率，需要对点云进行预处理.

基于统计学原理利用统计滤波去除大部分噪声点^［11］. 通过构造多维二叉树（K-dimension tree，KD）的方式，建立点云间的拓扑关系，并由此找到每一个点的k_{\mathrm{1}}个近邻点，即距离每一个点最近的k_{\mathrm{1}}个点. 计算每一个点与其k_{\mathrm{1}}个近邻点之间的平均距离，并计算这些平均距离的均值与标准差，计算方法如下：

{\overline{S}}_j=\frac{\mathrm{1}}{k_{\mathrm{1}}}{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{1}}}}{S}_{ij}（1）

\mu =\frac{\mathrm{1}}{n}{\displaystyle \sum _{j=\mathrm{1}}^n}{\overline{S}}_j（2）

\sigma =\sqrt[]{\frac{\mathrm{1}}{n-\mathrm{1}}{{\displaystyle \sum _{j=\mathrm{1}}^n}({\overline{S}}_j-\mu )}^{\mathrm{2}}}（3）

式中：{\overline{S}}_j为点云中第j个点与其k_{\mathrm{1}}个近邻点之间的平均距离；S_{ij}为点云中第j个点到其第i个近邻点的距离；\mu为平均距离的均值；\sigma为标准差；n为点云中点的个数. 设s为标准差倍数，当某个点与其k_{\mathrm{1}}个近邻点之间的平均距离在区间\left[\mu -s\sigma ,\mu +\right.\left.s\sigma \right]内时，将此点保留，否则，将其视为噪声点，予以剔除.

统计滤波示意图如图3所示，当k_{\mathrm{1}}=3时，以图中A点与B点为例，A点到其距离最近的3个点的平均距离满足阈值区间要求，将其保留，而B点到其距离最近的3个点的平均距离超出要求范围，将其视作噪声点，予以剔除.

通过创建三维体素栅格的方法，简化点云数据^［12］. 根据点云数据的坐标，分别求出X、Y、Z 3个坐标轴上的最大值x_{\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{、}{y}_{\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{、}{z}_{\mathrm{m}\mathrm{a}\mathrm{x}}及最小值x_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{、}{y}_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{、}{z}_{\mathrm{m}\mathrm{i}\mathrm{n}}，由此计算3个坐标轴上的点云最小包围盒边长l_x\mathrm{、}{l}_y\mathrm{、}{l}_z，计算过程如下：

\left\{\begin{array}{l}{l}_x=x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}\\ l_y=y_{\mathrm{m}\mathrm{a}\mathrm{x}}-y_{\mathrm{m}\mathrm{i}\mathrm{n}}\\ l_z=z_{\mathrm{m}\mathrm{a}\mathrm{x}}-z_{\mathrm{m}\mathrm{i}\mathrm{n}}\end{array}\right.（4）

由所得的点云最小包围盒边长计算体素栅格的尺寸，计算过程如下：

\left\{\begin{array}{l}{D}_x=⌊\frac{l_x}{r}⌋\\ D_y=⌊\frac{l_y}{r}⌋\\ D_z=⌊\frac{l_z}{r}⌋\end{array}\right.（5）

式中：r为本文所设置的体素单元边长；D_x\mathrm{、}{D}_y\mathrm{、}{D}_z分别为3个坐标轴方向上体素栅格的尺寸；符号⌊⌋表示向下取整. 计算点云中每个体素单元内的索引编号h，计算过程如下：

\left\{\begin{array}{l}{h}_x=⌊\frac{x-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{r}⌋\\ h_y=⌊\frac{y-y_{\mathrm{m}\mathrm{i}\mathrm{n}}}{r}⌋\\ h_z=⌊\frac{z-z_{\mathrm{m}\mathrm{i}\mathrm{n}}}{r}⌋\\ h=h_x+h_y{D}_x+h_z{D}_x{D}_y\end{array}\right.（6）

对索引h里的元素按从小到大排列，求出每一个体素单元的重心点，用重心点代替体素单元内的所有点，若体素单元的重心点不存在，则使用体素单元内距离重心最近的数据点代替重心点. 通过使用体素降采样，达到简化缩小点云数据且保留点云形状结构特点的目的. 体素降采样原理图如图4所示.

2.2 点云法向量估计

为了估计点云的法向量，采用基于邻域的近似计算方法，该方法的主要思想如下.

对于点云中的每一个点，通过KD，再次搜索与其距离最近的k_{\mathrm{2}}个近邻点，用最小二乘法拟合这些最近邻域点的局部平面，表示如下：

P(\mathit{n},d)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{(\mathit{n},d)}{\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}(\mathit{n}\cdot {\mathit{p}}_i{-d)}^{\mathrm{2}}（7）

式中：\mathit{n}为局部平面P的单位法向量；d为坐标轴原点到局部平面P的距离；{\mathit{p}}_i为当前点局部邻域中的第i个点，坐标为(x_i,y_i,z_i).

把由此拟合出的局部平面的单位法向量称作当前点的法向量^［13-15］，而局部平面的单位法向量可以使用主成分分析法求得. 主成分分析法通过降低维度的思想，利用线性变换，将一组变量转化为另一组不相关的变量，转化后，变量总方差不变，且最大的方差在第1个分量上，称为第1主成分，第2大的方差在第2个分量上，称为第2主成分，以此类推. 点云数据的变量为坐标的集合，通过变换，得到3个主成分. 对于点云中的平面，由于垂直于平面的方向，点云分布集中，方差最小，所以平面的第3主成分即为平面的单位法向量，如图5所示，求解平面单位法向量的问题转为求解平面点构成的协方差矩阵最小特征值对应的特征向量的问题.

由分析可知，局部平面P经过当前点k_{\mathrm{2}}个最近邻域点的三维质心{\mathit{p}}_{\mathrm{0}}(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}})，其坐标值如下：

\left\{\begin{array}{l}{x}_{\mathrm{0}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}{x}_i}{k_{\mathrm{2}}}\\ y_{\mathrm{0}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}{y}_i}{k_{\mathrm{2}}}\\ z_{\mathrm{0}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}{z}_i}{k_{\mathrm{2}}}\end{array}\right.（8）

通过当前点的k_{\mathrm{2}}个最近邻域点构建协方差矩阵 M ：

\mathit{M}=\frac{\mathrm{1}}{k_{\mathrm{2}}}{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}({\mathit{p}}_i-{\mathit{p}}_{\mathrm{0}})({\mathit{p}}_i-{\mathit{p}}_{\mathrm{0}}{)}^{\mathrm{T}}（9）

式中：上标T代表转置运算. 展开可得：

\mathit{M}=\left[\begin{array}{ccc}{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{x}_i\mathrm{\Delta }{x}_i& {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{x}_i\mathrm{\Delta }{y}_i& {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{x}_i\mathrm{\Delta }{z}_i\\ {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{x}_i\mathrm{\Delta }{y}_i& {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{y}_i\mathrm{\Delta }{y}_i& {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{y}_i\mathrm{\Delta }{z}_i\\ {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{x}_i\mathrm{\Delta }{z}_i& {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{y}_i\mathrm{\Delta }{z}_i& {\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{2}}}}\mathrm{\Delta }{z}_i\mathrm{\Delta }{z}_i\end{array}\right]（10）

式中：\mathrm{\Delta }{x}_i=x_i-x_{\mathrm{0}};\mathrm{\Delta }{y}_i=y_i-y_{\mathrm{0}};\mathrm{\Delta }{z}_i=z_i-z_{\mathrm{0}}. 对协方差矩阵\mathit{M}进行特征值分析，即：

\mathit{M}\cdot {\mathit{v}}_j={\lambda }_j\cdot {\mathit{v}}_j,j=\mathrm{0,1},\mathrm{2}（11）

式中：{\mathit{v}}_j为协方差矩阵的第j个特征向量；{\lambda }_j为协方差矩阵的第j个特征值. 协方差矩阵\mathit{{\rm M}}的最小特征值所对应的特征向量即为局部平面P的单位法向量，即当前点的法向量. 由于所求点的法向量通常具有二义性，所以需要指定视点{\mathit{p}}_{\mathrm{v}}对法向量进行定向，过程如下：

{\mathit{n}}_{\mathrm{c}\mathrm{u}\mathrm{r}}=\left\{\begin{array}{cc}{\mathit{n}}_{\mathrm{c}\mathrm{u}\mathrm{r}},& {\mathit{n}}_{\mathrm{c}\mathrm{u}\mathrm{r}}\cdot ({\mathit{p}}_{\mathrm{v}}-{\mathit{p}}_{\mathrm{c}\mathrm{u}\mathrm{r}})>\mathrm{0}\\ -{\mathit{n}}_{\mathrm{c}\mathrm{u}\mathrm{r}},& {\mathit{n}}_{\mathrm{c}\mathrm{u}\mathrm{r}}\cdot ({\mathit{p}}_{\mathrm{v}}-{\mathit{p}}_{\mathrm{c}\mathrm{u}\mathrm{r}})\le \mathrm{0}\end{array}\right.（12）

式中：{\mathit{n}}_{\mathrm{c}\mathrm{u}\mathrm{r}}为当前点的法向量；{\mathit{p}}_{\mathrm{c}\mathrm{u}\mathrm{r}}为当前点坐标.

2.3 DBSCAN求初始点集合

DBSCAN是一种基于密度的空间聚类算法，通过分析邻域进而分析样本集的紧密程度^［16-19］，参数(r_{\mathrm{e}\mathrm{p}\mathrm{s}},k_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{P}\mathrm{t}\mathrm{s}})用来描述样本数据点分布的紧密情况. 其中，r_{\mathrm{e}\mathrm{p}\mathrm{s}}为邻域半径，即当前点的邻域距离阈值，k_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{P}\mathrm{t}\mathrm{s}}为邻域半径范围内数据点的最小个数. 关于DBSCAN，有

1）核心对象：对于点云数据中的任意一点p，如果在其以r_{\mathrm{e}\mathrm{p}\mathrm{s}}为邻域半径的邻域范围内，至少包含k_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{P}\mathrm{t}\mathrm{s}}个样本点（包括当前点自身），则称点p为核心对象.

2）密度直达：如果点p是核心对象，当点q位于点p以r_{\mathrm{e}\mathrm{p}\mathrm{s}}为邻域半径的邻域范围内时，则称点q由点p密度直达.

3）密度可达：如果在点云样本中，存在序列点p_{\mathrm{1}},p_{\mathrm{2}},p_{\mathrm{3}},\cdots ,p_m，满足点p_{i+\mathrm{1}}由点p_i密度直达，且i=\mathrm{1,2},\mathrm{3},\cdots ,m-\mathrm{1}，则称点p_m由点p_{\mathrm{1}}密度可达.

4）密度相连：如果对于点p和点q，存在核心对象点o，使得点p和点q均可由点o密度可达，则称点p和点q密度相连. 将密度相连的点的最大集合定义为簇.

DBSCAN示意图如图6所示，当k_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{P}\mathrm{t}\mathrm{s}}=5时，红色圆点的邻域半径范围内的样本点都至少为5个，所以红色圆点都是核心对象，位于它们各自邻域半径范围内的点（除去它们自身），即图6中的绿色圆点，都可由它们密度直达. 图6中，使用蓝色箭头连接起来的核心对象，构成一组密度可达的核心对象样本序列，这些样本序列的邻域半径范围内的所有样本点相互都是密度相连的，所有红色圆点与绿色圆点构成一个簇.

运用DBSCAN，寻找初始点集合的步骤如下：在平面的密度紧密区域，根据事先设置的r_{\mathrm{e}\mathrm{p}\mathrm{s}}\mathrm{、}{k}_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{P}\mathrm{t}\mathrm{s}}，随机选择一个核心对象点，将其标记为已处理，通过KD树，搜索可由其密度直达的点，将这些点加入队列，在队列中，继续寻找未处理的核心对象点，并重复上述步骤直到队列中找不到未处理的核心对象点为止，由此得到的队列集合，便构成一个点云簇. 经过DBSCAN后，得到的这个点云簇，由于是在平面密度紧密区域计算的密度相连点的最大集合，簇内点为平面内点的概率将大大增加，所以，这个点云簇便是本文所求的初始点集合.

2.4 平面模型拟合及内点判定

在所求的初始点集合内，随机抽取3个点，根据这3个点的坐标，计算求得平面模型参数. 空间平面方程如下：

Ax+By+Cz+D=\mathrm{0}（13）

式中：A、B、C为平面单位法向量的3个坐标分量值，即满足A^{\mathrm{2}}+B^{\mathrm{2}}+C^{\mathrm{2}}=\mathrm{1}；D为坐标原点到平面的距离. 使用点云数据集中余下的所有点对平面模型进行验证.

计算所有剩余点到拟合平面的距离：

\begin{array}{l}{d}_i=\frac{\left|A{x}_i+B{y}_i+C{z}_i+D\right|}{\sqrt[]{A^{\mathrm{2}}+B^{\mathrm{2}}+C^{\mathrm{2}}}}=\\ \left|A{x}_i+B{y}_i+C{z}_i+D\right|\end{array}（14）

设置距离阈值d_{\mathrm{t}}，将d_i>d_{\mathrm{t}}的点直接归为外点，对d_i\le d_{\mathrm{t}}的点进行进一步判断.

计算d_i\le d_{\mathrm{t}}的点法向量与拟合平面的法向量夹角：

{\theta }_j=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{{\mathit{n}}_{\mathrm{p}}\cdot {\mathit{n}}_{\mathrm{j}}}{\left|{\mathit{n}}_{\mathrm{p}}\right|\left|{\mathit{n}}_{\mathrm{j}}\right|}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}({\mathit{n}}_{\mathrm{p}}\cdot {\mathit{n}}_{\mathrm{j}})（15）

式中：{\mathit{n}}_{\mathrm{p}}为拟合平面的法向量；{\mathit{n}}_{\mathrm{j}}为满足距离要求的点的法向量；\left|{\mathit{n}}_{\mathrm{p}}\right|,\left|{\mathit{n}}_{\mathrm{j}}\right|为法向量的模，其值均为1.

点与拟合平面法向量夹角如图7所示. 理论上，平面上的点，其法向量与平面法向量相互平行，由此设置角度阈值\alpha，将{\theta }_j\le \alpha的点视为内点. 本文算法的内点判定准则为：将同时满足点到平面距离阈值要求和点与拟合平面法向量夹角阈值要求的点视为内点，以此降低误识别概率.

记录内点个数满足要求的模型，重复上述步骤，直到采样次数达到迭代次数阈值，结束迭代. 迭代结束后，从记录的模型中选出内点个数最多的平面模型作为最终平面拟合的结果.

RANSAC算法的迭代次数阈值与抽取到合格样本概率的关系为：

\omega =\mathrm{1}-(\mathrm{1}-{\epsilon }^{\mathrm{3}}{)}^N（16）

式中：\omega为采样过程中选取的3个初始点均为平面点的概率；\epsilon为从初始点集合中选取到一个平面点的概率；N为迭代次数阈值. 由式（16）可知，本文算法相对传统RANSAC算法，在经过DBSCAN得到初始点集合之后，\epsilon将得以提高，在迭代次数阈值N相同的情况下，选取到的3个初始点均为平面点的概率\omega也将提高，即N相同时，本文算法的点云平面拟合准确度优于传统RANSAC算法.

3 实 验

为了验证本文算法的效果，基于\mathrm{W}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{s}\mathrm{10}系统和\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{l}(\mathrm{R})\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{e}(\mathrm{T}\mathrm{M})\mathrm{i}\mathrm{7}-\mathrm{7700}\mathrm{H}\mathrm{Q}\mathrm{C}\mathrm{P}\mathrm{U}@\mathrm{2.80}\mathrm{G}\mathrm{H}\mathrm{z}处理器，在Visual Studio2017上用C++语言进行算法程序的编写.

3.1 仿真分析

使用仿真数据进行试验分析. 令待拟合的点云空间平面方程为：

-\frac{\mathrm{1}}{\mathrm{5}}x+\frac{\mathrm{2}\sqrt[]{\mathrm{6}}}{\mathrm{5}}z=\frac{\sqrt[]{\mathrm{6}}}{\mathrm{5}}（17）

含不同数量异常噪声的仿真数据如图8所示.使用MATLAB软件随机选取待拟合点云平面上的 1 000个点，如图8（a）所示. 分别在1 000个点的仿真数据中任意加入300个和500个异常值作为噪声点，如图8（b）和图8（c）所示，其中红色点为待拟合平面上的点. 分别用传统RANSAC算法、最小二乘法、结合特征值的改进RANSAC算法（以下简称特征RANSAC算法）和本文算法对图8中的各仿真数据进行平面拟合，得到的仿真结果分别如表1、表2和表3所示，各算法对图8（b）和图8（c）仿真数据的平面拟合效果分别如图9和图10所示.

表1、表2和表3中的{\sigma }_{\mathrm{d}}为拟合平面的内点到拟合平面的距离标准差，计算方法如下：

{\sigma }_{\mathrm{d}}=\sqrt[]{\frac{{\displaystyle \sum _{j=\mathrm{1}}^{n_{\mathrm{i}\mathrm{n}}}}{\left(d_j-\frac{\mathrm{1}}{n_{\mathrm{i}\mathrm{n}}}{\displaystyle \sum _{j=\mathrm{1}}^{n_{\mathrm{i}\mathrm{n}}}}{d}_j\right)}^{\mathrm{2}}}{n_{\mathrm{i}\mathrm{n}}-\mathrm{1}}}（18）

式中：n_{\mathrm{i}\mathrm{n}}为内点个数；d_j为第j个内点到拟合平面的距离. 由表1可知，当点云数据中没有异常噪声点干扰时，4种算法所得的平面系数结果与设定的理论值近似相同，拟合效果都较为理想，但本文算法所得平面系数与设定的理论值最为接近，且距离标准差{\sigma }_{\mathrm{d}}最小，所拟合的平面离散程度最低，拟合精度高于另外3种算法. 由图9和图10可知，当点云数据中分别包含300个和500个异常噪声点时，最小二乘法拟合效果均表现最差，特征RANSAC算法所拟合平面与理论平面上的点均存在一定偏离，但相比传统RANSAC算法具有改进；本文算法所拟合平面均与理论平面上的点最为贴合，拟合效果均表现最佳. 结合表2和表3可知，当点云数据中含有异常噪声点时，最小二乘法受干扰最为严重，所得平面系数与设定的理论值具有较大偏差，且噪声点越多，偏差越大. 由于传统RANSAC算法具有一定鲁棒性，故所得结果略好于最小二乘法，但拟合结果与预测精度同样欠佳. 特征RANSAC算法拟合效果优于传统RANSAC算法，但所得平面系数与设定的理论值也存在一定偏差. 本文算法所得平面系数，在两种噪声情况下，均接近理论设定值，且距离标准差相较其他3种算法最小，拟合效果与预测精度以及算法稳定性均优于其他3种算法.

3.2 实例验证

为了更好地验证本文算法的效果，对实际工件点云进行实验分析. 通过RVCX3D结构光相机，在水平放置工件以及竖直放置工件两种工况下，扫描并获取带V形槽矩形工件的点云数据.带V形槽的矩形工件如图11所示.

RVCX3D结构光相机的测量工作距离为200~ 3 000 mm，Z向精度最高可达0.005 mm，点云合成频率最高可达15 Hz. 从获取到的点云数据中提取工件点云，并进行滤波与体素降采样，由此得到的矩形工件点云如图12所示.

分别采用传统RANSAC算法、最小二乘法、特征RANSAC算法和本文算法，在两种工况下，对工件上表面的平面点云进行拟合，并利用拟合结果，分割提取平面点云. 工件上不属于上表面的点以及工件外的点均为噪声点. 水平放置和竖直放置工件时各算法的拟合效果分别如图13和图14所示，分割工件上表面点云效果分别如图15和图16所示，工件上表面的拟合结果分别如表4和表5所示.

由图13可知，在水平放置工件时，本文算法所拟合平面相较于另外3种算法，与矩形工件上表面最为贴合. 观察图15可知，使用4种算法对水平放置的工件上表面点云进行分割时，最小二乘法丢失较多平面信息［图15（b）］，对比图15（a）与图15（b）可知，传统RANSAC算法分割效果优于最小二乘法，但提取到的平面点云完整度欠佳. 对比图15（a）与 图15（c）可知，特征RANSAC算法分割效果较传统RANSAC算法有所提升，但提取到的平面点云仍不完整，且边缘处存在误识别的噪声点. 由图15可知，本文算法分割所得平面完整度较高，且不存在明显误识别噪声点，分割效果优于另外3种算法. 由表4可知，在对水平放置的带V形槽的矩形工件上表面的平面点云进行拟合时，使用本文算法得到的距离标准差为0.007，小于其他3种算法，并且内点比率高于另外3种算法，相较传统RANSAC算法提高了24.7%，相较特征RANSAC算法提高了10.9%，且点云离散程度最低，平面拟合精度高于其他3种算法，所得结果与图示结果吻合. 分析图14和图16以及表5可知，当工件竖直摆放时，本文算法对工件上表面的拟合效果与分割效果和工件水平放置时所得结论一致，同样优于另外3种算法.内点比率相较传统RANSAC算法提高了24.6%，相较特征RANSAC算法提高了9.4%.

4 结 论

本文基于DBSCAN，改变传统RANSAC算法完全随机的初始点选点策略，并使用PCA法计算点云各点法向量，增加法向量夹角这一内点判定条件，来提高RANSAC算法拟合点云平面的准确性. 仿真试验与实例分析结果均表明，本文所提算法相较传统RANSAC算法以及结合特征值的改进RANSAC算法，具有更为可靠的拟合效果，在点云平面拟合方面具有一定的实际意义. 但本文算法也存在局限性，需要人为设置法向量夹角阈值以及点到平面距离阈值，这将是未来需要完善以及深入研究的地方.

参考文献

[1] 李明磊，李广云，王力，等 ．采用八叉树体素生长的点云平面提取［J］．光学精密工程，2018，26（1）： 172-183．

[2] LI M L， LI G Y， WANG L，et al ．Planar feature extraction from unorganized point clouds using octree voxel-based region growing［J］. Optics and Precision Engineering，2018，26（1）：172-183．（in Chinese）

[3] 官云兰， 程效军， 施贵刚. 一种稳健的点云数据平面拟合方 法［J］．同济大学学报（自然科学版），2008，36（7）：981-984．

[4] GUAN Y L， CHENG X J， SHI G G ．A robust method for fitting a plane to point clouds［J］．Journal of Tongji University （Natural Science），2008，36（7）： 981-984．（in Chinese）

[5] 李宗民，姚纯纯，刘玉杰，等 ．点云场景下基于结构感知的车辆检测［J］．计算机辅助设计与图形学学报，2021，33（3）：405-412．

[6] LI Z M， YAO C C， LIU Y J，et al ．Vehicle detection based on structure perception in point cloud［J］．Journal of Computer-Aided Design & Computer Graphics，2021，33（3）：405-412．（in Chinese）

[7] ZHENG L M， WANG R D， WANG S Y，et al.Point cloud plane fitting based on RANSAC and robust eigenvalue method［C］//2022 IEEE 8th International Conference on Computer and Communica- tions （ICCC）， December 9-12，2022.Chengdu，China： IEEE，2022：1368-1372．

[8] KANG Z Z， LI Z. Primitive fitting based on the efficient multiBaySAC algorithm［J］. PLoS One， 2015， 10（3）： e0117341．

[9] CHEN D， ZHANG L Q， MATHIOPOULOS P T，et al ．A methodology for automated segmentation and reconstruction of urban 3-D buildings from ALS point clouds［J］．IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 2014， 7（10）： 4199-4217．

[10] OEHLER B， STUECKLER J， WELLE J，et al ．Efficient multi-resolution plane segmentation of 3D point clouds［M］//Lecture Notes in Computer Science．Berlin，Heidelberg：Springer Berlin Heidelberg，2011：145-156．

[11] LI L， YANG F， ZHU H H，et al ．An improved RANSAC for 3D point cloud plane segmentation based on normal distribution transformation cells［J］．Remote Sensing，2017，9（5）：433．

[12] YANG L N， LI Y C， LI X C，et al ．Efficient plane extraction using normal estimation and RANSAC from 3D point cloud［J］．Computer Standards & Interfaces，2022，82：103608．

[13] 李茂月，赵伟翔，马康盛，等 ．结构光检测点云精简与重构参数自动调节方法［J］．仪器仪表学报，2022，43（8）：122-130．

[14] LI M Y， ZHAO W X， MA K S，et al ．Point cloud simplification and reconstruction parameters’ automatic adjustment method of structured light detection［J］．Chinese Journal of Scientific Instrument，2022，43（8）：122-130．（in Chinese）

[15] OUYANG Q， LIN Y H， ZHANG X L，et al ．Application of 3D reconstruction technology based on an improved MC algorithm in a shotcreting robot［J］．Applied Optics，2022，61（29）：8649-8656．

[16] SHI W J， XU J W， ZHU D C，et al ．RGB-D semantic segmentation and label-oriented voxelgrid fusion for accurate 3D semantic mapping［J］. IEEE Transactions on Circuits and Systems for Video Technology，2022，32（1）：183-197．

[17] 尹佳琪，王世勇，李范鸣 ．基于改进主成分分析的分焦平面偏振图像去噪算法［J］．光学学报，2021，41（7）： 64-73．

[18] YIN J Q， WANG S Y， LI F M ．Division-of-focal-plane polarization image denoising algorithm based on improved principal component analysis［J］．Acta Optica Sinica，2021， 41（7）： 64-73．（in Chinese）

[19] 冯林，李斌兵 ．一种基于最小广义方差估计的TLS点云抗差法向量求解方法［J］．武汉大学学报（信息科学版），2018， 43（11）： 1647-1653．

[20] FENG L， LI B B ．A robust normal estimation method for terrestrial laser scanning point cloud based on minimum covariance determinant［J］．Geomatics and Information Science of Wuhan University，2018，43（11）：1647-1653．（in Chinese）

[21] 白创，陈立，闫昱 ．一种基于PVDAC描述子的RGB-D三维点云配准方法［J］．湖南大学学报（自然科学版），2023，50（2）：95-101．

[22] BAI C， CHEN L， YAN Y ．A RGB-D 3D point cloud registration method based on PVDAC descriptor［J］．Journal of Hunan University （Natural Sciences），2023，50（2）：95-101．（in Chinese）

[23] QIU Z B， MA Y， FAN F，et al ．Improved DBSCAN for infrared cluster small target detection［J］．IEEE Geoscience and Remote Sensing Letters，2023，20：5511905．

[24] CHOWDHURY S， NA H L， CORDEIRO DE AMORIM R ．Feature weighting in DBSCAN using reverse nearest neighbours［J］.Pattern Recognition， 2023，137：109314．

[25] ROS F， GUILLAUME S， RIAD R，et al ．Detection of natural clusters via S-DBSCAN a self-tuning version of DBSCAN［J］．Knowledge-Based Systems，2022，241：108288．

[26] LIU Y T， ZHANG L， LI P J，et al ．Laser radar data registration algorithm based on DBSCAN clustering［J］．Electronics，2023， 12（6）：1373．

基金资助

国家重点研发计划资助项目(2018YFB1308603)

National Key Research and Development Program of China(2018YFB 1308603)

福建省高校产学合作项目(2022H6016)

Fujian Province Industry Academia Cooperation Project(2022H6016)

教育部产学合作协同育人项目(231003084265411)

Collaborative Education Project of the Ministry of Education(231003084265411)