正演模拟作业

2h 2 4v 2 t 2 k k 1 ( ) u u i, j i , j s ( k ) * (i i0 ) * ( j j0 ) 2 h
震源位置控制函数
如何得到更高阶精度的有限差分格式？
2 2 2u u u 2 v ( 2 ) S(t) 2 2 t x z
������+1 ������+1 需将������������,������ 、������������,������ 展开到四阶精度才可
时间二阶、空间四阶精度的差分格式为：
2 2 v t 1 k 4 k 5 k k 1 k k 1 k k ui , j 2ui , j ui , j { [ ui 2, j ui 2, j ] [ui 1, j ui 1, j ] ui , j } 2 h 12 3 2 v 2 t 2 1 k 4 k 5 k k k { [ u i 2, j u i 2, j ] [u i 1, j u i 1, j ] u i , j 2 12 3 2 h
u t
1 2u t k *t * t 2 t 2
2 2 * t o ( t ) t k *t
将上两式相加，略去高阶小量，整理得(i, j )点k时刻的二阶时间微商为:
u 2 t
2 1 k k 1 u ik, 2 u u j i, j i, j
根据单程波的近似程度，吸收边界条件方程具有不同的阶数 一阶吸收边界条件方程： 二阶吸收边界条件方程： 三阶吸收边界条件方程： 四阶吸收边界条件方程：
( v )u 0 t x
2 2 1 2 2 ( 2 v v )u 0 tx 2 z 2 t
3 3 1 3 3 3 2 3 ( 3 v 2 v v )u 0 t t x 4 xz 2 4 tz 2
h vmin /(Gf N )
G取2或4(四阶精度差分格式G取2)，vmin为最小速度， f N 为奈奎斯特频率，一般取主频的两倍。
时间采样步长
t
太小意味着极大的计算和内存消耗，太大 意味着精度太低，更为关键的是，太大会 产生误差累计，出现不稳定现象。 稳定性条件为：
1 vmax t 2 2 x z 1

uik, j m uik, j m z
2
)
v t c0 (
2 2
u
x
k i, j 2

u
z
k i, j 2
) s (k )* (i i0 )*( j j0 )
震源位置控制函数
二阶、四阶精度差分格式模拟结果
s ( k )为 震 源 函 数 ，
一般将一个理论子波离散化：
i 2, j 1 i 2, j2
i, j 1 i, j 2
（3）连续函数离散化 网格间隔长度 时间采样步长
，
h
x ih z jh
t
t kt
u
k i, j
表示(i,j)点k时刻的波场值
（4）微商化为差商
2 2 2u u u 2 v ( 2 ) S(t) 2 2 t x z
t 2
2 2 2u u u 2 v ( 2 ) S(t) 2 2 t x z
同理可得(i, j )点k时刻的二阶空间微商分别为:
u 2 x
2
u ik1, j 2u ik, j u ik1, j h 2
u 2 t
2 1 k k 1 u ik, 2 u u j i, j i, j
将u 、u 分别在(i, j )点k时刻应用Taylor展式展开：
1 k u ik, u j i, j
k 1 i, j
k 1 i, j
u t
1 2u t k *t * t 2 t 2
2 2 * t o ( t ) t k *t
二阶 精度
1 k u ik, u j i, j
s (t ) * (i i0 ) * ( j j0 )
N
进一步可得时间二阶、空间2m阶精度的差分格式为：
1 k k 1 2 2 uik, 2 u u v t cm ( j i, j i, j m 1
cm为差分系数
2
uik m, j uik m , j x
4 4 4u 4u 4u 3 u 4 u 4v 3 6v 2 2 4v v 4 0 4 3 t t x t x t x x
完全匹配层方法： ◆在波场外部构建一系列的吸收层（匹配层）； ◆吸收层内设置衰减因子，地震波在吸收层内传播按指 数规律衰减； ◆中心波场与边界方程协调计算。
◆通过用网格点上函数的差商代替微商，将含连续变量
的波动方程求解问题化成差分方程的求解； ◆求解差分方程实现地震波正演模拟。 ◆同样要考虑初、边值条件（一般都设为零）。
声波方程有限差分数值模拟： （1）二阶标量声波方程：
2 2 2u u u 2 v ( 2 ) S(t) 2 2 t x z
地震记录（seismic record）: 某一时刻所有网格点的波场值
模拟实验（1）：各种波的认识
直达波、透射波、反射波
多次波
折射波
绕射波
模拟实验（2）： 各种模型的地震波 传播规律
水平层状介质模拟：
断层模型模拟：
盆地模型模拟：
本次大作业的总体要求为： （1）C程序或Fortran程序实现波动方程的有限差分数值模拟， 图件应用Matlab绘制。 （2）以科技论文的形式撰写，包括题目 、摘要、关键词、正 文等。
程序实现流程
（1）开辟子波数组（一维）、速度模型数组（二维）， 波场数组（三维） （2）模型和子波离散化（选定合适的空间、时间间隔） （3）子波赋值，速度赋值，波场初始化（包括边界）； （4）三重循环，计算波场。
（5）输出成果（波前快照和地震记录）
波前面（wavefront）: 某一时刻所有网格点的波场值
雷克型子波波形示例
雷克子波（零相位） 1
0.5
0
-0.5
0
10
20
30
40
50
60
70
80
90
100
网格间隔长度 时间采样步长
，
h
t
空间网格和时间间隔能否 无限大或无限小？
网格间隔长度
h
太小意味着极大的计算和内存消耗，太大 意味着精度太低，精度通常通过频散关系 式进行评估。 频散关系式为：
i 1, j2 i 1, j 1
i 2, j2 i 2, j 1
i, j 1
i 2, j
h
i 1, j
i 1, j 1
i 1, j2
h
i, j
i 1, j
i 1, j 1 i 1, j2
i 2, j
i 2, j 1 i 2, j2
◆适应于正常矩形网格模拟； ◆能够描述且只能描述纵波的传播规律，包括直 达波、反射波、透射波、折射波等，但不能描 述横波传播规律。
◆需要的已知条件包括： 1）震源函数 2）地层速度 3）边界条件
（2）空间模型网格化
i 2, j2 i 2, j 1 i 1, j 1
i, j 2
(2 f / ) t cos 2 ft s (t ) e
2 2
雷克型 子波
s (k ) e (2 f / ) ( k t *k t ) cos 2 f (k t )
2 2
f 为中心频率，一般取为15 25HZ
时间采 样间隔
为控制频带宽度的参数，一般取25
地震正演模拟需给出波场边值条件，一般分两种情况： ①边界上的位移为零（狄利克雷边界）：
u
k i, j
0
②吸收边界： 可有效消除边界反射
Hale Waihona Puke 利克雷边界会产生强烈的边界反射：吸收边界方法（主要有Higdon边界和PML边界等）， 可有效吸收边界反射：
当前主要的边界处理方法：
◆吸收边界条件方法（ABC）； ◆完全匹配层方法（PML）。 吸收边界条件方法： ◆从双程波动方程中分离单程波方程作为边界条 件方程； ◆基于边界条件方程构建边界条件差分方程； ◆中心波场与边界方程协调计算。
vx
k 1 2 i, j 1 2
1 k 1 1 t N k k 2 1 0 . 5 t v a P P 1 x m i , j m i , j m1 k i, j 1 0.5 t 1 x m1 2 i, j 2
k k k 2 u u i , j 1 2u i , j u i , j 1 2 z h 2
t 2
（5）形成差分方程，并可递推求解
将差商代入声波方程可得：
1 2 k k k k uik, ( v * t / h ) * ( u u u u j i 1, j i 1, j i , j 1 i , j 1 )
本次大作业的具体要求为： 1.绘制出雷克型子波
雷克子波（零相位） 1
0.5
0
-0.5
0
10
20
30
40
50
60
70
80
90
100
2、设置一个均匀介质模型，整个区域的速度值可设为常数，即 只有一种介质，将震源点放在模型中间，分别记录两个时刻的 波前快照（即该时刻区域内所有网格点的波场值）。第一时刻 为地震波还未传播到边界上的某时刻，第二时刻为地震波已经 传播到边界上的某时刻，体会其人工边界反射；



v
1 2 z 1 i , j 2 k