空间直角坐标与大地坐标转换的拉格朗日反演方法

过家春,赵秀侠,吴艳兰

1.安徽农业大学理学院,安徽合肥 230036;2.安徽大学资源与环境工程学院,安徽合肥 230601; 3.安徽省农业科学院水产研究所,安徽合肥 230051

空间直角坐标与大地坐标转换的拉格朗日反演方法

过家春1,2,赵秀侠3,吴艳兰2

1.安徽农业大学理学院,安徽合肥 230036;2.安徽大学资源与环境工程学院,安徽合肥 230601; 3.安徽省农业科学院水产研究所,安徽合肥 230051

以拉格朗日反演理论为基础,导出空间直角坐标向大地坐标转换的一种新的直接解法。该方法将归化纬度的正弦函数sinμ表达为以空间直角坐标(X,Y,Z)为基础的相关变量的多项式。为验证公式的精度水平和实用性,将WGS-84椭球参数代入验算,结果表明,新解法在－2×106m≤H≤＋1010m范围内,展开至b4项纬度反解误差不超过9.71×10－6″,展开至b5项的误差不超过9.67× 10－8″,有效范围比Bowring公式广,在H≥＋105m的区域精度明显优于Bowring公式。与迭代算法相比,新方法在H≥－2×106m的区域精度与迭代算法精度相当,但计算效率明显优于迭代算法,展开至b4项的CPU执行时间分别约为迭代4次的1/2、展开至b5项的CPU执行时间约为迭代5次的1/10。兼顾精度和效率,本文算法优于Bowring公式和迭代算法。

空间直角坐标;大地坐标;坐标转换;拉格朗日反演定理

1 引 言

空间直角坐标系与大地坐标系之间的转换是空间大地测量技术中的基本问题,亦是近几十年来国内外大地测量学者讨论的热点问题之一[1-7]。众所周知,由大地坐标解算空间直角坐标可直接得到,而由空间直角坐标反解大地坐标则比较复杂。国内外学者对此提出许多不同的解法,总体上可归纳为间接解法和直接解法两大类。间接解法以各种迭代算法为主[8-13],直接解法则较为多样,包括构造并求解不同形式的一元四次方程[14-17]、各类近似公式、闭合解析解[18-23]等。各类直接解法中,Bowring公式较为著名,但精度不够理想,若要提高精度,则需进一步迭代;而其他各类直接解法,或因公式过于冗长,或因需要进行复杂的求根讨论,以致求解过程复杂、推广困难。因此,对于这一问题,国内外大都仍采用经典的迭代算法求解。

本文以拉格朗日反演理论为基础,导出空间直角坐标向大地坐标转换的泰勒级数展开新方法。

2 理论基础及基本思路

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

theorem)[24-26]

若函数y＝f(x)在点x0的某一邻域U(x0)内可展开为泰勒级数

且f′(x0)≠0,则其反函数x＝f－1(y)在相应的y0点的某一邻域V(y0)内也可展开为泰勒级数

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

拉格朗日反演定理表明,在一阶导数不等于零的情况下,解析函数的反函数也可展开为泰勒级数。

2.2 基本思路

旋转椭球下,大地坐标(L,B,H)与空间直角坐标(X,Y,Z)之间的关系为

图1 空间直接坐标与大地坐标之关系Fig.1 Cartesian coordinates and geodetic coordinates

图1中,OI＝Ne2sin B＝be′2sinμ,OQ＝Ne2cos B＝ae2cosμ,C为P0点的曲率中心,坐标为(ae2cos3μ,－be′2sin3μ),对应的JK弧为椭圆的渐屈线。以P、C两点坐标求解tan B即为Bowring公式的思路。

由式(4)及图1关系,可得

式(8)、式(9)两式中an、bn的关系为[25-26]

由归化纬度与大地纬度的关系式(5)可得

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

上式椭球面上及椭球外取“＋”,椭球内取“－”。

3 系数求解及简化

由上述讨论可首先求得式(8)的系数

再由关系式(10)可得

该系数较为复杂,为简化其结构,令

综合以上各式可得

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

以上为使推导过程方便,引入参数w、z,实际计算时,r、t1—t4可由空间直角坐标(X,Y,Z)直接得出

4 验算分析

尽管式(17)本质上是以拉格朗日反演定理为基础的泰勒级数展开,但由于该问题的复杂性,按泰勒级数的余项理论给出其误差估计仍十分困难。为验证公式的正确性及其精度水平,下面以WGS-84椭球参数代入验算,并与经典的直接解法Bowring公式及迭代算法进行对比。

4.1 盲区对比

式(17)与许多已有的其他解法一样,在接近椭球中心的区域存在盲区。

首先,由以上讨论可知,该解算方法要求(w2＋z2)3/2－e2w2≠0。同时,经验算分析,当(w2＋z2)3/2－e2w2＜0、z≠0时,式(17)右端的值域超出sinμ的值域范围。由此可见,(w2＋z2)3/2－e2w2≤0、z≠0所含区域为该方法的盲区。

