有限差分纵波波场模拟

1、弹性波课程报告 有限差分纵波波场模拟姓名： 邓 健 学号： 1201410215 课程名称：弹性波理论 任课老师：顾 汉 明 专业：地球探测与信息技术 时间： 2015年5月 评 语对课程论文的评语:平时成绩：课程论文成绩：总 成 绩：评阅人签名：注：1、无评阅人签名成绩无效；2、必须用钢笔或圆珠笔批阅，用铅笔阅卷无效；3、如有平时成绩，必须在上面评分表中标出，并计算入总成绩。摘要地震波场数值模拟是勘探地震学的重要研究课题之一,也是研究波动现象的重要手段。地震波场数值模拟常用的方法有伪谱法、有限差分法、有限元法等。有限差分法具有计算速度快、占用内存小等优点，该方法对于近远场及复杂边界都有广泛的

2、适用性，能够准确地模拟波在各种介质及复杂结构地层中的传播规律，是勘探地震学中应用最广泛的数值计算方法。本文以波动理论为理论基础，利用泰勒级数展开式推导出波动方程的有限差分格式及其离散表达式。针对传统二阶精度差分方法模拟精度较低的问题，推导出时间域二阶空间域四阶精度的有限差分方程，并综合分析了初始条件、震源项、算法稳定性及数值频散等因素在有限差分法数值模拟中的影响。有限差分数值模拟计算是在一个有限的区域内进行的，当波传播到边界时就会产生边界反射，这是进行数值模拟计算时所不期望出现的。本报告通过调研大量文献资料，对效果较好的透明边界条件进行了详细分析。在此基础上,分别推导得到这种边界条件方法的离散

3、表达式，并在数值模拟过程中进行验证。通过建立均匀模型、层状模型及高速透镜体模型，基于Matlab语言编程，论文实现了波动方程有限差分算法的数值模拟。对不同地质理论模型进行模拟，从模拟结果的对比分析中可以看出，模型参数的选择及各参数间的相互关系对模拟结果有着显著影响。无论是模拟精度，还是模拟计算效率，有限差分算法都具有一定优势。通过综合研究，论文认为波动方程有限差分算法具有算法简单、计算效率高、模拟精度较高等特点，理论研究意义大，应用前景广阔。关键词：数值模拟，波动方程，有限差分，边界条件目录1 前言11.1 地震波场数值模拟概述11.2 有限差分法及其与几种常用数值模拟方法的比较12 波动方程

4、有限差分法数值模拟22.1 有限差分原理22. 2 波动方程的建立33 二维声波方程有限差分格式的建立63.1 声波方程63.2 波动方程有限差分格式的建立73.3 高阶差分方程94 波动方程有限差分的几个问题124.1 初始条件124.2 震源函数124.3 差分方程的稳定性及收敛性134.4 频散问题135 有限差分算法中边界条件的处理145.1 一维透明边界条件155.2 二维透明边界条件165.3 透明边界条件的差分格式186 简单模型试算206.1均匀模型206.2 水平层状模型和高速透镜体模型227结论与建议25参考文献26附录281 前言1.1 地震波场数值模拟概述地震波场数值模

5、拟简单说来，就是已知地下介质构造及其参数，再用理论计算方法来研究地震波在地下介质中的传播规律，合成地震记录的一种技术方法。随着地震勘探技术的发展，数值模拟方法已经贯穿于地震数据的采集、处理和解释的全过程，而且在确定观测的合理性、检验处理和解释的正确性等方面都有了广泛的应用1-3。地震波数值模拟问题的研究主要包括三个方面：(1)地震波场数值模拟原理；(2)地震波场数值模拟算法；(3)数值模拟的计算机实现。地震波场数值模拟方法的理论基础主要是波动理论和射线理论。射线理论方法是求出质点运动方程的渐进近似解，这种方法计算速度快，但是精度较低。而波动理论方法则是用数值计算方法直接求出波动方程的解，其模拟

