根据地球上任意两地的经纬度,可以计算它们在球面上的最短距离(Great-circle Distance / Orthodromic Distance) 及相对始末位置的方位角(Bearing)。本文基于地球是标准球体的假设,使用解三角形法(平面三角形和球面三角形)推导上述公式。

基本概念

纬度是地球上某点与球心的连线与地球赤道面所成的线面角(如图所示$\varphi$),北纬为正数,南纬为负数。

经度是地球上某点经线平面与本初子午面所成的二面角(如图所示的$\lambda$),东经为正数,西经为负数。

球面上两点P、Q的最短距离为过P、Q的大圆对应两点之间的劣弧的长度(如图实线所示$\overset{\frown}{PQ}$)。

方位角是从某点的指北方向线起,依顺时针方向到目标方向线之间的水平夹角(如图所示$\theta$)。

将上述物理量等效到单位球中,如下图所示:

球面距离公式及其推导

球面距离$d$的常用计算公式为Haversine公式[1,2]:

\begin{align*} a &= \sin^2\frac{\Delta\varphi}{2} + \cos\varphi_1 \cos\varphi_2 \, sin^2\frac{\Delta\lambda}{2}\\
\delta &= 2\,atan2\frac{\sqrt{a}}{\sqrt{1-a}} \\
d &= R\,\delta\\
\tag{I} \end{align*}

其中,$\varphi$是纬度,$\lambda$是经度,$R$是地球平均半径($R=6371km$)。$\Delta\varphi=\varphi_2-\varphi_1$,$\Delta\lambda=\lambda_2-\lambda_1$ 分别为纬度、经度的差值。

注意,代入计算中的经、纬度为换算后的弧度值。

推导

思路:求出$\overset{\frown}{PQ}$对应的球心角$\delta$后使用$d = R\,\delta$求出球面距离。

在球面三角形NPQ中直接使用球面三角形的余弦定理(Spherical law of cosines)[3,4]得到

\begin{equation*} \cos\delta = \sin\varphi_1 \sin\varphi_2 + \cos\varphi_1 \cos\varphi_2 \cos\Delta\lambda \tag{*}\label{eq:1} \end{equation*}

上式即可求出$\overset{\frown}{PQ}$对应的球心角,但是,当P、Q相邻很近(例如,数米之间)即$\delta$趋近于0时,$\cos\delta$将趋近于1(0.99999999)。 其后果是对于低浮点精度(low floating-point precision)的计算设备,可能由于舍入误差而造成错误的结果。当然,对于现代的64位浮点数运算 的计算机,一般不会出现上述问题。

为了使上述公式具有更好的计算兼容性,在$\Delta{OPQ}$中使用平面三角形余弦定理

\begin{align*} {PQ}^2 &= {OP}^2 + {OQ}^2 - 2\,{OP}\,{OQ}\,\cos\delta \\
&= 2 - 2\,( \sin\varphi_1 \sin\varphi_2 + \cos\varphi_1 \cos\varphi_2 \cos\Delta\lambda )\\
&= 2 - 2\,\left[ \sin\varphi_1 \sin\varphi_2 + \cos\varphi_1 \cos\varphi_2 \left( 1-2\,\sin^2\frac{\Delta\lambda}{2} \right) \right ]\\
&= 2 - 2(\sin\varphi_1 \sin\varphi_2 + \cos\varphi_1 \cos\varphi_2) + 4\,\cos\varphi_1 \cos\varphi_2\,\sin^2\frac{\Delta\lambda}{2} \\
&= 2\,\left( 1-\cos\Delta\varphi \right) + 4\,\cos\varphi_1 \cos\varphi_2 \,\sin^2\frac{\Delta\lambda}{2} \\
&= 4\,\sin^2\frac{\Delta\varphi}{2} + 4\,\cos\varphi_1 \cos\varphi_2 \,\sin^2\frac{\Delta\lambda}{2} \\
\end{align*}

