地震资料数字处理程序

程序

1. 褶积滤波

# include "stdio.h" # include "math.h" # include "conv.c" # define pi 3.1415926 # define N 100 # define dt 0.002 main() {

float x[100], h[101], h1[101], y_low[200], y_band[200]; 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,"%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]); //output lowpass filter } fp2=fopen("synthesisdata_lowpass.DAT","w"); conv(x,m,h,n,y_low,l); for(i=50;i<l-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]); // output bandpass filter } fclose(fp5); conv(x,m,h1,n,y_band,l); fp3=fopen("synthesisdata_bandpass.DAT","w"); for(i=50;i<l-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"

# define pi 3.1415926

main() { 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); 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); xr[i]=z; } fclose(fp1); //补零至128位 for(i=100;i<nfft;i++) { xr[i]=0.0; } for(i=0;i<nfft;i++) { xi[i]=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)H[i]=1.0; else H[i]=0.0; } //滤波器对称 for(i=nfft/2;i<nfft;i++) {

f=i*df; H[i]=H[nfft-i]; } //filtering in frequency domain for(i=0;i<nfft;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<nfft;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<nfft;i++) { xr[i]=0.0; } for(i=0;i<nfft;i++) { xi[i]=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)H[i]=1.0; else H[i]=0.0; }

for(i=nfft/2+1;i<nfft;i++) { f=i*df; H[i]=H[nfft-i]; } //filtering in frequency domain for(i=0;i<nfft;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<nfft;i++) { fprintf(fp3,"%f\n",xr[i]); }

}

3. 褶积滤波与递归滤波 褶积滤波

# include "stdio.h" # include "math.h" # include "stdlib.h" # include "conv.c" # include "fft.c"

# define PI 3.1415926 main() { void conv(); float x[50],h[20],y[69],hreverse[20],hzero[39],yreverse[69]; 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; //read x(n) fp1=fopen("INPUT3.DAT","r"); for(i=0;i<50;i++) { fscanf(fp1,"%f",&x[i]); } 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",&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<l;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<p;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<q;i++) { fprintf(fp4,"%f\n",yreverse[i]); } fclose(fp4); }

递归滤波

#include<stdio.h> #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(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++)a[i]=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++)b[i]=0.0; //正向递归滤波 recur1(x,a,b,fil1); fp4=fopen("正向递归结果.DAT","wb"); for(i=0;i<np;i++) {

fprintf(fp4,"%12.4f\n",fil1[i]); }

fclose(fp4);

for(i=0;i<np;i++) {

printf("%12.4f\n",fil1[i]); }

printf("\n"); //反向递归滤波

recur2(fil1,a,b,fil2); fp5=fopen("反向递归结果.DAT","wb"); for(i=0;i<np;i++) {

fprintf(fp5,"%12.4f\n",fil2[i]); }

fclose(fp5);

for(i=0;i<np;i++) {

printf("%12.4f\n",fil2[i]); } }

void recur1(float x[],float a[],float b[],float fil1[]) { int i,j,k;

float y1[np],y2[np]; for(i=0;i<np;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]; } }

