基于BFGS算法的空間圓擬合方法

秦鋒，鹿松，張振虎

(1.中國核工業(yè)集團(tuán)三門核電有限公司，浙江 臺(tái)州 317112； 2.上?？辈煸O(shè)計(jì)研究院(集團(tuán))有限公司，上海 200093)

1 引 言

工程應(yīng)用中經(jīng)常需要對(duì)圓形物體進(jìn)行檢測(cè)或建模，大多數(shù)情況下圓形物體并非水平放置，傳統(tǒng)的平面圓擬合方法并不適用，而需要進(jìn)行空間圓的擬合?？臻g圓無直接公式描述，基于公式擬合的方法難以實(shí)施?？臻g圓通?？赏ㄟ^空間平面與球面相截得到，故大多數(shù)空間圓擬合采用平面與球面相截法，如文獻(xiàn)[1～4]先擬合空間平面求解空間圓所在平面法向量，然后利用球面中垂面的性質(zhì)構(gòu)造中垂面方程聯(lián)立空間平面方程解算空間圓圓心及半徑；文獻(xiàn)[5～7]通過擬合空間平面和一個(gè)球心在空間平面上的球面并解出空間圓所在平面法向量、空間圓圓心及半徑；文獻(xiàn)[8]分別擬合空間平面及球面，然后將球心投影到空間平面計(jì)算圓心及半徑。另有部分文獻(xiàn)將空間圓擬合問題變?yōu)槠矫鎴A擬合問題，如文獻(xiàn)[9-11]在擬合空間平面得到平面法向量后構(gòu)造旋轉(zhuǎn)矩陣通過坐標(biāo)轉(zhuǎn)換方式將空間圓變?yōu)槠矫鎴A，然后擬合平面圓得到平面圓圓心及半徑，最后再通過坐標(biāo)轉(zhuǎn)換得到空間圓圓心；文獻(xiàn)[12]先擬合空間平面，再將擬合點(diǎn)投影到空間平面上，建立誤差方程并聯(lián)立空間平面方程按附有限制條件的間接平差迭代求解。考慮到空間圓擬合數(shù)據(jù)中可能含有粗差，文獻(xiàn)[13～15]在上述算法的基礎(chǔ)上進(jìn)行了改進(jìn)，通過RANSAC(Random Sample Consensus，簡稱RANSAC)算法剔除粗差點(diǎn)，然后再采用上述算法求解空間圓擬合參數(shù)。

上述算法均采用兩步求解，即第一步計(jì)算出空間圓所在平面法向量，第二步求出空間圓圓心及半徑。本文采用一步求解，即通過Broyden-Fletcher-Goldfarb-Shanno(簡稱BFGS)算法一次性解算出全部空間圓擬合參數(shù)。為確保優(yōu)化函數(shù)收斂至極值，給出了一種簡便的空間圓擬合參數(shù)初始值計(jì)算方法。然后結(jié)合案例，與文獻(xiàn)[1]、文獻(xiàn)[9]中方法進(jìn)行了對(duì)比，驗(yàn)證了本文算法的準(zhǔn)確性和有效性。

2 數(shù)學(xué)模型及算法

2.1 空間圓擬合模型

空間圓所在平面方程可用下式表示：

ax+by+cz+d=0，s.t.a2+b2+c2=1

(1)

式中s.t.是subject to的縮寫，表示受限制于某條件。

為減少參數(shù)的計(jì)算，參考球面坐標(biāo)系表示方法，令：

a=cosαsinβ，b=sinαsinβ，c=cosβ

(2)

空間圓所在平面方程又可表示為：

cosαsinβx+sinαsinβy+cosβz+d=0

(3)

空間圓外一點(diǎn)(xi，yi，zi)，(1≤i≤n)到空間圓所在平面的法向偏差dvi可表示為：

dvi=cosαsinβxi+sinαsinβyi+cosβzi+d

(4)

假設(shè)空間圓的圓心坐標(biāo)為(x0，y0，z0)，其應(yīng)位于空間圓所在平面上，則有：

