地震波数值模拟方法研究综述

地震波数值模拟方法研究综述 在地学领域，对于许多地球物理问题，人们已经得到了它应 遵循的基本方程(常微分方程或偏微分方程)和相应的定解条 件，但能用解析方法求得精确解的只是少数方程性质比较简 单，且几何形状相当规则的问题。对于大多数问题，由于方 程的非线性性质，或由于求解区域的几何形状比较复杂，贝y 不能得到解析解。这类问题的解决通常有两种途径。一是引 入简化假设，将方程和几何边界简化为能够处理的情况，从 而得到问题在简化状态下的解答。 但这种方法只是在有限的 情况下是可行的，过多的简化可能导致很大的误差甚至错误 的解答。因此人们多年来寻找和发展了另一种求解方法 一- 数值模拟方法。 地震数值模拟(SeismicNumericalModeling) 是地震勘探和 地震学的基础，同时也是地震反演的基础。所谓地震数值模 拟，就是在假定地下介质结构模型和相应的物理参数已知的 情况下，模拟研究地震波在地下各种介质中的传播规律，并 计算在地面或地下各观测点所观测到的数值地震记录的一 种地震模拟方法。地震波场数值模拟是研究复杂地区地震资 料采集、处理和解释的有效辅助手段，这种地震数值模拟方 法已经在地震勘探和天然地震领域中得到广泛应用。 地震数值模拟的发展非常迅速， 现在已经有各种各样的 地震数值模拟方法在地震勘探和地震学中得到广泛而有效 的应用。这些地震波场数值模拟方法可以归纳为三大类，即 几何射线法、积分方程法和波动方程法。波动方程数值模拟 方法实质上是求解地震波动方程， 因此模拟的地震波场包含 了地震波传播的所有信息， 但其计算速度相对于几何射线法 要慢。几何射线法也就是射线追踪法， 属于几何地震学方法， 由于它将地震波波动理论简化为射线理论， 主要考虑的是地 震波传播的运动学特征，缺少地震波的动力学信息，因此该 方法计算速度快。因为波动方程模拟包含了丰富的波动信 息，为研究地震波的传播机理和复杂地层的解释提供了更多 的佐证，所以波动方程数值模拟方法一直在地震模拟中占有 重要地位。

1地震波数值模拟的理论基础 地震波数值模拟是在已知地下介质结构的情况下， 研究地震 波在地下各种介质中传播规律的一种地震模拟方法， 其理论 基础就是表征地震波在地下各种介质中传播的地震波传播 理论。上述三类地震波数值模拟方法相应的地震波传播理论 的数学物理表达方式不尽相同。 射线追踪法是建立在以射线 理论为基础的波动方程高频近似理论基础上的， 其数学表形 式为程函方程和传输方程。积分方程法是建立在以惠更斯原 理为基础的波叠加原理基础上的， 其数学表达形式为波动方 程的格林函数域积分方程表达式和边界积分方程表达式 动方程数值解法是建立在以弹性或粘弹性理论和牛顿力学 为基础的双曲型偏微分方程一波动方程的理论基础上的。 由 于地下介质性质不同，其相应的地震波传播方程也不同。 由于地震波波动方程在复杂介质中地震波传播研究的 广泛适应性及地震波方程数值解法在地震波数值模拟中应 用的广泛性和有效性， 本文将重点研究地震波波动方程数值 解法，同时对几何射线法和积分方程法也作适当的讨论。

