一種基于最小二乘法的離散點螺旋線式擬合算法

尤中桐 王太勇 劉清建 辛全琦

1.天津大學機械工程學院，天津，300350

2.天津職業大學機械工程與自動化學院，天津，300410

0 引言

在對非圓二次曲線與樣條曲線輪廓進行機械加工時，一般做法是在特定的容差范圍內用一系列折線替代原始曲線［1］，或用圓弧、雙圓弧逼近原始曲線［2?3］。但是這種方法存在擬合精度不夠高或程序段多等缺點，造成機床頻繁加減速，進而影響到實際的加工精度與效率。所以探索一種既能保證較高精度，又能減少數據量的曲線離散點擬合算法顯得十分重要。一些學者采用阿基米德螺線擬合曲線離散點，進行了數控加工的研究工作。王可等［4］在分析了阿基米德螺線特性與誤差計算方法后，用阿基米德螺線逼近渦旋盤輪廓，并將其與直線、圓弧擬合對比，指出阿基米德螺線的逼近特性優于直線和圓弧。張彥博［5］提出了一種用阿基米德螺線逼近橢圓來獲得曲線節點的算法，并通過一個實例證明，在精度滿足要求的前提下，該算法比等誤差直線逼近的程序段少得多。他們在求解阿基米德螺線參數時只用到了首末兩點，沒有兼顧到中間點對參數計算的影響。李澤蓉等［6］用平面二維優化方法求解阿基米德螺線參數來逼近大直徑圓弧，但是其計算涉及迭代，較為復雜，計算量也較大。

本文分析了阿基米德螺線上某一點處切矢、徑矢夾角與螺線參數的關系，為近似確定阿基米德螺線段的中心提供依據。擬合算法基于最小二乘原理，把中間數據點的誤差平方和作為目標函數，引入約束條件，考慮了中間數據點對阿基米德螺線參數計算的影響。通過對離散點集進行前尋、回溯，在保證精度的前提下，用阿基米德螺線與小線段完成離散點集的擬合。最后通過一個實例驗證所提算法，并將其與文獻［4?5］的方法和文獻［2］的圓弧擬合方法相對比，證明本文方法的優勢與應用價值。

1 基于不完全約束最小二乘原理的擬合算法

1.1 阿基米德螺線方程

如圖1所示，阿基米德螺線的極坐標方程可表示為

式中，ρ為螺線上的P點到螺線中心O的距離；θ為射線OP與X軸正方向的夾角；ρ0、v0為常數，ρ0為零度角對應的螺旋線上的點到螺線中心的距離；v0為單位角度變化量對應的距離增量。

圖1 平面阿基米德螺線Fig.1 Archimedes spiral in the plane

由式（1）可得到對應的直角坐標參數方程：

1.2 阿基米德螺線中心坐標的確定

如圖1所示，假設P點處螺線切矢T=（Tx，Ty），徑矢R=（Rx，Ry），二者的夾角為α，則有

依次取 θ=0o，60o，120o；以|ρ0/v0|為自變量研究其與夾角α的關系，在MATLAB中繪制二者的關系曲線，如圖2所示。

圖2 |ρ0/v0|與α的關系Fig.2 The relationship between|ρ0/v0|and α

就擬合計算得到的阿基米德螺線而言，|ρ0/v0|一般大于50，由圖2可知此時α接近90°，即T與R兩向量近似垂直。所以在確定螺線中心時，可以把阿基米德螺線當作圓弧處理，將首末兩端點的法線交點作為螺線中心。如圖3所示，實際計算時，以通過該段首點且垂直于首點與第二個點連線的直線為該段首點處的法線，同理，可確定該段尾點的法線。聯立兩直線方程，即可求得該段螺線中心點O的坐標。

圖3 計算螺線中心點坐標Fig.3 Calculation of center point of the spiral

1.3 ρ0與v0的求解

