何帥,陳富民,李建華,張廳方

(1.西安交通大學新能源與質(zhì)量工程研究所,710049,西安; 2.東方汽輪機有限公司,618000,四川德陽)

目前對于精密零件的測量,一般采用在三坐標機上測得坐標,再通過相關(guān)算法計算所需要的零件質(zhì)量特性,這種方法屬于坐標法。采用坐標法測量零件坐標時,由于無法直接利用被測零件本身基準為測量基準,測量數(shù)據(jù)里包含測量坐標系與設(shè)計坐標系不重合造成的位置誤差。葉片屬于精密零件,采用坐標法對葉片型線進行質(zhì)量評定,精度受基準不重合影響比較大,因此需消除測量基準與設(shè)計基準不重合對偏差評定的影響[1-2]。輪廓度評定時需要測量坐標與理論坐標的位置一致,按照最小區(qū)域原則[3]進行葉身型線的輪廓度誤差評定,但實現(xiàn)葉身型線最小條件比較困難[4]。

在評定過程中有兩個核心問題需要解決:實測點所對應(yīng)理論輪廓的最近點計算;如何統(tǒng)一實測基準與設(shè)計基準。對于第一個問題,一般對理論型線采用樣條函數(shù)進行擬合,對擬合的理論型線采用分割逼近[5]的方法轉(zhuǎn)換為點與點的搜索。對于第二個即基準統(tǒng)一的問題,本質(zhì)是剛性變換求解的問題,目前求解此問題主要有兩種方法:第一種方法是采用進化算法求解變換參數(shù),第二種方法是采用二維最近點迭代算法(ICP)進行求解。采用進化算法求解變換參數(shù)需要構(gòu)建不同的適應(yīng)度函數(shù),如遺傳算法[5-8]或者粒子群算法[9]求解基準變換,理論上可以求到變換參數(shù)的準確解。進化算法求解齊次變換,需要生成大量初始解,計算量大導致求解配準參數(shù)耗費時間較多。由于進化算法的隨機性,導致配準計算的穩(wěn)定性很差,可能陷入局部最優(yōu)。依據(jù)形位誤差評定原則最小區(qū)域法構(gòu)建適應(yīng)度函數(shù),計算出的旋轉(zhuǎn)量與理論值差異較大,而平移量較為準確。

許多學者依據(jù)文獻[10]提出的最近點迭代算法,提出用二維形式的最近點迭代算法計算實測輪廓與理論輪廓的配準參數(shù)。文獻[11]根據(jù)不同區(qū)域的公差不同,引入點對權(quán)值約束,實現(xiàn)輪廓的高精度配準。文獻[12]引入熵的概念,通過最大化熵準則的最近點迭代算法,實現(xiàn)有噪聲和離群點的二維點集配準,但其只適用粗配準精度。文獻[13]針對平面零件二維輪廓與矢量圖形的配準問題,使用二維最近點迭代算法的解析形式進行輪廓間的精配準,可以較好地實現(xiàn)輪廓的位置調(diào)整,但是平移量計算誤差受實測輪廓影響較大。二維最近點迭代算法求解變換參數(shù)速度快,依據(jù)原則為最小二乘原則,旋轉(zhuǎn)量求解非常精確,但由于實測點與理論型線總是存在偏差,采用質(zhì)心平移的方法獲得的平移量誤差較大。

基于以上原因,本文提出一種最近點迭代結(jié)合粒子群算法的葉身型線基準誤差求解方法。首先采用非均勻有理B樣條(NURBS)曲線的矩陣形式對葉身型線進行擬合,利用區(qū)間搜索的方法尋找實測點所對應(yīng)理論輪廓上的最近點及偏差;其次運用二維最近點迭代的解析形式計算基準變換的旋轉(zhuǎn)量,然后對最近點迭代計算得到的平移量加微小擾動,以最短距離均值為優(yōu)化目標,使用自適應(yīng)權(quán)重粒子群算法重新計算平移量。本文對實測點進行基準變換后計算偏差,得到實測輪廓的輪廓度,設(shè)定理論變換向量構(gòu)造仿真數(shù)據(jù),以計算結(jié)果與理論值的差異驗證方法的有效性,并與目前自由曲線配準的方法進行對比。

1 葉身型線數(shù)學模型

