规则网格有限差分解声波方程个人总结报告

地球探测科学与技术学院 总 结 报 告学 校：吉林大学学 院：地球探测科学与技术学院专 业：勘查技术与工程（应用地球物理）科 目：科学计算方法-有限差分解声波方程姓 名： 学 号： 目录一 相关理论基础31. 地震波场模拟32. 波动方程类型及其局限性33. 数值算法类型及其优缺点4二有限差分解声波方程基础理论知识61.需要的已知条件包括：62弹性波方程63.声波方程的有限差分法数值模拟64. 稳定性条件75. 频散关系式86. 有限差分参数8三程序及结果成图8四通过实验所发现的问题和认识12五他人所做的有限差分解波动方程程序及结果成图12 参考文献及资料19有限差分解声波方程总结报告1 相关理论基础1. 地震波场模拟地震波场模拟即地震正演，是指已知模型结构，通过物理或数值计算的方法模拟该地质结构下的地震波的传播，最终合成地震记录，也可以认为其是野外数据采集过程的室内再现。物理模拟花费昂贵，人们一般采用比较经济的数值模拟技术。地震波场数值模拟是在给定数学模型（如弹性波方程，声波方程等）、震源和地下几何界面、物性参数（岩层密度、速度等）情况下，研究弹性波或声波的传播规律。2. 波动方程类型及其局限性（1） 声波方程：二阶标量声波方程：一阶压力-速度方程组：波动方程能够描述且只能描述纵波的传播规律，包括直达波、反射波、透射波、折射波等，但不能描述转换波传播规律。需要的已知条件包括：震源函数、地层速度、密度边界条件（2） 弹性波方程：弹性波方程能够描述纵、横波的传播规律，包括直达波、反射波、透射波、折射波以及转换波等。需要的已知条件包括：震源函数、地层速度或根据方程的类型需要提供的地层的其它弹性参数、边界条件。（3） 粘声波/弹性波方程前面讨论的是理想弹性介质，波在其中传播时，没有能量的损耗，介质中应力和应变关系严格遵循胡克定律（这种理想介质称虎克固体），但波在实际介质中传播时，是有能量损耗的，这就是所谓的弹性波吸收。波在传播过程中，实际介质的不同部位之间会出现某种摩擦力，称为内摩擦力或粘滞力。这种力导致机械能向其他形式能量转换，最终转化为热能消耗掉。在地震勘探中，地震波传播的实际介质是十分复杂的。在一定条件下，即震源作用时间短，作用力微小，地球介质可以看作完全弹性模型，但随着地震勘探技术的发展，勘探精度要求提高，面临复杂地质目标时，要求地震勘探采用更加符合实际的介质模型进行研究。粘弹性介质模型更符合实际。但是到目前为止，在地震资料反演处理中应用最多的还是声波方程，弹性波以及粘弹性波方程的应用还只是停留在模拟层次上。3. 数值算法类型及其优缺点地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。（1） 克希霍夫积分法 本质是波动方程积分解的一个数值计算，某种程度上相当于绕射叠加。 特点该方法计算速度较快，但由于射线追踪中存在着焦散、多重路径等问题，故其一般只能适合于较简单的模型，难以模拟复杂地层的波场信息。（2） 傅里叶变换法 原理利用空间的全部信息对波场函数进行三角函数插值，能更加精确地模拟地震波的传播规律，同时，利用快速傅里叶变换(FFT)进行计算，还可以提高运算效率 特点精度高，占用内存小，但缺点是计算速度较慢，对模型的适用性差，尤其是不适应于速度横向变化剧烈的模型。（3） 波动方程有限元法 原理将变分法用于单元分析，得到单元矩阵，然后将单元矩阵总体求和得到总体矩阵，最后求解总体矩阵得到波动方程的数值解 优点理论上可适宜于任意地质体形态的模型，保证复杂地层形态模拟的逼真性，达到很高的计算精度 缺点占用内存和运算量均较大，不适用于大规模模拟，因此该方法在地震波勘探中尚未得到广泛地应用。（4） 有限差分法 基本原理该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分方法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的函数值为未知数的代数方程组。 特点该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 主要内容v 建立地球物理问题的离散有限差分模型* 根据问题的特点将定解区域做网格划分* 在所有网格节点上用有限差分格式对导数求近似，对函数、初始条件和边界条件求近似* 把原方程离散化为代数方程组v 保证计算过程的可行性和计算结果的正确性：解的相容性、稳定性、收敛性v 数值求解差分方程组二有限差分解声波方程基础理论知识1.需要的已知条件包括：（1）震源函数（2）地层速度（波速）（3）边界条件2弹性波方程：3.声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述： （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)同理，对于空间微分，采用同样的二阶精度差分格式，即有： (4-5) (4-6)这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得： (4-7)式中为介质速度的空间离散值，是空间离散步长，为时间离散步长，4. 稳定性条件： (4-8)这里表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-8)式时，第推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。5. 频散关系式： (4-9)式中为最小速度，为Nyquist频率。一般取震源子波中的主频的2倍值参与计算，G为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分取8，而对于空间为四阶差分的情况则取4方能有效减少频散。6. 有限差分参数：（1）震源：雷克子波， f=25Hz t=0.001s (4-10)（2）模型参数：各向同性介质，v=2000ms，大小 100*100m2（3）空间离散步长：dx=dz=1m（4）时间离散步长：dt=0.001s三程序及结果成图1.震源函数为雷克子波，离散绘图如下：clear;fm=25;t=0.001;for n=1:200f(n)=(1-2*pi*pi*fm*fm*n*n*t*t)*exp(-pi*pi*fm*fm*n*t*n*t);%雷克子波 endfigure(1),plot(w); 图1 震源雷克子波（延续时间200ms）2.设定模型为各向均匀介质，将震源点放在模型中间，震源位置，分别记录两个时刻的波前快照（即区域内所有网格点的波场值）：（1） 第一时刻为地震波还未传播到边界上的某时刻：eg t=500、2000ms（2） 第二时刻为地震波传播到人工边界上的某时刻：eg t=3000ms程序代码：clc;clear;v=2000; %速度dt=0.0001; %时间步长tm=0.8; %峰值，单位秒 fm=25; %主频，单位Hzdx=1;dz=1; %方向步长Sx=50;Sz=50;%震源x=100;z=100;%网格tmp=4.0000e-004; %tmp=v2*dt2/dx2p1=zeros(100,100);p1(:,:)=1;p2=zeros(100,100); p2(:,:)=1; for nt=1:1:7000 t=nt*dt; f=(1-2*pi*pi*fm*fm*t*t)*exp(-pi*pi*fm*fm*t*t);%雷克子波 for i=2:dx:99 for j=2:dz:99 p3(i,j)=tmp*(p2(i+1,j)-4*p2(i,j)+p2(i-1,j)+p2(i,j+1)+p2(i,j-1)+2*p2(i,j)-p1(i,j); end end if nt=500 %可换为t=1000/ imagesc(p3); figure(gcf); end for i=2:dx:99 for j=2:dz:99 p1(i,j)=p2(i,j); p2(i,j)=p3(i,j); p2(Sx,Sz)=p2(Sx,Sz)+f; end end end波场图件： 图2 t=500ms时波场图 图3 t=2000ms时波场图 图4 t=3000ms时波场图 图5 t=5000ms时波场图 图件分析：从实验的四张图中可看出波场向外传播的方式，但要注意超过一定时间后就会产生虚假反射波，是由于人工边界造成的，解决方案就是设计吸收边界。注意频散现象，这是由于所选择的差分公式阶数太低造成的，本组采用空间和时间都为二阶的差分公式，所以计算精度不够，要解决这个问题，可以采用高阶差分提高精度并且选用衰减较快的子波震源。 四通过实验所发现的问题和认识（1） 实验过程中深刻体会到了编程能力的匮乏，采用c语言编程时，很多小问题不容忽视:大括号、变量申明、数据读写、输入法等等，还有借助于matlab成图时要注意引用数据文件设置成全英文等等。（2） 今后的学习中一定要注意积累编程经验，学习编程语言（3） 相比之下，matlab编程要更容易、快捷五他人所做的有限差分解波动方程程序及结果成图1. 波动方程：2. 参数 震源参数：雷克子波，f=30,r=3 模型1参数：大小100*100m2，各向同性介质，速度2000ms 有限差分参数：dx=dz=5m，dt=0.001s3. 结果成图及程序代码（1）不考虑吸收边界的波场试验 图 6 t=100ms 时的波场图 图7 t=160ms 时的波场图程序代码#include #include #include #define PI 3.#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 u1XNZN,u2XNZN,u3XNZN,u4XNZN,fXNZN; /不能直接初值为0 float u5XNZN,vXNZN,wKN,uu0,uu1,uu2; 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); for(i=0;iXN;i+) /定义f函数，当且仅当i,j同时为50时，f为1，其余为0 for(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; /速度相同表示同一介质 for(k=0;kKN;k+) 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; 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(fp=fopen(wavefront.txt,w)!=NULL) for(i=0;iXN;i+) for(j=0;jZN;j+) fprintf(fp,%ft,u4ij); fprintf(fp,n); fclose(fp); （2） 水平层状介质的波前快照图8 水平双层介质的波场快照图程序代码：#include #include #include #define PI 3.#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 u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNZN; /不能直接初值为0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2; 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); 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+) 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=u1i10; /地震记录 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(fp=fopen(wavefront.dat,w)!=NULL) fprintf(fp,%dn,XN); fprintf(fp,%dn,ZN); for(i=0;iXN;i+) for(j=0;jZN;j+) fprintf(fp,%fn,u4ij); fclose(fp); （3） 水平层状介质的合成地震记录 图9 水平双层介质的拟合地震记录图程序代码：#include #include #include #define PI 3.#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 u1XNZN,u2XNZN,u3XNZN,u4XNZN,u5XNZN; /不能直接初值为0 float fXNZN,vXNZN,wKN,uu0,uu1,uu2; 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); 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