read_segy_file.m

基于Matlab的地震数据处理显示和测井数据显示于处理的小程序

function [seismic,ebcdic_header,binary_header]=read_segy_file(filename,varargin)
% Function reads a SEG-Y file and outputs a seismic structure
% Written by: E. R.: March 5, 2000
% Last updated: April 13, 2006: Add field "from" with the full filename of the 
%             dataset read
%
%             [seismic,ebcdic_header,binary_header]=read_segy_file(filename,varargin)
%
% INPUT       (all input arguments are optional)
% filename    name of the  file to read
%             if the name is omitted, empty, or a file with this name is not found
%             a file selector box will pop up to allow interactive file selection
% varargin    Variable number of arguments. Each argument is a cell array whose first
%             element is a keyword and whose other elements can be strings, numeric
%             values, or other cell arrays
%             Possible keywords are:
%     format  floating point format to use; this cell array has the form
%             {'format',fp-format} where fp_format is one of the three strings:
%             'ibm'      IBM floating point format (standard and default)
%             'ieee'     IEEE format, big endian  (Sun, SGI, etc.)
%             'header'   The binary header bytes 25-26 are used to determine the format
%             Default: {'format','header'}; 
%     headers  header values to be read from binary trace header; this cell array has
%             the form {'headers',{mnem1,first,bytes,units,dimension,description}, ...
%             {mnem2,first,bytes,units,dimension,description},...} where
%             "mnem1", "mnem2", ... denote header mnemonics (such as CDP, OFFSET), 
%             "first" denotes the first byte in the binary header, 
%             "bytes" denotes the number of bytes occupied by the mnemonic (2 or 4), 
%             "units" denotes the units of measurement for the header value, and 
%             "description" is a description of the header value. 
%             Example:
%             {'headers',{'ILINE_NO',181,4,'n/a','CDP number'}, ...
%                        {'OFFSET',37,4,'m','Source-receiver distance'}}
%             See below for a list of headers retrieved by default.
%             Default: {'headers',{}}
%     times   times to output; this cell array has the form
%             {'times',first,last}; all samples with times between (and including)
%             first and last (in ms) are output.
%             Default: {'times',[]}  i.e. all times
%     traces  select traces to output; this cell array has the form
%             {'traces',criterion}; the variable "criterion" is used to 
%             determine which traces to output. "criterion" can be an index
%             vector specifying the traces to output
%             Examples: {'traces',1:2:100}
%                       {'traces',[2,5,7:10,22]}
%             Alternatively, it can be a string with a logical expression involving 
%             trace headers such as '10 <= cdp & 100 >= cdp'
%             Examples: {'traces','cdp == 100 & offset > 100'}
%                       {'traces','14000 < cdp  & (14660 >= cdp | 14680 <= cdp)'};
%             The variables in the logical relationships must be headers of the 
%             data set; Use of functions "fix", "mod", and "round" are permitted; all
%             other function names will be interpreted as headers and likely cause an
%             error; the case of the headers in an expression does not matter.
%             Default: {'traces',[]}
%    ignoreshift  By default this function reads byte locations 109-110 (see 
%             header "lag" below) and applies the shift to the seismic data;
%             This behavior can be overwritten by setting this parameter to true;
%             Default: {'ignoreshift',logical(0)}
%
% Headers retrieved by default are (any one of these headers is removed if it turns out to 
% be identically zero):
%        ds_seqno      Trace sequence number within line (1-4)
%        ffid          Original Field record number (9-12)
%        o_trace_no    Trace sequence number within original field record (13-16)
%        source        Energy source point number (17-20)
%        cdp           CDP ensemble number (21-24)
%        seq_cdp       Trace sequence number within CDP ensemble (25-28)
%        trc_type      Trace ID (1=live,2=dead,3=dummy,4=time break,...) (29-30)
%        offset        Distance from source point to receiver group (37-40)
%        rec_elev      Receiver elevation (41-44);
%        sou_elev      Surface elevation at source (45-48)
%        depth         Source depth below surface (49-52)
%        sou_h2od      Water depth at source (61-64)
%        rec_h2od      Water depth at receiver group (65-68)
%        sou_x         X coordinate of source (73-76)
%        sou_y         Y coordinate of source (77-80)
%        rec_x         X coordinate of receiver (81-84)
%        rec_y         Y coordinate of receiver (85-88)
%        lag           Lag time between shot and recording start in ms (109-110)
%                      (the value of lag is added to the start time of the 
%                      seismic; hence it can be used to simulate non-zero start
%                      time of the data)
%                      see also parameter "ignoreshift", above.
%        cdp_x         X coordinate of CDP (181-184)
%        cdp_y         Y coordinate of CDP (185-189)
%        iline_no      In-line number (189-192)
%        xline_no      Cross-line number (193-196)
%               The numbers in parentheses at the end of the line denote the location 
%               of the corresponding bytes in the SEG-Y trace header
%    
% OUTPUT
% seismic       Seismic structure
%       seismic.type               'seismic' (type of structure)
%       seismic.traces             Array of seismic traces
%       seismic.first              Start time of seismic (in ms)
%       seismic.last               End time of seismic (in ms)
%       seismic.step               Sample interval of seismic (in ms)
%       seismic.units              Time units used (ms)
%       seismic headers            Matrix with header mnemonics (one row 
%                                  per header)
%       seismic.header_info        Three-column cell array with header info 
%                                  (one row per header)
%       seismic.line_number        Line number (5-8)
%       seismic.reel_number        Reel number (9-12)
%       seismic.cdp_fold           CDP fold
%       seismic.traces_per_record  Data traces per record (13-14)
%       seismic.aux_per_record     Auxiliary traces per record (15-16)
%       seismic.history            A four element cell array. The first element
%                                  is the start date/time of the program that 
%                                  invoked this function; the second element is
%                                  the start date/time this function was executed;
%                                  and the last cell contains the name if the file 
%                                  that was read
%            Example:
%               seismic.offset     contains the offset for each trace
%               seismic.header_info.cdp two-element cell array {'m','Offset'}
%                                  the first element represents the units of measurement, the 
%                                  second is a description of the header
%                 
% ebcdic_header         EBCDIC reel header converted to ASCII
% binary_header         Binary reel header
%              
% 	[seismic,ebcdic_header,binary_header]=read_segy_file(filename,varargin)

