鄧華波，王瑋琪，梁務(wù)宇，劉思遠(yuǎn)

(沈陽(yáng)建筑大學(xué)機(jī)械工程學(xué)院，遼寧 沈陽(yáng) 110168)

航空、航天、航海領(lǐng)域中所用的葉輪葉片的加工要求加工精度高、加工表面光順等特點(diǎn)以實(shí)現(xiàn)設(shè)備行駛過程中的噪聲減小[1-2]，但由于此類零件形狀通過空氣動(dòng)力學(xué)等復(fù)雜公式計(jì)算得到，零件表面極為復(fù)雜，一旦其數(shù)控加工刀具路徑規(guī)劃不合理，將導(dǎo)致機(jī)床加工出的零件誤差較大、表面不光滑。因此，很多學(xué)者在上位機(jī)進(jìn)行數(shù)控加工軌跡的規(guī)劃及刀具速度、刀軸矢量的優(yōu)化[3-4]，然后后處理生成G代碼或采用NURBS插補(bǔ)器加工，進(jìn)而解決加工時(shí)的加工精度和表面光順問題。也有學(xué)者通過對(duì)曲面擬合的方式求解復(fù)雜曲面數(shù)控加工的刀具路徑規(guī)劃問題，在曲面造型中擬合精度要求較高[5-6]，精加工葉片時(shí)生成的刀具路徑是的空間螺旋曲線，加工時(shí)無退刀、無抬刀、省時(shí)[7]。以往學(xué)者大多在平面上進(jìn)行曲線擬合[8-10]，而空間復(fù)雜螺旋曲線投影至平面會(huì)產(chǎn)生變形導(dǎo)致擬合結(jié)果較差，從而影響刀具路徑規(guī)劃效果。故筆者在數(shù)控加工路徑規(guī)劃中直接對(duì)復(fù)雜空間螺旋曲線進(jìn)行三維空間擬合。

空間曲線擬合方法中一般通過最小二乘法[11-13]選擇部分?jǐn)?shù)據(jù)點(diǎn)經(jīng)過插值反求得到擬合曲線并計(jì)算擬合誤差，若擬合精度沒達(dá)到要求，則增加數(shù)據(jù)點(diǎn)的數(shù)量重新插值并計(jì)算擬合誤差，直到達(dá)到擬合誤差要求為止。很多學(xué)者采用遺傳算法、粒子群算法等對(duì)曲線最小二乘法逼近擬合進(jìn)行優(yōu)化[14-15]，但在空間中，逼近算法容易陷入局部最優(yōu)解，導(dǎo)致運(yùn)算量龐大[16-17]。S.Tulsyan等[18]提出了一種拼接方法，將相鄰的線性導(dǎo)軌段之間插值同步的五次B樣條用來平滑刀軸矢量位置和方向。該方法中，公式參數(shù)的連續(xù)性對(duì)曲線的連續(xù)性有很大的影響，拼接后需要降階，并且五次B樣條的公式過于復(fù)雜，計(jì)算基函數(shù)時(shí)計(jì)算量大，且B樣條在描述非拋物線時(shí)會(huì)產(chǎn)生一定的誤差。SHEN Liyong等[19]給出了一種利用三次Bezier樣條曲線在給定參數(shù)區(qū)間逼近的方法，由于Bezier樣條靈活特性，改動(dòng)一個(gè)控制點(diǎn)，全部控制點(diǎn)會(huì)發(fā)生變化，會(huì)導(dǎo)致在誤差計(jì)算時(shí)的計(jì)算量大，效率低。NURBS技術(shù)是目前國(guó)際標(biāo)準(zhǔn)組織(ISO)規(guī)定的表示自由型曲線曲面的唯一國(guó)際標(biāo)準(zhǔn)。工程中一般使用C2連續(xù)的三次NURBS曲線進(jìn)行擬合，擬合結(jié)果光順、計(jì)算量小、計(jì)算快、精度高。但C2連續(xù)的三次樣條曲線在光順時(shí)使一些數(shù)據(jù)點(diǎn)程度不同的偏離原始位置，有時(shí)會(huì)出現(xiàn)超出規(guī)定的誤差[20]。LIU Jiyuan[21]采用調(diào)整權(quán)因子的方式來改變局部曲線形狀，但在空間復(fù)雜螺旋曲線插值擬合過程中誤差較大點(diǎn)處調(diào)整的效果不明顯，并且權(quán)因子的調(diào)整要通過計(jì)算擬合曲線與原始數(shù)據(jù)點(diǎn)的最短距離，計(jì)算量龐大，且權(quán)因子的調(diào)節(jié)規(guī)律性不強(qiáng)。全局插值法雖然使擬合后的曲線盡可能穿過原始數(shù)據(jù)點(diǎn)，但仍有部分?jǐn)M合誤差較大。由于全局逼近或其他方法擬合后的曲線穿過原始數(shù)據(jù)點(diǎn)的數(shù)量不如全局插值擬合方法多，但全局插值在部分點(diǎn)的擬合誤差較大無法達(dá)到工況要求的擬合精度。由于NURBS曲線的局部可調(diào)性，筆者提出一種基于非均勻有理B樣條擬合的空間螺旋刀具路徑整體插值逼近方法，通過插入虛擬數(shù)據(jù)點(diǎn)，并計(jì)算插入數(shù)據(jù)點(diǎn)后的曲線扎人合誤差，擬合后的曲線僅用一個(gè)公式對(duì)空間復(fù)雜螺旋軌跡進(jìn)行表達(dá)，為優(yōu)化刀具的進(jìn)給速度或刀軸矢量提供數(shù)學(xué)表達(dá)式，擬合誤差能達(dá)到精加工的要求、誤差變化幅度明顯減小情況下，計(jì)算量較小。

