改进的粒子群算法求解飞机位姿评估问题*

陈远志，黄 杰，章易镰，朱利民

（1.上海交通大学机械工程学院，上海 200240；2.上海拓璞数控科技股份有限公司，上海 201111）

飞机数字化装配技术相较于传统的固定型架装配方法，拥有高效、准确、经济和柔性等特点，因此被国内外飞机制造厂商广泛采用于制造各类先进的军民机[1]。飞机数字化装配主要包含3 个阶段：位姿测量、位姿评估、位姿调整。位姿测量是指采用激光跟踪仪、iGPS 等数字化测量系统测量固定在飞机上的调姿基准点在装配坐标系中的三维坐标值[2–3]；位姿评估是指比对调姿基准点的测量坐标值（测量系统提供）与理论坐标值（理论数模提供），从而计算出飞机部件位姿移动参数；位姿调整是指根据运动学关系，将飞机部件的位姿参数转化为数控定位器的运动参数，由数控定位器驱动飞机部件到达目标位姿。

位姿评估是求解点匹配问题在飞机装配领域的应用[4]，大量学者针对此问题进行了研究。奇异值分解法[5]、正交矩阵法[6]、四元组法[7]等是位姿评估常见的算法。上述方法是非迭代的，因为使用方便，已广泛应用于飞机装配现场。

此外，几种以非迭代算法获得位姿参数初值的迭代的位姿评估算法获得了更高的匹配精度。例如Li 等[8]提出基于迭代最近邻算法（Iterative Closest Point，ICP）的位姿评估算法，将姿态评估问题转化为空间点与曲面的配准问题、寻找点与点的对应关系，其内部采用奇异值分解（SVD）法求解；俞慈君等[4]提出一种考虑对称、共面、共线等工程约束的刚体匹配算法，利用SVD 法获得位姿参数初值，带入加权最小二乘牛顿迭代法求解非线性最优化目标函数，以获得刚体匹配最优位姿参数；张斌等[9]提出基于鞍点规划的位姿评估算法，综合利用激光跟踪仪和位移传感器进行为此测量，运用SVD 法和单纯形法求解飞机位姿评估问题。

随着人工智能技术的发展，群体智能优化算法被广泛研究并应用到各类最优化问题中。群体智能优化算法是一种生物启发式优化方法，主要包括蚁群算法、粒子群优化算法（Particle Swarm Optimization，PSO）、人工蜂群算法等[10]。位姿评估问题可以转化为最小二乘的优化问题，因此可以采用群体智能优化算法进行求解。PSO 算法具有结构简单、易于理解、计算量小的优点，因此已有学者将其应用到位姿评估的研究中。朱绪胜等[11]提出了一种综合考虑卫星舱段装配过程中的多个控制目标的多目标位姿参数优化算法，利用PSO 算法求解多目标优化函数，仿真结果证实了此算法的可行性；Hou 等[12]设计了一种简易机械测量装置测量飞机位姿，并采用PSO 和ICP 算法进行位姿参数求解。

上述针对位姿评估的研究均取得了理想的效果，但是当安装调姿测量点的位置的飞机部件制造误差或安装误差较大时，在按照位姿评估算法计算出的位姿参数调整飞机位姿后，会存在某一点或某些点的误差较大甚至远远超过其他的点。这是因为测量系统测量出的该调姿基准点在飞机部件局部坐标系中的实际坐标与数模定义的局部坐标系中的理论坐标相差较大，而上述文献主要考虑的是调姿测量点误差总和或者误差平方和最小，当所有调姿测量点均在合理坐标范围内时，上述文献的算法效果较好，但当调姿测量点不在合理坐标范围内时，会存在某一点或某些点的误差较大甚至远远超过其他的点。综合考虑调姿测量点最大误差以及误差总和的研究未见报道。本文针对这种情况，在基本PSO算法的基础上，提出了一种综合考虑最大误差以及误差总和的改进粒子群算法，首先以误差总和最小为优化条件，利用基本PSO 算法求解位姿参数，再通过动态约束优化的PSO 算法调整位姿参数，以使调姿基准点最大误差满足条件。

1 飞机数字化装配系统

飞机数字化装配系统是典型的测量辅助装配（Measurement Assisted Assembly，MAA）系统[13]，为了便于阐述此类系统，以较为常见的某型商务客机的机头中机身数字化装配系统为例加以说明，如图1所示。数字化装配系统首先通过激光跟踪仪测量安装在前机头部件和中机身部件上的调姿测量点的三维坐标，其次通过比对调姿测量点的当前坐标与理论坐标计算出飞机部件位姿移动参数，然后通过数控定位器驱动飞机部件到达目标位姿，最后再次通过激光跟踪仪测量调姿基准点是否符合容差要求，若不符合则重复上述步骤，否则调姿完成进入后续装配工作。

