聞洪峰　和　會

1　武漢大學(xué)測繪學(xué)院，武漢市珞喻路129號，430079　2　河北省第一測繪院，石家莊市中山東路495號，050032　3　河北省基礎(chǔ)地理信息中心，石家莊市中山東路495號，050032

?

基于Levenberg-Marquarat算法的非線性三維直角坐標(biāo)轉(zhuǎn)換方法

聞洪峰1,2和會3

1武漢大學(xué)測繪學(xué)院，武漢市珞喻路129號，4300792河北省第一測繪院，石家莊市中山東路495號，0500323河北省基礎(chǔ)地理信息中心，石家莊市中山東路495號，050032

摘要：提出一種基于Levenberg-Marquarat算法的非線性三維直角坐標(biāo)轉(zhuǎn)換方法,在法矩陣病態(tài)或者奇異時依然有效，并通過修正旋轉(zhuǎn)角參數(shù)的方法，有效解決了平移量與旋轉(zhuǎn)角量綱不同造成的迭代發(fā)散問題。設(shè)計出簡潔有效的迭代求解模式,獲得了穩(wěn)定的參數(shù)解。最后通過模擬數(shù)據(jù)對比分析,證明該方法的有效性和正確性。

關(guān)鍵詞：Levenberg-Marquarat算法;坐標(biāo)轉(zhuǎn)換;七參數(shù);病態(tài);奇異

傳統(tǒng)的三維直角坐標(biāo)轉(zhuǎn)換通常采用線性化的七參數(shù)模型，該模型只適用于小旋轉(zhuǎn)角，不適合于大旋轉(zhuǎn)角的三維直角坐標(biāo)轉(zhuǎn)換問題[1-4]。許多學(xué)者對大旋轉(zhuǎn)角的三維直角坐標(biāo)轉(zhuǎn)換問題進(jìn)行研究。文獻(xiàn)[1]分析局部區(qū)域三維七參數(shù)轉(zhuǎn)換中平移參數(shù)與旋轉(zhuǎn)參數(shù)的強(qiáng)相關(guān)性，提出通過正則化方法改善參數(shù)解的精度；文獻(xiàn)[2]提出基于改進(jìn)的高斯-牛頓法的非線性三維直角坐標(biāo)轉(zhuǎn)換方法；文獻(xiàn)[3]在文獻(xiàn)[2]的基礎(chǔ)上通過對坐標(biāo)進(jìn)行中心化，省卻了計算搜索步長的過程；文獻(xiàn)[4]通過變更求解參數(shù)的方式解決了7參數(shù)量綱差異的問題，提高了收斂性和解的穩(wěn)定性。文獻(xiàn)[1]提出七參數(shù)轉(zhuǎn)換中法方程病態(tài)的一種處理思路，但卻是針對傳統(tǒng)的七參數(shù)線性模型的。文獻(xiàn)[2-4]都是針對七參數(shù)非線性模型的，但都要求法矩陣為良態(tài)，在法矩陣呈病態(tài)或奇異時計算結(jié)果將很不穩(wěn)定。本文提出一種基于Levenberg-Marquarat算法的非線性三維直角坐標(biāo)轉(zhuǎn)換方法。該方法不僅適合小旋轉(zhuǎn)角，也適用于大旋轉(zhuǎn)角的三維坐標(biāo)轉(zhuǎn)換，在法矩陣病態(tài)或奇異時仍能正確解算，并通過旋轉(zhuǎn)角修正保證解的收斂性及穩(wěn)定性。

1　Levenberg-Marquarat算法

Levenberg-Marquardt算法[5]既具有最速下降法的收斂速度，又克服了Gauss-Newton法對初始值的苛刻要求，同時在法方程病態(tài)或奇異時也仍然有效。

設(shè)性能優(yōu)化指數(shù)F(x)的牛頓方法是：

(1)

式中，A(X(k))=2F(X(k)),g(X(k))=F(X(k))，上標(biāo)表示迭代次數(shù)，X為輸入量(待求參數(shù))向量矩陣，x為向量矩陣的元素。

如果設(shè)F(x)是殘差平方函數(shù)之和，即

(2)

(3)

其中，

N為輸出量個數(shù)，n為輸入量個數(shù)。

赫森矩陣可表示為：

(4)

(5)

將式(3)和式(5)代入式(1)，可以得到高斯-牛頓算法的迭代算式：

(6)

問題是，H=JT(X(k))J(X(k))可能不可逆，這時可以采用改進(jìn)的近似赫森矩陣：

(7)

設(shè)H的特征值和特征向量分別為{λ1,λ2,…,λn}和{z1,z2,…,zn}，則有：

