基于PSO与WSVD的飞机部件位姿拟合方法

何晓煦，雷沛，潘登，杨阳，李现坤，邓珍波

1.航空工业成都飞机工业（集团）有限责任公司，成都 610092

2.四川省航空智能制造装备工程技术研究中心，成都 610092

现代飞机制造已广泛采用数字化调姿对合技术［1］以实现飞机部件姿态的自动测量和调整［2］，解决了传统人工作业方式存在的调姿精度差、效率低下、劳动强度大等问题［3］。飞机部件数字化调姿包含3个主要步骤：位姿测量、位姿拟合、位姿调整。位姿测量［4］是利用激光跟踪仪等大尺寸测量设备测量飞机部件上的调姿基准点在现场装配坐标系下的坐标［5］；位姿拟合是将测量数据和理论数据进行配准，拟合出部件在全机坐标系下的位姿参数；位姿调整是根据调姿机构运动学关系，将部件的位姿参数转化为数控定位器的运动量，驱动飞机部件达到目标位姿［6］。

部件位姿拟合是数字化调姿的核心环节之一，所用的算法直接影响调姿精度和效率［7］。当前较常用的方法可以分为非迭代算法和迭代算法两类。常用的非迭代算法包括奇异值分解法（Singular Value Decomposition， SVD）［8-9］、三点法、四元数法等。迭代算法包括最小二乘法（Least Squares， LS）［10-11］、随机抽样一致性算法（Random Sample Consensus， RANSAC）、迭代最近点法（Iterative Closest Point， ICP）［12］及其各类改进算法［13］等。

由于飞机部件装配环节存在各类累积误差，实际产品与理论数模不可避免地存在差异。为了控制部件的装配精度及对接后的整体外形精度，不同区域的基准点通常会设计不同的公差要求。传统的SVD算法虽然控制了所有基准点的综合转换残差和，但经常导致精度要求高的基准点超差而精度要求低的点还有较大的调整余量。Chen等［14］提出了一种加权位姿拟合算法，利用有限元分析不同区域基准点理论值的偏移量，偏移越大的点其权重设置越小，但权重的计算公式是固定的，适用范围会受到影响。朱绪胜和郑联语［15］提出了一种基于关键装配特性的部件位姿多目标优化算法，并结合粒子群优化（Particle Swarm Optimization， PSO）算法［16］对位姿参数进行求解，将位姿参数作为六维粒子进行迭代计算，但未考虑权重因素影响。陈远志等［17］提出了一种改进粒子群算法，引入惩罚函数处理基准点的容差约束，减少拟合残差，但未考虑权重影响。蒋睿嵩等［18］研究了一种基于权值约束的点云精确配准模型，通过对配准模型中较为重要的区域赋予较高的权重来提高该区域的配准精度，但未考虑其他约束条件。

针对以上算法的不足，结合飞机大部件调姿的工艺特点，提出一种基于粒子群优化与加权奇异值分解（Weighted Singular Value Decomposition， WSVD）的位姿拟合算法，根据基准点的精度要求赋予相应的权重，将所有基准点的转化残差均满足精度要求作为约束条件，将点的权重作为粒子进行迭代优化，求解大部件的位姿参数，并在不同型号的飞机大部件调姿对合系统中进行应用验证。

1 加权奇异值分解算法

1.1 飞机部件位姿描述

在飞机坐标系下，飞机部件的基准点理论位置坐标为Qi(Qi，x，Qi，y，Qi，z)，i∈(1，2，…，n)，n为基准点个数，由激光跟踪仪测量的实测位置坐标为Pi(Pi，x，Pi，y，Pi，z)。假设部件经过位姿调整量M变换后到达理论位置，则有

式中：Tx、Ty、Tz分别为x、y、z方向的空间平移量；α、β、γ分别为绕x、y、z方向的空间角度；R为3×3的空间旋转矩阵；T为3×1的空间平移向量。

飞机部件位姿调整量M可基于误差平方和最小化原则进行拟合，其中，LS算法求解的旋转矩阵由于没有正交性的约束，可能造成过分追求数据逼近而引入较大误差。相比之下，SVD算法在求解过程中利用了旋转矩阵的正交性而更符合实际情况［19］。

由于两组点具有相同的质心，可通过质心化简化计算过程，即先将平移向量T从式（2）中分离，优先求解旋转矩阵R。分别求出理论位置和实测位置的质心坐标分别为

式中：n≥3。质心化后的坐标为

根据式（1）～式（6），有

1.2 奇异值分解法

1）SVD优化配准算法

