適用于高速高精加工的橢圓弧平滑壓縮插補算法

李 浩，吳文江，陳淥萍，韓文業，郭 安

1(中國科學院大學，北京 100049) 2(中國科學院沈陽計算技術研究所 高檔數控國家工程研究中心， 沈陽 110168) 3(沈陽高精數控智能技術股份有限公司， 沈陽 110168) 4(工業和信息化部計算機與微電子發展研究中心 中國軟件評測中心，北京 100048)

1 引 言

橢圓是圓錐曲線中的一種，其擁有幾何形式，能夠被精確地表示.由于橢圓弧長的計算涉及到橢圓積分問題，且其原函數不能用有限形式表示，所以無法直接對橢圓進行精確插補.又因其具有較強的工程實用性，國內外的學者對此進行了廣泛的研究.

蘇[1]等提出了映射法，將橢圓映射到某一平面，使其投影為圓，通過對圓進行插補，然后利用坐標變換得到橢圓上的插補點.胡[2]等推導出，在弧長增量相對橢圓周長非常小的情況下，下一插補點與當前插補點和離心角增量的關系.上述方法雖簡單可行，但都存在輪廓誤差和速度波動無法控制的問題.于[3]等提出了等誤差方法，雖然提高了橢圓的加工精度，但仍無法有效地控制速度波動.因此大量學者對數值方法進行了研究，通過較為精確地推導出插補點與弧長的關系，來有效地降低速度波動，如胡[4]提出的冪級數展開法、Toshio[5]提出的半參數和雙參數轉換法、John[6]提出的經驗參數迭代法和Toshio[7]提出的逆橢圓積分法，這些方法因涉及到多次迭代與復雜數值計算，實時性較差，所以均無法應用于數控系統.

隨著設計和制造技術的發展，越來越多的人使用計算機輔助設計(CAD)系統來進行復雜零部件的設計，然后利用計算機輔助制造(CAM)系統將CAD設計的自由曲線或曲面，在設定的最大輪廓誤差范圍內，用微小直線段去逼近，從而生成由大量指令點構成的數控加工程序[8].

通過對微小直線段進行數控加工[9,10]，可以不用考慮復雜的橢圓弧長計算問題.但由于加速度和加加速度的頻繁變化，會引起機床的震動，降低加工的效率和質量.而現在對于微小線段高速加工的研究主要分為兩種方法.一種是在相鄰微小線段的拐角處插入過渡樣條曲線[11].如Zhang等[12]提出通過在拐角處插入三次參數樣條曲線的方法，提高拐角處的速度，從而提高加工效率.但由于插補點在直線段和過渡樣條曲線段上的循環出現，使得插補步長不一致，從而導致加工速度的波動，如果指令點越密集，速度的波動就會越頻繁.另一種則是通過插值或者逼近的方式將離散的指令點擬合成平滑的加工路徑[13].例如Zhang[14]提出利用二次Bézier曲線將連續微小線段擬合成光滑的樣條曲線，通過樣條插補來實現對自由曲線的高速高精加工.雖然擬合方式能較好地逼近原設計曲線，但由于擬合曲線較復雜，無法精確計算插補步長對應的插補參數，導致加工速度波動，降低了加工精度.

橢圓弧是一種不能被直接精確快速插補的圖形.在離散小線段的形式下采用樣條插補方式，可以避免復雜的橢圓弧長計算，但會存在計算精度和效率的問題.所以，基于以上研究和分析，本文提出了一種適用于高速高精加工的橢圓弧平滑壓縮插補算法.

2 橢圓弧平滑壓縮插補算法

本文算法的基本思想是：從由離散小線段構成的加工路徑中識別出連續微小線段加工區域，然后尋找到型值點，通過控制擬合誤差，將指令點指定的折線加工路徑轉化為平滑的二次有理Bézier曲線加工路徑.然后，根據樣條曲線的特征識別出橢圓弧段，將有理形式轉換為幾何形式，同時合并相鄰的橢圓弧段.最后，通過在幾何形式的橢圓弧上進行插補，實現高速高精加工.

