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

1、错误 !未定义书签。错误 !未定义书签。第八章 走时数据的反演8.9.8.1球层介质的速度反演.9.0.9.0.9.1.8.2单台地震定位及 Inglada 法定位.1.0.0.8 。1.1 拐点法（古登堡法 ） 8.1。2 Hergloz -Wiechert 方法1.0.0.1.0.6.8.38.2。 3 采用三台站求解地震位置的计算机算法地震定位结果与台站分布的关系1.0.7.1.1.0.8.2.1 单台定位法8.2.2 三站求取震中位置1.1.4.1.1.4.1.1.5.1.3.2.8 。 4 迭代定位方法8.4.1 迭代定位算法的基本思路8。4。 2 迭代方法的具体实现8 。 5 相对

2、定位方法第八章 走时数据的反演在前两章我们根据已知的速度结构对射线追踪和走时曲线的计算问题进行 探讨，推导了波速只随深度变化的一维速度模型的射线追踪的表达式。 三维结构 的射线追踪虽然较复杂， 但它遵循的原则是类似的。 现在我们来研究根据观测到 时数据得到速度结构和地震位置的问题。前面两章是已知速度结构和地震位置， 如何计算地震波走时的问题， 而本章是根据地震波走时得到地震结构和震源位置 的问题 ,是前面两章的反问题 .在这个问题的研究中，地震学家通常把问题进行简 化：第一种简化是对于震源位置已经清楚的地震， 根据地震波走时数据求解速度 结构，第二种简化是已经得出了该地区的速度结构， 根据地震

3、波到时求解地震位 置 .限于篇幅，对于第一种简化 ,本书仅讲授球层介质一维速度结构的反演 .对于第 二种简化，我们讲授地震定位的基本方法。8。1球层介质的速度反演8。1.1拐点法（古登堡法）这个方法是利用走时曲线的拐点处与从震源处水平方向射出的射线相对应，据此可 求出震源处的波速.根据第七章球层介质中的 Snell定律有:rhSi nihrsinipvVh（错误!未定义书签。-错误!未定义书签。一1）其中h代表震源处的参数，rhR h,vh vrh都是常数从震源向不同方向射出的射线，i h值不同，相应的射线参数p也不同，当ih=90时，rhS!nh达到极大值，即射线参数也达到极大值Pmaxrh

4、Vh，因此有dpih2错误!未定义书签。（错误!未定义书签。-1-1 ）根据球层介质中的Bendoff定律（7-5-3）式,pd2t伙MERGEFORMAT错误!未定义书签。（错误!未定义书签。扌签。d2t曲线的拐点，并确定该点的走时曲线的斜率dtd，则M0对应于走时曲线的拐点因此走时曲线的拐点与从震源处水平射出的射线相对应.这样，我们只要求得某地震的震源深度及有观察分析归纳得到其走时曲线，找出走时SSI nihVhd£ d m错误!未定义书签(错误!未定义书签。一1错误!未定义书签。)而对应的M点有ih ，因此有Vhdtd mdtRd伙 MERGEFORMAT扌签。Vm错误!未定义

5、书签。一错误!未定义书签。)其中，h为震源深度，VM为走时曲线拐点 M点的视速度(参看第七章球层介质的Benndorf定律)，R为地球半径此方法原理清楚、方法简单、计算方便。但由于地震发生的深度大约在0 700km的范围内，采用这种方法只能求出0-700km处的波速。由于拐点不易找准，得到的精度较差,它又要求对应每个深度的地震就要总结出一条走时曲线，因而资料分析工作量很 大。8.1.2 Hergloz -Wiechert 方法考虑球面介质，速度随半径r变化v( r ),根据7.2节中得出震中距的表达式有R2pr1rdr-错误!未定义书签伙MERGEFORMAT (0 错误！未定义书签rr si

6、n i其中 二，r1为射线转折点距地心的距离，p= , R为地球半径。改变积分变数vv(r)则c 01dr ,2pd2 2 d p r ip d(0-错误!未定义书签。4)其中在此介绍一积分式,错误!未定义书签其中i匕p,，为转折半径为*的射线参数，所以Vi伙MERGEFORMAT （错误!未定义书签。-1-4）式表示从转折半径为 0至到转折半径为r,的所有射线的积分 将 伙MERGEFORMAT （错误!未定义书签。-1-4）式作用于震中距表 达式 伙MERGEFORMAT （0-错误!未定义书签。一4）式的两边，得到dp.P2 i20dp12p021dr d伙 MERGEFORMAT签。5

7、)方程的左边变为0 dpd卫12"，考虑至U cosh 1 xX2X'-，可以用分1部积分的形式给出cosh 1() |p1cosh1 dp1 p()dp1伙 MERGEFORMAT扌签。扌签。（错误!未定义书签。-1-错误!未定义书签。）式的第一项代入积分上下限，由于 cosh ( 1)=0 ,且在p= 0，半径为R,A =0，因此第一项消失仅剩下第二项，可写成:1cosh0MERGEFORMAT错误!未定义书签。（错误!未定义书签。一错误!未定义书签。一6）（错误!未定义书签。-错误！未定义书签。一5 ）的右侧为一双重积分式，改变积分变数 得到：R dr pr1 r pp

