渦輪葉片密集點(diǎn)云數(shù)據(jù)與CAD模型配準(zhǔn)方法

黃勝利 卜 昆 程云勇 周麗敏

西北工業(yè)大學(xué)現(xiàn)代設(shè)計(jì)與集成制造技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,西安,710072

0 引言

在渦輪葉片模具型腔反變形優(yōu)化設(shè)計(jì)中,計(jì)算分析葉片測(cè)量數(shù)據(jù)相對(duì)于CAD模型的三維偏差至關(guān)重要。葉片測(cè)量數(shù)據(jù)的主要獲取方法為三坐標(biāo)測(cè)量和光學(xué)掃描測(cè)量。近年來,隨著激光技術(shù)的發(fā)展,光學(xué)掃描儀獲取數(shù)據(jù)的速度和精度都有顯著的提高,應(yīng)用越來越廣泛,但其獲得的海量數(shù)據(jù)給后續(xù)的數(shù)據(jù)處理、分析、計(jì)算帶來很大的不便。由于測(cè)量數(shù)據(jù)是在測(cè)量設(shè)備坐標(biāo)系下獲得的,故在對(duì)測(cè)量數(shù)據(jù)進(jìn)行偏差分析前,需把測(cè)量數(shù)據(jù)與CAD模型的坐標(biāo)系統(tǒng)一起來,也就是實(shí)現(xiàn)測(cè)量數(shù)據(jù)與CAD模型的配準(zhǔn)定位。

目前應(yīng)用比較廣泛的配準(zhǔn)算法是Besl等[1]提出的迭代最近點(diǎn)(iterative closest point,ICP)算法。該方法通過最小二乘法構(gòu)造配準(zhǔn)的目標(biāo)函數(shù),具有收斂速度快、配準(zhǔn)精度高的優(yōu)點(diǎn)。基于ICP的改進(jìn)算法[2]在含有自由曲面的工件配準(zhǔn)中得到了廣泛的應(yīng)用。ICP算法要求待配準(zhǔn)兩模型的初始位置要足夠接近,否則不能得到好的收斂結(jié)果,因此,使用ICP算法進(jìn)行精確配準(zhǔn)前,需要對(duì)兩待配準(zhǔn)模型進(jìn)行預(yù)配準(zhǔn)。武殿梁等[3]采用遺傳算法實(shí)現(xiàn)了曲面之間的預(yù)配準(zhǔn),此方法比較耗時(shí)。吳鋒等[4]基于輪廓的力矩主軸實(shí)現(xiàn)了醫(yī)學(xué)圖像的預(yù)配準(zhǔn),該方法的缺點(diǎn)是對(duì)數(shù)據(jù)的缺失比較敏感,需要整個(gè)物體的全部信息。劉晶[5]通過建立模型的包圍盒實(shí)現(xiàn)三維測(cè)量模型與CAD模型之間的預(yù)配準(zhǔn),該方法對(duì)測(cè)量模型的完整性要求較高,且包圍盒的建立有一定的難度。

光學(xué)測(cè)量的密集點(diǎn)云數(shù)據(jù)的配準(zhǔn)計(jì)算非常耗時(shí),為了縮短配準(zhǔn)計(jì)算的時(shí)間,本文提出一種基于簡(jiǎn)化模型曲率計(jì)算的配準(zhǔn)方法,首先對(duì)密集點(diǎn)云模型進(jìn)行簡(jiǎn)化,估算簡(jiǎn)化模型的曲率,根據(jù)曲率提取對(duì)應(yīng)特征面;然后尋找匹配點(diǎn)對(duì)從而計(jì)算出變換矩陣,把求出的變換矩陣應(yīng)用于簡(jiǎn)化前的原始點(diǎn)云模型以實(shí)現(xiàn)模型的預(yù)配準(zhǔn);再通過奇異值分解和最近點(diǎn)迭代相結(jié)合(singular value decompo-sition and iterative closest point,SVD-ICP)算法[6-7]進(jìn)行精確配準(zhǔn);最后以某型渦輪葉片的蠟?zāi)槔?驗(yàn)證方法的有效性和實(shí)用性。

