当前位置:首页 >> 信息与通信 >>

地震声波数值模拟中的吸收边界条件


第 23 卷 第 4 期 2006 年 12 月

上海第二工业大学学报 JOURNAL OF SHANGHAI SECOND POLYTECHNIC UNIVERSITY

Vol.23 No.4 Dec.2006

文章编号: 1001-4543(2006)04-0272-07

地震声波数值模拟

中的吸收边界条件
邢 丽
(上海第二工业大学理学院,上海 201209)
摘 要:在地震波传播的数值模拟中,在有限区域建立吸收边界条件是一个很重要的问题。在不同类型的吸收边界 条件中,比较经典的是 Clayton 和 Engquist 提出的 Clayton-Engquist 吸收边界条件;PML 吸收边界条件是 Berenger 在 1994 年提出的一种新的人工边界条件。这种边界条件几乎达到零反射。在地震波传播数值模拟中,PML 已经被 用于速度应力形式的一阶声波方程, 而尚未被应用于位移形式的二阶声波方程。 推导了二阶位移形式的 PML 公式, 并应用于声波方程的数值模拟中;通过数值例子,与被传统应用的 CE 吸收边界条件作了比较。 关键词:地震声波数值模拟;有限差分;吸收边界条件;CE;PML 中图分类号:O24 文献标识码:A

0 引言
在地球物理、土木工程和电磁学等领域内,只要涉及到对局部区域内某种波动现象的解释,并在此解 释中需要考虑环境介质的影响,就会涉及到近场波动问题。为了实现波动问题的数值模拟,需要引入一个 包围所考虑区域的局部区域的有限计算区。在计算区域内用有限元或有限差分技术将对应的方程和物理边 界条件转化为离散的运动方程,从而实现区内波动现象的数值模拟。而计算区的边界是人为引入的,其上 的边界条件可称为吸收边界条件。 在地震波传播的数值模拟中,在有限区域建立吸收边界条件是一个很重要的问题。1977 年,Clayton 和 Engquist[1]对二维空间的声波方程提出了一系列不同阶数的吸收边界条件(简称 CE) 。PML(Perfectly 这种边界条件几乎达到零反射。 PML Matched Layer) Berenger[2]在 1994 年提出的一种新的吸收边界条件, 是 已经被用于速度应力形式的一阶声波方程[3],而没有被应用于位移形式的二阶声波方程。下面的目标是推 导二阶位移形式的 PML 公式。

1 吸收边界条件
1.1 Clayton-Engquist ABC 考虑直角坐标系 ( x, y ) 中的二维标准声波方程
2 ? 2u ? 2u ? 2?? u =c ? 2 + 2 ? ?t 2 ?y ? ? ?x

Clayton 和 Engquist 提出了一系列具有不同精度的吸收边界条件: 一阶近似

?u 1 ?u + =0 ?x c ?t

收稿日期:2006-04-00;修回日期:2006-09-26 作者简介:邢丽(1978—),女, 山东人, 硕士,主要研究领域为偏微分方程数值解。 基金项目:上海高校选拔培养优秀青年教师科研专项基金(No.LX306003)

第 23 卷

邢 丽:地震声波数值模拟中的吸收边界条件

273

? 2u 1 ? 2u c ? 2u 二阶近似 =0 ? + ?x?t c ?t 2 2 ?y 2
三阶近似

(1)

1 ? 3u 3 ? 3u ? 3u c 2 ? 3u =0 ? c + ? ?x?t 2 4 ?x?y 2 c ?t 3 4 ?t?y 2

这就是著名的 Clayton-Engquist 旁轴近似(paraxial approximation)吸收边界条件,简称 CE。 1.2 PML 吸收边界条件 对于二维标准声波方程, 其频率域形式的方程为

? ? 2u ? 2u ? ?ω u = c ? 2 + 2 ? ?y ? ? ?x
2 2

(2)

? ? 其中 ω 为频率,c 为波速。方程(2)解的形式为 A exp ?i (k ? x - ωt )? ,这里 A 为振幅, k = k x x + k y y, k x , k y 是 ? ?

? ? 笛卡儿坐标分量, x = xx + yy 。PML 方法的主要目标是建立一个新的波动方程,使其解的形式为
A exp [ i ( k ? x - ω t ) ? γ x ] , γ > 0 。这样,波沿着 x 轴方向呈指数递减状态,因此确保了介质和 PML 区域之 间的反射系数趋于零。为此,引入一个新的变量 x = x ? iγ ,这样就把方程(2)转化为以 x, y 为变量的方程

