毕设论文--粘声波正演模拟研究

《毕设论文--粘声波正演模拟研究》由会员分享，可在线阅读，更多相关《毕设论文--粘声波正演模拟研究（36页珍藏版）》请在金锄头文库上搜索。

1、-本本题科科 毕毕 业业 设设 计（论文）计（论文）学生姓名：xxx学号：xxx专业班级：xxx指导教师：xxx2015 年 6 月 20 日-目：粘声波正演模拟方法研究-粘声波正演模拟方法研究摘要地球上介质的黏滞性会引起大地的吸收效应 ,它会影响波场所有的频率成分 ,尤其对于高频的影响最大,导致地震分辨率降低。黏滞吸收作用会影响地震波波形、频带、振幅等因素。一个高效的粘声波正演模拟方法，可以考虑到由于实际介质造成的地震波的吸收衰减作用。可以更加准确模拟地震波在非完全弹性实际地层中的传播，在这里，本文通过编程建立不同的粘声波方程数值模拟模型跟正常的声波方程数值模拟模型进行对比分析，从而了解粘声

2、波正演模拟方法的优越性。关键词：粘声波；正演模拟；有限差分；-Study on the forward modeling of viscoelastic acousticStudy on the forward modeling of viscoelastic acousticwaveswavesAbstractAbstractThe absorption effect is mainly caused by the viscosity of the earth mediaitself.The viscous stagnation can affect all the frequency co

3、mponents of the wavefield.And the effect of the high frequency components is bigger,which leads to thedecrease of seismic resolution.The absorption of the absorption has a great influence onthe wave, frequency and amplitude of the seismic wave. A highly effective viscoelasticforward modeling method

4、can take into account the absorption and attenuation ofseismic waves by real media. Accurate simulation of the propagation of seismic wavesin the actual strata of the imperfect elasticity. Here. In this paper, the program, establishdifferent visco acoustic wave equation numerical simulation model wi

5、th normal acousticwave equation numerical simulation model for comparative analysis, to understand thevisco acoustic forward modeling method of superiority.KeywordsKeywords：Viscoelastic acoustic wave；Viscoelastic acoustic wave；Finitedifference;-目录第 1 章引言.1 1.1 研究意义.1 1.2 正演模拟方法.1 1.3 国内外研究现状.3 1.4 本

6、论文研究内容.3第 2 章粘性介质基本理论.52.1 粘性介质的基本特点.52.2 粘性介质模型的构建.62.2.1 开尔芬固体模型.62.2.2 标准线性体.72.3 品质因子.7 2.3.1 定义.8 2.3.2 与吸收系数的关系.8 2.3.3 与介质速度 v 的关系.8 2.3.4 与震源频率的关系.8第 3 章二维各向同性介质粘声波方程数值模拟.103.1 粘声波方程及其交错网格高阶差分格式.10 3.1.1 粘声波方程的推导.10 3.1.2速度-应力方程的推导.11 3.1.3交错网格有限差分.123.2 模拟震源.16 3.2.1 震源的选择.16 3.2.2 震源的类型.17

7、3.3 边界条件.183.3.1 衰减边界条件.193.3.2 PML 边界条件 .19 3.4 稳定性.20- 3.5 数值频散问题.22 3.6 模型试算.23 3.6.1 地震波的传播规律研究.23 3.6.2粘声波的衰减规律研究.27第 4 章结论.28致谢.30参考文献.31-第 1 章引言1.1 研究意义我国的石油天然气等资源类工业发展的非常快，随着油气资源的不断发现并且开采，现阶段所存在的油气资源越来越少，因此油气勘探工作变的越来越困难。勘探的目标更加复杂。而现阶段主要的勘探手段地震勘探，它效果的好坏取决于理论模型与实际地层的符合程度。在地震勘探中，我们通常把地下介质为看成弹性各

8、向同性介质。但是在实际中，地下介质存在各向异性和吸收衰减现象，并着这种现象非常的普遍。它会产生大地吸收效应1。同时，频率成分（特别是高频成分）会受到粘滞性影响，致使地震记录、地震成像的垂向分辨率降低。造成理论上的研究结果跟实际中观测到的记录在频率、振幅等方面存在很大的不同。因此，以前的那种将地下介质当作均匀各向同性完全弹性的理论是与实际情况差别很大的。实际的介质更趋向于粘弹性，对粘弹性介质进行理论研究，更加的符合实际情况。充分的顾及到了实际介质存在的吸收效应跟衰减效应，更加准确的还原地下介质中的实际情况2。综上所述， 粘弹介质正演模拟方法的研究对于还原地下介质的真实情况具有较强的理论意义和实际

9、意义。但是，粘弹性介质的正演模拟计算成本较高，计算繁琐，所以。为了降低计算成本，本文进行了声学近似下的粘声波正演模拟方法研究。1.2 正演模拟方法正演方法主要分为：数值模拟、物理模拟。数值模拟是最常用的。它通过给出的结构模型和物理参数，模拟地震波的传播轨迹，了解其规律以及过程，然后通过计算来推断观测点的地震记录。在勘探过程中，数值模拟的作用很大。例如： 1、采集上，可用于设计或者优化野外观测系统；2、处理上，可以通过数值模拟来检验是否采用了正确的反演方法。将正演反演不断的逼近，从而使结果更加准确； 3、解释上，还可以检测一下解释的资料是否正确。地震波正演模拟方法基于理论基础和表达形式可以分为三

