基于最小二乘原理的圓擬合及誤差評定算法研究*

黃力峰，汪 偉，吳南星

(景德鎮(zhèn)陶瓷大學 機械電子工程學院，江西 景德鎮(zhèn) 333403)

0 引言

在機械產(chǎn)品的生產(chǎn)加工中，往往會遇到對圓形零件特征參數(shù)的檢測和圓度誤差的評定，以提高零件加工和裝配的精度。在圓的參數(shù)求解中需要通過采集圓周上若干離散點的坐標，利用不同的方法求解出擬合圓的圓心坐標和半徑，其方法有近似算法、非線性最小二乘法、最小區(qū)域法優(yōu)化評定等[1]，其中最小二乘法分別由高斯和勒讓德獨立提出，是最早采用的擬合圓及圓度誤差評定的方法，因其評定模型簡單、易被工程人員掌握等優(yōu)點而得到較為廣泛的應用[2]。

在求解最小二乘圓的數(shù)學模型時，因其無法精確求解，需要對非線性方程組進行數(shù)值計算，國內(nèi)外學者對最小二乘方法的適應性和精度提高等方面進行了研究。Kasa[3]提出了基于最小二乘原理擬合圓的方法，并研究了采集點沿圓周分布的情況對圓心和半徑隨機誤差的影響。田社平等[4]給出了最小二乘通用算法，并與近似算法相比較，討論兩者的適用范圍和計算精度。Chan Y T等[5]通過加權(quán)平均所有組合的計算中心作為求解方法，避免了最小二乘法中近似或優(yōu)化帶來的問題。王秀梅等[6]針對最小二乘方法求解提出了改進的混沌優(yōu)化算法，并通過數(shù)據(jù)驗證了算法的有效性。

在基于最小二乘法對圓參數(shù)求解時，有基于非線性最小二乘問題的最優(yōu)化解法以及改變模型殘差項的代數(shù)解法[7]兩種方法。其中，最優(yōu)化解法存在如何選擇合理的初始值及迭代算法計算效率低等問題；而采用代數(shù)解法時，樣本中噪聲點即異常值的數(shù)量或采樣集中等因素會導致擬合精度下降。本文研究了代數(shù)解法的推導過程和最優(yōu)化解法中各類優(yōu)化函數(shù)的最佳選擇，以及噪聲點數(shù)目及分布情況的不同等影響因素，利用MATLAB進行數(shù)據(jù)仿真計算，分析兩種方法的擬合效果及計算精度，明確各方法的適用范圍及優(yōu)缺點，以實現(xiàn)最小二乘圓快速準確的求解。

1 最小二乘法模型的建立

圖1 最小二乘算法模型

使得各點到最小二乘圓上的距離的平方和最小，可得到模型一：

(1)

設(shè)最小二乘圓的方程式為：(xi-a)2+(yi-b)2=R2，可令殘差項為：ei=(xi-a)2+(yi-b)2-R2，則模型二可表示為：

.

(2)

由于模型一存在非線性方程組求解困難等問題，只能采用近似或優(yōu)化的方法求解[8]。針對模型二，式(2)中F的偏導數(shù)必須滿足最小化的條件為：

令

D=n∑xiyi-∑xi∑yi.

(3)

(4)

(5)

(6)

(7)

聯(lián)立式(3)～式(6)可解得模型二中A、B、C：

(8)

由式(8)可得到最小二乘意義下的擬合圓的圓心坐標(a,b)和半徑R為：

(9)

2 優(yōu)化算法的實現(xiàn)

MATLAB軟件的優(yōu)化工具箱(Optimization Toolbox)中含有一系列的優(yōu)化算法函數(shù)，當使用最小二乘法對最小二乘圓參數(shù)求解時，可以利用求解多變量函數(shù)極小值問題的fmincon、fminsearch函數(shù)和非線性最小二乘問題的lsqnonlin函數(shù)進行優(yōu)化求解[9]。而這些方法的收斂速度、求解效率和精度各不相同。引用文獻[10]一組圓度測量數(shù)據(jù)，分別對三種函數(shù)進行測試，其計算結(jié)果如表1所示。

表1 各優(yōu)化函數(shù)測試結(jié)果

此外，為測試三種優(yōu)化函數(shù)對初始值的依賴程度，利用MATLAB中rand命令產(chǎn)生2 000組隨機數(shù)作為初始值，三種函數(shù)達到全局最優(yōu)的組數(shù)占比分別為17.7%、14.3%、44.55%，故模型一直接利用MATLAB中解決非線性最小二乘問題的lsqnonlin優(yōu)化函數(shù)進行求解較好。該方法不涉及具體的算法，程序編制簡單，求解精度較高，可通過如下設(shè)置求解：

[x，resnorm]=lsqnonlin(fun，x0，lb，ub，options).

其中，x為使用迭代法搜索的最優(yōu)參數(shù)，為最小二乘圓的圓心坐標(a,b)和半徑r；resnorm為各項誤差，即目標函數(shù)的函數(shù)值；fun為目標函數(shù)；x0為變量的迭代初值；lb和ub為x值的上、下界，即x值的范圍設(shè)定；options定為默認設(shè)置。

