刘冠男　张卫东　舒亚祥

摘 要：数值模拟是一种重要的地震正演的主要手段。它能够解决复杂地质模型中地震波的传播问题。在许多数值模拟方法中，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.