制作合成地震记录，并与原始地震记录进行对比

clear,clc
%t is time，V is velocity，d is density;
%Z is wave impedance
%con is convlution，syn is synthetic seismogram
%all is international unit
%set terminal time is 2.2 second
t=1:1:2200;
for i=1:2200
    %The TERTIARY formation,base at 4550 fts(1386m)，time is 0.782s
    for i=1:782
        V(i)=3545.738;
        d(i)=999.547;
        Z(i)=V(i)*d(i);
    end
    %The PALEOCENE formation，base at 4932 fts(1503m)，time is 0.889s
    for i=783:889
        V(i)=2190.902;
        d(i)=1441.654;
        Z(i)=V(i)*d(i);
    end
    %The CRETACEOUS formation，base at 5537 fts(1687m)，time is 1.037s
    for i=890:1037
        V(i)=2494.483;
        d(i)=1089.250;
        Z(i)=V(i)*d(i);
    end
    %The J-UNCONFORMITY formation，base at 5754 fts(1753m)，time is 1.092s
    for i=1038:1092
        V(i)=2397.862;
        d(i)=1217.397;
        Z(i)=V(i)*d(i);
    end
    %The BRENT formation，base at 6082 fts(1853m)，time is 1.183s
    for i=1093:1183
        V(i)=2190.902;
        d(i)=1441.654;
        Z(i)=V(i)*d(i);
    end
    %The DUNLIN formation，base at 6512 fts(1984m)，time is 1.293s
    for i=1184:1293
        V(i)=2397.862;
        d(i)=1217.397;
        Z(i)=V(i)*d(i);
    end
    %The STATFJORD formation，base at 7212 fts(2918m)，time is 1.487s
    for i=1293:1487
        V(i)=2190.902;
        d(i)=1441.654;
        Z(i)=V(i)*d(i);
    end
    %The TRIASSIC formation，time is 2.2s
    for i=1488:2200
        V(i)=2397.862;
        d(i)=1217.397;
        Z(i)=V(i)*d(i);
    end
end
%Calculate ricker wavelet 
fm=30;          %main frequency is 30 Hz
dt1=0.001;      %sample interval of time domain
number=217;     %shot number
t1=-number/2+1:number/2;
a=(1-2*(pi*fm*t1*dt1).^2).*exp(-(pi*fm*t1*dt1).^2);
%reflection factor
for i=1:2199
    R(i)=(Z(i+1)-Z(i))/(Z(i+1)+Z(i));
end
R(1,2200)=0;
%calculate convlution
con=conv(a,R);
for i=1:2200
    syn(i)=con(i+108);
end
%load dataset
m=load('line90.asc');
k=1;
for i=1:240
    for j=1:751
        data(i,j)=m(k);
        k=k+1;
    end
end

%Figure one: Making synthetic seismogram
figure(1)
%Subgraph one: Dependence of both velocity and time
subplot(1,4,1)
plot(V,t./1000,'r','LineWidth',3);
title('V-T');
xlabel('Velocity');
ylabel('Time');
set(gca,'ydir','reverse')
%Subgraph two: Dependence of both density and time
subplot(1,4,2)
plot(d,t./1000,'r','LineWidth',3);
title('D-T');
xlabel('Density');
ylabel('Time');
set(gca,'ydir','reverse')
%Subgraph three: Dependence of both impedance and time
subplot(1,4,3)
plot(R,t./1000,'g','LineWidth',3);
title('I-T');
xlabel('Impedance');
ylabel('Time');
set(gca,'ydir','reverse')
%Subgraph four: Dependence of both ricker's amplitude and time
subplot(1,4,4)
plot(a,t1,'b','LineWidth',3);
axis([-0.5 1.0 -100 100]);
xlabel('Amplitude');
ylabel('Time');
title('Ricker-wavelet');
set(gca,'ydir','reverse')
%Figure two: Comparision Synthetic seismogram with Origin seismic channel record
figure(2)
%Subgraph one: Synthetic seismogram
subplot(1,2,1)
plot(syn,t./1000,'k','LineWidth',3);
axis([-0.08 0.1 0.7 2.2]);
title('Synthetic seismogram');
set(gca,'ydir','reverse')
%Subgraph two: Origin seismic channel record
subplot(1,2,2)
plot(data(60,:)./1000,0.7:0.002:2.2,'k','LineWidth',3)
axis([-0.08 0.1 0.7 2.2]);
title('Origin seismic channel record');
set(gca,'ydir','reverse')