秦鋒，鹿松，張振虎

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

1 引 言

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

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

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

2.1 空間圓擬合模型

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

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

(1)

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

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

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

(2)

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

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

(3)

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

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

(4)

假設空間圓的圓心坐標為(x0，y0，z0)，其應位于空間圓所在平面上，則有：

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

(5)

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

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

+cosβ(zi-z0)

(6)

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

(7)

則空間偏差dsi可表示為：

(8)

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

(9)

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

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

(10)

2.2 BFGS優(yōu)化算法

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

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

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

(11)

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

算法1-1：

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

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

步驟2，計算dk=-Hkgk，得到dk。

步驟3，設mk是滿足下列不等式的最小非負整數(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 初始值計算

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

將平面方程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)

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

(13)

得到空間圓圓心初始值后，就可計算出圓心到各個點的距離，然后對這些點取平均值就可以得到空間圓的初始半徑R0。

(14)

2.4 空間圓擬合效果評價

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

(15)

3 工程算例

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

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

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

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

4 結(jié) 論

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