10、种：1积分方程法(基于惠更斯原理）；2地震波方程数值解法（基于波动方程）；3射线追踪法（基于射线理论）。见图 1-1。-根据求解方法，地震波方程数值解法可分为有限元法、伪谱法、有限差分法。具体情况见表 1-2.由于有限差分法具有内存占用较小，精度较高等优点，因此本文主要采用这种方法进行模拟6。-1.3国内外研究现状地层的吸收作地震波在地下介质中传播时会不断衰减。为了更加准确的模拟地震波在地层中的传播，国内外相关专家学者做了大量的研究工作。在这里，我主要介绍一下国内外关于粘弹或粘声介质正演模拟方法（主要是有限差分法）的一些发展历程。在国外，1845 年，Stocks 在研究过程中，最早提出粘弹性

11、理论，认为地层中非弹性机械能转化成热能的能量损失是因为地层的内摩擦。到了 20 世纪，N.H.Ricker 于四十年代提出粘弹介质更加接近实际地层，推动了粘弹理论的发展。在 1968 年，为了研究研究层状介质弹性波数值模拟，卡拉尔和奥尔特曼将有限差分法应用其中。1974 年，Alford 又验证了声波方程有限差分法模拟的精确性。在1988 年，Carcione 提出了一种方法，用于模拟粘滞声波在地层中传播。到1990 年，Crase E 将高阶差分与交错网格相结合。这个也是这篇论文所研究的一部分 在 1995 年 RKnight 等人提出了在二维随机介质中粘声波的传播以及吸收和波散衰减。199

12、5 年，How-Wei Chen 等人提出了一种并行算法，用来进行粘滞声波数值模拟，并在二维介质上得到了验证。1996 年，Jose M.Carcione 等人进行了有关流体饱和粘声波孔隙介质中纵波的全频率瞬时解的研究。 今年来相关的研究有所减少。在国内，粘弹的研究很多，但是关于粘声的研究相对少一些，起步比国外晚。2010年，宫同举在孙成禹导师的指导下，进行了近地表地震波衰减规律及品质因子提取方法研究。2011 年，廖建平等人进行了二维频率空间域粘声波正演模拟研究，采用了嵌套剖分网格排序法，以便于节省计算机内存。通过数值试验得到结论：频率空间域对于模拟吸收与衰减效应非常简单，更加有利于后面的保

13、振幅处理工作。2013 年，孙天真等人进行了基于粘声介质的反滤波叠前深度偏移方法研究。2014 年，郭成锋、杜启振进行了有关二维 VTI 介质粘声波正演模拟的研究。同年，吴玉等人进行了基于分数阶拉普拉斯算子解耦的粘声介质地震正演模拟与逆时偏移的研究， 逼近时间导数利用的是交错网格有限差分，解分数阶拉普拉斯算子利用的是改进的伪谱法。李金丽等人进行了粘介质交错网格有限差分数值模拟研究。1.4 本论文研究内容本文主要在前人研究成果的基础上， 来研究粘声波的正演模拟。首先推导了非均匀各向同性介质情况下粘声波的一阶应力-速度方程， 然后推导构建二维粘声波方程正演模拟的交错网格有限差分格式，之后利用编程，

14、获得粘声波正演模拟的结果，最后，建-立不同的模型， 对粘声波方程数值模拟模型跟正常的声波方程数值模拟模型进行对比分析，从而了解粘声波正演模拟方法的优越性。按章节来讲：第一章为引言，主要介绍粘声波正演模拟研究的意义，以及国内外的主要研究现状。并且简单的介绍了正演模拟方法的基础内容。第二章主要介绍粘性介质的基础内容。包含其基本特点、构建模型的方法等。并且提到品质因子的概念。第三章主要是推导一阶应力-速度方程，构建交错网格有限差分格式，进行模拟。并对粘声波进行模拟得到的模型与声波模型进行对比，分析其特点。第四章主要分析得到的结论。第五章主要介绍缺点以及展望。-第第 2 2 章章粘性介质基本理论粘性介

15、质基本理论实际地层并不是完全弹性各向同性的，而是粘性介质，地震波在地层中传播时，能量逐渐衰减，不断被吸收。粘性介质理论及实践，在实际生产和生活中发挥着巨大的作用。在这一章，主要介绍一些关于粘性介质的基本概念、特点、模型的构建、品质因子简介等等。2.1 粘性介质的基本特点完全弹性介质没有考虑到压力、温度等因素引起的波形的改变。完全弹性介质理论认为，应力应变关系满足虎克定律，即应力应变成正比关系。完全弹性体因外力引起的应变是可逆的，当外力消失的时候，应变完全恢复。也就是说弹性体不会脱离一种准平衡的状态14。其内部每一处的温度都是完全一样的。认为热运动跟机械运动之间不会产生能量的转换，因此地震波在其

16、中传播时没有衰减。粘性流体介质当存在外力时，它随着时间的进行，发生的形变是永久的，当外力消失不会完全恢复原状，这体现的是流体性质的完全非弹性。粘弹性介质的性质表现为介于粘性流体与弹性固体之间的所有性质， 实际的介质相比于完全弹性介质与粘性流体，更加近似粘性介质13。粘性介质的行为随时间而变化。有以下几个特点：（1）蠕变在恒定应力作用下，粘性介质会发生连续的缓慢的应变，叫做蠕变。能力随蠕变的增大而增强。（2）应力松弛应力松弛，指的是粘性介质的应变持续不变的时候，应力会逐渐减弱的现象。可以简单的称之“松弛”。（3）松弛模量、蠕变柔量“蠕变模量”为时间域常数应力产生的应变，而“松弛模量”为对应的应变

17、产生的应力。（4）应变率敏感“应变敏感率”指的是，像“杨氏模量”等某些与介质的应变速率相关的介质参数。（5）恢复所谓的“恢复”指的是当其承受的外力消失的时候，粘性介质的应变也会慢慢的消失的现象。-（6）迟滞在加载或去除外力的过程中，它的应变响应相对于应力要滞后，导致介质的应变和应力形成“迟滞回线”，它下面的面积表示能量损失的多少。（7）记忆变量粘性介质中，应力相当于是不同阶段的应力相互叠加的结果。应力与应变的卷积的关系让在时间域的相关的处理变得非常繁琐，故引入概念“记忆变量”，用于避免这种卷积关系6。2.2粘性介质模型的构建在实际地层（非完全弹性）中，由于介质内摩擦力的作用，应变、应力不在成正

18、比关系。关于粘性介质中的应力与应变之间的关系，相关专家学着进行了很多探索，提出了很多不同的粘性介质模型的建立方法。其中，主要是两种方法：累计积分法跟器件组合法。累计积分法的思想是指：总应力 （应变） 是各阶段应变（应力） 所引起的的应力（应变）的线性叠加。波尔兹曼粘弹性介质以波尔兹曼叠加原理为基础，它在累计积分法的模型中非常典型。器件组合法是指：性质不同的单元体通过混联、并联、串联的方法组合起来。其中，有两种最基本的方法，即： Kelvin 模型、Maxwell 模型。由于本文在进行模型的构建时候，用的是器件组合法之中的Kelvin固体模型，所以，此处主要介绍一下关于开尔芬模型的相关内容。2.

19、2.1开尔芬固体模型开尔芬固体模型是由一个弹性的原件跟一个粘性的原件通过并联组合在一起， 如图2-1：图2-1 Kelvin（开尔芬）固体模型由单元体之间的应变应力跟整个组合的总的应变应力的关系（,其中，为弹性原件的应力，为粘性原件的应力)，我们能够得到Kelvin（开尔芬）固-体模型的本构方程： = =q0 + +q1在时间域上， 均匀弹性各向同性介质中的物性参数跟分别对应开尔芬粘固体模型中的+lt和+lt。其中 ，各种参数的含义为：、（ 与 ）表示均匀弹性各向同性介质（粘性介质）中的拉梅系数（粘性拉梅系数）；lt为时间偏导数算符；表示的是材料的松弛时间（ F /）7。2.2.2标准线性体由

