基于Kriging模型插值的孔位修正策略

石循磊，張繼文，劉順濤，陳懇,*

1. 清華大學(xué) 機械工程系，北京 100084 2. 航空工業(yè)成都飛機工業(yè)(集團)有限責(zé)任公司，成都 610091

飛機普遍采用鉚接或螺栓連接的形式實現(xiàn)零部件的裝配，疊層制孔和連接任務(wù)占據(jù)了大部分裝配任務(wù)[1]?，F(xiàn)代飛機一般具有復(fù)雜的氣動外形，且零件剛度較低，裝配時容易產(chǎn)生較大的變形和位置偏移[2-3]，因此難以像普通加工機床一樣，通過離線編程和在線對刀的形式一次性完成制孔作業(yè)任務(wù)。為了適應(yīng)零部件的變形和位置移動，保證制孔孔位的邊距準(zhǔn)確性，在機器制孔之前通過人工方式或特殊方式制備基準(zhǔn)孔(或稱預(yù)連接孔)，再通過激光跟蹤儀或者其他視覺測量系統(tǒng)檢測這些基準(zhǔn)孔的實際位置，將其與理論位置相比較，計算出基準(zhǔn)孔的孔位偏差[4-5]。依據(jù)已制的基準(zhǔn)孔位偏差，通過插值算法推算出待制孔的位置偏差，再對自動制孔程序進行校正，調(diào)整實際制孔位置，實現(xiàn)孔位補償?shù)墓δ堋?/p>

在自動孔位尋找校正過程中，通過已制基準(zhǔn)孔位推斷待制孔孔位偏差的補償算法顯得尤為重要。然而，飛機零件的弱剛性以及復(fù)雜外形特性為孔位推算帶來了多個難點：其一，裝配過程中，工件的變形并非簡單的剛性位置移動或旋轉(zhuǎn)，通常還伴隨有較大幅度的彈性扭轉(zhuǎn)和彎曲[6]，因此無法基于一系列基準(zhǔn)孔采用裝配體的剛體變換形式獲得待制孔孔位;其二，基準(zhǔn)孔通常位于不同的零件裝配件上，例如飛機的隔框縱梁長桁結(jié)構(gòu)交替出現(xiàn)，基準(zhǔn)孔也隨之交替出現(xiàn)，若基準(zhǔn)孔位與零件綁定，則為鑒別和記錄帶來繁雜的標(biāo)記過程，因而極易出現(xiàn)錯誤，進而導(dǎo)致整個孔位校正失效。

目前大多數(shù)工程化解決方案采用線性或三次樣條插值的方法校正兩個或多個基準(zhǔn)孔之間的孔位偏差。這是一種近似的解決方案，飛機的復(fù)雜曲面及變形特征很難服從這一假設(shè)，特別在一些曲率變化較大的位置，發(fā)生孔位校正錯誤的概率大大增加[7]。為了保證校正的局部特性，不得不增加基準(zhǔn)孔的數(shù)量，并嚴(yán)重依賴工藝實驗和人工經(jīng)驗選擇基準(zhǔn)孔位[8]。同時線性插值、三次樣條插值需要選定插值節(jié)點，然而，對于三維空間位置而言，很難區(qū)分樣值點的內(nèi)外關(guān)系，也難以確定插值節(jié)點，而且在外插值情況下容易出現(xiàn)較大誤差[9]。

在理論研究方面，Zhu等[10]通過基準(zhǔn)孔構(gòu)建了雙線性誤差平面，實現(xiàn)了對制孔區(qū)域的孔位進行線性插值補償。周叔陽等[11]將模擬退火算法應(yīng)用于預(yù)連接孔的排布問題，并抽象為數(shù)學(xué)模型，給出了預(yù)連接孔排布優(yōu)化算法及流程。董輝躍等[12]從飛機曲面構(gòu)造原理角度，提出了一種基于誤差Coons曲面函數(shù)的補償方法。王龍飛等[13]通過視覺輔助和手眼標(biāo)定，計算出基準(zhǔn)孔對應(yīng)的機器人驅(qū)動坐標(biāo)誤差變換矩陣，進而補償待制孔誤差。王青等[14]提出一種基于孔邊距偏差修正約束孔位和利用Shepard插值修正其余孔位的策略，實現(xiàn)了較好的孔位修正效果，但是缺少對孔位預(yù)測誤差的分析以及基準(zhǔn)孔的補制策略。本文提出一種采用Kriging模型插值算法的孔位偏差校正算法。該方法利用Kriging插值模型，依據(jù)基準(zhǔn)孔位偏差推算待制孔偏差和預(yù)測標(biāo)準(zhǔn)差。算法依賴孔位估計的均方誤差評估基準(zhǔn)孔的數(shù)量和分布質(zhì)量，進而可以在基準(zhǔn)孔實際選擇和布置時，依賴反復(fù)的迭代和計算，指導(dǎo)增添新的基準(zhǔn)孔，從而提升孔位偏差估計的精度。

