合成地震记录

% 地震合成记录 % 日期： 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]; %速度和密度

n

v和

den

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

R

%

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 :','s'); if isempty(reply)

f = 40; %子波频率 dt = ; else

clear a;

a = sscanf(reply,'%f',[2 1]); f = a(1); dt = a(2); end %

% 计算各反射界面所对应的采样点数%

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); %只有在有界面的地方反射系数才

f和采样间隔dt

nR

nsample/2

有值

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('地震合成记录');