声波方程数值模拟实验报告

声波方程数值模拟实验报告

一．基础理论知识

需要的已知条件包括：

1．1）震源函数

2）地层速度（波速） 3）边界条件

2．弹性波方程：?????????+??=??+??+??=??)

()()(222222

22222

222z w x w v t w t S z u x u

v t u s p 声波方程的有限差分法数值模拟

对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：

)()(2222

222t S z

u 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 展式展开：

z

?,i j

1,i j

+2,i j

+1,i j

-2,i j

-,2

i j -,2

i j +,1i j +,1

i j -1,1i j -+1,2

i j -+2,1i j -+2,2

i j -+1,2

i j ++2,2

i j ++1,1

i j +-2,1i j +-1,1i j ++2,1i j ++1,1i j --1,2i j +-2,2i j +-2,

2i j --2,1

i j --1,2i j --x

?

)(*21*22*2

2*,1,t o t t u

t t

u u

u

t

k t t k t k j

i k j

i ?+???+

???+=?=?=+

(4-2)

将1

,-k j i u 在(i,j)点k 时刻用Taylor 展式展开：

)(*21*22*2

2*,1,t o t t u t t

u u

u

t

k t t

k t k j

i k j

i ?+???+

???-=?=?=-

(4-3)

将上两式相加，略去高阶小量，整理得(i,j)点k 时刻的二阶时间微商为：

2

1

,,1,222t

u u u t u k j

i 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 u

1,1++k j i u 分别在(i,j)点k 时刻展开到四阶小量，消除四阶小量并解出二阶微分得：

}25][34][121{1,,1,1,2,22

22k j i k

j i k j i k j i k j i u u u u u x

x u -+++-?=??+-+- (4-5) 同理可得：

}25][34][12

1{1,1

,1,2,2,222k j i k

j 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,22

221,,1,k j i k j

i k j i k j i k j i k j

i k j

i k j

i u u u u u h t v u

u

u

-+++-??+-=+-+--+

k j i k

j i k j i k j i k j i u u u u u h

t v ,,1,1,2,22

2225][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()(2

2-= （4-8）

上式中，t 为时间， f 为中心频率，一般取为20-40HZ ，γ为控制频带宽度的参数，

一般取3-5。在实际计算过程中，需把此震源函数离散，参与波场计算。)

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

稳定性条件：8

3

/*max <

??h t v (4-8)

这里max v 表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-8)式时，第推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。 频散关系式：

)/(min N Gf v h ≤? (4-9)

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

二．实验步骤

1、 应用声波方程作为正演模拟的波动方程，忽略转换波的产生、传播；

2、 将所提供震源函数离散后绘图；

震源函数为雷克子波，离散绘图如下： fm=30; r=3; t=0.002; for n=1:200

w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n); end

figure(1),plot(w);

3、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，

分别记录两个时刻的波前快照（即区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻。

由稳定条件，设v为2000，可以令DT=0.001,DH=5,此时精度较高，且满足频散关系

程序，图形如下：

#include

#include

#include

#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.001

void 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

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

for(j=0;j

{if(i==50&&j==50)

f[i][j]=1;

else f[i][j]=0;

}

for(i=0;i

for(j=0;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

{for(i=2;i

for(j=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

for(n=0;n

{u1[m][n]=u2[m][n];

u2[m][n]=u3[m][n];

}

if(k==100)

for(m=0;m

for(n=0;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

for(j=0;j

fprintf(fp,"%f\n",u4[i][j]);

fclose(fp);

}

}

第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；

程序图形如下：

#include

#include

#include

#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.001

void 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

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

for(j=0;j

{if(i==50&&j==50)

f[i][j]=1;

else f[i][j]=0;

}

for(i=0;i

for(j=0;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

{for(i=2;i

for(j=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

for(n=0;n

{u1[m][n]=u2[m][n];

u2[m][n]=u3[m][n];

}

if(k==160) //只需改动K值，即显示边界反射

for(m=0;m

for(n=0;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

for(j=0;j

fprintf(fp,"%f\n",u4[i][j]);

fclose(fp);

}

}

4、对于大模型，定义为水平层状速度模型；做两个实验，一是将震源点放在区域表层任一

点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律，指出哪是反射波，哪是透射波；

这时取小模型的常量，为减少频散，速度v至少为2400

程序图形如下：

#include

#include

#include

#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.001

void 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

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

for(j=0;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; //第一层速度为2400

else

v[i][j]=3000; //第二层速度为3000

}

f[100][10]=1; //定义f函数，在100，10为1

for(k=0;k

{for(i=2;i

for(j=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

for(n=0;n

{u1[m][n]=u2[m][n];

u2[m][n]=u3[m][n];

}

if(k==100)

for(m=0;m

for(n=0;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

for(j=0;j

fprintf(fp,"%f\n",u4[i][j]);

fclose(fp);

}

}

反射波

透射波

二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波？

这时取小模型的常量，为减少频散，速度v至少为2400

程序图像如下：

#include

#include

#include

#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.001

void 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

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

for(j=0;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; //第一层速度为2400

else

v[i][j]=3100; //第二层速度为3100

}

f[100][10]=1; //定义f函数，在100，10为1

for(k=0;k

{for(i=2;i

for(j=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

for(n=0;n

{u1[m][n]=u2[m][n];

u2[m][n]=u3[m][n];

}

if(k==100)

for(m=0;m

for(n=0;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

fprintf(fp,"%f\n",u5[i][j]); fclose(fp); } }

从上向下看的两条直线为与震源同深度点的各点所有时刻的波场值，在向下与两直线紧邻的双曲线为第二层界面反射的记录。

三． 通过实验所发现的问题和认识

(1) 在用C 语言写程序时，很多小问题不注意比如将i=i+2写成i+2,在m ，n 循环中出现i ，

j 等等，导致无限循环，数据过大，不能出现dat 文件等问题。经过这次试验对C 语言和matlab 的使用有了更加深刻的认识。

（2）程序成图时，图形不合适，调试数据时出现The data file length error! 按照

8

3/*m a x <

??h t v ，)/(min N Gf v h ≤? 对数据进行多次修改，得到较合理的图形。 （3）通过这次试验对地下介质中地震波的传播规律和声波方程有了深刻的认识。