在飞机数字化系统中通常存在3 种坐标系：全局坐标系（Global Coordinate System，GCS）、测量坐标系（Measurement Coordinate System，MCS）以及与移动装配体固联的局部坐标系（Local Coordinate System，LCS），具体定义可参阅文献[11]。通过增强的系统参考点（Enhanced Reference System Point，ERS Point）和调姿基准点的测量，3 种坐标系可以互相转换，这种转换可以通过专业的大尺寸测量软件Spatial Analyses（SA）自动完成。为简化模型推导过程，省略各点的坐标值和各部件的位姿参数在不同坐标系之间的转换关系的推导，将它们都定义在GCS 中。

调姿基准点固定安装在待装配的飞机部件上，坐标值为调姿基准点的可以表示为：

式中，P为调姿基准点；Type 表示坐标值的类型，m代表调姿基准点的个数。5 种坐标值类型的关系在下文详述。需要特别说明的是，因为激光跟踪仪的测量精度所限，调姿基准点的测量值与真实值存在微小误差（Ux，y，z= 15μm+16μm/m），远小于飞机部件的制造误差和调姿基准点安装误差，所以下文中不特殊说明时调姿基准点的测量值即为真实值。

图1 飞机数字化装配系统示意图Fig.1 Schematic of aircraft digital assembly system

2 位姿评估模型的构造

在飞机数字化装配过程中，常把飞机部件作为刚体进行处理。位置和姿态是确定刚体之间关系的关键几何特征，它们描述了刚体在GCS 中的几何状态，在进行模型构建之前需要预先确定。刚体在GCS 中的位姿移动参数可以用六维矢量X=（α，β，γ，tx，ty，tz）T表示，其中α、β、γ表示固联在刚体上的LCS 在GCS 下的旋转角R、P、Y（Roll，Pitch，Yaw），tx、ty、tz表示LCS 的原点在GCS 下的位移量。

刚体运动可以由3×3 的矩阵R表示旋转

式中，cx=cos（x），sx=sin（x）。

位移矩阵T为3×1。

调姿基准点的估计值可以表示为

因为R和T是X的函数，所以位姿评估问题即是求解最优的X，使得估计值PiEs与理论值PiTh的误差总和最小且误差值PiEr小于容差值PiTo。其中，Th 代表理论值（Theoretical Value）由CAD 模型定义；M 代表测量值（Measured Value），由激光跟踪仪测量得到；Es 代表估计值（Estimated Value），是理论计算得到的调姿完成后的调姿基准点坐标值；Er 代表误差值（Error Value），是估计值与理论值的欧式距离；To 代表容差值（Tolerance），根据飞机设计要求确定。误差值PiEr可以表示为

则位姿评估问题的优化目标可以表示为

式中，μi是归一化参数，保证容差值较小同时Pi=[xi，yi，zi]T位置精度更高的调姿基准点不被轻视。

在飞机调姿问题中，其中最常见的约束条件是点位约束，可以表示为

即调姿完成后，所有调姿测量点的坐标误差值均符合容差要求。不满足约束条件意味着飞机部件调姿完成后无法达到工艺要求，所以要极力避免超差的存在。

3 基于PSO 算法求解位姿参数

3.1 基本PSO算法

PSO 算法是受鸟群觅食行为的启发，进而建立的群体智能优化模型。式（6）构成的位姿评估模型是一个典型的最优化问题，可以采用PSO 算法求解。

求解时将位姿移动矢量X作为粒子，其独立参数α、β、γ、tx、ty、tz作为粒子的维度，式（6）作为适应度函数。粒子群由N个位置为Xi、速度为Vi的粒子构成，Xi和Vi是在可行域范围内的六维向量。首先随机设置初始位置Xi0和初始速度Vi0；随后根据适应度函数（式（6））计算所有粒子的适应度f（Xi0）；再通过比较适应度，记录单个粒子最佳位置Pibest和群体最佳位置Pgbest；接着根据式（9）和式（10）更新粒子的速度和位置。

式中，Xik和Vik分别是第i个粒子第k代时的位置向量和速度向量；c1和c2是加速因子；r1和r2是0 到1 之间的随机数；w是惯性权重。

3.2 求解约束优化问题的PSO算法

