《地震数值模拟》实验报告一、实验题目叠加地震记录的相移波动方程正演模拟二、实验目的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<stdio.h>#include<math.h>#include<conio.h>#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.参数设置定义符号--- --#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 of2\n",Nt);exit(0); }if(Po2Judge(Nx) !=1) { printf("Nx=%d is not the Power of2\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;elsepi[i] = -90.0;elsepi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;for(i=0;i<n;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<Nx;Ix++)Abs[Ix]=1;//*****for(Ix=0;Ix<Lx-1;Ix++)Abs[Ix]=sqrt(sin((pai/2)*(Ix/(Lx-1))));Abs[Nx-Ix-1]=Abs[Ix];//*****for(Ix=0;Ix<Nx;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<Nx;Ix++)for(Iz=0;Iz<Nz;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<Nx;Ix++)for(Iz=0;Iz<Nz;Iz++)if(Iz<=3*Nz/4-1)//*****Vlcty[Iz]=5000;//*****elseVlcty[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<Nx;Ix++)printf("Wavefield0_FFT: Ix=%d\n",Ix);for(Iz=0;Iz<Nz;Iz++) //赋值，爆炸反射界面t=0时刻起爆fread(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct);for(It=0;It<Nt;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<Nt/2+1;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 Functionif ( Frqcy2Time( ) !=1 ) {printf("Frqcy2Time is error\n"); exit(0); }// Call Functionreturn 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<Nz;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<Nx;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<Nw/2+1;Iw++)// 4.1初始化当前波场for(Ix=0;Ix<Nx;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<Nx;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延拓到Izfor(Ikx=0;Ikx<Nx/2+1;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<Nx;Ix++)fwrite(&Wfldr[Ix],sizeof(Wfldr[Ix]),1,fp_Wfldr);fwrite(&Wfldi[Ix],sizeof(Wfldi[Ix]),1,fp_Wfldi);// 5.关闭文件，删除中间文件。