基于Levenberg-Marquardt算法的圓柱面擬合研究

秦 鋒，張振虎，張 威

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

在AP1000(Advanced Passive Pressurized Water Reactor，先進(jìn)壓水反應(yīng)堆)核電項(xiàng)目核島安裝工程中，一回路主設(shè)備安裝精度要求很高，例如主泵其安裝精度要求優(yōu)于0.3 mm，且精密對接區(qū)域空間狹小、無法通視，傳統(tǒng)的高精度全站儀，如徠卡TM50全站儀，無法滿足測量精度及通視要求，因此，需要引入更高精度的測量儀器及測量新方法輔助這些設(shè)備的安裝。激光跟蹤儀作為當(dāng)前精度最高的大范圍移動式坐標(biāo)測量系統(tǒng)，已在飛機(jī)制造、高能物理等領(lǐng)域得到了應(yīng)用。為了確保核島主設(shè)備安裝精度，三門核電項(xiàng)目在核電項(xiàng)目中率先采用了激光跟蹤儀測量整體建模技術(shù)輔助設(shè)備精密安裝。因核島主設(shè)備外形多為圓柱形，故采用激光跟蹤儀測量建模時需要進(jìn)行圓柱面的精確擬合。圓柱面擬合的關(guān)鍵在于確定圓柱面的軸線和半徑。文獻(xiàn)[1]通過高斯圖法計(jì)算圓柱面擬合參數(shù)，原理較復(fù)雜；文獻(xiàn)[2]在文獻(xiàn)[1]基礎(chǔ)上進(jìn)行了改進(jìn)，將高斯圖法得到的圓柱面參數(shù)作為初始值，并采用非線性最小二乘進(jìn)行優(yōu)化計(jì)算；文獻(xiàn)[3]根據(jù)矩陣特征分解、測量坐標(biāo)系和標(biāo)準(zhǔn)坐標(biāo)系轉(zhuǎn)換關(guān)系提出了按特征值、旋轉(zhuǎn)角和平移量為參數(shù)的二次曲面擬合方法；文獻(xiàn)[4]基于遺傳算法進(jìn)行圓柱面的擬合；文獻(xiàn)[5-8]通過圓柱面在平面上的投影圓度判別估算圓柱面參數(shù)初值，然后，采用優(yōu)化算法或平差方法進(jìn)行圓柱面擬合；文獻(xiàn)[9]先以點(diǎn)到圓柱中軸線距離平方與半徑平方之差為對象構(gòu)造優(yōu)化函數(shù)得到圓柱面參數(shù)初始值，再以點(diǎn)到圓柱面的距離為對象構(gòu)造誤差函數(shù)進(jìn)行精確求解；文獻(xiàn)[10-11]基于坐標(biāo)轉(zhuǎn)換的方式進(jìn)行圓柱面的擬合；文獻(xiàn)[12]綜合采用了文獻(xiàn)[9]和文獻(xiàn)[10-11]所述的兩種方法；文獻(xiàn)[13]采用主成分分析法和線性最小二乘法確定擬合初值，再采用非線性最小二乘法計(jì)算圓柱模型參數(shù)；文獻(xiàn)[14]采用遍歷算法通過三點(diǎn)共線法確定圓柱軸線向量初始值，再進(jìn)行圓柱面參數(shù)的精確求解。從以上算法可以看出，擬合圓柱面一般需采用合適的方法估算圓柱面擬合參數(shù)初始值，然后再進(jìn)行參數(shù)的精確求解。本文提出了一種基于Levenberg-Marquardt算法(以下簡稱LM算法)的新方法，該方法不用估算初始值，基于指定的初始值就可以進(jìn)行圓柱參數(shù)的計(jì)算。經(jīng)數(shù)據(jù)模擬及案例驗(yàn)證，取得了良好效果。

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

