合成地震记录

1、% 地震合成记录% 日期：07.07.19% clcclearreply = input(请输入层数n(Default=5):,s); %层数为nif isempty(reply) n = 5;else n = sscanf(reply,%f,1 1);endreply = 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

2、 2000 2500 3000; %速度和密度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 R(ilayer) = (z2-z1) / (z2+z1);end% % 计算各反射界面所对应

3、的时间tlength% tlength(1) = 2*h(1)/V(1);for ilayer = 2:n-1 tlength(ilayer) = tlength(ilayer-1) + 2*h(ilayer)/V(ilayer);endreply = 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); f = a(1); dt = a(2);end% % 计算各反射

4、界面所对应的采样点数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/2for ilayer = 1:n-1 RR(nR(ilayer) = R(ilayer); %只有在有界面的地方反射系数才有值end%subplot(2,2,1);stem(RR);title(反射系数序列);% % 形成一个Ricker子波wavelet%

5、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(单道记录);% % 计算多层反射界面

6、是的平均速度Vav,等效深度h0以及获得偏移距deltX% for i = 1:n-1 v(i) = V(i);endVav = 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(deltX2 + 4*h02) / Vav;del

7、tT(1) = 0;for i = 2:ntrace t(i) = sqrt(deltX*i)2 + 4*h02) / 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; endend% 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% endfor i = 1:length(SS) Tss(i) = dt*i;endfor i = 1:nsample if abs(S(i)0 bg=i; break; endendfor i = 1:ntrace Qt(i) = dt*(delay(i)+bg);end for i = 1:ntrace Q(i) = SS(delay(i)+bg