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

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

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

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

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

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

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

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

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

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

(1)

式中:C(u)為NURBS曲線的參數(shù)形式;u為參數(shù);wi為權(quán)重因子;pi為控制點坐標(biāo)向量;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ù)理論坐標(biāo)數(shù)據(jù)計算得到n-2組系數(shù)

(6)

2 葉身型線配準(zhǔn)

在葉身型線測量過程中,由于測量基準(zhǔn)與理論設(shè)計基準(zhǔn)存在不重合誤差,對實測點偏差的計算精度影響較大,所以需要減小基準(zhǔn)不重合對質(zhì)量評定的影響。對于基準(zhǔn)間的剛性變換參數(shù)的求解問題,可以視為實測點集與理論輪廓對應(yīng)的最近點集間的配準(zhǔn)問題。最近點迭代(ICP)算法的優(yōu)化目標(biāo)符合形位誤差評定的最小二乘原則,算法求取旋轉(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ū)間劃分搜索最近點,如黃金分割或二分法。

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

E=max(pd)-min(pd)

(8)

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

2.2 改進的葉身型線配準(zhǔn)算法

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

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

(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。二維點集配準(zhǔn)旋轉(zhuǎn)矩陣R由一個旋轉(zhuǎn)角度θ構(gòu)成,公式如下

(12)

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

(13)

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

(14)

求得f(θ)的極值點

(15)

Step4:求取平移向量

(16)

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

P=PR+T0

(17)

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

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

(18)

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

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

Step1:設(shè)置算法參數(shù),種群數(shù)量N為30,自我學(xué)習(xí)因子c1為0.9,群體學(xué)習(xí)因子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:重復(fù)Step3～Step4,直到收斂或者固定迭代次數(shù),輸出最佳個體zg,即重新計算的平移量T。

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

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)建基準(zhǔn)變換向量T1(0.02,0.1,-0.1),其中0.02為旋轉(zhuǎn)角度,0.1為x方向平移量,-0.1為y方向平移量。將點集P0進行T1的基準(zhǔn)變換,變換后的點集P2作為實測數(shù)據(jù)?；鶞?zhǔn)變換T1為理論基準(zhǔn)到實測基準(zhǔn),而配準(zhǔn)是求實測基準(zhǔn)到理論基準(zhǔn)的變換向量,所以配準(zhǔn)參數(shù)的理論值應(yīng)為基準(zhǔn)變換向量T1構(gòu)成的齊次剛性變換矩陣的逆矩陣對應(yīng)的基準(zhǔn)變換向量。配準(zhǔn)的理論基準(zhǔn)變換向量為(-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ù),計算實測基準(zhǔn)到理論基準(zhǔn)的變換向量。以求出的變換向量與理論變換向量的歐式距離度量基準(zhǔn)轉(zhuǎn)換的誤差,樣本數(shù)據(jù)及配準(zhǔn)結(jié)果如圖1所示,計算結(jié)果如表1所示。

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

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

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

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

綜上所求解的基準(zhǔn)變換向量理論值為(-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)化目標(biāo),求得的基準(zhǔn)變換向量更接近理論值,所求的輪廓度也更接近理論值。因此,本文方法選擇最短距離均值為適應(yīng)度函數(shù),重新求取平移量。

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

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

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

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

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

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

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

4 結(jié) 論

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