半正矢公式（Haversine公式）

前两天机缘巧合之下，需要批量求两个已知经纬度的点之间的距离。

心想着这年头，电脑辅助的计算已经惠及生活各个角落，Excel里说不定都自带这个函数。

甚至天真地想用勾股定理直接套进去，后来才发现并没有简单到那个程度。

在维基百科里，关于相关的计算有非常清晰、简练的成果，那就是半正矢公式（Haversine公式）：

Haversine formula

公式非常简洁有力，但没有写明推导过程，网上能搜到不少相关的说明，但不少都是不同的来源拼凑而成，图片也不太清楚。

于是自己再过一遍，作为学习笔记：

haversin函数，针对任意角的一种函数。

然后定义Θ=d/R，其中d是两地间距离，R为地球半径。

就可以得到如下的公式：

这就是半正矢公式（Haversine公式）了。

公式里：

φ1和φ2表示两点的纬度；

λ1和λ2表示两点的经度。

简洁，有力。

这个公式想直接放进Excel里硬算距离的话，得做一些变形：

这个就可以直接带进Excel里了。

接下来记录一下学习过程中得到的推导过程：

这是一个没有颜色的地球模型。

已知A(φ1，λ1)，B(φ2，λ2)，地球半径R。

在图中，A，B为地球表面已知经纬度的两点。N为北极点，S为南极点。弧NHS为本初子午线。弧HEF为赤道。弧NADES为经过A点的经线，弧AC为经过A点的纬线。弧NCBFS为经过B点的经线，弧DB为经过B点的纬线。

线段OA’是过O点的线段AD的垂线，又因为AO=DO，所以线段OA’是线段AD的中垂线；同理线段A’’O’’是线段AC的中垂线；线段B’O’是线段BD的中垂线。另，线段AG是过A点的线段BD的垂线，这根垂线是等腰梯形ACBD的高。

φ1=∠AOE=∠O’’AO，λ1=∠EOH，φ2=∠BOF=∠DOE=∠O’DO，λ2=∠FOH。

除了红线，其他线可以忽略。

这个算法的思路，是计算线段AB的长度LAB，再结合线段AO和BO的长度都等于地球半径R，可以反向计算∠AOB的大小，并计算出弧AB的长度。

加一句题外话，在AB距离较小，∠AOB趋近于零的时候，弧AB的长度近似等于线段AB也不是不可以。

还是只看红线。

为了计算线段AB的长度LAB，就得借助一套辅助线，就是等腰梯形ACBD，线段AB是等腰梯形ACBD的对角线。

而用一系列的三角函数和勾股定理，可以根据A、B点的坐标计算出线段AD、BD、AC、BC的长度，其中AD=BC。

继续只看红线。

先求线段AD也就是线段BC的长度LAD，这里是用三角函数来求的。

∠AOE=φ1，∠DOE=∠BOF=φ2。

所以∠AOD=(φ1-φ2)。因为OA’是线段AD的中垂线，所以∠AOA’=∠AOD/2。LAA’=R*Sin(∠AOA’)=R*Sin(∠AOD/2)=R*Sin((φ1-φ2)/2)。

LAD=2* LAA’=2R*Sin(∠AOD/2)=2R*Sin((φ1-φ2)/2)。

只看红线。

然后来计算线段BD的长度LBD，要算BD还是和前面算LAD一样的，用线段DO’的长度乘以Sin(∠DO’B’)，即LBD=2*LDO’*Sin(∠DO’B/2)。

这里关于∠DO’B要多说一句：因为弧NAEDS和弧NCDFS都是经线，而NS是穿南北极的地轴，所以∠DO’B=∠EOF=∠AO’’C=λ2-λ1。

那么先算线段DO’的长度LDO’，LDO’=R*Cos(∠DOE)=R*Cos(φ2)。

结合以上两步，得出LBD=2R*Cos(φ2)*Sin((λ2-λ1)/2)。

红线。

这时再来求线段AC的长度LAC，思路和求LBD是一模一样的，就不再写一遍了。

LAC=2R*Cos(φ1)*Sin((λ2-λ1)/2)。

红。

好了，这个时候安心回到等腰梯形ACBE上。

经过前面的一系列计算，已经得出了以下数据：

LAC=2R*Cos(φ1)*Sin((λ2-λ1)/2)。

LBD=2R*Cos(φ2)*Sin((λ2-λ1)/2)。

LAD=2R*Sin((φ1-φ2)/2)。

一个已知四条边长度的等腰梯形的对角线长度，应该有非常成熟的公式，但既然来写笔记，就一并过一遍。

先算线段AG的长度LAG，而要算LAG，就在直角三角形ADG里面用勾股定理，LAG^2=LAD^2-LDG^2。那就得先算线段DG的长度LDG，LDG=(LDB-LAC)/2。

于是LAG^2=(2R*Sin((φ1-φ2)/2))^2-(R*(Cos(φ2)-Cos(φ1))*Sin((λ2-λ1)/2)^2。

这里就先不忙开方了，因为为了算线段AB的长度LAB，又会用到一个勾股定理，LAB^2=LAG^2+LBG^2。

先线段BG的长度，LBG=LBD-LDG=(LDB+LAC)/2，把这个带进上面的勾股定理里面，就可以得到：

LAB^2=(2R*Sin((φ1-φ2)/2))^2-(R*(Cos(φ2)-Cos(φ1))*Sin((λ2-λ1)/2)^2+(R*(Cos(φ2)+Cos(φ1))*Sin((λ2-λ1)/2)^2=4R^2*(Sin((φ1-φ2)/2)^2+Cos(φ1)*Cos(φ2)*Sin((λ2-λ1)/2))

LAB=2R*(Sin((φ1-φ2)/2)^2+Cos(φ1)*Cos(φ2)*(Sin((λ2-λ1)/2))^2)^0.5

到这步，对于比较小角度的AB两点直线距离已经算是得出来了，这一路上都没有考虑正负的问题，因为一直都是二次方开路，一直到了最后这一步才有一个开方，而开方的括号里，虽然有可能出现负值，但这些负值都有一个二次方进行洗礼，所以都不用管正负问题了。当然这些全是建立在没有考虑南北纬东西经的基础上的，如果要考虑，那么在计算之前还得换算一下。

因为有计算器的存在，LAB的公式虽然长，但只要确定它是对的，就可以硬算。然后要算地球表面弧长的话，也不是那么费劲。这是一个已知扇形半径和弦长，求弧长的公式：首先算∠AOB的大小：∠AOB=2*ArcSin(LAB/2R)=2*ArcSin((Sin((φ1-φ2)/2)^2+4*Cos(φ1)*Cos(φ2)*(Sin((λ2-λ1)/2))^2)^0.5)。

弧AB的长度=∠AOB*R。

就是前面那个可以带进Excel硬算的半正矢公式(Haversine公式)了。

【完】