基于PSO與WSVD的飛機部件位姿擬合方法

何曉煦，雷沛，潘登，楊陽，李現(xiàn)坤，鄧珍波

1.航空工業(yè)成都飛機工業(yè)（集團）有限責(zé)任公司，成都 610092

2.四川省航空智能制造裝備工程技術(shù)研究中心，成都 610092

現(xiàn)代飛機制造已廣泛采用數(shù)字化調(diào)姿對合技術(shù)［1］以實現(xiàn)飛機部件姿態(tài)的自動測量和調(diào)整［2］，解決了傳統(tǒng)人工作業(yè)方式存在的調(diào)姿精度差、效率低下、勞動強度大等問題［3］。飛機部件數(shù)字化調(diào)姿包含3個主要步驟：位姿測量、位姿擬合、位姿調(diào)整。位姿測量［4］是利用激光跟蹤儀等大尺寸測量設(shè)備測量飛機部件上的調(diào)姿基準(zhǔn)點在現(xiàn)場裝配坐標(biāo)系下的坐標(biāo)［5］；位姿擬合是將測量數(shù)據(jù)和理論數(shù)據(jù)進行配準(zhǔn)，擬合出部件在全機坐標(biāo)系下的位姿參數(shù)；位姿調(diào)整是根據(jù)調(diào)姿機構(gòu)運動學(xué)關(guān)系，將部件的位姿參數(shù)轉(zhuǎn)化為數(shù)控定位器的運動量，驅(qū)動飛機部件達到目標(biāo)位姿［6］。

部件位姿擬合是數(shù)字化調(diào)姿的核心環(huán)節(jié)之一，所用的算法直接影響調(diào)姿精度和效率［7］。當(dāng)前較常用的方法可以分為非迭代算法和迭代算法兩類。常用的非迭代算法包括奇異值分解法（Singular Value Decomposition， SVD）［8-9］、三點法、四元數(shù)法等。迭代算法包括最小二乘法（Least Squares， LS）［10-11］、隨機抽樣一致性算法（Random Sample Consensus， RANSAC）、迭代最近點法（Iterative Closest Point， ICP）［12］及其各類改進算法［13］等。

由于飛機部件裝配環(huán)節(jié)存在各類累積誤差，實際產(chǎn)品與理論數(shù)模不可避免地存在差異。為了控制部件的裝配精度及對接后的整體外形精度，不同區(qū)域的基準(zhǔn)點通常會設(shè)計不同的公差要求。傳統(tǒng)的SVD算法雖然控制了所有基準(zhǔn)點的綜合轉(zhuǎn)換殘差和，但經(jīng)常導(dǎo)致精度要求高的基準(zhǔn)點超差而精度要求低的點還有較大的調(diào)整余量。Chen等［14］提出了一種加權(quán)位姿擬合算法，利用有限元分析不同區(qū)域基準(zhǔn)點理論值的偏移量，偏移越大的點其權(quán)重設(shè)置越小，但權(quán)重的計算公式是固定的，適用范圍會受到影響。朱緒勝和鄭聯(lián)語［15］提出了一種基于關(guān)鍵裝配特性的部件位姿多目標(biāo)優(yōu)化算法，并結(jié)合粒子群優(yōu)化（Particle Swarm Optimization， PSO）算法［16］對位姿參數(shù)進行求解，將位姿參數(shù)作為六維粒子進行迭代計算，但未考慮權(quán)重因素影響。陳遠志等［17］提出了一種改進粒子群算法，引入懲罰函數(shù)處理基準(zhǔn)點的容差約束，減少擬合殘差，但未考慮權(quán)重影響。蔣睿嵩等［18］研究了一種基于權(quán)值約束的點云精確配準(zhǔn)模型，通過對配準(zhǔn)模型中較為重要的區(qū)域賦予較高的權(quán)重來提高該區(qū)域的配準(zhǔn)精度，但未考慮其他約束條件。

