近震定位和應(yīng)用本科畢業(yè)論文

1、. . . . 目 次摘要Abstract1 前言 12 近震的理論基礎(chǔ)22.1 地震學(xué)的基本名詞和概念2 2.2 地震的分類(lèi)22.3 近震的主要震相3 3 正演計(jì)算（模型試算） 53.1 基本原理53.1.1 坐標(biāo)變換53.1.2 正演計(jì)算63.2 正演模型的建立74 反演計(jì)算84.1 基本原理84.1.1 當(dāng)速度V未知的初定方法84.1.2 當(dāng)速度V已知的初定方法104.2 模型反演結(jié)果114.3 誤差分析115 實(shí)例剖析136 結(jié)束語(yǔ)17參考文獻(xiàn) 18附錄 近震定位程序1924 / 251 前 言汶川地震和地震發(fā)生時(shí)的地動(dòng)山搖、房倒屋塌、生死分離的那種凄涼、悲慘的場(chǎng)面，又一次的讓我們感受

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

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

4、法等十幾種方法。計(jì)算機(jī)的近震定位方法主要分為速度未知的初定為方法和速度已知的初定為方法。目前定位精度已達(dá)到較高的水平。如果有合理的臺(tái)網(wǎng)分布和適用的走時(shí)表，對(duì)于<100km的區(qū)域性事件，定位精度為1km左右；對(duì)于100-1000km的近震，可達(dá)2至5km；遠(yuǎn)震為510km。一般淺源地震（h<100km）的震源深度誤差為深度值的10%左右。震源愈深，相對(duì)誤差愈小；h>300km時(shí)，誤差小于其深度的5%。發(fā)震時(shí)刻的誤差為1/101/2s.2 近震的理論基礎(chǔ)地震學(xué)主要研究地震的發(fā)生、地震波的傳播與地球部構(gòu)造，地震波能穿透到地球部，并能把地球部的信息帶回地面，具體說(shuō)來(lái)，它主要是根據(jù)天然

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

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

7、為研究方便，按震動(dòng)的性質(zhì)，可分為天然地震、人工地震與脈動(dòng)三類(lèi)。對(duì)于天然地震，有下述分類(lèi)：（一）按震源深度分：1淺源地震：震源深度小于60公里的天然地震稱(chēng)為淺震；也稱(chēng)正常深度地震。2中源地震：震源深度在60公里至300公里之間的地震稱(chēng)為中源地震。3深源地震：震源深度大于300公里的地震稱(chēng)為深震。已記錄到的最深地震的震源深度約700公里。有時(shí)也將中源地震和深源地震統(tǒng)稱(chē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 近震震相與其射線(xiàn)路徑示意圖 （c、c分別為康拉德界面與莫霍面的臨界角）1. Pg/Sg震相由于地球淺部存在強(qiáng)烈的橫向與縱向上的速度變化，尤其是有沉積該層地區(qū)的垂向梯度的影響，地震波射線(xiàn)在地殼淺部10公里的圍發(fā)生廻折而到達(dá)地表，并在震中附近數(shù)公里的圍，這種波束與直達(dá)波，被稱(chēng)之為廻折或潛波（Diving Wave），對(duì)于縱、橫波來(lái)說(shuō)，分別命名為Pg和Sg。2. P*/S*震相P*/S*分別為地震波在地殼部的康拉德界面產(chǎn)生的折射縱、橫波，它們的速度

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

10、震時(shí)刻最常用的方法是將非線(xiàn)性問(wèn)題化為線(xiàn)性問(wèn)題，然后利用最小二乘法原理迭代，全主元消元法等求解。通過(guò)反復(fù)修定，以得到震源參數(shù)的最佳值。此方法對(duì)近震經(jīng)緯度的確定比較準(zhǔn)確的。近震的定位要確定四個(gè)參數(shù)，即震源的空間坐標(biāo)（，h）與發(fā)震時(shí)刻t。為此，我們需要知道臺(tái)網(wǎng)各個(gè)臺(tái)站的坐標(biāo)，以與至少五個(gè)臺(tái)的地震波到時(shí)。由臺(tái)網(wǎng)提供的到時(shí)數(shù)據(jù)主要是指P 波的到時(shí)，因?yàn)楦鞣N續(xù)至波的相位通常是難以識(shí)別的。如果地震發(fā)生在臺(tái)網(wǎng)，那么只利用P波到時(shí)定位并不困難。如果地震發(fā)生在臺(tái)網(wǎng)以外，這種答案是不確定的。為了進(jìn)行參數(shù)的求解，一般要對(duì)地殼平均速度做出假設(shè)。為了對(duì)此問(wèn)題進(jìn)行研究，有必要建立近震定位的正反演模型。對(duì)于正演模型，其基本思

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

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

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