cosαsinβx0+sinαsinβy0+cosβz0+d=0

(5)

將式(4)與式(5)相減，得到空間圓外一點(diǎn)(xi，yi，zi)到空間圓所在平面的法向偏差dvi的另一種表述形式：

dvi=cosαsinβ(xi-x0)+sinαsinβ(yi-y0)

+cosβ(zi-z0)

(6)

設(shè)空間圓半徑為R，空間圓外一點(diǎn)(xi，yi，zi)到空間圓的徑向偏差dri可表示為：

(7)

則空間偏差dsi可表示為：

(8)

根據(jù)最小二乘法，以空間偏差平方和為對(duì)象構(gòu)建優(yōu)化函數(shù)fs：

(9)

在優(yōu)化函數(shù)fs取得最小值時(shí)得到的空間圓擬合參數(shù)為最佳參數(shù)，即求解以下數(shù)學(xué)問題：

minfs(w)，w=[x0，y0，z0，α，β，R]T

(10)

2.2 BFGS優(yōu)化算法

式(10)為無約束非線性最小二乘問題，根據(jù)文獻(xiàn)[16]，解決非線性最小二乘問題常用算法有Gauss-Newton(簡稱GN)算法、Levenberg-Marquardt(簡稱LM)算法、牛頓法及擬牛頓法。GN算法對(duì)于小殘量問題有較快的局部收斂速度，對(duì)于不是很嚴(yán)重的大殘量問題，有較慢的局部收斂速度，對(duì)于殘量很大的問題或非線性程度很大的問題往往不收斂。GN算法要求雅可比(Jacobi)矩陣滿秩，如果不滿秩，方法沒有定義，而且GN算法不一定總體收斂。LM算法最早由Levenberg提出，在GN算法加入正則化參數(shù)改善矩陣條件，后由Marquardt進(jìn)行了改進(jìn)。該方法是一種介于梯度下降法和GN算法之間的算法，常用于求解非線性最小二乘問題，在正則化參數(shù)大時(shí)相當(dāng)于梯度下降法，參數(shù)小時(shí)相當(dāng)于GN算法。在GN算法中，一旦雅可比矩陣奇異的情況發(fā)生，使得算法常常收斂到一個(gè)非駐點(diǎn)，LM算法通過添加正則化參數(shù)顯著改善了這一情況，但其并不能完全避免矩陣奇異的情況。牛頓法需要計(jì)算海森(Hessian)矩陣，但有的目標(biāo)函數(shù)海森矩陣很難計(jì)算或者計(jì)算量過大，導(dǎo)致牛頓法應(yīng)用受限。擬牛頓法通過利用目標(biāo)函數(shù)和一階導(dǎo)數(shù)的信息構(gòu)造出目標(biāo)函數(shù)的曲線近似，而不需要明顯形成海森矩陣，其計(jì)算量較牛頓法小并兼具收斂速度快的優(yōu)點(diǎn)，導(dǎo)致其比牛頓法更受歡迎、應(yīng)用范圍更廣。擬牛頓法中最著名的算法為Broyden-Fletcher-Goldfarb-Shanno(簡稱BFGS)算法，其由Broyden[17]、Fletcher[18]、Goldfarb[19]、Shanno[20]等學(xué)者于20世紀(jì)60年代末70年代初提出，是當(dāng)前最常用的擬牛頓算法。BFGS算法中校正矩陣始終保持正定性，并且算法具有超線性收斂速度，當(dāng)采用不精確線性搜索時(shí)，BFGS算法還具有總體收斂性質(zhì)。因此本文采用BFGS算法和Armijo不精確線性搜索策略[21]解決式(10)中的無約束非線性最小二乘問題。

首先對(duì)函數(shù)fs進(jìn)行偏微分，為便于展示，令：

對(duì)各參數(shù)的偏微分如下：

(11)

根據(jù)文獻(xiàn)[16]，算法步驟詳見以下算法1-1。

算法1-1：

