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

1、錯誤 !未定義書簽。錯誤 !未定義書簽。第八章 走時數(shù)據(jù)的反演8.9.8.1球?qū)咏橘|(zhì)的速度反演.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 采用三臺站求解地震位置的計算機算法地震定位結(jié)果與臺站分布的關系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 迭代方法的具體實現(xiàn)8 。 5 相對

2、定位方法第八章 走時數(shù)據(jù)的反演在前兩章我們根據(jù)已知的速度結(jié)構(gòu)對射線追蹤和走時曲線的計算問題進行 探討，推導了波速只隨深度變化的一維速度模型的射線追蹤的表達式。 三維結(jié)構(gòu) 的射線追蹤雖然較復雜， 但它遵循的原則是類似的。 現(xiàn)在我們來研究根據(jù)觀測到 時數(shù)據(jù)得到速度結(jié)構(gòu)和地震位置的問題。前面兩章是已知速度結(jié)構(gòu)和地震位置， 如何計算地震波走時的問題， 而本章是根據(jù)地震波走時得到地震結(jié)構(gòu)和震源位置 的問題 ,是前面兩章的反問題 .在這個問題的研究中，地震學家通常把問題進行簡 化：第一種簡化是對于震源位置已經(jīng)清楚的地震， 根據(jù)地震波走時數(shù)據(jù)求解速度 結(jié)構(gòu)，第二種簡化是已經(jīng)得出了該地區(qū)的速度結(jié)構(gòu)， 根據(jù)地震

3、波到時求解地震位 置 .限于篇幅，對于第一種簡化 ,本書僅講授球?qū)咏橘|(zhì)一維速度結(jié)構(gòu)的反演 .對于第 二種簡化，我們講授地震定位的基本方法。8。1球?qū)咏橘|(zhì)的速度反演8。1.1拐點法（古登堡法）這個方法是利用走時曲線的拐點處與從震源處水平方向射出的射線相對應，據(jù)此可 求出震源處的波速.根據(jù)第七章球?qū)咏橘|(zhì)中的 Snell定律有:rhSi nihrsinipvVh（錯誤!未定義書簽。-錯誤!未定義書簽。一1）其中h代表震源處的參數(shù)，rhR h,vh vrh都是常數(shù)從震源向不同方向射出的射線，i h值不同，相應的射線參數(shù)p也不同，當ih=90時，rhS!nh達到極大值，即射線參數(shù)也達到極大值Pmaxrh

4、Vh，因此有dpih2錯誤!未定義書簽。（錯誤!未定義書簽。-1-1 ）根據(jù)球?qū)咏橘|(zhì)中的Bendoff定律（7-5-3）式,pd2t伙MERGEFORMAT錯誤!未定義書簽。（錯誤!未定義書簽。扌簽。d2t曲線的拐點，并確定該點的走時曲線的斜率dtd，則M0對應于走時曲線的拐點因此走時曲線的拐點與從震源處水平射出的射線相對應.這樣，我們只要求得某地震的震源深度及有觀察分析歸納得到其走時曲線，找出走時SSI nihVhd£ d m錯誤!未定義書簽(錯誤!未定義書簽。一1錯誤!未定義書簽。)而對應的M點有ih ，因此有Vhdtd mdtRd伙 MERGEFORMAT扌簽。Vm錯誤!未定義

