記錄數值模擬的褶積模型法實驗報告

1、地震記錄數值模擬的褶積模型法實驗報告專 業(yè)： 勘查技術以工程學 號： XXXXXXXXXXXXXXX姓 名： XXXXXXXXXXXXX實驗題目：地震記錄數值模擬的褶積模型法實驗報告一、 實驗目的掌握褶積模型基本理論;實現方法與程序編制;由褶積模型初步分析地震信號的分辨率問題。二、 實驗原理和方法：1、褶積原理地震勘探的震源往往是帶寬很寬的脈沖，在地下傳播、反射、繞射到測線，傳播經過中高頻衰減，能量被吸收。吸收過程可以看成濾波的過程，濾波可以用褶積完成。在濾波中，反射系數與震源強弱關聯，吸收作用與子波關聯。最簡單的地震記錄數值模擬，可以看成反射系數與子波的褶積。通常，反射系數是脈沖，子波取雷克

2、子波。（1） 雷克子波:（2） 反射系數:（3） 褶積公式:數值模擬地震記錄trace(t):trace(t) =rflct(t)*wave(t)反射系數的參數由z變成了t，怎么實現？在簡單水平層介質，分垂直和非垂直入射兩種實現，分別如圖1和圖2所示。1） 垂直入射：t=2*h/v2）非垂直入射：2、褶積方法（1） 離散化（數值化） 計算機數值模擬要求首先必須針對連續(xù)信號離散化處理。反射系數在空間模型中存在，不同深度反射系數不同，是深度的函數。子波是在時間記錄上一延續(xù)定時間的信號，是時間的概念。在離散化時，通過深度采樣完成反射系數的離散化，通過時間采樣完成子波的離散化。如果記錄是Trace(t

3、)，則記錄是時間的函數，以時間采樣離散化。時間采樣間距以Dt表示，深度采樣間距以Dz表示。在做多道的數值模擬時，還有橫向x的概念，橫向采樣間隔以Dx表示。離散化的實現：或：D（2）離散序列的褶積å¥ -¥ =- 三、實驗內容1、 垂直入射地震記錄數值模擬的褶積模型；2、 非垂直入射地震記錄數值模擬的褶積模型；3、 點繞射的地震記錄數值模擬的褶積模型；三、 實驗步驟（一） 實驗一：1、根據垂直入射褶積模型理論算法，填充程序（附后）的下劃線部分，使程序完整，調試程序，算出結果，用“Fimage”顯示軟件顯示褶積結果；2、 根據非零偏移距算法，編制非零偏移距褶積模型程序

4、，算出結果，用“Fimage”顯示軟件顯示褶積結果。（參考垂直入射褶積模型理論算法和程序，子波與反射層不變）；3、變換子波的長度Nw和主頻fm，重復1和2；（1）Nw=80的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（二）實驗二：兩層反射界面如圖3所示為兩層反射界面的模型，對比不同間距H12情況下，不同長子波長度Nw和主頻fm的褶積模型。1、H12=20m時：（1）Nw=80的情況下：1）fm=5，2）fm=10，3）fm=15，4）fm=20，

5、5）fm=55；（2）Nw=160的情況下：1）fm=5，2）fm=10，3）fm=15，4）fm=20，5）fm=55；2、H12=40m時：（1）Nw=80的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；3、H12=60m時：（1）Nw=80的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；4、H12=8

6、0m時：（1）Nw=80的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情況下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（三）點繞射的褶積模型（1）地質模型與參數：如圖4所示為點繞射模型（下圖）及其地震記錄（上圖）。參數為：Nx=128，Nz=200，V=3500m/s,h=1000m。如圖4圖下圖所示，地下繞射點位置為：(Ix0=NX/2,Iz0=1000/20)。其他參數與前述同。（2）走時計算如圖4所示，繞射點正上方檢波器的自激自收世間為t0：t0=2h/v任意位置檢波器的自激自收世間為ti

