(完整版)近震定位及其應(yīng)用本科畢業(yè)論文.doc

1、目次摘要 ,Abstract,1前言 ,12近震的理論基礎(chǔ) ,22.1地震學(xué)的基本名詞和概念 ,22.2地震的分類 ,22.3近震的主要震相 ,33正演計算（模型試算） ,53.1基本原理 ,53.1.1坐標(biāo)變換 ,53.1.2正演計算 ,63.2正演模型的建立 ,74反演計算 ,84.1基本原理 ,84.1.1當(dāng)速度 V 未知的初定方法 ,84.1.2當(dāng)速度 V 已知的初定方法 ,104.2模型反演結(jié)果 ,114.3誤差分析 ,115實例剖析 ,136結(jié)束語 ,17參考文獻 ,18附錄 近震定位程序 ,191 前 言四川汶川地震和青海玉樹地震發(fā)生時的地動山搖、房倒屋塌、生死分離的那種凄涼、悲

2、慘的場面，又一次的讓我們感受到自然災(zāi)害的威力，更是讓災(zāi)區(qū)的人們心驚肉跳，恐慌難安，甚至談“地震”色變，這就是地震給人們的最直觀的印象。一次強烈地震的發(fā)生，常伴隨著地面變形和地層錯動，其破壞力是相當(dāng)大的。主要表現(xiàn)為：大型建筑物破壞， 普通民房破壞， 山崩地裂，人畜傷亡。 如今地震發(fā)生越來越頻繁， 2010 年也被“尊稱”為國際地震年。任何事情都具有雙面性，地震對我們的生命安全造成了很大的威脅的同時也為專家學(xué)者們的研究地球提供了豐富的資料。 目前人類對地球內(nèi)部結(jié)構(gòu)的認(rèn)識主要是來自地震學(xué)研究， 因此地震學(xué)也就成為地球科學(xué)中的一個重要學(xué)科。由于天然地震具有很大的能量，它所產(chǎn)生的地震波可以穿透很大的深度

3、，傳播很遠(yuǎn)的距離，而且當(dāng)遇到地球深部的波阻抗差異界面時產(chǎn)生反射波，其中，天然地震的中的近震資料為研究地球深部結(jié)構(gòu)提供了一個重要的信息資源和途徑。利用地震記錄進行定位始源于歐洲和日本。最初使用方位角方法，隨后是幾何作圖法和地球投影法。國際地震匯編（ISS）最先采用最小二乘法進行計算修訂震中。 1961年，博爾特和威爾莫合作，改進了計算方案，并首先在ISS使用，隨后，國際地震中央局與美國海岸和大地測量局先后使用。在我國，李善邦先生最初使用方位角和最小二乘進行地震的觀測和定位。1953年，我國開始采用大量觀測數(shù)據(jù)修訂震中。目前我國的地震的定位方法兼作圖法和計算法。其中近震的定位方法主要有石川法、和達

4、法、高橋法、外心方位角法、假定發(fā)震時刻定位、等時量板法等十幾種方法。計算機的近震定位方法主要分為速度未知的初定為方法和速度已知的初定為方法。目前定位精度已達到較高的水平。如果有合理的臺網(wǎng)分布和適用的走時表，對于5km；遠(yuǎn)震為 510km。一般淺源地震 （h<100km）的震源深度誤差為深度值的10%左右。震源愈深，相對誤差愈小；h>300km時，誤差小于其深度的5%。發(fā)震時刻的誤差為1/10 1/2s.2 近震的理論基礎(chǔ)地震學(xué)主要研究地震的發(fā)生、地震波的傳播及地球內(nèi)部構(gòu)造，地震波能穿透到地球內(nèi)部，并能把地球內(nèi)部的信息帶回地面，具體說來，它主要是根據(jù)天然地震或人工地震的資料，運用物理

5、學(xué)、數(shù)學(xué)及地質(zhì)學(xué)的知識，來研究地震發(fā)生的狀況及地震波傳播的規(guī)律，以求進一步達到研究地震和預(yù)報地震的目的。同時地震是地球表層的振動，是地殼構(gòu)造運動的一種形式，所以還利用地震波的傳播特征來研究地殼和地球內(nèi)部的構(gòu)造。在彈性介質(zhì)內(nèi)傳播的縱波和橫波，由于存在整個彈性空間，因而這些波統(tǒng)稱為體波。相對于體波而言，在彈性界面附近還存在著另一類波動，從能量上而說他們只分布在彈性分界面附近，稱為面波。面波又分為瑞雷波和勒夫波(如右圖 2.1 所示)。為便于對地震的了解研究，現(xiàn)將其基本知識介紹如下。2.1地震學(xué)的基本名詞和概念1震源：地球內(nèi)部發(fā)生地震的地方稱為震源（或稱震源區(qū)）。理論上將震源看成一個點 , 而實際上

