有限差分法地震波传播数值模拟

)
⎪⎫ ⎬
2 ⎪⎭
v v τ β z
=
0.5⎪⎨⎧ ⎪⎩
(old ) +
z
ρ
(old
zz
)
1
⎪⎫ ⎬ ⎪⎭
τ CC τ Cβ v xx
=
0.5⎪⎨⎧ ⎪⎩
+ 13 (old ) zz 33
13 1
(old )⎪⎫
Z
⎬ ⎪⎭
τ {τ β v } = 0.5 (old ) + ρ
(old )
zz
zz
1z
∂U ∂x
+
B1
∂U ∂z
∂U ∂t
=
A2
∂U ∂x
+
B2
∂U ∂z
其中：
A1 =
1 2
⎜⎛ c ⎜c ⎜⎝ 0
c c 0
0 ⎟⎞ 0⎟ 0 ⎟⎠
B1
=
1 2
⎜⎛ c ⎜0 ⎜⎝ c
0 0 0
c ⎟⎞ 0⎟ c ⎟⎠
A2
=
1 2
⎜⎛ ⎜ ⎜⎝
−c c 0
c −c 0
0⎟⎞ 0⎟ 0⎟⎠
B2
=
1 2
⎜⎛ ⎜ ⎜⎝
•Direct methods (Grid methods or Full-wave equation methods) ------没有限制，网格小时精度很高，效率是问题
声波方程：
∇⋅(1
ρ
ΔP)
=
1 K
∂2P ∂t 2
+
f
(r,t)δ (rr
−
rr0 )
弹性波方程：
形变位移方 程： 本构方 程： 应力运动方程：
离散 z 边界条件的实施是稳定的 问题 z 吸收效率要高
z 计算效率要高，不能因边界条件而引入过量的计算，失 去采用边界条件的意义
吸收边界条件类型：
z ABCs Based on Paraxial Approximation z ABCs Based on Damping Method z PML（Perfect Matched Layer） absorbing boundary conditions
M
C m (M ) 2M +2 m m=1
= 1 + const ⋅φ 2M
其
φ = kΔx = 2π Δx
中，
λ
z 当2M=2时， z 当2M=4时， z 当2M=6时， z 当2M=8时， z …...
c2 c02
≈ 1− φ2 12
c2 c02
≈ 1− φ4 90
c2 c02
≈ 1− φ6 560
上左前角A1B1C1 上左后角A1B2C1 上右前角A2B1C1 上右后角A2B2C1 下左前角 A1B1C2 下左后角A1B2C2 下右前角A2B1C2 下右后角A2B2C2
前人所用的一种弹性波吸收边界条件
τ β v v 底边界吸收边界条件：
x
=
0.5⎪⎨⎧ ⎪⎩
(old ) +
x
ρ
(old
xz
2
MM 4
⎤ ⎥ ⎥
⎢⎢⎡CC12((MM
) )
⎤ ⎥ ⎥
M = M
6
⎥ ⎥ ⎥
⎢⎢C3(M
)
⎥ ⎥
⎢M⎥
2M ⎥
M⎦
⎢⎣CM(M ) ⎥⎦
⎡1⎤ ⎢⎢0⎥⎥ ⎢0⎥ ⎢⎢ M ⎥⎥ ⎢⎣0⎥⎦
⎡1
⎢ ⎢
13
⎢ 15
⎢ ⎢
M
⎢⎣12N −1
3 33 35
M 32 N −1
5 53 55
M 52 N −1
一般吸收边界条件存在的问题：
z 吸收程度不高 z 边角（3D）及四个角（2D）处的人为边界反射难以处理 z 存在计算稳定性问题
ABCs Based on Paraxial Approximation
Acoustic Wave: (Bottom Boundary)
Pxx
+ Pzz
−
1 v2
Ptt
ቤተ መጻሕፍቲ ባይዱ=0
kz = ±
ω2
v2
−
k
2 x
kz
=
ω
v
1
−
v
2
k
2 x
ω2
kz
=
ω
v
⎪⎨⎧1 − ⎪⎩
1 2
v
2
k
2 x
ω2
−
1 8
⎜⎜⎝⎛
v
2
k
2 x
ω2
⎟⎟⎠⎞2
+ LL⎪⎬⎫ ⎪⎭
First-order Approximation:
kz
≈
ω
v
∂P + 1 ∂P = 0 ∂z v ∂t
Second-order Approximation:
∑C ∂f = 1 N
∂x Δx n=1
( N )⎧
n
⎨ ⎩
f
⎢⎣⎡ x
+
Δx 2
(2n
−1)⎥⎦⎤
−
f
⎡ ⎢⎣
x
−
Δx 2
(2n
−
1)⎥⎦⎤
⎫ ⎬ ⎭
+
O(Δx
2
N
)
高阶差分系数的确定：
⎡ ⎢ ⎢
12 14
⎢ ⎢
16
⎢M
1⎢ 2M
⎣
22 24 26
M
22M
32 34 36
M
32M
L L L O L
-----------J.M. Carcione
地震波传播数值模拟应用领域
地震波传播理论
数据采集
理论指导 物性参数
研究传播规律
正演模拟
指导设计 观测系统
验证
地震解释
提供理论数据 试验处理流程
数据处理
提供正演方法
岩石物理
参数反演
断层下覆界面反射能量强
炮点
T=2000ms
炮点 T=2300ms
炮点位于11km处的单炮记录
有限差分法 地震波传播数值模拟
主要内容
一、有限差分法地震波传播数值模拟方法 二、网格频散问题 三、算法稳定性问题 四、边界问题（自由边界、吸收边界） 五、震源模拟问题
一、地震波传播数值模拟技术
地震波模拟就是模拟地震波在地下的传播，模 拟的目的就是在给定地质模型情况下，预测可能 接收的地震记录。
“Seismic modeling is one of the cornerstones of geophysical data processing.”
1992，1994年，Tessmer et al.在模拟二维以及三维不规则 地表面波时，同样也是使用的如上吸收边界。
AA
+
A
M −1 (+)
AA
) ∂U + ( A ∂x
M −1 (−)
BB
+
B
M −1 (+)
BB
) ∂U B ∂z
L L M −1 (+) AA
∂U A ∂x
L L M −1 (+) BB
∂U B ∂z
L L M −1 (+) AA
∂U A ∂x
我们采用的吸收边界条件如下（2D）：
左边界
∂U ∂t
=
kΔx
=
2π λ
Δx ≤ 1
，即只要一个波长包含几个空间步
长，随着差分精度2M的提高，上述高阶差分解法产生
的数值频散会逐渐减小。
不同差分精度空间频散曲线
不同差分精度时间频散曲线
五、边界问题
自由边界条件
内部边界条件 吸收边界条件
设计吸收边界条件的目标：
z 方程+边界条件数学上是非病态的
连续
问题 z 方程+边界条件可以近似描述无限介质中的物理过程 z 边界条件和内部点的计算方式是相容、不冲突的
L L L O L
2N −1
C (2N −1)3 C (2N −1)5
C M C (2N −1)2N −1
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣
(N
1 (N
2 (N
3
M
(N
N
)⎤
)
⎥ ⎥
)⎥
⎥
⎥
)⎥ ⎥⎦
=
⎡1⎤
⎢⎢0⎥⎥
⎢0⎥
⎢ ⎢
M
⎥ ⎥
⎢⎣0⎥⎦
o(Δx2 )
o(Δx10 )
界处理困难，效率不高。
3、有限体法（95年） 介于有限元和有限差分之间，具有有限元的剖分灵活
性，又具有差分的效率。但处理边界困难。
4、有限差分法（68年）----兼顾精度和效率，使用方便灵活，
Marmousi和盐丘模型均使用FD方法
四、有限差分法地震波传播数值模拟方法 常规差分法存在的问题
z 数值频散影响模拟精度，严重时造成计算溢出 z 人为边界反射处理不当影响模拟效果及模拟效率
=
A1
∂U ∂x
+B
∂U ∂z
右边界
上边界 左上角 左下角
∂U ∂t
=
A
∂U ∂x
+
B1
∂U ∂z
∂U ∂t
=
A1
∂U ∂x
+
B1
∂U ∂z
∂U ∂t
=
A1
∂U ∂x
+
B2
∂U ∂z
下边界 右上角 右下角
∂U ∂t
=
A2
∂U ∂x
+B
∂U ∂z
∂U ∂t
=
A
∂U ∂x
+
B2
∂U ∂z
∂U ∂t
=
A2
三、地震波场数值模拟方法综述
•Ray-tracing methods (Asymptotic methods) ------高频近似，适应简单模型和特定波型，振幅计算困难， 效率高
•Integral-equation methods (Based Huygen’s principle) •Volume integral equations •Boundary integral equations ------具有解析特征，只适用特殊模型，如除局部外介质均匀
∂x2 ∂y2 ∂z2 c2 ∂t2
组 ∂U = A ∂U + B ∂U + C ∂U
∂t
∂x
∂y
∂z
( 简写为ABC）。
左面A1BC 右面A2BC 前面AB1C 后面AB2C 底面ABC1 顶面ABC2
上左棱A1BC1 上右棱A2BC1 上前棱AB1C1 上后棱AB2C1 下左棱A1BC2 下右左棱 A2BC2 下前棱AB1C2 下后棱AB2C2 左前棱A1B1C 左后棱A1B2C 右前棱A2B1C 右后棱A2B2C
τ {τ β v } = 0.5 (old ) + ρ
(old )
xz
xz
2x
其中，1986年，Bayliss et al.在使用Mac Cormack分裂算法 模拟二维各向同性介质中的弹性波时就是使用的如上吸收边界。
1990年，D.Kosloff and D.Kessler在用Chebychev谱方法模 拟二维各向同性介质中的弹性波时也是使用的如上吸收边界，由 于仍存在边界反射，只好利用阻尼衰减法（Cerjan et al,1985） 再次吸收。
2k
2 x
ω2
⎟⎟⎠⎞2
⎪⎫ ⎬ ⎪⎭
∂3P ∂z∂t 2
−
v2 4
∂3P ∂x2∂z
+
3v 4
∂3P ∂x2∂t
−
1 v
∂3P ∂t 3
=
0
Elastic Wave: (Bottom Boundary)
Utt = α 2U xx + β 2U zz + (α 2 − β 2 )Wxz Wtt = β 2Wxx + α 2Wzz + (α 2 − β 2 )U xz
(15 Degree)
kz
≈
ω
v
⎧ ⎨1
−
⎩
1 2
v
2
k
2 x
ω2
⎫ ⎬ ⎭
∂2P + 1 ∂2P − v ∂2P = 0 ∂z∂t v ∂t 2 2 ∂x2
Third-order Approximation:
(45 Degree)
kz
=
ω
v
⎪⎨⎧1 − ⎪⎩
1 2
v
2
k
2 x
ω2
−
1 8
⎜⎜⎝⎛
v
断层下部界面 反射较右边弱
炮点位于17km处的单炮记录
断层下部界面 反射与右边相当
二、地质模型的描述
建模方法要和模拟方法相适应。
1、介质用层状或块状描述，输入各块边界信息 （折线顶点）以及介质参数。三维时按片输 入。------FE方法、射线方法。
2、建立模型的物性参数场。------FD方法。
eij
=
1 2
(ui,
j
+ u j,i )
σ = c e ij
ijkl kl
σ ij, j + fi = ρu&&i
Direct methods：
1、有限元法（80年代初） 不受边界几何形状的限制，具有灵活的剖分方法，效率低， 模拟精度有待提高。
2、傅立叶变换法（80年代中期） 利用傅立叶变换求波场的空间微分，空间精度高，边
−c 0 c
0 0 0
c ⎟⎞ 0⎟ − c⎟⎠
边界吸收前后波前快照
吸
收
边
界 条
边界吸收前
件
试
验
（
各
向
同
性 边界吸收后
介
质
）
X 分量
Z 分量
TI
吸 收 边 边界吸收前 界 条 件 试 验 （
介 质
边界吸收后
）
X 分量
Z 分量
3D吸收边界条件
3D吸收边界条件对三维地震波数值模拟显得尤为重要！
三维声波方程 ∂2u + ∂2u + ∂2u = 1 ∂2u 转化为一阶方程
采取的措施
z 高阶差分解法--提高计算精度，减小数值频散 z 采用基于特征分析方法得到的吸收边界条件
数值频散问题------高阶差分解法
声波：
∑ ∂2 f
∂x 2
=1 Δx 2
M
Cm(M )[ f (x + mΔx) − 2 f (x) +
m=1
f (x − mΔx)] + o(Δx2M )
弹性波：
数值频散试验dx=dz=10m, dt=1ms
高阶差分为何会消除数值频散？
∑ ∂2 f
∂x 2
=1 Δx 2
M
Cm(M )[ f (x + mΔx) − 2 f (x) +
m=1
f (x − mΔx)] + o(Δx2M )
2M阶空间差分数值频散关系：
∑ c2
c02
=1+
2(−1)M φ 2M
(2M + 2)!
ABCs Based on Damping Method
Pnew (i, j) = Pold (i, j) * exp(−σd 2 )
•d
ABCs Based on Characteristic Analysis
如2D弹性波方程可以分解为：
L L L L L L L L ∂U = (
∂t
M −1 (−)
c2 c02
≈1−
φ8 3150
其中，φ = kΔx = 2π Δx λ
数值频散理论分析小结
z c < c0，即空间离散造成的数值频散出现在正常波形之
后，作为一个“尾巴”出现。
z 对于相同频率子波，空间网格间距越大，数值频散越 严重；对于相同的空间网格间距，子波频率越高，数 值频散越严重。
z
只要
φ