global S4M ABORTED

ABORTED=logical(0);
seismic=[];
ebcdic_header=[];
binary_header=[];

run_presets_if_needed

%	Set default parameters
param.format='header';
param.headers={};
param.ignoreshift=logical(0);
param.times=[];
param.traces=[];

if ~isempty(varargin)
   param=assign_input(param,varargin);   % Interpret/check input data
end

if strcmpi(param.format,'ieee')
   param.format='ieee-be';
end
o_format='ieee-be';

if nargin < 1
   filename='';
end
if nargout < 2
   parameters=initiate_segy_file_reading_no1(filename,param,o_format);
elseif nargout == 2
   [parameters,ebcdic_header]=initiate_segy_file_reading_no1(filename,param,o_format);
else
   [parameters,ebcdic_header,binary_header]=initiate_segy_file_reading_no1(filename,param,o_format);
end

if ABORTED
   return
end

param.format=parameters.format;
fid=parameters.fid;
seismic=parameters.seismic;
indices2=parameters.indices2;
indices4=parameters.indices4;
idx_time=parameters.idx_time;
filename=parameters.filename;
no_samples=parameters.no_samples;
n_default_headers=parameters.n_default_headers;
nindices2=parameters.nindices2;
nindices4=parameters.nindices4;

[header,nh_select,max_traces,idx_trace,h_index,trace_select,param]=select_traces(seismic,fid,param);

ier=0;

%       Reserve room for header and seismic
if max_traces <= 1000
  grab=max_traces;
else
  grab=1000;
  grab_inc=grab;
end
txt=['Make room for ',num2str(grab),' traces'];
ltext=length(txt);
%fprintf('\n')

headers=zeros(size(seismic.header_info,1),grab);
seismic.traces=zeros(length(idx_time),grab);
ntraces=0;
ntraces_read=0;

%       Read seismic trace header
four_bytes=zeros(60,1);
two_bytes=zeros(90,1);

test=fread(fid,1,'int32');

                try  % Catch errors due to premature end of SEG-Y file

          while ~isempty(test) & ntraces < max_traces   % read trace
ntraces_read=ntraces_read+1;
% trace_no=ntraces_read;          % For use in trace selection

if ntraces >= grab
   grab=grab+grab_inc;
   headers=[headers,zeros(size(seismic.header_info,1),grab_inc)];
   seismic.traces=[seismic.traces,zeros(length(idx_time),grab_inc)];
   if S4M.alert
      txt=['Make room for ',num2str(grab),' traces'];
      if grab == 2*grab_inc
         fprintf(txt)
      else
         fprintf([char(8*ones(1,ltext)),txt])
      end
      ltext=length(txt);
   end
end 
four_bytes(1)=test;
four_bytes(2:7)=fread(fid,6,'int32');
two_bytes(15:18)=fread(fid,4,'int16');
four_bytes(10:17)=fread(fid,8,'int32');
two_bytes(35:36)=fread(fid,2,'int16');
four_bytes(19:22)=fread(fid,4,'int32');
two_bytes(45:90)=fread(fid,46,'int16');
four_bytes(46:60)=fread(fid,15,'int32');

depth_scalar=two_bytes(35);
if depth_scalar > 0
  four_bytes(11:17)=four_bytes(11:17)*depth_scalar;
elseif depth_scalar < 0
  four_bytes(11:17)=four_bytes(11:17)/abs(depth_scalar);
else
  if ier < 2, disp('Alert from "read_segy_file": Scalar for elevations and depths (bytes 69-70) is zero'); end
  ier=ier+1;
end

coord_scalar=two_bytes(36);
if coord_scalar > 0
   four_bytes(19:22)=four_bytes(19:22)*coord_scalar;
