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

1、目錄摘要 摘要 1 前言 12 近震理論基礎(chǔ)22.1 地震學(xué)基本術(shù)語(yǔ)和概念 22.2 地震分類 . 22.3 近期地震的主要階段33 正演計(jì)算（模型試算） 53.1 基本原理 . 53.1.1坐標(biāo)變換53.1.2正向計(jì)算63.2 正演模型的建立74 反演計(jì)算. 84.1 基本原理 . 84.1.1速度 V 未知時(shí)的初始確定方法 84.1.2當(dāng)速度V已知時(shí)，初始確定方法104.2 模型反演結(jié)果 . 114.3 誤差分析115 案例研究. 136 結(jié)束語(yǔ)17參考文獻(xiàn). 18附錄 近震定位程序. 191 簡(jiǎn)介汶川大地震，大地震動(dòng)、房屋倒塌、生死分離的凄涼慘烈場(chǎng)面，再次讓我們感受到了自然災(zāi)害的威力，也

2、讓災(zāi)區(qū)人民心驚肉跳。安，甚至說(shuō)起“地震”變色，這是地震給人最直觀的印象。強(qiáng)地震的發(fā)生往往伴隨著地表變形和地層錯(cuò)位，其破壞力相當(dāng)大。主要表現(xiàn)形式有：大型建筑物毀壞、普通房屋毀壞、山體滑坡和地面裂縫、人畜傷亡等。這些天地震發(fā)生的頻率越來(lái)越高，2010年也被“尊重”為國(guó)際地震年。一切都有兩個(gè)方面。地震對(duì)我們的生命安全構(gòu)成極大威脅，為專家學(xué)者研究地球提供了豐富的資料。目前，人類對(duì)地球結(jié)構(gòu)的認(rèn)識(shí)主要來(lái)自于對(duì)地震學(xué)的研究，因此地震學(xué)已成為地球科學(xué)的一門重要學(xué)科。由于天然地震能量大，其產(chǎn)生的地震波可以穿透大深度，行進(jìn)距離長(zhǎng)，在地球深處遇到波阻抗差界面時(shí)，會(huì)產(chǎn)生反射波。地球近震資料為研究地球深部結(jié)構(gòu)提供了重要

3、的信息資源和途徑。使用地震記錄的定位起源于歐洲和日本。最初使用方位角方法，然后是幾何映射和地球投影。國(guó)際地震收集中心（ISS）是第一個(gè)使用最小二乘法計(jì)算修正震中的。 1961 年，博爾特和威爾莫爾合作改進(jìn)了計(jì)算方案，該方案首先在國(guó)際空間站使用，后來(lái)被國(guó)際地震中心和美國(guó)海岸與測(cè)量局使用。在我國(guó)，山幫先生最初使用方位角和最小二乘法來(lái)觀察和定位地震。 1953年，我國(guó)開(kāi)始利用大量觀測(cè)資料對(duì)震中進(jìn)行修正。目前我國(guó)對(duì)地震的定位方法既有測(cè)繪方法，也有計(jì)算方法。其中，局部地震的定位方法有石川法、和田法、高橋法、外中心方位角法、推定地震發(fā)生時(shí)間定位、等時(shí)板塊法等十幾種方法?；谟?jì)算機(jī)的局地地震定位方法主要分為

4、速度未知的初始確定方法和速度已知的初始確定方法。目前，定位精度已達(dá)到較高水平。如果有合理的臺(tái)網(wǎng)分布和合適的走時(shí)表，對(duì)于100km的區(qū)域性事件，定位精度1km在25km左右；近 100 級(jí) 的 地震 ,-1000km可達(dá)2 到 5公里;10km一般來(lái)說(shuō)，淺層地震的震源深度誤差（h 時(shí)，誤差小于其深度的 5%。地震時(shí)間誤差為1 /10-1/2s。2 近期地震的理論基礎(chǔ)地震學(xué)主要研究地震的發(fā)生、地震波的傳播和地球的結(jié)構(gòu)。地震波可以穿透地球，將地球的信息帶回地面。具體來(lái)說(shuō)，主要是基于自然地震或人工地震的數(shù)據(jù)。 ，利用物理、數(shù)學(xué)和地質(zhì)學(xué)的知識(shí)，研究地震的發(fā)生情況和地震波的傳播規(guī)律，以進(jìn)一步達(dá)到地震研究和

