GPS经纬度坐标WGS84到东北天坐标系ENU的转换

  • 常用坐标系介绍
    • 地理坐标系 (Geographic Coordinate System, GCS)
    • 地心地固坐标系 (ECEF)
    • 当地东、北、上 (ENU) 坐标
  • 基坐标相互转化
    • 地理坐标系到地心地固坐标系 (GCS to ECEF)
    • 地心地固坐标系到东北天、站心坐标系 ECEF to ENU
    • 地理坐标系到东北天、站心坐标系(GCS to ENU)

常用坐标系介绍

地理坐标系 (Geographic Coordinate System, GCS)

可以说是最为广泛应用的一个地球坐标系,它给出一点的大地纬度、大地经度和大地高程而更加直观地告诉我们该点在地球中的位置,故又被称作纬经高坐标系。WGS-84坐标系的X轴指向BIH(国际时间服务机构)1984.0定义的零子午面(Greenwich)和协议地球极(CTP)赤道的交点。Z轴指向CTP方向。Y轴与X、Z轴构成右手坐标系。

一句话解释就是:把前面提到的ECEF坐标系用在GPS中,就是WGS-84坐标系。

其中:
(1):大地纬度是过用户点P的基准椭球面法线与赤道面的夹角。纬度值在-90°到+90°之间。北半球为正,南半球为负。
(2):大地经度是过用户点P的子午面与本初子午线之间的夹角。经度值在-180°到+180°之间。
(3):大地高度h是过用户点P到基准椭球面的法线距离,基准椭球面以内为负,以外为正。

GPS经纬度坐标WGS84到东北天坐标系ENU的转换-编程知识网
参考资料:地理信息系统中常用坐标系

地心地固坐标系 (ECEF)

在GPS中使用的被称为地球中心,地球固定(ECEF)。ECEF使用三维XYZ坐标(以米为单位)来描述GPS用户或卫星。“地球”这个词来自于坐标轴的原点(0,0,0)位于质心处(通过多年的跟踪卫星确定轨迹)。“地球固定”这个词意味着坐标轴相对于地球是固定的(也就是说,它们随地球旋转)。其 x 轴延伸通过本初子午线(经度 0 度)和赤道(0度纬度)。z 轴延伸穿过真正的北极(即与地球自转轴重合)。y 轴完成右手坐标系,穿过赤道和 90 度经度,如图 1所示

GPS经纬度坐标WGS84到东北天坐标系ENU的转换-编程知识网
图1 ECEF直角坐标系

当地东、北、上 (ENU) 坐标

站心坐标系也叫东北天坐标系 ENU,用于表示以观察者为中心的其他物体
的运动规律。以站心为坐标系原点,z 轴与椭球法线重合,向上为正,y 与椭球
短半轴重合(北向),x 轴与椭球长半轴重合(东向)
GPS经纬度坐标WGS84到东北天坐标系ENU的转换-编程知识网

基坐标相互转化

地理坐标系到地心地固坐标系 (GCS to ECEF)

这里因为地球近似为一个椭球,所以沿着椭球表面的法线方向不在OeO_{e}OeMMM的连线方向,首先我们假设它的椭球表面的法线方向的反向延长线交ZZZ轴与于点M′M{'}M,设∣PQ∣=N,∣PM∣=h|PQ|= N,|PM| = hPQ=N,PM=h,在XECEF−ECEF−ECEF在X_{ECEF}-_{ECEF}-_{ECEF}XECEFECEFECEF基坐标系下,根据几何关系,M的坐标为
((N+h)cos⁡(φ)cos⁡(λ),(N+h)cos⁡(φ)sin⁡(λ),(N+h)sin⁡(φ))−OQ)\left((N+h)\cos \left(\varphi\right) \cos (\lambda), (N+h) \cos \left(\varphi\right) \sin (\lambda), (N+h) \sin \left(\varphi\right)\right)-OQ)((N+h)cos(φ)cos(λ),(N+h)cos(φ)sin(λ),(N+h)sin(φ))OQ)

