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

1、錯(cuò)誤 !未定義書簽。錯(cuò)誤 !未定義書簽。第八章 走時(shí)數(shù)據(jù)的反演8.9.8.1球?qū)咏橘|(zhì)的速度反演.9.0.9.0.9.1.8.2單臺(tái)地震定位及 Inglada 法定位.1.0.0.8 。1.1 拐點(diǎn)法（古登堡法 ） 8.1。2 Hergloz -Wiechert 方法1.0.0.1.0.6.8.38.2。 3 采用三臺(tái)站求解地震位置的計(jì)算機(jī)算法地震定位結(jié)果與臺(tái)站分布的關(guān)系1.0.7.1.1.0.8.2.1 單臺(tái)定位法8.2.2 三站求取震中位置1.1.4.1.1.4.1.1.5.1.3.2.8 。 4 迭代定位方法8.4.1 迭代定位算法的基本思路8。4。 2 迭代方法的具體實(shí)現(xiàn)8 。 5 相對(duì)

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

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

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

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

6、n i其中 二，r1為射線轉(zhuǎn)折點(diǎn)距地心的距離，p= , R為地球半徑。改變積分變數(shù)vv(r)則c 01dr ,2pd2 2 d p r ip d(0-錯(cuò)誤!未定義書簽。4)其中在此介紹一積分式,錯(cuò)誤!未定義書簽其中i匕p,，為轉(zhuǎn)折半徑為*的射線參數(shù)，所以Vi伙MERGEFORMAT （錯(cuò)誤!未定義書簽。-1-4）式表示從轉(zhuǎn)折半徑為 0至到轉(zhuǎn)折半徑為r,的所有射線的積分 將 伙MERGEFORMAT （錯(cuò)誤!未定義書簽。-1-4）式作用于震中距表 達(dá)式 伙MERGEFORMAT （0-錯(cuò)誤!未定義書簽。一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扌簽。扌簽。（錯(cuò)誤!未定義書簽。-1-錯(cuò)誤!未定義書簽。）式的第一項(xiàng)代入積分上下限，由于 cosh ( 1)=0 ,且在p= 0，半徑為R,A =0，因此第一項(xiàng)消失僅剩下第二項(xiàng)，可寫成:1cosh0MERGEFORMAT錯(cuò)誤!未定義書簽。（錯(cuò)誤!未定義書簽。一錯(cuò)誤!未定義書簽。一6）（錯(cuò)誤!未定義書簽。-錯(cuò)誤！未定義書簽。一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錯(cuò)誤!未定義書簽（錯(cuò)誤！未定義書簽。一1-6）與伙MERGEFORMAT （錯(cuò)誤！未定義書簽錯(cuò)誤!未定義書簽。一6）式結(jié)合，得出ln（R） 01 cosh 1 ）dA1伙MER

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

10、即其最深轉(zhuǎn)折點(diǎn)。根據(jù)轉(zhuǎn)折點(diǎn)為r1的1 土，得出其對(duì)應(yīng)此深度的速度V1V1j1伙MERGEFORMAT錯(cuò)誤!未定義書簽。(錯(cuò)誤味定義書簽。-1-7 )依此類推，得出速度與深度的分布圖F面歸納一下求速度分布的具體步驟1、做出0-范圍內(nèi)的走時(shí)曲線，即t-曲線，震源要取在地表上，否則 要進(jìn)行修正；在該步驟中，通常在計(jì)算機(jī)上利用三次自然樣條函數(shù) 對(duì)一元函數(shù)進(jìn)行成組插值及微分，從而使得走時(shí)曲線具有較好的平 滑性能，采用該方法能保證所插值的函數(shù)及其一階導(dǎo)數(shù)、二階導(dǎo)數(shù) 連續(xù)。2、由t-曲線求曲線上各點(diǎn)的斜率，從而做出 空 曲線；該步驟求出pd隨震中距的分布。3、在0-范圍內(nèi)任取1，取出1點(diǎn)對(duì)應(yīng)的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 ;在計(jì)算積分的過 0dd i0i程中，通常采用Simpson積分?jǐn)?shù)值公式進(jìn)行。6、根據(jù) 伙MERGEFORMAT （錯(cuò)誤 未定義書簽。-1-錯(cuò)誤!未定義書簽。）式求 出r1，然后根據(jù) 衣MERGEFORMAT （錯(cuò)誤!未定義書簽。-1-7 ）式得到r1所 對(duì)應(yīng)的速度v。7、取不同的1值，重復(fù)上述