5、地震預(yù)測(cè)的目的。同時(shí)，地震是地球表面的振動(dòng)，是地殼構(gòu)造運(yùn)動(dòng)的一種形式，因此地震波的傳播特性也被用來(lái)研究地殼和地球各部分的結(jié)構(gòu)。圖2.1 地震波的分類由于整個(gè)彈性空間的存在，在彈性介質(zhì)中傳播的縱波和橫波統(tǒng)稱為體波。與體波相比，在彈性界面附近還有另一種波，它們僅在能量上分布在彈性界面附近，稱為面波。面波又分為瑞利波和洛夫波（如右圖圖2.1 地震波的分類2.1 地震學(xué)的基本術(shù)語(yǔ)和概念1、震中：地球上發(fā)生地震的地方稱為震中（或震中）。理論上，震源被視為一個(gè)點(diǎn)，但實(shí)際上它是一個(gè)區(qū)域。2、焦深：將震中看成一個(gè)點(diǎn)，該點(diǎn)到地面的垂直距離稱為焦深，一般用字母h表示。3、震中：震中在地面上的投影點(diǎn)稱為震中（或震中

6、區(qū)域）。同時(shí)，地面上受災(zāi)最嚴(yán)重的區(qū)域稱為震中區(qū)。理論上，震中區(qū)域與震中區(qū)域相同。事實(shí)上，由于地表當(dāng)?shù)貝毫訔l件的影響，震中區(qū)不一定是震中區(qū)。與震中相對(duì)的地球直徑的另一端稱為震中；或震中到點(diǎn)。4、震中距：在地面上，沿大圓弧從震中到任意一點(diǎn)的距離稱為震中距。通常由字母表示?？梢杂镁€距或地心角來(lái)表示。5、地震矩：地震發(fā)生的時(shí)刻，一般用字母O或T 0表示。標(biāo)有時(shí)間。它比格林威治標(biāo)準(zhǔn)時(shí)間早 8 小時(shí)。2.2 地震分類為方便研究，根據(jù)振動(dòng)的性質(zhì)，可分為自然地震、人工地震和脈動(dòng)三類。對(duì)于自然地震，有以下分類：(1)按焦深：1、淺層地震：震源深度小于60公里的自然地震稱為淺層地震；也稱為正常深度地震。2.中度地

7、震：震源深度在60公里到300公里之間的地震稱為中度地震。3、深地震：震源深度大于300公里的地震稱為深地震。有記錄以來(lái)最深的地震焦深約700公里（2）按震中距分： 1、局地地震：震中距小于100公里的地震。2、近震：震中在1000公里以內(nèi)的地震。3、遠(yuǎn)震：震中大于1000公里的地震。2.3 近期地震的主要階段由于近期地震震中距小，地震波傳播主要局限于地殼和莫霍面附近，主要震相為Pg/Sg( )、P*/S*、Pn/Sn , P*P/S*S , P*S/S*P, PmP/SmS 和 PmS/SmP 等，如圖 2.2 所示。圖 2.2 局地震相及其射線路徑示意圖( c , c分別是康拉德界面和莫霍

8、面的臨界角)1. Pg/Sg震相由于地球淺部強(qiáng)烈的橫向和垂直速度變化，特別是沉積區(qū)垂直梯度的影響，地震波射線10公里在地殼淺部折疊并到達(dá)地表，到達(dá)地表震中附近數(shù)公里的地表。這種波束和直達(dá)波稱為彎曲波或潛水波（Diving Wave）。對(duì)于縱波和橫波，它們分別命名為 Pg 和 Sg2. P*/S*震相P*/S*是地震波在地殼Conrad界面產(chǎn)生的折射縱波和橫波，它們的速度分別為6.3 7.1km/s和3.6 3.9km/s。需要注意的是，由于康拉德接口不是通用接口，因此在某些地區(qū)無(wú)法觀察到。3. Pn/Sn震相Pn/Sn是地震波在地殼底部莫霍面產(chǎn)生的折射縱波和橫波，它們的速度分別為7.9 8.2

