基于动态精度评估机制的自适应定位方法

谭高山，张 涛，张丽艳，李 博

(1.安徽工业大学 数理科学与工程学院，安徽 马鞍山 243032; 2.南京航空航天大学 机电学院，江苏 南京 210016；3.河南工程学院 计算机学院, 河南 郑州 451191)

0 引言

为了提高产品性能，复杂曲面的应用日益广泛。复杂曲面是一种极难定义和加工的几何元素，具有复杂曲面的零部件的加工定位、几何精度评估等问题也对复杂曲面的几何定位提出了更高的要求。复杂曲面定位是保证复杂零部件几何精度的关键技术之一。

近年来，三维数据定位研究方面的成果较多[1-2]。因为区域设计精度和制造误差存在差异，基准定位不再适合复杂曲面精度评估问题，所以满足不同拟合准则的定位数学模型及其求解是复杂曲面评估定位的核心问题。基于最小二乘准则的定位方法，因其计算效率高且鲁棒性好，应用最为广泛[2-5],求解算法也较多，如四元数法求解的迭代最近点(Iterative Closest Point，ICP)算法[3]、最小二乘目标函数的伪逆法[4]、拟粒子群算法(Quantum-behaved Particle Swarm Optimization，QPSO)[5]等。由于受测点数量和分布等测量问题的影响，最小二乘定位的结果倾向于大偏差数据点和稠密数据点，从而过高估计了高精度区域的误差[2,6]。事实上，有些复杂曲面不同区域的数据本身具有不同的精度，考虑整体数据的最小二乘拟合并不适合复杂曲面误差评估问题。学者们也研究了基于其他数据拟合准则的定位方法，如：利用向量法和指数罚函数法求解的极差极小化原则的定位方法[7-8]、转化为带约束非线性优化问题并用内点法求解的最小区域定位方法[9]、用线性Tayler展式逼近的非线性定位方法[10]等。基于这类准则的定位目标函数是非光滑的，一般通过转化为近似光滑模型进行逼近求解，求解精度不高，实验部分也只给出了直线、平面、圆、圆柱等简单几何特征的定位结果[9]。显然，非光滑定位模型并不适合复杂曲面评估问题，特别是测量数据量较大时，计算问题更加突出。

基准定位中被选做基准的高精度曲面具有绝对的定位优势，非基准区域在定位中不被考虑，因此不适合复杂曲面定位。除基准定位,目前大部分定位算法和商业软件对测量数据都是整体进行处理，即不同区域的数据点对定位结果的影响是相同的。由于在工艺、效率、成本和产品性能等综合因素下，复杂曲面各区域的精度并不相同，无区别地对待所有测量点的定位方法与实际情况不符。

为满足特殊工程要求，学者们建立了约束最小二乘定位模型[11-13]，并分别用序列线性规划法、有限存储的BFGS算法(Limited memory Broyden-Fletcher-Goldfarb-Shanno，LBFGS)广义乘子法和遗传法-单纯形法求解。理论上，可以通过约束不同区域的精度进行误差评估问题的定位，但约束优化问题的迭代求解对约束边界比较敏感，受设计精度和制造精度的影响，工程上很难找到复杂曲面的合适精度约束[14]。另外，约束最小二乘方法需要求解带约束的非线性最优化问题，这类优化问题求解效率并不理想。基于几何特征的定位[15-16]与基准定位类似，这类方法考虑精度差异对定位结果的影响，按照特征面精度高低依次进行定位，先选取的特征构成了欧氏群的一个齐次空间，并对模型进行定位，后增加的特征只能在补空间中进行定位迭代，以不影响已有空间中的定位结果。若先选取的特征完全决定了6个自由度，后面的特征将不再参与定位。虽然这类方法区别对待不同特征，但过分强调了精度差异性，没有考虑量的影响，显然也不适用于区域精度不同的复杂曲面定位问题。

基于简单几何特征的计算几何理论，利用高效且鲁棒的加权最小二乘方法[17]可得到平面度和圆柱度的理想定位结果。较早的加权最小二乘思想[18]利用离散点代表的面积表示定位权因子。不同加权方式[19-21]的主要目的在于处理粗差点对计算速度和定位精度的影响。加权最小二乘模型可通过推广最小二乘定位求解方法来实现，如求解加权最小二乘问题的四元数法[22]。如何根据工程需要给出合适的加权策略一直是加权最小二乘法的核心问题。

