方 興,曾文憲,劉經(jīng)南,2,姚宜斌

1.武漢大學(xué)測繪學(xué)院,湖北 武漢 430079;2.武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心,湖北 武漢 430079

三維坐標(biāo)轉(zhuǎn)換的通用整體最小二乘算法

方 興1,曾文憲1,劉經(jīng)南1,2,姚宜斌1

1.武漢大學(xué)測繪學(xué)院,湖北 武漢 430079;2.武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心,湖北 武漢 430079

三維坐標(biāo)轉(zhuǎn)換模型屬于非線性EIV(errors-in-variables)模型,現(xiàn)有整體最小二乘算法均設(shè)定了某些特殊的假設(shè)條件,如僅適用于小角度或者屬于非統(tǒng)計(jì)意義上的數(shù)值解,并且不能用于結(jié)構(gòu)性的系數(shù)矩陣等,算法適用性受到極大限制。本文提出三維坐標(biāo)轉(zhuǎn)換模型的通用加權(quán)整體最小二乘算法,該算法適用于任意旋轉(zhuǎn)角度以及一般性權(quán)矩陣的三維坐標(biāo)轉(zhuǎn)換,并且將結(jié)構(gòu)性系數(shù)矩陣,同時(shí)包含隨機(jī)和非隨機(jī)元素的系數(shù)矩陣等情況納入統(tǒng)一的算法。實(shí)例計(jì)算表明,本文提出的算法具有通用性,適用于實(shí)際應(yīng)用中的各類三維坐標(biāo)轉(zhuǎn)換模型。

整體最小二乘估計(jì);三維坐標(biāo)轉(zhuǎn)換;變量含誤差模型;非線性算法

1 引 言

三維坐標(biāo)轉(zhuǎn)換模型是大地測量、攝影測量、地圖投影、計(jì)算機(jī)視覺等專業(yè)領(lǐng)域常用的變量含誤差(errors-in-variables,EIV)模型之一,模型描述了源坐標(biāo)系統(tǒng)[X Y Z]至目標(biāo)坐標(biāo)系統(tǒng)[x y z]的函數(shù)映射關(guān)系[1]

式中,n表示坐標(biāo)點(diǎn)個(gè)數(shù);1n表示n×1的元素為1的向量;

其中,M表示坐標(biāo)旋轉(zhuǎn)矩陣;α1、α2和α3表示坐標(biāo)系繞x、y和z軸旋轉(zhuǎn)角度參數(shù);M1、M2和M3表示相應(yīng)的旋轉(zhuǎn)矩陣;[ΔxΔyΔz]T表示平移參數(shù)向量;μ表示尺度參數(shù);式(1)為七參數(shù)三維坐標(biāo)轉(zhuǎn)換模型。

坐標(biāo)轉(zhuǎn)換模型的整體最小二乘估計(jì)(total least squares method,TLS)具有漸進(jìn)無偏性[2]?？梢钥吹?模型(1)的系數(shù)矩陣為[A 1n],包含隨機(jī)量與非隨機(jī)量,并且向量化后的系數(shù)矩陣中的隨機(jī)量呈結(jié)構(gòu)性特征[3];其次,旋轉(zhuǎn)參數(shù)[α1α2α3]與系數(shù)矩陣呈復(fù)雜的非線性關(guān)系。因此,三維坐標(biāo)轉(zhuǎn)換模型無法直接采用TLS算法求解。文獻(xiàn)[4]將式(1)轉(zhuǎn)換為線性化的Gauss-Helmert模型(Gauss-Helmert model,GHM)求解,由于算法的迭代更新項(xiàng)存在問題[5],得到的并不是最優(yōu)的TLS解。文獻(xiàn)[6]提出了一種僅適用于小旋轉(zhuǎn)角度的TLS解。文獻(xiàn)[7]基于Gauss-Markov模型和EIV模型的正交回歸解的等價(jià)性提出了加權(quán)TLS算法,但對權(quán)陣有嚴(yán)格的限制條件。文獻(xiàn)[1,8—9]采用迭代GHM方法推導(dǎo)了一般權(quán)陣條件下的加權(quán)TLS算法,迭代GHM方法盡管被認(rèn)為是解決非線性問題的有效方法,但到目前為止,該方法并沒有從理論上予以證明[10]。近年來,國內(nèi)大地測量領(lǐng)域?qū)θS坐標(biāo)轉(zhuǎn)換模型的TLS算法展開了研究[11-13]。文獻(xiàn)[14]討論了小角度、特殊權(quán)矩陣條件下的基于奇異值分解的TLS算法。文獻(xiàn)[15]在單位權(quán)陣條件下,將旋轉(zhuǎn)矩陣全部9個(gè)元素視為待估參數(shù),構(gòu)造了附有限制條件的間接平差模型,采用奇異值分解方法,得到了任意旋轉(zhuǎn)角度的TLS解法。由于基于奇異值分解的TLS算法本質(zhì)上屬于數(shù)值解法,無法得到統(tǒng)計(jì)意義上的估計(jì)值,在測量中的應(yīng)用受到了限制[16-17]。文獻(xiàn)[18]通過將式(1)轉(zhuǎn)換為附有限制條件的平差模型導(dǎo)出了特殊權(quán)矩陣條件下的加權(quán)TLS算法。

