陶　浩　何改云　王太勇　董甲甲　張永賓

天津大學(xué)機械工程學(xué)院，天津，300072

0　引言

在復(fù)雜曲線曲面的數(shù)控加工中，傳統(tǒng)的刀具軌跡生成方法是通過CAM系統(tǒng)生成G01代碼，再將其輸入到數(shù)控系統(tǒng)中進行直線插補加工。該方法原理簡單，但相鄰線段的銜接處易出現(xiàn)速度和加速度的波動，從而顯著降低零件的加工質(zhì)量。為了減小速度和加速度的波動，呂強等［1］提出了通過限制相鄰線段的銜接速度以實現(xiàn)速度平滑過渡的方法。近年來，隨著直接參數(shù)曲線插補方法的提出和發(fā)展，已有許多學(xué)者證明參數(shù)曲線插補方法能滿足目前數(shù)控加工行業(yè)高速高精的要求［2-5］。另外，相比G01代碼而言，參數(shù)曲線的數(shù)據(jù)量非常小，節(jié)省了數(shù)控系統(tǒng)的內(nèi)存空間。對于復(fù)雜曲線曲面的加工路徑而言，現(xiàn)有的CAM系統(tǒng)仍以產(chǎn)生G01代碼為主，因此，人們利用參數(shù)曲線對G01代碼進行擬合［6-8］，不但可壓縮數(shù)據(jù)量，還可從本質(zhì)上平滑加工路徑。無論是參數(shù)曲線插補還是擬合，NURBS和B樣條曲線的應(yīng)用最為廣泛。

數(shù)據(jù)點的擬合方法一般分為兩大類：插值擬合和逼近擬合。插值擬合精度較高，但不能縮減數(shù)據(jù)量；逼近擬合可大幅度縮減數(shù)據(jù)量，但為了在指定的精度內(nèi)逼近給定的數(shù)據(jù)點，必須通過迭代的方法不斷增加控制點（一般以最少的控制點個數(shù)開始擬合），因此逼近擬合效率極低。且當(dāng)控制點和節(jié)點的數(shù)量增加后，某些節(jié)點區(qū)間可能為空，會導(dǎo)致在逼近擬合求解時出現(xiàn)奇異方程組的情況，需進一步特殊處理。此外，關(guān)于NURBS曲線擬合時如何設(shè)定權(quán)值的研究很少，大多數(shù)情況下，簡單地將所有權(quán)值都取為1，即進行B樣條擬合（B樣條曲線實質(zhì)上是控制點權(quán)因子均為1的NURBS曲線）。且實際上對逼近擬合來說，允許其選取任意的權(quán)因子可能會生成具有較少控制點的曲線；但對插值擬合來說，其控制點個數(shù)固定，沒有合適的理由選取任意的權(quán)因子［9］。

YAU等［10］利用最小二乘法，將數(shù)據(jù)點逼近擬合成NURBS曲線。為了減少擬合時的數(shù)據(jù)量，PARK等［11］提出了對主導(dǎo)點進行擬合的方法。ZHANG等［12］同樣提出了幾何特征點的自適應(yīng)查尋方法。ZHAO等［13］提出了根據(jù)離散點曲率和弦高誤差等條件進行主導(dǎo)點選取的方法。TSAI等［14］提出了一種可實時具有C2連續(xù)性的局部Bezier曲線插值擬合方法，在逼近擬合時需計算數(shù)據(jù)點的擬合精度。傳統(tǒng)的牛頓迭代法雖精度較高，但效率較低。ZHU等［15］推導(dǎo)出了加工過程中實時計算輪廓跟隨誤差的Taylor二階展開算法，該算法可移植到數(shù)據(jù)點擬合時的精度校驗過程中，其計算速度較快且能達到合理的控制精度。

