read_segy_file.m

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

if isempty(header_ebcdic)
   errordlg('EBCDIC header is empty; requested file is either empty or not an SEG-Y file.', ...
      S4M.name)
   ABORTED=logical(1);
end

if nargout > 1
   ebcdic_header=char(ebcdic2ascii(reshape(header_ebcdic,80,40)'));
end

% Read binary real header
bh=fread(fid,400,'uchar');

two_bytes=bh(1:2:399)*256+bh(2:2:400);
four_bytes=((bh(1:4:9)*256+bh(2:4:10))*256+bh(3:4:11))*256+bh(4:4:12);

if nargout > 2
   binary_header=two_bytes;
end

seismic.type='seismic';
seismic.tag='unspecified';
[dummy,name]=fileparts(S4M.filename);
seismic.name=name;
seismic.from=fullfile(S4M.pathname,S4M.filename);
seismic.line_number=four_bytes(2);
seismic.reel_number=four_bytes(3);
seismic.traces_per_record=two_bytes(7);
seismic.aux_per_record=two_bytes(8);
seismic.first=0;
seismic.step=two_bytes(9)/1000;
no_samples=two_bytes(11);
seismic.last=seismic.step*(no_samples-1);
seismic.units='ms';

if two_bytes(13) == 1 
   if strcmpi(param.format,'header')  | strcmpi(param.format,'ibm')
      param.format='ibm';
 
   else
      display('IEEE format requested')
      warning(['Data apparently stored as 32-bit IBM floating point numbers; ' ...
         'data sample format code = ',num2str(two_bytes(13))]);
      
   end

elseif two_bytes(13) == 5
   if strcmpi(param.format,'header')  | strcmpi(param.format,'ieee-be')
      param.format='ieee-be';
 
   else
      display('IBM format requested')
      warning(['Data apparently stored as 32-bit IEEE big-endian floating point numbers; ' ...
         'data sample format code = ',num2str(two_bytes(13))]);
      
   end

else

   myerror(['Data in unsupported format; ' ...
         'data sample format code = ',num2str(two_bytes(13))]);
   if ABORTED
      return
   end
end


if two_bytes(14) ~= 0
  seismic.cdp_fold=two_bytes(14);
end

if two_bytes(28) == 1
  units='m';
elseif two_bytes(28) == 2
  units='ft';
else
  units='unknown';
end

[seismic,idx_time]=select_times(seismic,no_samples,fid,param);

%       Compute indices of two-byte and four-byte header words which are to be saved
%       by default
indices2=([29;109]-1)/2+1;
indices4=([1;9;13;17;21;25;37;49;61;65;73;77;45;81;85;41;181;185;189;193]-1)/4+1;
nindices2=length(indices2);
nindices4=length(indices4);
i_lag=nindices4+2;       % Index of lag
[seismic,indices2,indices4,n_default_headers] = ...
             select_headers(seismic,units,indices2,indices4,i_lag,param);


parameters.fid=fid;                     % File ID
parameters.seismic=seismic;             % Seismic structure without traces and headers
parameters.indices2=indices2;           % Indices of 2-byte headers
parameters.indices4=indices4;           % Indices of 4-byte headers
parameters.idx_time=idx_time;           % indices of trace samples to output
parameters.format=param.format;         % Data format to read
parameters.no_samples=no_samples;       % Number of samples per trace input
parameters.filename=filename;           % Name of SEG-Y file
parameters.traces=param.traces;         % Parameters controlling trace selection
parameters.n_default_headers=n_default_headers;   % Number of headers read by default
parameters.i_lag=i_lag;                 % Description index of lag
parameters.nindices2=nindices2;
parameters.nindices4=nindices4;        


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [header,nh,max_trace,idx_trace,h_index,trace_select,param] = ...
   select_traces(seismic,fid,param)
% Function creates parameters which  control which traces are output

if isempty(param.traces)
  max_trace=inf;
  idx_trace=[];
  h_index=[];
  trace_select=0;
  header=[];
  nh=[];

elseif ischar(param.traces)  % Compute criteria for outputting/discarding traces
  trace_select=1;
  max_trace=inf;
  idx_trace=[];
  param.traces=lower(param.traces);
  header=unique(extract_words(param.traces));

  nh=length(header);
  if nh == 0
    fclose(fid);
    error([' No header in trace constraint ',param.traces])
  end

  h_index=zeros(nh,1);
  counter=0;
  for ii=1:nh
    temp=find(ismember(seismic.header_info(:,1),lower(header(ii))));
    if ~isempty(temp)
      if length(temp) > 1
         disp([' More than one header with mnemonic "',header{ii},'" found when '])
         disp(' analyzing the trace-selection constraints.')
         disp(' Probably caused by explicitely requesting a header that is already ')
         disp(' read by default and then using it in the constraint.')
         pause(0)
         error(' Abnormal termination')
      end
      counter=counter+1;
      header(counter)=header(ii);
      h_index(counter)=temp;
    end
  end
  nh=counter;

else
  trace_select=2;
  idx_trace=param.traces;
  h_index=[];

     try
  max_trace=max(idx_trace);
     catch
  disp(' Apparently wrong syntax using keyword "traces"')
  disp(max_trace)
  error(' Abnormal termination')
     end

  header=[];
  nh=[];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [seismic,idx_time]=select_times(seismic,no_samples,fid,param)
% Function sets parameter that determine the time range to output

if iscell(param.times)
   param.times=cat(2,param.times{:});
end

%       Selection criteria for traces and times
if ~isempty(param.times)    % Compute index vector for time samples

   param.times(1)=ceil(param.times(1)/seismic.step)*seismic.step;
   if param.times(1) > 0
      seismic.first=param.times(1);
   end
   if param.times(2) < seismic.first
      fclose(fid);
      error([' Selected end time (',num2str(param.times(2)), ...
          ') smaller than start time (',num2str(seismic.first),')']);
   end
   if param.times(2) < seismic.last
      seismic.last=fix(param.times(2)/seismic.step)*seismic.step;
   end
    
   ita=ceil(seismic.first/seismic.step)+1;
   ite=floor(seismic.last/seismic.step)+1;
   idx_time=(max([1,ita]):min([ite,no_samples]))';
   if isempty(idx_time)
      fclose(fid);
      error(' No time samples selected. Check parameter "times"')
   end
else
  idx_time=(1:no_samples)';   % Keep all time samples
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [seismic,indices2,indices4,n_default_headers]= ...
       select_headers(seismic,units,indices2,indices4,i_lag,param)
% Function sets parameters which contol what headers to output

%global param

%         Select standard headers and add user-defined headers
% Descriptions for standard binary trace headers
seismic.header_info(:,3)=[...
      {'Trace sequence number within line'}; ...        % 4-byte headers 
      {'Original Field record number'}; ...
      {'Trace sequence number within original field record'}; ...
      {'Energy source point number'}; ...
      {'CDP number'}; ...
      {'Trace sequence number within CDP ensemble'}; ...
      {'Offset'}; ...
      {'Source depth below surface'}; ...
      {'Water depth at source'}; ...
      {'Water depth at receiver group'}; ...
      {'X coordinate of source'}; ...
      {'Y coordinate of source'}; ...
      {'Sourface elevation at source'}; ...
      {'X coordinate of receiver'}; ...
      {'Y coordinate of receiver'}; ...
      {'Receiver elevation'}; ...
      {'X-coordinate of CDP'}; ...
      {'Y-coordinate of CDP'}; ...
      {'In-line number'}; ...
      {'Cross-line number'}; ...
      {'Trace type (1=live,2=dead,3=dummy,...)'}; ...   % 2-byte headers
      {'Lag time between shot and recording start'}; ...
];

seismic.header_info(:,1)=[ ...  % 4-byte headers
      {'ds_seqno'}; ... 
      {'ffid'}; ...
      {'o_trace_no'}; ...  
      {'source'}; ...   
      {'cdp'}; ...      
      {'seq_cdp'}; ... 
      {'offset'}; ...     
      {'depth'}; ...  
      {'sou_h2od'}; ...  
      {'rec_h2od'}; ... 
      {'sou_x'}; ...     
      {'sou_y'}; ...      
      {'sou_elev'}; ...      
      {'rec_x'}; ...       
      {'rec_y'}; ...     
      {'rec_elev'}; ...
      {'cdp_x'};...
      {'cdp_y'}; ...
      {'iline_no'}; ...
      {'xline_no'}; ...
      {'trc_type'}; ...         % 2-byte headers
      {'lag'}; ...
];

n_default_headers=size(seismic.header_info,1);

seismic.header_info(1:n_default_headers,2)=deal({'n/a'});
seismic.header_info(7:16,2)={units};
seismic.header_info(i_lag,2)={'ms'};

%      Defaults for some commonly requested headers
defaults.iline_no={'n/a','In-line number'};
defaults.xline_no={'n/a','Cross-line number'};
defaults.cdp_x={units,'X-coordinate of CDP'}; 
defaults.cdp_y={units,'Y-coordinate of CDP'}; 
nh=length(param.headers);
add=cell(nh,5);
for ii=1:nh
   nh1=length(param.headers{ii});
   add(ii,1:nh1)=param.headers{ii};
   if nh1 == 5
     if isempty(add{ii,4}); add(ii,4)={'n/a'}; end
     if isempty(add{ii,5}); add(ii,5)={'not available'}; end
   else
     if sum(ismember(fieldnames(defaults),lower(add{ii,1})))
        add(ii,4:5)=getfield(defaults,lower(add{ii,1}));
%        add(ii,4:5)=defaults.(lower(add{ii,1}));
     else
        add(ii,4:5)={'n/a','not available'};
     end
   end
end
get.mnem=add(:,1);
get.first=add(:,2);
get.bytes=add(:,3);
get.units=add(:,4);
get.descriptions=add(:,5);
nheader=size(get.mnem,1); 

for jj=1:nheader
   if get.bytes{jj} == 2
      strt=(get.first{jj}-1)/2+1;
      if fix(strt) ~= strt
         error(['First byte for header ',get.mnem{jj}, ' is wrong'])
      end
      indices2=[indices2;strt];
   else
      strt=(get.first{jj}-1)/4+1;
      if fix(strt) ~= strt
         error(['First byte for header ',get.mnem{jj}, ' is wrong'])
      end
      indices4=[indices4;strt];
   end
end

seismic.header_info=cat(1,seismic.header_info, ...
   cat(2,lower(get.mnem(:)),get.units(:),get.descriptions(:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ascii=ebcdic2ascii(ebcdic)
% Function converts EBCDIC string to ASCII
% see http://www.room42.com/store/computer_center/code_tables.shtml
% Date Feb. 20, 2000;  written by E. R.
% INPUT
% ebcdic	EBCDIC string
% OUTPUT
% ascii		ASCII string
%		  ascii=ebcdic2ascii(ebcdic)

pointer= ...
[ 0    16    32    46    32    38    45    46    46    46    46    46   123   125    92    48
  1    17    33    46    46    46    47    46    97   106   126    46    65    74    46    49
  2    18    34    50    46    46    46    46    98   107   115    46    66    75    83    50
  3    19    35    51    46    46    46    46    99   108   116    46    67    76    84    51
  4    20    36    52    46    46    46    46   100   109   117    46    68    77    85    52
  5    21    37    53    46    46    46    46   101   110   118    46    69    78    86    53
  6    22    38    54    46    46    46    46   102   111   119    46    70    79    87    54
  7    23    39    55    46    46    46    46   103   112   120    46    71    80    88    55
  8    24    40    56    46    46    46    46   104   113   121    46    72    81    89    56
  9    25    41    57    46    46    46    46   105   114   122    46    73    82    90    57
 10    26    42    58    46    33   124    58    46    46    46    46    46    46    46    46
 11    27    43    59    46    36    44    35    46    46    46    46    46    46    46    46
 12    28    44    60    60    42    37    64    46    46    46    46    46    46    46    46
 13    29    45    61    40    41    95    39    46    46    91    93    46    46    46    46
 14    30    46    46    43    59    62    61    46    46    46    46    46    46    46    46
 15    31    47    63   124    94    63    34    46    46    46    46    46    46    46    46];
pointer=reshape(pointer,1,256);

ascii=pointer(ebcdic+1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = ibm2ieee (ibmf)

% Name:         ibm2ieee
% Abstract:     convert a matrix of IBM/360 32-bit floats
%               to IEEE doubles.
%
%               IBMF is the matrix of IBM/360 32-bit floats each
%               stored as a 32 bit unsigned big-endian integer
%               in a MATLAB double.
%
%               The format of a IBM/360 32-bit float is:
%
%                  sign 7-bit exponent  24 bit fraction
%                  The base is 16. The decimal point is to
%                  the left of the fraction. The exponent is
%                  biased by +64.
%
%               The basic idea is that you use floating point on
%               the various fields.
%
%               ieee = sign * 16 ^ (exponent - 64) * fraction / 16 ^ 6
%
% By:           Martin Knapp-Cordes
%               The MathWorks, Inc.
%
% Date(s):      Jun 95 - 28, 29

% $Revision: 1.2 $  $Date: 1995/06/29 14:50:03 $
% This M-file is not officially supported by The MathWorks, Inc.  Please use
% it as is, or modify it to your specific need
%
%Assuming you have a file, which contains IBM float 32 format binary data, called
%5702.seg, then you must use the following FOPEN and FREAD call the read the
%file:
%
%   fid = fopen('5702.seg','r','b');
%
%%
%% Read first data record - IBM/360 32-bit floating format
%%                          Read them as unsigned (32-bit) integers.
%% Convert to IEEE doubles using ibm2ieee.
%%
%% size - number of elements to read
%
%    ibm1 = fread(fid,size,'uint');
%    ieee1 = ibm2ieee(ibm1);
%----------------------------------------------------------------------------
%
    if (nargin ~= 1)
        error ('Wrong number of arguments.');
    elseif (isempty(ibmf))
        error ('Argument is an empty matrix');
    end
%
    aibmf = sprintf('%08x',ibmf);
%
% hexd(1) - 1st hex digit - first bit is sign, next 3 bits high order exponent
% hexd(2) - 2nd hex digit - bits of low order exponent
% hexd(3) - 3rd-8th hex digit - bits of fraction
%
    hexd = sscanf(aibmf,'%1x%1x%6x',[3,inf]);
    d = (1 - (hexd(1,:) >= 8) .* 2) .* ...
        16 .^ ((hexd(1,:) - (hexd(1,:) >= 8) .* 8) .* 16 + hexd(2,:) ...
        - 70).* hexd(3,:);
    d = reshape(d,size(ibmf));