reffind.m

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

function [reflog] = reffind(structure)
% This program calculates a reflectivity time series for a simple
% geology produced by the student.  This geology is simply in
% terms of depths of interfaces and interval velocities.
% In this version of the program you must edit the geology below
% prior to running.  The geology is contained in 'structure'
% The output reflog consists of a two column matrix with the first
% column equal to the time base (presently set equal to that of our
% data of 512 ms) and the second is the corresponding reflectivity
% time series.  These can be used as inputs to the program synth
% but be sure to do help synth before running it.
% To run  [reflecitvity] = reffind(structure) where an example structure
% is:
%structure =  [  0        1500; % Geologic structure defined in terms of
%           8        2000;      % the depth of an geologic contact and the
%         30        1400;      % seismic velocity below that contact
%         36        1000;      % for example the velocity for the first layer from
%         40        1400;      % the surface to 8 meters depth is 1500 m/s
%	  300	    1800;
%        1000        2000 ] This was a simple model over the university farm but
% you could put in sonic log information if you so desired.
%
% Some typical velocities that have been observed (see Waters, Reflection Seismology, for example)
% Glacial Till - unsaturated (0.43 - 1.04 km/s), saturated (1.73 km/s)
% Wet Clay - 1.50-1.65 km/s
% Clay - 1.1 - 2.5 km/s
% Impermeable Argillacous Clay - 2.0 km/s
% Weathered Layer - 0.3-0.9 km/s
% Sand and Gravel - unsaturated (0.38-0.5 km/s), saturated (1.67 km/s)
% Cretaceous and Tertiary Sandstone (depths > 0.3 km) 2.1 to 3.5 km/s
% Coal - 1.2 km/s (this is a quess on my part).
close all
[m,n] = size(structure);
structure(:,3:4) = zeros(m,2);
structure(2:m,3) = diff(structure(:,1))./structure(1:m-1,2);% Transit Time to interface in seconds
structure(:,3)
structure(:,4) = 2*cumsum(structure(:,3)); % Two way reflection time to each interface.
structure(:,4)

reflectivity = diff(structure(:,2))./(structure(1:m-1,2)+structure(2:m,2))

l = length(reflectivity);
delt = 0.001;  % Sampling rate

reflog = [1:512]'*delt;  % log of reflectivity - first column is timebase
reflog(:,2) = zeros(512,1);


for i = 1:l
 temp = round(structure(i+1,4)/delt);  % index of the reflection coefficient in time
 if temp < length(reflog(:,1))
        reflog(temp,2) = reflectivity(i);
 end % if
end % for i

figure
subplot(2,1,1)
plot(reflog(:,1),reflog(:,2),'r')  % Note - plot is only for a rough comparison
ylabel('Reflection Coefficient'), xlabel('Time (seconds)'), title('Reflectivity - No Multiples')
axis([0 0.512 -1.1*max(abs(reflog(:,2))) 1.1*max(abs(reflog(:,2)))]);
subplot(2,1,2)
stairs(structure(:,4),structure(:,2))
junk = axis; junk(2) = 0.512; axis(junk)
ylabel('Interval Velocity m/s'), xlabel('Time (seconds)'), title('Interval Velocity')