定位约束分析是讨论定位信息对定位结果的影响。用灵敏度矩阵[4]来度量测量点集在扰动下的定位稳定性，并用以估算定位误差；用信息矩阵[23]的行列式来定义定位精度；从测量点到理论曲面的梯度变化的角度提出定位约束矩阵[24]也可进行类似分析。

针对区域精度不同和测量点自身对定位的影响问题，本文提出一种基于动态精度评估机制的自适应定位算法，迭代过程中根据区域精度预估实时更新不同区域的定位权因子，因为定位迭代是基于曲面精度预估的，所以本算法选取相对精确的最小二乘定位作为初始位姿。通过不同权因子控制定位信息阵来约束定位方向，从而控制不同曲面定位精度。本算法中，每步迭代优化问题通过求解四阶线性方程组实现，求解速度快且具有较好的鲁棒性。

1 问题分析与模型思想

假设复杂曲面测量点集为P={Pi|i=1,2,…,n}，其对应理论曲面为S，集合P中的点是在测量坐标系下获得的，而理论曲面是相对于设计坐标系定义的。曲面定位就是要找到一个刚体变换，对齐设计坐标系和测量坐标系，使得测量点Pi(i=1,2,…,n)与其在理论曲面上对应点Qi之间的距离最小，这一距离度量可以定义为如下范数：

di(R,T)=‖RPi+T-Qi‖p。

(1)

其中,R∈3×3为空间三维旋转矩阵，T∈3为平移向量；‖‖p表示两数据模型之间的距离，当p=1时，为Chebyshev距离，当p=2时，即为欧氏距离，这里p可根据工程需要定义不同的度量。

(2)

T=(tx,ty,tz)T。

(3)

式中：s和c分别表示正弦函数和余弦函数,θx,θy,θz和tx,ty,tz分别表示沿3个坐标轴的旋转角度和平移大小。

为了使测量点定位到理论曲面，最小二乘定位要求所有测量点到理论曲面距离的平方和最小，这会过高估计高精度区域的误差，不能满足区域精度存在差异的复杂曲面定位要求。本文算法的基本思想是：首先按照精度要求将曲面进行分割，并利用最小二乘定位作为初始位置；然后利用区域距离的统计特征预估区域精度，构造相应的区域定位影响因子；最后，在每一步迭代中优化调节不同区域的定位姿态，从而自适应地获得测量数据与理论曲面的定位结果。

基于动态精度预估机制，建立定位问题的自适应加权数学模型如下：

(4)

式中：dj(R,T)表示数据模型对应点间的距离;m为不同精度区域的分类数;li为第i类曲面对应的测量点个数;ωi为第i类曲面对应定位影响因子，由预估精度来确定，它表明该类曲面在定位过程中的作用。

2 基于动态精度评估机制的定位方法

2.1 不同区域对定位结果的影响

定位的目的是获得测量数据相对理论曲面满足工程要求的空间位姿，具体来说就是获得测量点与其在理论模型上对应点之间的最佳距离分布。本文通过给不同测量点对之间的距离赋予相应权因子来控制不同精度数据点在定位中的影响。由自适应加权优化模型(式(4))可获知权因子与点对间距离变换的关系：权因子与距离变化方向相反，即权因子变大，绝对距离减小；权因子减小，则距离变大。可见，调节权因子可控制不同区域曲面与目标曲面之间的距离分布。本文通过动态地统计分析两数据模型之间的精确位姿，预估区域精度，获得下一步定位的权因子。该方法适用于具有不同区域精度的曲面的定位。为了清楚地说明不同区域曲面对整个定位结果的影响，用如图1所示仿真工件进行自适应加权定位实验，数据模型均采用均匀的三角网格描述。该仿真工件为立方体的3个相邻表面构成的组合曲面，3个不同精度的区域平面更能直观显示定位结果。图1a为理论模型，3个平面分别为Q1,Q2,Q3。图1b是3个面加上从小到大不同程度噪声后的模型，模拟实际测量数据中各区域从高到低的精度。图1c显示了模拟数据的理论精度情况。图1d和图1e中蓝色连续曲面表示理论模型，法矢向外；黄色网格表示定位变换后的测量模型。基于同一初始状态，图1d和图1e分别是经典的ICP算法和本文算法的定位结果。图1d中，精度最高的面Q1的误差估计过高，整体发生了偏移，误差分布与实际不符；而精度最低的面Q3的误差评估结果却得到了改善。图1e中定位结果更接近理论定位，特别是精度要求较高的面Q1获得了较好的定位结果(第3章将分析算法的工作原理，揭示两种定位的作用机理，并给出具体数据分析)。可见，区域精度对定位结果的影响不容忽视，是提高定位精度和适应性的关键。