6、是一個區(qū)。2震源深度：將震源看作一個點，此點到地面的垂直距離稱為震源深度，一般用字母 h 表示。3.震中：震源在地面上的投影點稱為震中（或稱震中區(qū)）。同時，地面上受破壞最嚴(yán)重的地區(qū)叫極震區(qū)，理論上震中區(qū)和極震區(qū)是相同的，實際上由于地表局部對質(zhì)條件的影響，極震區(qū)不一定是震中區(qū)。與震中相對的地球直徑的另一端稱為對震中；或稱震中對點。4震中距離：在地面上，從震中到任一點沿大圓弧測量的距離稱為震中距離。一般用字母表示。它可用線距離表示，也可用地心所張的角來表示。5發(fā)震時刻：發(fā)生地震的時刻，一般用字母 或 T0 表示。以北京時間標(biāo)出。它比格林威治時間早8 小時。2.2地震的分類為研究方便，按震動的性質(zhì)，

7、可分為天然地震、人工地震及脈動三類。對于天然地震，有下述分類：（一）按震源深度分：1淺源地震：震源深度小于60 公里的天然地震稱為淺震；也稱正常深度地震。2中源地震：震源深度在60 公里至 300 公里之間的地震稱為中源地震。3深源地震：震源深度大于300 公里的地震稱為深震。已記錄到的最深地震的震源深度約700 公里。有時也將中源地震和深源地震統(tǒng)稱為深震。（二）按震中距分：1. 地方震：震中距小于 100 公里的地震。2. 近震：震中距小于 1000 公里的地震。3. 遠(yuǎn)震：震中距大于 1000 公里的地震。2.3近震的主要震相由于近震的震中距很小，其地震波的傳播主要局限于地殼及莫霍面附近，

8、主要震相有 Pg/Sg（）、 P*/S* 、 Pn/Sn、 P*P/S*S 、P*S/S*P 、PmP/SmS及 PmS/SmP等，如圖2.2 所示。圖 2.2 近震震相及其射線路徑示意圖（ c 、 c 分別為康拉德界面及莫霍面的臨界角）1. Pg/Sg震相由于地球淺部存在強烈的橫向及縱向上的速度變化，尤其是有沉積該層地區(qū)的垂向梯度的影響，地震波射線在地殼淺部10 公里的范圍內(nèi)發(fā)生廻折而到達地表，并在震中附近數(shù)公里的范圍，這種波束與直達波，被稱之為廻折或潛波（Diving Wave ），對于縱、橫波來說，分別命名為Pg 和 Sg。2. P*/S*震相P*/S* 分別為地震波在地殼內(nèi)部的康拉德界

9、面產(chǎn)生的折射縱、橫波，它們的速度分別為 6.37.1km/s及 3.63.9km/s 。應(yīng)當(dāng)指出，由于康拉德界面并不是全球普遍存在的界面，在一些地區(qū)觀測不到。3. Pn/Sn震相Pn/Sn 分別為地震波在地殼底部的莫霍界面產(chǎn)生的折射縱、橫波，它們的速度分別為 7.98.2km/s及 4.44.8km/s 。一般來說， Pn 震相在莫霍面的臨界反射以遠(yuǎn)地區(qū)總是已初至波的形式出現(xiàn)，與其它折射波一樣，波至呈現(xiàn)為線性分布。4. P*P/S*S及 P*S/S*P 震相P*P/S*S分別為地震波在地殼內(nèi)部的康拉德界面產(chǎn)生的反射縱、橫波，P*S/S*P分別為地震縱 / 橫波入射地殼內(nèi)部的康拉德界面后發(fā)生轉(zhuǎn)換

