空间直角坐标系与大地坐标系互转(C++实现)

在空间大地测量数据过程中,经常要进行坐标系系统转换以及基准转换。其中空间直角坐标转换(x,y,z)与大地坐标(B,L,H)之间的转换就是测量数据处理过程中一个最基本和最重要的问题。由大地坐标转空间直角坐标比较简单,而空间直角坐标转大地坐标则相对复杂(本质上就是数学问题)。

目前由空间直角坐标转大地坐标的方法主要三种1:

  • Bowring 直接解法
  • 迭代解法
  • 数值导数法

其中Bowring直接解法在椭球表面附近(1000km内)精度较高,基本适合大多数情况,而数值导数法精度虽然,计算复杂度也较高。而迭代法计算复杂度较低,是常用的转换计算方法。不过迭代次数影响计算速度,所以最好采用直接解法的值作为初始值,然后利用迭代法得到更高精度的转换结果。

一、计算必要条件

实现转换过程只需知道参考椭球的长短半轴的长度(a,b)即可。目前最为常用的是美国国防部制图局在1984年构建的WGS94。下面是最为常见的几种参考椭球2:

椭球名称 长半轴(m) 短半轴(m) 使用国家和地区
克拉克(Clarke) 1866 6 378 206.4 6 356 583.8 北美
白塞尔(Bessel)1841 6 377 397.155 6 356 078.965 日本及台湾
International 1924 6 378 388 6 356 911.9 欧洲、北美及中东
克拉索夫斯基1940 6 378 245 6 356 863 俄罗斯、中国
1975年国际会议推荐 6 378 140 6 356 755 中国
GRS 1980 6 378 137 6 356 752.3141
WGS 1984 6 378 137 6 356 752.3142 全球
球形(6371 km) 6 371 000 6 371 000

本次计算采用了WGS 1984参考椭球。

二、BLH转xyz

由大地测量学基础知识可知,xyz转BLH的基本公式如下:
x =(N + H) cosB cosL \\ y =(N + H) cosB sinL \\ z =[N(1-e^2) + H] sinB
其中:
N = \frac{a}{\sqrt{1-e^2sin^2B}}, e^2 = \frac{ a^2 - b^2 }{a^2}

三、xyz转BLH

直角坐标转大地坐标系比较麻烦的就是计算纬度,下面是两种计算方式:

3.1 直接算法计算纬度

B = arctan(\frac{z+e'^2 b sin^3\theta}{\sqrt{x^2+y^2} - e^2 a cos^3 \theta})
其中:
\theta = arctan(\frac{az}{b\sqrt{x^2 + y^2}}), e'^2 = \frac{a^2 - b^2}{b^2}

3.2 使用迭代计算纬度

B_{n+1} = arctan\left [ \frac{z}{\sqrt{x^2 + y^2}} \left ( 1+ \frac{ae^2sinB_n}{z \sqrt{1-e^2sin^2B_n}} \right  )\right ]
为了减少迭代次数,可以使用直接算法的计算结果作为初始值来迭代计算纬度。

3.3 经度L与大地高H的计算

由于经度分东经和西经以及会出现tan 90度极值情况,所以计算经度的时候分类讨论:
L = arctan \frac{y}{x} \qquad (x>0) \\ L =  \pi /2 \quad (x = 0, y > 0)\\ L =  3\pi/2 \quad (x = 0, y > 0) \\ L =  arctan \frac{y}{x} + \pi \qquad (x>0)
而大地高则可以通过下式计算:

H = \sqrt{x^2 + y^2} cosB + zsinB - N ( 1 - e^2 sin^2B )

四、精度验证

为了保证计算过程的正确性和确认方法的精度,需要对上述转换过程进行测试。测试方法为:
将空间直角坐标转换成大地坐标,然后再转换为空间直角坐标,然后将计算结果进行比较
经过测试表明,基本精度是能够达到10^{-6} \sim 10^{-8} 米的水准的。

实现c++源码地址已托管到github
https://github.com/yinflying/BlogSource/blob/master/xyz2blh.cpp

Reference


  1. 史海锋.空间直角坐标与大地坐标转换算法研究.大地测量学与地球动力学.2012/10 
  2. 参考椭球体.维基百科.2017/09 
此条目发表在CODING, GNSS分类目录,贴了, 标签。将固定链接加入收藏夹。

发表评论

电子邮件地址不会被公开。