地理测绘基础知识(1) 坐标系经纬度与ECEF直角坐标的基本换算

经纬度与ECEF直角坐标的基本换算

我们目前最常用的全球坐标系是WGS-84坐标系，各种手机、地图基本用经纬度来标记位置。然而，经纬度对于空间的计算是很复杂的，需要很多三角函数操作。平面直角坐标系利用向量的运算，可以非常方便的计算角度、距离等参数，在实际应用中往往作为中间计算的工具。

目前用到的很多GIS、遥感与测绘工具里都有这种功能，比如利用 libproj、RTK等工具，直接进行转换。

1. WGS-84 坐标系

WGS-84 坐标系是目前使用最为广泛的全球坐标系统。其示意图如下：

WGS-84坐标系里，地球等效为一个椭圆，直角坐标系原点O位于地球的中心。

• 坐标轴 $\vec{Z}=\overrightarrow{ON}$ 指向北极点。
• 坐标轴 $\vec{X}=\overrightarrow{OP}$ 指向格林威治0度经线与赤道面的交点P。
• 坐标轴 $\vec{Y}=\vec{X}\otimes \vec{Z}$ 呈右手螺旋关系。

当前位置 T 的经纬度定义：

从高度 H 开始定义：

• T是一个距离参考海平面（椭球面）高度为 $\vec{H}=\overrightarrow{tT}$ 的点。
• t是点T在椭球面上的正射投影，矢量 $\overrightarrow{tT}$ 是椭球面上t点切面的法向量。
• 需要注意的是，和正球面不同，法向量$\overrightarrow{tT}$ 的反向延长线与Z轴的交点是 O’，和原点O并不重合。

对于经度 $\theta$ 的定义，是这样的：

• 0度经度大圆$\stackrel{⌢}{\mathrm{NOP}}$是由北极点 N, 地心 O, 格林威治经线与赤道面的交点 P 组成的平面。
• 当前经度大圆 $\stackrel{⌢}{N{O}^{\prime }T}$是由北极点 N, 交点O’, 当前位置T组成的平面。
• O、 O’ 、t、 T 都在当前经度面上，直线 O’T 与赤道面的交点为u，矢量 $\overrightarrow{Ou}$ 与赤道的交点是Q。
• O’、u、t、T四点共线。
• 经度 $\theta =\measuredangle POQ$

对于纬度 $\phi$ 的定义，是这样的：

• 椭球的切平面的法线与Z的交点 O’ 和原点并不重合。
• 纬度定义为法向量 $\overrightarrow{uT}$ 与赤道面的夹角。在北半球为正，南半球为负
• 纬度 $\phi =\measuredangle TuQ=\measuredangle Ou{O}^{\prime }$

2 从经纬度转换到ECEF

• ECEF 坐标系有利于利用矢量计算简化角度、距离的求取。

已知:

• 赤道半径 a = 6378137 米，对应上图$\overrightarrow{OP}$的长度;
• 椭球扁率 $f=(a-b)/a$= 1.0/298.257223563。扁率就是衡量椭球是不是很“扁”，正球的扁率为0.
• 椭球偏心率 $e=\sqrt{a^2-b^2}/a$= 0.0818191908426
• e和f的关系为 $e^2=f(2-f)$
• 极地半径 b = 6356752.3142 米，对应上图$\overrightarrow{ON}$的长度;

设：经度 $\theta$， 纬度 $\phi$， 海拔 H，则有：

$$r=\left|\vec{R}\right|=\left|\overrightarrow{O^{\prime }t}\right|=\frac{a}{\sqrt{1-e^2{\sin }^2\phi }}=\frac{a}{\sqrt{1-f(2-f){\sin }^2\phi }}$$

经纬度极坐标可以利用矢量 $\vec{R}$ 构造一个Z轴平移了的正球，进行极坐标运算。这个正球的球心是 O’。由于x,y没有平移，故而两个坐标系的x、y是重合的。Z坐标z’与z的关系是平移，平移量为OO’的长度：

$$\left|\overrightarrow{O{O}^{\prime }}\right|=r{e}^2\sin \phi$$

下图是把经度大圆$\stackrel{⌢}{\mathrm{NOQ}}$切割出来，观察：

在正球模型下，根据极坐标的基本定义，可以直接计算x,y:

$$x=\left(r+H\right)\cdot \cos \phi \cdot cos\theta$$

$$y=\left(r+H\right)\cdot \cos \phi \cdot sin\theta$$

Z的坐标，要先计算后，再平移即可得到:

$$z^{\prime }=\left(r+H\right)\cdot sin\phi$$
$$z=z^{\prime }-r{e}^2\sin \phi$$

带入后，可写成如下等效形式：

$$z=(r(1-e^2)+H)\sin \phi$$
$$z=(r(1-f(2-f))+H)\sin \phi =(r(1-f{)}^2+H)\sin \phi$$

