声波方程有限差分正演

1、(3)题目：使用Ricker子波，刚性边界条件，并且初值为零，在均匀各向同性介质条件下，利用交错网格法求解一阶二维声波方程数值解。解：一阶二维声波方程(1)1d2Pd2Pd2P=+c2dt20X2dz2将其分解为：10P0V0Vc20V0Pdtdx(2)0V0PzTOC o 1-5 h z0t0z对分解后的声波方程进行离散，可得到Vn+2=Vn-2+2艺cPn-Pn.丄丄.hmi+m,ji-m-1,jXi,/xi-,/22m=1+1丄AtVVn+2=Vn-2+乙cPn-Pnz-丄zi,j-1h丄mi，j+m+1i，j-m22m=1Atc2丄丄丄丄Pn+1=Pn+cVn+2-Vn+2+V+2-V

2、+2i,ji,jhmxi+m+1,jxi-m,jzi,j+mzi,j-m-1m=1Ax=Az=h针对公式(1)，使用二阶中心差商公式：P(i,j,n+1)2P(i,j,n)+P(i,j,n1)At2”P(i+1,j,n)2P(i,j,n)+P(i-1,j,n)十、Ax2=c2P(i,j+1n)2P(i,j,n)+P(i,j丄n)、Az2一变形：P(i,j,n+1)=2P(i,j,n)-P(i,j,n-1)P(i+1,j,n)-2P(i,j,n)+P(i-1,j,n)十Mt、Az2对离散格式作时间和空间三重Fourier变换P(i,j,n)吕P(k,k,w：0 xz)，P(i,j,n+1)吕P(

3、k,k,w)*exp(iwAt)0 xzP(i+1,j,n)吕P(k,k,w)*exp(-ikAx),P(i,j+1,n)吕P(k,k,w)*exp(-ikAz)0 xzx0 xzz对公式(4)进行Fourier变换：exp(iwAt)=2一(-iwAt)+At2c2exp(-ikAx)-2+exp(ikAx)+h2exp(-ikAz)-2+exp(ikAz)zzh2exp(iwAt)-2+(-iwAt)=At2c2exp(-ikAx)-2+exp(ikAx)XX+h2exp(-ikAz)-2+exp(ikAz)h2kAxkAzsm2+sm2v(5)TOC o 1-5 h zwAt.22、sm