10、而產(chǎn)生的反射橫/ 縱波。3 正演計算（模型試算）測定震源位置及發(fā)震時刻最常用的方法是將非線性問題化為線性問題，然后利用最小二乘法原理迭代，全主元消元法等求解。通過反復(fù)修定，以得到震源參數(shù)的最佳值。此方法對近震經(jīng)緯度的確定比較準(zhǔn)確的。近震的定位要確定四個參數(shù)，即震源的空間坐標(biāo)（，h）及發(fā)震時刻t 。為此，我們需要知道臺網(wǎng)各個臺站的坐標(biāo)，以及至少五個臺的地震波到時。由臺網(wǎng)提供的到時數(shù)據(jù)主要是指P 波的到時，因為各種續(xù)至波的相位通常是難以識別的。如果地震發(fā)生在臺網(wǎng)內(nèi)，那么只利用P 波到時定位并不困難。如果地震發(fā)生在臺網(wǎng)以外，這種答案是不確定的。為了進行參數(shù)的求解，一般要對地殼平均速度做出假設(shè)。為了對

11、此問題進行研究，有必要建立近震定位的正反演模型。對于正演模型，其基本思路是通過設(shè)定震源初始參數(shù)（經(jīng)度、緯度、高程、發(fā)震時刻及均勻介質(zhì)中的波速）計算出地震波從震源到各臺站的走時，進而求出各臺站接收到地震波的時刻；這些數(shù)據(jù)就作為反演的數(shù)據(jù)來源。此模型的建立主要涉及兩個關(guān)鍵問題：坐標(biāo)的變換、走時計算。3.1基本原理3.1.1坐標(biāo)變換對于近震及地方震的計算，采用平面坐標(biāo)比球面坐標(biāo)方便。為此需要將經(jīng)緯度換算成局部的直角坐標(biāo)系。這里需要指出，為了處理數(shù)據(jù)方便，一般是將國內(nèi)基準(zhǔn)臺網(wǎng)及區(qū)域臺網(wǎng)的全部臺站的經(jīng)緯度及高程均輸入計算機內(nèi)。但是，在作具體計算時，只用到其中一部分臺站，也就是只需要把這部分臺站的經(jīng)緯度（

12、球面坐標(biāo)）化為平面坐標(biāo)。因而，直角坐標(biāo)系的原點是浮動的，每次都選來盡可能的靠近震中區(qū)，以減少變換過程中帶來的誤差。在計算程序中，坐標(biāo)原點確定在接收到該次地震信號各臺站經(jīng)緯度的平均值0 及 0 上， x 軸向東，y 軸向北，如圖3.1 所示。圖 3.1坐標(biāo)轉(zhuǎn)換示意圖于是，地面上任一點（其經(jīng)緯度為、）相對于坐標(biāo)原點O（ 0, 0 ）的直角坐標(biāo)、的換算公式為：（3.1 ）其中參數(shù) :而參數(shù)和中的常數(shù)又為：=6378.160 km(地球長半徑 )= 6356.778km(地球短半徑 )經(jīng)換算為平面坐標(biāo)系后，就可以用計算機確定震源參數(shù)。但這樣得到的震中位置仍是局部的、坐標(biāo)。因此，必須化為經(jīng)緯度，輸出。將

13、直角坐標(biāo)值，化為經(jīng)緯度，的換算公式為：(3.2)3.1.2正演計算如果已知震源及臺站坐標(biāo)；設(shè)震源坐標(biāo), 發(fā)生時刻，均勻介質(zhì)中波速；有n 個臺站，其中任一臺站坐標(biāo)，相應(yīng)直達波到時如圖3.2 所示。圖 3.2正演計算示意圖由物理知識可知：Ti( xix)2( yiy) 2z 2 / vT0，然而實際中，臺站及震中坐標(biāo)是以經(jīng)、緯度（地理坐標(biāo)）表示的，在計算前首先要把地理坐標(biāo)（經(jīng)、緯度）轉(zhuǎn)換成大地坐標(biāo)（直角坐標(biāo)）。3.2正演模型建立假定某地震事件發(fā)生后被8 個臺站檢測到 ( 圖 3.3)，這 8 個臺站的經(jīng)緯度分別是：（ 96.000000 ， 28.000000 ）、（ 95.000000 ， 29

14、.000000 ）、（ 97.000000 ， 29.000000 ）（ 94.000000 ， 30.000000 ）、（ 98.000000 ， 30.000000 ）、（ 95.000000 ， 31.000000 ）（ 97.000000 ， 31.000000 ）、（ 96.000000 ， 32.000000 ）通過給定兩個模型：Mod1的震中位置為：（96.000000 ， 30.000000 ）Mod2的震中位置為：（95.500000 ， 29.500000 ），見圖 3.3圖 3.3臺站分布圖正演的結(jié)果：Mod1：t= 44.386244 、29.397192 、29.39

