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

2023-06-23 20:41| 来源: 网络整理| 查看: 265

目录写在前面方法1: taup方法2: ObsPy方法3: Mapping Toolbox的distance函数方法4: 自己写的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两点的坐标：

参数名意义AwA点纬度AjA点经度BwB点纬度BjB点经度

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

公式

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

d = R × arccos ⁡ ( sin ⁡ ( A w ) × sin ⁡ ( B w ) + cos ⁡ ( A w ) × cos ⁡ ( B w ) × cos ⁡ ( A j − B j ) ) d=R\times \arccos(\sin(Aw)\times\sin(Bw)+\cos(Aw)\times\cos(Bw)\times\cos(Aj-Bj)) d=R×arccos(sin(Aw)×sin(Bw)+cos(Aw)×cos(Bw)×cos(Aj−Bj))

2.Haversine公式 Haversine公式是球面余弦公式的一个变换，即使距离很小，也不会有较大误差。维基百科推荐使用Haversine公式。 d = 2 × R × arcsin ⁡ ( sin ⁡ 2 ( B w − A w 2 ) + cos ⁡ ( A w ) × cos ⁡ ( B w ) × sin ⁡ 2 ( B j − A j 2 ) ) 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) d=2×R×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 % [email protected] % 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

【本文地址】

公司简介

联系我们