14、00000， 30.000000）、（95.000000， 31.000000）（97.000000， 31.000000）、（96.000000， 32.000000）通過(guò)給定兩個(gè)模型：Mod1的震中位置為：（96.000000， 30.000000）Mod2的震中位置為：（95.500000， 29.500000），見(jiàn)圖3.3圖3.3 臺(tái)站分布圖正演的結(jié)果：Mod1：t=44.386244、29.397192、29.397192、38.647902、38.647902、29.52392329.523923、44.386244Mod2：t= 34.676012、14.846849、31.12

15、8296、31.261550、49.893035、34.715948、44.367994、56.320175對(duì)于mod1而言，震中在其中心，它與（96.000000，28.000000）臺(tái)站相距20，有換算關(guān)系10 111.19，且速度v=5.0可知此正演結(jié)果與實(shí)際情況相符合。Modd2也是如此。4 反演計(jì)算對(duì)于反演模型，主要目的就是利用地震波在各接收臺(tái)站的到時(shí)，結(jié)合臺(tái)站經(jīng)緯度，通過(guò)數(shù)學(xué)推導(dǎo)而確定震源參數(shù)。實(shí)際計(jì)算表明，只要輸入數(shù)據(jù)正確，初定的震源位置與發(fā)震時(shí)刻足夠精確，就可供地震速報(bào)之用。其具體方法如下。在此主要利用前面正演模型所得數(shù)據(jù)，作為反演之用，進(jìn)而推導(dǎo)出震源參數(shù)。4.1 基本原理4.

16、1.1 波速V未知的初定為方法此方法通常用在臺(tái)數(shù)n>4，震中距小于170km，已知直達(dá)縱波P到時(shí)，而且臺(tái)站方位和震中距分布合理的情況下。設(shè)接收到某次地震直達(dá)波到時(shí)的臺(tái)站有n個(gè)（n4）。其中任一臺(tái)站的坐標(biāo)為（xi,yi）相應(yīng)的直達(dá)波到時(shí)為T(mén)i。待求的震中位置為x,y,震源深度z,發(fā)震時(shí)刻為T(mén)（圖4.1）；按照均勻地殼模型，可列出下列方程： 圖4.1 震源計(jì)算機(jī)定位原理圖 （4.1） i=1,2,，n.式中為地震波在地殼中的平均傳播速度，因區(qū)域不同而異。式（4.1）為非線(xiàn)性方程組，為便于直接求解，首先應(yīng)將其線(xiàn)性化，變?yōu)榫€(xiàn)性方程組。為此將（4.1）式展開(kāi)得 (4.2) i=1,2, ,n.對(duì)于

17、第一個(gè)臺(tái)，（i=1）可以寫(xiě)出 (4.3)將(4.3)式中i分別代以2，3，n與(4.3)式相減得 （4.4）共有n-1個(gè)線(xiàn)性方程式。(4.4)式以消去了深度參數(shù)z，只有x、y、T、V四個(gè)未知數(shù)。由于n-1一般都大于4，故(4.4)式所表示的線(xiàn)性方程個(gè)數(shù)多于未知數(shù)個(gè)數(shù)，即構(gòu)成所謂的超定方程組 （4.5）為了表達(dá)簡(jiǎn)潔起見(jiàn)，將此超定方程組用矩陣形式表達(dá) (4.6)其中，對(duì)于n個(gè)臺(tái)站的數(shù)據(jù),其中U=V2， W=V2T。當(dāng)?shù)卣鹋_(tái)站數(shù)目n小于5時(shí)，上述方程組為不定方程，有多解；大于5時(shí)，方程組為超定方程，通過(guò)變換（A為A的轉(zhuǎn)置）使得超定方程組變?yōu)檎ǚ匠探M，從而計(jì)算出震源參數(shù)，這就是最小二乘法原理。只要解