1 空間螺旋刀具路徑的三次NURBS曲線擬合

1.1 NURBS曲線表示與節(jié)點(diǎn)計(jì)算

三維空間的NURBS曲線定義為對(duì)于給定的控制點(diǎn)di=[xiyizi]及其權(quán)因子ωi(i=0,1，…,n)，其帶權(quán)控制點(diǎn)Di=[ωidiωi]=[ωixiωiyiωiziωi],定義一條四維的k次非有理B樣條曲線P(u)，取它在第四坐標(biāo)ω=1的超平面上的中心投影，能夠得到三維空間的一條k次NURBS曲線p(u)。步驟如下：

(1)確定所給控制頂點(diǎn)di=[xiyizi]的帶權(quán)控制點(diǎn)Di=[ωidiωi]=[ωixiωiyi

ωiziωi]。

(2)定義一條四維k次非有理B樣條曲線。

(1)

(3)將其投影到ω=1超平面上并獲得一條k次NURBS曲線p(u)。

(2)

其中，Ni,k(u)為由節(jié)點(diǎn)矢量。

U=[u0,u1,…,un+k,un+k+1],為按de Boor-Cox遞推公式?jīng)Q定的k次規(guī)范B樣條基函數(shù)。

當(dāng)k=0時(shí)，

(3)

當(dāng)k∈[1,+∞)時(shí)，

(4)

通過對(duì)型值點(diǎn)的參數(shù)化得到節(jié)點(diǎn)矢量。NURBS曲線一般采用弦長(zhǎng)參數(shù)化法效果較好，其中福利參數(shù)化法下生成的插值曲線顯現(xiàn)出最好的光順性，擬合誤差較小。其計(jì)算公式為

(5)

規(guī)范化使節(jié)點(diǎn)矢量在[0，1]內(nèi)，進(jìn)而得到(u0=…=uk=0,uk+1,…,un+k=…=un+2k=1)節(jié)點(diǎn)矢量。

1.2 邊界條件定義

由于節(jié)點(diǎn)矢量中存在重復(fù)的節(jié)點(diǎn)[19]，會(huì)導(dǎo)致首(末)頂點(diǎn)和首(末)數(shù)據(jù)點(diǎn)重合，首末端點(diǎn)處的條件不足。因此，需要增加由邊界條件確定的附加方程求出帶權(quán)控制頂點(diǎn)Di，然后將帶權(quán)控制頂點(diǎn)在ω=1的超平面上進(jìn)行投影，得到待求控制頂點(diǎn)di。當(dāng)ωi(i=0,1，…,n)都為1時(shí)，非均勻有理B樣條曲線就轉(zhuǎn)化為非均勻B樣條曲線。因此只要定義非均勻B樣條的邊界條件亦可求出非均勻有理B樣條曲線的邊界條件方程。采用連續(xù)的C2三次NURBS曲線對(duì)空間螺旋曲線進(jìn)行擬合，要求曲率均勻變化，首末端點(diǎn)采用拋物線條件，有二階導(dǎo)且為實(shí)數(shù)。

