零件曲线拟合不确定度的评估与降低方法

韦进文

(广西大学机械工程学院,530004,南宁)

曲线拟合广泛应用于逆向工程、产品检测与加工质量控制等许多领域[1]。在应用传统的最小二乘法(LSM)拟合时,需要给出所用直线、圆弧等几何元素类型(以下简称元素)的目标函数,以及分割点(边界)处的约束条件,力求使元素尽可能逼近待拟合数据序列[2]。但是,数据序列所包含的元素种类、数量等通常是不可知的,也就无法预先给出分割点处的约束条件,甚至正确的目标函数。很多情况下,元素类型、数量及各自的数据范围是主观指定的,这显然不合理。在很多情况下,需要揭示数据对象的内在信息比单纯的数据逼近更有用[3]。所以,待拟合数据序列的辨识、分割是一种重要技术[4]。

以霍夫变换为代表的一类传统数据分割方法主要用于直线与圆等简单元素识别[5]。把数据序列的点映射到某元素参数空间,若多个映射交汇于一点,则对应数据点属于相同元素。这种用点映射“投票”来决定其所属元素的方法受噪声或误差影响时,在参数空间中可能由于“票数分散”引起识别误差[6],表明其鉴别能力有限。此外,霍夫变换用二维参数空间一个点来表示数据序列的参数恒定性,不适于其他高维参数元素的辨识与分割;简单元素如直线或圆弧的数据分割方法还包括遗传算法[7]、曲率分析[8]等。曲率分析涉及微分运算,受数据噪声影响较大;而遗传算法的分割精度不高。较复杂的其他元素如样条曲线,可以用模板匹配(广义霍夫变换)[9]、粒子群优化或杜鹃搜索算法[10]等分割方法,但精度有限;也有人采用曲线二阶导数符号变化来判断元素分割点[11],但受噪声影响很大。

拟合中的数据分割是一个随机事件。传统拟合只考虑减小数据噪声与误差的影响,使元素尽量逼近所拟合的数据序列,而对正确元素选用以及数据分割准确性对拟合的影响重视不够;分割方法一般存在适用特征维数低、受噪声影响大且分割点的随机偏移量大等不足。本文探讨数据分割误差与用错元素等因素对拟合的影响,据此提出一种几何元素辨识及其数据分割方法。

1 不确定度评估

根据JJF1059.1—2012,测量不确定度评估方法有指南法(GUM)与蒙特卡罗法(MCM)[12]等,其中GUM用于线性、对称或正态测量模型,而MCM用于非线性模型[13]。曲线拟合不确定度来源主要有:数据噪声与误差、数据分割误差以及用错元素(例如把圆弧误当做直线进行拟合)。数据噪声与误差在测量领域普遍存在,相应的不确定度评估已有成熟方法例如GUM[12];分割误差不仅由数据误差引起,还常常源于分割方法;而用错元素的原因与后果则更为复杂。因此,曲线拟合一般是一个强非线性模型。所以本文考虑MCM评估,以一个凸轮为对象,对其轮廓拟合实施MCM数字仿真实验。

1.1 实验对象设计

某凸轮轮廓曲线包含4个元素,名称与参数分别为:基圆弧b1与b2,圆心均为(0,0)、半径均为18;顶圆弧T,圆心为(20,0)、半径为5;右圆弧R,圆心为(-10.475 0,22.699 7)、半径为43;Bezier曲线段Bez,顶点为Q0(30.253 2,-8.957 3)、Q1(22.214 9,13.278 1)与Q2(-1.213 9,19.458),各参数的单位均为mm。圆弧R与曲线段Bez的两端分别关于x轴对称;各元素在连接点处满足C1连续;由于Q0与Q2均在轮廓外,曲线段Bez仅部分在轮廓上,两端分别对应伯恩斯坦基函数的参数t=0.3与t=0.8。