5、書簽。一錯誤!未定義書簽。)其中，h為震源深度，VM為走時曲線拐點 M點的視速度(參看第七章球?qū)咏橘|(zhì)的Benndorf定律)，R為地球半徑此方法原理清楚、方法簡單、計算方便。但由于地震發(fā)生的深度大約在0 700km的范圍內(nèi)，采用這種方法只能求出0-700km處的波速。由于拐點不易找準，得到的精度較差,它又要求對應每個深度的地震就要總結(jié)出一條走時曲線，因而資料分析工作量很 大。8.1.2 Hergloz -Wiechert 方法考慮球面介質(zhì)，速度隨半徑r變化v( r ),根據(jù)7.2節(jié)中得出震中距的表達式有R2pr1rdr-錯誤!未定義書簽伙MERGEFORMAT (0 錯誤！未定義書簽rr si

6、n i其中 二，r1為射線轉(zhuǎn)折點距地心的距離，p= , R為地球半徑。改變積分變數(shù)vv(r)則c 01dr ,2pd2 2 d p r ip d(0-錯誤!未定義書簽。4)其中在此介紹一積分式,錯誤!未定義書簽其中i匕p,，為轉(zhuǎn)折半徑為*的射線參數(shù)，所以Vi伙MERGEFORMAT （錯誤!未定義書簽。-1-4）式表示從轉(zhuǎn)折半徑為 0至到轉(zhuǎn)折半徑為r,的所有射線的積分 將 伙MERGEFORMAT （錯誤!未定義書簽。-1-4）式作用于震中距表 達式 伙MERGEFORMAT （0-錯誤!未定義書簽。一4）式的兩邊，得到dp.P2 i20dp12p021dr d伙 MERGEFORMAT簽。5

7、)方程的左邊變?yōu)? dpd衛(wèi)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 ）的右側(cè)為一雙重積分式，改變積分變數(shù) 得到：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整理變?yōu)閎u2cu1 .arcs in - cb2 4acR drr1 r伙MERGEFORMAT錯誤!未定義書簽（錯誤！未定義書簽。一1-6）與伙MERGEFORMAT （錯誤！未定義書簽錯誤!未定義書簽。一6）式結(jié)合，得出ln（R） 01 cosh 1 ）dA1伙MER

9、GEFORMAT錯誤!未定義書簽。（錯誤！未定義書簽。一錯誤!未定義書簽。一6）令3 o'cosh1（衛(wèi)）d ，可以得到:1r1 R exp（s ）衣MERGEFORMA錯誤!未定義書簽。-1-錯誤！未定義書簽。）這就是射線最深點的半徑的表達式。此表達式可以根據(jù)觀測走時曲線 t （）得出。如果t（）為一平滑曲線，我們可由其切線斜率得到p（）（由于p=也）。dd伙MERGEFORMAT （錯誤!未定義書簽。一錯誤!未定義書簽。一6）式中的-JX1Pi()1,為在震中距1的走時曲線斜率.S1的積分式由震中距A從 0到A 1d的p值，得到S1后，根據(jù)伙 MERGEFORMAT簽。)式求出A，

10、即其最深轉(zhuǎn)折點。根據(jù)轉(zhuǎn)折點為r1的1 土，得出其對應此深度的速度V1V1j1伙MERGEFORMAT錯誤!未定義書簽。(錯誤味定義書簽。-1-7 )依此類推，得出速度與深度的分布圖F面歸納一下求速度分布的具體步驟1、做出0-范圍內(nèi)的走時曲線，即t-曲線，震源要取在地表上，否則 要進行修正；在該步驟中，通常在計算機上利用三次自然樣條函數(shù) 對一元函數(shù)進行成組插值及微分，從而使得走時曲線具有較好的平 滑性能，采用該方法能保證所插值的函數(shù)及其一階導數(shù)、二階導數(shù) 連續(xù)。2、由t-曲線求曲線上各點的斜率，從而做出 空 曲線；該步驟求出pd隨震中距的分布。3、在0-范圍內(nèi)任取1，取出1點對應的dtd 1 ?

11、這就是1dTP1 ()1 °d4、將dt曲線中d0-1段除以常數(shù)理d即11，得到，即衛(wèi)曲線，從而求出cosh1 p 1 曲線;這就是cosh 1(-p)隨震中距的變化。15、在0 1范圍內(nèi)求出cosh 1 p 1 曲線下的面積,即求cosh 1 蘭 蘭 d ，這就是1 cosh 1 （ ）d ;在計算積分的過 0dd i0i程中，通常采用Simpson積分數(shù)值公式進行。6、根據(jù) 伙MERGEFORMAT （錯誤 未定義書簽。-1-錯誤!未定義書簽。）式求 出r1，然后根據(jù) 衣MERGEFORMAT （錯誤!未定義書簽。-1-7 ）式得到r1所 對應的速度v。7、取不同的1值，重復上述

