高精度地图及定位常用坐标系简介、坐标系转换

1. 常用坐标系

高精度地图及定位常涉及到的坐标系有：地心地固坐标系、大地坐标系、站心坐标系。
1、地心地固坐标系：Earth-Centered, Earth-Fixed，简称$ECEF$。$Z$轴指向协议地球北极方向，$X$轴指向零子午面和赤道的交点，$Y$轴与$Z$轴、$X$轴垂直构成右手坐标系。

2、大地坐标系：地面点的位置用大地经度、大地纬度和大地高度表示。亦称为经纬高坐标系，即$LLA$坐标系。

3、站心坐标系：站心坐标系以用户所在位置P为坐标原点，三个轴分别指向东向，北向和天向，也叫东北天坐标系（$ENU$坐标系）。站心坐标系的天向方向和大地坐标系的高度方向是一致的。站心坐标系用在惯性导航和卫星俯仰角计算中较多。

地心地固坐标系	大地坐标系	站心坐标系
$ECEF$坐标系	$LLA$坐标系	$ENU$坐标系
$XYZ$	经纬高	$xyz$

有一点特别说明：地理坐标系是指用经纬度表示地面点位的球面坐标系。在大地测量学中，对于地理坐标系统中的经纬度有三种描述：即天文经纬度、大地经纬度和地心经纬度。地理研究和小比例尺地图制图对精度要求不高，故常把椭球体当作正球体看待，地理坐标采用地球球面坐标，经纬度均用地心经纬度。地图学中常采用大地经纬度。

2. WGS84 & CGCS2000

$CGCS2000$与$WGS84$坐标系都属于地心地固坐标系。这点可能会有人疑惑，可以理解为$ECEF$坐标系用于$GPS$中，用在美国$GPS$是$WGS84$坐标系，用在中国北斗是$CGCS2000$。

在定义上，$CGCS2000$与$WGS84$是一致的，即关于坐标系原点、尺度、定向及定向演变的定义都是相同的。两个坐标系使用的参考椭球也非常相近，具体地说，在4个椭球常数a、f、GM、ω中，唯有扁率f有微小差异。

$WGS84$和$CGCS2000$两者参数对比：

参数	$WGS84$	$CGCS2000$
长半轴a	6378137.0m	6378137.0m
扁率$f$	1/298.257222101	1/298.257223563
地球自转角速度$\omega$	7.292115*1e-5	7.292115*1e-5
地心引力常数GM	3.986004418*1e14	3.986004418*1e14

椭球参数的区别而引起的同一点经纬度的差异有多大呢？程鹏飞院长的论文给出的研究数据是：“给定点位在某一框架和某一历元下的空间直角坐标，投影到$CGCS2000$椭球和$WGS84$椭球上所得的纬度的最大差异相当于0.11mm。”

关于为什么已经有了$WGS84$我国还要自己搞一套$CGCS2000$，简单地说，$WGS84$是全球的，准确说是为美国服务，顺便服务全球。相对应的$CGCS2000$就是要专心为我国服务。中国$CGCS2000$战略意义是深远的，不具体展开了，有一篇文章写的很好，感兴趣的可以看下：中国使用2000坐标系而不使用$WGS84$坐标系的原因和意义是什么？

3. 坐标系转换

3.1 LLA坐标系转ECEF坐标系

