過(guò)家春,趙秀俠,吳艷蘭

1.安徽農(nóng)業(yè)大學(xué)理學(xué)院,安徽合肥 230036;2.安徽大學(xué)資源與環(huán)境工程學(xué)院,安徽合肥 230601; 3.安徽省農(nóng)業(yè)科學(xué)院水產(chǎn)研究所,安徽合肥 230051

空間直角坐標(biāo)與大地坐標(biāo)轉(zhuǎn)換的拉格朗日反演方法

過(guò)家春1,2,趙秀俠3,吳艷蘭2

1.安徽農(nóng)業(yè)大學(xué)理學(xué)院,安徽合肥 230036;2.安徽大學(xué)資源與環(huán)境工程學(xué)院,安徽合肥 230601; 3.安徽省農(nóng)業(yè)科學(xué)院水產(chǎn)研究所,安徽合肥 230051

以拉格朗日反演理論為基礎(chǔ),導(dǎo)出空間直角坐標(biāo)向大地坐標(biāo)轉(zhuǎn)換的一種新的直接解法。該方法將歸化緯度的正弦函數(shù)sinμ表達(dá)為以空間直角坐標(biāo)(X,Y,Z)為基礎(chǔ)的相關(guān)變量的多項(xiàng)式。為驗(yàn)證公式的精度水平和實(shí)用性,將WGS-84橢球參數(shù)代入驗(yàn)算,結(jié)果表明,新解法在－2×106m≤H≤＋1010m范圍內(nèi),展開(kāi)至b4項(xiàng)緯度反解誤差不超過(guò)9.71×10－6″,展開(kāi)至b5項(xiàng)的誤差不超過(guò)9.67× 10－8″,有效范圍比Bowring公式廣,在H≥＋105m的區(qū)域精度明顯優(yōu)于Bowring公式。與迭代算法相比,新方法在H≥－2×106m的區(qū)域精度與迭代算法精度相當(dāng),但計(jì)算效率明顯優(yōu)于迭代算法,展開(kāi)至b4項(xiàng)的CPU執(zhí)行時(shí)間分別約為迭代4次的1/2、展開(kāi)至b5項(xiàng)的CPU執(zhí)行時(shí)間約為迭代5次的1/10。兼顧精度和效率,本文算法優(yōu)于Bowring公式和迭代算法。

空間直角坐標(biāo);大地坐標(biāo);坐標(biāo)轉(zhuǎn)換;拉格朗日反演定理

1 引 言

空間直角坐標(biāo)系與大地坐標(biāo)系之間的轉(zhuǎn)換是空間大地測(cè)量技術(shù)中的基本問(wèn)題,亦是近幾十年來(lái)國(guó)內(nèi)外大地測(cè)量學(xué)者討論的熱點(diǎn)問(wèn)題之一[1-7]。眾所周知,由大地坐標(biāo)解算空間直角坐標(biāo)可直接得到,而由空間直角坐標(biāo)反解大地坐標(biāo)則比較復(fù)雜。國(guó)內(nèi)外學(xué)者對(duì)此提出許多不同的解法,總體上可歸納為間接解法和直接解法兩大類。間接解法以各種迭代算法為主[8-13],直接解法則較為多樣,包括構(gòu)造并求解不同形式的一元四次方程[14-17]、各類近似公式、閉合解析解[18-23]等。各類直接解法中,Bowring公式較為著名,但精度不夠理想,若要提高精度,則需進(jìn)一步迭代;而其他各類直接解法,或因公式過(guò)于冗長(zhǎng),或因需要進(jìn)行復(fù)雜的求根討論,以致求解過(guò)程復(fù)雜、推廣困難。因此,對(duì)于這一問(wèn)題,國(guó)內(nèi)外大都仍采用經(jīng)典的迭代算法求解。

本文以拉格朗日反演理論為基礎(chǔ),導(dǎo)出空間直角坐標(biāo)向大地坐標(biāo)轉(zhuǎn)換的泰勒級(jí)數(shù)展開(kāi)新方法。

2 理論基礎(chǔ)及基本思路

2.1 拉格朗日反演定理(Lagrange inversion

theorem)[24-26]

若函數(shù)y＝f(x)在點(diǎn)x0的某一鄰域U(x0)內(nèi)可展開(kāi)為泰勒級(jí)數(shù)

