基于三维点云螺纹中轴线的拟合方法研究

谢张宁， 姜志华， 朱惠臣， 李智玮， 雷李华,，孔 明， 王道档， 张 波

(1. 中国计量大学计量测试工程学院，浙江杭州310018；2. 上海计量测试技术研究院，上海201203)

1 引 言

高精度螺纹规是精密机械、航天航空、石油化工、核电等领域中极为重要的基础零部件，其几何精度对装配质量具有决定性影响[1]。现有的高精度测量方法众多，但均是基于一维、二维的螺纹测量技术，无法获取三维综合参数，检测结果难以反映实际精度，导致装配质量不能准确控制。要改变“粗糙螺纹”生产的现状，需要对三维螺纹测量仪器设备进行研究开发。

当前的三维测量方法在实际测量前，首先要实现机器坐标系MCS与工件坐标系WCS转换[2]，在螺纹面的测量中，螺纹中轴线的准确测定对螺纹的后续测量至关重要。然而，在对螺纹规测量时，仍旧是依靠测量人员的经验,手动调整螺纹规位置，然后进行后续测量。这种测量方法不仅增加了测量难度，也增大了测量结果的不确定度。

本文研究了一种能够快速精准拟合螺纹中轴线的方法。在测量时，螺纹规放置在工作台上进行调整，但是工件坐标系不可能与机器坐标系完全重合，所以要使用螺纹规的基准来建立工件坐标系。将测量仪器对准螺纹规上下底面的顶针孔实现定位，在机器坐标系MCS下，非接触扫描得到的螺纹规外表面点云;经过随机增量算法计算点云凸包，制成螺纹凸包点集;基于最小二乘的概念建立螺纹中轴线的拟合数学模型，得出中轴线的单位向量v(a,b,c);再基于单位向量得到螺纹底面中心A，坐标旋转矩阵R，坐标偏移矩阵S;螺纹测量仪器根据上述参数得出工件坐标系。

由于螺纹复杂的三维形体结构，在建立数学模型对螺纹面进行直接计算时，牙槽表面、螺纹表面缺陷点云数据均会影响螺纹中轴线的拟合，导致结果不准确。为防止干扰，本文用随机增量算法计算螺纹表面三维点云凸包(convex hull)[3]；将螺纹牙槽表面、螺纹表面缺陷点云数据筛去，对于空间离散点集，求解最小的凸多面体，其中把空间中点集所有的点包含在里面的凸多面体就是凸包[4]。基于三维点云凸包建立中轴线拟合数学模型，能有效滤除掉螺纹牙槽、硬件缺陷、扫描误差等形成的干扰点云，减少了后续算法的运算量，加快了运行时间，提升了运算结果。

经过计算凸包的螺纹表面点云实质上是以公称直径为上下底直径的圆柱面。要拟合该圆柱面的轴线就要以面上的点与轴线的位置关系建立数学模型，并采用优化算法求解。常用的优化算法有最小二乘法、爬山算法、模拟退火算法、遗传算法等[5～8]。用爬山算法、模拟退火算法、遗传算法等随机优化算法确定唯一解的方法虽然可行，但随机优化算法大多处理高阶次线性问题；目前通用、认可度高的算法仍为最小二乘法，在工程测量上，用最小二乘法解决超定方程[9]、线性拟合[10～12]等数据量大的低阶次线性问题速度快、结果准确；且随机算法优化过程具有随机性，容易出现早熟或过拟合等问题。因此，拟合螺纹中轴线采用最小二次算法具有较强的实用性。

2 螺纹中轴线数学模型

在空间机器坐标系中，螺纹中轴线l可以表示为经过点A(x0,y0,z0)，单位向量为v(a,b,c)的直线l，其点斜式为[13]：

根据最小二乘原理，使得经过三维凸包计算的螺纹表面点云P(xi,yi,zi)，i=1,…,N到拟合直线距离d的平方和最小：

(1)

(2)

假设：

(3)

将式(2)、式(3)代入式(1)，可得最小二乘意义上的螺纹中轴线数学模型：

(bδz-cδy)2]+f(a,b,c)

(4)

由式(4)可以看出，对于任意的a,b,c均有：

(5)

当δx=δy=δz=0时，有：

min{f(a,b,c)}
a2+b2+c2=1

其中：

f(a,b,c)=(1-a2)B11+(1-b2)B22+

(1-c2)B33-2abB12-2acB13-2bcB23

(6)

3 基于三维点云的最小二乘拟合

3.1 制成三维凸包点集

由于螺纹表面三维结构复杂以及可能存在缺陷的原因，要对点云数据计算凸包。计算凸包是根据算法在点云中筛选出能构造包含所有点的最小多面体的点集。