15、7192 、38.647902 、38.647902 、29.52392329.523923 、44.386244Mod2：t= 34.676012 、14.846849 、31.128296 、31.261550 、49.893035 、34.715948 、44.367994 、56.320175對于 mod1而言，震中在其中心，它與（96.000000 ， 28.000000 ）臺站相距20，有換算關(guān)系10 111.19 ，且速度 v=5.0 可知此正演結(jié)果與實際情況相符合。Modd2也是如此。4 反演計算對于反演模型，主要目的就是利用地震波在各接收臺站的到時，結(jié)合臺站經(jīng)緯度，通過數(shù)學(xué)推

16、導(dǎo)而確定震源參數(shù)。實際計算表明，只要輸入數(shù)據(jù)正確，初定的震源位置及發(fā)震時刻足夠精確，就可供地震速報之用。其具體方法如下。在此主要利用前面正演模型所得數(shù)據(jù)，作為反演之用，進而推導(dǎo)出震源參數(shù)。4.1基本原理4.1.1波速 V 未知的初定為方法此方法通常用在臺數(shù)n>4，震中距小于170km，已知直達縱波P 到時，而且臺站方位和震中距分布合理的情況下。設(shè)接收到某次地震直達波到時的臺站有n 個（ n 4）。其中任一臺站的坐標(biāo)為（ xi ,y i ）相應(yīng)的直達波到時為Ti 。待求的震中位置為x,y, 震源深度z, 發(fā)震時刻為T（圖 4.1 ）；按照均勻地殼模型，可列出下列方程：圖 4.1震源計算機定

17、位原理圖( x xi ) 2( y yi ) 2z 2V 2 (T Ti ) 2（4.1 ）i=1,2, ， n.式中為地震波在地殼中的平均傳播速度，因區(qū)域不同而異。式（ 4.1 ）為非線性方程組，為便于直接求解，首先應(yīng)將其線性化，變?yōu)榫€性方程組。為此將（ 4.1 ）式展開得x 22xi x x2y 22 yi y yi2z2(4.2)V 2 (Ti22TTiT 2 )i=1,2,n.對于第一個臺，（i=1 ）可以寫出x 22x1 x x 2y22 y1 y y12z2(4.3)V 2(T122TT1 T2)將 (4.3)式中 i 分別代以2， 3，, ， n 與 (4.3)式相減得( xix

18、1 ) x ( yiy1 ) y (t i2t12 )V 2 / 2 (tit1 )TV 2( x2x2y 2y2 ) / 2（4.4 ）i1i1共有 n-1 個線性方程式。 (4.4)式以消去了深度參數(shù)z，只有 x 、y、T、V 四個未知數(shù)。由于 n-1 一般都大于4，故 (4.4) 式所表示的線性方程個數(shù)多于未知數(shù)個數(shù)，即構(gòu)成所謂的超定方程組(x2x1 )x ( y2y1 ) y (t 22t12 )V 2 / 2 (t2t1 )TV 2( x22x12y22y12 ) / 2(x3x1 ) x ( y3y1 ) y (t 32t12 )V 2 / 2(t3t1 )TV 2( x32x12

19、y32y12 ) / 2,( xn x1 ) x ( y ny1 ) y (t n2t12 )V 2 / 2(t nt1 )TV 2( xn2x12yn2y12 ) / 2（4.5 ）為了表達簡潔起見，將此超定方程組用矩陣形式表達(4.6)其中，對于n 個臺站的數(shù)據(jù) ,x2x1y2y1(t22t12 ) / 2t1t 2x3x1y3y1(t32t12 ) / 2t1t 3Axnx1yny1(tn2t12 ) / 2t1t n22T。其中 U=V， W=V當(dāng)?shù)卣鹋_站數(shù)目 n小于 5時，上述方程組為不定方程，有多解；大于5時，方程組為超定方程，通過變換（ A為 A的轉(zhuǎn)置）使得超定方程組變?yōu)檎ǚ匠?/p>

