长江大学matlab程序设计报告.doc

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 xibm2num(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
frac=bitand(b,uint32(hex2dec('00ffffff'))); % get mantissa from last 3 bytes
%format long
frac=double(frac);
frac=frac/2^24;
x=(1-2*sign).*16.^exp .* frac;
err = frac==0 & (exp~=-64 | sign~=0); % bias removal is incorrect for zero
if any(err)
% TMH 19/06/2003
disp(['WARNING : ',mfilename,' Invalid zero input --> Sure data are IBM FLOAT formatted ?'])
return;
end
err = frac~=0 & (frac<1/16 | frac>=1);
if any(err)
% TMH 19/06/2003
disp(['WARNING : ',mfilename,' Invalid mantissa input --> Sure data are IBM FLOAT formatted ?'])
return;
end

Wiggle函数代码：
Function
wiggle(x,t,Data,style,dmax,showmax,plImage,imageax,example_plot);

if (nargin==9);
np=3;
subplot(np,np,1); wiggle(Data);
subplot(np,np,2); wiggle(Data,dmax);
subplot(np,np,3); wiggle(x,t,Data);
subplot(np,np,4); wiggle(x,t,Data,style,dmax);
subplot(np,np,5); wiggle(x,t,Data,style,dmax,showmax);
subplot(np,np,6); wiggle(x,t,Data,style,dmax,showmax,plImage);
if isempty(dmax),
dmax=max(abs(Data(:)));
end
subplot(np,np,7); wiggle(x,t,Data,style,dmax,showmax,plImage,dmax./10);
return
end
showmax_def=100;
style_def='wiggle';
if nargin==1,
Data=x;
t=[1:1:size(Data,1)];
x=[1:1:size(Data,2)];
dmax=max(Data(:));
style=style_def;
showmax=showmax_def;
end
if nargin==2,
Data=x;
dmax=t;
t=[1:1:size(Data,1)];
x=[1:1:size(Data,2)];
style=style_def;
showmax=showmax_def;
end
if nargin==3,
style=style_def;
dmax=max(abs(Data(:)));
showmax=showmax_def;
end
if nargin==4,
dmax=max(abs(Data(:)));
showmax=showmax_def;
end
if nargin==5,
showmax=showmax_def;
end
if nargin<7
plImage=0;
end
if isempty(dmax),
% Set scaling factor dmax if empty
dmax=max(abs(Data(:)));
end
if isempty(showmax),
showmax=100;
end
if nargin==7,
imageax=[-1 1].*dmax;
end
if plImage==1,
imagesc(x,t,Data);
if (length(imageax)==1)
imageax=[-1 1].*abs(imageax);
end
caxis(imageax);
hold on
end
if (showmax>0)
if length(x)>1, dx=x(2)-x(1); end
ntraces=length(x);
d=ntraces/showmax;
if d<=1; d=1; end
d=round(d);
dmax=dmax/d;
LineWidth=0.0001;
EdgeColor=[0 0 0];
for i=1:d:length(x)
xt=dx*Data(:,i)'./dmax;
if (strmatch('VA',style)==1)
xt1=xt;xt1(find(xt1>0))=0;
f1=fill(x(i)+[xt,fliplr(xt1)],[t,fliplr(t)],[0 0 0]);
set(f1,'LineWidth',LineWidth)
set(f1,'EdgeColor',EdgeColor)
%set(f1,'EdgeAlpha',[0]); % GIVES ROCKY IMAGES
hold on
if (strmatch('VA2',style,'exact')==1)
xt2=xt;xt2(find(xt2<0))=0;
f2=fill(x(i)+[xt,fliplr(xt2)],[t,fliplr(t)],[1 0 0]);
set(f2,'LineWidth',LineWidth)

set(f2,'EdgeColor',EdgeColor)
%set(f2,'EdgeAlpha',[0])
end
else
plot(xt+x(i),t,'k-','linewidth',.05);
end
if i==1, hold on;end
end
end
hold off;
try
axis([min(x)-(x(2)-x(1)) max(x)+(x(2)-x(1)) min(t) max(t)])
catch
axis([min(x) max(x) min(t) max(t)])
end
set(gca,'Ydir','revers')

%getFrequencySpectrum函数得到频谱：
function [a,df]=getFrequencySpectrum(data,dt)
N=length(data);
t=[1:N]*dt;
nfft=2^nextpow2(N);
f=1/dt*(0:(nfft/2))/nfft;
W=fft(data,nfft);
a=abs(W);
a=a(1:(nfft/2+1));
a=a/max(a);
p=angle(W);
p=p(1:(nfft/2+1));
NN=length(a);
df=1/dt/(NN-1)/2;
figure;
subplot(3,1,1);
plot(dt*[0: N-1],data);
set(gca,'ygrid','on');
xlabel('时间/s');
ylabel('振幅/m');
title('第50道地震波波形图');
subplot(3,1,2);
plot(df*[0: NN-1],a);
set(gca,'ygrid','on');
xlabel('频率/Hz');
subplot(3,13);
plot(df*[0: NN-1],p);
set(gca,'ygrid','on');

滤波前波形及其频谱图：

滤波后波形及其频谱图：

干扰波和有效波对比图：
图形界面设计：
%modifysegy函数代码：
function modifysegy
if 0, callback; end
h=figure('name','Seismic_data Modify','units','pixels','NumberTitle','off',...
'menubar','none','toolbar','none','DoubleBuffer','on',...
'Units','pixels','resize','off',...
'WindowButtonDownFcn','callback(''buttondown'')',...
'WindowButtonMotionFcn','callback(''buttonmove'')',...
'WindowButtonUpFcn','callback(''buttonup'')',...
'CloseRequestFcn','callback(''close'')');
height=700; width=1200;
set(h,'position',[0 0 width height]);
movegui(h,'center');
uicontrol('style','pushbutton','Units','pixels', 'Position',[115 height-40 80 40],'string','wiggle','callback','callback(''wiggle'')');
uicontrol('style','pushbutton','Units','pixels', 'Position',[195 height-40 80 40],'string','imagesc','callback','callback(''imagesc'')');
uicontrol('style','pushbutton','Units','pixels', 'Position',[275 height-40 80 40],'string','close','callback','callback(''close'')');
c=get(gcf,'color');

%callback函数代码：
function callback(flag) % seismic data modify
global data head k jj filename h saveflag

switch flag
case 'wiggle'
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);
wiggle(seismic_data);
Xlabel('炮间距/m');
Ylabel('采样点时间/s');
Title('地震剖面图') ;

