MATLAB_SWPU第七章_地震处理相关知识PPT参考课件.ppt

1、第七章 地震处理相关知识,7.1 地震数据格式 7.2 常用处理 7.3 地震属性提取,1,7.1 地震数据格式,常见的地震数据体有SEG格式和NOSEG 格式， SEG是（ The Society of Exploration Geophysicists）的缩写， 而SEG格式有多种，其中最常用的就是SEGY格式地震数据体文件。它有3200个字节的图头（也叫卷头），400字节的文件头，240字节的道头，道数据，.道头，道数据，很多个。图头信息我们一般不关心，不太常用。文件头里有些信息可以用到，例如数据格式了，样品数了，采样间隔了。而道头在实际加载的过程中最常用，其中最常用的是X、Y、道号和线

2、号。但在实际的数据中它们存放的位置不一定相同，要具体数据体具体分析，可以用一个SEGY格式分析器进行转换,2,7.1 地震数据格式,SEGY数据由文件头和数据体组成 文件头总长度为3600字节，分两部分。 文件头第一部分 长度： 3200bytes；组成： 80bytes*40；特性：EBCDIC字符集，参数卡，需要转换为ASCII码后才能显示。 文件头第二部分 长度：400bytes；数据类型：32位、16位的整型；特性：二进制头，记录数据体信息。 数据体由多个数据道组成。每道数据分两部分：道头、采样数据，数据类型：32位的浮点型 ，占4个字节。,3,7.1 地震数据格式,工作站SEGY数据

