疊加地震記錄的相移波動方程正演模擬數(shù)值模擬實(shí)驗(yàn)書

1、地震數(shù)值模擬實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)時(shí)間2013年6月1-30日開課單位地球物理學(xué)院指導(dǎo)教師XXXXXXXXXXXXXX實(shí)驗(yàn)題目： 疊加地震記錄的相移波動方程正演模擬姓名學(xué)號班級勘查四班專業(yè)勘查技術(shù)與工程院系地球物理學(xué)院單項(xiàng)成績內(nèi)容理解寫作結(jié)構(gòu)程序設(shè)計(jì)模型設(shè)計(jì)計(jì)算結(jié)果結(jié)果分析總 成 績一、 實(shí)驗(yàn)題目疊加地震記錄的相移波動方程正演模擬二、 實(shí)驗(yàn)?zāi)康?.掌握各向同性介質(zhì)任意構(gòu)造、水平層狀速度結(jié)構(gòu)地質(zhì)模型的相移波動方程正演模擬基本理論2.實(shí)現(xiàn)方法與程序編制3.由正演記錄初步分析地震信號的分辨率。三、 實(shí)驗(yàn)原理1、 地震波傳播的波動方程設(shè)（x,z）為空間坐標(biāo)，t為時(shí)間，地震波傳播速度為v(x,z)，則二位介質(zhì)中任

2、意位置、任意時(shí)刻的地震波場為p(z,x,t)：壓縮波縱波。則二維各向同性均勻介質(zhì)中地震波傳播的遵循聲波方程為2、 傅里葉變換的微分性質(zhì)p(t)與其傅里葉變換的P(w)的關(guān)系：ïïîï ï íì = = ò ò¥¥ - ¥3、 地震波傳播的相移外推公式 令速度v不隨x變化，只隨z變化，則利用傅里葉變換微分性質(zhì)把波動方程(變換到頻率-波數(shù)域，得：4、 初始條件和邊界條件按照爆炸界面理論，反射界面震源在t=0時(shí)刻同時(shí)起爆，此時(shí)刻的波場就是震源。根據(jù)不同情況，可直接使用反射系數(shù)脈沖或

3、子波作震源。如果直接使用反射系數(shù)作震源脈沖，則初始條件可表示為：5、邊界處理（1）邊界反射問題把實(shí)際無窮空間區(qū)域中求解波場的問題化為有窮區(qū)域求解時(shí)，左右兩邊使用零邊界條件。物理上假設(shè)探區(qū)距Xmin與Xmax兩個(gè)端點(diǎn)很遠(yuǎn)，在兩個(gè)端點(diǎn)上收到的反射波很弱。但是，上述條件在實(shí)際中不能成立，造成零邊界條件反而成為絕對阻止波通過的強(qiáng)反射面。在正演模擬的剖面上出現(xiàn)了邊界假反射干涉正常界面的反射。（2）邊界強(qiáng)反射的處理鑲邊法、削波法、吸收邊界都能有效消除邊界強(qiáng)反射。削波法就是在波場延拓過程中，沒延拓一次，在其兩側(cè)均勻衰減到零，從而消除邊界強(qiáng)反射的影響。假設(shè)橫向總長度為NX，以兩邊Lx道吸波為例，有以下吸波公式

4、：四、 實(shí)驗(yàn)內(nèi)容1、 基本要求（1） 點(diǎn)繞射構(gòu)造和水平層狀速度模型（參數(shù)如圖1所示）的正演數(shù)值模擬；五、 實(shí)驗(yàn)步驟1、參數(shù)初始化；2、形成邊界削波數(shù)據(jù)；3、波場初始化；4、Zmax層波場延拓到深度Zmax-1；5、Zi+1層波場延拓到深度Zi；6、重復(fù)5，從Iz=Nz-1開始，直到Iz=1，得測線上的頻率空間域波場；7、頻率-空間域波場對頻率做反傅里葉變換，得時(shí)間-空間波場；8、使用Fimage軟件顯示所得結(jié)果。六、實(shí)驗(yàn)編程文檔可自由編輯打印 /*1.頭文件*/#include<stdio.h>#include<math.h>#include<conio.h>

5、;#include<stdlib.h>#include<string.h>/-2.Function Request功能要求函數(shù)說明-int kkfft(float *,float *,int,int);int Absorb(int); /削波函數(shù)int Rflct(); /反射函數(shù)int Vlcty(); /速度函數(shù)int WvFld0(); /波長函數(shù)int exp_ikzDz(float *,int,float,int,float,float);/int PsFrwd();/int Po2Judge(int);/-#define Nx 128 /-3.參數(shù)設(shè)置定義符

6、號- -#define Nt 256#define Nz 100#define Dx 20#define Dt 0.004#define Dz 20#define pai 3.1415926int main()int Labs=10; /定義削波邊界范圍 if(Po2Judge(Nt) !=1) printf("Nt=%d is not the Power of 2n",Nt);exit(0); if(Po2Judge(Nx) !=1) printf("Nx=%d is not the Power of 2n",Nt);exit(0); if(Absor

7、b(Labs) !=1) printf("Absorb is errorn"); exit(0); if(Rflct() !=1) printf("Rflction is errorn"); exit(0); if(Vlcty() !=1) printf("Vlcty is errorn"); exit(0); if(WvFld0() !=1) printf("WvFld is errorn"); exit(0); if(PsFrwd() !=1) printf("PsFrwd is errorn&quo

8、t;); exit(0); return 1;int Po2Judge(int N) /判斷是否是2的倍數(shù)/int k=0;long Ln=0;for(k=0;N-Ln>0;k+)Ln=(long)pow(2,k);Ln=(long)pow(2,k-1);if (fabs(Ln-N)>=1)return (0);return 1;/定義快速傅立葉函數(shù)/int kkfft(float pr,float pi,int n,int l) int it,m,is,i,j,nv,l0,il=0; float p,q,s,vr,vi,poddr,poddi;float fr4096,fi409