20、組，從而計算出震源參數(shù)，這就是最小二乘法原理。只要解出此方程組 (4.6) ，相應(yīng)震源參數(shù)也可求得。即：震中位置 :（x， y）介質(zhì)的波速 :V發(fā)震時刻 :T=W/U.又由：(x xi ) 2( y yi ) 2z2V 2 (T Ti ) 2任取一臺站 (i)數(shù)據(jù)代入此式便得震源深度：zV 2 (TTi ) 2( xxi ) 2( yyi ) 24.1.2當(dāng)波速 V 已知的初定為方法因為所研究區(qū)域位于青藏高原，其地殼的厚度較大，所以近震的深度基本在上地殼的范圍之內(nèi)，其速度大約為v=5.8km/s. 當(dāng)速度為已知值后，在上述方法中的線性方程組就會減少一個未知數(shù)，所以（3.5 ）變化之后，( xn

21、x1 )x ( yny1 ) y (t nt1 )TV 2（4.7 ）(x2x2y 2y 2 ) / 2 (t2t2)V 2/ 2n1n1n1其中x2x1y2y1t1t 2x3x1y3y1t1t 3Axnx1yny1t1t nx22x12y22y12t 22t121 x32x12y32y12t32t12B2xn2x12yn2y12t n2t122W=VT。只要解出此方程組(4.7) ，相應(yīng)震源參數(shù)也可求得。即：震中位置 :（x， y）發(fā)震時刻 : T=W/ V2.同理，再由zV 2 (TTi ) 2( xxi )2( yyi )2 求出震源的深度。這兩種方法的主要流程為：開始臺站經(jīng)緯度轉(zhuǎn)成直角

22、坐標(biāo)值計算方程組系數(shù)矩陣A和常數(shù)矩陣B求解方程組AX=B依 方 程 組 的 解 進 一 步 求 得 震 源 參 數(shù) (x,y,z,T)作出震源與臺站相對位置圖形直角坐標(biāo)值轉(zhuǎn)回經(jīng)緯度表示輸出結(jié)束圖 4.2震源計算機定位流程圖4.2模型反演結(jié)果為了更好的檢驗程序，在正演t 的結(jié)果上均加了6s，下面為各個反演方法所得到的結(jié)果：1. 當(dāng)介質(zhì)的速度 V 未知時，兩個模型反演出的地震參數(shù)分別為：2. 當(dāng)介質(zhì)的速度 V 作為已知條件（ V=5.0 ）時，兩個模型反演出的地震參數(shù)分別為：4.3誤差分析從數(shù)據(jù)的反演來看，深度的誤差很大，這主要是因為震源深度對地震波到時的讀取誤差反應(yīng)極其敏感，是地震學(xué)中最難準(zhǔn)確測

23、定的參數(shù)之一，各種方法對于震源深度的估計都具相當(dāng)程度的不確定性，影響著人們對震源過程的認(rèn)識。具體來說，震源深度的精度誤差受震中距、到時誤差和速度模型（地殼模型）三個因素制約，而且這些因素對震源深度的影響是非線性的。1. 當(dāng)?shù)卣鸩▊鞑ニ俣纫欢〞r，震源深度的誤差隨著震中距或臺站位置的增大和走時殘差的增大而增大。2. 走時殘差一定時，震源深度誤差隨著震中距的增大和地震波速度的增大而增大。3. 當(dāng)速度已知，走時殘差一定時，越淺的地震，定位誤差可能越大。這一結(jié)論與許多深源地震的定位結(jié)果相一致。盡管震源深度是地震學(xué)中很難精確測定的參數(shù)之一，各種方法測定震源深度的結(jié)果不同，但可以將不同的測定震源深度的結(jié)果進

24、行對比分析，可能會改善對震源深度的測定精度。地震發(fā)生后，立即在震源區(qū)布設(shè)流動觀測地震臺站（網(wǎng)）是修正主震震源深度精度的有效方法。5 實例剖析所研究區(qū)域位于青藏高原南迦巴瓦地段，該處臺站的分布如圖5.1 。5.1區(qū)域臺站及示例震源分布圖1臺網(wǎng)記錄的地震事件里該近震記錄根據(jù)美國國家高級地震系統(tǒng)（ANSS）提供的全球地震事件目錄，可知該事件的38 道記錄，震中距范圍在0.5 0-3 0。經(jīng)過上述方法反演后得到的位置為：圖分別是利用 geotool軟件在選定地震事件前后，通過屏幕抓圖得到的圖像。其中，紅色的p標(biāo)注的時刻是geotool軟件算出來的該地震記錄道的起跳時刻，是理論到時。黑色p 標(biāo)注處的時刻