螺纹表面点云是三维数据点，且数据量庞大，用随机增量算法计算凸包效果好、速度最快[15]。计算流程为：

(1)首先用筛选函数找到x,y,z轴上的边界点，如x,y,z值最小的点，设其为P1。

(2)随机找到与P1不共线的2个点P2,P3，再找到不共面的第4个点P4，构成四面体。

(3)若新点Pi在凸包内部或边界，则忽略。若新点Pi在凸包外部，见图1(a)，就要计算原凸包相对于Pi的面，见图1(b)所示;将Pi点能看到的所有面删除，同时将边上的点与Pi相连组成新的凸包，见图1(c)。

(4)依次遍历余下的点，制成三维凸包点集。

图1 随机增量算法计算凸包Fig.1 Random incremental algorithm to calculate convex hull

3.2 单位向量的求解

将式(6)引入矩阵B=[Bij]，整理成矩阵式，则可得：

由于B是对称方阵，不难看出，当[a,b,c]T是矩阵B的最大特征值λmax(B)对应的单位特征向量时，上式可以取得最小[16]。此时，[a,b,c]T为螺纹中轴线单位方向向量：

3.3 底面中心点的求解

在转换后的坐标系o-x′y′z′中，取圆环最外层，设其拟合曲线为：

式中M′为螺纹外径拟合结果。将上式展开，令：

整理得：

最小二乘法求解上述超定方程解得[18]：

其中：

3.4 测量坐标系到工件坐标系的转换

测量坐标系到工件坐标系的转换分两步：

(1)先将测量坐标系进行原点偏移。即将A(x0,y0,z0)作为工件坐标系的新原点o′，只需在原坐标系各轴减去相应的偏移量，可得偏移矩阵S：

(2)将坐标轴旋转，使z轴与对中轴线单位向量重合。坐标旋转矩阵R在求解中轴线底面中心时已经求出。

4 仿真实验

4.1 有效性验证

用MATLAB软件编写确定螺纹中轴线的程序。用MATLAB自带的[v,λ]=eig(X)函数求解特征向量λ对应矩阵的特征值，v中3个列向量为λ对应的特征向量。求解单位向量程序框图见图2。

图2 求解单位向量程序框图Fig.2 Block diagram for solving unit vector

用MATLAB编程求解底面中心，将螺纹点云对应的三个轴数据向量经过矩阵R后得到处于o-x′y′z′坐标系的点云数据，再由MATLAB自带的X-1=inv(X)函数进行坐标反变换，得到原坐标系的底面中心数据。

图3 求解底面中心坐标程序框图Fig.3 Block diagram of the bottom center coordinates

用圆柱点云图检验本文中轴线拟合算法的有效性。已知该圆柱点云中轴线参数为：

点云数据由scatter3函数绘制，在计算完成的中轴线上取若干点用MATLAB自带函数plot3绘制，如图4所示。

图4 基于三维凸包的最小二乘拟合结果Fig.4 Results of least square fitting based on 3D convex hull

中轴线拟合算法解得单位向量vt：

圆柱底面中心点At：

坐标旋转矩阵R：

坐标偏移矩阵S：

在有限值域区间内，以一条直线上的点到另一条直线的距离的最大值作为指针，定量评定不同直线与已知直线之间的符合程度。本文以两直线在螺纹旋合长度内，某一直线两端点到另一条直线的距离的最大者为评定参数，评价两条直线的符合程度。其算法过程为：

(1)首先以任意一条直线l的单位向量为法向量，取点云边界点Pα生成上截面α；下截面法向量相同。

(2)将上截面方程转换为参数方程，将参数方程代入两直线lt、lz的点斜式，可解到拟合直线lt在上截面α的端点坐标(xlt，α,ylt，α,zlt，α)，已知直线lz在上截面α的端点坐标(xlz，α,ylz，α,zlz，α)。

(3)同理按步骤(1)生成下截面β，拟合直线lt在下界面β的端点坐标(xlt，β,ylt，β,zlt，β)，已知直线lz在下截面β的端点坐标(xlz，β,ylz，β,zlz，β)。

(4)分别计算截面α,β上，两直线端点距离Eα、Eβ，用MATLAB自带的Max函数取两者的较大值即为评定参数E。

运算所得E=0.21 μm，所得中轴线参数与lz直线基本重合，算法可行。

4.2 中轴线拟合仿真实验

为了比较研究中轴线拟合算法的可行性，本文用MATLAB软件编译基于三维凸包的最小二乘算法拟合螺纹中轴线，并与投影法比较验证其可行性。根据国家标准《普通螺纹基本尺寸》表1[19]里的螺纹小径、中径数据建立了公称直径为80 mm,螺距为4 mm的普通外螺纹面的三维模型，并在模型表面随机成块删除点云形成螺纹自身表面缺陷、扫描失真等造成的干扰因素，如图5所示。

