飞机等大型机械产品的生产制造，常常采用分段制造最后再大部件组合装配的模式。大型产品常以预设的全局坐标系作为统一的装配基准，各个大部件也在此统一的坐标系下进行生产制造。为了提高飞机装配质量，现代飞机装配技术大量使用如激光跟踪仪、iGPS等高精度的大尺寸测量设备来提高测量精度。为了建立测量设备测量坐标系与装配坐标系的关系，通常需要通过测量固定在装配现场的增强参考系统点(enhanced reference system, ERS)进行转换参数计算，这个过程称为转站测量。根据ERS点已知的名义坐标计算当前测量坐标系在全局坐标系中的位姿，将测量设备的坐标系转换到全局坐标系下。测量坐标系与装配坐标系的坐标转换参数称作转站参数，包含旋转矩阵以及平移矩阵。由于工件的尺寸较大，环境温度变化引起的热膨胀变形明显，这将对飞机工件姿态评价及对接精度产生较大影响。为了减少这种影响，研究工件三维膨胀变形下激光跟踪仪转站参数优化，从而提高跟踪仪设备的转站精度以及工件的测量精度。

求解转站参数通常采用点匹配算法，包括奇异值分解法^[1]、正交矩阵法^[2]、双四组元法^[3]等。在大型产品工件的装配过程中，整个测量系统包含跟踪仪测量坐标系、工装坐标系以及产品工件坐标系等。需要经常使用激光跟踪仪进行转站，从而获得测量坐标系与被测坐标系的转站参数(包含旋转矩阵以及平移矩阵)。大型产品工件由于装配的材料众多、结构复杂以及受工装约束变形等影响，产品在x、y、z 3个方向上的热膨胀变形系数会不一样。而这3个不同方向上的热膨胀变形系数也会改变转站参数的结果。

管官等^[4]提出了一种考虑船舶工程约束的船体分段测量点集自动匹配算法，利用权值向量实现对不同方向上精度要求的误差分配，求解包含水平度、垂直度、平面度、直线度等工程约束的目标函数优化模型进而得出最佳匹配结果。俞慈君等^[5]把两点对称、多点在同一平面上、多点在同一直线上等工程约束引入到点匹配的优化目标函数中，利用Newton非线性最优化方法求解6个姿态参数，利用权值矢量实现对多个约束进行同时控制的误差分配求解。Predmore^[6]根据测量点不确定度椭球误差模型，采用基于参考点加权的马氏距离(Mahalanobis distance)最小的转站方法来提高转站精度。王青等^[7]提出了一种基于容差约束的机翼位姿评价及优化算法，通过奇异值分解法获得机翼位姿的初值，再采用Lagrange法结合拟Newton法迭代求解，得到满足机翼容差约束且综合装配误差最小的最优位姿。Xu等^[8]以测量不确定度为目标函数，采用粒子群优化算法迭代优化跟踪仪站位布局。

在使用点匹配算法时，由于被测对象在环境温度等变化下引起膨胀变形，测量值与理论值不能完全匹配，造成转站误差。俞慈君等^[9]建立了ERS点线性热变形模型，利用线性热变形补偿系数矩阵对ERS点热变形进行补偿，采用基于Levenberg-Marquardt (LM) 算法优化方法求解ERS理论值与测量值之间最小加权距离误差函数的优化模型。Pérez等^[10]研究了激光跟踪仪在热稳定以及空气扰动情况下的测量精度，热稳定误差是系统误差且可以被重现，因而可以采取相应的补偿措施。黄鹏等^[11]提出了一种加权点匹配计算方法来解决转站误差主要与公共点的测量误差的问题，以匹配残差加权平方和最小作为解的最优化条件，按照公共点的点位误差确定权重，并分析了测量误差和权重对匹配结果的影响。

在航空制造中广泛使用铝合金等热膨胀变形敏感材料，且工件的尺寸大，这导致工件容易受到外界温度变化的影响而发生热膨胀变形。工件等效热膨胀中心对跟踪仪转站参数具有重要影响。传统的奇异值分解(singular value decomposition, SVD)^[1]方法可以迅速地求得初步的转站参数。这种计算方法没有考虑产品工件不同尺度的三维热变形。Horn等^[2]转站计算方法考虑了产品工件在3个方向上有相同的热变形系数，默认工件的膨胀中心就是坐标原点。对于大型产品工件而言，膨胀中心不一定在工件的坐标原点，而膨胀中心会影响到转站参数的求解结果。