20、于应力下的应变产生突变或者突然失去应力的情况开尔芬模型考虑不到， 因此对于大多数粘性介质它描述得不够全面。很多科学家基于这一点，提出了一种“标准线性体”的思想。可以很好的改善这一点，模型如 2-2 所示。图 2-2 标准线性粘弹性介质本构方程形式为：其中，松弛弹性模量为 MR;松弛时间为跟。在时间域中， 均匀弹性各向同性介质与标准线性粘弹性固体模型的物性参数的对应原则为：121lt12ltlt-2.3 品质因子在研究粘性介质的时候，人们为了补偿地震波在传播过程中造成的能量的衰减，引入了一个物理量：品质因子 Q。品质因子在粘性介质研究过程中扮演着非常重要的角色3。2.3.1 定义品质因子的表达式

21、为：1EQ2E其中，假设E为一个周期内的最大能量，E为一个周期内消耗的能量。2.3.2 与吸收系数的关系Q 与的关系为：Q 从关系式可以看出，两者成反比例关系。当品质因子增大时，吸收系数减小，也就是介质的吸收作用减小，能量消耗减小，介质更加趋向于完全弹性介质，当品质因子无穷大时，则吸收系数趋向于零，介质成为完全弹性介质。同时，相比于吸收系数，品质因子受频率的影响接近为零，这也是人们在研究过程中经常使用品质因子而不是吸收系数的原因。2.3.3 与介质速度 v 的关系根据经验，Q 与 V 的关系为Q aVb院士李庆忠根据实际观测的资料计算 a、b，给出了一个Q与V的经验关系式：Q 14V2.2(2

22、-15)其中，速度V的单位是km s，Q 为无量纲。2.3.4 与震源频率的关系对于开尔文粘弹性的介质：-其中，Qp 为纵波的品质因子；Qs 为横拨的品质因子；Vp 为纵波速度；Vs 为横波速度； 与 为粘性拉梅系数；为圆频率。通过上式可以看出，两者成反比例关系，也就是说，品质因子随着震源频率的增大而减小，震源频率越大，衰减越强。-第第 3 3 章章二维各向同性介质粘声波方程数值模拟二维各向同性介质粘声波方程数值模拟在这一章，我主要介绍一下二维粘声波方程正演模拟的交错网格有限差分格式构建过程，并进行相关的模拟与模型试算。3.1 粘声波方程及其交错网格高阶差分格式3.1.1 粘声波方程的推导对于

