地震资料处理实习报告.docx
《地震资料处理实习报告.docx》由会员分享，可在线阅读，更多相关《地震资料处理实习报告.docx（43页珍藏版）》请在冰点文库上搜索。
地震资料处理实习报告
《地震资料数据处理》课程设计
总结报告
专业班级：
地球物理学1302班
姓名：
学号：
成绩：
2016年12月31日
一、设计内容………………………………………………………………
（1）褶积滤波……………………………………………………………
（2）快变滤波……………………………………………………………
（3）褶积滤波与快变滤波的比较………………………………………
（4）设计高通滤波因子…………………………………………………
（5）频谱分析……………………………………………………………
（6）分析补零对振幅谱的影响…………………………………………
（7）线性褶积与循环褶积………………………………………………
（8）最小平方反滤波……………………………………………………
（9）零相位转换…………………………………………………………
（11）静校正………………………………………………………………
二、程序…………………………………………………………………………
一、设计内容
1、褶积滤波
理想低通滤波因子
理想带通滤波因子
原始信号
低通滤波结果
带通滤波结果
2、快变滤波
原始数据
快变低通滤波结果
快变带通滤波结果
3．褶积滤波与递归滤波的比较
原始数据
零相位褶积滤波结果
非零相位褶积滤波结果
零相位递归滤波
非零相位递归滤波
4．设计高通滤波因子
原始数据
高通滤波因子
频率域高通滤波因子
5．频谱分析
正弦函数频谱
尖脉冲频谱
地震波频谱
6．分析补零对振幅谱的影响
正弦函数n=60
正弦函数n=128
\
非周期波形n=60
非周期波形n=64
非周期波形n=128
7．线性褶积与圆周褶积
线性褶积模型
圆周褶积结果
长度与圆周褶积相等的线性褶积
8．最小平方反滤波
原始反射系数序列
求出的反射系数序列
9．零相位转换
非零相位子波
零相位子波
11、地震记录
原始数据
选择第一道为参考到静校正结果
二、程序
1.褶积滤波
#include"stdio.h"
#include"math.h"
#include"conv.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;
else
h[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;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;
else
h1[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;i{
fprintf（fp3,"%f\n",y_band[i]）;
}
fclose（fp1）;fclose（fp2）;fclose（fp3）;
}
2.快变滤波
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#include"fft.c"
#definepi3.1415926
main（）
{
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（np>pow（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;i<100;i++）
{
fscanf（fp1,"%f",&z）;xr[i]=z;
}
fclose（fp1）;
//补零至128位
for（i=100;i{
xr[i]=0.0;
}
for（i=0;i{
xi[i]=0.0;
}
//transfertofrequencydoaminfft（xr,xi,k,1）;
//generatelowpassfilter（zero-phase）for（i=0;i{
f=i*df;
if（f<=fc2&&f>=fc1）H[i]=1.0;elseH[i]=0.0;
}
//滤波器对称for（i=nfft/2;i{
f=i*df;H[i]=H[nfft-i];
}
//filteringinfrequencydomainfor（i=0;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;i{
fprintf（fp2,"%f\n",xr[i]）;
}
fclose（fp2）;
//获取高通截至频率fpar=fopen（"bandpass.par","r"）;fscanf（fpar,"%f%f",&fc1,&fc2）;fp1=fopen（"INPUT1.DAT","r"）;
for（i=0;i<100;i++）
{
fscanf（fp1,"%f",&z）;xr[i]=z;
}
for（i=100;i{
xr[i]=0.0;
}
for（i=0;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;else
H[i]=0.0;
}
for（i=nfft/2+1;i{
f=i*df;H[i]=H[nfft-i];
}
//filteringinfrequencydomainfor（i=0;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;i{
fprintf（fp3,"%f\n",xr[i]）;
}
}
3.褶积滤波与递归滤波
褶积滤波
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#include"conv.c"
#include"fft.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;i<50;i++）
{
fscanf（fp1,"%f",&x[i]）;
}
fclose（fp1）;
//readfilterfactorh（n）fp2=fopen（"hn.dat","r"）;fp5=fopen（"h_reverse.dat","w"）;for（i=0;i<20;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;i{
fprintf（fp3,"%f\n",y[i]）;
}
fclose（fp3）;
p=n+n-1;q=m+p-1;
//构造零相位滤波因子
conv（h,n,hreverse,n,hzero,p）;fp6=fopen（"zerophasefilterfactor.dat","w"）;for（i=0;i{
fprintf（fp6,"%f\n",hzero[i]）;
}
fclose（fp6）;
//零相位滤波conv（x,m,hzero,p,yreverse,q）;fp4=fopen（"convfilterreverse.dat","w"）;for（i=0;i{
fprintf（fp4,"%f\n",yreverse[i]）;
}
fclose（fp4）;
}
递归滤波
#include
#include
#include
#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;i//输入a数组值fp2=fopen（"a（n）.txt","r"）;for（i=0;i<5;i++）fscanf（fp2,"%f",a+i）;
fclose（fp2）;for（i=5;i//输入b数组值fp3=fopen（"b（n）.txt","r"）;for（i=0;i<5;i++）fscanf（fp3,"%f",b+i）;fclose（fp3）;for（i=5;i//正向递归滤波
recur1（x,a,b,fil1）;fp4=fopen（"正向递归结果.DAT","wb"）;for（i=0;i{
fprintf（fp4,"%12.4f\n",fil1[i]）;
}
fclose（fp4）;for（i=0;i{
printf（"%12.4f\n",fil1[i]）;
}
printf（"\n"）;
//反向递归滤波
recur2（fil1,a,b,fil2）;fp5=fopen（"反向递归结果.DAT","wb"）;for（i=0;i{
fprintf（fp5,"%12.4f\n",fil2[i]）;
}
fclose（fp5）;for（i=0;i{
printf（"%12.4f\n",fil2[i]）;
}
}
voidrecur1（floatx[],floata[],floatb[],floatfil1[]）
{
inti,j,k;
floaty1[np],y2[np];for（i=0;i{y1[i]=0.0;
y2[i]=0.0;
for（j=0;j<=i;j++）y1[i]=y1[i]+a[j]*x[i-j];
if（i==0）y2[i]=0.0;
else{for（k=1;k<=i;k++）y2[i]=y2[i]+b[k]*fil1[i-k];
}
fil1[i]=y1[i]-y2[i];
}
}
voidrecur2（floatfil1[],floata[],floatb[],floatfil2[]）
{
inti,j,k;
floaty3[np],y4[np];
for（i=np-1;i>=0;i--）
{
else{
}
y3[i]=0.0;
y4[i]=0.0;
for（j=0;j<=np-1-i;j++）y3[i]=y3[i]+a[j]*fil1[i+j];if（i==np-1）y4[i]=0.0;
for（k=1;k<=np-1-i;k++）y4[i]=y4[i]+b[k]*fil2[i+k];
}
fil2[i]=y3[i]-y4[i];
}
4.设计高通滤波因子
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#include"fft.c"
#defineN101
#definedt0.004
#definePI3.1415926
main（）
{
double*h,*hi;inti,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）
h[i+50]=1.0/dt-60;
else
}
h[i+50]=-sin（2*PI*30.0*i*dt）/（PI*i*dt）;
for（i=100;i<128;i++）
{
h[i]=0.0;
}
for（i=0;i<128;i++）
{
hi[i]=0.0;
}
fp1=fopen（"timedomain.dat","w"）;for（i=0;i<128;i++）
{
fprintf（fp1,"%f\n",h[i]）;
}
fclose（fp1）;fft（h,hi,k,1）;
fp2=fopen（"frequencydoamin.dat","w"）;for（i=0;i<128;i++）
{
fprintf（fp2,"%f\n",sqrt（h[i]*h[i]+hi[i]*hi[i]））;
}
fclose（fp2）;printf（"itisok!
\n"）;
}
5.分析补零对振幅谱的影响
1、不补零
#include"stdio.h"
#include"math.h"
#include"dft.c"
#defineN60
#definePI3.1415926
#definedt0.004main（）
{
floatx[N],xr[N],xi[N],w[N],wr[N],wi[N],z;inti,k;
FILE*fp,*fp1,*fp2,*fp3;
fp3=fopen（"sin.dat","w"）;for（i=0;i{
x[i]=sin（2.0*PI*（i+1）/30）;
fprintf（fp3,"%f\n",x[i]）;
}
fclose（fp3）;
fp=fopen（"WAVE.DAT","r"）;for（i=0;i{
fscanf（fp,"%f",&z）;w[i]=z;
}
fclose（fp）;
printf（"itisok!
\n"）;dft（x,xr,xi,N）;
dft（w,wr,wi,N）;fp1=fopen（"frequencydomain1_60.dat","w"）;fp2=fopen（"frequencydomain2_60.dat","w"）;for（i=0;i{
fprintf（fp1,"%f\n",xr[i]）;
fprintf（fp2,"%f\n",wr[i]）;
}
fclose（fp1）;fclose（fp2）;
}
2、补到64位
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#include"fft.c"
#defineN60
#definePI3.1415926
#definedt0.004main（）
{
voidfft（）;
double*x1,*xi1,*x2,*xi2;floatz;
inti,k,nfft;
FILE*fp,*fp1,*fp2,*fp3,*fp4;k=log（N）/log
（2）;
if（N>pow（2,k））k=k+1;nfft=pow（2,k）;
x1=（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{
x1[i]=sin（2*PI*（i+1）/30）;
}
fp=fopen（"WAVE.DAT","r"）;for（i=0;i{
fscanf（fp,"%f",&z）;x2[i]=z;
}
for（i=0;i{
xi1[i]=0.0;
xi2[i]=0.0;
}
//补到64位
for（i=N;i<64;i++）
{
x1[i]=0.0;
x2[i]=0.0;
xi1[i]=0.0;
xi2[i]=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{
fprintf（fp1,"%f\n",sqrt（x1[i]*x1[i]+xi1[i]*xi1[i]）*dt）;
fprintf（fp2,"%f\n",sqrt（x2[i]*x2[i]+xi2[i]*xi2[i]）*dt）;
}
}
3、补到128位
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#include"fft.c"
#defineN60
#definePI3.1415926
#definedt0.004main（）
{
voidfft（）;
double*x1,*xi1,*x2,*xi2;floatz;
inti,nfft;
FILE*fp,*fp1,*fp2;nfft=128;
x1=（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{
x1[i]=sin（2*PI*（i+1）/30）;
}
fp=fopen（"WAVE.DAT","r"）;for（i=0;i{
fscanf（fp,"%f",&z）;x2[i]=z;
}
for（i=0;i{
xi1[i]=0.0;
xi2[i]=0.0;
}
//补到128位
for（i=N;i<128;i++）
{
x1[i]=0.0;
x2[i]=0.0;
xi1[i]=0.0;
xi2[i]=0.0;
}
fft（x1,xi1,7,1）;
升级会员 