因此，本文建立考虑多变量的转站优化目标函数，并采用粒子群等现代智能优化算法进行最优求解。粒子群算法^[12](particle swarm optimization, PSO)是一种随机搜索方法。该优化算法受到群体动物捕食行动的启发，如鱼类和飞鸟等。种群中每个粒子的状态调整受3个因素^[13]的影响，包括随机因素、粒子自身运动惯性因素以及社会群体中最优个体因素。PSO算法首先在可行的搜索空间内随机初始化人工粒子并赋初值。粒子的维度代表设计变量的个数。建立目标函数，通过对每个粒子的适应度值进行计算，以搜索该粒子自身的最优化位置为止。每个粒子的运动由粒子自身历史最优值和种群全局最优位置引导，同时加入一定的随机影响因素，最后迭代更新每个粒子的位置。当粒子迭代寻优到更好的位置，这些新的优化结果将被选择来引导粒子群的运动。不停的重复迭代这个过程，直到最终找出满足约束条件的解。

1 采用SVD方法计算获得初步的转站参数

在飞机的装配过程中，整个测量系统包含跟踪仪测量坐标系、工装坐标系以及飞机坐标系等，如图 1所示。在靶标点的布置过程中，参照以下基本原则：1) 冗余设计原则：为了提高转站精度以及方便激光跟踪仪能测量足够的ERS点进行转站，ERS广泛分布在整个装配空间内。2) 空间包络原则：为了保证测量的精度，所有测量点必须位于ERS点所构成的包络空间之内。因为处于ERS点包络空间内的测量点，其测量误差必然小于所有ERS点中的最大测量误差。结合本项目自身的特点，即工件在4个支撑立柱中间，故靶标点主要分布在4个支撑立柱上，需要经常使用激光跟踪仪进行转站，即将ERS点的实际测量坐标值与其理论坐标值进行匹配，从而获得测量坐标系{B}与被测坐标系{A} 的转站参数(包含旋转矩阵以及平移矩阵)。

在不考虑工件热变形以及忽略工件膨胀中心点的条件下，设转站参数为(θ, ϕ, ϑ, P_x, P_y, P_z)，其中(θ, ϕ, ϑ)是测量坐标系{B}在被测坐标系{A}下的旋转参数，用X-Y-Z固定角表示，分别为回转角、俯仰角和偏转角；(P_x, P_y, P_z)是测量坐标系{B}原点在被测坐标系{A}下的位置值。那么测量坐标系{B}到被测坐标系{A}下的转站参数为

$ _B^A{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {\cos \theta \cos \phi }&{\cos \theta \sin \phi \sin \vartheta - \sin \theta \cos \vartheta }&{\cos \theta \sin \phi \cos \vartheta + \sin \theta \sin \vartheta }\\ {\sin \theta \cos \phi }&{\sin \theta \sin \phi \sin \vartheta + \cos \theta \cos \vartheta }&{\sin \theta \sin \phi \cos \vartheta - \cos \theta \sin \vartheta }\\ { - \sin \phi }&{\cos \phi \sin \vartheta }&{\cos \phi \cos \vartheta } \end{array}} \right], $

(1)

$ _B^A{\boldsymbol{T}} = {\left[ {\begin{array}{*{20}{l}} {{P_x}}&{{P_y}}&{{P_z}} \end{array}} \right]^{\rm{T}}}. $

(2)

其中，_B^AR是测量坐标系{B}到被测坐标系{A}的旋转矩阵，_B^AT是测量坐标系{B}到被测坐标系{A}的平移矩阵。

设在被测坐标系下布置了N个ERS标志点，第i个ERS点在测量坐标系{B}原点在被测坐标系{A}下的测量值和理论值分别为P_Ai和P_Bi。建立转站误差最小二乘函数：

$ {\sum ^2} = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left\| {{P_{Ai}} - \left( {_B^A{\boldsymbol{R}}{P_{Bi}} + _B^A{\boldsymbol{T}}} \right)} \right\|}^2}} . $

(3)

通过SVD^[1]分解方法，可以求得测量坐标系{B}与被测坐标系{A}的旋转矩阵_B^AR和平移矩阵_B^AT。

2 建立工件三维热膨胀变形的跟踪仪转站参数目标优化函数

转站参数包括采用X-Y-Z固定角表示的旋转矩阵参数(θ, ϕ, ϑ)，平移矩阵参数(P_x, P_y, P_z)以及在x、y、z 3个方向上的热膨胀变形系数(s_x, s_y, s_z)等。以转站后的坐标值与理论值之间误差的平方和最小为目标函数，建立考虑产品工件不同尺度因子的三维热膨胀变形转站的目标优化函数。