elseif coord_scalar < 0
   four_bytes(19:22)=four_bytes(19:22)/abs(coord_scalar);
else
   if ier < 2, disp('Alert from "read_segy_file": Scalar for coordinates (bytes 71-72) is zero'); end
   ier=ier+1;
end

headers(:,ntraces+1)=[four_bytes(indices4(1:nindices4)); ...
                               two_bytes(indices2(1:nindices2)); ...
                              four_bytes(indices4(nindices4+1:end)); ...
                               two_bytes(indices2(nindices2+1:end))];

if trace_select == 0
   boolean=ntraces < max_traces;
elseif trace_select == 1 
   for jj=1:nh_select
      eval([header{jj},'=headers(h_index(jj),ntraces+1);']);
   end
                         try
   boolean=ntraces < max_traces & eval(param.traces);
                        catch
   error([' The input argument  ',param.traces, ...
         '  of read_segy_file is probably not a valid logical expression.'])
                         end
else
                        try
   boolean=ntraces < max_traces & sum(ismember(idx_trace,ntraces_read));
                       catch
   error([' The input argument  ',param.traces, ...
         '  of read_segy_file is probably not a valid logical expression'])
                        end
end

%       Read seismic trace
if strcmpi(param.format,'ibm')    % IBM format (requires conversion)
   %if ntraces < 10; disp(param.format); end %test
   temp=fread(fid,no_samples,'uint');
   if boolean
      ntraces=ntraces+1;
      seismic.traces(:,ntraces)=ibm2ieee(temp(idx_time ));
   end
else                                       % Other format (no conversion required)
   % if ntraces < 10; disp('ieee'); end%test
   temp=fread(fid,no_samples,'float32');
   if boolean
      ntraces=ntraces+1;
      try
      seismic.traces(:,ntraces)=temp(idx_time);
      catch
      % keyboard
      end
   end
end

test=fread(fid,1,'int32');                 % Check if there is one more trace

             end

                                        catch
 msgdlg('Prematude end of SEG-Y file or other file reading error.')

                                        end

fclose(fid);
fprintf('\n')

if S4M.history
   seismic=s_history(seismic,'add',filename);
end

%       Remove unused reserved space
%seismic.traces=seismic.traces(:,1:ntraces);
seismic.traces(:,ntraces+1:end)=[];
headers(:,ntraces+1:end)=[];

%       Remove headers with duplicate header mnemonics
if isfield(seismic,'header_info')
   mnems=seismic.header_info(:,1);
   lmnems=length(mnems);
   bool=logical(zeros(lmnems,1));
   if lmnems > 1
      try
      if length(unique(mnems)) < lmnems
         for ii=lmnems-1:-1:2
            bool1=ismember(mnems(1:ii),mnems(ii+1));
            bool(1:ii)= bool1(1:ii) | bool(1:ii);
         end
         seismic.header_info(bool,:)=[];
         headers(bool,:)=[];
      end
      catch
         dbstack
         ple
         alert('An error of unknown orign has been cought.')
         keyboard
      end
   end
end


%       Remove header mnemonics that contain only zeros
nh=size(headers,1);
index=1:nh;
for ii=n_default_headers:-1:1;
  idx=find(headers(ii,:) ~= 0);
  if isempty(idx)
     index(ii)=[];
  end
end

if ~isempty(index)
   seismic.header_info=seismic.header_info(index,:);
   seismic.headers=headers(index,:);

%       Check if Header "lag" still exists
   [index,ier]=header_index1(seismic,'lag');


   if ier == 0  &  ~param.ignoreshift
      disp('Seismic data shifted since header "lag" is not identically zero.')
      disp(['Lag varies from ',num2str(min(seismic.headers(index,:))),' to ',num2str(max(seismic.headers(index,:)))])
      seismic=s_shift(seismic,{'shifts',seismic.headers(index,:)});
   end
else
   seismic=rmfield(seismic,'header_info');
end

if seismic.aux_per_record == 0
   seismic=rmfield(seismic,'aux_per_record');
end

if isempty(seismic.traces)
   msgdlg([' Alert from "read_segy_file": No seismic traces read from file ',filename])
end

seismic.fp_format_of_segy_file=param.format;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [parameters,ebcdic_header,binary_header]=initiate_segy_file_reading_no1(filename,param,o_format)

global S4M ABORTED

parameters=[];

%       Open the file
if ~isempty(filename)
   fid=fopen(filename,'r',o_format);
   if fid == -1 
      if ~S4M.compiled
         fprintf(['... unable to find file ',filename,'\n']);
      end
   else
      filename2S4M(filename)
   end
else
   fid=-1;
end

if fid==-1     % Open file selector window
   [filename,ierr]=get_filename4r('sgy');
   if ierr
      return
   end
   fid=fopen(filename,'r',o_format);
   if fid < 0
      error(['File "',filename,'" could not be opened.'])
   end
end


%  Read EBCDIC reel header 1
header_ebcdic=fread(fid,3200,'uchar');