首端點(diǎn)拋物條件：

(6)

末端點(diǎn)拋物條件：

(7)

1.3 控制頂點(diǎn)求解

通過式(1)對(duì)3次非有理B樣條曲線的定義，可知當(dāng)ωi都為1時(shí)，NURBS曲線即為3次非均勻B樣條曲線。

根據(jù)反算C2連續(xù)的三次B樣條插值曲線的控制頂點(diǎn)矩陣方程：

(8)

曲線首端點(diǎn)與第一個(gè)型值點(diǎn)重合，末端點(diǎn)與最后一個(gè)型值點(diǎn)重合q0=d0，qn=dn+2。矩陣方程(8)中首行的元素b1、c1、a1與e1表示首端點(diǎn)邊界條件。末行的bn+1、cn+1、an+1與en+1表示末端點(diǎn)邊界條件。

采用追趕法解方程(8)，求出全部控制頂點(diǎn)di，對(duì)于非均勻有理B樣條，只需將方程(8)中的di改成帶權(quán)的控制頂點(diǎn)Di=[ωidiωi]，再投影到ω=1超平面上，即用ωidi除以ωi就得到控制頂點(diǎn)di。

1.4 擬合誤差計(jì)算

原數(shù)據(jù)點(diǎn)到空間曲線的最短距離為曲線擬合誤差，通過設(shè)定允許誤差界ε，若得到的擬合曲線最大誤差在允許的誤差界ε內(nèi)，則可作為最終曲線使用。求解曲線的最大擬合誤差，一般通過曲線p(u)到型值點(diǎn)qi的距離，取得qi到曲線點(diǎn)最短距離。豪斯多夫距離是曲線擬合偏差檢查通常檢查最大范數(shù)距離，其公式：

(9)

其中，qi為原始型值點(diǎn)；p(u)為擬合得到的曲線上的點(diǎn)。實(shí)際擬合過程中，型值點(diǎn)qi到擬合曲線p(u)的最短距離可通過二分法快速得到一個(gè)允許誤差范圍的近似解。

2 插值逼近擬合方法的求解

2.1 擬合方法的求解流程

筆者綜合運(yùn)用插值、逼近算法實(shí)現(xiàn)曲線擬合算式的求解，插值逼近擬合方法的算法流程如圖1所示。通過導(dǎo)入數(shù)據(jù)點(diǎn)即原始數(shù)據(jù)點(diǎn)，選擇曲線曲率變化較大處的特征點(diǎn)qm。確定特征的位置點(diǎn)qm和要插入的虛擬數(shù)據(jù)點(diǎn)的數(shù)量和位置。插入虛擬數(shù)據(jù)點(diǎn)后，對(duì)所有點(diǎn)進(jìn)行NURBS全局插值，對(duì)u在區(qū)間[m+2-r,m+4+r]的曲線采用式(9)進(jìn)行擬合誤差計(jì)算，采用二分法計(jì)算。當(dāng)計(jì)算得到的擬合誤差大于指定誤差時(shí)，停止對(duì)剩余曲線誤差的計(jì)算，進(jìn)入下一步增加虛擬數(shù)據(jù)點(diǎn)的插入數(shù)量。

圖1 插值逼近方法求解流程Fig.1 The flow of interpolation approximation

插入虛擬數(shù)據(jù)點(diǎn)后，再對(duì)所有點(diǎn)進(jìn)行NURBS全局插值，對(duì)u在區(qū)間(其中r為特征點(diǎn)兩側(cè)各選取的原始型值點(diǎn)數(shù)量)，曲線采用式(8)進(jìn)行擬合誤差計(jì)算，采用二分法或牛頓迭代法計(jì)算。當(dāng)計(jì)算得到的擬合誤差大于指定誤差時(shí)，停止對(duì)剩余曲線誤差的計(jì)算，繼續(xù)增加虛擬數(shù)據(jù)點(diǎn)的插入數(shù)量g。然后重復(fù)計(jì)算過程，直至計(jì)算得到的擬合誤差全者小于指定誤差為止。

