波动方程逆时偏移算法
分享下近期学习的逆时偏移相关算法，以下内容来源于网页地址，网站所提供的代码很全面（可直接下载），可以修改代码用于地震勘探数据处理。
以下是计算结果。
part1
vel=[repmat(1000,[1,30]), repmat(1200,[1,30]), repmat(1500,[1,21])];
vel=repmat(vel',[1 201]);
[nz,nx]=size(vel); dx=5; x = (0:nx-1)*dx; z = (0:nz-1)*dx;
sx = (nx-1)/2*dx; sz = 0;
gx=(0:2:(nx-1))*dx; gz=zeros(size(gx));
nbc=40; nt=2001; dt=0.0005;isFS=false;
freq=25; s=ricker(freq,dt);
figure(1);set(gcf,'position',[0 0 800 400]);subplot(221);imagesc(x,z,vel);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
figure(1);subplot(222);plot((0:numel(s)-1)*dt,s);
xlabel('Time (s)'); ylabel('Amplitude');title('wavelet');
tic;seis=a2d_mod_abc28_snapshot(vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
part2
clc
close all
clear
vel=[repmat(1000,[1,30]), repmat(1200,[1,30]), repmat(1500,[1,21])];
vel=repmat(vel',[1 201]);[nz,nx]=size(vel); dx=5; x = (0:nx-1)*dx; z = (0:nz-1)*dx;
sx = (nx-1)/2*dx; sz = 0;
gx=(0:2:(nx-1))*dx; gz=zeros(size(gx)); ng=numel(gx); g=1:ng;
nbc=40; nt=2001; dt=0.0005; t=(0:nt-1)*dt;isFS=false;freq=25; s=ricker(freq,dt);
tic;seis=a2d_mod_abc28(vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
[vel_ss,refl_ss]=vel_smooth(vel,3,3,1);
figure(3);set(gcf,'position',[0 0 1000 400]);subplot(231);imagesc(x,z,vel);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
figure(3); subplot(232);imagesc(g,t,seis);colormap(gray);
title('Seismic Profile');ylabel('Time (s)');xlabel('g #');caxis([-0.25 0.25]);
tic;[img,illum]=a2d_rtm_abc28_snapshot(seis,vel_ss,nbc,dx,nt,dt,s,sx,sz,gx,gz);toc;
figure(4); set(gcf,'position',[0 0 400 600]);colormap(gray);
subplot(311);imagesc(x,z,img);caxis([-10 10]);
xlabel('X (m)'); ylabel('Z (m)'); title('rtm image');
figure(4);subplot(312);imagesc(x,z,illum);caxis([-100 100]);
xlabel('X (m)'); ylabel('Z (m)'); title('illumination compensation');
figure(4);subplot(313);imagesc(x,z,img./illum);caxis([-1 1]);
xlabel('X (m)'); ylabel('Z (m)'); title('rtm image after compensation');
vel_homo=zeros(size(vel))+min(vel(:));
tic;seis_homo=a2d_mod_abc28(vel_homo,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
figure(3);set(gcf,'position',[0 0 1000 400]);subplot(231);imagesc(x,z,vel);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
seis=seis-seis_homo;
figure(3); subplot(232);imagesc(g,t,seis);figure_title='Seismic Profile';
title(figure_title);ylabel('Time (s)');xlabel('g #');caxis([-0.25 0.25]);
[img,illum]=a2d_rtm_abc28_snapshot(seis,vel_ss,nbc,dx,nt,dt,s,sx,sz,gx,gz);
figure(5); set(gcf,'position',[0 0 400 600]);colormap(gray);
subplot(311);imagesc(x,z,img);caxis([-10 10]);xlabel('X (m)'); ylabel('Z (m)'); title('rtm image');
figure(5);subplot(312);imagesc(x,z,illum);caxis([-100 100]);xlabel('X (m)'); ylabel('Z (m)'); title('illumination compensation');
figure(5);subplot(313);imagesc(x,z,img./illum);caxis([-1 1]);xlabel('X (m)'); ylabel('Z (m)'); title('rtm image after compensation');
part3
clc
close all
clear
vel=[repmat(1000,[1,30]), repmat(1200,[1,30]), repmat(1500,[1,21])];
vel=repmat(vel',[1 201]); [nz,nx]=size(vel); dx=5; x = (0:nx-1)*dx; z = (0:nz-1)*dx;
sx = (nx-1)/2*dx; sz = 0; gx=(0:2:(nx-1))*dx; gz=zeros(size(gx));
nbc=40; nt=2001; dt=0.0005; t=(0:nt-1)*dt;isFS=false;freq=25; s=ricker(freq,dt);
[vel_ss,refl_ss]=vel_smooth(vel,3,3,1);
figure(2);set(gcf,'position',[0 0 1000 500]);colormap(gray);
subplot(231);imagesc(x,z,vel_ss);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('smooth velocity');
figure(2);subplot(232);imagesc(x,z,refl_ss);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('reflectivity');
figure(2);subplot(233);plot((0:numel(s)-1)*dt,s);
xlabel('Time (s)'); ylabel('Amplitude');title('wavelet');
tic;seis=a2d_bmod_abc28_snapshot(vel_ss,refl_ss,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
% Forward Modeling to save BC
tic; [~,bc_top,bc_bottom,bc_left,bc_right,bc_p_nt,bc_p_nt_1]=...
a2d_mod_abc28(vel_ss,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
% RTM
img=a2d_rtm_abc28(seis,vel_ss,nbc,dx,nt,dt,s,sx,sz,gx,gz,...
bc_top,bc_bottom,bc_left,bc_right,bc_p_nt,bc_p_nt_1);
figure(3);set(gcf,'position',[0 0 600 300]);colormap(gray);
imagesc(x,z,img);caxis([-100 100]);
xlabel('X (m)'); ylabel('Z (m)'); title('RTM Image with Born Data');
tic;seis_new=a2d_bmod_abc28(vel_ss,img,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
a=sum(refl_ss(:).*img(:));
b=sum(seis(:).^2);
c=sum(img(:).^2);
d=sum(seis(:).*seis_new(:));
display(['<Lm,Lm> = ',num2str(b)]);
display(['<m,LTLm> = ',num2str(a)]);
display(['<LTLm,LTLm> = ',num2str(c)]);
display(['<Lm,LLTLm> = ',num2str(d)]);
vel(:)=min(vel(:)); refl=zeros(size(vel)); refl((nz+1)/2,(nx+1)/2)=1;
[seis,bc_top,bc_bottom,bc_left,bc_right,bc_p_nt,bc_p_nt_1]...
=a2d_bmod_abc28(vel,refl,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);
img=a2d_rtm_abc28(seis,vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,...
bc_top,bc_bottom,bc_left,bc_right,bc_p_nt,bc_p_nt_1);
figure(4);set(gcf,'position',[0 0 600 800]);colormap(gray);
subplot(311);imagesc(x,z,vel);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('Velocity');
figure(4);subplot(312);imagesc(x,z,refl);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('Reflectivity');
figure(4);subplot(313);imagesc(x,z,img);caxis([-10 10]);
xlabel('X (m)'); ylabel('Z (m)'); title('Migration Greens Function');
part4
clc
close all
clear
nz=81;nx=201;
vel=zeros(nz,nx)+1000; dx=5;
x = (0:nx-1)*dx; z = (0:nz-1)*dx;
sx = (nx-1)/4*dx; sz = 0; gx=(0:2:(nx-1))*dx; gz=zeros(size(gx)); g=1:numel(gx);
nbc=40; nt=2001; dt=0.0005; t=(0:nt-1)*dt;isFS=false;freq=25; s=ricker(freq,dt);
tic; seis=a2d_mod_abc28(vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
figure(1);set(gcf,'position',[0 0 1000 500]);
subplot(221);imagesc(x,z,vel);colorbar;xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
figure(1);subplot(222);imagesc(g,t,seis);colormap(gray);
xlabel('Time (t)'); title('seismic gather');
figure(1);subplot(223);plot(t,seis(:,76));xlabel('Time (t)'); title('seismic trace');
tic; wp=a2d_wavepath_abc28(seis(:,76),vel,nbc,dx,nt,dt,s,sx,sz,gx(76),gz(76)); toc;
figure(1);subplot(224);imagesc(x,z,wp);colormap(gray);caxis([-10 10]);
xlabel('X (m)'); ylabel('Z (m)'); title('wave path for direct wave');
dx = 5;
vel=(1000:1000/150:2000);
vel=repmat(vel',[1,401]);
% vel(40:80,160:240) = 1500;
[nz,nx]=size(vel);
x = (0:nx-1)*dx;
z = (0:nz-1)*dx;
sx = (nx-1)/8*dx;
sz = 0;
gx=(0:2:nx-1)*dx;
gz=zeros(size(gx));
g=0:numel(gx);
nbc=4; nt=4001;
dt=0.0005;
t=(0:nt-1)*dt;
isFS=false;
freq=25;
s=ricker(freq,dt);
tic; seis=a2d_mod_abc28(vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
figure(2);set(gcf,'position',[0 0 1000 500]);
subplot(221);imagesc(x,z,vel);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
figure(2);subplot(222);imagesc(g,t,seis);colormap(gray);caxis([-1 1]);
xlabel('Time (t)'); title('seismic gather');
gx=gx(176);gz=gz(176); trace=seis(:,176); trace(2901:end)=0;
figure(2);subplot(223);plot(t,trace);
xlabel('Time (t)'); title('seismic trace with direct wave muted');
tic; wp=a2d_wavepath_abc28(trace,vel,nbc,dx,nt,dt,s,sx,sz,gx,gz); toc;
figure(2);subplot(224);imagesc(x,z,wp);colormap(gray);caxis([-0.1 0.1]);
xlabel('X (m)'); ylabel('Z (m)'); title('wave path for diving wave');
vel=[repmat(1000,[1 41]),repmat(2000,[1,80])];vel=repmat(vel',[1,401]);
[nz,nx]=size(vel);x = (0:nx-1)*dx; z = (0:nz-1)*dx;
sx = (nx-1)/8*dx; sz = 0; gx=(0:2:nx-1)*dx; gz=zeros(size(gx)); g=0:numel(gx);
nbc=80; nt=4001; dt=0.0005; t=(0:nt-1)*dt;isFS=false;freq=25; s=ricker(freq,dt);
tic; seis=a2d_mod_abc28(vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
trace=seis(:,176); trace(2801:end)=0;
figure(3);set(gcf,'position',[0 0 1000 500]);
subplot(221);imagesc(x,z,vel);colorbar;xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
figure(3);subplot(222);imagesc(g,t,seis);colormap(gray);caxis([-0.025 0.025]);xlabel('Time (t)'); title('seismic data');
figure(3);subplot(223);plot(t,trace);xlabel('Time (t)'); title('refraction wave');
tic; wp=a2d_wavepath_abc28(trace,vel,nbc,dx,nt,dt,s,sx,sz,gx(176),gz(176)); toc;
figure(3);subplot(224);imagesc(x,z,wp);colormap(gray);caxis([-0.0002 0.0002]);
xlabel('X (m)'); ylabel('Z (m)'); title('wave path for refraction wave');
vel=zeros(101,201)+1000; vel(81,101)=1200;[nz,nx]=size(vel);x = (0:nx-1)*dx; z = (0:nz-1)*dx;
sx = (nx-1)/8*dx; sz = 0; gx=(0:2:nx-1)*dx; gz=zeros(size(gx)); g=0:numel(gx);
nbc=80; nt=2601; dt=0.0005; t=(0:nt-1)*dt;isFS=false;freq=20; s=ricker(freq,dt);
[vel_ss,refl_ss]=vel_smooth(vel,3,3,1);
tic; seis=a2d_mod_abc28(vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
trace=seis(:,51); trace(1:1800)=0;trace(2301:end)=0;
figure(4);set(gcf,'position',[0 0 1000 500]);
subplot(221);imagesc(x,z,vel);colorbar;xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
figure(4);subplot(222);imagesc(g,t,seis(1:2601,:));colormap(gray);caxis([-0.025 0.025]);
xlabel('Time (t)'); title('seismic data for one trace');
figure(4);subplot(223);plot(t,trace);xlabel('Time (t)'); title('diffraction wave');
tic; wp=a2d_wavepath_abc28_refl(trace,vel_ss,refl_ss,nbc,dx,nt,dt,s,sx,sz,gx(51),gz(51)); toc;
figure(4);subplot(224);imagesc(x,z,wp);colormap(gray);caxis([-0.002 0.002]);
xlabel('X (m)'); ylabel('Z (m)'); title('wave path for diffraction wave');
vel=zeros(101,201)+1000; vel(81:end,:)=1500;[nz,nx]=size(vel);x = (0:nx-1)*dx; z = (0:nz-1)*dx;
sx = (nx-1)/8*dx; sz = 0; gx=(0:2:nx-1)*dx; gz=zeros(size(gx)); g=0:numel(gx);
nbc=80; nt=2601; dt=0.0005; t=(0:nt-1)*dt;isFS=false;freq=20; s=ricker(freq,dt);
[vel_ss,refl_ss]=vel_smooth(vel,3,3,1);
tic; seis=a2d_mod_abc28(vel,nbc,dx,nt,dt,s,sx,sz,gx,gz,isFS);toc;
trace=seis(:,76); trace(1:1800)=0;trace(2301:end)=0;
figure(5);set(gcf,'position',[0 0 1000 500]);
subplot(221);imagesc(x,z,vel);colorbar;
xlabel('X (m)'); ylabel('Z (m)'); title('velocity');
figure(5);subplot(222);imagesc(g,t,seis(1:2601,:));colormap(gray);caxis([-0.1 0.1]);
xlabel('Time (t)'); title('seismic data for one trace');
figure(5);subplot(223);plot(t,trace);xlabel('Time (t)'); title('reflection wave');
tic; wp=a2d_wavepath_abc28_refl(trace,vel,refl_ss,nbc,dx,nt,dt,s,sx,sz,gx(76),gz(76)); toc;
figure(5);subplot(224);imagesc(x,z,wp);colormap(gray);caxis([-2 2]);
xlabel('X (m)'); ylabel('Z (m)'); title('wave path for reflection wave');
参考资料[1].https://csim.kaust.edu.sa/files/ErSE328.2013/LAB/WE_LABS/index.html
142
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