求解三维热膨胀变形转站参数算法包含9个求解变量(θ, ϕ, ϑ, P_x, P_y, P_z, s_x, s_y, s_z)。根据产品工件ERS点的理论值以及通过激光跟踪仪实际测量的对应ERS点位置结果，采用SVD方法获得初步的跟踪仪转站参数(θ, ϕ, ϑ, P_x, P_y, P_z)；工件等效热膨胀中心对跟踪仪转站参数具有重要影响。根据产品工件的三维数模，通过ANSYS有限元热分析计算，并结合实际激光跟踪仪测量的同一个ERS点在不同温度下坐标值变化，有限元仿真分析结果与跟踪仪实际测量结果对比分析工件初步热膨胀变形的中心点P₀=(P_x0, P_y0, P_z0)。

以SVD方法获得的初步转站参数以及有限元计算的初步膨胀中心位置为依据，适当扩大9个求解变量的取值范围。根据各个变量的转站参数变化范围，对种群个体进行随机初始化，并根据建立的优化目标函数计算每个种群个体的适应值。

采用PSO对种群个体进行迭代优化，通过对每个粒子的适应度值进行计算，以搜索该粒子自身的最优化位置为止。每个粒子的运动由粒子自身历史最优值和种群全局最优位置引导，同时加入一定的随机影响因素，最后迭代更新每个粒子的位置。当粒子迭代寻优到更好的位置，这些新的优化结果将被选择来引导粒子群的运动。对每个粒子个体进行判断，如果没有达到循环结束条件，则更新所有粒子的速度以及位置。如果达到循环结束条件，则返回种群的最优值，不停的重复迭代这个过程，直到最终找出满足约束条件的解。

2.1 建立优化的目标函数

复杂工件的热变形各向异性，通过引入三维热变形系数(s_x, s_y, s_z)，工件的等效膨胀中心为P₀=(P_x0, P_y0, P_z0)^T。某ERS点在温度为T₁时标定的理论位置为P₁=(P_x1, P_y1, P_z1)^T，在温度为T₂时位置值为P₂=(P_x2, P_y2, P_z2)^T，那么整理后可得

$ {{\boldsymbol{P}}_2} = {{\boldsymbol{P}}_0} + {\boldsymbol{S}} \cdot \Delta P = {{\boldsymbol{P}}_0} + {\boldsymbol{S}} \cdot \left( {{{\boldsymbol{P}}_1} - {{\boldsymbol{P}}_0}} \right). $

(4)

其中：S表示三维热变形比例系数矩阵，且S=diag(s_x, s_y, s_z); ΔP表示在标定温度下被测ERS点与等效膨胀中心原点的距离差值。

当从测量坐标系{B}中的某点P_Bi转换到被测坐标系{A}时，当工件温度发生改变，即工件三维热变形系数为(s_x, s_y, s_z)时，新的坐标位置为

$ {\boldsymbol{P}}{'_{Ai}} = {{\boldsymbol{P}}_0} + {\boldsymbol{S}} \cdot \left( {_B^A{\boldsymbol{R}} \cdot {P_{Bi}} + _B^A{\boldsymbol{T}} - {{\boldsymbol{P}}_0}} \right). $

(5)

采用最小二乘法计算的转站误差表达式为

$ \begin{array}{l} \;{\sum ^2} = \frac{1}{N}\mathop \sum \limits_{i = 1}^N \left\| {{P_{Ai}} - \left( {{{\boldsymbol{P}}_0} + } \right.} \right.\\ {\left. {\left. {{\boldsymbol{S}} \cdot \left( {_B^A{\boldsymbol{R}} \cdot {P_{Bi}} + _B^A{\boldsymbol{T}} - {{\boldsymbol{P}}_0}} \right)} \right)} \right\|^2}. \end{array} $

(6)

令${\mu _A} = \frac{1}{N}\sum\limits_{i = 1}^N {{P_{Ai}}} $，${\mu _B} = \frac{1}{N}\sum\limits_{i = 1}^N {{P_{Bi}}} $，${q_{Ai}} = {P_{Ai}} - {\mu _A}$，${q_{Bi}} = {P_{Bi}} - {\mu _B}$，$\Delta D = {\mu _A} - {\boldsymbol{S}} \cdot _B^A{\boldsymbol{R}} \cdot {\mu _B} - \left( {{{\boldsymbol{P}}_0} + S\left( {_B^A{\boldsymbol{T}}} \right.} \right.\left. {\left. { - {{\boldsymbol{P}}_0}} \right)} \right).$.

上式可以化简为

$ \begin{array}{l} {\sum ^2} = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left\| {{q_{Ai}} - {\boldsymbol{S}} \cdot _B^A{\boldsymbol{R}} \cdot {q_{Bi}} - \Delta D} \right\|}^2}} = \\ \;\;\;\;\;\;\;\;\;\frac{1}{N}\left( {\sum\limits_{i = 1}^N {{{\left\| {{q_{Ai}} - {\boldsymbol{S}} \cdot _B^A{\boldsymbol{R}} \cdot {q_{Bi}}} \right\|}^2}} - } \right.\\ \left. {2\Delta D\sum\limits_{i = 1}^N {{{\left( {{q_{Ai}} - {\boldsymbol{S}} \cdot _B^A{\boldsymbol{R}} \cdot {q_{Bi}}} \right)}^2}} + N{{\left\| {\Delta D} \right\|}^2}} \right) \end{array} $

