地震波正演模拟技术及其应用

复杂模型， Offset=2000 fresnel_volume
复杂模型， Offset=4000 fresnel_volume
50hz单频菲涅尔体
有限频菲涅尔体(主频50hz)
30hz单频菲涅尔体
有限频菲涅尔体(主频30hz)
有限频菲涅尔体(主频50hz) 有限频菲涅尔体(主频30hz)
R
S
S
R
R
S
S
R
Offset=2km
Offset=4km
不同炮检距的炮检对初至波所对应的菲涅尔体（50hz）
理论模型
射线层析 结果
表层精度显著提高 近似菲涅尔 体层析结果
层析结果与理论模型的比较
地表以下60米深度处的速度剖面图 理论（红）、射线层析（绿）、近似菲涅尔体层析（蓝）
利用不同层析结果进行波动方程基准面延拓后的炮道集
j
+
1)
+
D
r U
(i
+
1,
j)
−
r U
(i,
j
+ 1)
其中，
⎛⎜1− Vs2
D
=
C −1B1C
=
⎜ ⎜ ⎜⎜⎝
Vp2 − Vs2
Vp2
Vs2 ⎞⎟
Vp2 ⎟
Vs2 Vp2
⎟ −1⎟⎟⎠
x' x
z'
Vertical Component(500ms) Horizontal Component(500ms)
• Schuster，Vasco，Watanabe等
– 层析在某一方向上所能分辩的最小空间尺度为 经过该点的所有菲涅尔体中在该方向上的最小 宽度。
s1
g1
S
s2
g2
s3
g3
R
f=30hz,offset=2000m
X方向分辨率 Z方向分辨率
w
=
⎧ ⎪ ⎨
1− 2 rs ⋅ rg
fΔt , +1
0 ≤ Δt ≤1 2f ,rsg >>1
⎪⎩
0 , (Δt >1 2f )
无限高频时，退化为射线走时方程：τ
≈
∫ c (r )−1 L
dr
初至波菲涅尔体层析 在反演表层结构中的初步应用
与传统射线层析的差别与优势：
1、将Raypath上一点的投影值根据权值弥散到Fresnel带上; 2、First peak 取代 First break。 3、更符合地震波的传播实际，能更好地反演低速异常体，反演精度提高; 2、层析矩阵稀疏程度降低，地下各点的覆盖次数大幅度提高，反演的稳 定性提高。
Horizontal Component
问题： 地形的一阶和二阶坡度数值计算误差。
（b） FD+FE区域分裂法
（Moczo,1997; Huang and Dong,2004）
网格剖分示意图
0.4s垂直分量 波场快照
有限元方法 耦合方法
（c）交错网格FD方法
（Graves,1996；Robertsson,1996；王秀明,2004）
Vertical Component(900ms) Horizontal Component(900ms)
（e）零速度填充方法
（Zahradnik and Moczo,1993 ）
500m/s 800m/s
准噶尔南缘起伏地表模型
震源在16km处模拟记录
声米波泉方地程区模某拟实记际录记录
弹性波方程模拟的垂直分量记录
+
cos2 φσ zz
= σ zz − 2 f 'σ xz + ( f ' )2σ xx
1+ ( f ')2
σ nt
=
σ
' xz
=
− sin φ
cosφσ xx
+ cos(2φ )σ xz
+ sinφ cosφσ zz
[ ] = f ' (σ zz − σ xx ) + 1 − ( f ' )2 σ xz 1+ ( f ')2
等梯度模型50hz单频Rytov近似波路径
层状模型与透射波射线、30hz菲涅尔体
层状模型50hz单频透射波Rytov近似波路径
起伏地表等梯度模型示意图
炮检对SR对应菲涅尔体(50hz)
起伏地表层状模型示意图
炮检对SR对应菲涅尔体
复杂模型，Offset=1000 fresnel_volume
在倾斜坐标系下，自由边界条件为：
r U
'
z
+
r B1U
'
x
=
0
(i, j)•
•O
在O点时，差分形式为：
z
{ } r
U
' (i,
j)
=
r U
' (i
+1,
j
+1)
+
B1
r U
'
(i
+1,
j)
−Ur
' (i,
j
+1)
由于Ur '(i,
j)
=
r CU (i,
j)
，所以：
{ } r
U
(i,
j)
=
r U
(i
+ 1,
例如，第19类网格点：
r
r
水平自由边界条件： U z + B1U x = 0
r
r
A
垂直自由边界条件： U x + B2U z = 0
B
••(i, j)
O CD
{ } r
UH (i,
j)
=
r U
(i,
j
+
2)
+
dz dx
B1
r U
(i
+
1,
j +1) −Ur(i −1,
j +1)
=0
{ } r
UV (i,
j)
=
r U
(i
−
2,
j)
−
dx dz
B2
r U (i
−1,
j +1) −Ur(i −1,
j −1)
=0
求A、B点时采用
r UV
(i,
j)
，求C、D点时采用
r UH
(i,
j)。
( ) 求O点时，用
r Uxz
=
1 2
D+xD+z + D−xD−z
替Urxz = D0xD0z ，以免使用（i,j）点
例如，第21类网格点：（这里假定dx=dz）
射线层析 近似菲涅尔体层析
利用不同层析结果自起伏地表 进行波动方程叠前深度偏移结果
射线层析
菲涅尔体层析
实际资料处理
(塔里木沙漠资料) 绿山射线层析
同济近似菲涅尔体层析
实际资料处理
(四川山地资料) 绿山射线层析
同济近似菲涅尔体层析
四、层析分辨率
为什么研究层析分辨率？
1. 了解层析的反演能力 2. 优化观测系统设计 3. 辅助模型离散化 4. 辅助平滑窗设计
垂直镜像
水平镜像
水平+垂直镜像
水平+垂直镜像
Z分量
Vertical Component
R
R
X分量
Horizontal Component
R
R
Vertical Component(1000ms) Vertical Component(1200ms)
FE
FD Horizontal Component(500ms)
σ zz
=
(
f
σ )' 2 xx
σ xz = f 'σ xx
（2）水平自由边界条件实现方法
即使是水平地表，要较好地实现自由边界条件也十分 不易。其中，
FE、BE具有优越性； PS几乎难以实现； FD为实现自由边界，必然降低整体模拟精度，且易造 成算法的不稳定；
• 隐式自由边界（Vidal and Clayton,1986）
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝
lij
⎛
⎞⎛ ⎞ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟ ⎟
⎜ ⎜
Δs
j
⎟ ⎟
⇒
⎜ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟⎠ ⎜⎝
⎟⎠
⎜ ⎜⎝
w lij1 ij1
∑ wij
L
w lij2 ij2
∑ wij
L
w lij3 ij3
∑ wij
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠
⎪⎧
⎪⎩⎨η
ξ
=
=x zmax f (x)
z
(0 ≤ x ≤ xmax )
• 自由边界条件
σ zz
=
(f
σ )' 2 xx
σ xz = f 'σ xx
优点：统一形式隐式实现自由边界条件。
模型
水平分量
垂直分量
Vertical Component
Horizontal Component
Vertical Component
2、模拟方法
• 起伏地形对地震波传播的影响的研究，50年代 开始，至今仍是研究的热点，特别是在天然地 震和山地地震勘探研究中。
• 起伏地表地震波模拟方法仍未成熟！
（1）主要困难：
----自由边界条件数值实现如何更精确、更稳定？
σ nn
=
σ
' zz
=
sin 2 φσ xx
− sin(2φ )σ xz
fΔt,
,
(0 ≤ (Δt
Δt ≤ 1 2
>1 2f )
f
)
（Watanaba,1999）
均匀介质20hz单频Born近似波路径 均匀介质20hz单频第一菲涅尔体
均匀介质50hz单频Born近似波路径 均匀介质50hz单频第一菲涅尔体
实部
虚部
均匀介质Born近似波路径（有限频）
等梯度模型以及对应的射线、50hz菲涅尔体
⎧ ⎪ ⎨ ⎪(1
−
∂U ∂z
2β 2
+ ∂W ∂x
) ∂U +
=0 ∂W
=0
⎩ α 2 ∂x ∂z
O
x
z
• 镜像方法（Levander,1988;Crase,1990;Robertson,1996）
−Δz
O
x
Δz
z
自由面以上波场充零不正确！
τ
zz
(0)
=
τ
zz
(−Δz)ຫໍສະໝຸດ Baidu
+
Δz
∂τ zz
∂z
= τ zz (Δz) + τ zz (−Δz)
FE
FD Horizontal Component(1000ms)
Vertical Component
Horizontal Component
（d）常规FD方法
（Masahiko Fuyuki，1980；董良国，2007）
二维空间网格点划分成22类（空间2阶差分精度），分别构造 不同的差分格式。
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝⎜⎜
Δs j1 M
Δs j2 M
Δs j3
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠⎟⎟
2000m/s
模型试验
3500m/s
4000m/s 4200m/s
5500m/s
4700m/s
2500m/s
模 拟 炮 记 录
模 拟 炮 记 录
ray
Fresnel
初至波菲涅尔体层析反演表层速度
ζ ⎡⎣rr s, g⎤⎦ = 2k02(rr)G0 [rr, g]⋅ψ0(rr s)
Rytov近似 波动方程层析
ΔΦ(g s) = ∫ ΔVV((rrrr)) ⋅ζ [rr s, g]drr
[ ] ζ
rr s, g
=
2k02
(rr)
G0
[rr, g]⋅ψ0(rr
ψ0(g s)
s)
将剩余波场反投影至波路径（Wavepath）上。
理论模型 走时层析 波形层析
2500 1958 1416 873 331 (m/s)
Jianming Sheng （2004）
菲涅尔体走时层析
∫ ∫ 菲涅尔体走时： τ ≈ d (r)−1 w(r,σ )c(r,σ )−1dσdr
L
D(r )
d (r) = ∫D(r) w(r,σ )dσ
( ) 权值的选择,如：
除纵向坐标变换外，其它方法均遇到过稳定性 问题！
二、地震波菲涅尔体与波路径的模拟
走时层析 Δt(g s) = ∫ Δw(rr)L[rr s, g]drr
将时差反投影至射线路径（Raypath）上。
Born近似 波动方程层析
Δψ
(g
s)
=
∫
ΔV (rr) V (rr)
⋅ξ
⎡⎣rr
s,
g,V
(rr)⎤⎦drr
（Kravtsov and Orlov,1980;Woodward，1989,1992）
菲涅尔体：
τ SP
+ τ PR
− τ SR
≤
1 2f
分别计算炮点S、检波点R在P点的旅
行时τ SP 、τ RP ，如果 τ SP、τ RP满足上
式，则P为菲涅尔体中的一点。
灵敏度函数选择：
w
=
⎧1− 2
⎨ ⎩
0
三、菲涅尔体层析成像
走时层析
基于高频近似, 沿高速区域优势采 样，反演精度不高,只能反演速度空 间变化的低波数成分。走时层析矩 阵稀疏，病态严重。 。
菲涅体层析
波动方程层析
（Wu and Toksoz,1987; Woodward,1992）
• 子波问题 • 低信噪比问题 • 实际地震波传播的准确描述问题 • 反演的强非线性问题 • 初始模型问题
地震波正演模拟技术及其应用
董良国 同济大学地震波传播与成像学科组
2007.10.27厦门
内容
• 起伏地表弹性波传播数值模拟 • 地震波菲涅尔体与波路径的模拟 • 菲涅尔体层析成像 • 层析反演分辨率
一、起伏地表弹性波传播数值模拟
1、应用领域
（1）近地表地震波传播机制研究 （2）研究地形对地面震动的影响 （3）山地军事应用 （4）山地资料低信噪比问题 （5）起伏地表成像研究
2
=0
τ zz (−Δz) = −τ zz (Δz)
• 密度镶边法（Frankel,1992）
−2Δz
0.16ρ
−Δz
0.4ρ
O
ρ
x
Δz
z
没有物理依据。
（3）起伏地表地震波模拟方法简介
声波方程无法模拟地震波在近地表的各种复杂传播现象！
选择何种地震波传播模拟方法，取决于研究的目的。 例如：研究地震波在近地表传播的各种复杂现象、利用初至波 旅行时或波形反演近地表速度结构、进行信噪比问题的研究，等 等，在模拟中均需要正确实施自由边界条件，这是模拟中的最关 键问题。
几种起伏地表弹性波传播模拟方法
（a）纵向坐标变换方法 （b）FD+FE区域分裂法 （c）交错网格FD方法
（自由边界实现采用镜像法） （d）常规FD方法 （e）零速度填充方法 （f）边界元法
（a）纵向坐标变换方法
(Tessmer and Kosloff,1992; Hestholm and Bent Ruud,1994)