有限差分法地震波场数值模拟中的吸收边界条件

图 3 无吸收边界条件波场快照图 a 为 520ms，b 为 580ms，c 为 680ms，d 为 800ms 时刻 图 3 所示的就是无吸收边界条件时二维声波波场快照图。a 是 520ms 时刻的波场快照， 此时波未传播到边界时的情形。b 是 580ms 时刻的波场快照，此时波刚好传播到人工边界。 c 是 680ms 时刻的波场快照，从图中可以看出四周都有很明显的因人工边界而引起的虚假反 射波。d 是 800ms 时刻的波场快照，此时边界反射就更加的明显了，能量几乎全部被反射回 来了。
1. 引言
在当今的中国乃至整个世界，经济的发展相当部分由工业决定，工业离不开石油业的发 展，石油业离不开地球物理勘探仪器与算法的进步，算法中最主要的是地震资料数据处理。 为了能更好地处理地震数据，掌握地震波传播规律，我们这时就要做好最根本的工作——地 震波场数值模拟。波场数值模拟（正演）方法有有限元法、克希霍夫积分法、射线追踪法、 波动方程有限差分法及以上方法的改进和综合等方法。由于有限差分法编程简单易实现，计 算速度快，占用内存少，模拟精度高，适用条件广而备受青睐。但是有限差分法模拟地震波 场时，计算机的内存不可能模拟无限大的无边界的地质模型，这时就考虑到人工边界问题。 下面就以二维声波波动方程为例，阐述人工边界反射问题（虚假反射）及其克服方法。
⎪⎪
− A2P(M − 2, j, k) + (2A −1)P(M , j, k −1) − 2AP(M −1, j, k −1)
⎨ ⎪P(i,
0,
k
+ 1)
=
(2
−
2
A
−
A2
)
P(i,
0,
k
)
+
2
A(1 +
A)P(i,1,
k
)
⎪ ⎪
− A2P(i, 2, k) + (2A −1)P(i, 0, k −1) − 2AP(i,1, k −1)
1 ∂2P + 1 ∂2P − 1 ∂2P = 0 v ∂n∂t v2 ∂t2 2 ∂τ 2
（6）
其中，n
为边界的外法线方向，τ
为边界的切线方向。我们用
1 v2
∂2P ∂t 2
−
∂2P ∂n2
替代
∂2P ∂τ 2
上式
可改写为
其差分格式为：
1 ∂2P + 2 ∂2P + ∂2P = 0 v2 ∂t2 v ∂t∂n ∂n2
USA: SEG, 1997.
-6-
http://www.paper.edu.cn
Overview on Absorbing Boundary Conditions of the Finite Difference Method for Seismic Wave Numerical Modeling
（7）
⎧P(0, j, k +1) = (2 − 2A − A2 )P(0, j, k) + 2A(1+ A)P(1, j, k)
⎪ ⎪
− A2P(2, j, k) + (2A −1)P(0, j, k −1) − 2AP(1, j, k −1)
⎪⎪P(M , j, k +1) = (2 − 2A − A2 )P(M , j, k) + 2 A(1+ A)P(M −1, j, k)
作者简介：陈敬国，男，1981 年 6 月生，硕士研究生，现主要从事地震波数值模拟及地震 资料偏移成像方法的学习与研究。
-7-
-3-
http://www.paper.edu.cn
4. 一阶吸收边界条件
从图 3 可以看出，如果没有吸收边界条件则会有很强的因人工边界造成的虚假反射，这
给地震数据的采集、处理及解释工作带来了诸多不便。由此我们急需解决因人工边界引起的
强虚假反射。这里我们采用Clayton_Engquist_majda一阶吸收边界条件[1]，其一阶形式为
可得
Байду номын сангаас
Pk +1 i, j
=
2Pi,kj
−
Pk −1 i, j
+
A2
(
Pk i +1,
j
+
Pk i−1, j
+
Pk i, j+1
+
Pk i, j−1
−
4 Pi ,kj+1 )
+
v2
(Δt)2
Src(k
)
（2）
式中，
Pk i, j
=
P(iΔx,
jΔz, kΔt)
，
Δx, Δz, Δt
为别为空间、时间离散步长，
-4-
http://www.paper.edu.cn
5. 二阶吸收边界条件
从图 4 中我们可以看出，尽管采用了一阶吸收边界条件，但还是有明显的边界虚假反射。 我们可以继续用更高一级的吸收边界条件，这样吸收效果就会更加的明显了。这里我们采用 的是Clayton_Engquist_majda二阶吸收边界条件[1]，二阶形式为
图 5 二阶吸收边界条件波场快照图 a 为 520ms，b 为 580ms，c 为 680ms，d 为 800ms 时刻 从图 5 可以看出：由于使用了二阶吸收边界条件，当能量传播到边界时，很好的被吸收 了，故 c、d 看不到因人工边界而产生的虚假反射。
6. 结论
通过亲自动手，独立编程模拟出了各种条件下的地震波场数值模拟，为地震资料的有限 差分法波动方程逆时偏移处理打下了坚实的基础。在用有限差分法解波动方程，然后模拟地 震波场的过程中，吸收边界条件处理是不可或缺的一个环节。吸收边界处理的本质就是当地 震波传到因人工造成的波阻抗时，不会产生回射波（backward wave）即人工边界变成无反 射的、透明的或者是吸收的。
参考文献
[1] 邵秀民,刘臻.带吸收边界条件的声波方程显式差分格式的稳定性分析.计算数学，2001，23（2）. [2] 孙若昧.地震波传播有限差分模拟的人工边界问题.地球科学进展，1996，11（3）：53-58. [3] R. P. Bording and L. R. Lines, Seismic modeling and imaging with the complete wave equation,
http://www.paper.edu.cn
有限差分法地震波场数值模拟中的吸收边界条件
陈敬国
中国地质大学（北京）地球物理与信息技术学院（100083）
E-mail: chenjg_cugb@163.com
摘要：有限差分法是地震波场数值模拟中最重要的方法之一，但是有限差分法模拟地震波场 时，计算机的内存不可能模拟无限大的无边界的地质模型，这时就考虑到吸收边界问题。本 文就以二维声波波动方程为例，通过对无吸收边界条件，一阶、二阶吸收边界条件效用的对 比，阐述人工边界反射问题（虚假反射）及吸收边界条件的重要性。 关键词：有限差分法，地震波场，数值模拟，吸收边界条件，伪反射
（5）
图 4 一阶吸收边界条件波场快照图 a 为 520ms，b 为 580ms，c 为 680ms，d 为 800ms 时刻 为了与图 3 相对比，图 4 中 a、b、c、d 均分别是 520ms、580ms、680ms、800ms 时刻 波场快照（下同）。从中可以看出：a、b 时刻波形未到达、刚到达边界，没有发生反射；c、 d 时刻波形正在通过边界，已经吸收了大部分能量，但仍有少部分能量发生了虚假反射。
当然，吸收边界有很多种方法如多次入射法、旁轴近似法、高吸收区法、完全匹配层法 及多种方法联立边界条件等[2]。在数值算法领域中，总存在这样一对矛盾：数值计算精度越 高，计算量就越大，计算效率就低。所以我们在处理实际问题时，根据实际问题所要求的目 标来确定计算精度，这样既提高了计算效率，又能完美地解决问题。
2. 模拟参数
二维标量声波波动方程（含震源）：
∂2P ∂x2
+
∂2P ∂z 2
+
f
(x, z)
=
1 v2
∂2P ∂t 2
（1）
其中： P 是声波波场 P(x, z, t) ， v 是声波速度 v(x, z) ， f (x, z) 是震源。
对（1）式进行时间、空间均 2 阶精度有限差分离散即 7 点差分格式（见图 1），整理后
A=
vΔt h
，
h = Δx = Δz ， Src(k) 为震源函数。
-1-
http://www.paper.edu.cn
震源函数[1]（波形见图 2）：
Src
=
⎧⎪sin(50t) ⎨
×
e−188(t −3π
/100)2
,
⎪⎩
0
,
0 ≤ t ≤ 6π /100 t > 6π /100
（3）
由于是计算机模拟，为了能说明问题且便于计算，我们设地质模型边界为 1，具体详细
∂P − 1 ∂P = 0 ∂n v ∂t
（4）
其中 n 为边界的外法线方向，差分格式为
⎧P(0, j, k +1) = AP(1, j, k) + (1− A)P(0, j, k) ⎪⎪P(M , j, k +1) = AP(M −1, j, k) + (1− A)P(M , j, k) ⎨⎪P(i, 0, k +1) = AP(i,1, k) + (1− A)P(i, 0, k) ⎪⎩P(i, N , k +1) = AP(i, N −1, k) + (1− A)P(i, N , k)
参数如下见表 1：
表 1 波场模拟参数一览表
横向采 模 型 边 模型边 空间步 空间步 时间步
样点数 界（X） 界（Z） 长（dx） 长（dz） 长（dt）
（Nx）
纵向采 样点数 （Nz）
时间采 样点数 （Nt）
声波波 速（v）
-2-
http://www.paper.edu.cn
1
1
0.01
0.01 0.002
⎪P(i, N , k +1) = (2 − 2A − A2 )P(i, N, k) + 2A(1+ A)P(i, N −1, k)
⎪
⎪⎩
− A2P(i, N − 2, k) + (2A −1)P(i, N , k −1) − 2AP(i, N −1, k −1) （8）
-5-
http://www.paper.edu.cn
Chen Jingguo
China University of Geosciences (Beijing) (100083) E-mail: chenjg_cugb@163.com Abstract
The Finite Difference Method (FDM) is one of the most important methods of seismic wave field numerical modeling. But when we model seismic wave field with FDM, absorption boundary problem should be considered because the machine can not model the geological body of infinite boundary. This paper presents 2-D acoustic wave equation sample, through contrasting effectiveness of non-absorbing with the 1st order and 2nd order absorbing boundary conditions, illuminate the artificial boundary reflection (pseudo-reflection) problem and show the importance of the absorbing boundary conditions. Key words: finite difference method, seismic wave field, numerical modeling, absorbing boundary conditions, pseudo-reflection
101
101
500
1
（表 1 中的参数均用于以下的无吸收边界、吸收边界条件的波场数值模拟）
3. 无吸收边界条件条件
由于地质模型的有限性，当地震波传播到地质边界时会产生因边界内外速度的巨大反差 （边界外的速度为零），从而产生了不是由目标体界面产生的反射，这就是虚假反射（伪反 射）。通过编程实现了地震声波波场快照。