长江大学matlab程序设计报告
- 格式：doc
- 大小：900.03 KB
- 文档页数：14
MATLAB程序设计报告 学 院： 地球物理与石油资源学院 班 级： 物探 姓 名： 学 号： 班内编号： 指导教师： 完成日期： 2015年6月30日 2
一、课程设计题目和要求 题目：地震数据（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
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')); 4
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 5
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, 6
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);
