關(guān)于地震問題的數(shù)學(xué)模型探究

關(guān)于地震問題的數(shù)學(xué)模型探究吳瑞劉思謝磊摘要地震定位是地震學(xué)中最經(jīng)典、最基本的問題之一，對(duì)于研究諸如地震活動(dòng)構(gòu)造、地球內(nèi)部結(jié)構(gòu)、震源的幾何構(gòu)造等此類地震學(xué)中的基本問題有重要意義。此外，基于快速準(zhǔn)確的地震定位的地震速報(bào)，對(duì)于震后的減災(zāi)、救災(zāi)工作也是至關(guān)重要的。因此，地震學(xué)家一直在不斷改進(jìn)或提出新的定位方法。地震定位問題的提法如下:根據(jù)臺(tái)站對(duì)地震到時(shí)的觀測(cè)資料，來(lái)確定震源的空間坐標(biāo)和發(fā)震時(shí)刻，有時(shí)還給出對(duì)解的評(píng)價(jià)。早期的地震定位方法以幾何作圖法為主。近三十年來(lái)由于計(jì)算機(jī)技術(shù)的飛速發(fā)展和廣泛應(yīng)用，基于科學(xué)計(jì)算和計(jì)算機(jī)技術(shù)的智能化數(shù)值自動(dòng)定位方法也得到了迅速發(fā)展，并業(yè)己成為當(dāng)前地震定位的主流方法。我國(guó)最初的地震定位工作由李善邦先生于1930年在北京鶩峰地震臺(tái)開創(chuàng)，1953年開始采用多臺(tái)站大規(guī)模觀測(cè)數(shù)據(jù)確定震中，現(xiàn)在大多使用國(guó)際流行的定位方法。針對(duì)本文的數(shù)據(jù)處理方面，主要運(yùn)用MATLAB軟件進(jìn)行數(shù)據(jù)處理。有關(guān)的圖像處理，也是運(yùn)用MATLAB軟件進(jìn)行處理即根據(jù)局部地形、地貌等高線圖，繪制相應(yīng)的三維地形、地貌曲面圖。關(guān)鍵字:最小二乘法Geiger方法差值擬合P波S波前言1.問題的背景中國(guó)位于世界兩大地震帶一環(huán)太平洋地震帶與歐亞地震帶的交匯部位，受太平洋板塊、印度板塊和菲律賓海板塊的擠壓，地震斷裂帶十分發(fā)育。大地構(gòu)造位置決定，地震頻繁震災(zāi)嚴(yán)重。在20世紀(jì)里，全球共發(fā)生3次8.5級(jí)以上的強(qiáng)烈地震，其中兩次發(fā)生在我國(guó);全球發(fā)生兩次導(dǎo)致20萬(wàn)人死亡的強(qiáng)烈地震也都發(fā)生在我國(guó)，一次是1920年寧夏海原地震，造成23萬(wàn)多人死亡;一次是1976年河北唐山地震，造成24萬(wàn)多人死亡。這兩次地震死亡人數(shù)之多，在全世界也是絕無(wú)僅有的。上世紀(jì)以來(lái)，中國(guó)共發(fā)生6級(jí)以上地震近800次，遍布除貴州、浙江兩省和香港、澳門特別行政區(qū)以外所有的省、自治區(qū)、直轄市。中國(guó)地震活動(dòng)頻度高、強(qiáng)度大、震源淺，分布廣，是一個(gè)震災(zāi)嚴(yán)重的國(guó)家。1900年以來(lái)，中國(guó)死于地震的人數(shù)達(dá)55萬(wàn)之多，占全球地震死亡人數(shù)的53%;1949年以來(lái)，100多次破壞性地震襲擊了22個(gè)省(自治區(qū)、直轄市)，其中涉及東部地區(qū)14個(gè)省份，造成27萬(wàn)余人喪生，占全國(guó)各類災(zāi)害死亡人數(shù)的54%，地震成災(zāi)面積達(dá)30多萬(wàn)平方公里，房屋倒塌達(dá)700萬(wàn)間。地震及其他自然災(zāi)害的嚴(yán)重性構(gòu)成中國(guó)的基本國(guó)情之一。統(tǒng)計(jì)數(shù)字表明，中國(guó)的陸地面積占全球陸地面積的十五分之一，即百分之六左右;中國(guó)的人口占全球人口的五分之一左右，即百分之二十左右，都不到百分之二十，然而中國(guó)的陸地地震竟占全球陸地地震的三分之一，即百分之三十三左右，而造成地震死亡的人數(shù)竟達(dá)到全球的1/2以上。當(dāng)然這也有特殊原因，一是中國(guó)的人口密、人口多;中國(guó)的經(jīng)濟(jì)落后，房屋不堅(jiān)固，容易倒塌，容易壞;第三與中國(guó)的地震活動(dòng)強(qiáng)烈頻繁有密切關(guān)系。統(tǒng)計(jì)，20世紀(jì)以來(lái)，中國(guó)因地震造成死亡的人數(shù)，占國(guó)內(nèi)所有自然災(zāi)害包括洪水、山火、泥石流、滑坡等總?cè)藬?shù)的54%，超過1/2。從人員的死亡來(lái)看，地震是群害之首;而在經(jīng)濟(jì)上所造成的損失，最大的主要是氣象災(zāi)害(洪澇)，氣象我國(guó)的地震活動(dòng)主要分布在五個(gè)地區(qū)的23條地震帶上。這五個(gè)地區(qū)是:①臺(tái)灣省及其附近海域;②西南地區(qū)，主要是西藏、四川西部和云南中西部;③西北地區(qū)，主要在甘肅河西走廊、青海、寧夏、天山南北麓;④華北地區(qū)，主要在太行山兩側(cè)、汾渭河谷、陰山一燕山一帶、山東中部和渤海灣;⑤東南沿海的廣東、福建等地。我國(guó)的臺(tái)灣省位于環(huán)太平洋地震帶上，西藏、新疆、云南、四川、青海等省區(qū)位于喜馬拉雅一地中海地震帶上，其他省區(qū)處于相關(guān)的地震帶上。災(zāi)害所造成的經(jīng)濟(jì)損失要比地震大的多。圖1我國(guó)地震帶分布圖圖2中國(guó)強(qiáng)震及地震分布圖圖3中國(guó)地震和火山圖4中國(guó)地震震中分布圖2.概念闡述2.1最小二乘法最小二乘法(又稱最小平方法)是一種數(shù)學(xué)優(yōu)化技術(shù)，它通過最小化誤差的平方和找到一組數(shù)據(jù)的最佳函數(shù)匹配。最小二乘法是用最簡(jiǎn)的方法求得一些絕對(duì)不可知的真值，而令誤差平方之和為最小。最小二乘法通常用于曲線擬合。很多其他的優(yōu)化問題也可通過最小化能量或最大化嫡用最小二乘形式表達(dá)。人們往往對(duì)由某一變量七或多個(gè)變量f......f構(gòu)成的相關(guān)變量Y感興趣。如彈簧的形變與所用的力相關(guān)，一個(gè)企業(yè)的盈利與其營(yíng)業(yè)額，投資收益和原始資本有關(guān)。為了得到這些變量同Y之間的關(guān)系，便用不相關(guān)變量去構(gòu)建Y。嘗試用如下函數(shù)模型(1)最小化問題的精度，依賴于所選擇的函數(shù)模型。q個(gè)相關(guān)變量或p個(gè)附加的相關(guān)變量去擬和。通常人們將一個(gè)可能的、對(duì)不相關(guān)變量七的構(gòu)成都無(wú)困難的函數(shù)類型充作函數(shù)模型(如拋物線函數(shù)或指數(shù)函數(shù))。參數(shù)x是為了使所選擇的函數(shù)模型同觀測(cè)值Y相匹配。(如在測(cè)量彈簧形變時(shí)，必須將所用的力與彈簧的膨脹系數(shù)聯(lián)系起來(lái))。其目標(biāo)是合適地選擇參數(shù)，使函數(shù)模型最好的擬合觀測(cè)值。一般情況下，觀測(cè)值遠(yuǎn)多于所選擇的參數(shù)。其次的問題是怎樣判斷不同擬合的質(zhì)量。高斯和勒讓德的方法是，假設(shè)測(cè)量誤差的平均值為0。令每一個(gè)測(cè)量誤差對(duì)應(yīng)一個(gè)變量并與其它測(cè)量誤差不相關(guān)(隨機(jī)無(wú)關(guān))。人們假設(shè)，在測(cè)量誤差中絕對(duì)不含系統(tǒng)誤差，它們應(yīng)該是純偶然誤差，圍繞真值波動(dòng)。除此之外，測(cè)量誤差符合正態(tài)分布，這保證了偏差值在最后的結(jié)果Y上忽略不計(jì)。確定擬合的標(biāo)準(zhǔn)應(yīng)該被重視，并小心選擇，較大誤差的測(cè)量值應(yīng)被賦予較小的權(quán)。并建立如下規(guī)則:被選擇的參數(shù)，應(yīng)該使算出的函數(shù)曲線與觀測(cè)值之差的平方和最小。用函數(shù)表示為:(2)用歐幾里德規(guī)則表達(dá)為:（3）最小化問題的精度，依賴于所選擇的函數(shù)模型。2.2P波在測(cè)震學(xué)中，震中距在1000公里以上的稱為遠(yuǎn)震。如果把測(cè)量點(diǎn)作為地球的一極，穿過地心和測(cè)量點(diǎn)相對(duì)應(yīng)的為另一極，中間的距離按照和地心與測(cè)量點(diǎn)之間的張角劃分為180度，也可以說震中距在10度以內(nèi)的稱為近震，10度一一105度的為遠(yuǎn)震，105度以上的稱為極遠(yuǎn)震。對(duì)于遠(yuǎn)震和極遠(yuǎn)震，不能再象近震一樣認(rèn)為地震波是直線傳播的，必須考慮地震波在傳播時(shí)的折射。從震源發(fā)出，一路不斷折射傳播過來(lái)的縱波叫做P波。它是遠(yuǎn)震測(cè)量中的首波。P波(P-waveorprimarywave)是二種體波(體波的命名是因?yàn)榇瞬ù┰降厍騼?nèi)部，相對(duì)于體波的是表面波，另一種體波是S波(secondarywave)中的一種。P波意指(primarywave)或是壓力波(pressurewave)。在所有地震波中，P波擁有最快的傳遞度速，因此地震發(fā)生時(shí)，P波是最早抵達(dá)測(cè)站，并被地震儀紀(jì)錄下來(lái)的地震波，這也是P波名稱的由來(lái)。P波的P也能代表壓力(pressure)，來(lái)自于其震動(dòng)傳遞類似聲波，屬于縱波的一種(或疏密波)，傳遞時(shí)介質(zhì)的震動(dòng)方向與震波能量的傳播方向平行。對(duì)于地球內(nèi)部構(gòu)造的了解和推論，大部分是藉由觀測(cè)地震波中的體波。地震波在不同介質(zhì)有不同傳播時(shí)間和路徑，在介質(zhì)交界面時(shí)會(huì)產(chǎn)生反射、折射，以及相位的改變，地震學(xué)家利用這些特性來(lái)獲得地球內(nèi)部資訊。當(dāng)體波穿越地球液態(tài)層時(shí)，P波在經(jīng)過下部地函與外地核時(shí)會(huì)稍許折射。造成P波在104。至1400間會(huì)有陰影區(qū)，即地震儀記錄不到。2.3S波S波(S-waveorsecondarywave)是二種體波(體波的命名是因?yàn)椴ù┰降厍騼?nèi)部，相對(duì)于體波的是表面波)中之一。它是因地震而產(chǎn)生的，被地震儀記錄下來(lái)。命名為S波(secondarywave)是因?yàn)樗乃俣葍H次于P波(最快的球內(nèi)部粒子的震動(dòng)方向與震波能量傳遞方向是垂直的。s波與P波不同的是，S波無(wú)法穿越外地核。所以S波的陰影區(qū)正對(duì)著地震的震源。2.4走時(shí)地震走時(shí)是指地震波從震源傳到觀測(cè)點(diǎn)所經(jīng)過的時(shí)間。地震縱波和橫波傳播的速度不同，同一距離傳播時(shí)間也不同。傳播時(shí)間的長(zhǎng)短，還受到所經(jīng)之處的物質(zhì)組成情況的影響，準(zhǔn)確地測(cè)出走時(shí)，對(duì)研究地震有重要意義。3.模型的建立與求解3.1問題一的模型的建立與求解3.1.1問題的重述假定地面是一個(gè)平面，在這個(gè)平面上建立坐標(biāo)系見圖5。圖5中給出了10個(gè)地震觀測(cè)站點(diǎn)(A-J)的坐標(biāo)位置。圖5地震觀測(cè)站點(diǎn)示意圖2009年4月1口某時(shí)在某一地點(diǎn)發(fā)生了一次地震，圖5中10個(gè)地震觀測(cè)站點(diǎn)均接收到了地震波，觀測(cè)數(shù)據(jù)見表1。地震觀測(cè)站橫坐標(biāo)x（千米）縱坐標(biāo)y（千米）接收地震波時(shí)間A50033004月1日9時(shí)21分9秒B3002004月1日9時(shí)19分29秒C80016004月1日9時(shí)14分51秒D140022004月1日9時(shí)13分17秒E17007004月1日9時(shí)11分46秒F230028004月1日9時(shí)14分47秒G250019004月1日9時(shí)10分14秒H29009004月1日9時(shí)11分46秒I320031004月1日9時(shí)17分57秒J34001004月1日9時(shí)16分49秒表1地震觀測(cè)站坐標(biāo)及接收地震波時(shí)間假定地震波在各種介質(zhì)和各個(gè)方向的傳播速度均相等，并且在傳播過程中保持不變?，F(xiàn)在需要根據(jù)表1中的數(shù)據(jù)確定這次地震的震中位置、震源深度以及地震發(fā)生的時(shí)間。3.1.2模型一的確定思路由于P波和S波傳播的速度不同，我們采用均勻介質(zhì)模型，后引入虛波，利用公式“路程=距離*時(shí)間”建立模型，提出一種線性單事件模型。3.1.3模型一的建立設(shè)震源坐標(biāo)為，發(fā)震時(shí)刻為，臺(tái)站i的坐標(biāo)為P波、S波到時(shí)分別為，震源到臺(tái)站i的距離為。采用均勻介質(zhì)模型，設(shè)P波、S波速度分別為，引入虛波速度:則有:，從而（4）對(duì)站臺(tái)i有：，（5）對(duì)站臺(tái)j有：，（6）上面兩式展開后相減，即可消去二次項(xiàng)，得線性方程：，（7）其中。將線性方程用于3對(duì)以上的站臺(tái)，即得到一關(guān)于的線性方程組，從而反演得到震源位置。由于z:表示臺(tái)站高度，當(dāng)臺(tái)站不具有深度時(shí)，遠(yuǎn)小于和，使得方程組的系數(shù)矩陣具有奇異性，導(dǎo)致z。解發(fā)散?？梢酝ㄟ^其它方法單獨(dú)求解：取震中解，將其代入(5)式求得，將各個(gè)臺(tái)站求得的的平均值作為震源深度的解。這樣的誤差僅決定于震中的定位誤差和中的到時(shí)讀數(shù)誤差(二者均是可控的)，從根本上解決了無(wú)法得到震源深度有效信息的問題。關(guān)于發(fā)震時(shí)刻，有方程：，（8）將方程(8)用于各個(gè)臺(tái)站，求各t。的平均值作為發(fā)震時(shí)刻的解。由于完全取決于到時(shí)(見((4)式)，所以和震源的求解完全獨(dú)立，其誤差的唯一來(lái)源是到時(shí)讀數(shù)的誤差，只要到時(shí)讀數(shù)足夠精確，即使震源定位誤差很大，也可求得很精確的。3.1.4模型一的求解采用奇異值分解(SVD)的廣義逆求解。將方程組表示為:，（9）其中：，，.對(duì)矩陣A進(jìn)行奇異值分解:,(10)其中，，。滿足：，，奇異值是矩陣的p個(gè)非零特征值的正平方根。則方程組(9)的廣義逆解為:（11）因?yàn)樵挤匠探M即為線性，故求解過程不需要迭代，也不存在初值選取問題。并且矩陣A的計(jì)算只需要簡(jiǎn)單的算術(shù)加減，非常容易，計(jì)算量小。當(dāng)A的條件數(shù)很大時(shí)，可以采用阻尼廣義逆解:（12）其中I為單位矩陣，少為阻尼系數(shù)，取值隨情況而變，對(duì)于震中解，取值一般在0~1之間。此時(shí)直接采用（11）式求震中解。前己述及，當(dāng)臺(tái)站不具有深度時(shí)，(7)式中的遠(yuǎn)小于和，這表現(xiàn)為矩陣A的3個(gè)奇異值中，最小值遠(yuǎn)小、于其它兩個(gè)奇異值，即A的條件數(shù)很大。而對(duì)應(yīng)模型空間中z方向的信息，故深度誤差己將到時(shí)讀數(shù)誤差(包含在b中)放大了許多倍，用廣義逆及阻尼廣義逆均無(wú)法得到深度的有效信息。本模型充分利用臺(tái)站與震源的距離可由P波和S波的到時(shí)差事先確定的特點(diǎn)，直接由簡(jiǎn)單的距離關(guān)系求得震源深度，避免了矩陣奇異問題。但是臺(tái)站一旦具有深度，增大，矩陣A的條件數(shù)將很快減小，此時(shí)廣義逆解中的深度解將獲得很大改善。再次利用，可事先確定的特點(diǎn)，使與震源參數(shù)完全解耦，發(fā)震時(shí)刻解的精度完全不受震源定位精度的影響。3.1.5模型一的討論選取哪兩個(gè)臺(tái)站的方程相減構(gòu)成方程組，對(duì)定位結(jié)果有著很大的影響。將臺(tái)站按照所有可能進(jìn)行兩兩組合，這是一種操作最簡(jiǎn)單的“盲組合”，并且含有的物理信息也最豐富。實(shí)驗(yàn)證明這樣的定位結(jié)果是最理想的圖(6)。本文分別用表示震中、震源深度和發(fā)震時(shí)刻的定位誤差。圖6’o’表示臺(tái)站,’+’表示震中。兩臺(tái)站間的連線表示這兩個(gè)臺(tái)站的方程相減所得到的方程,第4種情況含有所有可能的臺(tái)站連線,即所有可能的兩兩組合,所得誤差最小。圖6臺(tái)站組合圖如果對(duì)一個(gè)臺(tái)站的到時(shí)在其讀數(shù)誤差范圍內(nèi)進(jìn)行隨機(jī)擾動(dòng),得到一系列不同的到時(shí),相當(dāng)于在該臺(tái)站附近增加了一系列虛擬臺(tái)站,對(duì)增加后的臺(tái)站進(jìn)行所有可能的兩兩組合,(從而增加了方程個(gè)數(shù)),可以進(jìn)一步提高定位結(jié)果的精度,(圖7)。實(shí)際上,到時(shí)讀數(shù)的確切誤差并不知道,只知道一個(gè)誤差范圍,增加虛擬臺(tái)站相當(dāng)于在該誤差范圍內(nèi)取平均效果,一般說來(lái),這比單個(gè)的到時(shí)讀數(shù)更接近準(zhǔn)確到時(shí),在數(shù)學(xué)形式上則表現(xiàn)為方程個(gè)數(shù)增加,從而所含信息量增大,故結(jié)果的精確度提高。圖7誤差曲線分析圖考慮到表1提供的數(shù)據(jù)只有10組,我們無(wú)法利用模型一得到較精確地結(jié)果。于是我們又提出了基于幾組數(shù)據(jù)即可得到較準(zhǔn)確地?cái)?shù)學(xué)模型。3.1.6模型二的確定思路假設(shè)采用均勻介質(zhì)模型,利用空間坐標(biāo)公式,求出臺(tái)站到震源的距離,利用“路程=距離×?xí)r間”算出走時(shí),再根據(jù)地震發(fā)生時(shí)間與地震臺(tái)站接收時(shí)間的關(guān)系列出函數(shù)表達(dá)式。3.1.7模型二的建立假設(shè)震源三維坐標(biāo)為,單位km。設(shè)發(fā)生地震的時(shí)間為2008年4月1日9時(shí)分，地震波傳播速度為km/s，,i=1,2,...分別表示地震觀測(cè)站點(diǎn)A—J的三維坐標(biāo)，用表示9時(shí)分接收到地震波。根據(jù)題設(shè)條件和以上假設(shè)可得如下模型:走時(shí)=震源到觀測(cè)臺(tái)站的距離/地震波傳播速度，觀測(cè)臺(tái)站接收時(shí)間=地震發(fā)生時(shí)間+走時(shí)。于是我們得到（13）其中，，，為未知參數(shù)。3.1.8模型二的求解令，，由(13)式可以看出t是x和y的函數(shù)，函數(shù)表達(dá)式中的未知參數(shù)，，，可由最小二乘法估計(jì)得到，相應(yīng)的matlab程序如下：function[b,fval]=dizhenxyz=[500,3300,0300,200,0800,1600,01400,2200,01700,700,02300,2800,02900,900,03200,3100,03400,100,0];min=[21191413111410111716];sec=[9295117464714465749];time=min+sec/60;b=nlinfit(xyz,time,@dizhenmoxing,[1000100111]);functiony=dizhenmoxing(beta,x)xb=x(:,1);yb=x(:,2);zb=x(:,3);y=sqrt((xb-beta(1)).^2+(yb-beta(2)).^2+(zb-beta(3)).^2/(60*beta(4))+beta(5);運(yùn)行以上程序得到：=2200.5=1399.9=35.1（14）=3=7將以上結(jié)果代入模型(13)，算得各觀測(cè)站點(diǎn)理論上接收到地震波的時(shí)間見表2。觀測(cè)站接收地震波實(shí)際時(shí)間(分)接收地震波理論時(shí)間(分)時(shí)間差(分)A21.1521.1566-0.0066B19.483319.4770.0063C14.8514.84990.0001D13.283313.27830.005E11.766711.7761-0.0049F14.783314.7881-0.0048G10.233310.23120.0021H11.766711.7677-0.001I17.9517.94640.0036J16.816716.81660.0001表2各觀測(cè)站點(diǎn)理論上接收到地震波的時(shí)間【注:時(shí)間差=接收地震波實(shí)際時(shí)間-接收地震波理論時(shí)間】3.1.9結(jié)果的探究由上表可知各觀測(cè)站接收到地震波的實(shí)際時(shí)間與理論時(shí)間之差比較小，結(jié)果(14)可以接受，即地震震中位置為原題圖1中點(diǎn)(2200.5,1399.9)處，震源深度為35.1千米，地震發(fā)生時(shí)刻為2008年4月1日9時(shí)7分。圖8為模型（13）求解后的立體曲線圖。即：的曲線圖圖8立體曲線圖附：matlab相應(yīng)程序x=linspace(0,5000,50);y=linspace(0,5000,50);[xx,yy]=meshgrid(x,y);zz=7+sqrt((xx-2200.5).^2+(yy-1399.9).^2+35.1*35.1)/180;plot3(xx,yy,zz)3.2問題二的模型的建立與求解3.2.1問題的重述地震災(zāi)區(qū)的地形、地貌對(duì)抗震救災(zāi)的進(jìn)展會(huì)有很大影響，根據(jù)衛(wèi)星遙感和飛機(jī)航拍得到的照片可以構(gòu)建災(zāi)區(qū)地形、地貌圖。所構(gòu)建的災(zāi)區(qū)局部地形、地貌等高線圖見圖9。圖9災(zāi)區(qū)局部地形、地貌等高線圖現(xiàn)在需要根據(jù)圖9中的局部地形、地貌等高線圖，建立數(shù)學(xué)模型，繪制出相應(yīng)的三維地形、地貌曲面圖。3.2.2圖像的處理首先對(duì)圖9進(jìn)行預(yù)處理，處理效果如圖10。圖10災(zāi)區(qū)局部地形、地貌等高線圖預(yù)處理效果圖將圖10保存成文件dx2.jpg，然后利用imread命令導(dǎo)入matlab進(jìn)行自動(dòng)識(shí)別，得出各條等高線所對(duì)應(yīng)點(diǎn)的三維坐標(biāo)，然后在圖10中人為選定一些參照點(diǎn)，現(xiàn)在選取如下參照點(diǎn)：(1,1,0)(1,2,0)(1,3,0)(1,4,0)(1,5,0)(2,1,0)(3,1,0)(4,1,0)(5,1,0)(5,2,0)(5,3,0)(5,4,0)(5,5,0)(2,5,0)(3,5,0)(4,5,0)對(duì)自動(dòng)識(shí)別出的各條等高線上的點(diǎn)和人為選定的參照點(diǎn)的三維坐標(biāo)進(jìn)行二維插值擬合，可以得到插值擬合曲面。相應(yīng)的matlab程序如下：I=imread('dx2.jpg');BW=im2bw(I);[m,n]=size(BW);x_xishu=4/(n-1);y_xishu=4/(m-1);gao=[0.1:0.1:0.7];[BW,n]=bwlabel(BW);xyz=[];fori=1:n;BW1=BW;BW1(BW1~=i)=0;[y,x]=find(BW1);xyz=[xyz;110;120;130;140;150;210;310;410;510;520;530;540;550;250;350;450];x=xyz(:1);y=xyz(:2);z=xyz(:3);[x1,y1]=meshgrid(1:0.1:5);z1=griddata(x,y,z,x1,y1,'v4');figuresurf(x1,y1,z1)zlim([00.8])插值擬合曲面，即與圖9中的局部地形、地貌等高線圖相應(yīng)的三維地形、地貌曲面圖如下：圖10三維地形、地貌曲面圖3.3問題三的模型的建立與求解3.3.1問題的重述現(xiàn)在需要查閱有關(guān)資料，了解地震波在各種介質(zhì)和各個(gè)方向的傳播速度問題，給出合理假設(shè)，根據(jù)表1中的數(shù)據(jù)確定這次地震的震中位置、震源深度以及地震發(fā)生的時(shí)間。3.3.2模型的確定思路在假設(shè)模型的過程中，我們查閱了大量的關(guān)于地震問題的文獻(xiàn)，我們發(fā)現(xiàn)Geiger模型對(duì)我們的問題具有指導(dǎo)性意義。我們將問題一的解作為Geiger方法的初值，因?yàn)槌踔惦x準(zhǔn)確解已經(jīng)相當(dāng)近，所以能夠保證迭代過程收斂至全局極小值，有效地提高了定位精度。3.3.3模型的建立設(shè)個(gè)臺(tái)站的觀測(cè)到時(shí)為，求震源及發(fā)震時(shí)刻使得目標(biāo)函數(shù)（15）最小。其中為到時(shí)殘差（16）為震源到第i個(gè)臺(tái)站的計(jì)算走時(shí)。3.3.4模型的求解使目標(biāo)函數(shù)取極小值也即（17）其中。為方便，記（18）則由(18)式,在真解附近任意試探解及其校正矢量滿足(19)也即（20）由的定義可得公式(20)的具體表達(dá)式（21）若偏離真解不大，則和較小，可忽略二階導(dǎo)數(shù)項(xiàng)，（21）式被簡(jiǎn)化為線性最小二乘解：（22）以矩陣形式表示，上式為（23）其中，.若二階導(dǎo)數(shù)項(xiàng)不可忽略，則(12)式給出非線性最小二乘解（24）通常各臺(tái)站的到時(shí)數(shù)據(jù)具有不同的精度，如果不加以區(qū)別，則具有較低精度的數(shù)據(jù)將嚴(yán)重干擾結(jié)果的精度，這一問題可以通過引入加權(quán)目標(biāo)函數(shù)來(lái)解決。設(shè)各臺(tái)站到時(shí)殘差的方差為，引入加權(quán)目標(biāo)函數(shù)，（25）按照上述同樣的步驟，通過求(25)式的極小值，得到如下加權(quán)線性最小二乘解，（26）其中為加權(quán)方差矩陣：。由方程(24),(25),或(26)求得后，以作為新的嘗試點(diǎn)，再求解相應(yīng)方程。如此反復(fù)迭代，直至或足夠?。ɑ驖M足一定的循環(huán)結(jié)束條件），此時(shí)即得估計(jì)解。然而，Geiger方法也有它的不足之處。Geiger方法強(qiáng)烈地依賴于初值的選取，如果初值不夠好，會(huì)導(dǎo)致迭代過程的發(fā)散，根本無(wú)法求得正確解，當(dāng)初值離準(zhǔn)確解較近時(shí)，能保證迭代收斂，且迭代次數(shù)較少。所以本題我們選取了模型二的解值作為我們的初值。3.3.5模型的進(jìn)一步探究我們想到能否用另一個(gè)模型可以補(bǔ)充Geiger方法的不足之處，考慮到地震波在各方向傳播的問題，我們采用震源位置與速度結(jié)構(gòu)的聯(lián)合反演法（SSH），SSH法不需要對(duì)波速進(jìn)行校準(zhǔn)，同時(shí)還可以獲得有關(guān)速度結(jié)構(gòu)的很多信息，是目前被廣泛使用的一