12、步驟，從而求出不同的R值所對應的速度值 w ;以一個具體地區(qū)的走時曲線為例子：采用的數(shù)據(jù)是北京臺網(wǎng)、大連、長春、牡丹江等26個臺站的記錄，并選用從北京地區(qū)到阿留申群島西端的300多個地震，經(jīng)過震源深度及剝殼校正后 （地球平均厚度35公 里）。得到的北京一薩哈林一線的P波走時曲線。數(shù)據(jù)如表 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 % 調(diào)用數(shù)據(jù) 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) ;%對數(shù)值進行樣條插值，得到插值點的y值，一階導數(shù)dy1,二階導數(shù)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);%視為常數(shù)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 )；%采用球?qū)覵nell定律得到最深點的速度；deg2rad()函數(shù)的作用是：將角度轉(zhuǎn)換為弧度，利用公式v -1%序號+1lndx=lndx+1 ；%pause endplot(v , r0 r1 , '.' )%繪制速度和深度剖面圖set(gca,'ydir'， reverse' ) %將y軸的方向進行翻轉(zhuǎn)利用上述程序調(diào)用表81-1中的北京一薩哈林地區(qū)的P波走時曲線的數(shù)據(jù)。對

17、 觀測曲線進行上述程序來反演即可得到速度剖面如圖8-1-2北京一薩哈林地區(qū)P波走時觀測反演速度剖面.8 1-1北京一薩哈林地區(qū)P波走時觀測反演速度剖面對反演數(shù)據(jù)進行重新插值處理，選擇出合適的數(shù)據(jù)，分別選擇如下的數(shù)據(jù)(km): 20,45, 85,125,145,165 265, 365, 405, 445,495,565, 705,875,965,1025,1165在上述MATLA函數(shù)后加入如：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'并利用得到的數(shù)據(jù)進行繪圖得到圖8 1-2定點插值的反演數(shù)據(jù)圖圖81-2定點插值的反演數(shù)據(jù)圖當有低速區(qū)時，這些公式是有問題的。在這種情況下，（P）不連續(xù)，沒有唯一的解。對這種情況，傳統(tǒng)的做法是設想在低速區(qū)存在若干有不同速度的均勻 層。由于沒有射線在這些層里轉(zhuǎn)折，因此低速區(qū)的這些層可以任意移來移去，穿 過低速區(qū)的射線的積分的走時和震中距不會改變。不論解析多么完全，在地震學中很少用赫格勞茲一維歇爾 （HW）公式，這至 少有兩個原因：首先，HW假定我們有已知的連續(xù)的t（）曲線，而實際上我們總是只有有限的走時點。這意味著

19、走時曲線必須在這些數(shù)據(jù)之間進行內(nèi)插，這些內(nèi)插方案的不同導致于有不同的速度剖面實際上，會有無數(shù)個稍有不同的速度模型 它們都適合于有限數(shù)目的t()點.然而，更嚴重的問題是實際的地震數(shù)據(jù)一般多 少都有噪聲，自身相互矛盾。圖8 1 3展示了實際數(shù)據(jù)的典型的例子。在左 邊的例子中，小的時間變動導致t()點的波動，不可能把這些點與能得到期望的1 D速度模型的光滑的、物理上可靠的t()曲線相聯(lián)系。在右邊的例子中，我們 注意到可以把幾次地震的數(shù)據(jù)結(jié)合在一起，用這些數(shù)據(jù)確定“平均”的速度剖面.圖8 1 3呈現(xiàn)離散的走時觀測結(jié)果往往使速度剖面的反演復雜化由于這兩個原因，HW公式不能直接應用.下面我們對地震學家反演

20、這些不完 善的數(shù)據(jù)的一些方法進行討論。HW公式的主要用途是闡明精確的t(X)曲線可給 出速度剖面的唯一的解.因此，許多反演的策略是把尋找最佳速度模型的問題簡 化為尋找最佳的t()曲線的問題。8.2錯誤!未定義書簽。單臺地震定位及INGLADA法定位根據(jù)走時數(shù)據(jù)確定地震的位置是地震學中的一個“古老"的，富有挑戰(zhàn)性的問題，一直是地震學研究的一個重要方面.地震通常按發(fā)震時間和震源位置來定 義.震源是地震的位置(x,y,z),而震中是震源正上方地表上的一個點(x, y)。在地 震定位中，通常把地震作為點源來對待。對于破裂幾十到幾百公里的大地震，震源 不一定是地震的“中心”，寧可把它看成是地震

21、發(fā)生時，地震能量開始輻射的那 個點。由于破裂的速度小于S波速度，可以不管地震最終的尺度和持續(xù)時間，只 根據(jù)初至到時來確定震源.標準目錄給出的地震資料，例如震中的初步確定(PDE) 依據(jù)的是高頻體波的走時，不能把這些發(fā)震時間和震源位置與長周期反演結(jié)果相 混淆，后者往往給出地震的矩心時間和位置，表示整個地震的平均發(fā)震時間和位 置本節(jié)我們首先討論單臺定位法和多臺定位的 Inglada方法。821單臺定位法單臺定位的原理就是根據(jù)縱波初動確定出震中方位角，然后根據(jù)震相到時（走時 表等）確定出震源到臺站的距離，根據(jù)水平向和垂直向記錄的數(shù)值，可以求得視入射角， 從而根據(jù)地震波速度估計出真入射角。根據(jù)震源到臺

22、站的距離和真入射角求得地震的震 中距和震源深度。注意震中方位角是指通過臺站子午線與地震臺站到震中連線間的夾 角，沿順時針量取為正。由圖8-2 1所示,設地震臺站的P波的北向初動半周期位移為 uN，東向初動半周期位移為uE，則水平向位移矢量與北向的夾角可以由下式來計算uEarctanUn錯誤!未定義書簽。（錯誤！未定義書簽。-2-錯誤!未定義書簽。）在MATLAB中可以用atan2（Ue,un）來求得，因為atan2函數(shù)根據(jù)的正負號給出到 的范圍，將其轉(zhuǎn)換為度，并且對于小于0度的角度加上360度，就轉(zhuǎn)換為0360度范圍。知道了水平向矢量方向，地震震中應該位于該矢量方向上，但究竟是沿著矢量 方向還

23、是與矢量方向相反？這需要采用垂直向記錄。我們定義垂直向記錄山向上為正。垂直方向的初動半周期決定地震射線的方向當垂直向初動向上時，地震射線的方向背向震中；當垂直初動向下時，地震波射線的方向指向震中。圖8-2-2為一斷層錯動圖，臺1觀測到的是壓縮波（+），其P波垂直向記 錄向上，即射線是“離源”（即背向震中）運動；此時震中位于水平矢量的相反方向上。 臺2觀測到的是膨脹波（-），其P波垂直向初動向下，即波射線為“向源”（指向震 中）的運動.此時震中位于水平矢量方向上。根據(jù)這種分析，我們可以得到結(jié)論：垂直方向初動向上的臺站，震中位于水平矢量方向的相反方向上；反之垂直方向初動向上的臺站，震中位于水平矢量

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顯示文本”不能橫跨多行！該物理量具有速度的量綱,通常被稱為虛波速度。設若可以根據(jù)當?shù)氐钠骄貧?S波速度和P波速度估計該參數(shù)Vp = 5.0 km/s , Vs = 3。0 km/s 則虛波速度 k=7.5km/s。根