本文綜合主導(dǎo)點和局部Bezier快速插值擬合的方法，給出了主導(dǎo)點的確定方法，并在離線數(shù)據(jù)預(yù)處理階段，通過局部Bezier插值擬合和精度校驗不斷增加新的主導(dǎo)點，最后對主導(dǎo)點進行在線插值擬合以及誤差校驗，擬合后的曲線只需對非主導(dǎo)點進行1次誤差校驗循環(huán)，就能達到預(yù)期的擬合精度要求，從而加快了擬合速度。

1　主導(dǎo)點的種類與選擇

1.1　曲率極大值點和超過曲率閾值的點

很多文獻指出，曲線的曲率值是用來判斷曲線平滑性的一種有力工具［16］。但對于離散數(shù)據(jù)點而言，無法根據(jù)曲線的導(dǎo)數(shù)求取曲率的準(zhǔn)確值。依據(jù)文獻［17］提供的方法，離散數(shù)據(jù)點的曲率估算公式如下：

式中，Qi-1、Qi和 Qi+1為三個連續(xù)的數(shù)據(jù)點；KQi為點 Qi處的曲率估算值；A(Qi-1,Qi,Qi+1)表示Qi-1、Q和Qi+1三點構(gòu)成的三角形的面積。

文獻［18］提出了離散數(shù)據(jù)點曲率的另一種估算公式，通過推導(dǎo)發(fā)現(xiàn)，兩套計算公式的本質(zhì)一致。但文獻［18］中的公式涉及反三角函數(shù)的求解，增加了計算的復(fù)雜性。

圖1所示為對半蝴蝶狀復(fù)雜曲線除去首末2點后的269個離散數(shù)據(jù)點的曲率估算結(jié)果。由圖1可以看出，將連續(xù)曲線離散后進行曲率的估值計算，不可避免地會引起計算值的波動，在曲率較大處的效果尤為明顯。為了更加準(zhǔn)確地提取離散數(shù)據(jù)點的曲率極大值點，本文通過將當(dāng)前曲率值點分別與前后5個連續(xù)相鄰的曲率值點相比較，進而篩選出曲率極大值點（部分）。對于前后比較點個數(shù)的確定并沒有嚴(yán)格的要求，但如果選取個數(shù)過少（前后至少各1個），則無法很好地排除數(shù)據(jù)波動點；而如果選取個數(shù)過多，可能會漏選一部分曲率極大值點（但在曲率值較大點處可通過曲率閾值點補選）。本文選定此數(shù)目為5，并將曲率極大值處的數(shù)據(jù)點標(biāo)記為主導(dǎo)點。為了保證擬合曲線能夠較好地進行插值或逼近曲率值較大處的數(shù)據(jù)點，除了曲率極大值點外，本文還設(shè)定了1個曲率閾值Km，若離散數(shù)據(jù)點的曲率值大于設(shè)定的曲率閾值Km，則同樣將其標(biāo)記為曲線擬合時的主導(dǎo)點。這也是本文不單純使用當(dāng)前曲率值點與前后較少數(shù)量（如2個）曲率值點作比較的原因，因為單純應(yīng)用曲率極大值點選取法且前后比較數(shù)量為2時，會漏選圖1中的P1和P2點。

圖1　半蝴蝶狀復(fù)雜曲線離散數(shù)據(jù)點的曲率估算值Fig.1　Curvature estimation of discrete points of semi-butterfly curve path

1.2　曲線的拐點或反曲點

在數(shù)學(xué)上，曲線的拐點或反曲點是定義曲線凹凸性發(fā)生變化的轉(zhuǎn)折點。對于3次Bezier曲線而言，其二階導(dǎo)函數(shù)為一條直線，依據(jù)曲線拐點的判斷原則，3次Bezier曲線內(nèi)部最多只能有1個拐點，因此，在利用3次Bezier曲線對離散數(shù)據(jù)點進行逼近擬合時，離散數(shù)據(jù)點的拐點是決定曲線幾何形狀的重要主導(dǎo)點之一。

圖2中，設(shè)Ti為向量QiQi+1與QiQi+2的叉積，Ti+1為Qi+1Qi+2與Qi+1Qi+3的叉積，αT為向量 Ti和Ti+1的夾角。離散數(shù)據(jù)點拐點的計算方法滿足：

