地震资料处理实习报告

1《地震资料数据处理》课程设计总结报告专业班级：地球物理学1302班姓名：学号：成绩：2016年12月31日2目录一、设计内容………………………………………………………………（1）褶积滤波……………………………………………………………（2）快变滤波……………………………………………………………（3）褶积滤波与快变滤波的比较………………………………………（4）设计高通滤波因子…………………………………………………（5）频谱分析……………………………………………………………（6）分析补零对振幅谱的影响…………………………………………（7）线性褶积与循环褶积………………………………………………（8）最小平方反滤波……………………………………………………（9）零相位转换…………………………………………………………（11）静校正………………………………………………………………二、程序…………………………………………………………………………3一、设计内容1、褶积滤波理想低通滤波因子理想带通滤波因子原始信号4低通滤波结果带通滤波结果2、快变滤波原始数据5快变低通滤波结果快变带通滤波结果3．褶积滤波与递归滤波的比较6原始数据零相位褶积滤波结果非零相位褶积滤波结果零相位递归滤波7非零相位递归滤波4．设计高通滤波因子原始数据高通滤波因子8频率域高通滤波因子5．频谱分析正弦函数频谱尖脉冲频谱9地震波频谱6．分析补零对振幅谱的影响正弦函数n=60正弦函数n=128\10非周期波形n=60非周期波形n=6411非周期波形n=1287．线性褶积与圆周褶积线性褶积模型圆周褶积结果12长度与圆周褶积相等的线性褶积8．最小平方反滤波原始反射系数序列求出的反射系数序列9．零相位转换13非零相位子波零相位子波11、地震记录14原始数据选择第一道为参考到静校正结果15二、程序1.褶积滤波#includestdio.h#includemath.h#includeconv.c#definepi3.1415926#defineN100#definedt0.002main(){floatx[100],h[101],h1[101],y_low[200],y_band[200];floatdf;inti,m=100,n=101,l=m+n-1;floatf=70.0;floatf1=10.0;floatf2=80.0;FILE*fp1,*fp2,*fp3,*fp4,*fp5;fp1=fopen(INPUT1.DAT,r);for(i=0;i=N;i++){fscanf(fp1,%f,&x[i]);}fp4=fopen(h_low.dat,w);//低通滤波设计for(i=-50;i=50;i++){if(i==0)h[i+50]=2*pi*f/pi;elseh[i+50]=sin(2*pi*f*i*dt)/(pi*i*dt);fprintf(fp4,%f,h[i+50]);//outputlowpassfilter}fp2=fopen(synthesisdata_lowpass.DAT,w);conv(x,m,h,n,y_low,l);for(i=50;il-50;i++){fprintf(fp2,%f\n,y_low[i]);}//带通滤波器fp5=fopen(h_band.dat,w);for(i=-50;i=50;i++){if(i==0)h1[i+50]=140;16elseh1[i+50]=sin(2*pi*f2*i*dt)/(pi*i*dt)-sin(2*pi*f1*i*dt)/(pi*i*dt);fprintf(fp5,%f\n,h1[i+50]);//outputbandpassfilter}fclose(fp5);conv(x,m,h1,n,y_band,l);fp3=fopen(synthesisdata_bandpass.DAT,w);for(i=50;il-50;i++){fprintf(fp3,%f\n,y_band[i]);}fclose(fp1);fclose(fp2);fclose(fp3);}2.快变滤波#includestdio.h#includemath.h#includestdlib.h#includefft.c#definepi3.1415926main(){double*xr,*xi;float*H;inti,np,nfft,k;floatt,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(nppow(2,k))k=k+1;nfft=pow(2,k);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));//readx(n)fp1=fopen(INPUT1.DAT,r);for(i=0;i100;i++){17fscanf(fp1,%f,&z);xr[i]=z;}fclose(fp1);//补零至128位for(i=100;infft;i++){xr[i]=0.0;}for(i=0;infft;i++){xi[i]=0.0;}//transfertofrequencydoaminfft(xr,xi,k,1);//generatelowpassfilter(zero-phase)for(i=0;infft/2;i++){f=i*df;if(f=fc2&&f=fc1)H[i]=1.0;elseH[i]=0.0;}//滤波器对称for(i=nfft/2;infft;i++){f=i*df;H[i]=H[nfft-i];}//filteringinfrequencydomainfor(i=0;infft;i++){xr[i]=xr[i]*H[i];xi[i]=xi[i]*H[i];}fft(xr,xi,k,-1);fp2=fopen(lowpass2.dat,w);for(i=0;infft;i++){fprintf(fp2,%f\n,xr[i]);}fclose(fp2);18//获取高通截至频率fpar=fopen(bandpass.par,r);fscanf(fpar,%f%f,&fc1,&fc2);fp1=fopen(INPUT1.DAT,r);for(i=0;i100;i++){fscanf(fp1,%f,&z);xr[i]=z;}for(i=100;infft;i++){xr[i]=0.0;}for(i=0;infft;i++){xi[i]=0.0;}//transfertofrequencydoaminfft(xr,xi,k,1);//generatelowpassfilter(zero-phase)for(i=0;i=nfft/2;i++){f=i*df;if(f=fc2&&f=fc1)H[i]=1.0;elseH[i]=0.0;}for(i=nfft/2+1;infft;i++){f=i*df;H[i]=H[nfft-i];}//filteringinfrequencydomainfor(i=0;infft;i++){xr[i]=xr[i]*H[i];xi[i]=xi[i]*H[i];}fft(xr,xi,k,-1);fp3=fopen(bandpass2.dat,w);for(i=0;infft;i++){fprintf(fp3,%f\n,xr[i]);}19}3.褶积滤波与递归滤波褶积滤波#includestdio.h#includemath.h#includestdlib.h#includeconv.c#includefft.c#definePI3.1415926main(){voidconv();floatx[50],h[20],y[69],hreverse[20],hzero[39],yreverse[69];floatdt=0.002;inti,m,n,l,p,q;FILE*fp1,*fp2,*fp3,*fp4,*fp5,*fp6;m=50;n=20;l=m+n-1;//readx(n)fp1=fopen(INPUT3.DAT,r);for(i=0;i50;i++){fscanf(fp1,%f,&x[i]);}fclose(fp1);//readfilterfactorh(n)fp2=fopen(hn.dat,r);fp5=fopen(h_reverse.dat,w);for(i=0;i20;i++){fscanf(fp2,%f,&h[i]);hreverse[i]=h[19-i];fprintf(fp5,%f\n,hreverse[i]);}fclose(fp2);fclose(fp5);conv(x,m,h,n,y,l);//非零相位褶积滤波fp3=fopen(con_filter.dat,w);for(i=0;il;i++){fprintf(fp3,%f\n,y[i]);}fclose(fp3);20p=n+n-1;q=m+p-1;//构造零相位滤波因子conv(h,n,hreverse,n,hzero,p);fp6=fopen(zerophasefilterfactor.dat,w);for(i=0;ip;i++){fprintf(fp6,%f\n,hzero[i]);}fclose(fp6);//零相位滤波conv(x,m,hzero,p,yreverse,q);fp4=fopen(convfilterreverse.dat,w);for(i=0;iq;i++){fprintf(fp4,%f\n,yreverse[i]);}fclose(fp4);}递归滤波#includestdio.h#includemath.h#includestdlib.h#definenp50voidmain(){float*x,*a,*b,*fil1,*fil2;inti;voidrecur1();voidrecur2();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(float));fil2=(float*)malloc(np*sizeof(float));//输入地震数据fp1=fopen(INPUT3.DAT,r);for(i=0;inp;i++)fscanf(fp1,%f,x+i);fclose(fp1);//输入a数组值fp2=fopen(a(n).txt,r);for(i=0;i5;i++)fscanf(fp2,%f,a+i);21fclose(fp2);for(i=5;inp;i++)a[i]=0.0;//输入b数组值fp3=fopen(b(n).txt,r);for(i=0;i5;i++)fscanf(fp3,%f,b+i);fclose(fp3);for(i=5;inp;i++)b[i]=0.0;//正向递归滤波recur1(x,a,b,fil1);fp4=fopen(正向递归结果.DAT,wb);for(i=0;inp;i++){fprintf(fp4,%12.4f\n,fil1[i]);}fclose(fp4);for(i=0;inp;i++){printf(%12