write_segy_file.m

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

     two_bytes(22)=0;
    end
  end
end

fwrite(fid,[jobid lineid reelid],'int32');
fwrite(fid,two_bytes,'int16');
if param.print, disp('Binary reel header written'), end;

%   Write headers and traces
nh=size(seismic.header_info,1);
start=zeros(nh,1);
bytes=zeros(nh,1);
index=zeros(nh,1);

kk=1;
[start,bytes,index,kk]=set_parameters(seismic,'ds_seqno',1,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'ffid',9,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'o_trace_no',13,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'source',17,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'cdp',21,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'seq_cdp',25,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'trc_type',29,2,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'offset',37,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'depth',49,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'sou_h2od',61,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'rec_h2od',65,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'sou_x',73,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'sou_y',77,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'sou_elev',45,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'rec_x',81,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'rec_y',85,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'rec_elev',41,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'gen_scale_factor',69,2,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'gen_scale_factor',71,2,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'gen_scale_factor',89,2,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'number_of_samples',115,2,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'samp_int',117,2,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'lag',109,2,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'cdp_x',181,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'cdp_y',185,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'iline_no',189,4,start,bytes,index,kk,param.print);
[start,bytes,index,kk]=set_parameters(seismic,'xline_no',193,4,start,bytes,index,kk,param.print);

%       Add user-specified headers
for ii=1:length(param.headers)
   if length(param.headers{ii}) ~= 3
      disp([char(13),'  Error in specification of headers to be written to file: ', ...
          cell2str(param.headers{ii}(1))])
      error(' Probably insufficient number of parameters')
   end
   [start,bytes,index,kk]=set_parameters(seismic,param.headers{ii}{1}, ...
        param.headers{ii}{2},param.headers{ii}{3},start,bytes,index,kk,param.print);
end

start(kk:end)=[];    % Remove unneeded, previously reserved array elements
bytes(kk:end)=[];    % Remove unneeded, previously reserved array elements
index(kk:end)=[];    % Remove unneeded, previously reserved array elements

for ii=1:ntr
  write_trace(fid,seismic.traces(:,ii),seismic.headers(index,ii),start,bytes,fidx);
end
fclose(fid);

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

function [start,bytes,index,kk]=set_parameters(seismic,mnem,sb,nb,start,bytes,index,kk,iprint)
% Function sets values start(index),bytes(index) to the starting location
% of the header mnemonic and increments index
% seismic    seismic structure
% mnem    header mnemonic whose starting location in the four-byte representation 
%         (nb=4) or the two-byte representation (nb=2) needs to be set
% sb      starting byte
% nb      number of bytes
% start   array where the index for the two-byte and four-byte header is stored
% bytes   array where the number of bytes required by this header is stored
% index   index into row of header which contains the header values of 
%         the mnemonic mnem
% kk      next location in arrays start and bytes
% OUTPUT
% start   updated start array
% bytes   updated bytes array
% index   updated index
% kk      input index incremented by 1
% iprint  print-out control parameter: 1 printout, 0 no printout

global S4M

if S4M.case_sensitive
   idx=find(ismember(seismic.header_info(:,1),mnem));
else
   idx=find(ismember(lower(seismic.header_info(:,1)),mnem));
end
if isempty(idx)
  if iprint, disp(['Header ',mnem,' not found and not written to file ']), end
  return, 
end
index(kk)=idx;
bytes(kk)=nb;

if nb == 2
  temp=fix((sb-1)/2);
  if temp*2 ~= sb-1
    error(['Starting byte location for header ',mnem,'(',num2str(sb),' is not odd'])
  else
    start(kk)=temp+1;
  end
else
    temp=fix((sb-1)/4);
  if temp*4 ~= sb-1
    error(['Starting byte location for header ',mnem,'(',num2str(sb),' is not 1 plus a multiple of 4'])
  else
    start(kk)=temp+1;
  end  
end
kk=kk+1;

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

function write_trace(fid,trace,headers,start,bytes,fidx)
% Function writes one seismic trace together with its header to file
% fid     file identification number
% trace   seismic trace
% headers header values associated with the trace
% start   starting four-byte or two-byte index in the 240-byte header 
%         for each header value
% bytes   number of bytes for each header value
% fidx    format index;
%         fidx=0; no conversion
%         fidx=1; conversion to IMB floating point format

% nsamp=length(trace);
nh=length(headers);

fbytes=zeros(60,1);
tbytes=zeros(120,1);
for ii=1:nh
   if bytes(ii) == 4
      fbytes(start(ii))=headers(ii);
   else
      tbytes(start(ii))=headers(ii);
   end