圖2　離散數(shù)據(jù)點的拐點計算Fig.2　The inflection point of discrete path points

對于平面數(shù)據(jù)點而言，若Qi+2為曲線的拐點，那么αT的值為π，反之，αT的值為0；如果是空間曲線，對于由CAM系統(tǒng)生成的連續(xù)復(fù)雜曲線的微小直線段路徑點而言，連續(xù)4點可近似看作位于同一平面內(nèi)，所以，若Qi+2為曲線的拐點，則αT為一個較大的值（接近于π），反之，αT為一個較小的值（接近于0）。由此可設(shè)定一個正參考值αthre，（如αthre=2/3π），若αT≥ αthre成立，則將Qi+2作為曲線的拐點處理，并標(biāo)記為主導(dǎo)點。

1.3　長度突變點

實驗驗證具有C2連續(xù)性的局部插值算法的要求之一是，相鄰微小直線段間的長度比值必須控制在一定范圍內(nèi)，否則插值生成的曲線會出現(xiàn)扭曲或尖點（圖3），這是由算法本身連續(xù)性的限制條件造成的。為了解決此矛盾，需通過長度均分策略增加數(shù)據(jù)主導(dǎo)點，以保證相鄰主導(dǎo)點線段之間的長度比值變化較均勻。本文定義了長度突變系數(shù)λ=，以及長度突變閾值λ（本實驗選取limλlim=3）。在長度均分過程中會出現(xiàn)圖4中的兩種情況，當(dāng)線段的長度突變系數(shù)λ大于設(shè)定閾值λlim時，需要對長線段進行均分處理。

由圖4a可以看出，第1段長度（第1個主導(dǎo)點與第2個主導(dǎo)點連線段長度）與第2段長度（第2個主導(dǎo)點與第3個主導(dǎo)點連線段長度）的比值大于λlim，預(yù)計將原始數(shù)據(jù)點下標(biāo)為istart,new=(istart+imiddle)/2的點設(shè)置為新的主導(dǎo)點，并將istart,new的值賦給istart。由于選取數(shù)據(jù)點下標(biāo)的中點并不能保證分割后的相鄰兩線段的長度接近相等，因此需對istart,new的值作進一步調(diào)整，以保證QistartQistart,new的長度和Qistart,newQiend的長度比值達到最合理。若istart,new=istart，則令 istart,new=imiddle。

圖3　長度不均勻造成的擬合尖角Fig.3　The sharp fitting angle caused by unevenness of the length between points

圖4　長度均分策略Fig.4　Length equalization strategy

由圖4b可以看出，第2段長度與第1段長度的比值大于λlim，預(yù)計將下標(biāo)為iend,new=(imiddle+iend)/2的點設(shè)置為新的主導(dǎo)點，并將iend,new的值賦給iend，同樣可通過進一步調(diào)整，保證QimiddleQiend,new的長度與Qiend,newQiend的長度比值達到最合理。若iend,new=imiddle成立，則將 imiddle賦值給 istart，否則，istart值不變。

通過不斷地循環(huán)迭代，可實現(xiàn)主導(dǎo)點間的長度比值達到預(yù)設(shè)范圍之內(nèi)的要求。由于迭代運算的原因，本文提出的查找主導(dǎo)點方法比較耗時，但是在離線數(shù)據(jù)預(yù)處理階段完成，并不會對在線擬合階段的效率造成影響。

2　C2連續(xù)的3次Bezier曲線局部插值算法

n次Bezier曲線C(u)的表達式如下：

式中，Pi為曲線控制點；Bi,n為Bernstein基函數(shù)。

n次Bezier曲線C(u)的一階求導(dǎo)公式為

根據(jù)式（6）可知，對于兩條相連的3次Bezier曲線，公共端點D處的一階導(dǎo)矢計算公式為

式中，D為相鄰兩條Bezier曲線C1(u)和C2(u)的連接點；為C1(u)的第3個控制點；為C2(u)的第2個控制點。