(8)

因此，G與H的特征向量相同，且G的特征值為λi+μ,對所有的λi，增加μ以保證λi+μ>0，則G必定正定可逆。從而得到LM算法的迭代算式：

(9)

2　數(shù)學(xué)模型

三維直角坐標(biāo)轉(zhuǎn)換模型[2]可表示為：

(10)

式中，

式(10)可表示為如下形式：

(11)

根據(jù)Levenberg-Marquardt算法，針對每個坐標(biāo)轉(zhuǎn)換用的公共點，可列出如下迭代公式：

(12)

從參數(shù)初值X0出發(fā)，按照上式估計出近似值X(k),然后以該近似值作為下一次迭代的初值計算X(k+1)。反復(fù)迭代，逐次逼近真正的極小點。

通常，μ的初值μ0可取(0,1)之間的任意小數(shù)，然后根據(jù)V(k)和V(k+1)來調(diào)整μ(k+1)，即

(13)

一般情況下,μ會隨著改正量的減小而逐步減小，有時會減小至非常趨近于0的小數(shù)，這時如果式(12)中方程系數(shù)陣JT(X(k))J(X(k))病態(tài)或奇異，則μ不能起到很好的修正作用。因此，在迭代過程中，如果μ小于某一閹值τ0(人為設(shè)定μ的最小值，本文取0.000 000 01)，應(yīng)該將其固定在閾值處，這樣既能保證系數(shù)矩陣可逆，同時由于μ值非常小，在系數(shù)矩陣非奇異的情況下也幾乎不會對最優(yōu)解產(chǎn)生影響。

另外，由于三角函數(shù)具有周期性，而平移參數(shù)的平移量又可能非常大，在局部區(qū)域，二者具有較強(qiáng)的相關(guān)性。在迭代過程中,3個旋轉(zhuǎn)角有可能會跟隨平移量逐漸增大，這時如果不加以處理，旋轉(zhuǎn)參數(shù)有時會由于改正過于劇烈而無法正確收斂。因此，筆者在迭代過程中不斷對旋轉(zhuǎn)參數(shù)進(jìn)行修正，即每次迭代改正后，將旋轉(zhuǎn)參數(shù)加減360°的整數(shù)倍，將其調(diào)整至0°～360°之間。

基于Levenberg-Marquarat算法的非線性三維直角坐標(biāo)轉(zhuǎn)換方法步驟如下：

2)根據(jù)式(12)列出誤差方程，若有n個基準(zhǔn)點,可組成3n個誤差方程，并計算參數(shù)X(k+1)；

3)根據(jù)式(11)計算VTV(k)和VTV(k+1)，如果VTV(k)和VTV(k+1)的差值絕對值小于給定的閾值ξ0(表征期望的收斂精度，本文取ξ0=0.000 000 001),則轉(zhuǎn)至步驟8)；

4)對X(k+1)中旋轉(zhuǎn)參數(shù)進(jìn)行修正，如果εX<0，則εX=εX+(int(εX/360)+1)×360；如果εX≥0，則εX=εX-int(εX/360)×360，εY、εZ作相同處理；

5)根據(jù)式(13)計算μ(k+1)；

6) 判斷μ(k+1)是否小于給定閾值τ0(本文取τ0=0.000 000 01)，若小于閾值τ0，則取μ(k+1)=τ0;

7)轉(zhuǎn)至步驟2)；

8) 輸出最終的轉(zhuǎn)換參數(shù)Xk+1,迭代結(jié)束。

3　算例和分析

算例1模擬算例來自文獻(xiàn)[2]和文獻(xiàn)[3]，變換前后的坐標(biāo)如表1。

根據(jù)本文推導(dǎo)的基于Levenberg-Marquarat算法的七參數(shù)嚴(yán)密解算公式和計算步驟,在VisualC++6.0中編程,計算結(jié)果如下:旋轉(zhuǎn)角εX=29°59′59.999 941″,εY=44°59′59.999 996″,εZ=60°00′00.000 068″;尺度比例 K=1.000 000 000 1;平移參數(shù)ΔX=500.000 0m,ΔY=1 000.000 0m,ΔZ=2 000.000 0m。結(jié)果與實際采用的變換參數(shù)相一致[2-3]，表明本文的方法是有效可靠的。

算例2模擬算例來自文獻(xiàn)[4]，變換前后的坐標(biāo)如表2。

表1　算例1示例數(shù)據(jù)

表2　算例2示例數(shù)據(jù)