已有算法的盲区(或多值区)多为椭球的渐屈线附近及其以内区域,如Bowring公式、经典迭代算法以及文献[20]所给的算法,与本文算法的盲区(w2＋z2)3/2－e2w2≤0、z≠0,即

的区域关系如图2所示。

图2 盲区对比Fig.2 Comparison of the“dead zone”between Bowring’s formula and new method

4.2 精度对比

以上盲区范围距地心距离不超过42.84 km,处于地核内核内部,卫星大地测量一般不涉及这一区域。因此,从实用的角度看,各国学者一般最大范围取大地高－6×106～＋1010m的范围进行分析,本文亦取此范围进行分析。将本文算法与Bowring公式、经典迭代算法对比,其纬度解算误差对比如表1和图3所示。因高程解算精度主要取决于纬度解算,以下仅对纬度解算误差进行分析。

4.3 效率对比

为对比3种算法的计算效率,笔者在同一台计算机上,用Mathematica软件以纬度1°、大地高108m为步长,分别利用3种方法计算B∈[0°, 90°],H∈[0,1010]范围内91×101个点的总CPU执行时间,进而得到单个离散点的平均CPU执行时间,计算结果对比如表2所示。

表1 本文算法与Bowring公式、迭代算法的误差对比Tab.1 Errors of Cartesian to geodetic coordinates between Bowring’s formula,iterative algorithm and new method

图3 本文算法与Bowring公式、迭代算法的误差对比Fig.3 Comparison of the errors of geodetic latitude between Bowring’s formula,iterative algorithm and our algorithm

表2 3种算法的CPU执行时间对比Tab.2 Comparison of the CPU times between Bowring’s formula,iterative algorithm and new method

4.4 综合分析

通过以上对比分析可知:

(1)新算法与两种经典算法在接近椭球中心区域均存在盲区,需进一步研究,但实际应用一般不涉及这一区域。

(2)精度上,若按照对应的点位误差不超过0.1 mm的原则,Bowring公式的有效范围约为－105m≤H≤＋105m,在H≥＋105m以外区域解算精度迅速下降,最弱精度位于H＝2a、B＝π/4处,约＋1.84×10－3″,相当于经线方向的0.11 m的点位误差。

新解法展开至b4项、b5项,以及迭代法迭代4次、5次在H≥－2×10－6m的区域范围能够满足要求,有效范围比Bowring公式大。

(3)计算效率上,Bowring公式最高,迭代法计算效率明显低于新解法计算效率。

综上可知,在－105m≤H≤＋105m的范围内,Bowring公式精度最优,计算效率上占有显著优势,对于地面上的定位与导航应用,Bowring公式占据优势。

但实际上导航卫星大都是中地球轨道卫星,卫星高度在107m级,如GPS、GLONASS、北斗等卫星导航系统的卫星高度分别约为20 200 km、19 100 km、21 500 km,此时Bowring公式精度明显较弱,如表1所示。兼顾精度和效率,本文算法较之于Bowring公式和迭代算法有一定优势。

此外,与其他直接解法相比,本文算法不必进行复杂的讨论,公式相对简洁,可根据精度需要决定末项取舍。

5 结 论

本文以拉格朗日反演理论为基础,得到空间直角坐标向大地坐标转换的泰勒级数展开新方法,是一种新的直接解法,且经验算分析表明其精度可靠,能够满足实际应用精度要求。新解法由式(17)、式(18)两式构成,与已有直接解法相比,公式相对简洁,无须复杂讨论,可根据精度需要决定末项取舍,可扩展性好。同时,该解法随大地高的增加精度提高显著,这对现代空间技术是十分有利的,如航天工程、探月工程等。该方法的缺点是接近椭球中心区域存在盲区,这亦是Bowring公式、经典迭代算法等解法的共同问题,需要进一步研究。

[1] CHEN Junyong.Direct Solution of Transforming Cartesian to Geodetic Coordinates[C]∥Special Issue on Geodesy.Beijing:Surveying and Mapping Press,1979:67-72.(陈俊勇.空间直角坐标与大地坐标换算的直接解法[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.(曾启雄.空间直角坐标直接解算大地坐标的闭合公式[J].测绘学报,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.(孔祥元,郭际明,刘宗泉.大地测量学基础[M].武汉:武汉大学出版社,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.(曾启雄.空间直角坐标直接计算高斯平面直角坐标公式[J].测绘学报,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].测绘学报,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网三维平差的数学模型[J].测绘学报,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.(束蝉方,李斐,沈飞.空间直角坐标向大地坐标转换的新算法[J].武汉大学学报:信息科学版,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/.

(责任编辑:陈品馨)

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

国家自然科学基金(41271443);江西省数字国土重点实验室开放研究基金(DLLJ201211);安徽农业大学学科学位点建设项目(XKXWD2013022)

2013-07-25

过家春(1981—),男,博士生,讲师,研究方向为大地测量学和生态学。

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.(过家春,赵秀侠,吴艳兰.空间直角坐标与大地坐标转换的拉格朗日反演方法[J].测绘学报,2014,43(10):998-1004.)

10.13485/j.cnki.11-2089.2014.0152

修回日期:2014-08-21