rmssimp.m

这是matlab在地球物理数据处理方面的源码

function [traceoff,timebase,csg] = rmssimp(h,v)

% Program to calculate a set of seismograms using simple short offset approximation
% through a stack of 
% n flat-lying layers defined by their velocity and thickness 
% The program also plots the expected traces and supeimposes on this the 
% traveltimes expected for the model using the short offset approximation.
%
% To run    [traceoff,timebase,csg] = raysimp(h,v);
%
% h = thickness of n layers in metres
% v = corresponding velocity of the layers in m/s.
% csg = resulting output set of seismograms observed, the 'pulse' is a hanning window
% traceoff is a vector containing the corresponding trace offsets.
% timebase is a vector which is the corresponding traces sample times.
% Current implementation set to 2 ms sampling, 1 second record, offsets @ 10 m spacing
% from 0 to 1 km.

numlayers = length(h);

maxoffset = 200; delx = 5; delt = 0.002; maxtime = 1; 
traceoff = 0:delx:maxoffset;
timebase = 0:delt:maxtime;
csg = zeros(length(timebase),length(traceoff));  % Common Shot Gather.
[numsamps,numtraces] = size(csg);

twoway = 2*h./v;
vrms = sqrt(cumsum(v.^2.*twoway)./cumsum(twoway))
[X,Vrms] = meshgrid(traceoff,vrms); 
[X,Twoway] = meshgrid(traceoff,cumsum(twoway)); 
times = sqrt(Twoway.^2 + X.^2./Vrms.^2);
arrtimeind = round(times/delt);
save junk arrtimeind times

for i = 1:numlayers  
   for k = 1:numtraces
      if arrtimeind(i,k) < numsamps 
         csg(arrtimeind(i,k),k) = csg(arrtimeind(i,k),k)+1;
      end
   end
end

waveletsamps = 21;
for i = 1:numtraces
   temp(:,i) = conv(hanning(waveletsamps),csg(:,i));
end



csg = temp(round(waveletsamps/2):numsamps+round(waveletsamps/2)-1,:);
imagesc(traceoff,timebase,csg), xlabel('Offset metres'), ylabel('Time seconds')