一种有效吸收边界条件的MATLAB实现

一种有效吸收边界条件的MATLAB 实现

陈敬国

中国地质大学（北京） 地球物理与信息技术学院 （100083）

E-mail: chenjg_cugb@

摘要：用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。但我们在实验室进行波场数值模拟时有限差分网格是限制在人工边界里面，即引入了人工边界条件。本文采用Clayton_Engquist_Majda 二阶吸收边界条件，通过MATLAB 编程实现了这一算法。依靠MATLAB 具有更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境，方便同行们交流。

关键词：有限差分法，地震波场，数值模拟，吸收边界条件，MATLAB

1. 引言

用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法[1]

。但我们在实验室进行波场数值模拟时，只能在有限的空间进行，所以有限差分网格是限制在人工边界里面，即引入了人为的边界条件。这种人为边界条件的引入将对有限区域内的波场值的计算带来严重影响，所以必须进行特殊的边界处理。边界条件处理的好坏直接影响地震正演模拟的最终效果。本文中我们采用Clayton_Engquist_Majda 二阶吸收边界条件[2]

。

被称作是第四代计算机语言的MATLAB 语言，利用其丰富的函数资源把编程工作者从繁琐的程序代码中解放出来。MATLAB 用更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境。本文介绍运用MATLAB 实现带有吸收边界条件的地震波场数值模拟方法和步骤，便于同行们交流，亦可用于本科地震理论的教学中，让学生们在程序演示中理解地震波的传播规律。

2. Clayton_Engquist_Majda 二阶吸收边界条件

我们给定二维标量声波波动方程（含震源）：

2222221(,)P P f x z 2P

x z v ∂∂∂++=∂∂∂t

（1） 式中：是声波波场，是声波速度，P (,,)P x z t v (,)v x z (,)f x z 是震源。

对（1）式进行时间和空间2阶精度有限差分离散（见图1），整理后可得

112,,,1,1,122,1,2( 4)()()

k k k k k

i j i j i j i j i j i j k k i j i j P P P A P P P P P v t Src k +−,1

k +−++−=−++++−+Δ （2）

式中，，,(,,k

i j P P i x j z k t =ΔΔΔ),,x z t ΔΔΔ为别为空间、时间离散步长，v t

A h

Δ=

，

- 1 -

h x =Δ=Δz /100

，为震源函数。

()Src k 震源函数：

2

188(3/100)sin(50),06/100

0,

6t t e t Src t πππ−−⎧×≤≤⎪=⎨

>⎪⎩

（3）

Clayton_Engquist_Majda 二阶吸收边界条件的微分表达式可参见文献[2]，其左、右、

上、下边界的差分格式分别为：

22

2(0,,1)(22)(0,,)2(1)(1,,) (2,,)(21)(0,,1)2(1,,1) (4-1)(,,1)(22)(,,)2(1)(1,,) P j k A A P j k A A P j k A P j k A P j k AP j k P M j k A A P M j k A A P M j k +=−−++−+−−−−+=−−++−2

2

2 (2,,)(21)(,,1)2(1,,1) (4-2)(,0,1)(22)(,0,)2(1)(,1,)

(,2,)(21)(,0,1)2(,1,1) (4-3)(,,A P M j k A P M j k AP M j k P i k A A P i k A A P i k A P i k A P i k AP i k P i N k −−+−−−−−+=−−++−+−−−−221)(22)(,,)2(1)(,1,) (,2,)(21)(,,1)2(,1,1) (4-4)A A P i N k A A P i N k A P i N k A P i N k AP i N k ⎧⎪⎪⎪⎪

⎪⎪⎨⎪⎪⎪

⎪+=−−++−⎪−−+−−−−−⎪⎩

3. 基本算法步骤

从图1可以看出，k +1时刻的波场值由k 时刻和k -1时刻的波场值决定。所以在MATLAB 里实现的基本算法步骤如下：

(1) 初始时刻的全波场值均为零，P (i , j , d t )=0（在MATLAB 中初始从d t 开始，不能从0开始）；

(2) 时刻2d t 时，在炮点S (m , n )给出一个脉冲震源Src (t )（见式（3）），其它点波场P =0；

- 2 -

(3) 时刻t大于或等于3d t时，P(i, j, k+1)根据式（2）计算，其它点波场P=0；

(4) 在波传播到四周边界时，左、右、上和下边界的波场值分别由式（4-1）、（4-2）、（4-3）和（4-4）计算出来。

4.数值模拟

由于是计算机模拟，为了能说明问题且便于计算，我们设地质模型边界为1，具体详细参数如下见表1：

表1 波场模拟参数一览表

模型边界 采 样 间 隔 总 采 样 点 数 声波波速 X Z dx dz dt Nx Nz Nt V1 V2

1 1 0.01 0.01 0.002101 101 500 1 2

4.1地质模型的构造

本文选取了两个较为简单的地质模型，分别是均匀介质模型（见图2）和层状均匀介质模型（见图3）。

4.2程序及相关说明

根据上述我们建立的地质模型，生成相应的速度文件，再联合表1中的模拟参数和吸收边界条件，就可以编程实现波场模拟。平时表示波场的习惯u(x,z,t)对应在matlab程序中，为了便于成图则被表示为u(z,x,t)，即u(i,j,k)变为u(j,i,k)。详细过程如下：

第一步：速度文件的载入及相关整理（替换）

clc; clear; %清除工作空间及显示屏幕

load vm_0.mat; % 载入速度文件，里面包含v(j, i)

Nx=101; Nz=101; Nt=800;hx=0.01;hz=0.01;dt=0.002; % 模拟参数见表1

for i=1:Nx

for j=1:Nz

r(j,i)=v(j,i)*dt/hx;

r2(j,i)=(v(j,i)*dt/hx)^2;

- 3 -