針對以上算法的不足，結(jié)合飛機大部件調(diào)姿的工藝特點，提出一種基于粒子群優(yōu)化與加權(quán)奇異值分解（Weighted Singular Value Decomposition， WSVD）的位姿擬合算法，根據(jù)基準(zhǔn)點的精度要求賦予相應(yīng)的權(quán)重，將所有基準(zhǔn)點的轉(zhuǎn)化殘差均滿足精度要求作為約束條件，將點的權(quán)重作為粒子進行迭代優(yōu)化，求解大部件的位姿參數(shù)，并在不同型號的飛機大部件調(diào)姿對合系統(tǒng)中進行應(yīng)用驗證。

1 加權(quán)奇異值分解算法

1.1 飛機部件位姿描述

在飛機坐標(biāo)系下，飛機部件的基準(zhǔn)點理論位置坐標(biāo)為Qi(Qi，x，Qi，y，Qi，z)，i∈(1，2，…，n)，n為基準(zhǔn)點個數(shù)，由激光跟蹤儀測量的實測位置坐標(biāo)為Pi(Pi，x，Pi，y，Pi，z)。假設(shè)部件經(jīng)過位姿調(diào)整量M變換后到達理論位置，則有

式中：Tx、Ty、Tz分別為x、y、z方向的空間平移量；α、β、γ分別為繞x、y、z方向的空間角度；R為3×3的空間旋轉(zhuǎn)矩陣；T為3×1的空間平移向量。

飛機部件位姿調(diào)整量M可基于誤差平方和最小化原則進行擬合，其中，LS算法求解的旋轉(zhuǎn)矩陣由于沒有正交性的約束，可能造成過分追求數(shù)據(jù)逼近而引入較大誤差。相比之下，SVD算法在求解過程中利用了旋轉(zhuǎn)矩陣的正交性而更符合實際情況［19］。

由于兩組點具有相同的質(zhì)心，可通過質(zhì)心化簡化計算過程，即先將平移向量T從式（2）中分離，優(yōu)先求解旋轉(zhuǎn)矩陣R。分別求出理論位置和實測位置的質(zhì)心坐標(biāo)分別為

式中：n≥3。質(zhì)心化后的坐標(biāo)為

根據(jù)式（1）～式（6），有

1.2 奇異值分解法

1）SVD優(yōu)化配準(zhǔn)算法

采用SVD進行優(yōu)化配準(zhǔn)計算時，質(zhì)心化后誤差平方和最小的目標(biāo)函數(shù)為

式中：sR為SVD求得的旋轉(zhuǎn)矩陣。則有

令

則上述問題等價于求max(trace(sRsH))，trace(·)為矩陣的跡。將sH進行奇異值分解得

式中：sU和sV為正交矩陣；sΛ為對角矩陣。則由SVD解算的sR的最優(yōu)值為

基準(zhǔn)點的平移向量為

由SVD解算的位姿調(diào)整量中的平移向量為

2）SVD的位姿擬合誤差

部件經(jīng)過位姿調(diào)整后，基準(zhǔn)點坐標(biāo)值為

調(diào)姿后基準(zhǔn)點的殘差為

sεi應(yīng)滿足飛機大部件基準(zhǔn)點容差要求εi，d，d∈(x，y，z)。

SVD的目標(biāo)函數(shù)是基準(zhǔn)點擬合誤差的平方和（見式（9）），但飛機部件裝配過程中，對不同基準(zhǔn)點的容差要求不同，這導(dǎo)致基于SVD的姿態(tài)擬合算法雖然能夠解算得到總殘差最小的部件位姿，但卻存在部分基準(zhǔn)點超出容差范圍的問題，而且往往是那些精度要求高、容差范圍小的關(guān)鍵基準(zhǔn)點。

1.3 加權(quán)的奇異值分解法

1）WSVD優(yōu)化配準(zhǔn)算法

針對SVD在飛機部件調(diào)姿方面的缺陷，提出一種能夠自適應(yīng)基準(zhǔn)點容差要求的姿態(tài)擬合算法。在位姿調(diào)整量的解算過程中引入權(quán)重，利用權(quán)值控制關(guān)鍵基準(zhǔn)點誤差，可以提高位姿調(diào)整量解算的準(zhǔn)確度，從而提高飛機部件調(diào)姿的合格率。