為了保證曲線C1(u)和C2(u)在連接點D處一階連續(xù)，聯(lián)立式（7）、式（8）可得

同理，兩條曲線在連接點D處的二階導(dǎo)矢計算公式為

同理，為了保證兩條曲線在連接點D處二階連續(xù)，聯(lián)立式（10）、式（11）可得

圖5　具有C2連續(xù)性的Bezier曲線局部插值過程Fig.5　Local Bezier curves interpolation fitting with C2 continuity

根據(jù)以上分析，為了保證所有分段Bezier曲線之間的C2連續(xù)，以下方程組必須成立：

通過進一步簡化，方程組（式（15））可改寫為

當(dāng)i=0和i=N時，令Q0=S0，QN=SN，則式（16）為一線性方程組，可將其改寫為矩陣形式。通過矩陣求解得到所有Si后，將其代入式（15），進而可求出所有Bezier曲線段的控制點Pi1和Pi2，且有 P(i)0=Di-1及 P(i)3=Di。

從以上求解Bezier曲線段控制點的過程中可發(fā)現(xiàn)，此方法不涉及Bernstein基函數(shù)的求解，大大減少了運算量，很適合計算機的快速求解。除了具有控制點求解快速、方便的優(yōu)點外，Bezier曲線還可同樣在不求解Bernstein基函數(shù)的情況下，快速計算曲線上對應(yīng)的參數(shù)點。

3　Bezier曲線擬合后的非主導(dǎo)點誤差檢測與輪廓誤差跟隨法

3.1　de Castejau算法簡述

本文提出的在數(shù)據(jù)預(yù)處理階段對主導(dǎo)點進行具有C2連續(xù)性的插值擬合算法的另一個優(yōu)點就是，在非主導(dǎo)點的誤差檢測過程中可充分利用de Castejau算法求解曲線及其導(dǎo)數(shù)曲線上的參數(shù)點，以提高非主導(dǎo)點誤差的計算速度。由于de Castejau算法是一種線性插值算法，在數(shù)值計算時受計算機浮點舍入誤差影響較小，因此也有利于提高計算精度。此外，由Bezier曲線的一階求導(dǎo)公式（式（6））可知，其導(dǎo)數(shù)曲線也是一條Bezier曲線，所以曲線及其導(dǎo)曲線上點的求取可利用同一套算法，從而精簡了程序的代碼量。

當(dāng)n=3時，3次Bezier曲線上參數(shù)點的求解展開式如下：

由式（17）可知u=0.5時的de Castejau計算過程（圖6）。

3.2　Bezier曲線擬合后的非主導(dǎo)點誤差計算

點到參數(shù)曲線的誤差計算通常使用牛頓迭代法，且誤差定義如下：

牛頓迭代法的初始值u，可通過對Bezier曲線段內(nèi)非主導(dǎo)數(shù)據(jù)點進行弦長參數(shù)化［9］求出。計算第i段Bezier曲線內(nèi)所有非主導(dǎo)點到曲線的距離誤差，并將超過設(shè)定閾值的誤差最大點設(shè)為新的主導(dǎo)點（圖7）；重新利用長度均分策略計算長度均分點，并再次將Bezier曲線擬合后計算出的非主導(dǎo)點中誤差超限且最大的點作為新的主導(dǎo)點，如此循環(huán)，直到所有非主導(dǎo)點的誤差達到預(yù)處理階段擬合誤差允許值的范圍內(nèi)。

圖6　de Castejau算法線性插值過程Fig.6　de Castejau linear interpolation process

圖7　非主導(dǎo)點誤差計算與新增主導(dǎo)點Fig.7　Non-dominant points error calculation and new added dominant point

3.3　輪廓誤差跟隨法計算B樣條擬合誤差