? ? 2u ? 2u ? ?ω 2u = c 2 ? 2 + 2 ? ?y ? ? ?x

(3)

把 x = 0 定义在规则区域与 PML 区域的交界处。在 PML 区域, x > 0 ,不失一般性。假设计算区域

? ? 与 PML 的区域交界线的法向量为 n ,定义 n 的增长方向为坐标轴 n(见图 1) 。

计算区域

PML 区域

0

n

? n
图 1 PML 区域示意图 Fig.1 The sketch of PML area

方程(3)的解的形式为 A exp ? i(k x x + k y y ? ωt ) ? = A exp [i( k ? x - ωt ) ? γ x ] ,其中 γ 的选择是关键的。首 ? ? 先定义一种衰减函数 d,使其在计算区域 d=0,在 PML 区域 d > 0,取 γ = 标n
1

ω



n

0

d ( s ) ds 。定义新的复数坐

274

上海第二工业大学学报

2006 年 第 4 期

n( n) = n ?

ω∫

i

n

0

d ( s ) ds

(4)

则有

?n iω = ?n iω + d (n)
以 x 方向为例,新的复数坐标为 ~ ,有 x

(5)

?x iω = 。 ?x iω + d ( x)

而式 (3) 中,

? ? ?x ? 2u ? ? ?u ? ? ? ?u ?x ? ? ? ?u iω iω ? 2u ?x ?u ? ? iω = ? ?= ? = ? = + ? ? ? ? 2 2 ?x ?x ? ?x ? ?x ? ?x ?x ? ?x ? ?x iω + d ( x) ? iω + d ( x) ?x ?x ?x ?x ? iω + d ( x) ? ?x
则式(3)变为
2 ?? ? ? 2u iω ω 2 d ′( x) ?u ? 2u ? ? ?ω u = c ? ? + + ? ? iω + d ( x) ? ?x 2 (iω + d ( x))3 ?x ?y 2 ? ? ? ? 2 2

(6)

令 u = u1 + u 2 + u3 ,则
2 2 ?? ? ? 2u ? ? 2u iω ω 2 d ′( x) ?u ? 2u ? iω 2 2? ? ? ω u1 = c ? ?ω (u1 + u2 + u3 ) = c ? ? + + ? 2 ? ? iω + d ( x) ? ?x 2 (iω + d ( x))3 ?x ?y 2 ? ? ? iω + d ( x) ? ?x ? ? 2 2



?ω 2u2 = c 2

ω 2 d ′( x) ?u (iω + d ( x))3 ?x
? 2u ?y 2

(7)

?ω 2u3 = c 2
把以上三式分别作拉普拉斯逆变换得

(? t + d ( x)) 2 u1 = c 2

? 2u ?x 2 ?u ?x
(8)

(? t + d ( x))3 u2 = ?c 2 d ′( x) ? t2u3 = c 2 ? 2u ?y 2

方程(8)就是我们前面期望得到的新的微分方程,其平面波解呈指数递减。对于其他的方向,可采用同样的 做法。在上式中取

第 23 卷

邢 丽:地震声波数值模拟中的吸收边界条件

275

? 1 ?? x ? d ( x) = log ? ? ? ? 2δ ? R ?? δ ?
3c
其中 δ 是 PML 区域的宽度,R 是反射系数,通常取 R= 10 导数,在数值计算中定义中间变量
?3[3]

2

。注意到 (8) 式中第二式左边是时间 t 的三阶

p = (? t + d ) u2
则 (8) 中第二式变为

(9)

(? t + d ( x)) 2 p = ?c 2 d ′( x)
式 (9) 可用 Euler 或 Runge-Kutta 方法求解。

?u ?x

2 数值解法和数值例子
2.1 数值解法 在本文中,对内部计算区域 ? ,有限元用双线性四边形元或双二次四边形元进行计算,差分用一般的 显格式。在对边界条件进行处理时,分别对上面讲述的两种人工边界条件用适当的数值方法进行离散,然 后通过数值例子对两种不同边界条件的计算结果加以分析和比较。 1)在数值例子中,由于二阶的 Clayton-Engquist ABC 有时会出现不稳定状态,所以在此采用一阶 Clayton-Engquist ABC

?u 1 ?u + =0 ?n c ?t
形式为

(10)

对边界 Γ 进行处理。式中 n 为边界法向量,用有限差分进行离散,如对于右边界,式 (10) 的差分离散