25、據(jù)臺站上讀取的該地震的S波和P波到時差和該地區(qū)的虛波速度，可求得震源到臺站的距離r，公式為r k ts tp（錯誤!未定義書簽-錯誤!未定義書簽。-錯誤!未定義書簽。）第三步：地震記錄的垂直向和水平向的記錄，求得視入射角，從而根據(jù)地球均勻速 度模型估計真入射角。地震在該臺站的水平向位移由南北向和東西向位移記錄合成:伙MERGEFORMAT錯誤!未定義書簽。（錯誤!未定義書簽。-2 錯誤！未定義書簽。）則視入射角為：rparctan 比Uu伙MERGEFORMAT （0錯誤!未定義書簽。-10）根據(jù)第四章的真入射角和視入射角的轉(zhuǎn)換公式（4243）,即sin ip sin得p 2到真入射角。第四步

26、根據(jù)真入射角和震源到臺站的距離r求得震中距 和震源深度h。r si nip伙MERGEFORMAT （錯誤!未定義書簽錯誤!11)h r cosi p伙MERGEFORMAT （錯誤!未定義書簽簽。震中方位角及震中距求得以后，便可決定出震中位置。這樣應用一個臺站的三分向記錄就可以估計地震的大概位置了。在計算機高速發(fā)展的今天，計算機可以代替人工操作。我們通常給出臺站的經(jīng)緯度， 需要求解地震的經(jīng)緯度和深度在近震的情況下，通常采用直角坐標來解決問題，因此需 要進行經(jīng)緯度轉(zhuǎn)換為直角坐標系的轉(zhuǎn)換。將觀測點的經(jīng)緯度i, i轉(zhuǎn)換為直角坐標 x , yi有很多成熟的公式可以使用(如，朱介壽等，1988 )。由