void recur2(float fil1[],float a[],float b[],float fil2[]) { int i,j,k; float y3[np],y4[np]; for(i=np-1;i>=0;i--) { 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; else{ 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" # define N 101 # define dt 0.004

# define PI 3.1415926

main() { 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) 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("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 x[N],xr[N],xi[N],w[N],wr[N],wi[N],z; int i,k; FILE *fp,*fp1,*fp2,*fp3;

fp3=fopen("sin.dat","w"); for(i=0;i<N;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<N;i++) { fscanf(fp,"%f",&z); w[i]=z; } fclose(fp);

printf("it is ok !\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<N;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" # 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);

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<N;i++) { x1[i]=sin(2*PI*(i+1)/30); } fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i++) { fscanf(fp,"%f",&z); x2[i]=z; } for(i=0;i<N;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<nfft;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" # 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=(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++) { x1[i]=sin(2*PI*(i+1)/30); } fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i++) { fscanf(fp,"%f",&z); x2[i]=z; } for(i=0;i<N;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); 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,"%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); } }

6. 线性褶积与循环褶积 # include "stdio.h" # include "math.h" # include "fft.c" # define dt 0.002

# define PI 3.1415926 # define L 101

main()

{ void conv(); void cir_conv(); float x[100],h[101],x1[L],h1[L],y[200],y1[100],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",&x[i]); }

fclose(fp);

for(i=-50;i<=50;i++) { if(i==0)h[i+50]=140.0; else h[i+50]=sin(2*PI*70*i*dt)/(PI*i*dt); } //linear convolution

conv(x,m,h,n,y,l); fp1=fopen("linearconv.dat","w"); for(i=0;i<l;i++) { fprintf(fp1,"%f\n",y[i]); } //circus convolution fp=fopen("INPUT1.DAT","r"); for(i=0;i<100;i++) { fscanf(fp,"%f",&x1[i]); }

fclose(fp); df=1.0/(dt*100); kc=70.0/df; for(i=0;i<=50;i++) { if(i>=0 && i<=kc)h1[i]=1.0; else h1[i]=0.0; } for(i=50;i<=100;i++) { h1[i]=h1[100-i]; } if(L>=100) { for(i=100;i<L;i++){ x1[i]=0.0; h1[i]=0.0; } } cir_conv(x1,h1,y1,L); fp2=fopen("circusconv.dat","w"); for(i=0;i<L;i++) { fprintf(fp2,"%f\n",y1[i]); printf("%f\n",y1[i]); } fclose(fp2); }

//线性褶积子程序

void conv(float x[],int m,float h[],int n,float y[],int l) { int k,i;

for(k=0;k<l;k++) { y[k]=0.0;

for(i=0;i<m;i++) if(k-i>=0&&k-i<=n-1) y[k]=y[k]+x[i]*h[k-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=h[k]; h[k]=h[n-k-1]; h[n-k-1]=temp; } for(k=0;k<n;k++) { temp=h[n-1]; for(j=n-2;j>=0;j--) h[j+1]=h[j]; h[0]=temp; y[k]=0.0; for(j=0;j<n;j++) y[k]=x[j]*h[j]*dt+y[k]; } return; }

7. 最小平方反滤波 # include "stdio.h" # include "math.h" # include "tlvs.c" # define N 200

//bn是地震子波序列 a是反射系数系列 main() {

void conv(); void autocorr(); float bn[60],a[200],x[211],rxx[211],a_cal[270]; double t[60],b[60],d[60]; 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",&a[i]); } fclose(fp); fp1=fopen("bn.dat","r"); for(i=0;i<12;i++) { fscanf(fp1,"%f",&bn[i]); } fclose(fp1); //synthetize seismic records conv(bn,n,a,m,x,l); fp11=fopen("synthseismicdata.dat","w"); for(i=0;i<211;i++) { fprintf(fp11,"%f\n",x[i]); }

autocorr(x,rxx,l);

fp12=fopen("autocorrdata.dat","w"); for(i=0;i<l;i++) { fprintf(fp12,"%f\n",rxx[i]); } for(i=0;i<60;i++) { if(i==0)d[i]=1.0; else d[i]=0.0; t[i]=rxx[i]; } tlvs(t,60,d,b); //输出反子波

fp13=fopen("waveletreverse.dat","w"); for(i=0;i<60;i++) { fprintf(fp13,"%f\n",b[i]);

} conv(x,l,b,60,a_cal,l+60-1);

fp2=fopen("sigma.dat","w"); for(i=0;i<270;i++) { fprintf(fp2,"%f\n",a_cal[i]); } }

//褶积子程序

void conv(float x[],int m,float h[],int n,float y[],int l) { int k,i; for(k=0;k<l;k++) { y[k]=0.0; for(i=0;i<m;i++) if(k-i>=0&&k-i<=n-1) y[k]=y[k]+x[i]*h[k-i]*0.004; } }

//自相关子程序

void autocorr(float x[],float y[],int n) { int i,j; for(i=0;i<n;i++) { y[i]=0.0; for(j=i;j<n;j++) y[i]=x[j]*x[j-i]+y[i]; } return; }

8. 零相位转换 # include "stdio.h" # include "math.h" # include "stdlib.h" # include "fft.c"

# define PI 3.1415926

main()

{ void fft(); double *xr,*xi; float z,w[25]; int i,nfft=32; FILE *fp,*fp1,*fp2; xr=(double*)malloc(nfft*sizeof(double)); xi=(double*)malloc(nfft*sizeof(double)); fp=fopen("wavelet.dat","r"); for(i=0;i<25;i++) { fscanf(fp,"%f",&z); xr[i]=z; } for(i=25;i<32;i++) { xr[i]=0.0; } for(i=0;i<32;i++) { xi[i]=0.0; } fft(xr,xi,5,1); fp1=fopen("amplitude.dat","w"); for(i=0;i<32;i++) { fprintf(fp1,"%f\n",sqrt(xr[i]*xr[i]+xi[i]*xi[i])); } for(i=0;i<32;i++){ xr[i]=sqrt(xr[i]*xr[i]+xi[i]*xi[i]); xi[i]=0.0; } fft(xr,xi,5,-1); fp2=fopen("changedwavelet.dat","w"); for(i=0;i<32;i++) { fprintf(fp2,"%f\n",xr[i]); } }

10. 最小相位转换 //最小相位转换 # include "stdio.h" # include "math.h"

# include "tlvs.c"

main() { void conv(); void autocorr(); double temp,bo[73],ao[25],g[49],gr[49]; float z; double b[25],rbb[25],d[25],y[25],c[49]; int i; FILE *fp,*fp1,*fp2; fp=fopen("wavelet.dat","r"); for(i=0;i<25;i++) { fscanf(fp,"%f",&z); b[i]=z;

//printf("%f\n",b[i]); } autocorr(b,25,b,25,rbb); for(i=0;i<25;i++) { if(i==0)d[i]=1.0; else d[i]=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,"%f\n",y[i]); }

//求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+c[i]*c[i]; printf("%f\n",c[i]); }

bo[0]=1.0/sqrt(temp);

本文来源：https://www.bwwdw.com/article/r4ui.html