1 點(diǎn)云模型簡(jiǎn)化

由于光學(xué)測(cè)量獲取的數(shù)據(jù)量巨大,為了提高計(jì)算的速度需對(duì)點(diǎn)云模型進(jìn)行適當(dāng)?shù)暮?jiǎn)化。圖1所示分別為某型葉片及其蠟?zāi)５墓鈱W(xué)測(cè)量模型。目前,常用的點(diǎn)云簡(jiǎn)化方法有包圍盒法、均勻網(wǎng)格法、聚類法、隨機(jī)采樣法、曲率采樣法、均勻采樣法等[8],它們各有優(yōu)缺點(diǎn)和不同的適用范圍。本文對(duì)點(diǎn)云模型進(jìn)行簡(jiǎn)化的主要目的是提高計(jì)算的速度,因此采用均勻采樣法。

圖1 光學(xué)測(cè)量密集點(diǎn)云模型

2 曲率估算

曲率是曲面的重要幾何特征,具有平移、旋轉(zhuǎn)和縮放不變性,可以作為曲面特征識(shí)別的重要依據(jù)。本文采用鄰域內(nèi)擬合二次參數(shù)曲面的方法估算散亂點(diǎn)云的曲率[9]。

2.1 散亂數(shù)據(jù)參數(shù)化

對(duì)原始點(diǎn)云進(jìn)行空間劃分,計(jì)算點(diǎn)云中每個(gè)點(diǎn)pi的K 鄰域,K鄰域記為Knbhd(pi)。每個(gè)Knbhd(pi)可以擬合成一曲面且此曲面可以被一個(gè)局部基面逼近,該局部基面記為L(zhǎng)BS(Knbhd(pi)),本文選擇點(diǎn)pi處的切平面為局部基面。切平面可通過各點(diǎn)的鄰近關(guān)系和最小二乘原理[10]求得。對(duì)于不同的模型,合理選取K的值(K一般取10～30),以保證鄰域內(nèi)凹凸一致,這樣得到的切平面才能更好地逼近原始點(diǎn)集,使投影點(diǎn)在局部基面的參數(shù)化更好地反映點(diǎn)云的參數(shù)化。

把Knbhd(pi)中的各點(diǎn)在相應(yīng)局部基面上的投影點(diǎn)記為PRO(Knbhd(pi)),投影點(diǎn)一般是隨機(jī)分布的。為了參數(shù)化PRO(Knbhd(pi)),首先需要確定 u、v兩個(gè)方向。在 PRO(Knbhd(pi))中 ,p′i是pi在局部基面上的投影點(diǎn),求出距離p′i最遠(yuǎn)的點(diǎn)p′,將連接p′i、p′兩點(diǎn)的直線方向作為u 方向,垂直于u方向的直線方向作為v方向。這里v方向可由u矢量和局部基面的法矢作向量積求得。然后將PRO(Knbhd(pi))中的每一點(diǎn)p′j(j=1,2,…,K)都與p′i相連,得到K 條有向線段,把這些有向線段分別與u矢量作點(diǎn)積,記為Dj,把得到的點(diǎn)積進(jìn)行排序,最大值為Dmax,最小值為Dmin。由式(1)就可求得PRO(Knbhd(pi))中各點(diǎn)的u參數(shù)值:

其中,j=2,3,…,K+1且U1=0。同理可以求得各點(diǎn)的v參數(shù)值。把所求的PRO(Knbhd(pi))中各點(diǎn)的參數(shù)值作為Knbhd(pi)中對(duì)應(yīng)點(diǎn)的參數(shù),就實(shí)現(xiàn)了局部散亂數(shù)據(jù)的參數(shù)化。

2.2 擬合二次參數(shù)曲面

由2.1節(jié)中方法得到每個(gè)點(diǎn)pi的Knbhd(pi)的參數(shù)u、v,就可以建立二次參數(shù)曲面:

其中,A為3×3的系數(shù)矩陣。運(yùn)用Yang[11]提出的二次參數(shù)曲面逼近法求得散亂數(shù)據(jù)點(diǎn)的二次參數(shù)曲面。首先引入新的矢量表達(dá)式:

則二次參數(shù)曲面可表示為

其中,系數(shù)a、b、c的下標(biāo)為式(2)中A的行號(hào)和列號(hào)。采用最小二乘法使Knbhd(pi)上各點(diǎn)到二次參數(shù)曲面的歐幾里德距離之和最小,可得

式中,X為鄰域Knbhd(pi)的矩陣。

2.3 曲率計(jì)算

由于計(jì)算得到的曲面法矢方向一般是不協(xié)調(diào)的,即并不是所有的法矢都指向曲面的同一側(cè),因此有必要對(duì)求得的法矢進(jìn)行歸一化調(diào)整,具體方法參見文獻(xiàn)[9]。

根據(jù)曲面的曲率性質(zhì)可以得出高斯曲率:

平均曲率:

上述物理量可由曲面的第一、第二基本形式推出。

3 提取特征面

渦輪葉片主要由自由曲面和一些平面組成,而自由曲面的曲率特征沒有一定的規(guī)律性,因此,曲面的提取有一定的難度。本文中,將測(cè)量模型和CAD模型上某一組對(duì)應(yīng)平面作為尋找匹配點(diǎn)對(duì)的依據(jù)。高斯曲率I和平均曲率H均分別反映了曲面的凸凹形狀,I和H的不同組合表現(xiàn)的是不同的幾何形狀,表1給出了幾種常見的曲面形狀與I、H組合的關(guān)系。

表1 點(diǎn)的局部曲面類型

由表1可以看出,當(dāng)I、H 都為0時(shí),曲面就是一平面,根據(jù)此條件可提取模型的某一平面為特征面,為尋找對(duì)應(yīng)匹配點(diǎn)作好準(zhǔn)備。

4 模型配準(zhǔn)

4.1 預(yù)配準(zhǔn)

為了分析鑄件相對(duì)于其CAD設(shè)計(jì)模型的彎扭變形,首先需要實(shí)現(xiàn)測(cè)量模型與CAD模型的配準(zhǔn)。運(yùn)用SVD-ICP算法進(jìn)行精確配準(zhǔn),要想獲得較好的效果,首先需要實(shí)現(xiàn)測(cè)量模型與CAD模型的預(yù)配準(zhǔn)。本文通過曲率提取特征面的方法實(shí)現(xiàn)預(yù)配準(zhǔn),具體方法如下:

(1)運(yùn)用均勻采樣法對(duì)原始測(cè)量點(diǎn)云{P}進(jìn)行簡(jiǎn)化得到{P′}。

(2)計(jì)算簡(jiǎn)化后點(diǎn)云{P′}和CAD模型的重心(分別記為O1、O2),平移兩模型,使兩模型的重心與坐標(biāo)原點(diǎn)O重合。

(3)運(yùn)用第2節(jié)中所述方法估算簡(jiǎn)化后測(cè)量模型和CAD模型的高斯曲率K與平均曲率H。

(4)根據(jù)所求得的K和H提取兩模型中相對(duì)應(yīng)的特征面,分別記為L(zhǎng)1、L2,并求取兩特征面的

形心 P1、Q1。

(5)把P1、Q1沿所在特征面的法線方向移動(dòng)一定距離 d(d值可適當(dāng)選取,本文實(shí)驗(yàn)中 d=10mm),得點(diǎn)P2、Q2的計(jì)算公式:

式中,P1、P2、Q1、Q2為點(diǎn) P1、P2、Q1、Q2的位置矢量;N1、N2分別為特征面的單位法矢。

(6)由三組對(duì)應(yīng)點(diǎn){O,P1,P2}和{O,Q1,Q2}求取旋轉(zhuǎn)矩陣R和平移矩陣T,方法如下:以O(shè)點(diǎn)為坐標(biāo)系原點(diǎn),由單位矢量組 e1 、e2 、e3 與e′1 、e′2 、e′3按右手法則構(gòu)建兩個(gè)局部坐標(biāo)系;設(shè)存在旋轉(zhuǎn)變換矩陣R與平移變換矩陣T使兩個(gè)局部坐標(biāo)系重合,于是有