葉身型線屬于自由曲線,理論形式以離散點坐標給出。非均勻有理B樣條的出現(xiàn)為標準解析形狀和自由曲線提供了統(tǒng)一的數(shù)學表示[14]。國際標準化組織(ISO)在1991年頒布的工業(yè)產(chǎn)品幾何定義的STEP標準中,NURBS被定義為唯一的自由型曲線曲面表示方法[15],因此本文采用NURBS對葉身型線進行描述。

常用的NURBS曲線的表達形式有兩種,一種是德布爾算法的遞推形式,另一種是矩陣形式。雖然德布爾算法概念清晰,但算法復雜,計算量大。曲線軌跡上任意點的求取速度受基函數(shù)計算影響較大,為提高曲線的計算速度,采用NURBS曲線矩陣形式計算曲線的系數(shù)[16]。NURBS曲線的非均勻性是指節(jié)點矢量與間距為任意值,可在不同區(qū)間上得到不同的混合函數(shù)形狀[17],曲線隨參數(shù)變化而變化,可以很好地對自由曲線進行描述[18]。

由n+1個控制頂點定義的k次NURBS曲線的數(shù)學定義[14]如下

(1)

式中:C(u)為NURBS曲線的參數(shù)形式;u為參數(shù);wi為權(quán)重因子;pi為控制點坐標向量;Ni,k(u)為k次B樣條基函數(shù),表達式如下

(2)

Ni,k(u)=

(3)

(4)

式中:Mi(t)為系數(shù)矩陣,表達式為

(5)

式中:m11、m13、m23、m33、m43、m44分別代替矩陣Mi中的元素。

NURBS曲線最終可以化為以下形式,根據(jù)理論坐標數(shù)據(jù)計算得到n-2組系數(shù)

(6)

2 葉身型線配準

在葉身型線測量過程中,由于測量基準與理論設(shè)計基準存在不重合誤差,對實測點偏差的計算精度影響較大,所以需要減小基準不重合對質(zhì)量評定的影響。對于基準間的剛性變換參數(shù)的求解問題,可以視為實測點集與理論輪廓對應(yīng)的最近點集間的配準問題。最近點迭代(ICP)算法的優(yōu)化目標符合形位誤差評定的最小二乘原則,算法求取旋轉(zhuǎn)量非常精確,由于存在加工誤差,求取的平移量與真值存在較大差異。因此,以ICP算法計算旋轉(zhuǎn)量,以所有實測點偏差的絕對值均值作為適應(yīng)度函數(shù),對ICP算法獲得的初始平移量加微小擾動,采用自適應(yīng)權(quán)重粒子群算法重新求取平移量。

2.1 實測點與理論輪廓的輪廓度計算

理論點集為有序點集,對理論點集進行NURBS矩陣形式擬合,計算實測點的偏差,公式如下

pd=sdmin

(7)

式中:pd為實測點的偏差;s為位置符號;dmin為實測點與理論輪廓對應(yīng)的最近點間的距離。

實測點的位置符號s可以用實測點與理論輪廓上對應(yīng)的兩個最近點組成的叉積進行判斷。實測點可能位于理論輪廓內(nèi),也可能位于理論輪廓外,實測點在理論輪廓外側(cè)時,s為1;在理論輪廓內(nèi)側(cè)時,s為-1。

實測點對應(yīng)的理論輪廓最近點計算步驟如下。

(1)求出實測點到理論輪廓的最近距離點所在段的索引。

(2)傳入最近點所在理論輪廓段的參數(shù),對參數(shù)進行區(qū)間劃分搜索最近點,如黃金分割或二分法。

當測量基準與理論設(shè)計基準統(tǒng)一時,輪廓度E定義為

E=max(pd)-min(pd)

(8)

式中:max(pd)為所有實測點偏差的最大值;min(pd)為所有實測點偏差的最小值。

2.2 改進的葉身型線配準算法

實測點集P,數(shù)目為m,理論點集Q0,數(shù)目為mQ,對理論點集Q0進行NURBS矩陣形式擬合。改進的葉身型線配準算法原理如下。

2.2.1 以二維最近點迭代的解析形式求解旋轉(zhuǎn)矩陣R和初始平移向量T0ICP算法依據(jù)的原則是最小二乘,優(yōu)化的目標為最短距離平方和最小,數(shù)學形式如下

(9)

式中:Qi為實測點集P對應(yīng)理論型線的最近點集Q中的第i個點;Pi為實測點集P的第i個點。

二維ICP算法的解析算法步驟如下。

Step1:對實測點集P搜索對應(yīng)理論型線的最近點集Q。

Step2:對點集P和點集Q去中心化,分別得到點集P1和點集Q1,去中心化公式如下

