基于Matlab实现的地震波场边界处理软件姓名：姚嘉德学号：2015301130007院系：资源与环境科学学院摘要：用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。
但我们在实验室进行波场数值模拟时有限差分网格是限制在人工边界里面，即引入了人工边界条件。
本文采用Clayton_Engquist_Majda二阶吸收边界条件，通过MATLAB编程实现了这一算法。
依靠MATLAB具有更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境，方便同行们交流。
利用Matlab本身所具有可视化功能以及像素识别功能，可以将生成的动画电影进行识别，用于地震局实时分析有着深远意义。
关键词：有限差分法，地震波场，吸收边界条件，MATLAB矢量帧,像素识别Abstract：Modeling seismic wave field with the Finite Difference Method (FDM) is an effective method to study theseismic wave propagation in the earth medium. When we model seismic wave field in the laboratory, the finitedifference grids are restricted in the artificial boundary. So it should introduce the artificial boundary conditions. This paper adopts Clayton_Engquist_Majda second absorbing boundary conditions and realizes the arithmetic with MATLAB.The MATLAB codes are direct and accord with our thinking custom. So it can provide the friendlyand succinct programming environment and is easy to communicate with ing the functions of Matlab that make visualization come true and identify the pixel,we can identify the earthquake wave field.Key words: finite difference method, seismic wave field, numerical modeling, absorbing boundary conditions,MATLAB一、引言用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。
但我们在实验室进行波场数值模拟时，只能在有限的空间进行，所以有限差分网格是限制在人工边界里面，即引入了人为的边界条件。
这种人为边界条件的引入将对有限区域内的波场值的计算带来严重影响，所以必须进行特殊的边界处理。
边界条件处理的好坏直接影响地震正演模拟的最终效果。
本文中我们采用Clayton_Engquist_Majda二阶吸收边界条件[2]。
被称作是第四代计算机语言的MATLAB语言，利用其丰富的函数资源把编程工作者从繁琐的程序代码中解放出来。
MATLAB用更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境。
本文介绍运用MATLAB实现带有吸收边界条件的地震波场数值模拟方法和步骤，便于同行们交流，亦可用于本科地震理论的教学中，让学生们在程序演示中理解地震波的传播规律。
二.、Clayton_Engquist_Majda二阶吸收边界条件我们给定二维标量声波波动方程（含震源）：（1）式中：是声波波场，是声波速度，是震源。
对（1）式进行时间和空间2阶精度有限差分离散（见图1），整理后可得（2）式中，，为别为空间、时间离散步长，，，为震源函数。
震源函数：（3）Clayton_Engquist_Majda 二阶吸收边界条件的微分表达式可参见文献[2]，其左、右、上、下 边界的差分格式分别为：三、基本算法步骤从图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；(3) 时刻t 大于或等于3d t 时，P (i , j , k+1)根据式（2）计算，其它点波场P =0；(4) 在波传播到四周边界时，左、右、上和下边界的波场值分别由式（4-1）、（4-2）、（4-3）和（4-4）计算出来。
四、数值模拟由于是计算机模拟，为了能说明问题且便于计算，我们设地质模型边界为1，具体详细参数如下见表1：第一步：速度文件的载入及相关整理（替换）clc; clear; %清除工作空间及显示屏幕load vm_0.mat; % 载入速度文件，里面包含v(j, i)Nx=101; Nz=101; Nt=800;hx=0.01;hz=0.01;dt=0.002; % 模拟参数见表1for i=1:Nxfor j=1:Nzr(j,i)=v(j,i)*dt/hx;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); % 分布空间，预值充零第二步：赋初值for k=1:2for j=1:Nzfor i=1:Nxu(j,i,k)=0;endendend第三步：边界条件处理及7点差分计算波场延拓for k=2:Nt-1for j=1:Nzfor i=2:Nx-1if j==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,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 j==Nzu(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,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 j==30&i==51&(k-1)*dt <= 6*pi/100 %炮点S(30,51)，子波持续时间条件u(j,i,k+1)=c(j,i)^2*dt^2*sin(50*(k-1)*dt)*exp(-188*((k-1)*dt-3*pi/100)^2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+…u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);elseu(j,i,k+1)=r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);endendend ；for i=1:Nx 4.1 地质模型的构造本文选取了两个较为简单的地质模型，分别是均匀介质模型（见图2）和层状均匀介质模型（见图3）。
4.2 程序及相关说明根据上述我们建立的地质模型，生成相应的速度文件，再联合表 1中的模拟参数和吸收边界条件，就可以编程实现波场模拟。
平时表示波场的习惯u(x,z,t)对应在matlab 程序中，为了便于成图则被表示为u(z,x,t)，即u(i,j,k)变为u(j,i,k)。
详细过程如下：for j=2:Nz-1if i==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,i+1,k)-r2(j,i)*u(j,i+2,k)+(2*r(j,i)-1)*u(j,i,k-1)-…2*r(j,i)*u(j,i+1,k-1); % 左边界吸收elseif i==Nxu(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,i-1,k)-r2(j,i)*u(j,i-2,k)+(2*r(j,i)-1)*u(j,i,k-1)-…2*r(j,i)*u(j,i-1,k-1); % 右边界吸收elseif j==30&i==51&(k-1)*dt <= 6*pi/100 % 同上u(j,i,k+1)=c(j,i)^2*dt^2*sin(50*(k-1)*dt)*exp(-188*((k-1)*dt-3*pi/100)^2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+…u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);elseu(j,i,k+1)=r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);endendendend第四步：四个角点的处理for k=1:Ntu(1,1,k)=1/2*(u(1,2,k)+u(2,1,k));u(1,Nx,k)=1/2*(u(1,Nx-1,k)+u(2,Nx,k));u(Nz,1,k)=1/2*(u(Nz-1,1,k)+u(Nz,2,k));u(Nz,Nx,k)=1/2*(u(Nz,Nx-1,k)+u(Nz-1,Nx,k));end第五步：保存结果及成图u % 在主窗口显示数值结果，按照时间切片save U u ; % 保存数值结果% 波场快照图显示for k=1:Nt % 表示所有时刻，也可以是等间隔时间如k=1:10:Ntfigure(k)imagesc(u(:,:,k)) ; colormap(gray); % 同理可用surf(u(:,:,k))xlabel('x'); ylabel('z’');% zlabel('Amplitude') 在surf(u(:,:,k))中用到title('Wave Field Snapshot');axis squareend% 二维电影动画放映，除了imagesc,还有contour、surf等绘图命令clf; shg;M=moviein(Nt); %预设画面矩阵for k=1:Ntimagesc(u(:,:,k))colormap(gray)xlabel('x');ylabel('z');title('Wave Field Snapshot');axis([0 Nx 0 Nz]);M(:,k)=getframe; % 捕获画面end% movie(M,5,10) %%电影回放，每秒10帧，重复播放5次xin_xi_ti_shi=('*****************该程序已经运行结束********************')4.3 运行结果显示运行以上程序，可以得到所有时刻的波场快照图，在这里只选择部分图件进行展示。