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

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

1.天津大學(xué)機(jī)械工程學(xué)院，天津，300350

2.天津職業(yè)大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，天津，300410

0 引言

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

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

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

1.1 阿基米德螺線方程

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

式中，ρ為螺線上的P點(diǎn)到螺線中心O的距離；θ為射線OP與X軸正方向的夾角；ρ0、v0為常數(shù)，ρ0為零度角對(duì)應(yīng)的螺旋線上的點(diǎn)到螺線中心的距離；v0為單位角度變化量對(duì)應(yīng)的距離增量。

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

由式（1）可得到對(duì)應(yīng)的直角坐標(biāo)參數(shù)方程：

1.2 阿基米德螺線中心坐標(biāo)的確定

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

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

圖2 |ρ0/v0|與α的關(guān)系Fig.2 The relationship between|ρ0/v0|and α

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

圖3 計(jì)算螺線中心點(diǎn)坐標(biāo)Fig.3 Calculation of center point of the spiral

1.3 ρ0與v0的求解

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

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

假設(shè)某段螺線要擬合的數(shù)據(jù)點(diǎn)為P1(ρ1,θ1)，P2(ρ2,θ2)，…，Pn(ρn,θn)。構(gòu)造目標(biāo)函數(shù)S：

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

引入螺線通過首點(diǎn)這一約束條件：

這樣對(duì)參數(shù)ρ0與v0的求解就轉(zhuǎn)化為約束條件下的優(yōu)化問題。根據(jù)拉格朗日乘子法［7］，定義新函數(shù)

這樣關(guān)于參數(shù)ρ0與v0的約束條件下的優(yōu)化問題就轉(zhuǎn)化成了關(guān)于F的無約束優(yōu)化問題，確定F(ρ0,v0,k)的極值時(shí)，可利用其偏導(dǎo)數(shù)為零進(jìn)行求解：

化簡(jiǎn)式（8）得到線性方程組：

即可求解ρ0與 v0。

1.4 擬合誤差的計(jì)算

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

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

設(shè)ε為該段螺線的擬合誤差，則

1.5 數(shù)據(jù)點(diǎn)的前尋與回溯

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

圖6 前尋回溯數(shù)據(jù)點(diǎn)Fig.6 Search forward and backward the data points

1.6 擬合算法的流程

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

（1）輸入待擬合的離散點(diǎn)集。

（2）讀入第一個(gè)數(shù)據(jù)點(diǎn)，將其作為首點(diǎn)。

（3）判斷待讀取數(shù)據(jù)點(diǎn)的個(gè)數(shù)：若為0，則結(jié)束擬合算法；若為1，則輸出首點(diǎn)與待讀取的最后一個(gè)點(diǎn)確定的小線段，結(jié)束擬合算法；若大于等于2，則讀入2個(gè)數(shù)據(jù)點(diǎn)，開始擬合計(jì)算。

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

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

（6）如果擬合的數(shù)據(jù)點(diǎn)數(shù)超過3，則回溯一個(gè)數(shù)據(jù)點(diǎn)，根據(jù)這時(shí)的數(shù)據(jù)點(diǎn)集計(jì)算螺線參數(shù)并輸出，同時(shí)根據(jù)這組參數(shù)更新尾點(diǎn)坐標(biāo)，作為下段擬合的首點(diǎn)；如果擬合的數(shù)據(jù)點(diǎn)數(shù)是3，則輸出前2個(gè)點(diǎn)確定的小線段，同時(shí)尾點(diǎn)回溯一個(gè)數(shù)據(jù)點(diǎn)，作為下段擬合的首點(diǎn)。轉(zhuǎn)到步驟（3）。

（7）判斷是否還有數(shù)據(jù)點(diǎn)待讀?。喝粲?，則前尋一個(gè)數(shù)據(jù)點(diǎn)，然后根據(jù)這時(shí)的數(shù)據(jù)點(diǎn)集重新計(jì)算螺線參數(shù)與擬合誤差ε，轉(zhuǎn)到步驟（5）；若沒有，則輸出上次計(jì)算出的螺線參數(shù)，結(jié)束擬合算法。

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

2 算法仿真與分析

上述算法已經(jīng)通過Visual C++編程實(shí)現(xiàn)，現(xiàn)以一含有255個(gè)數(shù)據(jù)點(diǎn)的點(diǎn)集為例，在計(jì)算機(jī)上進(jìn)行該算法的仿真測(cè)試分析。測(cè)試環(huán)境如下：操作系統(tǒng)為Windows 7旗艦版；內(nèi)存為8 GB；CPU為Inter Core i5?4570 3.20 GHz。

給定擬合精度5 μm，采用本文算法計(jì)算得到所需阿基米德螺線段數(shù)為27，如圖8所示，實(shí)心點(diǎn)為離散點(diǎn)，空心點(diǎn)為螺線端點(diǎn)，實(shí)線為螺線，表明這些螺線段很好地?cái)M合了離散點(diǎn)。

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

給定擬合精度5 μm，采用文獻(xiàn)［2］的圓弧擬合算法去擬合這些數(shù)據(jù)點(diǎn)，計(jì)算得到所需圓弧段數(shù)為56，如圖9所示。

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

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

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

圖10 基于文獻(xiàn)［4-5］方法的阿基米德螺線擬合結(jié)果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 基于文獻(xiàn)［4-5］方法的阿基米德螺線擬合誤差Fig.12 Spiral fitting error based on the method from reference［4-5］

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

為了比較本文所提算法與前人算法的時(shí)間花費(fèi)，在給定相同的數(shù)據(jù)點(diǎn)集（上述255個(gè)數(shù)據(jù)點(diǎn)）與擬合精度（10 μm、5 μm、1 μm）的前提下，通過調(diào)用Windows系統(tǒng)的API高精度計(jì)時(shí)函數(shù)Query?PerformanceCounter（）可得到3種算法的運(yùn)行時(shí)間，如表1所示?？芍疚奶岢龅臄M合算法在計(jì)算效率上介于文獻(xiàn)［4?5］的擬合算法與文獻(xiàn)［2］的擬合算法之間。

表1 三種算法的運(yùn)行時(shí)間Tab.1 The running time of the three algorithms

3 結(jié)論

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

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

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