声波方程正演模拟

k

式中

k 1

v(i, j ) 为介质速度的空间离散值

(4-7)

h 是空间离散步长 t 为时间采样率
30

同理可得空间四阶精度的差分格式为：

 u
2

x
2

2



1 x

{ 2

1 12

[u

k i  2, j

u

k i  2, j

]

4 3

[u

k i 1, j

u

k i 1, j

]

5 2

ui , j }
21

相对于上述几种方法，有限差分法是一种更为快速
有效的方法。虽然其精度比不上有限元法，但因其具

有计算速度快，占用内存较小的优点，在地震学界受
到广泛的重视与应用。

22

内容提纲
一、地震勘探基本原理
二、波动方程类型及其 局限性

三、数值算法类型及其局限性
四、声波方程的有限差分法数值模拟

23

四、声波方程的有限差分法数值模拟

t

2

(4-4)

28

同理可得(i,j)点k时刻的二阶空间微商分别为：

 u
2

x
2

2



u i 1, j  2u i , j  u i 1, j
k k k

h
k k

2

(4-5)

 u z
2



u i , j 1  2u i , j  u i , j 1
k

h

2

(4-6)

29

这就实现了用网个点波场值的差商代替了偏微 分方程的微商，将上三个式子代入(4-1)式中得：
胡克定律（这种理想介质称虎克固体），但波在实

际介质中传播时，是有能量损耗的，这就是所谓的
弹性波吸收。波在传播过程中，实际介质的不同部 位之间会出现某种摩擦力，称为内摩擦力或粘滞力。 这种力导致机械能向其他形式能量转换，最终转化 为热能消耗掉。
16

在地震勘探中，地震波传播的实际介质是十分

复杂的。在一定条件下，即震源作用时间短，作用
2

2

h
1 12

2

{
k i  2, j

1 12

[u

k i  2, j

u
4 3

k i  2, j

]

4 3

[u

k i 1, j

u
5 2

k i 1, j

]

5 2

ui , j }

k



v t h
2

