半正矢公式

半正矢公式是一种根据两点的经度和纬度来确定大圆上两点之间距离的计算方法，在导航有着重要地位。它是球面三角学中“半正矢定理”公式的特例，该定理涉及了球面三角形的边和角。

尽管第一份英文版的半正矢表由詹姆斯·安德鲁在1805年印刷出版，但弗洛里安·卡乔里相信José de Mendoza y Ríos在1801年就使用过类似的术语。“半正矢”这个名字由詹姆斯·英曼在1835年创造。

这个公式正如其名，用半正矢函数表达，通过 $\operatorname{haversin}(\theta) = \sin^2\left(\frac {\theta}{2}\right)$ 而来。该公式可以用半正矢函数的任意倍数表达，如正矢函数（半正矢函数的两倍）。在计算机出现之前，为了计算简便，人们会利用对数来计算乘积和利用半正矢函数计算距离，所以在十九和二十世纪初的导航和三角测量书中包含了半正矢值表和对数表。现在，将该公式用半正矢函数表达也很方便，因为它能避免 $\sin^2$ 的系数。

## 半正矢公式

对于任何球面上的两点，圆心角的半正矢值可以通过如下公式计算：

$$\operatorname{hav}\left(\frac {d}{r}\right)=\operatorname{hav}(\varphi_2-\varphi_1)+\cos(\varphi_1) \cos(\varphi_2) \operatorname{hav}(\lambda_2-\lambda_1)$$

$$\operatorname{hav}(\theta) = \sin^2\left(\frac{\theta}{2}\right) = \frac{1 - \cos(\theta)}{2}$$

- $d$ 是两点之间的距离（沿大圆，见球面距离）。
- $r$ 是球的半径。
- $\varphi_1 \varphi_2$：点 1 的纬度和点 2 的纬度，以弧度制度量。
- $\lambda_1 \lambda_2$：点 1 的经度和点 2 的经度，以弧度制度量。

左边的等号 $\frac {d}{r}$ 是圆心角，以弧度来度量。

可以通过应用反半正矢函数（如果可以查到值）或通过使用反正弦函数来解出 $d$：

$$d = r\operatorname{archav}(h) = 2r\arcsin\left(\sqrt{h}\right)$$

$h = \operatorname{hav} \frac {d}{r}$，代入可得：

$$\begin{align} <pre> d &= 2r \arcsin\left(\sqrt{\operatorname{hav}(\varphi_2 - \varphi_1) + \cos(\varphi_1) \cos(\varphi_2)\operatorname{hav}(\lambda_2 - \lambda_1)}\right) \\ &= 2r \arcsin\left(\sqrt{\sin^2\left(\frac{\varphi_2 - \varphi_1}{2}\right) + \cos(\varphi_1) \cos(\varphi_2)\sin^2\left(\frac{\lambda_2 - \lambda_1}{2}\right)}\right) </pre> \end{align}$$

在使用这个公式时，必须确保 $h$ 不超过 1 （$d$ 只当 $h$ 的值在 0 到 1 间才有意义）。$h$ 仅在球面上两点连线过球心时才为 1，当测量精度有限时，会产生较大的计算误差。$d$ 如果大到接近圆周的一半长度 $\pi R$ 时（这个情况并不常见），一个小错误通常是一个大问题在。虽然还有其他的大圆距离公式，能够避免这个问题。上述公式有时候用反正切函数表达，但是 $h$ 接近 1 时这个问题仍然存在。

但是仍有对策来避免上述问题。类似的公式可以使用余弦函数表达（有时亦称为球面余弦定理，但是并非平面余弦定理）而不是半正矢，但如果这两点十分靠近（例如在地球上的一公里外）你可能得到的结果为 $\cos{\frac {d}{r}}=0.99999999$，导致一系列的误差。而由于半正矢式使用正弦函数，避免了这一问题。所以两个公式应该互补使用。

因为地球是一个不完美的球体，故当其应用于地球时，无论哪一个公式只是做一个近似测算，“地球半径” $R$ 在极点地区是 $6356.752$ 公里，在赤道地区为 $6378.137$ 公里。另外，半径曲率在极点处（$\approx 6399.594$ 公里）比在赤道处($\approx 6335.439$ 公里)大 1%，所以半正矢公式或者球面余弦定理是不能保证 0.5% 以内的误差。更准确的方法，应该是使用考虑地球离心率的Vincenty的公式或其他有关地理距离的论文所给出方法。

## 半正矢定理

<image>

给出一个单位球，一个在表面的球面三角形三个过三点 $(u,v,w)$ 的大圆所围出来的区域。如图，这个球面三角形的三边分别是 $a$（$\mathbf u$ 至 $\mathbf v$）， $b$ （$\mathbf u$ 到 $\mathbf w$）和 $c$ （$\mathbf v$ 至 $\mathbf w$）并且角 $C$ 对边 $c$ 那么有如下关系

$\operatorname{hav}(c)=\operatorname{hav}(a-b)+\sin(a) \sin(b) \operatorname{hav}(C)$

由于这是一个单位球，弧长 $a$、 $b$ 和 $c$ 与两端点同球心构成的圆心角（弧度制）相等。

为了从该定理中得到半正矢公式，只需要考虑一种特殊情况，即 $\mathbf u$ 是北极，而 $\mathbf v$ 和 $\mathbf w$ 是需要确定距离 $\mathbf d$ 两个点。在这个情况下， $a$ 和 $b$ 由 $\frac{\pi}{2} - \varphi_{1,2}$ 来计算（即 90° - 维度）， $C$ 为经度分隔 $\Delta \lambda$ 和 $c$ 为求得的$\frac{d}{R}$。结合 $\sin\left(\frac{\pi}{2} - \varphi\right) = \cos(\varphi)$，半正矢公式即可推导得出。

或者可以从球面余弦定理中得出

$$\cos(c)=\cos(a)\cos(b)+\sin(a)\sin(b)\cos(C)$$

当 $c$ 很小时，这个公式精度不高。所以我们可以用 $\cos(\theta) = 1 - 2\operatorname{\operatorname{hav}}(\theta)$ 来替换，并引用三角恒等式 $\cos(a - b) = \cos(a)\cos(b) + \sin(a)\sin(b)$，半正矢定理就能推出。

## 延伸阅读

- 美国人口普查局的 地理信息系统的问题(内容已被转移到 什么是最好的方式来计算两点间的距离?)
- R.W.Sinnott的"半正矢的优点", Sky and Telescope 68 (2), 159 (1984).
- 推导半正矢公式，Ask Dr.Math(Apr. 20–21,1999)。
- Romuald Ireneus'Scibor-Marchocki, 球三角形，初等几何(1997年)。
- W.Gellert,S.Gottwald,M.Hellwich，H.Kästner和H.Küstner的,The VNR Concise Encyclopedia of Mathematics,2nd ed., ch. 12(Van Nostrand Reinhold：纽约，1989年)。

## 外部链接

- 用的91种编程语言rosettacode.org实现半正矢公式 或 17种编程语言codecodex.com
- 其它在 C++、 C(mac os)、 Pascal、 Python、 Ruby、 JavaScript、 PHP、 Matlab、 MySQL中的实现。
- 由任何两者之间的纬度和经度计算距离、方位等。http://www.movable-type.co.uk