就一段阿基米德螺線的方程而言，其未知參數為ρ0、v0、中心坐標（x0，y0），上文已經確定了中心坐標，所以還剩ρ0、v0待求解。前人往往僅利用首末兩點構建方程組來求解ρ0與v0，忽視了中間點的影響。本文采用不完全約束最小二乘法來保證擬合精度，同時妥善處理擬合后曲線連續性的問題，考慮中間數據點對參數計算的影響，以較少段螺線擬合所有數據點。不完全約束最小二乘法是指，擬合計算后曲線的起點與終點是原有數據點列的起點與終點，螺線段之間的連接節點不一定是原數據點列中的點。如圖4所示，對以一段螺線擬合的數據點列而言，螺線終點一般并不與該數據點列的最后一個點重合。把該數據點列最后一個點對應的角度代入擬合計算出的螺線方程，可以得到螺線終點，用得到的螺線終點替代該數據點列的最后一個點，同時作為下一段螺線的起點，如圖4中的連接節點G。

圖4 不完全約束最小二乘螺線擬合Fig.4 Spiral fitting based on incomplete constraint least square method

假設某段螺線要擬合的數據點為P1(ρ1,θ1)，P2(ρ2,θ2)，…，Pn(ρn,θn)。構造目標函數S：

其中，ρi為點 Pi到螺線中心的距離，θi為點 Pi與螺線中心的連線同X軸正方向的夾角，二者均可通過點Pi的直角坐標求得。

引入螺線通過首點這一約束條件：

這樣對參數ρ0與v0的求解就轉化為約束條件下的優化問題。根據拉格朗日乘子法［7］，定義新函數

這樣關于參數ρ0與v0的約束條件下的優化問題就轉化成了關于F的無約束優化問題，確定F(ρ0,v0,k)的極值時，可利用其偏導數為零進行求解：

化簡式（8）得到線性方程組：

即可求解ρ0與 v0。

1.4 擬合誤差的計算

在圖5中，根據點集{ps,ps+1,ps+2,…,pi，…，pe-2,pe-1,pe}，利用前述計算方法得到的螺線方程ρ = ρ0+v0θ，中心點為 O′。除首點 ps外，點集中某點Pi的擬合誤差δpi的計算可采用以下方法：①計算點Pi在坐標系X′O′Y′中的坐標，進而得到其對應的極徑ρp與轉角α；②把α代入螺線方程，得到點對應的極徑；③pi與兩點極徑差值的絕對值即為點pi處的擬合誤差。即

圖5 擬合誤差計算Fig.5 The calculation of the fitting error

設ε為該段螺線的擬合誤差，則

1.5 數據點的前尋與回溯

對于離散點的擬合計算，由于無法預知其轉折點、尖點等特征信息，所以在擬合計算過程中，最好能遍歷所有數據點，以保證不會丟失某些關鍵數據點，從而不喪失原輪廓的特征信息。本文采用了一種前尋、回溯方法［8］來搜尋、遍歷所有數據點。如圖6所示，若點集的當前尾點為pi且擬合誤差滿足要求時，則向前搜尋點pi+1作為新的尾點，重新進行擬合計算，判斷是否滿足誤差條件，這稱為數據點前尋；若點集的當前尾點為pi且擬合誤差不滿足要求，則向后搜尋點pi?1作為新尾點，這時的點集即為某段待擬合的點集，這稱為數據點回溯。該前尋、回溯方法最大限度地擬合一段螺線較多的數據點。

圖6 前尋回溯數據點Fig.6 Search forward and backward the data points

1.6 擬合算法的流程

擬合算法流程見圖7，具體步驟如下：

（1）輸入待擬合的離散點集。

（2）讀入第一個數據點，將其作為首點。

（3）判斷待讀取數據點的個數：若為0，則結束擬合算法；若為1，則輸出首點與待讀取的最后一個點確定的小線段，結束擬合算法；若大于等于2，則讀入2個數據點，開始擬合計算。

（4）判斷該段擬合的首點與步驟（3）讀取的兩點是否共線：若共線，則輸出第一點與第三點確定的小線段，同時把第三點作為下段擬合的首點，轉到步驟（3）；若不共線，則根據上文的計算方法用這三點計算螺線參數，然后計算擬合誤差ε。

