长江大学matlab程序设计报告

MATLAB程序设计报告

学院：地球物理与石油资源学院

班级：物探

姓名：

学号：

班内编号：

指导教师：

完成日期： 2015年6月30日

一、课程设计题目和要求

题目：地震数据（segy格式数据）读写

日期：2014年3月27日至2014年4月27日

主要内容：

1、理解segy格式

2、二维动态数组

3、将读取的数据显示出来

要求：

1、源代码格式规范

2、运行结果以SURFER的图片形式提交

3、独立完成，严禁抄袭，逐一验收

二、编程思路分析：

1、segy格式的了解：文件前3600字节为格式属性等的描述，跳过3216字节可读取采样率，跳过3220字节可读取采样数等。3600字节之后为每一道检波器检测的信息，每道前240字节为每一道的描述，之后数据才是实际测得数据！所以运用fopen,freed,fseek,ftell,fclose可以获得采样率，采样数等信息。

2、读取检波器测的数据并保存在一个矩阵[SI Sp]中，其中SI为每道采样数，Sp为道数。因为读到的数据是ibm格式，所以我们还要将ibm格式的数据转换成数值格式（调用ibm2num函数）。这里就有两种选择：第一，将数据读完后再进行格式转换，第二，边读数据边转换格式。经过测试第二种方案速度更快所以选择第二种方案。

3、数据的显示，这里也有两种方案：第一，如果直接用读到的矩阵[SI Sp]进行画图会出现一个问题，那就是各道的波形就会叠加在一起，根本得不到地震剖面。所以我们要让各道波形分开，通过n*100*ones[SI Sp],其中n为道号，由于没有读出道间距我就假设道间距为100，然后n*100*ones[SI Sp]+ [SI Sp]就相当于每道的波形错开一个道间距，再显示。第二，直接调用wiggle或imagesc显示。

三、源代码

主程序：

segyid=fopen('101.seg','rb');

if ~segyid

{

disp('can''t open file!');

exit;

} ;

end

fseek(segyid,3220,'bof'); %读取采样点数

SI=fread(segyid,1,'int16','b') ;

fseek(segyid,3216,'bof'); %读取采样率

Sp=fread(segyid,1,'int16','b') ;

fseek(segyid,0,'eof'); %计算总文件字节数

BytesInFile=ftell(segyid);

Tn=(BytesInFile-3600)/(SI*4+240); %计算道数

fclose(segyid);

seismic_data=zeros(SI,Tn);

segyid=fopen('101.seg','r','b');

if ~segyid

{

disp('can''t open file!');

exit;

} ;

end

%读取地震数据

for n=1:Tn

fseek(segyid,3600+n*240+(n-1)*SI*4,'bof');

seismic=fread(segyid ,[1 SI],'uint32','b');

seismic_data(:,n)=ibm2num(uint32(seismic'));

end

fclose(segyid);

% h=0:Sp/1000:Sp/1000*(SI-1);

% for n=1:Tn

% q=100*n*ones(SI,Tn);

% c=q(:,n);

% seismic_data11=seismic_data1(:,n);

% seismic_data2=seismic_data11+c;

% end;

% plot(seismic_data2,h);

% hold on;

wiggle(seismic_data);

Xlabel('炮间距');

Ylabel('采样点时间/s');

Title('地震剖面图') ;

%选择第50道得到它的波形和频谱

%零相位滤波filter一般滤波

data=seismic_data(:,50)';

t = 0:1/250:4.496;

[a,df]=getFrequencySpectrum(data,1/250);

[b,a]=butter(4,20/(250/2),'low'); %低通

yy = filtfilt(b,a,data); %零相位滤波filter一般滤波

figure;

plot(t,data-500,'r',t,yy+500,'g');

xlabel('时间/s');

ylabel('振幅/m');

legend('干绕波波形','有效波波形',1);

title('第50道地震波干绕波和有效波波形图');

getFrequencySpectrum(yy,1/250);

Ibm2num函数代码：

function x=ibm2num(b)

x=repmat(NaN,size(b));

sign=bitget(b,32); % get sign from first bit sign=double(sign);

% format hex

exp=bitand(b,uint32(hex2dec('7f000000'))); % get exponent from first byte, last 7 bits

exp=bitshift(exp,-24);

%format long

exp=double(exp)- 64; % remove bias from exponent

%format hex