7、：（3）褶積模型爆炸反射界面?zhèn)鞑サ綑z波點的脈沖值游俠是計算：î（四）觀察并記下實驗結果圖四、實驗結果分析1、 結果顯示（1）源程序段：實驗一.垂直入射地震記錄數值模擬的褶積模型單層反射界面/ 預處理部分 /#include<math.h>#include<stdio.h>#include<string.h>int Cnltn(float,float);int Rflct(float,float,float);int Wave(float,float);#define Nx 128#define Nt 256#define Nw 160/ 主程序波分

8、 /void main() float dt=0.004,dx=20,fm=20,h=1000,v=3000; int iflag_Co,iflag_Re,iflag_Wv; if(iflag_Wv=Wave(fm,dt)!=1)printf("Wave is error"); if(iflag_Re=Rflct(dt,h,v)!=1)printf("Reflection is error"); if(iflag_Co=Cnltn(dt,dx)!=1)printf("Convosion is error");/ 函數實現部分 /Wav

9、e Formaingfunction /int Wave(float fm,float dt) FILE *fpw; int It; float WaNw,t; double pai=3.1415926; if(fpw=fopen("wave.dat","wb")=NULL)printf("Connot open file ""wave"""); for(It=0;It<Nw;It+) t=(float)(It+1)*dt; / 雷克子波的計算 / WaIt=(float)(1-2*pai

10、*pai*fm*fm*t*t)*(float)exp(-2*pai*pai*fm*fm*t); fwrite(&WaIt,sizeof(WaIt),1,fpw); fclose(fpw); return(1);/ Reflect Formaing function /int Rflct(float dt,float h,float v) FILE *fpr; int It,Ix,Ltdpth; float t; float ReNt; if(fpr=fopen("Reflect.dat","wb")=NULL)printf("Conno

11、t open file ""Reflect"""); for(Ix=0;Ix<Nx;Ix+) for(It=0;It<Nt;It+) ReIt=0; t=(float)2*h/v; Ltdpth=(int)(t/dt); ReLtdpth=1; for(It=0;It<Nt;It+) fwrite(&ReIt,sizeof(ReIt),1,fpr); fclose(fpr); return(1);/ Convolution function /int Cnltn(float dt,float dx) FILE *fpc

12、,*fpw,*fpr; int It,Ix,Itao; float Wa1Nw,WaNw,ReNt+Nw+Nw,Re1Nt; float ConNt+Nw; if(fpc=fopen("Convosion.dat","wb")=NULL)printf("Connot open file ""Convosion"""); if(fpw =fopen("wave.dat","rb")=NULL)printf("Connot open file &qu

13、ot;"wave"""); if(fpr =fopen("Reflect.dat","rb")=NULL)printf("Connot open file ""Reflect"""); for(Ix=1;Ix<2;Ix+) for(It=0;It<Nw;It+) fread(&Wa1It,sizeof(Wa1It),1,fpw);/ 子波數據反褶 / WaNw-1-It=Wa1It; fclose(fpw); for(Ix=0;Ix&l

14、t;Nx;Ix+) for(It=0;It<Nt;It+) fread(&Re1It,sizeof(&Re1It),1,fpr); for(It=0;It<Nt+2*Nw;It+) ReIt=0.; for(It=0;It<Nt;It+) ReIt+Nw=Re1It; /反射系數數據移動: "0 到 Nt" 移為 "Nw 到Nw+Nt" / for(It=0;It<Nt+Nw;It+) ConIt=0; for(Itao=0;Itao<Nw;Itao+) ConIt=ConIt+WaItao*ReIt+Ita

15、o; / 褶積運算 / for(It=Nw/2;It<Nt+Nw/2;It+) fwrite(&ConIt,sizeof(ConIt),1,fpc); fclose(fpw); fclose(fpr); fclose(fpc); return(1);雙層反射界面/ 預處理部分 /#include<math.h>#include<stdio.h>#include<string.h>int Cnltn(float,float);int Rflct(float,float,float,float);int Wave(float,float);#def