該算法采用預處理方式，其由7個部分組成：加工區域的識別、型值點的選取、型值點的擬合、橢圓弧的識別、幾何形式的轉換、橢圓弧段的合并和橢圓弧的插補，其中加工區域的識別主要是從加工路徑中尋找可以擬合的連續區；型值點的選取則是尋找連續區中的曲率極值點和拐點；型值點的擬合用于將連續區內的折線加工路徑轉化成光滑的樣條曲線；橢圓弧的識別用于識別樣條曲線中的橢圓弧段；幾何形式的轉換是將橢圓弧段的有理表示形式轉換為幾何表示形式；橢圓弧段的合并用于合并已識別出的相鄰橢圓弧段；橢圓弧的插補用于對幾何形式的橢圓弧進行插補計算.

圖1 連續微小線段驗證Fig.1 Continuous micro-line verification

2.1 加工區域的識別

平滑壓縮算法并非針對所有的G01段，只有連續的微小線段才適合于擬合為Bézier曲線.以雙弓高誤差為判斷依據，從由大量微小線段構成的加工路徑中識別出連續微小線段加工區域.如圖1所示Pi-1、Pi和Pi+1為順序的三個指令點l1、l2，為小線段的段長，θ為小線段間的拐角，雙弓高誤差可由式(1)和(2)得到.

δ1=R(1-cosφ1)

(1)

δ2=R(1-cosφ2)=R(1-cos(π-θ-φ1))

(2)

(3)

(4)

如果δ1或δ2大于設定的最大弓高誤差值δmax，則Pi為斷點；若兩個相鄰的斷點間存在指令點，那么兩個斷點連同它們之間的指令點就構成了一個連續微小線段加工區域.

2.2 型值點的選取

為了減少對連續微小線段加工區域的擬合次數、增加程序段的壓縮量，通過如下三個步驟來選取型值點.

1)將連續加工區域的開始點和結束點，即斷點，標記為型值點.

2)將局部曲率最大值點標記為型值點.

如圖1所示，三個指令點的坐標為Pi-1(xi-1，yi-1)、Pi(xi，yi)和Pi+1(xi+1，yi+1)，由式(5)可以計算出離散指令點Pi的曲率值[15]，其中ΔPi-1PiPi+1為帶符號的三角型面積.

(5)

(6)

假設kl為Pi左邊的局部曲率最小值，kr為Pi右邊的局部曲率最小值，如果滿足下面兩個條件，則Pi為局部曲率最大值點，如圖2所示.

1)|ki|>|kl|并且|ki|>|kr|

2)|ki|-|kl|≥δf或者|ki|-|kr|≥δf，δf為設定的最大曲率差值.

圖2 局部曲率最大值Fig.2 Local curvature maximum

3)將加工路徑彎曲方向改變的點，即拐點，標記為型值點.

利用第(2)步計算出的指令點曲率值進行判斷，如果ki-1ki>0并且kiki+1<0，則Pi為拐點.

2.3 型值點的擬合

假設連續微小線段加工區域中有n個型值點，為了達到平滑壓縮加工路徑的目的，使用n-1段二次有理Bézier曲線對其進行擬合.

二次有理Bézier曲線的定義為

(7)

其中，Pi=(xi，yi，zi)為控制點，wi是標量，為權因子，u為參變量，且u∈[0，1].P0P1P2構成了曲線的控制多邊形，如圖3所示.

圖3 二次有理Bézier曲線Fig.3 Quadric rational Bézier curve

本文采用標準型的二次有理Bézier曲線，取w0=w1=1，w1>0，標準型的二次有理Bézier曲線為式(8)所示

(8)

由圖3可知，當給定首末端點P0和P2，以及這兩點處的切線方向T0和T2時，可以很容易地通過直線[P0，T0]和[P2，T2]的交點得到P1，再給定一點P就能唯一確定該曲線，因而確定了w1.

將所求的曲線看成是由點P0，P1和P2確定的拋物線的投影，P1為投影中心.如圖3所示，將直線段[P0，P2]投影到要求的曲線上，則點P和Q為相對應的投影點.令w1=0，得到直線段L(u)=[P0，P2]，即

(9)

L(u)是點P0和P2的凸組合，因此|P0Q|和|QP2|的比值為u2：(1-u)2，從而推出

(10)

(11)

將u和P帶入式(8)，得到w1，從而得到所求曲線.

(12)

假定型值點Qi與Qj之間由指令點Qi+1，Qi+2，…，Qj-2與Qj-1指定的折線加工路徑組成，且型值點切矢已知，則型值點Qi與Qj之間的擬合詳細步驟如下所述.