把该轮廓线从基圆弧的中部开始,绕原点进行以1°为步长的等角采样离散,可以得到360个采样点,分为5个子集,即圆弧b1:{P1,…,P115},圆弧R:{P116,…,P171},圆弧T:{P172,…,P187},曲线段Bez:{P188,…,P245},及圆弧b2:{P246,…,P360},即各元素数据段末端数据点下标为{ieb,ieR,ieT,ieB}={115,171,187,245}。

对采样数据施加均方差σ=0.1 μm的高斯噪声。

1.2 曲线拟合算法

1.2.1 圆弧拟合 待拟合的参数为半径r与圆心(xc,yc)。设数据序列(xi,yi),i=1,…,mc,由最小二乘法(LSM)得目标函数为

(1)

1.2.2 Bezier曲线拟合 二阶Bezier曲线具有足够曲率与平滑性,可以很好拟合一般的样条曲线。由控制点Q0～Q2生成的Bezier曲线中,将其参数方程中的参数消掉,得到抛物线形式的Bezier曲线

(Cyx-Cxy)2+[2(ByCx-BxCy)Cy+

Ay(AxCy-AyCx)]x+[2(BxCy-ByCx)Cx+

Ax(AyCx-AxCy)]y+(AxBy-AyBx)(AxCy-

AyCx)+(ByCx-BxCy)2=0

式中:A=2(Q1-Q0);B=Q0;C=(Q0+Q2-2Q1);x、y可互换。

设Cx≠0,Bezier曲线LSM目标函数的一般形式为

(2)

作为一般的圆锥曲线,其LSM目标函数则为

传统的样条曲线参数化(如累加弦长参数化CCL)拟合通常需要曲线的控制端点Q0或Q2,用抛物线拟合则不受此限制,可以很好地拟合曲线片段。

1.2.3 含多元素曲线的LSM拟合 把前述实验对象分割后的数据(xi,yi)按整体拟合,其LSM目标函数为

s.t.C1连续

在圆弧R两端R1与R2满足

(xb-xR)2+(yb-yR)2-(rb-rR)2=0,(xT-xR)2+(yT-yR)2-(rT-rR)2=0;在曲线段Bez两端(R1x,-R1y)与(R2x,-R2y)满足

而

Z={xb,yb,rb,xR,yR,rR,xT,yT,rT,b1,b2,b3,b4}

(3)

为待拟合的13个参数,可用诸如高斯-牛顿法解之。为了简单方便,以下用拟合误差表示相应的参数

δZ=

{δxb,δyb,δrb,δxR,δyR,δrR,δxT,δyT,δrT,δ1,δ2,δ3,δ4}

1.3 不确定度模型

1.3.1 数据分割误差造成的拟合不确定度 分割误差可视为随机变量。考虑到前述实验对象4个分割点关于x轴对称,分割点下标{ieb,ieR,ieT,ieB}的误差设为{d1,d2,-d2,-d1},其中{d1,d2}服从{μ1,μ2}=0,Σ=diag(σ1,σ2)的二维正态分布。各元素分割系依据其各自局部参数特征,其独立性是有保证的,所以d1与d2不相关,即

d1～N(0,σ12);d2～N(0,σ22)

(4)

在离散测量数据序列中,{d1,d2}应离散为整数,即把小数点后的数字4舍5入。考虑分割误差后,上述凸轮轮廓4个元素数据段的下标邻域依次变为

[1+ieB-d1,ieb+d1];[1+ieb+d1,ieR+d2]

[1+ieR+d2,ieT-d2];[1+ieT-d2,ieB-d1]

以下设定几组分割误差的组合{σ1,σ2},每对组合实施5×104次数字仿真实验,各组实验数据仍用前述的噪声水平。由于概率密度函数(PDF)与直方图的相似性,图1给出xb、yb、b1、b4等几个参数用直方图表示的实验结果,δzi(σ1,σ2)～F,F为实验出现频次。参数期望与不确定度见表1。

表1 各种分割误差的拟合期望与不确定度Table 1 Some fitting results under various segmentation errors