由式(13)可求得旋轉(zhuǎn)變換矩陣:

待配準(zhǔn)模型中的任一點(diǎn)Pi變換到對(duì)應(yīng)點(diǎn)Qi的關(guān)系式為

把P1、Q1代入式(15)可得到平移變換公式:

(7)把求得的旋轉(zhuǎn)矩陣和平移矩陣應(yīng)用到簡(jiǎn)化前的原始測(cè)量點(diǎn)云模型上,得到變換后的測(cè)量模型,再按向量OO2進(jìn)行整體平移就實(shí)現(xiàn)了測(cè)量模型與CAD模型的預(yù)配準(zhǔn)。

4.2 精確配準(zhǔn)

實(shí)現(xiàn)了測(cè)量模型與CAD模型的預(yù)配準(zhǔn)后,就可以運(yùn)用SVD-ICP算法進(jìn)行精確配準(zhǔn)。記測(cè)量點(diǎn)云數(shù)據(jù)為{pi},與其對(duì)應(yīng)的CAD模型上的對(duì)應(yīng)點(diǎn)集為{p′i}。配準(zhǔn)目標(biāo)函數(shù)表示為

對(duì)配準(zhǔn)點(diǎn)集{pi}和{p′i}運(yùn)用文獻(xiàn)[7]中奇異值分解法求取變換矩陣R和T,并將變換矩陣作用到測(cè)量點(diǎn)云數(shù)據(jù)上,迭代地進(jìn)行這個(gè)操作,直到使目標(biāo)函數(shù)值滿足設(shè)定的閾值為止。這樣就實(shí)現(xiàn)了測(cè)量點(diǎn)云模型與CAD模型坐標(biāo)系的統(tǒng)一,即實(shí)現(xiàn)了測(cè)量模型到CAD模型的精確配準(zhǔn)。

5 實(shí)例驗(yàn)證

為了驗(yàn)證本文方法的有效性,以VC++6.0為平臺(tái),以VTK可視化類庫(kù)為輔助工具編程實(shí)現(xiàn)了上述方法,并以某型渦輪葉片的蠟?zāi)槔M(jìn)行了驗(yàn)證。試驗(yàn)中所用的測(cè)量點(diǎn)云模型是由ATOS流動(dòng)式光學(xué)掃描儀獲取的,其數(shù)據(jù)點(diǎn)數(shù)為102 908。試驗(yàn)所用的電腦配置為Intel(R)Core(TM)2 Duo CPU E6550@2.33GHz;內(nèi)存為2.00GB。圖2所示為蠟?zāi)y(cè)量模型與CAD模型的初始位置,圖3所示為經(jīng)上述方法計(jì)算得到的匹配點(diǎn)對(duì)。圖4和圖5分別為預(yù)配準(zhǔn)和精配準(zhǔn)的結(jié)果圖。

圖2 蠟?zāi)Ｅ錅?zhǔn)前初始位置

圖3 對(duì)應(yīng)匹配點(diǎn)對(duì)

圖4 預(yù)配準(zhǔn)結(jié)果

圖5 精配準(zhǔn)結(jié)果

運(yùn)用本文方法對(duì)原始蠟?zāi)y(cè)量模型進(jìn)行不同程度的簡(jiǎn)化,實(shí)現(xiàn)預(yù)配準(zhǔn)和精確配準(zhǔn),并輸出配準(zhǔn)后的數(shù)據(jù),通過Geomagic Studio 7.0對(duì)配準(zhǔn)后的測(cè)量數(shù)據(jù)和CAD模型進(jìn)行三維偏差計(jì)算,結(jié)果如圖6所示。表2記錄了不同簡(jiǎn)化程度下,配準(zhǔn)計(jì)算所消耗的時(shí)間;表3記錄了不同簡(jiǎn)化程度下對(duì)應(yīng)的三維偏差結(jié)果。

圖6 三維偏差結(jié)果