1)構造以P0=Qi，P1和P2=Qj為控制點的擬合曲線，使其分別插值于Qk(k=i+1，…，j-1).這將產生中間的權值wk及相應的肩點坐標sk=wk/(1+wk)；

2)通過對sk求平均值得到逼近曲線的肩點坐標，即

(13)

則中間的權值為w=s/(1-s)；

3)根據P0，P1，P2和w即可確定型值點Qi與Qj之間的擬合曲線.

2.3.1 型值點處切矢的計算

由于數控加工程序中并不提供指令點處的切矢量，但可以通過型值點以及它周圍的四個指令點來計算型值點處的切矢量.如圖4所示，Qi為型值點，Qi-2、Qi-1、Qi+1、Qi+2為它周圍的四個指令點，Qi的單位切矢表達式如下

圖4 型值點處切矢的計算Fig.4 Tangent vector of reference point

(14)

Vi=(1-αi)qi+αiqi+1

(15)

qi=Qi-Qi-1

(16)

(17)

由于計算公式中未使用參數u，所求的單位切矢只能被看作切矢的方向.將下式帶入式(14)～(17)，即可得到T0，T1和Tn-1，Tn.

q0=2q1-q2

(18)

q-1=2q0-q1

(19)

qn+2=2qn+1-qn

(20)

qn+1=2qn-qn-1

(21)

2.3.2 加工精度控制

雖然曲線段經過型值點Qi與Qj，但不能保證它到Qi與Qj之間的所有指令點的距離都滿足最大擬合誤差.因此，將Qi與Qj之間的所有指令點投影到擬合的曲線上，檢查它們的偏差是否都在允許的誤差范圍內，即小于等于系統設定的最大輪廓誤差δc，如果滿足誤差要求，進行下一段的曲線擬合；否則，將偏差最大的指令點設置為新的型值點，用新的型值點重新進行曲線擬合，再次檢驗擬合精度，重復這個過程，直至滿足誤差要求.

2.4 橢圓弧的識別

由于標準型的二次有理Bézier曲線只有一個權因子，其表達能力比非有理形式的Bézier曲線強，可以表達多種曲線，如圖5所示，當01時，曲線為雙曲線.

圖5 不同w1定義的曲線Fig.5 Curves defined by w1

假設Ci(u)為加工路徑中某段的標準型二次有理Bézier曲線，如果其權值，0

2.5 幾何形式的轉換

如果Ci(u)為加工路徑中某段的標準型二次有理Bézier曲線，且其對應于一個橢圓弧段，如圖6所示.

根據Ci(u)可得到該橢圓弧段的起點坐標p_start和終點坐標p_end.由Lee[16]的推導可知，該橢圓弧段的中心坐標為

圖6 橢圓弧段幾何信息的計算Fig.6 Calculation of geometric information of elliptical arcs

P1+ε(S+T)

(22)

S=P0-P1

(23)

T=P2-P1

(24)

(25)

(26)

其中，P0、P1和P2為Ci(u)的控制點，w1為權值.

該橢圓弧段的長、短半徑分別為

(27)

(28)

假設λ2>λ1>0，且是以下二次方程的根

2δλ2-(kη+4β)λ+2(k-1)=0

(29)

(30)

(31)

β=S·T

(32)

該橢圓弧段長軸上的兩個點為

q1=P1+(ε+r1x0)S+(ε+r1y0)T

(33)

q2=P1+(ε-r1x0)S+(ε-r1y0)T

(34)

根據q1、q2，可以得到長軸的斜率kl，和長軸與x軸正半軸的夾角d_kl；以長軸為x′軸，短軸為y′軸建立局部坐標系，利用q1、q2、中心坐標、起點坐標和終點坐標，求得相對于局部坐標系x′正半軸的起始角度d_start、終點角度d_end和橢圓弧的方向turn.

當連續微小線段加工區域完成橢圓弧的識別和幾何形式轉換后，會得到一個橢圓弧段數組Ellipse_Arcs[ ]，其中橢圓弧段的數據結構如下：

structEllipse_Arc_INF//橢圓弧段幾何信息

{

doublecenter[2]； //中心坐標

doubler_l； //長半徑

doubler_s； //短半徑

doublep_start[2]； //起點坐標

doublep_end[2]； //終點坐標

doubled_start； //起點角度

doubled_end； //終點角度

doubleturn； //橢圓弧的方向

intkl； //長軸的斜率

doubled_kl； //長軸與x正半軸的夾角

}；