图1中,若分割误差{σ1,σ2}=0,拟合误差直方图在噪声作用下大致呈现正态分布的形态δZ～N(μ,Σ),几乎所有的拟合误差期望均在μm级以下,不难从中建立其PDF及各参数的置信区间。然而,随着分割误差的出现,拟合误差直方图呈现出复杂的多峰形态(d1,d2离散化导致的栏栅效应)。在MCM实验中发现,拟合误差状况与元素及其参数的对称性、数据段大小、拟合算法鲁棒性等有关。

(a)圆弧b段参数xb拟合误差随分割误差的直方图

(b)圆弧b段参数yb拟合误差随分割误差的直方图

(c)曲线段Bez参数a1拟合误差随分割误差的直方图

(d)曲线段Bez参数a4拟合误差随分割误差的直方图图1 部分变量拟合的MCM实验结果Fig.1 The results of some MCM experiments

(1)元素/参数对称性及数据段大小。圆弧b与圆弧T均关于x轴对称,而圆弧R与曲线段Bez既不对称也不互相对称。圆弧b与圆弧T的对称性会随着混入圆弧R与曲线段Bez的数据点错误而变差。因此,δyb的正态形状随着分割误差而逐渐变坏;而圆弧T则由于数据量过少,在分割误差很小时δyT即出现栏栅效应。

(2)算法鲁棒性。式(3)中曲线段Bez的参数a4阶次最低,鲁棒性最差。故δ4对分割误差敏感性远超过其对数据噪声的敏感性,会随着分割误差而快速增大甚至发散。曲线段Bez的其他参数相对圆弧参数的鲁棒性也比较差。若各元素分别用式(1)或(2)单独拟合,圆弧b与圆弧T元素只分别受σ1与σ2影响,而圆弧R与曲线段Bez则都受到二者影响;但若用式(3)整体拟合,任一分割误差均会影响所有参数。

1.3.2 用错元素类型造成的拟合不确定度 最可能出现的用错元素情况是曲线段Bez被当作圆弧处理。不考虑曲线段Bez与圆弧R的对称性,另设一圆弧B及其参数{xB,yB,rB},用圆弧B的取代Bezier的LSM目标函数以仿真用错元素。这样式(3)中共有12个参数。图2与表2给出部分MCM实验结果。

(a)用圆弧拟合曲线段Bez时参数xB拟合误差随分割误差的直方图

(b)用圆弧拟合曲线段Bez时参数yB拟合误差随分割误差的直方图图2 用圆弧拟合曲线段Bez的部分MCM实验结果Fig.2 The results of the MCM experiments on the Bez fitted by an arc

表2 用圆弧拟合曲线段Bez的参数期望与不确定度Table 2 The fitting results of the Bez by an arc

由实验结果看出其特点:没有分割误差时,圆弧B拟合误差直方图在一定程度上呈现期望不为0的正态性,但拟合准确度不如Bezier曲线,也不及与其相对称的圆弧R,且不确定度也比较大;有分割误差时,出现了栏栅效应,不确定度增大但不显著。由于样条曲线拟合算法鲁棒性较弱[14],圆弧B对分割误差的敏感性远小于Bezier,用式(1)拟合曲线段Bez的准确度可能比式(2)高。

2 降低曲线拟合不确定度

上述实验表明,对曲线拟合影响最大的是元素选择与数据分割,而非数据噪声。为了降低不确定度,以下提出一种从零件曲线中正确识别所含元素类型并准确分割其数据段的方法。

2.1 基本原理

待拟合数据序列包含多个元素时,每个数据段所属元素的参数特征恒定不变。由此,本文的元素识别与数据分割原理为:将待拟合数据序列按序划分为多个包含少量局部数据点的子集(称为数据子集),用于计算某元素的一组参数特征;顺序扫描各数据子集,如果参数特征恒定则相应子集属于该元素,不恒定则无该元素;逐个测试各个可能元素,通过观察相应参数特征恒定性,可以判断是否存在其数据段;以此类推,实现整个序列的元素识别与数据分割。

