油气勘探方法程序设计课程设计(论文)-最小平方反滤波

1、油气勘探方法程序设计课程设计(论文)设计（论文）题目 最小平方反滤波 学院名称 地球物理学院 专业名称 勘查技术与工程（石油物探） 学生姓名 学生学号 任课教师 田仁飞 设计（论文）成绩 教务处 制2015年7月1日最小平方反滤波一、方法原理由地震波的传播理论可知，在粘弹性介质中地震波是以地震子波的形式在地下传播，地面接收到的反射波地震记录是地层反射系数与地震子波的褶积，因此，地层相当于一个滤波器，是反射系数序列变成了由子波组成的地震记录，降低了地震勘探的纵向分辨率。反滤波的目的就是要设计一个反滤波器，来对地震记录滤波，消除地层滤波的作用，提高地震记录的纵向分辨率。最小平方反滤波的基本思想在于

2、设计一个滤波算子，用它把已知的输入信号转换为与给定的期望输出信号在最小平方误差的意义下是最佳输出。1）用最小平方法求反滤波因子对输入子波b(t)反滤波后的期望输出为d(t),实际输出为y(t)，按最小平方原理，使二者的误差平方和Q为最小时求得的反滤波因子称为最小平方反滤波因子，用它对地震记录x(t)进行的反滤波为最小平方反滤波。设输入离散信号为地震子波b(n)=b(1),b(2), ,b(m),待求的反滤波因子a(n)=a(0),a(1),a(2), ,a(m), a(t)的起始时间为0，（m+1）为a(t)的延续长度，b(n)与a(n)的褶积为实际输出y(n),即 yn=an*bn=0mab

3、n- 1实际输出与期望输出的误差平方和为Q=n=0myn-d(n)2 =n=0m=0mabn-d(n)2 （2）要使Q为最小，数学上就是求Q的极值问题，即求满足 Qa(l)=0 l=0,1,m （3）的滤波因子a(t)。即满足方程 =0marbbl-=rbdl l=0,1,m （4）其中rbbl-=n=0mbn-b(n-l)为地震子波的自相关函数，rbdl=n=0md(n)b(n-l),写成矩阵形式 rbb(0)rbb(1)rbb(m)rbb(1)rbb(m)rbb(0)rbb(m-1)rbb(m-1)rbb(0)a0a1am=adb0adb1adbm （5）上式系数矩阵称为托布利兹方程。当期

4、望输出是单位脉冲t时，即 dt=t=1,t=00,t0 （6）则 rbdl=b-l （7） 由于b(t)是地震子波，b(t)=0,当t<0时，上述矩阵可写成： rbb(0)rbb(1)rbb(m)rbb(1)rbb(m)rbb(0)rbb(m-1)rbb(m-1)rbb(0)a0a1am=b000=100 （8）用地震记录自相关系数rxx()代替ybb()2）已知地震记录x(t) x(t)=b(t)* (t) （9）其中x(t)为地震记录，(t)是反射系数。根据自相关函数的性质，可证明：地震记录自相关rxx(n)是地震子波自相关rbb(n)和反射系数自相关r(n)的褶积,即： rxxn=

5、rbbn*rn （10）假定反射系数是随机白噪音，则 rn=N0,n=00,n0 （11）故可得： rxxn=t=0+rbbtr(n-t)=N0rbbn （12）故得到 rxx(1)rxx(m)rxx(0)rxx(m-1)rxx(m-1)rxx(0)a0a1am=100 （13）其中 ai=1N0b0ai （14）解出方程，即可求得反滤波因子a(t)。二、方法步骤在地震记录上选取自相关段，要选择在波形单一，能量强的地方。自相关段的长度应大于或等于子波的长度；求选取段的自相关函数： rxxn=t=1m-nxt×xt+n, n=0,1,2,M-1 （15） 其中M为自相关长度；把rxxn

6、代入式中求解at，并由式得反滤波因子a(t)，式中，通常取N0=1；将地震记录x(t)与反滤波因子a(t)褶积，完成反滤波； yt=0m-1a×xt-, t=m,m+1,L （16）其中m为反滤波因子长度，L为地震记录长度。三、源程序1）雷克子波模型源程序#include "stdio.h"#include "math.h" #include "stdlib.h" #include "malloc.h" #include "string.h" #define PI 3.1415926

7、#define fm 25 /主频 #define dt 0.001 /采样间隔 #define dz 10 /深度采样间隔#define N 100 /反射系数序列长度#define Nw 100 /子波长度#define P 199 /合成地震记录长度/=建立速度模型=/void velocity(float v,int n)int i,z;float h=1000;for(i=0;i<n;i+)z=i*dz;if(z>=0)&&(z<h/5)vi=4000; if(z>=h/5)&&(z<h*2/5)vi=5000;if(z&