25、是經(jīng)人工判斷后標(biāo)記的該地震道記錄的起跳時刻，是實際到時。圖 5.2 中， a 為定位前的地震記錄圖像， b 為定位后的圖像。由圖可看出反演后的結(jié)果。反演后理論到時和實際的起跳時刻仍然存在幾秒的誤差，因為模型畢竟是理想狀態(tài)下的情況，而實際的資料并不是那么完美，在實際操作中很難得到極其精確的臺站到時，一方面存在人為因素，如時間拾取的精度；另一方面也存在機器設(shè)備的物理因素。而且由于地表的覆蓋層厚度，松散度等等的不同等等都會引起波至到時的誤差，所以利用地震資料到時時，還要經(jīng)過判斷和人工調(diào)整。在該地震事件中，各個記錄道起跳刻的最小誤差為2s，最大的誤差為5s。a： 近震定位前b： 近震定位后圖5.2實例

26、一近震定位前（a）和定位后（b）2臺站網(wǎng)記錄的地震事件里沒有近震記錄時當(dāng)臺站網(wǎng)記錄的地震事件里沒有近震記錄時，利用上述反演方法進行地震的定位，可確定該地震事件發(fā)生的經(jīng)度和緯度為：事件就行處理后，結(jié)果如下圖5.3 所示，其中最小的誤差為0. 1s ，最大的誤差為1.3s 。a: 近震定位前b：近震定位后圖 5.3 實例二近震定位前（ a）和定位后（ b）為了進一步驗證此種定位的正確與否，再加一例輔證。同樣采用上述反演方法定位后，可確定該地震事件發(fā)生的位置為a: 近震定位前b：近震定位后圖5.4實例二近震定位前（a）和定位后（b）6 結(jié)束語通過本次畢業(yè)設(shè)計，使我了解了近震震源參數(shù)定位的方法和原理，

27、從正演和反演的理論知識落腳到實際資料的處理。解決此問題的基本思想就是將非線性問題化為線性問題，然后利用最小二乘法原理迭代，全主元消元法等求解，通過反復(fù)修定，以得到震源參數(shù)的最佳值。此次設(shè)計的目的在于借助震源的計算機初步定位這個問題，進一步了解一般科學(xué)計算研究所采用的正、反演計算方法及實例應(yīng)用。本次畢業(yè)設(shè)計是在王教授悉心指導(dǎo)下，經(jīng)過不斷的學(xué)習(xí)和修改完成的。王教授淵博的學(xué)識，豐富的實踐經(jīng)驗，高瞻遠(yuǎn)矚、敏銳的科學(xué)眼光，將是我永遠(yuǎn)學(xué)習(xí)的楷模；他樂觀的心態(tài)、謙遜的為人、樸實的生活態(tài)度，令我深深敬佩；他嚴(yán)謹(jǐn)細(xì)致、一絲不茍的作風(fēng)將是我以后學(xué)習(xí)、工作中的榜樣，他循循善誘的教導(dǎo)和不拘一格的思路給予我無盡的啟迪。

28、謹(jǐn)致以衷心的感謝和崇高的敬意。本次畢業(yè)論文的撰寫也得到眾多老師和同學(xué)的關(guān)心和幫助。在此一并表示衷心的感謝。祝愿他們身體健康，工作順利，事業(yè)上取得更大成功。由于本人水平有限，本篇論文及程序有待進一步改進，請各位老師批評指正。參考文獻1SethStein,MichaelWysession.AnIntroductiontoSeismologyEarthquakes,and Earth Structure.Blackwell Publishing, 2003.2 李善邦，中國地震 . 地震出版社 ,1981.3 單娜琳 , 程志平 , 劉云禎 . 工程地震勘探 . 冶金工業(yè)出版社， 2006.4 姜枚

29、 , 王有學(xué) , 錢輝 . 造山的高原青藏高原及其鄰區(qū)的寬頻地震探測與地殼上地幔結(jié)構(gòu) . 地質(zhì)出版社， 2009.5 熊章強 , 方根顯 . 淺層地震勘探 . 地震出版社， 2002.6 大港油田科技叢書編委會 . 地震勘探資料處理和解釋技術(shù) . 石油工業(yè)出版社， 1999.7 時振梁 , 張少泉 , 趙榮國等 . 地震工作手冊 . 地震出版社， 1990.8 傅淑芳 , 劉寶城 , 李文藝 . 地震學(xué)教程（上、下冊） . 地震出版社， 1980.9 彭國倫 .Fortran 95 程序設(shè)計 . 中國電力出版社， 2002.附錄程序program testimplicit noneintege

