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

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

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

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

前言

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

1 FDTD方法

1.1 FDTD算法的基本格式

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語

FDTD方法是解決地震波場數(shù)值模擬問題的一種有效方法。理論上，選取的空間步長與時間步長越短，模擬的效果越好。但在另一方面，這樣做會增加算法的復(fù)雜度。為了解決這一問題，我們可以通過計算數(shù)值色散因子選取合適的步長。事實上，在一些更深入、復(fù)雜的研究中采用的FDTD算法，具有更高的維度、多變的差分方向式與網(wǎng)格形狀。這些變化都是為了適應(yīng)研究對象的特點。關(guān)于地震波場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ì)大學(xué)學(xué)報，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]朱金明，王麗燕.地震波走時有限差分計算[J].地球物理學(xué)報，1992，35（1）：86-92.