6、结果中包含丰富的波动信息，模拟结果较为精确。数值模拟又可分为波动方程的解析解和数值解。对简单的地质模型就可以得到其精确解，但对于复杂的地质模型，只有通过数值解来说明波在底下介质中的传播。通常解析解只是作为数值解的一种检验4-6。1.2 有限差分法及其与几种常用数值模拟方法的比较在地震勘探中，为了研究地震波在地下各种介质中的传播规律，就需要对波场进行数值模拟。常用的数值模拟方法有伪谱法、有限差分法，有限元法等。有限元法依据变分原理，通过灵活的网格剖分，用简单形态逼近实际的地质体，能处理多种介质和自然边界条件。适用于非规则网格的计算，方便有效，模拟的精度高。但是有限元法的问题是不适用于大规模的模型

7、的计算，而且计算量大7。伪谱法在20世纪七十年代引入数值模拟计算领域的，它是利用傅立叶变换来计算波场的空间导数，用差分方法来计算时间的导数，可以看成是一种无限阶的有限差分法，是传统有限差分法的一个推广。伪谱法在粗网格上也能实现高精度的计算，相对有限元法实现起来较容易，在非线性波动问题及气象预测等领域有着广泛的应用8-10。有限差分法是一种最常用的数值模拟方法，它是将波动方程中波场函数的空间导数和时间导数用相应的空间、时间的差分代替。有限差分法具有计算速度快、占用内存小等优点，该方法对于近远场及复杂边界都有广泛的适用性，能够准确地模拟波在各种介质及复杂结构地层中的传播规律。有限差分法方法简单、高

8、效的优点是其他方法难以比拟的，因此有限差分法目前仍然是勘探地震学中应用最广泛的数值计算方法11。2 波动方程有限差分法数值模拟2.1 有限差分原理有限差分法就是把求解区域划分为差分网格，然后用有限的网格节点代替连续的求解域，利用微商与差商的近似关系将描述介质传播的微分方程转化为差分方程进行求解。差分离散的方法有两种，一种是将单变量的二阶波动方程直接转化为时间空间的二阶中心差分进行离散求解；第二种方法是把用位移表示的二阶波动方程转化为用应力及质点速度表示的一阶方程组，然后用应力和速度的交错网格求解,8。构造差分的方法也有很多种，一般采用泰勒级数展开法。常见的差分格式有显示差分、隐式差分和显隐交替

9、格式；按照差分的精度又可以分为一阶差分、二阶差分和高阶差分等，目前通常见到的差分格式主要是几种差分格式的组合。我们知道，差分方程的建立首先选择网格布局和差分形式，然后以有限差分代替无限微分，以差商代替微商，并以差分方程代替微分方程及其边界条件，最后建立差分方程12。在建立差分方程的时候要注意到两个方面。一是合理选择网格布局与步长。我们将离散后各相邻离散点之间的距离或者离散化单元的长度称为步长。在所选定的区域内进行网格划分是差分方程建立的第一步，其方法比较灵活，但是实际应用中往往遵循误差最小原则。常用的典型的网格剖分方法有矩形网格、三角形网格等，网格样式的选择和区域的形状有关。其次是将微分方程转

10、化为差分方程。这个过程就是选择一种差分格式代替微分形式的过程。构造差分的方法有多种形式，本文采用的是泰勒级数展开法2,5,13。地震波场数值模拟以地震波动理论为基础，用有限差分法解波动方程时，对变量离散化，也即对连续的物理量只考虑其在离散的空间位置和离散时刻的值，然后把方程中的导数用这些离散的采样值表示14。对于一个单变量的函数f(x)，将其离散化，那么在采样点x =lx的采样值就是f(lx),其中x表示步长，l为整数。则有限差分法中f(x)在采样点x=lx的导数就可以近似表示为： 其中an是系数，N表示差分格式的长度，差分的格式是由这两个系数来决定。在实际应用中常用的差分格式有向前差分、向后