$LLA$坐标系下的$(lon,lat,alt)$转换为$ECEF$坐标系下点$(X,Y,Z)$
$$\{\begin{array}{l}X=(N+alt)cos(lat)cos(lon)\\ Y=(N+alt)cos(lat)sin(lon)\\ Z=(N(1-e^2)+alt)sin(lat)\end{array}$$

其中 $e$为椭球偏心率，$N$ 为基准椭体的曲率半径
$$\{\begin{array}{l}{e}^2=\frac{a^2-b^2}{a^2}\\ N=\frac{a}{\sqrt{1-e^{2}si{n}^{2}lat}}\end{array}$$
由于$WGS84$坐标系扁率$f=(a-b)/a$，偏心率 $e$ 和扁率 $f$ 之间的关系:
$$e^2=f(2-f)$$
坐标转换公式也可以为
$$\{\begin{array}{l}X=(N+alt)cos(lat)cos(lon)\\ Y=(N+alt)cos(lat)sin(lon)\\ Z=(N(1-e^2)+alt)sin(lat)\\ N=\frac{a}{\sqrt{1-f(2-f)si{n}^{2}lat}}\end{array}$$

3.2 ECEF坐标系转LLA坐标系

$ECEF$坐标系下点$(X,Y,Z)$转换为$LLA$坐标系下的$(lon,lat,alt)$
$$\begin{array}{rl}lon& =arctan\frac{y}{x}\\ alt& =\frac{p}{cos(lat)-N}\\ lat& =arctan[\frac{z}{p}(1-e^2\frac{N}{N+alt}{)}^{-1}]\\ p& =\sqrt{x^2+y^2}\end{array}$$
公式3在计算纬度时还会用到纬度值，所以不能直接用上述公式计算纬度角，可以先将大地纬度近似为地心纬度处理，经过几次迭代计算就能收敛。具体计算逻辑详见参考[5]。

3.3 ECEF坐标系转ENU坐标系

用户所在坐标点$P0(x0,y0,z0)$，计算点$P=(x,y,z)$在以$P_0$为坐标系原点的$ENU$坐标$(e,n,u)$，这里需要用到$LLA$坐标系的数据，$P_0$的$LLA$坐标为$LL{A}_0=(lo{n}_0,la{t}_0,al{t}_0)$
$$\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]-\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]$$

$$\left[\begin{array}{c}e\\ n\\ u\end{array}\right]=S\cdot \left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]=\left[\begin{array}{ccc}-sin(lo{n}_0)& cos(lo{n}_0)& 0\\ -sin(la{t}_0)cos(lo{n}_0)& -sin(la{t}_0)sin(lo{n}_0)& cos(la{t}_0)\\ cos(la{t}_0)cos(lo{n}_0)& cos(la{t}_0)sin(lo{n}_0)& sin(la{t}_0)\end{array}\right]\cdot \left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]$$

即坐标变换矩阵：
$$S=\left[\begin{array}{ccc}-sin(lo{n}_0)& cos(lo{n}_0)& 0\\ -sin(la{t}_0)cos(lo{n}_0)& -sin(la{t}_0)sin(lo{n}_0)& cos(la{t}_0)\\ cos(la{t}_0)cos(lo{n}_0)& cos(la{t}_0)sin(lo{n}_0)& sin(la{t}_0)\end{array}\right]$$
上述：$enu$是$P$在$P_0$(东北天)坐标系下坐标

3.4 ENU坐标系转ECEF坐标系

S为单位正交矩阵
$$\begin{array}{rl}& S^{-1}=S^T\\ & \left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]=S^{-1}\cdot \left[\begin{array}{c}e\\ n\\ u\end{array}\right]=S^T\cdot \left[\begin{array}{c}e\\ n\\ u\end{array}\right]\end{array}$$

3.5 LLA坐标系转ENU坐标系

上述可以看到，从$LLA$坐标系转到$ENU$坐标系有较多计算量，在考虑地球偏心率$e$很小的前提下，可以做一定的近似公式计算：
$$\left[\begin{array}{c}\Delta e\\ \Delta n\\ \Delta u\end{array}\right]=\left[\begin{array}{ccc}a\cdot cos(lat)\cdot \Delta lon& 0& 0\\ 0& a\cdot \Delta lat& 0\\ 0& 0& \Delta alt\end{array}\right]$$

参考文献及博客链接：
[1]魏子卿.2000中国大地坐标系及其与WGS84的比较[J].大地测量与地球动力学,2008(05):1-5.
[2]程鹏飞,文汉江,成英燕,王华.2000国家大地坐标系椭球参数与GRS 80和WGS 84的比较[J].测绘学报,2009,38(03):189-194.
[3]https://www.zhihu.com/question/35775670/answer/362106176
[4]https://www.cnblogs.com/langzou/p/11388520.html
[5]https://www.cnblogs.com/taqikema/p/8678596.html