case '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);
imagesc(seismic_data);
Colormap（gary）；
Xlabel('炮间距/m');
Ylabel('采样点时间/s');
Title('地震剖面灰度图') ;

case 'close'
if saveflag==1
button = questdlg('save?','quit','yes','no','cancel','yes');
if strcmp(button,'yes')
callback('save');
delete(gcf);
return;
elseif strcmp(button,'no')
delete(gcf);
return;
elseif strcmp(button,'cancel')
return;
end
end
delete(gcf);
end

四、 总结和体会
MATLAB是一款商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。Matlab具有很强的实用性，在实际中发挥着很大的作用。在上次地震勘探原理实验中，通过编写运行matlab程序，实现对“地震波传播规律”，“地震波时距曲线”等规律的绘图成像，对课本的知识又有了更为直观，深入的认知。
对于这个题目，当我刚看到这个题目时，我感到很茫然，完全无从下手。由于自己的编程能力比较差，再加上自己上课没有紧跟老师的步伐，没有认真听讲，感到无从下手。还有就是c语言文件读写那部分当时就没有学扎实，基础比较薄弱。我通过网络视频，同时也看了老师给的大量实例，对文件读写程序有了一定的了解。通过参考其他程序，修改运行，还是有一定收货的。
1、通过这次实验，我初步了解了地震的数据格式——segy格式，初步了解了ibm2num函数、wiggle函数等函数的调用与功能。
2、通过这次实验以及课外自己的练习，我基本了解了一些图形界面制作的程序，对制作界面有了一些理解，在参考了一些程序后，对一些程序代码进行修改，最后也作出了自己的图形界面。
3、要想提高自己的matlab水平，不仅要在实验课上认真完成实验，在实验后，也要认真思考，多做练习，这样才能有所提高，有所进步。更应学会与人讨论交流，共同进步。还要学会充分利用网络视频资源，多学习，多实践。
总的来说自己在这次作业中学到了很多，自己的编程动手能力得到了，加深了对地震数据格式的认知，以后更应该多动手。多操作。
最后，感谢老师开设这门课程，给我了我们更多的机会，让我们从中收益。
15