RTM.rar_RTM matlab_Reverse time _migration_seismic migration

% Pavel Sinitsyn % I.D. 321254898 % migration: ex2- migration and demigration clear all tcur=0.8; %%%%%%%%%%model parameters%%%%%%%%%%%%%%%%%%%%% nt=1001;%number of sampels dt=0.004;%sample interval nx=256;%number of traces dx=10;%trace interval dz=8; nz=501; pml_size=300;%size of absorption aria in terms of elements %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%read binary files %%%%%%%%%%%%%%% fid1=fopen('STK_2SYNC'); fid2=fopen('VEL_2SYNC'); vel=fread(fid2,[nz nx],'float'); seis=fread(fid1,[nt nx],'float'); disz=0:(nz-1).*dz; disx=0:(nx-1).*dx; tim=0:(nt-1).*dt; % imagesc(disz,disx,vel); % figure % imagesc(disx,tim,seis); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% FD parameters %%%%%%%%%%%%%%%%%%%%%%%%% nx2=round(nx*dx/dz); totx=2*pml_size+4+nx2;%total number of elements at x direction totz=2*pml_size+4+nz;%total number of elements at z direction cmax=max(max(vel));%naximal velocity dt2=dt/round(dt/(dz/(2*cmax)));%dt for time FD propogation nt2=round(nt*dt/dt2); ddt=dt/dt2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%pml%%%%%%%%%%%%%%%%%%%%%% alfa=0.003; gammamax=0.5*alfa; gamma=zeros(totz,totx); lx=pml_size*dz; for ii=1:totz for jj=1:totx if (ii<=pml_size) gamma(ii+2,jj)=(1+cos(pi*(ii-1)*dz/lx))*gammamax/2; end if (jj<=pml_size) ezer=(1+cos(pi*(jj-1)*dz/lx))*gammamax/2; if (ezer>=gamma(ii,jj+2)) gamma(ii,jj+2)=ezer; end end if ((totz-ii+1)<=pml_size) ezer=(gammamax/2+cos(pi*(-totz+(ii))*dz/lx)*gammamax/2); if (ezer>=gamma(ii-2,jj)) gamma(ii-2,jj)=ezer; end end if ((totx-jj+1)<=pml_size) ezer=(gammamax/2+cos(pi*(-totx+(jj))*dz/lx)*gammamax/2); if (ezer>=gamma(ii,jj-2)) gamma(ii,jj-2)=ezer; end end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 %%%%% new velocities matrix %%%%%%%%%%%%%%%%% c=zeros(totz,totx); for ii=(pml_size+3):(totz-pml_size-2) for jj=(pml_size+3):(totx-pml_size-2) c(ii,jj)=vel(ii-pml_size-2,round((jj-pml_size-3)*dz/dx)+1); end end for ii=1:(pml_size+2) for jj=(pml_size+3):(totx-pml_size-2) c(ii,jj)=c(pml_size+3,jj); c(totz-ii+1,jj)=c(pml_size+2+nz,jj); end end for ii=1:totz for jj=1:(pml_size+2) c(ii,jj)=c(ii,pml_size+3); c(ii,totx-jj+1)=c(ii,pml_size+2+nx2); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%new stack matrix (interpolated)%%%%%%%%% stk=zeros(nt,nx2); for ii=1:nt for jj=1:nx2 xnew=(jj-1)*dz; %%%%%%%% sinc interpolation val=0; for kk=1:nx val=val+seis(ii,kk)*(sinc(xnew/dx-(kk-1))); end stk(ii,jj)=val; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f1=zeros(totz,totx);%future time step field f2=zeros(totz,totx);%current time step field f3=zeros(totz,totx);%past time step field %%%%%%%%%% finite difference wave propagation %%%%%%%%% %%%Implicit in space (4rth order) and explicit in time (2nd order)%%%%% %%%%% starting the system at the beginning %%%%%%%%%%%%%% for ii=1:nx2 f3(pml_size+3,ii+pml_size+2)=stk(nt,ii); end for ii=3:(totz-2) for jj=3:(totx-2) f2(ii,jj)=(c(ii,jj)/2)^2*dt2^2/(24*dz^2)*(16*f3(ii+1,jj)+16*f3(ii-1,jj)+16*f3(ii,jj+1)+16*f3(ii,jj-1)-f3(ii+2,jj)-f3(ii-2,jj)-f3(ii,jj+2)-f3(ii,jj-2)-60*f3(ii,jj))+f3(ii,jj); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%next time step calculation %%%%%%%%%%%%%%%%% count1=1; count2=1; figure; imagesc(f2((pml_size+3):(totz-pml_size-2),(pml_size+3):(totx-pml_size-2))); axis tight caxis([-1.6096e-004 1.8627e-004]); set(gca, 'nextplot','replacechildren', 'Visible','off'); mov(1:nt) = struct('cdata',[], 'colormap',[]);%open movie file for t=2:(nt2-3) %t*dt2 count2=count2+1; for ii=3:(totz-2) for jj=3:(totx-2) f1(ii,jj)=(c(ii,jj)/2)^2*dt2^2/(12*dz^2)*(16*f2(ii+1,jj)+16*f2(ii-1,jj)+16*f2(ii,jj+1)+16*f2(ii,jj-1)-f2(ii+2,jj)-f2(ii-2,jj)-f2(ii,jj+2)-f2(ii,jj-2)-60*f2(ii,jj))+2*f2(ii,jj)-f3(ii,jj); end end f3=f2; f2=f1.*(1-gamma); if count2==ddt for ii=1:nx2 f2(pml_size+3,ii+pml_size+2)=f2(pml_size+3,ii+pml_size+2)+stk(nt-count1,ii); end mov(count1) = getframe(gcf);%make a movie imagesc(f1((pml_size+3):(totz-pml_size-2),(pml_size+3):(totx-pml_size-2))); caxis([-1.6096e-004 1.8627e-004]); count1=count1+1 count2=0; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% movie2avi(mov, 'RTM.avi', 'compression', 'None');%make avi movie close(gcf); fclose(fid1); fclose(fid2);