地震學(xué)百科知識(shí)（八）
——計(jì)算地震學(xué)＊

地震科普

地震學(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ù)值近似解。

1 歷史

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］。

2 基本思想

計(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é)果。

3 常用方法

計(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ū)域的解。

4 應(yīng)用

計(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

5 前景

計(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。