地震数值模拟实验报告

1、冷孙强z去專本科生实验报告实验课程数值模型模拟学院名称地球物理学院专业名称勘查技术与工程学生姓名ZRY学生学号指导教师实验地点624实验成绩二O五年4月二O五年5月成都理工大学地震数值模拟实验报告实验时间2015年5月31日开课单位地球物理学院指导教师实验题目：叠加地震记录的相移波动方程正演模拟实验姓名学号班级专业勘查技术与工程院（系）地球物理学院单内容理解写作结构项程序设计成模型设计计算结果绩结果分析总成绩实验二叠加地震记录的相移波动模拟实方程正演验摘要利用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）边界反射问题把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左

5、右两边使用零边界条件物理上假设探区距和两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。（2）边界强反射的处理镶边法、削波法、吸收边界都能有效消除边界强反射。削波法就是在波场延拓过程中，没延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX,以两边Lx道吸波为例，有以下吸波公式:其他6、数字化根据数字信号处理的采样定理，把连续的信号变为计算机能处理的数字信号，使相移法正演模拟得以实现。频域抽样定理：一个频谱受限信号,如果时间只占据-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