(7)

由上式可知，转站参数的旋转矩阵_B^AR，三维热变形比例系数矩阵S都与膨胀中心P₀无关。而平移矩阵_B^AT与等效膨胀中心P₀有关。

当已知被测坐标系{A}的N个ERS点在标定温度为T₁时标定的理论位置P_Ai(i=1, 2, …, N)，以及在测量坐标系下对应的ERS点坐标值P_Bi(i=1, 2, …, N), 待求解的转站参数包含旋转矩阵参数(θ, ϕ, ϑ)、平移矩阵参数(P_x, P_y, P_z)、三维缩放矩阵参数(s_x, s_y, s_z)以及膨胀中心原点(P_x0, P_y0, P_z0)等。

2.2 采用粒子群优化算法迭代求解转站参数

目标函数式(7)包含9个变量，且结构复杂，难以采用传统的方法求解最优值。粒子群优化算法能够对多变量目标函数进行寻优，采用粒子群优化算法流程如图 2所示。以转站参数(θ, ϕ, ϑ, P_x, P_y, P_z, s_x, s_y, s_z)作为粒子群优化的变量，通过SVD分解方法获得初步的转站参数。设粒子种群中有m个粒子，每个粒子包含9个变量，第i个粒子的第j个变量在第k次运算时的速度与位置分别为v_ij^k、x_ij^k。根据粒子群算法，他们更新的规则为

$ v_{ij}^{k + 1} = wr_{ij}^k + {c_1}{r_1}\left( {p_{ij, {\rm{ best }}}^k - x_{ij}^k} \right) + {c_2}{r_2}\left( {g_{ij, {\rm{ best }}}^k - x_{ij}^k} \right), $

(8)

$ \begin{array}{l} x_{ij}^{k + 1} = x_{ij}^k + v_{ij}^{k + 1}, i = 1, 2, \cdots , m, \\ \;\;j = 1, 2, \cdots , 9, k = 1, 2, \cdots , n. \end{array} $

(9)

其中：c₁、c₂是权重系数，r₁和r₂是范围[0, 1] 间的随机均匀分布，p_{ij, best}^k、g_{ij, best}^k是在第k次迭代中个体历史最优值以及群体全局最优值。

w是权重系数, 取值范围一般在0.1~0.9。为了提高粒子群算法的收敛性能，将w随算法迭代的次数进行线性减小。设w_max为最大权重系数, w_min为最小权重系数, iter为当前迭代次数, iter_sum为算法迭代总次数, 则有:

$ w = {w_{\min }} + \frac{{{\rm{ iter}}{{\rm{ }}_{{\rm{sum }}}} - {\rm{ iter }}}}{{{\rm{ iter}}{{\rm{ }}_{{\rm{sum }}}}}}\left( {{w_{\max }} - {w_{\min }}} \right). $

(10)

在每个粒子位置更新过程中, 粒子每一维的最大速率限制设为v_max, 粒子每一维的坐标也被限制在允许范围之内。同时, 根据迭代计算结果，不断更新每个粒子的历史最优值以及种群全局最优值，即p_{ij, best}^k与g_{ij, best}^k, 最后输出的(θ, ϕ, ϑ, P_x, P_y, P_z, s_x, s_y, s_z, P_x0, P_y0, P_z0)就是粒子群优化算法得到的考虑三维热变形因素下转站参数的最优解，如图 3所示。

3 实验和数据分析

