何帅,陈富民,李建华,张厅方

(1.西安交通大学新能源与质量工程研究所,710049,西安; 2.东方汽轮机有限公司,618000,四川德阳)

目前对于精密零件的测量,一般采用在三坐标机上测得坐标,再通过相关算法计算所需要的零件质量特性,这种方法属于坐标法。采用坐标法测量零件坐标时,由于无法直接利用被测零件本身基准为测量基准,测量数据里包含测量坐标系与设计坐标系不重合造成的位置误差。叶片属于精密零件,采用坐标法对叶片型线进行质量评定,精度受基准不重合影响比较大,因此需消除测量基准与设计基准不重合对偏差评定的影响[1-2]。轮廓度评定时需要测量坐标与理论坐标的位置一致,按照最小区域原则[3]进行叶身型线的轮廓度误差评定,但实现叶身型线最小条件比较困难[4]。

在评定过程中有两个核心问题需要解决:实测点所对应理论轮廓的最近点计算;如何统一实测基准与设计基准。对于第一个问题,一般对理论型线采用样条函数进行拟合,对拟合的理论型线采用分割逼近[5]的方法转换为点与点的搜索。对于第二个即基准统一的问题,本质是刚性变换求解的问题,目前求解此问题主要有两种方法:第一种方法是采用进化算法求解变换参数,第二种方法是采用二维最近点迭代算法(ICP)进行求解。采用进化算法求解变换参数需要构建不同的适应度函数,如遗传算法[5-8]或者粒子群算法[9]求解基准变换,理论上可以求到变换参数的准确解。进化算法求解齐次变换,需要生成大量初始解,计算量大导致求解配准参数耗费时间较多。由于进化算法的随机性,导致配准计算的稳定性很差,可能陷入局部最优。依据形位误差评定原则最小区域法构建适应度函数,计算出的旋转量与理论值差异较大,而平移量较为准确。

许多学者依据文献[10]提出的最近点迭代算法,提出用二维形式的最近点迭代算法计算实测轮廓与理论轮廓的配准参数。文献[11]根据不同区域的公差不同,引入点对权值约束,实现轮廓的高精度配准。文献[12]引入熵的概念,通过最大化熵准则的最近点迭代算法,实现有噪声和离群点的二维点集配准,但其只适用粗配准精度。文献[13]针对平面零件二维轮廓与矢量图形的配准问题,使用二维最近点迭代算法的解析形式进行轮廓间的精配准,可以较好地实现轮廓的位置调整,但是平移量计算误差受实测轮廓影响较大。二维最近点迭代算法求解变换参数速度快,依据原则为最小二乘原则,旋转量求解非常精确,但由于实测点与理论型线总是存在偏差,采用质心平移的方法获得的平移量误差较大。

基于以上原因,本文提出一种最近点迭代结合粒子群算法的叶身型线基准误差求解方法。首先采用非均匀有理B样条(NURBS)曲线的矩阵形式对叶身型线进行拟合,利用区间搜索的方法寻找实测点所对应理论轮廓上的最近点及偏差;其次运用二维最近点迭代的解析形式计算基准变换的旋转量,然后对最近点迭代计算得到的平移量加微小扰动,以最短距离均值为优化目标,使用自适应权重粒子群算法重新计算平移量。本文对实测点进行基准变换后计算偏差,得到实测轮廓的轮廓度,设定理论变换向量构造仿真数据,以计算结果与理论值的差异验证方法的有效性,并与目前自由曲线配准的方法进行对比。

1 叶身型线数学模型

叶身型线属于自由曲线,理论形式以离散点坐标给出。非均匀有理B样条的出现为标准解析形状和自由曲线提供了统一的数学表示[14]。国际标准化组织(ISO)在1991年颁布的工业产品几何定义的STEP标准中,NURBS被定义为唯一的自由型曲线曲面表示方法[15],因此本文采用NURBS对叶身型线进行描述。

常用的NURBS曲线的表达形式有两种,一种是德布尔算法的递推形式,另一种是矩阵形式。虽然德布尔算法概念清晰,但算法复杂,计算量大。曲线轨迹上任意点的求取速度受基函数计算影响较大,为提高曲线的计算速度,采用NURBS曲线矩阵形式计算曲线的系数[16]。NURBS曲线的非均匀性是指节点矢量与间距为任意值,可在不同区间上得到不同的混合函数形状[17],曲线随参数变化而变化,可以很好地对自由曲线进行描述[18]。

由n+1个控制顶点定义的k次NURBS曲线的数学定义[14]如下

(1)

式中:C(u)为NURBS曲线的参数形式;u为参数;wi为权重因子;pi为控制点坐标向量;Ni,k(u)为k次B样条基函数,表达式如下

(2)

Ni,k(u)=

(3)

(4)

式中:Mi(t)为系数矩阵,表达式为

(5)

式中:m11、m13、m23、m33、m43、m44分别代替矩阵Mi中的元素。