12、步驟，從而求出不同的R值所對(duì)應(yīng)的速度值 w ;以一個(gè)具體地區(qū)的走時(shí)曲線為例子：采用的數(shù)據(jù)是北京臺(tái)網(wǎng)、大連、長春、牡丹江等26個(gè)臺(tái)站的記錄，并選用從北京地區(qū)到阿留申群島西端的300多個(gè)地震，經(jīng)過震源深度及剝殼校正后 （地球平均厚度35公 里）。得到的北京一薩哈林一線的P波走時(shí)曲線。數(shù)據(jù)如表 8 1 1為北京-薩哈林剖面P波走時(shí)觀測(cè)值。表8-1-1北京-薩哈林剖面P波走時(shí)觀測(cè)值? (°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然后將走時(shí)曲線用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) ' %走時(shí) dh=0。 2;x=0.2:dh:47.8;n=length(x )；y,dy1, dy2 =splineinsert (xO,yO, x) ;%對(duì)數(shù)值進(jìn)行樣條插值，得到插值點(diǎn)的y值，一階導(dǎo)數(shù)dy1,二階導(dǎo)數(shù)dy2figure ( 1)%實(shí)際觀測(cè)曲線plot ( x,y, ' r* 'grid ontitle ('實(shí)際觀測(cè)曲線)xla

15、bel('震中距)；ylabel('走時(shí)');R=637135；In dx=1 ；r1=zeros(1，(length (3: 2: n) )；%射線能夠達(dá)到的最低點(diǎn)的半徑v=zeros(1，(length(3:2:n ) ;%射線最低點(diǎ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 公式進(jìn)行積分 :f(x)dx f(a) 4f(

16、) f(b)%采用公式求解射線最深點(diǎn)的半徑r1 (Indx)=R/exp (Int/pi)；r1 R exp(3)v(lndx ) =deg2rad (r1(Indx)/ksi1 )；%采用球?qū)覵nell定律得到最深點(diǎn)的速度；deg2rad()函數(shù)的作用是：將角度轉(zhuǎn)換為弧度，利用公式v -1%序號(hào)+1lndx=lndx+1 ；%pause endplot(v , r0 r1 , '.' )%繪制速度和深度剖面圖set(gca,'ydir'， reverse' ) %將y軸的方向進(jìn)行翻轉(zhuǎn)利用上述程序調(diào)用表81-1中的北京一薩哈林地區(qū)的P波走時(shí)曲線的數(shù)據(jù)。對(duì)

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

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

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

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

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

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

24、方向上。圖8-2-2兩個(gè)臺(tái)站接收到地震 P波初動(dòng)符號(hào)，判斷地震所在方位的示意圖第二步，確定震源到臺(tái)站的距離VpVsr-設(shè)ts, tp分別為該臺(tái)站S波及P波的到時(shí)，則地震到達(dá)該臺(tái)站的 S波和P波的走時(shí) 之差等于S波及P波的到時(shí)差，則有：, r r 11 ts tprVsVpVs Vp伙MERGEFORMAT （0-錯(cuò)誤!未定義書簽錯(cuò)誤!通常定義1,11VpVskVs VpVp Vs顯示文本”不能橫跨多行！該物理量具有速度的量綱,通常被稱為虛波速度。設(shè)若可以根據(jù)當(dāng)?shù)氐钠骄貧?S波速度和P波速度估計(jì)該參數(shù)Vp = 5.0 km/s , Vs = 3。0 km/s 則虛波速度 k=7.5km/s。根

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

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