假設(shè)第i個基準(zhǔn)點x、y、z坐標(biāo)對應(yīng)的權(quán)重值相等，均為wi。將權(quán)重代入誤差平方和最小的目標(biāo)函數(shù)中，此時目標(biāo)函數(shù)表示為

式中：wR為WSVD求得的旋轉(zhuǎn)矩陣。與SVD求解方法類似，將式（18）展開得

將求解最小化問題minwΕ轉(zhuǎn)換為求解

則上述問題等價于求max(trace(wRwH))。對wH進行奇異值分解得

式中：wU和wV為正交矩陣；wΛ為對角矩陣。則wR的最優(yōu)解為

先求出基準(zhǔn)點加權(quán)后的平移向量為

則WSVD位姿調(diào)整量中的平移向量為

2）WSVD的位姿擬合誤差

飛機部件調(diào)姿后，基準(zhǔn)點的三維坐標(biāo)為

用wεi表示調(diào)姿后的殘差為

wεi應(yīng)滿足飛機大部件基準(zhǔn)點容差要求εi，d。

引入權(quán)值進行加權(quán)的奇異值分解，能較好地分配飛機部件基準(zhǔn)點的殘差在最小化問題的權(quán)值，使結(jié)果趨近于權(quán)重值更大的基準(zhǔn)點，即偏向裝配精度要求更高、容差范圍更小的基準(zhǔn)點，從而避免這些基準(zhǔn)點超出容差范圍的問題，使求解出的部件位姿調(diào)整量更貼近實際調(diào)姿需求。

2 基于粒子群優(yōu)化的權(quán)重分配算法

2.1 粒子群優(yōu)化算法

1）PSO目標(biāo)函數(shù)

PSO具有較強的魯棒性和全局尋優(yōu)能力，能保持個體間持續(xù)的信息交互，并且易與其他算法進行結(jié)合［20］。利用PSO自動搜尋能使擬合后所有基準(zhǔn)點滿足容差要求的權(quán)重值，優(yōu)化目標(biāo)是基準(zhǔn)點在x、y、z方向的超差項的殘差盡可能小。因此定義目標(biāo)函數(shù)適應(yīng)度值E為擬合誤差超差的基準(zhǔn)點，其殘差的絕對值之和，即

式中：為擬合誤差超差的基準(zhǔn)點在d方向的殘差的絕對值。

利用PSO迭代求解基準(zhǔn)點權(quán)值時，將待求解的權(quán)值作為粒子。粒子群每次更新后，權(quán)值的改變將引起WSVD擬合結(jié)果中殘差的變化，從而使目標(biāo)函數(shù)E中的超差項隨粒子群的更新而動態(tài)調(diào)整，引起PSO適應(yīng)度值更新，進而指導(dǎo)粒子群的搜索方向，因此具有自適應(yīng)的尋優(yōu)能力。

2）多維粒子群算法

每個粒子包含一個位置向量和一個速度向量?；鶞?zhǔn)點對應(yīng)的權(quán)值向量為w，則粒子的維度為n。假設(shè)選取m個粒子，則每次迭代時粒子群存在n×m維的位置矩陣x，以及n×m維的速度矩陣v。

粒子群初始化時，將位置矩陣x初始化為每項均為1的單位矩陣。則初始位置因子=1，速度矩陣v初始化為每項為[0， 1]的隨機數(shù)矩陣。

粒子搜索解空間時，每個粒子根據(jù)自身的適應(yīng)度值保存各自搜索到的個體最優(yōu)經(jīng)歷位置pp，并在個體最優(yōu)經(jīng)歷位置中選取群體最優(yōu)經(jīng)歷位置pg。粒子群按式（28）調(diào)整位置矩陣：

式中：c1和c2為加速因子，取c1=c2=2；ξ和η為[0， 1]中均勻分布的隨機數(shù)；λ為慣性權(quán)重因子，采用線性遞減策略對其進行賦值，有

每個粒子的權(quán)值總和