且f′(x0)≠0,則其反函數(shù)x＝f－1(y)在相應(yīng)的y0點(diǎn)的某一鄰域V(y0)內(nèi)也可展開(kāi)為泰勒級(jí)數(shù)

式中,記y0＝f(x0)＝f(0)(x0)。

拉格朗日反演定理表明,在一階導(dǎo)數(shù)不等于零的情況下,解析函數(shù)的反函數(shù)也可展開(kāi)為泰勒級(jí)數(shù)。

2.2 基本思路

旋轉(zhuǎn)橢球下,大地坐標(biāo)(L,B,H)與空間直角坐標(biāo)(X,Y,Z)之間的關(guān)系為

圖1 空間直接坐標(biāo)與大地坐標(biāo)之關(guān)系Fig.1 Cartesian coordinates and geodetic coordinates

圖1中,OI＝Ne2sin B＝be′2sinμ,OQ＝Ne2cos B＝ae2cosμ,C為P0點(diǎn)的曲率中心,坐標(biāo)為(ae2cos3μ,－be′2sin3μ),對(duì)應(yīng)的JK弧為橢圓的漸屈線。以P、C兩點(diǎn)坐標(biāo)求解tan B即為Bowring公式的思路。

由式(4)及圖1關(guān)系,可得

式(8)、式(9)兩式中an、bn的關(guān)系為[25-26]

由歸化緯度與大地緯度的關(guān)系式(5)可得

然后按式(12)求得大地高H

上式橢球面上及橢球外取“＋”,橢球內(nèi)取“－”。

3 系數(shù)求解及簡(jiǎn)化

由上述討論可首先求得式(8)的系數(shù)

再由關(guān)系式(10)可得

該系數(shù)較為復(fù)雜,為簡(jiǎn)化其結(jié)構(gòu),令

綜合以上各式可得

由此即可按式(11)、式(12)解得大地緯度B和大地高H。

以上為使推導(dǎo)過(guò)程方便,引入?yún)?shù)w、z,實(shí)際計(jì)算時(shí),r、t1—t4可由空間直角坐標(biāo)(X,Y,Z)直接得出

4 驗(yàn)算分析

盡管式(17)本質(zhì)上是以拉格朗日反演定理為基礎(chǔ)的泰勒級(jí)數(shù)展開(kāi),但由于該問(wèn)題的復(fù)雜性,按泰勒級(jí)數(shù)的余項(xiàng)理論給出其誤差估計(jì)仍十分困難。為驗(yàn)證公式的正確性及其精度水平,下面以WGS-84橢球參數(shù)代入驗(yàn)算,并與經(jīng)典的直接解法Bowring公式及迭代算法進(jìn)行對(duì)比。

4.1 盲區(qū)對(duì)比

式(17)與許多已有的其他解法一樣,在接近橢球中心的區(qū)域存在盲區(qū)。

首先,由以上討論可知,該解算方法要求(w2＋z2)3/2－e2w2≠0。同時(shí),經(jīng)驗(yàn)算分析,當(dāng)(w2＋z2)3/2－e2w2＜0、z≠0時(shí),式(17)右端的值域超出sinμ的值域范圍。由此可見(jiàn),(w2＋z2)3/2－e2w2≤0、z≠0所含區(qū)域?yàn)樵摲椒ǖ拿^(qū)。

已有算法的盲區(qū)(或多值區(qū))多為橢球的漸屈線附近及其以內(nèi)區(qū)域,如Bowring公式、經(jīng)典迭代算法以及文獻(xiàn)[20]所給的算法,與本文算法的盲區(qū)(w2＋z2)3/2－e2w2≤0、z≠0,即

的區(qū)域關(guān)系如圖2所示。

圖2 盲區(qū)對(duì)比Fig.2 Comparison of the“dead zone”between Bowring’s formula and new method

4.2 精度對(duì)比

以上盲區(qū)范圍距地心距離不超過(guò)42.84 km,處于地核內(nèi)核內(nèi)部,衛(wèi)星大地測(cè)量一般不涉及這一區(qū)域。因此,從實(shí)用的角度看,各國(guó)學(xué)者一般最大范圍取大地高－6×106～＋1010m的范圍進(jìn)行分析,本文亦取此范圍進(jìn)行分析。將本文算法與Bowring公式、經(jīng)典迭代算法對(duì)比,其緯度解算誤差對(duì)比如表1和圖3所示。因高程解算精度主要取決于緯度解算,以下僅對(duì)緯度解算誤差進(jìn)行分析。

