地震波交错网格高阶差分数值模拟研究摘要: 地震波数值模拟技术是勘探地球物理学中的重要组成部分，研究通过弹性波一阶速度——应力方程，采用交错网格高阶有限差分法实现了地震波在各向同性介质中的高精度的数值模拟，并采用完全匹配层( PML) 吸收边界来消除边界反射，可取得较好的效果。
通过模型的正演计算和复杂模型的处理结果表明，交错网格高阶有限差分法数值模拟是一种快速有效的地震波数值模拟方法。
关键词: 地震勘探; 交错网格; 有限差分; 数值模拟引言地震数值模拟是模拟地震波在介质中传播的一种数值模拟技术，随着地震波理论在天然地震和地震勘探中的应用，地震模拟技术便应运而生，并随着地震波理论和计算机技术的发展，地震数值模拟技术自20世纪60年代以来也得到了飞速发展，形成了目前具有有限差分法、有限元法、虚谱法和积分方程法等各种数值模拟方法的现代地震数值模拟技术。
有限差分法是偏微分方程的主要数值解法之一。
在各种地震数值模拟方法中，最早出现的数值模拟方法是有限差分法。
Ａlterman和Karal(1968)首先将有限差分法应用于层状介质弹性波传播的数值模拟中。
此后，Boore(1972)又将有限差分法用于非均匀介质地震波传播的模拟。
Alford等（1974）研究了声波方程有限差分法模拟的精确性。
Kelly等（1976）研究了用有限差分法制作人工合成地震记录的方法。
Virieux(1986)提出了应用速度——应力一阶方程交错网格有限差分法模拟P——SV波在非均匀介质中的传播。
交错网格方法提高了地震模拟的精度和稳定性，并消除了部分假想。
有限元法也是偏微分方程的数值解法之一。
Lysmer和Drake（1972）最早将有限元法应用于地震数值模拟。
Marfurt(1984)研究对比了模拟弹性波传播的有限差分法和有限元法的精度。
Seron等（1990,1996）给出了弹性波传播有限元模拟方法。
Padovani等（1994）研究了地震波模拟的低阶和高阶有限元法。
Sarma等（1998）给出了三维声波模拟的虚谱法。
积分方程法是建立在波动方程的积分表达式的基础上的，其理论基础是惠更斯原理。
积分方程法也是有限元法之后发展起来的一种地震数值模拟方法。
Pao 和Varatharajulu（1976）提出了弹性波散射的积分表达式。
Bennett和Mieras(1981)给出了流体目标声波散射的时间域积分方程解。
Bouchon(1987)给出了裂隙或孔洞弹性波绕射的离散波数法模拟方法。
Bouchon等（1989）研究了具有不规则界面的多层介质中波传播的边界积分方程——离散波数法。
Bakamjian（1992）给出了三维地震波传播模拟的边界积分方程法。
符力耘和牟永光（1994）提出了弹性波正演模拟的边界元法。
符力耘等（1997）提出了非线性Fredholm积分方程的正演问题。
符力耘（2003）给出了含起伏地表的广义Lipmann—Schwinger积分方程的数值模拟方法。
射线追踪方法是建立在波动方程的高频近似基础上的一种地震数值模拟方法（cerveny等，1977）。
这种方法实际只计算了最奇异部分的解，即旅行时和振幅函数的特征曲线，它们分别是程函方程和传播方程的解。
这种方法计算效率高。
但是，一些复杂的本构方程由于积分方程法和射线追踪法不满足假设条件而限制了这些方法的应用。
上述这些地震数值模拟方法各有优缺点。
对于复杂构造、复杂地质体和复杂岩性地震模拟而言，交错网格高阶有限差分法其综合性能（占内存大小、模拟精度、计算效率和并行算法实现）最好，是实用性最好的方法。
在声波方程正演的数值模拟中,由于有限的计算空间区域无形之中引入了人为边界,不可避免的需要对在数值网格边界上产生的反射或回绕能量进行合适的处理,否则这些人为产生的反射或回绕能量会在很大程度上扭曲真正的波动传播信号,使模拟剖面变得模糊不清,不利于对地层构造信息进行解释。
为了消除这些人为产生的边界反射或回绕能量,人们发展了多种方法,其中最常见的主要有以下几种:一是最简单的扩展边界法,即在需要计算的数值网格外增加一些额外的网格数目,这样可以使人为边界反射效应远离所需要的计算网格,但是这种方法带来的负面效应是所需要的网格数目大大增加,因而也大大增加了计算量,对计算机的计算速度和存储能力提出了更高的要求,所以这种方法并没有得到很好的推广;二是海绵吸收法,即在计算网格的边界区域设置一定宽度的阻尼带,利用某些衰减函数对数值模拟波场进行逐步衰减;三是反周期扩展法,即利用正反周期函数极性相反的特点消除回绕波场;四是傍轴近似法,即利用波动方程的傍轴近似条件来消除计算边界上的反射。
如何选择合适的边界吸收是一个值得研究的问题。
数值模拟基本原理各向同性介质是最基本的一种介质模型，目前地震勘探中大多都是基于这种介质模型，根据弹性介质位移，应力和应变之间的关系，可以推导出各向同性介质中的弹性波方程，在二维介质情况下为:（1）式中x V ———质点位移速度的水平分量;z V ———质点位移速度的垂直分量; xx σ———X 方向正应力; zz σ———Z 方向正应力; xz σ———切应力分量; λ，μ———拉梅系数。
一阶声波方程交错网格差分方程的建立1）一阶声波方程时间导数的2M 阶差分精度算法在计算机中进行数值计算时，需要对连续函数离散化，即对方程(1)中的微(2)(2)()x xx xzxz z zz xx x z x zz z xz x z v t x z v t x zv v t x z v v t z x v v t z x σσσσσσσρρλμλλμλμ∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂⎧=+⎪=+⎪⎪=++⎨⎪=++⎪⎪=+⎩分用有限差分来代替。
设函数f(t)单值连续且存在任意阶导数，则其Taylor 展开式为：232323111222!23!2!2()()()()...+()+()m mf f f f mm t t t t t t m t t t f t f t O t ∂∂∂∂∆∆∆∆∆∂∂∂∂+=++++∆232323111222!23!2!2-(-)()-()()...+(-)+()m mf f ff mm tt t t t t m t t t f t f t O t ∂∂∂∂∆∆∆∆∆∂∂∂∂=++∆将以上两式相减，得到2M 阶精度的时间差分近似式为：2-12-12-11(2m-1)!=12222*()*+()=(-)+2*()m m Mfm t m M tt t f t f t O t ∂∂∆∆∆+∆∑（2）Δt 为时间步长，在采用交错网格技术解声波方程时，速度场在2+tt ∆ 时刻计算，而应力场在t +Δt 时刻计算。
当M ＝1时，对方程(1)中的时间微分进行数值离散，可得2阶时间差分近似为：22(+)=(-)+[]xx xz t t t x x x z v t v t σσρ∂∂∆∆∆∂∂+22(+)=(-)+[]xz zz tt tz z xzv t v t σσρ∂∂∆∆∆∂∂+(+)=()+[(2)]x zv v xx xx xz t t t t σσλμλ∂∂∂∂∆∆++ (+)=()+[(2)]xzv v zz zz x z t t t t σσλλμ∂∂∂∂∆∆++(+)=()+[]x z v v xz xz xzt t t t σσμμ∂∂∂∂∆∆+对于时间高阶差分，计算2-12-1m m f t ∂∂时将涉及较多的时间层，需要大量的存储空间因此利用一阶弹性波方程中的速度－应力关系式，将速度对时间的任意奇数阶高阶导数用应力对空间的导数代替，而应力对时间的任意奇数阶高阶导数用速度对空间的导数代替，这样，在计算一个时间层上的速度或应力场时，只需要前一个时刻的速度或应力场以及两时间层之间的应力或速度场，不需要过多的时间层，从而节省了内存(董良国，2000)。
2）一阶弹性波方程空间导数的2N 阶差分精度算法为了提高计算精度，空间导数也要采用高阶差分近似。
在交错网格技术中，对空间变量的导数是在相应的空间变量网格点之间的半程上计算，对于空间函数f (x)，假设其存在2N ＋1阶导数，则f (x)在2-102=n x x x ±∆处的2N ＋1阶Taylor展开式为：2-122+1()(i)2+22-10002!=1()=()+()+O(),n=1,2,...,N i n N N n i i f x x f x f x x ±±∆∆∑（3）由于交错网格一阶导数2N 阶精度差分近似式可表示为：()(2N+1)2+12(+1)+12-12-100022==1={[+]-[-]}+e ()+()Nf x N N n n n N x x xn xc f x x f x x f x x O x ∂∂∆∆∆∆∆∑ （4）其中，n c 为差分系数，N e 为误差系数； 将(3)式代入(4)式中，化简得到：2+12+12+1-1(2n-1)(2n-1)(1)(1)(2i+1)2+1(2N+1)0000(2+1)!(2+1)!=1=1=1=12(+1)+1()=(2n-1)()+()+()+() i i N NN N Nx N n n ni N n n i n N xf x c xf x c fx c x f x O x ∆∆∆∆∆∑∑∑∑5（）根据上式，可得差分系数由下面的方程确定：1113322-12-1113(2N-1)013(2N-1)=013(2N-1)N N N c c c ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 求出差分系数n c 之后，则可得到空间导数的2N 阶精度差分近似：()22-12-1122=1={[+]-[-]}+()Nf x N n n n xx n c f x x f x x O x ∂∂∆∆∆∆∑ （6）3）一阶弹性波方程差分格式图1 精度为O(24+tx ∆∆)的交错网格示意图(据Levander 1988)如图1所示，水平速度x V 定义在网格点(i,j)，垂直速度z V 定义在网格点1122(i+,j+)正应力xx σ,zz σ都定义在网格点12(i+,j)，剪切应力xz σ定义在网格点12(i,j+)。
速度分量x V 、z V 均定义在时间层12-k ，12+k ，应力分量xx σ,zz σ,xz σ均定义在时间层k ,k+1。
设+12+12,+12,+12+12,+12,,+12,,,,k k k k ki j i j i j i j i j U V P Q S 分别为质点速度分量x V 、z V 和应力分量xx σ,zz σ,xz σ的离散值，ρ,λ,μ的离散值分别为,,,,,i j i j i j l m ρ，水平与垂直方向网格间距均为Δx 。