9、km/s和4.4 4.8km/s。一般來(lái)說(shuō)，Pn相在莫霍面上的臨界反射在遠(yuǎn)區(qū)總是以初至波的形式出現(xiàn)。像其他折射波一樣，波的到達(dá)是線性分布的。4. P*P/S*S 和 P*S/S*P 階段P*P/S*S 分別是地震波在地殼康拉德界面產(chǎn)生的反射縱波和橫波，P*S/S*P 是入射到地殼后轉(zhuǎn)換的地震縱波/橫波分別為地殼的康拉德界面。產(chǎn)生的反射橫波/縱波。3 正演計(jì)算（模型試算）確定震源位置和地震發(fā)生時(shí)刻最常用的方法是將非線性問(wèn)題轉(zhuǎn)化為線性問(wèn)題，然后利用最小二乘原理迭代求解，采用全樞軸法。通過(guò)反復(fù)修改，可以得到源參數(shù)的最優(yōu)值。這種方法對(duì)于確定本次地震的經(jīng)緯度較為準(zhǔn)確。局部地震定位需要確定四個(gè)參數(shù)，即震中

10、空間坐標(biāo)（ , , h ）和地震發(fā)生時(shí)間 t。為此，我們需要知道網(wǎng)絡(luò)中每個(gè)臺(tái)站的坐標(biāo)，以對(duì)應(yīng)來(lái)自至少五個(gè)臺(tái)站的地震波。臺(tái)網(wǎng)提供的到站數(shù)據(jù)主要是指P波的到站，因?yàn)楦鞣N連續(xù)到站波的相位通常難以識(shí)別。如果地震發(fā)生在網(wǎng)絡(luò)中，僅使用P 波到達(dá)不難定位。如果地震發(fā)生在網(wǎng)絡(luò)之外，這個(gè)答案是不確定的。為了求解參數(shù)，通常需要對(duì)地殼的平均速度做出假設(shè)。為了研究這個(gè)問(wèn)題，有必要建立局部地震位置的正反演模型。對(duì)于正向模型，其其基本思想是通過(guò)設(shè)置震源的初始參數(shù)（經(jīng)度、緯度、高程、地震發(fā)生時(shí)間和波速在均質(zhì)介質(zhì)中），計(jì)算地震波從震源到各臺(tái)站的傳播時(shí)間，然后找出每個(gè)臺(tái)站接收地震波的時(shí)間；這些數(shù)據(jù)作為反演的數(shù)據(jù)源。該模型的建立

11、主要涉及兩個(gè)關(guān)鍵問(wèn)題：坐標(biāo)變換和走時(shí)計(jì)算。3.1 基本原則3.1.1坐標(biāo)變換對(duì)于局地地震和局地地震的計(jì)算，使用平面坐標(biāo)比球坐標(biāo)更方便。為此，需要將緯度和經(jīng)度轉(zhuǎn)換為本地笛卡爾坐標(biāo)系。這里需要指出的是，為方便數(shù)據(jù)處理，一般將國(guó)家基準(zhǔn)網(wǎng)和區(qū)域網(wǎng)中各臺(tái)站的經(jīng)度、緯度和高程輸入計(jì)算機(jī)。但在具體計(jì)算中，只用到了部分臺(tái)站，即這些臺(tái)站的經(jīng)緯度（球坐標(biāo)）只需轉(zhuǎn)換成平面坐標(biāo)即可。因此，笛卡爾坐標(biāo)系的原點(diǎn)是浮動(dòng)的，每次都選擇盡可能靠近震中，以減少變換過(guò)程帶來(lái)的誤差。在計(jì)算程序中，坐標(biāo)原點(diǎn)由接收到地震信號(hào)的臺(tái)站經(jīng)緯度0和0的平均值確定，x軸為東，y軸為北，如圖所示圖 3.1。圖 3.1 坐標(biāo)變換示意圖，地面上任意一點(diǎn)