11、差分以及中心差分：（1）向前差分：（2）向后差分：（3）中心差分：2. 2 波动方程的建立建立波动方程是在已知物体形状、位置、弹性常数及外力分布等参数的情况下求出物体的位移、应力与应变的分布。这个问题具体包括以下内容：（1）应力部分的平衡方程（2）应变部分（3）应力与应变的关系将式（2-7）代入（2-6），推导可以得下式：方程（2-8）是物体在平衡状态下的平衡方程。当物体处于不平衡状态时，式（2-8）则变为：在二维情况下，我们只需要考虑x，z方向的位移分量，这时式（2-9）就可转换为关于x,z的方程：当fx=fz=0时，上面的这两个方程就变为：（2-11）就是二维非均匀介质弹性波方程。而在均匀

12、介质中，、和均为常数，则二维均匀介质弹性波方程就为：把（2-12）中的x和z分量合并就得到了二维均匀介质弹性波方程的矢量形式表达式在没有弹性横波只有弹性纵波存在时，对（2-13）两边取散度：利用旋度与散度的关系，交换上式的微分次序并化简可得：其中表示纵波速度。根据赫姆霍兹在他著名的有关涡流运动的著作中证明了下属定理：任何向量点函数，若它的散度和旋度具有位，则它可以表示为一个无旋部分和一个旋转部分之和。亦即对于任何一种向量场，如果在其定义域内有散度和旋度，则该向量场可表示为一个标量位的梯度场和一个向量位的旋度场之和，而空间传播的波动正是无旋运动和旋转运动这两种运动之和15-20 令代入式（2-8

13、）中，取其中的无旋部分则有：然后再消去式中的梯度,就得到了用位移位表示的纵波方程：3 二维声波方程有限差分格式的建立3.1 声波方程一般地，二维非均匀介质的声波波动方程可表示为：式中 U = U ( x , z , t)表示声压，V表示声波在介质中的传播速度，是密度，它是随空间各点而变的。其中 s ( x , z , t )为震源函数。对于均匀介质，密度为常数，则二维声波波动方程就可以表示为：在实际问题中，我们总是用一个有限宽频带的时间函数来代替 s ( x , z , t )函数，以便能够真实地反映地震波的传播21。3.2 波动方程有限差分格式的建立上一节中讨论了差分的原理和几种常见的差分格

14、式，在这个基础上，我们结合波动方程，对纵波波动方程的有限差分格式进行推导。令，其中x，z是空间间隔，t是时间间隔。用k表示时间方向的离散网格，m表示x方向的离散网格，n表示z方向的离散网格。利用泰勒级数展开式将 在 展开可得：同理，在处的展开式为：用（3-3）式减去（3-4）式，就可以推出关于t的一阶中心差分：如果将（3-3）式与（3-4）式相加可得进一步推导就可以得到关于t的二阶中心差分：同理也可以推导出关于x，z的中心差分格式：关于x的一阶中心差分关于x的二阶中心差分： 关于z的一阶中心差分：关于z的二阶中心差分：若令x =z=h，利用上面关于x,z,t的中心差分方程就可以得到二维波动方程

15、的有限差分方程：对上式继续推导可得：其中，上式就是二阶波动方程的有限差分格式。如果我们继续利用泰勒级数展开式，则可以得到：这样我们就能够得到时间域二阶空间域四阶的波动方程有限差分格式：上式中k表示时间方向的离散网格，m表示x方向的离散网格，n表示z方向的离散网格，。需要注意的是，这里同样也假设x=z=h，即网格步长相等5 ,22-24。3.3 高阶差分方程将采用高阶差分会有很多时间层，计算较复杂。为此时间上采用2阶，空间上可采用高阶。为了方便起见，本文空间步长在x，z轴上相等，则通过上述类似推倒，可以得到不同精度的声波差分方程。时间2阶，空间6阶时间2阶，空间8阶时间2阶，空间10阶时间2阶，

