第八章地震定位与地球速度反演.doc

个人收集整理 勿做商业用途
Error! Bookmark not defined.Error! Bookmark not defined. 第八章 走时数据的反演 89
8.1 球层介质的速度反演 90
8。1.1 拐点法（古登堡法) 90
8.1。2 Hergloz – Wiechert方法 91
8.2 单台地震定位及Inglada法定位 100
8.2.1单台定位法 100
8.2.2三站求取震中位置 105
8.2。3 采用三台站求解地震位置的计算机算法 106
8.3 地震定位结果与台站分布的关系 109
8。4 迭代定位方法 113
8.4.1 迭代定位算法的基本思路 113
8。4。2迭代方法的具体实现 114
8。5 相对定位方法 131
第八章 走时数据的反演

在前两章我们根据已知的速度结构对射线追踪和走时曲线的计算问题进行探讨，推导了波速只随深度变化的一维速度模型的射线追踪的表达式。三维结构的射线追踪虽然较复杂，但它遵循的原则是类似的。现在我们来研究根据观测到时数据得到速度结构和地震位置的问题。前面两章是已知速度结构和地震位置，如何计算地震波走时的问题，而本章是根据地震波走时得到地震结构和震源位置的问题,是前面两章的反问题.在这个问题的研究中，地震学家通常把问题进行简化：第一种简化是对于震源位置已经清楚的地震，根据地震波走时数据求解速度结构，第二种简化是已经得出了该地区的速度结构，根据地震波到时求解地震位置.限于篇幅，对于第一种简化,本书仅讲授球层介质一维速度结构的反演.对于第二种简化，我们讲授地震定位的基本方法。
8。1 球层介质的速度反演
8。1.1 拐点法(古登堡法）
这个方法是利用走时曲线的拐点处与从震源处水平方向射出的射线相对应，据此可求出震源处的波速.
根据第七章球层介质中的Snell定律有：
(Error! Bookmark not defined.-Error! Bookmark not defined.—1）
其中h代表震源处的参数，都是常数,从震源向不同方向射出的射线，值不同，相应的射线参数p也不同，当=90°时，达到极大值，即射线参数也达到极大值，因此有
Error! Bookmark not defined.（Error! Bookmark not defined.-1-1）
根据球层介质中的Bendorff定律(7-5-3)式,，有
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.-1-Error! Bookmark not defined.）
而对应于走时曲线的拐点.因此走时曲线的拐点与从震源处水平射出的射线相对应.
这样，我们只要求得某地震的震源深度及有观察分析归纳得到其走时曲线,找出走时曲线的拐点，并确定该点的走时曲线的斜率，则
Error! Bookmark not defined.（Error! Bookmark not defined.—1—Error! Bookmark not defined.)
而对应的M点有，因此有
\＊ MERGEFORMAT (Error! Bookmark not defined.—Error! Bookmark not defined.—Error! Bookmark not defined.)
其中，h为震源深度,为走时曲线拐点M点的视速度（参看第七章球层介质的Benndorf定律），R为地球半径。
此方法原理清楚、方法简单、计算方便。但由于地震发生的深度大约在0—700km的范围内，采用这种方法只能求出0-700km处的波速。由于拐点不易找准,得到的精度较差,它又要求对应每个深度的地震就要总结出一条走时曲线，因而资料分析工作量很大。
8.1.2 Hergloz – Wiechert方法
考虑球面介质，速度随半径r 变化[v( r ）],根据7.2节中得出震中距的表达式有
\＊ MERGEFORMAT （0—Error! Bookmark not defined.-Error! Bookmark not defined.）
其中=，为射线转折点距地心的距离，p=，R为地球半径。改变积分变数,则
, （0-Error! Bookmark not defined.—4）
其中。
在此介绍一积分式，
， Error! Bookmark not defined.(Error! Bookmark not defined.-1-4)
其中，为转折半径为的射线参数,所以式表示从转折半径为0至到转折半径为的所有射线的积分.将式作用于震中距表达式式的两边，得到
\＊ MERGEFORMAT (Error! Bookmark not defined.-Error! Bookmark not defined.—5）
方程的左边变为,考虑到，可以用分部积分的形式给出
\＊ MERGEFORMAT (Error! Bookmark not defined.-1-Error! Bookmark not defined.)
式的第一项代入积分上下限,由于cosh—1（1)=0，且在p=，半径为R，Δ=0，因此第一项消失.仅剩下第二项,可写成：
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.—Error! Bookmark not defined.—6)
的右侧为一双重积分式，改变积分变数,得到：
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.-Error! Bookmark not defined.—Error! Bookmark not defined.)
这里采用了不定积分公式：
整理变为
\＊ MERGEFORMAT Error! Bookmark not defined.(Error! Bookmark not defined.—1-6)
与式结合，得出