12、（其經(jīng)緯度為, ）的直角坐標(biāo)相對(duì)于坐標(biāo)原點(diǎn)O ( 0 , 0 )的轉(zhuǎn)換公式為：(3.1)其中參數(shù)：參數(shù)sum中的常數(shù)為：= =0.0066935165（地球偏心率）= 6378.160 km（地球長(zhǎng)半徑）= 6356.778km（短地球半徑）轉(zhuǎn)換成平面坐標(biāo)系后，可以通過(guò)計(jì)算機(jī)確定源參數(shù)。但是，這種方式得到的震中位置仍然是局部的、坐標(biāo)的。因此，必須將其轉(zhuǎn)換為經(jīng)緯度，并輸出。將笛卡爾坐標(biāo)轉(zhuǎn)換為經(jīng)緯度，轉(zhuǎn)換公式為：(3.2)3.1.2前向計(jì)算如果源和站點(diǎn)的坐標(biāo)是已知的；設(shè)置了震源坐標(biāo)、發(fā)生時(shí)間、均質(zhì)介質(zhì)中的波速；有n個(gè)站點(diǎn)，其中任意一個(gè)站點(diǎn)的坐標(biāo)，對(duì)應(yīng)的直達(dá)波到達(dá)時(shí)間如圖3.2所示。圖 3.2 正向

13、計(jì)算示意圖從物理知識(shí)上可以知道：但在實(shí)踐中，臺(tái)站和震中的坐標(biāo)是用經(jīng)度和緯度（地理坐標(biāo)）來(lái)表示的。在計(jì)算之前，應(yīng)將地理坐標(biāo)（經(jīng)緯度）轉(zhuǎn)換為坐標(biāo)（直角坐標(biāo)）。3.2 正演模型的建立假設(shè)有 8 個(gè)臺(tái)站檢測(cè)到一次地震事件（圖 3.3），這 8 個(gè)臺(tái)站的經(jīng)緯度分別為：（ 96.000000，28.00000000 ） ， （ 95.000000，29.000000 ） ， （ 9 7.000000，29.000000 ） （ 94.000000，30.000000 ） ， （ 9 8.000000，30.000000 ）( 97.000000 , 31.000000), (96.000000 , 32.

14、000000 ) _ _通過(guò)給出兩個(gè)模型：Mod1 的震中是： ( 9 6 .000000 , 30 .000000 )Mod2的震中為：( 9 5.5 00000 , 29. 5 00000 )，見(jiàn)圖3.3圖 3.3 車站分布圖正在進(jìn)行的結(jié)果：Mod1 : t= 44.386244 , 29.397192 , 29.397192 , 38.647902 , 38.647902 , 29.52392329.523923 , 44.386244 _Mod2 : t = 34.676012 , 14.846849 , 31.128296 , 31.261550 , 49.893035 , 34.7

15、15948 , _ _44.367994 , 56.320175 _，距（96.000000，28.000000）站20 ，有10111.19的轉(zhuǎn)換關(guān)系，速度v= 5.0，可以看出正演結(jié)果與實(shí)際情況。 Modd2也是如此。4 反演計(jì)算對(duì)于反演模型，主要目的是利用地震波到達(dá)各接收站的時(shí)間，結(jié)合該站的經(jīng)緯度，通過(guò)數(shù)學(xué)推導(dǎo)確定震源參數(shù)。實(shí)際計(jì)算表明，只要輸入數(shù)據(jù)正確，震中初始位置和地震發(fā)生時(shí)間就足夠準(zhǔn)確，可用于地震速報(bào)。具體方法如下。在這里，從前面的正演模型得到的數(shù)據(jù)主要用于反演，然后推導(dǎo)出源參數(shù)。4.1 基本原則4.1.1未知波速V的初步確定方法臺(tái)站n4，震中距小于170km當(dāng)方位角和震中距分布合

