⎧ ⎪ ⎨ ⎪(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)
+
Δz
∂τ zz
∂z
= τ zz (Δz) + τ zz (−Δz)
地震波正演模拟技术及其应用
董良国 同济大学地震波传播与成像学科组
2007.10.27厦门
内容
• 起伏地表弹性波传播数值模拟 • 地震波菲涅尔体与波路径的模拟 • 菲涅尔体层析成像 • 层析反演分辨率
一、起伏地表弹性波传播数值模拟
1、应用领域
（1）近地表地震波传播机制研究 （2）研究地形对地面震动的影响 （3）山地军事应用 （4）山地资料低信噪比问题 （5）起伏地表成像研究
σ zz
=
(
f
σ )' 2 xx
σ xz = f 'σ xx
（2）水平自由边界条件实现方法
即使是水平地表，要较好地实现自由边界条件也十分 不易。其中，
FE、BE具有优越性； PS几乎难以实现； FD为实现自由边界，必然降低整体模拟精度，且易造 成算法的不稳定；
• 隐式自由边界（Vidal and Clayton,1986）
FE
FD Horizontal Component(1000ms)
Vertical Component
Horizontal Component
（d）常规FD方法
（Masahiko Fuyuki，1980；董良国，2007）
二维空间网格点划分成22类（空间2阶差分精度），分别构造 不同的差分格式。
2
=0
τ zz (−Δz) = −τ zz (Δz)
• 密度镶边法（Frankel,1992）
−2Δz
0.16ρ
−Δz
0.4ρ
O
ρ
x
Δz
z
没有物理依据。
（3）起伏地表地震波模拟方法简介
声波方程无法模拟地震波在近地表的各种复杂传播现象！
选择何种地震波传播模拟方法，取决于研究的目的。 例如：研究地震波在近地表传播的各种复杂现象、利用初至波 旅行时或波形反演近地表速度结构、进行信噪比问题的研究，等 等，在模拟中均需要正确实施自由边界条件，这是模拟中的最关 键问题。
除纵向坐标变换外，其它方法均遇到过稳定性 问题！
二、地震波菲涅尔体与波路径的模拟
走时层析 Δt(g s) = ∫ Δw(rr)L[rr s, g]drr
将时差反投影至射线路径（Raypath）上。
Born近似 波动方程层析
Δψ
(g
s)
=
∫
ΔV (rr) V (rr)
⋅ξ
⎡⎣rr
s,
g,V
(rr)⎤⎦drr
⎪⎧
⎪⎩⎨η
ξ
=
=x zmax f (x)
z
(0 ≤ x ≤ xmax )
• 自由边界条件
σ zz
=
(f
σ )' 2 xx
σ xz = f 'σ xx
优点：统一形式隐式实现自由边界条件。
模型
水平分量
垂直分量
Vertical Component
Horizontal Component
Vertical Component
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝
lij
⎛
⎞⎛ ⎞ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟ ⎟
⎜ ⎜
Δs
j
⎟ ⎟
⇒
⎜ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟⎜ ⎟ ⎜
⎟⎠ ⎜⎝
⎟⎠
⎜ ⎜⎝
w lij1 ij1
∑ wij
L
w lij2 ij2
∑ wij
L
w lij3 ij3
∑ wij
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠
fΔt,
,
(0 ≤ (Δt
Δt ≤ 1 2aba,1999）
均匀介质20hz单频Born近似波路径 均匀介质20hz单频第一菲涅尔体
均匀介质50hz单频Born近似波路径 均匀介质50hz单频第一菲涅尔体
实部
虚部
均匀介质Born近似波路径（有限频）
等梯度模型以及对应的射线、50hz菲涅尔体
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）
等梯度模型50hz单频Rytov近似波路径
层状模型与透射波射线、30hz菲涅尔体
层状模型50hz单频透射波Rytov近似波路径
起伏地表等梯度模型示意图
炮检对SR对应菲涅尔体(50hz)
起伏地表层状模型示意图
炮检对SR对应菲涅尔体
复杂模型，Offset=1000 fresnel_volume
例如，第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,
ζ ⎡⎣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）上。
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、层析矩阵稀疏程度降低，地下各点的覆盖次数大幅度提高，反演的稳 定性提高。
+
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
• Schuster，Vasco，Watanabe等
– 层析在某一方向上所能分辩的最小空间尺度为 经过该点的所有菲涅尔体中在该方向上的最小 宽度。
s1
g1
S
s2
g2
s3
g3
R
f=30hz,offset=2000m
X方向分辨率 Z方向分辨率
在倾斜坐标系下，自由边界条件为：
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,
（Kravtsov and Orlov,1980;Woodward，1989,1992）