27、于本研究的區(qū)域范圍較小，我們采用下述簡單公式(時振 梁等,1990 )。將觀測點的經(jīng)緯度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 )%艮據(jù)坐標系原點lat0和long0，將地理坐標(lat , long )轉(zhuǎn)換為直角坐標x,y%輸入lat , long分別緯度和經(jīng)度%

28、lat0 和long0為坐標原點%輸岀x和y為得到的直角坐標，x沿經(jīng)向方向的坐標(東向 )，y沿緯向方向(北向)c1 = 111。199;% 1 度弧長y=c1* (lat-lat0) ;%北向坐標(8-2 9)式第一式rlatP=deg2rad (lat+lat0 ) /2 ；%緯度的平均值作為計算經(jīng)度的大圓的半徑x=c1 *( long-long0)*cos(rlatP);%經(jīng)向坐標(8 2 9)式第二式return采用上面的程序得到震源位置在直角坐標系的表示后，還需將直角坐標系中的坐標Xi,yi轉(zhuǎn)換為經(jīng)緯度i, i ,公式為:ii111.1990X0111.199 cos2顯示文本”不能

29、橫跨多行！其中的符號及意義與上式相同。采用MATLAB進行轉(zhuǎn)換的程序如下:functionlatP,longP =dist2geo(x , y,lat0,long0)%根據(jù)坐標系原點lat0和long0，將直角坐標x,y轉(zhuǎn)換為地理坐標 %輸入x , y分別為緯度方向和經(jīng)度方向的直角坐標距離% lat0和long0為坐標原點% 輸出 latp 和 longp 為得到的直角坐標點的緯度和經(jīng)度c1 = 111.199;%1度弧長latP=y/c1+lat0 ； (x,y )點的緯度，( 8-2 10)式第一式rlatP=deg2rad (latP+lat0 )/2 ；%緯度的平均值作為計算經(jīng)度的大圓

30、的半徑IongP=x/cos(rlatP)/c1+longO；%(x, y)點的經(jīng)度，(8210 )式第二式return為驗證上面的程序， 我們以北緯 25°,東經(jīng) 120° 作為坐標原點， 求北緯 26° ，東經(jīng) 121°的直角坐標， 運行“ x,y=geo2dist(26,121 ，25，120)", 則可以得到 x= 100.37km，y=111。2km的結(jié)果，然后用此結(jié)果求解其經(jīng)緯度，運行"latp , Iongp =dist2geo(x ,y,25,120 )”的語句 , 可以得到北緯 26°,東經(jīng) 121

