基于随机测量点的机翼精加工位姿计算方法

中航工业成都飞机工业（集团）有限责任公司 苗 勇 何 凯 陈国强

西北工业大学现代设计与集成制造技术教育部重点实验室 余剑峰

机翼作为具有气动外形要求的关键部件，在其接头精加工及翼身对接过程中，为提高精加工质量，需要对机翼的空间位姿参数进行计算，以便指导机翼空间位姿的调整。文献[1-5]中，通过使用激光跟踪仪进行测量，获得部件上用于表达部件位姿的测量点的测量值，并结合这些测量点在产品理论姿态下的理论值，采用奇异分解法或最小二乘法，计算测量点与其理论值间平方和最小的转换关系，从而得到部件的空间位姿。这种方法已成功应用于机身的姿态计算，具有计算简单、精度高的特点，但由于这种方法在实际使用时，需要根据产品数模在产品上预先布置用于表达位姿的点，并对这些点在空间中的实际位置进行测量，使得其可能会受到来自产品结构以及加工平台结构等因素的限制而无法使用，同时，要实现对产品上确定点的测量，通常要使用到激光跟踪仪等先进测量设备，对测量设备的要求较高。例如，就机翼精加工而言，由于其产品对象的曲面封闭结构以及精加工型架的结构紧凑性，使得在机翼上不易布置位姿表达点，并使用激光跟踪仪对这些空间点进行测量，从而造成该方法难以使用。

翼面水平测量是机翼精加工工艺过程中必不可少的环节，能够起到检验位姿精度以及定性指导位姿调整的双重作用，但由于这种测量所得的点并非翼面上的预先布置的确定点，而是依赖于机翼实际位姿的随机点，此时，无法预知测量点的理论坐标值，因此，前文所述方法不能适用于该种情况下的机翼位姿精确计算。

为此，本文通过对翼面水平测量原理的分析，提出一种基于随机测量点的机翼精加工位姿计算方法。该方法旨在利用翼面水平测量数据，实现位姿参数的求解，以起到简化精加工工艺流程、提高精加工质量的作用。

1 原理分析

翼面水平测量数据是本文计算机翼位姿的基础，下面将通过翼面水平测量原理的分析，对本文机翼位姿参数计算的基本思想进行阐述。如图 1所示，Ω、Ω'分别表示具有理论位姿与实际位姿的翼面，OXYZ为全局坐标系。为描述机翼的空间位姿变化，在机翼部件上建立一个随动坐标系，该坐标系依机翼位姿改变而改变，如O1X1Y1Z1与O'1X'1Y'1Z'1分别对应于理论位姿和实际位姿。在翼面的设计模型上分布有一定数量的关键点，这些点称为水平测量点，如点Ai。点Ai的正下方设有直线位移传感器，完成对这些点的测量，其过程为：由直线位移传感器沿Z向运动，当接触到蒙皮后，可通过传感器读数得出测点的Z向高度值，结合直线位移传感器的位置坐标值（x，y）便可得到测点的空间坐标。由于误差的存在，机翼上架后，其实际状态并不处于理论姿态，如Ω'所示，因此，直线位移传感器所测得的高度值并非理论测点Ai的高度值，而是实际位姿翼面（即Ω'）上的另一个点B'i，将B'i称为“伪水平测量点”。考虑到误差影响的随机性，B'i可认为是翼面上随机测量的1个点。以表示测量点集（i=1，2，…，N），其中N为测量点数目，在不考虑翼面变形及测量误差的情况下，实测点集B'实际上是理论位姿翼面上的某一对应点集B={B1，B2，…，BN}经过旋转与平移得到，将B称为B'的匹配点集。点集B与B'间的旋转与平移关系能够反映出机翼的位姿信息，若能求得点集B，便能采用现有的位姿计算方法，求出机翼的实际位姿。