27、于本研究的區(qū)域范圍較小，我們采用下述簡(jiǎn)單公式(時(shí)振 梁等,1990 )。將觀測(cè)點(diǎn)的經(jīng)緯度i, i換算為直角坐標(biāo) xi, yi的公式如下：y 111.199x 111.199i00 cos -2伙MERGEFORMAT錯(cuò)誤!未定義書簽(0錯(cuò)誤!未定義書簽式中0, 0為直角坐標(biāo)的原點(diǎn)。0, 0, i, i以度為單位，Xi,yi的單位為km。其MATLAB 程序如下：functionx , y =geo2dist (lat , Iong,latO, longO )%艮據(jù)坐標(biāo)系原點(diǎn)lat0和long0，將地理坐標(biāo)(lat , long )轉(zhuǎn)換為直角坐標(biāo)x,y%輸入lat , long分別緯度和經(jīng)度%

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

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

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

31、6;的結(jié)果。讀者也可以采用其他的數(shù)據(jù)進(jìn)行檢驗(yàn)。有了上面的坐標(biāo)轉(zhuǎn)換程序，我們可以編寫根據(jù)單臺(tái)記錄求解地震位置的計(jì)算機(jī)程序?yàn)椋篠lat=35 ； Slon=115 ；%臺(tái)站的經(jīng)緯度ue=3； un=3;uu=-2; 地震臺(tái)記錄的南北分向、東西分向和上下分向Parrival=5 ； Sarrival=8 ；% P波和S波到時(shí)，單位為秒vp=5;vs=3 ；%平均的P波速度和S波速度，單位為km/sk=(vp vs)/(vp vs)；%虛波速度azi0=atan2(ue ， un) ；%按( 82-1 )式求解方位角if (uu<0 )azi=azi0 ；disp ('地震相對(duì)于臺(tái)站的方

32、位角為：' ,num2str(rad2deg(azi) ) );elseazi=2 pi-azi0;disp( '地震相對(duì)于臺(tái)站的方位角為：' ,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 )式求震源到臺(tái)站的距離h=r*cos ( ip ) ;%求地震的

33、深度delta=r*sin ( ip )；%求地震的震中距x=delta*s in (azi);%得到以臺(tái)站為坐標(biāo)原點(diǎn)的 x坐標(biāo)y=delta*cos(azi ) ;%得到以臺(tái)站為坐標(biāo)原點(diǎn)的y坐標(biāo)Elat ， Elon=dist2geo(x,y ， Slat,Slon);disp(sprintf('地震的經(jīng)度東經(jīng) f度，緯度為北緯 f度，深度 f km' ,Elon,Elat , h)；disp ( sprintf( '地震到臺(tái)站的震中距為% fkm' ,delta) )disp(sprintf('P 波的走時(shí)為 1秒'，r/vp )disp(s

34、printf(' S波的走時(shí)為 %f秒' , r/vs)8。2.2三站求取震中位置的石川法以三站為一組的作法，稱為虛波速度法，亦稱石川法。虛波速度在近震范圍內(nèi)是相當(dāng)穩(wěn)定的，可由本地經(jīng)驗(yàn)得到，一般都在8公里/秒左右設(shè)選定的三個(gè)觀測(cè)站為 Si，S2，S3,先將它們按照地理位置，標(biāo)在地圖上 ，如圖82-2(a )所示，再按各站 記錄到的(ts tp)到達(dá)時(shí)間差，用公式(錯(cuò)誤！未定義書簽。-錯(cuò)誤!未定義書簽。-錯(cuò)誤!未定義書簽。)求得各站相應(yīng)的震 源距離：Di, D2，D3。然后以Si為中心，以Di為半徑，在操作地圖上作第一地震 圓，其實(shí)應(yīng)為一個(gè)隱藏在下的半球面震源應(yīng)當(dāng)是此球面上的一點(diǎn)

35、同樣以S2及S3為中心，以及為半徑作第二、第三地震圖，各有其下隱的半球面。這三個(gè)半球面互 相交切，每兩個(gè)構(gòu)成一個(gè)弧形交切痕，共有三條，投影到圖上為：aa， bb， cc。這 三條弧相交于一點(diǎn)，因它是三個(gè)下隱半球面的共同點(diǎn)，當(dāng)然就是震源，投影到圖上為E,就是震中，便可在操作地圖上定出其地理位置。過震中E,垂直于觀測(cè)點(diǎn)Si (S2或S3亦可)與震中的連線，作一直線與第一地 震圓交于A，則EA之長就是震源深度h。這個(gè)從圖822的右圖(b )清楚的看 至V，bi是第一地震圓的半球面，今若過 E并垂直SiE于作一垂直面切下，則這個(gè)垂 面必通過震源F,并與地震圓的半球面交割成一半圓形割痕，為b2,很顯然A

36、E=EF=h，都是這半圓形的半徑最后用操作地圖的比例尺，將AE之長折合成公里數(shù)，便得到震源深度的數(shù)據(jù)，如圖8-2-3.(a)(b)圖8 2-3三站交切法測(cè)定震源位置示意圖四個(gè)臺(tái)站（或以上）求解地震位置的INGLADA算法的坐標(biāo)為（x,yi,z）,P波、S采用均勻介質(zhì)模型，設(shè)P波、設(shè)震源坐標(biāo)為（X0,y°,z0）,發(fā)震時(shí)刻為t°，臺(tái)站i波到時(shí)分別為tpi和tsi,震源到臺(tái)站i的距離為RS波速度分別為Vp和Vs ,引入虛波速度VsVp1VpVSVp Vs，則有t R R Rtpi錯(cuò)誤!未定義書簽。（0-錯(cuò)誤!未定義書簽錯(cuò)誤!未定義書簽。）從而V tSitPiMERGEFORMA

37、T錯(cuò)誤!未定義書簽。（0-錯(cuò)誤!未定義書簽。-12）對(duì)臺(tái)站i有2 2 2 2(xo Xi)(yo y)(z。zj R錯(cuò)誤!未定義書簽。（0 錯(cuò)誤!未定義書簽。-錯(cuò)誤!未定義書簽。）對(duì)臺(tái)站j有2 2 2 2(Xo Xj)(yo yj)(zo Zj)Rj衣MERGEFORMAT （錯(cuò)誤!未定義書簽。-錯(cuò)誤!未定義書簽。-錯(cuò)誤!未定義書簽。）將方程 衣MERGEFORMAT1錯(cuò)誤！未定義書簽。-錯(cuò)誤！未定義書簽。）、衣MERGEFORMA錯(cuò)誤！未定義書簽。-錯(cuò)誤！未定義書簽。-錯(cuò)誤！未定義書簽。） 展開后相減，得到線性方程R21 2 2 2 (Xi Xj)Xo (yi yjy° (w zj

38、z°r RjMERGEFORMAT (錯(cuò)誤味定義書簽。一2-錯(cuò)誤!未定義書簽其中ri2x22 2yiz。將方程（錯(cuò)誤!未定義書簽。一2-錯(cuò)誤!未定義書簽。）用于3對(duì)以上的臺(tái)站，即得到一關(guān)于（xo,yo,zo）的線性方程組，方程組寫成矩陣形式為：XiX2yiy2zi Z2XoXiX3yiy3ziZ3yo*4*Zo2 ri222 riRi2Ri2顯示文本”不能橫跨多行衣MERGEFORMAT （錯(cuò)誤!未定義書簽。錯(cuò)誤!未定義書簽。-錯(cuò)誤!未定義書簽。）可以寫為Ax B的形式，從而方程的解就變成了求解 A的逆矩陣的形式，方程的解可以表示為：x A iB。這樣就得到了震源位置由于Zi表示臺(tái)站

39、高度，當(dāng)臺(tái)站海拔高度相差不大時(shí)，乙Zj遠(yuǎn)小于Xi Xj和y yj ，使得方程組* MERGEFORMAT （錯(cuò)誤!未定義書簽。一錯(cuò)誤!未定義書簽。-錯(cuò)誤!未定義書簽。） 的系數(shù)矩陣具有奇異性，導(dǎo)致Zo解發(fā)散?？梢酝ㄟ^其它方法單獨(dú)求解 Zo：取震中解&, y°，將其代入* MERGEFORMAT（0-錯(cuò)誤!未定義書簽。-錯(cuò)誤!未定義書簽。）式求得 將各個(gè)臺(tái)站求得的Zo的平均值作為震源深度的解。這樣Zo的誤差僅決定于震中的定位誤差和R的到時(shí)讀數(shù)誤差（二者均是可控的），從根本上解決了無法得到震源深度有效信息的問題。關(guān)于發(fā)震時(shí)刻，有方程t0tpi旦Vp錯(cuò)誤!未定義書簽。（0 2 14

40、）將* MERGEFORMAT214）用于各個(gè)臺(tái)站，然后求各臺(tái)站得到的t。的平均值 即可作為發(fā)震時(shí)刻。由于R完全取決于到時(shí)（0-錯(cuò)誤！未定義書簽。-i2 ）式）， 所以to和震源的求解完全獨(dú)立，其誤差的唯一來源是到時(shí)讀數(shù)的誤差 ，只要到時(shí)讀數(shù)足夠精確，即使震源定位誤差很大，也可求得很精確的to。根據(jù)上述方法,可以初步求解出地震發(fā)生的經(jīng)度、緯度、深度和發(fā)震時(shí)間的估計(jì)解下面為求解上述問題的計(jì)算機(jī)程序:function FirLocal =ingladawan(Vp, Vs,Pg , Sg,Stx,Sty , Stz )%虛波速度VSpeed= (Vp*Vs)/(Vp-Vs);%計(jì)算地震位置l =

41、1; %方程序號(hào)從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ā)震時(shí)刻和地震深度沒有采用前面的計(jì)算結(jié)果，下面分別計(jì)算for

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

44、cal（ 3 ） > 100）%如果地震深度超出可信范圍，則強(qiáng)制給定地震深度FirLocal(3) = 12;endFirLocal（ 4） = mean （ SeiTime） ； %估計(jì)地震發(fā)震時(shí)刻的平均值return一個(gè)例子Stx=50 0 50 100 100 100 50 0 0;%臺(tái)站位置的 x坐標(biāo)Sty= : 50 100 100 100 50 0 0 0 50；%臺(tái)站位置的 y坐標(biāo)Stz=zeros（1,9）; %臺(tái)站位置的z坐標(biāo)，在此設(shè)置為零Pg= : 6。4 18.5 11。9 11.9 6 。4 11。9 11。9 18.5 15 。5 ; % 各個(gè)臺(tái)站的直達(dá) P波到

45、時(shí)Sg=10.7 30。8 19.8 19 。8 10.7 19.8 19.8 30.8 25。9 ; %各個(gè)臺(tái)站的直達(dá) S波到時(shí)FirLocal =ingladawan （5， 3,Pg ， Sg,Stx,Sty,Stz ） % 調(diào)用和達(dá)法求解地震位置在上面的程序中， 我們均采用直角坐標(biāo) . 如果已知的臺(tái)站坐標(biāo)的經(jīng)緯度， 求解地震的經(jīng)緯度， 則需要采用（ 8-29）式首先轉(zhuǎn)換為直角坐標(biāo),待求得地震位置后再采用（82-10 ）式轉(zhuǎn)換為地震的經(jīng)緯度地理坐標(biāo) .錯(cuò)誤 !未定義書簽。 8.3 地震定位結(jié)果與臺(tái)站分布的關(guān)系 地震定位問題需要求解地震的空間位置和發(fā)震時(shí)刻 .我們把這些參數(shù)稱為模型， 定義

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

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

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

49、近似相等。在數(shù)理統(tǒng)計(jì)中2分布于自由度個(gè)數(shù)及置信度的關(guān)系如表8-3 1所示。表8-3-12分布與自由度和置信度的關(guān)系ndf2(95 %)22 (50%2(5%51。154.3511。07103。949。3418。312010。8519。3431.415034。7649.3367.5010077。9399。33124。 34最小二乘法常被用來求解這類問題， 這是因?yàn)樵谇笞钚≈档膯栴}中，它使方 程有簡(jiǎn)單的解析形式。如果t與tp之間的擬合差是由t中非相關(guān)的隨機(jī)噪聲所造 成的，那么最小二乘法將有助于給出正確的答案。在本門課程的討論中，如不特殊說明，我們假定誤差是由服從高斯正態(tài)分布的隨機(jī)噪聲導(dǎo)致的?，F(xiàn)在我

50、們以一組臺(tái)站，假設(shè)為均勻的地殼模型，以1km和0.1s為網(wǎng)格點(diǎn)搜索地震的位置。MATLAB程序如下：clcx0=30 ； y0=29;%假定的震源位置，發(fā)震時(shí)刻假定為零S1%x0=50 ; y0=20; %假定的震源位置，發(fā)震時(shí)刻假定為零S2% x0=10;y0=60; % 假定的震源位置，發(fā)震時(shí)刻假定為零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 ； %假定的臺(tái)站位置N=length(st )； 臺(tái)站個(gè)數(shù)Pt=sqrt ( st(:,1 ) -x0 )。A2+ (st(: , 2) y0)。A2)./vel+ ( 1-rand(N , 1) *0.5;%計(jì)算各臺(tái)站的預(yù)測(cè)走時(shí)，并加上幅度為 0.5s 的誤差ndf=N3；%該問題需要求解地震的發(fā)震時(shí)刻和 x,y 坐標(biāo)，因此有三個(gè)未知數(shù)resmin=1 。 0E50;%給出最小值是一個(gè)較大的數(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)式求得預(yù)測(cè)到時(shí)for orig=origt-10:0 。 1:origtres=sum(Pt-Pred_torig).A2 )/ndf ； %按照 (8-3 5 )式計(jì)算平均殘差if (res<resmin)resmin=res; xmin=x;ymin=y;origmin=orig；endendendendxmin ， ymin ， origmin ， resmin%輸出找到的最優(yōu)結(jié)果orig=origmin;%將最優(yōu)的發(fā)震時(shí)刻保存figure ( 1)for x=1:1 ：101for y=1 ：1 ： 101Pred_t=sqrt( xst (：,1) ).A2+(y-st( ：,2 )。 A2)/vel ；%按( 8-3 4)式求得預(yù)測(cè)到時(shí)chisq(x,y)=sum(Pt Pred_t orig) 。 A2); %按照 (8-3 5)式計(jì)算 Chisqendendchisqlim=interp1(5,10,20,50,100, : 1。15, 3.94 , 10。85,34.76,7