seisrectime[1].m

一个地震正演程序,matlab源码,对相关专业的初学者有帮助。

function varargout =seisrectime[1](varargin)

%z denotes shot number
%x denotes time t  （x1的时间大，看作多次波）
z=-3:0.025:3;
deltx=0.1;
x=0:deltx:30;
v1=1;
v2=1.4;
v3=1.8;
m=length(x);
n=length(z);
%反射系数
%r1temp=(v2-v1)/(v2+v1);
%r2temp=((v3-v2)/(v3+v2))*(1-r1temp^2);
r1temp=0.15;
r2temp=0.3;
%时距曲线方程（双曲线）
x1=abs(sqrt(28.74*z.^2+297.29));
x2=abs(sqrt(33.13*z.^2+132.89));
%把时间分配到网格上，与deltx有关
for i=1:n
    if (x1(i)/deltx-fix(x1(i)/deltx))>=0.5
        x1(i)=fix(x1(i)/deltx)*deltx+deltx;
    else x1(i)=fix(x1(i)/deltx)*deltx;
    end
end

for i=1:n
    if (x2(i)/deltx-fix(x2(i)/deltx))>=0.5
        x2(i)=fix(x2(i)/deltx)*deltx+deltx;
    else x2(i)=fix(x2(i)/deltx)*deltx;
    end
end

%初始化、转换反射系数矩阵
r1=zeros(m,n);
for i=1:n
    r1(fix(x1(i)*10),i)=r1temp;
end
r2=zeros(m,n);
for i=1:n
    r2(fix(x2(i)*10),i)=r2temp;
end

%打开子波文件，去掉文件头
fid=fopen('wavelet.dat','r');
wl=fscanf(fid,'%f');
fclose(fid);
wl=WKEEP(wl,241,'r');

r=r1+r2;
%用反褶积合成地震记录
sr=[];
for i=1:n
    srtemp=conv2(wl,wl,r(:,i),'same');
    sr=[sr srtemp];
end



%%%%%zyh
fid=fopen('D:\temp\seisrectime1.txt','w');
fprintf(fid,'%d  %d\r\n',m,n);
for j=1:n
    fprintf(fid,'%f  ',sr(1:m,j));
    fprintf(fid,'\r\n');
end
fclose(fid);

wigb(sr,1,z,x);