(10)

(11)

式中:Q1i為最近點集Q1中的第i個點;P1i為實測點集P1的第i個點。

Step3:計算旋轉(zhuǎn)矩陣R。二維點集配準旋轉(zhuǎn)矩陣R由一個旋轉(zhuǎn)角度θ構(gòu)成,公式如下

(12)

將式(12)代入式(11),得到

(13)

對式(13)所示的函數(shù)進行求導,得到

(14)

求得f(θ)的極值點

(15)

Step4:求取平移向量

(16)

Step5:更新實測點集P,更新公式如下

P=PR+T0

(17)

Step6:重復Step1～Step5,直到算法收斂,得到基準變換的旋轉(zhuǎn)矩陣R和初始平移向量T0。

2.2.2 采用自適應(yīng)權(quán)重粒子群算法重新計算基準變換的平移量T由于存在偏差,ICP算法難以精確獲得基準變換的平移量,建立如下適應(yīng)度函數(shù)

(18)

式中:pdi為第i個實測點的偏差。

運用自適應(yīng)權(quán)重粒子群算法[20],以T0加微小擾動生成初始解空間,重新計算平移向量,計算步驟如下。

Step1:設(shè)置算法參數(shù),種群數(shù)量N為30,自我學習因子c1為0.9,群體學習因子c2為0.9,最大慣性權(quán)重wmax為0.8,最小慣性權(quán)重wmin為0.2。

Step2:初始化解空間和速度,計算所有粒子的適應(yīng)度值,并且找出解空間的最佳個體zg和最佳適應(yīng)度值fz

(19)

式中:xi為解空間的第i個解向量;vi為解空間的第i個解的速度向量;N1為均值為0、方差為0.04的正態(tài)分布生成的2維度的隨機向量;N2為均值為0、方差為1的正態(tài)分布生成的2維度的隨機向量。

Step3:進入主循環(huán),計算種群的適應(yīng)度,并且記錄種群最小適應(yīng)度fmin及種群適應(yīng)度的均值favg。計算個體的慣性權(quán)重w,更新速度及個體,更新規(guī)則如下

(20)

式中:f為第i個粒子的適應(yīng)度值;zi為第i個粒子在迭代過程中的歷史最佳個體。

Step4:最后比較個體更新前后的適應(yīng)度值,如果更新后的個體適應(yīng)度更小,則保留更新后的個體,否則保留更新前的個體。計算更新后的種群適應(yīng)度值,尋找出最佳個體及其適應(yīng)度值代替zg和fz。

Step5:重復Step3～Step4,直到收斂或者固定迭代次數(shù),輸出最佳個體zg,即重新計算的平移量T。

計算出旋轉(zhuǎn)矩陣R和平移量T后,對實測點集進行基準變換,此時實測型線與理論型線達到最佳配準。

3 算法驗證及實例分析

3.1 樣本構(gòu)造

本文樣本的構(gòu)造方法如下。以理論數(shù)據(jù)Q0(有序點集,數(shù)目為200)進行NURBS擬合,得到199段的理論輪廓。以正態(tài)分布N(u,σ2)生成兩組199個偏差pd,其中u為0,σ為0.03。根據(jù)參數(shù)曲線的等距公式生成偏距點集P0,每個偏距點的位置取理論輪廓曲線段的中間位置。構(gòu)建基準變換向量T1(0.02,0.1,-0.1),其中0.02為旋轉(zhuǎn)角度,0.1為x方向平移量,-0.1為y方向平移量。將點集P0進行T1的基準變換,變換后的點集P2作為實測數(shù)據(jù)?；鶞首儞QT1為理論基準到實測基準,而配準是求實測基準到理論基準的變換向量,所以配準參數(shù)的理論值應(yīng)為基準變換向量T1構(gòu)成的齊次剛性變換矩陣的逆矩陣對應(yīng)的基準變換向量。配準的理論基準變換向量為(-0.02,-0.098,0.102),樣本1偏差極差為0.215 4,樣本2偏差極差為0.163 6。本文算法實驗在Matlab2017b中進行。

3.2 不同適應(yīng)度函數(shù)求解結(jié)果的對比

