地震數值模擬實驗報告

1、冷孫強z去專本科生實驗報告實驗課程數值模型模擬學院名稱地球物理學院專業(yè)名稱勘查技術與工程學生姓名ZRY學生學號指導教師實驗地點624實驗成績二O五年4月二O五年5月成都理工大學地震數值模擬實驗報告實驗時間2015年5月31日開課單位地球物理學院指導教師實驗題目：疊加地震記錄的相移波動方程正演模擬實驗姓名學號班級專業(yè)勘查技術與工程院（系）地球物理學院單內容理解寫作結構項程序設計成模型設計計算結果績結果分析總成績實驗二疊加地震記錄的相移波動模擬實方程正演驗摘要利用C語言編制地質模型的相移波動方程正演模擬，改變繞射點位置、速度，再做正演模擬。關鍵字：地震模型；正演記錄實驗目的掌握各向同性介質任意構造

2、、水平層狀速度結構地質模型的相移波動方程正演模擬基本理論、實現方法與程序編制，由正演記錄初步分析地震信號的分辨率。實驗內容1、基本要求：（1）點繞射構造和水平層狀速度模型（參數如圖1所示）的正演數值模擬；1）削波的正演；2）無削波的震正演；（2）計算中點和兩個邊界的信號位置，分析實驗結果的正確性；（3）做同樣模型的褶積模型數值模擬，對比分析分析兩者的異同。（4）改變繞射點位置、速度，再做正演模擬。2、較高要求：（1）使用雷克子波做爆炸源，對三個不同的主頻：25hz、50hz和75hz分別做點繞射模型的正演模擬；（2）設計復雜反射構造模型，再做正演模擬。實驗原理1、地震波傳播的波動方程設（x,z

3、）為空間坐標，t為時間，地震波傳播速度為v（x,z），則二位介質中任意位置、任意時刻的地震波場為p（z,x,t）：壓縮波一一波。則二維各向同性均勻介質中地震波傳播的遵循聲波方程為（）2、傅里葉變換的微分性質p（t）與其傅里葉變換的P（w）的關系:則有時間微分性質w為頻率，w=2k/T,T為周期。同理有空間微分性質:k為頻率，k=2兀/X，X為波長。3、地震波傳播的相移外推公式令速度v不隨x變化，只隨z變化，則利用傅里葉變換微分性質（3）和（4）式，把波動方程（1）式變換到頻率波數域，得:或:令:則（5）式的解為:包括上行波和下行波兩項。正演模擬取上行波：若和間隔為，速度v（z）在此間隔內不隨Z

4、變的常數,（7）式實現波場從到的延拓，即：在深度zj+i開始向上延拓到z，若延拓深度為零，即：az=z+i-Zj=o，貝y對于任意深度Zj+1到Zj的延拓，可得正演模擬中地震波的傳播方程（延拓公式）4、初始條件和邊界條件按照爆炸界面理論，反射界面震源在t=0時刻同時起爆，此時刻的波場就是震源。根據不同情況，可直接使用反射系數脈沖或子波作震源。如果直接使用反射系數作震源脈沖，貝初始條件可表示為：其他對時間Z和空間X做二維傅立葉變換，則得頻率-波數域的初始波場邊界條件：其他，其他，其他其他參數都是在范圍內定義的。5、邊界處理（1）邊界反射問題把實際無窮空間區(qū)域中求解波場的問題化為有窮區(qū)域求解時，左

5、右兩邊使用零邊界條件物理上假設探區(qū)距和兩個端點很遠，在兩個端點上收到的反射波很弱。但是，上述條件在實際中不能成立，造成零邊界條件反而成為絕對阻止波通過的強反射面。在正演模擬的剖面上出現了邊界假反射干涉正常界面的反射。（2）邊界強反射的處理鑲邊法、削波法、吸收邊界都能有效消除邊界強反射。削波法就是在波場延拓過程中，沒延拓一次，在其兩側均勻衰減到零，從而消除邊界強反射的影響。假設橫向總長度為NX,以兩邊Lx道吸波為例，有以下吸波公式:其他6、數字化根據數字信號處理的采樣定理，把連續(xù)的信號變?yōu)橛嬎銠C能處理的數字信號，使相移法正演模擬得以實現。頻域抽樣定理：一個頻譜受限信號,如果時間只占據-tm+tm