2.2 基于精度预估机制的定量影响权因子

为给出不同精度区域对整体定位结果的影响，首先要对理论曲面进行分类，依精度不同将理论曲面分为m类：S={Q(i)|i=1,2,…,m}; 然后根据区域精度预估，分别赋予不同曲面类相应的权因子ω1,ω2,…,ωm。设与曲面类Q(i)对应的测量点云为S(i)，定位目的是找到两对应数据之间的理想相对位姿，以满足误差或精度需要。均方根误差和平均误差是常用的两个衡量曲面精度的指标，定义

(5)

rmsi和avdi分别为第i类曲面的均方根误差和平均误差，其中均方根误差最小与该类曲面的最小二乘定位目标表示一致。不考虑区域精度差异的最小二乘定位获得了整体的最小均方根误差，即各区域均方根误差均匀分布，从而导致高精度区域误差估计过大，因此，给出区域定量影响权因子，是区域精度不同的复杂曲面进行适应性定位的必要前提。

分析第i类测量点和对应理论曲面的当前位姿，根据2.1节中的分析，利用式(5)两个定位指标构造第i类曲面下一步定位变换的影响权因子如下：

(6)

式中dj(R(k),T(k))是该i类曲面当前第j个对应距离。统计当前位姿的均方根误差和平均误差，给出区域精度预估，并度量下一步迭代定位中区域精度的定量影响程度，用定位指标rmsi和avdi优化控制定位位姿。每一步迭代前的区域精度预估，给出了式(4)定位方法中不同精度区域在定位中的影响程度，不但有效防止了低精度曲面导致的整体位姿发生偏离，而且对不同区域的定位影响进行了相应度量，使得定位结果具有更高的适应性。

2.3 自适应加权最小二乘问题的线性化求解

定位问题中的距离函数di(R,T)=‖RPi+T-Qi‖是旋转变量θx,θy,θz和平移变量tx,ty,tz的非线性函数，式(4)所示加权定位模型为高度非线性优化问题。考虑空间三维旋转变换和平移变换的齐次表示

(7)

式中：g为变换的线性表达矩阵，因此距离函数转化为向量ξ的线性表示，令

ξ=vec(g)。

f(ξ)=ξTAξ-2bTξ+c。

(8)

其中:

(9)

(10)

(11)

Aξ=b。

(12)

式中A，b如式 (9)和式(10)所示。

利用矩阵分块思想，可通过求解3个系数矩阵相同的四阶线性方程组得到方程组式(12)的解，并由此得变换矩阵g，进而获得旋转矩阵和平移变量。由于每次迭代中变换矩阵g(k)可以通过线性求解ξ(k)得到，不需要迭代求解非线性问题，因此求解速度快，鲁棒性高。用经典的奇异值分解法或四元素法求解加权最小二乘问题时，需要分别利用奇异矩阵和特征向量构造变换矩阵，计算量比本算法大，实验部分将给出具体效率比较。

2.4 算法步骤

本文所提基于精度预估机制的自适应加权最小二乘定位方法通过动态更新权因子，预估不同曲面区域在定位中的定量影响。由于算法对初始位置要求较高，本文选取经典的ICP算法作为初始定位。在每次定位迭代过程中，先预估不同区域的制造精度，并由式 (6) 构造区域定位影响权因子，因为高精度区域在定位优化求解中得到优先考虑，所以得到更好的定位。算法流程如图2所示，具体算法步骤如下：

步骤1初始化旋转矩阵R为单位矩阵，平移向量T为零向量，令初始迭代次数k=0，给定各曲面类初始权因子为1/li(i=1,…,m)，旋转变量和平移变量的停止条件为ε1,ε2，定位指标rmsi和avdi的停止条件为ε3,ε4，并设定最大迭代次数kmax。