4、2=At2c2(2h2公式(5)右端必须满足下列条件：kAxkAzsm2+sm20At2c2(2)1h2取k和k最大值，即k=,k=，则有：0At2c221xzxAxzAzh2因此AtcV即为所求得的稳定性条件。在程序试算中，选中相关参数如下：dt=0.001s（时间网格）Ax=Az=10m（空间网格）v=3000m/s（速度）x=z=200m（模型范围）t=ls（采样长度）sx=sz=100m（震源位置）f=30Hz（主频）mN=15（空间精度）满足稳定性条件。因为vAt=3000*0.001=3-L=巴，22波长长=f=警=100m空间采样间隔h=10m，一个波长内的空间采样点数m=-=1

5、00=10，基本满足网格色散条件。h10程序运行输出P（x,z,t）波场快照：(a)200ms(b)300ms(c)400ms(d)500ms(f)700ms(g)800ms(h)900ms图1波场快照主要程序代码附录：震源Ricker子波函数floatf(floatt1,floatf0)floatt00=1/f0,y;y=(1.0-2.0*pow(pi*f0*(t1-t00),2)*exp(-pow(pi*f0*(t1-t00),2);returny;计算交错网格有限差分系数voidCal_1D_FdCoff(float*c1,intn)float*A=newfloat*n;for(inti

6、=0;in;i+)Ai=newfloatn;float*x=newfloatn;float*b=newfloatn;/for(inti=0;in;i+)for(intj=0;jn;j+)Aij=pow(2.0*(j+1)-1,2.0*(i+1)-1);for(inti=0;in;i+)if(i1)bi=1.0;elsebi=0;/Gauss(A,n,b,x);for(inti=1;i=n;i+)c1i=xi-1;for(inti=0;in;i+)deleteAi;deleteA;deletex;deleteb;主函数：voidmain()FILE*fp;charname1000;floatdt

7、,h;floatTs,f0;intNX,NZ,NT,s_x,s_z,nPoints;fp=fopen(2D_Parameters.txt,r);fscanf(fp,%nn,name);fscanf(fp,%fn,&dt);/TimeIntervalfscanf(fp,%nn,name);fscanf(fp,%dn,&NX);/XGridDimensionfscanf(fp,%nn,name);fscanf(fp,%dn,&NZ);/ZGridDimensionfscanf(fp,%nn,name);fscanf(fp,%dn,&s_x);/SourceXfscanf(fp,%nn,name);

8、fscanf(fp,%dn,&s_z);/SourceZfscanf(fp,%nn,name);fscanf(fp,%fn,&h);/SpaceIntervalfscanf(fp,%nn,name);fscanf(fp,%fn,&Ts);/TotalTimefscanf(fp,%nn,name);fscanf(fp,%fn,&f0);/DominantFrequencyfscanf(fp,%nn,name);fscanf(fp,%dn,&nPoints);/AccuracyofFDmethodintmodel;fscanf(fp,%nn,name);fscanf(fp,%dn,&model);

9、/modelfclose(fp);NT=(int)(Ts+1e-6)/dt);/Theprocessofopeningthearrayfloat*U1,*U2;float*Txx1,*Txx2;float*Tzz1,*Tzz2;float*Vp;U1=Create2DActivityArray(NZ,NX);U2=Create2DActivityArray(NZ,NX);Txx1=Create2DActivityArray(NZ,NX);Txx2=Create2DActivityArray(NZ,NX);Tzz1=Create2DActivityArray(NZ,NX);Tzz2=Create

10、2DActivityArray(NZ,NX);Vp=Create2DActivityArray(NZ,NX);/SetModelVelocity_Model(model,Vp,NX,NZ);float*c1;c1=(float*)malloc(sizeof(float)*(nPoints+1);Cal_1D_FdCoff(c1,nPoints);/ComputeFDCoefficient/PrintCoefficientfor(inti=1;i=nPoints;i+)printf(%en,c1i);/printf(*ForwardCalculationstarts!*n);floatsum1,

11、sum3;for(intk=0;kNT;k+)for(intj=0;jNZ;j+)for(inti=0;iNX;i+)sum1=0.0;sum3=0.0;for(intm=0;mnPoints;m+)/U/Xif(i-m-1NX-1)sum1-=(float)c1m+1*Txx1ji-m-1;elsesum1+=(float)c1m+1*(Txx1ji+m-Txx1ji-m-1);/U/Zif(j-mNZ-1)sum3-=(float)c1m+1*Tzz1j-mi;elsesum3+=(float)c1m+1*(Tzz1j+m+1i-Tzz1j-mi);/U2ji=-Vpji*Vpji*(dt

12、/h)*(sum1+sum3)+U1ji;/fori/forj/for(intj=0;jNZ;j+)for(inti=0;iNX;i+)sum1=0.0;sum3=0.0;for(intm=0;mnPoints;m+)/V/Xif(i-mNX-1)sum1-=(float)c1m+1*U2ji-m;elsesum1+=(float)c1m+1*(U2ji+m+1-U2ji-m);/V/Zif(j-m-1NZ-1)sum3-=(float)c1m+1*U2j-m-1i;elsesum3+=(float)c1m+1*(U2j+mi-U2j-m-1i);/Txx2ji=-(dt/h)*sum1+Txx1ji;Tzz2ji=-(dt/h)*sum3+Tzz1ji;/添加震源项if(i=s_x&j=s_z)Txx2ji+=f(k*dt,f0);/fori/forj/Tranferthevalue./for(intj=0;jNZ;j+)for(inti=0;i0)Fwrite_Snapshot_U(U2,NX,NZ,k,50);printf(The%dmsForwardfilehasfinished!n,k);/FinishthecirculationofK.Release2DActivit