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

1、陈敬国中国地质大学（北京）地球物理与信息技术学院（ 100083 ）cheiyg_cugb地震波场边界处理的MATLAB实现摘要：用有限差分法模拟地震波场是研究地変波在地球介质中传播的有效方法.但我们在实脸室进 行波场数值模拟时有限差分网格是限制在人工边界里面，即引入了人工边界条件.本文采用 Clayton_Eiigquist_Majda二阶吸收边界条件，通过MATLAB编程实现了这一算法.依靠MATLAB具有 更加直观的、符合大仗思维习惯的代码，为用户提供了友好、简洁的程序开发环境，方便同行们交 流.关虢词：有限差分法，地震波场，数值模拟，吸收边界条件，MATLAB1.引言用有限差分法模拟地

2、震波场是研究地爲波在地球介质中传播的何效方法I但我们在实验电进 仃波场数值模拟时，只能在仃限的空间进行，所以仃限差分网格是限制在人工边界里面，即引入了 人为的边界条件。这种人为边界条件的引入将对冇限区域内的波场值的计算带来严匝影响，所以必 须进行特殊的边界处理。边界条件处理的飓坏口接影响地震正演模拟的赧终效果。本文中我们釆用 C layton_Engquist_Maj da 二阶吸收边界条件口。被祐作是第皿代计算机语言的MATLAB语言，利用其丰富的函数资源把编程匚作者从繁琐的程 序代码中解放出来.MATLAB用更加直观的、符合人众思维习惯的代码，为用户提供了友好、简洁的 程序开发坏境。本文介

