有限差分纵波波场模拟(7)

2018-12-12 22:32

[19]奚先，姚姚.二维随机介质及波动方程正演模拟[J].石油地球物理勘探，2001，36（5）：546～552 [20]董良国等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报，2000，43（3）:411～419 [21]奚先，姚姚.二维弹性随机介质中的波场特征[J].石油地球物理勘探，2004，39（6）：679～685 [22]邵秀民，蓝志凌．非均匀各向同性介质中地震波传播的数值模拟[J].地球物理学报，1935；38(1)：39～55

[23]孙若昧.地震波传播有限差分模拟的人工边界问题[J].地球物理学进展，1996，11（3）：53 [24] Cerjan C, Kosloff D, Kosloff R and ReshefM. A nonreflecting boundary condit -ion for discrete acoustic and elastic wave equation. Geophysics, 1985, 50 (4) : 705～708

[25]Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equation.Bulletin of the Seismological Society of America, 1977;66(11):1529~1540

[26]Reynolds. Boundary conditions for the numerical solution of wave propagation problem.Geophysics, 1978;43(6):1099~1110

[27]崔兴福，张关泉.地震波方程人工边界的两种处理方法[J].石油物探，2003，42（4）：452～455 [28]张毅.声波正演模拟中的人工边界问题[J].工程地球物理学报，2007，4（4）

[29]李文杰，魏修成，刘洋.声波正演中一种新的边界条件-双重吸收边界条件[J].石油物探，2004，43（6）：528～531

[30]邢丽.地震波数值模拟中的吸收边界条件[J].上海第二工业大学学报，2006，23（4）

[31]Dahlain M A. The application of high difference to the scalar waveequation. Geophysics, 1986;51(1): 54~66

[32] Higdon R L. Absorbing boundary condition for elastic waves. Geophysics,1991, 56 (2) :231～241 [33]孙成禹，张吉辉．完全纵波方程有限差分数值模拟[J].石油地球物理勘探，2005；40(3)：289～294

[34]黄自萍，徐伶俐，周继顺.起伏地表声波方程的数值模拟[J].石油地球物理勘探，2006，41（3）：275～280

[35]宛书金，董敏煌等.横向各向同性介质中的地震波传播特性[J].石油大学学报，2002,26(1)：23~28

附录

本课程报告所编matlab程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 2015.5.20 纵波波场模拟

%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 初始化 mode=12; %选择差分介数，最高支持16介 Model=load('Model.mat');

V_model=Model.UniformModel;%速度模型 c=load('c.mat');

c=c.c; %系数矩阵 cellnum=size(V_model); %网格数目 dt=100*10^(-6); %时间采样间隔 100us dx=1; dz=1; %空间采样间隔 x0=250; z0=20; %震源坐标 f0=30; %震源子波频率

S_t=2*10^(-3); %震源持续时间 2ms T=140*10^(-3); %波场快照时间 %%%%%%%%%%%%%%%%%%%%%%%%%%%% U0=zeros(cellnum(1),cellnum(2)); %波场初始值 U1=zeros(cellnum(1),cellnum(2)); U2=zeros(cellnum(1),cellnum(2)); for k=1:T/dt; if k

U1(z0+1,x0+1)=sin(2*pi*f0*k*dt)/(k*dt);%震源函数。 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if mode==2

for m=mode/2+1:cellnum(2)-mode/2 % 2介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+... 2*U1(m,n)-U0(m,n); end end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if mode==4

for m=mode/2+1:cellnum(2)-mode/2 % 4介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+... c(mode/2,2)*a^2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n))+... 2*U1(m,n)-U0(m,n); end end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if mode==6

for m=mode/2+1:cellnum(2)-mode/2 % 6介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if mode==8

for m=mode/2+1:cellnum(2)-mode/2 %8介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

c(mode/2,4)*a^2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n))+... 2*U1(m,n)-U0(m,n); end end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if mode==10

for m=mode/2+1:cellnum(2)-mode/2 % 10介差分方程

for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if mode==12

for m=mode/2+1:cellnum(2)-mode/2 % 12介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

U2(m,n)=c(mode/2,1)*a^2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n))+... c(mode/2,2)*a^2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n))+... c(mode/2,3)*a^2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n))+... c(mode/2,4)*a^2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n))+... c(mode/2,5)*a^2*(U1(m+5,n)+U1(m-5,n)+U1(m,n-5)+U1(m,n+5)-4*U1(m,n))+... c(mode/2,6)*a^2*(U1(m+6,n)+U1(m-6,n)+U1(m,n-6)+U1(m,n+6)-4*U1(m,n))+... 2*U1(m,n)-U0(m,n); end end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%if mode==14

for m=mode/2+1:cellnum(2)-mode/2 % 14介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if mode==16

for m=mode/2+1:cellnum(2)-mode/2 % 16介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx;

2*U1(m,n)-U0(m,n); end end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fori=mode/2:-1:1 % 透明边界条件

j=i-1;

U2(:,cellnum(2)-j)=U1(:,cellnum(2)-j)+U1(:,cellnum(2)-j-1)-U0(:,cellnum(2)-j-1)-... (V_model(:,cellnum(2)-j)*dt/dx).*((U1(:,cellnum(2)-j)-U1(:,cellnum(2)-j-1))-... (U0(:,cellnum(2)-j-1)-U0(:,cellnum(2)-j-2)));%右边界条件

U2(:,1+j)=U1(:,1+j)+U1(:,2+j)-U0(:,2+j)-(V_model(:,1+j)*dt/dx).*((U1(:,1+j)-... U1(:,2+j))-(U0(:,2+j)-U0(:,3+j)));%左边界条件

U2(cellnum(1)-j,:)=U1(cellnum(1)-j,:)+U1(cellnum(1)-j-1,:)-U0(cellnum(1)-j-1,:)-... (V_model(cellnum(1)-j,:)*dt/dz).*((U1(cellnum(1)-j,:)-U1(cellnum(1)-j-1,:))-... (U0(cellnum(1)-j-1,:)-U0(cellnum(1)-j-2,:)));%下边界条件

U2(1+j,:)=U1(1+j,:)+U1(2+j,:)-U0(2+j,:)-(V_model(1+j,:)*dt/dz).*((U1(1+j,:)-... U1(2+j,:))-(U0(2+j,:)-U0(3+j,:)));%上边界条件 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% U0=U1; U1=U2; end

colormap(gray) imagesc(U2)

共8页:

有限差分纵波波场模拟(7).doc 将本文的Word文档下载到电脑下载失败或者文档不完整，请联系客服人员解决！

下载这篇word文档