9、6;int k=0;long Ln=0;for(k=0;n-Ln>0;k+)Ln=(long)pow(2,k);k=k-1;for (it=0; it<=n-1; it+)m = it;is = 0;for(i=0; i<=k-1; i+)j = m/2;is = 2*is+(m-2*j);m = j;frit =pris;fiit = piis;pr0 = 1.0; pi0 = 0.0;p = 6.283185306/(1.0*n);pr1 = (float) cos(p); pi1 = -(float)sin(p);if (l!=0) pi1=-pi1;for (i=2;

10、 i<=n-1; i+) p = pri-1*pr1; q = pii-1*pi1;s = (pri-1+pii-1)*(pr1+pi1);pri = p-q; pii = s-p-q;for (it=0; it<=n-2; it=it+2) vr = frit; vi = fiit;frit = vr+frit+1; fiit = vi+fiit+1;frit+1 = vr-frit+1; fiit+1 = vi-fiit+1;m = n/2; nv = 2;for (l0=k-2; l0>=0; l0-) m = m/2; nv = 2*nv;for(it=0; it&l

11、t;=(m-1)*nv; it=it+nv)for (j=0; j<=(nv/2)-1; j+) p = prm*j*frit+j+nv/2;q = pim*j*fiit+j+nv/2;s = prm*j+pim*j;s = s*(frit+j+nv/2+fiit+j+nv/2);poddr = p-q; poddi = s-p-q;frit+j+nv/2 = frit+j-poddr;fiit+j+nv/2 = fiit+j-poddi; frit+j = frit+j+poddr;fiit+j = fiit+j+poddi;if(l!=0)for(i=0; i<=n-1; i+

12、) fri = fri/(1.0*n);fii = fii/(1.0*n); if(il!=0)for(i=0; i<=n-1; i+) pri = sqrt(fri*fri+fii*fii);if(fabs(fri)<0.000001*fabs(fii) if (fii*fri)>0) pii = 90.0;else pii = -90.0;elsepii = atan(fii/fri)*360.0/6.283185306;for(i=0;i<n;i+)pri=fri;pii=fii;return(1); /*調(diào)用削波函數(shù)，形成削波數(shù)據(jù)存入一個(gè)文件*/int Abso