通过上述基本PSO 算法或者奇异值分解法、牛顿欧拉法等求解飞机位姿评估问题，将求解出的位姿移动参数应用到飞机数模上进行仿真，时常出现两种问题：（1）调姿基准点最大误差（Max Error of OTP，MEO）不符合容差要求；（2）MEO 符合容差要求，但远超其余调姿基准点的误差。

为了解决问题（1），本文提出一种求解约束优化问题的PSO 算法（PSO for constrained optimization problems，PSO-COPs），用于求解带点位约束的飞机调姿问题。采用惩罚函数作为约束处理机制，将式（8）的点位约束条件转化为惩罚函数：

式中，σ=αk是指数型惩罚因子，α是大于1 的加速参数，k是优化代数，σ随着优化的进行而逐渐增大，使粒子Xi受到更强的约束力；φ是放大因子取10，因为调姿基准点的通常安装在刚度和精度较好的飞机部件上，所以PiEr–PiTo通常为负数，只有极少数情况为较小的正数，因此采用φ放大两者的差值，增强惩罚效果。

通过上述约束处理机制，可以将由式（6）和式（8）组成的约束优化问题转化为由式（12）描述的无约束优化问题。

在基本PSO 算法基础上，将式（12）作为新的适应度函数，用于求解带点位约束的飞机调姿问题。

3.3 改进的PSO算法

为解决3.2 节所述问题（2），本文提出一种改进的PSO 算法（Improved Particle Swarm Optimization,IPSO），以调姿基准点误差总和以及最大误差为优化目标，动态改变点位约束条件，形成嵌套循环优化模型。将式（8）中固定容差值的PiTo改为动态约束参数（i表示调姿基准点编号，j表示外层循环次数），惩罚函数改为：

式中，PiDC的初值为PiTo，动态更新规则下文详述。IPSO 算法的优化目标改为

IPSO 算法流程如图2所示，具体步骤如下。

Step1：输入CAD 模型定义的调姿基准点理论值PiTh、根据工艺要求确定的容差值PiTo以及由激光跟踪仪测量得到的测量值PM。

Step2：生成N个离散的粒子Xi，每个粒子拥有6 个自由度（α、β、γ、tx、ty、tz）。

Step3：在粒子位置和速度的可行域范围内，初始化粒子的位置Xi0和速度Vi0。

Step4：根据式（14）计算Xi的适应度，记录单个粒子最佳位置Pibest和群体最佳位置Pgbest。

Step5：以式（9）和式（10）更新粒子的速度和位置，若粒子速度超出可行域，则速度更新为

若粒子位置超出可行域，则位置更新为

同时将粒子位置超出可行域的维度相对应的速度维度取反，例如Xik中α的数值超出可行域，则Vik中α的速度取反，以此类推。

Step6：判断内层循环终止条件，即群体最佳适应度在T=50 代内变化小于预定义的最小变化MINVAR或者迭代次数超过最大迭代次数MAXGEN=1000。满足条件则进入步骤7，反之重复步骤4、5、6。

Step7：判断外层循环终止条件，即外层循环运行次数大于等于两次并且本次循环优化结果的MEO 大于前次的MEO。满足条件则进入步骤9，反之进入步骤8。

式中，jd是第j次外层循环的衰减系数

式中，dmin=0.8，dmax=0.99。

Step9：输出移动位姿参数R，T以及评估结果。

图2 改进的PSO算法流程图Fig.2 IPSO algorithm diagram

4 算例分析与蒙特卡洛仿真

以图1中的中机身调姿定位为例，验证算法的有效性。为了简化计算模型和便于理解，将中机身上左右两边共8 个调姿基准点简化为一个正方体模型上的8 个顶点，如图3所示。

4.1 算例分析

仿真算例中的模拟测量坐标值由式（19）生成

图3 算例模型示意图Fig.3 Calculation model diagram

表1 算例模型调姿基准点坐标值Table 1 Example model appearance benchmark coordinate values

采用本文提出的PSO-COPs 和IPSO 求解飞机位姿评估问题，并利用求得的位姿移动参数计算调姿基准点评估值PiEs，然后计算各点到相应的理论值PiTh的距离误差和以及调姿基准点最大误差。根据杨维等[14]描述的粒子群算法参数选取方法，为了避免优化过程中陷入局部极小点，取粒子群算法的w=1，c1=1.5，c2=2。将本文提出的两种方法与基本PSO 以及SVD 方法比较，如表2所示。