16、理時(shí)。假設(shè)接收到某次地震的直波時(shí)有n個(gè)臺(tái)站（ n4 ）。任一站的坐標(biāo)為( xi , y i ) ，對(duì)應(yīng)的直達(dá)波到達(dá)時(shí)間為T i 。待確定的震中位置為x，y，震源深度為z，地震發(fā)生時(shí)間為T（圖4.1）。根據(jù)均質(zhì)地殼模型，可列出以下方程：圖 4.1 源計(jì)算機(jī)位置示意圖(4.1 )i=1,2,n。其中是地震波在地殼中的平均傳播速度，因地區(qū)而異。方程（4.1 ）是一個(gè)非線性方程系統(tǒng)。為了便于直接求解，應(yīng)先線性化，成為線性方程組。因此，(4.1 )展開(kāi)得到(4.2)i=1,2, ,n。對(duì)于第一站， (i=1 )可以寫(xiě)為(4.3)i in (4.3) with 2, 3, ., n 和減去 (4.3) 得

17、到(4.4 )總共有n-1個(gè)線性方程。式（4.4）去掉了深度參數(shù)z，只有四個(gè)未知數(shù)x、y、T、V。由于n-1一般大于4，所以式（4.4）表示的線性方程的個(gè)數(shù)比未知數(shù)的個(gè)數(shù)，構(gòu)成所謂的超定方程組(4.5)為簡(jiǎn)潔起見(jiàn)，這個(gè)超定方程組以矩陣形式表示(4.6)其中，對(duì)于n個(gè)站的數(shù)據(jù)，其中U=V 2 ，W=V 2 T。當(dāng)?shù)卣鹋_(tái)站數(shù)n小于5時(shí)，上述方程為多解不定方程；（A 是A的轉(zhuǎn)置）將超定方程轉(zhuǎn)化為正定方程，計(jì)算源參數(shù)，這就是最小二乘法的原理。只要解出這個(gè)方程組（4.6），也可以得到相應(yīng)的源參數(shù)。這是：震中位置：(x, y )介質(zhì)的波速：V=地震矩：T=W/U。還通過(guò)：將任意臺(tái)站（i）的數(shù)據(jù)代入該公式即

18、可得到焦深：4.1.2當(dāng)波速V已知時(shí)，初步確定的方法是由于研究區(qū)地處青藏高原，地殼厚度較大，所以近期地震的深度基本在上地層。在外殼周圍，它的速度約為 v= 5.8km/s。當(dāng)速度已知時(shí)，上述方法中的線性方程組會(huì)減少一個(gè)未知數(shù)，所以經(jīng)過(guò)（3.5）式變化，(4.7)在W=V 2 T。只要求解這個(gè)方程組（4.7），也可以得到相應(yīng)的源參數(shù)。這是：震中位置：(x, y )地震時(shí)間：T=W/V 2 。同理，得到震源深度。這兩種方法的主要過(guò)程是：圖 4.2 震源計(jì)算機(jī)定位流程圖4.2 模型反演結(jié)果為了更好地檢查程序，前向 t 的結(jié)果中增加了 6s。以下是每種反演方法得到的結(jié)果：1.當(dāng)介質(zhì)速度V未知時(shí)，兩種模

19、型反演的地震參數(shù)為：: 震中 ( 96.000000 , 29.999999 ) t0= 6.0000 25 v=5.000002 h=9.998836mod2：震中（ 95.490894，29.495549 ）t0= 5.797093 v=4.994111 h=20.3332112.當(dāng)介質(zhì)速度V為已知條件（ V =5.0）時(shí)，兩種模型反演的地震參數(shù)為：: 震中 ( 96.000000 , 29.999999 ) t0= 6.0000 00 h=9.999972mod2：震中（ 95.490701,29.495359）t0= 5.874814 h=18.7946024.3 誤差分析從數(shù)據(jù)反演來(lái)