18、出此方程組(4.6)，相應(yīng)震源參數(shù)也可求得。即：震中位置: （x，y）介質(zhì)的波速: V發(fā)震時(shí)刻: T=W/U.又由： 任取一臺(tái)站(i)數(shù)據(jù)代入此式便得震源深度：4.1.2 當(dāng)波速V已知的初定為方法因?yàn)樗芯繀^(qū)域位于青藏高原，其地殼的厚度較大，所以近震的深度基本在上地殼的圍之，其速度大約為v=5.8km/s.當(dāng)速度為已知值后，在上述方法中的線(xiàn)性方程組就會(huì)減少一個(gè)未知數(shù)，所以（3.5）變化之后 ， （4.7）其中 W=V2T。只要解出此方程組(4.7)，相應(yīng)震源參數(shù)也可求得。即：震中位置: （x，y）發(fā)震時(shí)刻: T=W/ V2.同理，再由求出震源的深度。這兩種方法的主要流程為： 圖4.2 震源計(jì)算

19、機(jī)定位流程圖4.2 模型反演結(jié)果 為了更好的檢驗(yàn)程序，在正演t的結(jié)果上均加了6s，下面為各個(gè)反演方法所得到的結(jié)果：1.當(dāng)介質(zhì)的速度V未知時(shí)，兩個(gè)模型反演出的地震參數(shù)分別為： mod1:震中(96.000000,29.999999) t0=6.000025 v=5.000002 h=9.998836 mod2:震中(95.490894,29.495549) t0=5.797093 v=4.994111 h=20.3332112.當(dāng)介質(zhì)的速度V作為已知條件（V=5.0）時(shí)，兩個(gè)模型反演出的地震參數(shù)分別為： mod1:震中(96.000000,29.999999) t0=6.000000 h=9.9

20、99972 mod2:震中(95.490701,29.495359) t0=5.874814 h=18.7946024.3 誤差分析從數(shù)據(jù)的反演來(lái)看，深度的誤差很大，這主要是因?yàn)檎鹪瓷疃葘?duì)地震波到時(shí)的讀取誤差反應(yīng)極其敏感，是地震學(xué)中最難準(zhǔn)確測(cè)定的參數(shù)之一，各種方法對(duì)于震源深度的估計(jì)都具相當(dāng)程度的不確定性，影響著人們對(duì)震源過(guò)程的認(rèn)識(shí)。具體來(lái)說(shuō)，震源深度的精度誤差受震中距、到時(shí)誤差和速度模型（地殼模型）三個(gè)因素制約，而且這些因素對(duì)震源深度的影響是非線(xiàn)性的。1.當(dāng)?shù)卣鸩▊鞑ニ俣纫欢〞r(shí)，震源深度的誤差隨著震中距或臺(tái)站位置的增大和走時(shí)殘差的增大而增大。2.走時(shí)殘差一定時(shí)，震源深度誤差隨著震中距的增大和地

21、震波速 度的增大而增大。3.當(dāng)速度已知，走時(shí)殘差一定時(shí)，越淺的地震，定位誤差可能越大。這一結(jié)論與許多深源地震的定位結(jié)果相一致。盡管震源深度是地震學(xué)中很難精確測(cè)定的參數(shù)之一，各種方法測(cè)定震源深度的結(jié)果不同，但可以將不同的測(cè)定震源深度的結(jié)果進(jìn)行對(duì)比分析，可能會(huì)改善對(duì)震源深度的測(cè)定精度。地震發(fā)生后，立即在震源區(qū)布設(shè)流動(dòng)觀測(cè)地震臺(tái)站（網(wǎng)）是修正主震震源深度精度的有效方法。 5 實(shí)例剖析所研究區(qū)域位于青藏高原南迦巴瓦地段，該處臺(tái)站的分布如圖5.1。5.1 區(qū)域臺(tái)站與示例震源分布圖1臺(tái)網(wǎng)記錄的地震事件里該近震記錄根據(jù)美國(guó)國(guó)家高級(jí)地震系統(tǒng)（ANSS）提供的全球地震事件目錄，可知該事件的位置（96.6050,