可以看到,現(xiàn)有三維坐標(biāo)轉(zhuǎn)換模型的TLS算法均受到特殊條件限制,實(shí)際應(yīng)用范圍有限,如攝影測量中物方坐標(biāo)系到像方坐標(biāo)系的轉(zhuǎn)換是大角度三維坐標(biāo)系統(tǒng)之間的轉(zhuǎn)換,不同大地基準(zhǔn)下點(diǎn)的坐標(biāo)測量技術(shù)手段可能不同,從而坐標(biāo)點(diǎn)的精度不等且相關(guān),因此,有必要研究一般情況下的TLS算法。本文通過提取系數(shù)矩陣的隨機(jī)量作為待估計(jì)參數(shù),研究了三維坐標(biāo)轉(zhuǎn)換模型的通用加權(quán)整體最小二乘算法,該算法適用于任意旋轉(zhuǎn)角度以及一般性的權(quán)矩陣,將結(jié)構(gòu)性系數(shù)矩陣,同時(shí)包含隨機(jī)和非隨機(jī)元素的系數(shù)矩陣等情況,納入統(tǒng)一的坐標(biāo)轉(zhuǎn)換模型算法,能夠很好地滿足現(xiàn)代大地測量領(lǐng)域?qū)Ω呔热S坐標(biāo)轉(zhuǎn)換模型估計(jì)的需要。

2 三維坐標(biāo)轉(zhuǎn)換的加權(quán)整體最小二乘算法

為采用統(tǒng)一的方法求解結(jié)構(gòu)性系數(shù)矩陣的EIV模型,將模型(1)中系數(shù)矩陣的隨機(jī)量與坐標(biāo)轉(zhuǎn)換參數(shù)一并作為待估計(jì)參數(shù),則三維坐標(biāo)轉(zhuǎn)換模型可描述為

式中,ˉA表示將系數(shù)矩陣的隨機(jī)量(即源坐標(biāo))作為模型的待估計(jì)參數(shù)所構(gòu)成的矩陣。

經(jīng)過以上變換,坐標(biāo)轉(zhuǎn)換的EIV模型可表示為非線性GM模型形式

在整體最小二乘準(zhǔn)則下,進(jìn)一步將式(3)的求解轉(zhuǎn)化為經(jīng)典的非線性模型的最優(yōu)化問題。令式(3)中觀測值l的一般性權(quán)陣為W,不限定觀測數(shù)據(jù)(源坐標(biāo)和目標(biāo)坐標(biāo)值)的精度或者相關(guān)性,式(3)的整體最小二乘最優(yōu)估計(jì)可表示為以下最優(yōu)化問題

導(dǎo)出式(4)的梯度即一階偏導(dǎo)項(xiàng)

式中,各項(xiàng)梯度公式分別為

導(dǎo)出三維坐標(biāo)轉(zhuǎn)換模型的梯度公式(5)—式(10)后,根據(jù)非線性最優(yōu)化理論,可以采用多種方法求解模型參數(shù)的最佳估值。如根據(jù)梯度公式,在參數(shù)的近似值處采用梯度下降法,即在反梯度方向搜索模型的最優(yōu)估計(jì)值,或者采用高斯牛頓法求解等。這些方法均為非線性估計(jì)的基本方法,具體請參閱文獻(xiàn)[19]。采用傳統(tǒng)的大地測量求解非線性Gauss-Markov模型的算法同樣可以求解模型(3)[20]。運(yùn)用各種迭代算法求解坐標(biāo)轉(zhuǎn)換參數(shù)時(shí),參數(shù)初始值可選取。以擬牛頓法為例,列出三維坐標(biāo)轉(zhuǎn)換模型式(2)的整體最小二乘算法的計(jì)算步驟如下:

(2)迭代計(jì)算擬牛頓法中的海森(Hessian)矩陣的近似表達(dá)式量,直到前后兩次參數(shù)估計(jì)之差Δβi小于設(shè)定的正微小量

(3)根據(jù)以下Broyden-Fletcher-Goldfarb-Shanno(BFGS)算法公式,得到三維坐標(biāo)轉(zhuǎn)換模型的整體最小二乘估計(jì)結(jié)果

3 實(shí) 例

筆者選用文獻(xiàn)[7]中大角度三維坐標(biāo)轉(zhuǎn)換模型實(shí)例數(shù)據(jù),采用本文提出的算法進(jìn)行計(jì)算。坐標(biāo)轉(zhuǎn)換數(shù)據(jù)來源于兩個(gè)表面模型坐標(biāo)轉(zhuǎn)換的實(shí)測觀測值,共包含4個(gè)公共點(diǎn),相應(yīng)的源坐標(biāo)值和目標(biāo)坐標(biāo)值見表1,模型估計(jì)結(jié)果見表2。

表1 公共點(diǎn)三維坐標(biāo)數(shù)據(jù)[7]Tab.1 Public points 3D coordinates

表2 三維坐標(biāo)轉(zhuǎn)換參數(shù)估計(jì)結(jié)果Tab.2 Results of the 3D datum transformation

當(dāng)坐標(biāo)數(shù)據(jù)等權(quán)即權(quán)陣為單位陣(W＝I24)的情況下,表2的第2列和第3列分別為轉(zhuǎn)換參數(shù)的最小二乘估計(jì)值和整體最小二乘估計(jì)值,第4列為設(shè)定坐標(biāo)數(shù)據(jù)非等權(quán)情況下(W＝I6?diag (1,4,6.25,16))的整體最小二乘估計(jì)值,表中的估計(jì)結(jié)果與文獻(xiàn)[1,7]一致。對比第2列和第3列數(shù)據(jù)可以看到,整體最小二乘估計(jì)結(jié)果由于同時(shí)顧及了目標(biāo)坐標(biāo)和源坐標(biāo)的隨機(jī)誤差,殘差平方和遠(yuǎn)小于最小二乘的殘差平方和,與文獻(xiàn)[21—23]實(shí)例中坐標(biāo)轉(zhuǎn)換模型的最小二乘和整體最小二乘殘差平方和對比結(jié)論一致。同時(shí),表中旋轉(zhuǎn)參數(shù)的最小二乘和整體最小二乘估計(jì)結(jié)果相等,與文獻(xiàn)[7]證明的當(dāng)觀測數(shù)據(jù)的權(quán)陣為一定特殊形式時(shí),三維坐標(biāo)轉(zhuǎn)換模型旋轉(zhuǎn)參數(shù)的最小二乘和整體最小二乘估計(jì)結(jié)果相等的結(jié)論一致。第3列和第4列數(shù)據(jù)表明,權(quán)陣對整體最小二乘估計(jì)結(jié)果影響顯著,若沒有采用正確的權(quán)陣,例如沒有考慮坐標(biāo)數(shù)據(jù)不等精度、相關(guān)性,或者權(quán)陣沒有考慮到系數(shù)矩陣的結(jié)構(gòu)性特點(diǎn)等,都會(huì)導(dǎo)致整體最小二乘的估計(jì)結(jié)果產(chǎn)生偏差。

4 結(jié) 論

