叠加地震记录的相移波动方程正演模拟数值模拟实验书

1、地震数值模拟实验报告实验时间2013年6月1-30日开课单位地球物理学院指导教师XXXXXXXXXXXXXX实验题目： 叠加地震记录的相移波动方程正演模拟姓名学号班级勘查四班专业勘查技术与工程院系地球物理学院单项成绩内容理解写作结构程序设计模型设计计算结果结果分析总 成 绩一、 实验题目叠加地震记录的相移波动方程正演模拟二、 实验目的1.掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论2.实现方法与程序编制3.由正演记录初步分析地震信号的分辨率。三、 实验原理1、 地震波传播的波动方程设（x,z）为空间坐标，t为时间，地震波传播速度为v(x,z)，则二位介质中任

2、意位置、任意时刻的地震波场为p(z,x,t)：压缩波纵波。则二维各向同性均匀介质中地震波传播的遵循声波方程为2、 傅里叶变换的微分性质p(t)与其傅里叶变换的P(w)的关系：ïïîï ï íì = = ò ò¥¥ - ¥3、 地震波传播的相移外推公式 令速度v不随x变化，只随z变化，则利用傅里叶变换微分性质把波动方程(变换到频率-波数域，得：4、 初始条件和边界条件按照爆炸界面理论，反射界面震源在t=0时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或

3、子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为：5、边界处理（1）边界反射问题把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距Xmin与Xmax两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。（2）边界强反射的处理镶边法、削波法、吸收边界都能有效消除边界强反射。削波法就是在波场延拓过程中，没延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx道吸波为例，有以下吸波公式

4、：四、 实验内容1、 基本要求（1） 点绕射构造和水平层状速度模型（参数如图1所示）的正演数值模拟；五、 实验步骤1、参数初始化；2、形成边界削波数据；3、波场初始化；4、Zmax层波场延拓到深度Zmax-1；5、Zi+1层波场延拓到深度Zi；6、重复5，从Iz=Nz-1开始，直到Iz=1，得测线上的频率空间域波场；7、频率-空间域波场对频率做反傅里叶变换，得时间-空间波场；8、使用Fimage软件显示所得结果。六、实验编程文档可自由编辑打印 /*1.头文件*/#include<stdio.h>#include<math.h>#include<conio.h>

5、;#include<stdlib.h>#include<string.h>/-2.Function Request功能要求函数说明-int kkfft(float *,float *,int,int);int Absorb(int); /削波函数int Rflct(); /反射函数int Vlcty(); /速度函数int WvFld0(); /波长函数int exp_ikzDz(float *,int,float,int,float,float);/int PsFrwd();/int Po2Judge(int);/-#define Nx 128 /-3.参数设置定义符

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的倍数/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;/定义快速傅立叶函数/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); /*调用削波函数，形成削波数据存入一个文件*/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; /*反射系数_构造形态模型的生成*/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;/*爆炸反射界面的生成，并调用FFT函数变换到频率域储存*/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时刻起爆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+) /利用付氏变换的对称性，保存一半的数据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调用波长函数相移延拓计算函数int Frqcy2Time( ); /调用波场做IFFT从频率域变换到时间域函数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预处理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速度与削波数据读入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. 每个频率的波场延拓 for(Iw=0;Iw<Nw/2+1;Iw+) / 4.1初始化当前波场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取出当前深度的初始波场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到波数域if( kkfft(Wfldr,Wfldi,Nx,0) !=1 ) printf("FFT is errorn");exit(0); / 4.2.3频率-波数域波场在从Iz+1延拓到Izfor(Ikx=0;Ikx<Nx/2+1;Ikx+) / 1.计算相移数据expikzdz（实部、虚部if( exp_

28、ikzDz(kz,Ix, (float)(VlctyIz/2.), Iw,Dw,Dkx) !=1) printf("exp_ikzDz is errorn");exit(0); ; / 2.波场延拓：波场=波场×相移数据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.关闭文件，删除中间文件。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+) /按道取出数据，按Ix循环fwrite(&WfldtrIt,sizeof(WfldtrI

35、t),1,fp_Record); /按道存入数据，实部数据存入fclose(fp_Wfldr);fclose(fp_Wfldi);remove("Wfldr.dat");remove("Wfldi.dat");fclose(fp_Record);return 1;七、 实验结果1、 结果显示改变绕射点位置观察截图： .改变速度观察截图 2、 对比分析削波后曲线更为平滑清晰改变绕射点位置时，曲线左右移动，移动趋势同绕射点移动方向一致改变速度时，由公式可知，速度越大曲线越缓，其物理意义为，速度大走时短八、 讨论建议1、 实验收获大量上机实践，编写褶积函数的C

36、语言源程序有很大程度提高，在与同学交流讨论过程中收获很大，不管是学习还是工作要扬长避短，但更要面对自己的不足，努力改善不足，完善各方面能力，尤其是编程能力。通过对比分析的实践，学会了从差异中学习，在差异中进行深入研究讨论学习。2、存在问题编程能力有待提高，对专业理论只是有待巩固加强。附：心得体会 其一，世上无难事只怕有心人，虽然对编程有惧，但是用心一步一步的进行其实不是很难。完成本次试验后感觉编程有很大提到，得益于老师提供提供的原程序段，让学生依葫芦画瓢，降低了学习难度。其二，功在平时，每次做好实验进行存档整理很重要，养成习惯不管是对学业、生活或是人生都是无形的财富。其三，学习可以和同学一起讨

37、论但是深入研究需要自己完成，其中的收获只有自己慢慢的领悟。毕业设计（论文）原创性声明和使用授权说明原创性声明本人郑重承诺：所呈交的毕业设计（论文），是我个人在指导教师的指导下进行的研究工作及取得的成果。尽我所知，除文中特别加以标注和致谢的地方外，不包含其他人或组织已经发表或公布过的研究成果，也不包含我为获得 及其它教育机构的学位或学历而使用过的材料。对本研究提供过帮助和做出过贡献的个人或集体，均已在文中作了明确的说明并表示了谢意。作 者 签 名： 日 期： 指导教师签名： 日期： 使用授权说明本人完全了解 大学关于收集、保存、使用毕业设计（论文）的规定，即：按照学校要求提交毕业设计（论文）的印

38、刷本和电子版本；学校有权保存毕业设计（论文）的印刷本和电子版，并提供目录检索与阅览服务；学校可以采用影印、缩印、数字化或其它复制手段保存论文；在不以赢利为目的前提下，学校可以公布论文的部分或全部内容。作者签名： 日 期： 学位论文原创性声明本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承担。作者签名： 日期： 年 月 日学位论文版权使用授权书本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并向