segy.m

利用Matlab语言编写的读取segy数据的小程序

%   union SI si;//采样点数
%   union SP sp;//定义采样率
%   union Data data;//定义数据
%   union Line line;//线号
%   union Trace trace;//道号
%   union X_cor x_cor;//X坐标
%   union Y_cor y_cor;//Y坐标
%   Line_num     /最大线号
%   Trace_num    /最大道号

Line_num=300;
Trace_num=1;
fid=fopen('d:\mywork\sgyread\dfd.sgy','r');
   if ~fid
   {
       disp('can''t open file!');
       exit;
   }
   end
  fseek(fid,3220,'bof');  %读取采样点数
  SI=fread(fid,1,'int16','b')
  fseek(fid,3216,'bof');  %读取采样率
  SP=fread(fid,1,'int16','b')
  fseek(fid, 0, 'eof'); %计算总文件字节数
  file_n=ftell(fid);
  Tn =(file_n-3600)/(SI*4+240)  %计算道数
  fclose(fid);
 
 % seismic_data=zeros(Line_num,Trace_num,SI);  
  seismic_data=zeros(Line_num,SI);  
  
 fid=fopen('d:\mywork\sgyread\dfd.sgy','r');
   if ~fid
   {
       disp('can''t open file!');
       exit;
   }
   end
   for j=1:Tn
%读取线号
      fseek(fid,3600+(j-1)*240+(j-1)*SI*4+8,'bof'); 
      Line=fread(fid,1,'int32','b');
     if j==1
     Line_first=Line;
     end
%读取道号
      fseek(fid,3600+(j-1)*240+(j-1)*SI*4+20,'bof'); 
      Trace=fread(fid,1,'int32','b');
     if j==1
     Trace_first=Trace;
     end
%读取X坐标	
      fseek(fid,3600+(j-1)*240+(j-1)*SI*4+72,'bof'); 
      X_cor=fread(fid,1,'int32','b');
%读取Y坐标	
      fseek(fid,3600+(j-1)*240+(j-1)*SI*4+76,'bof'); 
      Y_cor=fread(fid,1,'int32','b');
%读取地震数据
      fseek(fid,3600+j*240+(j-1)*SI*4,'bof'); 
      seismic=fread(fid,[1,SI],'float32','b');
 %    seismic_data(Line+1-Line_first,Trace+1-Trace_first,:)=seismic;
      seismic_data(Line+1-Line_first,:)=seismic;
   end
 fclose(fid)   
 
  %中值滤波、去除野值
%Length :滤波器长度（样点个数）

  Length=5;
  seismic_filter=seismic_data;
  filter=zeros(1,Length);
  for i=1:Tn
      for j=(Length+1)/2:(SI-(Length-1)/2)
          filter=seismic_data(i,(j-(Length-1)/2):(j+(Length-1)/2));
          seismic_filter(i,j)=median(sort(filter));            
      end
   end
   seismic_data=seismic_filter; %滤波后的数据
   clear seismic_filter
   clear filter
 
 seismic_max=max(max(seismic_data)); %最大值
 seismic_min=min(min(seismic_data)); %最小值
 
%  fid2=fopen('d:\mywork\sgyread\data.txt','w');
%  if ~fid2   %新建文本文件存储数据
%   {
%       disp('can''t open file!');
%       exit;
%   }
%  end 
%  fprintf(fid2,'%8.6f\n',seismic_data);
%  fclose(fid2);

 fid4=fopen('d:\mywork\sgyread\data.grd','w');
 if ~fid4   %新建文本文件存储数据
  {
      disp('can''t open file!');
      exit;
  }
 end 
 
%   plot(seismic_data)
% surf(x,y,z)
%  colormap(cool)
%  xlable('x'),ylable('y'),zlable('z')

 %输出grd格式
  Dx=25.0;
  fprintf(fid4,'DSAA\n');
  fprintf(fid4,'%d %d\n',Tn,SI);
  temp1=1;
  temp2=Tn;
  fprintf(fid4,'%f %f\n',temp1,temp2);
  Dt=2.0;
  temp1=0;
  temp2=(SI-1)*Dt;
  fprintf(fid4,'%f %f\n',-temp2,-temp1); 
%   temp1=0;
%   temp2=128;
  fprintf(fid4,'%f %f\n',seismic_min,seismic_max); 
  for j=1:SI
      for i=1:Tn   
          fprintf(fid4,'%f ',seismic_data(i,SI-j+1));
        end
       fprintf(fid4,'\n'); 
  end
  fclose(fid4);