由表2可见，4 种方法的距离误差和均小于随机误差和Σ‖Ei‖。因为PSO 方法直接以距离误差和最小为优化目标，所以PSO 方法的距离误差和最小，但最大误差值最大，超过了调姿基准点容差值PiTo。在PSO 基础上改进的PSO-COPs 的距离误差和略大于PSO 方法，但因为采用了惩罚函数作为约束处理机制，其最大误差值达到了容差要求，在在距离误差和与最大误差值方面的优化结果均好于SVD 方法。IPSO 因为不断强化约束条件，所以获得了最小的最大误差值，是4 种方法中唯一最大误差值小于max‖Ei‖的方法。

4.2 蒙特卡洛仿真

在实际飞机位姿评估问题中，‖Ei‖≤PiTo不是总能满足，偶尔会出现个别调姿基准点制造误差超过设计要求。采用蒙特卡洛方法仿真在不同预设最大误差情况下，比较上述4 种方法的优化效果，以对本文提出的两种方法的可行性进行验证。蒙特卡洛法是一种数值模拟方法，需要产生大量随机样本，通常样本量越大，模拟结果越准确。因此在δ=0.1，0.2，…，2.0 这20 种情况下，每种情况均根据式（19）生成1000 组模拟测量坐标值，然后用4 种方法求解。各种情况下4 种方法的1000 次仿真平均最大误差值和平均距离误差和如图4所示，其中图4（b）中的最大误差修正率可以表示为

类似的，图4（d）中的误差和修正率可以表示为

图例中的MAXERROR 和SUMERROR 分别表示最大随机误差max‖Ei‖和随机误差和Σ‖Ei‖。

表2 SVD、PSO、PSO-COPs及IPSO距离误差和以及最大误差比较Table 2 SVD,PSO and PSO-COPs and IPSO distance error and error comparison

图4 蒙特卡洛仿真结果Fig.4 Results of Monte Carlo simulation

由图4（a）和图4（b）可见：预设最大误差小于等于0.8mm 时，SVD 方法的平均最大误差值最小；当预设最大误差大于0.8mm 时，IPSO 方法的平均最大误差值最小，并且IPSO方法的最大误差修正率开始上升，直到预设最大误差达到1.2mm后开始稳定在20%附近；PSO-COPs 方法的最大误差修正率的趋势与PSO 方法类似，超越SVD 方法的转折点是预设最大误差1.3mm，最大误差修正率从上升到稳定的转折点是预设误差达到1.6mm；PSO 方法的最大误差值始终最大、最大误差修正率始终最小。4 种方法的蒙特卡洛仿真1000 次平均最大误差超过容差值PiTo的转折点为SVD 1.4mm，PSO 1.3mm，PSO-COPs 1.5mm，IPSO 1.7mm。

由图4（c）和图4（d）可见：4 种方法求解飞机位姿评估问题得到的距离误差和均小于随机误差和Σ‖Ei‖；PSO 方法的距离误差和在本文的仿真条件中始终有最小的距离误差和，且距离误差和修正率保持在17%附近；SVD 方法的距离误差和修正率保持在15%左右，并且与PSO 方法的距离误差和的变化趋势保持一致；PSO-COPs 方法和IPSO 方法在预设误差小于等于0.7mm 时，因为各调姿基准点随机误差较小，惩罚函数并未达到触发条件，所以两者的距离误差和修正率与PSO 方法一致，因为此时三者的优化目标是一致的；当预设最大误差大于0.7mm 后，惩罚函数被激活，PSO-COPs 方法和IPSO 方法的优化目标发生变化，所以距离误差和修正率开始下降，而最大误差修正率则开始上升。

5 结论

（1）本文采用惩罚函数方式将点位约束条件引入到飞机位姿评估问题中，因为常规的优化方法容易陷入局部极小点，而PSO 方面能获得更接近全局最优的解，所以本文采用求解约束优化问题的粒子群算法求解。将带点位约束的飞机评估问题转化为带惩罚函数的优化模型，其他复杂约束也可以采用此方法添加到优化模型中，采用粒子群算法求解。

（2）在本文求解约束优化问题的粒子群算法基础上提出了一种改进的粒子群算法，采用动态改变约束参数来加强约束条件，可以通过迭代方式进一步减小飞机调姿基准点最大距离误差。但是因为PSO 方法是迭代的，IPSO 方法是嵌套迭代的，所以在算法效率上，SVD 方法最快，PSO、PSO-COPs 次之，IPSO最慢。

（3）通过一个算例和蒙特卡洛仿真，可以看到相比于传统的SVD 方法和新近提出的PSO 方法，本文提出的两种方法在优化最大距离误差方面优于SVD 方法和PSO 方法，在调姿基准点本身的误差较大时仍可以达到较好的优化结果，使最大距离误差仍满足点位约束条件。