2.6 橢圓弧的合并

為了提高橢圓弧加工的效率和精度，將連續微小線段加工區域內的橢圓弧段合并，如果相鄰的橢圓弧段具有相同的中心center、長半徑r_l、短半徑r_s、方向turn、長軸與x正半軸的夾角d_kl，并且上一橢圓弧段的終點坐標p_endi-1與下一橢圓弧段的起點坐標p_starti相等，則它們屬于同一橢圓弧上連續的橢圓弧段，通過修改上一橢圓弧段的終點角度d_endi-1和終點坐標p_endi-1將它們合并為一個橢圓弧段.重復上述步驟，直至橢圓弧段不能再合并.

2.7 橢圓弧的插補

2.7.1 橢圓弧長計算

圖7 區域及變量劃分Fig.7 Region and variable partitions

(35)

根據高斯勒讓德求積公式[17]，可將上式轉換為

(36)

(37)

其它象限x或y為變量的情況，其弧長的推導方法相同.

2.7.2 弧長與變量的關系

2)根據公式(36)計算y0、y1、y2和y3到橢圓弧的起點y0的弧長，s0、s1、s2、s3.

3)將(s0，y0)、(s1，y1)、(s2，y2)和(s3，y3)代入如下三次多項式中進行插值，求得多項式的系數，即得到弧長與變量的關系.

y=a0+a1s+a2s2+a3s3，s∈[s0，s3]

(38)

其它象限x或y為變量的情況，其弧長與變量關系的推導方法相同.

2.7.3 橢圓弧實時插補

為了保證算法的實時性，對上述過程采用預處理方式.當加工程序讀入數控系統后，對整個G01段進行處理，提前將加工路徑中所有橢圓弧段的弧長與其參數關系進行計算并保存下來.在預處理結束后，再進行加工，直接利用計算好的弧長與參數關系實現實時插補.即在插補過程中，只需要7次浮點運算和5次加法運算，就能計算出第i個插補周期插補點坐標(xi，yi)，其具體步驟如下所示.

圖8 橢圓弧平滑壓縮算法流程圖Fig.8 Flow chart of the proposed algorithm

假定插補周期為T，P0(xi-1，yi-1)和si-1分別為第i-1個插補周期的插補點坐標與弧長，第i個插補周期的進給速度為vi-1，則弧長增量為Δs=vi-1T，第i個插補周期的弧長為si=si-1+Δs，將si代入弧長與變量的關系中，得到對應變量xi或yi的值，再將其代入橢圓方程中，從而求得第i個插補周期的插補點坐標P(xi，yi).

2.8 算法流程

整個橢圓弧平滑壓縮插補算法采用預處理方式，其流程如圖8所示.

首先，讀入一系列指令點，通過連續微小線段判斷條件，識別出連續加工區域.在連續區域內，根據指令點的曲率，尋找型值點.然后，將相鄰型值點間折線加工路徑擬合為平滑的二次有理Bézier曲線加工路徑，并檢驗每個指令點的擬合誤差，將誤差最大的點添加為新的型值點，重新進行擬合；如果都滿足誤差要求，從擬合曲線中識別出橢圓弧，并將其轉換為幾何形式；最后，將轉換后的橢圓弧段與相鄰的橢圓弧段合并.如果還有未擬合的型值點，則選取下一個型值點繼續擬合，否則，計算橢圓弧段的弧長與變量的關系，然后判斷是否還有未處理的指令點，如果有，重復上述操作直至最后一個指令點.

3 試驗驗證

為了驗證該橢圓弧平滑壓縮插補算法的正確性和有效性，本文使用了由沈陽高精數控智能技術股份有限公司自主研發的開放式四軸可控三軸聯動數控系統GJ301，如圖9所示，該系統采用工業級PC硬件平臺，以及RTLinux實時操作系統.在驗證過程中，為了與本文所提出的算法形成對比，將通過在小線段拐角處插入五次樣條曲線的插補方法[11]稱為傳統算法.

圖9 GJ301數控系統Fig.9 GJ301 NC System

