玉榮， 徐 金 亭， 孫 玉 文＊

（1.大連理工大學(xué) 機械工程學(xué)院，遼寧 大連 116024；2.大連理工大學(xué) 汽車工程學(xué)院，遼寧 大連 116024）

0 引 言

基于測量數(shù)據(jù)的截面輪廓重構(gòu)一直是計算機輔助幾何設(shè)計、計算機圖形學(xué)領(lǐng)域的研究熱點［1－3］，也是各類復(fù)雜曲面構(gòu)造方法如蒙皮、掃掠、回轉(zhuǎn)及放樣等操作的前提基礎(chǔ).傳統(tǒng)的截面輪廓曲線重構(gòu)方法可分為兩類：對測量數(shù)據(jù)的插值和對測量數(shù)據(jù)的逼近.對散亂截面數(shù)據(jù)的插值或逼近而言，散亂數(shù)據(jù)的分割及序化處理、參數(shù)化方法的選擇、節(jié)點矢量的確定、曲線控制點的反算等都是截面輪廓曲線重構(gòu)的關(guān)鍵步驟，任何一步處理不當(dāng)都會導(dǎo)致最終的截面曲線無法反映真實的截面輪廓［4－5］.目前，多數(shù)研究集中于上述某一關(guān)鍵步驟的處理，如散亂數(shù)據(jù)的參數(shù)化、節(jié)點矢量的優(yōu)化等方面［6－7］.影響因素的多變性在帶來了造型的靈活性的同時，也增加了很大的不確定性，往往需要手工參與，促使新的造型技術(shù)的發(fā)展.最近，基于設(shè)計模板的輪廓曲線構(gòu)造方法在醫(yī)療圖像及制造工程中逐漸受到重視并成為新的研究熱點［8］.如，葉片等復(fù)雜曲面零件因在負(fù)載及應(yīng)力等工況下長期運行，會產(chǎn)生不同程度的變形，導(dǎo)致其實際的幾何形狀與名義幾何模型存在較大誤差.在獲取實際葉型截面輪廓時，如能以設(shè)計模板為基礎(chǔ)，通過模板變形進(jìn)行截面輪廓重構(gòu)，不僅可以避免實測散亂數(shù)據(jù)的序化處理、參數(shù)化和節(jié)點矢量優(yōu)化等傳統(tǒng)步驟，在快速實現(xiàn)散亂數(shù)據(jù)的精確輪廓重建的同時，還能建立起設(shè)計模板與實際截形之間的解析關(guān)系，精確判斷零件的變形情況，從而用于指導(dǎo)零件形狀的校正工藝，使得用于葉片等零件破損區(qū)域修復(fù)的截形更為合理.

當(dāng)前，基于變形模板的曲線重構(gòu)中，變形模板主要采用離散數(shù)據(jù)表達(dá)形式，有些B樣條形式表達(dá)的設(shè)計模板在變形時往往不能定量給定最大變形量的幅值，這給變形模板的精確表達(dá)和對模板變形的定量控制帶來不便.為實現(xiàn)這一目的，本文提出一種基于變形模板的復(fù)雜截面輪廓重構(gòu)新方法.在建立的模型中，采用能夠?qū)崿F(xiàn)單點或多點變形量精確控制的B樣條曲線變形策略.首先，通過精確的剛性配準(zhǔn)使得設(shè)計模板與實測點集盡量貼合，進(jìn)而采用精配準(zhǔn)與模板變形迭代進(jìn)行的方式，使模板與實測點集最大限度地重合，在將重構(gòu)誤差控制在允許范圍內(nèi)的同時使模板的變形最小化.

1 變形截面輪廓重構(gòu)的數(shù)學(xué)模型

基于模板變形的截面輪廓重構(gòu)是基于給定的截面模板曲線D（u）和對應(yīng)的實測零件截形數(shù)據(jù)點集｛Qi｝n1，在模板變形最小的條件下使得實測點與變形后的模板曲線間的幾何誤差最小.具體的重構(gòu)過程主要包括數(shù)據(jù)點和設(shè)計模板曲線的給定、設(shè)計模板與實測點的剛性配準(zhǔn)以及模板變形與配準(zhǔn)迭代優(yōu)化等主要過程.在本章中，將首先提出基于變形模板的截面輪廓重構(gòu)的數(shù)學(xué)模型.

