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

1、/预处理/头文件/#include "ps.h"/需要的头文件和初始定义在此头文件#include "fft.c"/fft变换函数存储文件/主函数/int main()float fm=75, dt=0.004;int Labs=10;/削波的宽度if(Po2Judge(Nt)!=1)printf("Nt=%d is not the power of 2n",Nt);exit(0);/判断Nt是否为2的整数次方if(Po2Judge(Nx) !=1)printf("Nx=%d is not the power of 2n&q

2、uot;,Nx);exit(0);/判断Nx是否为2的整数次方if(Absorb(Labs)!=1)printf("Absorb is errorn");exit(0);/削波程序是否有错if(Rflct(fm,dt)!=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()!=

3、1)printf("PsFwrd is errorn");exit(0);return 0;/ 功能函数实现/1.判断某数据是否为2的次方/int Po2Judge(int N)int k=0;long Ln=0;for(k=0;N-Ln>0;k+)Ln=(long)pow(2,k);/Ln为小于N的最大2的整数次方数Ln=(long)pow(2,k-1);if(fabs(Ln-N)>=1)return(0);return(1);/2.快速傅里叶变换/3.削边界/int Absorb(int Labs)FILE *fp_Abs;int Ix;float AbsN

4、x;if(fp_Abs=fopen("Absb.dat","wb")=NULL)printf("Connot open file ""Absb""");for(Ix=0;Ix<Nx;Ix+)/削波前全初始化赋值为1AbsIx=1.;for(Ix=0;Ix<Labs;Ix+)/对左侧Ix长度进行削波AbsIx=(float)(pai/2.);AbsIx=(float)(AbsIx/(Labs-1.);AbsIx=(float)(Ix*AbsIx);AbsIx=(float)(sin(A

5、bsIx);AbsIx=(float)(sqrt(AbsIx);AbsNx-1-Ix=AbsIx;/由对称性，对右侧Ix长度进行削波（与左侧对称相等） for(Ix=0;Ix<Nx;Ix+)fwrite(&AbsIx,sizeof(AbsIx),Nx,fp_Abs);fclose(fp_Abs);return (1);/4.反射系数/int Rflct(float fm,float dt)FILE *fp_Rflct;int Ix,Iz;float t,RflctNz,WaNx;if(fp_Rflct=fopen("Rflct.dat","wb&qu

6、ot;)=NULL)printf("Connot open file ""Reflection""");for(Ix=0;Ix<Nx;Ix+) WaIx=0.; Wa63=1.; for(Ix=0;Ix<Nx;Ix+) t=(float)(Ix-Nx/2)*dt;WaIx=(float)(1-2*pai*pai*fm*fm*t*t)*(float)exp(-2*pai*pai*fm*fm*t*t);/雷克子波的计算？if(Ix<=25) WaIx=0.;WaNx-Ix-2=WaIx;for(Ix=0;Ix<Nx

7、;Ix+)for(Iz=0;Iz<Nz;Iz+)RflctIz=0;/区域点赋值0if(Iz=(Nz/2)RflctIz=WaIx;/中心红点处赋值=1fwrite(&RflctIz,sizeof(RflctIz),1,fp_Rflct);fclose(fp_Rflct);return (1);/5.速度模型/int Vlcty()FILE *fp_Vlcty;int Ix,Iz;float VlctyNz;if(fp_Vlcty=fopen("Vlcty.dat","wb")=NULL)printf("Connot open f

8、ile ""Vlcty""");for(Ix=0;Ix<Nx/4;Ix+)for(Iz=0;Iz<Nz;Iz+)VlctyIz=5000;/0到Nz深度 v=5000m/sfwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for(Iz=(int)(3*Nz/4);Iz<Nz;Iz+)VlctyIz=5500;/3*Nz/4到Nz 深度v=5500m/s,此时其他还是5000m/s fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for

9、(Ix=Nx/4;Ix<Nx/2;Ix+)for(Iz=0;Iz<Nz;Iz+)VlctyIz=5000;/0到Nz深度 v=5000m/sfwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for(Iz=(int)(Nz/2);Iz<Nz;Iz+)VlctyIz=5500;/3*Nz/4到Nz 深度v=5500m/s,此时其他还是5000m/s fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for(Ix=Nx/2;Ix<3*Nx/4;Ix+)for(Iz=0;Iz<Nz

10、;Iz+)VlctyIz=5000;/0到Nz深度 v=5000m/sfwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for(Iz=(int)(4*Nz/5);Iz<Nz;Iz+)VlctyIz=5500;/3*Nz/4到Nz 深度v=5500m/s,此时其他还是5000m/s fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for(Ix=3*Nx/4;Ix<Nx;Ix+)for(Iz=0;Iz<Nz;Iz+)VlctyIz=5000;/0到Nz深度 v=5000m/sfwrite

11、(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);for(Iz=(int)(3*Nz/4);Iz<Nz;Iz+)VlctyIz=5500;/3*Nz/4到Nz 深度v=5500m/s,此时其他还是5000m/s fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);fclose(fp_Vlcty);return (1);/6.反射界面生成/int WvFld0()FILE *fp_Rflct,*fp_Wfld0r,*fp_Wfld0i;int Ix,Iz,It;float RflctNz,Wfld0rNt,Wfl

12、d0iNt;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 o

13、pen Rflct.dat");for(Ix=0;Ix<Nx;Ix+)printf("Wavefield0_FFT: Ix=%dn",Ix);for(Iz=0;Iz<Nz;Iz+)fread(&RflctIz,sizeof(RflctIz),1,fp_Rflct);for(It=0;It<Nt;It+)/实部，虚部初始化为0,即反射界面T=0时刻起爆Wfld0rIt=0.;Wfld0iIt=0.;Wfld0r0=RflctIz;if(kkfft(Wfld0r,Wfld0i,Nt,0)!=1)printf("FFT is erro

14、r"); 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_Wfld0r);fclose(fp_Wfld0i);fclose(fp_Rflct);return(1);/7.波场相移延拓/int PsFrwd()/申明功能/int PhaseShift();/相移计算int Frqcy2Time();

15、/波场做FFT,频域转化到时域/功能调用,若出错则退出/if (PhaseShift()!=1) printf("PhaseShift is errorn");exit(0);if (Frqcy2Time()!=1) printf("Frqcy2Time is errorn");exit(0);return (1);/相移延拓/int PhaseShift()/1.预处理; FILE *fp_Wfldr,*fp_Wfldi,*fp_Wfld0r,*fp_Wfld0i,*fp_Vlcty,*fp_Absb;float kz2;int Ix,Iz,Iw,Nw

16、=Nt;long Mgrtn;float VlctyNz,AbsbNx,WfldrNx,WfldiNx,Wfld0iNx,Wfld0rNx;float Kxmax,Dkx,Wmax,Dw,Newdr,Newdi;Kxmax=(float)(2*pai/Dx);Dkx=(float)(Kxmax/Nx);/采样定理：波数_空间Wmax=(float)(2*pai/Dt);Dw =(float)(Wmax/Nt);/采样定理：频率_时间/2.速度与削波数据读入if(fp_Vlcty=fopen("Vlcty.dat","rb")=NULL)printf(&q

17、uot;connot open Rflct.dat");for(Iz=0;Iz<Nz;Iz+)fread(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);fclose(fp_Vlcty);if(fp_Absb=fopen("Absb.dat","rb")=NULL) printf("connot open Absb.dat");for(Ix=0;Ix<Nx;Ix+)fread(&AbsbIx,sizeof(AbsbIx),1,fp_Absb);fclose(fp_Absb

18、);remove("Absb.dat");/3.波场文件打开 if(fp_Wfld0r=fopen("Wfld0r.dat","rb")=NULL) printf("Connot open Wfld0r.dat");if(fp_Wfld0i=fopen("Wfld0i.dat","rb")=NULL) printf("Connot open Wfld0i.dat");if(fp_Wfldr =fopen("Wfldr.dat",&quo

19、t;wb")=NULL) printf("Connot open Wfldr.dat");if(fp_Wfldi =fopen("Wfldi.dat","wb")=NULL) printf("Connot open Wfldi.dat");/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最深处开始，延

20、拓到Iz=1测线深度for(Iz=Nz-1;Iz>0;Iz-)/4.2.1 形成新波场for(Ix=0;Ix<Nx;Ix+)/取出当前深度的初始波场实部和虚部 ,在什么位置？离文件起始距离或字节数是多少？Mgrtn=(Ix*Nz+Iz)*(Nt/2+1)+Iw;fseek(fp_Wfld0r, sizeof(Wfld0rIx)*Mgrtn,0);fread(&Wfld0rIx,sizeof(Wfld0rIx),1,fp_Wfld0r);fseek(fp_Wfld0i, sizeof(Wfld0iIx)*Mgrtn,0);fread(&Wfld0iIx,sizeof(

21、Wfld0iIx),1,fp_Wfld0i);/形成新波场准备向上延拓=原始爆炸界面波场+下一层延拓到此层的波场WfldrIx=WfldrIx+Wfld0rIx;WfldiIx=WfldiIx+Wfld0iIx;/边界削波：新波场=新波场×削波因子(实部和虚部分别作同样的计算)WfldrIx=WfldrIx*AbsbIx;WfldiIx=WfldiIx*AbsbIx;/4.2.2 新波场FFT到波数域 if(kkfft(Wfldr,Wfldi,Nx,0)!=1)printf("FFT is error");exit(0);/4.2.3 频率-波数域波场在从Iz+1

22、延拓到Izfor( Ix=0;Ix<Nx/2+1;Ix+)/4.2.3.1 计算相移数据expikzdz（实部、虚部） if(exp_ikzDz(kz,Ix,(float)(VlctyIz/2.),Iw,Dw,Dkx)!=1)printf("exp_ikzDz is error");exit(0);/4.2.3.2 波场延拓：波场=波场×相移数据,下面为基本复数运算 Newdr=WfldrIx*kz0-WfldiIx*kz1;Newdi=WfldiIx*kz0+WfldrIx*kz1;WfldrIx=Newdr; WfldiIx=Newdi;if(Ix!=0

23、&&Ix!=Nx/2) Newdr=WfldrNx-Ix*kz0-WfldiNx-Ix*kz1;Newdi=WfldiNx-Ix*kz0+WfldrNx-Ix*kz1;WfldrNx-Ix=Newdr;WfldiNx-Ix=Newdi;/4.2.4 波场反FFT到空间域if(kkfft(Wfldr, Wfldi,Nx,1) !=1 ) printf("FFT is error");exit(0);/4.3 存储延拓到了测线的波场for(Ix=0;Ix<Nx;Ix+)/波场从最深处延拓到侧线深度：Z=Z0fwrite(&WfldrIx,sizeo

24、f(WfldrIx),1,fp_Wfldr);fwrite(&WfldiIx,sizeof(WfldiIx),1,fp_Wfldi);/5.关闭文件，删除中间文件。fclose(fp_Wfld0r); remove("Wfld0r.dat"); fclose(fp_Wfld0i); remove("Wfld0i.dat");fclose(fp_Wfldr); fclose(fp_Wfldi);return (1);/*相移因子函数*/ int exp_ikzDz(float eikzdz,int Ix,float Vc,int Iw,float

25、Dw,float Dkx)float kz=0;eikzdz0=0.;eikzdz1=0.;kz=(float)(pow(Iw*Dw/Vc,2)-pow(Ix*Dkx,2);if(kz>0.)kz=(float)sqrt(kz);eikzdz0=(float)cos(kz*Dz);/实部eikzdz1=-(float)sin(kz*Dz);/虚部return(1);/int Frqcy2Time() /频率域到时间域FILE *fp_Wfldr,*fp_Wfldi;FILE *fp_Record;int Ix,It,Iw,Nw=Nt;float WfldtrNt,WfldtiNt;long AddFrmStrt;if(fp_Wfldr =fopen("Wfldr.dat","rb")=NULL)printf("connot open file Wfldr.dat""");exit(0);if(fp_Wfldi =fopen("Wfldi.dat","rb")=NULL)p