资源描述：
《地震波场边界处理的MATLAB实现》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、地震波场边界处理的MATLAB实现陈敬国中国地质大学（北京）地球物理与信息技术学院（100083）E-mail:chenjg_cugb@163.com摘要：用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。但我们在实验室进行波场数值模拟时有限差分网格是限制在人工边界里面，即引入了人工边界条件。本文采用Clayton_Engquist_Majda二阶吸收边界条件，通过MATLAB编程实现了这一算法。依靠MATLAB具有更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境，方便同行们交流。关键词：有限差分法，地震
2、波场，数值模拟，吸收边界条件，MATLAB1.引言[1]用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。但我们在实验室进行波场数值模拟时，只能在有限的空间进行，所以有限差分网格是限制在人工边界里面，即引入了人为的边界条件。这种人为边界条件的引入将对有限区域内的波场值的计算带来严重影响，所以必须进行特殊的边界处理。边界条件处理的好坏直接影响地震正演模拟的最终效果。本文中我们采用[2]Clayton_Engquist_Majda二阶吸收边界条件。被称作是第四代计算机语言的MATLAB语言，利用其丰富的函数资源把编程工作者从繁琐的程
3、序代码中解放出来。MATLAB用更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境。本文介绍运用MATLAB实现带有吸收边界条件的地震波场数值模拟方法和步骤，便于同行们交流，亦可用于本科地震理论的教学中，让学生们在程序演示中理解地震波的传播规律。2.Clayton_Engquist_Majda二阶吸收边界条件我们给定二维标量声波波动方程（含震源）：（1）式中：是声波波场，是声波速度，是震源。对（1）式进行时间和空间2阶精度有限差分离散（见图1），整理后可得（2）式中，，为别为空间、时间离散步长，，，为震源函数。震源函数：
4、（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
5、=0；(4)在波传播到四周边界时，左、右、上和下边界的波场值分别由式（4-1）、（4-2）、（4-3）和（4-4）计算出来。4.数值模拟由于是计算机模拟，为了能说明问题且便于计算，我们设地质模型边界为1，具体详细参数如下见表1：表1波场模拟参数一览表模型边界采样间隔总采样点数声波波速XZdxdzdtNxNzNtV1V2110.010.010.002101101500124.1地质模型的构造本文选取了两个较为简单的地质模型，分别是均匀介质模型（见图2）和层状均匀介质模型（见图3）。4.2程序及相关说明根据上述我们建立的地质模型，生成相应的速度文
6、件，再联合表1中的模拟参数和吸收边界条件，就可以编程实现波场模拟。平时表示波场的习惯u(x,z,t)对应在matlab程序中，为了便于成图则被表示为u(z,x,t)，即u(i,j,k)变为u(j,i,k)。详细过程如下：第一步：速度文件的载入及相关整理（替换）clc;clear;%清除工作空间及显示屏幕loadvm_0.mat;%载入速度文件，里面包含v(j,i)Nx=101;Nz=101;Nt=800;hx=0.01;hz=0.01;dt=0.002;%模拟参数见表1fori=1:Nxforj=1:Nzr(j,i)=v(j,i)*dt/hx
7、;r2(j,i)=(v(j,i)*dt/hx)^2;s(j,i)=2-4*(v(j,i)*dt/hx)^2;endend%简缩“常量”u=zeros(Nz,Nx,Nt);%分布空间，预值充零第二步：赋初值fork=1:2forj=1:Nzfori=1:Nxu(j,i,k)=0;endendend第三步：边界条件处理及7点差分计算波场延拓fork=2:Nt-1forj=1:Nzfori=2:Nx-1ifj==1u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j+1,i
8、,k)-r2(j,i)*u(j+2,i,k)+(2*r(j,i)-1)*u(j,i,k-1)-…2*r(j,i)*u(j+1,i,k-1);%上边界吸收elseif