6、的范圍，若在頻域以不大于1/2tm頻率間隔華1/2tm對信號f(t)的頻譜F()采樣，則抽樣到的離散信號F可以唯一表示原信號。時域抽樣定理：一個時間受限信號f(t)，如果頻譜只占據-m+%的范圍，則信號(t)可以用等間隔的抽樣值唯一表示出來，而時間At抽樣間隔必須不大于1%,m=2兀fm，g/2fm。實驗過程/#includePsFrwrdMdlParameter.h#include#include#include#include#includeintkkfft(floatpr,floatpi,intn,intl);intAbsorb(intLabs);intPo2Judge(int);int

7、Rflct();intVlcty();intWvFld0();intPsFrwd();intexp_ikzDz(floateikzdz,intIx,floatVc,intIw,floatDw,floatDkx);intPhaseShift();/相移intFrqcy2Time();#defineNx128/TraceNumber#defineNt256/RecordNumber#defineNz100/DepthNumber#defineDx20./TraceInterval#defineDt0.004/RecordInterval#defineDz20./DepthInterval#defi

8、nePI3.1415voidmain()intLabs=10;/LengthOfBoundaryAbsorbingif(Po2Judge(Nt)!=1)printf(Nt=%disnotthePowerof2n,Nt);exit(0);elseprintf(Nt=%disthePowerof2n,Nt);if(Po2Judge(Nx)!=1)printf(Nx=%disnotthePowerof2n,Nx);exit(0);elseprintf(Nx=%disthePowerof2n,Nx);if(Absorb(Labs)!=1)printf(Absorbiserrorn);exit(0);e

9、lseprintf(Absorbisokayn);if(Rflct()!=1)printf(Reflectioniserrorn);exit(0);elseprintf(Reflectionisokayn);if(Vlcty()!=1)printf(Vlctyiserrorn);exit(0);elseprintf(Vlctyisokayn);if(WvFld0()!=1)printf(WvFldiserrorn);exit(0);elseprintf(WvFldisokayn);if(PsFrwd()!=1)printf(PsFurdiserrorn);exit(0);elseprintf(

10、PsFurdisokayn);remove(Wfld0r.dat);remove(Wfld0i.dat);remove(PsFrwrdMdl.ncb);remove(PsFrwrdMdl.dsp);remove(PsFrwrdMdl);remove(Abs);/return();/intPo2Judge(intN)/JudgeNtorNxisPowerof2doublea=0;for(inti=0;a=1)return(0);elsereturn(1);/intAbsorb(intLabs)/FormingAbsorbingBoudaryFILE*fp_Absorb,*fp_Abs;intIx