步骤2对于测量数据点Pi，搜索理论曲面模型上的最近点Qi，并获得Qi所在曲面类的权因子ωi。

步骤3利用四阶线性方程组求解加权最小二乘问题式 (4)，得到齐次变换阵，利用式(13)计算旋转角θx,θy,θz和平移量tx,ty,tz，再利用式 (2) 计算旋转矩阵和平移向量Rk和Tk。

θy=arcsin(r13),θx=arctan(-r23/r33),θz=arctan(-r12/r11),tx=r14,ty=r24,tz=r14。

(13)

其中rij(i,j=1，…，4)为式(7)中变换矩阵元素。

R=RkR,T=RkT+Tk。

(14)

步骤5若式(15)中变换停止条件ΔR≤ε1,ΔT≤ε2成立，转步骤6，否则转步骤2。

(15)

步骤6若式(16)中位姿停止条件Δrmsi≤ε3,Δavdi≤ε4成立或者k=kmax，停止迭代，退出；否则转步骤2。

(16)

3 加权定位约束

所提加权最小二乘定位模型通过权因子控制定位结果，从而获得不同区域的误差分布。本章研究权因子在定位中的作用机理。考虑算法的优化目标函数f(X)随空间变量X=(θx,θy,θz,tx,ty,tz)T的变化，由全微分公式得到

(17)

E(f)=dXTAωdX。

(18)

其中二次型E(f)表示目标函数值随自变量的微小改变而发生的改变，其中二次型的系数矩阵Aω如下：

(19)

这是一个6×6的对称半正定矩阵，由于Aω中包含了加权定位中的权因子、测量数据点及其对应点，称Aω为加权定位信息矩阵。为研究空间变量的扰动对目标函数造成的影响，对加权信息矩阵进行奇异值分解，得

Aω=QΛQT=(q1,q2,q3,q4,q5,q6)diag(λ1,λ2,λ3,λ4,λ5,λ6)(q1,q2,q3,q4,q5,q6)T。

(20)

其中diag(·)表示对角型矩阵；λ1≥λ2≥λ3≥λ4≥λ5≥λ6为Aω的特征值；qi(i=1,2,…,6)表示对应于特征值λi(i=1,…,6)的单位特征向量，特征向量表示变换变量的微分变换。现将式(20)代回式(18)得到

(21)

式中q1是最大特征值对应的最大约束变换，微分变换q1引起的能量变化最大，说明该方向上定位约束最大；反之，最小特征值对应的特征方向上微分变换产生的能量最小，若最小特征值为零，则其方向上的定位处于自由状态，即该方向上的微分变换对定位目标函数无影响。由式(21)可知，λi与qi变换引起的函数变化成比例，因此大的特征值是主约束力，对应的特征向量是主约束方向。迭代定位中，加权定位信息矩阵Aω通过权因子调节主方向和主约束力，从而调节不同权对应区域的定位结果。

图1利用仿真数据模型的定位结果说明了动态定位权因子的重要作用。这一数据模型是立方体3个互相垂直的相邻面，分别平行于3个坐标面，这样不仅容易观察定位结果，还便于直观理解定位约束分析中的空间变换。本章利用这一数据模型进行定位约束分析，说明本文算法的作用机理。图3a～图3c是初始状态下3个面的正视图。对ICP算法和本文算法进行定位约束分析，表1给出了不同权因子的最大约束变换方向。ICP算法没有考虑曲面的不同区域在定位中的影响程度，因此可视为区域权因子均为1的加权定位，最大约束变换方向为绕x轴顺时针旋转0.45弧度，由图3a可见曲面Q2,Q3精度得到改善；绕y轴顺时针旋转0.25弧度，由图3b可见曲面Q1,Q3定位位姿得到改善；再沿x轴负方向平移0.85 mm，曲面Q1定位精度下降。本文算法中，曲面Q1的权因子相对变大，最大约束变换方向为绕z轴顺时针旋转0.095弧度，由图3c可见曲面Q1,Q2位姿得到改善；绕y轴顺时针旋转0.089弧度，由图3b可见曲面Q1,Q3定位位姿得到改善；再沿x轴平移0.216 mm，曲面Q1定位精度得到进一步改善。由最大约束变换方向可知：在ICP算法中，精度较低的面Q3是主要定位目标，这与最小二乘法获得均衡误差分布的理论是一致的；本文算法中，具有最高精度的曲面Q1通过区域精度预估被赋予了更大的权因子，空间变换优先照顾Q1的定位。可见,加权算法通过精度预估获得区域定位权因子，从而控制空间变换方向，能够获得区域适应性更高的定位结果。