一般情況下，截面模板曲線可由B樣條曲線描述，其數(shù)學(xué)表達(dá)式為

式中：di為第i個控制頂點；Ni，k（u）為B樣條曲線的基函數(shù)，k為樣條基的次數(shù).對于一測量點Qi，可定義其到模板曲線的距離函數(shù)如下式所示：

而變形模板與實測點之間重構(gòu)就是實現(xiàn)測量點到模板曲線間的誤差最小.根據(jù)式（2）及最小二乘原理，變形模板與實測點之間重構(gòu)的目標(biāo)函數(shù)可寫為

式中：f（Δd）為模型必須滿足模板變形的最小條件；Δd為模板變形后的控制點變動量；E為配準(zhǔn)時的變換矩陣為變形后的模板曲線，

利用上述模型進(jìn)行重構(gòu)時，首先采用剛性配準(zhǔn)使模板與實測點最大限度地貼合.在此基礎(chǔ)上，再采用模板變形與精確配準(zhǔn)迭代的策略實現(xiàn)模板與實測點的最終貼合.其中，在迭代中交替計算時限定模板變形量的約束點可按下式給出：

式中：Dj（u）為與測量數(shù)據(jù)Qj在模板曲線上的對應(yīng)點；Δt為一小量.一般情況下，模板曲線的變形可取一點或多點約束.

2 實測數(shù)據(jù)的精確配準(zhǔn)及優(yōu)化算法

2.1 實測數(shù)據(jù)與模板曲線間的配準(zhǔn)模型

在基于模板變形的曲線重構(gòu)中，實測數(shù)據(jù)與模板曲線間的精確配準(zhǔn)可以被描述為求解測量數(shù)據(jù)到模板曲線距離的平方和相對于剛體變換6個參數(shù)的極小值問題.因此，可以根據(jù)最小二乘原理構(gòu)造如下的目標(biāo)函數(shù)：

式中：D（ui）為Qi在給定曲線D（u）上的對應(yīng)點；R為由繞坐標(biāo)軸的3個旋轉(zhuǎn)量α、β、γ構(gòu)成的旋轉(zhuǎn)變換矩陣；T為沿坐標(biāo)軸的3個平移量tx、ty、tz構(gòu)成的平移矢量.目前所提出的優(yōu)化上述目標(biāo)函數(shù)的迭代方法，如ICP算法［9］、Menq算法［10］等大都采用輪換變量法，即計算點到曲線的最近點和求解旋轉(zhuǎn)、平移變換，對迭代過程進(jìn)行優(yōu)化處理，可在保證定位精度的同時加快迭代過程向全局最優(yōu)的收斂速度，特別適用于復(fù)雜曲線的配準(zhǔn)定位問題［11］，因此本文也采用輪換變量法.具體優(yōu)化過程如下：首先給定初始變換估計R0、T0，并將其作為迭代過程的初始值，進(jìn)而計算出點｛Qi｝n1在模板曲線上的最近點作為初始目標(biāo)點［11］；然后利用輪換變量法進(jìn)行迭代優(yōu)化.在每一次迭代過程中，都計算目標(biāo)函數(shù)ek，同時檢查迭代終止準(zhǔn)則

式中：εiter為給定計算精度.如果終止準(zhǔn)則成立，則結(jié)束迭代，得到最優(yōu)定位變換Ropt、Topt；否則迭代繼續(xù)進(jìn)行.經(jīng)過上述優(yōu)化處理，就實現(xiàn)了測量數(shù)據(jù)與模板曲線之間的精確定位.

2.2 點到B樣條曲線的最近點計算

在計算點到曲線的最近點時，通常采用的方法是數(shù)值迭代法.此類方法可達(dá)到很高的計算精度，但對初始值的要求卻比較苛刻，如果初始點選擇不當(dāng)，迭代過程可能會陷入局部極值甚至根本無法收斂.為解決這一問題，本文提出了一種新的基于伯恩斯坦多項式算術(shù)運算［12］和遞歸細(xì)分的最近點計算方法.