設定系統的插補周期為2ms，x、y軸的最大加速度為A= 2000mm/s2，進給速度為F= 3000mm/min，雙弓高誤差為δmax= 0.01mm，最大曲率差值為δf= 0.0005mm-1，輪廓誤差為δc= 0.01mm，試驗中輸入的初始加工路徑如圖10所示，該加工路徑是按照0.001mm的精度將三個順序相連的橢圓弧離散成211段小線段所生成的.

圖10 加工路徑Fig.10 Machining trajectory

由圖11和圖12的對比可以看出，傳統算法在小線段的每個拐角處插入過渡曲線，平滑加工路徑，提高加工速度，而本文的算法可以識別出所有的橢圓弧段并以三條曲線表示，加工路徑的壓縮比達到了70.3：1，同時具有較好的光滑性.

圖11 傳統算法拐角插入過渡曲線Fig．11 Inserttransitioncurveatthecorner圖12 本文算法橢圓弧識別Fig．12 Recognizeellipsearc

加工過程中的速度曲線如圖13、圖14所示，傳統算法和本文算法的速度曲線形狀相同，但傳統算法存在拐角，且加工時間分別為3.6s和3.464s.將勻速段放大后發現，傳統算法存在頻繁的波動，而本文算法的波動不明顯.

圖13 傳統算法加工的速度曲線Fig.13 Speed curve of traditional algorithm

圖14 本文算法加工的速度曲線Fig.14 Speed curve of proposed algorithm

計算輪廓誤差，即由相鄰插補點構成的小線段到原始(CAD)橢圓弧的最遠距離，如圖15和圖16所示，傳統算法和本文算法的輪廓誤差均小于系統設置的最大值，但采用傳統算法加工時，插補點位于小線段或過渡曲線上，產生的最大誤差值個數較多，且平均誤差值為0.0041mm，而本文算法生成的擬合曲線能夠精確地表示橢圓弧，并且與原始(CAD)橢圓弧重合，使得插補點位于原始(CAD)橢圓弧上，其平均誤差僅為5.6773x10-5mm，精度更高.對于系統中設置的雙弓高誤差δmax= 0.01mm和輪廓誤差δc= 0.01mm，僅用于連續微小線段加工區域的識別與曲線擬合，雖然CAM編程相對于原始(CAD)橢圓弧存在0.001mm的誤差，且擬合橢圓弧相對于CAM編程存在0.01mm的擬合誤差，但這些誤差主要是折線到曲線的最大距離，又因為擬合橢圓弧與原始(CAD)橢圓弧重合，且插補點位于擬合橢圓弧，而不是折線，所以本文算法能夠在進給速度為3000mm/min的條件下，達到E-5mm量級的加工精度.

圖15 傳統算法加工的輪廓誤差Fig.15 Contour error of traditional algorithm

圖16 本文算法加工的輪廓誤差Fig.16 Contour error of proposed algorithm

4 結 語

1)本文提出的橢圓弧平滑壓縮插補算法采用預處理方式，不僅能夠識別出數控加工程序中的橢圓弧，而且還能以橢圓弧的幾何表示形式進行插補.

2)與傳統的樣條插補方法相比，提出的橢圓弧平滑壓縮插補算法能夠精確地識別并表示加工軌跡中的橢圓弧，提高了加工路徑的壓縮量和橢圓弧的加工質量，又因為本文用幾何形式的橢圓弧進行插補，能夠精確地計算出弧長對應的插補參數，降低了計算的復雜度和加工中速度的頻繁波動，提高了加工效率.

3)本文提出的橢圓弧平滑壓縮插補算法主要適用于在離散小線段形式下，橢圓弧部分的高速高精加工.

[1] Su Ai-lin，Zhang Chun-lei，Zhao Xue-feng，et al.A new method of ellipse interpolation-mapping method [J].Journal of Sichuan Union University，1999，3(1)：110-113.

[2] Hu Chi-bing，Cao Liang-wei，Yan Lin.An improved algorithm of elliptic interpolation [J].Journal of Lanzhou University of Technology，2007，33(6)：33-35.

[3] Yu Guang-wei，Zhao Yu-gang，Zhang Jian.Study of algorithm and path simulation by equal-error method to approach elliptical arc with beeline [J].Journal of Shandong University of Technology，2010，24(1)：90-92.

[4] Hou Zhi-qun.Setting-up and practice of equal segmental arc of ellipse curve points [J].Surveying and Mapping of Geology and Mineral Resources，2001，1(4)：14-16.