令$a=(PQ/2)^2$,则上式写为:

\[a = \sin^2\frac{\Delta\varphi}{2} + \cos\varphi_1 \cos\varphi_2 \,\sin^2\frac{\Delta\lambda}{2}\]

进而得到球心角$\delta$为:

\[\delta = 2\,atan2\frac{\sqrt{a}}{\sqrt{1-a}}\]

方位角公式及其推导

方位角$\theta$的计算公式:

\[\theta = atan2\frac{\sin\Delta\lambda \cos\varphi_2}{\cos\varphi_1 \sin\varphi_2-\sin\varphi_1 \cos\varphi_2\,\cos\Delta\lambda} \tag{II}\]

推导

思路:球面三角形的正弦定理和余弦定理。

在球面三角形NPQ中分别使用正弦定理和余弦定理:

\begin{align*} \frac{\sin\theta}{\cos\varphi_2} &= \frac{\sin\Delta\lambda}{\sin\delta} \\
\sin\varphi_2 &= \sin\varphi_1 \cos\delta + \cos\varphi_1 \sin\delta \, \cos\theta \tag{**}\label{eq:2} \end{align*}

得到:

\begin{align*} \sin\theta &= \frac{\cos\varphi_2}{\sin\delta} \, \sin\Delta\lambda\\
\cos\theta &= \frac{\sin\varphi_2-\sin\varphi_1 \cos\delta}{\cos\varphi_1 \sin\delta} \end{align*}

两式相除,并将式\eqref{eq:1}代入得:

\begin{align*} \tan\theta &= \frac{\sin\Delta\lambda \cos\varphi_1 \cos\varphi_2}{\sin\varphi_2-\sin\varphi_1 \cos\delta} \\
&= \frac{\sin\Delta\lambda \cos\varphi_1 \cos\varphi_2}{\sin\varphi_2-\left( \sin^2\varphi_1 \sin\varphi_2 + \sin\varphi_1 \cos\varphi_1 \cos\varphi_2 \cos\Delta\lambda \right)}\\
&= \frac{\sin\Delta\lambda \cos\varphi_2}{ \cos\varphi_1 \sin\varphi_2 - \sin\varphi_1 \cos\varphi_2 \cos\Delta\lambda} \end{align*}

取反函数即可得到方位角$\theta$的表达式。

已知起始点、方位角及距离,求目标位置

已知起始点经纬度$P(\varphi_1,\lambda_1)$,方向角$\theta$及球面距离$d$,计算目标位置$Q(\varphi_2,\lambda_2)$的公式为:

\begin{align*} \varphi_2 &= asin \left( \sin\varphi_1 \cos\delta + \cos\varphi_1 \sin\delta \cos\theta \right)\\
\lambda_2 &= \lambda_1 + atan2 \frac{\sin\theta \sin\delta \cos\varphi_1}{\cos\delta-\sin\varphi_1 \sin\varphi_2} \tag{III} \end{align*}

其中,球心角$\delta = d/R$。

推导

思路:已知两边及夹角,解球面三角形问题。

根据球面三角形NPQ中的余弦定理,即式\eqref{eq:2}中的第二个式子,可直接得到Q点的纬度值$\varphi_2$。

根据球面三角形NPQ中的正弦定理,即式\eqref{eq:2}中的第一个式子,得到:

\[\cos\varphi_2 \sin\Delta\lambda = \sin\theta \sin\delta\]

根据式\eqref{eq:1}可知:

\[\cos\varphi_2 \cos\Delta\lambda = \frac{\cos\delta - \sin\varphi_1 \sin\varphi_2}{\cos\varphi_1}\]

以上二式相除:

\[\tan\Delta\lambda = \frac{\sin\theta \sin\delta \cos\varphi_1}{\cos\delta - \sin\varphi_1 \sin\varphi_2}\]

从上式即可直接解出Q点的经度$\lambda_2 = \lambda_1 + \Delta\lambda$。

参考资料