由于B樣條曲線可以采用節(jié)點插入的方式方便地轉(zhuǎn)化為Bézier曲線，不失一般性，下面將以Bézier曲線為例推導(dǎo)點到曲線最近點計算的數(shù)學(xué)模型，進(jìn)而再將其推廣到B樣條曲線.數(shù)學(xué)上，一條m次Bézier曲線可以表示為一段伯恩斯坦多項式函數(shù)曲線：

式中：bi為Bézier曲線r（u）的控制頂點；Bi，m（u）為伯恩斯坦基函數(shù)；u為曲線的參數(shù)；m為曲線的次數(shù).對于一點p0，它到曲線r（u）的平方距離函數(shù)可定義如下：

求解d（u）局部極小值的一般方法是，計算距離函數(shù)的導(dǎo)數(shù)函數(shù)（u）的所有零點，并檢查在這些零點處d（u）是否達(dá)到最小.對于非線性方程（u）＝0而言，可利用Newton－Raphson或其他數(shù)值迭代優(yōu)化方法進(jìn)行求解.然而，實際的計算表明此類迭代方法即便事先給定很好的初始值，仍有可能出現(xiàn)計算錯誤［13］.為此，本文將充分利用德卡斯特里奧算法和Bézier曲線的凸包性質(zhì)進(jìn)行求解，以避免求解過程對初始值的依賴.

利用伯恩斯坦多項式算術(shù)運算及基函數(shù)的規(guī)范性，可將（u）改寫為一伯恩斯坦多項式s（u）：

式中：gi∈R，為多項式s（u）的伯恩斯坦系數(shù).為了得到更為直觀的數(shù)學(xué)模型，利用伯恩斯坦基函數(shù)的線性精度性質(zhì)，函數(shù)s（u）在u參數(shù)軸上的圖形可由如下式所示的一條Bézier曲線表示：

式中：gi∈R2，是Bézier曲線rs（u）的控制頂點.

從上面的推導(dǎo)過程可以看到，如果方程（u）＝0成立，則曲線rs（u）必與u軸相交.這樣，點到曲線的最近點計算問題就轉(zhuǎn)化為直觀的曲線與參數(shù)軸之間交點的計算問題.對于B樣條曲線，可先將B樣條曲線轉(zhuǎn)化為一組Bézier曲線，然后對每條Bézier曲線應(yīng)用上述方法，通過比較計算確定點在B樣條曲線上的最近點.圖1給出了一個B樣條曲線Bézier曲線細(xì)分和最近點計算的算例.

圖1 B樣條曲線Bézier曲線細(xì)分及最近點計算的算例Fig.1 An example of B－spline curve subdivided by Bézier curve，and closest point calculation

2.3 旋轉(zhuǎn)和平移變換矩陣的計算

設(shè)Q＝ ｛Q1，Q2，…，Qn｝為測量點集，P＝｛P1，P2，…，Pn｝為Q在模板曲線上的對應(yīng)點集，其中Qi和Pi為一對匹配點.假設(shè)滿足式（6）的最小二乘解為Rr和Tr，則點集P′＝ ｛P′i｜P′i＝RrPi＋Tr｝和Q具有相同的質(zhì)心，即CP′＝CQ.其中

這樣，旋轉(zhuǎn)平移變換可分為兩步計算：（1）計算使式（13）取得最小值的旋轉(zhuǎn)變換R；（2）按T＝CQRCP計算平移矢量.本文將采用Arun等提出的采用奇異值分解法（SVD）［14］求解使方程式（13）取得最小值的旋轉(zhuǎn)變換矩陣R.首先計算P和Q之間的協(xié)方差矩陣

對H進(jìn)行奇異值分解可得H＝UΛVT，X＝VUT，計算行列式det（X），如果det（X）＝＋1，則旋轉(zhuǎn)矩陣R＝X［15］，求得旋轉(zhuǎn)矩陣后，進(jìn)而求解平移矢量T.

3 模板曲線的約束變形

B樣條曲線的單點約束、多點約束變形要通過計算曲線控制點的變動量來達(dá)到所要求的變形，這可由最小二乘法重新配置控制點來實現(xiàn)［15］，具體計算過程如下.將變形后的模板曲線寫成矩陣的形式，可以得到