對主導(dǎo)點進行插值擬合生成B樣條曲線后，需對非主導(dǎo)點的擬合誤差進行檢測。由于此過程需要在線完成，而傳統(tǒng)牛頓迭代法的誤差計算時間較長，算法效率較低，因此為提高算法的效率和減少誤差計算時間，本文利用文獻［15］中提出的輪廓誤差跟隨檢測方法，計算非主導(dǎo)點到擬合曲線的誤差。實驗結(jié)果表明：相比牛頓迭代法，輪廓誤差跟隨法對同樣數(shù)量點到復(fù)雜曲線的誤差檢測時間更短，且其計算所需的曲率值可提前算出，并用于后續(xù)B樣條曲線插補速度規(guī)劃預(yù)處理時速度突變點的檢測，從而提高了算法的效率。

圖8中，C(ui)為非主導(dǎo)點Qi在曲線上的估算投影點，C(ui+δ)為非主導(dǎo)點Qi在曲線上的實際投影點，δ為Qi對應(yīng)的ui估算誤差，e(Qi)為非主導(dǎo)點到曲線的實際誤差。根據(jù)文獻［15］，e(Qi)的二階泰勒展開表達式為

式中，K為B樣條曲線在C(ui)處的曲率；dt為矢量D在矢量T上的投影長度，o3(dt)表示截斷誤差；dc為矢量D在矢量C上的投影長度；T為曲線在C(ui)處的切向矢量；N為曲線在C(ui)處的單位主法矢量；C為D在C(ui)處主平面上的投影，且T、N、C均為單位向量。

圖8　誤差跟隨法計算B樣條擬合誤差Fig.8　B-spline fitting error calculation by error following method

由圖8可知，為了提高誤差跟隨法的誤差檢測精度，必須合理估算非主導(dǎo)點Qi對應(yīng)的參數(shù)值ui。通常，對主導(dǎo)點進行B樣條插值擬合時，其節(jié)點向量是通過對主導(dǎo)點進行弦長參數(shù)化得到，而非主導(dǎo)點的參數(shù)值可通過再次對兩主導(dǎo)點間的數(shù)據(jù)進行弦長參數(shù)化得到。為了統(tǒng)一節(jié)點向量和數(shù)據(jù)點對應(yīng)參數(shù)值的計算方法，本文采用原始數(shù)據(jù)，對主導(dǎo)點進行B樣條擬合時的節(jié)點向量進行合理規(guī)劃。設(shè)主導(dǎo)點Dj在原始數(shù)據(jù)點中的下標(biāo)存儲在數(shù)組mark[N]中，其中N為原始數(shù)據(jù)點的個數(shù)，則進行B樣條擬合時，主導(dǎo)點處的參數(shù)值計算方法如下：

4　算法流程與仿真分析

4.1　算法流程

為了驗證本文所提方法及算法的正確性，以Visual Studio 2008為平臺，編寫了主導(dǎo)點（數(shù)據(jù)預(yù)處理）選取、主導(dǎo)點B樣條插值擬合及非主導(dǎo)點誤差計算的C++源程序，其算法流程見圖9。其中，數(shù)據(jù)預(yù)處理階段的誤差檢測容差e1設(shè)置較大，本文取e1=0.1 mm；實時處理階段的誤差結(jié)果一定小于e1，其檢測容差e2直接設(shè)置為最終的擬合誤差，本文取e2=0.03 mm，且在線處理時只需1次誤差檢測循環(huán)，并將超過誤差的非主導(dǎo)點增加為新的主導(dǎo)點，對更新后的主導(dǎo)點再次進行插值擬合生成B樣條曲線，此時的擬合誤差一定在e2范圍內(nèi)。

圖9　算法流程圖Fig.9　Algorithm flowchart圖12　主導(dǎo)點B樣條插值擬合Fig.12　Dominant points and B-spline curve

4.2　仿真分析

本文的仿真案例為1條連續(xù)的半蝴蝶狀曲線，通過UG NX8.0軟件生成加工路徑，離散點內(nèi)外公差絕對值均設(shè)置為0.02 mm。通過數(shù)據(jù)預(yù)處理后，篩選出138個主導(dǎo)點，原始數(shù)據(jù)點和主導(dǎo)點的分布見圖10。