設(shè)i=1,2,…,n(n>3)為被測實際圓周上的測量采樣點，最小二乘圓圓心坐標為(a，b)、半徑為R，令目標函數(shù)為：

然后調(diào)用MATLAB優(yōu)化函數(shù)lsqnonlin求得最小二乘圓的圓心坐標和半徑。

3 數(shù)據(jù)仿真與分析

在仿真實驗的設(shè)計中，已知樣本(x1,y1),(x2,y2),…,(xn,yn),要求擬合出一個模型(函數(shù))F，其預測值F(x)與樣本實際值y的誤差最小。考慮到樣本數(shù)據(jù)并不是真實值本身，假設(shè)真實模型(函數(shù))是f，則采樣值y=f(x)+ε，其中ε代表噪聲，所以擬合出函數(shù)后，可以采用偏差和均方誤差函數(shù)來度量其擬合的好壞程度。故在圓心為a=0，b=0，半徑R=10圓周上不同弧度范圍內(nèi)任意取n個樣本點，分別加上一定比例均值為0、方差為σ2的高斯隨機變量ε，σ∈[0,1]，從45°<θ<360°變化，n可取6、10、20、30，對每個不同的θ、n和σ值，生成數(shù)據(jù)以進行仿真實驗，每組重復實驗30次，以式(10)和式(11)計算每種方法的偏差B和均方誤差M的值，其值的大小作為衡量優(yōu)化方法和代數(shù)解法擬合圓效果的優(yōu)劣評價指標。

(10)

(11)

此外為驗證各方法在計算圓度誤差時的評定精度，從圓周上采集點均勻分布(12點)以及少量點集中分布(6點)兩種情況，計算最小二乘圓的參數(shù)及圓度誤差，并與文獻[4]和文獻[5]中的通用算法和質(zhì)心法進行比較。實驗數(shù)據(jù)分別為：

均勻采樣12點：

x=[23.30，15.67，3.85，-9.02，-19.47，-24.70，-23.31，-15.68，-3.83，9.03，19.48，24.71]；

y=[9.05，19.47，24.70，23.31，15.68，3.83，-9.03，-19.46，-24.70，-23.30，-15.67，-3.81]。

離散采樣6點：

x=[1.0，2.0，5.0，7.0，9.0，3.0]；

y=[7.0，6.0，8.0，7.0，5.0，7.0]。

仿真實驗結(jié)果如表2所示，表2中，θ為所取數(shù)據(jù)的弧度范圍。比較表2中兩種方法計算所得的偏差值與均方誤差值的大小(在零噪聲條件下即理想模型的偏差和均方誤差為0)可知：在小弧度范圍內(nèi)，隨著噪聲即σ的增大，兩種方法的擬合精度和穩(wěn)定變差，但優(yōu)化算法優(yōu)于代數(shù)解法，具有更好的擬合精度；隨著弧度的增大至大部分圓周或整圓時，兩種方法擬合圓精度大致相同且效果較好，能包容較大的噪聲點，因此選擇何種圓擬合法，視弧度、采集點數(shù)和噪聲情況而定。

表3和表4分別為均勻采樣12點和非均勻采樣6點計算結(jié)果。表3結(jié)果表明：在采集點均勻分布情況下，各種算法圓度誤差的評定精度相差不大。由表4可得出：對于采集點分布較集中的情況，代數(shù)解法計算圓度誤差值較小，評定精度相對較高，但優(yōu)化算法擬合效果比其余三種算法的效果要稍好，擬合效果如圖2～圖5所示。

表2 兩種方法計算結(jié)果

表3 均勻采樣12點計算結(jié)果 mm

表4 非均勻采樣6點計算結(jié)果 mm

圖2通用算法擬合效果圖3質(zhì)心法擬合效果

4 結(jié)論

研究了基于最小二乘法對圓擬合及誤差的評定，對其中的優(yōu)化算法和代數(shù)解法分別計算并進行比較分析。代數(shù)解法可直接求得準確的解析解，優(yōu)化算法利用MATLAB優(yōu)化工具箱選擇了易達到全局最優(yōu)解且計算效率較高的lsqnonlin函數(shù)進行優(yōu)化求解。仿真實驗結(jié)果表明：當采集點較多地分布在圓周時，優(yōu)化算法和代數(shù)解法擬合圓的精度大致相同；當采集點較少且集中分布時，優(yōu)化算法擬合精度更高，結(jié)果更加穩(wěn)定；在采集點均勻分布在圓周上時，各方法計算圓度誤差的結(jié)果基本相同；當較少的采集點非均勻分布時，代數(shù)解法計算的圓度誤差值較小。本文提出的最小二乘算法能準確快速求解出圓的參數(shù)及圓度誤差值，對實際中的測量應用具有指導意義。

圖4代數(shù)解法擬合效果圖5優(yōu)化算法擬合效果