为了验证本文智能转站算法的有效性，以某部装配现场的转站数据为样本。图 4是实验室测量示意图(本文采用的测量数据是在工厂真实环境中的工件立柱上进行的测量。图 4是实验室进行原理验证时的照片，用于说明测量方法)。该部件的装配现场共有12个ERS点。采用跟踪仪测量某一温度下立柱上靶球座的坐标值。当环境温度变化且温度值恒定后，再次测量立柱上靶球座的坐标值。现场使用的Leica AT901跟踪仪，设备的标定精度为±15 μm+6 μm/m。部件的转站误差要求低于0.2 mm。12个ERS点的理论值以及对应的跟踪仪实测值如表 1所示。粒子群优化算法的最大迭代次数为400次，种群的粒子数为1 000个，学习因子c₁=2、c₂=2，惯性因子w_max=0.9、w_min=0.1。工件膨胀中心的等效原点采用有限元热膨胀分析以及采用激光跟踪仪测量ERS点在不同温度下的热膨胀变形，通过有限元仿真分析以及跟踪仪实测结果相互结合方法来初步确定工件等效膨胀中心，将此值作为粒子群个体膨胀中心的初始位置。

优化计算的转站参数_B^AR、_B^AT、S以及P₀的优化值为

$ _B^A{\boldsymbol{R}} = \left[ \begin{array}{l} \begin{array}{*{20}{c}} {\;\;0.171947}&{ - 0.985106}&{\;\;0.000098} \end{array}\\ \begin{array}{*{20}{c}} {\;\;0.000742}&{\;\;0.000229}&{\;\;0.999999} \end{array}\\ \begin{array}{*{20}{c}} { - 0.985106}&{ - 0.171947}&{\;\;0.000770} \end{array} \end{array} \right], $

$ _B^A{\boldsymbol{T}} = \left[ \begin{array}{l} \; - 79.405170\\ 12581.320114\\ \;3242.619615 \end{array} \right], {\boldsymbol{S}} = \left[ {\begin{array}{*{20}{l}} {0.999934}\\ {0.999932}\\ {0.999925} \end{array}} \right], $

$ {{\boldsymbol{P}}_0} = \left[ {\begin{array}{*{20}{c}} { - 10.9374}\\ {13873.9224}\\ { - 739.4043} \end{array}} \right]. $

图 5是全局最优目标函数值随迭代次数的变化，由图可知，目标函数随着迭代次数增加而迅速收敛，当迭代约230次时，全局最优目标函数值趋于稳定，最终的优化转站结果为0.026 778。而直接采用SVD转站的误差为0.190 375，采用Horn算法转站的平均误差为0.027 436。实验测试工件的温度变化最大差值约5 ℃左右，温度变化范围相对较小，热膨胀效应不是特别突出，这是导致本文优化结果虽然优于Horn算法，但差异不太明显的一个原因。但传统的SVD转站求解方法以及Horn算法等，适用于一些理想的变形情况。对于复杂情况的约束，采用本文的智能迭代求解方法，能够比较好地考虑这些额外的约束条件，为进一步的多约束条件求解提供了方法。

图 6是12个ERS点的理论值以及优化转站后的坐标值对应的误差向量图。为了能够比较直观的对比转站前后每个ERS点的误差，每个ERS点的误差放大了5 000倍。向量由每个ERS点的理论值指向改点转站后的点，向量的长度则是该ERS点误差的大小。如果某些点的转站误差太大，应该重点分析原因，如是否是测量误差等。通过误差向量的方向分布，可以判断工件变形的方向。如果误差都指向某个方向，则可能工件在这个方向上的变形比较明显，则应该采取相应的措施来抑制这类变形。图 6除了个别ERS点的误差较大以外，其余的ERS点的转站误差较小且比较均匀。图 7是12个ERS点的理论值以及优化转站后的坐标值对应的误差向量分别在XY、XZ、YZ 3个平面上的投影。

4 结论

根据工件的三维数模，通过有限元热分析计算工件热膨胀变形的中心点，建立考虑工件三维热膨胀变形的转站目标优化函数。采用粒子群优化算法对种群个体进行迭代优化，直到计算得到最优的跟踪仪转站参数。通过仿真分析与实验验证，得到了以下结论：

1) 工件的等效膨胀中心也是影响工件转站参数的重要因素。采用有限元仿真找出工件的等效膨胀中心，能够提高工件的转站参数精度。

2) 目标函数随着迭代次数增加而迅速收敛，当迭代约230次时，全局最优目标函数值趋于稳定，最终的优化转站结果为0.026 778。而直接采用SVD转站的误差为0.190 375，采用Horn算法转站的平均误差为0.027 436。

大型工件由于受温度变化而产生变形，从而影响跟踪仪的转站进度。复杂工件由于组成材料结构不同，因而在三维方向上的膨胀系数也略微不同。对于复杂情况的约束，采用本文的智能迭代求解方法，能够比较好地考虑这些额外的约束条件，为进一步的多约束条件求解提供了方法。