NURBS曲线最终可以化为以下形式,根据理论坐标数据计算得到n-2组系数

(6)

2 叶身型线配准

在叶身型线测量过程中,由于测量基准与理论设计基准存在不重合误差,对实测点偏差的计算精度影响较大,所以需要减小基准不重合对质量评定的影响。对于基准间的刚性变换参数的求解问题,可以视为实测点集与理论轮廓对应的最近点集间的配准问题。最近点迭代(ICP)算法的优化目标符合形位误差评定的最小二乘原则,算法求取旋转量非常精确,由于存在加工误差,求取的平移量与真值存在较大差异。因此,以ICP算法计算旋转量,以所有实测点偏差的绝对值均值作为适应度函数,对ICP算法获得的初始平移量加微小扰动,采用自适应权重粒子群算法重新求取平移量。

2.1 实测点与理论轮廓的轮廓度计算

理论点集为有序点集,对理论点集进行NURBS矩阵形式拟合,计算实测点的偏差,公式如下

pd=sdmin

(7)

式中:pd为实测点的偏差;s为位置符号;dmin为实测点与理论轮廓对应的最近点间的距离。

实测点的位置符号s可以用实测点与理论轮廓上对应的两个最近点组成的叉积进行判断。实测点可能位于理论轮廓内,也可能位于理论轮廓外,实测点在理论轮廓外侧时,s为1;在理论轮廓内侧时,s为-1。

实测点对应的理论轮廓最近点计算步骤如下。

(1)求出实测点到理论轮廓的最近距离点所在段的索引。

(2)传入最近点所在理论轮廓段的参数,对参数进行区间划分搜索最近点,如黄金分割或二分法。

当测量基准与理论设计基准统一时,轮廓度E定义为

E=max(pd)-min(pd)

(8)

式中:max(pd)为所有实测点偏差的最大值;min(pd)为所有实测点偏差的最小值。

2.2 改进的叶身型线配准算法

实测点集P,数目为m,理论点集Q0,数目为mQ,对理论点集Q0进行NURBS矩阵形式拟合。改进的叶身型线配准算法原理如下。

2.2.1 以二维最近点迭代的解析形式求解旋转矩阵R和初始平移向量T0ICP算法依据的原则是最小二乘,优化的目标为最短距离平方和最小,数学形式如下

(9)

式中:Qi为实测点集P对应理论型线的最近点集Q中的第i个点;Pi为实测点集P的第i个点。

二维ICP算法的解析算法步骤如下。

Step1:对实测点集P搜索对应理论型线的最近点集Q。

Step2:对点集P和点集Q去中心化,分别得到点集P1和点集Q1,去中心化公式如下

(10)

(11)

式中:Q1i为最近点集Q1中的第i个点;P1i为实测点集P1的第i个点。

Step3:计算旋转矩阵R。二维点集配准旋转矩阵R由一个旋转角度θ构成,公式如下

(12)

将式(12)代入式(11),得到

(13)

对式(13)所示的函数进行求导,得到

(14)

求得f(θ)的极值点

(15)

Step4:求取平移向量

(16)

Step5:更新实测点集P,更新公式如下

P=PR+T0

(17)

Step6:重复Step1～Step5,直到算法收敛,得到基准变换的旋转矩阵R和初始平移向量T0。

2.2.2 采用自适应权重粒子群算法重新计算基准变换的平移量T由于存在偏差,ICP算法难以精确获得基准变换的平移量,建立如下适应度函数

(18)

式中:pdi为第i个实测点的偏差。

运用自适应权重粒子群算法[20],以T0加微小扰动生成初始解空间,重新计算平移向量,计算步骤如下。

Step1:设置算法参数,种群数量N为30,自我学习因子c1为0.9,群体学习因子c2为0.9,最大惯性权重wmax为0.8,最小惯性权重wmin为0.2。

Step2:初始化解空间和速度,计算所有粒子的适应度值,并且找出解空间的最佳个体zg和最佳适应度值fz

(19)

式中:xi为解空间的第i个解向量;vi为解空间的第i个解的速度向量;N1为均值为0、方差为0.04的正态分布生成的2维度的随机向量;N2为均值为0、方差为1的正态分布生成的2维度的随机向量。

Step3:进入主循环,计算种群的适应度,并且记录种群最小适应度fmin及种群适应度的均值favg。计算个体的惯性权重w,更新速度及个体,更新规则如下

(20)

式中:f为第i个粒子的适应度值;zi为第i个粒子在迭代过程中的历史最佳个体。

Step4:最后比较个体更新前后的适应度值,如果更新后的个体适应度更小,则保留更新后的个体,否则保留更新前的个体。计算更新后的种群适应度值,寻找出最佳个体及其适应度值代替zg和fz。

Step5:重复Step3～Step4,直到收敛或者固定迭代次数,输出最佳个体zg,即重新计算的平移量T。

计算出旋转矩阵R和平移量T后,对实测点集进行基准变换,此时实测型线与理论型线达到最佳配准。