20、看，深度誤差很大，主要是因?yàn)檎鹪瓷疃葘?duì)地震波到達(dá)的讀數(shù)誤差極為敏感，是地震學(xué)中最難準(zhǔn)確測(cè)量的參數(shù)之一。深度估計(jì)具有相當(dāng)大的不確定性，影響人們對(duì)源過(guò)程的理解。具體來(lái)說(shuō)，震源深度的精度誤差受震中距、到時(shí)誤差和速度模型（地殼模型）三個(gè)因素的制約，這些因素對(duì)震源深度的影響是非線性的。1.當(dāng)?shù)卣鸩▊鞑ニ俣纫欢〞r(shí)，震源深度誤差隨著震中距或臺(tái)站位置的增加以及走時(shí)殘差的增加而增加。2 .當(dāng)走時(shí)殘差一定時(shí)，震源深度誤差隨著震中距的增大和地震波速的增大而增大。3.當(dāng)速度已知且走時(shí)殘差不變時(shí)，地震越淺，定位誤差可能越大。這一結(jié)論與許多深部地震的定位結(jié)果是一致的。雖然震源深度是地震學(xué)中難以準(zhǔn)確測(cè)量的參數(shù)之一，而且各種

21、方法測(cè)量震源深度的結(jié)果也不盡相同，但可以對(duì)不同震源深度的結(jié)果進(jìn)行比較分析，從而提高震源深度的準(zhǔn)確度。焦深。在地震發(fā)生后立即通過(guò)在震源區(qū)設(shè)置移動(dòng)觀測(cè)地震臺(tái)站（臺(tái)網(wǎng)）來(lái)修正主震震源深度的準(zhǔn)確性是一種有效的方法。5 實(shí)例分析研究區(qū)位于青藏高原南迦巴瓦段，站點(diǎn)分布如圖5.1所示。5.1 區(qū)域臺(tái)站分布圖及震源示例1. 臺(tái)網(wǎng)記錄的地震事件中，近震記錄根據(jù)國(guó)家先進(jìn)地震系統(tǒng)（ANSS）提供的全球地震事件目錄，事件發(fā)生地點(diǎn)（96.6050，29.5730）已知，地震發(fā)生在2003年第230天09:03?？傆?8個(gè)這次地震事件的記錄。震中距離約為 0.5 0 -3 0 。上述方法反演后得到的位置為：(95.610

22、2, 29.4248)，反演結(jié)果與實(shí)際經(jīng)緯度分別相差(0.0052, 0.1382)。下面的圖片是在所選地震事件前后使用geotool軟件獲得的屏幕截圖。其中，紅色p標(biāo)注的時(shí)間為geotool軟件計(jì)算的地震記錄的跳躍時(shí)間，即理論到達(dá)時(shí)間。黑p標(biāo)記的時(shí)間是人工判斷標(biāo)記的地震道記錄的起飛時(shí)間，即實(shí)際到達(dá)時(shí)間。圖 5.2 中，a 為定位前的地震記錄圖像，b 為定位后的圖像。從圖中可以看出反演后的結(jié)果。反演后，理論到站時(shí)間和實(shí)際起飛時(shí)間還有幾秒的誤差，因?yàn)槟Ｐ彤吘故抢硐肭闆r，實(shí)際數(shù)據(jù)不是那么完善，很難得到極其準(zhǔn)確的實(shí)際運(yùn)行中的車站到達(dá)時(shí)間。一方面是人為因素，比如時(shí)間揀選的準(zhǔn)確性；另一方面，還有機(jī)器設(shè)備

23、的物理因素。而且，由于表層覆蓋層的厚度和松散度等不同，波浪的到來(lái)會(huì)造成誤差。因此，當(dāng)?shù)卣饠?shù)據(jù)到達(dá)時(shí)，需要人工判斷和調(diào)整。在本次地震事件中，每條記錄道的最小誤差為2s，最大誤差為5s。a : 當(dāng)?shù)氐卣鹞恢弥?b: 當(dāng)?shù)氐卣鹞恢弥髨D 5.2 例 1 局地地震定位前后 (a) 和 (b)2.臺(tái)網(wǎng)記錄的地震事件中沒(méi)有近震記錄時(shí)當(dāng)臺(tái)網(wǎng)記錄的地震事件中沒(méi)有近震記錄時(shí)，采用上述反演方法進(jìn)行地震定位，可確定地震事件的經(jīng)緯度為：(96.7124, 25.7474)。它將獲取其經(jīng)度、緯度、地震時(shí)間、深度等信息，按照標(biāo)準(zhǔn)格式編寫(xiě)起源事件，并使用geotool軟件對(duì)事件進(jìn)行處理。結(jié)果如下圖5.3所示，其中最小誤差

