地理坐标转换
发布时间：2024-06-29
归档分类：
标签：

近些日子的工作频繁接触到了大地测量中的各种知识，在编写程序的时候更是涉及各种地理坐标的转换，虽然没有深入去学习其核心的变化原理，但是简单地使用现存的公式配合程序进行计算还是可以做到的，以下是有关知识的记录。

1.名词解释

1.1 LLA 坐标系

LLA 坐标是地理坐标系中的一种表示方法，LLA 代表经纬度和高度（Latitude, Longitude, and Altitude），它通常用于描述地球表面上的位置。

  • Latitude (纬度)：表示从地球赤道向北或向南的角度，范围是 -90° 到 +90°。正值表示北纬，负值表示南纬。
  • Longitude (经度)：表示从本初子午线（通常是通过格林尼治的子午线）向东或向西的角度，范围是 -180° 到 +180°。正值表示东经，负值表示西经。
  • Altitude (高度)：表示相对于平均海平面的高度，可以是正值（高于海平面）或负值（低于海平面）。

LLA 坐标系也叫全球地理坐标系、大地坐标系、WGS-84坐标系。纬度和经度的数值可以以多种不同的单位或格式出现：

  • 六十进制度：度、分、秒：40° 26′ 46“ N 79° 58′ 56” W
  • 度和十进制分：40° 26.767′ N 79° 58.933′ W
  • 十进制度： +40.446 -79.982

横纬竖经，在计算的过程中，主要也是将之转换为十进制来计算的，以下是计算公式： $$ decimal = degress + \frac{minutes}{60} + \frac{seconds}{3600} $$

关于 WGS-84 （ World Geodetic System 1984，1984年世界大地测量系统），它是目前全球范围内使用最广泛的地理坐标系统和地球模型，由美国国防部制定和维护，主要用于全球定位系统（GPS）和各种地理信息系统（GIS）。

WGS-84 有几个关键常量用于定义参考椭球体的形状和尺寸，这些常量包括半长轴、半短轴、扁率和离心率等：

  • 半长轴 (A)：椭球体的赤道半径
    • 其值为 $ 6378137.0 $ 米。
  • 半短轴 (B)：椭球体的极半径
    • 其计算公式为： $ B = A*(1-F) $
    • 其值为 $ 6356752.3142 $ 米
  • 扁率 (F)：描述椭球体扁平程度的参数
    • 其计算公式为：$ F = \frac{A-B}{A}$
    • 其值为：1/298.257223563
  • 第一离心率 (E)：描述椭球体形状的一种参数，反映了椭球体的偏离程度
    • 其计算公式：$ E = \sqrt{1 - (\frac{B}{A}) ^2} $
    • 其值约为： 0.0818191908426
  • 第一离心率的平方 (E²)：减少计算所用
    • 其中值约为： 0.00669437999014

1.2 ECEF坐标系

ECEF（Earth-Centered, Earth-Fixed）坐标系是一种三维笛卡尔坐标系，用于表示地球上的位置。ECEF 坐标系也称为地心地固坐标系，它的原点位于地球质心，并且随着地球的自转而旋转。

  • X 轴：指向穿过地球赤道与本初子午线（通过格林尼治的子午线）交点的方向。
  • Y 轴：指向穿过地球赤道与东经90度子午线交点的方向。
  • Z 轴：指向北极方向，与地球自转轴一致。

ECEF 坐标系提供了一个统一的三维坐标框架，可以精确地表示地球表面和近地空间的任何位置。

1.3 ENU坐标系

ENU（East-North-Up）坐标系是一种局部笛卡尔坐标系，用于表示相对于某个参考点的三维位置。

  • ENU 坐标系的原点通常位于地球表面的某个参考点，该点的地理坐标为 (Latitude, Longitude, Altitude)。
  • E 轴（东向轴）：指向地平线的东方。
  • N 轴（北向轴）：指向地平线的北方。
  • U 轴（上向轴）：垂直向上，指向天空。

2. 坐标系转换

2.1 从 LLA 坐标到 ECEF 坐标

这里约定LLA的经度为 $\phi$，纬度为 $\lambda $，海拔为 $ h $, 选取 WGS-84 坐标系参数，$a$和$b$分别是是赤道半径（半长轴）和极半径（半短轴），$e^2 = 1 - \frac{b^2}{a^2}$是偏心率的平方, $f=1-\frac{b}{a} $ 是基准椭球体的极扁率。 $$ \begin{align} X& = (N(\phi) + h) cos \phi cos\lambda \\ Y& = (N(\phi) + h) cos \phi sin\lambda \\ Z& = (\frac{b^2}{a^2}N(\phi) + h)sin\phi \\ &=((1-e^2)N(\phi) +h)sin\phi \\ &=((1-f)^2N(\phi)+h)sin\phi \end{align} $$ 其中 $$ \begin{align} N(\phi) = \frac{a^2}{\sqrt{a^2cos^2\phi + b^2sin^2\phi}} = \frac{a}{\sqrt{1-e^2sin^2\phi}} \end{align} $$

2.2 从 ECEF 到 LLA坐标

将 ECEF 坐标 (X, Y, Z) 转换为经纬度和高度 (Latitude, Longitude, Altitude) 需要迭代计算。 $$ \begin{align} \lambda = atan2(Y,X) \end{align} $$

其中，atan2 是反正切函数

$\lambda$ 是唯一能直接算出的，其余的纬度和高度的转换所要涉及 N 的循环关系 $$ \begin{align} \frac{Z}{p}cot\phi = 1 - \frac{e^2N}{N+h} \\ h = \frac{p}{cos\phi} -N \end{align} $$ 其中，如公式$(6)$所示，N 的变化取决于 $\phi$ 的值

纬度和高度需要迭代求解， 例如，从第一个猜测 h≈0 开始，然后更新 N。

其流程是这样的：

  1. 猜测 h = 0，通过公式$(8) $ 可以推算 $cot\phi = \frac{(1-e^2)p}{Z} $ 求解出一个$\phi_1$
  2. 将 $\phi_1$ 带入公式 $(6)$，求解出 $N$
  3. 将 $N$ 带入公式$(9)$ 求解出 $h$
  4. 将 $h$ 带入公式$(8)$ 求解出一个新的 $\phi_2$
  5. 如果 $\phi_1$和$\phi_2$足够接近，则说明迭代出了一个正确的值，如果不是，则需要回到第2步继续迭代。

2.3 从 ECEF 到 ENU 坐标

ENU 坐标又称为站心坐标，6也就是局部坐标，这个坐标需要一个局部参考点。在实际例子中，这个参考点通常是雷达、基站这些地方。

如果基站位于 ${X_0,Y_0,Z_0}$，测站位于 ${X_1,Y_1,Z_1}$

$$ \begin{bmatrix} e\\ n\\ u \end{bmatrix}= \begin{bmatrix} -sin\lambda_0& cos\lambda_0 & 0\\ -sin\phi_0cos\lambda_0& -sin\phi_0sin\lambda_0 &cos\phi_0 \\ cos\phi_0cos\lambda_0& cos\phi_0sin\lambda_0 & sin\phi_0 \end{bmatrix} \begin{bmatrix} X_1 -X_0\\ Y_1-Y_0\\ Z_1-Z_0 \end{bmatrix} $$

(完)