FDTD算法在模擬地震波傳播方面的應(yīng)用研究

劉冠男　張衛(wèi)東　舒亞祥

摘 要：數(shù)值模擬是一種重要的地震正演的主要手段。它能夠解決復(fù)雜地質(zhì)模型中地震波的傳播問(wèn)題。在許多數(shù)值模擬方法中，F(xiàn)DTD方法是一種非常有效的方法。文章從數(shù)學(xué)與物理學(xué)的角度討論了FDTD方法的基本原理，包括差分格式、吸收邊界條件、算法穩(wěn)定性，又利用MATLAB軟件對(duì)簡(jiǎn)單地質(zhì)模型中的地震波場(chǎng)進(jìn)行了模擬。結(jié)果顯示，利用FDTD算法模擬的地震波場(chǎng)能夠體現(xiàn)出實(shí)際地震波傳播的基本規(guī)律。本研究對(duì)地震波場(chǎng)的時(shí)域有限差分正演問(wèn)題提供了基本的思路與參考。

關(guān)鍵詞：FDTD；地震波場(chǎng)；數(shù)值模擬

前言

地質(zhì)構(gòu)造的復(fù)雜程度非常高，我們不可能用解析的形式來(lái)描述地震波的傳播問(wèn)題。為了解決這個(gè)難題，學(xué)者們引入數(shù)值模擬的方法來(lái)解決地震波在復(fù)雜介質(zhì)中的傳播問(wèn)題。地震波數(shù)值模擬方法有很多種，如有限元法（FE）、時(shí)域有限差分法（FDTD），以及傳輸矩陣法（TM）等。其中，時(shí)域有限差分法是應(yīng)用最廣泛的方法。時(shí)域有限差分法是科學(xué)家Yee在1966年提出的。Yee將含時(shí)的Maxwell方程離散化并轉(zhuǎn)化為差分格式，這就是最初的FDTD算法。在此之后，科學(xué)家們針對(duì)FDTD算法的穩(wěn)定性、邊界條件的處理方法、以及高維與高階FDTD算法進(jìn)行了研究并取得了豐富的成果。隨著數(shù)學(xué)與物理學(xué)進(jìn)一步發(fā)展，F(xiàn)DTD算法已經(jīng)突破了傳統(tǒng)的二維正方形網(wǎng)格的局限，針對(duì)不同的坐標(biāo)和區(qū)域形狀，差分網(wǎng)格的大小和形狀可以做相應(yīng)的改變。自適網(wǎng)格的差分格式逐漸成為研究的熱門。文章主要討論了時(shí)域有限差分法在模擬地震波傳播方面的應(yīng)用，包括數(shù)值解法、邊界條件與數(shù)值色散，在理論上詳細(xì)描述了FDTD方法的原理與過(guò)程，并通過(guò)MATLAB軟件編程實(shí)現(xiàn)。

1 FDTD方法

1.1 FDTD算法的基本格式

FDTD方法的原理是將微分方程離散化，再利用遞推關(guān)系求得數(shù)值解。函數(shù)的一階微分與二階微分分別可以表示為

根據(jù)式（1），我們可以將描述地震波傳播問(wèn)題的偏微分方程完全轉(zhuǎn)化為離散的形式。在深度-水平距離二維空間中，地震波傳播方程的差分格式如圖1所示。

圖1（a）給出的是地震波傳播方程的“五點(diǎn)式”差分格式。“五點(diǎn)式”差分格式僅考慮了位于兩個(gè)主要坐標(biāo)軸上的元胞，算法精度較低。為提高精度，我們將對(duì)角線方向上的元胞考慮到“五點(diǎn)式”差分格式中，得到了“九點(diǎn)式”差分格式，如圖1（b）。函數(shù)在對(duì)角線方向上的一階微分形式與二階微分形式為：

由于對(duì)角線方向上元胞之間的距離較大，它們產(chǎn)生的影響較小。我們可以引入一個(gè)影響系數(shù)a（0-0.5），用于描述不同方向上的影響比重。這樣，Laplace算符就可以表示為：

1.2 邊界條件與算法穩(wěn)定性