確定一個圓柱面一般需要圓柱面的中軸線及圓柱的半徑，而確定圓柱面的中軸線需要知道中軸線上兩點(diǎn)坐標(biāo)或中軸線方向向量和軸線上的一個點(diǎn)坐標(biāo)等6個參數(shù)，故擬合一個圓柱面需要求解圓柱面中軸線及半徑等7個參數(shù)。設(shè)空間圓柱面其軸線方向向量為(a,b,c)，軸線上一點(diǎn)坐標(biāo)為(X0，Y0，Z0)，圓柱面上測量數(shù)據(jù)為(Xi，Yi，Zi)。

1.1 圓柱面數(shù)學(xué)模型

圓柱面可以用下式表示：

(1)

其中，

A=c(Yi-Y0)-b(Zi-Z0)

B=a(Zi-Z0)-c(Xi-X0)

C=b(Xi-X0)-a(Yi-Y0)

為便于計(jì)算，令(a,b,c)為單位向量，即：a2+b2+c2=1。

(2)

上述問題為約束優(yōu)化問題，求取最佳參數(shù)較繁瑣。為簡化求解，可將約束優(yōu)化問題轉(zhuǎn)化為無約束優(yōu)化問題。參考球面坐標(biāo)系表示方法，令：

(3)

此外，圓柱面中軸線上有無數(shù)點(diǎn)，其中軸線上一點(diǎn)可以任意指定，故其三維坐標(biāo)的某一坐標(biāo)可以取任意值，為便于說明，本文中指定Z0為常數(shù)。

綜上，擬合空間圓柱面原需要計(jì)算7個參數(shù)(a,b,c,X0，Y0，Z0,R)，經(jīng)優(yōu)化后，只需計(jì)算5個參數(shù)，分別為(α,β,X0,Y0,R)，且將問題轉(zhuǎn)化為以下無條件極值問題：

(4)

1.2 參數(shù)求解算法

式(4)為非線性最小二乘問題，根據(jù)文獻(xiàn)[15]，解決非線性最小二乘問題常用算法有最速下降法、Gauss-Newton算法(GN算法)及LM算法。最速下降法以負(fù)梯度方向作為極小化算法的下降方向，是無約束優(yōu)化中最簡單的方法。研究發(fā)現(xiàn)，最速下降方向僅是算法的局部性質(zhì)，對于許多問題，最速下降法并非“最速下降”，而是下降緩慢，特別當(dāng)靠近極小值時，下降變得非常緩慢。GN算法使用泰勒級數(shù)一階展開式去近似地代替非線性回歸模型，然后通過多次迭代，多次修正待求參數(shù)值，直至限差達(dá)到閾值時終止。GN算法不一定總體收斂，其對于小殘量問題有較快的局部收斂速度，對于不是很嚴(yán)重的大殘量問題有較慢的局部收斂速度，對于殘量很大的問題或非線性程度很大的問題則往往不收斂。LM算法最早由Levenberg[16]提出，在GN算法加入正則化參數(shù)改善矩陣條件，后由Marquardt[17]進(jìn)行了改進(jìn)。該方法是一種介于梯度下降法和GN算法之間的算法，常用于求解非線性最小二乘問題，在正則化參數(shù)大時相當(dāng)于梯度下降法，參數(shù)小時相當(dāng)于GN算法。在GN算法中，通常要求雅可比(Jacobi)矩陣是滿秩的，一旦雅可比矩陣奇異的情況發(fā)生，使得算法常常收斂到一個非駐點(diǎn)，LM算法顯著地改善了這一情況。在一定條件下，LM算法具有全局收斂性，其具有超線性收斂速度。在解決非線性最小二乘問題上，LM法是最常用的一種算法，故本文選用LM法求解圓柱面擬合參數(shù)，其主要分為數(shù)據(jù)重心化處理、構(gòu)造函數(shù)及求一階偏導(dǎo)數(shù)、指定初值及迭代求解等步驟。

1.2.1 數(shù)據(jù)重心化處理

為降低數(shù)據(jù)量級及減少迭代過程中可能出現(xiàn)的矩陣奇異情況，將圓柱面上的數(shù)據(jù)進(jìn)行重心化處理。重心化公式如下：

(5)

重心化處理后，圓柱面軸線向量及半徑不變，但圓柱面中軸線位置發(fā)生了改變，故中軸線上一點(diǎn)變?yōu)?x0,y0,z0)，其與重心化處理前圓柱面中軸線上的一點(diǎn)(X0，Y0，Z0)關(guān)系為：