2.2 尋找特征點(diǎn)

在復(fù)雜曲面刀具路徑生成過程中，常采用導(dǎo)入原始型值點(diǎn)，同時(shí)考慮曲線的撓率和曲率。為避免計(jì)算時(shí)受空間維度限制，曲線可看成由連續(xù)原始數(shù)據(jù)點(diǎn)qi、qi-1、qi+1所在平面上的弧線組成，根據(jù)圓弧估算法求離散曲率κi的公式[16]：

(10)

其中,夾角∠qi-1qiqi+1的區(qū)間為[90°,180°]，|qi-1qi+1|為弧線的弦長(zhǎng)。

正弦函數(shù)sin(∠qi-1qiqi+1)在[90°,180°]內(nèi)單調(diào)遞減，為方便計(jì)算可視為與夾角角度∠qi-1qiqi+1成反比。因此，結(jié)合式(10)可知曲率半徑、夾角角度∠qi-1qiqi+1和弦長(zhǎng)|qi-1qi+1|成正比。通過觀察發(fā)現(xiàn)，在刀具路徑彎曲程度改變較大時(shí)，弦長(zhǎng)和夾角通常同時(shí)增大或減小。因此，在尋找特征點(diǎn)的過程中，僅需尋找夾角角度變化和弦長(zhǎng)這兩個(gè)自變量變化均較大的數(shù)據(jù)點(diǎn)即可，即計(jì)算弦長(zhǎng)之差的變化量‖qi+1qi+3|-|qiqi+2‖與夾角角度之差變化量|∠qi+1qi+2qi+3-∠qiqi+1qi+2|的乘積Am+1。

Am+1=‖qi+1qi+3|-|qiqi+2‖×|∠qi+1qi+2qi+3-∠qiqi+1qi+2|.

(11)

其中，|∠qi+1qi+2qi+3-∠qiqi+1qi+2|為qi、qi+2小線段和qi+1、qi+3小線段的夾角與qi、qi+1小線段和qi+1、qi+1小線段夾角之差的絕對(duì)值，單位為rad;‖qi+1qi+3|-|qiqi+2‖為弦長(zhǎng)之差的絕對(duì)值，mm；m的取值為1，2，…，n-2；|qi+1qi+3|為qi+3到qi+1弧線的弦長(zhǎng)；|qiqi+2|為qi+2到qi弧線的弦長(zhǎng)；乘積A可看成以兩個(gè)夾角角度之差為弧度并且以弦長(zhǎng)之差為半徑的弧長(zhǎng)，其單位為mm?？筛鶕?jù)乘積A簡(jiǎn)單判斷曲線曲率半徑變化大小，A越大，曲線彎曲程度變化越大，進(jìn)而導(dǎo)致曲線擬合誤差越大。計(jì)算A的平均值，將A中大于平均值所對(duì)應(yīng)的m值找出，如圖2所示，其中m為初始數(shù)據(jù)點(diǎn)序列號(hào)。

圖2 A與m值的關(guān)系圖Fig.2 The relationship between A and m

2.3 虛擬數(shù)據(jù)點(diǎn)插入

將虛擬型值點(diǎn)添加在相鄰兩個(gè)原始型值點(diǎn)之間的直線上，在特征點(diǎn)qm兩側(cè)各選取原始型值點(diǎn)數(shù)為r。在特征點(diǎn)qm前后先各加0個(gè)虛擬型值點(diǎn)。通過計(jì)算分析當(dāng)特征點(diǎn)兩側(cè)所選擇的原始型值點(diǎn)的數(shù)量小于等于3時(shí)，插入虛擬型值點(diǎn)的位置在相鄰兩原始型值點(diǎn)線段上，距離前一點(diǎn)的5/8處擬合誤差較小，如圖3中qm-1與qm之間線段上距離qm點(diǎn)5/8處的點(diǎn)J。當(dāng)添加虛擬型值點(diǎn)前后各添加點(diǎn)數(shù)大于3時(shí)，根據(jù)計(jì)算在相鄰兩原始型值點(diǎn)線段中點(diǎn)位置添加虛擬數(shù)據(jù)點(diǎn)，擬合后的最大誤差相較其他位置最小。