end
fwrite(fid,fbytes( 1: 7),'int32');
fwrite(fid,tbytes(15:18),'int16');
fwrite(fid,fbytes(10:17),'int32');
fwrite(fid,tbytes(35:36),'int16');
fwrite(fid,fbytes(19:22),'int32');
fwrite(fid,tbytes(45:90),'int16');
fwrite(fid,fbytes(46:60),'int32');

% Write trace
if fidx == 0
   fwrite(fid,trace,'float32');
else
   error('IBM floating point not yet implemented')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ascii_header=make_header
% Function creates ASCII version of standard EBCIC header of SEG-Y format

ascii_header=char(...
'C 1 CLIENT                        COMPANY                       CREW NO', ...         
'C 2 LINE            AREA                        MAP ID ', ...                          
'C 3 REEL NO           DAY-START OF REEL     YEAR      OBSERVER', ...                   
'C 4 INSTRUMENT: MFG            MODEL            SERIAL NO', ...                       
'C 5 DATA TRACES/RECORD        AUXILIARY TRACES/RECORD         CDP FOLD', ...           
'C 6 SAMPLE INTERNAL         SAMPLES/TRACE       BITS/IN      BYTES/SAMPLE', ...        
'C 7 RECORDING FORMAT        FORMAT THIS REEL        MEASUREMENT SYSTEM', ...           
'C 8 SAMPLE CODE: FLOATING PT     FIXED PT     FIXED PT-GAIN     CORRELATED ', ...      
'C 9 GAIN  TYPE: FIXED     BINARY     FLOATING POINT     OTHER ', ...                   
'C10 FILTERS: ALIAS     HZ  NOTCH     HZ  BAND    -     HZ  SLOPE    -    DB/OCT ', ...  
'C11 SOURCE: TYPE            NUMBER/POINT        POINT INTERVAL', ...                   
'C12     PATTERN:                           LENGTH        WIDTH', ...                   
'C13 SWEEP: START     HZ  END     HZ  LENGTH      MS  CHANNEL NO     TYPE', ...         
'C14 TAPER: START LENGTH       MS  END LENGTH       MS  TYPE', ...                      
'C15 SPREAD: OFFSET        MAX DISTANCE        GROUP INTERVAL', ...                    
'C16 GEOPHONES: PER GROUP     SPACING     FREQUENCY     MFG          MODEL', ...        
'C17     PATTERN:                           LENGTH        WIDTH', ...                   
'C18 TRACES SORTED BY: RECORD     CDP     OTHER', ...                                  
'C19 AMPLITUDE RECOVEY: NONE      SPHERICAL DIV       AGC    OTHER', ...                
'C20 MAP PROJECTION                      ZONE ID       COORDINATE UNITS', ...           
'C21 PROCESSING:', ...                                                                  
'C22 PROCESSING:', ...                                                                  
'C23 ', ...                                                                             
'C24 ', ...                                                                             
'C25 ', ...                                                                             
'C26 ', ...                                                                             
'C27 ', ...                                                                             
'C28 ', ...                                                                             
'C29 ', ...                                                                             
'C30 ', ...                                                                             
'C31 ', ...                                                                             
'C32 ', ...                                                                             
'C33 ', ...                                                                             
'C34 ', ...                                                                             
'C35 ', ...                                                                             
'C36 ', ...                                                                             
'C37 ', ...                                                                             
'C38 ', ...                                                                             
'C39 ', ...                                                                             
'C40 END EBCDIC')';

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

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

pointer =  ...
  [0  16  64 240 124 215 125 151  75  75  75  75  75  75  75  75
   1  17  90 241 193 216 129 152  75  75  75  75  75  75  75  75
   2  18 127 242 194 217 130 153  75  75  75  75  75  75  75  75
   3  19 123 243 195 226 131 162  75  75  75  75  75  75  75  75
   4  20  91 244 196 227 132 163  75  75  75  75  75  75  75  75
   5  21 108 245 197 228 133 164  75  75  75  75  75  75  75  75
   6  22  80 246 198 229 134 165  75  75  75  75  75  75  75  75
   7  23 125 247 199 230 135 166  75  75  75  75  75  75  75  75
   8  24  77 248 200 231 136 167  75  75  75  75  75  75  75  75
   9  25  93 249 201 232 137 168  75  75  75  75  75  75  75  75
  10  26  92 122 209 233 145 169  75  75  75  75  75  75  75  75
  11  27  78  94 210 173 146 192  75  75  75  75  75  75  75  75
  12  28 107  76 211 224 147 106  75  75  75  75  75  75  75  75
  13  29  96 126 212 189 148 208  75  75  75  75  75  75  75  75
  14  30  75 110 213  95 149 161  75  75  75  75  75  75  75  75
  15  31  97 111 214 109 150  75  75  75  75  75  75  75  75  75];
pointer=pointer(:);

ebcdic=pointer(ascii+1);