1 孔位補償原理

根據(jù)基準(zhǔn)孔位偏差推算待制孔位置偏差的原理如圖1所示。

圖1 由基準(zhǔn)孔位偏差推算待制孔位置偏差的原理Fig.1 Principle of estimating position deviation of hole to be made from reference hole

已經(jīng)依據(jù)理論位置及工件裝配實際狀況制備出N個基準(zhǔn)孔?；鶞?zhǔn)孔實際位置與CAD模型中的理論位置存在一定偏差。

(1)

設(shè)CAD模型中的一個待制孔理論位置坐標(biāo)為(u*,v*,w*)T，則補償目標(biāo)是推算出輸出實際制孔位置與待制孔理論位置的偏差值：

(2)

上述過程本質(zhì)是分別建立了孔位偏差的3個空間坐標(biāo)分量Δu、Δv、Δw與孔位理論位置[u,v,w]T之間的響應(yīng)面模型。建立這種響應(yīng)面關(guān)系的一種方法是通過材料力學(xué)、彈性力學(xué)等復(fù)雜的理論推導(dǎo)出這一關(guān)系，但難度過大，而本文方法通過多次實驗的點以及一種特殊的函數(shù)關(guān)系模擬出這種關(guān)系，建立Kriging模型。

2 Kriging模型建立

Kriging模型屬于響應(yīng)面模型，其本質(zhì)是一種數(shù)據(jù)擬合的過程，將函數(shù)視作一種隨機過程，通過有限的數(shù)據(jù)擬合出這一隨機過程的概率分布規(guī)律，從而得到函數(shù)，計算出在空間坐標(biāo)變化時函數(shù)數(shù)值的變化規(guī)律，進而計算出未知點的數(shù)值。Kriging模型具有良好的非線性函數(shù)擬合特性，不僅提供了未知點的預(yù)測值，還提供了預(yù)測標(biāo)準(zhǔn)差[15]。

Kriging模型的基本表達式為

y(x)=fT(x)β+μ(x)

(3)

式中：f(x)為q維基函數(shù)向量；β為q維系數(shù)向量；μ(x)為隨機變量，每給出一個采樣點(或稱自變量)x，在f(x)、β和μ(x)已知的條件下，即可計算出其響應(yīng)函數(shù)值y(x)[16]。β和μ(x)未知時，當(dāng)具有一系列采樣點x(1),x(2),…,x(N)及對應(yīng)的y(1),y(2),…,y(N)，即可估計出其參數(shù)，進而求得其他位置x*的函數(shù)值y(x*)。

如果將μ(x)當(dāng)做滿足正態(tài)分布、與位置無關(guān)的隨機變量μ，則可以通過多次實驗消除。然而，在Kriging模型內(nèi)，μ(x)是一個依賴于x的隨機變量，且兩個樣值點x(i)、x(j)距離越近，其相關(guān)性越強。針對全局范圍內(nèi)的采樣點x，滿足數(shù)學(xué)期望E[μ(x)]=0，同時，對于任意兩個采樣點x(i)、x(j)，滿足協(xié)方差關(guān)系：

(4)

式中：σ2為采樣點對應(yīng)的隨機變量的方差；R(θ,p,x(i),x(j))為關(guān)聯(lián)函數(shù)，體現(xiàn)了x(i)、x(j)兩個采樣點間的相關(guān)系數(shù)：

(5)

式中：θh為x向量中h分量的相關(guān)強度；ph為x向量中h分量的相關(guān)形式[17]。與擬合回歸不同，Kriging模型主要通過μ(x)描述y(x)的響應(yīng)值，并通常將fT(x)設(shè)置為簡單的1T常數(shù)。該模型假設(shè)μ(x)和y(x)僅與狀態(tài)x有關(guān)，因而不存在多項式擬合中的擬合誤差。

在fT(x)=1T和p預(yù)先指定的條件下，當(dāng)給出多組采樣點時，通過最大似然估計法，估算參數(shù){β,σ2,θ}，從而得到模型估計結(jié)果：

