基于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 結 論

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