基于變形模板的復雜截面輪廓重構方法研究

玉榮， 徐 金 亭， 孫 玉 文＊

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

0 引 言

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

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

1 變形截面輪廓重構的數學模型

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

一般情況下，截面模板曲線可由B樣條曲線描述，其數學表達式為

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

而變形模板與實測點之間重構就是實現測量點到模板曲線間的誤差最小.根據式（2）及最小二乘原理，變形模板與實測點之間重構的目標函數可寫為

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

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

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

2 實測數據的精確配準及優化算法

2.1 實測數據與模板曲線間的配準模型

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

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

式中：εiter為給定計算精度.如果終止準則成立，則結束迭代，得到最優定位變換Ropt、Topt；否則迭代繼續進行.經過上述優化處理，就實現了測量數據與模板曲線之間的精確定位.

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

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

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

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

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

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

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

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

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

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

2.3 旋轉和平移變換矩陣的計算

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

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

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

3 模板曲線的約束變形

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

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

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

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

可以通過計算矩陣N（u）的廣義逆即N（u）的方式求解上述最小二乘問題，即Δd＝N（u）ΔD，從而使模板曲線控制點擾動量Δd達到最小［16］.對于秩為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）的秩分解并不是唯一的，但任意兩種分解所得的廣義逆矩陣卻是相同的，上述求解方法既可用于欠定線性方程組，也可用于超定線性方程組.

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

式中：D0為模板曲線D（u）上參數為u0的點為變形后模板曲線上相應的目標點；ΔD0為目標點與初始點的位移.對于單點約束，由于矩陣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 實例驗證

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

圖3 初始截形配準Fig.3 Initial match between template curve and measured sectional points

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

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

5 結 論

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

［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］朱心雄.自由曲線曲面造型技術 ［M］.北京：科學出版社，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.