表2是理论定位、ICP定位和本文算法定位后各区域的误差数据。相对ICP定位，本文算法结果中，曲面Q1,Q2的均方根误差的精度分别提高42.5% 和35.5%，相应地，曲面Q3的精度下降12.7%；本文算法更接近理论结果(如图1)。从计算效率来看，本文算法迭代4次，计算效率是ICP算法(迭代1次)的71.74%。

与理论结果比较，ICP定位算法和本文定位算法的算法误差如表3所示，表中：“+”表示曲面误差被过高估计，“-”表示曲面误差被过低估计。ICP定位的算法误差与前文分析一致，高精度曲面Q1,Q2的误差分别被过高估计93.75%和47.06%，而精度较低的曲面Q3的误差相应地被过低估计了；3类曲面的本文算法误差分别是+12.50%、-5.88%和9.02%，本文优势明显。

表1 最大约束变换

表2 立方体3个面的定位结果比较 mm

表3 定位算法误差比较

4 整体叶轮精度评估的定位实验

本文实验在CPU 2.10 GHz和RAM 2.5 GB的电脑上用VS 2008软件编程实现。

如图4所示为整体叶轮精度评估的定位实例，叶片的压力面和吸力面的制造精度较高，而靠近叶根部位由于对气流影响不大而精度要求相对较低，轮毂轴端面、叶片前后缘等曲率变换大的地方精度最差。现根据工程要求将整体曲面按照精度从高到低分为3类Qj(j=1,2,3)，如图4a所示。图4b是叶轮的初始位置，蓝色曲面为理论数模，黄色为测量三角网格数据，可见叶片位置偏差较大。将17 068个测量数据点与数模进行直接ICP配准，获得如图4c所示的结果，易见叶片定位误差仍较大；采用本文提出的基于精度预估机制的自适应动态加权定位方法，结果如图4d，在制造精度要求较高的叶身处得到了较好的定位精度。如表4所示为ICP方法和本文方法定位结果中测量数据与理论模型的偏差值。比较可知：ICP算法得到的3类曲面的定位精度较均衡，低精度部位的数据使定位位姿偏离，在最小二乘意义下，整体上均方根误差达到最小0.53 mm；本文算法在精度较高的曲面Q1,Q2处定位精度分别比ICP算法相应提高了3.22%和5.66%；ICP定位算法效率稍高，综合考虑，本文方法更符合实际工程要求。

表4 整体叶轮定位结果比较

图5是叶轮3类精度曲面的定位优化过程，实时采用均方根误差和平均误差构造权因子，位姿不断得到优化。图5a和图5b分别是均方根误差和平均误差随迭代次数的收敛变换过程，可见所提算法可以逐步得到测量数据和理论模型的逼近，迭代超过40次后，算法趋于收敛。

5 结束语

考虑复杂曲面不同区域的精度差异性，基于区域精度预估机制进行的动态加权最小二乘法能够灵活地对测量数据进行自适应定位，在直接最小二乘定位下不满足要求的高精度区域可获得更优的定位结果。动态加权定位模型是连续的无约束非线性优化问题，不需要额外的精度约束，本文利用线性化方法将问题转化为四阶线性方程组求解，计算高效且具有鲁棒性。该方法能够根据误差分布自适应地控制测量曲面的定位位姿且计算高效，因此适用于复杂曲面大数据定位问题。本算法通过预估区域精度来构造该区域在定位中的影响程度。虽然加权定位方法取得了不错的定位效果，但区域精度与定位权系数之间的函数关系需要进行更细致地优化。另外，通过最近点算法进行预定位，为本文算法提供收敛的初始位姿，这是一种理想状态。在预定位算法中应该考虑区域数据量对均衡误差分布的影响，为本文算法提供可靠的收敛初始位姿。