uin++1 j = uin+1, j ? c 1,

?t n (ui+1, j ?u i ,nj ) ?x

2)在此采用有限差分法对 PML 公式 (8) 进行离散。本文采用通常所用的五点差分格式。 2.2 数值例子 图 2 是均匀半空间点震源模型。Gauss 震源(*表示)加在计算区域 ? 的中间,在四条边界 Γ 上我们分 别使用 CE 吸收边界条件和 PML 吸收边界条件。在 R1 处放一接收器。取计算区域为正方形,区域的大小 取为 1 000 m×1 000 m,共划分处 200×200 个节点,取 x 方向和 y 方向的步长 ?x = ?y =5 m,时间步长为

?t = 0.000 125 s。介质参数 c = 2 000 m / s 。

276

上海第二工业大学学报

2006 年 第 4 期

图 2 震源模型 Fig.2 The model of hypocenter

经过计算机数值模拟计算,得到下面波场快照,见下列一系列图形。

图 2 a 波场快照图 Fig.2 a The wave snapshot

图 2 b 波场快照图 Fig.2 b The wave snapshot

图 2 c 波场快照图 Fig.2 c The wave snapshot

图 2 d 波场快照图 Fig.2 d The wave snapshot

第 23 卷

邢 丽:地震声波数值模拟中的吸收边界条件

277

图 2 e 波场快照图 Fig. 2 e The wave snapshot

图 2 f 波场快照图 Fig. 2 f The wave snapshot

图 2 g 波场快照图 Fig. 2 g The wave snapshot

图 2 h 波场快照图 Fig. 2 h The wave snapshot

图 2 a 至图 2 h 是分别利用上述两种吸收边界条件有限元合成的记录。其中图 2 a 至图 2 d 是波向外传播 时的波场快照图;图 2 e 至图 2 h 是波到达边界时分别加两种边界条件的情况;图 2 e 与 2 f,2 g 与 2 h 分别 是同一时刻加 CE 条件和 PML 条件。从图中可以看出,加 CE 边界条件后,波到达边界时会出现明显的反 射;而利用 PML 边界条件,波到达边界时,几乎没有发生反射。 为了更清楚详细地比较两种边界条件差异,我们分别画出接收器的时程图,并与精确解作比较。这里 用细网格高阶差分并通过远置人工边界的方法求出精确解。计算结果见下面的图 3。

a CE ABC

b PML ABC

图 3 接收器时程图(实线:精确解;虚线:数值解) Fig.3 The time-distance graph of receiver (real line:accurate solution, broken line: numerical solution)

278

上海第二工业大学学报

2006 年 第 4 期

图 3 中,图 a 是加 CE 边界条件,图 b 是加 PML 边界条件。实线是计算出的精确解,虚线是数值解。 从图中我们可以明显地看出,图 a 中两条曲线在后半部分相差很大,而图 b 中两条曲线几乎重叠。也就是 说,PML 边界条件明显地好于 CE 和 MTF 边界条件。

3 结论
本文重点讨论了两种吸收边界条件。针对地震波数值模拟中的声波方程,重点推导了位移形式的 PML 公式,并把 PML 和 CE 边界条件进行了比较。通过数值例子,我们可以概括出以下结论: 1. CE 条件是以微分形式给出的连续型人工边界条件,易于实现。但在实际运算中,两阶 CE 条件有时 会出现不稳定现象。一般使用一阶 CE 条件,但一阶 CE 条件的精度很低,对人工边界上的反射波不能很好 地吸收。另外,CE 条件对计算区域角点的处理也不是很理想。 2. PML 是一种连续型人工边界条件,优点是由所计算的波动方程导出。如果选择了适当宽度的 PML 区 域和阻尼函数,PML 对人工边界的反射波几乎全部吸收。并且以位移形式给出的 PML 公式易于程序实现。 而且对角点区域的处理也很容易,整个计算过程是稳定的。但 PML 也有缺点,比如对于声波方程来说,用 声波方程导出的 PML 区域上的方程共有三个,这样共有四个微分方程,因此,在程序实现时不免有些繁琐。 本文用对 PML 区域对应的三个微分方程用有限差分进行离散。其他的方法,比如有限元、谱方法也是可行 的。虽然程序实现比较困难,但 PML 人工边界条件的优越性还是大于它的缺点,还是值得我们应用的。