8、dpR drr1 r1、p 1 . pp2 pdpp 1 . p212 . 2 IR drr1dp2221R dr .arcs inr1 r2 2 2、2p (1 )2 21R dr(arcs in 1 rarcsinR drR drririln r |RInri显示文本”不能横跨多行！这里采用了不定积分公式:du2cu b整理变为bu2cu1 .arcs in - cb2 4acR drr1 r伙MERGEFORMAT错误!未定义书签（错误！未定义书签。一1-6）与伙MERGEFORMAT （错误！未定义书签错误!未定义书签。一6）式结合，得出ln（R） 01 cosh 1 ）dA1伙MER

9、GEFORMAT错误!未定义书签。（错误！未定义书签。一错误!未定义书签。一6）令3 o'cosh1（卫）d ，可以得到:1r1 R exp（s ）衣MERGEFORMA错误!未定义书签。-1-错误！未定义书签。）这就是射线最深点的半径的表达式。此表达式可以根据观测走时曲线 t （）得出。如果t（）为一平滑曲线，我们可由其切线斜率得到p（）（由于p=也）。dd伙MERGEFORMAT （错误!未定义书签。一错误!未定义书签。一6）式中的-JX1Pi()1,为在震中距1的走时曲线斜率.S1的积分式由震中距A从 0到A 1d的p值，得到S1后，根据伙 MERGEFORMAT签。)式求出A，

10、即其最深转折点。根据转折点为r1的1 土，得出其对应此深度的速度V1V1j1伙MERGEFORMAT错误!未定义书签。(错误味定义书签。-1-7 )依此类推，得出速度与深度的分布图F面归纳一下求速度分布的具体步骤1、做出0-范围内的走时曲线，即t-曲线，震源要取在地表上，否则 要进行修正；在该步骤中，通常在计算机上利用三次自然样条函数 对一元函数进行成组插值及微分，从而使得走时曲线具有较好的平 滑性能，采用该方法能保证所插值的函数及其一阶导数、二阶导数 连续。2、由t-曲线求曲线上各点的斜率，从而做出 空 曲线；该步骤求出pd随震中距的分布。3、在0-范围内任取1，取出1点对应的dtd 1 ?

11、这就是1dTP1 ()1 °d4、将dt曲线中d0-1段除以常数理d即11，得到，即卫曲线，从而求出cosh1 p 1 曲线;这就是cosh 1(-p)随震中距的变化。15、在0 1范围内求出cosh 1 p 1 曲线下的面积,即求cosh 1 兰 兰 d ，这就是1 cosh 1 （ ）d ;在计算积分的过 0dd i0i程中，通常采用Simpson积分数值公式进行。6、根据 伙MERGEFORMAT （错误 未定义书签。-1-错误!未定义书签。）式求 出r1，然后根据 衣MERGEFORMAT （错误!未定义书签。-1-7 ）式得到r1所 对应的速度v。7、取不同的1值，重复上述

12、步骤，从而求出不同的R值所对应的速度值 w ;以一个具体地区的走时曲线为例子：采用的数据是北京台网、大连、长春、牡丹江等26个台站的记录，并选用从北京地区到阿留申群岛西端的300多个地震，经过震源深度及剥壳校正后 （地球平均厚度35公 里）。得到的北京一萨哈林一线的P波走时曲线。数据如表 8 1 1为北京-萨哈林剖面P波走时观测值。表8-1-1北京-萨哈林剖面P波走时观测值? (°t(s)? (°t (s)0。0020.0265。21.014。421.02762。028.722。0286。83.042.823。0297。64。056.824。0308。35。070。825.

13、0317。56.084.826。0326.57.099。027.0335.58。0113.428。0344。59。0127。829。0353.310。0141.830。0362。111。0154。832。0379.712.0167。834。0397.313.0180.836。0414.114。0193.838.0430.915.0206.840。0447.716。0219。842。0464.217。0232.844.0479.518.0243。646.0494.519.0254.448。0509。5然后将走时曲线用Herglotz-Wiechert方法利用MATLAB中的三次样条插值和辛普森积

14、分方法，编制成 MATLAB程序。程序如下： clcclfclear allload bj_shl_zs.txt % 调用数据 x0=bj_shl_zs(:,1)'%震中距 yO=bj_shl_zs(:, 2) ' %走时 dh=0。 2;x=0.2:dh:47.8;n=length(x )；y,dy1, dy2 =splineinsert (xO,yO, x) ;%对数值进行样条插值，得到插值点的y值，一阶导数dy1,二阶导数dy2figure ( 1)%实际观测曲线plot ( x,y, ' r* 'grid ontitle ('实际观测曲线)xla

15、bel('震中距)；ylabel('走时');R=637135；In dx=1 ；r1=zeros(1，(length (3: 2: n) )；%射线能够达到的最低点的半径v=zeros(1，(length(3:2:n ) ;%射线最低点的速度for ii=3 : 2:nksi仁dy1 ( ii);%视为常数f=acosh(dy1(1:ii ) /ksi1 )；Int=dh/3 *( f+f(ii)+2 衣 sum(f(3 : 2： ii) +4*sum (f(2:2 : ii) )*pi/180 ；%采用 Simpson 公式进行积分 :f(x)dx f(a) 4f(

16、) f(b)%采用公式求解射线最深点的半径r1 (Indx)=R/exp (Int/pi)；r1 R exp(3)v(lndx ) =deg2rad (r1(Indx)/ksi1 )；%采用球层Snell定律得到最深点的速度；deg2rad()函数的作用是：将角度转换为弧度，利用公式v -1%序号+1lndx=lndx+1 ；%pause endplot(v , r0 r1 , '.' )%绘制速度和深度剖面图set(gca,'ydir'， reverse' ) %将y轴的方向进行翻转利用上述程序调用表81-1中的北京一萨哈林地区的P波走时曲线的数据。对

17、 观测曲线进行上述程序来反演即可得到速度剖面如图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在上述MATLA函数后加入如：x= 20,45, 85,125,145,165,265,365,405 445, 495, 565,705,875,965,1025,1165；y, dy1, dy2 =splineinsert (r

18、0 r1, v,x);r0-r1' v'x ' y'并利用得到的数据进行绘图得到图8 1-2定点插值的反演数据图图81-2定点插值的反演数据图当有低速区时，这些公式是有问题的。在这种情况下，（P）不连续，没有唯一的解。对这种情况，传统的做法是设想在低速区存在若干有不同速度的均匀 层。由于没有射线在这些层里转折，因此低速区的这些层可以任意移来移去，穿 过低速区的射线的积分的走时和震中距不会改变。不论解析多么完全，在地震学中很少用赫格劳兹一维歇尔 （HW）公式，这至 少有两个原因：首先，HW假定我们有已知的连续的t（）曲线，而实际上我们总是只有有限的走时点。这意味着

19、走时曲线必须在这些数据之间进行内插，这些内插方案的不同导致于有不同的速度剖面实际上，会有无数个稍有不同的速度模型 它们都适合于有限数目的t()点.然而，更严重的问题是实际的地震数据一般多 少都有噪声，自身相互矛盾。图8 1 3展示了实际数据的典型的例子。在左 边的例子中，小的时间变动导致t()点的波动，不可能把这些点与能得到期望的1 D速度模型的光滑的、物理上可靠的t()曲线相联系。在右边的例子中，我们 注意到可以把几次地震的数据结合在一起，用这些数据确定“平均”的速度剖面.图8 1 3呈现离散的走时观测结果往往使速度剖面的反演复杂化由于这两个原因，HW公式不能直接应用.下面我们对地震学家反演

20、这些不完 善的数据的一些方法进行讨论。HW公式的主要用途是阐明精确的t(X)曲线可给 出速度剖面的唯一的解.因此，许多反演的策略是把寻找最佳速度模型的问题简 化为寻找最佳的t()曲线的问题。8.2错误!未定义书签。单台地震定位及INGLADA法定位根据走时数据确定地震的位置是地震学中的一个“古老"的，富有挑战性的问题，一直是地震学研究的一个重要方面.地震通常按发震时间和震源位置来定 义.震源是地震的位置(x,y,z),而震中是震源正上方地表上的一个点(x, y)。在地 震定位中，通常把地震作为点源来对待。对于破裂几十到几百公里的大地震，震源 不一定是地震的“中心”，宁可把它看成是地震

21、发生时，地震能量开始辐射的那 个点。由于破裂的速度小于S波速度，可以不管地震最终的尺度和持续时间，只 根据初至到时来确定震源.标准目录给出的地震资料，例如震中的初步确定(PDE) 依据的是高频体波的走时，不能把这些发震时间和震源位置与长周期反演结果相 混淆，后者往往给出地震的矩心时间和位置，表示整个地震的平均发震时间和位 置本节我们首先讨论单台定位法和多台定位的 Inglada方法。821单台定位法单台定位的原理就是根据纵波初动确定出震中方位角，然后根据震相到时（走时 表等）确定出震源到台站的距离，根据水平向和垂直向记录的数值，可以求得视入射角， 从而根据地震波速度估计出真入射角。根据震源到台

22、站的距离和真入射角求得地震的震 中距和震源深度。注意震中方位角是指通过台站子午线与地震台站到震中连线间的夹 角，沿顺时针量取为正。由图8-2 1所示,设地震台站的P波的北向初动半周期位移为 uN，东向初动半周期位移为uE，则水平向位移矢量与北向的夹角可以由下式来计算uEarctanUn错误!未定义书签。（错误！未定义书签。-2-错误!未定义书签。）在MATLAB中可以用atan2（Ue,un）来求得，因为atan2函数根据的正负号给出到 的范围，将其转换为度，并且对于小于0度的角度加上360度，就转换为0360度范围。知道了水平向矢量方向，地震震中应该位于该矢量方向上，但究竟是沿着矢量 方向还

23、是与矢量方向相反？这需要采用垂直向记录。我们定义垂直向记录山向上为正。垂直方向的初动半周期决定地震射线的方向当垂直向初动向上时，地震射线的方向背向震中；当垂直初动向下时，地震波射线的方向指向震中。图8-2-2为一断层错动图，台1观测到的是压缩波（+），其P波垂直向记 录向上，即射线是“离源”（即背向震中）运动；此时震中位于水平矢量的相反方向上。 台2观测到的是膨胀波（-），其P波垂直向初动向下，即波射线为“向源”（指向震 中）的运动.此时震中位于水平矢量方向上。根据这种分析，我们可以得到结论：垂直方向初动向上的台站，震中位于水平矢量方向的相反方向上；反之垂直方向初动向上的台站，震中位于水平矢量

24、方向上。图8-2-2两个台站接收到地震 P波初动符号，判断地震所在方位的示意图第二步，确定震源到台站的距离VpVsr-设ts, tp分别为该台站S波及P波的到时，则地震到达该台站的 S波和P波的走时 之差等于S波及P波的到时差，则有：, r r 11 ts tprVsVpVs Vp伙MERGEFORMAT （0-错误!未定义书签错误!通常定义1,11VpVskVs VpVp Vs显示文本”不能横跨多行！该物理量具有速度的量纲,通常被称为虚波速度。设若可以根据当地的平均地壳 S波速度和P波速度估计该参数Vp = 5.0 km/s , Vs = 3。0 km/s 则虚波速度 k=7.5km/s。根

25、据台站上读取的该地震的S波和P波到时差和该地区的虚波速度，可求得震源到台站的距离r，公式为r k ts tp（错误!未定义书签-错误!未定义书签。-错误!未定义书签。）第三步：地震记录的垂直向和水平向的记录，求得视入射角，从而根据地球均匀速 度模型估计真入射角。地震在该台站的水平向位移由南北向和东西向位移记录合成:伙MERGEFORMAT错误!未定义书签。（错误!未定义书签。-2 错误！未定义书签。）则视入射角为：rparctan 比Uu伙MERGEFORMAT （0错误!未定义书签。-10）根据第四章的真入射角和视入射角的转换公式（4243）,即sin ip sin得p 2到真入射角。第四步

26、根据真入射角和震源到台站的距离r求得震中距 和震源深度h。r si nip伙MERGEFORMAT （错误!未定义书签错误!11)h r cosi p伙MERGEFORMAT （错误!未定义书签签。震中方位角及震中距求得以后，便可决定出震中位置。这样应用一个台站的三分向记录就可以估计地震的大概位置了。在计算机高速发展的今天，计算机可以代替人工操作。我们通常给出台站的经纬度， 需要求解地震的经纬度和深度在近震的情况下，通常采用直角坐标来解决问题，因此需 要进行经纬度转换为直角坐标系的转换。将观测点的经纬度i, i转换为直角坐标 x , yi有很多成熟的公式可以使用(如，朱介寿等，1988 )。由

27、于本研究的区域范围较小，我们采用下述简单公式(时振 梁等,1990 )。将观测点的经纬度i, i换算为直角坐标 xi, yi的公式如下：y 111.199x 111.199i00 cos -2伙MERGEFORMAT错误!未定义书签(0错误!未定义书签式中0, 0为直角坐标的原点。0, 0, i, i以度为单位，Xi,yi的单位为km。其MATLAB 程序如下：functionx , y =geo2dist (lat , Iong,latO, longO )%艮据坐标系原点lat0和long0，将地理坐标(lat , long )转换为直角坐标x,y%输入lat , long分别纬度和经度%

28、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采用上面的程序得到震源位置在直角坐标系的表示后，还需将直角坐标系中的坐标Xi,yi转换为经纬度i, i ,公式为:ii111.1990X0111.199 cos2显示文本”不能

29、横跨多行！其中的符号及意义与上式相同。采用MATLAB进行转换的程序如下:functionlatP,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 ；%纬度的平均值作为计算经度的大圆

30、的半径IongP=x/cos(rlatP)/c1+longO；%(x, y)点的经度，(8210 )式第二式return为验证上面的程序， 我们以北纬 25°,东经 120° 作为坐标原点， 求北纬 26° ，东经 121°的直角坐标， 运行“ x,y=geo2dist(26,121 ，25，120)", 则可以得到 x= 100.37km，y=111。2km的结果，然后用此结果求解其经纬度，运行"latp , Iongp =dist2geo(x ,y,25,120 )”的语句 , 可以得到北纬 26°,东经 121

31、6;的结果。读者也可以采用其他的数据进行检验。有了上面的坐标转换程序，我们可以编写根据单台记录求解地震位置的计算机程序为：Slat=35 ； Slon=115 ；%台站的经纬度ue=3； un=3;uu=-2; 地震台记录的南北分向、东西分向和上下分向Parrival=5 ； Sarrival=8 ；% P波和S波到时，单位为秒vp=5;vs=3 ；%平均的P波速度和S波速度，单位为km/sk=(vp vs)/(vp vs)；%虚波速度azi0=atan2(ue ， un) ；%按( 82-1 )式求解方位角if (uu<0 )azi=azi0 ；disp ('地震相对于台站的方

32、位角为：' ,num2str(rad2deg(azi) ) );elseazi=2 pi-azi0;disp( '地震相对于台站的方位角为：' ,num2str(rad2deg(azi);enduh=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 ) ;%求地震的

33、深度delta=r*sin ( ip )；%求地震的震中距x=delta*s in (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 波的走时为 1秒'，r/vp )disp(s

34、printf(' S波的走时为 %f秒' , r/vs)8。2.2三站求取震中位置的石川法以三站为一组的作法，称为虛波速度法，亦称石川法。虛波速度在近震范围內是相当稳定的，可由本地经验得到，一般都在8公里/秒左右设选定的三个观测站为 Si，S2，S3,先将它们按照地理位置，标在地图上 ，如图82-2(a )所示，再按各站 记录到的(ts tp)到达时间差，用公式(错误！未定义书签。-错误!未定义书签。-错误!未定义书签。)求得各站相应的震 源距离：Di, D2，D3。然后以Si为中心，以Di为半径，在操作地图上作第一地震 圆，其实应为一个隐藏在下的半球面震源应当是此球面上的一点

35、同样以S2及S3为中心，以及为半径作第二、第三地震图，各有其下隐的半球面。这三个半球面互 相交切，每两个构成一个弧形交切痕，共有三条，投影到图上为：aa， bb， cc。这 三条弧相交于一点，因它是三个下隐半球面的共同点，当然就是震源，投影到图上为E,就是震中，便可在操作地图上定出其地理位置。过震中E,垂直于观测点Si (S2或S3亦可)与震中的连线，作一直线与第一地 震圆交于A，则EA之长就是震源深度h。这个从图822的右图(b )清楚的看 至V，bi是第一地震圆的半球面，今若过 E并垂直SiE于作一垂直面切下，则这个垂 面必通过震源F,并与地震圆的半球面交割成一半圆形割痕，为b2,很显然A

36、E=EF=h，都是这半圆形的半径最后用操作地图的比例尺，将AE之长折合成公里数，便得到震源深度的数据，如图8-2-3.(a)(b)图8 2-3三站交切法测定震源位置示意图四个台站（或以上）求解地震位置的INGLADA算法的坐标为（x,yi,z）,P波、S采用均匀介质模型，设P波、设震源坐标为（X0,y°,z0）,发震时刻为t°，台站i波到时分别为tpi和tsi,震源到台站i的距离为RS波速度分别为Vp和Vs ,引入虚波速度VsVp1VpVSVp Vs，则有t R R Rtpi错误!未定义书签。（0-错误!未定义书签错误!未定义书签。）从而V tSitPiMERGEFORMA

37、T错误!未定义书签。（0-错误!未定义书签。-12）对台站i有2 2 2 2(xo Xi)(yo y)(z。zj R错误!未定义书签。（0 错误!未定义书签。-错误!未定义书签。）对台站j有2 2 2 2(Xo Xj)(yo yj)(zo Zj)Rj衣MERGEFORMAT （错误!未定义书签。-错误!未定义书签。-错误!未定义书签。）将方程 衣MERGEFORMAT1错误！未定义书签。-错误！未定义书签。）、衣MERGEFORMA错误！未定义书签。-错误！未定义书签。-错误！未定义书签。） 展开后相减，得到线性方程R21 2 2 2 (Xi Xj)Xo (yi yjy° (w zj

38、z°r RjMERGEFORMAT (错误味定义书签。一2-错误!未定义书签其中ri2x22 2yiz。将方程（错误!未定义书签。一2-错误!未定义书签。）用于3对以上的台站，即得到一关于（xo,yo,zo）的线性方程组，方程组写成矩阵形式为：XiX2yiy2zi Z2XoXiX3yiy3ziZ3yo*4*Zo2 ri222 riRi2Ri2显示文本”不能横跨多行衣MERGEFORMAT （错误!未定义书签。错误!未定义书签。-错误!未定义书签。）可以写为Ax B的形式，从而方程的解就变成了求解 A的逆矩阵的形式，方程的解可以表示为：x A iB。这样就得到了震源位置由于Zi表示台站

39、高度，当台站海拔高度相差不大时，乙Zj远小于Xi Xj和y yj ，使得方程组* MERGEFORMAT （错误!未定义书签。一错误!未定义书签。-错误!未定义书签。） 的系数矩阵具有奇异性，导致Zo解发散。可以通过其它方法单独求解 Zo：取震中解&, y°，将其代入* MERGEFORMAT（0-错误!未定义书签。-错误!未定义书签。）式求得 将各个台站求得的Zo的平均值作为震源深度的解。这样Zo的误差仅决定于震中的定位误差和R的到时读数误差（二者均是可控的），从根本上解决了无法得到震源深度有效信息的问题。关于发震时刻，有方程t0tpi旦Vp错误!未定义书签。（0 2 14

40、）将* MERGEFORMAT214）用于各个台站，然后求各台站得到的t。的平均值 即可作为发震时刻。由于R完全取决于到时（0-错误！未定义书签。-i2 ）式）， 所以to和震源的求解完全独立，其误差的唯一来源是到时读数的误差 ，只要到时读数足够精确，即使震源定位误差很大，也可求得很精确的to。根据上述方法,可以初步求解出地震发生的经度、纬度、深度和发震时间的估计解下面为求解上述问题的计算机程序:function FirLocal =ingladawan(Vp, Vs,Pg , Sg,Stx,Sty , Stz )%虚波速度VSpeed= (Vp*Vs)/(Vp-Vs);%计算地震位置l =

41、1; %方程序号从1开始for ii=1 ： length(Pg) 1Ri2 = (VSpeed*(Sg(ii ) - Pg (ii)A2；% 震中距的平方，(8-2-12)式ri2 = Stx(ii )卜2 + Sty(ii ) Q + Stz(叩2 ；%公式 x?yi2Z?for jj=ii+1: length ( Pg)% ii 和jj 不能相同Rj2 = (VSpeed* (Sg (jj) - Pg(jj ) A2； % 震中距的平方，(8-2 12)式rj2 = Stx(jj)A2 + Sty (jj)A2 + Stz(jj)A2 ；%采用 ri2xi2y2 zi22 222ri r

42、j Rj R LMat(l , 1) = (ri2 rj2 + Rj2 Ri2 ) /2; %- jj-2RMat(l,1) = Stx (ii) - Stx(jj)；%(xi -x j)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

43、ii=1length (Pg)SeiDist =VSpeed*(Sg(ii) Pg (ii) ； % 震中距SeiTime（ii） = Pg（ii ） SeiDist/Vp ；%计算发震时刻SeiDeeps(ii)=Stz (ii ) + sqrt(SeiDistA2- (FirLocal (1) - Stx(ii ) )A2 - (FirLocal(2) Sty ( ii）A2）;%根据（8-2 7）式计算地震深度end %for ii循环结束FirLocal（3） = real（ mean（SeiDeeps ） ; 估计地震深度的平均值if （ FirLocal（3 ） 0 | FirLo

44、cal（ 3 ） > 100）%如果地震深度超出可信范围，则强制给定地震深度FirLocal(3) = 12;endFirLocal（ 4） = mean （ SeiTime） ； %估计地震发震时刻的平均值return一个例子Stx=50 0 50 100 100 100 50 0 0;%台站位置的 x坐标Sty= : 50 100 100 100 50 0 0 0 50；%台站位置的 y坐标Stz=zeros（1,9）; %台站位置的z坐标，在此设置为零Pg= : 6。4 18.5 11。9 11.9 6 。4 11。9 11。9 18.5 15 。5 ; % 各个台站的直达 P波到

45、时Sg=10.7 30。8 19.8 19 。8 10.7 19.8 19.8 30.8 25。9 ; %各个台站的直达 S波到时FirLocal =ingladawan （5， 3,Pg ， Sg,Stx,Sty,Stz ） % 调用和达法求解地震位置在上面的程序中， 我们均采用直角坐标 . 如果已知的台站坐标的经纬度， 求解地震的经纬度， 则需要采用（ 8-29）式首先转换为直角坐标,待求得地震位置后再采用（82-10 ）式转换为地震的经纬度地理坐标 .错误 !未定义书签。 8.3 地震定位结果与台站分布的关系 地震定位问题需要求解地震的空间位置和发震时刻 .我们把这些参数称为模型， 定义

46、模型矢量：m (g,m2,m3,m4)(T,x, y,z)衣MERGEFORMAT (错误!未定义书签。一错误!未定义书签。一15)现在假定我们有地震台观测的n个震相观测到时.为了得到参数m，我们首先必 须假定一个参考的地球模型，如一维平层介质模型。这样，按照第6 7章的知识, 根据m值和台站位置可以计算第i个台站的预测到时：tip Fi (m)错误!未定义书签。(错误!未定义书签。一错误!未定义书签。一15)这里F是算子，其复杂程度由假定地球模型的复杂程度和地震波传播的射线路径决定。观测到时与预测的理论到时之差为：ri ti 丁 ti Fi (m)衣MERGEFORMAT (错误!未定义书签

47、。-错误!未定义书签。一错误!未定义书签。)这里ri是第i个台站的残差.在某种意义上来说，我们希望得到使观测到时与预测 到时的残差为最小的m.注意，F是地球模型和每个台站位置的函数，根据第 6, 7章的知识就可以精确计算。然而，走时与地震位置参数的非线关系使反演最佳 地震模型的工作大大复杂化。这种非线性的影响基至在均匀速度进行定位的简单 情况下，也是明显的，在这种情况下，座标为(xi, yi)的台站记录震源点(x, y)发震 时刻为T的到时为：ti(x x)2 (y yj2衣MERGEFORMAT错误!未定义书签。(错误!未定义书签。一错误!未定义书签。-16)这里 是速度。在此方程中，到时t

48、不论是与x还是y都不成线性关系.该结果表 明，我们不能用解线性方程的方法来得到震源位置。现在有了功能强大的计算机，我们可在所有可能的地震位置和发震时间范围 里作网格搜索，计算每个台站预测的到时。那么，我们可以求得使预测的时间tip 与观测的时间ti最一致的特定的m。怎样规定“最”一致呢？流行的选择是最小 二乘法，即求使下式达最小的 m值：tii 1ti显示文本”不能横跨多行这里n是台站数目。把平均的平方残差力叫做方差。在数理统计上更普遍的形式是用ndf来定义方差这里ndf是自由度的数目（ndf为n减去拟合中的独立参数的数目）对通常需要解决的问题，拟合参数的数目比数据的数目少得多 所以n和ndf

49、近似相等。在数理统计中2分布于自由度个数及置信度的关系如表8-3 1所示。表8-3-12分布与自由度和置信度的关系ndf2(95 %)22 (50%2(5%51。154.3511。07103。949。3418。312010。8519。3431.415034。7649.3367.5010077。9399。33124。 34最小二乘法常被用来求解这类问题， 这是因为在求最小值的问题中，它使方 程有简单的解析形式。如果t与tp之间的拟合差是由t中非相关的随机噪声所造 成的，那么最小二乘法将有助于给出正确的答案。在本门课程的讨论中，如不特殊说明，我们假定误差是由服从高斯正态分布的随机噪声导致的。现在我

50、们以一组台站，假设为均匀的地壳模型，以1km和0.1s为网格点搜索地震的位置。MATLAB程序如下：clcx0=30 ； y0=29;%假定的震源位置，发震时刻假定为零S1%x0=50 ; y0=20; %假定的震源位置，发震时刻假定为零S2% x0=10;y0=60; % 假定的震源位置，发震时刻假定为零S3vel=6 ;%假定的地震波速度St=9.0,24.0;24。0，13。2; 33。0,4.8;45.0，10。8; 39。0,27。0;54。0，30。0;15.0 ，39。 0；36。0，42。0;27.0,48 。0；48.0,48 。0；15。0，42.0 ；18。0，15.0

51、；33。 0,36.0 ； %假定的台站位置N=length(st )； 台站个数Pt=sqrt ( st(:,1 ) -x0 )。A2+ (st(: , 2) y0)。A2)./vel+ ( 1-rand(N , 1) *0.5;%计算各台站的预测走时，并加上幅度为 0.5s 的误差ndf=N3；%该问题需要求解地震的发震时刻和 x,y 坐标，因此有三个未知数resmin=1 。 0E50;%给出最小值是一个较大的数xmin=0;ymin=0 ；origt=min(Pt )；%初始的位置for x=0 ： 1 ： 100for y=0:1 ：100Pred_t=sqrt(xones(N， 1

52、)-st (：,1) 。 A2+(y*ones (N， 1)st( ：， 2).A2 )/vel ；%按 (8-3 4)式求得预测到时for orig=origt-10:0 。 1:origtres=sum(Pt-Pred_torig).A2 )/ndf ； %按照 (8-3 5 )式计算平均残差if (res<resmin)resmin=res; xmin=x;ymin=y;origmin=orig；endendendendxmin ， ymin ， origmin ， resmin%输出找到的最优结果orig=origmin;%将最优的发震时刻保存figure ( 1)for x=1:1 ：101for y=1 ：1 ： 101Pred_t=sqrt( xst (：,1) ).A2+(y-st( ：,2 )。 A2)/vel ；%按( 8-3 4)式求得预测到时chisq(x,y)=sum(Pt Pred_t orig) 。 A2); %按照 (8-3 5)式计算 Chisqendendchisqlim=interp1(5,10,20,50,100, : 1。15, 3.94 , 10。85,34.76,7