采用SVD进行优化配准计算时，质心化后误差平方和最小的目标函数为

式中：sR为SVD求得的旋转矩阵。则有

令

则上述问题等价于求max(trace(sRsH))，trace(·)为矩阵的迹。将sH进行奇异值分解得

式中：sU和sV为正交矩阵；sΛ为对角矩阵。则由SVD解算的sR的最优值为

基准点的平移向量为

由SVD解算的位姿调整量中的平移向量为

2）SVD的位姿拟合误差

部件经过位姿调整后，基准点坐标值为

调姿后基准点的残差为

sεi应满足飞机大部件基准点容差要求εi，d，d∈(x，y，z)。

SVD的目标函数是基准点拟合误差的平方和（见式（9）），但飞机部件装配过程中，对不同基准点的容差要求不同，这导致基于SVD的姿态拟合算法虽然能够解算得到总残差最小的部件位姿，但却存在部分基准点超出容差范围的问题，而且往往是那些精度要求高、容差范围小的关键基准点。

1.3 加权的奇异值分解法

1）WSVD优化配准算法

针对SVD在飞机部件调姿方面的缺陷，提出一种能够自适应基准点容差要求的姿态拟合算法。在位姿调整量的解算过程中引入权重，利用权值控制关键基准点误差，可以提高位姿调整量解算的准确度，从而提高飞机部件调姿的合格率。

假设第i个基准点x、y、z坐标对应的权重值相等，均为wi。将权重代入误差平方和最小的目标函数中，此时目标函数表示为

式中：wR为WSVD求得的旋转矩阵。与SVD求解方法类似，将式（18）展开得

将求解最小化问题minwΕ转换为求解

则上述问题等价于求max(trace(wRwH))。对wH进行奇异值分解得

式中：wU和wV为正交矩阵；wΛ为对角矩阵。则wR的最优解为

先求出基准点加权后的平移向量为

则WSVD位姿调整量中的平移向量为

2）WSVD的位姿拟合误差

飞机部件调姿后，基准点的三维坐标为

用wεi表示调姿后的残差为

wεi应满足飞机大部件基准点容差要求εi，d。

引入权值进行加权的奇异值分解，能较好地分配飞机部件基准点的残差在最小化问题的权值，使结果趋近于权重值更大的基准点，即偏向装配精度要求更高、容差范围更小的基准点，从而避免这些基准点超出容差范围的问题，使求解出的部件位姿调整量更贴近实际调姿需求。

2 基于粒子群优化的权重分配算法

2.1 粒子群优化算法

1）PSO目标函数

PSO具有较强的鲁棒性和全局寻优能力，能保持个体间持续的信息交互，并且易与其他算法进行结合［20］。利用PSO自动搜寻能使拟合后所有基准点满足容差要求的权重值，优化目标是基准点在x、y、z方向的超差项的残差尽可能小。因此定义目标函数适应度值E为拟合误差超差的基准点，其残差的绝对值之和，即

式中：为拟合误差超差的基准点在d方向的残差的绝对值。

利用PSO迭代求解基准点权值时，将待求解的权值作为粒子。粒子群每次更新后，权值的改变将引起WSVD拟合结果中残差的变化，从而使目标函数E中的超差项随粒子群的更新而动态调整，引起PSO适应度值更新，进而指导粒子群的搜索方向，因此具有自适应的寻优能力。

2）多维粒子群算法

每个粒子包含一个位置向量和一个速度向量。基准点对应的权值向量为w，则粒子的维度为n。假设选取m个粒子，则每次迭代时粒子群存在n×m维的位置矩阵x，以及n×m维的速度矩阵v。

粒子群初始化时，将位置矩阵x初始化为每项均为1的单位矩阵。则初始位置因子=1，速度矩阵v初始化为每项为[0， 1]的随机数矩阵。

粒子搜索解空间时，每个粒子根据自身的适应度值保存各自搜索到的个体最优经历位置pp，并在个体最优经历位置中选取群体最优经历位置pg。粒子群按式（28）调整位置矩阵：

式中：c1和c2为加速因子，取c1=c2=2；ξ和η为[0， 1]中均匀分布的随机数；λ为惯性权重因子，采用线性递减策略对其进行赋值，有

每个粒子的权值总和

式中：λmax和λmin分别为λ值的上下限，取λmax=0.95，λmin=0.4；a为当前迭代次数；A为总迭代次数，A=50。

粒子群的位置矩阵更新为

式中：δ为位置更新时的约束因子，取δ=1。

