八心圓弧擬合橢圓誤差的理論解析及最優解

王國順， 唐立波

(大連交通大學機械工程學院，遼寧 大連 116028)

針對用圓弧擬合橢圓的誤差分析沒有精確理論解析這一現狀，根據橢圓圖形本身的特性，利用圖形學理論，提出了等差擬合弧的概念，在給定橢圓長短軸的情況下，就可以推導出八圓弧擬合橢圓誤差的理論解析表達式，該解析式是關于法向誤差變量的超越方程，為求解超越方程，首先確定了定解區間，然后利用數學中的二分法對其求解。對于一個具體給定的橢圓，可以求出使用八心圓弧擬合的法向最小誤差，即八心圓弧擬合的最高精度，如果要求更高的擬合精度，只能通過增加擬合圓弧的數量來實現；計算過程采用AutoCAD軟件內嵌的Visual LISP語言編程實現。在數控加工中，采用擬合圓弧作為刀具的運行控制軌跡，加工出來的曲線總是存在一定的誤差，為此，只能將橢圓誤差控制在允許的變動范圍內，該范圍定義為橢圓公差帶[1]。彭云柯[2]采用約束非線性切比雪夫逼近論，給出了一種用四段圓弧逼近橢圓的方法，可使最大逼近誤差達到最小，從而提高了數控機床和機器人的橢圓加工作業的精度，并且提供了最佳半徑的相當準確的近似計算公式可作為最后結果使用，從而避免編程迭代。周曉銘等[3]發現人們對于圓弧插補的研究，針對大多數集中在圓弧逼近算法本身而沒有考慮機床的最大切割半徑，一旦擬合圓弧半徑大于機床半徑時，機床不能處理的情況，提出了一種新的基于二分法的曲線圓弧擬合方法，該方法繼承單圓弧等弦高誤差控制法和傳統二分法取中點的習慣，以最大弦高點為節點分割曲線，以相鄰三節點作圓弧擬合曲線，在擬合半徑大于數控系統能處理的最大半徑時，用兩直線擬合該段圓弧，方便數控自動編程。樊琳[4]研究了列表曲線擬合問題，根據已給出曲線上某些坐標點，由曲線坐標離散的型值點來描述工件表面形狀曲線的大致走向，但沒有給出方程。韓鴻鸞[5]、丁克會和席平原[6]分析了四心圓弧擬合橢圓的偏差變化規律，給出的誤差是沿橢圓的擬合圓弧徑向。這與橢圓實際要求控制的法向誤差不完全一致，因而不能準確描述四心圓弧擬合橢圓的精度，擬合的近似程度較模糊。沈紀桂等[7]介紹了傳統的曲線擬合方法，即單圓弧一次逼近法、兩次逼近法和雙圓弧一次逼近法，上述三種方法中，兩次逼近法和雙圓弧一次逼近法中的雙圓弧均是在兩相鄰型值點間進行擬合，其擬合圓弧段數多，計算、輸出的數據量大，而單圓弧一次逼近法雖采用大跨距擬合，但具有一定的盲目性，擬合精度差。王林鴻和陳永輝[8]研究了用傳統的四心圓法作近似橢圓與實際橢圓的誤差分布規律，結果表明，除橢圓的四個交點以及與坐標軸的四個交點外，處處存在偏差，其中大圓弧段偏差較小，小圓弧段偏差最大，近似橢圓比精確橢圓面積大，橢圓越扁平，偏差越大，在極端扁平的情況下，最大偏差超過長半軸的3%、短半軸的15%。文獻[9-11]研究了數控加工領域中的關于橢圓及曲線擬合的很多方法，普遍認為，根據曲線段斜率變化規律確定圓心方位、以凸弧代替凸曲線的等弧長圓弧擬合方法是可行的，并在該方法中使用了級數展開來簡化超越三角方程，避免了迭代求解的麻煩，這種方法的不足是沒能從理論上給出擬合誤差的解析式，因而不能更好的描述擬合精度。文獻[12-14]認為雙圓弧方法適用于對數量較大的實驗數據作近似擬合，該方法在數控加工及圖形繪制中可以得到較好的應用，雙圓弧曲線是一階連續的，數控加工后零件的光滑度好，從而，它可有效地對任意列表平面曲線和空間曲面作數控加工。唐立波等[15]給出了四心圓弧擬合橢圓的解析算法及其最優解，但是擬合精度更高時無法滿足。