（5）判斷擬合誤差ε是否滿足要求：若滿足要求，則轉到步驟（7）；若不滿足要求，則轉到步驟（6）。

（6）如果擬合的數據點數超過3，則回溯一個數據點，根據這時的數據點集計算螺線參數并輸出，同時根據這組參數更新尾點坐標，作為下段擬合的首點；如果擬合的數據點數是3，則輸出前2個點確定的小線段，同時尾點回溯一個數據點，作為下段擬合的首點。轉到步驟（3）。

（7）判斷是否還有數據點待讀取：若有，則前尋一個數據點，然后根據這時的數據點集重新計算螺線參數與擬合誤差ε，轉到步驟（5）；若沒有，則輸出上次計算出的螺線參數，結束擬合算法。

圖7 擬合算法流程圖Fig.7 The flow chart of fitting algorithm

2 算法仿真與分析

上述算法已經通過Visual C++編程實現，現以一含有255個數據點的點集為例，在計算機上進行該算法的仿真測試分析。測試環境如下：操作系統為Windows 7旗艦版；內存為8 GB；CPU為Inter Core i5?4570 3.20 GHz。

給定擬合精度5 μm，采用本文算法計算得到所需阿基米德螺線段數為27，如圖8所示，實心點為離散點，空心點為螺線端點，實線為螺線，表明這些螺線段很好地擬合了離散點。

圖8 基于最小二乘法的阿基米德螺線擬合結果Fig.8 Spiral fitting result based on least square method

給定擬合精度5 μm，采用文獻［2］的圓弧擬合算法去擬合這些數據點，計算得到所需圓弧段數為56，如圖9所示。

圖9 基于文獻［2］最小二乘法的圓弧擬合結果Fig.9 Circle fitting result based on the least square method from reference［2］

給定擬合精度5 μm，采用文獻［4?5］中僅用首末兩點的方法計算得到所需螺線段數為31，如圖10所示，比本文擬合算法多出4段螺線，且增加的4段螺線均出現在曲率變化劇烈的部分，可見本文提出的算法要優于文獻［4?5］方法。對比以上3種結果可知，相同擬合精度下，基于最小二乘法的阿基米德螺線擬合算法能使數據量得到縮減。

使用53段阿基米德螺線，分別采用本文提出的算法與文獻［4?5］的方法去擬合該離散點集，再使用53段圓弧采用文獻［2］的方法去擬合該離散點集，得到各自的擬合誤差如圖11～圖13所示。段數相同時，基于最小二乘法的阿基米德螺線擬合算法的精度更高，再將其與已有的阿基米德螺線插補算法［9?10］結合，便可實現離散點集基于離散點螺旋線式擬合算法的數控加工。

圖10 基于文獻［4-5］方法的阿基米德螺線擬合結果Fig.10 Spiral fitting result based on the method from reference［4-5］

圖11 基于最小二乘法的阿基米德螺線擬合誤差Fig.11 Sspiral fitting error based on the least square method

圖12 基于文獻［4-5］方法的阿基米德螺線擬合誤差Fig.12 Spiral fitting error based on the method from reference［4-5］

圖13 基于文獻［2］最小二乘法的圓弧擬合誤差Fig.13 Circle fitting error based on the least square method from reference［2］

為了比較本文所提算法與前人算法的時間花費，在給定相同的數據點集（上述255個數據點）與擬合精度（10 μm、5 μm、1 μm）的前提下，通過調用Windows系統的API高精度計時函數Query?PerformanceCounter（）可得到3種算法的運行時間，如表1所示。可知本文提出的擬合算法在計算效率上介于文獻［4?5］的擬合算法與文獻［2］的擬合算法之間。

表1 三種算法的運行時間Tab.1 The running time of the three algorithms

3 結論

（1）阿基米德螺線上某點處的切矢、徑矢夾角隨 |ρ0/v0|的增大而趨近于90°。

（2）通過最小二乘法引入了中間數據點對螺線參數計算的影響。與前人算法相比，本文算法具有更好的數據壓縮效果。

（3）在相同段數的前提下，使用基于最小二乘法的阿基米德螺線擬合離散點集，具有更高的擬合精度。