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

声波方程数值模拟实验报告学号：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;for n=1:50w(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n);endsubplot(211)plot(w)t1=0.002;for n=1:50w1(n)=(1-2*(pi*fm*(t1*n-1/fm)2)*exp(-(pi*fm*(t1*n-1/fm)2);endsubplot(212)plot(w1)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.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /不能直接初值为0 float vXNZN,wKN,uu0,uu1,uu2;for(i=0;iXN;i+) /定义f函数，当且仅当i,j同时为50时，f为1，其余为0for(j=0;jZN;j+)if(i=50&j=50)fij=1;else fij=0; for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0;u2ij=0.0;u3ij=0.0;u4ij=0.0;vij=2000; /速度相同表示同一介质 if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(k=0;kKN;k+) if(0=flag2) /时间二阶差分格式、空间四阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; else if(1=flag2) /时间二阶差分格式、空间二阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=u2i+1j-2*u2ij+u2i-1j; uu2=u2ij+1-2*u2ij+u2ij-1; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(m=0;mXN;m+) for(n=0;nZN;n+) u1mn=u2mn; u2mn=u3mn; if(k=100) for(m=0;mXN;m+) for(n=0;nZN;n+) u4mn=u3mn;/记录波前快照,中间点 if(0=flag1 & 0=flag2)fp=fopen(wavefront14.dat,w);else if(1=flag1 & 0=flag2)fp=fopen(wavefront24.dat,w);else if(0=flag1 & 1=flag2)fp=fopen(wavefront12.dat,w);else if(1=flag1 & 1=flag2)fp=fopen(wavefront22.dat,w); if(fp!=NULL) for(i=0;iXN;i+)for(j=0;jZN;j+)fprintf(fp,%d %d %fn,i,j,u4ij);fclose(fp); void main()function(0,0); /时间二阶差分格式、空间四阶差分格式、震源1function(0,1); /时间二阶差分格式、空间二阶差分格式、震源1function(1,0); /时间二阶差分格式、空间四阶差分格式、震源2function(1,1); /时间二阶差分格式、空间二阶差分格式、震源2（1）时间二阶差分格式、空间四阶差分格式、震源1（2）时间二阶差分格式、空间二阶差分格式、震源1（3）时间二阶差分格式、空间四阶差分格式、震源2（4）时间二阶差分格式、空间二阶差分格式、震源2第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射（此处用Reynolds边界条件）；程序图形如下：#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.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /不能直接初值为0 float vXNZN,wKN,uu0,uu1,uu2,uu3; for(i=0;iXN;i+) /定义f函数，当且仅当i,j同时为50时，f为1，其余为0for(j=0;jZN;j+)if(i=50&j=50)fij=1;else fij=0; for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0;u2ij=0.0;u3ij=0.0;u4ij=0.0;vij=2000; /速度相同表示同一介质 if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(k=0;kKN;k+)if(0=flag2) /时间二阶差分格式、空间四阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; else if(1=flag2) /时间二阶差分格式、空间二阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=u2i+1j-2*u2ij+u2i-1j; uu2=u2ij+1-2*u2ij+u2ij-1; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(m=0;mXN;m+) /波向前传播记录连续三个时刻的值for(n=0;nZN;n+)u1mn=u2mn; u2mn=u3mn; /Reynolds边界条件uu3=vij*(DT/DH);u30j=u20j+u21j-u11j+uu3*(u21j-u20j-u12j+u11j);/左边界u3XNj=u2XNj+u2XN-1j-u1XN-1j+uu3*(u2XN-1j-u2XNj-u1XN-2j+u1XN-1j);/右边界u3i0=u2i0+u2i1-u1i1+uu3*(u2i1-u2i0-u1i2+u1i1);/下边界u3iZN=u2iZN+u2iZN-1-u1iZN-1+uu3*(u2iZN-1-u2iZN-u1iZN-2+u1iZN-1);/上边界if(k=180) /只需改动K值，即显示边界反射for(m=0;mXN;m+)for(n=0;nZN;n+)u4mn=u3mn;/记录波前快照，边界 if(0=flag1 & 0=flag2)fp=fopen(wavefront14.dat,w);else if(1=flag1 & 0=flag2)fp=fopen(wavefront24.dat,w);else if(0=flag1 & 1=flag2)fp=fopen(wavefront12.dat,w);else if(1=flag1 & 1=flag2)fp=fopen(wavefront22.dat,w); if(fp!=NULL) for(i=0;iXN;i+)for(j=0;jZN;j+)fprintf(fp,%d %d %fn,i,j,u4ij);fclose(fp); void main()function(0,0); /时间二阶差分格式、空间四阶差分格式、震源1function(0,1); /时间二阶差分格式、空间二阶差分格式、震源1function(1,0); /时间二阶差分格式、空间四阶差分格式、震源2function(1,1); /时间二阶差分格式、空间二阶差分格式、震源2（1）时间二阶差分格式、空间四阶差分格式、震源1（2）时间二阶差分格式、空间二阶差分格式、震源1（3）时间二阶差分格式、空间四阶差分格式、震源2（4）时间二阶差分格式、空间二阶差分格式、震源24、 对于大模型，定义为水平层状速度模型；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律，指出哪是反射波，哪是透射波；这时取小模型的常量，为减少频散，速度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.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN; /不能直接初值为0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2;if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0; u2ij=0.0; u3ij=0.0; u4ij=0.0; fij=0.0; if (j=30) vij=2400; /第一层速度为2400 else vij=3000; /第二层速度为3000 f10010=1; /定义f函数，在100，10为1 for(k=0;kKN;k+) if(0=flag2) /时间二阶差分格式、空间四阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; else if(1=flag2) /时间二阶差分格式、空间二阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=u2i+1j-2*u2ij+u2i-1j; uu2=u2ij+1-2*u2ij+u2ij-1; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; for(m=0;mXN;m+) for(n=0;nZN;n+) u1mn=u2mn; u2mn=u3mn; if(k=100)for(m=0;mXN;m+)for(n=0;nZN;n+)u4mn=u3mn;/记录波前快照 if(0=flag1 & 0=flag2)fp=fopen(wavefront14.dat,w);else if(1=flag1 & 0=flag2)fp=fopen(wavefront24.dat,w);else if(0=flag1 & 1=flag2)fp=fopen(wavefront12.dat,w);else if(1=flag1 & 1=flag2)fp=fopen(wavefront22.dat,w); if(fp!=NULL) for(i=0;iXN;i+)for(j=0;jZN;j+)fprintf(fp,%d %d %fn,i,j,u4ij);fclose(fp); void main()function(0,0); /时间二阶差分格式、空间四阶差分格式、震源1function(0,1); /时间二阶差分格式、空间二阶差分格式、震源1function(1,0); /时间二阶差分格式、空间四阶差分格式、震源2function(1,1); /时间二阶差分格式、空间二阶差分格式、震源2（1）时间二阶差分格式、空间四阶差分格式、震源1（2）时间二阶差分格式、空间二阶差分格式、震源1（3）时间二阶差分格式、空间四阶差分格式、震源2（4）时间二阶差分格式、空间二阶差分格式、震源2二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波。这时取小模型的常量，为减少频散，速度v至少为2400程序图像如下：#include #include #include #define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 200#define ZN 180#define DH 5#define DT 0.001void function(const int flag1,const int flag2)FILE *fp;int i,j,k,m,n; float u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNKN; /不能直接初值为0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2;if(0=flag1)for(k=0;kKN;k+)wk=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT)*cos(2*PI*FM*k*DT);else if(1=flag1)for(k=0;kKN;k+)wk=(1-2*pow(PI*FM*(k*DT-1./FM),2)*exp(-2*pow(PI*FM*(k*DT-1./FM),2); for(i=0;iXN;i+) for(j=0;jZN;j+) u1ij=0.0; u2ij=0.0; u3ij=0.0; u4ij=0.0; fij=0.0; if (j=30) vij=2400; /第一层速度为2400 else vij=3000; /第二层速度为3000 for(i=0;iXN;i+) for(k=0;kKN;k+) u5ik=0.0; /速度相同表示同一介质 f10010=1; /定义f函数，在100，10为1 for(k=0;kKN;k+) if(0=flag2) /时间二阶差分格式、空间四阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)*(DT/DH); uu1=-1.0/12*(u2i-2j+u2i+2j)+4.0/3*(u2i-1j+u2i+1j)-5.0/2*u2ij; uu2=-1.0/12*(u2ij-2+u2ij+2)+4.0/3*(u2ij-1+u2ij+1)-5.0/2*u2ij; u3ij=2*u2ij-u1ij+uu0*uu1+uu0*uu2+wk*fij; u5ik=u3i10; else if(1=flag2) /时间二阶差分格式、空间二阶差分格式for(i=2;iXN-2;i+)for(j=2;jZN-2;j+) uu0=(vij)*(vij)*(DT/DH)