测地计算

介绍

考虑一个具有赤道半径的旋转椭球 \(a\) ，极半轴 \(b\) ，并变平 \(f=(a-b)/a\) . 椭球表面上的点以它们的纬度为特征 \(\phi\) 经度 \(\lambda\) . （注意这里的纬度是指 地理纬度 ，椭球体法线与赤道面之间的角度）。

椭球体上两点间的最短路径 \((\phi_1,\lambda_1)\) 和 \((\phi_2,\lambda_2)\) 称为测地线。它的长度是 \(s_{{12}}\) 从点1到点2的测地线具有正方位角 \(\alpha_1\) 和 \(\alpha_2\) 在两个端点。在这个数字中，我们有 \(\lambda_{{12}}=\lambda_2-\lambda_1\) .

通过要求任何足够小的线段是最短路径，测地线可以无限延伸；测地线也是曲面上最直的曲线。

测地线问题的求解

传统上考虑两个测地线问题：

• 直接问题-给定 \(\phi_1\) ， \(\lambda_1\) ， \(\alpha_1\) ， \(s_{{12}}\) ，确定 \(\phi_2\) ， \(\lambda_2\) ， \(\alpha_2\) .

• 反问题-给定 \(\phi_1\) ， \(\lambda_1\) ， \(\phi_2\) ， \(\lambda_2\) ，确定 \(s_{{12}}\) ， \(\alpha_1\) ， \(\alpha_2\) .

Proj合并了 C library for Geodesics 。这个库提供了解决测地线正问题和反问题的例程。保持完全双精度精度，前提是 \(\lvert f\rvert<\frac1{{50}}\) 。

测地例程的接口在两个方面与项目的其余部分不同：

• 角度（纬度、经度和方位角）以度（而不是弧度）为单位；

• 椭球体的形状由压扁度指定 \(f\) ；这可以是负数以表示长椭球；设置 \(f=0\) 对应于一个球体，在这种情况下，测地线变成一个大圆。

PROJ还包括一个命令行工具， 格奥德 （1），用于执行简单的测地计算。

其他属性

这些例程还计算其他几个感兴趣的量

• \(S_{{12}}\) 是从点1到点2的测地线与赤道之间的面积；即逆时针测量的带角四边形的面积 \((\phi_1,\lambda_1)\) ， \((0,\lambda_1)\) ， \((0,\lambda_2)\) 和 \((\phi_2,\lambda_2)\) . 以米为单位² .

• \(m_{{12}}\) 当初始方位角受到 \(d\alpha_1\) （弧度）然后第二个点被 \(m_{{12}}\,d\alpha_1\) 在垂直于测地线的方向上。 \(m_{{12}}\) 以米为单位。在曲面上，缩减长度服从对称关系， \(m_{{12}}+m_{{21}}=0\) . 在平面上，我们有 \(m_{{12}}=s_{{12}}\) .

• \(M_{{12}}\) 和 \(M_{{21}}\) 是测地尺度。如果两条测地线在点1处平行并相隔一小段距离 \(dt\) ，然后它们被一段距离隔开 \(M_{{12}}\,dt\) 在第2点。 \(M_{{21}}\) 定义类似（测地线在点2处彼此平行）。 \(M_{{12}}\) 和 \(M_{{21}}\) 是无量纲量。在平面上，我们有 \(M_{{12}}=M_{{21}}=1\) .

• \(\sigma_{{12}}\) 辅助球面上的弧长。这是一个将问题转化为球面三角问题的构造。从一个赤道交叉点到下一个赤道交叉点的球面弧长度总是 \(180^\circ\) .

如果点1、2和3位于单个测地线上，则以下加法规则适用：

• \(s_{{13}}=s_{{12}}+s_{{23}}\) ，

• \(\sigma_{{13}}=\sigma_{{12}}+\sigma_{{23}}\) ，

• \(S_{{13}}=S_{{12}}+S_{{23}}\) ，

• \(m_{{13}}=m_{{12}}M_{{23}}+m_{{23}}M_{{21}}\) ，

• \(M_{{13}}=M_{{12}}M_{{23}}-(1-M_{{12}}M_{{21}})m_{{23}}/m_{{12}}\) ，

• \(M_{{31}}=M_{{32}}M_{{21}}-(1-M_{{23}}M_{{32}})m_{{12}}/m_{{23}}\) .

多条最短测地线

通过求解反问题找到的最短距离（显然）是唯一定义的。然而，在一些特殊情况下，有多个方位角产生相同的最短距离。以下是这些案例的目录：

• \(\phi_1=-\phi_2\) （两个点都不在杆子上）。如果 \(\alpha_1=\alpha_2\) ，测地线是唯一的。否则有两条测地线，第二条测地线通过设置 \([\alpha_1,\alpha_2]\leftarrow[\alpha_2,\alpha_1]\) ， \([M_{{12}},M_{{21}}]\leftarrow[M_{{21}},M_{{12}}]\) ， \(S_{{12}}\leftarrow-S_{{12}}\) . （这发生在经度差接近时 \(\pm180^\circ\) 对于扁椭球体。）

• \(\lambda_2=\lambda_1\pm180^\circ\) （两个点都不在杆子上）。如果 \(\alpha_1=0^\circ\) 或 \(\pm180^\circ\) ，测地线是唯一的。否则有两条测地线，第二条测地线通过设置 \([\alpha_1,\alpha_2]\leftarrow[-\alpha_1,-\alpha_2]\) ， \(S_{{12}}\leftarrow-S_{{12}}\) . （当 \(\phi_2\) 就在附近 \(-\phi_1\) 对于长椭球体。）

• 点1和点2在相反的极点。通过设置可以生成无限多条测地线 \([\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,-\delta]\) ，用于任意 \(\delta\) . （对于球体，当第1点和第2点是反足的时，此处方适用。）

• \(s_{{12}}=0\) （重合点）。通过设置可以生成无限多条测地线 \([\alpha_1,\alpha_2]\leftarrow[\alpha_1,\alpha_2]+[\delta,\delta]\) ，用于任意 \(\delta\) .

多边形的面积

测地线面的面积可以通过求和来确定 \(-S_{{12}}\) 对于多边形的连续边 (\(S_{{12}}\) 被否定，从而使多边形的顺时针遍历产生正区域)。但是，如果面环绕极点，则总和必须调整 \(\pm A/2\) ，在哪里 \(A\) 是整个椭球体的面积，其中选择了放置结果的符号 \((-A/2, A/2]\) 。

背景

The algorithms implemented by this package are given in [Karney2013] (addenda) and are based on [Bessel1825] and [Helmert1880]; the algorithm for areas is based on [Danielsen1989]. These improve on the work of [Vincenty1975] in the following respects:

• 结果精确到地球椭球体的四舍五入（距离误差小于15纳米，而文森蒂的误差为0.1毫米）。

• 反问题的解总是可以找到的。（Vincenty的方法无法收敛到几乎相反的点。）

• 这些例程计算测地线的微分和积分性质。例如，这允许计算测地多边形的面积。

附加的背景资料提供了地理的CLIB的 geodesic bibliography ，维基百科的文章“ Geodesics on an ellipsoid", and [Karney2011] (errata) .