R(ilayer) = (z2-z1) / (z2+z1);
end
%
%计算各反射界面所对应的时间tlength
%
tlength(1) = 2*h(1)/V(1);
for ilayer = 2:n-1
tlength(ilayer) = tlength(ilayer-1) + 2*h(ilayer)/V(ilayer);
end
reply = input('请输入Ricker子波的频率f和采样间隔dt(Defalt=40 0.004):','s');
if isempty(reply)
f = 40;%子波频率f和采样间隔dt
dt = 0.004;
else
clear a;
a = sscanf(reply,'%f',[2 1]);
deltT(i) = t(i) - t(1);
delay(i) = floor(deltT(i)/dt);
end
%
%形成地震记录SS不明白！！！！
%
for i = 1:ntrace
SS(1:2*nsample,i) = i*1;
for j = delay(i)+1:2*nsample
SS(j,i) = S(j-delay(i))+1*i;
plot(SS,Tss,'b');
hold on;
plot(Q,Qt,'r-');
set(gca,'XAxisLocation','top');
set(gca,'YDir','reverse');
title('地震合成记录');
for i = 1:length(wavelet);
Tw(i) = i*dt;
end
%subplot(2,2,2);
plot(Tw,wavelet);%应该是能画出雷克子波
title('Ricker子波');
%
%褶积后得到单道记录S
%
S = conv(RR,wavelet);
for i = 1:length(S)
Ts(i) = i*dt;
end
%subplot(2,2,3);
plot(S,Ts);
set(gca,'XAxisLocation','top');
set(gca,'YDir','reverse');
title('单道记录');
%
%计算多层反射界面是的平均速度Vav,等效深度h0以及获得偏移距deltX
f = a(1);
dt = a(2);
end
%
%计算各反射界面所对应的采样点数nR
%
nsample = floor(tlength(n-1)/dt);
for ilayer = 1:n-1
nR(ilayer) = floor(tlength(ilayer)/dt);
end
%
%形成反射系数序列RR
%
RR(1:2*nsample) = 0;%？这个地方反射系数的长度应该是nsample/2
a = sscanf(reply,'%f',[1 1]);
deltX = a;
end
%
%计算延迟时间delay完全不知道算的啥意思！！！！！@
%
ntrace = 10;
t(1) = sqrt(deltX^2 + 4*h0^2) / Vav;
deltT(1) = 0;
for i = 2:ntrace
t(i) = sqrt((deltX*i)^2 + 4*h0^2) / Vav;
for ilayer = 1:n-1
RR(nR(ilayer)) = R(ilayer); %只有在有界面的地方反射系数才有值
end
%subplot(2,2,1);
stem(RR);
title('反射系数序列');
%
%形成一个Ricker子波wavelet
%
wavelet = ricker(f,dt);
%
for i = 1:n-1
v(i) = V(i);
end
Vav = sum(h) / sum(h/v);
h0 = sum(h);
reply = input('请输入偏移距deltX(Defalt = 400):','s');
if isempty(reply)
deltX = 400;
else
clear a;
V = a(1,:);
dens = a(2,:);
h = a(3,:);
end
%
%计算反射系数R
%
for ilayer = 1:n-1
z1(ilayer) = V(ilayer) * dens(ilayer);
z2(ilayer) = V(ilayer+1) * dens(ilayer+1);%各层反射系数R
%地震合成记录
%日期：07.07.19
%
clc
clear
reply = input('请输入层数n(Default=5):','s');%层数为n
if isempty(reply)
n = 5;
else
n = sscanf(reply,'%f',[1 1]);
end
reply = input...
('请输入各层速度、密度及层厚(Defaul=[600 1000 1500 2000 2500;1500 1800 2000 2500 3000;500 700 400 300]):','s');
end
end
% for i = 1:ntrace
% S(1:nsample,i) = i*1;
% delay = 10*(i-1);
% for j = delay+1:nsample
% S(j,i) = SS(j-delay)+1*i;
% end
% end
for i = 1:length(SS)
Tss(i) = dt*i;
if isemp00];
dens = [1500 1800 2000 2500 3000];%速度和密度v和den
h = [500 700 400 300];
else
clear a;
a = sscanf(reply,'%f',[3 n]);
end
for i = 1:nsample
if abs(S(i))>0
bg=i;
break;
end
end
for i = 1:ntrace
Qt(i) = dt*(delay(i)+bg);
end
for i = 1:ntrace
Q(i) = SS(delay(i)+bg,i);
end
%subplot(2,2,4);