寻求“伪水平测量点”的匹配点的过程可简化为一个测量数据点集与理论曲面匹配的问题，这个匹配问题可描述为：将测量点集看作一个刚体，相对于曲面的CAD模型（即具有理论位姿时的翼面模型），经过坐标变换使得测量数据尽可能地包容该模型，并且使得测量数据与曲面间的距离误差最小，此时，理论曲面上与经过旋转的测量点最近的点就是原测量点的匹配点。在现有寻求匹配点算法中，以迭代最近点（Iterative Closest Point, ICP）[6-10]最为常用，该算法基于最小二乘法目标函数，其基本思想是假定测量点与曲面充分贴合时，点集中的所有点到曲面的最近距离平方和最小。国内外学者对ICP算法开展了较多应用研究[10-14]，用于自由曲面的误差评价及曲面加工余量的优化。本文结合机翼精加工工艺过程，通过采用ICP算法，寻求翼面水平测量实测点在翼面设计模型上的匹配点，以便为位姿计算提供完整的信息。

图1 水平测量原理图Fig.1 Principle diagram of level measuring

在求得匹配点后，便可进行位姿的求解。机翼位姿可采用6个参数表示，将其记为其中x1、x2、x3为刚体的欧拉角，表达刚体由当前姿态到理论姿态的一个变换序列，x4、x5、x6为刚体旋转后的平移矢量。分别表示实测点与其匹配点的空间坐标值，则这2个点与位姿参数间的关系可表示为：

式中，R为旋转矩阵，见式（2），其中有

2 位姿计算模型

2.1 位姿参数的简化

采用6个参数的刚体姿态表达，是一种通用的表达方法。然而，由于机翼精加工有其特定的工艺规程，这会导致机翼处于某种约束下，例如，机翼在位姿调整或进行精加工前，有一个预定位的过程，会采用精加工型架上的交点定位器对机翼进行定位，此时，机翼的实际位姿与理论位姿相比，将会很接近，只在某些位姿参数上与理论位姿参数存在差异。因此，通过结合机翼精加工工艺特点，可以有效简化位姿参数的求解。

以实际加工为基础，本文假定在机翼精加工预定位时，优先保证机翼部件上的某一交点孔（图2中红色圆圈），并采用短销进行精确定位。采用这种定位方式，既可以保证机翼的预定位精度，也可有效减少位姿求解参数，其定位效果等价于在机翼与精加工型架之间施加了一个固定点约束。将固定坐标系设在该固定点处，并假定具有理论位姿时，固连在机翼上的动坐标系与固定坐标系重合，那么，在该固定点约束下，机翼的自由度及位姿参数见图 2。由于存在固定点约束，3个平移自由度被限制，机翼实际位姿中只有3个旋转参数与理论位姿存在差异，因此，在机翼实际位姿求解时，只需考虑3个旋转参数，其实际位姿可表示为以下内容将以水平测量数据及翼面设计模型为基础，采用ICP算法与Levenberg-Marquardt算法，实现匹配点的寻求及3个参数的求解。

图2 固定点约束下的机翼自由度Fig.2 Freedom of aircraft wing under fixed point constraint

2.2 算法原理及步骤

设具有理论位姿的标准曲面（设计翼面）为S，测量点为Pi,0（0表示初始值，即由直线位移传感器测得的值），它到曲面S的最近距离点为Bi,0（i=1，2，…,N），其中，N为测量点的数量。第k次匹配后的点为Pi,k，为原测量点经过k次旋转所到的新位置，Bi,k为曲面上对应的最近点。由Pi,k、Bi,k所组成的2个点集间的距离可表示为

当dk值最小时，可认为达到最佳匹配，此时，所得的点Bi,k即为原测量点Pi,0的匹配点。根据假设，不考虑平移位姿参数，只需对测量点集进行旋转变换，其旋转矩阵为R（X）或R（x1，x2，x3）。在已知测量点集及标准曲面的情况下，求测量点集的匹配点集的算法步骤如下。

（1）算法 1（ICP）。

步骤 1:令k= 0，计算Bi,0和d0；给定精度ε＞0；

步骤2:k=k+1；

步骤3:求旋转矩阵Rk（X），使得如下函数

达到最小值，采用算法2求解；

步骤4:计算

步骤5:求解到曲面S的最近距离点

步骤6:计算目标函数

步骤7:令则停止计算；若则转步骤2，否则转步骤8；

步骤8:即为测量点的匹配点（i=1，2，…，N），为求得位姿参数X,再次调用算法2，使目标函数最小，停止。