30、r : i, j, k, nstn, istnreal*8 : x0, y0, v0, t0, xx, yy, coef(1000,5)real*8 : delt(1000), solut(4,5), x(1000), y(1000), t(1000)real*8 : ph0, la0, phi, lai ,v,h, MinLat, MaxLat, MinLon, MaxLoncharacter : filename*100write(*,*) "input filename"read(*,'(a)') filename!read the picks for

31、 all stationsopen(1,file=filename)read(1,*) nstnx0=0;y0=0do istn=1,nstnread(1,*) x(istn),y(istn), t(istn)x(istn)=x(istn)*atan(1.00)/45y(istn)=y(istn)*atan(1.00)/45x0=x0+x(istn)y0=y0+y(istn)end dola0=x0/nstnph0=y0/nstnclose(1)goto 10!calculating sythetic datala0=99*atan(1.00)/45ph0=27*atan(1.00)/45op

32、en(2,file='mod13.dat')write(2,'(I2)') nstnv=5.0h=10.0do istn=1,nstnlai=x(istn)phi=y(istn)call lp2xy(lai,phi,la0,ph0,xx,yy)t(istn)=sqrt(xx*2.00+yy*2.00+h*2.00)/vwrite(2,'(6f15.6)') x(istn)*45.00/atan(1.00),y(istn)*45.00/atan(1.00), t(istn)+6end dowrite(2,'(4f20.6)') la

33、0*45.00/atan(1.00),ph0*45.00/atan(1.00),v,hclose(2)！ stop10 continue !(lai,phi)->(x,y) do istn=1,nstnlai=x(istn)phi=y(istn)call lp2xy(lai,phi,la0,ph0,xx,yy) x(istn)=xxy(istn)=yy end do!calculating coefficents for the linear equations do istn=2,nstncoef(istn,1)=x(istn)-x(1) coef(istn,2)=y(istn)-y(

34、1) coef(istn,3)=(t(istn)*2-t(1)*2)/2.00 coef(istn,4)=-t(istn)+t(1)delt(istn) =(x(istn)*2-x(1)*2+y(istn)*2-y(1)*2)/2.00 end do!least squar-root method to get the solution! dt(1)de(1,1)de(1,2) de(1,3) de(1,4) de(1,5)! dt(2)de(2,1)de(2,2) de(2,3) de(2,4) de(2,5)! .! .! .! dt(n)de(n,1)de(n,2) de(n,3) de

35、(n,4) de(n,5)! de(1,1) de(2,1) . de(n,1)! de(1,2) de(2,2) . de(n,2)! de(1,3) de(2,3) . de(n,3)! de(1,4) de(2,4) . de(n,4)! de(1,5) de(2,5) . de(n,5) do i=1,4do j=1,4 solut(i,j)=0 do k=1,nstnsolut(i,j)=solut(i,j)+coef(k,i)*coef(k,j)end doend doend dodo i=1,4solut(i,5)=0do k=1,nstnsolut(i,5)=solut(i,5

36、)+coef(k,i)*delt(k)end doend docall gaussian(solut,4,5)xx=solut(1,5)yy=solut(2,5)v0=sqrt(solut(3,5)t0=solut(4,5)/solut(3,5)h=0！ print *,xx,yy do istn=1,nstnh=h+sqrt(v0*v0*(t(7)-t0)*(t(7)-t0)-(x(7)-xx)*2-(y(7)-yy)*2) end doh=h/nstncall xy2lp(xx,yy,la0,ph0,lai,phi)print *, lai*45.00/atan(1.00),phi*45.00/atan(1.00),v0,t0, hend program testSUBROUTINE Gaussian(A,N1,M1)implicit noneinteger :n1,m1,k1,k2,ik,i,j ,l, i1REAL*8 : A(n1,m1)REAL*8 : BMAX,T,EPSEPS=0.DO K1=1,N1BMAX=0.DO I=K1,N1IF(BMAX-ABS(A(I,K1).LT.0) THENBMAX=ABS(A(I,K1)L=IEND IFEND DOIF(BMAX.LT.EPS) STOP 4444IF(L.NE.K1) THENDO J