圖3 虛擬數(shù)據(jù)點(diǎn)添加位置Fig.3 The add location for virtual data point

若擬合誤差小于指定誤差，則確定NURBS曲線公式符合工況要求。若大于指定誤差，則將添加型值點(diǎn)的個(gè)數(shù)前后各增加1，在原始型值點(diǎn)qm前后先各插入1個(gè)虛擬型值點(diǎn)，圖3中J點(diǎn)和K點(diǎn)。若大于指定誤差，插入圖3中的H、L點(diǎn)。本次添加的虛擬型值點(diǎn)位置在圖3中的H、J、K、L點(diǎn)，擬合后計(jì)算與原始型值點(diǎn)的最短距離中的最大值與指定誤差比較，若大于指定誤差則繼續(xù)增加型值點(diǎn)個(gè)數(shù)，即圖3中的G、H、J、K、L、M。當(dāng)添加虛擬型值點(diǎn)前后各添加點(diǎn)數(shù)大于3時(shí)，根據(jù)計(jì)算在相鄰兩點(diǎn)線段中點(diǎn)位置添加型值點(diǎn)，擬合后的最大誤差相較其他位置最小。因此，當(dāng)添加虛擬型值點(diǎn)前后各添加點(diǎn)數(shù)大于3時(shí)，添加型值點(diǎn)的位置變?yōu)橄噜弮牲c(diǎn)線段中點(diǎn)位置。通過不斷添加虛擬型值點(diǎn)，并計(jì)算相應(yīng)的擬合誤差。以上為插入虛擬型值點(diǎn)總數(shù)g小于n時(shí)的情況。

當(dāng)曲線最大擬合誤差大于指定誤差ε且插入虛擬型值點(diǎn)總數(shù)g等于n×t時(shí)，h值增加1，當(dāng)本次擬合誤差大于指定誤差時(shí)，執(zhí)行步數(shù)為h+1步，重復(fù)上述過程，并以此類推。

2.4 擬合誤差計(jì)算

根據(jù)式(9)得到擬合誤差計(jì)算公式:

(12)

其中，ε為指定誤差。通過二分法或牛頓迭代法快速計(jì)算插入虛擬型值點(diǎn)部分的曲線u的范圍為[m+2-r，m+4+r]上的擬合點(diǎn)與原始型值點(diǎn)之間的最短直線距離得到一個(gè)允許誤差范圍的近似解，m按從小到大的順序取，只要有一個(gè)m所在區(qū)間的擬合誤差值大于指定誤差C2就停止擬合誤差的計(jì)算，進(jìn)入下一步。若區(qū)間上所有曲線最大擬合誤差小于或等于指定誤差ε，則得到所求NURBS曲線公式。

3 局部逼近擬合方法算法算例

3.1 插值逼近方法算例

算例采用如圖4所示的葉片模型，其工件坐標(biāo)系下刀具路徑的刀位點(diǎn)為筆者算例的數(shù)據(jù)點(diǎn)。葉輪毛坯直徑為100 mm，選用球頭銑刀，采用葉片精加工驅(qū)動(dòng)方法，切削模式選用螺旋銑削模式。行距參數(shù)采用等高殘留法并設(shè)置0.2 mm的殘留高度，設(shè)置加工時(shí)內(nèi)外公差均為0.02 mm，通過UG軟件實(shí)現(xiàn)刀具路徑的初步規(guī)劃如圖4所示。

圖4 刀具路徑Fig.4 The tool path