(6)

1.2.2 函數(shù)構(gòu)造及求導(dǎo)

以重心化后的圓柱面上的點(diǎn)(xi,yi,zi)到圓柱中軸線的距離與圓柱面半徑之差構(gòu)造函數(shù)：

(7)

其中,

A′=cosβ(yi-y0)-sinαsinβ(zi-z0)

B′=cosαsinβ(zi-z0)-cosβ(xi-x0)

C′=sinαsinβ(xi-x0)-cosαsinβ(yi-y0)

W=[α,β,x0,y0,R]T

因圓柱面中軸線上有無數(shù)點(diǎn)，其上的一點(diǎn)可隨意指定，為方便計(jì)算，本文指定z0=0。

令：r(W)=[r1(W)r2(W) …rn(W)]T

對(α,β,x0,y0,R)分別求偏導(dǎo)數(shù)，構(gòu)造雅克比矩陣：

(8)

1.2.3 算法步驟

LM算法如下：

步驟1：選取參數(shù)μ0(一般取0.001)，初始點(diǎn)W0∈Rn，容許誤差0≤ε?1，令k:=0；

步驟2：令Jk=D(Wk)，計(jì)算g(Wk)=Jkr(Wk)，若‖g(Wk)‖≤ε，停止迭代，輸出Wk作為近似極小點(diǎn)；

1.2.4 初始值及終值

(9)

(10)

1.3 結(jié)果評價(jià)

圓柱面擬合結(jié)果主要通過計(jì)算圓柱面圓度及根據(jù)圓柱面圓度計(jì)算出的均方根誤差(Root Mean Squared Error, RMSE)來評價(jià)?？臻g圓柱面的圓度ΔRi為圓柱面表面點(diǎn)到圓柱中軸線的距離與圓柱面半徑之差，可以通過式(7)計(jì)算。RMSE通過下式計(jì)算：

(11)

2 案例分析

2.1 案例1

以激光跟蹤儀測量的某設(shè)備圓柱面為例，進(jìn)行圓柱面的擬合，觀測的設(shè)備圓柱面三維坐標(biāo)數(shù)據(jù)如表1所示。

表1 觀測點(diǎn)三維坐標(biāo)/mm

采用本文LM算法計(jì)算圓柱面擬合參數(shù)，容許誤差ε=10-6。經(jīng)迭代12次后，算法收斂至迭代終止。迭代結(jié)果及根據(jù)迭代結(jié)果計(jì)算出的圓柱面參數(shù)如表2所示。

表2 迭代結(jié)果和圓柱面參數(shù)

根據(jù)迭代結(jié)果計(jì)算出圓度值，圓度值統(tǒng)計(jì)如表3所示。

表3 圓度值統(tǒng)計(jì) /mm

為驗(yàn)證本文算法，采用了行業(yè)內(nèi)著名的SA(Spatial Analyzer)測量軟件進(jìn)行計(jì)算，計(jì)算結(jié)果與本文一致。

2.2 案例2

采用文獻(xiàn)[4]案例中數(shù)據(jù)，根據(jù)本文算法計(jì)算圓柱面擬合參數(shù)，容許誤差ε=10-6，并與文獻(xiàn)[4]和文獻(xiàn)[14]計(jì)算結(jié)果進(jìn)行對比，相關(guān)結(jié)果如表4所示。

表4 文獻(xiàn)[4]、文獻(xiàn)[14] 算法和本文算法結(jié)果

經(jīng)過7次迭代，算法迭代終止。因圓柱面中軸線上一點(diǎn)可以隨意指定，為便于與文獻(xiàn)[4]和文獻(xiàn)[14]中數(shù)據(jù)對比，將本文計(jì)算結(jié)果中軸線上一點(diǎn)的X值移動至與文獻(xiàn)[4]和文獻(xiàn)[14]一致，得到的中軸線上的一點(diǎn)變?yōu)?10, 22.649 6, 4.484 3)。對比本文算法和文獻(xiàn)[14]算法可以發(fā)現(xiàn)，本文算法計(jì)算結(jié)果和文獻(xiàn)[14]中計(jì)算結(jié)果一致，說明兩種算法精度相當(dāng)。另對比本文算法和文獻(xiàn)[4]算法可以看出，本文算法迭代較快，RMSE值較小，這表明本文算法略優(yōu)于文獻(xiàn)[4]算法。