[5] Fukushima T.Precise and fast computation of a general incomplete elliptic integral of second kind by half and double argument transformations [J].Journal of Computational and Applied Mathematics，2011，235(14)：4140-4148.

[6] Boyd J P.Numerical，perturbative and chebyshev inversion of the incomplete elliptic integral of the second kind [J].Applied Mathematics and Computation，2012，218(13)：7005-7013.

[7] Fukushima T.Numerical inversion of a geVneral incomplete elliptic integral [J].Journal of Computational and Applied Mathematics，2013，237(1)：43-61.

[8] Gai Rong-li，Wang Yun-sen，Sun Yi-lan，et al.Survey on spline curve interpolation methods [J].Journal of Chinese Computer System，2012，33(12)：2744-2749.

[9] Tian Jun-feng，Lin Hu，Yao Zhuang，et al.Study on S-shape curve acceleration and deceleration control fast planning on CNC system [J].Journal of Chinese Computer System，2013，34(1)：168-172.

[10] Sun Shu-jie，Lin Hu，Zheng Liao-mo，et al.Algorithm of adaptive NURBS interpolation points calculation [J].Journal of Chinese Computer System，2014，35(4)：895-899.

[11] Zhang Xiao-hui，Yu Dong，Yang Dong-sheng，et al.Corner curve transition interpolation algorithm for high speed machining of micro-line segment [J].Journal of Mechanical Engineering，2010，46(19)：183-191.

[12] Zhang L B，You Y P，He J，et al.The transition algorithm based on parametric spline curve for high-speed machining of continuous short line segments [J].International Journal of Advanced Manufacturing Technology，2011，52(1)：245-254.

[13] Zhang Xiao-hui，Yu Dong，Hong Hai-tao，et al.Research on smooth compression interpolation algorithm in CNC machining [J].Journal of Mechanical Engineering，2011，47(5)：156-169.

[14] Zhang M，Yan W，Yuan C M，et al.Curve fitting and optimal interpolation on CNC machines based on quadratic B-splines [J].Science China Information Sciences，2011，54(7)：1407-1418.

[15] Yang Z Y，Shen L Y，Yuan C M，et al.Curve fitting and optimal interpolation for CNC machining under confined error using quadratic B-splines [J].Computer-Aided Design，2015，66 (7)：62-72.

[16] Farin G E.Geometric modeling：algorithms and new trends [M].Society for Industrial and Applied Mathematics(SIAM)，1987.

[17] Yin Chen-chen，Wu Wen-jiang.A new ellipse interpolation method based on arc length [J].Modular Machine Tool & Automatic Manufacturing Technique，2012，4(1)：1-4.

附中文參考文獻：

[1] 蘇愛林，張春雷，趙雪峰，等.橢圓插補的新方法-映射法[J].四川聯合大學學報，1999，3(1)：110-113.

[2] 胡赤兵，操良偉，閆 琳.一種橢圓插補的改進算法[J].蘭州理工大學學報，2007，33(6)：33-35.

[3] 于光偉，趙玉剛，張 鍵.橢圓弧等誤差直線擬合的算法研究與軌跡仿真[J].山東理工大學學報，2010，24(1)：90-92.

[4] 侯至群.等弧長橢圓曲線點的放樣與實踐[J].地礦測繪，2001，1(4)：14-16.

[8] 蓋榮麗，王允森，孫一蘭，等.樣條曲線插補方法綜述[J].小型微型計算機系統，2012，33(12)：2744-2749.

[9] 田軍鋒，林 滸，姚 壯，等.數控系統S型曲線加減速快速規劃研究[J].小型微型計算機系統，2013，34(1)：168-172.

[10] 孫樹杰，林 滸，鄭飂默，等.計算復雜度自適應的NURBS曲線插補算法[J].小型微型計算機系統，2014，35(4)：895-899.

[11] 張曉輝，于 東，楊東升，等.面向微線段高速加工的拐角曲線過渡插補算法[J].機械工程學報，2010，46(19)：183-191.

[13] 張曉輝，于 東，洪海濤，等.數控加工中的平滑壓縮插補算法研究[J].機械工程學報，2011，47(5)：156-169.

[17] 殷晨晨，吳文江.基于弧長的橢圓插補新算法[J].組合機床與自動化加工技術，2012，4(1)：1-4.