以上算法中有2个判断停止计算的条件，如步骤7，分别为：（1）时，说明经过k次迭代后，Pi,k有远离曲面的趋势，此时算法不收敛，应该停止计算需要重新对机翼进行预算定位，提高机翼的预定位精度，并再次进行水平测量和调用ICP算法进行位姿评估，（2）说明经过k次迭代后，Pi,k已离曲面足够近，实现匹配点的搜索，经过步骤8的计算可求得机翼位姿。

在算法1中，每进行一次循环都需要调用算法2。给出算法2之前，参考文献[5]先对其目标函数进行前处理。算法2的目的是实现一个旋转变换的求解，使得2个点集的距离平方和最小，其目标函数为：

将式（4）展开，得到：

式中，的第j行所组成的行向量，式（5）表达的是一个非线性最小二乘优化问题，可采用Levenberg-Marquardt方法[15]，即L-M法进行求解。L-M算法中需要计算的梯度，为此，给出式（5）的梯度计算公式：

其中，

（2）算法2（L-M法）。

步骤1:取给定

步骤2:停算；

步骤3:求解方程组

步骤4:搜索步长，令是满足下面不等式的最小非负整数

步骤5:令转到步骤2。

在求解匹配点时，算法1可能会收敛于某一局部解，算法求解的初值X0对结果有着很重要的影响；此外，计算过程中，还需要反复调用算法2，X0对该算法的收敛速度也有较大影响。然而，通过2.1节的分析可知，机翼精加工之前存在一个预定位工艺过程，在此之后，机翼的位姿精度已经很贴近理论值，在这种实际工程背景下，本文取作为迭代初值，既能有效保证算法收敛于全局最优解，也可使算法效率不受大的影响。

3 测量点精度补偿

采用以最小二乘法为目标函数的ICP法寻求测量点的匹配点，实现位姿参数计算，具有目标函数简单、收敛速度快的特点。然而，实际加工过程中，由于变形的存在，测量点的实际位置与其理论位置间有偏差，导致测量时存在测量误差，并且各测量点的测量精度不一，例如，由于机翼沿翼展方向的变形累积及刚性变化，机翼外侧测量点的测量精度较内侧相比会有所降低。若采用存在误差的测量点与理论位姿曲面进行匹配，势必会导致误差，因为这些测量点并非理论位姿翼面上的点。通过分析算法原理及工艺流程可以得出，影响机翼位姿计算的误差主要包含2个部分（不考虑测量系统自身的测量误差）： （1）部件总装时引入的误差，这类误差包括零件制造误差、定位误差、铆接变形等，最终会反映到翼面外形上，称为局部外形偏差； （2）部件在精加工型架上完成预定位，并释放部分定位器，由于部件自重变形而引起的误差。如图3所示，Δi1、Δi2分别表示第一部分误差与第二部分误差，为理论上的实测值，为实际测量值，i表示第i个测点。从图3中可以看出，由于Δi1及Δi2的存在，使得并不重合。

Δi1可通过翼面外形测量预知其值大小。翼面外形测量在机翼完成部件总装后进行，主要针对翼面水平测量点，用于检验翼面的变形情况。Δi2与部件刚性有关，很难进行精确预测，可根据产品结构及以往产品的检验历史数据，大致分析。为减少Δi1及Δi2对翼面位姿计算的影响，提出测量点精度补偿量：

图3 测点误差示意图Fig.3 Diagram of measurement point errors

式中，Δi表示第i个测量点的精度补偿量调整系数，为第i个测量点的公差要求，用于反映Δi2。由测量原理分析可知，测点的x,y值取决于直线位移传感器的安装精度，其误差可以忽略，误差较大是其z向值，因此，只需对测点的z进行精度补偿。

4 算例分析

为验证算法的有效性，本文给出一个仿真算例，以一个圆柱面表示机翼翼面，固定点约束为坐标系原点，翼面上共有8个水平测量点，沿y轴对称分布，见图 4，各点坐标值其公差要求见表1。预先给翼面设定一个位姿单位为弧度。在不考虑误差，即Δ= 0的情况下，其实测点坐标值见表2。

得到水平测量点的实测点后，采用ICP算法可求得其对应的匹配点，见表2中的匹配点。利用表2中具有对应关系的两个点集，采用L-M算法，可求得部件位姿。Δ时，主要包括Δi1、Δi2两部分，其中Δi1由实际测量得到，可直接用于对测量数据的精度补偿，在此不予考虑。为仿真Δi2，在各测点上z向增加一个误差扰动，模拟机翼因自重产生的变形，并通过结合各测点公差要求及调整系数示意测点精度补偿方法的使用。且未进行补偿时，其测量点及其匹配点见表3；取进行补偿后，其测量点及其匹配点见表4，位姿计算结果见表5。

