地震數(shù)值模擬實(shí)驗(yàn)報(bào)告

1、冷孫強(qiáng)z去專本科生實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)課程數(shù)值模型模擬學(xué)院名稱地球物理學(xué)院專業(yè)名稱勘查技術(shù)與工程學(xué)生姓名ZRY學(xué)生學(xué)號(hào)指導(dǎo)教師實(shí)驗(yàn)地點(diǎn)624實(shí)驗(yàn)成績(jī)二O五年4月二O五年5月成都理工大學(xué)地震數(shù)值模擬實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)時(shí)間2015年5月31日開(kāi)課單位地球物理學(xué)院指導(dǎo)教師實(shí)驗(yàn)題目：疊加地震記錄的相移波動(dòng)方程正演模擬實(shí)驗(yàn)姓名學(xué)號(hào)班級(jí)專業(yè)勘查技術(shù)與工程院（系）地球物理學(xué)院?jiǎn)蝺?nèi)容理解寫作結(jié)構(gòu)項(xiàng)程序設(shè)計(jì)成模型設(shè)計(jì)計(jì)算結(jié)果績(jī)結(jié)果分析總成績(jī)實(shí)驗(yàn)二疊加地震記錄的相移波動(dòng)模擬實(shí)方程正演驗(yàn)摘要利用C語(yǔ)言編制地質(zhì)模型的相移波動(dòng)方程正演模擬，改變繞射點(diǎn)位置、速度，再做正演模擬。關(guān)鍵字：地震模型；正演記錄實(shí)驗(yàn)?zāi)康恼莆崭飨蛲越橘|(zhì)任意構(gòu)造

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

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

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

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

6、的范圍，若在頻域以不大于1/2tm頻率間隔華1/2tm對(duì)信號(hào)f(t)的頻譜F()采樣，則抽樣到的離散信號(hào)F可以唯一表示原信號(hào)。時(shí)域抽樣定理：一個(gè)時(shí)間受限信號(hào)f(t)，如果頻譜只占據(jù)-m+%的范圍，則信號(hào)(t)可以用等間隔的抽樣值唯一表示出來(lái)，而時(shí)間At抽樣間隔必須不大于1%,m=2兀fm，g/2fm。實(shí)驗(yàn)過(guò)程/#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);/原為時(shí)間域函數(shù)，變后為頻率域。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預(yù)處理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速度與削波數(shù)據(jù)讀入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波場(chǎng)文件打開(kāi)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、每個(gè)頻率的波場(chǎng)延拓for(Iw=0;IwNw/2+1;Iw+)printf(PhaseShift:Iw=%dn,Iw);/4.1InitializingWavefieldofEveryFrequency:Allget0初始化當(dāng)前波場(chǎng)for(Ix=0;Ix0;Iz-)/4.2.1FormingNewWavefieldforExtrapolatingonSpaceDomain形成新波長(zhǎng)for(Ix=0;IxNx;Ix+)Mgrtn=(Ix*Nz+Iz)*(Nw/2+1)+Iw;/TakeoutInitialWaveFieldDataWitheveryDepth取出當(dāng)前深度的初始波場(chǎng)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新波場(chǎng)二初始波場(chǎng)+從下面延拓到此處的波場(chǎng)WfldrIx=WfldrIx+Wfld0rIx;WfldiIx=WfldiIx+Wfld0iIx;/Boundaryabsorbin

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

26、utingPhaseshiftdata計(jì)算相移數(shù)據(jù)expikzdz(實(shí)部、虛部)if(exp_ikzDz(kz,Ikx,(float)(VlctyIz/2.),Iw,Dw,Dkx)!=1)printf(exp_ikzDziserror);exit(0);/4.2.3.2WaveFieldmultiplyPhaseshiftdata波場(chǎng)延拓：波場(chǎng)二波場(chǎng)X相移數(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)

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波場(chǎng)反FFT到空間域if(kkfft(Wfldr,Wfldi,Nkx,1)!=1)printf(FFTiserror);exit(0);/4.3StoringWavefieldDataonSurveyLine:Z二Zmin存儲(chǔ)延拓到了測(cè)線的波場(chǎng)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