由于模擬區(qū)域存在邊界，模擬的地震波在到達(dá)邊界時(shí)會(huì)產(chǎn)生反射。為了去掉反射波，我們可以采用吸收邊界條件。文章在處理邊界條件方面采用了一階近似的Engquist-Majda吸收邊界條件。設(shè)任意方向l上波函數(shù)？漬滿足：

式（4）是一階近似條件下Engquist-Majda吸收邊界條件的一般形式。根據(jù)l的方向，我們可以將算符 +jkl或 -jkl作用于波函數(shù)？漬消去反射波。除吸收邊界條件外，我們還需要考慮數(shù)值色散。空間步長(zhǎng)與時(shí)間步長(zhǎng)都會(huì)影響數(shù)值色散。較大的時(shí)間或空間步長(zhǎng)會(huì)使模擬結(jié)果產(chǎn)生嚴(yán)重錯(cuò)誤。但如果步長(zhǎng)選取得過(guò)小，算法的復(fù)雜度就會(huì)過(guò)高。

2 數(shù)值模擬結(jié)果與分析

如圖3（a）所示，作者建立了一個(gè)速度場(chǎng)模型。其中，深藍(lán)色區(qū)域中地震波的傳播速度為2000m/s，青色區(qū)域中地震波的傳播速度為2400m/s，黃色區(qū)域中地震波的傳播速度為2800m/s，褐色區(qū)域中地震波的傳播速度為3200m/s。

作者通過(guò)MATLAB軟件進(jìn)行編程計(jì)算得到了地震波場(chǎng)的數(shù)值模擬結(jié)果，輸出150-550ms時(shí)刻的波場(chǎng)，如圖3（b）所示。地震波在左上角被激發(fā)后向右下方向傳播。在到達(dá)兩層的交界處時(shí)，波形發(fā)生了變化。如果進(jìn)一步仔細(xì)觀察，我們還能發(fā)現(xiàn)地層交界處的反射波與透射波。圖3（b）所示的結(jié)果說(shuō)明，F(xiàn)DTD方法在模擬地震波傳播的過(guò)程中能夠體現(xiàn)地震波傳播的基本規(guī)律。文章只是針對(duì)簡(jiǎn)單的彈性波方程簡(jiǎn)單討論了FDTD方法在地震波場(chǎng)模擬方面的應(yīng)用。在實(shí)際研究中，地質(zhì)模型的復(fù)雜度較高。這就對(duì)FDTD方法提出了更高的要求。

3 結(jié)束語(yǔ)

FDTD方法是解決地震波場(chǎng)數(shù)值模擬問(wèn)題的一種有效方法。理論上，選取的空間步長(zhǎng)與時(shí)間步長(zhǎng)越短，模擬的效果越好。但在另一方面，這樣做會(huì)增加算法的復(fù)雜度。為了解決這一問(wèn)題，我們可以通過(guò)計(jì)算數(shù)值色散因子選取合適的步長(zhǎng)。事實(shí)上，在一些更深入、復(fù)雜的研究中采用的FDTD算法，具有更高的維度、多變的差分方向式與網(wǎng)格形狀。這些變化都是為了適應(yīng)研究對(duì)象的特點(diǎn)。關(guān)于地震波場(chǎng)FDTD數(shù)值模擬的此方面研究尚在進(jìn)行中。

參考文獻(xiàn)

[1]C·B·Zhao. Computational simulation of wave propagation problems in infinite domains [J].SCIENCE CHINA Physics， Mechanics & Astronomy， 2010，53（8）：1397-1407.

[2]馮英杰，等.地震波有限差分法綜述[J].地球物理學(xué)進(jìn)展，2007，22（2）：487-491.

[3]陸偉民，許哲明.地震波的計(jì)算機(jī)生成方法[J].同濟(jì)大學(xué)學(xué)報(bào)，1984，2：110-124.

[4]裴正林，牟永光.地震波傳播數(shù)值模擬[J].地球物理學(xué)進(jìn)展，2004，19（4）：933-941.

[5]孫書榮，凡友華.地震波數(shù)值模擬技術(shù)發(fā)展現(xiàn)狀[J].油氣地球物理， 2009，7（1）：18-22.

[6]朱金明，王麗燕.地震波走時(shí)有限差分計(jì)算[J].地球物理學(xué)報(bào)，1992，35（1）：86-92.