为提高算法的寻优能力和收敛速度，在粒子每次更新后，位置矩阵x应满足非负性条件。否则可能导致粒子朝相反方向搜索造成非全局最优解，或搜索时间被无限放大。当位置因子xi，k不满足非负性条件时，将更新后的第k个粒子的第i个基准点的权值xi，k重新初始化为1，即

同时，在PSO迭代过程中，位置矩阵x应满足归一性条件。通过归一化处理来保证每个粒子的权值总和为n不变，有：

此，每个基准点权值的空间范围为[0，n]。

2.2 本文方法计算流程

本文方法的计算流程如图1所示，首先用基准点进行SVD拟合，当不满足算法终止条件时，进入PSO部分，粒子根据目标函数适应度值，在约束空间内按相应的规律搜索，不断调整自身的位置矩阵和速度矩阵，每次得到更新的权值后，再带入WSVD部分进行加权奇异值分解，从而更新PSO目标函数适应度值，继而指导粒子群的搜索方向，最后在满足算法终止条件的情况下停止计算。具体步骤为

图1 本文方法计算流程Fig.1 Calculation process of this method

步骤1将理论位置Qi和实测位置Pi用SVD拟合，求得位姿调整后的残差sεi（见式（3）～式（17）），判断是否满足容差要求。满足进入步骤10，反之进入步骤2。

步骤2初始化粒子群。利用步骤1中的残差计算PSO适应度值E（见式（27）），得到个体最优位置pp和群体最优位置pg。进入步骤3。

步骤3更新速度矩阵v和位置矩阵x（见式（28）～式（32））。进入步骤4。

步骤4将更新后的位置矩阵x的每个粒子分别作为权值w代入WSVD拟合，求得位姿调整后的残差wεi（见式（18）～式（26））。进入步骤5。

步骤5根据步骤4的残差wεi求得每个粒子的适应度值E（见式（27）），更新粒子个体最优位置pp和群体最优位置pg。进入步骤6。

步骤6将更新后的群体最优位置pg作为权值，判断对应的WSVD分解后的残差wεi是否满足容差要求。满足进入步骤9，反之进入步骤7。

步骤7判断当前是否达到最大迭代次数A。满足进入步骤8，反之进入步骤3。

步骤8输出步骤6中的群体最优位置pg作为权值w，以及对应的WSVD的位姿矩阵，该结果为当前迭代条件下，最接近容差要求的解。如需要继续寻找满足容差要求的解，可以通过增加粒子数或迭代次数，或者调整迭代参数后，从步骤1重新开始计算。

步骤9输出步骤6中的群体最优位置pg作为权值w，以及对应的WSVD的位姿矩阵。找到了满足容差要求的解，计算结束。

步骤10输出SVD的位姿矩阵，权值w每项均为1。找到了满足容差要求的解，计算结束。

可见，本文方法在SVD无法满足容差的情况下才启用PSO+WSVD自动寻优，因此既具有常规SVD数据处理量小的优点，又具备PSO+WSVD全局搜索和自适应权重分配的能力。

3 案例验证与分析

3.1 C919机头上下部调姿对合

为了验证本文方法的可行性，以C919机头上部和下部的调姿对合过程为算例进行分析，以机头下部作为基准部件，机头上部作为调姿部件与机头下部对接。分别使用SVD方法和本文方法计算出机头上部的位姿调整量为M1和M2，若机头上部经过M1调姿后的残差超差，而经过M2调姿后的残差在容差范围内，则说明本文方法具有可行性。C919机头上部有7个基准点，理论值和实测值如表1所示，根据实际装配中的工艺要求，基准点对应的三维容差约束为± 1 mm。

表1 基准点位置信息Table 1 Location information of datum points

用1.2节的SVD方法进行计算。计算时间为0.005 s，位姿调整量为M1=［4.056 mm，5.279 mm， 8.286 mm， 0.011°， 0.006°， -0.013°］。机头上部经过M1调姿后，基准点的残差如表2所示。调姿后仍有1个超差项：2号基准点的y方向，残差为-1.235 mm（超差-0.235 mm）。因此，SVD方法无法使所有基准点的残差都满足容差要求。

表2 基准点残差比较Table 2 Datum points error comparison

在实际工程应用环境下，在使用M1自动调姿后，仍需要人工进行干预，工人重新测量基准点的实际位置，凭经验进行手工微调。据现场统计，为使部件调整到符合对接精度要求的姿态，需要人工手动微调1～5次，花费时间长达数十分钟到数小时不等。人工干预环节对工人技能和经验要求高，较大程度地制约了部件调姿对合的效率。并且，仍存在无法调整合格的风险。