上述原理是在理想状况下的。实际上数据序列都包含噪声、错误等,各顺序数据子集即使属于同一元素,它们的特征一般也不相同。这就需要一种方法尽量抑制噪声、错误的影响,凸显元素几何特征,以利于其恒定性判断。高维叉积运算能满足该需求。

2.1.1 高维叉积的定义及其鲁棒性质[15]

定义对于n×(n+1)叉矩阵X={X1;…;Xn},对应的叉积运算定义为Xn+1=cross(X)={Xn,1,…,Xn,n+1},且

Xn+1,j=

式中:Xi∈Rn+1;i,j=1,…,n+1。

积向量的模为

‖Xn+1‖=

‖X1‖…‖Xn‖|B|0.5

式中:βij为向量Xi与Xj的夹角。

或

即‖Xn+1‖包含了Xi的模与方位信息,且与‖Xi‖之比为可变的比例系数cosαi/cosαr。由于αr越接近π,在Xr附近‖Xn+1‖对‖Xi‖的变化率越小,即Xn+1的鲁棒性越强;反之,若αr→π/2,则Xn+1的鲁棒性变差。故αr可用作积向量对Xi的鲁棒性指标。

2.1.2 基于高维叉积的元素识别与数据分割 对于待拟合序列Pi,i=1,…,m,顺序构造某测试元素的数据子集序列(Pi,Pi+s,Pi+2s,…),从中计算各相应的参数特征向量Xi∈Rn+1。为了观察该参数特征序列是否恒定,把它们映射到一个标量即积向量模

‖Xn+1(X)‖=‖cross[X;E]‖

(5)

式中:E={E1;…;En-1}∈R(n-1)×(n+1)为常数映射矩阵;设X的正交角为α。

当扫描某个数据段时,假设该数据段属于所测试的元素,其对应的特征向量为常向量Xr,对应的正交角为αr,有α,αr∈[π/2,π]。以Xr为参照,若噪声或数据误差使特征向量从Xr变为X,对应参考正交角变化为αr→α。选用不同的αr导致的比例系数cosα/cosαr随X的方位变化情况:若αr接近π/2,比例系数随X正交角而急剧变化,‖Xn+1(X)‖亦将随X而剧烈变化;若αr接近π,比例系数只随正交角在1附近稍有平缓变化,噪声或错误引起X的波动在‖Xn+1(X)‖中被抑制。如果数据属于被测元素,则‖Xn+1(X)‖保持恒定直至数据不再属于被测元素而出现转折或拐点,这样通过检测‖Xn+1(X)‖的恒定性及其范围而实现被测元素的识别与数据分割。

参考正交角αr取决于映射矩阵E的构造。

2.1.3 根据特征恒定性的要求构造映射矩阵 映射矩阵E决定于所选用的参考向量Xr。为了了解各数据子集的特征向量X是否与Xr相同,应该抑制X中的噪声干扰,以便于观察积向量模恒定性。所以,针对Xr构造映射矩阵E使αr→π。以被测试序列的第1个特征向量X1为参照,若其他特征向量X与X1相同,则式(5)中的积向量模能最大限度克服数据噪声与误差的影响而展现出相同性;反之,若X与X1不同,由于积向量模鲁棒性好,也能够最大限度克服数据噪声与错误的影响而展现出特征向量序列X与X1的差异,利于判断测试元素的不成立。由式(5)可得

(6)

式中:βr,i为Xr与Ei的夹角;βi,j为Ei与Ej的夹角。由于要求αr=π,可得cosβref,i=0,即Xr对E的约束为

2.1.4 根据准确定位分割点的要求构造映射矩阵 为了精确定位两个相邻元素的边界,应该抑制测试数据段在Xn+1的噪声,同时放大其相邻数据段在Xn+1的噪声,使它们之间形成强烈反差,以便于鉴别它们的分割点。假设相邻元素的参考特征向量为XN,则映射矩阵应该满足αN→π/2。XN对E的约束为

由于要求在XN附近放大噪声对‖Xn+1(X)‖的影响,所以cosβN,i→1。建立以E为变量的方程组

(7)