4.3 效率對(duì)比

為對(duì)比3種算法的計(jì)算效率,筆者在同一臺(tái)計(jì)算機(jī)上,用Mathematica軟件以緯度1°、大地高108m為步長(zhǎng),分別利用3種方法計(jì)算B∈[0°, 90°],H∈[0,1010]范圍內(nèi)91×101個(gè)點(diǎn)的總CPU執(zhí)行時(shí)間,進(jìn)而得到單個(gè)離散點(diǎn)的平均CPU執(zhí)行時(shí)間,計(jì)算結(jié)果對(duì)比如表2所示。

表1 本文算法與Bowring公式、迭代算法的誤差對(duì)比Tab.1 Errors of Cartesian to geodetic coordinates between Bowring’s formula,iterative algorithm and new method

圖3 本文算法與Bowring公式、迭代算法的誤差對(duì)比Fig.3 Comparison of the errors of geodetic latitude between Bowring’s formula,iterative algorithm and our algorithm

表2 3種算法的CPU執(zhí)行時(shí)間對(duì)比Tab.2 Comparison of the CPU times between Bowring’s formula,iterative algorithm and new method

4.4 綜合分析

通過(guò)以上對(duì)比分析可知:

(1)新算法與兩種經(jīng)典算法在接近橢球中心區(qū)域均存在盲區(qū),需進(jìn)一步研究,但實(shí)際應(yīng)用一般不涉及這一區(qū)域。

(2)精度上,若按照對(duì)應(yīng)的點(diǎn)位誤差不超過(guò)0.1 mm的原則,Bowring公式的有效范圍約為－105m≤H≤＋105m,在H≥＋105m以外區(qū)域解算精度迅速下降,最弱精度位于H＝2a、B＝π/4處,約＋1.84×10－3″,相當(dāng)于經(jīng)線方向的0.11 m的點(diǎn)位誤差。

新解法展開(kāi)至b4項(xiàng)、b5項(xiàng),以及迭代法迭代4次、5次在H≥－2×10－6m的區(qū)域范圍能夠滿足要求,有效范圍比Bowring公式大。

(3)計(jì)算效率上,Bowring公式最高,迭代法計(jì)算效率明顯低于新解法計(jì)算效率。

綜上可知,在－105m≤H≤＋105m的范圍內(nèi),Bowring公式精度最優(yōu),計(jì)算效率上占有顯著優(yōu)勢(shì),對(duì)于地面上的定位與導(dǎo)航應(yīng)用,Bowring公式占據(jù)優(yōu)勢(shì)。

但實(shí)際上導(dǎo)航衛(wèi)星大都是中地球軌道衛(wèi)星,衛(wèi)星高度在107m級(jí),如GPS、GLONASS、北斗等衛(wèi)星導(dǎo)航系統(tǒng)的衛(wèi)星高度分別約為20 200 km、19 100 km、21 500 km,此時(shí)Bowring公式精度明顯較弱,如表1所示。兼顧精度和效率,本文算法較之于Bowring公式和迭代算法有一定優(yōu)勢(shì)。

此外,與其他直接解法相比,本文算法不必進(jìn)行復(fù)雜的討論,公式相對(duì)簡(jiǎn)潔,可根據(jù)精度需要決定末項(xiàng)取舍。

5 結(jié) 論

本文以拉格朗日反演理論為基礎(chǔ),得到空間直角坐標(biāo)向大地坐標(biāo)轉(zhuǎn)換的泰勒級(jí)數(shù)展開(kāi)新方法,是一種新的直接解法,且經(jīng)驗(yàn)算分析表明其精度可靠,能夠滿足實(shí)際應(yīng)用精度要求。新解法由式(17)、式(18)兩式構(gòu)成,與已有直接解法相比,公式相對(duì)簡(jiǎn)潔,無(wú)須復(fù)雜討論,可根據(jù)精度需要決定末項(xiàng)取舍,可擴(kuò)展性好。同時(shí),該解法隨大地高的增加精度提高顯著,這對(duì)現(xiàn)代空間技術(shù)是十分有利的,如航天工程、探月工程等。該方法的缺點(diǎn)是接近橢球中心區(qū)域存在盲區(qū),這亦是Bowring公式、經(jīng)典迭代算法等解法的共同問(wèn)題,需要進(jìn)一步研究。