3、绍运用MATLAB实现带何吸收边界条件的地雇波场数值模拟方法和步骤，便 同行们交流，亦可用本科地震理论的教学屮，让学生们在程序演示屮理解地震波的传播规律。2. Clayton_Engquist_Majda二阶吸收边界条件我们给定二维标杲声波去动方程（含震源）：涙 dz1 八一（加式中：是声波波场，是声波速度，*vvx,z）/（x?z）l是震源。对（1）式进行时间和空间2阶精度冇限差分离散（见图1）,整理后可得+昭厂4用门+円&几氐伙）(2)vAr式中，巴产PQMjgkW，共为别为空间、时间离散步长，h力=& = & ,为靈源隨数。皿如震源函数:Src =sin(50f)xe-188(|-3jr

4、/loo),0r 6zr/100图1波动方程差分1&式示意图Clayton_Engquist_Majda 阶吸收边界条件的微分表达式可参见文献2,其左、右、上、卜 边界的差分格式分别为：卩(0,护 + 1) (2 _ 2 A -j, A) + 24(1 + 4)如 j. A)-QPQ, Z灼 + (21)P(0,2 -1)_2AP(1,(4-1)P(MtJ,k+)=(2-2A- A2)P(M, 7，i)+ 24(1 +k)-A2P(M- 2, J,町 + (2 A J： & -1) _ 2 AP(M 1, J,上一1) (4 - 2)PQ, 0, + l) = (2-24-夕讯y)亠 24(1

5、 + 4)P(?, 1, Q-A2P(i, 2tk)(2A- 1)卩伉0北 一 1)一 2APQJ上-1)(4 3)卩皿上十1)=(2-2虫一护)7川比)+2人(1十A)PQ,N-g.-PFQ,N 一 2,町十(2 卫一 )P(i, N,k-l)-2AP(i9N-l,k-T)(4 -4)3. 基本算法步骤从图1可以看出,R+1时刻的波场值由R时刻和hl时刻的波场值决定.所以在MATLAB里实现的 某本算法步骤如卜 I(1) 初始时刻的全波场值均为零,P(iJ, d/尸0 (在MATLAB中初始从d/开始,不能从0开始)；(2) 时刻2d川寸，在炮点S(m,“)给山-个脉冲震源Src(t)(见式

6、(3),其它点波场P =0：(3) 时刻f大F或等丁3df时,P(i,j,k+1)根据式(2)计算，其它点波场P=0：(4) 在波传播到四周边界时，左、右、上和卜边界的波场值分别由式(4-1) . (4-2) . (4-3) 和(4-4)计算出來。4. 数值模拟由于是计算机模拟,为了能说明问题且便于计算,我们设地质模型边界为1,具体详细参数如 卜见表1:表1波场模拟参数一览表模型边界采样间隔总采样点数声波 波速XZdxdzdtNxNzNtVIV2110.010.010. 002101101500124.1地质模型的构造木文选取了两个较为简单的地质模熨，分别是均匀介质模熨(见图2)和层状均匀介质

7、模粮(见 图3) 图2 层均匀介质模型图3二层均匀介质梗型4.2程序及相关说明根据上述我们建立的地质模型，生成相应的速度文件，再联介农1屮的模拟参数和吸收边界条 件，就可以编程实现波场模拟。平时表示波场的习惯u(x,z,t)对丿在niatlab程序屮，为了便/成图则 被表示为U(zxt), UPu(ijJc)变为u(j,ijc)o详细过程如F：笫-步：速度文件的载入及相关整理(替换)clc; clear; %淸除工作空间及显示屏胖load vm_0 mat; % 入速度文件.里而包禽v(j.i)Nx-101, Nz-101. Nt-800;hx-0 01Jiz-0.01,dt-0.002; %

8、 模拟参数见农 1for i-l :Nxfdij-l:NzKM 尸vO，Ldhx:r204MvU4)*dt/hx)A2; s(j.i)-2-4w(v(ja)dt.hx)A2; endend %简缩喘费u-zeros(Nz.Nx.Nt); %分布空间.预值充零第二步：赋初值for k-1 2fdij-l:Nz for i-l:Nxu(j4.k)-0.endendend第三步：边界条件处理及7点差分计算波场延拓fork-2Nt-lfoij-l:Nzfor i-2:Nx-lif j1u(j.tkl)-(2-2*r04)-r2(j.i)*u0j.k)+2*rg4)*(Rr(j.i)*u(j+l j.k

9、H2Qd)uQ+2Jjc-l)-.2*r(j4)*uO+l.iJc-l);% 上边界吸收elseif jNzu0xklM2-2*iU4br204)*u0a.k)+2*r0j)*(l+r04)*uu44JcH20j)*u(j-2,iJc)n2*r(ja)-l)u(jaJc-l)-.%下边界吸收elseif j30&i51&(k-l)dt l); % 右边界吸收elseifj30&i51&(k-ldt -6*pvlOO % 同上 u(jxk-t-l)-c(j4)A2*dtA2*siii(50*(k-l)*dt)*exp(-188(k-l),dt-3*pL100)A2)-t-r2(j4)*(u(jj

10、+Lkrt-u(j 上 l.k)+ uglj,k)-ni(j-l,i.k)+s(j4)*u(j,iJc)-u(j,iJc-l);elseu(j.hkl)-r2(j4)*(u(j.计 lk)+u(jiJk)+u(j+liJc)rUJik)+s(j)5(jkHi0jkl)：endendendend第四步:四个角点的处理for k-1 Ntu(l,l .k)-l/2 (u( 12 Jpr(2,1u(l ,Nx.k)-l/2*(u( 1 Nx 1 .k)十 u(工 Nx 上);u(Nz,l ,k)-l (u(Nz-l ,1 .k)u(N 乙 2、k); u(Nz.Nx.k)-l/2*(u(N z.Nx

11、-1 ,k)+u(Nz-1 NxJ:);end第五步：保存结采及成图u%在左窗口显示数值结果.按照时间切片save U u ;%保存数值结果%波场快照图显示for k-1 Nt %表示所有时刻.也可以是等间隔时间jOI k-l:10:Ntfiguie(k)inngesc(u(:,:,k); colorinap(gray); % 同理可用surf(u(:,：k) xlabel(*x*); ylabef); % zlabel(fAmplitude1)在surf(u(:,:J0)中用到 utle(*Wave Field Snapshot*);axis squareend% 维电蹈动ilHj放映除门m

12、agesc,还/(contour, surf绘图命令elf. shg.M-inoviein(Nt); % 预设刚 |fli 矩阵for k-1 Ntimagesc(u(:colormap(giay)xlabel(,x,);yhbelCz,);titAVave Field Snapshot*);axis(0 Nx 0 Nz);M(: ,k)-gecfraine; % 捕获 R |fijend% movie(M.5J0)%电影|1放，毎秒10帧.巫复播放5次xmXi-LshLr” 该程序己经运行结束 4.3运行结果显示运行以I:程序町以得到所有时刻的波场快照图.在这里只选择部分图件进行展示。20-

13、80-Wave Field SnapshotWave Field Snapshot204060801001.20 40 60 E0 100x1001_!.I20 40 60 80 100x目4均匀介质波场快毀因图4足同速度均匀介质模型的波场快照图。a、b、cflk/Jii分别200ms、400ms、600ms和800ms 时刻的波场快照，震源位于模型中心，波速为VI。从四个小图匕可以看出在波逐步向外传播的过程 中，当波传播到边界时逐渐彼吸收从而没有人为的反射波。图5是二层均匀介质模型的波场快照图。a. b. c和d是分别200ms、400ms. 600ms和800ms时刻 的波场快照，震源位J

14、模型上屮部，波速分别为VI和V2。从四个小图上可以看出在波逐步向外传 播的过程中，半波传播到地层边界时有很强的反射，而当波传播到人匸边界时逐渐被吸收没冇人为 的反射波。5.结束语木文侧朿通过数值模拟介绍运用MATLAB实现带吸收边界条件的波场模拟的简要方法和基木 步骤，因而文中数值模型和其计算差分节点设置都仅考虑十分简单的情形。在实际应用中,除计算 节点数増多外.人们不断通过引入差分格式差分精度等一系列措施改进正演模拟.从而提高计算 效率，减少计算谋差。以上数值模拟程序在MATLAB7.0中运行通过，可供同彳J：们交流。对J：二层均匀介质模怕 运行的结果以P(S)形式保存成dat文件，然后用地

15、震成图软204060Wave Field Snapshotx80inn Iiiiit20 40 60 80 10005二层均匀介质波歧快照图件打开，可以看到介成地震记录(以厉另文再述)o对J：时间切片的波场快照图件可以用flash软 件做成swf文件，这样不需平台就可以单独显示动画了。在相关问题的探讨上，李晓光硕士提供了莫人的於助，作者表示感谢。参考文献1孙若昧地艇波传播有限基分模拟的人工边界问题地球物理学报.1996, 11 (3) 邵秀氏刘臻带吸收边界条件的声波方程显式足分格式的稳定性分析计算数学.2001. 23 (2) 3 R Clayton and B Enquist. Absorb

16、ing boundaiy conditions for acoustic aiid elastic wave equations. Bull Seis. Soc. Am . 67 6(1977), pp.1529-1540.R. P. Boidmg and L. R Lines, Seisimc modeling and imaging with the complete wave equation, USA: SEG, 1997.5张志涌徐彦琴.MATLAB教程搭于6.x版本北京：北京触空航天大学出版社.2001, 4.Realization of the Effectual Absorbi

17、ng Boundary Conditions withMATLABChen JingguoChina Umversitv of Geosciences (Beijmg) (100083)E-mail: chenjg_cugb AbstractModeling 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 seisnuc wave fie

18、ld in the laboratory, the finite difference grids are restricted in the artificial boundary So it should introduce the artificial boundary conditions, ndThis paper adopts Clayton_Engquist_Majda 2 absorbing boundaiy conditions and realizes the arithmetic with MATLAB The MATLAB codes are direct and accord with our thinking custom So it can provide