三維坐標(biāo)轉(zhuǎn)換模型是大地測量、攝影測量等眾多專業(yè)領(lǐng)域常用的基本EIV模型之一,現(xiàn)代大地測量中高精度、三維動(dòng)態(tài)大地坐標(biāo)系統(tǒng)的建立對三維坐標(biāo)轉(zhuǎn)換模型的估計(jì)算法提出了更高的要求?，F(xiàn)有TLS算法均受到一定限制條件的制約,如僅適用于小角度坐標(biāo)轉(zhuǎn)換模型、屬于非統(tǒng)計(jì)意義上的數(shù)值解以及需要對坐標(biāo)轉(zhuǎn)換模型的結(jié)構(gòu)性系數(shù)矩陣的權(quán)陣進(jìn)行特殊處理等。本文通過將模型中的源坐標(biāo)與坐標(biāo)轉(zhuǎn)換參數(shù)一并作為待估計(jì)參數(shù),推導(dǎo)了三維坐標(biāo)轉(zhuǎn)換模型的TLS通用算法,該算法適用于任意旋轉(zhuǎn)角的坐標(biāo)轉(zhuǎn)換以及一般性的權(quán)矩陣,并可直接應(yīng)用于坐標(biāo)轉(zhuǎn)換模型結(jié)構(gòu)性系數(shù)矩陣的情況。本文通過大角度三維坐標(biāo)轉(zhuǎn)換的實(shí)例對該算法進(jìn)行了說明和驗(yàn)證,結(jié)果表明新算法突破了現(xiàn)有算法的限制條件,能有效地應(yīng)用于坐標(biāo)轉(zhuǎn)換模型的一般情況,其通用性有利于應(yīng)用中的編程實(shí)現(xiàn)。

[1] FANG Xing.Weighted Total Least Squares Solution for Application in Geodesy[D].Hanover:Leibniz University of Hanover,2011.

[2] HUFFEL S V,VANDEWALLE J.The Total Least-squares Problem:Computational Aspects and Analysis[M].Philadelphia:Society for Industrial and Applied Mathematics,1991.

[3] FANG Xing.A Structured and Constrained Total Leastsquares Solution with Cross-covariances[J].Studia Geophysica et Geodaetica,2014,58(1):1-16.

[4] BLEICH P,ILLNER M.Strenge L?sung der R?umlichen Koordinatentransformation Durch Iiterative Berechnung[J].Allgemeine Vermessungs-Nachrichten,1989,96(4):133-144.

[5] POPE A.Some Pitfalls to be Avoided in the Iterative Adjustment of Nonlinear Problems[C]∥Proceedings of the 38th Annual Meeting of American Society Photogrammetry.Washington DC:[s.n.],1972:449-473.

[6] ACAR A,?YLüDEMIR MT,AKYILM AZ O,et al.Deformation Analysis with Total Least Squares[J].Natural Hazards and Earth System Sciences,2006,6: 663-669.

[7] FELUS F,BURTCH R.On Symmetrical Three-dimensional Datum Conversion[J].GPS Solutions,2009,13 (1):65-74.

[8] AKYILMAZ O.Solution of the Heteroscedastic Datum Transformation Problem[R].Munich:International Association of Geodesy,2012.

[9] LU Jue,CHEN Yi,FANG Xing,et al.Performing 3-D Similarity Transformation Using the Weighted Total Least-squares Method[R].Munich:International Association of Geodesy,2012.

[10] NEITZEL F.Generalization of Total Least-squares on Example of Unweighted and Weighted 2D Similarity Transformation[J].Journal of Geodesy,2010,84(12):751-762.

[11] LI Bofeng,SHEN Yunzhong,LI Weixiao.The Seamless Model for Three-dimensional Datum Transformation[J].Science China:Earth Science,2012,55(12):2099-2108.

[12] LI Bofeng,SHEN Yunzhong,ZHANG Xingfu,et al.Seamless Multivariate Affine Error-in-variables Transformation and Its Application to Map Rectification[J].International Journal of Information Science,2013,27(8):1572-1592.[13] CHEN Yi,LU Jue.Performing 3D Similarity Transformation by Robust Total Least Squares[J].Acta Geodaetica et Cartographica Sinica,2012,41(5):715-722.(陳義,陸玨.以三維坐標(biāo)轉(zhuǎn)換為例解算穩(wěn)健總體最小二乘方法[J].測繪學(xué)報(bào),2012,41(5):715-722.)

[14] LU Jue,CHEN Yi,ZHENG Bo.Applying Total Least Squares to the Three-dimensional Datum Transformation [J].Journal of Geodesy and Geodynamics,2008,28(5): 77-81.(陸玨,陳義,鄭波.總體最小二乘方法在三維坐標(biāo)轉(zhuǎn)換中的應(yīng)用[J].大地測量與地球動(dòng)力學(xué),2008,28(5): 77-81.)