16、空间12阶时间2阶，空间14阶图3-1为均匀介质中各阶精度波场快照对比图。模型参数为网格500×500，震源位置（250，250），空间步长x和z方向均为1m,采样时间间隔100us。所用的震源是雷克子波，子波频率30Hz。从下图可以看出，时间2阶，空间2阶精度波场出现明显频散，随着空间精度的增高，频散逐渐减弱。所以提高差分精度可以有效减轻频散，但计算量也随之增加。图3-1均匀介质中各阶精度波场快照图4 波动方程有限差分的几个问题4.1 初始条件在进行波动方程有限差分数值模拟时，初始时刻的波场值是不知道的，而在计算差分方程时，只有知道时间间隔前一时刻的值才能递推计算出后面时刻的值。因

17、此在计算开始时就必须要考虑到初始条件。一般假设在震源发生震动前，各个质点均处于静止状态：初始时刻波场的速度也同样为零，即4.2 震源函数震源函数的给出通常有两种方式：一种是初值法，就是将理论结果当作初始值给出；另外一种是力源法，就是以力源的方式给出。初值法的震源除了自由表面或内界面附近，可以在模型的其他任意处进行定义；而力源法可以将震源定义在自由表面的附近，但必须是在网格点上。两种方法各有其优缺点1,26。在进行波动方程数值模拟计算的过程中，震源问题的处理对模拟的结果有着重要的影响。子波是描述震源的时间延续特征的时间函数，对于地震子波而言，子波延续的时间越短，频带越宽，子波的垂直分辨率也就越高

18、。由于做有限差分计算时会出现数值频散，尤其当空间采样不足时，子波的高频成分频散就会更严重。所以要根据模型的速度及网格间距合理选择子波主频。本文采用了雷克子波和高斯子波作为震源函数27。4.3 差分方程的稳定性及收敛性一个有实际应用意义的数值模拟算法必须是稳定的，也就是说对算法的数值解与解析解的差值是有限定的。对于有限差分法，一般来说如果差分方程的解的误差不随时间的推进而增加，那么就说该差分方程的解是稳定的，这也就是说差分方程的解是有界的。差分的稳定条件根据差分阶数的不同而有所不同，差分的阶数越高，稳定性越差。对于不同的差分格式，参数的选取一般也是不同的5,28,29。对于均匀介质，差分方程的解

19、的稳定性是确定的，稳定性条件如下式所示：其中V是均匀介质中波的传播速度，应当注意，这里的x是差分格点中x与z的最小值。而对于非均匀介质应满足：其中Vmax是各种均匀介质中的最大波速值。4.4 频散问题在进行地震波数值模拟时，我们希望尽可能同时达到较高的模拟精度和较快的计算速度。但是在用微分方程数值模拟的方法求解波动方程时，就会产生数值频散或网格频散。这种频散是用微分方程求解波动方程时固有的本质特征，是避免不了的，只能尽量减弱这种频散2。通过对有限差分法数值频散的理论分析研究可以知道，影响频散的因素有很多。在有限差分的计算过程中，不仅需要对空间进行网格化离散，同时还需要对时间进行离散。波动方程经

20、过空间离散化之后就引入了差分算子误差；而时间间隔的选择也对数值模拟的精度和数值频散有一定影响。当子波频率越高，则一个波场内的网格点数就越少，则频散现象越严重；对于一定频率的子波，介质的速度越低，则一个波场内的网格点数也少，数值频散现象也越严重。因此，从理论的角度来说，在进行数值模拟时，对于给定的子波频率，提高时间和空间的差分阶数，减小网格步长和时间间隔都可以提高数值计算的精度和改善数值频散情况。但是过细的网格剖分会增大计算所需的存储量，从而增加计算的时间，所以就需要选取合适的参数，在保证精度和计算速度的同时尽量减少频散29-33。5 有限差分算法中边界条件的处理我们在计算机上进行有限差分数值模

21、拟计算时，由于受到计算机内存和计算速度等因素的制约，只能在一个有限的区域内进行计算。这样就会在两侧与底界上产生边界，这个边界就是人工边界。当波传播到人工边界时就会产生反射，反射必然会造成一定的干扰，产生边界效应。从图5-1（左）中可以看出，在没有加入边界条件时，当波传播到边界时，在边界处产生了很强的反射，这是我们在进行数值模拟计算时所不期望出现的。为了消除或减弱这种边界反射效应，得到地质地层真实的反射信息，就需要对人工边界进行处理，从而得到更接近于实际空间中波的传播规律。在处理边界反射时常用的边界条件方法有高衰减区法、傍轴近似法、吸收边界法和反周期扩展法1,2,25。高衰减区法是在靠近边界处引

