油氣勘探方法程序設(shè)計(jì)課程設(shè)計(jì)(論文)-最小平方反濾波

1、油氣勘探方法程序設(shè)計(jì)課程設(shè)計(jì)(論文)設(shè)計(jì)（論文）題目 最小平方反濾波 學(xué)院名稱 地球物理學(xué)院 專業(yè)名稱 勘查技術(shù)與工程（石油物探） 學(xué)生姓名 學(xué)生學(xué)號(hào) 任課教師 田仁飛 設(shè)計(jì)（論文）成績(jī) 教務(wù)處 制2015年7月1日最小平方反濾波一、方法原理由地震波的傳播理論可知，在粘彈性介質(zhì)中地震波是以地震子波的形式在地下傳播，地面接收到的反射波地震記錄是地層反射系數(shù)與地震子波的褶積，因此，地層相當(dāng)于一個(gè)濾波器，是反射系數(shù)序列變成了由子波組成的地震記錄，降低了地震勘探的縱向分辨率。反濾波的目的就是要設(shè)計(jì)一個(gè)反濾波器，來對(duì)地震記錄濾波，消除地層濾波的作用，提高地震記錄的縱向分辨率。最小平方反濾波的基本思想在于

2、設(shè)計(jì)一個(gè)濾波算子，用它把已知的輸入信號(hào)轉(zhuǎn)換為與給定的期望輸出信號(hào)在最小平方誤差的意義下是最佳輸出。1）用最小平方法求反濾波因子對(duì)輸入子波b(t)反濾波后的期望輸出為d(t),實(shí)際輸出為y(t)，按最小平方原理，使二者的誤差平方和Q為最小時(shí)求得的反濾波因子稱為最小平方反濾波因子，用它對(duì)地震記錄x(t)進(jìn)行的反濾波為最小平方反濾波。設(shè)輸入離散信號(hào)為地震子波b(n)=b(1),b(2), ,b(m),待求的反濾波因子a(n)=a(0),a(1),a(2), ,a(m), a(t)的起始時(shí)間為0，（m+1）為a(t)的延續(xù)長(zhǎng)度，b(n)與a(n)的褶積為實(shí)際輸出y(n),即 yn=an*bn=0mab

3、n- 1實(shí)際輸出與期望輸出的誤差平方和為Q=n=0myn-d(n)2 =n=0m=0mabn-d(n)2 （2）要使Q為最小，數(shù)學(xué)上就是求Q的極值問題，即求滿足 Qa(l)=0 l=0,1,m （3）的濾波因子a(t)。即滿足方程 =0marbbl-=rbdl l=0,1,m （4）其中rbbl-=n=0mbn-b(n-l)為地震子波的自相關(guān)函數(shù)，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）上式系數(shù)矩陣稱為托布利茲方程。當(dāng)期

4、望輸出是單位脈沖t時(shí)，即 dt=t=1,t=00,t0 （6）則 rbdl=b-l （7） 由于b(t)是地震子波，b(t)=0,當(dāng)t<0時(shí)，上述矩陣可寫成： rbb(0)rbb(1)rbb(m)rbb(1)rbb(m)rbb(0)rbb(m-1)rbb(m-1)rbb(0)a0a1am=b000=100 （8）用地震記錄自相關(guān)系數(shù)rxx()代替ybb()2）已知地震記錄x(t) x(t)=b(t)* (t) （9）其中x(t)為地震記錄，(t)是反射系數(shù)。根據(jù)自相關(guān)函數(shù)的性質(zhì)，可證明：地震記錄自相關(guān)rxx(n)是地震子波自相關(guān)rbb(n)和反射系數(shù)自相關(guān)r(n)的褶積,即： rxxn=

5、rbbn*rn （10）假定反射系數(shù)是隨機(jī)白噪音，則 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)。二、方法步驟在地震記錄上選取自相關(guān)段，要選擇在波形單一，能量強(qiáng)的地方。自相關(guān)段的長(zhǎng)度應(yīng)大于或等于子波的長(zhǎng)度；求選取段的自相關(guān)函數(shù)： rxxn=t=1m-nxt×xt+n, n=0,1,2,M-1 （15） 其中M為自相關(guān)長(zhǎng)度；把rxxn

6、代入式中求解at，并由式得反濾波因子a(t)，式中，通常取N0=1；將地震記錄x(t)與反濾波因子a(t)褶積，完成反濾波； yt=0m-1a×xt-, t=m,m+1,L （16）其中m為反濾波因子長(zhǎng)度，L為地震記錄長(zhǎng)度。三、源程序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 /反射系數(shù)序列長(zhǎng)度#define Nw 100 /子波長(zhǎng)度#define P 199 /合成地震記錄長(zhǎng)度/=建立速度模型=/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;/=形成反射系數(shù)序列=/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);/=求褶積函數(shù)=/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; /=自相關(guān)函數(shù)=/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為矩陣的第一行或者第一列數(shù)據(jù)，b為方程組右側(cè)數(shù)據(jù)，x為計(jì)算的反濾波因子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); /=形成反射系數(shù)序列=/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); /=子波自相關(guān)=/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); /=形成方程組右側(cè)數(shù)據(jù)=/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 /=自相關(guān)函數(shù)=/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; /=求褶積函數(shù)=/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為矩陣的第一行或者第一列數(shù)據(jù)，b為方程組右側(cè)數(shù)據(jù)，x為計(jì)算的反濾波因子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("原始數(shù)據(jù).csv","w"); if(fp=NULL) prin

25、tf("can't open output filen"); exit(0); fpr=fopen("處理后數(shù)據(jù).csv","w"); if(fpr=NULL) printf("can't open output filen"); exit(0); fread(&dat.buf0,1,LEN,fp1);/文件頭部3200字節(jié)特殊編碼部分的讀取 fwrite(&dat.buf0,1,LEN,fp2);/把文件頭部3200字節(jié)特殊編碼部分寫入 fread(&dat.buf0,1,

26、400,fp1);/文件頭400字節(jié)二進(jìn)制部分的讀取 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字節(jié)二進(jìn)制部分寫入 nflg=dat.shd12;/數(shù)據(jù)采樣格式碼 if(nflg=1) le

27、n=len1*4+240;/浮點(diǎn)數(shù)4字節(jié) if(nflg=2) len=len1*4+240;/定點(diǎn)數(shù)4字節(jié) if(nflg=3) len=len1*2+240;/定點(diǎn)數(shù)2字節(jié) n=0; jhd=60;/道頭240字節(jié)=60*4/=seg_y數(shù)據(jù)的讀取=/ /*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;/形成方程組右側(cè)數(shù)據(jù) 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);/=對(duì)所有道地震記錄做濾波=/ /*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;