%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % chap3b.m performs forward modeling (via 2-2 staggered grid FD 1st-order % acoustic eqns) to generate line-src response of velocity model for a line % src at A and receivers along rectangle boundary. Then it % backpropagates these data (via 2-2 FD algorithm for 2nd-order acoustic % wave eqn) back to receivers at A and B using the reciprocity thm of % correlation type. % Theory: The reciprocity equation of correlation type for positive t>0 values is given by % % G(B,t|A,0) = sum_x [dG(A,t|x,0)/dn @ G(x,-t|B,0) - G(A,t|x,0) @ dG(x,-t|B,0)/dn]dx % % where @ denotes temporal convolution. The code is structured so that % 1). dG(A,t|x,0)/dn and G(A,t|x,0) are computed by a FD code % for a source at A and receivers at x along a square boundary. % 2). The dG(A,t|x,0)/dn and G(A,t|x,0) values are time reversed to give dG(A,-t|x,0)/dn % and G(A,-t|x,0). These time reversed values will be % used by another FD code as the temporal history of sources on the square boundary, to % generate the pressure field at B. This new pressure field at B according to the above % formula should give us G(B,t|A,0). % % (nt,he,wi,c,fr) - input- (# time steps, height model, width model, velocity, peak frequency) % (dx,dt,nx,nz) - input- (gridpoint interval, time sample interval, # model grid pts in x and z directions) % (nxsrc,nzsrc) - input- (x,z) coordinates of source at A % (nxsrc1,nzsrc1) - input- (x,z) coordinates of source at B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; % % Input modeling parameters % [nt,he,wi,c,fr,vmax,vmin,invrho,kappa,dx,dt,nx,nz,u,w,px,pz,p,nxsrc,nzsrc,seismo,alpha]=indata1; nxsrc1=nxsrc;nzsrc1=nzsrc+100; % Locations of A and B in the model interior indata2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % FORWARD MODELING STEP TO GET G and dG/dn along square boundary % % 2-2 FD of the 1st-order acoustic equations of motion to generate data at % the recatngular boundary for sources at A=(nxsrc,nzsrc) and at % B=(nxsrc1,nzsrc1) (staggered grid code by Mark Fan of China). Data % created are G(A|x), dG(A|x)/dn, G(B|x), dG(B|x)/dn %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for is=1:2; % Loop over sources at A and B %for is=1:1; % Loop over sources at A and B xs=nxsrc;zs=nzsrc;xr=nxsrc1;zr=nzsrc1; if is==2; xs=nxsrc1;zs=nzsrc1;xr=nxsrc;zr=nzsrc;u=u.*0;w=u;px=u;pz=u;p=u; end; for it=2:nt % Loop over time steps for forward modeling u(1:nx-1,:) = (u(1:nx-1,:).* dampx(1:nx-1,:) - kappa1.*diff(p,1,1))./dampxa(1:nx-1,:); w(:,1:nz-1) = (w(:,1:nz-1).* dampz(:,1:nz-1) - kappa2.*diff(p,1,2))./dampza(:,1:nz-1); px(2:nx,:) = (px(2:nx,:).* dampx(2:nx,:) - invrho1.*diff(u(:,:),1,1) )./dampxa(2:nx,:); pz(:,2:nz) = (pz(:,2:nz).* dampz(:,2:nz) - invrho2.*diff(w(:,:),1,2) )./dampza(:,2:nz); p = px + pz; %Adding source p(xs,zs)=p(xs,zs)+rick(it); store1;plot6; % Store field values dG(A|x)/dn and G(A|x) along square boundary at x pAB(it,is)=p(xr,zr); end plot5; end p0=zeros(nx,nz);p1=p0;p2=p0;bc=zeros(nx,1); cns=4*(c*dt/dx).^2; UU=diff(UU);WW=diff(WW);seismo1=seismo(nt,2); %--------------------------------------------------------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % BACKWARD MODELING STEP TO BACKPROPAGATE FIELD FROM BOUNDARY TO A and B % % 2-2 FD of the 2nd-order acoustic wave equation to backpropagate data % from the rectangular boundary to receivers at A=(nxsrc,nzsrc) and at % B=(nxsrc1,nzsrc1). Reciprocity Eqn is % G(B|A)=int_B [dG(A,t|x)/dn@G(x,-t|B) - G(x,t|A)@dG(B,-t|x)/dn]dx %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for is=1:2 x=nxsrc1;z=nzsrc1; if(is==2) x=nxsrc;z=nzsrc;p1=p1.*0.0;p2=p1.*0.0;p0=p1.*0.0; end for it=nt-2:-1:2 % Loop over reverse time for backpropagation backward into medium for ns=1:nsp;p1(xsr(ns),zsr(ns))=UU(it,ns,is)*nxx(ns)+WW(it,ns,is)*nzz(ns);end; % dG/dn at t for ns=1:nsp;p0(xsr(ns),zsr(ns))=UU(it+1,ns,is)*nxx(ns)+WW(it+1,ns,is)*nzz(ns);end;% dG/dn at t-1 p2 = 2*p1 - p0 + cns.*del2(p1); if it==round(it/30)*30;imagesc([1:nx]*dx,[1:nz]*dx,p2');title(num2str(it));pause(0.1);end; p0=p1;p1=p2; %xs=nxsrc;zs=nzsrc;xr=nxsrc1;zr=nzsrc1; seismo1(it,is) = p2(x,z); % seismo1(it,2) = p2(x1,z1); end; end %--------------------------------------------------------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PLOT RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [nt1]=length(pAB(:,1)); subplot(111);a=diff(pAB(:,1));amx=max(abs(a));plot([0:nt1-1]*dt,[0; a(1:nt1-1)/amx],'-'); %hold on;nt1=length(seismo1(:,2));aa2=seismo1(:,2);mx=max(abs(aa2(:))); hold on;nt1=length(seismo1(:,2));aa1=seismo1(:,1);aa2=seismo1(:,2);aa12=aa1+aa2;mx=max(abs(aa12(:))); plot([0:nt1-1]*dt,aa12/mx,'--');hold off xlabel('Time (s)');ylabel('Normalized Amplitude');title('G(B,t|A,0)'); text(.05,.6,'Solid = Analytic');text(.05,.7,'Dashed = Interferometric'); %---------------------------------------------------------------------------