(6)

式(6)顯示Kriging模型估計的結(jié)果是一系列關(guān)聯(lián)函數(shù)的線性組合所表達的連續(xù)可微函數(shù)式：

(7)

式中：ai為最大似然假設(shè)所估計的N個常數(shù)[19]。

依據(jù)R(θ,x*,x(i))的定義(設(shè)ph=2時)：

(8)

3 孔位偏差估計及改進策略

針對孔位偏差問題，距離待制孔距離越大的基準(zhǔn)孔位偏差對當(dāng)前待制孔偏差的影響越小，而距離待制孔距離越小的基準(zhǔn)孔孔位偏差則對當(dāng)前待制孔的影響較大，這種特性能夠有效描述待制孔與基準(zhǔn)孔之間的偏差關(guān)聯(lián)關(guān)系。如圖2所示。

根據(jù)Kriging模型的定義及參數(shù)估計方法，當(dāng)x*=x(i)時，即取已知實驗點時，由響應(yīng)面函數(shù)計算可得：

(9)

式(9)表明由式(6)計算出的基準(zhǔn)孔x(i)處的偏差值仍然是實測值y(i)，且其估計偏差為0。因此，Kriging模型所建立的是對已知實驗點的一種特殊的插值，而相比于線性插值或樣條插值，孔位偏差不再限于內(nèi)插值，而是可以實現(xiàn)一定程度的外插值。同時對于一個待制孔的誤差，其不僅取決于兩端的基準(zhǔn)孔，還取決于與其臨近的若干個孔的偏差，這與一般裝配的誤差特性一致。

圖2 待制孔與基準(zhǔn)孔之間的偏差關(guān)聯(lián)關(guān)系Fig.2 Deviation correlation between hole to be made and reference hole

(10)

(11)

因此，采用Kriging模型插值推算孔位偏差，能夠依賴孔位估計的均方誤差評估基準(zhǔn)孔的數(shù)量和分布質(zhì)量，進而可以在基準(zhǔn)孔實際選擇和布置時，依賴反復(fù)的迭代和計算，指導(dǎo)增添新的基準(zhǔn)孔，從而提升孔位偏差估計的精度。此外，可以通過式(7)采用一種連續(xù)可微的函數(shù)形式模擬孔位偏差，從而繪制出孔位誤差分布圖，有利于裝配質(zhì)量以及基準(zhǔn)孔制備的迭代改進。算法的總體結(jié)構(gòu)如圖3所示。

圖3 算法總體結(jié)構(gòu)Fig.3 Overall structure of algorithm

4 實驗與分析

飛機孔位偏差來源可分為待制孔零件的變形和平移旋轉(zhuǎn)兩種，實驗中分別驗證兩種偏差來源形式下本文算法的有效性。

針對零件變形帶來的孔位偏差，采用有限元仿真實驗來驗證算法的有效性。利用有限元軟件Abaqus建模，如圖4所示，將蒙皮簡化為邊長約4.8 m、厚度3 mm的梯形曲面工件，連接于骨架上，并施加載荷。坐標(biāo)系y軸定義為與曲面蒙皮3個相距較遠頂點組成的平面垂直的方向，即y軸近似與蒙皮表面垂直，而x和z軸在該平面內(nèi)。

圖4 有限元仿真載荷示意圖Fig.4 Load diagram of finite element simulation

有限元仿真結(jié)果如圖5所示，圖中標(biāo)出了25個基準(zhǔn)孔，在蒙皮上另取50個點作為待制孔的位置。將有限元仿真結(jié)果中基準(zhǔn)孔的位置作為實際孔位，可得基準(zhǔn)孔的位置偏差(位移)Δu、Δv、Δw作為Kriging模型的因變量，以基準(zhǔn)孔理論位置為自變量，可建立3個Kriging模型。利用建立好的模型估計待制孔的應(yīng)制孔位置，結(jié)果如圖6所示。

圖6中，eu、ev、ew是Kriging模型估計的孔位偏差值(位移)的誤差，定義為

(12)

從圖6中可以看出，變形導(dǎo)致的x和z方向位移量較小，孔位估計誤差均小于0.1 mm。y方向位移達到10 mm量級，可認(rèn)為是實際裝配時壁板的極限變形，因此主要分析y方向?qū)嶒灲Y(jié)果。

針對零件整體平移及旋轉(zhuǎn)帶來的孔位偏差，偏移前后孔位坐標(biāo)有Q=R0Q0+T0

