地震波数值模拟方法研究综述
地震波数值模拟方法研究综述
在地学领域，对于许多地球物理问题，人们已经得到了它应遵循的基本方程(常微分方程或偏微分方程)和相应的定解条件，但能用解析方法求得精确解的只是少数方程性质比较简单，且几何形状相当规则的问题。对于大多数问题，由于方程的非线性性质，或由于求解区域的几何形状比较复杂，贝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］运用高阶差分法，
通过交错网格技术，对一阶速度一应力弹性波方程进行了数值求解，并进行了算法稳定性分析。孙若昧等［1们采用多阶振型和有限差分联立的混合法，模拟了1976年唐山地震所
引起的北京西集一郎府地区的剪切波运动。Yang等［11］提出一种基于矢量和矩阵在二维各向异性介质中应用的快速有限差分方法，并碍到了稳定性方程。殷文等n胡采用25点
优化差分算子，再根据最优化理论求取的优化系数，建立了频率空间域中弹性波方程的差分格式，有效地克服了常规差
分算子的数值频散。王一博等［1朝采用紧集正交小波基对空间域进行多尺度离散，采用二阶精度有限差分算子对时域离散，推导得到了多尺度有限差分方法正演模拟的递推公式，并实现了相应的波传播过程的数值模拟。
有限元法也是偏微分方程的数值解法之一。Lysmer和Drake［14］最早将有限元法应用于地震波数值模拟。Seron 等［15d63 给出了弹性波传播的有限元数值模拟方法；Padovani等‘”］研究了地震波数值模拟的低阶和高阶有限元方法；Sarma等n81给出了弹性波传播有限元数值模拟的无反射边界条件。张美根等口93研究了各向异性弹性波
有限元正演系统的精度和效率问题，提出了一种透射加衰减
的组合人工边界方案（吸收边界条件）。杨顶辉等口叩基于双
相各向异性介质模型，推导了双相各向异性介质中弹性波传
播的动力学方程及其Galerkin变分方程和有限元运动方程，对双
相PTL介质和双相各向同性介质中的弹性波传播进行了数值模拟。
伪谱法是偏微分方程的另一种数值解法，它最早由
Kreiss和Oliger［z妇提出。进入上世纪90年代之后，该方法有了飞跃式的发展。石玉梅［223给出了流体饱和多孔隙介质中弹性波传播数值模拟的伪谱法。张文生等‘ 231用伪谱法进行了二维横向各向同性介质波动方程的正演模拟，特别是对边界吸收问题作了有效的处理。Takashi等阻］首次针对
伪谱法提出了反周期扩展边界方法；Takenaka和王彦宾等
瞳“ 263利用伪谱法分别计算了球对称全球模型和具有垂向速度梯度的沉积盆地模型中地震波的传播问题。之后不久，王彦宾和Takenaka等乜71利用不连续网格傅立叶伪谱多域方法模拟了区域地球模型中弹性波的传播，接下来王彦宾等利用伪谱法模拟了二维柱坐标下全地球模型中地震波的传播。赵志新等［2盯给出了非均匀介质中地震波传播数值模拟的错格实数傅立叶伪谱法；赵景霞等口们利用多块映射和
超限插值技术将直角坐标系下的曲界面变换成曲坐标系下的规则界面，在曲坐标系下利用伪谱法模拟波场；张丽琴等给出了三维离散介质中地震波传播数值模拟的伪谱法。
谱元法是一种新的波动方程数值解法，它首先由Patera[3胡在流体动力学研究中提出。PriOlO和Seriani[3 刀首先将谱元法用于地震波传播的数值模拟。
Komatitsch[341 采用拉格朗日插值多项式的谱元法对大尺
度地质模型进行了三维地震波场数值模拟，首次得到了地震
应用中的对角质量矩阵，从而得到非常简单的显式时间差分
方案和有效的并行实现方案。Chaljub[353 首先将谱元法用
于全球地震波传播的数值模拟；Capdeville 等[36删将谱元
法与基于球谐函数展开后正常振型求和的形式解方法相结合数值模拟了全球尺度的弹性波传播；Komatitsch 等‘ 38]
将谱元法推广到了二维介质中三角形网格单元的弹性波数值模拟；Komatitsch 和Tromp[391 在谱元法中采用了适用于弹性动力学方程变分公式的完美匹配层吸收边界条件。谱元法的这些最新发展趋势表明该方法是地震波数值模拟的重要工具，在将来的应用中具有巨大的潜力，但在方法论上还不十分完善。
积分方程法是建立在波动方程的积分表达基础上的，其理论基础是惠更斯原理。Bakamjian[4 阳给出了三维地震波传播数值模拟的边界积分方程法；符力耘和牟永光[413提出
了弹性波正演模拟的边界元法。之后，符力耘等n2]提出了非线性Fredh01m 积分方程的正演问题。接着符力耘[43] 给出了含起伏地表的广义Lipmann —Schwinger 积分方程
的数值模拟方法。
射线追踪法是建立在波动方程的高频渐近基础上的地震波数值模拟方法。这种方法实际只是计算了最奇异部分的解，即走时
和振幅函数的特征曲线，它们分别是程函方程和传输方程的解。这种方法计算效率高，但是一些复杂的本构方程由于积分方程法和射线追踪法不满足假设条件而限制了它们的应用。陈九辉和刘启元[4幻提出了一种基于Maslov理论的横向非均匀三维介质接收函数的计算方法，不仅在保证追踪精度的同时提高了射线追踪的效率，而且可完全避免传统射线追踪过程中的射线编码问题。王辉和常旭[4明为了克服最小走时射线追踪方法存在的问题，对该方法计算过程中的关键步骤进行了改进。改进的最小走时射线追
踪方法具有精度高、速度快的特点。赵爱华等[463提出了一种基予惠更斯原理和费马原理求取地震波走时及其反射波射线路径的新方法；高亮等H73以及秦孟兆和陈景波[48]
把辛几何算法引入射线追踪问题，不仅提高了效率，而且克服了焦散问题。
除了以上常规的解析和数值算法以外，一些研究者还在
把科学和工程的其它领域中出现的新技术引入复杂介质中的地震波数值模拟研究方面作了有益的探索。王真理等[49] 把基于区域分解的格子自动机地震波传播模拟并行化算法，应用于实际地质模型中平面波反射地震响应的计算，极大地提高了计算效率。刘劲松等[5叼提出一种新的用于地震波场
模拟的变网格边长声格固体模型；罗明秋等［513给出了哈密
顿体系下地震波传播的表述和用于波场计算的一些辛格式。数值
结果表明，辛算法是一类非常好的波场模拟数值算法。
陈凌等［523从Kaiser的声子波理论出发，给出了通过分别引入点源波形的复时间函数和点源虚时间坐标来构成声子波的两种解释，并对点源模型的合成地震图和实际复杂模型的地震波资料进行了时一空域的声子波变换。刘礼农等口33 利用单程波方法对复杂构造中的地震波传播进行了数值模拟。谢桂生等口41导出了一种由单程波方程利用反射/透射算子的可分表示方法模拟复杂介质中一次反射地震波的数值算法，该算法具有足够高的模拟精度，且计算效率成倍地提高。孙建国等D53利用有限正弦和有限余弦变换模拟地震波场；贾晓峰等‘驰1利用无单元法进行了地震波动方程数值模拟和成像的研究。
4地震波数值模拟方法的评述
基于高频近似的射线追踪是一类研究均匀介质中地震波传播，特别是在不规则界面上的反射波行为的重要工具。随着地震勘探技术的发展，薪的射线追踪技术也不断涌现，以满足大的数据处理（如三维数据）和较高精度要求下对复杂地质体研究的需要。这些技术主要研究焦点是如何精确地划分地质体，如何实现旅行时的快速准确计算以及对已有方法
的改良等。射线法的主要优点是概念明确，显示直观，运算方便，适应性强；其缺陷是应用有一定限制条件，对于复杂构造进
行两点三维射线追踪往往比较麻烦，为了计算波沿射线的走时和波沿射线的振幅变化，首先都必须知道波的传播
路径。近几年射线法取得了很多新的发展，最近提出的高斯射线束法、傍轴射线法、马斯罗夫渐进理论等就是明显的例证。另外，射线法和波动方程理论相结合解决地震学和地震勘探中的各种正、反演问题，也是目前发展的一种趋势。
从惠更斯原理出发，根据指定的边界条件或震源分布建立波动方程积分形式解的方法称为积分解法。这时惠更斯原理一般采用两种形式，即体积分方程和边界积分方程，两者分别有各自的高频渐近或者说射线追踪方法，这是一种在地震波模拟和成像中应用广泛的方法。高频近似方法一般不考
虑全波场，但从另一方面讲，它计算效率比较高，所以在地震成像中占有很重要的地位。
波动方程法是研究地震波场特点的最根本方法，虽然费时，易引进干扰波，但是其波场齐全，信息丰富，因此对于研究复杂条件下的各种波场最为有效，具有广阔的发展前景。随着计算机运算能力的提高和各种新方法的诞生，波动方程法会越来越显示出它的优越性。有限差分法的主要优点
是计算速度快，占用内存小；缺点是精度低，难以克服频散效应，而要解决频散问题，须加密数值计算的网格，这势必
导致计算量增加、效率下降，所以仅适合于相对较简单的地质模型。伪谱法由于利用空闯的全部信息对波场函数进行三角函数插
值，所以能更加精确地模拟地震波的传播规律。同时，利用快速傅立叶变换(FFT)进行计算，还可以提高运算效率。有限元法的主要优点是适宜予模拟任意形态地质体，可以任意三角形单元逼近地层界面，保证复杂地层形态模拟
的逼真性；缺点是占用内存和运算量均较大，使得计算速度过慢。谱元法本质上是一种基于Galerkin 变分法的近似求
解以积分形式表示的偏微分方法的离散方法，它结合了伪谱法和有限元法的优点，既具有伪谱法的高精度和快速收敛性，又具有有限元法的几何灵活性。
21世纪出现的新方法，如单程波方法、波传播的哈密顿系统描述、可分表示方法和无单元法等具有比传统方法更高的模拟精度，昭示着地震波数值模拟方法发展的美好未来。
5地震波数值模拟方法应用前景及展望
地震波传播数值模拟就是模拟波在地下介质中的传播过程；并且获得给定介质模型在地表接收器记录的理论地震图。地震波数值模拟和合成地震图主要用于：
(1)研究地震
波在各种复杂介质中的传播行为，无论对勘
探地震剖面还是天然地震台站记录，都可以作为数据解释的
强有力的辅助工具
(2)为检验各种地震成像方法的效果而提供模型合成记录，
比如在工业界已经成为检验成像及速度估计标准的Marmousi和SEG /EAGE盐丘模型记录。
(3)为大多数反演方法提供理论基础。
(4)为勘探地震数据采集的方案设计及评价提供依据。
地震波理论是固体地球物理学研究的重要基础，地震波
研究领域的任何实质性进展都会促进固体地球物理学的发展。在过去的几十年里，中国地球物理学家在该领域做了很多有意义的研究工作，其中不乏创新性的理论工作。当前地震波数值模拟研究领域的重要课题包括：
(1)复杂地球介质中地震波激发与传播理论。
(2)高效而精确地计算三维全球和区域尺度介质中地震波传
播的数值方法。
(3)利用先进的地震波数值模拟方法，开展地震与强地面
运动的数值模拟研究，为精细的地震危险分析与预测奠定基
础。