22、入一个高吸收区，通过耗散因子使入射波仔这个吸收区域内逐渐衰减，从而抑制边界附近的人工反射。旁轴近似法是将边界区的波动方程用单向的傍轴近似波动方程代替，只模拟波能量由差分网络向外边界的单向传播，从而消除边界反射的问题，不过这个方法只适用于小角度入射的情况。吸收边界法是先导出吸收边界条件方程，然后让方程与计算区域内的波动方程联立求解，从而使从人工边界向计算区域反射回来的波能够全部或部分被吸收。反周期扩展法，即利用正反周期函数极性相反的特点消除回绕波场，这个方法在理论上是能够完美的消除回绕波场而不产生任何的反射，但是在实际应用时，由于这种方法需要进行反周期波场的扩展，大大的增加了计算量，而且为了得到

23、最后叠加的平均波场，正反周期波场都必须存储，这样同时也会占用很大内存空间25,34,35。下面将就本文将要用到的透明边界条件进行讨论。5.1 一维透明边界条件对于一维均匀介质波动方程：式中 V 表示声波的传播速度，U 表示声压。我们引入一个向右行波（k =/V），将这个平面波的方程代入上面的均匀介质波动方程（5-1）中可以解得：其中r表示反射系数。在右边界x=a上反射波与入射波的振幅相等，即|r|=1。当引入一个左行波时，左边界x=a上同样也会产生反射，且|r|=1。由此可以看出，如果边界条件选取不当会使左右两侧的边界上产生强反射，因此就需要在左右边界处选取的边界条件能使r=0。则右行波就满足

24、方程左行波就满足方程由此可以看出，如果在右行波上加上一个左行波项，则只有当r = 0时方程（5-5）才成立，即当选式（5-5）作为右边界条件时，则在右边界上不会产生反射，因此可以将式上式作为右边界条件。同理可得左边界条件为通常用右行波和左行波之和来表示波函数，因此也只有当 r = 0时，波函数才满足式（5-4），所以在引入边界条件以后左边界和右边界上都不会产生反射。5.2 二维透明边界条件设在平面区域为，对于已知某初始条件下的二维波动方程如果给定一个边界条件则平面波在边界x = a、 x = a及 z = b处会产生反射。考虑一个向右传播的平面波其中表示入射角，即平面波前与 x 轴的夹角。那么

25、这时反射波为其中r为边界产生的反射波的反射系数。则总波场可以表示为以x=a处为例，把式（5-11）代入U(a,z,t)=0可得|r|=1，即在右边界上产了反射。由于对于波动方程（5-6）当波只在水平方向传播时，即=0时在右边界x = a，如果则由式（5-12）可得把式（5-11）代入式（5-15）可得r=0，这时在右边界x=a上就不会产生反射。同理可得在左边界x=a处，要使r=0，则取如果取并在边界上令当x=a时取当x=a时取以上这两种边界条件都只有很小的边界反射，但还需要考虑到有限差分法的稳定性条件。如果把项的系数看成 当的特殊情况下，的系数就为1/2。则（5-19）式可以写为由波动方程式（

26、5-6）和上式联立推导可得对式（5-22）进行整理化简就可得，x = a时的右边界条件为同理可得x=-a时的左边界条件为当z=b时，底边界条件为以上得到的式（5-23）、式（5-24）及式（5-25）就分别为二维均匀介质波动方程的右边界条件、左边界条件和底边界条件。5.3 透明边界条件的差分格式在推导出二维均匀介质波动方程的透明边界条件基础上，下面在进一步推导一下这些边界条件的差分格式。对右边界条件（5-24）两边同乘以h可得下式其中h=x=z，s=Vt/x。由于通过泰勒级数展开式有因此把式（5-28）和式（5-27）代入式（5-26）可得上式中的微分项可以表示为然后把式（5-30）各式代入式

