做建模或者研究空间数据，可能会遇到“根据经纬度计算两点间的球面距离”的问题，网上的资料很多，都是各种公式推导，但是一旦按公式编程计算，很可能得不到正确的距离。根本原因是在“角度-弧度的转化”与“软件语法”结合上出错了。
为此，我梳理了两种常用的计算方法，并用 Matlab 编程实现，保证能计算出正确的距离。
方法一：Great-Circle距离(基于球面余弦公式)
计算公式为：
其中，
为地球半径，
分别表示
点的经度、纬度；
分别表示
点的经度、纬度。
公式的具体推导过程略(可参考文献[1][2])
Matlab实现：
先定义一个角度转弧度的函数：
functionrad =D2R(theta)rad = theta*pi/180;
再定义计算Great-Circle距离的函数：
functiond =SphereDist(x,y,R)%根据两点的经纬度计算大圆距离(基于球面余弦公式)
%x为A点[经度, 纬度], y为B点[经度, 纬度]
if nargin < 3
R = 6378.137;
end
x = D2R(x);
y = D2R(y);
DeltaS = acos(cos(x(2))*cos(y(2))*cos(x(1)-y(1))+sin(x(2))*sin(y(2)));
d = R*DeltaS;
注：R 设为了默认参数，这里要吐槽一下 Matlab 在这方面的啰嗦，其他语言直接 SphereDist(x,y,R=6378.137) 就行了。
方法二：基于Haversine公式
前面的球面余弦公式中有
项，当两点间距离很短时(比如相距几百米的两点)，余弦函数会得出
的结果，会导致较大的舍入误差。Haversine 方法进行了某种变换消除了
项，能够避免该问题。
实际上，以现在电脑软件的一般计算精度，这已经不是问题。用计算器算可能有这个问题。
Haversine公式：
其中，
这里，
为地球半径；
表示两点的纬度；
表示两点经度的差值。
公式的具体推导过程略(可参考文献[2])
Matlab实现：
functiond =SphereDist2(x,y,R)%根据两点的经纬度计算大圆距离(基于Haversine公式)
%x为A点[经度, 纬度], y为B点[经度, 纬度]
if nargin < 3
R = 6378.137;
end
x = D2R(x);
y = D2R(y);
h = HaverSin(abs(x(2)-y(2)))+cos(x(2))*cos(y(2))*HaverSin(abs(x(1)-y(1)));
d = 2 * R * asin(sqrt(h));
functionh =HaverSin(theta)h=sin(theta/2)^2;
end
end
注：嵌套了一个 HaverSin 函数。
测试：
%计算北京到天津的球面距离
SphereDist([116.4,39.9], [117.2, 39.1])
SphereDist2([116.4,39.9], [117.2, 39.1])
运行结果：
参考文献：
[1] 根据两点的经纬度求方位角和距离，等
原创作品，转载请注明。