2905伪谱法地震波正演模拟.zip

function Forward_modeling_PS() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% This is a program for calculating elastic wavefield in isotropic media using %%% %%% pseudo-spectral method. %%% Reference: D. Kosloff, M. Reshef, and D. Loewenthal, 1984. BSSA, 74(3): 875-891.%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; %%% setting parameters %%% %%%% number of grids of working area Nx=301; %% x direction Ny=201; %% y direction %% spatial interval (meter) dx=10; dy=10; af_t=50; % central frequency of Ricker wavelet %%% end of setting parameters %%% % set location of source S_loco_x=(Nx-1)/2; % x S_loco_y=(Ny-1)/2+3; % y % location of acquisition array detct_y=S_loco_y; %%%%%% quality contral (don't change!!!!!) %%%%%%%%%%% Xmax=(Nx-1)*dx; Ymax=(Ny-1)*dy; x=0:dx:Xmax; y=0:dy:Ymax; edge_num=20; % number of grids occupied by boundary layer [rou,Vp,Vs,t_max]=elastic_model(Nx,Ny); dd=max(dx,dy); %%% DT=0.4*dd/max(max(Vp)); %%% Nt=t_max/DT-mod(t_max/DT,1); %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% re_rou=rou.^-1; % give 1/rou % build the constitutive relation, unit: kg/(m*s*s) lamda=rou.*(Vp.^2-2*Vs.^2); miu=rou.*Vs.^2; M=lamda+2*miu; c11=M; c12=lamda; c16=0; c21=c12; c22=M; c26=0; c61=c16; c62=c26; c66=2*miu; %% calculate the wavenumber x1=x; for ii=1:Ny-1 x1=vertcat(x1,x); end y1=y'; for jj=1:Nx-1 y1=horzcat(y1,y'); end fx1=(-(Nx-1)/2:-1)/Nx/dx; fx2=(0:(Nx-1)/2)/Nx/dx; fx=horzcat(fx2,fx1); Kx=2*pi*fx; Kx=complex(0,Kx); fy1=(-(Ny-1)/2:-1)/Ny/dy; fy2=(0:(Ny-1)/2)/Ny/dy; fy=horzcat(fy2,fy1); Ky=2*pi*fy; Ky=complex(0,Ky); Kx1=Kx; for ii=1:Ny-1 Kx1=vertcat(Kx1,Kx); end Ky1=Ky'; for jj=1:Nx-1 Ky1=horzcat(Ky1,Ky'); end Kx=Kx1; Ky=Ky1; %% end of calculation of wavenumber %% %% build the attenuation boundary %%% for ii=1:edge_num absorb(ii,1)=exp(-(0.015*(edge_num-ii))^2); absorb_re(edge_num-ii+1,1)=exp(-(0.015*(edge_num-ii))^2); end absorb1=absorb; absorb11=absorb_re; absorb2=absorb'; absorb22=absorb_re'; for jj=1:Nx-1 absorb1=horzcat(absorb1,absorb); % up absorb11=horzcat(absorb11,absorb_re); % down end for jj=1:Ny-1 absorb2=vertcat(absorb2,absorb'); % left absorb22=vertcat(absorb22,absorb_re'); % right end temp1=ones(Ny-2*edge_num,Nx); temp1=vertcat(absorb1,temp1); absorb1=vertcat(temp1,absorb11); temp2=ones(Ny,Nx-2*edge_num); temp2=horzcat(absorb2,temp2); absorb2=horzcat(temp2,absorb22); %% end of building the attenuation boundary %%% %% the spatial distribution of source energy %% x0=(dx*S_loco_x)*ones(Ny,Nx); y0=(dx*S_loco_y)*ones(Ny,Nx); f_y=exp(-((x1-x0).^2+(y1-y0).^2)); %%%%% time_max=(Nt-1)*DT; time=0:DT:time_max; ht=(1-2*pi^2*af_t^2*((time).^2)).*exp(-pi^2*af_t^2*(time).^2); %% Ricker wavelet for ii=1:Nt if ht(ii)==min(ht) jjj=ii; end end for ii=jjj:Nt if abs(ht(ii))<1e-3 T0=DT*ii; break; end end if Nt<50 T0=0.03; end T0=T0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% initial values f_x=zeros(Ny,Nx); ux=zeros(Ny,Nx); uy=zeros(Ny,Nx); uxnsemi_t=zeros(Ny,Nx); uynsemi_t=zeros(Ny,Nx); P_xx=zeros(Ny,Nx); P_yy=P_xx; P_xy=P_xx; % ***************************************** % for ii=1:Nt uxn_semi_t=uxnsemi_t; uyn_semi_t=uynsemi_t; uxn_1=ux; uyn_1=uy; iii=ii; Nt; %% source function (Ricker wavelet) h=(1-2*pi^2*af_t^2*((ii*DT-T0)^2))*exp(-pi^2*af_t^2*(ii*DT-T0)^2); % Ricker wavelet f_y1=h*f_y; fai=rou.*f_y1; F_P_xx=fft2(P_xx+fai); F_P_yy=fft2(P_yy+fai); F_P_xy=fft2(P_xy); F_U=Kx.*F_P_xx+Ky.*F_P_xy; F_V=Kx.*F_P_xy+Ky.*F_P_yy; P=re_rou.*ifft2(F_U); Q=re_rou.*ifft2(F_V); uxnsemi_t=uxn_semi_t+DT*P; ux=uxn_1+DT*uxnsemi_t; uynsemi_t=uyn_semi_t+DT*Q; uy=uyn_1+DT*uynsemi_t; ux=real(ux); uy=real(uy); %% apply the attenuation boundary ux=ux.*absorb1; ux=ux.*absorb2; uy=uy.*absorb1; uy=uy.*absorb2; uxnsemi_t=uxnsemi_t.*absorb1; uxnsemi_t=uxnsemi_t.*absorb2; uynsemi_t=uynsemi_t.*absorb1; uynsemi_t=uynsemi_t.*absorb2; F_P=fft2(ux); F_Q=fft2(uy); F_exx=Kx.*F_P; F_eyy=Ky.*F_Q; F_exy=(Kx.*F_Q+Ky.*F_P)/2; exx=ifft2(F_exx); eyy=ifft2(F_eyy); exy=ifft2(F_exy); exx=real(exx); eyy=real(eyy); exy=real(exy); P_xx=c11.*exx+c12.*eyy+c16.*exy; P_yy=c21.*exx+c22.*eyy+c26.*exy; P_xy=c61.*exx+c62.*eyy+c66.*exy; %% record the seismogram %% ux_graph(ii,:)=ux(detct_y,:); uy_graph(ii,:)=uy(detct_y,:); end ux123=ux(edge_num+1:Ny-edge_num,edge_num+1:Nx-edge_num); uy123=uy(edge_num+1:Ny-edge_num,edge_num+1:Nx-edge_num); ux_graph123=ux_graph(:,edge_num+1:Nx-edge_num); uy_graph123=uy_graph(:,edge_num+1:Nx-edge_num); %% save snapshots and seismograms save ux ux123; save uy uy123; save ux_graph ux_graph123 Nx Ny dx dy DT; save uy_graph uy_graph123 Nx Ny dx dy DT; %% plot the vertical component snapshot figure(1); imagesc(x,y,real(uy)); figure(gcf); axis([0 Xmax 0 Ymax]); axis equal tight; %% plot the horizontal component snapshot figure(2); imagesc(x,y,real(ux)); figure(gcf); axis([0 Xmax 0 Ymax]); axis equal tight;