式中：λmax和λmin分別為λ值的上下限，取λmax=0.95，λmin=0.4；a為當(dāng)前迭代次數(shù)；A為總迭代次數(shù)，A=50。

粒子群的位置矩陣更新為

式中：δ為位置更新時的約束因子，取δ=1。

為提高算法的尋優(yōu)能力和收斂速度，在粒子每次更新后，位置矩陣x應(yīng)滿足非負性條件。否則可能導(dǎo)致粒子朝相反方向搜索造成非全局最優(yōu)解，或搜索時間被無限放大。當(dāng)位置因子xi，k不滿足非負性條件時，將更新后的第k個粒子的第i個基準(zhǔn)點的權(quán)值xi，k重新初始化為1，即

同時，在PSO迭代過程中，位置矩陣x應(yīng)滿足歸一性條件。通過歸一化處理來保證每個粒子的權(quán)值總和為n不變，有：

此，每個基準(zhǔn)點權(quán)值的空間范圍為[0，n]。

2.2 本文方法計算流程

本文方法的計算流程如圖1所示，首先用基準(zhǔn)點進行SVD擬合，當(dāng)不滿足算法終止條件時，進入PSO部分，粒子根據(jù)目標(biāo)函數(shù)適應(yīng)度值，在約束空間內(nèi)按相應(yīng)的規(guī)律搜索，不斷調(diào)整自身的位置矩陣和速度矩陣，每次得到更新的權(quán)值后，再帶入WSVD部分進行加權(quán)奇異值分解，從而更新PSO目標(biāo)函數(shù)適應(yīng)度值，繼而指導(dǎo)粒子群的搜索方向，最后在滿足算法終止條件的情況下停止計算。具體步驟為

圖1 本文方法計算流程Fig.1 Calculation process of this method

步驟1將理論位置Qi和實測位置Pi用SVD擬合，求得位姿調(diào)整后的殘差sεi（見式（3）～式（17）），判斷是否滿足容差要求。滿足進入步驟10，反之進入步驟2。

步驟2初始化粒子群。利用步驟1中的殘差計算PSO適應(yīng)度值E（見式（27）），得到個體最優(yōu)位置pp和群體最優(yōu)位置pg。進入步驟3。

步驟3更新速度矩陣v和位置矩陣x（見式（28）～式（32））。進入步驟4。

步驟4將更新后的位置矩陣x的每個粒子分別作為權(quán)值w代入WSVD擬合，求得位姿調(diào)整后的殘差wεi（見式（18）～式（26））。進入步驟5。

步驟5根據(jù)步驟4的殘差wεi求得每個粒子的適應(yīng)度值E（見式（27）），更新粒子個體最優(yōu)位置pp和群體最優(yōu)位置pg。進入步驟6。

步驟6將更新后的群體最優(yōu)位置pg作為權(quán)值，判斷對應(yīng)的WSVD分解后的殘差wεi是否滿足容差要求。滿足進入步驟9，反之進入步驟7。

步驟7判斷當(dāng)前是否達到最大迭代次數(shù)A。滿足進入步驟8，反之進入步驟3。

步驟8輸出步驟6中的群體最優(yōu)位置pg作為權(quán)值w，以及對應(yīng)的WSVD的位姿矩陣，該結(jié)果為當(dāng)前迭代條件下，最接近容差要求的解。如需要繼續(xù)尋找滿足容差要求的解，可以通過增加粒子數(shù)或迭代次數(shù)，或者調(diào)整迭代參數(shù)后，從步驟1重新開始計算。

步驟9輸出步驟6中的群體最優(yōu)位置pg作為權(quán)值w，以及對應(yīng)的WSVD的位姿矩陣。找到了滿足容差要求的解，計算結(jié)束。

步驟10輸出SVD的位姿矩陣，權(quán)值w每項均為1。找到了滿足容差要求的解，計算結(jié)束。

可見，本文方法在SVD無法滿足容差的情況下才啟用PSO+WSVD自動尋優(yōu)，因此既具有常規(guī)SVD數(shù)據(jù)處理量小的優(yōu)點，又具備PSO+WSVD全局搜索和自適應(yīng)權(quán)重分配的能力。