24、為0.1s，最大誤差為1.3s。a: 局部地震位置前 b: 局部地震位置后圖 5.3 例 2 局地地震定位前后 (a) 和 (b)為了進(jìn)一步驗(yàn)證這個(gè)定位的正確性，增加了一個(gè)例子。用同樣的反演方法定位后，可以確定地震事件的位置為（96.7124，25.7474）。使用geotool軟件進(jìn)行p波到達(dá)標(biāo)定結(jié)果如圖5.4所示，最小誤差為0.09s，最大誤差為3s。a: 局部地震位置前 b: 局部地震位置后圖 5.4 第二次地震定位前后 (a) 和 (b)六，結(jié)論通過(guò)這個(gè)畢業(yè)設(shè)計(jì)，我了解了近震源參數(shù)定位的方法和原理，從正演反演的理論知識(shí)到實(shí)際數(shù)據(jù)的處理。求解該問(wèn)題的基本思路是將非線性問(wèn)題轉(zhuǎn)化為線性問(wèn)題，

25、然后利用最小二乘原理迭代求解，采用全主元消元法求解。本次設(shè)計(jì)的目的是借助源計(jì)算機(jī)初步定位問(wèn)題，進(jìn)一步了解通用科學(xué)計(jì)算研究所采用的正反演計(jì)算方法和實(shí)際應(yīng)用。這個(gè)畢業(yè)設(shè)計(jì)是在王教授的悉心指導(dǎo)下，不斷學(xué)習(xí)和修改完成的。王教授淵博的學(xué)識(shí)、豐富的實(shí)踐經(jīng)驗(yàn)、遠(yuǎn)見(jiàn)卓識(shí)和敏銳的科學(xué)眼光，將是我永遠(yuǎn)學(xué)習(xí)的榜樣。我深深佩服他樂(lè)觀的態(tài)度，謙遜的人，簡(jiǎn)單的生活態(tài)度。他一絲不茍的工作作風(fēng)將成為我今后學(xué)習(xí)和工作的榜樣。他有說(shuō)服力的教義和不拘一格的思想給了我無(wú)窮的靈感。我由衷的感謝和崇高的敬意。畢業(yè)論文的撰寫(xiě)也得到了眾多師生的關(guān)心和幫助。在這里我要表達(dá)我誠(chéng)摯的謝意。祝愿他們身體健康，工作順利，事業(yè)更上一層樓。由于本人水平有

26、限，本文及程序有待進(jìn)一步完善，敬請(qǐng)批評(píng)指正。參考1 賽斯斯坦，邁克爾威塞辛。地震學(xué)和地球結(jié)構(gòu)概論。布萊克威爾出版社，2003 年。2 山幫，中國(guó)地震。地震，1981 年。3 單娜琳，程志平，云珍工程地震勘探。冶金工業(yè)，2006。4 梅，王有雪，錢慧造山高原-青藏高原及鄰區(qū)寬帶地震探測(cè)與地殼上地幔結(jié)構(gòu)J.地質(zhì)學(xué)，2009 年。5 熊章強(qiáng)，方根先淺層地震勘探。地震，2002 年。6 大港油田科技叢書(shū)編委。地震勘探數(shù)據(jù)處理與解釋技術(shù)。石油工業(yè)，1999 年。7 石振良，紹全，榮國(guó)等。地震工作手冊(cè)。地震，1990 年。8 付淑芳，寶成，文藝地震學(xué)課程（第 1 卷和第 2 卷）。地震，1980 年。9