27、（5-29）即得到右边界条件的差分格式同理可以得到左边界条件x=a的差分格式为底边界条件z = b的差分格式为顶边界条件z=0的差分格式为上面得到的式（5-34）、式（5-33）、式（5-32）及式（5-31）就是二维均匀介质波动方程透明边界条件的差分格式。如下图5-1（右）所示，可以看出在加入了边界条件后，边界反射明显减少，说明透明边界条件对边界反射有很好的吸收效果。图5-1 未加边界条件（左）与加入边界条件（右）波场快照对比图（所加边界条件为透明边界，差分精度为时间2阶空间2阶）6 简单模型试算6.1均匀模型模型速度为2000m/s，密度为常数，网格大小为500×500，网格空间

28、采样步长为1m，采样时间间隔为100s，震源坐标为（250m，20m），震源震动时间5ms，子波为雷克子波，频率为30Hz，模型见图6-1。图6-2中包含了20-120ms的波场快照，差分精度为时间2阶空间12阶，边界条件为透明边界条件。从图中可以看出，波以震源为中心，呈圆弧向外传播，当波传播至顶边界时有微弱的反射波，说明利用透明边界条件并不能无完全消除人工边界产生的反射波。另外，从图中还可以看出有轻微的频散现象。 图6-1均匀模型示意图图6-2 均匀模型波场快照图（差分精度为时间2阶，空间12阶，边界条件为透明边界条件）6.2 水平层状模型和高速透镜体模型图6-3（左）所示的水平层状模型，模

29、拟网格大小为500×500，网格间距1m；表层速度为2000m/s，厚度为100m，第二层速度为3000m/s，厚度为200m，第三层速度为3500m/s，厚度为200m。采样时间间隔为100s，震源坐标为（250m，20m），持续时间为5ms，子波选用雷克子波，子波频率为30Hz。图6-3（右）所示的高速透镜体模型，模拟网格大小为500×500，网格间距1m。透镜体为一个半长轴为100m，半短轴为50m的椭圆，中心坐标为：z=250m，x=250m；高速透镜体的速度为3500m/s，周围介质速度为2000m/s。采样时间间隔为100s，震源坐标为（250m，20m），持续

30、时间为5ms，子波选用雷克子波，子波频率为30Hz。图6-4中T=40ms时，波在第一层与第二层的界面上开始传播，从图中T=60,80ms中可以清晰的看到第一个界面的反射波和透射波，T=100ms时波还没有传播到第二层与第三层的界面，T=120ms时波传播第三层，从图中可以看到产生了一个反射波。顶边界反射波比较明显，同时，随着波向外传播，频散现象也逐渐明显。图6-5中T=80ms时，波时波传播到透镜体内，在透镜体的顶界面产生了反射波，T=140ms波场快照中，出现了波从透镜体内向外传播时产生的反射波。T=100ms的波场快照中频散现象较为明显。图6-3 层状模型（左）和高速透镜体模型（右）示意

31、图图6-4 水平层状模型波场快照图（差分精度为时间2阶，空间12阶，边界条件为透明边界条件）图6-5 高速透镜体模型波场快照图（差分精度为时间2阶，空间12阶，边界条件为透明边界条件）7结论与建议地震波场的数值模拟是地震勘探学的重要研究内容之一，是进一步研究地震资料的采集、处理和解释的有效辅助手段。论文以波动理论为基础，对有限差分模拟的一些关键性问题进行了探讨和研究，主要包括波动方程有限差分格式的建立、边界条件、初始条件、震源问题、有限差分的稳定性及频散问题等。经过本次课程报告，我对有限差分波场模拟有以下认识：（1）有限差分算法做波动方程数值模拟时，具有计算速度快、算法简单等特征。（2）有限差