11、;floatAbsNx;if(fp_Absorb=fopen(Absb.dat,wb)=NULL)printf(CannotopenfileAbsorb);elseprintf(Absbisokayn);if(fp_Abs=fopen(Abs.txt,w)=NULL)printf(CannotopenfileAbs);elseprintf(Absisokayn);for(Ix=0;IxNx;Ix+)AbsIx=1.;for(Ix=0;IxLabs;Ix+)AbsIx=sqrt(sin(PI*Ix/(2*(Labs-1);AbsNx-Ix-1=AbsIx;for(Ix=0;IxNx;Ix+)fw

12、rite(&AbsIx,sizeof(AbsIx),Nx,fp_Absorb);/printf(%fn,AbsIx);fprintf(fp_Abs,%fn,AbsIx);fclose(fp_Absorb);fclose(fp_Abs);return(1);/intRflct()/FormingReflectStructureModelFILE*fp_Rflct;intIx,Iz;floatRflctNz;if(fp_Rflct=fopen(Rflct.dat,wb)=NULL)printf(CannotopenfileReflection);for(Ix=0;IxNx;Ix+)for(Iz=0

13、;IzNz;Iz+)RflctIz=0.0;if(Ix=(int)(Nx/2)&Iz=(int)(Nz/2)RflctIz=1.;fwrite(&RflctIz,sizeof(RflctIz),1,fp_Rflct);fclose(fp_Rflct);return(1);/intVlcty()/FormingVelociyModelFILE*fp_Vlcty;intIx,Iz;floatVlctyNz;if(fp_Vlcty=fopen(Vlcty.dat,wb)=NULL)printf(CannotopenfileVlcty);for(Ix=0;IxNx;Ix+)for(Iz=0;Iz2*N

14、z/3;Iz+)VlctyIz=3000.;/fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for(Iz=(2*Nz/3);IzNz;Iz+)VlctyIz=6000.;/fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);fwrite(&Vlcty0,sizeof(VlctyIz),Nz,fp_Vlcty);fclose(fp_Vlcty);return(1);/intWvFld0()/WaveFieldInitializationFILE*fp_Rflct,*fp_Wfld0r,*fp_Wfld0i,*fp_Wf

15、ld0isee,*fp_Wfld0rsee,*fp_Wfld0r0see;intIx,Iz,It;floatRflctNz;floatWfld0rNt,Wfld0iNt;if(fp_Wfld0r=fopen(Wfld0r.dat,wb)=NULL)printf(CannotopenWfld0r.dat);if(fp_Wfld0i=fopen(Wfld0i.dat,wb)=NULL)printf(CannotopenWfld0i.dat);if(fp_Rflct=fopen(Rflct.dat,rb)=NULL)printf(CannotopenRflct.dat);/if(fp_Wfld0is

16、ee=fopen(fp_Wfld0isee.txt,wb)=NULL)printf(Cannotopenfp_Wfld0isee.dat);/if(fp_Wfld0rsee=fopen(fp_Wfld0rsee.txt,wb)=NULL)printf(Cannotopenfp_Wfld0rsee.dat);/if(fp_Wfld0r0see=fopen(fp_Wfld0r0see.txt,wb)=NULL)printf(Cannotopenfp_Wfld0r0see.dat);for(Ix=0;IxNx;Ix+)/printf(Wavefield0_FFT:Ix=%dn,Ix);for(Iz=

17、0;IzNz;Iz+)fread(&RflctIz,sizeof(RflctIz),1,fp_Rflct);for(It=0;ItNt;It+)Wfld0rIt=0.;Wfld0iIt=0.;if(It=0)Wfld0rIt=RflctIz;/fprintf(fp_Wfld0r0see,%ft,Wfld0rIz);if(kkfft(Wfld0r,Wfld0i,Nt,0)!=1)printf(FFTiserror);exit(0);/原為時間域函數，變后為頻率域。for(It=0;It(int)(Nt/2)+1;It+)fwrite(&Wfld0rIt,sizeof(Wfld0rIt),1,fp

18、_Wfld0r);fwrite(&Wfld0iIt,sizeof(Wfld0iIt),1,fp_Wfld0i);/fprintf(fp_Wfld0isee,%ft,Wfld0iIt);/fprintf(fp_Wfld0rsee,%ft,Wfld0rIt);/fprintf(fp_Wfld0isee,n);/fprintf(fp_Wfld0rsee,n);/fprintf(fp_Wfld0r0see,n);fclose(fp_Rflct);fclose(fp_Wfld0r);fclose(fp_Wfld0i);return(1);/intPsFrwd()/WaveFieldPhaseShifti

19、ntPhaseShift();intFrqcy2Time();if(PhaseShift()!=1)printf(PhaseShiftiserrorn);exit(0);if(Frqcy2Time()!=1)printf(Frqcy2Timeiserrorn);exit(0);return(1);/intPhaseShift()/相移/1.Prepprocedure預處理FILE*fp_Wfldr,*fp_Wfldi,*fp_Wfld0r,*fp_Wfld0i,*fp_Vlcty,*fp_Absb;floatkz2;intIx,Iz,Iw,Nw=Nt,Ikx,Nkx=Nx;longMgrtn;

20、floatVlctyNz;floatWfld_r,Wfld_i;floatAbsbNx,Wfld0rNx,Wfld0iNx,WfldrNx,WfldiNx;floatKxmax,Dkx,Wmax,Dw;Wmax=PI/Dt;Dw=Wmax/(float)Nt;Kxmax=PI/Dx;Dkx=Kxmax/(float)Nx;open/2.ReadinVelocityandAbsorbingBoundaryData速度與削波數據讀入if(fp_Vlcty=fopen(Vlcty.dat,rb)=NULL)printf(CannotRflct.dat);for(Iz=0;IzNz;Iz+)fread

21、(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);fclose(fp_Vlcty);if(fp_Absb=fopen(Absb.dat,rb)=NULL)printf(CannotopenAbsb.dat);for(Ix=0;IxNx;Ix+)fread(&AbsbIx,sizeof(AbsbIx),1,fp_Absb);fclose(fp_Absb);/remove(Absb.dat);/3.OpenInitialWaveFieldFileandCurrentWaveFieldFileusingInWaveFieldExtrapolating波場文件打開if(fp_

22、Wfld0r=fopen(Wfld0r.dat,rb)=NULL)printf(CannotopenWfld0r.dat);if(fp_Wfld0i=fopen(Wfld0i.dat,rb)=NULL)printf(CannotopenWfld0i.dat);if(fp_Wfldr=fopen(Wfldr.dat,wb)=NULL)printf(CannotopenWfldr.dat);if(fp_Wfldi=fopen(Wfldi.dat,wb)=NULL)printf(CannotopenWfldi.dat);/4.WaveFiedExtrapolateWithEveryFrequency

23、每個頻率的波場延拓for(Iw=0;IwNw/2+1;Iw+)printf(PhaseShift:Iw=%dn,Iw);/4.1InitializingWavefieldofEveryFrequency:Allget0初始化當前波場for(Ix=0;Ix0;Iz-)/4.2.1FormingNewWavefieldforExtrapolatingonSpaceDomain形成新波長for(Ix=0;IxNx;Ix+)Mgrtn=(Ix*Nz+Iz)*(Nw/2+1)+Iw;/TakeoutInitialWaveFieldDataWitheveryDepth取出當前深度的初始波場fseek(fp

24、_Wfld0r,sizeof(Wfld0rIx)*Mgrtn,0);fread(&Wfld0rIx,sizeof(Wfld0rIx),1,fp_Wfld0r);fseek(fp_Wfld0i,sizeof(Wfld0iIx)*Mgrtn,0);fread(&Wfld0iIx,sizeof(Wfld0iIx),1,fp_Wfld0i);/NewWaveFiedEqualtotheInitialonPlusCurrentone新波場二初始波場+從下面延拓到此處的波場WfldrIx=WfldrIx+Wfld0rIx;WfldiIx=WfldiIx+Wfld0iIx;/Boundaryabsorbin

25、gofNewWaveFied邊界削波：新波場二新波場X削波因子WfldrIx=WfldrIx*AbsbIx;WfldiIx=WfldiIx*AbsbIx;/4.2.2FFTofNewWavefieldFromSpacetoWaveNumber新波場FFT到波數域if(kkfft(Wfldr,Wfldi,Nx,0)!=1)printf(FFTiserror);exit(0);/4.2.3WavefieldExtrapolatingonFrequency-WaveNumberDomain頻率-波數域波場在從lz+1延拓到Izfor(Ikx=0;IkxNx/2+1;Ikx+)/4.2.3.1Comp

26、utingPhaseshiftdata計算相移數據expikzdz(實部、虛部)if(exp_ikzDz(kz,Ikx,(float)(VlctyIz/2.),Iw,Dw,Dkx)!=1)printf(exp_ikzDziserror);exit(0);/4.2.3.2WaveFieldmultiplyPhaseshiftdata波場延拓：波場二波場X相移數據Wfld_r=WfldrIkx*kz0-WfldiIkx*kz1;Wfld_i=WfldiIkx*kz0+WfldrIkx*kz1;WfldrIkx=Wfld_r;WfldiIkx=Wfld_i;if(Ikx!=0&Ikx!=Nkx/2)

27、Wfld_r=kz0*WfldrNkx-Ikx-kz1*WfldiNkx-Ikx;Wfld_i=kz1*WfldrNkx-Ikx+kz0*WfldiNkx-Ikx;WfldrNkx-Ikx=Wfld_r;WfldiNkx-Ikx=Wfld_i;/4.2.4IFFTofNewWavefieldFromWaveNumbertoSpaceDomain波場反FFT到空間域if(kkfft(Wfldr,Wfldi,Nkx,1)!=1)printf(FFTiserror);exit(0);/4.3StoringWavefieldDataonSurveyLine:Z二Zmin存儲延拓到了測線的波場for(I

28、x=0;Ix0.)eikzdz0=(float)cos(kz*Dz);eikzdz1=(float)-sin(kz*Dz);return(1);/intFrqcy2Time()FILE*fp_Wfldr,*fp_Wfldi,*fp_Record;intIx,It,Iw,Nw=Nt;floatWfldtrNt,WfldtiNt;longAddFrmStrt;if(fp_Wfldr=fopen(Wfldr.dat,rb)=NULL)printf(CanntWfldr.datn);if(fp_Wfldi=fopen(Wfldi.dat,rb)=NULL)printf(CanntWfldi.datn)

29、;if(fp_Record=fopen(Record.dat,wb)=NULL)printf(CanntRecord.datn);for(Ix=0;IxNx;Ix+)/printf(Frqcy2Time:Ix=%dn,Ix);for(Iw=0;IwNw/2+1;Iw+)AddFrmStrt=Iw*Nx+Ix;fseek(fp_Wfldr,sizeof(WfldtrIw)*AddFrmStrt,0);fseek(fp_Wfldi,sizeof(WfldtiIw)*AddFrmStrt,0);fread(&WfldtrIw,sizeof(WfldtrIw),1,fp_Wfldr);fread(&W

30、fldtiIw,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(FFTiserrorn);exit(0);for(It=0;It0;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.283