根据式(6),若要求βN,1=…=βN,n-1,则cosβN,i≈1/(n-1)0.5。

由此构造的映射矩阵满足了映射要求。

2.1.5 测试元素的特征向量 把数据子集序列转化为测试元素的特征序列是本方法的一个重要步骤。对于曲面,常用法向量用作为特征向量[16],而对于简单一维元素,则可用其主要参数做其特征向量。

(1)圆弧特征:用圆弧所在的圆半径r与圆心(xc,yc)作为圆弧元素的特征向量,即X={xc,yc,r},其数据子集为(xi,yi),i=1,…,mc,用式(1)将该数据子集转化为特征向量。

(2)Bezier曲线特征:数据子集{Pi,Pj,Pk}可用参数化方法转化为Bezier曲线元素的特征向量,但这样的特征向量没有独立性。用抛物线参数Xi={b1,…,b4}作为Bezier曲线特征,其对应的数据子集为{Pi,Pi+s,Pi+2s,Pi+3s},用式(2)进行数据子集到特征向量的转化,则特征向量有独立性或局部性。

2.2 凸轮轮廓元素识别与数据分割实验

考虑到每个元素数据段的起始点实际上是未知的,任取一点如P1作为识别、分割扫描的起始点。通过观察特征向量的恒定性确定,一个数据段是否属于测试元素,如果属于则找出其终点。如此反复直至完成整个数据序列的扫描。

2.2.1 霍夫变换的圆弧数据识别分割 传统方法在数据点Pi处构造圆弧法线做为霍夫变换,容易受噪声影响[5]。因此,改用弦PiPi+s的平分线作为霍夫变换,即把凸轮轮廓数据点Pi(xi,yi)映射为

其中可取数据步长s≡常数,i=1,2,…;也可取s=1,2,…,i≡1。如果以各连续点的霍夫变换交点O作为其所在圆弧的圆心,即若Pi辨识为一段圆弧的点,则Pi+1,…,Pi+s也都位于该圆弧上,并可求出圆弧半径,找到圆弧的最后一个点后,即完成了该圆弧的数据分割。

这种圆弧数据分割精度受霍夫矩阵维数及所设的圆心与弦平分线距离容差影响比较大。经过反复尝试,取霍夫矩阵维数5 400×5 520与s≡3,得到弧b的圆心为(3.52,-4.17) μm,相应数据段下标为

i∈[1,116]∪[179]∪[245,360]

多出了下标为i=116,245,179的3个点。维数低容易误分,降低容差则容易丢失一些点。

2.2.2 高维叉积映射的数据分割 先以圆锥曲线为测试元素,其数据子集为{Pi,Pi+s,Pi+2s,Pi+3s,Pi+4s,Pi+5s},数据步长s=6。对应的特征向量序列Xi={b1,…,b5},以Xr=X1设计映射矩阵,用式(5)进行映射,结果见图3a。图中显示在下标i∈[1,85]范围内的圆锥曲线特征向量是恒定的,从i=86开始为其他特征向量。由于X86对应的数据子集为{P86,P92,P98,P104,P110,P116},可知P116为下一个数据段的第1个点。所以下标i∈[1,115]的数据段为一段圆锥曲线。

以圆弧为测试元素,其数据子集为{Pi,Pi+s,Pi+2s,Pi+3s},s=1,对应的特征向量为Xi={xc,yc,r}。以Xr=X1,XN=X116设计映射矩阵,再进行特征向量序列的映射,见图3b,下标为i∈[1,112]范围的特征向量是恒定的;i=113处出现拐点,而对应的数据子集为{P113,P114,P115,P116},可知P116属于下一元素的点。同样可得知,下标i∈[1,115]范围内的数据段为一段圆弧。另外,图3a中i=86的拐点由特征向量的变化产生的,其后的特征向量相对于参考向量变化较大,缓慢过度则可能无法形成清晰的拐点;图3b中虽然水平噪声相同,但是两个参考向量的比例系数不同,使i≤113与i≥115分别形成了两种噪声形态;两个邻域之间过度区出现的阶跃是由两边特征向量的差异造成。因此,噪声形态差异与特征向量差异均显示出数据分割点。