8、gt;=h*2/5)&&(z<h*3/5)vi=4700;if(z>=h*3/5)&&(z<h*4/5)vi=5600;if(z>=h*4/5)&&(z<h)vi=6700;/=形成反射系数序列=/void reflection(float r,float v,int n)int i;for(i=0;i<n-1;i+)ri=(float)(vi+1-vi)/(vi+1+vi);/=形成地震子波=/void wave(float w,int n)int i;double t;for(i=0;i<n;i+)t

9、=i*dt; wi=exp(-2*fm*fm*t*t*log(2)*sin(2*PI*fm*t);/=求褶积函数=/void convolution(float w,float r,float c,int M,int L) int i,j,k; float tp=0; k=M+L-1; for(i=0;i<k;i+) for(j=0;j<M;j+) if(i-j)>=0&&(i-j)<L) tp+=wj*r(i-j); ci=tp; tp=0; /=自相关函数=/void autocorrelation(float *a,float *r,int n)

10、int i,j; float s=0; for(i=0;i<n;i+) for(j=0;j<=i;j+) s=s+aj*an-1-i+j;/从最左边开始，移到两者重合 rn-1-i=s; s=0; /=莱文森递推解托普利茨方程组=/ /t为矩阵的第一行或者第一列数据，b为方程组右侧数据，x为计算的反滤波因子int seq_Toeplitz(float t,int n,float b,float x) int i,j,k; float a,beta,q,c,h,*y,*s; s=(float*)calloc(n,sizeof(double); y=(float*)calloc(n,s

11、izeof(double); a=t0; if (fabs(a)+1.0=1.0) free(s); free(y); printf("failn"); return(-1); y0=1.0;x0=b0/a; for(k=1;k<=n-1;k+) beta=0.0;q=0.0; for(j=0;j<=k-1;j+) beta=beta+tj+1*yj; q=q+tk-j*xj; if(fabs(a)+1.0=1.0) free(s);free(y);printf("failn");return(-1); c=-beta/a;s0=c*yk-1

12、;yk=yk-1;if(k!=1) for(i=1;i<=k-1;i+)si=yi-1+c*yk-i-1; a=a+c*beta; if(fabs(a)+1.0=1.0) free(s);free(y);printf("failn");return(-1); h=(bk-q)/a; for(i=0;i<=k-1;i+) xi=xi+h*si;yi=si; xk=h*yk; free(s); free(y); return(1); void main() int i; float vN,rN=0.,wNw,c199;float rww100,b100,a100,y

13、c298; FILE *fpv,*fpr,*fpw,*fpc,*fpaw,*fpaf,*fps; /=建立速度模型=/fpv=fopen("velocity.dat","wb");if(fpv=NULL)printf("can't open velocity filen");exit(0); velocity(v,N);for(i=0;i<N;i+)fwrite(&vi,sizeof(vi),1,fpv); /=形成反射系数序列=/fpr=fopen("reflection.csv",&quo

14、t;w");if(fpr=NULL)printf("can't open reflection filen");exit(0);reflection(r,v,N);for(i=0;i<N;i+)fprintf(fpr,"%fn",ri); /=形成子波=/fpw=fopen("wave.csv","w");if(fpw=NULL)printf("can't open wave filen");exit(0);wave(w,Nw);for(i=0;i<Nw;i

15、+) fprintf(fpw,"%fn",wi); /=褶积合成地震记录=/fpc=fopen("convolution.csv","w");if(fpc=NULL)printf("can't open convolution filen");exit(0);convolution(w,r,c,N,Nw); for(i=0;i<P;i+) fprintf(fpc,"%fn",ci); /=子波自相关=/fpaw=fopen("autocorrelation.csv&quo

16、t;,"w");if(fpaw=NULL)printf("can't open autocorrelation filen");exit(0);autocorrelation(w,rww,Nw); for(i=0;i<Nw;i+) fprintf(fpaw,"%fn",rwwi); /=形成方程组右侧数据=/for(i=0;i<Nw;i+) bi=0; b0=1; /=解Toeplitz方程组，求反滤波因子=/fpaf=fopen("antifilter.csv","w");

17、if(fpaf=NULL)printf("can't open antifilter filen");exit(0);seq_Toeplitz(rww,Nw,b,a); for(i=0;i<Nw;i+) fprintf(fpaf,"%fn",ai); /=反滤波因子与合成地震记录褶积=/ fps=fopen("deconvolution.csv","w");if(fps=NULL)printf("can't open deconvolution filen");exit(0

18、);convolution(c,a,yc,P,N); for(i=0;i<N;i+) fprintf(fps,"%fn",yci); fclose(fpv);fclose(fpr);fclose(fpw);fclose(fpc);fclose(fpaw);fclose(fpaf);fclose(fps);2) 地震记录最小平方反滤波源程序#include <stdio.h>#include <math.h>#include <stdlib.h> #include <malloc.h>#include <string

19、.h>#include <ctype.h>#define LEN 3200#define NQ2 30000#define Ns 1251 /=自相关函数=/void autocorrelation(float *a,float *r,int n) int i,j; float s=0; for(i=0;i<n;i+) for(j=0;j<=i;j+) s=s+aj*an-1-i+j;/从最左边开始，移到两者重合 rn-1-i=s; s=0; /=求褶积函数=/void convolution(float w,float r,float c,int M,int L

20、) int i,j,k; float tp=0; k=M+L-1; for(i=0;i<k;i+) for(j=0;j<M;j+) if(i-j)>=0&&(i-j)<L) tp+=wj*r(i-j); ci=tp; tp=0; /=莱文森递推解托普利茨方程组=/ /t为矩阵的第一行或者第一列数据，b为方程组右侧数据，x为计算的反滤波因子int seq_Toeplitz(float t,int n,float b,float x) int i,j,k; float a,beta,q,c,h,*y,*s; s=(float*)calloc(n,sizeof

21、(double); y=(float*)calloc(n,sizeof(double); a=t0; if (fabs(a)+1.0=1.0) free(s); free(y); printf("failn"); return(-1); y0=1.0;x0=b0/a; for(k=1;k<=n-1;k+) beta=0.0;q=0.0; for(j=0;j<=k-1;j+) beta=beta+tj+1*yj; q=q+tk-j*xj; if(fabs(a)+1.0=1.0) free(s);free(y);printf("failn");r

22、eturn(-1); c=-beta/a;s0=c*yk-1;yk=yk-1;if(k!=1) for(i=1;i<=k-1;i+)si=yi-1+c*yk-i-1; a=a+c*beta; if(fabs(a)+1.0=1.0) free(s);free(y);printf("failn");return(-1); h=(bk-q)/a; for(i=0;i<=k-1;i+) xi=xi+h*si;yi=si; xk=h*yk; free(s); free(y); return(1); void main () int nflg,n,l,i;int len,l

23、en1,jhd;float kiNQ2,koNQ2;float xNs,rssNs,bNs,aNs,r2*Ns-1;char namei50,nameo50; FILE *fp,*fpr,*fp1,*fp2;/ union int hdNQ2; float dataNQ2; char bufNQ2; short int shdNQ2; dat; printf("input file name: "); scanf("%s",&namei); printf("output file name: "); scanf("%

24、s",&nameo); fp1=fopen(namei,"rb"); if(fp1=NULL) printf("can't open input segy filen"); exit(0); fp2=fopen(nameo,"wb"); if(fp2=NULL) printf("can't open output segy filen"); exit(0); fp=fopen("原始数据.csv","w"); if(fp=NULL) prin

25、tf("can't open output filen"); exit(0); fpr=fopen("处理后数据.csv","w"); if(fpr=NULL) printf("can't open output filen"); exit(0); fread(&dat.buf0,1,LEN,fp1);/文件头部3200字节特殊编码部分的读取 fwrite(&dat.buf0,1,LEN,fp2);/把文件头部3200字节特殊编码部分写入 fread(&dat.buf0,1,

26、400,fp1);/文件头400字节二进制部分的读取 len1=dat.shd10;/sample numbers printf ("sample numbers: %d n",len1); printf ("sample rate(ms): %u n",dat.shd8/1000); printf ("seg_y flag: %d n",dat.shd12); fwrite(&dat.buf0,1,400,fp2);/把文件头400字节二进制部分写入 nflg=dat.shd12;/数据采样格式码 if(nflg=1) le

27、n=len1*4+240;/浮点数4字节 if(nflg=2) len=len1*4+240;/定点数4字节 if(nflg=3) len=len1*2+240;/定点数2字节 n=0; jhd=60;/道头240字节=60*4/=seg_y数据的读取=/ /*while (l=fread(&dat.buf0,1,len,fp1)>0) for(i=0;i<len1;i+)kii=dat.datai+jhd;if(n=1)for(i=0;i<len1;i+)xi=kii;fprintf(fp,"%fn",xi);for(i=0;i<len1;

28、i+)dat.datai+jhd=kii; fwrite(&dat.buf0,1,len,fp2);n+; */ for(i=0;i<Ns;i+) bi=0; b0=1;/形成方程组右侧数据 while (l=fread(&dat.buf0,1,len,fp1)>0) for(i=0;i<len1;i+)kii=dat.datai+jhd; /=第n道地震记录做反滤波=/if(n=1)for(i=0;i<len1;i+)xi=kii;fprintf(fp,"%fn",xi);autocorrelation(x,rss,Ns);seq_Toeplitz(rss,Ns,b,a);convolution(x,a,r,Ns,Ns);for(i=0;i<Ns;i+)fprintf(fpr,"%fn",ri);/=对所有道地震记录做滤波=/ /*for(i=0;i<len1;i+)xi=kii; autocorrelation(x,rss,Ns);seq_Toeplitz(rss,Ns,b,a);convolution(x,a,r,Ns,Ns);for(i=0;i<Ns;