式中：N（u）為基函數(shù)向量，N（u）＝（N0，k（u） …Nn，k（u））；d為控制點向量，d＝（d0…dn）；Δd為控制點擾動向量，Δd＝（Δd0… Δdn）T.

對于復(fù)雜的物體，一般需移動多點才能達(dá)到所要求的變形.假設(shè)曲線有m＋1個初始點，它們對應(yīng)的參數(shù)分別為uj（j＝0，1，…，m），則模板曲線多點約束方程可寫為

根據(jù)式 （15）中各變量的含義，并令 ΔD＝（ΔD0… ΔDm）T，則式（16）可改寫為如下的矩陣形式：

可以通過計算矩陣N（u）的廣義逆即N（u）的方式求解上述最小二乘問題，即Δd＝N（u）ΔD，從而使模板曲線控制點擾動量Δd達(dá)到最?。?6］.對于秩為r，（m＋1）×（n＋1）的矩陣N（u）而言，其秩分解為

式中：R是秩為r的（m＋1）×r矩陣；S是秩為r的r×（n＋1）矩陣.則矩陣N（u）的廣義逆可寫為

將式（19）代入Δd＝N（u）ΔD即可得到最小二乘意義下的控制點擾動量Δd.需要注意的是，矩陣N（u）的秩分解并不是唯一的，但任意兩種分解所得的廣義逆矩陣卻是相同的，上述求解方法既可用于欠定線性方程組，也可用于超定線性方程組.

特別地，當(dāng)是單點約束時，模板曲線約束方程可直接寫為

式中：D0為模板曲線D（u）上參數(shù)為u0的點為變形后模板曲線上相應(yīng)的目標(biāo)點；ΔD0為目標(biāo)點與初始點的位移.對于單點約束，由于矩陣N（u）的秩為1，它的廣義逆為

令ΔD＝ΔD0，N（u）＝N（u0），則式（17）的解可表示為

基于上面的計算公式，圖2（a）給出了單點約束時B樣條曲線變形，圖2（b）為多點約束時B樣條曲線變形的情況.

圖2 兩種B樣條曲線變形的算例Fig.2 Two types of B－spline deformation′s examples

4 實例驗證

為驗證所提出的基于變形模板的截面輪廓重構(gòu)模型及方法的有效性，本文針對葉片曲面進(jìn)行了重構(gòu)實驗.如圖3所示，算例中的截面輪廓為葉片截形.從圖中可以看到，配準(zhǔn)前設(shè)計模板與實際截形的空間位置相對誤差較大，設(shè)計模板和實際截形間的最大距離為9.78mm，截形弦長500 mm，重構(gòu)中實際截形的采樣點為100個.圖3給出了經(jīng)過配準(zhǔn)后的設(shè)計模板和實際截形的相對位置.圖4分別給出了基于模板變形的截面輪廓重構(gòu)過程中經(jīng)過第1、2、6、26、36次迭代時設(shè)計模板與實際截形的相對位置.從表1也可以看出，設(shè)計模板的形狀在不斷演化，且同實際截形的誤差（誤差以設(shè)計模板與實際截形間的最大距離表示）在逐步減小，這表明所提出方法具有充分的形狀變化適應(yīng)性.同時，由于模板變形和精確配準(zhǔn)的交替迭代使用，最大限度保證了模板在盡可能小的變形條件下實現(xiàn)與實際截形的貼合.所提算法已經(jīng)在Visual C＋＋中編碼實現(xiàn)，經(jīng)過36次迭代計算時，設(shè)計模板與實際截形之間最大誤差僅為0.001 mm.此外，該算法采用迭代求解的方式是為了減少設(shè)計模板的變形程度，同常規(guī)曲線重構(gòu)的迭代求解以減少擬合誤差一樣，會增加程序的運行時間，但運行時間（CPU 2.33GHz，內(nèi)存4GB）在毫秒級，可以滿足工程計算需求.

圖3 初始截形配準(zhǔn)Fig.3 Initial match between template curve and measured sectional points

圖4 截形配準(zhǔn)過程Fig.4 Different matching stages of template curve and measured sectional points

表1 不同迭代次數(shù)時與實際截形的誤差Tab.1 Error with practical sectional points on different iterations

5 結(jié) 論

