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

《地震数值模拟》实验报告

一、实验题目

叠加地震记录的相移波动方程正演模拟

二、实验目的

1.掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论

2.实现方法与程序编制

3.由正演记录初步分析地震信号的分辨率。

三、实验原理

1、地震波传播的波动方程

设（x,z）为空间坐标，t为时间，地震波传播速度为v(x,z)，则二位介质中任意位置、任意时刻的地震波场为p(z,x,t)：压缩波——纵波。则二维各向同性均匀介质中地震波传播的遵循声波方程为

2、傅里叶变换的微分性质

p(t)与其傅里叶变换的P(w)的关系：

3、地震波传播的相移外推公式

令速度v不随x变化，只随z变化，则利用傅里叶变换微分性质把波动方程(变换到频率-波数域，得：

4、初始条件和边界条件

按照爆炸界面理论，反射界面震源在t=0时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为：

5、边界处理

（1）边界反射问题

把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距Xmin与Xmax两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。

（2）边界强反射的处理

镶边法、削波法、吸收边界都能有效消除边界强反射。

削波法就是在波场延拓过程中，没延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx道吸波为例，有以下吸波公式：

四、实验内容

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

#include

#include

#include

#include

//---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.参数设置定义符号--- --

#define Nt 256

#define Nz 100

#define Dx 20

#define Dt 0.004

#define Dz 20

#define pai 3.1415926

int main()

int Labs=10; //定义削波边界范围

if(Po2Judge(Nt) !=1) { printf("Nt=%d is not the Power of

2\n",Nt);exit(0); }

if(Po2Judge(Nx) !=1) { printf("Nx=%d is not the Power of

2\n",Nt);exit(0); }

if(Absorb(Labs) !=1) { printf("Absorb is error\n"); exit(0); }

if(Rflct() !=1) { printf("Rflction is error\n"); exit(0); }

if(Vlcty() !=1) { printf("Vlcty is error\n"); exit(0); }

if(WvFld0() !=1) { printf("WvFld is error\n"); exit(0); }

if(PsFrwd() !=1) { printf("PsFrwd is error\n"); 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 fr[4096],fi[4096];

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;

fr[it] =pr[is];

fi[it] = pi[is];

pr[0] = 1.0;

pi[0] = 0.0;

p = 6.283185306/(1.0*n);

pr[1] = (float) cos(p);

pi[1] = -(float)sin(p);

if (l!=0)

pi[1]=-pi[1];

for (i=2; i<=n-1; i++)

p = pr[i-1]*pr[1];

q = pi[i-1]*pi[1];

s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1]);

pr[i] = p-q;

pi[i] = s-p-q;

for (it=0; it<=n-2; it=it+2)

vr = fr[it];

vi = fi[it];

fr[it] = vr+fr[it+1];

fi[it] = vi+fi[it+1];

fr[it+1] = vr-fr[it+1];

fi[it+1] = vi-fi[it+1];

m = n/2;

nv = 2;

for (l0=k-2; l0>=0; l0--)

m = m/2;

nv = 2*nv;

for(it=0; it<=(m-1)*nv; it=it+nv)

for (j=0; j<=(nv/2)-1; j++)

p = pr[m*j]*fr[it+j+nv/2];

q = pi[m*j]*fi[it+j+nv/2];

s = pr[m*j]+pi[m*j];

s = s*(fr[it+j+nv/2]+fi[it+j+nv/2]);

poddr = p-q;

poddi = s-p-q;

fr[it+j+nv/2] = fr[it+j]-poddr;

fi[it+j+nv/2] = fi[it+j]-poddi;

fr[it+j] = fr[it+j]+poddr;

fi[it+j] = fi[it+j]+poddi;

if(l!=0)

for(i=0; i<=n-1; i++)

fr[i] = fr[i]/(1.0*n);

fi[i] = fi[i]/(1.0*n);

if(il!=0)

for(i=0; i<=n-1; i++)

pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);

if(fabs(fr[i])<0.000001*fabs(fi[i]))

if ((fi[i]*fr[i])>0)

pi[i] = 90.0;

else

pi[i] = -90.0;

else

pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;

for(i=0;i

pr[i]=fr[i];

pi[i]=fi[i];

return(1);

//***调用削波函数，形成削波数据存入一个文件**/

int Absorb(int Lx) //Nx已为全局变量不参与传递

FILE *fp_Abs;

int Ix;

float Abs[Nx];

if((fp_Abs=fopen("Absb.dat","wb"))==NULL) {printf("Connot open file""Absb"""); }

for(Ix=0;Ix

Abs[Ix]=1;//*****

for(Ix=0;Ix

Abs[Ix]=sqrt(sin((pai/2)*(Ix/(Lx-1))));

Abs[Nx-Ix-1]=Abs[Ix];//*****

for(Ix=0;Ix

fwrite(&Abs[Ix],sizeof(Abs[Ix]),1,fp_Abs);

fclose(fp_Abs);

return 1;

/******反射系数_构造形态模型的生成*****/

int Rflct()

FILE *fp_Rflct;

int Ix,Iz;

float Rflct[Nz];

if((fp_Rflct=fopen("Rflct.dat","wb"))==NULL)

{printf("Connot open file""Reflection""");}

for(Ix=0;Ix

for(Iz=0;Iz

Rflct[Iz]=0;//*****

if(Ix==Nx/2-1&&Iz==Nz/2-1)//*****

Rflct[Iz]=1;//*****

fwrite(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct);

fclose(fp_Rflct);

return 1;

int Vlcty() /********速度模型的生成*********/

FILE *fp_Vlcty;

int Ix,Iz;

float Vlcty[Nz];

if((fp_Vlcty=fopen("Vlcty.dat","wb"))==NULL)

{printf("Connot open file""Vlcty""");}

for(Ix=0;Ix

for(Iz=0;Iz

if(Iz<=3*Nz/4-1)//*****

Vlcty[Iz]=5000;//*****

else

Vlcty[Iz]=5500;//*****

fwrite(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);

fclose(fp_Vlcty);

return 1;

/********爆炸反射界面的生成，并调用FFT函数变换到频率域储存*******/

int WvFld0()

FILE *fp_Rflct,*fp_Wfld0r,*fp_Wfld0i;

int Ix,Iz,It;

float Rflct[Nz],Wfld0r[Nt],Wfld0i[Nt];

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" ,"rb"))==NULL)

{printf("Connot open Rflct.dat");}

for(Ix=0;Ix

printf("Wavefield0_FFT: Ix=%d\n",Ix);

for(Iz=0;Iz

fread(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct);

for(It=0;It

Wfld0r[It]=0;

Wfld0i[It]=0;

if(It==0) Wfld0r[It]=Rflct[Iz];//*****

if(kkfft(Wfld0r,Wfld0i,Nt,0)!=1) {printf("FFT is error");exit(0);}

for(It=0;It

fwrite(&Wfld0r[It],sizeof(Wfld0r[It]),1,fp_Wfld0r);//*****

fwrite(&Wfld0i[It],sizeof(Wfld0i[It]),1,fp_Wfld0i);//*****

fclose(fp_Rflct);

fclose(fp_Wfld0r);

fclose(fp_Wfld0i);

return 1;

/////**********PhaseShift Forward Modeling **********/////

int PsFrwd() //--波场相移延拓int PhaseShift( ); // Requset Function:PhaseShift

调用波长函数相移延拓计算函

数

int Frqcy2Time( ); //调用波场做IFFT从频率域变换到时间域函数

if ( PhaseShift( ) !=1 ) {printf("PhaseShift is error\n"); exit(0); }// Call Function

if ( Frqcy2Time( ) !=1 ) {printf("Frqcy2Time is error\n"); exit(0); }// Call Function

return 1;

int PhaseShift()

// 1. Prepprocedure预处理

FILE

*fp_Wfldr,*fp_Wfldi,*fp_Wfld0r,*fp_Wfld0i,*fp_Vlcty,*fp_Absb;

float kz[2];

int Ix,Ikx,Nkx=Nx,Iz,Iw,Nw=Nt;

long Mgrtn;

float Vlcty[Nz];

float Absb[Nx],Wfld0r[Nx],Wfld0i[Nx],Wfldr[Nx],Wfldi[Nx];

float Wfld_r,Wfld_i;

float Kxmax,Dkx,Wmax,Dw;

Wmax = 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 file Vlcty.dat\n");}

for(Iz=0;Iz

fread(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);

fclose(fp_Vlcty);

if((fp_Absb = fopen("Absb.dat","rb"))==NULL) {printf("Cann't open file Absb.dat\n");}

for(Ix=0;Ix

fread(&Absb[Ix],sizeof(Absb[Ix]),1,fp_Absb);

fclose(fp_Absb);

// 3. Open Initial Wave Field File and Current Wave Field File using In Wave Fied Extrapolating波场文件打开

if((fp_Wfld0r = fopen("Wfld0r.dat","rb"))==NULL)

{printf("Cann't open file Wfld0r.dat\n");}

if((fp_Wfld0i = fopen("Wfld0i.dat","rb"))==NULL)

{printf("Cann't open file Wfld0i.dat\n");}

if((fp_Wfldr = fopen("Wfldr.dat","wb")) ==NULL)

{printf("Cann't open file Wfldr.dat \n");}

if((fp_Wfldi = fopen("Wfldi.dat","wb")) ==NULL)

{printf("Cann't open file Wfldi.dat \n");}

// 4. 每个频率的波场延拓

for(Iw=0;Iw

// 4.1初始化当前波场

for(Ix=0;Ix

Wfldr[Ix]=0.;

Wfldi[Ix]=0.;

// 4.2波场从Iz=Nz-1最深处开始，延拓到

Iz=1测线深度

for(Iz=Nz-1;Iz>0;Iz--)

// 4.2.1形成新波场

for(Ix=0;Ix

// 1. Take out Initial Wave Field Data With every Depth取出当前深度的初始波场

Mgrtn=(Ix*Nz+1+Iz)*(Nt/2+1)+Iw;

fseek(fp_Wfld0r,sizeof(Wfld0r[Ix])*Mgrtn,SEEK_SET);

fseek(fp_Wfld0i,sizeof(Wfld0i[Ix])*Mgrtn,SEEK_SET);

fread(&Wfld0r[Ix],sizeof(Wfld0r[Ix]),1,fp_Wfld0r);

fread(&Wfld0i[Ix],sizeof(Wfld0i[Ix]),1,fp_Wfld0i);

// 2.新波场=初始波场+从下面延拓到此处的波场

Wfldr[Ix] = Wfldr[Ix]+Wfld0r[Ix];

Wfldi[Ix] = Wfldi[Ix]+Wfld0i[Ix];

// 3.边界削波：新波场=新波场×削波因子

Wfldr[Ix] = Wfldr[Ix]*Absb[Ix];

Wfldi[Ix] = Wfldi[Ix]*Absb[Ix];

// 4.2.2 新波场FFT到波数

域

if( kkfft(Wfldr,Wfldi,Nx,0) !=1 ) { printf("FFT is error\n");exit(0); }

// 4.2.3频率-波数域波场在从Iz+1延拓到Iz

for(Ikx=0;Ikx

// 1.计算相移数据expikzdz（实部、虚部

if( exp_ikzDz(kz,Ix, (float)(Vlcty[Iz]/2.),

Iw,Dw,Dkx) !=1) { printf("exp_ikzDz is error\n");exit(0); };

// 2.波场延拓：波场=波场×相移数据

Wfld_r = Wfldr[Ikx]*kz[0]-Wfldi[Ikx]*kz[1];

Wfld_i = Wfldi[Ikx]*kz[0]+Wfldr[Ikx]*kz[1];

Wfldr[Ikx] = Wfld_r;

Wfldi[Ikx] = Wfld_i;

if(Ikx!=0&&Ikx!=Nkx/2)

Wfld_r =

kz[0]*Wfldr[Nkx-Ikx]-kz[1]*Wfldi[Nkx-Ikx];

Wfld_i =

kz[1]*Wfldr[Nkx-Ikx]+kz[0]*Wfldi[Nkx-Ikx];

Wfldr[Nkx-Ikx] = Wfld_r;

Wfldi[Nkx-Ikx] = Wfld_i;

// 4.2.4 波场反FFT到空间域

if( kkfft(Wfldr,Wfldi,Nkx,1) !=1 ) { printf("FFT is error\n");exit(0); }

// 4.3 存储延拓到了测线的波场

for(Ix=0;Ix

fwrite(&Wfldr[Ix],sizeof(Wfldr[Ix]),1,fp_Wfldr);

fwrite(&Wfldi[Ix],sizeof(Wfldi[Ix]),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 Wfld0r.dat

//remove("Wfld0i.dat"); //Delete Wfld0i.dat

return(1);

int exp_ikzDz(float eikzdz[],int Ix,float Vc,int Iw,float Dw,float Dkx)

float kz=0;

eikzdz[0]=0;

eikzdz[1]=0;

kz=sqrt(pow(Iw*Dw/Vc,2)-pow(Ix*Dkx,2));

if(kz>0)

eikzdz[0]=(float)cos(kz*Dz);

eikzdz[1]=(float)-sin(kz*Dz);

return 1;

FILE *fp_Wfldr,*fp_Wfldi;

FILE *fp_Record;

int Ix,It,Iw,Nw=Nt;

float Wfldtr[Nt],Wfldti[Nt];

long AddFrmStrt;

if((fp_Wfldr = fopen("Wfldr.dat","rb")) ==NULL)

{printf("Cann't open file Wfldr.dat ");exit(0);}

if((fp_Wfldi = fopen("Wfldi.dat","rb")) ==NULL)

{printf("Cann't open file Wfldi.dat ");exit(0);}

if((fp_Record = fopen("Record.dat","wb")) ==NULL)

{printf("Cann't open file Record.dat ");exit(0);}

for(Ix=0;Ix

for(Iw=0;Iw

AddFrmStrt=Iw*Nx+Ix;

fseek(fp_Wfldr,sizeof(Wfldtr[Iw])*AddFrmStrt,SEEK_SET);

fseek(fp_Wfldi,sizeof(Wfldti[Iw])*AddFrmStrt,SEEK_SET);

fread(&Wfldtr[Iw],sizeof(Wfldtr[Iw]),1,fp_Wfldr);

fread(&Wfldti[Iw],sizeof(Wfldti[Iw]),1,fp_Wfldi);

if(Iw!=0&&Iw!=Nw/2)

Wfldtr[Nw-Iw] = Wfldtr[Iw];

Wfldti[Nw-Iw] = -Wfldti[Iw];

if(kkfft(Wfldtr,Wfldti,Nw,1)!=1){printf("FFT is error");exit(0);}

for(It=0;It

fwrite(&Wfldtr[It],sizeof(Wfldtr[It]),1,fp_Record);

} //按道存入数据，实部数据存入

fclose(fp_Wfldr);

fclose(fp_Wfldi);

remove("Wfldr.dat");

remove("Wfldi.dat");

fclose(fp_Record);

return 1;

七、实验结果

1、结果显示

①．改变绕射点位置观察截图：

②.改变速度观察截图

2、对比分析

①削波后曲线更为平滑清晰

②改变绕射点位置时，曲线左右移动，移动趋势同绕射点移动方向一致

③改变速度时，由公式可知，速度越大曲线越缓，其物理意义为，速度大走时短

八、讨论建议

1、实验收获

大量上机实践，编写褶积函数的C语言源程序有很大程度提高，在与同学交流讨论过程中收获很大，不管是学习还是工作要扬长避短，但更要面对自己的不足，努力改善不足，完善各方面能力，尤其是编程能力。通过对比分析的实践，学会了从差异中学习，在差异中进行深入研究讨论学习。

2、存在问题

编程能力有待提高，对专业理论只是有待巩固加强。

附：心得体会

其一，世上无难事只怕有心人，虽然对编程有惧，但是用心一步一步的进行其实不是很难。完成本次试验后感觉编程有很大提到，得益于老师提供提供的原程序段，让学生依葫芦画瓢，降低了学习难度。

其二，功在平时，每次做好实验进行存档整理很重要，养成习惯不管是对学业、生活或是人生都是无形的财富。

其三，学习可以和同学一起讨论但是深入研究需要自己完成，其中的收获只有自己慢慢的领悟。

毕业设计（论文）原创性声明和使用授权说明

原创性声明

本人郑重承诺：所呈交的毕业设计（论文），是我个人在指导教师的指导下进行的研究工作及取得的成果。尽我所知，除文中特别加以标注和致谢的地方外，不包含其他人或组织已经发表或公布过的研究成果，也不包含我为获得及其它教育机构的学位或学历而使用过的材料。对本研究提供过帮助和做出过贡献的个人或集体，均已在文中作了明确的说明并表示了谢意。

作者签名：日期：

指导教师签名：日期：

使用授权说明

本人完全了解大学关于收集、保存、使用毕业设计（论文）的规定，即：按照学校要求提交毕业设计（论文）的印刷本和电子版本；学校有权保存毕业设计（论文）的印刷本和电子版，并提供目录检索与阅览服务；学校可以采用影印、缩印、数字化或其它复制手段保存论文；在不以赢利为目的前提下，学校可以公布论文的部分或全部内容。

作者签名：日期：

学位论文原创性声明

本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承担。

作者签名：日期：年月日

学位论文版权使用授权书

本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授权大学可以将本学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。

涉密论文按学校规定处理。

作者签名：日期：年月日

导师签名：日期：年月日

致谢

时间飞逝，大学的学习生活很快就要过去，在这四年的学习生活中，收获了很多，而这些成绩的取得是和一直关心帮助我的人分不开的。

首先非常感谢学校开设这个课题，为本人日后从事计算机方面的工作提供了经验，奠定了基础。本次毕业设计大概持续了半年，现在终于到结尾了。本次毕业设计是对我大学四年学习下来最好的检验。经过这次毕业设计，我的能力有了很大的提高，比如操作能力、分析问题的能力、合作精神、严谨的工作作风等方方面面都有很大的进步。这期间凝聚了很多人的心血，在此我表示由衷的感谢。没有他们的帮助，我将无法顺利完成这次设计。

首先，我要特别感谢我的知道郭谦功老师对我的悉心指导，在我的论文书写及设计过程中给了我大量的帮助和指导，为我理清了设计思路和操作方法，并对我所做的课题提出了有效的改进方案。郭谦功老师渊博的知识、严谨的作风和诲人不倦的态度给我留下了深刻的印象。从他身上，我学到了许多能受益终生的东西。再次对周巍老师表示衷心的感谢。

其次，我要感谢大学四年中所有的任课老师和辅导员在学习期间对我的严格要求，感谢他们对我学习上和生活上的帮助，使我了解了许多专业知识和为人的道理，能够在今后的生活道路上有继续奋斗的力量。

另外，我还要感谢大学四年和我一起走过的同学朋友对我的关心与支持，与他们一起学习、生活，让我在大学期间生活的很充实，给我留下了很多难忘的回忆。

最后，我要感谢我的父母对我的关系和理解，如果没有他们在我的学