GPS经纬度坐标WGS84到东北天坐标系ENU的转换-编程知识网
然后我们将MMM点所在的子午圈椭圆投影到一个二维的平面上,如下图所示,(gps的高度h不是与地心在一条直线上)
GPS经纬度坐标WGS84到东北天坐标系ENU的转换-编程知识网

首先设p点所在的子午圈的方程设为
x2Re2+z2Rp2=1①\tag*{①}\frac{x^{2}}{R_{e}^{2}}+\frac{z^{2}}{R_{p}^{2}}=1Re2x2+Rp2z2=1 椭圆的扁率定义为
f=Re−RpRef=\frac{R_{e}-R_{p}}{R_{e}} f=ReReRp 椭圆的离心率定义为
e=Re2−Rp2Re②e=\tag*{②}\frac{\sqrt{R_{e}^{2}-R_{p}^{2}}}{R_{e}}e=ReRe2Rp2由②式可得Rp=Re1−e2③\tag*{③}\ R_{p}=Re \sqrt{1-e^{2}}\, Rp=Re1e2
将椭圆方程①两边同时对xx\,x求导,并结合③式,可得2xRe2+2zdz/dxRe2(1−e2)=0\frac{2 x}{R_{e}^{2}}+\frac{2 z \mathrm{~d} z / \mathrm{d} x}{R_{e}^{2}\left(1-e^{2}\right)}=0Re22x+Re2(1e2)2z dz/dx=0移项整理得dzdx=−(1−e2)xz\frac{\mathrm{d} z}{\mathrm{~d} x}=-\left(1-e^{2}\right) \frac{x}{z} dxdz=(1e2)zx式中dzdx\frac{\mathrm{d} z}{\mathrm{~d} x} dxdz表示椭圆在P点的切线的斜率,显然切线PE⊥\bot \,法线PQ,所以KPEKPQ=−1K_{PE} K_{PQ}=-1\,KPEKPQ=1,因为PQPQPQ的斜率为tan⁡φ\tan \varphi \,tanφ,则
dzdxtan⁡φ=−(1−e2)xztan⁡φ=−1\frac{\mathrm{d} z}{\mathrm{~d} x} \tan \varphi =-\left(1-e^{2}\right) \frac{x}{z} \tan \varphi =-1\, dxdztanφ=(1e2)zxtanφ=1
整理可得z=x(1−e2)tan⁡φ④\tag*{④}z=x\left(1-e^{2}\right) \tan \varphiz=x(1e2)tanφ
结合③④带入椭圆方程①中,可得到以地理纬度φ\varphi \,φ的椭圆的参数方程x=Re1−e2sin⁡2φcos⁡φz=Re(1−e2)1−e2sin⁡2φsin⁡φ}⑤\tag*{⑤}\left.\begin{array}{l} x=\frac{R_{e}}{\sqrt{1-e^{2} \sin ^{2} \varphi}} \cos \varphi \\ z=\frac{R_{e}\left(1-e^{2}\right)}{\sqrt{1-e^{2} \sin ^{2} \varphi}} \sin \varphi \end{array}\right\}x=1e2sin2φRecosφz=1e2sin2φRe(1e2)sinφ
根据子午圈投影的图可得x=Ncosφ⑥\tag*{⑥} x =N cos \varphix=Ncosφ 将⑤中的xx\,x带入⑥式中N=Re1−e2sin⁡2φN =\frac{R_{e}}{\sqrt{1-e^{2} \sin ^{2} \varphi}} N=1e2sin2φRe
因此参数方程可以简写为x=Ncos⁡φz=N(1−e2)sin⁡φ}⑥\tag*{⑥}\left.\begin{array}{l} x=N\cos \varphi \\ z=N\left(1-e^{2}\right) \sin \varphi \end{array}\right\}x=Ncosφz=N(1e2)sinφ}
通过上面的椭圆参数方程⑥,结合投影到2D平面上的Z轴与ECEF的Z轴2D平面上的Z轴与ECEF的Z轴2DZECEFZ是同一个可得MMM点的坐标系zM=N(1−e2)sin⁡φ+hsinφ=(z_{M}=N\left(1-e^{2}\right) \sin \varphi + hsin \varphi =(zM=N(1e2)sinφ+hsinφ=(,结合之前所说的{ECEF坐标系下PPP的坐标系为((N+h)cos⁡(φ)cos⁡(λ),(N+h)cos⁡(φ)sin⁡(λ),(N+h)sin⁡(φ))−OQ)\left((N+h)\cos \left(\varphi\right) \cos (\lambda), (N+h) \cos \left(\varphi\right) \sin (\lambda), (N+h) \sin \left(\varphi\right)\right)-OQ)((N+h)cos(φ)cos(λ),(N+h)cos(φ)sin(λ),(N+h)sin(φ))OQ)
可以得出
xM=(N+h)cos(φ)cos(λ)x_{M} = (N+h)cos(φ)cos(λ)xM=(N+h)cos(φ)cos(λ)
yM=(N+h)cos(φ)sin(λ)y_{M} = (N+h)cos(φ)sin(λ)yM=(N+h)cos(φ)sin(λ)
zM=(b2a2N+h)sin⁡φz_{M} =\left(\frac{b^{2}}{a^{2}} N+h\right) \sin \varphizM=(a2b2N+h)sinφ

附上各参数的含义
φ=latitude\varphi = latitudeφ=latitude
λ=longitude\lambda = longitudeλ=longitude
h=heightaboveellipsoid(meters)h = height above ellipsoid (meters)h=heightaboveellipsoid(meters)
N=RadiusofCurvature(meters),definedas:a1−e2sin⁡2φN = Radius of Curvature (meters), defined as:\frac{a}{\sqrt{1-e^{2} \sin ^{2} \varphi}}N=RadiusofCurvature(meters),definedas:1e2sin2φa,这里的a就是我们上面推导公式的ReR_{e}Re,详细的参考⑥⑥下面的推导.

WGS84 Parameters
a=6378137a = 6378137a=6378137(椭圆的长半轴)这里表示赤道椭圆的长半轴
b=a(1−f)=6356752.31424518b = a (1 -f ) = 6356752.31424518b=a(1f)=6356752.31424518
f=1298.257223563f =\frac{1}{298.257223563}\,f=298.2572235631
e=a2−b2a2e=\frac{\sqrt{a^{2}-b^{2}}}{a^{2}}e=a2a2b2详细含义参照②

我们把MMM点放大,来说说gps的高度
GPS经纬度坐标WGS84到东北天坐标系ENU的转换-编程知识网
高度h,h \,,h,gps ,的的altitude,(从上面图中可知gps的测量高度是GPS天线与参考椭球面的垂直距离,WGS84椭球面在世界范围内近似大地水准面),地球的长半轴aaa与短半轴bbb是已知的,然后地理纬度φ\varphi \,φ也是已知的,所以M点在ECEFECEFECEF坐标系下已知.

地心地固坐标系到东北天、站心坐标系 ECEF to ENU

这里采用几何意义的方式来讲这个转化关系,这里为了直观,我们在p点画ENU坐标系
GPS经纬度坐标WGS84到东北天坐标系ENU的转换-编程知识网因为East(xLocal)East(x Local)East(xLocal)的方向与M点所在的子午线(P点子午线)垂直,所以我们首先把绕着ZECEFZ_{ECEF}ZECEF旋转至XECEFX_{ECEF}XECEFEast(xLocal)East(x Local)East(xLocal)的方向一致①ECEF坐标系ECEF坐标系ECEF绕着ZECEFZ_{ECEF}ZECEF逆时针旋转(π2+λ)(\frac{π}{2} + \lambda)(2π+λ),然后绕着XECEFX_{ECEF}XECEF旋转至ZECEFZ_{ECEF}ZECEFUp(ZLocal)Up(Z Local)Up(ZLocal)方向一致②绕着XECEFX_{ECEF}XECEF逆时针旋转(π2−φ)(\frac{π}{2} – \varphi)(2πφ)③平移的部分就是之前LLA转到ECEF坐标系下的MMM的坐标
Rzecef[(π2+λ)]=(cos⁡[(π2+λ)]sin⁡[(π2+λ)]0−sin⁡[(π2+λ)]cos⁡[(π2+λ)]0001)R_{z_{ecef}}[(\frac{π}{2} + \lambda)]=\left(\begin{array}{ccc}\cos [(\frac{π}{2} + \lambda)] & \sin [(\frac{π}{2} + \lambda)] & 0 \\ -\sin [(\frac{π}{2} + \lambda)] & \cos [(\frac{π}{2} + \lambda)] & 0 \\ 0 & 0 & 1 \end{array}\right)Rzecef[(2π+λ)]=cos[(2π+λ)]sin[(2π+λ)]0sin[(2π+λ)]cos[(2π+λ)]0001

Rxecef[(π2−φ)]=[1000cos⁡[(π2−φ)]sin⁡[(π2−φ)]0−sin⁡[(π2−φ)]cos⁡[(π2−φ)]]R_{x_{ecef}}[(\frac{π}{2} – \varphi)]=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos [(\frac{π}{2} – \varphi)] & \sin [(\frac{π}{2} – \varphi)] \\ 0 & -\sin [(\frac{π}{2} – \varphi)] & \cos [(\frac{π}{2} – \varphi)] \end{array}\right]Rxecef[(2πφ)]=1000cos[(2πφ)]sin[(2πφ)]0sin[(2πφ)]cos[(2πφ)]

