地震资料处理实习报告

1、 地震资料数据处理课程设计总结报告专业班级： 地球物理学1302班 姓 名： 学 号： 成 绩： 2016年12月31日目 录一 、 设计内容（1）褶积滤波（2）快变滤波（3）褶积滤波与快变滤波的比较（4）设计高通滤波因子（5）频谱分析（6）分析补零对振幅谱的影响（7）线性褶积与循环褶积（8）最小平方反滤波（9）零相位转换（11）静校正二、程序一、设计内容 1、褶积滤波 理想低通滤波因子 理想带通滤波因子原始信号低通滤波结果带通滤波结果2、快变滤波原始数据快变低通滤波结果快变带通滤波结果3褶积滤波与递归滤波的比较原始数据零相位褶积滤波结果非零相位褶积滤波结果零相位递归滤波非零相位递归滤波4设计

2、高通滤波因子原始数据高通滤波因子频率域高通滤波因子5频谱分析正弦函数频谱尖脉冲频谱地震波频谱6分析补零对振幅谱的影响正弦函数n=60正弦函数n=128非周期波形n=60非周期波形n=64非周期波形n=1287线性褶积与圆周褶积线性褶积模型圆周褶积结果长度与圆周褶积相等的线性褶积8最小平方反滤波原始反射系数序列求出的反射系数序列9零相位转换非零相位子波零相位子波11、地震记录原始数据选择第一道为参考到静校正结果二、程序1. 褶积滤波# include "stdio.h"# include "math.h"# include "conv.c&quo

3、t;# define pi 3.1415926# define N 100# define dt 0.002 main()float x100, h101, h1101, y_low200,y_band200; float df;int i,m=100,n=101,l=m+n-1; float f=70.0;float f1=10.0; float f2=80.0;FILE *fp1,*fp2,*fp3,*fp4,*fp5; fp1=fopen("INPUT1.DAT","r");for(i=0;i<=N;i+)fscanf(fp1,"%

4、f",&xi);fp4=fopen("h_low.dat","w");/低通滤波设计 for(i=-50;i<=50;i+)if(i=0)hi+50=2*pi*f/pi;elsehi+50=sin(2*pi*f*i*dt)/(pi*i*dt);fprintf(fp4,"%f",hi+50);/output lowpass filterfp2=fopen("synthesisdata_lowpass.DAT","w"); conv(x,m,h,n,y_low,l);for(

5、i=50;i<l-50;i+)fprintf(fp2,"%fn",y_lowi);/带通滤波器 fp5=fopen("h_band.dat","w"); for(i=-50;i<=50;i+)if(i=0)h1i+50=140;23elseh1i+50=sin(2*pi*f2*i*dt)/(pi*i*dt)-sin(2*pi*f1*i*dt)/(pi*i*dt);fprintf(fp5,"%fn",h1i+50); / output bandpass filterfclose(fp5); conv(x,

6、m,h1,n,y_band,l);fp3=fopen("synthesisdata_bandpass.DAT","w"); for(i=50;i<l-50;i+)fprintf(fp3,"%fn",y_bandi);fclose(fp1);fclose(fp2);fclose(fp3);2. 快变滤波# include "stdio.h"# include "math.h"# include "stdlib.h"# include "fft.c"#