將工件坐標(biāo)系下的刀具路徑導(dǎo)出，文件格式為CLS。刀具路徑文件數(shù)據(jù)中，前三列為刀位點(diǎn)數(shù)據(jù)，分別對(duì)應(yīng)工件坐標(biāo)系X、Y、Z軸。算例生成的刀具路徑共2 637個(gè)數(shù)據(jù)點(diǎn)。精加工葉輪葉片允許誤差設(shè)為0.02 mm，擬合誤差為加工誤差的1/5時(shí)對(duì)葉輪葉片加工影響較小，故算例設(shè)定NURBS擬合最大誤差為 4 μm。

導(dǎo)入數(shù)據(jù)點(diǎn)，通過式(11)計(jì)算乘積A值，m為刀具路徑文件中的刀位點(diǎn)序列號(hào)。A與m關(guān)系如圖5所示。計(jì)算乘積A的平均值，并找出大于乘積A平均值的對(duì)應(yīng)點(diǎn)分別為第m個(gè)。

圖5 算例中A與m值的關(guān)系圖Fig.5 The relationship between A and m

超過A平均值所對(duì)應(yīng)的m值有343個(gè)。按局部逼近擬合方法步驟在qm前后各插入0個(gè)虛擬型值點(diǎn)，此時(shí)為全局插值擬合，此時(shí)計(jì)算通過式(9)計(jì)算曲線擬合誤差，可得到最大擬合誤差為0.025 0 mm，擬合誤差大于4 μm，則根據(jù)插值逼近法流程添加插入虛擬型值點(diǎn)的個(gè)數(shù)并根據(jù)式(11)計(jì)算相應(yīng)曲線段的擬合誤差。當(dāng)r為3時(shí)得到刀具路徑曲線擬合誤差，如圖6所示。

圖6 插值逼近法擬合誤差Fig.6 The fitting error of interpolation approximation

3.2 兩種方法擬合誤差參數(shù)對(duì)比

通過全局插值方法擬合曲線，并采用式(9)計(jì)算得到2 637個(gè)數(shù)據(jù)點(diǎn)擬合后誤差，全局插值方法擬合誤差如圖7所示。

圖7 全局插值擬合方法擬合誤差Fig.7 The fitting error of global interpolation fitting method

將插值逼近的擬合誤差結(jié)果如圖6所示，對(duì)比兩種擬合方法的最大誤差、控制頂點(diǎn)以及方差(見表1)。

表1 兩種方法擬合2 637個(gè)數(shù)據(jù)點(diǎn)后曲線參數(shù)對(duì)比Table 1 Comparison of curve parameters after fitting 2 637 data points by two methods

對(duì)兩種方法擬合后的誤差、控制頂點(diǎn)數(shù)以及方差進(jìn)行比較，在提供數(shù)據(jù)點(diǎn)相同的情況下，筆者方法的最大誤差和方差遠(yuǎn)小于全局插值方法。當(dāng)要求最大擬合誤差達(dá)到4 μm以下時(shí)，采用最小二乘法在擬合時(shí)需要重復(fù)計(jì)算整段曲線的擬合誤差，而筆者方法僅需重復(fù)計(jì)算部分曲線段的擬合誤差，且需要計(jì)算擬合誤差的曲線段長(zhǎng)度小于整體曲線段的36.7%，并且擬合得到的曲線穿過大部分的型值點(diǎn)，在C2參數(shù)連續(xù)的情況下，能夠更準(zhǔn)確的表達(dá)刀具路徑的特點(diǎn)。

4 結(jié) 論

(1)采用筆者方法對(duì)空間復(fù)雜螺旋曲線擬合，得到了擬合軌跡既符合C2連續(xù)的三次NURBS曲線光順的條件，最大擬合誤差也達(dá)到了精密加工的工況要求，并且僅用一個(gè)公式表達(dá)出空間螺旋刀具路徑。

(2)三次NURBS曲線連續(xù)性好，擬合后曲線最大擬合誤差比采用原始全局插值方法擬合后的誤差要低，且擬合誤差的波動(dòng)范圍減小。

(3)與最小二乘法對(duì)比，筆者方法在計(jì)算擬合后曲線的擬合誤差方面的計(jì)算量小，得出曲線穿過大部分的數(shù)據(jù)點(diǎn)，能夠更準(zhǔn)確地表達(dá)刀具路徑的形狀特點(diǎn)，在零件加工時(shí)保證刀具路徑參數(shù)的準(zhǔn)確性。