23、完全弹性的物体，它的动量方程是：i,j;i 表示应力i,j（i,j=1，2，3）对 i 求导。为密度，i 表示 ui 对时间的二阶导数，Fi 表示体力。本构关系：体积相对形变：；拉梅系数：，；应变张量：ei,j;单位张量：i,j；并且：当形变很小的时候，实际介质近似为Voigt 黏弹性体。可以得到Voight 体的本构方程：结合方程（1）和方程（5），不考虑纵横波分离，体力为0 的情况下，可得到声波方程：P 表示压力或速度位,如果把地层近似为各向同性粘滞介质，则：-为圆频率，Qp 为纵波品质因子，令：则：3.1.2速度-应力方程的推导根据牛顿第二定律可得vx1 p (12-1)txvz1 p

24、(12- 2)tz将 12 式代入 11 式可得：vx2 vxvz2vx2vz pttx2t2ttt13交换积分顺序并对两边按时间进行积分得：vz2vx2vzu2vx v()()txzxtzt14此处，为了方便美观，将替换为v2p，v 替换为vp，可得最终的应力-速度方程为：2uvz2vz2vx2vx vp() vp()xzxtzttvx1 u xtvz1 u zt(15)在上式中，各字母所代表的含义为：法线应力u；质点速度vz，vx；密度；纵波速度 v=vp，=vp2vp2Qp，其中，为圆频率，Qp为品质因子。-3.1.33.1.3 交错网格有限差分交错网格有限差分(1)交错网格简介网格是用

25、来实现介质模型离散化的一种常用手段。它的分类如下（图 3-1）：三角形网格等贴体网格等网格结构网格非结构网格规则网格常规网格曲线网格交错网格图 3-1 网格分类用高阶差分方法模拟波的传播，一般用两种方法：1、利用常规网格差分方程；2、利用交错网格高阶差分方程。交错网格相比于常规网格存在很多优点。例如，计算量相同时，交错网格数值模比常规网格更精确，数值计算更稳定，收敛速度更快，压制数值频散更彻底。固再此我使用交错网格进行模拟。图 3-2 是二维交错网格的应力、速度分布。图 3-2 交错网格（二维）的应力-速度空间展布(2) 计算差分系数5如图 3-3 所示，在 x 方向上，关于x0对称分布的 2

26、N 个网格节点的坐标分别为x0qNxxxxxx，x0qN1，x0q1，x0q1，x0qN1，x0qN。222222其中，x表示节点间的最小间距；qi表示任意正整数。2N 个网格节点所对应的函数值-图 3-3 以x0为中心对称的 2N 个网格节点分布图已知，分别为fx0qNx/2，fx0qN1x/2，fx0q1x/2，fx0q1x/2，fx0qN1x/2，fx0qNx/2。 利用 Taylor 级数展开求解fx在点x0处的一阶导数近似值。fx0qix/2 fx0qix/2f1x01qix/22f2x021qix/22Nf2Nx0Oqix/22N12N!fx0qix/2 fx0qix/2f1x01