16、ine Nx 128#define Nt 256#define Nw 160/ 主程序波分 /void main() float dt=0.004,dx=20,fm=15,h=1000,dh=80,v=4000; int iflag_Co,iflag_Re,iflag_Wv; if(iflag_Wv=Wave(fm,dt)!=1)printf("Wave is error"); if(iflag_Re=Rflct(dt,h,dh,v)!=1)printf("Reflection is error"); if(iflag_Co=Cnltn(dt,dx)!=

17、1)printf("Convosion is error");/ 函數實現部分 / Wave Formaing function /int Wave(float fm,float dt) FILE *fpw; int It; float WaNw,t; double pai=3.1415926; if(fpw=fopen("wave.dat","wb")=NULL)printf("Connot open file ""wave"""); for(It=0;It<Nw;I

18、t+) t=(float)(It+1)*dt; / 雷克子波的計算 / WaIt=(float)(1-2*pai*pai*fm*fm*t*t)*(float)exp(-2*pai*pai*fm*fm*t); fwrite(&WaIt,sizeof(WaIt),1,fpw); fclose(fpw); return(1);/ Reflect Formaing function /int Rflct(float dt,float h,float dh,float v) FILE *fpr; int It,Ix,Ltdpth,Ltdpth1; float t,t1; float ReNt;

19、if(fpr=fopen("Reflect.dat","wb")=NULL)printf("Connot open file ""Reflect"""); for(Ix=0;Ix<Nx;Ix+) for(It=0;It<Nt;It+) ReIt=0; t=(float)2*h/v; Ltdpth=(int)(t/dt); ReLtdpth=1; t1=(float)2*(h+dh)/v; Ltdpth1=(int)(t1/dt); ReLtdpth1=1; for(It=0;It&l

20、t;Nt;It+) fwrite(&ReIt,sizeof(ReIt),1,fpr); fclose(fpr); return(1);/ Convolution function /int Cnltn(float dt,float dx) FILE *fpc,*fpw,*fpr; int It,Ix,Itao; float Wa1Nw,WaNw,ReNt+Nw+Nw,Re1Nt; float ConNt+Nw; if(fpc=fopen("Convosion.dat","wb")=NULL)printf("Connot open fil

21、e ""Convosion"""); if(fpw =fopen("wave.dat","rb")=NULL)printf("Connot open file ""wave"""); if(fpr =fopen("Reflect.dat","rb")=NULL)printf("Connot open file ""Reflect"""); for(

22、Ix=1;Ix<2;Ix+) for(It=0;It<Nw;It+) fread(&Wa1It,sizeof(Wa1It),1,fpw);/子波數據反褶 / WaNw-1-It=Wa1It; fclose(fpw); for(Ix=0;Ix<Nx;Ix+) for(It=0;It<Nt;It+) fread(&Re1It,sizeof(&Re1It),1,fpr); for(It=0;It<Nt+2*Nw;It+) ReIt=0.; for(It=0;It<Nt;It+) ReIt+Nw=Re1It; /反射系數數據移動: "

23、;0 到 Nt" 移為 "Nw 到Nw+Nt" / for(It=0;It<Nt+Nw;It+) ConIt=0; for(Itao=0;Itao<Nw;Itao+) ConIt=ConIt+WaItao*ReIt+Itao; / 褶積運算 / for(It=Nw/2;It<Nt+Nw/2;It+) fwrite(&ConIt,sizeof(ConIt),1,fpc); fclose(fpw); fclose(fpr); fclose(fpc); return(1);實驗二.非垂直入射地震記錄數值模擬的褶積模型/預處理部分/#includ

24、e<math.h> #include<stdio.h>#include<string.h>int Cnltn(float,float);int Rflct(float,float,float,float);int Wave(float,float);#define Nx 128#define Nt 256#define Nw 60/主程序波分/void main() float dt=0.004,dx=20,fm=5,h=1000,v=3000;int iflag_Co,iflag_Re,iflag_Wv;if(iflag_Wv=Wave(fm,dt)!=1