[15] XU Chaoqian,YAO Yibin,XIONG Siting,et al.3 D Rectangular Coordinate Transformation Adapted to Arbitrary Rotation Angle Based on Total Least-SquaresRegression[J].Journal of Geomatics,2010,35(5):46-48.(許超鈐,姚宜斌,熊思婷,等.三維任意旋轉(zhuǎn)角度坐標(biāo)轉(zhuǎn)換的整體最小二乘回歸解法[J].測繪信息與工程, 2010,35(5):46-48.)

[16] XU Peiliang,LIU Jingnan,SHI Chuang.Total Least Squares Adjustment in Partial Errors-in-variables Models:Algorithm and Statistical Analysis[J].Journal of Geodesy,2012,86(8): 661-675.

[17] FANG Xing.Weighted Total Least Squares:Necessary and Sufficient Conditions,Fixed and Random Parameters[J].Journal of Geodesy,2013,87(8):733-749.

[18] GE Xuming,WU Jicang.Iterative Method of Weight Constraint Total Least Squares for Three Dimensional Datum Transformation[J].Science of Wuhan University:Geomatics and Information,2012,37(2):178-182.(葛旭明,伍吉倉.三維基準(zhǔn)轉(zhuǎn)換的約束加權(quán)混合整體最小二乘的迭代解法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2012,37(2):178-182.)

[19] NOCEDAL J,WRIGHT S.Numerical Optimization[M].Berlin:Springer,2006.

[20] LENZMANN L,LENZMANN E.Zur L?sung des Nichtlinearen Gauss-Markov-Modells[J].Zeitschrift für Geod?sie, Geoinformation und Landmanagement,2007,132:108-120.

[21] KONG Jian,YAO Yibin,WU Han.Iterative Method for Total Least-squares[J].Science of Wuhan University: Geomatics and Information,2010,35(6):711-714.(孔建,姚宜斌,吳寒.整體最小二乘的迭代解法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2010,35(6):711-714.)

[22] SCHAFFRIN B,FELUS Y A.On the Multivariate Total Least-squares Approach to Empirical Coordinate Transformations:Three Algorithms[J].Journal of Geodesy, 2008,82(6):373-383.

[23] CAI Jianqing,GRAFAREND EW.Systematical Analysis of the Transformation between Gauss-Krueger-coordinate/ DHDN and UTM-coordinate/ETRS89 in Baden-Württemberg with Different Estimation Methods[C]∥Geodetic Reference Frame International Association of Geodesy Symposia.Munich:International Association of Geodesy,2009,134: 205-211.

(責(zé)任編輯:陳品馨)

A General Total Least Squares Algorithm for Three-dimensional Coordinate Transformations

FANG Xing1,ZENG Wenxian1,LIU Jingnan1,2,YAO Yibin1
1.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China;2.Research Center of GNSS, Wuhan University,Wuhan 430079,China

The model of three-dimensional coordinate transformation is a nonlinear errors-in-variables model.The methods proposed in published literatures are always restricted in practice for their special assumptions,such as size of rotation angles,structured design matrix and special weight matrix.A general weighted TLS algorithm of a three-dimensional coordinate transformation is investigated in this paper.The new algorithm can be applied in transformations with an arbitrary rotation angles and any applicable weights of the observations,as well as the structured design matrix or the design matrix with both random and fixed elements.The example given in this paper illustrates that this algorithm is general and suits to all kinds of three-dimensional coordinate transformations in practice.

total least-squares method;three-dimensional coordinate transformations;errors-in-variables model;nonlinear program

FANG Xing(1981—),male,PhD,majors in the theory and method of surveying data processing.E-mail:xfang＠sgg.whu.edu.cn

ZENG Wenxian

P207

A

1001-1595(2014)11-1139-05

國家自然科學(xué)基金(41474006;41404005;41231174;41174012);中央高?；究蒲谢?2042014kf053)

2014-01-13

方興(1981—),男,博士,現(xiàn)主要從事測量數(shù)據(jù)處理理論與應(yīng)用研究。

曾文憲

E-mail:wxzeng＠sgg.whu.edu.cn

FANG Xing,ZENG Wenxian,LIU Jingnan,et al.A General Total Least Squares Algorithm for Three-dimensional Coordinate Transformations[J].Acta Geodaetica et Cartographica Sinica,2014,43(11):1139-1143.(方興,曾文憲,劉經(jīng)南,等.三維坐標(biāo)轉(zhuǎn)換的通用整體最小二乘算法[J].測繪學(xué)報(bào),2014,43(11):1139-1143.)

10.13485/j.cnki.11-2089.2014.0193

修回日期:2014-08-29