forward.m

井中地震正演程序

function forward()
xx=100;%%检波器排列长度
zz=100;%%向下延托深度
dx=1;%%道间距
dz=1;%%深度采样间隔
dt=0.0005;%%时间采样间隔
t=0.035;%%总采样时间
nx=xx/dx;
nz=zz/dz;
nt=t/dt;%%计算总采样时间
a=0.1;%%定义脉冲宽度
x0=70;%%震源位置
z0=0;%%震源位置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%构建速度模型
vmode(1:nx,1:nz)=600;
vmode(80,60)=1200;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pf(1:nx,1:nz)=0;
ps(1:nx,1:nz)=0;
pt(1:nx,1:nz)=0;
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始时刻0的值
    for ix=1:nx
        for iz=1:nz
            pf(ix,iz)=exp(-a*((ix-x0/dx)^2+(iz-z0/dz)^2));
        end
    end
% pf(nx/2,1)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始时刻1的值
    for ix=2:nx-1
        for iz=2:nz-1
            alfa=dt*vmode(ix,iz)/dx;
            beita=dt*vmode(ix,iz)/dz;
            ps(ix,iz)=(1-alfa^2-beita^2)*pf(ix,iz)+0.5*(alfa^2*(pf(ix+1,iz)+pf(ix-1,iz))+beita^2*(pf(ix,iz+1)+pf(ix,iz-1)));
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%开始延托第二个时刻的值
for it=1:nt
    for ix=2:nx-1
        for iz=2:nz-1
            alfa=dt*vmode(ix,iz)/dx;
            beita=dt*vmode(ix,iz)/dz;
            zyh=2*(1-alfa^2-beita^2)*ps(ix,iz)-pf(ix,iz)+alfa^2*(pf(ix+1,iz)+pf(ix-1,iz))+beita^2*(pf(ix,iz+1)+pf(ix,iz-1));
            pt(ix,iz)=zyh;
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下是确定右边边界条件
% for iz=1:nz
%     pt(nx,iz)=ps(nx,iz)+ps(nx-1,iz)-pf(nx-1,iz)+alfa*(ps(nx-1,iz)-ps(nx,iz)-pf(nx-2,iz)+pf(nx-1,iz));
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上是确定右边边界条件
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下是确定左边边界条件
% for iz=1:nz
%     pt(1,iz)=ps(1,iz)+ps(2,iz)-pf(2,iz)+alfa*(ps(2,iz)-ps(1,iz)-pf(3,iz)+pf(2,iz));
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上是确定左边边界条件
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下是确定下边边界条件
% for ix=1:nx
%     pt(ix,nz)=ps(ix,nz)+ps(ix,nz-1)-pf(ix,nz-1)+beita*(ps(ix,nz-1)-ps(ix,nz)-pf(ix,nz-2)+pf(ix,nz-1));
% end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上是确定下边边界条件
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%临时输出数据
    if(it==nt)
        fid=fopen('d:\VSP.grd','w');
        fprintf(fid,'%s\r\n','DSAA');
        fprintf(fid,'%d  %d\r\n',nx,nz);
        fprintf(fid,'%d  %d\r\n',0,nx-1);
        fprintf(fid,'%d  %d\r\n',0,nz-1);
        fprintf(fid,'%f  %f\r\n',min(min(pt)),max(max(pt)));
        for iz=1:nz
            fprintf(fid,'%f  ',pt(1:nx,iz));
            fprintf(fid,'\r\n');
        end
        fclose(fid);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    pf=ps;
    ps=pt;
end