步驟0，選取參數(shù)β∈(0，1)，σ∈(0，0.5)，初始點(diǎn)w0∈Rn，終止誤差0≤ε?1，初始對(duì)稱正定陣H0(通常取單位矩陣In)，令k：=0

步驟1，計(jì)算gk=?fs(wk)，若||gk||≤ε，停止計(jì)算，輸出wk作為近似極小點(diǎn)。

步驟2，計(jì)算dk=-Hkgk，得到dk。

步驟3，設(shè)mk是滿足下列不等式的最小非負(fù)整數(shù)m：

令αk=βmk，wk+1：=wk+αkdk

步驟4，由以下校正公式確定Hk+1

其中，sk=wk+1-wk，yk=?fs(wk+1)-?fs(wk)

步驟5，令k：=k+1，轉(zhuǎn)步驟1。

2.3 初始值計(jì)算

函數(shù)fs為非凸函數(shù)，存在多個(gè)駐點(diǎn)，因此，為確保函數(shù)在迭代過程中收斂到極小值，迭代初始值需要合理選取。以下是一種較簡單的擬合空間圓初始值計(jì)算方法。

將平面方程ax+by+cz+d=0兩側(cè)同除以d(d≠0)，令a′=a/d，b′=b/d，c′=c/d方程可以簡化為：

a′x+b′y+c′z=-1

(12)

因方程的三個(gè)參數(shù)之間存在著線性關(guān)系，故可采用線性最小二乘解算(a′，b′，c′)T，再將解算出的(a′，b′，c′)T單位化，就可根據(jù)式(2)反算出角度(α，β)的初始值(α0，β0)。

(13)

得到空間圓圓心初始值后，就可計(jì)算出圓心到各個(gè)點(diǎn)的距離，然后對(duì)這些點(diǎn)取平均值就可以得到空間圓的初始半徑R0。

(14)

2.4 空間圓擬合效果評(píng)價(jià)

空間圓擬合效果可通過均方根誤差(Root Mean Squared Error，RMSE)評(píng)價(jià)。各點(diǎn)擬合后的法向偏差dvi、徑向偏差dri及空間偏差dsi可分別用式(6)、式(7)及式(8)進(jìn)行計(jì)算?？臻g圓擬合法向偏差均方根誤差mdv、法向偏差均方根誤差mdr及空間偏差均方根誤差mds可根據(jù)下式計(jì)算。

(15)

3 工程算例

以激光跟蹤儀測(cè)量的某設(shè)備端面圓為例，進(jìn)行空間圓的擬合。采集的原始數(shù)據(jù)如表1所示。

表1 某設(shè)備端面圓測(cè)量原始數(shù)據(jù)

表2 各種算法計(jì)算結(jié)果(單位/mm)

通過表2可發(fā)現(xiàn)，本文算法法向誤差RMSE小于文獻(xiàn)[1]改進(jìn)及文獻(xiàn)[9]算法，說明前者算法空間平面擬合精度優(yōu)于后者。對(duì)比徑向誤差RMSE及空間誤差RMSE可發(fā)現(xiàn)，本文算法及文獻(xiàn)[9]算法在徑向誤差、空間誤差等方面均優(yōu)于文獻(xiàn)[1]改進(jìn)算法。從整體比較，本文算法和文獻(xiàn)[9]算法優(yōu)于文獻(xiàn)[1]改進(jìn)算法，本文算法略優(yōu)于文獻(xiàn)[9]算法，說明本文算法可行，精度較高。

4 結(jié) 論

本文借助球面坐標(biāo)表示方式構(gòu)造了空間圓擬合誤差函數(shù)，通過采用BFGS數(shù)值優(yōu)化算法，可一步解算出所有空間圓擬合參數(shù)，改變了傳統(tǒng)的兩步擬合空間圓方法。為確保優(yōu)化函數(shù)快速收斂至極值，本文還給出了一種簡便的空間圓擬合參數(shù)初始值計(jì)算方法。通過案例分析與驗(yàn)證，本文方法精度高，迭代快，可適用于高精度的空間圓擬合。