7、define pi 3.1415926main()double *xr,*xi; float *H;int i,np,nfft,k;float t,dt,df,f,z,fc1,fc2; FILE *fpar,*fp1,*fp2,*fp3;/从参数文件中获得截至频率fpar=fopen("lowpassfilter.par","r"); fscanf(fpar,"%f%f",&fc1,&fc2); np=100;k=log(np)/log(2); if(np>pow(2,k)k=k+1; nfft=pow(2,k)

8、; dt=0.002;df=1.0/(nfft*dt); xr=(double*)calloc(nfft,sizeof(double); xi=(double*)calloc(nfft,sizeof(double); H=(float*)calloc(nfft,sizeof(float);/ read x(n) fp1=fopen("INPUT1.DAT","r");for(i=0;i<100;i+)fscanf(fp1,"%f",&z); xri=z;fclose(fp1);/补零至 128 位for(i=100;i&

9、lt;nfft;i+)xri=0.0;for(i=0;i<nfft;i+)xii=0.0;/transfer to frequency doamin fft(xr,xi,k,1);/generate lowpass filter(zero-phase) for(i=0;i<nfft/2;i+)f=i*df;if(f<=fc2 && f>=fc1)Hi=1.0; else Hi=0.0;/滤波器对称 for(i=nfft/2;i<nfft;i+)f=i*df; Hi=Hnfft-i;/filtering in frequency domain for

10、(i=0;i<nfft;i+)xri=xri*Hi;xii=xii*Hi;fft(xr,xi,k,-1); fp2=fopen("lowpass2.dat","w"); for(i=0;i<nfft;i+)fprintf(fp2,"%fn",xri);fclose(fp2);/获取高通截至频率 fpar=fopen("bandpass.par","r"); fscanf(fpar,"%f%f",&fc1,&fc2); fp1=fopen("

11、;INPUT1.DAT","r");for(i=0;i<100;i+)fscanf(fp1,"%f",&z); xri=z;for(i=100;i<nfft;i+)xri=0.0;for(i=0;i<nfft;i+)xii=0.0;/transfer to frequency doamin fft(xr,xi,k,1);/generate lowpass filter(zero-phase) for(i=0;i<=nfft/2;i+)f=i*df;if(f<=fc2 && f>=fc1

12、)Hi=1.0; elseHi=0.0;for(i=nfft/2+1;i<nfft;i+)f=i*df; Hi=Hnfft-i;/filtering in frequency domain for(i=0;i<nfft;i+)xri=xri*Hi;xii=xii*Hi;fft(xr,xi,k,-1); fp3=fopen("bandpass2.dat","w"); for(i=0;i<nfft;i+)fprintf(fp3,"%fn",xri);3. 褶积滤波与递归滤波褶积滤波# include "stdi

13、o.h"# include "math.h"# include "stdlib.h"# include "conv.c"# include "fft.c"# define PI 3.1415926 main()void conv();float x50,h20,y69,hreverse20,hzero39,yreverse69; float dt=0.002;int i,m,n,l,p,q;FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6; m=50;n=20;l=m+n-1;/rea

14、d x(n) fp1=fopen("INPUT3.DAT","r");for(i=0;i<50;i+)fscanf(fp1,"%f",&xi);fclose(fp1);/read filterfactor h(n) fp2=fopen("hn.dat","r"); fp5=fopen("h_reverse.dat","w"); for(i=0;i<20;i+)fscanf(fp2,"%f",&hi); hre

15、versei=h19-i; fprintf(fp5,"%fn",hreversei);fclose(fp2);fclose(fp5);conv(x,m,h,n,y,l);/非零相位褶积滤波 fp3=fopen("con_filter.dat","w"); for(i=0;i<l;i+)fprintf(fp3,"%fn",yi);fclose(fp3);p=n+n-1; q=m+p-1;/构造零相位滤波因子conv(h,n,hreverse,n,hzero,p); fp6=fopen("zeropha

16、sefilterfactor.dat","w"); for(i=0;i<p;i+)fprintf(fp6,"%fn",hzeroi);fclose(fp6);/零相位滤波 conv(x,m,hzero,p,yreverse,q); fp4=fopen("convfilterreverse.dat","w"); for(i=0;i<q;i+)fprintf(fp4,"%fn",yreversei);fclose(fp4);递归滤波#include<stdio.h>

17、#include<math.h>#include<stdlib.h>#define np 50 void main()float *x,*a,*b,*fil1,*fil2; int i;void recur1(); void recur2();FILE *fp1,*fp2,*fp3,*fp4,*fp5; x=(float*)malloc(np*sizeof(float); a=(float*)malloc(np*sizeof(float); b=(float*)malloc(np*sizeof(float); fil1=(float*)malloc(np*sizeof(

18、float); fil2=(float*)malloc(np*sizeof(float);/输入地震数据fp1=fopen("INPUT3.DAT","r");for(i=0;i<np;i+)fscanf(fp1,"%f",x+i); fclose(fp1);/输入 a 数组值 fp2=fopen("a(n).txt","r"); for(i=0;i<5;i+)fscanf(fp2,"%f",a+i);fclose(fp2); for(i=5;i<np;i+

19、)ai=0.0;/输入 b 数组值 fp3=fopen("b(n).txt","r"); for(i=0;i<5;i+)fscanf(fp3,"%f",b+i); fclose(fp3); for(i=5;i<np;i+)bi=0.0;/正向递归滤波recur1(x,a,b,fil1); fp4=fopen("正向递归结果.DAT","wb"); for(i=0;i<np;i+)fprintf(fp4,"%12.4fn",fil1i);fclose(fp4)

20、; for(i=0;i<np;i+)printf("%12.4fn",fil1i);printf("n");/反向递归滤波recur2(fil1,a,b,fil2); fp5=fopen("反向递归结果.DAT","wb"); for(i=0;i<np;i+)fprintf(fp5,"%12.4fn",fil2i);fclose(fp5); for(i=0;i<np;i+)printf("%12.4fn",fil2i);void recur1(float x

21、,float a,float b,float fil1)int i,j,k;float y1np,y2np; for(i=0;i<np;i+) y1i=0.0;y2i=0.0;for(j=0;j<=i;j+)y1i=y1i+aj*xi-j;if(i=0)y2i=0.0;else for(k=1;k<=i;k+)y2i=y2i+bk*fil1i-k;fil1i=y1i-y2i;void recur2(float fil1,float a,float b,float fil2)int i,j,k;float y3np,y4np;for(i=np-1;i>=0;i-)else

22、y3i=0.0;y4i=0.0;for(j=0;j<=np-1-i;j+)y3i=y3i+aj*fil1i+j; if(i=np-1)y4i=0.0;for(k=1;k<=np-1-i;k+)y4i=y4i+bk*fil2i+k;fil2i=y3i-y4i;4. 设计高通滤波因子# include "stdio.h"# include "math.h"# include "stdlib.h"# include "fft.c"# define N 101# define dt 0.004# define

23、PI 3.1415926main()double *h,*hi; int i,k,nfft; FILE *fp1,*fp2;k=log(N)/log(2); if(N>pow(2,k)k=k+1; nfft=pow(2,k);h=(double*)malloc(nfft*sizeof(double); hi=(double*)malloc(nfft*sizeof(double); for(i=-50;i<=50;i+)if(i=0)hi+50=1.0/dt-60;elsehi+50=-sin(2*PI*30.0*i*dt)/(PI*i*dt);for(i=100;i<128;

24、i+)hi=0.0;for(i=0;i<128;i+)hii=0.0;fp1=fopen("timedomain.dat","w"); for(i=0;i<128;i+)fprintf(fp1,"%fn",hi);fclose(fp1); fft(h,hi,k,1);fp2=fopen("frequencydoamin.dat","w"); for(i=0;i<128;i+)fprintf(fp2,"%fn",sqrt(hi*hi+hii*hii);fclo

25、se(fp2); printf("it is ok!n");5. 分析补零对振幅谱的影响1、 不补零# include "stdio.h"# include "math.h"# include "dft.c"# define N 60# define PI 3.1415926# define dt 0.004 main()float xN,xrN,xiN,wN,wrN,wiN,z; int i,k;FILE *fp,*fp1,*fp2,*fp3;fp3=fopen("sin.dat","

26、;w"); for(i=0;i<N;i+)xi=sin(2.0*PI*(i+1)/30);fprintf(fp3,"%fn",xi);fclose(fp3);fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i+)fscanf(fp,"%f",&z); wi=z;fclose(fp);printf("it is ok !n"); dft(x,xr,xi,N);dft(w,wr,wi,N); fp1=fopen("frequency

27、domain1_60.dat","w"); fp2=fopen("frequencydomain2_60.dat","w"); for(i=0;i<N;i+)fprintf(fp1,"%fn",xri);fprintf(fp2,"%fn",wri);fclose(fp1);fclose(fp2);2、 补到 64 位# include "stdio.h"# include "math.h"# include "stdlib.h&q

28、uot;# include "fft.c"# define N 60# define PI 3.1415926# define dt 0.004 main()void fft();double *x1,*xi1,*x2,*xi2; float z;int i,k,nfft;FILE *fp,*fp1,*fp2,*fp3,*fp4; k=log(N)/log(2);第 19 页if(N>pow(2,k)k=k+1; nfft=pow(2,k);x1=(double*)malloc(nfft*sizeof(double); xi1=(double*)malloc(nfft

29、*sizeof(double); x2=(double*)malloc(nfft*sizeof(double)xi2=(double*)malloc(nfft*sizeof(double); for(i=0;i<N;i+)x1i=sin(2*PI*(i+1)/30);fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i+)fscanf(fp,"%f",&z); x2i=z;for(i=0;i<N;i+)xi1i=0.0;xi2i=0.0;/补到 64 位for(i=N;i<64

30、;i+)x1i=0.0;x2i=0.0;xi1i=0.0;xi2i=0.0;fft(x1,xi1,6,1);fft(x2,xi2,6,1); fp1=fopen("frequencydomain1_64.dat","w"); fp2=fopen("frequencydomain2_64.dat","w"); for(i=0;i<nfft;i+)fprintf(fp1,"%fn",sqrt(x1i*x1i+xi1i*xi1i)*dt);fprintf(fp2,"%fn",

31、sqrt(x2i*x2i+xi2i*xi2i)*dt);3、 补到 128 位# include "stdio.h"# include "math.h"# include "stdlib.h"# include "fft.c"# define N 60# define PI 3.1415926# define dt 0.004 main()void fft();double *x1,*xi1,*x2,*xi2; float z;int i,nfft;FILE *fp,*fp1,*fp2; nfft=128;x1=(

32、double*)malloc(nfft*sizeof(double); xi1=(double*)malloc(nfft*sizeof(double); x2=(double*)malloc(nfft*sizeof(double); xi2=(double*)malloc(nfft*sizeof(double);for(i=0;i<N;i+)x1i=sin(2*PI*(i+1)/30);fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i+)fscanf(fp,"%f",&z); x2i=

33、z;for(i=0;i<N;i+)xi1i=0.0;xi2i=0.0;/补到 128 位for(i=N;i<128;i+)x1i=0.0;x2i=0.0;xi1i=0.0;xi2i=0.0;fft(x1,xi1,7,1);fft(x2,xi2,7,1); fp1=fopen("frequencydomain1_128.dat","w"); fp2=fopen("frequencydomain2_128.dat","w"); for(i=0;i<nfft;i+)fprintf(fp1,"%

34、fn",sqrt(x1i*x1i+xi1i*xi1i)*dt);fprintf(fp2,"%fn",sqrt(x2i*x2i+xi2i*xi2i)*dt);6. 线性褶积与循环褶积# include "stdio.h"# include "math.h"# include "fft.c"# define dt 0.002# define PI 3.1415926# define L 101main()void conv(); void cir_conv();float x100,h101,x1L,h1L,

35、y200,y1100,df;int i,m,n,l,kc;FILE*fp,*fp1,*fp2; m=100;n=101;l=m+n-1; fp=fopen("INPUT1.DAT","r");for(i=0;i<100;i+)fscanf(fp,"%f",&xi);fclose(fp); for(i=-50;i<=50;i+)if(i=0)hi+50=140.0;elsehi+50=sin(2*PI*70*i*dt)/(PI*i*dt);/linear convolutionconv(x,m,h,n,y,l); f

36、p1=fopen("linearconv.dat","w"); for(i=0;i<l;i+)fprintf(fp1,"%fn",yi);/circus convolution fp=fopen("INPUT1.DAT","r");for(i=0;i<100;i+)fscanf(fp,"%f",&x1i);fclose(fp);df=1.0/(dt*100); kc=70.0/df; for(i=0;i<=50;i+)if(i>=0 &

37、& i<=kc)h1i=1.0; elseh1i=0.0;for(i=50;i<=100;i+)h1i=h1100-i;if(L>=100)for(i=100;i<L;i+) x1i=0.0;h1i=0.0;cir_conv(x1,h1,y1,L); fp2=fopen("circusconv.dat","w"); for(i=0;i<L;i+)fprintf(fp2,"%fn",y1i);printf("%fn",y1i);fclose(fp2);/线性褶积子程序void c

38、onv(float x,int m,float h,int n,float y,int l)int k,i; for(k=0;k<l;k+)yk=0.0;for(i=0;i<m;i+)if(k-i>=0&&k-i<=n-1) yk=yk+xi*hk-i*dt;/圆周褶积子程序void cir_conv(float x,float h,float y,int n)int k,j; int temp;for(k=0;k<n/2;k+)temp=hk; hk=hn-k-1;hn-k-1=temp;for(k=0;k<n;k+)temp=hn-1;

39、for(j=n-2;j>=0;j-)hj+1=hj; h0=temp; yk=0.0;for(j=0;j<n;j+)yk=xj*hj*dt+yk;return;7. 最小平方反滤波# include "stdio.h"# include "math.h"# include "tlvs.c"# define N 200/bn 是地震子波序列 a 是反射系数系列main()void conv(); void autocorr();float bn60,a200,x211,rxx211,a_cal270; double t60,

40、b60,d60;int i,m,n,l;FILE *fp,*fp1,*fp2,*fp11,*fp12,*fp13; n=12;m=N;l=m+n-1; fp=fopen("INPUT8.DAT","r");for(i=0;i<200;i+)fscanf(fp,"%f",&ai);fclose(fp); fp1=fopen("bn.dat","r"); for(i=0;i<12;i+)fscanf(fp1,"%f",&bni);fclose(fp1)

41、;/synthetize seismic records conv(bn,n,a,m,x,l); fp11=fopen("synthseismicdata.dat","w"); for(i=0;i<211;i+)fprintf(fp11,"%fn",xi);autocorr(x,rxx,l); fp12=fopen("autocorrdata.dat","w"); for(i=0;i<l;i+)fprintf(fp12,"%fn",rxxi);for(i=0;i&

42、lt;60;i+) if(i=0)di=1.0;elsedi=0.0;ti=rxxi;tlvs(t,60,d,b);/输出反子波 fp13=fopen("waveletreverse.dat","w"); for(i=0;i<60;i+)fprintf(fp13,"%fn",bi);conv(x,l,b,60,a_cal,l+60-1);fp2=fopen("sigma.dat","w"); for(i=0;i<270;i+)fprintf(fp2,"%fn",a

43、_cali);/褶积子程序void conv(float x,int m,float h,int n,float y,int l)int k,i; for(k=0;k<l;k+)yk=0.0;for(i=0;i<m;i+)if(k-i>=0&&k-i<=n-1) yk=yk+xi*hk-i*0.004;/自相关子程序void autocorr(float x,float y,int n)int i,j; for(i=0;i<n;i+)yi=0.0;for(j=i;j<n;j+)yi=xj*xj-i+yi;return;8. 零相位转换# in

44、clude "stdio.h"# include "math.h"# include "stdlib.h"# include "fft.c"# define PI 3.1415926main()void fft(); double *xr,*xi; float z,w25; int i,nfft=32;FILE *fp,*fp1,*fp2;xr=(double*)malloc(nfft*sizeof(double); xi=(double*)malloc(nfft*sizeof(double); fp=fopen(

45、"wavelet.dat","r"); for(i=0;i<25;i+)fscanf(fp,"%f",&z); xri=z;for(i=25;i<32;i+)xri=0.0;for(i=0;i<32;i+)xii=0.0;fft(xr,xi,5,1); fp1=fopen("amplitude.dat","w"); for(i=0;i<32;i+)fprintf(fp1,"%fn",sqrt(xri*xri+xii*xii);for(i=0;i

46、<32;i+) xri=sqrt(xri*xri+xii*xii); xii=0.0;fft(xr,xi,5,-1); fp2=fopen("changedwavelet.dat","w"); for(i=0;i<32;i+)fprintf(fp2,"%fn",xri);10. 最小相位转换/最小相位转换# include "stdio.h"# include "math.h"# include "tlvs.c"main()void conv(); void au

47、tocorr();double temp,bo73,ao25,g49,gr49; float z;double b25,rbb25,d25,y25,c49; int i;FILE *fp,*fp1,*fp2; fp=fopen("wavelet.dat","r"); for(i=0;i<25;i+)fscanf(fp,"%f",&z); bi=z;/printf("%fn",bi);autocorr(b,25,b,25,rbb); for(i=0;i<25;i+)if(i=0)di=1.0;el

48、sedi=0.0;/求 y(t) tlvs(rbb,25,d,y);/printf("it is ok !n"); fp1=fopen("yt.dat","w");for(i=0;i<25;i+)fprintf(fp1,"%fn",yi);/求 bo(0) conv(b,25,y,25,c,49);/c(t)=b(t)*g(t) temp=0.0;for(i=0;i<49;i+)temp=temp+ci*ci; printf("%fn",ci);bo0=1.0/sqrt(temp);

49、printf("bo(0)=%fn",bo0);/求 ao(t) for(i=0;i<25;i+)aoi=bo0*yi;conv(b,25,ao,25,g,49);/g(t)=b(t)*ao(t)/求 g(-t) for(i=0;i<49;i+)gri=g48-i;/求 bo(t)=b(t)*g(-t)conv(b,25,gr,49,bo,73); fp2=fopen("minphasew.dat","w"); for(i=0;i<73;i+)fprintf(fp2,"%fn",boi);/褶积子

50、程序void conv(double x,int m,double h,int n,double y,int l)int k,i; for(k=0;k<l;k+)yk=0.0;for(i=0;i<m;i+)if(k-i>=0&&k-i<=n-1) yk=yk+xi*hk-i*0.004;/自相关子程序void autocorr(double x,int m,double h,int n,double y )int i,j; for(i=0;i<m;i+)yi=0.0;for(j=1;j<n;j+)if(i+j<m) yi=xj*hj+i

51、+yi;第 25 页return;11. 静校正处理# include "stdio.h"# include "math.h"main()void ccor();float x10100,y100,rxy10199;float z,max;int i,j,k,tao10; FILE *fp,*fp1,*fp2;/从文件中读取地震数据fp=fopen("seis.dat","r"); for(i=0;i<10;i+)for(j=0;j<100;j+)fscanf(fp,"%f",&am

52、p;z);xij=z;fclose(fp);/求参考道 for(j=0;j<100;j+) yj=0.0;for(i=0;i<10;i+)yj=yj+xij;yj=0.1*yj;/printf("%fn",yj);/互相关求取静校正量值 for(i=0;i<10;i+)ccor(xi,100,y,100,rxyi);fp1=fopen("rxy.dat","w"); for(i=0;i<10;i+)for(j=0;j<199;j+) fprintf(fp1,"%fn",rxyij);f

53、close(fp1);/求取 tao 值 for(i=0;i<10;i+)max=rxyi0; taoi=0; for(j=0;j<199;j+)if(rxyij>max)max=rxyij; taoi=j;taoi=taoi-99;printf("%dn",taoi);/静校正处理 for(i=0;i<10;i+)if(taoi>=0) for(j=0;j<100;j+)k=j+taoi; if(k<100)xij=xik;if(taoi<0)for(j=99;j>=0;j-)k=j+taoi; if(k>=0)

54、xij=xik;fp2=fopen("processeddata.DAt","w"); for(i=0;i<10;i+)for(j=0;j<100;j+) fprintf(fp2,"%fn",xij);fclose(fp2);/自相关程序void ccor(float x,int m,float h,int n,float y)int i,j;for(i=-n+1;i<=m-1;i+)yi+n-1=0.0;for(j=0;j<=m-1;j+)if(j-i>=0&&j-i<=n-1)y

55、i+n-1=yi+n-1+xj*hj-i;部分子程序：/程序名：fft.c#include "stdio.h"#include "math.h"#include "stdlib.h"#define PI 3.1415926/*sr,sx：双精度型一维数组,输入（输出）信号的实部和虚部 */*m0: 2 的次方数, 2*m0=nfft */*inv=1 forward transform; inv=-1 inverse transform*/ void fft(double sr,double sx,int m0,int inv)int i,j,lm,li,k,lmx,np,lix,mm2; double t1,t2,c,s,cv,st,ct;if(m0<0)return; lmx=1;for(i=1;i<=m0;+i)lmx+=lmx;/form 2*