matlab.txt

子波求取

%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);

合成地震记录初始速度的确定方法
1，VSP速度
2，经验公式
3，checkshot

雷克子波 源代码
function r=leikewavelet()
dt=10000;
r=zeros(1000,1);
for f=20:80
for t=-499:500
a=(1-2*(pi*f*t/dt).^2).*exp(-(pi*f*t/dt).^2);
r(t+500)=r(t+500)+a;
end
end