3 案例驗證與分析

3.1 C919機頭上下部調(diào)姿對合

為了驗證本文方法的可行性，以C919機頭上部和下部的調(diào)姿對合過程為算例進行分析，以機頭下部作為基準(zhǔn)部件，機頭上部作為調(diào)姿部件與機頭下部對接。分別使用SVD方法和本文方法計算出機頭上部的位姿調(diào)整量為M1和M2，若機頭上部經(jīng)過M1調(diào)姿后的殘差超差，而經(jīng)過M2調(diào)姿后的殘差在容差范圍內(nèi)，則說明本文方法具有可行性。C919機頭上部有7個基準(zhǔn)點，理論值和實測值如表1所示，根據(jù)實際裝配中的工藝要求，基準(zhǔn)點對應(yīng)的三維容差約束為± 1 mm。

表1 基準(zhǔn)點位置信息Table 1 Location information of datum points

用1.2節(jié)的SVD方法進行計算。計算時間為0.005 s，位姿調(diào)整量為M1=［4.056 mm，5.279 mm， 8.286 mm， 0.011°， 0.006°， -0.013°］。機頭上部經(jīng)過M1調(diào)姿后，基準(zhǔn)點的殘差如表2所示。調(diào)姿后仍有1個超差項：2號基準(zhǔn)點的y方向，殘差為-1.235 mm（超差-0.235 mm）。因此，SVD方法無法使所有基準(zhǔn)點的殘差都滿足容差要求。

表2 基準(zhǔn)點殘差比較Table 2 Datum points error comparison

在實際工程應(yīng)用環(huán)境下，在使用M1自動調(diào)姿后，仍需要人工進行干預(yù)，工人重新測量基準(zhǔn)點的實際位置，憑經(jīng)驗進行手工微調(diào)。據(jù)現(xiàn)場統(tǒng)計，為使部件調(diào)整到符合對接精度要求的姿態(tài)，需要人工手動微調(diào)1～5次，花費時間長達數(shù)十分鐘到數(shù)小時不等。人工干預(yù)環(huán)節(jié)對工人技能和經(jīng)驗要求高，較大程度地制約了部件調(diào)姿對合的效率。并且，仍存在無法調(diào)整合格的風(fēng)險。

用2.2節(jié)的本文方法進行計算。PSO的粒子數(shù)為30，機頭上部有7個基準(zhǔn)點，對應(yīng)粒子維度為7，因此每個基準(zhǔn)點權(quán)值的空間范圍為[0， 7]。PSO的初始適應(yīng)度值為1.235 mm，經(jīng)過第1次迭代，適應(yīng)度值降為1.043 mm，經(jīng)過第2次迭代降為0，計算時間為0.144 s?；鶞?zhǔn)點權(quán)值為w1=［0.779， 1.646， 0.858， 0.746， 0.627，1.265， 1.079］。位姿調(diào)整量為M2=［4.117 mm，5.686 mm， 8.334 mm， 0.009°， 0.007°， -0.015°］。

機頭上部使用M2調(diào)姿后，基準(zhǔn)點的殘差如表2所示，殘差均滿足容差要求。由表2可以看出，本文方法通過權(quán)值的自動尋優(yōu)，使2號基準(zhǔn)點y方向的對接誤差減小為-0.962 mm，對比SVD方法提高了0.273 mm。本文方法不僅使SVD方法的超差項滿足容差要求，并且保證其他基準(zhǔn)點殘差均在容差范圍內(nèi)，避免了人工微調(diào)環(huán)節(jié)帶來的效率損失和無法調(diào)整到位的風(fēng)險，實現(xiàn)了完全的自動調(diào)姿。

在SVD方法和本文方法計算時間都很短的情況下，本文方法能夠解決SVD方法局部容差超差的情況，證明了本文方法的可行性。

3.2 實例分析