表2 不同簡(jiǎn)化程度下配準(zhǔn)所耗時(shí)間對(duì)照表

表3 不同簡(jiǎn)化程度時(shí)三維偏差對(duì)照表

由表2可以看出,對(duì)原始測(cè)量點(diǎn)云數(shù)據(jù)進(jìn)行不同程度的簡(jiǎn)化可以減少預(yù)配準(zhǔn)所消耗的時(shí)間,精確配準(zhǔn)時(shí)由于仍然是對(duì)原始測(cè)量點(diǎn)云數(shù)據(jù)進(jìn)行計(jì)算,所以,此階段所耗時(shí)間是相同的。由表3可知,與由原始點(diǎn)云計(jì)算的偏差相比,簡(jiǎn)化到不同點(diǎn)數(shù)時(shí)得到的三維偏差結(jié)果相差很小,以簡(jiǎn)化到12 864個(gè)點(diǎn)時(shí)的結(jié)果來看,與以原始點(diǎn)云作為控制點(diǎn)集計(jì)算得出的結(jié)果相比,其配準(zhǔn)精度并沒有受到影響,而其預(yù)配準(zhǔn)階段的時(shí)間從30s降到了4s,可見速度有了很大的提高。綜合表2和表3的結(jié)果可知,對(duì)原始測(cè)量點(diǎn)云進(jìn)行不同程度的簡(jiǎn)化,可以提高配準(zhǔn)計(jì)算的速度,并且保證了配準(zhǔn)的精度。

6 結(jié)語(yǔ)

本文提出一種基于簡(jiǎn)化模型的曲率計(jì)算配準(zhǔn)方法。實(shí)例表明,該方法既保證了配準(zhǔn)的精度又提高了配準(zhǔn)的速度。但該方法只提高了預(yù)配準(zhǔn)階段的速度,并沒有縮短精確配準(zhǔn)階段所耗的時(shí)間。

[1]Besl P J,MaKay N D.A Method of Registration of 3D Shapes[J].IEEE Transaction on Pattern Analysis and Machine Intelligence,1992,14(2):239-255.

[2]Li Y D,Gu P H.Free-form Surface Inspection Techniques State of the Artreview[J].Computer-Aided Design,2004,36(3):1395-1417.

[3]武殿梁,黃海量,丁玉成,等.基于遺傳算法和最小二乘法的曲面匹配[J].航空學(xué)報(bào),2002,23(3):285-288.

[4]吳鋒,錢宗才,杭洽時(shí),等.基于輪廓的力矩主軸法在醫(yī)學(xué)圖像配準(zhǔn)中的應(yīng)用[J].第四軍醫(yī)大學(xué)學(xué)報(bào),2001,22(6):567-569.

[5]劉晶.葉片數(shù)字化檢測(cè)中的模型配準(zhǔn)技術(shù)及應(yīng)用研究[D].西安:西北工業(yè)大學(xué),2006.

[6]戴靜蘭,陳志楊,葉修梓.ICP算法在點(diǎn)云配準(zhǔn)中的應(yīng)用[J].中國(guó)圖像圖形學(xué)報(bào),2007,12(3):517-521.

[7]張二虎,卞正中,張燕,等.基于ICP和 SVD的視網(wǎng)膜圖像特征點(diǎn)配準(zhǔn)算法[J].小型微型計(jì)算機(jī)系統(tǒng),2004,25(10):1810-1813.

[8]戴靜蘭.海量點(diǎn)云預(yù)處理算法研究[D].杭州:浙江大學(xué),2006.

[9]賀美芳,周來水,神會(huì)存.散亂點(diǎn)云數(shù)據(jù)的曲率估算及應(yīng)用[J].南京航空航天大學(xué)學(xué)報(bào),2005,37(4):515-519.

[10]張麗艷,周儒榮,蔡煒斌,等.海量測(cè)量數(shù)據(jù)簡(jiǎn)化技術(shù)研究[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2001,13(11):1019-1023.

[11]Yang M,Lee E.Segmentation of Measured Datausing a Parametric Quadric Surface Approximation[J].Computer-Aided Design,1999,31(7):449-457.