[1] CHEN Junyong.Direct Solution of Transforming Cartesian to Geodetic Coordinates[C]∥Special Issue on Geodesy.Beijing:Surveying and Mapping Press,1979:67-72.(陳俊勇.空間直角坐標(biāo)與大地坐標(biāo)換算的直接解法[C]∥大地測(cè)量研究專輯.北京:測(cè)繪出版社,1979:67-72.)

[2] ZENG Qixiong.Closed Formula for Direct Solution of Geodetic Coordinates from Rectangular Space Coordinates[J].Acta Geodaetica et Cartographica Sinica,1981,10(2):83-95.(曾啟雄.空間直角坐標(biāo)直接解算大地坐標(biāo)的閉合公式[J].測(cè)繪學(xué)報(bào),1981,10(2):83-95.)

[3] DEAKIN R E,HUNTER M N.Geometric Geodesy:Part A[M].Melbourne:[s.n.],2010:89-102.

[4] KONG Xiangyuan,GUO Jiming,LIU Zongquan.Foundation of Geodesy[M].Wuhan:Wuhan University Press,2001: 36-39.(孔祥元,郭際明,劉宗泉.大地測(cè)量學(xué)基礎(chǔ)[M].武漢:武漢大學(xué)出版社,2001:36-39.)

[5] ZENG Qixiong.Direct Calculation Formula of Gauss Plane Rectangular Coordinates from Space Rectangular Coordinates [J].Acta Geodaetica et Cartographica Sinica,1993, 22(1):74-79.(曾啟雄.空間直角坐標(biāo)直接計(jì)算高斯平面直角坐標(biāo)公式[J].測(cè)繪學(xué)報(bào),1993,22(1):74-79.)

[6] LIU Jingnan,LIU Dajie.The Influence of the Accuracy in Geodetic and Geocentric Coordinates on the Accuracy in the Results of Simultaneous Adjustment[J].Acta Geodaetica et Cartographica Sinica,1985,14(2):133-144.(劉經(jīng)南,劉大杰.大地坐標(biāo)和地心坐標(biāo)精度對(duì)聯(lián)合平差的精度影響[J].測(cè)繪學(xué)報(bào),1985,14(2):133-144.)

[7] LIU Dajie,BAI Zhengdong.A Mathematic Model of the 3-dimensional Adjustment of GPS Network[J].Acta Geodaetica et Cartographica Sinica,1997,26(1):37-41.(劉大杰,白征東.一種GPS網(wǎng)三維平差的數(shù)學(xué)模型[J].測(cè)繪學(xué)報(bào),1997,26(1):37-41.)

[8] BOWRING B R.Transformation from Spatial to Geographical Coordinates[J].Survey Review,1976,181:323-327.

[9] LIN K C,WANG J.Transformation from Geocentric to Geodetic Coordinates Using Newton’s Iteration[J].Bulletin Géodésique,1995,69:300-303.

[10] FUKUSHIMA T.Fast Transform from Geocentric to Geodetic Coordinates[J].Journal of Geodesy,1999,73:603-610.

[11] POLLARD J.Iterative Vector Methods for Computing Geodetic Latitude and Height from Rectangular Coordinates[J].Journal of Geodesy,2002,76:36-40.

[12] FELTENS J.Vector Method to Compute the Cartesian(X, Y,Z)to Geodetic(φ,λ,h)Transformation on a Triaxial Ellipsoid[J].Journal of Geodesy,2009,83:129-137.

[13] SHU Chanfang,LI Fei,SHEN Fei.A New Algorithm for Transforming from Cartesian to Geodetic Coordinates[J].Geomatics and Information Science of Wuhan University, 2009,34(5):561-563.(束蟬方,李斐,沈飛.空間直角坐標(biāo)向大地坐標(biāo)轉(zhuǎn)換的新算法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2009,34(5):561-563.)

[14] BORKOWSKI K M.Accurate Algorithms to Transform Geocentric to Geodetic Coordinates[J].Bulletin Géodésique, 1989,63:50-56.

[15] VERMEILLE H.Direct Transformation from Geocentric to Geodetic Coordinates[J].Journal of Geodesy,2002, 76:451-454.

[16] VERMEILLE H.Computing Geodetic Coordinates from Geocentric Coordinates[J].Journal of Geodesy,2004,78: 94-95.

[17] VERMEILLE H.An Analytical Method to Transform Geocentric into Geodetic Coordinates[J].Journal of Geodesy, 2011,85:105-117.