Rxecef[(π2−φ)]Rzecef[(π2+λ)]=(−sin⁡λcos⁡λ0−cos⁡λsin⁡φ−sin⁡λsin⁡φcos⁡φcos⁡λcos⁡φsin⁡λcos⁡φsin⁡φ)R_{x_{ecef}}[(\frac{π}{2} – \varphi)] \, R_{z_{ecef}}[(\frac{π}{2} + \lambda)] = \left(\begin{array}{ccc} -\sin \lambda & \cos \lambda & 0 \\ -\cos \lambda \sin \varphi & -\sin \lambda \sin \varphi & \cos \varphi \\ \cos \lambda \cos \varphi & \sin \lambda \cos \varphi & \sin \varphi \end{array}\right)Rxecef[(2πφ)]Rzecef[(2π+λ)]=sinλcosλsinφcosλcosφcosλsinλsinφsinλcosφ0cosφsinφ

在车辆运动中,我们一般车辆启动假设gps天线的初始位置在{Xr,Yr,Zr}\left\{X_{r}, Y_{r}, Z_{r}\right\}{Xr,Yr,Zr},车辆运动的位置为{Xp,Yp,Zp}\left\{X_{p}, Y_{p}, Z_{p}\right\}{Xp,Yp,Zp},则基于启动为原点的ENU的坐标为:
[xyz]=[−sin⁡λcos⁡λ0−cos⁡λsin⁡φ−sin⁡λsin⁡φcos⁡φcos⁡λcos⁡φsin⁡λcos⁡φsin⁡φ][Xp−XrYp−YrZp−Zr]\left[\begin{array}{l} x \\ y \\ z \end{array}\right]=\left[\begin{array}{ccc} -\sin \lambda & \cos \lambda & 0 \\ -\cos \lambda \sin \varphi & -\sin \lambda \sin \varphi & \cos \varphi \\ \cos \lambda \cos \varphi & \sin \lambda \cos \varphi & \sin \varphi \end{array}\right]\left[\begin{array}{c} X_{p}-X_{r} \\ Y_{p}-Y_{r} \\ Z_{p}-Z_{r} \end{array}\right]xyz=sinλcosλsinφcosλcosφcosλsinλsinφsinλcosφ0cosφsinφXpXrYpYrZpZr
这里讲的也很好

地理坐标系到东北天、站心坐标系(GCS to ENU)

GCS和ENU之间的坐标系转换通过先转换为ECEF实现。
代码有时间再放上来吧