声波方程数值模拟实验报告一．基础理论知识需要的已知条件包括：1．1）震源函数2）地层速度（波速） 3）边界条件2．弹性波方程：⎪⎪⎩⎪⎪⎨⎧∂∂+∂∂=∂∂+∂∂+∂∂=∂∂)()()(22222222222222z w x w v t w t S zu x uv t u s p声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：)()(2222222t S zu x u v t u +∂∂+∂∂=∂∂ （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震源函数。

为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。

为此，先把空间模型网格化（如图4-1所示）。

设x 、z 方向的网格间隔长度为h ∆，t ∆为时间采样步长，则有：h i x ∆= （i 为正整数）h j z ∆= (j 为正整数)t n t =∆ (n 为正整数)k j i u , 表示在(i,j)点，k 时刻的波场值。

将1,+k j i u 在(i,j)点k 时刻用Taylor 展式展开：)(*21*22*22*,1,t o t t ut tu uutk t t k t k ji k ji ∆+∆∂∂+∆∂∂+=∆=∆=+(4-2)将1,-k j i u 在(i,j)点k 时刻用Taylor 展式展开：)(*21*22*22*,1,t o t t u t tu uutk t tk t k ji k ji ∆+∆∂∂+∆∂∂-=∆=∆=-(4-3)将上两式相加，略去高阶小量，整理得(i,j)点k 时刻的二阶时间微商为：21,,1,222t u u u t u k ji k j i k j i ∆+-=∂∂-+ (4-4)对于空间微分，采用四阶精度差分格式，（以X 方向为例）即将1,2++k j i u 、1,1++k j i u 、1,1+-k j i u1,1++k j i u 分别在(i,j)点k 时刻展开到四阶小量，消除四阶小量并解出二阶微分得：}25][34][121{1,,1,1,2,2222k j i kj i k j i k j i k j i u u u u u xx u -+++-∆=∂∂+-+- (4-5) 同理可得：}25][34][121{1,1,1,2,2,222k j i kj i k j i k j i k j i u u u u u z z u -+++-∆=∂∂+-+- (4-6)这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得：}25][34][121{2,,1,1,2,22221,,1,k j i k ji k j i k j i k j i k ji k ji k ji u u u u u h t v uuu-+++-∆∆+-=+-+--+k j i kj i k j i k j i k j i u u u u u h t v ,,1,1,2,222225][34][121{-+++-∆∆++-+-})(**)(*)(00j j i i t s --+δδ(4-7)式中),(j i v 为介质速度的空间离散值，h ∆是空间离散步长，t ∆为时间离散步长，)(k s 为震源函数，关于)(k s 一般使用一个理论的雷克型子波代替，即：ft t f t s e πγπ2cos )/2()(22-= （4-8）上式中，t 为时间， f 为中心频率，一般取为20-40HZ ，γ为控制频带宽度的参数，一般取3-5。

在实际计算过程中，需把此震源函数离散，参与波场计算。

)(*)(00j j i i --δδ确定震源位置。

稳定性条件：83/*max <∆∆h t v (4-8)这里m ax v 表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-8)式时，第推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。

频散关系式：)/(min N Gf v h ≤∆ (4-9)式中m in v 为最小速度，N f 为Nyquist 频率。

一般取震源子波中的主频f 的2倍值参与计算，G 为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分G 取8，而对于空间为四阶差分的情况则G 取4方能有效减少频散。

二．实验步骤1、 应用声波方程作为正演模拟的波动方程，忽略转换波的产生、传播；2、 将所提供震源函数离散后绘图；震源函数为雷克子波，离散绘图如下： fm=30; r=3; t=0.002; for n=1:200w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n); endfigure(1),plot(w);3、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即区域内所有网格点的波场值）。

第一时刻为地震波还未传播到边界上的某时刻。

由稳定条件，设v为2000，可以令DT=0.001,DH=5,此时精度较高，且满足频散关系程序，图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN]; //不能直接初值为0 float u5[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++) //定义f函数，当且仅当i,j同时为50时，f为1，其余为0for(j=0;j<ZN;j++){if(i==50&&j==50)f[i][j]=1;else f[i][j]=0;}for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;v[i][j]=2000; //速度相同表示同一介质}for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==100)for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照,中间点}if((fp=fopen("wavefront.dat","w"))!=NULL){fprintf(fp,"%d\n",XN);fprintf(fp,"%d\n",ZN);for(i=0;i<XN;i++)for(j=0;j<ZN;j++)fprintf(fp,"%f\n",u4[i][j]);fclose(fp);}}第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；程序图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN]; //不能直接初值为0 float u5[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++) //定义f函数，当且仅当i,j同时为50时，f为1，其余为0for(j=0;j<ZN;j++){if(i==50&&j==50)f[i][j]=1;else f[i][j]=0;}for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;v[i][j]=2000; //速度相同表示同一介质}for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==160) //只需改动K值，即显示边界反射for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照，边界}if((fp=fopen("wavefront.dat","w"))!=NULL){fprintf(fp,"%d\n",XN);fprintf(fp,"%d\n",ZN);for(i=0;i<XN;i++)for(j=0;j<ZN;j++)fprintf(fp,"%f\n",u4[i][j]);fclose(fp);}}4、对于大模型，定义为水平层状速度模型；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律，指出哪是反射波，哪是透射波；这时取小模型的常量，为减少频散，速度v至少为2400程序图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 200#define ZN 100#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],u5[XN][ZN]; //不能直接初值为0float f[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;f[i][j]=0.0;if (j<=30)v[i][j]=2400; //第一层速度为2400elsev[i][j]=3000; //第二层速度为3000}f[100][10]=1; //定义f函数，在100，10为1for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];u5[i][k]=u1[i][10]; //地震记录}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==100)for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照}if((fp=fopen("wavefront.dat","w"))!=NULL){fprintf(fp,"%d\n",XN);fprintf(fp,"%d\n",ZN);for(i=0;i<XN;i++)for(j=0;j<ZN;j++)fprintf(fp,"%f\n",u4[i][j]);fclose(fp);}}反射波透射波二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波？这时取小模型的常量，为减少频散，速度v至少为2400程序图像如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 200#define ZN 200#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],u5[XN][ZN]; //不能直接初值为0 float f[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;f[i][j]=0.0;if (j<=30)v[i][j]=2400; //第一层速度为2400elsev[i][j]=3100; //第二层速度为3100}f[100][10]=1; //定义f函数，在100，10为1for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];u5[i][k]=u1[i][10]; //地震记录}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==100)for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照}if((fp=fopen("sei_record.dat","w"))!=NULL){fprintf(fp,"%d\n",XN/2);fprintf(fp,"%d\n",KN);for(i=0;i<XN;i=i+2)for(j=0;j<KN;j++)fprintf(fp,"%f\n",u5[i][j]);fclose(fp);} }从上向下看的两条直线为与震源同深度点的各点所有时刻的波场值，在向下与两直线紧邻的双曲线为第二层界面反射的记录。