{

[u

u

k i  2, j

]

[u

k i 1, j

u

k i 1, j

]

ui, j
能够描述且只能描述纵波的传播规律，包括 直达波、反射波、透射波、折射波等，但不能描 述转换波传播规律。

需要的已知条件包括：
1）震源函数 2）地层速度/密度 3）边界条件

13

2、弹性波方程：
2 2   2u  u 2  u  vp (  )  S (t )  2 2 2  t x z  2 2 2  w  w  w 2   vs (  ) 2 2  t 2 x z  v x  xx  xz    t x z  x z v z  zz    t x z  xx  vx v z   (  2u )  t x z  vx  zz v z    (  2u ) t x z  xz v x v z  u u 14 t x z

方程等）、震源和地下几何界面、物性参数（岩层
密度、速度等）情况下，研究弹性波或声波的传播

规律。
10

内容提纲
一、地震勘探基本原理
二、波动方程类型及其 局限性

三、数值算法类型及其局限性
四、声波方程的有限差分法数值模拟

11

二、波动方程类型及其局限性
1、声波方程： 二阶标量声波方程：
 p
2

t

2

v (

i  2, j
i  2, j 1 i  2, j2

z

i  2, j 1 i  2, j2

i, j  1

i, j  2

x

图4-1 差分网格划分示意图
25

网格间隔长度 时间采样步长
，

h

t

x  ih
z  jh

t  kt

u

k i, j

u i, j

k

表示(i,j)点k时刻的波场值
边界的处理是地震波场正演数值模拟的一个重要课题。

本次作业不涉及边界条件的使用，可通过增大模型
来避免边界反射干扰，有精力的同学可通过查阅资料获

得解决边界问题的方法。

38

本次大作业的具体要求为：
1、应用声波方程作为正演模拟的波动方程；

2、将所提供震源函数离散后绘图；
2

 p
2

x
2

2



 p
2

z

2

)  S(t)

P t v x

一阶压力-速度方程组：

x 1 P  ( ) t  x 1 P  ( ) t  y

  C (

v x



v y y



v z z

)

v y

v z

1 P  ( ) t  z

12

二、波动方程类型及其局限性
33

 (i  i0 ) *  ( j  j 0 )

4.1 稳定性条件
对于特定的偏微分方程只有特定的几种有限差分
格式是无条件或有条件稳定的，（4-7）、（4-8）式即 是已被证明的有条件稳定格式，其稳定性条件分别为：

v max * t / h 

2/2

(4-9)

vmax * t / h  3 / 8
声波方程数值模拟
――地球物理场论

基础Ⅰ期末作业（1） 任课教师： 宋鹏
1

内容提纲
一、地震勘探基本原理
二、波动方程类型及其 局限性

三、数值算法类型及其局限性
四、声波方程的有限差分法数值模拟

2

一、地震勘探基本原理
＊ ▽ ▽ ▽ ▽ ▽

x

t

3

同相轴为双曲线，即反射波的时距曲线为双

曲线，反射波一个同相轴可带来一个地层的信息。
实际地下介质非常复杂，所得到的炮集记录也包含

更多的地下信息。实际的炮集记录见图1-1和1-2。

4

图1-1 陆上某区实际地震记录
5

图1-2 海上某区实际地震记录
6

广义的地震反演即是从地震炮集记录出发，经 过复杂的去噪、速度分析以及偏移成像处理等手段 得到反映地下的地质结构的地震剖面。实际的地震 剖面见图1-3和1-4。
ui , j  (v * t / h) * (ui 1, j  ui 1, j  ui , j 1  ui , j 1 )
2 k k k k k 1

(

2h  4v t
2 2

2

h

2

)ui , j  ui , j  s (k ) *  (i  i0 ) *  ( j  j0 )

k

 s (t ) *  (i  i0 ) *  ( j  j0 )

(4-8)
32

s (k为震源函数， )

一般使用一个理论的雷克型子波代替，即：
s (t ) 


e

( 2f /  ) t
2

2

cos 2ft

t


为时间

f 为中心频率，一般取为20-40HZ
为控制频带宽度的参数，一般取2-5 确定震源位置
k k 1

u t
t  k *t

* t 

1 u
2

2 t

2

t  k *t

* t  o(t )
2 2

(4-3)

27

将上两式相加，略去高阶小量，整理得(i,j)点k时刻 的二阶时间微商为：

 u
2

t

2



u i , j  2u i , j  u i , j
k

k 1

k 1

内容提纲
一、地震勘探基本原理
二、波动方程类型及其 局限性

三、数值算法类型及其局限性
四、声波方程的有限差分法数值模拟

18

三、数值算法类型及其优缺点
地震波波动方程数值模拟方法主要包括克希霍夫积

分法、傅里叶变换法、有限元法和有限差分法等。
克希霍夫积分法引入射线追踪过程，本质上是波动

方程积分解的一个数值计算，在某种程度上相当于绕射
及地下介质的波速满足稳定性条件及频散关系式时，就 可以应用(4-7)式递推求得所给的速度—深度模型内任意 时刻、任何采样点的波场值，并可最终得到一个人工合 成的地震记录。

37

*4.3 边界条件
在地震波场正演模拟中，必须引入人工边界来界定
计算区域。人工边界若不做特殊处理，就会随着波场的

递推计算在边界上产生虚假反射波从而扰乱波场，人工

7

图1-3 陆上某区地震剖面
8

图1-4 海上某区地震剖面
9

地震波场模拟即地震正演，是指已知模型结构， 通过物理或数值计算的方法模拟该地质结构下的地 震波的传播，最终合成地震记录，也可以认为其是 野外数据采集过程的室内再现。物理模拟花费昂贵，

人们一般采用比较经济的数值模拟技术。地震波场
数值模拟是在给定数学模型（如弹性波方程，声波
对于二维速度-深度模型，地下介质中地震波
的传播规律可以近似地用声波方程描述：
 u
2

t

2

v (
2

 u
2

x

2



 u
2

z

)  S (t ) 2

(4-1)

v ( x, z ) 是介质在点（x , z）处的纵波速度，

u
s (t )

为描述速度位或者压力的波场， 为震源函数。

24

空间模型网格化（如图4-1所示）：

(4-10)

34

4.2 频散关系式
同时，在差分计算过程中，如果空间和时间采样间隔不当， 就会导致波形畸变，甚至派生出多个同相轴，这种现象称为频散 现象。 偏微分方程本身没有频散，网格频散是由于差分方程近似替

代微分方程引起的。当波场按照波动方程所表示的微分方程传播
时，波场的传播速度就是波动方程中的速度，但当波场按照波动 方程离散化后的差分方程传播时，波场的传播速度就不再是波动 方程中的速度了，而是与波的频率和波数有关的函数，具有不同 频率和波数的波有不同的传播速度，因而在传播过程中会出现频

能够描述纵、横波的传播规律，包括直达波、反
射波、透射波、折射波以及转换波等。 需要的已知条件包括： 1）震源函数 2）地层速度或根据方程的类型需要提 供的地层的其它弹性参数 3）边界条件

15

3、粘声波/弹性波方程 前面讨论的是理想弹性介质，波在其中传播时，

没有能量的损耗，介质中应力和应变关系严格遵循

小，但缺点是计算速度较慢，对模型的适用性差，尤
其是不适应于速度横向变化剧烈的模型。

20

波动方程有限元法的做法是：将变分法用于单元分 析，得到单元矩阵，然后将单元矩阵总体求和得到总体 矩阵，最后求解总体矩阵得到波动方程的数值解；其主 要优点是理论上可适宜于任意地质体形态的模型，保证 复杂地层形态模拟的逼真性，达到很高的计算精度，但 有限元法的主要问题是占用内存和运算量均较大，不适 用于大规模模拟，因此该方法在地震波勘探中尚未得到 广泛地应用。

k

 u z
2



1 z

{ 2

1 12

[u

k i , j 2

u

k i, j2

]

4 3

[u

k i , j 1

u

k i , j 1

]

5 2

ui , j }

k

31

同理可得时间二阶、空间四阶精度的声波方程差分格式为：

u

k 1 i, j

 2u

k i, j

u
2

Hale Waihona Puke Baiduk 1 i, j
2



v t
26

时间二阶、空间二阶差分格式推导如下： 将
k 1

u 在(i,j)点k时刻用Taylor展式展开：
k

k 1 i, j

ui, j  ui, j 

u t
t  k *t

* t 

1 u
2

2 t 2

t  k *t

* t  o( t )
2 2

(4-2)
ui, j  ui, j 
i  2, j2 i  2, j 1 i  1, j 1 i  1, j2 i  1, j 1 i  2, j2 i  2, j 1

i, j  2

i, j  1

i  2, j

i  1, j
i  1, j 1
i  1, j2

i, j

i  1, j
i  1, j 1 i  1, j2
叠加。该方法计算速度较快，但由于射线追踪中存在着 诸如焦散、多重路径等问题，故其一般只能适合于较简 单的模型，难以模拟复杂地层的波场信息。

19

傅里叶变换法是利用空间的全部信息对波场函数
进行三角函数插值，能更加精确地模拟地震波的传播

规律，同时，利用快速傅里叶变换(FFT)进行计算，还
可以提高运算效率，其主要优点是精度高，占用内存
力微小，地球介质可以看作完全弹性模型，但随着

地震勘探技术的发展，勘探精度要求提高，面临复
杂地质目标时，要求地震勘探采用更加符合实际的 介质模型进行研究。粘弹性介质模型更符合实际。 但是到目前为止，在地震资料反演处理中应用 最多的还是声波方程，弹性波以及粘弹性波方程的 应用还只是停留在模拟层次上。
17

散，发生畸变，且随走时的增加而增加。
35

Dablain给出了一个能有效减少频散的经验公式为：

h  v min /( Gf N )
f N 为Nyquist频率，一般取为 主频的两倍

(4-9)

G为每个波长所占的网格点数，时间、空间为两阶差
分的情况G取8，而时间、空间为四阶差分的情况G

取4。

36

当所给震源函数、空间网格间隔、时间采样间隔以