根據(jù)本文推導(dǎo)的基于Levenberg-Marquarat算法的七參數(shù)嚴(yán)密解算公式和計算步驟,取2、3、4、5號點為已知點，1號點為檢查點，轉(zhuǎn)換參數(shù)計算結(jié)果如下:旋轉(zhuǎn)角εX=142°04′11.790 503″,εY=229°03′03.348 73″,εZ=12°54′30.581 561″;尺度比例K=0.999 923 072 0;平移參數(shù)ΔX=-2 831 698.766 1 m,ΔY=4 675 677.770 0 m,ΔZ=3 275 369.428 0 m。已知點轉(zhuǎn)換殘差及檢核點殘差如表3所示。

由表3的殘差統(tǒng)計可見，轉(zhuǎn)換精度在mm級，與文獻(xiàn)[4]的轉(zhuǎn)換精度相當(dāng)。而對于算例2，如果采用高斯-牛頓方法，則迭代發(fā)散、無法收斂。

表3　算例2殘差統(tǒng)計

4　結(jié)　語

本文提出的基于Levenberg-Marquarat算法的非線性三維直角坐標(biāo)轉(zhuǎn)換方法,不需要中心化、正則化、計算搜索步長及變更求解參數(shù)，相比文獻(xiàn)[2-4]的算法更加簡單、易行。由于Levenberg-Marquarat算法的特性，本文方法屬于全局收斂，即便系數(shù)矩陣病態(tài)或奇異，仍能穩(wěn)定收斂并得到可靠結(jié)果。由于旋轉(zhuǎn)角與平移量量綱不同，二者在數(shù)值上相差很大時可能會造成不收斂現(xiàn)象，通過迭代過程中對旋轉(zhuǎn)角的不斷修正，可以有效保證收斂性及解的穩(wěn)定性。

參考文獻(xiàn)

[1]沈云中,胡雷鳴,李博峰. Bursa模型用于局部區(qū)域坐標(biāo)變換的病態(tài)問題及其解法[J].測繪學(xué)報, 2006, 35 (2):95-98(Shen Yunzhong, Hu Leiming,Li Bofeng. Ill-Posed Problem in Determination of Coordinate Transformation Parameters with Small Area’s Data Based on Bursa Model [J].Acta Geodaetica et Cartographica Sinica,2006,35(2):95-98)

[2]羅長林,張正祿,鄧勇，等.基于改進(jìn)的高斯-牛頓法的非線性三維直角坐標(biāo)轉(zhuǎn)換方法研究[J].大地測量與地球動力學(xué),2007,27(1): 50-53(Luo Changlin,Zhang Zhenglu,Deng Yong,et al.On Nonlinear 3D Rectangular Coordinate Transformation Method Based on Improved Guss-Newton Method [J]. Journal of Geodesy and Geodynamics,2007,27(1):50-53)

[3]劉仁釗,段文貴，張玉堂.Bursa轉(zhuǎn)換模型七參數(shù)嚴(yán)密解算方法研究[J].資源環(huán)境與工程，2010,24(4)：416-423(Liu Renzhao,Duan Wengui,Zhang Yutang. Study on Close Computation of Parameters of Bursa Transformation Model[J].Resources Environment & Engineering, 2010,24(4)：416-423)

[4]李金嶺,劉鸝,喬書波,等. 關(guān)于三維直角坐標(biāo)七參數(shù)轉(zhuǎn)換模型求解的討論[J].測繪科學(xué)，2010,35(4)：76-78(Li Jinling,Liu Li,Qiao Shubo,et al. Discussion on the Determination of Transformation Parameters of 3D Cartesian Coordinates[J]. Science of Surveying and Mapping,2010,35(4)：76-78)

[5]王解先.七參數(shù)轉(zhuǎn)換中參數(shù)之間的相關(guān)性[J].大地測量與地球動力學(xué),2007,27(2):43-46(Wang Jiexian. Correlations Among Parameters in Seven-Parameter Transformation Model [J].Journal of Geodesy and Geodynamics, 2007,27(2):43-46)

[6]姚吉利.三維坐標(biāo)轉(zhuǎn)換參數(shù)直接計算的嚴(yán)密公式[J].測繪通報, 2006(5):7-10(Yao Jili.Rigorous Formula for Direct Calculating Parameter in 3D Transformation[J].Bulletin of Surveying and Mapping, 2006(5):7-10)

[7]Awange L J, Grafarend E W. Closed Form Solution of the Over Determined Nonlinear 7 Parameter Datum Transformation [J]. AVN, 2003(4): 130-148

[8]李瀟,尹暉.基于最小二乘配置的三維空間坐標(biāo)轉(zhuǎn)換[J]. 測繪工程, 2008,17(2):16-18(Li Xiao,Yin Hui. Three-Dimensional Spatial Coordinate Transformation Based on Least Squares Collocation[J].Engineering of Surveying and Mapping, 2008,17(2):16-18)

