声波方程数值模拟实验报告.docx
《声波方程数值模拟实验报告.docx》由会员分享，可在线阅读，更多相关《声波方程数值模拟实验报告.docx（34页珍藏版）》请在冰豆网上搜索。
声波方程数值模拟实验报告
声波方程数值模拟实验报告
学号：
201107010230
姓名：
包灵
指导老师：
程冰洁老师
完成时间：
2014年11月26日星期三
实验要求：
1、应用声波方程作为正演模拟的波动方程；
2、将所提供震源函数离散后绘图；
3、给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图形来；
4、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。
第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；
5、对于大模型，定义为水平层状速度模型（至少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波。
实验目的：
1.通过本次作业，加深对波动方程的理解，明白波动方程所代表的物理意义。
2.通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规律。
3.通过模拟水平层状速度模型，体会地震波在两种介质分界面的传播规律，并能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。
4.通过模拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波形变化趋势。
需要的已知条件包括：
1）震源函数
2）地层速度（波速）
3）边界条件
2．弹性波方程：
声波方程的有限差分法数值模拟
对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：
（4-1）
是介质在点（x,z）处的纵波速度，
为描述速度位或者压力的波场，
为震源函数。
为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。
为此，先把空间模型网格化（如图4-1所示）。
设x、z方向的网格间隔长度为
，
为时间采样步长，则有：
（i为正整数）
（j为正整数）
（n为正整数）
表示在（i,j）点，k时刻的波场值。
将
在（i,j）点k时刻用Taylor展式展开：
（4-2）
将
在（i,j）点k时刻用Taylor展式展开：
（4-3）
将上两式相加，略去高阶小量，整理得（i,j）点k时刻的二阶时间微商为：
（4-4）
对于空间微分，采用四阶精度差分格式，（以X方向为例）即将
、
、
分别在（i,j）点k时刻展开到四阶小量，消除四阶小量并解出二阶微分得：
（4-5）
同理可得：
（4-6）
这就实现了用网格点波场值的差商代替了偏微分方程的微商，将上三个式子代入（4-1）式中得：
}
（4-7）
式中
为介质速度的空间离散值，
是空间离散步长，
为时间离散步长，
为震源函数，关于
一般使用一个理论的雷克型子波代替，即：
（1）
（4-8）
（2）
（4-9）
上式中，
为时间，
为中心频率，一般取为20-40HZ，
为控制频带宽度的参数，一般取3-5。
在实际计算过程中，需把此震源函数离散，参与波场计算。
确定震源位置。
稳定性条件：
（4-10）
这里
表示的是地下介质的最大波速；若地下介质网格间隔、最大速度及时间采样间隔不符合（4-8）式时，递推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。
频散关系式：
（4-11）
式中
为最小速度，
为Nyquist（奈奎斯特）频率。
一般取震源子波中的主频
的2倍值参与计算，G为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分
取8，而对于空间为四阶差分的情况则
取4方能有效减少频散。
同理：
对于空间微分，采用二阶精度差分格式为
（4-12）
考虑Reynolds边界条件（详细推到见文献“基于Marmousi模型的声波方程有限差分正演算法”），差分格式如下：
（4-13）
注：
其中参数的设置：
借用以前实验的数据，当然可以根据限制条件4-10、4-11计算得到；至于震源放于[101*5,101*5]的矩形中心，根据速度与位移可以计算传播到边界时的时间，此处忽略。
二．实验步骤
1、应用声波方程作为正演模拟的波动方程，忽略转换波的产生、传播；
2、将所提供震源函数离散后绘图；
震源函数为雷克子波，离散绘图如下：
fm=30;
r=3;
t=0.002;
forn=1:
50
w（n）=exp（-（2*pi*fm/r）^2*（t*n）^2）*cos（2*pi*fm*t*n）;
end
subplot（211）
plot（w）
t1=0.002;
forn=1:
50
w1（n）=（1-2*（pi*fm*（t1*n-1/fm））^2）*exp（-（pi*fm*（t1*n-1/fm））^2）;
end
subplot（212）
plot（w1）
3、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即区域内所有网格点的波场值）。
第一时刻为地震波还未传播到边界上的某时刻。
由稳定条件，设v为2000，可以令DT=0.001,DH=5,此时精度较高，且满足频散关系
程序，图形如下：
#include
#include
#include
#definePI3.141593
#defineFM30
#defineR3
#defineKN200
#defineXN101
#defineZN101
#defineDH5
#defineDT0.001
voidfunction（constintflag1,constintflag2）
{
FILE*fp;
inti,j,k,m,n;
floatu1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN];//不能直接初值为0
floatv[XN][ZN],w[KN],uu0,uu1,uu2;
for（i=0;ifor（j=0;j{if（i==50&&j==50）
f[i][j]=1;
else
f[i][j]=0;
}
for（i=0;ifor（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;//速度相同表示同一介质
}
if（0==flag1）
{
for（k=0;kw[k]=exp（-（2*PI*FM/R）*（2*PI*FM/R）*（k*DT）*（k*DT））*cos（2*PI*FM*k*DT）;
}
elseif（1==flag1）
{
for（k=0;kw[k]=（1-2*pow（（PI*FM*（k*DT-1./FM））,2））*exp（-2*pow（（PI*FM*（k*DT-1./FM））,2））;
}
for（k=0;k{
if（0==flag2）//时间二阶差分格式、空间四阶差分格式
{
for（i=2;ifor（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];
}
}
elseif（1==flag2）//时间二阶差分格式、空间二阶差分格式
{
for（i=2;ifor（j=2;j{
uu0=（v[i][j]）*（v[i][j]）*（DT/DH）*（DT/DH）;
uu1=u2[i+1][j]-2*u2[i][j]+u2[i-1][j];
uu2=u2[i][j+1]-2*u2[i][j]+u2[i][j-1];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
}
}
for（m=0;mfor（n=0;n{
u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
if（k==100）
for（m=0;mfor（n=0;nu4[m][n]=u3[m][n];//记录波前快照,中间点
}
if（0==flag1&&0==flag2）
fp=fopen（"wavefront14.dat","w"）;
elseif（1==flag1&&0==flag2）
fp=fopen（"wavefront24.dat","w"）;
elseif（0==flag1&&1==flag2）
fp=fopen（"wavefront12.dat","w"）;
elseif（1==flag1&&1==flag2）
fp=fopen（"wavefront22.dat","w"）;
if（fp!
=NULL）
{
for（i=0;ifor（j=0;jfprintf（fp,"%d%d%f\n",i,j,u4[i][j]）;
fclose（fp）;
}
}
voidmain（）
{
function（0,0）;//时间二阶差分格式、空间四阶差分格式、震源1
function（0,1）;//时间二阶差分格式、空间二阶差分格式、震源1
function（1,0）;//时间二阶差分格式、空间四阶差分格式、震源2
function（1,1）;//时间二阶差分格式、空间二阶差分格式、震源2
}
（1）时间二阶差分格式、空间四阶差分格式、震源1
（2）时间二阶差分格式、空间二阶差分格式、震源1
（3）时间二阶差分格式、空间四阶差分格式、震源2
（4）时间二阶差分格式、空间二阶差分格式、震源2
第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射（此处用Reynolds边界条件）；
程序图形如下：
#include
#include
#include
#definePI3.141593
#defineFM30
#defineR3
#defineKN200
#defineXN101
#defineZN101
#defineDH5
#defineDT0.001
voidfunction（constintflag1,constintflag2）
{
FILE*fp;
inti,j,k,m,n;
floatu1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN];//不能直接初值为0
floatv[XN][ZN],w[KN],uu0,uu1,uu2,uu3;
for（i=0;ifor（j=0;j{if（i==50&&j==50）
f[i][j]=1;
else
f[i][j]=0;
}
for（i=0;ifor（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;//速度相同表示同一介质
}
if（0==flag1）
{
for（k=0;kw[k]=exp（-（2*PI*FM/R）*（2*PI*FM/R）*（k*DT）*（k*DT））*cos（2*PI*FM*k*DT）;
}
elseif（1==flag1）
{
for（k=0;kw[k]=（1-2*pow（（PI*FM*（k*DT-1./FM））,2））*exp（-2*pow（（PI*FM*（k*DT-1./FM））,2））;
}
for（k=0;k{
if（0==flag2）//时间二阶差分格式、空间四阶差分格式
{
for（i=2;ifor（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];
}
}
elseif（1==flag2）//时间二阶差分格式、空间二阶差分格式
{
for（i=2;ifor（j=2;j{
uu0=（v[i][j]）*（v[i][j]）*（DT/DH）*（DT/DH）;
uu1=u2[i+1][j]-2*u2[i][j]+u2[i-1][j];
uu2=u2[i][j+1]-2*u2[i][j]+u2[i][j-1];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
}
}
for（m=0;mfor（n=0;n{
u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
//Reynolds边界条件
uu3=v[i][j]*（DT/DH）;
u3[0][j]=u2[0][j]+u2[1][j]-u1[1][j]+uu3*（u2[1][j]-u2[0][j]-u1[2][j]+u1[1][j]）;//左边界
u3[XN][j]=u2[XN][j]+u2[XN-1][j]-u1[XN-1][j]+uu3*（u2[XN-1][j]-u2[XN][j]-u1[XN-2][j]+u1[XN-1][j]）;//右边界
u3[i][0]=u2[i][0]+u2[i][1]-u1[i][1]+uu3*（u2[i][1]-u2[i][0]-u1[i][2]+u1[i][1]）;//下边界
u3[i][ZN]=u2[i][ZN]+u2[i][ZN-1]-u1[i][ZN-1]+uu3*（u2[i][ZN-1]-u2[i][ZN]-u1[i][ZN-2]+u1[i][ZN-1]）;//上边界
if（k==180）//只需改动K值，即显示边界反射
for（m=0;mfor（n=0;nu4[m][n]=u3[m][n];//记录波前快照，边界
}
if（0==flag1&&0==flag2）
fp=fopen（"wavefront14.dat","w"）;
elseif（1==flag1&&0==flag2）
fp=fopen（"wavefront24.dat","w"）;
elseif（0==flag1&&1==flag2）
fp=fopen（"wavefront12.dat","w"）;
elseif（1==flag1&&1==flag2）
fp=fopen（"wavefront22.dat","w"）;
if（fp!
=NULL）
{
for（i=0;ifor（j=0;jfprintf（fp,"%d%d%f\n",i,j,u4[i][j]）;
fclose（fp）;
}
}
voidmain（）
{
function（0,0）;//时间二阶差分格式、空间四阶差分格式、震源1
function（0,1）;//时间二阶差分格式、空间二阶差分格式、震源1
function（1,0）;//时间二阶差分格式、空间四阶差分格式、震源2
function（1,1）;//时间二阶差分格式、空间二阶差分格式、震源2
}
（1）时间二阶差分格式、空间四阶差分格式、震源1
（2）时间二阶差分格式、空间二阶差分格式、震源1
（3）时间二阶差分格式、空间四阶差分格式、震源2
（4）时间二阶差分格式、空间二阶差分格式、震源2
4、对于大模型，定义为水平层状速度模型；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律，指出哪是反射波，哪是透射波；
这时取小模型的常量，为减少频散，速度v至少为2400
程序图形如下：
#include
#include
#include
#definePI3.141593
#defineFM30
#defineR3
#defineKN200
#defineXN200
#defineZN100
#defineDH5
#defineDT0.001
voidfunction（constintflag1,constintflag2）
{
FILE*fp;
inti,j,k,m,n;
floatu1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN];//不能直接初值为0
floatf[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;
if（0==flag1）
{
for（k=0;kw[k]=exp（-（2*PI*FM/R）*（2*PI*FM/R）*（k*DT）*（k*DT））*cos（2*PI*FM*k*DT）;
}
elseif（1==flag1）
{
for（k=0;kw[k]=（1-2*pow（（PI*FM*（k*DT-1./FM））,2））*exp（-2*pow（（PI*FM*（k*DT-1./FM））,2））;
}
for（i=0;ifor（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{
if（0==flag2）//时间二阶差分格式、空间四阶差分格式
{
for（i=2;ifor（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];
}
}
elseif（1==flag2）//时间二阶差分格式、空间二阶差分格式
{
for（i=2;ifor（j=2;j{
uu0=（v[i][j]）*（v[i][j]）*（DT/DH）*（DT/DH）;
uu1=u2[i+1][j]-2*u2[i][j]+u2[i-1][j];
uu2=u2[i][j+1]-2*u2[i][j]+u2[i][j-1];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
}
}
for（m=0;mfor（n=0;n{
u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
if（k==100）
for（m=0;mfor（n=0;nu4[m][n]=u3[m][n];//记录波前快照
}
if（0==flag1&&0==flag2）
fp=fopen（"wa
升级会员 