3 工程應(yīng)用

AP1000主泵為核1級設(shè)備，是垂直安裝、無密封的單級高慣性離心屏蔽泵，其為核島反應(yīng)堆廠房中核心主設(shè)備之一。每臺機(jī)組含4 臺主泵，每臺主泵可拆卸組件凈重約67.4 t，最大外形尺寸為φ2 044.8 mm×5 874.3 mm。在主泵安裝過程中，需要先將吸入適配器頂升進(jìn)入蒸汽發(fā)生器底部泵殼內(nèi)，再將主泵頂升進(jìn)入吸入適配器中。主泵安裝過程中精度要求很高，如主泵擴(kuò)散器與吸入適配器設(shè)計(jì)間隙為0.38～0.64 mm，相應(yīng)地要求測量精度優(yōu)于0.3 mm。

主泵安裝過程中采用了三維整體建模技術(shù)，該技術(shù)先采用激光跟蹤儀采集主泵外形數(shù)據(jù)以及與主泵匹配的泵殼及吸入適配器外形數(shù)據(jù)。以主泵測量所用的美國自動精密工程公司(Automatic Precision Industry, API) T3激光跟蹤測量系統(tǒng)為例，其角度測量分辨率為0.05″，干涉法距離測量的分辨率為0.1 μm，測量精度1 μm/m；絕對測距精度的分辨率為1 μm，在10 m之內(nèi)測量精度可以達(dá)到15 μm?？紤]到激光跟蹤儀測量精度很高，且主泵所在房間空間狹小，測量距離不超過5 m，根據(jù)文獻(xiàn)[18]研究，API公司T3激光跟蹤儀全量程點(diǎn)位精度為5 μm/m，整個測量過程實(shí)測數(shù)據(jù)的精度不會超過25 μm，與模型誤差相比，實(shí)測數(shù)據(jù)的誤差可忽略不計(jì)，因此，對主泵圓柱面建模主要考慮建模誤差，不考慮實(shí)測數(shù)據(jù)誤差。主泵等設(shè)備的外形和端面均為圓柱形和圓形，因此，在采集外形數(shù)據(jù)后需要進(jìn)行圓柱面和空間圓的建模。對圓柱面的建模采用本文方法進(jìn)行。根據(jù)測量數(shù)據(jù)，分別建立吸入適配器和泵殼組合模型如圖1所示、主泵模型如圖2所示。

圖1 吸入適配器和泵殼組合模型

圖2 主泵建模模型

建模完成后，將主泵模型與泵殼和吸入適配器組合模型進(jìn)行虛擬裝配，計(jì)算得出吸入適配器與主泵擴(kuò)散器之間間隙及位置偏差等，并判斷是否滿足設(shè)計(jì)要求。如不滿足設(shè)計(jì)要求，可根據(jù)偏差數(shù)據(jù)進(jìn)行調(diào)整。位置調(diào)整可分階段多次進(jìn)行，同時在主泵頂升過程中采用激光跟蹤儀監(jiān)測實(shí)時間隙。通過精密建模及實(shí)時調(diào)整，順利完成了核島主泵的精密安裝，安裝過程中未發(fā)生任何碰撞，各項(xiàng)偏差均滿足設(shè)計(jì)要求，這也驗(yàn)證了本文建模方法的準(zhǔn)確性。

4 結(jié) 論

本文基于LM算法，探討了一種圓柱面擬合方法，該方法參考了球坐標(biāo)系表示方式，將約束優(yōu)化問題變成了無約束優(yōu)化問題，方便了參數(shù)計(jì)算。該方法不用計(jì)算初始值，采用指定初始值，提高了計(jì)算效率。經(jīng)案例及工程驗(yàn)證，該方法精度高、收斂快，可以用于高精度的圓柱面擬合。