[9]陳宇,白征東,羅騰.基于改進(jìn)的布爾沙模型的坐標(biāo)轉(zhuǎn)換方法[J]. 大地測量與地球動力學(xué),2010,30(3):71-73(Chen Yu,Bai Zhengdong,Luo Teng.An Improved Bursa Model for Coordinate Transformation[J]. Journal of Geodesy and Geodynamics, 2010,30(3):71-73)

[10]沈云中,衛(wèi)剛.利用過渡坐標(biāo)系改進(jìn)3維坐標(biāo)變換模型[J]. 測 繪 學(xué) 報,1998,27(2):161-165(Shen Yunzhong,Wei Gang.Improvement of Three Dimensional Coordinate Transformation Model by Use of Interim Coordinate System[J].Acta Geodaetica et Cartographica Sinica,1998,27(2):161-165)

[11]陳義, 沈云中, 劉大杰. 適用于大旋轉(zhuǎn)角的三維基準(zhǔn)轉(zhuǎn)換的一種簡便模型[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2004,29(12):1 102-1 104(Chen Yi,Shen Yunzhong,Liu Dajie.A Simplified Model of Three Dimensional-Datum Transformation Adapted to Big Rotation Angle[J].Geomatics and Information Science of Wuhan University, 2004,29(12):1 102-1 104)

[12]胡亞江,楊曉梅,沙月進(jìn).大歐拉角的空間直角坐標(biāo)轉(zhuǎn)換方法探討[J].現(xiàn)代測繪,2006,29(6):10-12(Hu Yajiang,Yang Xiaomei,Sha Yuejin. Discussion on the Way of Space Right-Angle Coordinate Transform of the Big Euler Angle[J].Modern Surveying and Mapping, 2006,29(6):10-12)

[13]Grafarend E W, Awange L J. Nonlinear Analysis of the Three-dimensional Datum Transformation[J]. Journal of Geodesy, 2003,77:66-76

[14]王建梅,覃文忠. 基于LM算法的BP神經(jīng)網(wǎng)絡(luò)分類器[J]. 武漢大學(xué)學(xué)報:信息科學(xué)版, 2005, 30(10):928-931(Wang Jianmei,Qin Wenzhong.BP Neural Network Classifier Based on Levenberg-Marquardt Algorithm[J]. Geomatics and Information Science of Wuhan University, 2005, 30(10):928-931)

[15]趙雙明,郭秋燕,羅研,等. 基于四元數(shù)的三維空間相似變換解算[J].武漢大學(xué)學(xué)報：信息科學(xué)版, 2009,34(10):1 214-1 217(Zhao Shuangming,Guo Qiuyan,Luo Yan,et al.Quaternion-Based 3D Similarity Transformation Algorithm[J]. Geomatics and Information Science of Wuhan University, 2009,34(10):1 214-1 217)

About the first author:WEN Hongfeng, engineer, majors in GPS data process and regional surface subsidence monitoring by InSAR,E-mail：81350890@qq.com.

收稿日期：2015-08-30

第一作者簡介：聞洪峰，工程師，主要研究方向為GPS數(shù)據(jù)處理及InSAR區(qū)域變形監(jiān)測,E-mail：81350890@qq.com。

DOI：10.14075/j.jgg.2016.08.018

文章編號：1671-5942(2016)08-0737-04

中圖分類號：P226

文獻(xiàn)標(biāo)識碼：A

A Nonlinear 3D Rectangular Coordinate Transformation Method Based on Levenberg-Marquarat Algorithm

WENHongfeng1,2HEHui3

1School of Geodesy and Geomatics,Wuhan University，129 Luoyu Road,Wuhan 430079,China2The First Surveying and Mapping Institute of Hebei Province，495 East-Zhongshan Road,Shijiazhuang 050032，China3Basic Geographic Information Center of Hebei Province，495 East-Zhongshan Road, Shijiazhuang 050032，China

Abstract:A nonlinear 3D rectangular coordinate transformation method based on the Levenberg-Marquarat algorithm,which is effective when normal equation is ill-conditioned or singular,is proposed. The problem of divergencein the process of iteration,which arises due to different dimensions between translation amount and rotation angle,is effectively resolved.A neat and effective model of iterative solving is designed and a stable parametric solution is achieved.Finally, the effectiveness and correctness of the method is proven by comparative analysis of simulated data.

Key words:Levenberg-Marquarat algorithm；coordinate transformation;seven-parameter; ill-conditional; singular