綜上所述，研究側重在數控加工應用領域；在圓弧擬合橢圓或曲線方面的理論研究很多，已有的研究一般是采用雙圓弧或其他一些近似擬合方法，而沒有給出關于誤差的精確解析，致使擬合誤差大，精度低，導致數控加工效率低，時間長，從而增加了產品的成本。本文導出的解析超越方程，填補了八心圓弧擬合橢圓法向誤差理論解析的空白。

1 八心圓弧擬合橢圓的邊界條件與約定

橢圓的參數方程：

其中α為橢圓的參數角，0°≤α≤360°，a、b為橫、縱半軸長，b＞a＞0，且a≠b；如果a=b，則為圓的參數方程。由于橢圓具有對稱性，故以下只按a＜b時第一象限的橢圓分析推導解析式，其他象限可由第一象限的結果作對稱變換獲得，a＞b時，可交換a、b求解，然后將結果相對直線x=y作對稱變換。為了較好地逼近理想橢圓(如圖1)，依據橢圓的特點，八心圓弧擬合橢圓應同時滿足以下5個邊界條件：

(1) 過橢圓端點；

(2) 相對橢圓中心對稱；

(3) 保持一階導數連續；

(4) 在給定的橢圓公差帶內；

(5) 分段擬合圓弧的半徑變化趨勢與橢圓的曲率半徑變化趨勢一致。

圖1 橢圓(1/2)

給定橢圓法向誤差δ(為討論方便，以下取δ>0)時，可以確定橢圓內外兩條等距線，如圖2所示，兩條等距線確定的范圍稱橢圓的誤差范圍，若此范圍與橢圓的公差帶相同，則這個誤差范圍就是橢圓的公差帶。

圖2 橢圓法向誤差δ

八心圓弧擬合橢圓，滿足上述條件下在第一象限共分三段，如圖3所示。其中，過橫軸端點的圓弧半徑最大，稱為大圓弧；過縱軸端點的圓弧最小，稱為小圓弧，剩下一段圓弧稱為中間圓弧。

圖3 擬合圓弧名稱

為便于討論，本文中的字母O、R、tP分別代表圓心、半徑和切點，小圓弧對應的變量符號角標用1表示，大圓弧對應的變量符號角標用2表示，中間圓弧對應的變量符號角標用3表示，圓心對應的坐標分別為：xo1,yo1、xo2,yo2和xo3,yo3，其中xo1=yo2=0。

2 橢圓的等差擬合圓弧

2.1 定義

在一個橢圓給定法向誤差δ時，分別切于δ所確定的橢圓內外兩條等距線的圓弧定義為橢圓等差擬合圓弧(如圖4中Pt、Pt1之間的圓弧)。

圖4 等差擬合圓弧

在區間[0,δ0]上取任一誤差δ，則區間[0,π/2]上任一參數角α對應橢圓上一點P，則橢圓上過點P的法向斜率：

過點P的法線方程：

再取區間[0,α]上任意一參數角α1，α1對應橢圓上一點P1，則過點P1的法線方程：

式(3)和式(4)聯立求解，可得交點O，其坐標：

OP交橢圓內等距線于點tP，O1P交橢圓外等距線于點Pt1，若：

則Pt、Pt1是在以O為圓心的同一個圓上，根據上述等差擬合圓弧的定義可知，此時以Pt、Pt1為端點的圓弧是橢圓的一等差擬合圓弧。以下稱Pt為等差擬合圓弧的左端點；稱Pt1為等差擬合圓弧的右端點。

2.2 橢圓最小等差擬合圓弧

根據上述等差擬合圓弧的定義，當α=π/2時，等差擬合圓弧最小，其左端點是縱軸與橢圓內等距線的交點Pt(0,b-δ)，如圖5所示。

圖5 最小等差擬合圓弧

設右端點對應的橢圓參數角為αm，過O和OPt1與橢圓的交點P作輔助線交M，由圖示關系可知OPt1的斜率：

將式(2)代入式(8)可得：

故最小等差擬合圓弧的圓心縱坐標和半徑分別為：

Pt1的坐標為：

由式(9)～(12)可以看出，只要能求解mα，最小等差擬合圓弧即可確定。為此，由式(7)和圖5所示幾何關系可得：

將式(10)代入化簡可得：

