smc003.m

创建地震正演模型

%  地震合成记录
%  日期：07.07.19
% 
clc
clear


reply = input('请输入层数n(Default=5):','s');
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');
if isempty(reply)
    V = [600 1000 1500 2000 2500];
    dens = [1500 1800 2000 2500 3000];
    h = [500 700 400 300];
else
    clear a;
    a = sscanf(reply,'%f',[3 n]);
    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(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;
    dt = 0.004;
else
    clear a;
    a = sscanf(reply,'%f',[2 1]);
    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;
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: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
% 
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;
    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;
    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;
    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;
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);
plot(SS,Tss,'b');
hold on;
plot(Q,Qt,'r-');
set(gca,'XAxisLocation','top');
set(gca,'YDir','reverse');
title('地震合成记录');