[18] FEATHERSTONE W E,CLAESSENS S J.Closed-form Transformation between Geodetic and Ellipsoidal Coordinates [J].Studia Geophysica et Geodaetica,2008,52:1-18.

[19] TURNER J D.A Non-iterative and Non-singular Perturbation Solution for Transforming Cartesian to Geodetic Coordinates[J].Journal of Geodesy,2009,83:139-145.

[20] ZHANG C D,HSU H T,WU X P.An Alternative Algebraic Algorithm to Transform Cartesian to Geodetic Coordinates [J].Journal of Geodesy,2005,79:413-420.

[21] LIGAS M.Cartesian to Geodetic Coordinates Conversion on a Triaxial Ellipsoid[J].Journal of Geodesy,2012,86: 249-256.

[22] LAUREANO G V,IRENE P B.A Symbolic Analysis of Vermeille and Borkowski Polynomials for Transforming 3D Cartesian to Geodetic Coordinates[J].Journal of Geodesy, 2009,83:1071-1081.

[23] FUKUSHIMA T.Transformation from Cartesian to Geodetic Coordinates Accelerated by Halley’s Method[J].Journal of Geodesy,2006,79:689-693.

[24] ABRAMOWITZ M,STEGUN I A.Handbook of Mathematical Functions with Formulas,Graphs,and Mathematical Tables[M].New York:Dover Publications Incorporation, 1972:14-16.

[25] WEISSTEIN E W.“Series Reversion”from Math World: A Wolfram Web Resource.[EB/OL].[2013-07-07].http:∥mathworld.wolfram.com/seriesreversion.html.

[26] SLOANE N J A.The On-line Encyclopedia of Integer Sequences:The OEIS Foundation.[EB/OL].[2013-07-07].http:∥oeis.org/.

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

Transformation from Cartesian to Geodetic Coordinates Using Lagrange Inversion Theorem

GUO Jiachun1,2,ZHAO Xiuxia3,WU Yanlan2
1.School of Science,Anhui Agricultural University,Hefei 230036,China;2.School of Resources and Environmental Engineering,Anhui University,Hefei 230601,China;3.Fishery Institute of Anhui Academy of Agricultural Sciences, Hefei 230051,China

According to Lagrange Inversion Theorem,a Taylor-series expansion method for transforming Cartesian to geodetic coordinates is obtained,which expresses the sine function of reduce latitude sinμas a convergent polynomial of geodetic coordinates(X,Y,Z).Comparative computation of the new method and Bowring’s formula show that the new method is sufficiently precise in the sense that the maximum error of the latitude is less than 9.71×10－6″,9.67×10－8″for the range of－2×106≤H≤＋1010,truncating at b4and b5respectively,while Bowring’s formula only works well for the range of－105≤H≤＋105.Meanwhile,new algorithm is around 50%～90%faster than the iterative algorithm with the approximate accuracy.

Cartesian coordinates;geodetic coordinates;coordinate transformation;Lagrange inversion theorem

GUO Jiachun(1981－),male,PhD candidate,lecturer,majors in geodesy and ecology.

WU Yanlan

P226

A

1001-1595(2014)10-0998-07

國(guó)家自然科學(xué)基金(41271443);江西省數(shù)字國(guó)土重點(diǎn)實(shí)驗(yàn)室開(kāi)放研究基金(DLLJ201211);安徽農(nóng)業(yè)大學(xué)學(xué)科學(xué)位點(diǎn)建設(shè)項(xiàng)目(XKXWD2013022)

2013-07-25

過(guò)家春(1981—),男,博士生,講師,研究方向?yàn)榇蟮販y(cè)量學(xué)和生態(tài)學(xué)。

E-mail:guojiachun＠ahau.edu.cn

吳艷蘭

E-mail:wylmp＠sina.com

GUO Jiachun,ZHAO Xiuxia,WU Yanlan.Transformation from Cartesian to Geodetic Coordinates Using Lagrange Inversion Theorem[J].Acta Geodaetica et Cartographica Sinica,2014,43(10):998-1004.(過(guò)家春,趙秀俠,吳艷蘭.空間直角坐標(biāo)與大地坐標(biāo)轉(zhuǎn)換的拉格朗日反演方法[J].測(cè)繪學(xué)報(bào),2014,43(10):998-1004.)

10.13485/j.cnki.11-2089.2014.0152

修回日期:2014-08-21