弹性波理论

地震波交错网格高阶差分数值模拟研究
摘要: 地震波数值模拟技术是勘探地球物理学中的重要组成部分，研究通过弹性波一阶速度——应力方程，采用交错网格高阶有限差分法实现了地震波在各向同性介质中的高精度的数值模拟，并采用完全匹配层( 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 xz
xz z zz xx x z x zz z xz x z v t x z v t x z
v v t x z v v t z x v v t z x σσσσσσσρ
ρ
λμλ
λμλ
μ
∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂⎧=+⎪=+⎪⎪=++⎨⎪=++⎪⎪=+⎩
分用有限差分来代替。

设函数f(t)单值连续且存在任意阶导数，则其Taylor 展开式为：
232
3
231112
2
2!2
3!2
!2
()()()()...+
()+()m m
f f f f m
m t t t t t t m t t t f t f t O t ∂∂∂∂∆∆∆∆∆∂∂∂∂+=+
+
+
+∆23
2
3
231112
2
2!2
3!2
!2
-(-)()-()
()...+
(-)+()m m
f f f
f m
m t
t t t t t m t t t f t f t O t ∂∂∂∂∆∆∆∆∆∂∂∂∂=+
+∆将以上两式相减，得到2M 阶精度的时间差分近似式为：
2-12-12-1
1(2m-1)!
=1
22
2
2
*(
)
*+()=(-)+2*()m m M
f
m t m M t
t t f t f t O t ∂∂∆∆∆+
∆∑
（2）
Δt 为时间步长，在采用交错网格技术解声波方程时，速度场在2
+t
t ∆ 时刻计算，而应力场在t +Δt 时刻计算。

当M ＝1时，对方程(1)中的时间微分进行数值离散，可得2阶时间差分近似为：
22(+)=(-)+[]xx xz t t t x x x z v t v t σσρ∂∂∆∆∆∂∂+
2
2
(+)=(-)+[
]xz zz t
t t
z z x
z
v t v t σσρ∂∂∆∆∆∂∂+
(+)=()+[(2)
]
x z
v v xx xx x
z t t t t σσλμλ∂∂∂∂∆∆++ (+)=()+[(2)]
x
z
v v zz zz x z t t t t σσλλμ∂∂∂∂∆∆++
(+)=()+[]x z v v xz xz x
z
t t t t σσμ
μ
∂∂∂∂∆∆+
对于时间高阶差分，计算
2-12-1
m m f t ∂∂时将涉及较多的时间层，需要大量的存储空
间因此利用一阶弹性波方程中的速度－应力关系式，将速度对时间的任意奇数阶高阶导数用应力对空间的导数代替，而应力对时间的任意奇数阶高阶导数用速度对空间的导数代替，这样，在计算一个时间层上的速度或应力场时，只需要前一个时刻的速度或应力场以及两时间层之间的应力或速度场，不需要过多的时间层，从而节省了内存(董良国，2000)。

2）一阶弹性波方程空间导数的2N 阶差分精度算法
为了提高计算精度，空间导数也要采用高阶差分近似。

在交错网格技术中，对空间变量的导数是在相应的空间变量网格点之间的半程上计算，对于空间函数f (x)，假设其存在2N ＋1阶导数，则f (x)在2-1
02
=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 ()+()N
f x N N n n n N x x x
n x
c f x x f x x f x x O x ∂∂∆∆∆∆∆∑ （4）
其中，n c 为差分系数，N e 为误差系数； 将(3)式代入(4)式中，化简得到：
2+12+1
2+1
-1
(2n-1)(2n-1)(1)
(1)
(2i+1)
2+1(2N+1)0000(2+1)!
(2+1)!=1
=1=1
=1
2(+1)+1()=(2n-1)()+()+()
+() i i N N
N N N
x N n n n
i N n n i n N xf x c xf x c f
x 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-11
22
=1={[+]-[-]}+()N
f x N n n n x
x n c f x x f x x O x ∂∂∆∆∆∆∑ （6）
3）一阶弹性波方程差分格式
图1 精度为O(2
4+t
x ∆∆)的交错网格示意图(据Levander 1988)
如图1所示，水平速度x V 定义在网格点(i,j)，垂直速度z V 定义在网格点
1122(i+,j+)正应力xx σ,zz σ都定义在网格点1
2(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 k
i 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 。

则方程(1)的精度为O(2
2+N
t x
∆∆)差分格式为：
2-12-12-12-1,2
2
2
2
+12-12,,+,-,,+,-=1
=1
=+[c (-)+c (S -)]n n n n i j N
N
k k k
k k k t i j
i j
n
n x i j i j i j i j n n U
U
P P
S ρ∆∆∑∑
1
1+12,+122
2
+12-12+12,+12
+12,+12
+12,++12,+1-+,+,-=1=1=+[c (S
-)+c (-)]
i j N
N
k k k
k k
k
t
i j i j n n i j n
i j n
x i n j i j n n V
V
S Q
Q
ρ∆∆∑∑2-12-12
2
+1+12+12+12+12+12,+12,+12,+,+1-,+12,+12,++12,-=1=1
=+[(+2)c (-)+c (-)]
n n N
N
k k k k k k t i j i j i j n i n j i n j i j n x i j i j n n P P l m U U l V V ∆∆∑∑2-12-12
2
+1+12+12+12+12+12,+12,+12,+12,+,+1-,+12,++12,-=1
=1
=+[(+2)c (-)+c (-)]
n n N
N
k k k k k k t i j i j i j n i j n i n j i n j x i j i j n n Q Q l m V V l U U ∆∆∑∑2-12-12
2
+1+12+12+12+12,+12,+12,+12,+,+1-+,+12-,+12=1
=1
=+[c (-)+c (-)]
n n N N
k k k k k k t i j i j i j n i j n i j n n x i j i j n n S S m U U V V ∆∆∑∑同理，可以根据实际需要选择不同的M,N 来达到预定的精度，精度越高，计算量相应的越大。

完全匹配层吸收边界条件
采用有限差分方法进行数值模拟时，边界条件的处理是非常关键的，选择合适的边界处理条件，不仅可以得到效果较好的数值模拟结果，而且也能提高地震照明结果的准确性。

因此，对弹性波方程的数值模拟采用Collino(2001)提出的弹性波方程中的完全匹配层吸收边界条件进行边界处理。

1）完全匹配层原理
基于波场分裂的完全匹配层原理是在PML 层中，将每一个波场分量可分解为两部分(二维情况)或三部分(三维情况)，每一部分代表一个空间变量的偏导数对场分量的贡献，且每一部分只在其坐标轴方向衰减，然后将同一个场分量的各个分裂部分相加得到该场分量，并用于下一个循环迭代，求解其他场。

由于只是在PML 层对波场分量进行分解，而PML 层相对于整个模型计算空间还是很小，因此增加的内存量并不很大。

考虑一阶双曲型偏微分方程定解问题：
--=0
(=0)=U U U
t x z A B U t U ∂∂∂∂∂∂ （7）
现在寻找一组新的方程，使之在PML 模型区域(S1,S2,S3)(如图2)与内部波场
的方程互相耦合，在两个区域之间的边界x ＝0上不会产生反射，并且在PML 区域内，波是呈指数衰减传播的。

图2完全匹配层示意图
首先，将U 分解为垂直于边界x ＝0和平行于边界x ＝0的两个分量U ⊥和U ，
U ⊥表示只包含对x 的偏导数，U ⊥表示只包含对z 的偏导数。

-=0-=0U U t z U U
t
x B A ⊥∂∂∂∂∂∂∂∂⎧⎪⎨⎪⎩ （8） 其中：=+U U U
⊥
其次，构造一个新的方程组，使之在PML 区域内指数衰减。

在U ⊥分量上引入一个衰减因子d (x)，方程的解为一个新的波场u ：
0+d()-=0-=0(t=0)=U u U t x u u
t z x u A B u ⊥
⊥∂∂∂∂∂∂∂∂⎧⎪⎨⎪⎩
（9）
其中：=+U U U
⊥
对比(8)和(9)可知，在内部波场(x<0)，u ≡U ，将(8)和(9)经Fourier 变换到频率域，根据波函数在Fourier 变换中
t
∂
∂与i ω的互换性，得：
-=0-=0U x U z i U A i U B ωω⊥
∂∂∂∂⎧⎪⎨⎪⎩
（10）
(+())-=0-=0u x u
z i d x u A i u B ωω⊥∂∂∂∂⎧⎪⎨⎪⎩
（11）
,,u ,u ,,u U U U ⊥⊥分别为,,u ,u ,,u U U U ⊥⊥的Fourier 变换。

(11)也可以写为：
+()-=0-=0i u
i d x x u
z i u A i u B ωωωω⊥∂∂∂∂⎧⎪⎨⎪⎩
（12）
对比(10)和(12)可知：PML 模型是由内部波场经过简单代换得到的：
+()=,i x
x i d x x ωω∂∂∂
∂∂∂→
即相当于对变量x 进行复变换：
()=-(s)x
i
x x x d ds ω⎰ （13）
设方程(7)存在如下形式的特解：
0=exp[-i(k x+k z-t)]=(k ,k )
x z x z U U k ω （14）
上式表示沿着k 方向以相速度k
ω传播的平面波，U 满足以下关系式 000++=0x z i U iAU k iBU k ω（15） 同样，(9)式的特解为如下形式：
0=exp[-i(k x+k z-t)]x z u U ω （16） 将(13)式代入(16)式中，得到：
00
=exp[-i(k x+k z-t)]exp[-(s)]x
x
k x z u U d ds ω
ω⎰
（17）
综上所述，可以得出如下结论：
1.在x ≤0区域，u ≡U ，表明PML 模型与内部波场是完全匹配的。

2.在PML 区域(x>0)，u 按指数衰减，衰减系数为：
(x)(x)
=
=exp[-(s)]x
x
u k d U d ds ω
α⎰
衰减与波的传播方向有关，体现在x k 上，对正入射波衰减最快，但对接近于平行入射到该界面的波，衰减速度很慢。

2）完全匹配层差分格式
将PML 原理应用到一阶弹性波方程，进行波场分裂，得到一组新的方程式：
=+,=+v v v σσσ⊥⊥
(+d(x))=;=(+d(x))=
;=xx
x
xz
xz z zz v x t x t z
v z t x
t z
v v σσ
σσρρρρ
∂∂∂⊥∂∂∂∂∂∂∂∂⊥
∂
∂∂∂∂
(+d(x))=(+2)
;
=x xx
z
xx v v t
x
t
z σ
σλμλ∂∂∂⊥
∂∂∂∂∂ （18）
(+d(x))=;
=(+2)x zz
z v v zz t x
t
z σ
σλ
λμ∂∂∂⊥∂∂∂∂∂
(+d(x))=;
=xz
x z v v xz t x
t
z σ
σμ
μ
∂∂∂⊥∂∂∂∂∂
其中，d (x)，d (z)为衰减因子，在内部波场为0，s1区域d (x)>0，s2区域d (z)>0，在s3区域，d (x)>0，且d (z)>0
完全匹配层的离散方法仍然采用交错网格差分格式，内部波场和PML 的交界面在离散过程中作为内部弹性界面处理。

则得到PML 差分格式：
,,,()=()+()n x n z n
x i j x i j x i j v v v
+12+12
+1+1+12,-12,,,,,,()-()(v )-(v )(v )+(v )2
+=n n x n x n x n x n
xx xx i j i j
x i j x i j
x i j x i j i j x i t
h
d σσρ∆
+12+12+1+1,+12,-12
,,,,,()-()(v )-(v )(v )+(v )2
+=
n n z n z n z n z n
xz xz i j i j x i j x i j
x i j x i j i j z j t
h
d
σσρ∆
+12,+12+12,+12+12,+12()=()+()n x n z n z i j z i j z i j v v v
+12+12
+1+1+12,+12+12,+12
+12,+12+12,+12
+1,+12,+12
+12,+12()-()(v )-(v )(v )+(v )+122
+=n n x n x n
x n x n
xz xz z i j z i j z i j z i j i j i j i j x i t
h
d σσρ∆
+12+12+1+1+12,+12+12,+12
+12,+12+12,+12
+12,+1+12,+12,+12()-()(v )-(v )(v )+(v )j+122
+=
n n z n z n z n z n
zz zz z i j z i j z i j z i j i j i j
i j z t
h
d
σσρ∆
+12+12+12
+12,+12,+12,()=()+()n x n z n xx i j xx i j xx i j σσσ
+12-12+12-12
+12,+12,+12,+12,+1,,i+12,()-()()+()()-()+12i+12,2
+=(+2)
n n n n x x x x n n xx xx xx xx i j i j
i j i j x i j x i j
j v v x i j t
h
d σσσσλμ∆
+12-12+12-12+12,+12,+12,+12,+12,+12+12,-12
()-()()+()()-()i+12,2
+=n n n n z z z z n n xx xx xx xx i j i j
i j i j z i j z i j v v z j j
t
h
d
σσσσλ∆
+12+12+12
+12,+12,+12,()=()+()n x n z n zz i j zz i j zz i j σσσ
+12-12+12-12
+12,+12,+12,+12,+1,,()-()()+()()-()+12i+12,2
+=n n n n x x x x n n zz zz zz zz i j i j
i j i j x i j x i j
v v x i j
t
h
d
σσσσλ∆
+12-12
+12-12
+12,+12,+12,+12,+12,+12+12,-12
i+12,()-()()+()()-()i+12,2
+=(+2)
n n n n z z z z n n zz zz zz zz i j i j
i j i j z i j z i j j v v z j j t
h
d
σσσσλμ∆
+12+12+12
,+12,+12,+12()=()+()n x n z n xz i j xz i j xz i j σσσ
+12-12+12-12
+12,+12,,+12,+12+12,+12-12,+12
()-()()+()()-()i,+12
2
+=n n n n x x z z n n xz xz xz xz i j i j
i j i j z i j z i j v v x i j t
h
d
σσσσμ∆
+12-12
+12-12
,+12,+12
,+12,+12,+1,()-()()+()()-()+12i,+12
2
+=n n n n z z z z n n xz xz xz xz i j i j i j i j x i j x i j
v v z j j t
h
d
σσσσμ∆
其中：x x v 表示x v ⊥，z
x v 表示x v
+1,,(v )+(v ),2
==,(x)(v )=x n x n x i j x i j x n
x x i j
i h x z d d
∆∆
,(v )=(i x,j z,n t)n x i j x v ∆∆∆，其它与此雷同。

数值计算结果
下面采用精度为O (24+t x ∆∆)的交错网格有限差分方法，对各向同性弹性介质中的地震波传播进行数值模拟。

完全匹配层衰减因子采用Collino(2001)给出的
公式：231002(x)=(),=log()x v
R d d d δδ，v 是波传播中最快的速度，x 为计算点与PML
内层边界的距离，δ为PML 层的厚度，0d 为最大衰减系数，它与理论反射系数有关。

1）不同厚度的PML 吸收边界吸收效果比较
均匀介质模型大小为2000m*2000m ，空间网格步长为5m ，震源位(1000m,1000m)处，子波主频为30Hz.图3a,b,c,d,e,f 中PML 网格厚度依次为0、5、
10、15、20、30,均为600ms时刻的波场快照。

由a图知无吸收边界时，当波传播到模型边界时，产生了强烈的反射，随着PML吸收边界的厚度增加，吸收效果也逐渐变好，由d、e、f图知，当PML边界达到一定厚度之后，来自边界的人工反射波已经被很好地吸收了，图e 和f几乎无反射存在。

因此，在实际计算时，可以适当地减少PML边界的厚度来减少计算量。

图3.a PML厚度为0个网格点的波场快照图3.b PML厚度为5个网格点的波场快照
图3.c PML厚度为10个网格点的波场快照图3.d PML厚度为15个网格点的波场快照
图3.e PML厚度为20个网格点的波场快照图3.f PML厚度为30个网格点的波场快照
2）PML 边界采用不同吸收函数吸收效果之比较
完全匹配层（PML ）方法提供了一种依据指数率迅速衰减，适用于任意频率、任意入射角度，并将人工边界的虚假反射迅速衰减为零的外行波模拟技术。

但是，目前关于PML 人工边界条件中衰减函数的使用尚无详细论述。

基于弹性波方程给出了完全匹配层算法的原理，特定选了一下几个函数作比较,具体见下表 图号 4.a 4.b 4.c 4.d
吸收因子α 0(1-)ih d δ 20(1-)ih d δ 02[1-sin()]ML i P d π 02cos()ML i P d π 类型 一次型 二次型 正弦型 余弦型
其中，3102=log()v R d δ v 是传播中最快的速度，δ为PML 层的厚度（=ML P h δ），0
d 为最大衰减系数，它与理论反射系数有关，i 为计算点距离PML 内边界的网格点数（=0,1,2,ML i P ）。

图4.a 衰减函数为一次型的波场快照 图4.b 衰减函数为二次型的波场快照
图4.c 衰减函数为正弦型的波场快照 图4.d 衰减函数为余弦型的波场快照
通过对比、分析多个典型的衰减函数可知，正余弦型吸收衰减因子的效果优于指数型吸收衰减因子的结论。

由上面的两组分析可知，PML厚度对吸收效果的影响比吸收函数对吸收效果的影响更明显，如果选择两者的最佳匹配组合，可以达到更好的吸收效果。

3）PML 边界与衰减边界之比较
为了验证最佳匹配层（PML）的优越性，特地将PML边界吸收与指数衰减吸收进行比较，人工边界的宽度均为10个网格，由图5.a和5.b可知，对于相同的模型，PML的吸收效果明显优于指数衰减的吸收效果，因此，选用PML吸收边界既可以节约计算成本，又可以达到更理想的效果。

图5.a指数衰减边界吸收效果图图5.b PML边界吸收效果图
小结
地震波数值模拟技术在地震学中的地位相当重要，这次学习主要研究了采用交错网格有限差分法实现了地震波在各向同性介质中的高精度的数值模拟，并采用完全匹配层( PML) 吸收边界来消除边界反射，在编程过程中遇到不少问题，幸运的是有师兄和同学们帮忙，问题得到了及时的解决，在此衷心地感谢他们！
参考文献
[1] 皮红梅. 双程波动方程数值模拟和照明分析方法研究[D].吉林：吉林大学，2009.
[2] 李万万. 基于波动方程模拟的地震采集参数论证方法及应用研究[D].北京：中国地质大学，2008.
[3] 牟永光, 裴正林著.三维复杂介质地震数值模拟[M].北京：石油工业出版社，2004.
[4] 周学明, 李庆春,等.地震波交错网格高阶差分数值模拟研究[J].铁道工程学报,2011,8:1~6
[5] 董良国.弹性波数值模拟中的吸收边界条件[J].石油地球物理勘探，1999.34(1)：45-56.
[6] 董良国，马在田，曹景忠等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报，2000，43(3)：
411-419.
[7] 董良国，马在田，曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究[J].地球物理学报，2000，
43(6)：856-864.
[8] 常君勇.交错网格有限差分法地震数值模拟[D].北京：中国地质大学(北京)，2007.
[9] Alford R M,Kelly K R,Boore D M.Accuracy of finite-difference modeling of the acoustic wave
equation[J].Geophysics,1974,39(6):834-842.
[10] Dablain M A.The application of high-order differencing to the scalar wave equation[J].
Geophysics,1986,51(1):54-66.
[11]李卫志, 李正文,等复杂模型高精度有限差分正演模拟[J].内蒙古石油化工,2006,10:132-135
[12] 裴正林, 牟永光.非均匀介质地震波传播交错网格高阶有限差分法模拟[J].石油大学学报。

2003,27(6):17-21
[13] 熊章强,张大洲,等瑞雷波数值模拟中的边界条件及模拟实例分析[J].中南大学学报，
2008,39(4):824-829.。