地震科普
地震學(xué)百科知識(shí)(八)
——計(jì)算地震學(xué)*
王彥賓1)陳曉非2)
1)北京大學(xué)地球與空間科學(xué)學(xué)院,北京100871
2)中國科學(xué)技術(shù)大學(xué)地球與空間科學(xué)學(xué)院,合肥230026
研究理論地震學(xué)計(jì)算問題所涉及的數(shù)值算法以及計(jì)算機(jī)實(shí)現(xiàn)是一門分支學(xué)科。它的誕生得益于計(jì)算機(jī)的廣泛應(yīng)用和大量的高質(zhì)量地震波形觀測資料的獲得,是介于理論地震學(xué)、計(jì)算科學(xué)和觀測地震學(xué)之間的新興交叉學(xué)科。
基于基本的物理定律和嚴(yán)格的數(shù)學(xué)工具,理論地震學(xué)給出了定量描述地震震源自發(fā)破裂動(dòng)力學(xué)過程、地震波激發(fā)、地震波傳播等問題的控制方程。但是,這些方程只對(duì)極少數(shù)相對(duì)簡單的地球模型能夠給出完整的解析解。隨著全球范圍內(nèi)地震觀測資料數(shù)量的增多和覆蓋區(qū)域的擴(kuò)大,地震學(xué)家認(rèn)識(shí)到地震震源的復(fù)雜性和地球內(nèi)部結(jié)構(gòu)的非均勻性。對(duì)于更接近真實(shí)地球的復(fù)雜模型,控制方程不能給出解析解。隨著計(jì)算機(jī)的出現(xiàn)和數(shù)值算法的發(fā)展,各類控制方程都可以通過數(shù)值方法快速、精確地求解,對(duì)于一些相對(duì)簡單的復(fù)雜模型,可以結(jié)合數(shù)值方法給出半解析解,對(duì)于任意復(fù)雜模型,通常只能依賴數(shù)值方法給出數(shù)值近似解。
1950—1955年間數(shù)字計(jì)算機(jī)開始引入地震學(xué)研究,帶來地震學(xué)的“第二次革命”[1]。計(jì)算地震學(xué)作為一門交叉學(xué)科在這一背景下產(chǎn)生。Thomson[2]和Haskell[3]分別于1950年和1953年給出了水平多層介質(zhì)中面波頻散的數(shù)值算法,1955年P(guān)ekeris[4]首次利用計(jì)算機(jī)數(shù)值求解蘭姆(H.Lamb)問題(半無限均勻彈性空間內(nèi)由點(diǎn)源或線源激發(fā)的彈性波傳播問題),給出了泊松固體介質(zhì)中點(diǎn)源激發(fā)的理論地震圖。隨后,1959年實(shí)現(xiàn)了全球模型運(yùn)動(dòng)方程的數(shù)值積分計(jì)算,1959—1960年間,Dorman等[5]發(fā)表了第一個(gè)基于數(shù)值計(jì)算給出的包含低速層的海洋地幔模型,1960年利用數(shù)值方法進(jìn)行了震中定位,1964年進(jìn)行了復(fù)雜地球模型中面波的激發(fā)模式和地球介質(zhì)Q值的反演,1965年和1969年利用快速傅里葉變換進(jìn)行了傅立葉振幅譜和功率譜的數(shù)值計(jì)算。1962—1969年間,數(shù)值計(jì)算方法提高了地震數(shù)據(jù)處理能力,數(shù)字濾波、多維傅里葉變換、F-K變換(F指頻率,K指波數(shù))、時(shí)間域譜分析等技術(shù)應(yīng)用于地震數(shù)據(jù)處理。此后,理論地震學(xué)中的正演計(jì)算、反演計(jì)算、數(shù)據(jù)處理越來越多地依賴于數(shù)值方法[1,6-7]。
計(jì)算地震學(xué)最根本的思想是如何將連續(xù)的地球介質(zhì)轉(zhuǎn)化為離散的模型,然后在計(jì)算機(jī)上對(duì)離散模型實(shí)現(xiàn)數(shù)值求解。通常將真實(shí)地球介質(zhì)在空間上劃分為細(xì)小的單元,在單元的節(jié)點(diǎn)上通過一定的數(shù)值算法求解方程,得到連續(xù)模型的近似數(shù)值解。一般包括以下一些步驟:
(1)預(yù)處理。
①確定物理問題的控制方程;
②確定物理問題對(duì)應(yīng)的真實(shí)地球介質(zhì)模型的空間范圍;
③將介質(zhì)在空間上離散為單元;
④定義相應(yīng)的初始條件和邊界條件。
(2)數(shù)值計(jì)算開始,迭代(或遞推)求解。
(3)計(jì)算結(jié)束,利用后處理方法分析計(jì)算結(jié)果,并用可視化技術(shù)表示結(jié)果。
計(jì)算地震學(xué)中常用的數(shù)值方法可以分為域類方法(domain method)和邊界類方法(boundary method)兩大類。
域類方法是在整個(gè)計(jì)算模型空間區(qū)域上進(jìn)行離散,用分布在整個(gè)區(qū)域上的有限個(gè)節(jié)點(diǎn)上函數(shù)的近似值來代替連續(xù)問題的解。常用的方法包括:
(1)有限差分方法;
(2)有限單元方法;
(3)偽譜方法;
(4)譜元方法。
有限差分法和有限單元法是應(yīng)用于計(jì)算地震學(xué)較早的傳統(tǒng)數(shù)值方法。偽譜法和譜元法則是近些年引入的新方法。偽譜方法(pseudo-spectral method)是譜方法的一種,和有限差分法類似,它只在離散節(jié)點(diǎn)上近似函數(shù)值,節(jié)點(diǎn)上函數(shù)的微分計(jì)算通過快速傅里葉(Fourier)變換或切比雪夫(Chebyshev)變換實(shí)現(xiàn),因此,它具有精度較高的特點(diǎn)。譜元法(spectral-element method)是一種高階有限單元方法,它選取以正交多項(xiàng)式表示的基函數(shù),它吸取了譜方法精度高的優(yōu)點(diǎn),提高了計(jì)算精度和效率。
此外,還發(fā)展了基于以上基本方法的混合方法,以提高計(jì)算精度和效率。
不同的域類方法劃分網(wǎng)格的方式不同,但是都能夠較好地處理空間上變化較大的非均勻介質(zhì)模型。
邊界方法主要包括邊界元、邊界積分方程等方法,它將計(jì)算模型空間的邊界離散化,在邊界的節(jié)點(diǎn)上求解。它是針對(duì)域類方法占用計(jì)算機(jī)內(nèi)存資源過多的缺點(diǎn)而發(fā)展起來的一種求解偏微分方程的數(shù)值方法,它的最大優(yōu)點(diǎn)是降維,將三維問題降維為二維問題,二維問題降維為一維問題,只在求解區(qū)域的邊界進(jìn)行離散就能求得整個(gè)模型區(qū)域的解。
計(jì)算地震學(xué)主要應(yīng)用于理論地震學(xué)的正演和反演兩大領(lǐng)域。對(duì)于給定的地球介質(zhì)模型和特定的地震學(xué)問題,正演是通過數(shù)值方法求解方程,給出預(yù)測的理論解,有助于理解地震學(xué)問題的物理過程。反演則是建立在正演基礎(chǔ)之上的數(shù)值計(jì)算,通過正演結(jié)果和地震觀測數(shù)據(jù)之間的擬合,估計(jì)地球介質(zhì)的地震波速度結(jié)構(gòu)、介質(zhì)物性參數(shù)和震源參數(shù)等,反演過程中的正演計(jì)算和反演算法都需要利用數(shù)值方法實(shí)現(xiàn)。在實(shí)際應(yīng)用方面,目前計(jì)算地震學(xué)已經(jīng)廣泛用于全球、區(qū)域和局部不同尺度的地球介質(zhì)速度結(jié)構(gòu)和物性參數(shù)成像以及地震震源過程研究,為地球動(dòng)力學(xué)研究提供證據(jù),為勘探地震學(xué)和強(qiáng)震地面運(yùn)動(dòng)預(yù)測提供新的手段。此外,數(shù)值方法也廣泛應(yīng)用于地震觀測數(shù)據(jù)的分析和處理中[8-11]。
王彥賓等[12]基于傅里葉偽譜法(Fourier pseudospectral method,PSM),發(fā)展了模擬具有任意橫向非均勻結(jié)構(gòu)的全地球模型中SH波傳播的一種數(shù)值計(jì)算方法。圖1給出了他們計(jì)算的一個(gè)模擬的深源地震(震源深度600km)激發(fā)出的SH波在地球中的傳播情況。計(jì)算中采用的地球地震波速結(jié)構(gòu)模型是PREM模型(1981由美國A.M.Dziewonski和D.L.Anderson提出的初步參考地球模型)。圖1展示出震源激發(fā)地震波后4個(gè)時(shí)刻的SH波的波場快照?qǐng)D,從圖中可看出,不同時(shí)刻SH波在各個(gè)界面上的反射、透射和繞射情況,圖中還標(biāo)出了在各個(gè)時(shí)刻SH波的多個(gè)震相在地球內(nèi)部到達(dá)的位置。
圖1 用PREM模型計(jì)算得到的4個(gè)時(shí)刻的SH波位移波場快照?qǐng)D
計(jì)算地震學(xué)的發(fā)展是由以下幾個(gè)因素推動(dòng)的:地震觀測數(shù)據(jù)不斷增多,人們認(rèn)識(shí)到數(shù)據(jù)所反映的地球介質(zhì)以及地震過程的復(fù)雜性,希望通過復(fù)雜介質(zhì)模型以及地震震源過程的正演和反演修正地震學(xué)的理論,并得到更符合實(shí)際的地球模型??焖侔l(fā)展的計(jì)算機(jī)硬件條件和數(shù)值算法為此提供了必不可少的條件[13]。
目前,計(jì)算地震學(xué)已經(jīng)能夠處理三維復(fù)雜模型和模擬地震破裂的動(dòng)態(tài)過程??梢灶A(yù)見,隨著計(jì)算機(jī)性能的不斷提高和數(shù)值方法的發(fā)展,計(jì)算地震學(xué)將發(fā)揮更重要的作用,通過更接近實(shí)際的地球模型的數(shù)值計(jì)算,提高我們對(duì)地震這一自然現(xiàn)象本質(zhì)的認(rèn)識(shí),減輕地震災(zāi)害對(duì)人類的威脅,提高我們對(duì)地下自然資源的勘探能力。
作者電子信箱 王彥賓#ybwang@pku.edu.cn$
[1]Ben-Menahem A.A concise history of mainstream seismology:origins,legacy and perspectives.Bull.Seis.Soc.Amer.,1995,85(4):1202-1225
[2]Thomson W T.Transmission of elastic waves through a stratified solid.J.Appl.Phys.,1950,21:89-93
[3]Haskell N.The dispersion of surface waves in multilayered media.Bull.Seis.Soc.Amer.,1953,43:17-34
[4]Pekeris C L.The seismic surface pulse.Proc.Nat.Acad.Sci.USA,1955,41:469-480
[5]Dorman J,Ewing M,Oliver J.Study of shear velocity distribution in the upper mantle by mantle Rayleigh waves.Bull.Seis.Soc.Amer.,1960,50(13):87-115
[6]Engdahl E R,Taggart J,Lobdell J L,et al.Computational methods.Bull.Seis.Soc.Amer.,1968,58(4):1339-1344
[7]Harkrider D G.The early years of computational seismology at Caltech.Bull.Seis.Soc.Amer.,1988,78(6):2105-2109
[8]Herrmann R B.Digital processing of regional network data.Bull.Seis.Soc.Amer.,1982,72(6B):S261-S276
[9]Allen R.Strong motion prediction using mathematical modeling techniques.Bull.Seis.Soc.Amer.,1982,72(6B):S29-S41
[10]Engdahl E R,Peterson J,Orsini N A.Global digital networks:current status and future directions.Bull.Seis.Soc.Amer.,1982,72(6B):S243-S259
[11]Aki K.Strong motion prediction using mathematical modeling techniques.Bull.Seis.Soc.Amer.,1982,72(6B):S29-S41
[12]王彥賓,Takenaka H.利用偽譜法模擬橫向非均勻全球SH波場.中國科學(xué):地球科學(xué),2012,42(1):140-148
[13]Tromp J.A new era in computational global seismology.Seism.Res.Lett.,2001,72(6):640-641
P315.01;
: A;
10.3969/j.issn.0235-4975.2014.02.007
2013-12-18。