31、6;的結(jié)果。讀者也可以采用其他的數(shù)據(jù)進行檢驗。有了上面的坐標轉(zhuǎn)換程序，我們可以編寫根據(jù)單臺記錄求解地震位置的計算機程序為：Slat=35 ； Slon=115 ；%臺站的經(jīng)緯度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) ) ; %根據(jù)( 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('地震的經(jīng)度東經(jīng) 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三站求取震中位置的石川法以三站為一組的作法，稱為虛波速度法，亦稱石川法。虛波速度在近震范圍內(nèi)是相當穩(wěn)定的，可由本地經(jīng)驗得到，一般都在8公里/秒左右設選定的三個觀測站為 Si，S2，S3,先將它們按照地理位置，標在地圖上 ，如圖82-2(a )所示，再按各站 記錄到的(ts tp)到達時間差，用公式(錯誤！未定義書簽。-錯誤!未定義書簽。-錯誤!未定義書簽。)求得各站相應的震 源距離：Di, D2，D3。然后以Si為中心，以Di為半徑，在操作地圖上作第一地震 圓，其實應為一個隱藏在下的半球面震源應當是此球面上的一點

35、同樣以S2及S3為中心，以及為半徑作第二、第三地震圖，各有其下隱的半球面。這三個半球面互 相交切，每兩個構(gòu)成一個弧形交切痕，共有三條，投影到圖上為：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之長折合成公里數(shù)，便得到震源深度的數(shù)據(jù)，如圖8-2-3.(a)(b)圖8 2-3三站交切法測定震源位置示意圖四個臺站（或以上）求解地震位置的INGLADA算法的坐標為（x,yi,z）,P波、S采用均勻介質(zhì)模型，設P波、設震源坐標為（X0,y°,z0）,發(fā)震時刻為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 （錯誤!未定義書簽。一錯誤!未定義書簽。-錯誤!未定義書簽。） 的系數(shù)矩陣具有奇異性，導致Zo解發(fā)散?？梢酝ㄟ^其它方法單獨求解 Zo：取震中解&, y°，將其代入* MERGEFORMAT（0-錯誤!未定義書簽。-錯誤!未定義書簽。）式求得 將各個臺站求得的Zo的平均值作為震源深度的解。這樣Zo的誤差僅決定于震中的定位誤差和R的到時讀數(shù)誤差（二者均是可控的），從根本上解決了無法得到震源深度有效信息的問題。關于發(fā)震時刻，有方程t0tpi旦Vp錯誤!未定義書簽。（0 2 14

40、）將* MERGEFORMAT214）用于各個臺站，然后求各臺站得到的t。的平均值 即可作為發(fā)震時刻。由于R完全取決于到時（0-錯誤！未定義書簽。-i2 ）式）， 所以to和震源的求解完全獨立，其誤差的唯一來源是到時讀數(shù)的誤差 ，只要到時讀數(shù)足夠精確，即使震源定位誤差很大，也可求得很精確的to。根據(jù)上述方法,可以初步求解出地震發(fā)生的經(jīng)度、緯度、深度和發(fā)震時間的估計解下面為求解上述問題的計算機程序: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循環(huán)結(jié)束end % for ii循環(huán)結(jié)束FirLocal = pinv(RMat) *LMat; % 公式 x=A-1B, pinv (x)為求解矩陣 x的偽逆%發(fā)震時刻和地震深度沒有采用前面的計算結(jié)果，下面分別計算for

43、ii=1length (Pg)SeiDist =VSpeed*(Sg(ii) Pg (ii) ； % 震中距SeiTime（ii） = Pg（ii ） SeiDist/Vp ；%計算發(fā)震時刻SeiDeeps(ii)=Stz (ii ) + sqrt(SeiDistA2- (FirLocal (1) - Stx(ii ) )A2 - (FirLocal(2) Sty ( ii）A2）;%根據(jù)（8-2 7）式計算地震深度end %for ii循環(huán)結(jié)束FirLocal（3） = real（ mean（SeiDeeps ） ; 估計地震深度的平均值if （ FirLocal（3 ） 0 | FirLo

44、cal（ 3 ） > 100）%如果地震深度超出可信范圍，則強制給定地震深度FirLocal(3) = 12;endFirLocal（ 4） = mean （ SeiTime） ； %估計地震發(fā)震時刻的平均值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 ） % 調(diào)用和達法求解地震位置在上面的程序中， 我們均采用直角坐標 . 如果已知的臺站坐標的經(jīng)緯度， 求解地震的經(jīng)緯度， 則需要采用（ 8-29）式首先轉(zhuǎn)換為直角坐標,待求得地震位置后再采用（82-10 ）式轉(zhuǎn)換為地震的經(jīng)緯度地理坐標 .錯誤 !未定義書簽。 8.3 地震定位結(jié)果與臺站分布的關系 地震定位問題需要求解地震的空間位置和發(fā)震時刻 .我們把這些參數(shù)稱為模型， 定義

46、模型矢量：m (g,m2,m3,m4)(T,x, y,z)衣MERGEFORMAT (錯誤!未定義書簽。一錯誤!未定義書簽。一15)現(xiàn)在假定我們有地震臺觀測的n個震相觀測到時.為了得到參數(shù)m，我們首先必 須假定一個參考的地球模型，如一維平層介質(zhì)模型。這樣，按照第6 7章的知識, 根據(jù)m值和臺站位置可以計算第i個臺站的預測到時：tip Fi (m)錯誤!未定義書簽。(錯誤!未定義書簽。一錯誤!未定義書簽。一15)這里F是算子，其復雜程度由假定地球模型的復雜程度和地震波傳播的射線路徑?jīng)Q定。觀測到時與預測的理論到時之差為：ri ti 丁 ti Fi (m)衣MERGEFORMAT (錯誤!未定義書簽

47、。-錯誤!未定義書簽。一錯誤!未定義書簽。)這里ri是第i個臺站的殘差.在某種意義上來說，我們希望得到使觀測到時與預測 到時的殘差為最小的m.注意，F(xiàn)是地球模型和每個臺站位置的函數(shù)，根據(jù)第 6, 7章的知識就可以精確計算。然而，走時與地震位置參數(shù)的非線關系使反演最佳 地震模型的工作大大復雜化。這種非線性的影響基至在均勻速度進行定位的簡單 情況下，也是明顯的，在這種情況下，座標為(xi, yi)的臺站記錄震源點(x, y)發(fā)震 時刻為T的到時為：ti(x x)2 (y yj2衣MERGEFORMAT錯誤!未定義書簽。(錯誤!未定義書簽。一錯誤!未定義書簽。-16)這里 是速度。在此方程中，到時t

48、不論是與x還是y都不成線性關系.該結(jié)果表 明，我們不能用解線性方程的方法來得到震源位置?，F(xiàn)在有了功能強大的計算機，我們可在所有可能的地震位置和發(fā)震時間范圍 里作網(wǎng)格搜索，計算每個臺站預測的到時。那么，我們可以求得使預測的時間tip 與觀測的時間ti最一致的特定的m。怎樣規(guī)定“最”一致呢？流行的選擇是最小 二乘法，即求使下式達最小的 m值：tii 1ti顯示文本”不能橫跨多行這里n是臺站數(shù)目。把平均的平方殘差力叫做方差。在數(shù)理統(tǒng)計上更普遍的形式是用ndf來定義方差這里ndf是自由度的數(shù)目（ndf為n減去擬合中的獨立參數(shù)的數(shù)目）對通常需要解決的問題，擬合參數(shù)的數(shù)目比數(shù)據(jù)的數(shù)目少得多 所以n和ndf

49、近似相等。在數(shù)理統(tǒng)計中2分布于自由度個數(shù)及置信度的關系如表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中非相關的隨機噪聲所造 成的，那么最小二乘法將有助于給出正確的答案。在本門課程的討論中，如不特殊說明，我們假定誤差是由服從高斯正態(tài)分布的隨機噪聲導致的?，F(xiàn)在我

50、們以一組臺站，假設為均勻的地殼模型，以1km和0.1s為網(wǎng)格點搜索地震的位置。MATLAB程序如下：clcx0=30 ； y0=29;%假定的震源位置，發(fā)震時刻假定為零S1%x0=50 ; y0=20; %假定的震源位置，發(fā)震時刻假定為零S2% x0=10;y0=60; % 假定的震源位置，發(fā)震時刻假定為零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 )； 臺站個數(shù)Pt=sqrt ( st(:,1 ) -x0 )。A2+ (st(: , 2) y0)。A2)./vel+ ( 1-rand(N , 1) *0.5;%計算各臺站的預測走時，并加上幅度為 0.5s 的誤差ndf=N3；%該問題需要求解地震的發(fā)震時刻和 x,y 坐標，因此有三個未知數(shù)resmin=1 。 0E50;%給出最小值是一個較大的數(shù)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%輸出找到的最優(yōu)結(jié)果orig=origmin;%將最優(yōu)的發(fā)震時刻保存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