基于接触式轮廓仪测量数据的工件形状自动标注方法
摘 要
接触式轮廓仪在测量工件或产品的轮廓时，由于存在探针沾污、扫描位置不准等问题，会对工件形状的准确标注带来影响，本文基于接触式轮廓仪的测量数据，研究工件形状的自动标志方法。
针对问题一，本文首先采用一种基于一阶二阶导数的模式识别方法，以此作为理论基础区分水平线、斜线和圆。针对测试数据中的噪声问题，采用一种基于滑动平均的去噪方法。并比较了有重叠滑动窗口和无重叠滑动窗口的去噪效果。通过以上方法可以确定不同曲线段的形状，然后采用一种基于直线和圆的方程的拟合方法，对不同曲线段数据进行拟合，求得各个分段函数方程。最后计算得出直线夹角、半径等参数。与原始数据图形对比，基本吻合。
针对问题二，本文采用一种基于直线旋转量的图形校正方法。水平平移不改变直线的斜率，工件的倾斜旋转的角度就是直线夹角的旋转量。首先采用问题一的方法识别出直线的图形区域，再计算倾斜后的直线旋转量，并基于此进行图形校正。校正后的数据与问题一原始数据对比误差较小。
针对问题三，本文采用一种基于多组测量数据复原完整轮廓的方法。首先采用问题一方法识别出每组数据不同曲线段的形状。然后基于问题二水平直线的旋转量对各组数据进行校正，最后综合各组数据的公共部分和左右两侧的最大冗余数据进行拼接，复原完整轮廓。
针对问题四，本文采用一种基于多组局部数据修正全局图像的方法。首先采用问题一方法识别出每组局部数据不同曲线段的形状。然后通过拟合计算不同曲线段形状的方程和参数，最后等效替代全局图形中对应部分的方程和参数，实现基于局部数据对全局图形的修正。
模型的评价。以上工件自动标注方法，可以较好的替代人工计算，很大程度的解决了接触式轮廓仪工作中存在的问题，标注误差较小。
关键词：数据处理 拟合分析 去噪 线性回归 最小二乘法
1 问题重述
1.1 问题背景:
接触式轮廓仪的工作原理是，探针接触到被测工件表面并匀速滑行，传感器感受到被测表面的几何变化，在X和Z方向分别采样，并转换成电信号。
在理想状况下，轮廓曲线是光滑的，但由于接触式轮廓仪存在探针玷污、探针缺陷、扫描位置不准等问题，检测到的轮廓曲线呈现出粗糙不平的情况。
1.2具体问题：
假设被测工件的轮廓线是由直线和圆弧构成的平面曲线，建立数学模型，根据附件所提供的轮廓仪测量数据，研究下列问题：
1.2.1 问题一:
根据工件1在水平状态下的测量数据，标注出轮廓线的槽口宽度 (如x₁，x₃等)、圆弧半径(如R₁，R₂等) 、圆心之间的距离(如c₁，c₂等) 、圆弧的长度、水平线段的长度(如x₂，x₄等)、斜线线段的长度、斜线与水平线之间的夹角(如∠1,∠2等) 和人字形线的高度 (z₁) 。
1.2.2问题二:
计算工件1在测量时的倾斜角度，并作水平校正，比较理想下的水平测量与倾斜测量状态数据之间的差异。
1.2.2问题三:
根据工件2的10次测量数据完成：(1) 每次测量工件2的倾斜角度；(2) 标注出工件2轮廓线的各项参数值； (3) 画出工件2的完整轮廓线。
1.2.2问题四:
利用工件2关于圆和角的9次局部测量数据，修正问题3的结论，并对工件的完整轮廓线作进一步修正。
2 问题分析
针对问题一，根据水平放置的工件1的测量数据，计算工件自身轮廓线的参数，已知轮廓线是由直线和圆构成，那么问题的关键是我们用什么指标和区分圆和直线，注意到直线的二阶导数为零，而圆的二阶导数不为零，用二阶导数区分
圆和直线，而斜线和水平线可以用一阶导数来区分。利用区分好的数据进行圆和直线的拟合从而得到圆和直线的函数表示，再得到整个轮廓线的分段函数表达，进而可以得到各种参数。
针对问题二，已知倾斜放置工件1的测量数据和水平放置工件1的测量数据，计算倾斜角度和水平平移量。首先利用前面第一问的方法区分出来直线和圆，把水平测量和倾斜测量的多个圆的位置做对比，得到多个圆的平移量，取平均得到整体轮廓线的平移量； 注意到工件1的水平测量数据中有很多的水平线，倾斜之后变成斜线，计算多条斜线的倾斜角取平均得到工件1的倾斜角度； 最后利用倾斜角度和平移量通过矩阵相乘把倾斜放置的工件2的测量数据计算得到水平放置的测量数据，然后利用问题1的计算参数方法计算得到工件1水平校正后轮廓线的参数，最后和工件1真实水平测量的标注参数进行对比。
针对问题三，题目中给出对工件2的10次不同的测量数据，首先计算每次测量工件2的倾斜角度，我们假设工件2中本身有很多的水平线，仍然使用问题一的方法区分出直线和圆，那么斜率差不多相等的直线段比较多，然后倾斜角的平均，就为此次测量的倾斜角度； 其次，因为多次测量位置不同，每次测量的整个工件2的部分，所以我们把每次的测量数据乘以旋转矩阵水平校正后，通过对比，得到整体工件2水平校正后的数据，然后利用问题1 的方法标注参数； 最后，通过水平校正后的数据画出工件2的完整的轮廓线。
针对问题四，附件中给出了9次角和圆的9次局部测量数据，修正工件2的参数，并给出修正后的完整轮廓线。
3 模型假设
1) 假设被测工件的轮廓线是由直线和圆弧构成的平面曲线；
2) 假设接触式轮廓仪检测准确没有误差；
3) 假设工件检测的轮廓线无外界干扰。
4 符号说明
| 符号 | 说明 |
| xi | 第 i 个测量点在x轴上的坐标 |
| zi | 第 i 个测量点距离水平线的高度 |
| ki | 第i条直线的斜率 |
| ki' | 旋转后的第i条直线的斜率 |
| b | 直线在 z轴上的截距 |
| R | 圆弧的半径 |
| f₁(x) | 工件1轮廓线函数 |
| f₂(x) | 工件2轮廓线函数 |
5 模型建立与求解
5.1问题一模型建立和求解
根据问题限定条件，工件只包含水平直线、斜线、圆三种图形。我们要识别这三种模式，需要讨论三种图形数学上的本质区别是什么，为此本文首先提出以下模型：
5.1.1基于一阶二阶导数的模式识别方法：
建立x，z直角坐标系，设x轴的方向为水平方向，z轴的方向为垂直方向。设测量水平样本数据为( xizi,i=1⋯n,n
为数据个数，设工件1轮廓线的函数为f₁(x), 其中. f₁xᵢ=zᵢ;
若轮廓线为水平直线，则满足下列条件：
f₁'xᵢ=0
(1-1)
若轮廓线为斜线：则满足以下条件：
f₁'xᵢ=C
(1-2)
若轮廓线为圆，则满足以下条件
f₁''xᵢ≠0
(1-3)
其中， f1'xi=zi+1-zixi+1-xi,f1''xi=f'xi+1-f1xixi+1-xi
综合以上分析，本文提出可以基于一阶导数和二阶导数是否为0来综合判断该曲线的图形模式。如下图所示，为问题一数据的一阶导数和二阶导数图形：
从以上图形可以看出，原图形直线部分对应的一阶导数并不为0，这是因为测量数据中包含了很多噪声。实际应用上理论方法还需要解决去噪问题。为此本文提出以下模型：
5.1.2基于滑动平均的去噪方法
针对离散数据的一阶导数为
f1'xi=zi+1-zixi+1-xi
(1-4)
离散数据的噪声会导致以上公式的一阶导数剧烈波动，为消除噪声影响，本文提出基于滑动窗口内求平均的去噪方法，根据滑动窗口之间是否重叠，本文对比研究了有重叠区域的滑动窗口法和无重叠区域的滑动窗口法。
(1) 有重叠区域的滑动窗口法
选择滑动窗口长度为N，起始横坐标x₁，终止横坐标为. Xi+N。
有重叠区域的滑动窗口法，第一个窗口覆盖x₁到xN，第二个窗口覆盖x₂到 xN+1,
窗口内求平均，以此类推，则滑动窗口内的平均一阶导数记为 f1'xi,
f1'xi=1N∑ii+Nf1'xi
(1-6)
同理，滑动窗口内的平均二阶导数记为 f1"xi,
f1nxi=1N∑ii+Nf1"xi
(1-7)
(2) 无重叠区域的滑动窗口法
无重叠区域的滑动窗口法，第一个窗口覆盖x₁到xN，第二个窗口覆盖 xN+1
到x₂N，窗口内求平均，以此类推，则滑动窗口内的平均一阶导数记为 f1'xi,
f1'xi=1N∑i-1N+1iNf1'xi
(1-6)
同理，滑动窗口内的平均二阶导数记为 f1"xi,
f1nxi=1N∑i-1N+1iNf1"xi
(1-7)
考虑到滑窗法并不能将实际误差完全降低到0，选择阈值λ₁和λ₂分别表示一阶导数阈值和二阶导数阈值，并即滑窗内区域的图形模式为Y，平行直线、斜线、圆分别使用1、2、3表示。则Y满足以下条件：
lf₁" f1"xi≤λ2)
(1-8)
应用以上模式可以求得图形中哪些数据属于水平线，哪些属于斜线，哪些属于圆。接下来由于分段的带噪声的数据分别求直线、斜线、圆的方程还需要进行拟合，即以下模型：
5.1.3基于圆和直线方程的拟合方法[1]
(1) 圆的拟合:
设拟合曲线为 x-A²+z-B²=R²,
A，B，R是拟合参数
变形可得到圆的拟合曲线的另外形式
所以 x²+z²+ax+bz+c=0
设截取的关于圆的测量数据为 xizi,i=1⋯N
达到 minxi2+zi2+axi+bzi+c2
求得到a,b,c
得到圆心 A=a-2,B=a-2,
半径 R=12a2+b2-4c
(2) 直线的拟合：
设拟合曲线为 z=kx+b
其中k，b是拟合的参数
优化目标： minkxᵢ+b-zᵢ²
得到 k=n∑i=1nxizi-∑i=1nxininin∑i=1nxi2-∑i=1nxi
b=n∑i=1nxi2∑i=1nzi-∑i=1nxizi∑i=1nxin∑i=1nxi2-∑i=1nxin
根据该模型，带入实测数据，可以求得各个分段图像的方程表达式，进而求得夹角等参数。
5.1.4问题一模型的求解：
(1) 去噪:
分别使用有重叠滑窗和无重叠滑窗方法对一阶导数和二阶导数进行去噪，如下图所示。
通过上图对比，原始工件的第2个图形竖直向下斜线部分，在有重叠滑窗条件下，几乎不能识别出来。而图4中无重叠滑窗部分存在一阶导数基本不变的短暂区域，可以大体识别出来。因此下文主要采用无滑窗的去燥方法。
原始数据第一个曲线交点为下图标注点，查原始数据该点为第6052个数据。
如上图是采用滑动平均后的去燥一阶导数数据，该图的第一个拐点为61，考虑到滑窗长度设置为了100，该点对应第6100 个数据。
以上对比表明，该方法可以识别出不同曲线段的转折点。下图左是抓取的水平测量数据递减直线的位置，下图右是抓取的倾斜位置测量数据递减直线的位置，效果很好。(其中，一个点表示1000 个数据)
(2) 拟合结果:
下面是水平放置的轮廓线的第一条和第四条和第一个圆的拟合情况，其他的见附录。
(3) 工件1的参数标注：
轮廓线的函数表达：
1.77 46.595892371385 <x≤x₁
-2.723867741639066x+133.8002234310488 x₁<x≤x₂
-0.49172-x-51.21732-4.2610
x₂<x≤x₃
2.754359108225884x-146.8171658730294 x₃<x≤x₄
1.77 x₄<x≤x₅
-3.724845609376202x+213.0069757236736 x₅<x≤x₆
-0.52592-x-58.60192-4.1492
x₆<x≤x₇
2.980598439534949x-179.8507093519987 x₇<x≤x₈
1.77 x₈<x≤x₉
-3.723602901285295x+239.2514907860618 x₉<x≤x₁₀
-0.38022-x-65.73382-4.3746
x₁₀<x≤x₁₁
3.511651703755757x-236.2555247229892 x₁₁<x≤x₁₂
f(x)= 1.77 x₁₂<x≤x₁₃
0.194764774178566x-15.745625965295059 x₁₃<x≤x₁₄
-0.200744158613707x+14.645116449010247 x₁₄<x≤x₁₅
1.77 x₁₅<x≤x₁₆
-0.91122-x-85.72382-1.4363
x₁₆<x≤x₁₇
1.77 x₁₇<x≤x₁₈
0.52992-x-89.76852-2.3050
x₁₈<x≤x₁₉
1.77 x₁₉<x≤x₂₀
-3.99232-x-98.04052-0.0869
x₂₀<x≤x₂₁
-MWWE--------------- x₂₁<x≤x₂₂
x₂₂<x≤x₂₃
x₂₃<x≤118.123985877831
把拟合得到的函数与x轴得到相交求交点，即对函数求零点得到所有槽口宽度x₂，x₄，水平线段的长度x₁，x₃，通过圆心的参数得到圆心之间的距离c₁，c₃，最后两条斜线的交点求得人字形的高度z₁。
表 1 轮廓线的槽口宽度及圆弧的长度、水平线段的长度
| x₁ | x₂ | x₃ | x₄ | x₅ | x₀ | x₇ |
| 2.8914075 | 4.9982367 | 2.0873795 | 4.9800147 | 2.0468106 | 4.9940161 | 9.9916311 |
| 13568545 | 01987779 | 02859584 | 21154198 | 48264795 | 36801619 | 60546746 |
| x₈ | x₉ | x₁₀ | x₁₁ | x₁₂ | x₁₃ | |
| 2.6672212 | 0.9356328 | 1.1740463 | 12.268778 4.6218374 10.153688 | |||
| 42061245 | 84535536 | 57141904 44708819551824376 233493440 | ||||
| 表 2 轮廓线的斜线与水平线之间的夹角 | ||||||
| ∠1 | ∠2 | ∠3 | ∠4 | |||
| 110.1594348872009 | 109.9539638548886 | 105.0276884989337 | 108.5467473758782 | |||
| ∠5 | ∠6 | ∠7 | ∠8 | |||
| 105.0324768433333 | 105.8951506831928 | 168.9787730174710 | 168.6490737801 013 | |||
5.3.1直线旋转量的图形校正方法[3]
首先利用前面第一问的二阶导数是否为零识别圆和直线，把所有的直线数据提取出来，分别通过拟合计算其斜率，斜率相等的直线所在的直线的斜率即是工件2的倾斜斜率，并基于此进行图形校正。
1.旋转角的确定：
设k,,i=1…l为第i条直线的斜率，k₁',i=1…m其中差不多相等那些直线中的第 i条，所以所在直线斜率
k=1m∑ki'=-0.131
所以倾斜角为
θ= arctank
利用前面的直线拟合算法, 用 matlab求解得到θ=-0.13=-7.4°.
2.水平校正：
水平校正分两步：
(1) 倾斜测量数据的旋转
把倾斜的数据旋转：
xizi0=cosθsinθ-sinθcosθxi,i=1⋯n
其中 xizi0,i=1⋯n
为工件1倾斜的数据水平旋转之后得到的数据， xi'zi,i=1⋯n
为倾斜放置的工件1 测量的原始数据。
(2) 比较平移量
计算每一个点的平移量： xi0-xi,i=1⋯n
取平均得到整个轮廓的水平平移量 a=1n∑i=1nxi0-xi
通过几个特别的点，比如圆心验证。
(3) 检验:
根据上面求得旋转角度为θ，水平偏移量为a，
mximzi=cosθ-sinθasinθcosθ0001xizi1,i=1⋯n
xizi,i=1⋯n
为水平放置的工件1 测量的原始数据，
mximzi,i=1⋯n
为水平放置的旋转角度得到的数据。
通过误差平方和的平均检验上面得到的倾斜角度和水平偏移量的准确性，即
1n∑i=1nxi-mxi2+zi-mzi2
从旋转后的和水平放置的原始数据对比图看出效果很好：
3. 对校正后的数据进行标注参数：
利用问题一模型里的三步：第一，识别圆和直线，第二，拟合直线和圆，第三，得到路廓线函数表达，计算参数。
5.4问题三的模型建立和求解：
本问题考虑的是在对工件作多次检测过程中，工件每次放置的角度、测量的起点和终点都会有偏差，此时如何根据已测得的多组数据复原工件完整轮廓的问题。
5.4.1基于多组测量数据复原完整轮廓的方法
1.每次测量的倾斜角度：
我们假设工件2中水平放置时本身有很多的水平线，仍然使用问题一的方法区分出直线和圆，那么斜率差不多相等的直线段比较多，那么这倾斜角的平均为此次测量的倾斜角度，所以类似问题二中求解倾斜角的模型得到每次测量的倾斜角。
2.对比得到工件2轮廓的全部水平测量数据
(1) 先把每一次的数据通过进行旋转
设为第j次的测量数据(x₁³, z₁'),i=1…n, j=1,2…10,
第j次旋转后的测量数据 xij0zij0,i=1⋯n,j=1,2⋯10.
(2) 整合为一组数据
因为多次测量位置不同，每次测量的整个工件2的部分，对比第j次旋转后的测量数据 xij0zij0,i=1⋯n,j=1,2⋯10.
得到工件2的水平校正后的整体数据(x₁⁰, z₁⁰),i=1…n,
5.4.2基于多组测量数据复原完整轮廓的方法具体求解步骤
步骤1：基于在问题1中建立起来的基于一阶二阶导数的模式识别方法和基于滑动平均的去噪方法对附件2中给出的10组测量数据进行直线或圆的图形类型的识别，并抓取每个图形所在的位置信息。
步骤2：基于步骤1获得的图形类型及其位置信息，挖掘每组数据的图形结构信息。具体地，首先确定每组中的直线有几条，圆有几个； 其次，根据位置信息将这些直线和圆的前后位置关系搞清楚，从而对所测工件的结构有了框架性的认识。
步骤3：基于每条直线的位置信息，抓取附件2中直线的原始数据。然后通过拟
合的方法，得到每条直线的方程。
步骤4：根据直线的斜率相等与否，对每组数据中的直线进行分类。把每组中直线数量最多的类作为目标类。将目标类中直线的斜率作为本组数据的旋转角度θ。
步骤5：基于问题2中建立的直线旋转量的的图像校正方法，以旋转角度θ对每组数据进行图像矫正。
步骤6：对每组数据矫正后的图像进行比较，选用每组数据中从左边起第一个圆的位置作为基准，对每组数据进行水平方法的平移矫正。
步骤7：对步骤6得到的矫正后的数据以横坐标范围的大小对数据进行同一整合，复原工件的轮廓。
依据问题1建立的方法，对复原之后的工件2进行测量，数据如下：
圆的圆心坐标: (37.8866, -32.5784) 半径: 5.7886
水平直线(从左边看起) 的长度：
第一条: 3.8256; 第二条: 3.3156; 第三条: 0.8752; 第四条: 1.1463
垂直直线 (从上边看起) 的长度：
第一条: 2.7381; 第二条: 1.9542
斜线(从上边看起) 的长度：
第一条: 9.0421; 第二条: 0.3862
5.5问题四的模型建立和求解
该问提供了多组局部的测试数据，包含了圆和角的多次测量。局部的测试数据相比第三问的全局测试数据更加精细，目的是为了解决全局测试时个别局部区域标注偏差较大的问题。基于此，本文采用一种基于多组局部数据修正全局图像的方法。
5.5.1基于多组局部数据修正全局图像的方法
该方法分为以下几步：
步骤一：使用问题一的方法识别出每组局部数据的图形模式，即它是圆的一部分还是直线。
步骤二：使用问题二的方法，对该组局部数据进行校正。
步骤三：使用校正后的局部数据与第三问全局数据进行等效匹配，从而确认校正后的局部数据是全局数据中的某一区域的放大。
步骤四：根据校正后的数据计算各个参数，并修正全局图形
步骤五：对其他各组局部数据重复以上步骤，得到修改后的完整工件图形。
按照以上步骤，以第一组圆形数据为例进行处理，其余数据的处理结果见附录。
如下图所示，是第一组圆的原始局部数据图形：
然后使用问题一和问题二中的算法进行校正，得到下图：
代码部分
clc, clear;
close all;
%导入数据
gj1= xlsread('附件1 工件1 的测量数据');
x gj1=gj1(:,1);z gj1=gj1(:,2);
plot(x gj1,z gj1)
%设置窗口大小及存放一阶导数、二阶导数、均值、方差等的数组
size chuangkou=1000;
p= round( length(x gj1)/ size chuangkou)-2;
m= ones(p,1);
n= ones(p,1);
x new= ones(p,1);
z new= ones(p,1);
zhi m= ones(p,1);
zhi n= ones(p,1);
zhi x new= ones(p,1);
%求一阶导数k1、二阶导数k2以及根据窗口去噪的方法求窗口均值和方差 for d=1:p
v=(d-1)* size chuangkou;
k1= ones( size chuangkou);
for i=1: size chuangkou
k1(i)=(z gj1(i+v+1)-z gj1(i+v))/(x gj1(i+v+1)-x gj1(i+v));
end
k2= ones( size chuangkou);
for i=1: size chuangkou
k2(i)=(k1(i+1)-k1(i))/(x gj1(i+v+1)-x gj1(i+v));
end
m(d)= mean(k2(1: size chuangkou));
n(d)= var(k2(1: size chuangkou));
x new(d)=x gj1(d* size chuangkou+1);
z new(d)=z gj1(d* size chuangkou+1);
zhi m(d)= mean(k1(1: size chuangkou));
zhi n(d)= var(k1(1: size chuangkou));
zhi x new(d)=x gj1(d* size chuangkou+1);
end
% plot( zhi x new(1:p), zhi m(1:p),’.’)
% hold on
%设置二阶导数的窗口均值的阈值，抓取圆域所在的位置
x circle= ones(p,1);
for i=1:p
if m(i)>3.5
x circle(i)=x new(i);
end
end
% plot(x new(1:p),x circle(1:p),'.')
% hold on
%通过抓取到的圆域的位置，抓取圆域的原始数据以备拟合圆的方程之用
j=0;
x nihe= ones(20000,1);
z nihe= ones(20000,1);
for i=1:p
if x circle(i)>2
for k=1: size chuangkou
x nihe(j+k)=x gj1(i* size chuangkou+k);
z nihe(j+k)=z gj1(i* size chuangkou+k);
end
j=j+ size chuangkou;
end
end
%圆的拟合计算圆心，半径
% xx= xlsread('新建 Microsoft Excel 工作表');
x=x gj1(10001:13000);y=z gj1(10001:13000);
N= length(x);
C=N* sum(x.^2)- sum(x)* sum(x);
D=N* sum(x.*y)- sum(x)* sum(y);
E=N* sum(x.^3)+N* sum(x.*y.^2)- sum(x.^2+y.^2)* sum(x);
G=N* sum(y.^2)- sum(y)* sum(y);
H=N* sum(x.^2.*y)+N* sum(y.^3)- sum(x.^2+y.^2)* sum(y);
a=(H*D-E*G)/(C*G-D^2);
b=(H*C-E*D)/(D^2-G*C);
c=-sumx.-2+y.⇀2+a*sumx+b*sumy/N;
A=a/(-2);
B=b/(-2);
R=1/2*sqrtar2+b∨2-4*c;
%设置一阶导数窗口均值的阈值，抓取直线所在的位置
x zhi= ones(p);
for i=1:p
if zhi m(i)>0.15
x zhi(i)= zhi x new(i);
end
if x circle(i)>2
x zhi(i)=1;
end
end
plot( zhi x new(1:p),x zhi(1:p),'.')
hold on
title(’ level水平直线位置(窗口: 1000; ) ’)
%通过抓取到的直线的位置，抓取直线的原始数据以备拟合直线的方程之用s=0;
x nihe zhi= ones(20000,1);
z nihe zhi= ones(20000,1);
for i=1:p
if x zhi(i)>2
for k=1: size chuangkou
x nihe zhi(s+k)=x gj1(i* size chuangkou+k);
z nihe zhi(s+k)=z gj1(i* size chuangkou+k);
end
s=s+ size chuangkou;
end
end
%直线的拟合求出斜率和截距
x=x gj1(22301:23500)’;y=z gj1(22301:23500)’
扫一扫
579
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