含有扁率的化简，是基于扁率 $f=(a-b)/a=1-b/a$ 计算得到的。

通过上面的处理，极坐标的 $T(\theta ,\phi ,H)$ 便转换为平面直角坐标 $T(x,y,z)$ 。

3 从 ECEF 到经纬度

对上述计算过程而言，逆向运算可以立刻得到经度，却不能解析得到纬度。纬度需要进行迭代。

主要原因是：

$$\left|\overrightarrow{O{O}^{\prime }}\right|=r{e}^2\sin \phi$$

的具体取值不知道导致的。

(1) 求取经度

由于

$$x/y=\tan \theta$$

故而

$$\theta ={\tan }^{-1}(x/y)$$

注意处理好象限问题，即可得到经度。

(2) 纬度、海拔计算

对于纬度，则比较复杂。求解三角方程

$$(r+H)\sin \phi =z+\left|\overrightarrow{O{O}^{\prime }}\right|$$

由于 H 不知道，r 也不知道，加入第二条件（相切，H向量与r共线）后，得到的是二元四次三角方程，很难解析求取。

可以首先假设$\left|\overrightarrow{O{O}^{\prime }}\right|=0$, 求取一个粗略的纬度，并得到 d，利用d反过来更新其他数值，进行迭代。这种迭代收敛的前提，是因为每次计算出的纬度一定小于真实的纬度，这是由三角关系决定的。

迭代的效果是不断地把原点向真实的原点移动，直到z达到一个很小的误差。迭代方法：

设起始偏移 $d=\left|\overrightarrow{O{O}^{\prime }}\right|=0$，迭代开始：

$$R^{\prime }=r+H=\sqrt{x^2+y^2+(z+d{)}^2}$$

$${\phi }^{\prime }={\sin }^{-1}\frac{z+d}{R^{\prime }}$$

$$r^{\prime }=\frac{a}{\sqrt{1-e^2{\sin }^2{\phi }^{\prime }}}$$

$$d^{\prime }=r^{\prime }{e}^2\sin {\phi }^{\prime }$$

$$z^{\prime }=R^{\prime }sin{\phi }^{\prime }-d^{\prime }$$

$$d=d^{\prime }$$

当
$$err=\left|z-z^{\prime }\right|<\epsilon$$
时停止迭代。此时

$$\phi ={\sin }^{-1}\frac{z+d}{\sqrt{x^2+y^2+(z+d{)}^2}}$$

$$H=R^{\prime }-r^{\prime }$$

4 代码与工程

代码与工程参考

https://gitcode.net/coloreaglestdio/geocalc/-/blob/master/geocalc.h

核心转换逻辑如下：

