计算地球上两点距离（震中距）的Matlab函数（兼容度数和度分秒）及另外三种方法

写在前面

最近要计算震中距就想起自己之前写作业时写的这个代码，打开一看，写的太烂了，数字居然要以字符串格式输入，要输入好几个引号，赶紧修改一下。顺带介绍一下计算震中距的另外两种方法。

方法1: taup

其实也可以使用taup来计算震中距，不过还需要安装，也不支持度分秒格式，可以参阅Seisman的博客地震“学”软件-TauP，taup实际上是计算一维球状分层模型下地震震相的走时和路径的软件，不过也顺带输出了震中距。

使用如下，-evt和-sta选项分别指定震源和台站的纬度和经度。

taup time -evt 27.3397 -128.352 -sta 18.81 119.911

输出的第一列是震中距（度）。

Model: iasp91
Distance   Depth   Phase   Travel    Ray Param  Takeoff  Incident  Purist    Purist
  (deg)     (km)   Name    Time (s)  p (s/deg)   (deg)    (deg)   Distance   Name 
-----------------------------------------------------------------------------------
   99.40     0.0   Pdiff    824.08     4.439     13.39    13.39    99.40   = Pdiff
   99.40     0.0   PKiKP   1093.71     1.787      5.35     5.35    99.40   = PKiKP
   99.40     0.0   SKS     1463.80     4.972      8.64     8.64    99.40   = SKS  
   99.40     0.0   Sdiff   1517.52     8.323     14.57    14.57    99.40   = Sdiff
   99.40     0.0   SKiKS   1526.39     1.880      3.26     3.26    99.40   = SKiKS

需要指出的是，使用 -evt 和 -sta 选项时，taup time 会假设地球是完美球体，这会产生一定的误差，一般情况下可以接受，如果需要更精确的结果，可以使用ObsPy中的gps2dist_azimuth来计算。

方法2: ObsPy

ObsPy 的 gps2dist_azimuth 函数采用 WGS84 椭球：赤道半径 6378.1370 km、扁率约为 0.0033528106647474805。第一个输出为大圆距离（单位m），再使用kilometers2degrees 函数转为度数，不过此时又会遇到完美球体的假设。

使用gps2dist_azimuth(lat1, lon1, lat2, lon2)来计算大圆距离

from obspy import geodetics
distance = geodetics.gps2dist_azimuth(27.3397, -128.352, 18.81, 119.911)[0]
angle = geodetics.kilometers2degrees(distance*0.001)
print(angle)

输出

99.55074173891092

而使用taup和下文Matlab函数的计算结果都为 99.4

方法3: Mapping Toolbox的distance函数

distance为Matlab的Mapping Toolbox中的函数，可以直接计算两个坐标点的距离及方位角。用法如下：[arclen, az] = distance(Aw,Aj,Bw,Bj)

[arclen, az] = distance(27.3397, -128.352, 18.81, 119.911)

输出为：

arclen =

   99.4001


az =

  296.9690

返回的是震中距（°）及方位角。

方法4: 自己写的Matlab函数

参数

需要输入AB两点的坐标：

参数名	意义
Aw	A点纬度
Aj	A点经度
Bw	B点纬度
Bj	B点经度

R: 地球半径
d: 两点之间的距离（弧度）
latitude: 纬度
longitude: 经度
北纬为正 南纬为负 东经为正 西经为负

公式

1.球面余弦公式
因公式用到了大量余弦函数，当系统的浮点运算精度不高时，在求算较近的两点间的距离（几百米）时会有较大误差（64位系统一般无较大误差）。

$$d=R\times \arccos (\sin (Aw)\times \sin (Bw)+\cos (Aw)\times \cos (Bw)\times \cos (Aj-Bj))$$

2.Haversine公式
Haversine公式是球面余弦公式的一个变换，即使距离很小，也不会有较大误差。维基百科推荐使用Haversine公式。
$$d=2\times R\times \arcsin \left(\sqrt{{\sin }^2\left(\frac{Bw-Aw}{2}\right)+\cos \left(Aw\right)\times \cos (Bw)\times {\sin }^2\left(\frac{Bj-Aj}{2}\right)}\right)$$

函数

这里使用Haversine公式计算两个坐标间的距离，为了使函数更方便使用，这里经纬度坐标兼容度数格式和度分秒格式输入。Matlab中的sin、cos、asin都是弧度。

若输入度分秒格式的坐标，则需要以字符串类型输入。若输入度数格式的坐标，则以数字或字符串都输入都可以。例如：

[distace, angle] = calc_distance(27.3397, 128.352, 18.81, 119.911)
calc_distance(27.3397, -128.352, '18.48.36', '119.54.40')

function [distance, angle] = calc_distance(Aw,Aj,Bw,Bj)

% [distance (km), angle (degree)] = calc_distance(A_lat,A_lon,B_lat,B_lon)
% Calculating the epicentral distance of two points
% You need to enter the coordinates of two points A and B
% All coordinates can be entered in degrees or degrees, minutes and seconds
% The coordinates in the format of degrees, minutes and seconds should be input in string, 
% the coordinates in degrees can be input in string or number
% e.g., [distace, angle] = calc_distance(27.3397, 128.352, 18.81, 119.911)
% e.g., calc_distance(27.3397, -128.352, '18.48.36', '119.54.40'), return distance
% e.g., [~, angle] = calc_distance('-27.3397', 128.352, '18.48.36', 119.911), return angle
% Yuechu WU
% 12131066@mail.sustech.edu.cn
% updated 2022-05-30

Aw = num2str(Aw);
Aj = num2str(Aj);
Bw = num2str(Bw);
Bj = num2str(Bj);

% Radius of the earth (km)
R = 6378.1370;

if length(find(Aw=='.')) > 1
    DMS = find(Aw=='.');M = DMS(1);S = DMS(2);
    Aw = str2double(Aw(1:M-1))+str2double(Aw(M+1:S-1))/60+str2double(Aw(S+1:end))/3600;
else
    Aw = str2double(Aw);
end

if length(find(Aj=='.')) > 1
    DMS = find(Aj=='.');M = DMS(1);S = DMS(2);
    Aj = str2double(Aj(1:M-1))+str2double(Aj(M+1:S-1))/60+str2double(Aj(S+1:end))/3600;
else
    Aj = str2double(Aj);
end

if length(find(Bw=='.')) > 1
    DMS = find(Bw=='.');M = DMS(1);S = DMS(2);
    Bw = str2double(Bw(1:M-1))+str2double(Bw(M+1:S-1))/60+str2double(Bw(S+1:end))/3600;
else
    Bw = str2double(Bw);
end

if length(find(Bj=='.')) > 1
    DMS = find(Bj=='.');M = DMS(1);S = DMS(2);
    Bj = str2double(Bj(1:M-1))+str2double(Bj(M+1:S-1))/60+str2double(Bj(S+1:end))/3600;
else
    Bj = str2double(Bj);
end

Aw = Aw*pi/180;
Aj = Aj*pi/180;

Bw = Bw*pi/180;
Bj = Bj*pi/180;


distance = 2*R*asin(sqrt((sin((Bw-Aw)/2).^2)+cos(Aw)*cos(Bw)*(sin((Bj-Aj)/2)).^2));

angle = (distance/R)*180/pi;

end

输入两点坐标(Aw,Aj,Bw,Bj)调用该函数

[distace, angle] = calc_distance(27.3397, -128.352, 18.81, 119.911)

第一个返回值是震中距（km）,第二个返回值是震中距（°）。

distace =

   1.1065e+04


angle =

   99.4001