3、存储格式有两种：IEEE和IBM IEEE和IBM的整型存储与微机格式的存储不同之处，IEEE和IBM的高字节在前、低字节在后，即BigEndian，微机则是低字节在前、高字节在后，即LittleEndian 。,4,文件头格式,5,6,道头格式,7,8,9,10,例 读取地震数据并绘图,clear all indata=altreadsegy(mianbomodle.sgy);%读取数据 Tr=length(indata(1,:);%道数 T=length(indata(:,1); %采样点数 s_wplot(indata);%绘剖面,11,12,例 提取第10道单独显示，并且加入随机噪声,

4、clear all indata=altreadsegy(mianbomodle.sgy); Tr=length(indata(1,:); T=length(indata(:,1); data=indata(:,10);%提取第10道来显示 temp=randn(T,1);%生成随机数据 data=data+temp; s_wplot(data); %s_cplot(data);%变面积显示,13,14,例 在第10道中加入随机噪声,clear all indata=altreadsegy(mianbomodle.sgy); Tr=length(indata(1,:); T=length(in

5、data(:,1); temp=randn(T,1);%加入随机噪声 indata(:,10)=indata(:,10)+temp; s_wplot(indata); %s_cplot(data);%变面积显示,15,16,7.2 常用处理-滤波,clear all indata=altreadsegy(mianbomodle.sgy); Tr=length(indata(1,:); T=length(indata(:,1); data=indata(:,10);%提取第10道来显示 temp=randn(T,1);%生成随机数据 data=data+temp; fdata=fft(data)

6、; plot(abs(fdata);,17,18,例 低通滤波器实现, fdata=fft(data); fT=length(fdata); for j=100:fT %低通滤波器 fdata(j)=0; end plot(abs(fdata);,19,20,例 低通滤波器实现-反变换, fdata=fft(data); fT=length(fdata); for j=100:fT %低通滤波器 fdata(j)=0; end plot(abs(fdata); data2=ifft(fdata);%逆变换 s_wplot(data); s_wplot(data2);,21,22,例 高通滤波器

7、实现, fdata=fft(data); fT=length(fdata); for j=1:200 %高通滤波器 fdata(j)=0; end plot(abs(fdata);,23,24,25,例 带通滤波器实现, fdata=fft(data); fT=length(fdata); for j=1:fT %带通滤波器 if(j300|j40) fdata(j)=0; end end plot(abs(fdata);,26,27,28,常用处理-子波提取,clear all read_segy_file(mig.sgy)； seismic=s_data; sdata=s_seismic2

8、wavelet(seismic,color,blue,type,zero-phase,wlength,80); s_wplot(sdata);,zero-phase zero-phase wavelet min-phase minimum-phase wavelet max-phase maximum-phase wavelet,29,30,常用处理-子波提取,clear all read_segy_file(mig.sgy); seismic=s_data; bt=s_seismic2wavelet(seismic,color,blue,type,zero-phase,wlength,80)

9、; kesai=zeros(1,200); kesai(12)=-0.2; kesai(33)=0.4; kesai(37)=-0.15; kesai(52)=-0.5; kesai(72)=0.35; kesai(87)=-0.1; plot(bt.traces); st=conv(bt.traces,kesai); s_wplot(st);,31,常用处理-希尔伯特变换,clear all seismic=read_segy_file(mig.sgy); sdata=s_hilbert(seismic); s_cplot(sdata); s_cplot(seismic);,32,33,常用

10、处理-能量谱,clear all seismic=read_segy_file(mig.sgy); sdata=s_power_spectrum(seismic); s_cplot(sdata);,34,35,7.3 地震属性提取-能量,clear all tempdata=altreadsegy(mig.sgy);%读地震数据 indata=tempdata;%矩阵转置 T=length(indata(1,:);%采样率 Tr=length(indata(:,1); %道数 lTr,ltime=textread(layer.txt,%f %f,headerlines,0)；%读取层位时间 t

11、window=5;%设置时窗大小 tspace=25;%采样间隔 for i=1:Tr energy(i)=0; for j=1:T if(jltime(i)/tspace-twindow,36,37,属性提取-FFT分析,clear all tempdata=altreadsegy(mig.sgy);%读地震数据 indata=tempdata;%矩阵转置 T=length(indata(1,:);%采样率 Tr=length(indata(:,1); %道数 lTr,ltime=textread(layer.txt,%f %f,headerlines,0);%读取层位时间 twindow=

12、5;%设置时窗大小 tspace=25;%采样间隔 for i=1:Tr for j=1:T if(jltime(i)/tspace-twindow end end end,38,属性提取-FFT分析,fdata=fft(data); x1=length(data(1,:); y1=length(data(:,1); x=1:x1; y=1:y1; z=abs(fdata); mesh(y,x,z); colorbar; view(2);,39,40,7.4 生成理论地震模型,function wavelet f=30; tdeta=0.004; N=30; nt=0:tdeta:N*tdet

13、a; bt=(1-2*f*f*pi*pi.*nt.*nt).*exp(-pi*pi*f*f.*nt.*nt); kesai=zeros(1,200); kesai(12)=-0.2; kesai(33)=0.4; kesai(37)=-0.15; kesai(52)=-0.5; kesai(72)=0.35; kesai(87)=-0.1; st=conv(bt,kesai);,41,nRandRange=0.2; ntrand=-nRandRange+rand(1,length(st)*nRandRange*2;%随机数组 str=sprintf(随机信号 n(t) 范围 %.1g到%.1g

14、,-nRandRange,nRandRange); %xt=st+ntrand; xt=st; nfft=2ceil(log(length(xt)/log(2);%产生大于xt的最接近2的n次方的整数 xw=fft(xt,nfft); plot(nt,bt); title(地震子波 b(t); axis(0 0.12 -1 1);,42,subplot(2,2,1); plot(ntrand); title(str); subplot(2,2,2); stem(kesai); title(反射系数序列 xi(t); subplot(2,2,3); plot(linspace(0,N*tdeta

15、,length(xt),xt); title(x(t); axis(0 0.12 -1 1); subplot(2,2,4); plot(abs(xw); title(|X(W)|);,43,7.5 相关计算,function XiangGuan A1=1;tdetaX=0.004;f1=35;MX=150; A2=1;beta=100;tdetaY=0.004;f2=50;MY=150; ntX=0:tdetaX:MX*tdetaX; ntY=0:tdetaY:MY*tdetaY; xt=A1*sin(2*pi*f1.*ntX); yt=A2*exp(-beta.*ntY.*ntY).*co

16、s(2*pi*f2.*ntY); rxy=xcorr(xt,yt);%相关计算 ryx=xcorr(yt,xt); r_xy=conv(-xt,yt);%离散序列做卷积运算,44,subplot(2,1,1); plot(ntX,xt); title(x(n); subplot(2,1,2); plot(ntY,yt); title(y(n); figure(201); subplot(3,1,1); plot(linspace(0,MX*tdetaX,length(rxy),rxy); title(bffontsize12gamma_xy(tau); subplot(3,1,2); plot

17、(linspace(0,MX*tdetaX,length(ryx),ryx); title(bffontsize12gamma_yx(tau); subplot(3,1,3); plot(linspace(0,MX*tdetaX,length(ryx),r_xy); title(bffontsize12gamma_xy(-tau);,45,7.6 复杂的滤波程序,function filter A1=1;A2=1;A3=0.5;%振幅初值 f1=25;f2=40;f3=5;f4=80;%频率初值 beta=100;%衰减系数 tdeta=0.004;%时域间隔 N=200;%取样点个数 M=6

18、0;%时间域滤波因子长度 nt=0:tdeta:N*tdeta;%生成时间数组 xt=A1*exp(-beta*nt.*nt).*sin(2*pi*f1*nt)+A2*exp(-beta*nt.*nt).*sin(2*pi*f2*nt)+A3*(sin(2*pi*f3*nt)+cos(2*pi*f4*nt);%地震记录 st=A1*exp(-beta*nt.*nt).*sin(2*pi*f1*nt)+A2*exp(-beta*nt.*nt).*sin(2*pi*f2*nt);%无干扰信号的地震记录,46,计算时间域的信号,w1=20*2*pi;w2=45*2*pi; w0=w1/2+w2/2;

19、detaW=w2/2-w1/2; hdeta=tdeta:tdeta:(M-1)*tdeta;%M-1个点 h0=2*detaW/pi;滤波因子起始振幅 h1=2/pi./hdeta.*cos(w0*hdeta).*sin(detaW*hdeta);%滤波因子 h=h0,h1;%滤波因子 yt=conv(xt,h);%时域滤波,47,计算时间域的信号,figure(101);%图名 subplot(3,1,1); plot(nt,xt); title(时域：x(nDelta);xlabel(t);ylabel(A); subplot(3,1,2); plot(linspace(0,M*tdet

20、a,length(h),h); title(滤波因子：h(nDelta);xlabel(t);ylabel(A); subplot(3,1,3); plot(linspace(0,N*tdeta,length(yt),yt); title(时域滤波：y(nDelta);xlabel(t);ylabel(A);,48,频率域滤波后结果,nfft=2ceil(log(length(xt)/log(2);%产生大于xt的最接近2的n次方的整数 xw=fft(xt,nfft);%时域信号付氏变换 nf=1/nfft/tdeta*(0:nfft-1);%必须减1，因为从0开始，求频率样点 fstart=

21、20;fend=45;%低截频、高截频 h1w=(nf-fstart) = 0 ,49,频率域滤波后结果,figure(102); subplot(3,1,1); plot(nf,abs(xw); title(str,|X(W)|);xlabel(f);ylabel(A); axis(0 250 0 60); subplot(3,1,2); plot(nf,abs(y1w); title(单门频域滤波后：|Y(W)|);xlabel(f);ylabel(A); axis(0 250 0 60); subplot(3,1,3); plot(nf,abs(y2w); title(双门频域滤波后：|

22、Y(W)|);xlabel(f);ylabel(A); axis(0 250 0 60);,50,设计的滤波器频域和时间域结果,figure(103);%设计的单双门滤波因子 subplot(2,1,1); plot(nf,abs(h1w); title(str,单门：|H(W)|);xlabel(f);ylabel(A); subplot(2,1,2); plot(nf,abs(h2w); title(双门：|H(W)|);xlabel(f);ylabel(A);,51,设计的滤波器频域和时间域结果,figure(104);%单双门滤波效果与时域褶积滤波效果比较 subplot(3,1,1)

23、; plot(linspace(0,N*tdeta,length(y1w),real(ifft(y1w); title(str,单门频域滤波时域：y(t);xlabel(f);ylabel(A); axis(0 0.8 -2 2); subplot(3,1,2); plot(linspace(0,N*tdeta,length(y2w),real(ifft(y2w); title(双门频域滤波时域：y(t);xlabel(f);ylabel(A); axis(0 0.8 -2 2); subplot(3,1,3); plot(linspace(0,N*tdeta,length(yt),yt); title(时域滤波：y(nDelta);xlabel(t);ylabel(A);,52,时域褶积滤波与原始信号在频域进行比较,nfft=2ceil(log(length(yt)/log(2);%产生大于xt的最接近2