用2.2节的本文方法进行计算。PSO的粒子数为30，机头上部有7个基准点，对应粒子维度为7，因此每个基准点权值的空间范围为[0， 7]。PSO的初始适应度值为1.235 mm，经过第1次迭代，适应度值降为1.043 mm，经过第2次迭代降为0，计算时间为0.144 s。基准点权值为w1=［0.779， 1.646， 0.858， 0.746， 0.627，1.265， 1.079］。位姿调整量为M2=［4.117 mm，5.686 mm， 8.334 mm， 0.009°， 0.007°， -0.015°］。

机头上部使用M2调姿后，基准点的残差如表2所示，残差均满足容差要求。由表2可以看出，本文方法通过权值的自动寻优，使2号基准点y方向的对接误差减小为-0.962 mm，对比SVD方法提高了0.273 mm。本文方法不仅使SVD方法的超差项满足容差要求，并且保证其他基准点残差均在容差范围内，避免了人工微调环节带来的效率损失和无法调整到位的风险，实现了完全的自动调姿。

在SVD方法和本文方法计算时间都很短的情况下，本文方法能够解决SVD方法局部容差超差的情况，证明了本文方法的可行性。

3.2 实例分析

为说明本文方法在其他机型和部件调姿中的有效性，采用文献［21］的实测数据作为算例进行分析。只考虑文献［21］的后机身调姿后，基准点的残差是否满足容差要求，如果满足则认为后机身能够实现与前机身对接，本文方法有效。根据文献［21］，后机身8个基准点的理论值和实测值如表3所示，基准点对应的三维容差约束为±1 mm。

表3 文献［21］基准点位置信息Table 3 Location information of datum points for Ref.［21］

按照± 1 mm的公差要求分别采用SVD方法和本文方法进行计算，两种方法的残差如表4所示。从表4可以看出，两种方法计算结果一致，且基准点残差均在容差范围内。由于SVD方法就能使后机身调姿后的基准点残差满足± 1 mm的容差要求，因此本文方法在通过SVD求出位姿调整量后就终止了计算，没有启用PSO+WSVD。此时位姿调整量为M3= ［-246.554 mm，212.303 mm， 8.260 mm， -0.569°， 2.496°，-1.249° ］，最大的残差值为-0.871 mm，为4号点的x方向。

表4 文献［21］基准点残差比较Table 4 Datum points error comparison for Ref.［21］

假设公差要求提高到± 0.65 mm，分别采用SVD方法和本文方法进行计算。

使用SVD方法进行计算。计算时间为0.002 s，求得位姿调整量仍为M3，基准点残差如表4所示。由表4可以看出，SVD方法在容差要求为± 1 mm和± 0.65 mm两种情况下计算结果一致，这是由于SVD是以控制总残差和最小为目标进行计算的，容差的改变并不影响计算结果。然而，这样的计算结果，在容差要求提高后，就出现了超差项：4号点的x方向，残差为-0.871 mm（超差-0.221 mm）。此时，SVD方法无法使所有基准点满足± 0.65 mm的容差要求。

使用本文方法进行计算。PSO的粒子数为30，后机身8个基准点对应粒子维度为8，因此每个基准点权值的空间范围为［0， 8］。PSO适应度曲线收敛图如图2所示，初始适应度值为0.871 mm，经过第1次降迭代后，适应度值降为0.790 mm，经过第2次迭代后降为0.736 mm，经过第3次迭代后降为0，计算时间为0.164 s。基准点权值为w2=［0.201， 0.600， 1.300，2.091， 0.353， 1.546， 0.743， 1.165］。位姿调整量为M4=［-246.439 mm， 210.657 mm，8.936 mm， -0.561°， 2.499°， -1.245° ］，基准点的残差如表4所示。

图2 PSO适应度值曲线收敛图Fig.2 PSO fitness curve convergence diagram

由表4可以看出，本文方法通过权值的自动寻优，使4号点x方向的对接误差减小为-0.644 mm，对比SVD方法提高了0.227 mm。本文方法不仅使SVD方法的超差项满足容差要求，并且保证其他基准点残差均在满足± 0.65 mm的容差要求。据此，验证了本文方法在其他机型和部件调姿中应用的有效性。

4 结 论

针对飞机大部件调姿中的位姿拟合问题，提出了一种基于粒子群优化结合加权奇异值分解的方法，将不同精度要求的基准点赋予不同的权重，权值的分配由粒子群算法进行迭代优化，并将每个基准点的转换残差均满足精度要求作为优化的约束条件。通过飞机数字化装配中的应用分析验证了算法的有效性。