参考文献:
[1]CLAYTON, R, ENGQUIST, B. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bull. Seism. Soc. Am., 1977, 67: 1529–1540. [2]BERENGER, J P. Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[J]. J. Comput. Phys., 1994, 127: 363–379. [3]CHEW W C, LIU Q. Perfectly matched layers for the electrodynamics: a new absorbing boundary condition[J]. J. Comp. Acoustic, 1996, 4: 341–359. [4]邵秀民,蓝志凌. 非均匀各项同性介质中地震波传播的数值模拟[J]. 地球物理学报,1995, 38: 39–55. [5]SMITH W D. A non-reflecting plane boundary for wave propagation problems[J]. J. Comput. Phys., 1974, 15: 492–503. [6]KOMATITSCH D, TROMP J. A perfectly matched layer for the second-order seismic wave equation[J]. Geophys. J. Int., 2003, 154: 146–153.

Absorbing Boundary Conditions for the Numerical Simulation of Acoustic Waves XING Li
( School of Science, Shanghai Second Polytechnic University, 201209, P.R. China )
Abstract:In the numerical simulation of the wave equation, to develop an absorbing boundary condition in the finite domain is a very important problem. In different kinds of absorbing boundary conditions , a classical one is Clayton-Engquist (CE) ABC which developed by Clayton&Engquist 1977. Berenger(1994)introduced a new condition called the perfectly matched layer(PML),which has the remarkable property of having a zero reflection. Over the last decades, the PML has been used in the context of the acoustic wave equation formulated as a first-order system in velocity and stress, but hasn’t been used in the context of numerical schemes that are based on the wave equation written as a second-order system in displacement. The goal of this article is to reformulate PML condition in order to use it and to compare with CE condition. Keywords:The numerical simulation of acoustic waves;PFD;absorbing boundary condition;CE;PML


相关文章:
声波方程数值模拟实验报告
?t 2 声波方程的有限差分法数值模拟对于二维速度-深度模型,地下介质中地震波的...2uik, j ? uik, j ?1 (4-12) 考虑 Reynolds 边界条件(详细推到见文献...
声波方程数值模拟实验报告
理论知识需要的已知条件包括: 1.1)震源函数 2)地层速度(波速) 3)边界条件 2...声波方程的有限差分法数值模拟对于二维速度-深度模型,地下介质中地震波的传播规律...
王川-任意偶数阶精度交错网格粘声波方程数值模拟
任意偶数阶精度交错网格粘声波方程数值模拟 摘要 如何提高正演模拟的精度,使之更...差分格式, 结合完全匹配层(PML)吸收边界条件模拟地震波在各种模型中的传播...
二维有限差分法地震波场数值模拟研究
摘要: 有限差分法是地震波场数值模拟中最重要的方法之一, 但是有限差分法模拟...通过编程实现了地震声波波场快照。如果没有吸收边界条件则会有很强的因人工边界...
地震波数值模拟方法研究综述
性及地震波方程数值解法在地震数值模拟中应 用的...对于均匀各向同性介质中的声波方程而言,地球物 理...了一种透射加衰减 的组合人工边界方案(吸收边界条件...
声波加pml边界条件
因本文中只考虑各项同性介质中的地震波传播规律,故可做 ?? ? ?x ? ?z ...pml 吸收边界 暂无评价 13页 免费 声波方程数值模拟实验报... 暂无评价 28页...
OUC地震波场论数值模拟实验报告
OUC地震波场论数值模拟实验报告 隐藏>> 声波方程数值模拟试验报姓名: XXX 专业...时为没到达边界时的值,n 较大时为人工反射的波前值,第一部分中的两幅图为...
弹性波理论
地震波交错网格高阶差分数值模拟研究摘要: 地震波...Alford 等(1974)研究了声波方程 有限差分法模拟的...弹性波数值模拟中的吸收边界条件[J].石油地球物理...
第八章 边界条件
第八章 边界条件 任何数值模拟都可以认为仅仅是在物理区域或系统的一部分中进行...在人为规定的边界上 求解吸收流出波浪的方法有很多。 【10】—【15】 。我们...
数值模拟实验一
地震波场边界处理的MATL... 7页 免费 有限元分析...数值模型模拟实验报告 实验名称:地震记录数值模拟的...反射、绕射到 测线,传播经过中高频衰减,能量被吸收...
更多相关标签:
地震次声波 | 声波方程数值模拟 | 地震时会产生次声波 | 地震次声波前兆观测 | 声波吸收材料 声呐 | 声波吸收 | 地震发出什么声波 | 声波吸收材料 |