地震子波波形顯示及一維地震合成記錄

1、地震勘探原理實(shí)驗(yàn)一地震子波波形顯示及一維地震合成記錄姓名： 學(xué)號(hào)：專業(yè)：地球物理勘察技術(shù)與工程 級(jí) 一、實(shí)驗(yàn)?zāi)康?. 認(rèn)識(shí)子波，對(duì)子波的波形有直觀的認(rèn)識(shí)。（名詞：零相位子波，混合相位子波，最小相位子波；了解子波的分辨率與頻寬的關(guān)系；）2. 利用褶積公式合成一維地震記錄。二、實(shí)驗(yàn)步驟1. 雷克子波 零相位子波 （最小相位子波） n= m1/m2為最大波峰m1和最大波谷m2之比 鐘型子波 xw為初相m為時(shí)間域主波峰與次波峰之比w(t)=exp(-2*Fm2*t2*ln(n)*sin(T-2*pi*Fm*t) n=m1/m2 最大相位子波（最大相位子波請同學(xué)們自己查找相關(guān)文獻(xiàn)完成，非必須完成）其中代

2、表子波的中心頻率， t =i*dt，dt為時(shí)間采樣間隔，i為時(shí)間離散點(diǎn)序號(hào); 這里可以為 = 10，25，40，100 Hz等，采樣間隔dt=0.002秒，i為0256;2. 根據(jù)公式編程實(shí)現(xiàn)不同頻率的零相位子波的波形顯示；不同中心頻率的零相位子波圖 = 25: = 100:3. 其地質(zhì)模型為：震源檢波器第一層介質(zhì)r(100)=1.0第二層介質(zhì)r(200)=-0.7第三層介質(zhì)r(300)=0.5第四層介質(zhì)r(400)=0.4第五層介質(zhì)r(450)=0.6設(shè)計(jì)反射系數(shù)(n=512)，n為地層深度，其中，為第一層介質(zhì)深度；，為第二層介質(zhì)深度；，為第三層介質(zhì)深度；，為第四層介質(zhì)深度；，為第五層介質(zhì)深

3、度；其它為0。地震波在介質(zhì)中傳播，當(dāng)?shù)竭_(dá)介質(zhì)分界面時(shí)，發(fā)生反射和透射，反射波被檢波器接受，生成地震記錄。反射系數(shù)表示地震波在兩層介質(zhì)分界面的能量重新分配，如r(100)=1.0，表示地震波入射到分界面時(shí)，只有一種波，反射縱波（或反射橫波）。反射系數(shù)不為1.0時(shí)，表示當(dāng)?shù)卣鸩ㄈ肷涞椒纸缑鏁r(shí)，產(chǎn)生兩種反射波。反射系數(shù)為正，表示反射波相位與入射波相位相差2；反射系數(shù)為負(fù)，表示反射波相位與入射波相位相差，存在半波損失。4. 應(yīng)用褶積公式合成一維地震記錄，并圖形顯示；應(yīng)用褶積公式求f（n）的程序?yàn)椋?include#include#define PI 3.1415926#define FM 100voi

