FDTD算法在模擬地震波傳播方面的應用研究

劉冠男　張衛東　舒亞祥

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

關鍵詞：FDTD；地震波場；數值模擬

前言

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

1 FDTD方法

1.1 FDTD算法的基本格式

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

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

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

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

1.2 邊界條件與算法穩定性

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

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

2 數值模擬結果與分析

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

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

3 結束語

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

參考文獻

[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].地球物理學進展，2007，22（2）：487-491.

[3]陸偉民，許哲明.地震波的計算機生成方法[J].同濟大學學報，1984，2：110-124.

[4]裴正林，牟永光.地震波傳播數值模擬[J].地球物理學進展，2004，19（4）：933-941.

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

[6]朱金明，王麗燕.地震波走時有限差分計算[J].地球物理學報，1992，35（1）：86-92.