弹性波课程报告 有限差分纵波波场模拟姓名： 邓 健 学号： 1201410215 课程名称：弹性波理论 任课老师：顾 汉 明 专业：地球探测与信息技术 时间： 2015年5月 评 语对课程论文的评语:平时成绩：课程论文成绩：总 成 绩：评阅人签名：注：1、无评阅人签名成绩无效；2、必须用钢笔或圆珠笔批阅，用铅笔阅卷无效；3、如有平时成绩，必须在上面评分表中标出，并计算入总成绩。摘要地震波场数值模拟是勘探地震学的重要研究课题之一,也是研究波动现象的重要手段。地震波场数值模拟常用的方法有伪谱法、有限差分法、有限元法等。有限差分法具有计算速度快、占用内存小等优点，该方法对于近远场及复杂边界都有广泛的适用性，能够准确地模拟波在各种介质及复杂结构地层中的传播规律，是勘探地震学中应用最广泛的数值计算方法。本文以波动理论为理论基础，利用泰勒级数展开式推导出波动方程的有限差分格式及其离散表达式。针对传统二阶精度差分方法模拟精度较低的问题，推导出时间域二阶空间域四阶精度的有限差分方程，并综合分析了初始条件、震源项、算法稳定性及数值频散等因素在有限差分法数值模拟中的影响。有限差分数值模拟计算是在一个有限的区域内进行的，当波传播到边界时就会产生边界反射，这是进行数值模拟计算时所不期望出现的。本报告通过调研大量文献资料，对效果较好的透明边界条件进行了详细分析。在此基础上,分别推导得到这种边界条件方法的离散表达式，并在数值模拟过程中进行验证。通过建立均匀模型、层状模型及高速透镜体模型，基于Matlab语言编程，论文实现了波动方程有限差分算法的数值模拟。对不同地质理论模型进行模拟，从模拟结果的对比分析中可以看出，模型参数的选择及各参数间的相互关系对模拟结果有着显著影响。无论是模拟精度，还是模拟计算效率，有限差分算法都具有一定优势。通过综合研究，论文认为波动方程有限差分算法具有算法简单、计算效率高、模拟精度较高等特点，理论研究意义大，应用前景广阔。关键词：数值模拟，波动方程，有限差分，边界条件目录1 前言11.1 地震波场数值模拟概述11.2 有限差分法及其与几种常用数值模拟方法的比较12 波动方程有限差分法数值模拟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 地震波场数值模拟概述地震波场数值模拟简单说来，就是已知地下介质构造及其参数，再用理论计算方法来研究地震波在地下介质中的传播规律，合成地震记录的一种技术方法。随着地震勘探技术的发展，数值模拟方法已经贯穿于地震数据的采集、处理和解释的全过程，而且在确定观测的合理性、检验处理和解释的正确性等方面都有了广泛的应用1-3。地震波数值模拟问题的研究主要包括三个方面：(1)地震波场数值模拟原理；(2)地震波场数值模拟算法；(3)数值模拟的计算机实现。地震波场数值模拟方法的理论基础主要是波动理论和射线理论。射线理论方法是求出质点运动方程的渐进近似解，这种方法计算速度快，但是精度较低。而波动理论方法则是用数值计算方法直接求出波动方程的解，其模拟结果中包含丰富的波动信息，模拟结果较为精确。数值模拟又可分为波动方程的解析解和数值解。对简单的地质模型就可以得到其精确解，但对于复杂的地质模型，只有通过数值解来说明波在底下介质中的传播。通常解析解只是作为数值解的一种检验4-6。1.2 有限差分法及其与几种常用数值模拟方法的比较在地震勘探中，为了研究地震波在地下各种介质中的传播规律，就需要对波场进行数值模拟。常用的数值模拟方法有伪谱法、有限差分法，有限元法等。有限元法依据变分原理，通过灵活的网格剖分，用简单形态逼近实际的地质体，能处理多种介质和自然边界条件。适用于非规则网格的计算，方便有效，模拟的精度高。但是有限元法的问题是不适用于大规模的模型的计算，而且计算量大7。伪谱法在20世纪七十年代引入数值模拟计算领域的，它是利用傅立叶变换来计算波场的空间导数，用差分方法来计算时间的导数，可以看成是一种无限阶的有限差分法，是传统有限差分法的一个推广。伪谱法在粗网格上也能实现高精度的计算，相对有限元法实现起来较容易，在非线性波动问题及气象预测等领域有着广泛的应用8-10。有限差分法是一种最常用的数值模拟方法，它是将波动方程中波场函数的空间导数和时间导数用相应的空间、时间的差分代替。有限差分法具有计算速度快、占用内存小等优点，该方法对于近远场及复杂边界都有广泛的适用性，能够准确地模拟波在各种介质及复杂结构地层中的传播规律。有限差分法方法简单、高效的优点是其他方法难以比拟的，因此有限差分法目前仍然是勘探地震学中应用最广泛的数值计算方法11。2 波动方程有限差分法数值模拟2.1 有限差分原理有限差分法就是把求解区域划分为差分网格，然后用有限的网格节点代替连续的求解域，利用微商与差商的近似关系将描述介质传播的微分方程转化为差分方程进行求解。差分离散的方法有两种，一种是将单变量的二阶波动方程直接转化为时间空间的二阶中心差分进行离散求解；第二种方法是把用位移表示的二阶波动方程转化为用应力及质点速度表示的一阶方程组，然后用应力和速度的交错网格求解,8。构造差分的方法也有很多种，一般采用泰勒级数展开法。常见的差分格式有显示差分、隐式差分和显隐交替格式；按照差分的精度又可以分为一阶差分、二阶差分和高阶差分等，目前通常见到的差分格式主要是几种差分格式的组合。我们知道，差分方程的建立首先选择网格布局和差分形式，然后以有限差分代替无限微分，以差商代替微商，并以差分方程代替微分方程及其边界条件，最后建立差分方程12。在建立差分方程的时候要注意到两个方面。一是合理选择网格布局与步长。我们将离散后各相邻离散点之间的距离或者离散化单元的长度称为步长。在所选定的区域内进行网格划分是差分方程建立的第一步，其方法比较灵活，但是实际应用中往往遵循误差最小原则。常用的典型的网格剖分方法有矩形网格、三角形网格等，网格样式的选择和区域的形状有关。其次是将微分方程转化为差分方程。这个过程就是选择一种差分格式代替微分形式的过程。构造差分的方法有多种形式，本文采用的是泰勒级数展开法2,5,13。地震波场数值模拟以地震波动理论为基础，用有限差分法解波动方程时，对变量离散化，也即对连续的物理量只考虑其在离散的空间位置和离散时刻的值，然后把方程中的导数用这些离散的采样值表示14。对于一个单变量的函数f(x)，将其离散化，那么在采样点x =lx的采样值就是f(lx),其中x表示步长，l为整数。则有限差分法中f(x)在采样点x=lx的导数就可以近似表示为： 其中an是系数，N表示差分格式的长度，差分的格式是由这两个系数来决定。在实际应用中常用的差分格式有向前差分、向后差分以及中心差分：（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）就是二维非均匀介质弹性波方程。而在均匀介质中，、和均为常数，则二维均匀介质弹性波方程就为：把（2-12）中的x和z分量合并就得到了二维均匀介质弹性波方程的矢量形式表达式在没有弹性横波只有弹性纵波存在时，对（2-13）两边取散度：利用旋度与散度的关系，交换上式的微分次序并化简可得：其中表示纵波速度。根据赫姆霍兹在他著名的有关涡流运动的著作中证明了下属定理：任何向量点函数，若它的散度和旋度具有位，则它可以表示为一个无旋部分和一个旋转部分之和。亦即对于任何一种向量场，如果在其定义域内有散度和旋度，则该向量场可表示为一个标量位的梯度场和一个向量位的旋度场之和，而空间传播的波动正是无旋运动和旋转运动这两种运动之和15-20 令代入式（2-8）中，取其中的无旋部分则有：然后再消去式中的梯度,就得到了用位移位表示的纵波方程：3 二维声波方程有限差分格式的建立3.1 声波方程一般地，二维非均匀介质的声波波动方程可表示为：式中 U = U ( x , z , t)表示声压，V表示声波在介质中的传播速度，是密度，它是随空间各点而变的。其中 s ( x , z , t )为震源函数。对于均匀介质，密度为常数，则二维声波波动方程就可以表示为：在实际问题中，我们总是用一个有限宽频带的时间函数来代替 s ( x , z , t )函数，以便能够真实地反映地震波的传播21。3.2 波动方程有限差分格式的建立上一节中讨论了差分的原理和几种常见的差分格式，在这个基础上，我们结合波动方程，对纵波波动方程的有限差分格式进行推导。令，其中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的中心差分方程就可以得到二维波动方程的有限差分方程：对上式继续推导可得：其中，上式就是二阶波动方程的有限差分格式。如果我们继续利用泰勒级数展开式，则可以得到：这样我们就能够得到时间域二阶空间域四阶的波动方程有限差分格式：上式中k表示时间方向的离散网格，m表示x方向的离散网格，n表示z方向的离散网格，。需要注意的是，这里同样也假设x=z=h，即网格步长相等5 ,22-24。3.3 高阶差分方程将采用高阶差分会有很多时间层，计算较复杂。为此时间上采用2阶，空间上可采用高阶。为了方便起见，本文空间步长在x，z轴上相等，则通过上述类似推倒，可以得到不同精度的声波差分方程。时间2阶，空间6阶时间2阶，空间8阶时间2阶，空间10阶时间2阶，空间12阶时间2