图4 位姿计算仿真模型Fig.4 Simulation model of posture calculating

表1 理论水平测量点坐标值

表2 Δ = 0时实测点及其匹配点

表3 Δ ≠ 0、未进行补偿时实测点及其匹配点

表4 Δ ≠ 0、进行补偿后实测点及其匹配点

表5 仿真计算结果

从表5中可以看出，当不考虑误差，即认为Δ=0时，位姿计算值与其理论值间的误差很小，几乎完全一致；当且未进行精度补偿时，位姿计算值的误差有所增大；进行精度补偿后，可提高位姿计算精度。因此，本文所述的位姿计算方法是一种有效的方法。

5 结论

利用机翼精加工的水平测量，采用迭代最近点算法（ICP算法）与L-M算法，实现机翼精加工位姿的求解。仿真算例表明，该方法是一种有效的机翼位姿计算方法。在采用该方法计算位姿时，应当注意以下2点。

（1）本文结合机翼精加工工艺特点，将机翼的6个位姿表达参数简化为3个。在实际应用时，可根据加工的实际情况，对待求解的位姿参数进行适当的调整。

（2）为提高位姿计算精度，需要对测量值进行精度补偿，其具体补偿量在一定程度上依赖于工程经验，后期将进一步研究测量值的精度补偿方法。

[1] 俞慈君, 李江雄, 余锋杰, 等. 带工程约束的点匹配算法.机械工程学报, 2010, 46(5):183-190.

[2] 张斌, 姚宝国, 柯映林. 基于鞍点规划理论的机翼水平位姿评估方法. 浙江大学学报, 2009, 43(10):1761-1765.

[3] 邱宝贵, 蒋君侠, 毕运波, 等. 大型飞机机身调姿与对接试验系统. 航空学报, 2011, 32(5):908-919.

[4] 柯映林, 杨卫东, 李江雄, 等. 基于激光跟踪仪的飞机机身姿态计算方法: 中国, CN200810121358, 4. 2009-02-11.

[5] 罗芳, 邹方, 周万勇. 飞机大部件对接中的位姿计算方法.航空制造技术, 2011(3):90-94.

[6] Besl P J, Mckay H A. A method for registration of 3-D shapes.IEEE Trasactions on Pattern Analysis and Machine Intelligence, 1992,14(2):239-256.

[7] Ahmad A, Christophe L, Mohamed D, et al. 3D registration using a new implementation of the ICP algorithm based on a comprehensive lookup matrix: Application to medical imaging. Pattern Recognition Letters,2007, 28: 1523-1533.

[8] Chetverikov D, Stepanov D, Krsek P. Robust Euclidean alignment of 3D point sets: The trimmed iterative closest point algorithm.Image and Vision Computing, 2005, 23(3):299-309.

[9] Liu Y H. Automatic registration of overlapping 3D point clouds using closest points. Image and Vision Computing, 2006, 24(7):762-781.

[10] Yahia H M, Huot E G, Herlin I L, et al. Geodesic distance evolution of surfaces: a new method for matching surfaces: a new method for matching surfaces. Computer Vision and Pat tern Recognition, 2000, 1(1):663-668.

[11] 武殿梁, 黄海量, 丁玉成, 等. 基于遗传算法和最小二乘法的曲面匹配. 航空学报, 2002, 23(3): 285-288.

[12] 刘元朋, 刘晶, 张力宁, 等. 复杂曲面测量数据最佳匹配问题研究. 中国机械工程, 2005, 16(12):1080-1082.

[13] 严思杰, 周云飞, 彭芳瑜. 大型复杂曲面零件加工余量均布优化问题研究. 华中科技大学学报, 2002, 30(10):35-37.

[14] 高栋, 孔令豹, 姚英学, 等. 光学自由曲面面形误差评定算法的研究. 哈尔滨工业大学学报, 2006, 38(10):1630-1632.

[15] 马昌风. 最优化方法及其Matlab程序设计. 北京:科学出版社, 2010.