22、29.5730），發(fā)震時(shí)刻為2003年第230天的09:03,該地震事件總有38道記錄，震中距圍在0.50-30。經(jīng)過(guò)上述方法反演后得到的位置為：（95.6102,29.4248），反演的結(jié)果與實(shí)際經(jīng)緯度分別相差（0.0052,0.1382）。下圖分別是利用geotool軟件在選定地震事件前后，通過(guò)屏幕抓圖得到的圖像。其中，紅色的p標(biāo)注的時(shí)刻是geotool軟件算出來(lái)的該地震記錄道的起跳時(shí)刻，是理論到時(shí)。黑色p標(biāo)注處的時(shí)刻是經(jīng)人工判斷后標(biāo)記的該地震道記錄的起跳時(shí)刻，是實(shí)際到時(shí)。圖5.2中，a為定位前的地震記錄圖像，b為定位后的圖像。由圖可看出反演后的結(jié)果。反演后理論到時(shí)和實(shí)際的起跳時(shí)刻仍然存在

23、幾秒的誤差，因?yàn)槟Ｐ彤吘故抢硐霠顟B(tài)下的情況，而實(shí)際的資料并不是那么完美，在實(shí)際操作中很難得到極其精確的臺(tái)站到時(shí)，一方面存在人為因素，如時(shí)間拾取的精度；另一方面也存在機(jī)器設(shè)備的物理因素。而且由于地表的覆蓋層厚度，松散度等等的不同等等都會(huì)引起波至到時(shí)的誤差，所以利用地震資料到時(shí)時(shí)，還要經(jīng)過(guò)判斷和人工調(diào)整。在該地震事件中，各個(gè)記錄道起跳刻的最小誤差為2s，最大的誤差為5s。a：近震定位前 b： 近震定位后圖 5.2 實(shí)例一近震定位前（a）和定位后（b）2臺(tái)站網(wǎng)記錄的地震事件里沒(méi)有近震記錄時(shí)當(dāng)臺(tái)站網(wǎng)記錄的地震事件里沒(méi)有近震記錄時(shí)，利用上述反演方法進(jìn)行地震的定位，可確定該地震事件發(fā)生的經(jīng)度和緯度為：（9

24、6.7124,25.7474）。將得出其經(jīng)度、緯度、發(fā)震時(shí)刻、深度等信息，按照標(biāo)準(zhǔn)格式寫(xiě)入origin事件，利用geotool軟件對(duì)該事件就行處理后，結(jié)果如下圖5.3所示，其中最小的誤差為0. 1s，最大的誤差為1.3s。a:近震定位前 b：近震定位后圖5.3 實(shí)例二近震定位前（a）和定位后（b）為了進(jìn)一步驗(yàn)證此種定位的正確與否，再加一例輔證。同樣采用上述反演方法定位后，可確定該地震事件發(fā)生的位置為（96.7124,25.7474）。利用geotool軟件的p波波至標(biāo)定后的結(jié)果如下圖5.4，其最小的誤差為0. 09s，最大的誤差為3s。a:近震定位前 b：近震定位后圖5.4 實(shí)例二近震定位前（

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

26、的作風(fēng)將是我以后學(xué)習(xí)、工作中的榜樣，他循循善誘的教導(dǎo)和不拘一格的思路給予我無(wú)盡的啟迪。謹(jǐn)致以衷心的感和崇高的敬意。 本次畢業(yè)論文的撰寫(xiě)也得到眾多老師和同學(xué)的關(guān)心和幫助。在此一并表示衷心的感。祝愿他們身體健康，工作順利，事業(yè)上取得更大成功。由于本人水平有限，本篇論文與程序有待進(jìn)一步改進(jìn)，請(qǐng)各位老師批評(píng)指正。參考文獻(xiàn)1 Seth Stein,Michael Wysession. An Introduction to Seismology Earthquakes,and Earth Structure.Blackwell Publishing, 2003.2 善邦，中國(guó)地震. 地震,1981.3 單

