MATLAB_SWPU第七章_地震处理相关知识

第七章地震处理相关知识7.1地震数据格式7.2常用处理7.3地震属性提取7.1地震数据格式常见的地震数据体有SEG格式和NOSEG格式，SEG是（TheSocietyofExplorationGeophysicists）的缩写，而SEG格式有多种，其中最常用的就是SEGY格式地震数据体文件。它有3200个字节的图头（也叫卷头），400字节的文件头，240字节的道头，道数据，.............道头，道数据，很多个。图头信息我们一般不关心，不太常用。文件头里有些信息可以用到，例如数据格式了，样品数了，采样间隔了。而道头在实际加载的过程中最常用，其中最常用的是X、Y、道号和线号。但在实际的数据中它们存放的位置不一定相同，要具体数据体具体分析，可以用一个SEGY格式分析器进行转换7.1地震数据格式SEGY数据由文件头和数据体组成文件头总长度为3600字节，分两部分。文件头第一部分长度：3200bytes；组成：80bytes*40；特性：EBCDIC字符集，参数卡，需要转换为ASCII码后才能显示。文件头第二部分长度：400bytes；数据类型：32位、16位的整型；特性：二进制头，记录数据体信息。数据体由多个数据道组成。每道数据分两部分：道头、采样数据，数据类型：32位的浮点型，占4个字节。7.1地震数据格式工作站SEGY数据存储格式有两种：IEEE和IBMIEEE和IBM的整型存储与微机格式的存储不同之处，IEEE和IBM的高字节在前、低字节在后，即BigEndian，微机则是低字节在前、高字节在后，即LittleEndian。文件头格式道头格式例读取地震数据并绘图clearallindata=altreadsegy(‘mianbomodle.sgy’);%读取数据Tr=length(indata(1,:));%道数T=length(indata(:,1));%采样点数s_wplot(indata);%绘剖面例提取第10道单独显示，并且加入随机噪声clearall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);%变面积显示例在第10道中加入随机噪声clearallindata=altreadsegy('mianbomodle.sgy');Tr=length(indata(1,:));T=length(indata(:,1));temp=randn(T,1);%加入随机噪声indata(:,10)=indata(:,10)+temp;s_wplot(indata);%s_cplot(data);%变面积显示7.2常用处理-滤波clearall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);plot(abs(fdata));例低通滤波器实现………fdata=fft(data);fT=length(fdata);forj=100:fT%低通滤波器fdata(j)=0;endplot(abs(fdata));例低通滤波器实现-反变换………fdata=fft(data);fT=length(fdata);forj=100:fT%低通滤波器fdata(j)=0;endplot(abs(fdata));data2=ifft(fdata);%逆变换s_wplot(data);s_wplot(data2);例高通滤波器实现………fdata=fft(data);fT=length(fdata);forj=1:200%高通滤波器fdata(j)=0;endplot(abs(fdata));例带通滤波器实现………fdata=fft(data);fT=length(fdata);forj=1:fT%带通滤波器if(j300||j40)fdata(j)=0;endendplot(abs(fdata));常用处理-子波提取clearallread_segy_file(‘mig.sgy’)；seismic=s_data;sdata=s_seismic2wavelet(seismic,{'color','blue'},{'type','zero-phase'},{'wlength',80});s_wplot(sdata);'zero-phase'zero-phasewavelet'min-phase'minimum-phasewavelet'max-phase'maximum-phasewavelet常用处理-子波提取clearallread_segy_file('mig.sgy');seismic=s_data;bt=s_seismic2wavelet(seismic,{'color','blue'},{'type','zero-phase'},{'wlength',80});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');常用处理-希尔伯特变换clearallseismic=read_segy_file('mig.sgy');sdata=s_hilbert(seismic);s_cplot(sdata);s_cplot(seismic);常用处理-能量谱clearallseismic=read_segy_file('mig.sgy');sdata=s_power_spectrum(seismic);s_cplot(sdata);7.3地震属性提取-能量clearalltempdata=altreadsegy(‘mig.sgy’);%读地震数据indata=tempdata‘;%矩阵转置T=length(indata(1,:));%采样率Tr=length(indata(:,1));%道数[lTr,ltime]=textread(‘layer.txt’,‘%f%f’,‘headerlines’,0)；%读取层位时间twindow=5;%设置时窗大小tspace=25;%采样间隔fori=1:Trenergy(i)=0;forj=1:Tif(jltime(i)/tspace-twindow&&jltime(i)/tspace+twindow)energy(i)=energy(i)+abs(indata(i,j));%振幅绝对值属性endendendplot(energy);s_cplot(tempdata);属性提取-FFT分析clearalltempdata=altreadsegy(‘mig.sgy’);%读地震数据indata=tempdata‘;%矩阵转置T=length(indata(1,:));%采样率Tr=length(indata(:,1));%道数[lTr,ltime]=textread(‘layer.txt’,‘%f%f’,‘headerlines’,0);%读取层位时间twindow=5;%设置时窗大小tspace=25;%采样间隔fori=1:Trforj=1:Tif(jltime(i)/tspace-twindow&&jltime(i)/tspace+twindow)tt=j-ltime(i)/tspace+twindow;data(i,tt)=indata(i,j);endendend属性提取-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);7.4生成理论地震模型functionwaveletf=30;tdeta=0.004;N=30;nt=0:tdeta:N*tdeta;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);nRandRange=0.2;ntrand=-nRandRange+rand(1,length(st))*nRandRange*2;%随机数组str=sprintf('随机信号n(t)范围%.1g到%.1g',-nRandRange,nRandRange);%xt=st+ntrand;xt=st;nfft=2^ceil(log(length(xt))/log(2));%产生大于xt的最接近2的n次方的整数xw=fft(xt,nfft);plot(nt,bt);title('地震子波b(t)');axis([00.12-11]);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,length(xt)),xt);title('x(t)');axis([00.12-11]);subplot(2,2,4);plot(abs(xw));title('|X(W)|');7.5相关计算functionXiangGuan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).*cos(2*pi*f2.*ntY);rxy=xcorr(xt,yt);%相关计算ryx=xcorr(yt,xt);r_xy=conv(-xt,yt);%离散序列做卷积运算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('\bf\fontsize{12}{\gamma}_{xy}(\tau)');subplot(3,1,2);plot(linspace(0,MX*tdetaX,length(ryx)),ryx);title('\bf\fontsize{12}{\gamma}