2地震波数值模拟的内容及特点 2 . 1地震波数值模拟的内容 地震波数值模拟是以地震波传播理论为基础的。 描述地 震波在各种介质中传播的波动方程是一个变系数的偏微分 方程。地震波波动方程的定解问题 (即正演问题)包括微分算 子、算子系数、震源、初始条件和边界条件等。地震波正演 过程中，求解微分方程可以计算出系统中表示状态的参数随 时间的变化。地震波的正演过程数学上可描述为 d=A(m) 式中d为合成地震数据向量； 以为正演算子； m为模型向量。 d的精度受m的离散化精度和正演算子以计 算精度的影响。 地震波波动方程正演问题的内容主要包括： (1) 地震波数值模拟的基本原理； (2) 地震波数值模拟的算法； (3) 程序设计及其质量，它主要受计算精度、计 算效率和计算稳定性的影响。 2. 2模型的离散化 研究的目的不同，构成地球物理模型的物理量 也不同。对于均匀各向同性介质中的声波方程而言，地球物 理模型小可以表示为 m=(p,v P,)；而对于均匀各向同性介质 中的弹性波方程而言，其地球物理模型则可以表示为 m — m(p , vp,vs)。 地球物理模型的离散化是通过对模型的空间剖分实现 的。地球物理模型的空间剖分方法目前主要有两种，即正交 网格剖分和非正交网格剖分。 所谓的正交网格就是在平面上 是矩形网格，而非正交网格在平面上是三角形网格和不规则 四边形网格。对于地下介质进行非正交网格剖分可以差分考 虑地下介质分布的几何形状，并且不受边界几何形态的限 制。基于这一点，非正交网格数值模拟方法要优于正交网格 的数值模拟方法。为了准确刻画地下介质物理性质的空间变 化，网格剖分必须要足够精细，但是模型剖分得越细，空间 网格点的数目就越多，这必然会占用大量的计算内存，加大 计算量，降低效率，从而增加计算的成本，同时引起误差的 积累。因此，模型离散化时必须考虑数值模拟的分辨率 （或网 格大小）和计算成本。 2 . 3波动方程的离散化 模型空间的网格化必然带来波场的网格化。 由于这种网 格化把一个连续的地震波动问题转化成一个离散的地震波 动问题，因此必然涉及到波场逼近， 并且在空间网格化以后， 尽可能以较小的逼近误差表示离散波场的空间微分。 有限差 分法通过有限差分算子将波动方程离散化，以差分代替微 分，将微分方程问题转化为代数方程问题，然后求解相关的 线性代数方程组以获得微分方程问题的数值解。 差分算子是 一个空间局部的算子，在空间域具有较高的分辨率，可以较 好地适应剧烈变化的地下介质。但是在频率域中，有限差分 算子的分辨率就非常低了。 算法的稳定性和收敛性受空间采 样率和时间采样率的影响，但算法的速度较快。 基于变分原理和网格插值的有限元法比较适合几何条 件和物理条件都较复杂的问题。 但是算法复杂，计算速度慢， 一般要求插值基函数是分段线性函数，不具有正交性，算子 也是空间局部算子，空间分辨率高，但是频率域中分辨率却 很低。 另外一种逼近空间微分的方法是伪谱法， 它是利用傅立 叶变换将波场函数表示为傅立叶级数的展开形式， 将时间域 的波动方程在频率域中求解。 伪谱法对微分算子的逼近程度 可以达到尼奎斯特频率，并且收敛速度快。但由于傅立叶变 换是基于整个空间域的， 改变空间中任一点的值就会改变频 率域中的所有值，因此每一个点的微分结果都要受到计算域 中其它点的影响。实际上，求导运算应该是一种局部运算， 对于空间物性剧烈变化的情形显然有其局限性。 与有限差分 法和有限元法相比而言，伪谱法在频率域中的分辨率高，而 在时间域的分辨率却相对较低。 将有限元法和伪谱法相结合就产生了现在流行的地震 波数值模拟方法——谱元法。 如果地震波场具有规则的特征， 即地下介质是均匀分布 的，那么上面几种算法都是适合的。事实上，由于地下介质 的分布具有高度非均匀性， 且这种非均匀性发生在非常大的 尺度范围内（其尺度从岩石粒度到全球球谐函数的最低阶 ）， 地面接收到的地震波场不仅包含反射和折射信息， 还包含了 散射信息（介质的奇异性信息）。介质的这种多尺度非均匀性 通过地震波动方程可以映射到地面接收的地震记录中， 其中 介质的物理性质变化通过波动方程的系数体现出来。

3 .20世纪90年代以来地震波数值模拟新进展 随着地震波理论在天然地震和勘探地震中的应用，地震 模拟技术应运而生，并随着波动理论和计算机技术的发展， 地震数值模拟技术自20世纪90年代以来得到了飞速发展， 到目前为止形成了射线追踪法、有限差分法、有限兀法、伪 谱法、谱元法和积分方程法等各种现代数值模拟技术。 有限差分法是偏微分方程的主要数值解法之一。在各种 数值模拟方法中，最早出现的就是有限差分法。 Alterman

等［11首先将有限差分法应用于层状介质弹性波传播的数值 模拟中，之后许多研究人员对该方法作了深入的探讨。 Tal

—EzerH等［21研究了线性粘弹性介质中地震波传播的数值 模拟方法；Robertsson 等［3］给出了粘弹性波有限差分模拟 方法；Carcione和Helle［41提出了孔隙粘弹性介质中地震 波传播的交错网格有限差分模拟方法。 一般的有限差分地震模拟方法是基于笛卡儿坐标系中 的规则网格，在模拟复杂地质构造和复杂地质体的复杂界面 时，必然会出现弯曲边界，在这种边界上必然会引起人为的 虚假绕射波，为了减弱这种虚假绕射， 就必须采用精细网格， 而这不仅会导致存储量的增加和计算量的加大， 而且会带来 误差的积累。为此，人们发展了基于可变网格和不规则网格 的地震波数值模拟方法。 Jastram 等口 ］提出了垂直间距可变网格的弹性波模拟 方法；Oprsal等…］提出了非均匀介质弹性波的矩形不规则 网格有限差分模拟方法； Nordstrom 等［7］提出了曲线坐标 下变形网格高阶有限差分法地震波数值模拟方法。 随着地表 复杂区地震勘探的发展，起伏地表地震波数值模拟技术受到 了越来越多的关注和重视。 董良国等［8 ' 9］运用高阶差分法，