\＊ MERGEFORMAT Error! Bookmark not defined.(Error! Bookmark not defined.—Error! Bookmark not defined.—6)
令,可以得到：
\＊ MERGEFORMAT (Error! Bookmark not defined.-1-Error! Bookmark not defined.)
这就是射线最深点的半径的表达式。此表达式可以根据观测走时曲线t（）得出。如果t()为一平滑曲线,我们可由其切线斜率得到p()（由于p=)。式中的,为在震中距的走时曲线斜率.S1的积分式由震中距Δ从0到Δ1的p值，得到S1后，根据式求出，即其最深转折点。根据转折点为的，得出其对应此深度的速度
\＊ MERGEFORMAT Error! Bookmark not defined.(Error! Bookmark not defined.-1-7）
依此类推，得出速度与深度的分布图。
下面归纳一下求速度分布的具体步骤
1、 做出0-范围内的走时曲线，即t-曲线，震源要取在地表上,否则要进行修正；在该步骤中，通常在计算机上利用三次自然样条函数对一元函数进行成组插值及微分,从而使得走时曲线具有较好的平滑性能，采用该方法能保证所插值的函数及其一阶导数、二阶导数连续。
2、 由t-曲线求曲线上各点的斜率，从而做出曲线；该步骤求出p随震中距的分布。
3、 在0-范围内任取，取出点对应的；这就是。
4、 将曲线中0-段除以常数即，得到，即曲线，从而求出曲线；这就是随震中距的变化。
5、 在0—范围内求出曲线下的面积,即求,这就是；在计算积分的过程中，通常采用Simpson积分数值公式进行。
6、 根据式求出，然后根据式得到所对应的速度。
7、 取不同的值，重复上述步骤，从而求出不同的值所对应的速度值；
以一个具体地区的走时曲线为例子：
采用的数据是北京台网、大连、长春、牡丹江等26个台站的记录，并选用从北京地区到阿留申群岛西端的300多个地震，经过震源深度及剥壳校正后(地球平均厚度35公里）。得到的北京—萨哈林一线的P波走时曲线。数据如表8—1—1为北京-萨哈林剖面P波走时观测值。
表8-1-1北京-萨哈林剖面P波走时观测值
∆（°)
t(s）
∆（°）
t（s)
0。0
0
20.0
265。2
1.0
14。4
21.0
276
2。0
28.7
22。0
286。8
3.0
42.8
23。0
297。6
4。0
56.8
24。0
308。3
5。0
70。8
25.0
317。5
6.0
84.8
26。0
326.5
7.0
99。0
27.0
335.5
8。0
113.4
28。0
344。5
9。0
127。8
29。0
353.3
10。0
141.8
30。0
362。1
11。0
154。8
32。0
379.7
12.0
167。8
34。0
397.3
13.0
180.8
36。0
414.1
14。0
193.8
38.0
430.9
15.0
206.8
40。0
447.7
16。0
219。8
42。0
464.2
17。0
232.8
44.0
479.5
18.0
243。6
46.0
494.5
19.0
254.4
48。0
509。5

然后将走时曲线用Herglotz-Wiechert方法,利用MATLAB中的三次样条插值和辛普森积分方法，编制成MATLAB程序。程序如下：

clc
clf
clear all
load bj_shl_zs.txt %调用数据
x0=bj_shl_zs(:,1)'；％震中距
y0=bj_shl_zs(:，2）';％走时
dh=0。2；
x=0.2:dh:47.8;
n=length(x）；
［y,dy1，dy2］=splineinsert（x0,y0，x）; ％对数值进行样条插值,得到插值点的y值，一阶导数dy1,二阶导数dy2
figure（1) %实际观测曲线
plot（x,y,’r*’）;
grid on
title（’实际观测曲线’)
xlabel(’震中距’);
ylabel(’走时'）;

R=6371—35；
Indx=1；
r1=zeros(1，(length（3：2：n）)）； %射线能够达到的最低点的半径
v=zeros(1，(length(3:2:n）））; %射线最低点的速度
for ii=3：2:n
ksi1=dy1（ii); ％视为常数
f=acosh(dy1(1:ii）/ksi1）；
Int=dh/3＊（f(1)+f(ii)+2＊sum(f(3：2：ii)）+4*sum（f(2:2：ii)）)*pi/180；
％采用Simpson 公式进行积分
r1（Indx)=R/exp（Int/pi)； ％采用公式求解射线最深点的半径
v(Indx）=deg2rad（r1(Indx)/ksi1）；
%采用球层Snell定律得到最深点的速度；deg2rad(）函数的作用是:将角度转换为弧度 ，利用公式
Indx=Indx+1； ％序号+1
%pause
end
plot(v，r0—r1，'.’) ％绘制速度和深度剖面图
set(gca,'ydir’,’reverse’） ％将y轴的方向进行翻转

利用上述程序调用表8—1-1中的北京—萨哈林地区的P波走时曲线的数据。对观测曲线进行上述程序来反演即可得到速度剖面如图8-1-2北京—萨哈林地区P 波走时观测反演速度剖面.

8—1-1北京—萨哈林地区P 波走时观测反演速度剖面
对反演数据进行重新插值处理,选择出合适的数据，分别选择如下的数据(km):20,45，85,125,145,165，265，365，405，445,495,565，705,875,965,1025,1165
在上述MATLAB函数后加入如：
x=［20,45，85,125,145,165,265,365,405，445，495，565,705，875,965,1025,1165]；
[y，dy1，dy2］=splineinsert（r0—r1，v,x);
［r0-r1' v']
[x’，y']
并利用得到的数据进行绘图得到图8—1-2定点插值的反演数据图

图8—1-2定点插值的反演数据图

当有低速区时，这些公式是有问题的。在这种情况下，不连续，没有唯一的解。对这种情况,传统的做法是设想在低速区存在若干有不同速度的均匀层。由于没有射线在这些层里转折,因此低速区的这些层可以任意移来移去，穿过低速区的射线的积分的走时和震中距不会改变。
不论解析多么完全，在地震学中很少用赫格劳兹—维歇尔(HW）公式，这至少有两个原因：首先，HW假定我们有已知的连续的曲线，而实际上我们总是只有有限的走时点。这意味着走时曲线必须在这些数据之间进行内插，这些内插方案的不同导致于有不同的速度剖面.实际上,会有无数个稍有不同的速度模型,它们都适合于有限数目的点.然而，更严重的问题是实际的地震数据一般多少都有噪声，自身相互矛盾。图8—1—3展示了实际数据的典型的例子。在左边的例子中,小的时间变动导致点的波动，不可能把这些点与能得到期望的1—D速度模型的光滑的、物理上可靠的曲线相联系。在右边的例子中，我们注意到可以把几次地震的数据结合在一起，用这些数据确定“平均”的速度剖面.

图8—1—3 呈现离散的走时观测结果往往使速度剖面的反演复杂化

由于这两个原因，HW公式不能直接应用.下面我们对地震学家反演这些不完善的数据的一些方法进行讨论。HW公式的主要用途是阐明精确的曲线可给出速度剖面的唯一的解.因此，许多反演的策略是把寻找最佳速度模型的问题简化为寻找最佳的曲线的问题。
Error! Bookmark not defined.8.2 单台地震定位及Inglada法定位
根据走时数据确定地震的位置是地震学中的一个“古老"的，富有挑战性的问题，一直是地震学研究的一个重要方面.地震通常按发震时间和震源位置来定义.震源是地震的位置,而震中是震源正上方地表上的一个点。在地震定位中,通常把地震作为点源来对待。对于破裂几十到几百公里的大地震,震源不一定是地震的“中心”，宁可把它看成是地震发生时，地震能量开始辐射的那个点。由于破裂的速度小于波速度，可以不管地震最终的尺度和持续时间，只根据初至到时来确定震源.标准目录给出的地震资料，例如震中的初步确定（PDE）依据的是高频体波的走时，不能把这些发震时间和震源位置与长周期反演结果相混淆，后者往往给出地震的矩心时间和位置,表示整个地震的平均发震时间和位置.本节我们首先讨论单台定位法和多台定位的Inglada方法。
8.2.1单台定位法
单台定位的原理就是根据纵波初动确定出震中方位角，然后根据震相到时（走时表等）确定出震源到台站的距离,根据水平向和垂直向记录的数值，可以求得视入射角，从而根据地震波速度估计出真入射角。根据震源到台站的距离和真入射角求得地震的震中距和震源深度。注意震中方位角是指通过台站子午线与地震台站到震中连线间的夹角，沿顺时针量取为正。
由图8-2—1所示,设地震台站的P波的北向初动半周期位移为，东向初动半周期位移为,则水平向位移矢量与北向的夹角可以由下式来计算
Error! Bookmark not defined.(Error! Bookmark not defined.-2-Error! Bookmark not defined.)
在MATLAB中可以用atan2（uE，uN)来求得，因为atan2函数根据的正负号给出到的范围，将其转换为度，并且对于小于0度的角度加上360度，就转换为0～360度范围。知道了水平向矢量方向，地震震中应该位于该矢量方向上，但究竟是沿着矢量方向还是与矢量方向相反？这需要采用垂直向记录。
我们定义垂直向记录向上为正。 垂直方向的初动半周期决定地震射线的方向.当垂直向初动向上时，地震射线的方向背向震中;当垂直初动向下时，地震波射线的方向指向震中。图8-2-2为一断层错动图,台1观测到的是压缩波（+），其P波垂直向记录向上，即射线是“离源”（即背向震中）运动；此时震中位于水平矢量的相反方向上。台2观测到的是膨胀波（-），其P波垂直向初动向下，即波射线为“向源”（指向震中）的运动.此时震中位于水平矢量方向上。根据这种分析，我们可以得到结论:垂直方向初动向上的台站，震中位于水平矢量方向的相反方向上;反之垂直方向初动向上的台站，震中位于水平矢量方向上。

图 8-2-2 两个台站接收到地震P波初动符号，判断地震所在方位的示意图
第二步,确定震源到台站的距离
设ts，tp分别为该台站S波及P波的到时,则地震到达该台站的S波和P波的走时之差等于S波及P波的到时差,则有：
\＊ MERGEFORMAT （0-Error! Bookmark not defined.—Error! Bookmark not defined.)
通常定义
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.—Error! Bookmark not defined.-Error! Bookmark not defined.）
该物理量具有速度的量纲，可以根据当地的平均地壳S波速度和P波速度估计该参数。通常被称为虚波速度。设若Vp = 5.0 km/s，Vs = 3。0 km/s则虚波速度k=7.5km/s。
根据台站上读取的该地震的S波和P波到时差和该地区的虚波速度,可求得震源到台站的距离r，公式为
…………………………………………………。. （Error! Bookmark not defined.-Error! Bookmark not defined.-Error! Bookmark not defined.）
第三步：地震记录的垂直向和水平向的记录，求得视入射角，从而根据地球均匀速度模型估计真入射角。
地震在该台站的水平向位移由南北向和东西向位移记录合成：
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.-2—Error! Bookmark not defined.)
则视入射角为：
\＊ MERGEFORMAT (0—Error! Bookmark not defined.-10)
根据第四章的真入射角和视入射角的转换公式（4—2—43），即得到真入射角。
第四步,根据真入射角和震源到台站的距离r求得震中距和震源深度h。
\＊ MERGEFORMAT （Error! Bookmark not defined.—Error! Bookmark not defined.—11)
\＊ MERGEFORMAT （Error! Bookmark not defined.—2—Error! Bookmark not defined.)
震中方位角及震中距求得以后，便可决定出震中位置。这样应用一个台站的三分向记录就可以估计地震的大概位置了。
在计算机高速发展的今天，计算机可以代替人工操作。我们通常给出台站的经纬度，需要求解地震的经纬度和深度.在近震的情况下，通常采用直角坐标来解决问题，因此需要进行经纬度转换为直角坐标系的转换。
将观测点的经纬度转换为直角坐标有很多成熟的公式可以使用（如，朱介寿等,1988）。由于本研究的区域范围较小，我们采用下述简单公式（时振梁等,1990）。将观测点的经纬度换算为直角坐标的公式如下:
\＊ MERGEFORMAT Error! Bookmark not defined.（0—Error! Bookmark not defined.-Error! Bookmark not defined.）
式中为直角坐标的原点。，以度为单位，xi，yi的单位为km。其MATLAB程序如下:
function [x，y］=geo2dist（lat，long,lat0，long0）
%根据坐标系原点lat0和long0，将地理坐标(lat，long）转换为直角坐标x,y
%输入 lat，long分别纬度和经度
％ lat0和long0为坐标原点
％ 输出x和y为得到的直角坐标，x沿经向方向的坐标（东向)，y沿纬向方向（北向)
c1 = 111。199； ％1度弧长
y=c1*（lat-lat0)； ％北向坐标(8-2—9）式第一式
rlatP=deg2rad（lat+lat0）/2； ％纬度的平均值作为计算经度的大圆的半径
x=c1＊（long-long0)*cos（rlatP); ％经向坐标（8—2—9）式第二式
return

采用上面的程序得到震源位置在直角坐标系的表示后,还需将直角坐标系中的坐标转换为经纬度,公式为:
\＊ MERGEFORMAT Error! Bookmark not defined.(Error! Bookmark not defined.-Error! Bookmark not defined.—Error! Bookmark not defined.)
其中的符号及意义与上式相同。
采用MATLAB进行转换的程序如下：
function [latP,longP］=dist2geo(x，y,lat0,long0)
%根据坐标系原点lat0和long0，将直角坐标x,y转换为地理坐标
%输入 x，y分别为纬度方向和经度方向的直角坐标距离
％ lat0和long0为坐标原点
% 输出latp和longp为得到的直角坐标点的纬度和经度
c1 = 111.199; %1度弧长
latP=y/c1+lat0； ％(x,y）点的纬度，（8-2—10）式第一式
rlatP=deg2rad（latP+lat0）/2； %纬度的平均值作为计算经度的大圆的半径
longP=x/cos(rlatP)/c1+long0； %(x，y）点的经度，（8—2—10）式第二式
return
为验证上面的程序，我们以北纬25°,东经120°作为坐标原点，求北纬26°，东经121°的直角坐标，运行“［x,y]=geo2dist(26,121，25，120)",则可以得到x= 100.37km，y=111。2km的结果，然后用此结果求解其经纬度,运行“[latp，longp］=dist2geo(x，y,25,120）”的语句,可以得到北纬26°,东经121°的结果。读者也可以采用其他的数据进行检验。
有了上面的坐标转换程序，我们可以编写根据单台记录求解地震位置的计算机程序为：
Slat=35； Slon=115； %台站的经纬度
ue=3；un=3;uu=-2; ％地震台记录的南北分向、东西分向和上下分向
Parrival=5；Sarrival=8； ％P波和S波到时，单位为秒
vp=5;vs=3； %平均的P波速度和S波速度，单位为km/s
k=(vp＊vs)/(vp—vs）； %虚波速度
azi0=atan2(ue，un)； ％按（8—2-1）式求解方位角
if(uu<0）
azi=azi0；
disp（［’地震相对于台站的方位角为：’,num2str(rad2deg(azi)）］);
else
azi=2＊pi-azi0;
disp(［’地震相对于台站的方位角为：’,num2str(rad2deg(azi))］);
end
uh=sqrt（ue*ue+un＊un)； ％南北分向和东西分向在水平面的投影
ip1=atan2（uh,abs(uu)）; %求视入射角
ip=asin(vp/vs＊sin（ip1/2)）; %根据（4-2-43）给出真入射角
r=k*(Sarrival—Parrival)； %按(8—2-4）式求震源到台站的距离
h=r*cos（ip）; %求地震的深度
delta=r*sin（ip）； ％求地震的震中距
x=delta*sin（azi); %得到以台站为坐标原点的x坐标
y=delta*cos(azi）; ％得到以台站为坐标原点的y坐标
［Elat，Elon]=dist2geo(x,y，Slat,Slon);
disp(sprintf（’地震的经度东经％f度，纬度为北纬％f度，深度％f km',Elon,Elat，h)）；
disp（sprintf(’地震到台站的震中距为％fkm',delta)）
disp(sprintf('P波的走时为%f秒’,r/vp））
disp(sprintf(’S波的走时为%f秒’，r/vs)）

8。2.2三站求取震中位置的石川法
以三站为一组的作法,称为虛波速度法，亦称石川法。虛波速度在近震范围內是相当稳定的，可由本地经验得到,一般都在8公里/秒左右.设选定的三个观测站为S1，S2，S3，先将它们按照地理位置，标在地图上,如图8—2-2(a）所示,再按各站记录到的(ts—tp)到达时间差，用公式求得各站相应的震源距离：D1，D2，D3。然后以S1为中心，以D1为半径，在操作地图上作第一地震圆，其实应为一个隐藏在下的半球面.震源应当是此球面上的一点.同样以S2及S3为中心,以及为半径作第二、第三地震图，各有其下隐的半球面。这三个半球面互相交切,每两个构成一个弧形交切痕,共有三条，投影到图上为：aa， bb， cc。这三条弧相交于一点，因它是三个下隐半球面的共同点，当然就是震源,投影到图上为E,就是震中,便可在操作地图上定出其地理位置。
过震中E,垂直于观测点S1（S2或S3亦可）与震中的连线，作一直线与第一地震圆交于A，则EA之长就是震源深度h。这个从图8—2—2的右图(b）清楚的看到，b1是第一地震圆的半球面，今若过E并垂直S1E于作一垂直面切下，则这个垂面必通过震源F，并与地震圆的半球面交割成一半圆形割痕，为b2，很显然AE=EF=h，都是这半圆形的半径.最后用操作地图的比例尺，将AE之长折合成公里数，便得到震源深度的数据，如图8-2-3.
A
A
F
E
h
(2)
b
A
A
S1
c
E
c
S3
D3
S2
D2
D1
b
a
a

S1
F
A
A
E
(1)
(震源)

（a) （b）
图8—2-3 三站交切法测定震源位置示意图
8.2.3 四个台站(或以上）求解地震位置的Inglada算法
设震源坐标为,发震时刻为，台站i 的坐标为,P 波、S 波到时分别为和 ,震源到台站i 的距离为 。采用均匀介质模型,设P 波、S 波速度分别为和 ,引入虚波速度，则有
Error! Bookmark not defined.(0-Error! Bookmark not defined.—Error! Bookmark not defined.)
从而
\＊ MERGEFORMAT Error! Bookmark not defined.(0-Error! Bookmark not defined.-12）
对台站i有 　 Error! Bookmark not defined.(0—Error! Bookmark not defined.-Error! Bookmark not defined.）
对台站j有 　 \＊ MERGEFORMAT （Error! Bookmark not defined.-Error! Bookmark not defined.-Error! Bookmark not defined.)
将方程 、 展开后相减,得到线性方程　
\＊ MERGEFORMAT (Error! Bookmark not defined.—2-Error! Bookmark not defined.)
其中。将方程 用于3 对以上的台站,即得到一关于的线性方程组，方程组写成矩阵形式为:
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.—Error! Bookmark not defined.-Error! Bookmark not defined.）
可以写为的形式，从而方程的解就变成了求解A的逆矩阵的形式,方程的解可以表示为：。这样就得到了震源位置.
由于 表示台站高度,当台站海拔高度相差不大时， 远小于和 ,使得方程组的系数矩阵具有奇异性，导致解发散。可以通过其它方法单独求解:取震中解,,将其代入式求得,将各个台站求得的的平均值作为震源深度的解。这样的误差仅决定于震中的定位误差和的到时读数误差（二者均是可控的) ,从根本上解决了无法得到震源深度有效信息的问题。
关于发震时刻，有方程
Error! Bookmark not defined.(0—2—14)
将用于各个台站,然后求各台站得到的的平均值即可作为发震时刻。由于Ri 完全取决于到时( 式) ，所以和震源的求解完全独立，其误差的唯一来源是到时读数的误差,只要到时读数足够精确,即使震源定位误差很大，也可求得很精确的。根据上述方法,可以初步求解出地震发生的经度、纬度、深度和发震时间的估计解.下面为求解上述问题的计算机程序：
function ［FirLocal］=ingladawan（Vp，Vs,Pg，Sg,Stx,Sty，Stz）
%虚波速度
VSpeed=（Vp*Vs)/(Vp-Vs);
%计算地震位置
l = 1; %方程序号从1开始
for ii=1：length(Pg)—1
Ri2 = (VSpeed*（Sg(ii） - Pg（ii)))^2； %震中距的平方，（8-2-12)式
ri2 = Stx（ii）^2 + Sty(ii）^2 + Stz（ii)^2； ％公式
for jj=ii+1：length（Pg） ％ii和jj不能相同
Rj2 = (VSpeed＊(Sg（jj) - Pg（jj）））^2; %震中距的平方,（8-2—12）式
rj2 = Stx（jj)^2 + Sty（jj)^2 + Stz（jj)^2； ％采用
LMat（l，1） = (ri2 —rj2 + Rj2 — Ri2）/2; %
RMat(l，1） = Stx（ii) - Stx(jj）; ％（xi-xj)
RMat(l，2) = Sty(ii） - Sty(jj）； %(yi—yj)
RMat（l,3) = Stz（ii） — Stz（jj); ％（zi—zj）
l=l+1；
end ％ for jj 循环结束
end % for ii循环结束
FirLocal = pinv(RMat）*LMat; ％ 公式x=A-1B，pinv（x)为求解矩阵x的伪逆
%发震时刻和地震深度没有采用前面的计算结果，下面分别计算
for ii=1：length（Pg）
SeiDist =VSpeed*（Sg(ii) — Pg（ii))； %震中距
SeiTime（ii) = Pg（ii） — SeiDist/Vp； %计算发震时刻
SeiDeeps(ii) = Stz（ii） + sqrt(SeiDist^2 - (FirLocal（1） - Stx(ii）)^2 - (FirLocal（2) — Sty（ii））^2); ％根据（8-2—7）式计算地震深度
end %for ii 循环结束
FirLocal(3) = real（mean(SeiDeeps））; ％估计地震深度的平均值
if （ FirLocal(3） 〈