声波方程有限差分正演

题目：使用Ricker 子波，刚性边界条件，并且初值为零，在均匀各向同性介质条件下，利用交错网格法求解一阶二维声波方程数值解。

解：一阶二维声波方程：22222221zPx P t P c ∂∂+∂∂=∂∂ (1)将其分解为：21P c t Px P z x z x z V V x z V tV t ∂∂∂⎧=+⎪∂∂∂⎪∂∂⎪=⎨∂∂⎪∂∂⎪=⎪∂∂⎩(2)对分解后的声波方程进行离散，可得到：112211,-1,,,122[]N n n n n m i m j i m j xi j xi j m t VVc P P h +-+---=∆=+-∑ 112211,1,,,122[]Nn n n n m i j m i j m zi j zi j m t VV c P P h +-++---=∆=+-∑ 1111212222,,m 1,,,,11[]Nn n n n n n i ji jmxi j xi m j zi j m zi j m m tc PP cVVVVh+++++++-+--=∆=+-+-∑h z x =∆=∆针对公式(1)，使用二阶中心差商公式：2P(,,1)2(,,)(,,1)i j n P i j n P i j n t +-+-∆222(1,,)2(,,)(1,,)(,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆=⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(3)变形：P(,,1)=2(,,)(,,1)i j n P i j n P i j n +--2222(1,,)2(,,)(1,,)t (,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆+∆⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(4)对离散格式作时间和空间三重Fourier 变换：0P(,,)(,,)x z i j n P k k w ↔ ，0P(,,1)(,,)*exp()x z i j n P k k w iw t +↔∆0P(1,,)(,,)*exp(k )x z x i j n P k k w i x +↔-∆，0z P(,1,)(,,)*exp(k )x z i j n P k k w i z +↔-∆对公式(4)进行Fourier 变换：2222exp()2exp()h exp()2()exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆=--∆+∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦2222exp()2exp()h exp()2()=exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆-+-∆∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦222222sin sin 22sin (2x z k x k zw tt c h∆∆+∆=∆） (5) 公式(5)右端必须满足下列条件：22222sin sin 220(x z k x k zt c h∆∆+≤∆≤）1 取x k 和z k 最大值，即=x x k π∆，z =k z π∆，则有：22220t c h≤∆≤1因此tc ∆≤即为所求得的稳定性条件。

在程序试算中，选中相关参数如下：dt 0.001()10()3000/20010030N m s x z m v m s x z msx sz m f Hz=⎧⎪∆=∆=⎪⎪=⎪==⎪⎨⎪⎪==⎪=⎪⎪⎩时间网格空间网格（速度）（模型范围）t=1s(采样长度)（震源位置）（主频）=15(空间精度)因为v 3000*0.0013t ∆==<=，满足稳定性条件。

波长v 3000=100m 30f λ==，空间采样间隔h=10m ，一个波长内的空间采样点数 100m=1010hλ==，基本满足网格色散条件。

程序运行输出P(,,)x z t 波场快照：(a) 200ms (b) 300ms(c) 400ms (d) 500ms(e) 600ms (f) 700ms(g) 800ms (h) 900ms图1 波场快照主要程序代码附录：震源Ricker子波函数float f(float t1, float f0){float t00=1/f0,y;y=((1.0-2.0*pow(pi*f0*(t1-t00),2))*exp(-pow(pi*f0*(t1-t00),2)));return y;}计算交错网格有限差分系数void Cal_1D_FdCoff(float *c1,int n){float **A = new float * [n];for (int i = 0; i < n; i++){A[i] = new float [n];}float *x = new float [n];float *b = new float [n];//for (int i = 0; i < n; i++){for (int j = 0; j < n; j++){A[i][j] = pow( 2.0*(j+1)-1, 2.0*(i+1)-1 );}}for (int i = 0; i < n; i++){if (i < 1)b[i] = 1.0;elseb[i]=0;}//Gauss(A,n,b,x);for (int i = 1; i <= n; i++){c1[i] = x[i-1];}for (int i = 0; i < n; i++){delete[] A[i];}delete A;delete x;delete b;}主函数：void main(){FILE *fp;char name[1000];float dt, h;float Ts, f0;int NX, NZ, NT, s_x, s_z, nPoints;fp=fopen("2D_Parameters.txt", "r");fscanf(fp, "%[^\n]\n", name);fscanf(fp, "%f\n", &dt); //Time Intervalfscanf(fp, "%[^\n]\n", name);fscanf(fp, "%d\n", &NX); //X Grid Dimensionfscanf(fp, "%[^\n]\n", name);fscanf(fp, "%d\n", &NZ); //Z Grid Dimensionfscanf(fp, "%[^\n]\n", name);fscanf(fp,"%d\n",&s_x); //Source Xfscanf(fp, "%[^\n]\n", name);fscanf(fp, "%d\n", &s_z); //Source Zfscanf(fp, "%[^\n]\n", name);fscanf(fp,"%f\n",&h); //Space Intervalfscanf(fp, "%[^\n]\n", name);fscanf(fp, "%f\n", &Ts); //Total Timefscanf(fp, "%[^\n]\n", name);fscanf(fp, "%f\n", &f0); //Dominant Frequencyfscanf(fp, "%[^\n]\n", name);fscanf(fp, "%d\n", &nPoints); //Accuracy of FD methodint model;fscanf(fp, "%[^\n]\n", name);fscanf(fp, "%d\n", &model); //modelfclose(fp);NT = (int)((Ts + 1e-6)/dt);///////////////////////////////////////////////////////////////// //The process of opening the arrayfloat **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 = Create2DActivityArray(NZ,NX);Vp = Create2DActivityArray(NZ,NX);// Set ModelVelocity_Model(model, Vp, NX, NZ);float *c1;c1 = (float*)malloc( sizeof(float) * (nPoints+1) );Cal_1D_FdCoff(c1, nPoints); // Compute FD Coefficient// Print Coefficientfor (int i = 1; i <= nPoints; i++){printf("%e\n", c1[i]);}///////////////////////////////////////////////////////////////// printf("***********Forward Calculation starts!************\n");float sum1, sum3;for (int k = 0; k < NT; k++){for (int j = 0; j < NZ; j++){for (int i = 0; i < NX; i++){sum1 = 0.0;sum3 = 0.0;for (int m = 0; m < nPoints; m++){// U/Xif ( i-m-1 < 0 ){sum1 += (float)c1[m+1] * Txx1[j][i+m];}else if ( i+m > NX-1 ){sum1 -= (float)c1[m+1] * Txx1[j][i-m-1];}else{sum1 += (float)c1[m+1] * (Txx1[j][i+m] -Txx1[j][i-m-1]);}// U/Zif ( j-m < 0 ){sum3 += (float)c1[m+1] * Tzz1[j+m+1][i];}else if ( j+m+1 > NZ-1 ){sum3 -= (float)c1[m+1] * Tzz1[j-m][i];}else{sum3 += (float)c1[m+1] * (Tzz1[j+m+1][i] -Tzz1[j-m][i]);}}//U2[j][i] = -Vp[j][i]*Vp[j][i] * (dt/h) * (sum1 + sum3) + U1[j][i];}//for i}//for j///////////////////////////////////////////////////////////////// for (int j = 0; j < NZ; j++){for (int i = 0; i < NX; i++){sum1 = 0.0;sum3 = 0.0;for (int m = 0; m < nPoints; m++){// V/Xif ( i-m < 0 ){sum1 += (float)c1[m+1] * U2[j][i+m+1];}else if ( i+m+1 > NX-1 ){sum1 -= (float)c1[m+1] * U2[j][i-m];}else{sum1 += (float)c1[m+1] * (U2[j][i+m+1] -U2[j][i-m]);}// V/Zif ( j-m-1 < 0 ){sum3 += (float)c1[m+1] * U2[j+m][i];}else if ( j+m > NZ-1 ){sum3 -= (float)c1[m+1] * U2[j-m-1][i];}else{sum3 += (float)c1[m+1] * (U2[j+m][i] - U2[j-m-1][i]);}}//Txx2[j][i] = -(dt/h) * sum1 + Txx1[j][i];Tzz2[j][i] = -(dt/h) * sum3 + Tzz1[j][i];//添加震源项if ( i==s_x && j==s_z ){Txx2[j][i]+=f(k*dt, f0);}}//for i}//for j////////////////////////////////////////////////////////////////Tranfer the value.////for (int j = 0; j < NZ; j++) {for (int i = 0; i < NX; i++){U1[j][i] = U2[j][i];Txx1[j][i] = Txx2[j][i];Tzz1[j][i] = Tzz2[j][i];}}////////////////////////////////////////////////////////////if (k > 0){Fwrite_Snapshot_U(U2, NX, NZ, k, 50);}printf("The %d ms Forward file has finished!\n",k);//Finish the circulation of K.}Release2DActivityArray(U1,NZ);Release2DActivityArray(U2,NZ);Release2DActivityArray(Txx1,NZ);Release2DActivityArray(Txx2,NZ);Release2DActivityArray(Tzz1,NZ);Release2DActivityArray(Tzz2,NZ);Release2DActivityArray(Vp,NZ);}。