13、rb(int Lx) /Nx已為全局變量不參與傳遞 FILE *fp_Abs;int Ix;float AbsNx;if(fp_Abs=fopen("Absb.dat","wb")=NULL) printf("Connot open file""Absb"""); for(Ix=0;Ix<Nx;Ix+)AbsIx=1;/*for(Ix=0;Ix<Lx-1;Ix+)AbsIx=sqrt(sin(pai/2)*(Ix/(Lx-1);AbsNx-Ix-1=AbsIx;/*for(Ix=0;

14、Ix<Nx;Ix+)fwrite(&AbsIx,sizeof(AbsIx),1,fp_Abs);fclose(fp_Abs);return 1; /*反射系數(shù)_構(gòu)造形態(tài)模型的生成*/int Rflct() FILE *fp_Rflct;int Ix,Iz;float RflctNz;if(fp_Rflct=fopen("Rflct.dat","wb")=NULL) printf("Connot open file""Reflection""");for(Ix=0;Ix<Nx;

15、Ix+)for(Iz=0;Iz<Nz;Iz+)RflctIz=0;/*if(Ix=Nx/2-1&&Iz=Nz/2-1)/*RflctIz=1;/*fwrite(&RflctIz,sizeof(RflctIz),1,fp_Rflct);fclose(fp_Rflct);return 1;int Vlcty() /*速度模型的生成*/ FILE *fp_Vlcty;int Ix,Iz;float VlctyNz;if(fp_Vlcty=fopen("Vlcty.dat","wb")=NULL) printf("Conn

16、ot open file""Vlcty""");for(Ix=0;Ix<Nx;Ix+)for(Iz=0;Iz<Nz;Iz+) if(Iz<=3*Nz/4-1)/*VlctyIz=5000;/*else VlctyIz=5500;/*fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);/*fclose(fp_Vlcty);return 1;/*爆炸反射界面的生成，并調(diào)用FFT函數(shù)變換到頻率域儲存*/int WvFld0() FILE *fp_Rflct,*fp_Wfld0r,*fp_W

17、fld0i;int Ix,Iz,It;float RflctNz,Wfld0rNt,Wfld0iNt;if(fp_Wfld0r=fopen("Wfld0r.dat","wb")=NULL) printf("Connot open Wfld0r.dat");if(fp_Wfld0i=fopen("Wfld0i.dat","wb")=NULL) printf("Connot open Wfld0i.dat");if(fp_Rflct =fopen("Rflct.dat&

18、quot; ,"rb")=NULL) printf("Connot open Rflct.dat"); for(Ix=0;Ix<Nx;Ix+)printf("Wavefield0_FFT: Ix=%dn",Ix);for(Iz=0;Iz<Nz;Iz+) /賦值，爆炸反射界面t=0時(shí)刻起爆fread(&RflctIz,sizeof(RflctIz),1,fp_Rflct);for(It=0;It<Nt;It+)Wfld0rIt=0;Wfld0iIt=0;if(It=0) Wfld0rIt=RflctIz;/*i

19、f(kkfft(Wfld0r,Wfld0i,Nt,0)!=1) printf("FFT is error");exit(0);for(It=0;It<Nt/2+1;It+) /利用付氏變換的對稱性，保存一半的數(shù)據(jù)fwrite(&Wfld0rIt,sizeof(Wfld0rIt),1,fp_Wfld0r);/*fwrite(&Wfld0iIt,sizeof(Wfld0iIt),1,fp_Wfld0i);/*fclose(fp_Rflct);fclose(fp_Wfld0r);fclose(fp_Wfld0i);return 1;/*PhaseShift

20、Forward Modeling */int PsFrwd() /-波場相移延拓int PhaseShift( ); / Requset Function:PhaseShift調(diào)用波長函數(shù)相移延拓計(jì)算函數(shù)int Frqcy2Time( ); /調(diào)用波場做IFFT從頻率域變換到時(shí)間域函數(shù)if ( PhaseShift( ) !=1 ) printf("PhaseShift is errorn"); exit(0); / Call Functionif ( Frqcy2Time( ) !=1 ) printf("Frqcy2Time is errorn");

21、 exit(0); / Call Function return 1;int PhaseShift() / 1. Prepprocedure預(yù)處理FILE *fp_Wfldr,*fp_Wfldi,*fp_Wfld0r,*fp_Wfld0i,*fp_Vlcty,*fp_Absb;float kz2;int Ix,Ikx,Nkx=Nx,Iz,Iw,Nw=Nt;long Mgrtn;float VlctyNz;float AbsbNx,Wfld0rNx,Wfld0iNx,WfldrNx,WfldiNx;float Wfld_r,Wfld_i;float Kxmax,Dkx,Wmax,Dw;Wmax

22、= pai/0.004;Dw = Wmax/Nt;Kxmax = pai/20.;Dkx = Kxmax/Nx; / 2. Read in Velocity and Absorbing Boundary Date速度與削波數(shù)據(jù)讀入if(fp_Vlcty = fopen("Vlcty.dat","rb")=NULL)printf("Cann't open n");for(Iz=0;Iz<Nz;Iz+)fread(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);fclose(fp_Vlcty

23、);if(fp_Absb = fopen("Absb.dat","rb")=NULL)printf("Cann't open n");for(Ix=0;Ix<Nx;Ix+)fread(&AbsbIx,sizeof(AbsbIx),1,fp_Absb);fclose(fp_Absb);/ 3. Open Initial Wave Field Current Wave Field In Wave Fied Extrapolating波場文件打開 if(fp_Wfld0r = fopen("Wfld0r.da

24、t","rb")=NULL)printf("Cann't open n");if(fp_Wfld0i = fopen("Wfld0i.dat","rb")=NULL)printf("Cann't open n");if(fp_Wfldr = fopen("Wfldr.dat","wb") =NULL)printf("Cann't open n");if(fp_Wfldi = fopen("Wf

25、ldi.dat","wb") =NULL)printf("Cann't open n"); / 4. 每個(gè)頻率的波場延拓 for(Iw=0;Iw<Nw/2+1;Iw+) / 4.1初始化當(dāng)前波場for(Ix=0;Ix<Nx;Ix+)WfldrIx=0.;WfldiIx=0.; / 4.2波場從Iz=Nz-1最深處開始，延拓到Iz=1測線深度 for(Iz=Nz-1;Iz>0;Iz-) / 4.2.1形成新波場for(Ix=0;Ix<Nx;Ix+)/ 1. Take out Initial Wave Field D

26、ata With every Depth取出當(dāng)前深度的初始波場Mgrtn=(Ix*Nz+1+Iz)*(Nt/2+1)+Iw;fseek(fp_Wfld0r,sizeof(Wfld0rIx)*Mgrtn,SEEK_SET);fseek(fp_Wfld0i,sizeof(Wfld0iIx)*Mgrtn,SEEK_SET);fread(&Wfld0rIx,sizeof(Wfld0rIx),1,fp_Wfld0r);fread(&Wfld0iIx,sizeof(Wfld0iIx),1,fp_Wfld0i); / 2.新波場=初始波場+從下面延拓到此處的波場WfldrIx = Wfldr

27、Ix+Wfld0rIx;WfldiIx = WfldiIx+Wfld0iIx; / 3.邊界削波：新波場=新波場×削波因子WfldrIx = WfldrIx*AbsbIx;WfldiIx = WfldiIx*AbsbIx; / 4.2.2 新波場FFT到波數(shù)域if( kkfft(Wfldr,Wfldi,Nx,0) !=1 ) printf("FFT is errorn");exit(0); / 4.2.3頻率-波數(shù)域波場在從Iz+1延拓到Izfor(Ikx=0;Ikx<Nx/2+1;Ikx+) / 1.計(jì)算相移數(shù)據(jù)expikzdz（實(shí)部、虛部if( exp_

28、ikzDz(kz,Ix, (float)(VlctyIz/2.), Iw,Dw,Dkx) !=1) printf("exp_ikzDz is errorn");exit(0); ; / 2.波場延拓：波場=波場×相移數(shù)據(jù)Wfld_r = WfldrIkx*kz0-WfldiIkx*kz1;Wfld_i = WfldiIkx*kz0+WfldrIkx*kz1; WfldrIkx = Wfld_r;WfldiIkx = Wfld_i;if(Ikx!=0&&Ikx!=Nkx/2)Wfld_r = kz0*WfldrNkx-Ikx-kz1*WfldiNkx

29、-Ikx;Wfld_i = kz1*WfldrNkx-Ikx+kz0*WfldiNkx-Ikx;WfldrNkx-Ikx = Wfld_r;WfldiNkx-Ikx = Wfld_i; / 4.2.4 波場反FFT到空間域 if( kkfft(Wfldr,Wfldi,Nkx,1) !=1 ) printf("FFT is errorn");exit(0); / 4.3 存儲延拓到了測線的波場for(Ix=0;Ix<Nx;Ix+)fwrite(&WfldrIx,sizeof(WfldrIx),1,fp_Wfldr);fwrite(&WfldiIx,siz

30、eof(WfldiIx),1,fp_Wfldi); / 5.關(guān)閉文件，刪除中間文件。fclose(fp_Wfld0r);fclose(fp_Wfld0i);fclose(fp_Wfldr);fclose(fp_Wfldi);/remove("Absb.dat"); /Delete Absb.dat/remove("Rflct.dat"); /Delete Rflcy.dat/remove("Vlcty.dat"); /Delete Vlcty.dat /remove("Wfld0r.dat"); /Delete W

31、fld0r.dat/remove("Wfld0i.dat"); /Delete Wfld0i.datreturn(1);int exp_ikzDz(float eikzdz,int Ix,float Vc,int Iw,float Dw,float Dkx)float kz=0;eikzdz0=0;eikzdz1=0;kz=sqrt(pow(Iw*Dw/Vc,2)-pow(Ix*Dkx,2);if(kz>0)eikzdz0=(float)cos(kz*Dz);eikzdz1=(float)-sin(kz*Dz);return 1;FILE *fp_Wfldr,*fp_

32、Wfldi;FILE *fp_Record;int Ix,It,Iw,Nw=Nt;float WfldtrNt,WfldtiNt;long AddFrmStrt;if(fp_Wfldr = fopen("Wfldr.dat","rb") =NULL)printf("Cann't open ");exit(0);if(fp_Wfldi = fopen("Wfldi.dat","rb") =NULL)printf("Cann't open ");exit(0);i

33、f(fp_Record = fopen("Record.dat","wb") =NULL)printf("Cann't open ");exit(0);for(Ix=0;Ix<Nx;Ix+)for(Iw=0;Iw<Nw/2+1;Iw+)AddFrmStrt=Iw*Nx+Ix;fseek(fp_Wfldr,sizeof(WfldtrIw)*AddFrmStrt,SEEK_SET);fseek(fp_Wfldi,sizeof(WfldtiIw)*AddFrmStrt,SEEK_SET);fread(&Wfldt

34、rIw,sizeof(WfldtrIw),1,fp_Wfldr);fread(&WfldtiIw,sizeof(WfldtiIw),1,fp_Wfldi);if(Iw!=0&&Iw!=Nw/2)WfldtrNw-Iw = WfldtrIw;WfldtiNw-Iw = -WfldtiIw;if(kkfft(Wfldtr,Wfldti,Nw,1)!=1)printf("FFT is error");exit(0);for(It=0;It<Nt;It+) /按道取出數(shù)據(jù)，按Ix循環(huán)fwrite(&WfldtrIt,sizeof(WfldtrI

35、t),1,fp_Record); /按道存入數(shù)據(jù)，實(shí)部數(shù)據(jù)存入fclose(fp_Wfldr);fclose(fp_Wfldi);remove("Wfldr.dat");remove("Wfldi.dat");fclose(fp_Record);return 1;七、 實(shí)驗(yàn)結(jié)果1、 結(jié)果顯示改變繞射點(diǎn)位置觀察截圖： .改變速度觀察截圖 2、 對比分析削波后曲線更為平滑清晰改變繞射點(diǎn)位置時(shí)，曲線左右移動，移動趨勢同繞射點(diǎn)移動方向一致改變速度時(shí)，由公式可知，速度越大曲線越緩，其物理意義為，速度大走時(shí)短八、 討論建議1、 實(shí)驗(yàn)收獲大量上機(jī)實(shí)踐，編寫褶積函數(shù)的C

36、語言源程序有很大程度提高，在與同學(xué)交流討論過程中收獲很大，不管是學(xué)習(xí)還是工作要揚(yáng)長避短，但更要面對自己的不足，努力改善不足，完善各方面能力，尤其是編程能力。通過對比分析的實(shí)踐，學(xué)會了從差異中學(xué)習(xí)，在差異中進(jìn)行深入研究討論學(xué)習(xí)。2、存在問題編程能力有待提高，對專業(yè)理論只是有待鞏固加強(qiáng)。附：心得體會 其一，世上無難事只怕有心人，雖然對編程有懼，但是用心一步一步的進(jìn)行其實(shí)不是很難。完成本次試驗(yàn)后感覺編程有很大提到，得益于老師提供提供的原程序段，讓學(xué)生依葫蘆畫瓢，降低了學(xué)習(xí)難度。其二，功在平時(shí)，每次做好實(shí)驗(yàn)進(jìn)行存檔整理很重要，養(yǎng)成習(xí)慣不管是對學(xué)業(yè)、生活或是人生都是無形的財(cái)富。其三，學(xué)習(xí)可以和同學(xué)一起討

37、論但是深入研究需要自己完成，其中的收獲只有自己慢慢的領(lǐng)悟。畢業(yè)設(shè)計(jì)（論文）原創(chuàng)性聲明和使用授權(quán)說明原創(chuàng)性聲明本人鄭重承諾：所呈交的畢業(yè)設(shè)計(jì)（論文），是我個(gè)人在指導(dǎo)教師的指導(dǎo)下進(jìn)行的研究工作及取得的成果。盡我所知，除文中特別加以標(biāo)注和致謝的地方外，不包含其他人或組織已經(jīng)發(fā)表或公布過的研究成果，也不包含我為獲得 及其它教育機(jī)構(gòu)的學(xué)位或?qū)W歷而使用過的材料。對本研究提供過幫助和做出過貢獻(xiàn)的個(gè)人或集體，均已在文中作了明確的說明并表示了謝意。作 者 簽 名： 日 期： 指導(dǎo)教師簽名： 日期： 使用授權(quán)說明本人完全了解 大學(xué)關(guān)于收集、保存、使用畢業(yè)設(shè)計(jì)（論文）的規(guī)定，即：按照學(xué)校要求提交畢業(yè)設(shè)計(jì)（論文）的印

38、刷本和電子版本；學(xué)校有權(quán)保存畢業(yè)設(shè)計(jì)（論文）的印刷本和電子版，并提供目錄檢索與閱覽服務(wù)；學(xué)?？梢圆捎糜坝　⒖s印、數(shù)字化或其它復(fù)制手段保存論文；在不以贏利為目的前提下，學(xué)校可以公布論文的部分或全部內(nèi)容。作者簽名： 日 期： 學(xué)位論文原創(chuàng)性聲明本人鄭重聲明：所呈交的論文是本人在導(dǎo)師的指導(dǎo)下獨(dú)立進(jìn)行研究所取得的研究成果。除了文中特別加以標(biāo)注引用的內(nèi)容外，本論文不包含任何其他個(gè)人或集體已經(jīng)發(fā)表或撰寫的成果作品。對本文的研究做出重要貢獻(xiàn)的個(gè)人和集體，均已在文中以明確方式標(biāo)明。本人完全意識到本聲明的法律后果由本人承擔(dān)。作者簽名： 日期： 年 月 日學(xué)位論文版權(quán)使用授權(quán)書本學(xué)位論文作者完全了解學(xué)校有關(guān)保留、使用學(xué)位論文的規(guī)定，同意學(xué)校保留并向