球面距离$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$。