![855ed3ef080158f522164854b2b679e8.png](https://i-blog.csdnimg.cn/blog_migrate/70374515636711d3100388e6669733f3.jpeg)
维基百科推荐使用Haversine公式,理由是Great-circle distance公式用到了大量余弦函数, 而两点间距离很短时(比如地球表面上相距几百米的两点),余弦函数会得出0.999...的结果, 会导致较大的舍入误差。而Haversine公式采用了正弦函数,即使距离很小,也能保持足够的有效数字。 以前采用三角函数表计算时的确会有这个问题,但经过实际验证,采用计算机来计算时,两个公式的区别不大。 稳妥起见,这里还是采用Haversine公式。
![f46ddc29427161285b53b51e5a7aa4ab.png](https://i-blog.csdnimg.cn/blog_migrate/513d691ecb6612988056db926bc0eb32.jpeg)
其中
![83289704793679867238170ab8e31875.png](https://i-blog.csdnimg.cn/blog_migrate/23216d0c5019c2bd1c0a8830adea051b.jpeg)
- R为地球半径,可取平均值 6371km;
- φ1, φ2 表示两点的纬度;
- Δλ 表示两点经度的差值。
公式来历:
VERSINE(F)=1-cos(F)
![16d8e622a0877927640547dfb5b4bf50.png](https://i-blog.csdnimg.cn/blog_migrate/c4d86e586e89318d6f4f8ddf4603bc7c.png)
Haversine名字来历是Ha-VERSINE,即Half-Versine ,表示sin的一半的意思。
hav(A) = (1-cos(A))/2 = sin(A/2)* sin(A/2)
推倒过程:
如下一个半径为1 的圆,O是圆心,A、B是弦(chord)。角度AOB=theta。则角度AOC=theta/2。OC是垂直于AB的垂线(perpendicular)。AC长度是sin(theta/2),AB长度是2*sin(theta/2)。
![d091bb47f8078838bc7451b4ed79f64e.png](https://i-blog.csdnimg.cn/blog_migrate/239f0af0b928b6116cfaf0e99a02421b.png)
(图1)
如下地球图所示,假设半径R为1,O是球心,A (lat1,lon1) 和 B (lat2,lon2) 是我们感兴趣的2个点。2跟经度线 lon1,lon2相交于北极(north pole)N。EF所在的线是赤道(equator)。ACBD是平面上的等腰梯形的四个顶点(vertice)。AC和DB的弦(直线)在图上没有画出。CD的位置是:C (lat2,lon1) and D (lat1,lon2)。角度AOC是A点与C点的纬度差 dlat。角度EOF是经度E点和经度F点的差dlon。
![15610611e75bf83091d5a6ea8092d469.png](https://i-blog.csdnimg.cn/blog_migrate/47172ab979a93193e14c8f2f65e4edc4.jpeg)
(图2)
弦AC的长度,参照图1的方式,那么是AC=2*sin(dlat/2),弦BD也是一样的长度。
E、F 2个点是赤道上的2个点,它们的纬度是0。EF的距离是EF=2*sin(dlon/2)
A、D2个点所在的纬度是lat1。AD所在纬度的圆平面的半径是cos(lat1)。从A作一条垂线(perpendicular)到OE为AG,AO是球半径,则OG=cos(lat1),即A、D所在纬度圆圈的半径(AO`)。
这时候,AD的弦长AD= 2*sin(dlon/2)*cos(lat1),类似的可以推出CB的长度= CB=2*sin(dlon/2)*cos(lat2)
下面看一下如何求AB的长度,回到平面等腰梯形,如下图:
![b879c53d7bccb87768e931697c01a447.png](https://i-blog.csdnimg.cn/blog_migrate/6605e78f3dd71dabe41a7f3b7e73d3d6.png)
(图3)
AH是到CB的垂线(perpendicular),CH= (CB-AD)/2。
根据勾股定理(Pythagorean theorem): 【^2表示2的平方】
AH^2 = AC^2 - CH^2
= AC^2 - (CB-AD)^2/4
HB 的长度是HB=AD+CH = AD+(CB-AD)/2 = (CB+AD)/2,根据勾股定理得到:
AB^2 = AH^2 + HB^2
= AC^2 - (CB-AD)^2/4 + (CB+AD)^2/4
= AC^2 + CB*AD
根据前面球面上的求经纬距离的方式,我们已经得到 AC、AD和CB的长度,代入公式得到:
AB^2 = 4*(sin^2(dlat/2) + 4*cos(lat1)*cos(lat2)*sin^2(dlon/2))
假设中间值h 是AB长度一半的平方,如下
h = (AB/2)^2
= (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
(请参看代码里的h)
最后一步,是求得代表AB长度的角度AOB。参照图1的方式,我们可以知道
![8b30ff02ee0878487cec8c45b10cf0c1.png](https://i-blog.csdnimg.cn/blog_migrate/ab567982acdaca0a87e34dcf79aab98e.png)
(图4)
设AC=
![39d732fcc13603c0fddac16de1075000.png](https://i-blog.csdnimg.cn/blog_migrate/37ce083bcea6e3c6966571110b906b36.png)
,根据勾股定理(Pythagorean theorem)得到:
OC=
![57b250aeb978c50502907eb5200fd32a.png](https://i-blog.csdnimg.cn/blog_migrate/101740c5a1edded1bdb35597f0f56a48.png)
= sqrt(OA^2 - AC^2)
=
![564078d45b029fc9d692d6b091af9fa5.png](https://i-blog.csdnimg.cn/blog_migrate/0d174490d589b20520ffab56a351dc03.png)
= sqrt(1-a) // sqrt表示开根号
如果设c是角AOB的度数值。
tan(<AOC) = tan(c)= AC/OC = sqrt(a)/sqrt(1-a)
则:
c = 2 * arctan(sqrt(a)/sqrt(1-a)),
最后的AB真实距离,把地球半径带上就可以了。
distance = 2 * EARTH_RADIUS * c。
参考文献:https://www.cnblogs.com/softfair/p/distance_of_two_latitude_and_longitude_points.html