长江大学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*onesSI Sp,其中n为道号，由于没有读出道间距我就假设道间距为100，然后n*100*onesSI Sp+ SI Sp就相当于每道的波形错开一个道间距，再显示。第二，直接调用wiggle或imagesc显示。三、源代码主程序：segyid=fopen(101.seg,rb); if segyid disp(cant 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(cant 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); endfclose(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 bitsign=double(sign);% format hexexp=bitand(b,uint32(hex2dec(7f000000); % get exponent from first byte, last 7 bitsexp=bitshift(exp,-24);%format longexp=double(exp)- 64; % remove bias from exponent %format hexfrac=bitand(b,uint32(hex2dec(00ffffff); % get mantissa from last 3 bytes%format longfrac=double(frac);frac=frac/224;x=(1-2*sign).*16.exp .* frac;err = frac=0 & (exp=-64 | sign=0); % bias removal is incorrect for zeroif any(err) % TMH 19/06/2003 disp(WARNING : ,mfilename, Invalid zero input - Sure data are IBM FLOAT formatted ?) return; enderr = frac=0 & (frac=1);if any(err) % TMH 19/06/2003 disp(WARNING : ,mfilename, Invalid mantissa input - Sure data are IBM FLOAT formatted ?) return; end Wiggle函数代码： Functionwiggle(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); returnendshowmax_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;endif nargin=2, Data=x; dmax=t; t=1:1:size(Data,1); x=1:1:size(Data,2); style=style_def; showmax=showmax_def;endif nargin=3, style=style_def; dmax=max(abs(Data(:); showmax=showmax_def;endif nargin=4, dmax=max(abs(Data(:); showmax=showmax_def;endif nargin=5, showmax=showmax_def;endif nargin0) if length(x)1, dx=x(2)-x(1); end ntraces=length(x); d=ntraces/showmax; if d0)=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(xt20)=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 endendhold 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)endset(gca,Ydir,revers)%getFrequencySpectrum函数得到频谱：function a,df=getFrequencySpectrum(data,dt)N=length(data);t=1:N*dt;nfft=2nextpow2(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,1,3);plot(df*0: NN-1,p);set(gca,ygrid,on);滤波前波形及其频谱图：滤波后波形及其频谱图：干扰波和有效波对比图： 图形界面设计：%modifysegy函数代码：function modifysegyif 0, callback; endh=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 modifyglobal data head k jj filename h saveflagswitch flagcase wiggle segyid=fopen(101.seg,rb); if segyid disp(cant 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(cant 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(cant 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(cant 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