/*!* GEOCalc 实现了基本的WGS-84坐标系的计算。* by CEStdio* 1997-2023*
*/
#include <cassert>
#include <cmath>
namespace CES_GEOCALC {inline const double a = 6378137;
inline const double pi = 3.14159265358979323846;
inline const double f = 1.0 / 298.257223563;
inline const double e = sqrt(f * (2 - f));
inline const double es = f * (2 - f);
inline const double deg2rad = pi / 180.0;
inline const double rad2deg = 180.0 / pi;/*!* \brief lla2ecef 经纬度坐标到ECEF，* \param lla  纬经高(默认)/经纬高, 量纲是度(默认)/弧度、米* \param ecef xyz，量纲是米* \param rad  角度量纲开关，false 是度，true 是弧度* \param latfirst 经纬度顺序，false 是经度\纬度\高度，true 是纬度\经度\高度
*/
inline void lla2ecef(const double lla[/*3*/],double ecef[/*3*/],const bool rad = false,const bool latfirst = true)
{const double lat = latfirst ? (rad ? (lla[0]) : (lla[0] * deg2rad)): (rad ? (lla[1]) : (lla[1] * deg2rad));const double lon = latfirst ? (rad ? (lla[1]) : (lla[1] * deg2rad)): (rad ? (lla[0]) : (lla[0] * deg2rad));const double H = lla[2];const double r = a / sqrt(1 - es * sin(lat) * sin(lat));const double d = r * es * sin(lat);ecef[0] = (r + H) * cos(lat) * cos(lon);ecef[1] = (r + H) * cos(lat) * sin(lon);ecef[2] = (r + H) * sin(lat) - d;
}/*!* \brief ecef2lla ECEF到经纬度坐标* \param ecef xyz，量纲是米* \param lla  纬经高(默认)/经纬高, 量纲是度(默认)/弧度、米* \param maxiter  最大迭代次数* \param piter  迭代次数输出，可以为null* \param eps  Z误差门限* \param rad  角度量纲开关，false 是度，true 是弧度* \param latfirst 经纬度顺序，false 是经度\纬度\高度，true 是纬度\经度\高度* \return 迭代收敛标志
*/
inline bool ecef2lla(const double ecef[/*3*/],double lla[/*3*/],const int maxiter = 32,int * piter = nullptr,const double eps = 1e-10,const bool rad = false,const bool latfirst = true)
{const double x = ecef[0], y = ecef[1], z = ecef[2] < 0 ? -ecef[2] : ecef[2];const int south = ecef[2] < 0 ? -1 : 1;double d = 0;double err = 0;int iter = 0;double Rp = 0, sinphi = 0, rp = 0, dp = 0, zp = 0;do {Rp = sqrt(x * x + y * y + (z + d) * (z + d));sinphi = (z + d) / Rp;rp = a / sqrt(1 - es * sinphi * sinphi);dp = rp * es * sinphi;zp = Rp * sinphi - dp;d = dp;assert(z >= zp);err = z - zp; //z always > zp++iter;} while (err > eps && iter < maxiter);const double lat = asin(sinphi) * south;const double lon = atan2(y, x);const double H = Rp - rp;lla[latfirst?0:1] = lat * (rad?1:rad2deg);lla[latfirst?1:0] = lon * (rad?1:rad2deg);lla[2] = H;if (piter)*piter = iter;return err <= eps;
}}; // namespace CES_GEOCALC

对各种边界和迭代次数进行测试, 以确定迭代的收敛性:

using namespace CES_GEOCALC;for (double lon = -180; lon <=180; lon+=60){for (double lat = -90; lat <=90; lat += 15){double LLA[] = {lat,lon,10000};double ECEF[] = {0, 0, 0};double LLA2[] = {0, 0, 0};lla2ecef(LLA, ECEF);int iter = 0;ecef2lla(ECEF,LLA2,32,&iter);printf("LAT=%12.7lf,LON=%12.7lf,dLAT=%10.7lf,dLON=%10.7lf,dALT=%5.3lf, iter = %d\n",lat,lon,LLA[0]-LLA2[0], LLA[1]-LLA2[1], LLA[2]-LLA2[2],iter);}}

输出：

LAT= -90.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -75.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -60.0000000,LON=-180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -45.0000000,LON=-180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT= -30.0000000,LON=-180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT= -15.0000000,LON=-180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=   0.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 1
LAT=  15.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  30.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  45.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT=  60.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  75.0000000,LON=-180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  90.0000000,LON=-180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -90.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 2
LAT= -75.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT= -60.0000000,LON=-120.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT= -45.0000000,LON=-120.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 7
LAT= -30.0000000,LON=-120.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT= -15.0000000,LON=-120.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT=   0.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 1
LAT=  15.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT=  30.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT=  45.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 7
LAT=  60.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT=  75.0000000,LON=-120.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT=  90.0000000,LON=-120.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 2
LAT= -90.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 2
LAT= -75.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT= -60.0000000,LON= -60.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT= -45.0000000,LON= -60.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 7
LAT= -30.0000000,LON= -60.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT= -15.0000000,LON= -60.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT=   0.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 1
LAT=  15.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT=  30.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 8
LAT=  45.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 7
LAT=  60.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT=  75.0000000,LON= -60.0000000,dLAT=-0.0000000,dLON=-0.0000000,dALT=0.000, iter = 6
LAT=  90.0000000,LON= -60.0000000,dLAT= 0.0000000,dLON=-0.0000000,dALT=0.000, iter = 2
LAT= -90.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -75.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -60.0000000,LON=   0.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -45.0000000,LON=   0.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT= -30.0000000,LON=   0.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT= -15.0000000,LON=   0.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=   0.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 1
LAT=  15.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  30.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  45.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT=  60.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  75.0000000,LON=   0.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  90.0000000,LON=   0.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -90.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -75.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -60.0000000,LON=  60.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -45.0000000,LON=  60.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT= -30.0000000,LON=  60.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT= -15.0000000,LON=  60.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=   0.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 1
LAT=  15.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  30.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  45.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT=  60.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  75.0000000,LON=  60.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  90.0000000,LON=  60.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -90.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -75.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -60.0000000,LON= 120.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -45.0000000,LON= 120.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT= -30.0000000,LON= 120.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT= -15.0000000,LON= 120.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=   0.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 1
LAT=  15.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  30.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  45.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT=  60.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  75.0000000,LON= 120.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  90.0000000,LON= 120.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -90.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2
LAT= -75.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -60.0000000,LON= 180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT= -45.0000000,LON= 180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT= -30.0000000,LON= 180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT= -15.0000000,LON= 180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=   0.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 1
LAT=  15.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  30.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 8
LAT=  45.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 7
LAT=  60.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  75.0000000,LON= 180.0000000,dLAT=-0.0000000,dLON= 0.0000000,dALT=0.000, iter = 6
LAT=  90.0000000,LON= 180.0000000,dLAT= 0.0000000,dLON= 0.0000000,dALT=0.000, iter = 2