27、娜琳,程志平,云禎.工程地震勘探. 冶金工業(yè)，2006.4 枚,王有學(xué),錢(qián)輝.造山的高原青藏高原與其鄰區(qū)的寬頻地震探測(cè)與地殼上地幔結(jié)構(gòu).地質(zhì)，2009.5 熊章強(qiáng),方根顯.淺層地震勘探.地震，2002.6 大港油田科技叢書(shū)編委會(huì).地震勘探資料處理和解釋技術(shù).石油工業(yè)，1999.7 時(shí)振梁,少泉,榮國(guó)等.地震工作手冊(cè).地震，1990.8 傅淑芳,寶城,文藝.地震學(xué)教程（上、下冊(cè)）. 地震，1980.9 國(guó)倫.Fortran 95 程序設(shè)計(jì).中國(guó)電力， 2002.附錄 程序program test implicit noneinteger : i, j, k, nstn, istn real*8

28、: 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 all stationsopen(1,file=file

29、name)read(1,*) nstnx0=0;y0=0do istn=1,nstn read(1,*) x(istn),y(istn), t(istn) x(istn)=x(istn)*atan(1.00)/45 y(istn)=y(istn)*atan(1.00)/45 x0=x0+x(istn) y0=y0+y(istn)end dola0=x0/nstn ph0=y0/nstnclose(1) goto 10!calculating sythetic data la0=99*atan(1.00)/45ph0=27*atan(1.00)/45 open(2,file='mod13

30、.dat')write(2,'(I2)') nstnv=5.0 h=10.0do istn=1,nstn lai=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)/v write(2,'(6f15.6)') x(istn)*45.00/atan(1.00), y(istn)*45.00/atan(1.00), t(istn)+6end dowrite(2,'(4f20.6)') la0*45.00/atan(

31、1.00),ph0*45.00/atan(1.00),v,hclose(2)！stop10 continue!(lai,phi)->(x,y)do istn=1,nstn lai=x(istn) phi=y(istn) call lp2xy(lai,phi,la0,ph0,xx,yy) x(istn)=xx y(istn)=yyend do!calculating coefficents for the linear equationsdo istn=2,nstn coef(istn,1)=x(istn)-x(1) coef(istn,2)=y(istn)-y(1) coef(istn,

32、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.00end 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(n,4) d

33、e(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,4 do j=1,4 solut(i,j)=0 do k=1,nstn solut(i,j)=solut(i,j)+coef(k,i)*coef(k,j) end do end doend dodo i=1,4 solut(i,5)=0 do k=1,nstn solut(i,5)=solut(i,5)

34、+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 do h=h/nstncall xy2lp(xx,yy,la0,ph0,lai,phi) print *,'-lai phi v0 t0 h-'

35、; print *, lai*45.00/atan(1.00),phi*45.00/atan(1.00),v0,t0, hend program testSUBROUTINE Gaussian(A,N1,M1) implicit none integer :n1,m1,k1,k2,ik,i,j ,l, i1 REAL*8 : A(n1,m1) REAL*8 : BMAX,T,EPS EPS=0. DO K1=1,N1 BMAX=0. DO I=K1,N1 IF(BMAX-ABS(A(I,K1).LT.0) THEN BMAX=ABS(A(I,K1) L=I END IF END DO IF(BMAX.LT.EPS) STOP 4444 IF(L.NE.K1) THEN DO J=K1,M1 T=A(L,J) A(L,J)=A(K1,J) A(K1,J)=T END DO END IF T=1./A(K1,K1)K2=K1+1 DO J=K2,M1 A(K1,J)=A(K1,J)*T DO I=K2,N1 A(I,J)=A