圖10　主導(dǎo)點的組成部分Fig.10　Components of dominant points

主導(dǎo)點主要由離散數(shù)據(jù)點的曲率極大值點、曲率閾值點、曲線拐點、長度突變點以及誤差超限點組成。由圖10中的3處尖角放大圖可知，通過比較當(dāng)前曲率值點與前后5個連續(xù)相鄰曲率值點確定曲率極大值點，能很好地消除數(shù)據(jù)波動造成的偽曲率極大值點，但會漏選曲率閾值點處的曲率極大值點（圖11），不過，這些點已通過曲率閾值點進行補選。對于除去特殊形狀標(biāo)記覆蓋的叉形標(biāo)記主導(dǎo)點，其生成方法是，對非主導(dǎo)點進行長度均分策略處理，將長度突變點記為主導(dǎo)點，以及將Bezier曲線擬合后的誤差超限點補選為主導(dǎo)點，此過程通過迭代計算完成，且每次迭代循環(huán)時都要重新進行長度均分處理及Bezier曲線擬合。

圖11　2點與5點比較法生成的曲率極大值點Fig.11　Local curvature max points with 2 and 5 points method

對主導(dǎo)點進行B樣條插值擬合后生成的曲線見圖12。

圖9　算法流程圖Fig.9　Algorithm flowchart圖12　主導(dǎo)點B樣條插值擬合Fig.12　Dominant points and B-spline curve

將主導(dǎo)點插值擬合成B樣條曲線后，需要對非主導(dǎo)點進行誤差檢測，本文采用誤差跟隨法計算擬合誤差，以保證此過程的快速和高效性，并將其結(jié)果與牛頓迭代法的計算值相比較。結(jié)果表明：兩種方法的方法誤差小于0.01 mm。圖13為兩種計算方法的誤差計算值對比圖。

圖13　兩種方法的誤差計算值對比圖Fig.13　Comparison of error calculation values of the two methods

測試結(jié)果表明：對133個非主導(dǎo)點進行誤差計算，通過牛頓迭代法（其距離檢測精度為0.001 mm）需要循環(huán)計算275次，共耗時8.6 ms；而誤差跟隨法只需要循環(huán)計算133次，共耗時5.537 ms，誤差跟隨法節(jié)省時間約35.6%。對138個主導(dǎo)點進行1次B樣條插值擬合的時間為647.435 ms。由于在線處理最多只需進行2次B樣條插值擬合，因此，相比傳統(tǒng)的插值擬合方法，在允許時間內(nèi)，誤差跟隨法可一次性在線擬合更多的數(shù)據(jù)點，從而提高了擬合效率。

5　結(jié)論

（1）通過仿真實驗可知，本文提出的基于主導(dǎo)點的B樣條插值擬合方法可實現(xiàn)對復(fù)雜數(shù)控加工刀具軌跡（即G01代碼）的平滑壓縮，且利用輪廓跟隨誤差法的二階泰勒展開式計算擬合誤差，可在允許的誤差范圍內(nèi)提高計算效率。

（2）由于所提方法的數(shù)據(jù)預(yù)處理過程可在離線階段完成，而在線處理階段只包含1次誤差檢測循環(huán)和最多2次B樣條插值擬合，因此根據(jù)數(shù)控系統(tǒng)的處理能力，若合理調(diào)整處理的數(shù)據(jù)量，此方法可實現(xiàn)刀具軌跡的在線甚至實時平滑壓縮計算。

（3）為了驗證理論的正確性，本文選取了形狀較復(fù)雜的半蝴蝶狀曲線作為仿真實例，實驗證明，將原始的271個離散數(shù)據(jù)點壓縮成以B樣條曲線控制點形式存儲的138個數(shù)據(jù)點，其壓縮效率可提高近2倍。若曲線形狀不復(fù)雜，可在數(shù)據(jù)預(yù)處理階段適當(dāng)增大長度突變閾值，從而減少主導(dǎo)點數(shù)量，其壓縮效率將會進一步提高。