3 算法验证及实例分析

3.1 样本构造

本文样本的构造方法如下。以理论数据Q0(有序点集,数目为200)进行NURBS拟合,得到199段的理论轮廓。以正态分布N(u,σ2)生成两组199个偏差pd,其中u为0,σ为0.03。根据参数曲线的等距公式生成偏距点集P0,每个偏距点的位置取理论轮廓曲线段的中间位置。构建基准变换向量T1(0.02,0.1,-0.1),其中0.02为旋转角度,0.1为x方向平移量,-0.1为y方向平移量。将点集P0进行T1的基准变换,变换后的点集P2作为实测数据。基准变换T1为理论基准到实测基准,而配准是求实测基准到理论基准的变换向量,所以配准参数的理论值应为基准变换向量T1构成的齐次刚性变换矩阵的逆矩阵对应的基准变换向量。配准的理论基准变换向量为(-0.02,-0.098,0.102),样本1偏差极差为0.215 4,样本2偏差极差为0.163 6。本文算法实验在Matlab2017b中进行。

3.2 不同适应度函数求解结果的对比

遗传算法和粒子群算法同属于进化算法,自适应权重粒子群算法简单,收敛较快,是一种全局寻优的进化算法。本文分别以式(8)(18)作为适应度函数,以自适应权重粒子群算法优化上述适应度函数,计算实测基准到理论基准的变换向量。以求出的变换向量与理论变换向量的欧式距离度量基准转换的误差,样本数据及配准结果如图1所示,计算结果如表1所示。

(a)样本1数据(b)样本2数据

(c)轮廓度配准结果(样本1)(d)轮廓度配准结果(样本2)

(e)最短距离均值配准结果(样本1)(f)最短距离均值配准结果(样本2)

图1 不同适应度函数样本数据及配准结果

综上所求解的基准变换向量理论值为(-0.020,-0.098,0.102),样本1的偏差极差为0.215 4,即样本1的轮廓度理论值为0.215 4;样本2的偏差极差为0.163 6,即样本2的轮廓度理论值为0.163 6。样本1(轮廓度)的轮廓度计算值与理论值差异为-0.032 5,样本1(最短距离均值)的轮廓度计算值与理论值差异为-0.015 8。样本2(轮廓度)的轮廓度计算值与理论值差异为0.011 4,样本2(最短距离均值)的轮廓度计算值与理论值差异为0.000 6。

综合样本1和样本2的计算结果,经以最短距离均值最小为优化目标,求得的基准变换向量更接近理论值,所求的轮廓度也更接近理论值。因此,本文方法选择最短距离均值为适应度函数,重新求取平移量。

3.3 本文算法与ICP算法的对比

标准ICP算法采用SVD分解去中心化后的点集的协方差矩阵。本文推导了二维ICP算法的解析形式,避免出现协方差矩阵分解失败的情况。ICP算法依据最小二乘原则,求解刚性变换参数非常快,但由于存在加工误差,平移量计算误差较大。根据3.2小节,以最短距离均值作为优化目标,计算结果更加贴近真值,因此采用自适应权重粒子群算法重新求取平移量。由于减小了优化的维度并且得到了初始的平移量,算法计算速度更快,同时保证了精度。分别以ICP算法和本文算法对样本1和样本2进行配准,以计算所得变换向量与理论变换向量的欧式距离度量基准转换误差,样本数据集配准结果如图2所示,计算结果如表2所示。

由表2可以看出,采用ICP算法虽然计算速度快,但计算出的基准变换向量误差较大,计算所得的轮廓度明显偏离理论值,主要影响因素为平移量计算不准确。本文算法基准变换向量计算值与理论值的误差在0.006左右,样本1的轮廓度计算相对误差为2.8%,样本2的轮廓度计算相对误差为1.2%,计算结果已经非常接近理论值。与采用进化算法相比,本文所提算法直接计算所有变换参数的速度要快40%左右。

(a)ICP算法配准结果(样本1)(b)ICP算法配准结果(样本2)

(c)本文算法配准结果(样本1)(d)本文算法配准结果(样本2)

图2 ICP算法与本文算法配准结果的对比

表2 ICP算法与本文算法配准计算结果的比较

4 结 论

本文采用NURBS的矩阵形式对叶身型线进行拟合,一次性计算出拟合曲线的所有参数,避免了德布尔算法递归形式的大量计算,采用区间搜索求解点到曲线的最近点,计算速度快,精度高。此外,提出了一种叶身型线配准算法,首先通过二维ICP算法的解析形式精确获得配准的旋转量,其次以实测点集最短距离均值作为适应度函数,采用自适应权重粒子群算法重新优化平移量,基准变换向量误差在0.006左右,相对于采用进化算法直接计算配准参数的速度提高了40%左右。所提算法适用于三坐标测量叶身型线数据后对叶身型线进行评价,能够有效减小由于基准不重合对偏差计算的影响,避免叶身型线加工质量误判。