% 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);
评论0