(a)测试元素为圆锥曲线

(b)测试元素为圆弧图3 两种测试元素的映射局部放大图Fig.3 Detailed views of testing conic and arc

用相似方法实施整个数据序列的分割,得到4个元素对应的分割点与它们实际值完全一致。

数据分割受很多因素影响,其中包括数据中的噪声与误差。虽然增大数据子集或数据步长可在一定程度上减小噪声影响,但会加大过渡区而不易精确定位分割点。此外,数据子集转换为特征向量时,不同元素的转换精度也不同。Bezier曲线比圆弧的转换更易受噪声影响。

3 实物对象的元素识别与分割拟合

图4 凸轮原型Fig.4 A physical cam for the verification

图5为两种传统数据分割的拟合效果,图5a为近似按平均点数分割点集为11段,而图5b采用主观识别的点云分割。图5a中,由于采用弧长均匀采样,曲率大的曲线段偏差较大。图5b中,部分节点与分割点不重合,说明分割不合理。这两种拟合效果均不理想。

(a)均匀分割

(b)主观分割图5 两种传统分割的Bezier曲线拟合效果Fig.5 Two Bezier fittings with traditional segmentations

在图6a～6c中,取数据子集{Pi,Pj,Pk}并用CCL参数化法[19]获得Bezier曲线控制点特征

式中:lk=Σ‖Pj-Pj-1‖,j=1,…k,即Pk对应参数tk≈lk/L,L为所在曲线的累计弦长。利用式(5)映射为积向量模,辨识得到测量序列的前3个分割点。如此逐段辨识其他各分割点,得其下标构成的向量

i={13,24,50,83,98,109,129,146,173}

(a)从i=1开始CCL识别

(b)从i=13开始CCL识别

(c)从i=24开始CCL识别

(d)从i=1开始抛物线识别(局部)图6 用二阶Bezier曲线特征识别测量序列Fig.6 Testing the cam contour with second order Bezier curve

总共9个分割点把数据序列分成10段Bezier曲线。图6d取数据子集{Pi,Pi+1,Pi+2,Pi+3},用式(2)计算抛物线特征向量Xi={b1,b2,b3,b4},图中分割点基本与CCL法一致,但只需一次映射即可显示所有拐点和分割点。若只用一个数据段(例如第1段或第2段)的参照特征设计映射矩阵E,可以获得较清晰的积向量模拐点,而其他段的拐点则可能不够清晰,而且抛物线的特征鲁棒性比CCL及圆弧等的差。图7显示了识别分割的拟合效果,其中两个坐标的误差平方和均小于前2种分割拟合。

图7 辨识分割的Bezier拟合Fig.7 The Bezier fittings after recognizing segmentation

4 结 论

(1)对拟合不确定度影响最大的因素是拟合元素选择的正确与否以及相应数据分割准确程度。拟合算法与数据噪声等因素也影响拟合精度,但远不及前二者。

(2)由数据噪声等因素产生的拟合误差分布为正态形状;错用元素的拟合准确性降低而分割误差则会使拟合误差分布呈现多峰、栏栅效应等复杂形态,不确定度大大增加。

(3)基于GPC鲁棒性,将任意维特征向量映射到标量并观察其恒定性及范围,可有效识别元素类型与其边界,从而降低了拟合不确定度。

(4)与仿真不同,实际零件元素及其数量未知,其GPC映射除了会受测量精度、噪声,还受拟合算法的鲁棒性、数据子集步长等因素影响。多种测试元素的映射均有恒定性时,可取恒定性范围最大的作为识别结果。除了圆弧等简单元素,二阶Bezier曲线拟合具有通用性,可获得较好效果。

(5)Bezier的CCL法特征无局部性,其积向量模拐点也可能比较模糊,需逐段辨识各分割点;抛物线特征具有局部性,可以一次识别多个分割点,利于降低拟合不确定度。