27、qix/22f2x021qix/22Nf2Nx0Oqix/22N12N!32其中，i=1,2，N将上述两式相减，省略式中的误差项，得到1x3fx0qix/2 fx0qix/2 qif1x0qi fx0x3!2qi x2N 1!22N12N2 (16）f2N1x0将相减后得到的式子整理成矩阵形式，有q1q2qNqq3132qN3q xf1x0f x0122q11x32N1q x fx0q21f x023!22x2N212N1 x2N1qNq xx0f f x0N2N 1!22 2N1q x f x012q x (17)f x022q xf x0N2为了简化矩阵，可以记作-q1qA 2qNq13q

28、23qN3q1xf x 02N12q12N1q xq2f x02，D 12x2N1qNq x fx0N2 q x f x012q xf x022qNxf x02同时，构造两个简单矩阵，辅助计算1101，E I 10N1NN假设存在A1， 使得AA1 I， 也可得A1AT I； 即A1为AT的逆， 得到ATA1式子两边右乘向量E就可得 TT T T I。ATA1E E由式（17）可得 (18)f1x0 ETA1D同时，假设(19)AE C1TTc1,c2,cNT（20）将C c1,c2,cN带入式（19），得1Nfx0cnfx0qnx/2 fx0qnx/2xn11（21）可结合式（18）和式（2

29、0），可得到矩阵计算式： q1q312N1q1q2q2q232N1qNc113 qNc20 2N1qNcN0（22）当qi的值确定后，可根据式（22）来求解cn的值，从而计算出f1x0的值。利用式（22）可以求得对称任意节点间距的一阶导数差分系数。其中，当qi取值为2n1(n 1,2,N),则式（22）可表示为-3 11333 2N12N1312N 1c11c02N 132 2N12N 1cN0（23）此时，所求得的cnn 1,2,N就是等节点间距的一阶导数的交错网格不同差分精度的差分系数（表 3-3 所示）。2L24681012a1a2a3a4-6.97545E-4-1.76566E-3-2

30、.96729E-3a5a61.000001.12500-4.16667E-21.17187-6.51042E-24.68750E-31.19629-7.97526E-29.57031E-31.21124-8.97217E-21.38428E-21.22134-9.69315E-21.74477E-21.18680E-43.59005E-4-2.18478E-5表 3-3一阶导数交错网格各阶精度的差分系数(3) 构建二维粘声波方程正演模拟的交错网格有限差分格式设Pi,kj1/ 2为应力u的离散值，Uik1/2, j是速度vx的离散值，Vi,kj1/ 2速度vz的离散值，U1ik1/ 2, j，V

31、1ik, j1/ 2是加速度分格式为：Pk1/2i, jvxvz,的离散值。则公式(15)的O(t2,x2N)的交错网格差tt Pk1/2i, jtvp2x2N(N)CnUik(2n1)/2, jUik(2n1)/2, jn1Ntvpz(N)CnVi,kj(2n1)/2Vi,kj(2n1)/2n1tvp2xtvp2z(N)CnU1ik(2n1)/ 2, jU1ik(2n1)/ 2, jn1NN(N)CnV1ik, j(2n1)/ 2V1ik, j(2n1)/ 2(24)n1Uki, jUk1i, jtN(N)k1/ 2k1/ 2CnPi(2n1)/ 2, j Pi(2n1)/ 2, j(25)

32、xn1-k1/ 2i, jk1/ 2i, jVVtN(N)k1/ 2CnPi, j(2n1)/ 2 Pi,kj1(/22n1)/ 2(26)zn1U1ki, jtN(N)k1/ 2k1/ 2 CnPi(2n1)/ 2, j Pi(2n1)/ 2, j (27)xn1tN(N)k1/ 2 CnPi, j(2n1)/ 2 Pi,kj1(/22n1)/ 2 (28)zn1V1k1/ 2i, j在公式(24)-(28)中，时间步长是t，x,z分别表示x,z方向的网格间距，vz(x,z zxt,t)为Vi,kj1/ 2，vx(x ,z,t)为Uik1/2, j，u(x,z,t )为Pi,kj1/ 2。2

33、223.2 模拟震源我们在进行正演模拟的过程中，需要加载震源，震源的加载方式以及类型等因素，都能影响正演模拟的结果，因此，我们也应该注意这方面的问题。选择最合适的震源类型以及加载方式。3.2.1 震源的选择在此处,本文选择的是雷克子波作为加载的震源，来模拟激发地震波。表达式（时间域）为：f2222ffMs fe 3fM通过傅立叶变换为频率域：2222221fst2Mt t0expfMt t0其中，初始时刻是t0，主频是fM。我选用频率为 30hz，然后采用时间域的雷克子波的公式做成图，如图 3-4 所示：图 3-4 频率为 30hz 的雷克子波的特征图（时间域）-3.2.23.2.2震源的类型

34、震源的类型正演模拟中有一个重要的部分就是加载震源。它加载的方式，模拟震源的质量，将会对模拟结果产生直接的影响。震源的分类方式有很多，按激发方式，按加载方式，按埋藏深度等等， 方式具体分类不一一列举。 在此， 本文主要介绍一下按加载方式的分类.按加载方式的不同，震源可以分为剪切力源、炸药震源、定向力源、等能量源，下面进行简要的介绍：震源类型震源类型剪切力源剪切力源炸药震源炸药震源定向力源定向力源等能量源等能量源激发方式激发方式网格内四个模拟实际情沿某一方向网格内四个格点 （构成圆形）格点 （构成圆况， 外力施加的某个格点的网格线方向。形）的切向在某些格点） 上各各向向同同性性横波 （类似纯纵波

35、（类似纯纵波、横波介质介质横波震源）纵波震源）纵波、横波纵波和横波（能量相当）纵波和横波（能量相当）各各向向异异性性纵波、横波纵波、横波介质介质补充补充其施加的外力在垂直跟水平方向分量的大小相等外力施加的外力分解到切向或者径向的方向倾斜、 水分量大小一样，由于产生的平、 垂直都可纵横波能量差不多，所以在以。模拟中可以有一些特殊的用途。表 3-5震源按加载方式不同的分类不同的激发方式有不通过的效果， 各有各的有点， 应该根据目的的不同， 认真分析，然后选用最合适的激发方式。在本文中，经过分析我们使用的是纵波震源。3.3 边界条件我们在进行数值模拟时，所应用的模型都是有一定的范围，即存在一定的边界

36、，即人为边界。而在实际情况下，真正的地下介质类似于半无限大的空间，因此，所产生的人为边界对于数值模拟的结果会产生不利的影响。为了尽量减少人为边界的影响，我们可以加入恰当的边界条件来更加准确的模拟无限的真实空间， 尽最大可能去削弱边界反射，改善正演模拟的效果。在历史上，已经有很多科学家进行了相关的研究。科学家 Cerjan 在 1985 年，为了-消除边界反射，他采用了在人工边界引入衰减带的方法，取得了不错的效果。科学家Berenger（法国）在 1994 年提出了 PML(完全匹配层)吸收边界条件，证明利用 PML 边界条件可以吸收不相同的频率、方向的电磁波，且理论上投射到表面的反射系数都为

37、0，即不产生反射，大大的改善了吸收效果，提高了正演模拟的精度。到了二十一世纪，Collino 等人又将此种方法跟弹性波数值模拟联系到一起，效果很好。上面所说的处理方法比较典型，这些方法根据原理的不同可以分成两大类：衰减边界条件、吸收边界条件。3.3.1 衰减边界条件图 3-6 衰减边界条件如图 3-6， 直观的展示了衰减边界条件的基本情况。 上面所说的 Cerjan 提出的方法相比于其他方法更加简单，具有很强的适应性，比较常用。3.3.2 PML 边界条件由于本文应用的是 PML 边界条件，因此，这里主要介绍一下 PML 吸收边界条件。图 3-7 PML 吸收边界条件如图 3-7，比较直观的表

38、现了 PML 边界条件的基本情况。声波方程完全匹配层吸收边界， 它的基本原理是： 将吸收层引入到研究区域的边界上,地震波由研究区域边界出发，当传播到吸收层时不会发生反射,在随着传播距离越来越大，逐渐的衰减（指数衰减）。为了测试边界条件的作用， 本文在此利用一个均匀介质模型， 计算网格大小为 10m-10m，速度采用 3000m/s，采样 2 阶时间，10 阶空间粘声波交错网格有限差分算子对拥有 pml 边界条件跟无边界条件的情况分别进行模拟，得到 t=400ms 时的波长快照，如下图（3-8）：无边界条件有边界条件图 3-8 有无边界条件对比图通过两个图片的对比，我们可以很明显的看出，使用了边

39、界条件的波长快照在边界处没有产生反射波，而没有进行边界处理的反射波很明显。3.4 稳定性稳定性是波动方程求解所必须考虑的，有限差分法正演模拟不能忽略稳定性问题。这个问题的关键是合理的设置参数，因为如果参数设置的不准确，那么模拟就偏离实际值，造成结果不合理，甚至产生溢出现象，无法进行下一步的计算。因此，要研究稳定性问题，就是要使时间采样间隔尽量符合稳定性条件。即，进行数值模拟，必须要分析方法的稳定性，也就是必须要有了解稳定的离散参数区域【9】。下面，本文主要来分析一下二维声波高阶差分模拟方法的稳定性。对于二维声波方程，精度不同，稳定性也不同，不同的精度满足不同的公式，二维声波方程 2n 阶精度的

40、稳定性条件（简写）为：-对上式进行付式变换（时间、空间），得：此时：-20C(N)n为空间差分系数，它是正负号不断改变的，固，要使kx最大，则，kxx（Nyquist 波数）。所以，声波方程不同阶精度的稳定性的条件是：这是稳定性条件的一般形式， 展开后可以分别表示各个精度的稳定性条件 （2n 阶） ：当 2n=2 时，稳定性条件为：VtVt1 0.707。x2111。当x=z 时，左式能够简化为：x2z2当 2n=4 时，稳定性条件为：Vt化为：11 0.866。当x=z 时，左式能够简x2z2Vt 0.612。x11 0.813。当x=z 时，左式能够简22xz当 2n=6 时，稳定性条件为

41、：Vt化为：Vt 0.575。x11 0.784。当当x=z 时，左式能够22xz-当 2n=8 时，稳定性条件为：Vt-简化为：Vt 0.555。x11 0.765。当x=z 时，左式能够简22xz当 2n=10 时，稳定性条件为：Vt化为：Vt 0.541。x在此处，V 代表的是所用模型中的最大速度。从上面的描述中可以看出，随着 n的增加，也就是随着阶数的增加，这种交错网格解法对于稳定性的要求差别不是太大，只是略有提高。因此，我们在高阶差分中也完全可以应用低阶精度的差分网格的相关参数，差分精度对所选择参数的影响不大。3.5 数值频散问题数值频散是有限差分法最大的缺点。在用有限差分法进行数值

42、模拟时，有很多因素会造成数值的频散，例如：1）震源模拟精度不高，造成了波形的畸变甚至产生重影现象；2）差分的时间采样和空间网格间距之间不太合适；3）当使用角度很低的有限差分方程来进行陡倾角成像使波场发生畸变；4）横向速度的剧烈变化造成的波长畸变【10】。所谓的数值频散现象，是指波前的形状在传播过程中由于某些原因发生变化、分散，这种现象是由于波的群速度和相速度不一致造成的。我们引入误差项，用于差分方程向微分方程的逼近。有时候引入的误差项会产生很多不好的影响，例如使相速度发生变化或者使结果的振幅值衰减，结果类似于频散或者物理耗散，这种物理效应是虚假的，可以称之为数值频散或者数值耗散。此类频散是差分

43、方程原有的特点，它跟波动方程造成的物理频散不一样【11】。地震波数值模拟效果由数值频散程度所直接决定，正演模拟过程中，消除或者压制数值频散的方法可以从三个方面来考虑，一是可以通过高阶差分压制频散；二是可以通过对差分算子进行校正来压制数值频散；三是通过 FCT 方法压制数值频散【12】。本文中是利用的第一种方法，通过提高差分方程的阶数来压制数值频散。下面我用不同阶数的差分方程所产生的波长快照来验证一下这种方法的正确性。本文通过建立一个均匀介质模型（二维声波各向同性完全弹性介质），仅用于测试差分算子阶数跟数值频散的关系。利用不同阶数的高阶有限差分算法进行了数值模拟，得到波场快照如图 3-9 所示。

44、- 2 阶精度 4阶精度 10 阶精度 12阶精度图 3-9 不同阶数高阶有限差分模拟波场快照从图中可以看出，在阶数较低时，出现很多同相轴，说明数值频散现象严重；对于高阶差分算子来说，算子阶数越高，压制数值频散效果越好，精度越高。3.6 模型试算3.6.1 地震波的传播规律研究为了模拟粘声波在不同水平层状介质中的传播， 研究不同介质模型情况下地震波的传播规律，采用一个十阶空间、两阶时间的交错网格有限差分算法，建立了三个简单的三层分布的模型。然后根据得到的波长快照来进行分析研究，假定上述水平层状模型的大小是 5000m5000m，其中，第一层为2000m，第二层为1000m，第三层为 2000m

45、。在这个模型中，震源的频率为 15hz，三层的速度、密度、品质因子的设定如下表(表 3-10）所示：层序vp(m/s)(g /cm3)Qp5570第一层2000第二层30002.02.2-第三层40002.585模型 1：正序水平层层序vp(m/s)(g /cm3)Qp557085第一层2000第二层第三层400030002.02.52.2模型 2:含高速夹层的水平层层序第一层vp(m/s)(g /cm3)Qp55708530002.22.02.5第二层2000第三层4000模型 3：含低速夹层的水平层表 3-10 不同水平层状模型的设定参数波场快照是用来研究地震波传播特点的一种手段， 它表示

46、的是某一特定时间地震波波前面的真实位置。为了了解粘声波在介质中的传播情况，我首先利用模型 1（正序水平地层）来分别得到 t=200s，t=300s，t=400s 的波场快照（x 方向），如图 3-12 所示。震源位置以及地层厚度如图 3-11 所示。2000m震源10002000m5000m图 3-11 水平层状介质模型地层厚度以及震源位置-t=200mst=300mst=400ms图 3-12正序水平地层不同时刻的波长快照（x 方向）图中可以看出，在未遇到界面前，地震波在均匀介质中的波前面一个圆。当遇到地层界面之后，在界面处发生了反射、透射和折射现象。沿测线方向的地震记录如图3-12所示。记

47、录中存在两条直线状的同相轴和两条近似双曲线的同相轴。由于直达波的时距曲线是直线，因此两条直线同相轴对应直达波；由于反射波的时距曲线是近似双曲线，因此近似双曲线同相轴对应的是反射波。1233-12 二维粘滞声波方程数值模拟记录其中，【1】为直达波，【2】为第一层反射波，【3】为第二层反射波。为了进一步研究地震波的传播规律，对模型 2、模型 3 进行模拟，震源位置改为x=2500m,y=2500m 处 （见图 3-13） ， 当 t=300ms 时的波长快照 （沿测线方向） 见图 3-14。-2000m1000震源2000m5000m图 3-13 水平层状介质模型地层厚度以及震源位置模型 1 顺序

48、地层模型 2含高速夹层模型 3 含低速夹层图 3-14 各模型在 t=300ms 时的波场快照从图 3-14 中能够看到：（1）当地震波入射到两种不同介质的分界面时，会产生透射波跟反射波；（2）在波前面位置产生极性反转现象，沿测线方向（水平方向）极性反转的方向是左右反转；（3）从图 3-14 中的模型 3（含低速夹层）中可以看到，当波入射角为临界角的时候，二三层的分界面出产生了折射波，这是因为第二层介质的速度为 2000m/s，远小于第三层介质的速度（4000m/s），符合折射波产生的条件。3.6.2 粘声波的衰减规律研究为了了解地震波在粘声介质中的传播，我重建了一个二维粘滞声波四层介质模型，

49、选用的模型大小为 3000m3000m，计算的网格大小为 10m10m，采用两阶时间差分，十阶空间差分的交错网格粘声波正演模拟算法，震源的主频为 30hz,dt 为 1ms。具体的-参数设置如下图（3-15）。层序vp(m/s)(g /cm3)Qp90909090第一层2500第二层第三层第四层3000350040001.82.12.42.7表 3-15 四层介质模型的参数设定通过进行模拟，得到粘声波方程正演模拟地震记录。然后与普通的二维 声波方程的正演模拟记录进行对比，如图 3-16（纵坐标为时间，横坐标为道数）。粘声方程声波方程图 3-16 二维声波方程与二维粘声方程的模拟记录对比图然后选

50、取第 30 道的声波记录跟粘声记录，放到一个表格中进行对比，如图3-17:图 3-17第 30 到声波记录跟粘声记录对比图-从上图 3-16 与 3-17 中我们能够得到： （1）粘声波的地震记录在深层的反射波的振幅明与普通声波的相比，有明显的衰减现象，同相轴相对教弱；（ 2）在深层，粘声波的记录的波形变化比较严重。也就是说，粘性介质考虑到了地层对能量的衰减作用，反射波的能量比完全弹性的介质减弱很多，更加的符合实际的底层情况。-第 4 章结论本文在在前人研究的基础上，通过粘声波正演模拟的方法，了解了地震波的传播规律、衰减规律以及波在粘性介质中的一些特点等等。本文应用的模型为开尔芬模型，通过正演

51、模拟，得到了以下几个结论。（1）我们在进行正演模拟的过程中，震源的加载方式以及类型等因素，都能影响正演模拟的结果，我们应该注意这方面的问题。选择最合适的震源类型以及加载方式。（2）使用了边界条件的波长快照在边界处没有产生反射波，而没有进行边界处理的反射波很明显。（3）对于有限差分法，一个很大的问题是数值频散问题，在阶数较低时，出现很多同相轴，数值频散现象严重；对于高阶差分算子来说，算子阶数越高，压制数值频散效果越好，精度越高。（3）在粘性介质中，波的衰减情况跟品质因子密切相关，而品质因子与震源频率有关，震源频率越大，品质因子越小，吸收越严重，也就是衰减越厉害。（4）当地震波入射到两种不同介质的

52、分界面时，会产生透射波跟反射波；并且在波前面位置产生极性反转现象， 当含低速夹层， 波入射角为临界角的时候， 在低速层 （上层）跟高速层（下层）的分界面会产生折射波。（5）粘声波的地震记录在深层的反射波的振幅明与普通声波的相比，有明显的衰减现象，同相轴相对教弱；在深层，粘声波的记录的波形变化比较严重。也就是说，粘性介质考虑到了地层对能量的衰减作用，反射波的能量比完全弹性的介质减弱很多，更加的符合实际的底层情况。在本文所介绍的研究中，还存在很多的不足之处：（1）开尔文模型本身存在着某些缺点，仅仅应用这个模型进行数值模拟不够精确，高频成分吸收严重如果应用广义标准线性体，即将开尔文模型与麦克斯韦模型

53、结合起来，这样的话能更好地进行模拟，更加精确的模拟地震波在粘性介质中的传播规律。（2）在这里，对于粘性介质中地震波的衰减的研究比较简单，只是通过与声波进行了地震记录与某一时刻振幅的简单对比，研究的较少。另外，还缺少一些其他方面，例如自由表面处理等方法的研究，需要继续学习的地方还有不少。对未来的展望：（1）q 值的提取这一块对于地震波能量的补偿意义重大。（2）在进行边界处理的过程中所应用的广义标准线性体中，为了避免应变应力复-杂的卷积关系，可以引入记忆变量。在粘声正演这一块，我需要掌握更多的知识，还有很长的路要走。我要进一步查阅有关文献，了解相关的学术动态，透彻的理解相关的理论知识，不断的提升自

54、己-致谢致谢大学时光在不知不觉中就这么结束了，在这里，我真诚的感谢短暂的四年里给我巨大的鼓励和帮助的老师、同学、朋友。在进行毕业设计的期间，我首先要感谢的是xxx 老师的帮助。他在我毕设的过程中给我提供了宝贵的意见。他虽然远在美国，但是仍然定期的对我进行督促，检查我的进展情况。他严谨的工作态度，对知识的努力追求，以及严以律己宽以待人的处事方式，无时无刻的影响着我。是我以后学习的标杆。其次，我也衷心的感谢xxx 给我提供的帮助。在我毕设的过程中，从初期的相关文献的检索，到后期论文的编写，都给了我巨大的帮助，耐心的给我讲解相关的程序与知识，让我得以顺利的完成自己的毕设。最后，感谢各位查阅的老师，抽

55、出时间来评阅我们的论文。还有很多没有提到的、帮助过我的人，我在心里默默祝福你们，对你们表示最诚挚的谢意！-参考文献参考文献1 郭成锋、杜启振，二维 VTI 介质粘声波正演模拟，中国地球科学联合学术年会，2014 年，P12372 廖建平、 刘和秀等，二维频率空间域粘声波正演模拟研究，地球物理学进展，2011年 12 月，第 26 卷第 6 期 P2075-20813 孙天真， 谷玉田等， 基于粘声介质的反滤波叠前深度偏移方法研究， 石油物探，2013 年月，第 52 卷第期，P275-2794 廖建平、刘和秀等，精确的频率空间域粘声波正演数值模拟，中国地球物理，2011年，P7265 陈可详，

56、标量声波波动方程高阶交错网格有限差分法，中国海上油气， 2009 年 8月，第 21 卷第 4 期 P232-2366 贾立坤，粘弹性介质地震波数值模拟，中国石油大学（华东）硕士论文，2011 年 5月7 刘庆敏， 高阶差分数值模拟方法研究与应用， 中国石油大学 （华东） 硕士论文， 2004年 9 月8 孙成禹、李振春，地震波动力学基础，石油工业出版社，2011 年 4 月9 董良国，马在田等，一阶弹性波方程交错网格高阶差分解法稳定性问题，地球物理学报，2000 年，P85686410 Jean Virieux，P-SV wave propagation in heterogeneous m

57、edia: Velocity-stressfinite-difference method ，GEOPHYSICS ，April 1986 ，number4 P889-90111 Alan R. Levander，Fourth-order finite-difference P-W seismogramsGEOPHYSICS, November 1988，VOL. 53. NO. 11P1425-143612 Wang TL，Tang XM，Finite-difference modeling of elastic wave propagation：A nonsplitting perfectly matched layer approachGeophysics，2003，P1749-175513 JohanO.A.Robertsson.Viscoelastic finite-difference modeling J.Geophysics, 1994, V ol59 NO.9 P1444-145614 Carcione.Wave propagation simulation in a linear viscoelastic medium.GeophysicalJournal,1988,95,597-11-