图5 带有干扰因素的外螺纹三维模型面Fig.5 Three-dimensional model surface of external thread with interference factors

将三维模型面保存为obj格式，在Visual Studio 2017环境下用PCL18.1编程库将其生成外螺纹面三维点云。外螺纹面底面圆心坐标为(-40,150,-40)，顶端圆心坐标为(-40,0,-40)，如图6所示。

图6 带有干扰的外螺纹三维点云Fig.6 3D point cloud image of external thread with interference

4.2.1 投影法拟合中轴线

投影法是轴体测量中常见的方法,在测量时，其将螺纹器具投影至xOy面、yOz面和zOx面，基于平面测量[20]。

本文将三维点云投影至xOy面、yOz面经过计算边界坐标，将xOy面、yOz面的边界点y值取中值，分别拟合直线。求得投影面上的中线，过所得中线作面，根据2个不平行的平面相交确定1条直线的原则[13]确定空间螺纹中轴线，用MATLAB编程得中轴线一般方程：

4.2.2 最小二乘拟合螺纹中轴线

用随机增量算法对带有干扰因素的外螺纹三维点云作convex hull计算，计算结果如图7所示。

图7 经过convex hull处理后外螺纹三维点云Fig.7 Thread 3D point cloud after convex hull processed

经过convex hull处理后的外螺纹三维点云导入MATLAB，将参数区间限制为点云数值区域内，用投影法和基于三维凸包的最小二乘算法拟合中轴线。图8、图9分别为螺纹中轴线三维拟合结果以及其在xOy面、yOz面的投影图。

图8 基于三维凸包的最小二乘算法和投影法拟合结果Fig.8 Fitting results of least squares algorithm and projection method based on three-dimensional convex hull

图9 拟合结果的投影Fig.9 Projection of the fitting result

用基于三维凸包的最小二乘算法拟合的螺纹中轴线所确定的中轴线lls参数为：

单位向量vls

螺纹底面中心Als

坐标旋转矩阵R

坐标偏移矩阵S

4.3 结果验证

由图8、图9可知，投影法拟合的直线与基于三维凸包最小二乘法拟合的直线基本重合。为了验证基于三维凸包的最小二乘算法拟合中轴线的精度和符合程度，本文采用计算lty,lls与外螺纹点云距离的方差[21]来表示算法拟合精度：

符合程度评价标准与上节相同，在有限值域区间内，将1条线段的2个端点到另1条直线的距离的最大者作为评定标准E。以max函数、min函数对点云的x，y，z值取最大或最小，找到上下边界点Pb1，Pb2，再按照第4.1节E算法过程求解。Xdata，Ydata，Zdata分别为点云的x，y，z数据集，上、下边界点的生成代码如下：

%生成上边界点

max_X=max(Xdata); max_Y=max(Xdata);

max_Z=max(Xdata);

P_a=[max_X,max_Y,max_Z];

%生成下边界点

min_X=min(Xdata); min_Y=min(Xdata);

min_Z=min(Xdata);

P_b=[min_X,min_Y,max_Z];

为比较投影法和最小二乘算法的计算速度，在处理器为Intel(R)Core(TM)i5-3230M的计算机上运行，结果见表1。

表1 拟合确定的中轴线比较表Tab.1 Table of fitted central axis variances

从拟合结果来看，用二维投影法和最小二乘法拟合的中轴线与外螺纹点云距离的方差结果都为0.34,且两直线在旋合长度范围内端面点距离的最大者，值为0.15 μm，达到了三维高精度测量标准，说明最小二乘法能够获得投影法一致的效果;在计算机执行速度上来看，最小二乘法拟合中轴线程序在MATLAB上运行时间只有15.3 ms，在同样的硬件条件下,逐次求解的投影法计算速度为1 196.7 ms。由此可以得出基于三维点云的最小二乘拟合算法是一种快速、准确的螺纹中轴线的拟合方法。

综上，最小二乘拟合中轴线的方法是从三维的角度上拟合中轴线，其囊括的维度更高更具有说服性，且计算时间短，在高精度螺纹规的测量上更具有实用性和推广价值。

5 结 论

本文从螺纹中轴线特征入手，建立了螺纹中轴线数学模型，提出了基于三维凸包的最小二乘算法拟合螺纹中轴线的方法，通过MATLAB软件编程验证算法的可行性，最后结合三维模拟仿真实验，并经二维投影法验证，得出以下结论：

(1)经过三维凸包计算的最小二乘算法在拟合具有三维复杂结构的螺纹中轴线具有与投影法一致的结果。

(2)经过三维凸包处理的最小二乘法从三维的角度解决问题，计算速度快，在高精度螺纹测量上，具有实用性和推广价值。