本文提出了基于變形模板的復(fù)雜截面輪廓重構(gòu)模型和方法，相比于常規(guī)B樣條曲線重構(gòu)，其優(yōu)勢在于不必進(jìn)行實測數(shù)據(jù)點的分割、序化和數(shù)據(jù)參數(shù)化.所提出的方法通過精確配準(zhǔn)和模板變形迭代優(yōu)化策略，有效避免了并行優(yōu)化時的計算耗費過大等問題，并保證了模板在盡可能小的變形條件下實現(xiàn)與實際截形的貼合.采用的ICP精確配準(zhǔn)模型和奇異值分解算法保證了配準(zhǔn)的可靠性和精確性；單點和多點約束變形能夠有效控制設(shè)計模板變形的趨勢和幅值；提出的一種新的基于伯恩斯坦多項式算術(shù)運算和遞歸細(xì)分的最近點計算方法，能夠克服經(jīng)典迭代算法需要給定初始值且其精度依賴初始值的問題，計算精度高、速度快.下一步的工作是進(jìn)一步將所提出的重構(gòu)模型和方法應(yīng)用于復(fù)雜曲面零件破損區(qū)域的幾何修復(fù)中.

［1］Varady T，Martin R R，Cox J.Reverse engineering of geometric models－an introduction［J］.Computer－Aided Design，1997，29（4）：255－268.

［2］Weiss V，Andor L，Renner G，etal.Advanced surface fitting techniques ［J］.Computer Aided Geometric Design，2002，19（1）：19－42.

［3］LI Yong－qing，HUANG Xiao－ping，GONG Chenhe，etal.An engineering rule based parameterization approach for turbine blade reverse engineering ［C］／／ Proceedings of the Geometric Modeling and Processing.Beijing：IEEE Computer Society，2004：311－318.

［4］LI Yong－qing，HUANG Xiao－ping，GONG Chenhe，etal.Sketch template based parametric modeling in reverse engineering ［J］.Computer－Aided Design ＆ Applications，2005，2（1－4）：19－28.

［5］KE Ying－lin，ZHU Wei－dong，LIU feng－shan，etal.Constrained fitting for 2Dprofile－based reverse modeling［J］.CAD Computer Aided Design，2006，38（2）：101－114.

［6］XIE Hui， QIN Hong.A novel optimization approach to the effective computation of NURBS knots［J］.International Journal of Shape Modeling，2001，7（2）：199－227.

［7］ZHANG Cai－ming，HAN Hui－jian，Cheng F F.Determining knots by minimizing energy ［J］.Journal of Computer Science and Technology，2006，21（2）：261－264.

［8］LI Yong－qing，HUANG Xiao－ping，GONG Chenhe，etal.Sketch template based parametric modeling in reverse engineering ［J］.Computer－Aided Design ＆ Applications，2005，2（1－4）：19－28

［9］Besl P J，Mckay N D.A method for registration of 3－D shapes ［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1992，14（2）：239－256.

［10］Menq C，Yau H，Lai G.Automated precision measurement of surface profile in CAD－directed inspection［J］.IEEE Transactions on Robotics and Automation，1992，8（2）：268－278.

［11］Li Z，Gou J，Chu Y.Geometric algorithms for workpiece localization ［J］.IEEE Transactions on Robotics and Automation，1998，14（6）：864－878.

［12］Berchtold J，Bowyer A.Robust arithmetic for multivariate Bernstein－form polynomials ［J］.CAD Computer Aided Design，2000，32（11）：681－689.

［13］MA Ying－liang，Hewitt W T.Point inversion and projection for NURBS curve and surface：control polygon approach ［J］.Computer Aided Geometric Design，2003，20（2）：79－99.

［14］Arun K S，Huang T S，Blostein S D.Least－squares fitting of two 3－D point sets［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1987，9（5）：687－700.

［15］朱心雄.自由曲線曲面造型技術(shù) ［M］.北京：科學(xué)出版社，2000.ZHU Xin－xiong.Free－form Curves／Surface Modeling Technology［M］.Beijing：Science Press，2000.（in Chinese）

［16］Hus W M，Hughes J F，Kaufman H.Direct manipulation of freeform deformation ［J］.Computer Graphics，1992，26（2）：177－184.