為說明本文方法在其他機型和部件調(diào)姿中的有效性，采用文獻［21］的實測數(shù)據(jù)作為算例進行分析。只考慮文獻［21］的后機身調(diào)姿后，基準(zhǔn)點的殘差是否滿足容差要求，如果滿足則認(rèn)為后機身能夠?qū)崿F(xiàn)與前機身對接，本文方法有效。根據(jù)文獻［21］，后機身8個基準(zhǔn)點的理論值和實測值如表3所示，基準(zhǔn)點對應(yīng)的三維容差約束為±1 mm。

表3 文獻［21］基準(zhǔn)點位置信息Table 3 Location information of datum points for Ref.［21］

按照± 1 mm的公差要求分別采用SVD方法和本文方法進行計算，兩種方法的殘差如表4所示。從表4可以看出，兩種方法計算結(jié)果一致，且基準(zhǔn)點殘差均在容差范圍內(nèi)。由于SVD方法就能使后機身調(diào)姿后的基準(zhǔn)點殘差滿足± 1 mm的容差要求，因此本文方法在通過SVD求出位姿調(diào)整量后就終止了計算，沒有啟用PSO+WSVD。此時位姿調(diào)整量為M3= ［-246.554 mm，212.303 mm， 8.260 mm， -0.569°， 2.496°，-1.249° ］，最大的殘差值為-0.871 mm，為4號點的x方向。

表4 文獻［21］基準(zhǔn)點殘差比較Table 4 Datum points error comparison for Ref.［21］

假設(shè)公差要求提高到± 0.65 mm，分別采用SVD方法和本文方法進行計算。

使用SVD方法進行計算。計算時間為0.002 s，求得位姿調(diào)整量仍為M3，基準(zhǔn)點殘差如表4所示。由表4可以看出，SVD方法在容差要求為± 1 mm和± 0.65 mm兩種情況下計算結(jié)果一致，這是由于SVD是以控制總殘差和最小為目標(biāo)進行計算的，容差的改變并不影響計算結(jié)果。然而，這樣的計算結(jié)果，在容差要求提高后，就出現(xiàn)了超差項：4號點的x方向，殘差為-0.871 mm（超差-0.221 mm）。此時，SVD方法無法使所有基準(zhǔn)點滿足± 0.65 mm的容差要求。

使用本文方法進行計算。PSO的粒子數(shù)為30，后機身8個基準(zhǔn)點對應(yīng)粒子維度為8，因此每個基準(zhǔn)點權(quán)值的空間范圍為［0， 8］。PSO適應(yīng)度曲線收斂圖如圖2所示，初始適應(yīng)度值為0.871 mm，經(jīng)過第1次降迭代后，適應(yīng)度值降為0.790 mm，經(jīng)過第2次迭代后降為0.736 mm，經(jīng)過第3次迭代后降為0，計算時間為0.164 s?；鶞?zhǔn)點權(quán)值為w2=［0.201， 0.600， 1.300，2.091， 0.353， 1.546， 0.743， 1.165］。位姿調(diào)整量為M4=［-246.439 mm， 210.657 mm，8.936 mm， -0.561°， 2.499°， -1.245° ］，基準(zhǔn)點的殘差如表4所示。

圖2 PSO適應(yīng)度值曲線收斂圖Fig.2 PSO fitness curve convergence diagram

由表4可以看出，本文方法通過權(quán)值的自動尋優(yōu)，使4號點x方向的對接誤差減小為-0.644 mm，對比SVD方法提高了0.227 mm。本文方法不僅使SVD方法的超差項滿足容差要求，并且保證其他基準(zhǔn)點殘差均在滿足± 0.65 mm的容差要求。據(jù)此，驗證了本文方法在其他機型和部件調(diào)姿中應(yīng)用的有效性。

4 結(jié) 論

針對飛機大部件調(diào)姿中的位姿擬合問題，提出了一種基于粒子群優(yōu)化結(jié)合加權(quán)奇異值分解的方法，將不同精度要求的基準(zhǔn)點賦予不同的權(quán)重，權(quán)值的分配由粒子群算法進行迭代優(yōu)化，并將每個基準(zhǔn)點的轉(zhuǎn)換殘差均滿足精度要求作為優(yōu)化的約束條件。通過飛機數(shù)字化裝配中的應(yīng)用分析驗證了算法的有效性。