1、地震数值模拟实验报告本科生实验报告实验课程 数值模型模拟 学院名称 地球物理学院 专业名称 勘查技术与工程 学生姓名 ZRY 学生学号 指导教师 实验地点 624 实验成绩二一五年4月 二一五年5月成都理工大学地震数值模拟实验报告实验时间2015年5月31日开课单位地球物理学院指导教师实验题目：叠加地震记录的相移波动方程正演模拟实验姓名学号班级专业勘查技术与工程院（系）地球物理学院单项成绩内容理解写作结构程序设计模型设计计算结果结果分析总成绩实验二 叠加地震记录的相移波动模拟实方程正演验摘要利用C语言编制地质模型的相移波动方程正演模拟，改变绕射点位置、速度，再做正演模拟。关键字：地震模型；正演
2、记录1.1实验目的掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论、实现方法与程序编制，由正演记录初步分析地震信号的分辨率。1.2实验内容1、 基本要求：（1） 点绕射构造和水平层状速度模型（参数如图1 所示）的正演数值模拟；1） 削波的正演；2） 无削波的震正演；（2） 计算中点和两个边界的信号位置，分析实验结果的正确性；（3） 做同样模型的褶积模型数值模拟，对比分析分析两者的异同。（4） 改变绕射点位置、速度，再做正演模拟。2、 较高要求：（1）使用雷克子波做爆炸源，对三个不同的主频：25hz、50hz 和75hz 分别做点绕射模型的正演模拟；（2）设计复杂
3、反射构造模型，再做正演模拟。1.3实验原理1、地震波传播的波动方程设（x,z）为空间坐标，t 为时间，地震波传播速度为v(x,z)，则二位介质中任意位置、任意时刻的地震波场为 p(z,x,t)：压缩波纵波。则二维各向同性均匀介质中地震波传播的遵循声波方程为2、傅里叶变换的微分性质p(t)与其傅里叶变换的P()的关系：则有时间微分性质 为频率，=2/T，T 为周期。同理有空间微分性质：k 为频率，k=2/， 为波长。3、地震波传播的相移外推公式令速度v 不随x 变化，只随z 变化，则利用傅里叶变换微分性质（3）和（4）式，把波动方程(1)式变换到频率-波数域，得：或：令：则（5）式的解为：包括上
4、行波和下行波两项。正演模拟取上行波：若和间隔为 Z ，速度v(z) 在此间隔内不随Z 变的常数,(7)式实现波场从到的延拓，即：在深度Zj+1 开始向上延拓到Zj，若延拓深度为零，即：Z= Zj+1-Zj=0，则对于任意深度Zj+1 到Zj 的延拓，可得正演模拟中地震波的传播方程（延拓公式）4、初始条件和边界条件按照爆炸界面理论，反射界面震源在 t=0 时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为： 对时间t和空间x做二维傅立叶变换，则得频率-波数域的初始波场边界条件：其他参数都是在范围内定义的。5、边
5、界处理（1）边界反射问题把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。（2）边界强反射的处理镶边法、削波法、吸收边界都能有效消除边界强反射。削波法就是在波场延拓过程中，没延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx 道吸波为例，有以下吸波公式：6、数字化根据数字信号处理的采样定理，把连续的信号变为计算机能处理的数字信号，使相移
6、法正演模拟得以实现。频域抽样定理：一个频谱受限信号f(t)，如果时间只占据 -tm +tm 的范围，若在频域以不大于1/2tm 频率间隔 f1/2tm 对信号f(t)的频谱F()采样，则抽样到的离散信号F1()可以唯一表示原信号。时域抽样定理：一个时间受限信号f(t)，如果频谱只占据 -m+m 的范围，则信号f(t)可以用等间隔的抽样值唯一 表示出来，而时间t 抽样间隔必须不大于1/2fm，m =2 fm，t1/2fm。1.4实验过程/#include PsFrwrdMdlParameter.h#include #include #include #include #include int k
7、kfft(float pr, float pi, int n, int l);int Absorb(int Labs);int Po2Judge(int);int Rflct();int Vlcty();int WvFld0();int PsFrwd();int exp_ikzDz(float eikzdz,int Ix, float Vc,int Iw,float Dw,float Dkx);int PhaseShift();/相移int Frqcy2Time();#define Nx 128 /Trace Number#define Nt 256 /Record Number#define
8、 Nz 100 /Depth Number#define Dx 20. /Trace Interval#define Dt 0.004/Record Interval#define Dz 20. /Depth Interval#define PI 3.1415void main() int Labs=10;/Length Of Boundary Absorbing if(Po2Judge(Nt) !=1)printf(Nt=%d is not the Power of 2n,Nt);exit(0); else printf(Nt=%d is the Power of 2n,Nt); if(Po
9、2Judge(Nx) !=1)printf(Nx=%d is not the Power of 2n,Nx);exit(0); else printf(Nx=%d is the Power of 2n,Nx); if(Absorb(Labs) !=1)printf(Absorb is errorn); exit(0); else printf(Absorb is okayn); if(Rflct() !=1)printf(Reflection is errorn); exit(0); else printf(Reflection is okayn); if(Vlcty() !=1)printf
10、(Vlcty is errorn); exit(0); else printf(Vlcty is okayn); if(WvFld0() !=1)printf(WvFld is errorn); exit(0); else printf(WvFld is okayn); if(PsFrwd() !=1)printf(PsFurd is errorn); exit(0); else printf(PsFurd is okayn); remove(Wfld0r.dat); remove(Wfld0i.dat); remove(PsFrwrdMdl.ncb); remove(PsFrwrdMdl.d
11、sp); remove(PsFrwrdMdl); remove(Abs);/ return( );/int Po2Judge(int N)/Judge Nt or Nx is Power of 2 double a=0; for(int i=0;a=1)return(0); else return(1); /int Absorb(int Labs) /Forming Absorbing Boudary FILE *fp_Absorb,*fp_Abs; int Ix; float AbsNx; if(fp_Absorb=fopen(Absb.dat,wb)=NULL) printf(Cannot
12、 open file Absorb); else printf(Absb is okayn); if(fp_Abs=fopen(Abs.txt,w)=NULL) printf(Cannot open file Abs); else printf(Abs is okayn); 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+) fwrite(&AbsIx,sizeof(AbsIx),Nx,fp_Absorb)
13、; /printf(%fn,AbsIx); fprintf(fp_Abs,%fn,AbsIx); fclose(fp_Absorb); fclose(fp_Abs); return(1); /int Rflct()/Forming Reflect Structure Model FILE *fp_Rflct; int Ix,Iz; float RflctNz; if(fp_Rflct=fopen(Rflct.dat,wb)=NULL) printf(Cannot open fileReflection); for(Ix=0;IxNx;Ix+) for(Iz=0;IzNz;Iz+) RflctI
14、z=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);/int Vlcty() /Forming Velociy Model FILE *fp_Vlcty; int Ix,Iz; float VlctyNz; if(fp_Vlcty=fopen(Vlcty.dat,wb)=NULL)printf(Cannot open file Vlcty); for(Ix=0;IxNx;Ix+) for(Iz=0;
15、Iz2*Nz/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); /int WvFld0()/Wave Field Initialization FILE *fp_Rflct,*fp_Wfld0
16、r,*fp_Wfld0i,*fp_Wfld0isee,*fp_Wfld0rsee,*fp_Wfld0r0see; int Ix,Iz,It; float RflctNz; float Wfld0rNt,Wfld0iNt; if(fp_Wfld0r=fopen(Wfld0r.dat,wb)=NULL)printf(Cannot open Wfld0r.dat); if(fp_Wfld0i=fopen(Wfld0i.dat,wb)=NULL)printf(Cannot open Wfld0i.dat); if(fp_Rflct=fopen(Rflct.dat,rb)=NULL)printf(Can
17、not open Rflct.dat); /if(fp_Wfld0isee=fopen(fp_Wfld0isee.txt,wb)=NULL)printf(Cannot open fp_Wfld0isee.dat); /if(fp_Wfld0rsee=fopen(fp_Wfld0rsee.txt,wb)=NULL)printf(Cannot open fp_Wfld0rsee.dat); /if(fp_Wfld0r0see=fopen(fp_Wfld0r0see.txt,wb)=NULL)printf(Cannot open fp_Wfld0r0see.dat); for(Ix=0;IxNx;I
18、x+) /printf(Wavefield0_FFT:Ix=%dn,Ix); for(Iz=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,%f t,Wfld0rIz); if(kkfft(Wfld0r,Wfld0i,Nt,0)!=1)printf(FFT is error);exit(0);/原为时间域函数，变后为频率域。 for(It=0;
19、It(int)(Nt/2)+1;It+) fwrite(&Wfld0rIt,sizeof(Wfld0rIt),1,fp_Wfld0r); fwrite(&Wfld0iIt,sizeof(Wfld0iIt),1,fp_Wfld0i); /fprintf(fp_Wfld0isee,%f t,Wfld0iIt); /fprintf(fp_Wfld0rsee,%f t,Wfld0rIt); /fprintf(fp_Wfld0isee,n); /fprintf(fp_Wfld0rsee,n); /fprintf(fp_Wfld0r0see,n); fclose(fp_Rflct); fclose(fp_
20、Wfld0r); fclose(fp_Wfld0i); return(1);/int PsFrwd()/Wave Field PhaseShift int PhaseShift(); int Frqcy2Time(); if( PhaseShift()!=1 )printf(PhaseShift is errorn); exit(0); if( Frqcy2Time()!=1 )printf(Frqcy2Time is errorn); exit(0); return(1); /int PhaseShift()/相移 /1.Prepprocedure预处理 FILE *fp_Wfldr,*fp
21、_Wfldi,*fp_Wfld0r,*fp_Wfld0i,*fp_Vlcty,*fp_Absb; float kz2; int Ix,Iz,Iw,Nw=Nt,Ikx,Nkx=Nx; long Mgrtn; float VlctyNz; float Wfld_r,Wfld_i; float AbsbNx,Wfld0rNx,Wfld0iNx,WfldrNx,WfldiNx; float Kxmax,Dkx,Wmax,Dw; Wmax=PI/Dt; Dw=Wmax/(float)Nt; Kxmax=PI/Dx; Dkx=Kxmax/(float)Nx; /2.Read in Velocity and
22、 Absorbing Boundary Data速度与削波数据读入 if(fp_Vlcty=fopen(Vlcty.dat,rb)=NULL) printf(Cannot open Rflct.dat); for(Iz=0;IzNz;Iz+) fread(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty); fclose(fp_Vlcty); if(fp_Absb=fopen(Absb.dat,rb)=NULL)printf(Cannot open Absb.dat); for(Ix=0;IxNx;Ix+) fread(&AbsbIx,sizeof(AbsbIx),1,f
23、p_Absb); fclose(fp_Absb); /remove(Absb.dat); /3.Open Initial Wave Field File and Current Wave Field File using In Wave Field Extrapolating波场文件打开 if(fp_Wfld0r=fopen(Wfld0r.dat,rb)=NULL) printf(Cannot open Wfld0r.dat); if(fp_Wfld0i=fopen(Wfld0i.dat,rb)=NULL) printf(Cannot open Wfld0i.dat); if(fp_Wfldr
24、=fopen(Wfldr.dat,wb)=NULL) printf(Cannot open Wfldr.dat); if(fp_Wfldi=fopen(Wfldi.dat,wb)=NULL) printf(Cannot open Wfldi.dat); /4.Wave Fied Extrapolate With Every Frequency每个频率的波场延拓 for(Iw=0;IwNw/2+1;Iw+) printf(PhaseShift:Iw=%dn,Iw); /4.1 Initializing Wavefield of Every Frequency: All get 0初始化当前波场
25、for(Ix=0;Ix0;Iz-) /4.2.1 Forming New Wavefield for Extrapolatingon Space Domain形成新波长 for(Ix=0;IxNx;Ix+) Mgrtn=(Ix*Nz+Iz)*(Nw/2+1)+Iw; /Take out Initial Wave Field Data With every Depth取出当前深度的初始波场 fseek(fp_Wfld0r,sizeof(Wfld0rIx)*Mgrtn,0); fread(&Wfld0rIx,sizeof(Wfld0rIx),1,fp_Wfld0r); fseek(fp_Wfld0
26、i,sizeof(Wfld0iIx)*Mgrtn,0); fread(&Wfld0iIx,sizeof(Wfld0iIx),1,fp_Wfld0i); /New Wave Fied Equal to the Initial on Plus Current one新波场=初始波场+从下面延拓到此处的波场 WfldrIx=WfldrIx+Wfld0rIx; WfldiIx=WfldiIx+Wfld0iIx; /Boundary absorbing of New Wave Fied边界削波：新波场=新波场削波因子 WfldrIx=WfldrIx*AbsbIx; WfldiIx=WfldiIx*AbsbIx; /4.2.2 FFT of New Wave field From Space to WaveNumber新波场FFT到波数