25、)printf("Wave is error");if(iflag_Re=Rflct(dt,h,v,dx)!=1)printf("Reflection is error");if(iflag_Co=Cnltn(dt,dx)!=1)printf("Convosion is error");/ 3.函數實現部分/ 3.1 Wave Formaing function/int Wave(float fm,float dt)FILE *fpw;int It;float WaNw,t;double pai=3.1415926;if(fpw=fo

26、pen("wave.dat","wb")=NULL)printf("Connot open file ""wave""");for(It=0;It<Nw;It+)t=(float)(It+1)*dt;WaIt=(float)(1.-2*pow(pai*fm*t),2)*(float)exp(-2*t*pow(pai*fm),2);/雷克子波的計算/fwrite(&WaIt,sizeof(WaIt),1,fpw);fclose(fpw);return(1);/3.2 Reflect

27、Formaing function/int Rflct(float dt,float h,float v,float dx)FILE *fpr;int It,Ix,Ltdpth;float t;float ReNt;if(fpr=fopen("Reflect.dat","wb")=NULL)printf("Connot open file ""Reflect""");for(Ix=0;Ix<Nx;Ix+)for(It=0;It<Nt;It+)ReIt=0.;t=(float)sqrt

28、(4*h*h+Ix*Ix*dx*dx)/v;Ltdpth=(int)(t/dt);ReLtdpth=1.;for(It=0;It<Nt;It+)fwrite(&ReIt,sizeof(ReIt),1,fpr);fclose(fpr);return(1);/3.2 Convolution function/int Cnltn(float dt,float dx)FILE *fpc,*fpw,*fpr;int It,Ix,Itao;float Wa1Nw,WaNw,ReNt+Nw+Nw,Re1Nt;float ConNt+Nw; if(fpc=fopen("Convosio

29、n.dat","wb")=NULL)printf("Connot open file ""Convosion""");if(fpw =fopen("wave.dat","rb")=NULL)printf("Connot open file ""wave"""); if(fpr =fopen("Reflect.dat","rb")=NULL)printf("C

30、onnot open file ""Reflect""");for(Ix=1;Ix<2;Ix+)for(It=0;It<Nw;It+)fread(&Wa1It,sizeof(Wa1It),1,fpw);WaNw-1-It=Wa1It;/子波數據反褶/fclose(fpw);for(Ix=0;Ix<Nx;Ix+)for(It=0;It<Nt;It+)fread(&Re1It,sizeof(&Re1It),1,fpr);for(It=0;It<Nt+2*Nw;It+)ReIt=0.;for(It

31、=0;It<Nt;It+)ReNw+It=Re1It;/反射系數數據移動: "0 到Nt" 移為"Nw 到Nw+Nt"/for(It=0;It<Nt+Nw;It+)ConIt=0; for(Itao=0;Itao<Nw;Itao+)ConIt=ConIt+WaItao*ReIt+Itao;/褶積運算/for(It=Nw/2;It<Nt+Nw/2;It+)fwrite(&ConIt,sizeof(ConIt),1,fpc);fclose(fpw);fclose(fpr);fclose(fpc);return(1);實驗三.點

32、繞射的褶積模型/ 預處理部分 / #include <math.h>#include <stdio.h>#include <string.h>int Cnltn(float,float);int Rflct(float,float,float,float);int Wave(float,float);#define Nx 128#define Nt 256#define Nw 80/ 主程序部分 /int main()float dt=0.004,dx=20.,fm=15.,h=1000.,v=3500.;int iflag_Co,iflag_Re,ifla

33、g_Wv;if(iflag_Wv=Wave(fm,dt)!=1)printf("Wave is error");if(iflag_Re=Rflct(dt,dx,h,v)!=1)printf("Reflection is error");if(iflag_Co=Cnltn(dt,dx)!=1)printf("Convosion is error");return 0;/函數實現部分 / 3.1 Wave Formaing function / int Wave(float fm,float dt)FILE *fpw;int It;flo

34、at WaNw,t;double pai=3.1415926;if(fpw=fopen("wave.dat","wb")=NULL)printf("Connot open file Wave");for(It=0;It<Nw;It+)t=(float)(It+1)*dt;WaIt=(float)(1.-2*pai*pai*fm*fm*t*t)*(float)exp(-pai*pai*fm*fm*t);fwrite(&WaIt,sizeof(WaIt),1,fpw);fclose(fpw);return 1;/3.2 Re

35、flect Formaing function /int Rflct(float dt,float dx,float h,float v)FILE *fpr;int It,Ix,Ltdpth;float t;float ReNt;if(fpr=fopen("Reflect.dat","wb")=NULL)printf("Connot open file Reflect");for(Ix=0;Ix<Nx;Ix+)for(It=0;It<Nt;It+)ReIt=0.;t=(float)2*sqrt( h*h+pow(Ix-Nx

36、/2)*dx,2) )/v;Ltdpth=(int)(t/dt);ReLtdpth=1.;for(It=0;It<Nt;It+)fwrite(&ReIt,sizeof(ReIt),1,fpr);fclose(fpr);return 1;/ 3.3 Convolution function /int Cnltn(float dt,float dx)FILE *fpc,*fpw,*fpr;int It,Ix,Itao;float Wa1Nw,WaNw,ReNt+Nw+Nw,Re1Nt;float ConNt+Nw;if(fpc=fopen("Convosion.dat&qu

37、ot;,"wb")=NULL)printf("Connot open file Convosion");if(fpw=fopen("wave.dat","rb")=NULL)printf("Connot open file Wave");if(fpr=fopen("Reflect.dat","rb")=NULL)printf("Connot open file Reflect"); for(It=0;It<Nw;It+)fread

38、(&Wa1It,sizeof(Wa1It),1,fpw);WaNw-1-It=Wa1It;fclose(fpw);for(Ix=0;Ix<Nx;Ix+)for(It=0;It<Nt;It+)fread(&Re1It,sizeof(&Re1It),1,fpr);for(It=0;It<Nt+2*Nw;It+)ReIt=0.;for(It=0;It<Nt;It+)ReNw+It=Re1It;for(It=0;It<Nt+Nw;It+)ConIt=0;for(Itao=0;Itao<Nw;Itao+)ConIt=ConIt+WaItao*R

39、eIt+Itao;for(It=Nw/2;It<Nt+Nw/2;It+)fwrite(&ConIt,sizeof(ConIt),1,fpc);fclose(fpw);fclose(fpr);fclose(fpc);return 1;2、 對比分析實驗一：垂直入射地震記錄數值模擬的褶積模型A單層反射界面1、 當h=1000m，v=3000m/s,Nw=80時，fm=5 fm=10 fm=15 fm=20 圖件依次如下分析：當地層深度h、子波長度Nw和地震波在地層中的傳播速度v固定時，隨著子波頻率fm的增大，檢波器接收到地震記錄的時間不變，但持續(xù)時間變長。在“Fimage”顯示軟件中

40、表現為顯示的記錄隨著fm的增大，同相軸依次變細。2、 當h=1000m，v=3000m/s,fm=10Hz時，Nw=80 Nw=160 圖件依次如下：分析：可以看出，當地層深度h、子波頻率fm和地震波在地層中的傳播速度v固定時，隨著子波長度Nw的增大，檢波器接收到地震記錄的時間變大。在“Fimage”顯示軟件中表現為顯示的記錄整體向上移動；并且隨著的子波長度的均勻增加，地震記錄向上移動是均勻的。B.雙層反射界面1、H12=20m時，（1）當h=1000m，v=4000m/s，Nw=80時，fm=5 fm=10 fm=15 fm=20 圖件(2)當h=1000m，v=4000m/s，Nw=160時fm=5 fm=10 fm=15 fm=20圖件2、H12=40m時，（1）當h=1000m，v=4000m/