(13)

式中：Q0為孔位理論坐標(biāo)；Q為偏移后坐標(biāo)；R0為旋轉(zhuǎn)矩陣；T0為平移變換矩陣，因此孔位偏差量為Q-Q0。

在實際裝配過程中，蒙皮壁板放置在飛機骨架上，受骨架約束，可認(rèn)為壁板僅繞圖8中y軸旋轉(zhuǎn)。驗證實驗中取旋轉(zhuǎn)角度為1.15°，可得旋轉(zhuǎn)矩陣：

圖5 有限元仿真位移結(jié)果云圖Fig.5 Displacement nephogram of finite element simulation

(14)

同樣采用圖5中的25個基準(zhǔn)孔建立Kriging模型，50個待制孔的x方向和z方向最大位移均達10 mm量級，而孔位估計誤差如圖9所示，均小于0.05 mm。因此對于零件整體旋轉(zhuǎn)平移導(dǎo)致的孔位偏差，該算法能很好地估計偏移后的位置。

利用上述方法在某型飛機數(shù)字化裝配現(xiàn)場實驗的結(jié)果見表1，實物實驗中基準(zhǔn)孔的偏移量均小于1 mm,該方法估計出的待制孔的偏差量與測量值相差很小。

考慮該算法的計算效率，實驗使用的計算機CPU為i7-8700k，利用R語言編程，由25個基準(zhǔn)孔建立Kriging模型并估計50個待制孔位所需運算時間為0.214 5 s，50個基準(zhǔn)孔估計100個待制孔耗時0.306 2 s，90個基準(zhǔn)孔估計150個待制孔耗時0.480 7 s，而實際生產(chǎn)中孔位估計并無實時性要求，所以計算效率完全能達到要求。與線性插值等傳統(tǒng)方法相比較，該算法除了能給出待制孔的預(yù)測標(biāo)準(zhǔn)誤差外，另一個顯著優(yōu)勢體現(xiàn)在曲面區(qū)域外插值的情況。例如在上述有限元仿真實驗中，一個待制孔的坐標(biāo)為(-1 592.7, 75.251, -1 008.3)，該待制孔附近的兩個基準(zhǔn)孔坐標(biāo)為(-1 692.1, 66.179, -1 334.5)和(-1 690.7, 66.312, -814.05)，以這兩個點為插值節(jié)點，線性插值得待制孔的y方向偏差值為154.78 mm，而仿真結(jié)果為2.990 3 mm，本文算法估計結(jié)果為2.847 1 mm，此時線性插值出現(xiàn)很大誤差，而該算法仍有較好效果。

圖6 變形導(dǎo)致孔位偏差的實驗結(jié)果Fig.6 Experimental results of hole displacement caused by deformation

圖7 增添基準(zhǔn)孔后實驗結(jié)果Fig.7 Experimental results after adding reference holes

圖8 實驗件坐標(biāo)系定義Fig.8 Definition of experimental coordinate system

多組實驗表明本文方法能較準(zhǔn)確、高效地通過基準(zhǔn)孔信息估計待制孔的位移，并且給出有實際參考意義的待制孔位置的預(yù)測標(biāo)準(zhǔn)誤差，用以指導(dǎo)排布或增添基準(zhǔn)孔，最后通過在合適位置增添基準(zhǔn)孔能有效提高待制孔位移估計的準(zhǔn)確度。

圖9 平移旋轉(zhuǎn)導(dǎo)致孔位偏差的實驗結(jié)果Fig.9 Experimental results of hole displacement caused by translation and rotation

表1 實物實驗待制孔位估計結(jié)果

5 結(jié) 論

1) 提出了一種基于基準(zhǔn)孔的孔位校正方法，該方法利用Kriging模型，依據(jù)基準(zhǔn)孔位偏差推算出待制孔偏差和預(yù)測標(biāo)準(zhǔn)誤差。

2)相比于線性或三次樣條插值，該方法適用范圍更大，且能給出孔位估計的標(biāo)準(zhǔn)誤差，以評估基準(zhǔn)孔的數(shù)量和分布質(zhì)量，進而指導(dǎo)增添新的基準(zhǔn)孔。

3)利用有限元仿真實驗，驗證了在零件變形及平移旋轉(zhuǎn)導(dǎo)致孔位偏差情況下本文方法的有效性，在合理增添基準(zhǔn)孔后，孔位估計的誤差小于0.3 mm，滿足航空裝配要求。