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

地震波数值模拟方法研究综述在地学领域，对于许多地球物理问题，人们已经得到了它应遵循的基本方程(常微分方程或偏微分方程)和相应的定解条件，但能用解析方法求得精确解的只是少数方程性质比较简单，且几何形状相当规则的问题。对于大多数问题，由于方程的非线性性质，或由于求解区域的几何形状比较复杂，则不能得到解析解。这类问题的解决通常有两种途径。一是引入简化假设，将方程和几何边界简化为能够处理的情况，从而得到问题在简化状态下的解答。但这种方法只是在有限的情况下是可行的，过多的简化可能导致很大的误差甚至错误的解答。因此人们多年来寻找和发展了另一种求解方法——数值模拟方法。地震数值模拟(SeismicNumericalModeling)是地震勘探和地震学的基础，同时也是地震反演的基础。所谓地震数值模拟，就是在假定地下介质结构模型和相应的物理参数已知的情况下，模拟研究地震波在地下各种介质中的传播规律，并计算在地面或地下各观测点所观测到的数值地震记录的一种地震模拟方法。地震波场数值模拟是研究复杂地区地震资料采集、处理和解释的有效辅助手段，这种地震数值模拟方法已经在地震勘探和天然地震领域中得到广泛应用。地震数值模拟的发展非常迅速，现在已经有各种各样的地震数值模拟方法在地震勘探和地震学中得到广泛而有效的应用。这些地震波场数值模拟方法可以归纳为三大类，即几何射线法、积分方程法和波动方程法。波动方程数值模拟方法实质上是求解地震波动方程，因此模拟的地震波场包含了地震波传播的所有信息，但其计算速度相对于几何射线法要慢。几何射线法也就是射线追踪法，属于几何地震学方法，由于它将地震波波动理论简化为射线理论，主要考虑的是地震波传播的运动学特征，缺少地震波的动力学信息，因此该方法计算速度快。因为波动方程模拟包含了丰富的波动信息，为研究地震波的传播机理和复杂地层的解释提供了更多的佐证，所以波动方程数值模拟方法一直在地震模拟中占有重要地位。1地震波数值模拟的理论基础地震波数值模拟是在已知地下介质结构的情况下，研究地震波在地下各种介质中传播规律的一种地震模拟方法，其理论基础就是表征地震波在地下各种介质中传播的地震波传播理论。上述三类地震波数值模拟方法相应的地震波传播理论的数学物理表达方式不尽相同。射线追踪法是建立在以射线理论为基础的波动方程高频近似理论基础上的，其数学表形式为程函方程和传输方程。积分方程法是建立在以惠更斯原理为基础的波叠加原理基础上的，其数学表达形式为波动方程的格林函数域积分方程表达式和边界积分方程表达式。波动方程数值解法是建立在以弹性或粘弹性理论和牛顿力学为基础的双曲型偏微分方程一波动方程的理论基础上的。由于地下介质性质不同，其相应的地震波传播方程也不同。由于地震波波动方程在复杂介质中地震波传播研究的广泛适应性及地震波方程数值解法在地震波数值模拟中应用的广泛性和有效性，本文将重点研究地震波波动方程数值解法，同时对几何射线法和积分方程法也作适当的讨论。2地震波数值模拟的内容及特点2．1地震波数值模拟的内容地震波数值模拟是以地震波传播理论为基础的。描述地震波在各种介质中传播的波动方程是一个变系数的偏微分方程。地震波波动方程的定解问题(即正演问题)包括微分算子、算子系数、震源、初始条件和边界条件等。地震波正演过程中，求解微分方程可以计算出系统中表示状态的参数随时间的变化。地震波的正演过程数学上可描述为d=A(m)式中d为合成地震数据向量；以为正演算子；m为模型向量。d的精度受m的离散化精度和正演算子以计算精度的影响。地震波波动方程正演问题的内容主要包括：(1)地震波数值模拟的基本原理；(2)地震波数值模拟的算法；(3)程序设计及其质量，它主要受计算精度、计算效率和计算稳定性的影响。2．2模型的离散化研究的目的不同，构成地球物理模型的物理量也不同。对于均匀各向同性介质中的声波方程而言，地球物理模型小可以表示为m=(p,vp,)；而对于均匀各向同性介质中的弹性波方程而言，其地球物理模型则可以表示为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胡在流体动力学研究中提出。Pri010和Seriani[3刀首先将谱元法用于地震波传播的数值模拟。Komatitsch[341采用拉格朗日插值多项式的谱元法对大尺度地质模型进行了三维地震波场数值模拟，首次得到了地震应用中的对角质量矩阵，从而得到非常简单的显式时间差分方案和有效的并行实现方案。Chaljub[353首先将谱元法用于全球地震波传播的数值模拟；Capdeville等[36删将谱元法与基于球谐函数展开后正常振型求和的形式解方法相结合数值模拟了全球尺度的弹性波传播；Komatitsch等‘38]将谱元法推广到了二维介质中三角形网格单元的弹性波数值模拟；Komatitsch和Tromp[391在谱元法中采用了适用于弹性动力学方程变分公式