32、分波动方程是利用泰勒级数展开式展开推导得到的，差分的阶数越高，有限差分的误差就越小，而有限差分的解就越精确。（3）有限差分模拟时，网格空间步长的选择不仅对模拟的精确度有影响，而且对边界的处理效果也有着明显的影响。同时由于震源项是以离散形式加入的，所以网格步长也必须满足采样定理，这样才能取到完整的子波，以保证计算结果不会出现较大的误差。因此，有限差分数值模拟计算时，必须合理地选择网格步长。（4）数值模拟过程是在特定的区域内进行，因此就必然会产生人工边界，波传播到边界就产生了反射，通过对几种模型的模拟表明，本文所用的透明边界条件和吸收边界条件对边界反射均有较好的吸收效果。本课程报告的工作在主要是在

33、前人的研究基础上开展的，但是还存在一些问题需要进一步的完善和改进，建议主要从以下几个方面进行：（1）网格剖分是影响模拟结果的主要因素之一，本文采用的仍然是常规的规则网格剖分，规则网格剖分方法在模拟计算中仍然存在一些缺陷，当模型较复杂或要求的精度较高时，往往不能达到理想的模拟效果。如果减小网格间距，增加网格数量，又会导致计算速度和计算效率降低，因此为了获得较满意的模拟结果，需要采用更优的网格剖分方法，以克服存在的缺陷。（2）对于模拟中出现的数值频散现象，还需要更深入的研究，以便能够尽可能的减小数值频散。（3）本报告的算法主要针对简单地质模型，如果对于复杂的介质模型，算法还需要进一步的研究和完善。

34、参考文献1裴正林，牟永光.地震波传播数值模拟.地球物理学进展J，2004，19(4)：933941.2马在田.计算地球物理学概论M.上海：同济大学出版社，1997.3何樵登.地震波理论M .地质出版社，19884何樵登，熊维钢.应用地球物理教程地震勘探M .地质出版社，19915吴清岭，张平，施泽龙.波动方程正演模拟即应用J.大庆石油地质与开发，1998，17（3）6江玉乐，雷苑著.地球物理数据处理教程M.北京：地质大学出版社，2006，77姚姚著.地震波场与地震勘探.北京：地质出版社，2006，6.8张永刚.地震波数值模拟方法J.石油物探，2003，42（2）：143-147.9陈伟.起伏地

35、表条件下二维地震波场的数值模拟J.勘探地球物理进展，2005，28（1）10吴永刚，吴清岭.有限差分法弹性波场数值模拟J.大庆石油地质与开发，1994，13（3）11Alterman. Propagation of elastic wave in layered media by finite difference methods. Bulletin of the Seismological Society of America,1968;58(1):36739812Alford, R. M. Kelly ,K .R. and Booer, D. M. Accuracy of finite d

36、ifference modeling of the acoustic wave equation.Geophysics,1974;39(6):83484213周家纪，贺振华.模拟地震波传播的大网格快速差分算法J.地球物理学报，1994；37（增刊）：45045414董良国.一阶弹性波方程交错网格高阶差分解法稳定性研究J.地球物理学报，2000，43（6）：85615邵秀民，蓝志凌非均匀各向同性介质中地震波传播的数值模拟J.地球物理学报，1935；38(1)：395516刘洋，李承楚，牟永光.任意偶数阶精度有限差分数值模拟J.石油地球物理勘探，1998，53（1）17褚春雷，王修田.非规则三角网

37、格有限差分法地震正演模拟J.中国海洋大学学报，2005，35（1）：04304818董良国，李培明地震波传播数值模拟中的频散问题J.天然气工业，2004，24(6)：53-5619奚先，姚姚.二维随机介质及波动方程正演模拟J.石油地球物理勘探，2001，36（5）：54655220董良国等.一阶弹性波方程交错网格高阶差分解法J.地球物理学报，2000，43（3）:41141921奚先，姚姚.二维弹性随机介质中的波场特征J.石油地球物理勘探，2004，39（6）：67968522邵秀民，蓝志凌非均匀各向同性介质中地震波传播的数值模拟J.地球物理学报，1935；38(1)：395523孙若昧.地震

