地震波场边界处理的MATLAB实现.pdfVIP

地震波场边界处理的MATLAB实现 陈敬国 中国地质大学（北京） 地球物理与信息技术学院 （100083） E-mail: chenjg_cugb@163.com 摘要：用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。但我们在实验室进 行波场数值模拟时有限差分网格是限制在人工边界里面，即引入了人工边界条件。本文采用 Clayton_Engquist_Majda二阶吸收边界条件，通过MATLAB编程实现了这一算法。依靠MATLAB具有 更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境，方便同行们交 流。 关键词：有限差分法，地震波场，数值模拟，吸收边界条件，MATLAB 1. 引言 [1] 用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法 。但我们在实验室进 行波场数值模拟时，只能在有限的空间进行，所以有限差分网格是限制在人工边界里面，即引入了 人为的边界条件。这种人为边界条件的引入将对有限区域内的波场值的计算带来严重影响，所以必 须进行特殊的边界处理。边界条件处理的好坏直接影响地震正演模拟的最终效果。本文中我们采用 [2] Clayton_Engquist_Majda二阶吸收边界条件 。 被称作是第四代计算机语言的MATLAB语言，利用其丰富的函数资源把编程工作者从繁琐的程 序代码中解放出来。MATLAB用更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的 程序开发环境。本文介绍运用MATLAB实现带有吸收边界条件的地震波场数值模拟方法和步骤，便于 同行们交流，亦可用于本科地震理论的教学中，让学生们在程序演示中理解地震波的传播规律。 2. Clayton_Engquist_Majda二阶吸收边界条件 我们给定二维标量声波波动方程（含震源）： （1） 式中：是声波波场，是声波速度， 是震源。 对（1）式进行时间和空间2阶精度有限差分离散（见图1），整理后可得 （2） 式中，， 为别为空间、时间离散步长， ， ，为震源函数。 震源函数： （3） Clayton_Engquist_Majda二阶吸收边界条件的微分表达式可参见文献[2]，其左、右、上、下 边界的差分格式分别为： 3. 基本算法步骤 从图1可以看出，k+1 时刻的波场值由k 时刻和k-1时刻的波场值决定。所以在MATLAB里实现的 基本算法步骤如下： (1) 初始时刻的全波场值均为零，P(i,j , dt)=0 （在MATLAB 中初始从dt开始，不能从0开始）； (2) 时刻2dt时，在炮点S (m, n)给出一个脉冲震源Src(t) （见式（3 ）），其它点波场P =0； (3) 时刻t大于或等于3dt时，P(i,j , k+ 1)根据式（2 ）计算，其它点波场P=0 ； (4) 在波传播到四周边界时，左、右、上和下边界的波场值分别由式（4-1