写在前面

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

方法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   = Pdiff99.40     0.0   PKiKP   1093.71     1.787      5.35     5.35    99.40   = PKiKP99.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   = Sdiff99.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.4001az =296.9690

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

方法4: 自己写的Matlab函数

参数

需要输入AB两点的坐标：

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

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

公式

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

d=arccos⁡(sin⁡(Aw)×sin⁡(Bw)+cos⁡(Aw)×cos⁡(Bw)×cos⁡(Aj−Bj))d=arc\cos(\sin(Aw)\times\sin(Bw)+\cos(Aw)\times\cos(Bw)\times\cos(Aj-Bj)) d=arccos(sin(Aw)×sin(Bw)+cos(Aw)×cos(Bw)×cos(Aj−Bj))

2.Haversine公式
Haversine公式是球面余弦公式的一个变换，即使距离很小，也不会有较大误差。维基百科推荐使用Haversine公式。
d=2×arcsin⁡(sin⁡2(Bw−Aw2)+cos⁡(Aw)×cos⁡(Bw)×sin⁡2(Bj−Aj2))d=2\times arc\sin\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) d=2×arcsin(sin2(2Bw−Aw)+cos(Aw)×cos(Bw)×sin2(2Bj−Aj))

函数

这里使用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-30Aw = num2str(Aw);
Aj = num2str(Aj);
Bw = num2str(Bw);
Bj = num2str(Bj);% Radius of the earth (km)
R = 6378.1370;if length(find(Aw=='.')) > 1DMS = 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;
elseAw = str2double(Aw);
endif length(find(Aj=='.')) > 1DMS = 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;
elseAj = str2double(Aj);
endif length(find(Bw=='.')) > 1DMS = 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;
elseBw = str2double(Bw);
endif length(find(Bj=='.')) > 1DMS = 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;
elseBj = str2double(Bj);
endAw = 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+04angle =99.4001

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

php 计算两点时间距离,PHP计算地球上两点之间的距离（示例详解）
给定经度和纬度,求地球上两点之间的距离.首先我们需要了解该问题的解决思路,然后再用PHP代码来实现计算. 此问题可以用半正矢(haversine)公式求解: 大圆距离或正交距离是球面(或地球表面)上两 ...
根据经纬度计算地球上两点之间的距离——Haversine公式介绍及计算步骤
目录摘要 1.半正矢公式(Haversine Formula)介绍 2.半正矢公式应用 3.半正矢公式计算 3.1 主要思路 3.2 计算步骤 3.2.1 平面向量计算方法 3.2.2 空间向量计算 ...
给定经纬度计算距离_根据经纬度计算地球上两点之间的距离js实现代码
利用JS实现的根据经纬度计算地球上两点之间的距离最近用到了根据经纬度计算地球表面两点间距离的公式,然后就用JS实现了一下. 计算地球表面两点间的距离大概有两种办法. 第一种是默认地球是一个光滑的球面 ...
Android 计算地球上两点的距离
private static double EARTH_RADIUS = 6378.137;//地球半径 @Override protected boolean isRouteDisplayed() ...
计算地球上两点之间的俯仰角和方位角
源代码:https://gitee.com/gnoyuin/jiaodu https://github.com/niuyong/jiaodu
php 地图两点距离计算,计算地图上两点间的距离PHP类
计算地图上两点间的距离,使用的是谷歌地图 class GeoHelper { /** * @param int $lat1 * @param int $lon1 * @param int $lat2 ...
如何计算地球上两点的距离（附公式推导）
前段时间,看了一些电子围栏的算法,对其中一段计算球面上两点距离的代码有些不解,然后找了一下相关算法,在维基百科的大圆距离词条中记录了相关的计算公式,大致思路就是求出这两点间的弧长对应的圆心角的余弦或正 ...
利用CASS使用三种方法计算两期土方
软件及数据准备:南方CASS9.1:开挖前后的坐标高程数据dat格式. 计算方法:方格网网法.三角网法.断面法. 首先,绘出计算的范围线,注意开挖前后范围线一般相同.为了方便后续的三种方法计算,先生成 ...
调节e18-d80nk的测量距离_线缆太长负载太远，负载端电压难测量？三种方法帮你搞定...
关键词:电缆, 线缆阻抗, 负载电流, 开尔文检测在产品设计,研发过程中,我们常会碰到稳压器与负载分离的情况.此时,如果电缆线太长,线缆阻抗无法忽略.负载电流太大,这些都会使得配电线上的压降增大,从 ...

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

目录