該方程是關于sinmα的一元二次方程。由于設定在第一象限求解，sinmα大于0，故舍去負解，將所得sinmα的正解化簡，根據反三角函數關系可得：

綜上所述，最小等差擬合圓弧由a、b、δ確定。

2.3 橢圓最大等差擬合圓弧

當α1=0時，等差擬合圓弧最大，其右端點是橫軸與橢圓外等距線的交點Pt1(a+δ,0)，如圖6所示。

圖6 最大等差擬合圓弧

參照最小等差擬合圓弧的解析方法，求解出αn，從而最大等差擬合圓弧即可確定，αn的解析表達見式(17)。

設左端點Pt對應的橢圓參數角為αn，過O和OPt與橢圓的交點P作輔助線交N，由圖示關系可知OPt的斜率：

2.4 橢圓等差擬合中間圓弧的求解

根據上述分析，給定橢圓法向誤差δ時，等差擬合圓弧的左端點所對應的橢圓參數角的定解區間是：α∈[αn,π/2]；右端點所對應的橢圓參數角的定解區間應是：α1∈[0,αm]。若α1=0，α=π/2，則橢圓的最大、最小等差擬合圓弧重合，其圓心在原點，設此時對應的法向誤差為δ′，由式(13)可得：

顯然，法向誤差δ≤δ′時，才存在等差擬合圓弧，進而在區間[αn,π/2]上任取α，由式(3)～(6)聯立求解α1可確定一個等差擬合圓弧。但式(3)～(6)聯立得到的是關于α1的超越方程，為此，設：

根據該式繪制F(α1)的變化曲線得到的結果如圖7所示，可以看出F(α1)在區間[0,α)上關于α1是單調減函數，且F(0)＞0，α1趨向α時F(α1)＜0，故可在區間[0,α)上用二分法[16]求解得到α1，至此，橢圓的等差擬合圓弧隨給定的α完全確定。

圖7 函數F(α1)隨α1變化曲線

用同樣的方法，在區間[0,αm]上任取α1求解α也可確定對應的等差擬合圓弧。

綜上所述，在給定橢圓一個法向誤差δ時，橢圓的等差擬合圓弧只與α或α1獨立相關，將α、α1代入式(5)、(6)得出xo3、yo3。

3 八心圓弧擬合橢圓解析

3.1 八心圓弧擬合橢圓的定解分析

在區間(0,δ0)上給定一個法向誤差δ時，擬合橢圓的大、小圓弧可按2.2節和2.3節的方法進行求解；中間擬合圓弧可按2.4節的方法求解，其兩個端點所對應的橢圓參數角的定解區間為α∈[αn,π/2]，α1∈[0,αm]，在定解區間內，有無數個橢圓等差擬合圓弧，如果等差擬合圓弧分別與大、小圓弧內切，則這三個圓弧就可構成滿足1節所述5個邊界條件的八心擬合圓弧。

3.2 與小圓弧相切的橢圓等差擬合圓弧

區間[αn,π/2]上任意一參數角α對應的等差擬合圓弧可按2節所述求解，設其半徑為R3，其圓心O3與小圓弧的圓心O1之間的距離為L1，則該等差擬合圓弧與小圓弧內切的條件是：

其中，

由2節可知，O3、R3均由α確定，故L1也由α確定，因此式(17)是關于α的獨立方程，但導出此方程卻是關于α的超越方程，為此，取

根據該式繪制F(α)的變化曲線如圖8所示，可以看出F(α)在區間[αn,π/2]上關于α是單調減函數，且F(αn)＞0，F(π/2)＜0，故可在區間[αn,π/2]上用二分法[16]求解方程F(α)=0，得到α，其所對應的等差擬合圓弧即是與小圓弧相切的等差擬合圓弧。

圖8 函數F(α)隨α變化曲線

3.3 中間圓弧求解

由3.2節可知，給定橢圓法向誤差δ時即可求得與小圓弧內切的一等差擬合圓弧，設其半徑為R3，其圓心與大圓弧圓心的距離為L2，則其與大圓弧內切的條件是：

其中，

該方程中的3個變量都只與δ相關，所以此方程是關于δ的獨立方程，但導出此方程卻是關于δ的超越方程，為此，取：

