一种有效吸收边界条件的MATLAB实现.pdf-资料库

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

http://www.paper.edu.cn h x= Δ = Δz ， Src k ( ) 为震源函数。震源函数： Src ⎧ ⎪= ⎨ ⎪⎩ e t sin(50 ) × 0 − 188( t − 3 /100) π 2 , , 0 6 /100 π ≤ ≤ /100 6 > π t t （3） Clayton_Engquist_Majda 二阶吸收边界条件的微分表达式可参见文献[2]，其左、右、上、下边界的差分格式分别为： 2 , j k , A A P j k A j k A P P , ) 2 (1 ) (0, , ) ) (1, (2 2 (0, 1) − + + − + = P A j k j k AP A P (1, (0, (2, , 1) 1) 2 (2 , ) 2 − − − + − A P M A A P M j k A P M j k 1, (2 2 ( ( , ) 2 (1 ) 1) ) , ( , 2 + − − − + + = P M j k j k A P M A , 1) , ) 2, (2 ( ( 1) 2 , 2 + − − − − − k k A P i A A P i P i A ) ( ,0, ) 2 (1 (2 2 ) ( ,1, ) ( ,0, 1) 2 + + = − + − k A P i k k P i A AP i ( ,1, 1) ( ,2, ) 1) 2 ( ,0, (2 2 − + − − − P i N k A A A P i N k A P i N k , ( , 1) 1, ) ) ( , , ) 2 (1 ) ( , (2 2 2 + + + = − − − AP i N P i N k A k A P i N 2, ) 1) 2 , ( , 1) (2 ( , 2 − − + − − j k , − j k , ) AP M − k 1) (4 -1) ( − 1, j k , − 1) (4 - 2) − 1) (4 -3) ( , − 1, k − 1) (4 - 4) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ 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； - 2 -

(3) 时刻 t 大于或等于 3dt 时，P(i, j, k+1)根据式（2）计算，其它点波场 P=0； (4) 在波传播到四周边界时，左、右、上和下边界的波场值分别由式（4-1）、（4-2）、（4-3） http://www.paper.edu.cn 和（4-4）计算出来。 4. 数值模拟由于是计算机模拟，为了能说明问题且便于计算，我们设地质模型边界为 1，具体详细参数如下见表 1：表 1 波场模拟参数一览表模型边界采样间隔总采样点数声波波速 Z 1 X 1 dx dz dt Nx Nz Nt 0.01 0.01 0.002 101 101 500 V1 1 V2 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 -

http://www.paper.edu.cn s(j,i)=2-4*(v(j,i)*dt/hx)^2; end end % 简缩“常量” u=zeros(Nz,Nx,Nt); % 分布空间，预值充零第二步：赋初值 for k=1:2 for j=1:Nz for i=1:Nx u(j,i,k)=0; end end end 第三步：边界条件处理及 7 点差分计算波场延拓 for k=2:Nt-1 for j=1:Nz for i=2:Nx-1 if j==1 u(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==Nz u(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); else u(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); end end end for i=1:Nx for j=2:Nz-1 if i==1 u(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==Nx u(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)+… - 4 -

http://www.paper.edu.cn u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1); else u(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); end end end end 第四步：四个角点的处理 for k=1:Nt u(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:Nt figure(k) imagesc(u(:,:,k)) ; colormap(gray); % 同理可用 surf(u(:,:,k)) xlabel('x'); ylabel('z’'); % zlabel('Amplitude') 在 surf(u(:,:,k))中用到 title('Wave Field Snapshot'); axis square end % 二维电影动画放映，除了 imagesc,还有 contour、surf 等绘图命令 clf; shg; M=moviein(Nt); %预设画面矩阵 for k=1:Nt imagesc(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 运行结果显示运行以上程序，可以得到所有时刻的波场快照图，在这里只选择部分图件进行展示。 - 5 -

http://www.paper.edu.cn 图 4 是同一速度均匀介质模型的波场快照图。a、b、c 和 d 是分别 200ms、400ms、600ms 和 800ms 时刻的波场快照，震源位于模型中心，波速为 V1。从四个小图上可以看出在波逐步向外传播的过程中，当波传播到边界时逐渐被吸收从而没有人为的反射波。图 5 是二层均匀介质模型的波场快照图。a、b、c 和 d 是分别 200ms、400ms、600ms 和 800ms 时刻的波场快照，震源位于模型上中部，波速分别为 V1 和 V2。从四个小图上可以看出在波逐步向外传播的过程中，当波传播到地层边界时有很强的反射，而当波传播到人工边界时逐渐被吸收没有人为的反射波。 5. 结束语本文侧重于通过数值模拟介绍运用 MATLAB 实现带吸收边界条件的波场模拟的简要方法和基本步骤，因而文中数值模型和其计算差分节点设置都仅考虑十分简单的情形。在实际应用中，除计算节点数增多外，人们不断通过引入差分格式，差分精度等一系列措施改进正演模拟，从而提高计算效率，减少计算误差。以上数值模拟程序在 MATLAB7.0 中运行通过，可供同行们交流。对于二层均匀介质模型，运行的结果以 P(z, x)形式保存成 dat 文件，然后用地震成图软 - 6 -

http://www.paper.edu.cn 件打开，可以看到合成地震记录（以后另文再述）。对于时间切片的波场快照图件可以用 flash 软件做成 swf 文件，这样不需平台就可以单独显示动画了。在相关问题的探讨上，李晓光硕士提供了莫大的帮助，作者表示感谢。参考文献 [1] 孙若昧.地震波传播有限差分模拟的人工边界问题.地球物理学报，1996，11（3）. [2] 邵秀民,刘臻.带吸收边界条件的声波方程显式差分格式的稳定性分析.计算数学，2001，23（2）. [3] R.Clayton and B. Enquist, Absorbing boundary conditions for acoustic and elastic wave equations, Bull. Seis. Soc. Am., 67:6(1977), pp.1529-1540. [4] R. P. Bording and L. R. Lines, Seismic modeling and imaging with the complete wave equation, USA: SEG, 1997. [5] 张志涌,徐彦琴.MATLAB 教程——基于 6.x 版本.北京：北京航空航天大学出版社，2001，4. - 7 -

Realization of the Effectual Absorbing Boundary Conditions http://www.paper.edu.cn with MATLAB Chen Jingguo China University of Geosciences (Beijing) (100083) E-mail: chenjg_cugb@163.com Abstract Modeling seismic wave field with the Finite Difference Method (FDM) is an effective method to study the seismic wave propagation in the earth medium. When we model seismic wave field in the laboratory, the finite difference grids are restricted in the artificial boundary. So it should introduce the artificial boundary conditions. This paper adopts Clayton_Engquist_Majda 2nd 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 friendly and succinct programming environment and is easy to communicate with other. Key words: finite difference method, seismic wave field, numerical modeling, absorbing boundary conditions, MATLAB 作者简介：陈敬国，男，1981 年 6 月生。2004 年本科毕业于中国地质大学（北京），现在该校继续攻读硕士研究生学位，主要从事地震波数值模拟及地震资料处理方法的学习与研究。 - 8 -

资料库

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

相关推荐

开发技术

热门标签

最新资料