27、郭倫 Fortran 95 編程。中國(guó)電力，2002。附錄程序程序測(cè)試隱式無(wú)整數(shù) : i, j, k, nstn, istn真實(shí)*8 : x0, y0, v0, t0, xx, yy, coef(1000,5)實(shí)數(shù)*8 : delt(1000), solut(4,5), x(1000), y(1000), t(1000)real*8 : ph0, la0, phi, lai ,v,h, MinLat, MaxLat, MinLon, MaxLon字符 : 文件名*100write(*,*) 輸入文件名讀取(*,(a) 文件名!閱讀所有電臺(tái)的精選打開(kāi)（1，文件=文件名）讀取(1,*) nstnx

28、0=0;y0=0做 istn=1,nstn讀取(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)結(jié)束做la0=x0/nstnph0=y0/nstn關(guān)閉(1)轉(zhuǎn)到 10!計(jì)算綜合數(shù)據(jù)la0=99*atan(1.00)/45ph0=27*atan(1.00)/45打開(kāi)（2，文件=mod13.dat）寫(xiě)(2,(I2) nstnv=5.0h=10.0做 istn=1,nstn來(lái)=x(istn)phi=y(istn)調(diào)用 lp

29、2xy(lai,phi,la0,ph0,xx,yy)t(istn)=sqrt(xx*2.00+yy*2.00+h*2.00)/v寫(xiě)(2,( 6f15.6) x(istn)*45.00/atan(1.00), y(istn)*45.00/atan(1.00), t(istn)+6結(jié)束做寫(xiě)(2,( 4f20.6) la0*45.00/atan(1.00),ph0*45.00/atan(1.00),v,h關(guān)閉(2)！停止10 繼續(xù)!(lai,phi)-(x,y)做 istn=1,nstn來(lái)=x(istn)phi=y(istn)調(diào)用 lp2xy(lai,phi,la0,ph0,xx,yy)x(istn

30、)=xxy(istn)=yy結(jié)束做!計(jì)算線性方程的系數(shù)做 istn=2,nstncoef(istn,1)=x(istn)-x(1)coef(istn,2)=y(istn)-y(1)coef(istn,3)=(t(istn)*2-t(1)*2)/2.00coef(istn,4)=-t(istn)+t(1)delt(istn) =(x(istn)*2-x(1)*2+y(istn)*2-y(1)*2)/2.00結(jié)束做!最小平方根法得到解!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

31、) de(2,5)！ .！ .！ .!dt(n) de(n,1) de(n,2) de(n,3) de(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)做我=1,4做 j=1,4溶液(i,j)=0k=1,nstn溶液(i,j)=溶液(i,j)+coef(k,i)*coef(k,j) 結(jié)束做結(jié)束做結(jié)束做做我=1,4解決方案(i,5)=0k=1,nstnsolut(

32、i,5)=solut(i,5)+coef(k,i)*delt(k)結(jié)束做結(jié)束做調(diào)用高斯(solut,4,5)xx=解決方案（1,5）yy=解決方案（2,5）v0=sqrt（解決方案（3,5） t0=解決方案（4,5）/解決方案（3,5）h=0！打印 *,xx,yy做 istn=1,nstnh=h+sqrt(v0*v0*(t(7)-t0)*(t(7)-t0)-(x(7)-xx)*2-(y(7)-yy)* 2)結(jié)束做 h=h/nstn調(diào)用 xy2lp(xx,yy,la0,ph0,lai,phi)print *,-lai phi v0 t0 h-打印 *, lai*45.00/atan(1.00),phi*45.00/atan(1.00),v0,t0,h結(jié)束程序測(cè)試子程序高斯(A,N1,M1)隱式無(wú)整數(shù) :n1,m1,k1,k2,ik,i,j ,l,i1實(shí)數(shù)*8 : A(n1,m1)真實(shí)*8 : BMAX,T,EPS每股收益=0。DO K1=1,N1BMAX=0。DO I=K1,N1IF(BMAX-ABS(A(I,K1).LT.0) 那么BMAX=ABS(A(I,K1)L=我萬(wàn)一結(jié)束做IF(BMAX.LT.EPS) 停止 4444如果（L.NE.K1）那么做 J=K1,M1T=A(L,J)A(L,J)=A(K1,J