遺傳算法和粒子群算法同屬于進化算法,自適應(yīng)權(quán)重粒子群算法簡單,收斂較快,是一種全局尋優(yōu)的進化算法。本文分別以式(8)(18)作為適應(yīng)度函數(shù),以自適應(yīng)權(quán)重粒子群算法優(yōu)化上述適應(yīng)度函數(shù),計算實測基準到理論基準的變換向量。以求出的變換向量與理論變換向量的歐式距離度量基準轉(zhuǎn)換的誤差,樣本數(shù)據(jù)及配準結(jié)果如圖1所示,計算結(jié)果如表1所示。

(a)樣本1數(shù)據(jù)(b)樣本2數(shù)據(jù)

(c)輪廓度配準結(jié)果(樣本1)(d)輪廓度配準結(jié)果(樣本2)

(e)最短距離均值配準結(jié)果(樣本1)(f)最短距離均值配準結(jié)果(樣本2)

圖1 不同適應(yīng)度函數(shù)樣本數(shù)據(jù)及配準結(jié)果

綜上所求解的基準變換向量理論值為(-0.020,-0.098,0.102),樣本1的偏差極差為0.215 4,即樣本1的輪廓度理論值為0.215 4;樣本2的偏差極差為0.163 6,即樣本2的輪廓度理論值為0.163 6。樣本1(輪廓度)的輪廓度計算值與理論值差異為-0.032 5,樣本1(最短距離均值)的輪廓度計算值與理論值差異為-0.015 8。樣本2(輪廓度)的輪廓度計算值與理論值差異為0.011 4,樣本2(最短距離均值)的輪廓度計算值與理論值差異為0.000 6。

綜合樣本1和樣本2的計算結(jié)果,經(jīng)以最短距離均值最小為優(yōu)化目標,求得的基準變換向量更接近理論值,所求的輪廓度也更接近理論值。因此,本文方法選擇最短距離均值為適應(yīng)度函數(shù),重新求取平移量。

3.3 本文算法與ICP算法的對比

標準ICP算法采用SVD分解去中心化后的點集的協(xié)方差矩陣。本文推導了二維ICP算法的解析形式,避免出現(xiàn)協(xié)方差矩陣分解失敗的情況。ICP算法依據(jù)最小二乘原則,求解剛性變換參數(shù)非?？?但由于存在加工誤差,平移量計算誤差較大。根據(jù)3.2小節(jié),以最短距離均值作為優(yōu)化目標,計算結(jié)果更加貼近真值,因此采用自適應(yīng)權(quán)重粒子群算法重新求取平移量。由于減小了優(yōu)化的維度并且得到了初始的平移量,算法計算速度更快,同時保證了精度。分別以ICP算法和本文算法對樣本1和樣本2進行配準,以計算所得變換向量與理論變換向量的歐式距離度量基準轉(zhuǎn)換誤差,樣本數(shù)據(jù)集配準結(jié)果如圖2所示,計算結(jié)果如表2所示。

由表2可以看出,采用ICP算法雖然計算速度快,但計算出的基準變換向量誤差較大,計算所得的輪廓度明顯偏離理論值,主要影響因素為平移量計算不準確。本文算法基準變換向量計算值與理論值的誤差在0.006左右,樣本1的輪廓度計算相對誤差為2.8%,樣本2的輪廓度計算相對誤差為1.2%,計算結(jié)果已經(jīng)非常接近理論值。與采用進化算法相比,本文所提算法直接計算所有變換參數(shù)的速度要快40%左右。

(a)ICP算法配準結(jié)果(樣本1)(b)ICP算法配準結(jié)果(樣本2)

(c)本文算法配準結(jié)果(樣本1)(d)本文算法配準結(jié)果(樣本2)

圖2 ICP算法與本文算法配準結(jié)果的對比

表2 ICP算法與本文算法配準計算結(jié)果的比較

4 結(jié) 論

本文采用NURBS的矩陣形式對葉身型線進行擬合,一次性計算出擬合曲線的所有參數(shù),避免了德布爾算法遞歸形式的大量計算,采用區(qū)間搜索求解點到曲線的最近點,計算速度快,精度高。此外,提出了一種葉身型線配準算法,首先通過二維ICP算法的解析形式精確獲得配準的旋轉(zhuǎn)量,其次以實測點集最短距離均值作為適應(yīng)度函數(shù),采用自適應(yīng)權(quán)重粒子群算法重新優(yōu)化平移量,基準變換向量誤差在0.006左右,相對于采用進化算法直接計算配準參數(shù)的速度提高了40%左右。所提算法適用于三坐標測量葉身型線數(shù)據(jù)后對葉身型線進行評價,能夠有效減小由于基準不重合對偏差計算的影響,避免葉身型線加工質(zhì)量誤判。