38、波传播有限差分模拟的人工边界问题J.地球物理学进展，1996，11（3）：5324 Cerjan C, Kosloff D, Kosloff R and ReshefM. A nonreflecting boundary condit -ion for discrete acoustic and elastic wave equation. Geophysics, 1985, 50 (4) : 70570825Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equation.B

39、ulletin of the Seismological Society of America, 1977;66(11):1529154026Reynolds. Boundary conditions for the numerical solution of wave propagation problem.Geophysics, 1978;43(6):1099111027崔兴福，张关泉.地震波方程人工边界的两种处理方法J.石油物探，2003，42（4）：45245528张毅.声波正演模拟中的人工边界问题J.工程地球物理学报，2007，4（4）29李文杰，魏修成，刘洋.声波正演中一种新的边界

40、条件-双重吸收边界条件J.石油物探，2004，43（6）：52853130邢丽.地震波数值模拟中的吸收边界条件J.上海第二工业大学学报，2006，23（4）31Dahlain M A. The application of high difference to the scalar waveequation. Geophysics, 1986;51(1): 546632 Higdon R L. Absorbing boundary condition for elastic waves. Geophysics,1991, 56 (2) :23124133孙成禹，张吉辉完全纵波方程有限差分数值模

41、拟J.石油地球物理勘探，2005；40(3)：28929434黄自萍，徐伶俐，周继顺.起伏地表声波方程的数值模拟J.石油地球物理勘探，2006，41（3）：27528035宛书金，董敏煌等.横向各向同性介质中的地震波传播特性J.石油大学学报，2002,26(1)：2328附录本课程报告所编matlab程序% 2015.5.20 纵波波场模拟%clcclear% 初始化mode=12; %选择差分介数，最高支持16介Model=load('Model.mat');V_model=Model.UniformModel;%速度模型c=load('c.mat');c=c

42、.c; %系数矩阵cellnum=size(V_model); %网格数目dt=100*10(-6); %时间采样间隔 100usdx=1; dz=1; %空间采样间隔x0=250; z0=20; %震源坐标f0=30; %震源子波频率S_t=2*10(-3); %震源持续时间 2ms T=140*10(-3); %波场快照时间 %U0=zeros(cellnum(1),cellnum(2); %波场初始值U1=zeros(cellnum(1),cellnum(2);U2=zeros(cellnum(1),cellnum(2);for k=1:T/dt; if k<S_t/dt U1(z

43、0+1,x0+1)=sin(2*pi*f0*k*dt)/(k*dt);%震源函数。 end % if mode=2 for m=mode/2+1:cellnum(2)-mode/2 % 2介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end % if mode=4 for m=mode/2+1:cellnum(2)-

44、mode/2 % 4介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end % % if mode=6 for m=mode/2+1:cellnum(2)-mode/2 % 6介差分方

45、程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. c(mode/2,3)*a2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end

46、% if mode=8 for m=mode/2+1:cellnum(2)-mode/2 %8介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. c(mode/2,3)*a2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1

47、(m,n+3)-4*U1(m,n)+. c(mode/2,4)*a2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n)+. 2*U1(m,n)-U0(m,n); end end end % if mode=10 for m=mode/2+1:cellnum(2)-mode/2 % 10介差分方程 for n=mode/2+1:cellnum(1)-mode/2 a=V_model(m,n)*dt/dx; U2(m,n)=c(mode/2,1)*a2*(U1(m+1,n)+U1(m-1,n)+U1(m,n-1)+U1(m,n+1)-4*U1(m,n)+. c(mode/2,2)*a2*(U1(m+2,n)+U1(m-2,n)+U1(m,n-2)+U1(m,n+2)-4*U1(m,n)+. c(mode/2,3)*a2*(U1(m+3,n)+U1(m-3,n)+U1(m,n-3)+U1(m,n+3)-4*U1(m,n)+. c(mode/2,4)*a2*(U1(m+4,n)+U1(m-4,n)+U1(m,n-4)+U1(m,n+4)-4*U1(m,n)+. c(mode/2,5)*a2*(U1(m+5,n)+U1(m-5,n)+U1(