4、d main() double fac(double x,double y,double z,int m,int n); FILE *fp; int i,j,x; double W,dt=0.002,t,a256; double b512=0; double r512=0; r100=1.0; r200=-0.7; r300=0.5; r400=0.4; r450=0.6; fp=fopen(Date.txt,w+); printf(please input x:n); scanf(%d,&x); for(i=0;i256;i+)t=i*dt; if(x=1) W=exp(-2*FM*FM*t

5、*t*log(2)*sin(2*PI*FM*t); else if (x=2) W=(1-2*pow(PI*FM*t,2)*exp(-pow(PI*FM*t,2); else if (x=3) W=exp(-FM*FM*t*t*log(2)*cos(2*PI*FM*t+PI/4); ai=W;fac(r,a,b,512,256);for(j=0;j512;k+)fprintf(fp,%fn,bj);double fac(double x,double y,double z,int m,int n) int i,j; for(i=0;i=m+n-1;i+) double sum=0.0; for

6、(j=0;j0&i-j=256) sum+=xj*yi-j; zi=sum;三、實(shí)驗(yàn)結(jié)果一維反射系數(shù)序列圖形顯示為：零相位子波與反射系數(shù)褶積后的地震記錄圖形顯示： = 25: = 100:最小相位子波與反射系數(shù)褶積后的地震記錄圖形顯示： = 25: = 100:混合相位子波與反射系數(shù)褶積后的地震記錄圖形顯示： = 25: = 100:最大相位子波與反射系數(shù)褶積后的地震記錄圖形顯示： = 25: = 100:零相位振幅圖形顯示： = 25: = 100:零相位幅角圖形顯示： = 25: = 100:最小相位振幅圖形顯示： = 25: = 100:最小相位幅角圖形顯示： = 25: = 100:混

7、合相位幅角圖形顯示： = 25: = 100:混合相位振幅圖形顯示： = 25: = 100:最大相位幅角圖形顯示： = 25: = 100:最大相位振幅圖形顯示： = 25: = 100:四、實(shí)驗(yàn)分析根據(jù)所學(xué)知識(shí)對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行分析；地震子波由震源激發(fā)，在地層中傳播，因?yàn)樵诔练e地層中，每層介質(zhì)的物理性質(zhì)不相同，從而使得地震波的傳播速度也不相同。當(dāng)?shù)卣鸩▊鞑サ絻蓪咏橘|(zhì)的分界面時(shí)，會(huì)發(fā)生反射，由于每層介質(zhì)的反射系數(shù)不同，所以反射波的能量也不相同，檢波器接收到不同時(shí)刻的、不同能量的反射波，形成一個(gè)地震記錄。由合成地震記錄中可以看出，最小相位子波相對(duì)零相位子波來說是相位滯后的，能量延遲的，但兩者為同一

8、家族的子波。合成地震記錄中橫坐標(biāo)為時(shí)間，縱坐標(biāo)為振幅。每一時(shí)刻的值由m個(gè)值的和組成，m為反射系數(shù)r(n)的長度，整個(gè)地震記錄由m+n-1個(gè)時(shí)刻的值組成。對(duì)于零相位的地震記錄來說，當(dāng)r(m)=1.0時(shí)，即j=100時(shí)，i=100時(shí)，w(i-j)=1.0，是能量最大的，即w(0)=1.0。同理，當(dāng)n=200,、300、400、450時(shí)，能量也是最大的。對(duì)于最小相位的地震記錄來說，當(dāng)r(m)=1.0時(shí)，即j=100時(shí)，但i=100時(shí)，w(i-j)不是最大能量的，即最大能量不是在w(0)出現(xiàn)，而是延遲出現(xiàn)。同理，當(dāng)n=200、300、400、450時(shí)，能量也不是最大的，而是要延遲出現(xiàn)。由振幅圖及幅角圖

9、可知，零相位子波能量聚集在首部，開始時(shí)就具有最大能量，無積累過程，當(dāng)振幅最大時(shí)，相位為零，即此時(shí)波的振幅為實(shí)數(shù)，達(dá)到最大值；最小相位子波能量聚集在序列首部，是最小能量延遲的，信號(hào)隨時(shí)間增大而減小，當(dāng)振幅最大時(shí)，相位不為零，是非零相位的，相對(duì)零相位子波來說，最大能量是延遲的；混合相位子波的能量聚集在序列中部，是混合能量延遲的；最大相位子波能量聚集在后部。最大相位子波和混合相位子波的信號(hào)信號(hào)不隨時(shí)間增大而減小。五、附：源程序代碼 #include#include13KFFT.C#include#define PI 3.1415926#define FM 100void main() double

10、fac(double x,double y,double z,int m,int n); FILE *fp,*fpr,*fpre,*fpi,*fpamp,*fpha; int i,j,x; double W,dt=0.002,t,a256,pr512,pi512=0.0,fr512,fi512,amp512,pha512; double b512=0; double r512=0; r100=1.0; r200=-0.7; r300=0.5; r400=0.4; r450=0.6; fp=fopen(褶積結(jié)果.txt,w+);fpr=fopen(反射系數(shù).csv,w+);fpre=fopen(

11、實(shí)部.txt,w+);fpi=fopen(虛部.txt,w+);fpamp=fopen(振幅.csv,w+);fpha=fopen(相位.csv,w+); for(i=0;i512;i+) fprintf(fpr,%fn,ri); fclose(fpr); printf(please input x:n); scanf(%d,&x); for(i=0;i256;i+)t=i*dt; if(x=1) W=(1-2*pow(PI*FM*t,2)*exp(-pow(PI*FM*t,2); else if (x=2) W=exp(-2*FM*FM*t*t*log(2)*sin(2*PI*FM*t);

12、else if (x=3) W=exp(-2*pow(FM*t,2)*log(2)*sin(0.512-2*PI*FM*t);else if(x=4)W=exp(-FM*FM*t*t*log(2)*cos(2*PI*FM*t+PI/4); ai=W;fac(r,a,b,512,256);for(j=0;j512;j+)fprintf(fp,%fn,bj);for(i=0;i512;i+)pri=bj;for(i=0;i512;i+)pri=fri+127;kfft(pr,pi,512,9,fr,fi,0,1);for(i=0;i512;i+)fprintf(fpre,%en,fri);fprintf(fpi,%en,fii); fprintf(fpamp,%en,pri);fprintf(fpha,%fn,pii