根據該式繪制F(δ)的變化曲線如圖9所示，可以看出F(δ)在區間[0,δ0]上關于δ是單調減函數，且F(0)＞0，F(δ0)＜0，故可在區間[0,δ0]上用二分法[16]求解方程F(δ)=0，得到的解設為δ3，δ3所對應的等差擬合圓弧與大、小圓弧同時相切，是擬合橢圓的八心圓弧中的中間圓弧。

圖9 函數F(δ)隨δ變化曲線

上述解析采用Visual LISP語言編程，在AutoCAD環境下運行，可快速求得八心圓弧擬合橢圓的結果。

4 實例分析

按照上述方法，以長半軸b等于60，短半軸a分別取5、10、15、20、25、30、35、40、45、50、55的系列橢圓為例，編程計算的八心圓弧擬合橢圓的誤差，結果見表1，誤差分布見圖10。

圖10 誤差分布圖

從表1可以看出：對于某一固定的橢圓，δ是用八心圓弧擬合的最小誤差，也是八心圓弧擬合所能達到的最高精度，如果要求更高的擬合精度，則必須增加擬合圓弧的段數；從圖10可以看出：橢圓越飽滿(即長、短軸相差較小)或極端扁平(即長、短軸相差甚大)時擬合精度高。

表1 八心圓弧擬合橢圓誤差的計算結果

5 結 論

橢圓的八心擬合圓弧有無數解。本文提出了橢圓等差擬合圓弧的概念，定義了等差擬合弧，從理論上推導出八心圓弧擬合橢圓時法向誤差的超越方程，并使用二分法進行求解。上述方法能根據橢圓長、短半軸快速求出八心圓弧擬合橢圓的最大法向誤差，即擬合精度，若擬合精度的絕對值大于或等于最小法向誤差時均可用該八心圓弧擬合，若小于最大法向誤差的更高精度橢圓是不能用八心圓弧擬合的，必須增加擬合圓弧段數，這就為八心圓弧是否可以擬合給定法向誤差的橢圓提供了準確的判斷依據。

將來會在此研究基礎上，繼續研究更多段圓弧擬合橢圓，以提高擬合精度，滿足工程需求。

[1] 劉巽爾.形狀和位置公差[M].北京: 中國標準出版社,2004: 34-55.

[2] 彭云柯.橢圓加工作業中的四圓弧切比雪夫逼近[J].機械工程學報,2004,40(12): 168-171.

[3] 周曉銘,陶俊才,秦海鷗.基于二分法的曲線等誤差圓弧擬合算法及其應用[J].南昌大學學報(理科版),2008,32(1): 100-102.

[4] 樊 琳.列表曲線的圓弧擬合[J].機械工藝師,1998,(12): 7-8.

[5] 韓鴻鸞.橢圓插補原理存在的缺陷[J].機械制造與研究,2006,35(3): 65-67.

[6] 丁克會,席平原.基于MATLAB的最大誤差雙圓弧逼近曲線的算法及實現[J].機械傳動,2007,31(6):57-59.

[7] 沈紀桂,方志梅,查顯亮.優化雙圓弧擬合的曲面加工[J].寧夏工學院學報(自然科學版),1997,9(1):5-10.

[8] 王林鴻,陳永輝.用四心圓法作近似橢圓的偏差分析[J].機械工程師,2002,(8): 59-60.

[9] 朱 虹,王立軍,楊 海.基于圓弧擬合的非圓二次曲線加工技術研究[J].機械設計與制造,2012,(12):227-229.

[10] 劉吉安,張 平.非圓曲線的雙圓弧擬合及數控加工的程序編制[J].機電工程技術,2009,38(6): 44-46,91.

[11] 姚必強,姚 進.數控加工曲線的等弧長圓弧擬合方法[J].四川大學學報(工程科學版),2008,40(1):171-174.

[12] 虞銘財,楊勛年,汪國昭.整體最優雙圓弧擬合[J].高校應用數學學報A輯(中文版),2004,19(2): 225-232.

[13] 沈紀桂,方志梅,李凌豐.優化雙圓弧擬合[J].工程圖學學報,1995,16(1): 53-56.

[14] 閻 童,王 琦.雙圓弧擬合在輪廓仿真加工中的應用[J].小型微型計算機系統,1998,19(10): 57-60.

[15] 唐立波,張 旭,李寶良.四心圓弧擬合橢圓的解析算法及其最優解[J].大連交通大學學報,2009,30(4):18-21.

[16] 關 治,陸金甫.數值分析[M].北京: 清華大學出版社,2006: 55-65.