基于三阶Adams格式的求解声波方程的多步算法

1、 China University of Mining & Technology-Beijing创新项目论文一种基于三阶Adams格式的求解声波方程的多步算法 摘 要一个准确、高效、低数值频散的正演算法能够提高反演精度、加快反演收敛速度，因此研究地震波场正演模拟技术具有重要意义。区别于传统的空间离散方法，利用空间插值, 用网格点处的函数值及其梯度共同逼近空间高阶偏导数的方法称为近似解析离散化方法。声波方程通过变换，并采用近似解析离散化方法进行空间离散，从而转变成为一个半离散化的常微分方程组，再利用三阶显式Adams格式进行时间推进，求解半离散化的常微分方程组，从而得到了一个新的求解声波

2、方程的有限差分方法（AD-STEM）。对AD-STEM进行了理论误差和数值误差分析、计算效率比较和数值波场模拟。研究表明，与传统方法AD-LWC比较， AD-STEM方法数值精度更高，数值频散更低，更高效，且与解析解匹配更好。AD-STEM方法能够通过压制数值频散而提高计算效率。在无可见数值频散的条件下，AD-STEM的计算速度是AD-LWC的1.88倍，而存储量只有其72%，更适合在粗网格下进行大规模地震波场数值模拟。关键词：近似解析离散化方法；三阶Adams格式；数值频散；有限差分目 录1 绪论11.1选题背景和研究意义 1.2粘弹性介质国内外研究现状1.3有限差分国内外研究现状1.4本文

3、主要研究内容2 粘弹性介质的基本模型33 方法介绍33.1 Stereo-modeling方法简介33.2 Lax-Wendroff correction方法简介54 粘弹性介质中的波场数值模拟64.1 波场快照64.2 波形图64.3 SEG模型的地表地震记录105 结论136 参考文献15中国矿业大学（北京）2016届本科毕业设计（论文）1 绪论1.1 选题背景和研究意义 随着我国石油和天然气工业的迅速发展，油气勘探的精度和难度也在不断地加大。这就要求我们必须深入了解地震波在地下介质中的传播规律，以便我们能更加准确地进行地下油气的勘探。地震波数值模拟技术是目前地震勘探领域人们研究地震传播规

4、律的一个重要手段。传统地震波数值模拟技术，一般假设地层介质是均匀和完全弹性的，但是我们发现：在这种假设情况下，数值模拟的结果和我们在野外观察到的实际记录有很大差别。造成这种差别的主要原因是因为实际的地层介质并不是均匀和完全弹性的，我们应该用粘弹性介质来代替，这对于我们研究地震波波场的传播规律具有了非常重要的理论和实际意义。1.2 粘弹性介质国内外研究现状 为了对实际地震资料的解释以及地震波的偏移反演等处理提供可靠依据，我们需要建立不同的介质模型来模拟地下介质，发展相应的理论和勘探方法。从以往的模拟结果来看，以经典弹性波理论为基础得到的理论记录实际地震波在大地中传播时所接收到的记录之间存在很大差

5、别，为了消除这种差别，我们通常要对球面扩散、界面的透射损失等进行补偿，但是在补偿过后，理论记录上来自深层的地震波振幅仍然较小、低频成分较丰富，这说明地下介质对地震波有吸收衰减的作用，特别是对高频成分，吸收更为严重。地震波衰减受多种因素影响，事实上现在没有一种机理可以描述所有环境条件下产生的损耗，介质的粘滞性就是其中的一个主要原因，为了弄清楚地震波在地层中的衰减机制，很多研究人员开始对粘弹性介质的理论和应用方法进行研究。（1） 国外研究现状 早在1845年Stokes就开始研究粘弹性介质中的地震波，他提出要考虑具有耗损的固体就必须使这种固体具有类似于粘滞液体的性质，因此他给固体的剪切模量附加上一

6、个剪切粘滞，并建立了由于粘滞型内摩擦引起的能量耗损的Stokes粘弹性波动方程。但是粘弹介质理论的发展却是非常缓慢的，直到上个世纪40年代，美国地球物理学家N.K.Ricker首次在粘弹介质理论中完整详细的描述了地震波在地下介质传播中时的衰减问题。其后，人们纷纷开始研究粘弹介质中地震波的传播理论，提出了一系列的理论和研究方法。Aki和Richards(1980)建立了一种新的粘弹性介质模型标准线性固体模型，相对于开尔芬佛格特模型和麦克斯韦尔模型，标准线性固体模型更加适合描述地震波衰减的物理机制。粘弹介质地震波数值模拟技术的研究是从上世纪八十年代开始的，Day和Minster (1984)第一次

7、成功的在粘弹性介质中利用Pade近似法进行2-D时间域数值模拟。EmmeriCh和Korn (1987)发现这种方法存在明显的缺陷即质量低劣和计算效率低，于是他们提出了一种新的近似称之为“广义标准线性固体”(GSLS)，同时对粘弹性模量的有理式进行推导，并发展二维有限差分算法使之适合标量波的传播，提高了计算效率。CarCione(1988)等对粘滞声波在地层中传播的进行了模拟并推导出了模拟方程。Robertsson(1994)等设计了一种速度一应力有限交错网法，并利用它研究了二维粘弹性介质下纵波和横波的传播规律，随后他们又将该算法推广到三维介质情况。(二)国内研究现状跟国外关于粘弹性介质理论及

8、其数值模拟的迅速发展，国内对这方面研究起步较晚。近年来国内地球物理学家也开始对地震波在粘弹性介质中传播进行研究。张剑锋、李幼铭(1994)把地层介质假设为水平层状介质，利用混合变量粘弹性波方程，直接逐层求解位移、应力的方法进行数值模拟四。毕玉英、杨宝俊(1995)给出一种实现二维粘滞介质完全波场模型计算的方法。该方法的独到之处在于将传播时间分解成了传播水平时间和传播垂直时间两部分，平面源人射改为线源人射，无需繁杂的射线追踪，只考虑入射角随偏移距的变化情况便可获得包括反射纵波、转换波、多次波、折射波及面波等在内的多种波场的时间一空间道集的正演模拟记录，而且还能灵活地模拟井间和VSP地震剖面。宋守

9、根(1996)提出了一种提高地震剖面纵横向分辨率的粘弹性介质波场外推方法，利用该方法进行反演时，无需对积分方程作线性化近似，适应性较强，并通过实例对该理论推断予以验正。张建贵(1999)等针对塔里木盆地的沙漠覆盖面大，地层埋藏深，地质构造幅度小的特点以大地介质为粘弹性介质为前提，利用粘弹性介质波动方程聚焦成像技术得到了一套高分辨率、高信噪比和高保真的地震处理剖面。孟凡顺(2000)等人根据粘滞弹性理论，推导出了粘滞弹性波的有限差分运动方程并对任意复杂地质体粘滞弹性波进行了正演模拟；程昌钧(2001)等系统研究了粘弹性介质中波的逆散射问题及其求解过程；崔建军和何继善 (2001)以二维粘弹性波动

10、方程为基础，研究了粘弹性波动方程F-K域的正演模拟和偏移方法的思路和理论基础，并阐明了推导过程；杨午阳(2003)提出了一种利用粘弹性声波波动方程进行偏移的新方法；奚先和姚姚(2004)通过波动方程的交错网格有限差分数值模拟方法，对地震波在二维粘弹性随机介质中的传播进行了模拟并研究了其波场特征；朱慧卿(2004)利用新定义的各向异性网格一广义拟一致网格对粘弹性方程各向异性有限元方法的超收敛性进行了分析22；刘财、张智(2005)以积分本构方程为基础，应用对应原理建立波动方程，对线性粘弹体中地震波场进行伪谱法模拟；宋常瑜(2006)等采用二维粘弹性波方程交错网格高阶有限差分数值模拟方法进行了井间

11、地震粘弹性波场数值模拟，并对地震波的传播规律和衰减特征进行了研究；邵志刚(2007)将以往两种粘弹介质中地震波模拟方法的优点结合起来，以模型理论和积分本构方程为基础，从理论上分析了模型对地震波场的影响并采用交错网格有限差分法对粘弹介质中的地震波进行数值模拟；孙成禹(2007)针对以往粘弹性模型在描述介质品质因子对频率的依赖关系方面存在不足，与实际观测结果不符的问题，提出了一种新的粘弹介质模型的数值方法，该方法较好地描述了品质因子对频率的非依赖性，可以从理论上较准确地对实际介质的粘弹性对地震波场的影响进行研究；唐启军(2009)将Von Karman型的随机各向同性背景引入粘弹性单斜各向异性波动

12、方程，并应用交错网格技术进行模拟；Yong Yun-dong (2010)等人推导出了交错网格有限差分的并行算法，并将其用于三维粘弹性随机溶洞介质数值模拟中，取得了良好的效果。1.3有限差分国内外研究现状 地震数值模拟作为地震勘探基础性研究，是我们认识地震波传播规律，检验各种处理方法正确性的重要工具，地震波的数值模拟是是研究地震波传播规律的最有效工具和手段，在地震勘探和地震学各工作阶段都有重要的作用，贯穿于地震资料的采集、处理、解释的整个过程中，同时也是地震反演的基础。有限差分法是最常用的一种数值模拟方法，基本原理就是在解偏微分方程时用有限差分算子代替微分，将微分方程化为相关的线性代数方程，通

13、过求解代数方程，得到偏微分方程的数值解。有限差分法数值模拟技术相对于射线方法具有更高的精度，同时比有限元方法计算量小，因此在实际应用中占很重要的地位。(一)国外研究现状有限差分法数值模拟技术的发展只有几十年的历史，早期对于地震波传播有限差分方法的研究都是基于二阶位移方程。有限差分数值模拟方法最早由Alterman和Karal(1968)提出，他们首次将有限差分数值模拟方法应用于层状介质的模拟，并分析弹性波在层状均匀介质中的传播规律。他们为利用差分方法解决各种勘探地震学的实际问题奠定了基础，随后差分方法被广泛应用于地震勘探中并在应用中不断得到发展。Boore(1972)进行了非均匀介质地震波有限

14、差分数值模拟并研究了地震波非均匀介质中的传播规律。Alford (1974)对有限差分模拟的影响因素(网格大小和地震波传播方向)进行分析，网格间距控制地震波数值模拟的计算精度，相同精度下高阶差分与低阶差分对网格间距的要求不同型;Kelly(1976)利用有限差分法进行地震记录的合成，还对合成的地震记录进行一些常规的数据处理分析，进一步拓宽了数值模拟在实际地震处理中的应用；前面都是基于笛卡尔坐标系下常规网格的离散差分，Madariaga(1976)首次提出交错网格方法，并将此方法应用于弹性介质地震波的模拟，其模拟的差分精度为O(t2+x2)。virieux (1984)提出了应力一速度一阶方程交

15、错网格有限差分法，用一阶速度一应力方程来代替二阶位移方程，由于一阶速度一应力方程不需要对弹性参数进行求导，因此它更加适用于复杂介质的模拟。 Levander对Virieux (1988)的方法进行了发展，采用四阶空间有限差分算子模拟弹性波的地震记录，进一步提高模拟精度。在解决计算精度和运算效率这一矛盾方面，Dablain(1986)利用Taylor级数给出了声波方程高阶差分格式，成功消除了空间差分频散对数值模拟的影响，极大的提高了有限差分数值模拟的精度，同时还分析了时间和空间离散的频散特征。Crase (1990)利用该方法推导了二阶弹性波方程高阶差分形式，并给出了任意阶的高阶交错网格差分方法

16、。前面讲到的差分方法它们的网格步长都是固定的，Jastram和Tessmer(1994)提出了垂直间距可变网格差分方法，并将此方法应用于弹性介质地震波的模拟，Graves (1996)将交错网格二维模型模拟推广到三维模型的模拟，该方法模拟能力有了明显提高。Pitarka（1999)在三维交错网格基础上发展了矩形非规则交错网格，并进行了三维弹性波的模拟；在有限差分算法中除了对网格步长的控制外，还可以通过变时间步长来减少计算量提高计算效率，Tessmer(2000)提出了二阶运动方程的变时间算法；为了进一步模拟地震波在非完全弹性的实际地层中的传播，Carcione(1988)等利用交错网格模拟了粘

17、滞声波在地层中的传播；Saenge和Bohlen(2004)提出了旋转交错网格有限差分方法，进一步提高了交错网格对复杂介质模拟的精确度，并将旋转交错网格应用于三维各向异性介质中的研究，carcione和Helle(1999)提出孔隙介质中粘弹性波交错网格有限差分数值模拟。(二)国内研究现状 何樵登(1996等对横向各向同性介质中地震波数值模拟方法进行了系统介绍；董良国等(2000)详细阐述了一阶速度应力方程交错网格高阶有限差分的实现方法，从理论上证明了交错网格具有较高的计算精度和计算效率；裴正林等(2005)进一步利用任意偶数阶精度交错网格有限差分实现了对不同复杂介质的模拟；殷文、印兴耀、吴国

18、忱等(2006)提出了高精度频率域弹性波方程有限差分方法，从模拟的结果看，该方法能很好地消除频散和边界反射；李胜军、孙成禹等(2007)等在不做插值处理的前提下实现了垂向任意倍数网格步长比高阶有限差分算法，并解决了高速层中含有低速层在模拟时存在的数值频散和稳定性不能均衡的问题。1.4本文主要研究内容2 粘弹性介质的基本模型 1.粘弹性的模型表述 介质的线粘弹性质，主要表现为粘滞性和弹性两方面，这两方面的性质可以由离散的弹性元件弹簧与粘性原件阻尼器进行描述。弹簧和阻尼器以不同的方式进行组合可以得到具有不同性质的粘弹性模型。弹性元件弹簧服从胡克定律 =E或=G 其中为正应力、为剪应力、为正应变、为

19、剪应变；E为拉压弹性模量，G为剪切弹性模量。我们可以看出应力和应变的关系与时间无关，不随时间的变化而变化，弹性变形和弹性回复都是瞬间完成的(图 2-4)。 图 2-3 弹性元件模型 图 2-4 弹性元件应力、应变 粘性元件即阻尼器，有时也被称为粘壶(图 2-5)，它服从牛顿粘性定律： = 或 图 2-5 粘性元件模型 粘弹介质交错网格有限差分数值模拟 图 2-6 粘性元件应力、应变响应 （左边应变响应、右边应力响应） 其中， 或为粘性系数；为应变率 阻尼器的流变特性，我们可以通过观察等应力和等应变作用下它们的响应来了解在应力 =H(t) 作用下，应变响应为即呈稳态流动(图2-6)其中 H(t)

20、 是单位阶跃函数， 关于 t=0 时的数值，将在后面具体问题中加以说明。 在(t) =H(t)的条件下，由式(22)得应力 =(t) )，式中 (t) 为单位脉冲函数，它满足两个条件： 因此，阻尼器受阶跃应变作用时，应力为无限大而后瞬即为零(图 2-6)。由于不可能产生一数值为无限大的力，所以实际上不能瞬时地使粘性元件产生有限应变。 2 麦克斯韦尔介质模型 麦克斯韦尔（Mexwell）介质模型由一个弹簧和一个阻尼器串联而成(图 2-7)。设在应力 作用下，弹簧和阻尼器的应变分别为和，则总应变为： 图 2-7 麦克斯韦尔（Mexwell）介质模型 得： 或 式中，和均表示材料常数。式(2-9)即

21、为麦克斯韦尔（Mexwell）介质模型的应力应变微分关系式。 3 开尔芬介质模型 开尔芬（Kelvin）介质模型由一个弹簧和一个阻尼器并联而成如图 2-8 所示。两个元件的应变都等于模型的总应变，而模型的总应力为两元件应力之和，即: 代人上式，得: 或这就是开尔芬（Kelvin）介质模型的本构关系。式中，图 2-8 开尔芬（Kelvin）介质模型 4 标准线性粘弹性介质模型 如图 2-9 所示，标准线性粘弹性介质是由一个弹簧与一个开尔芬粘弹性单元体相互串联而成的三元固体。 图 2-9 标准线性粘弹性介质模型 对于两个组成元件，分别有： 和总应变为: 进行拉普拉斯变换及其逆变换就可得到标准线性粘

22、弹性介质微分形式的本构方程，如下式所示： 3 方法介绍3.1 Stereo-modeling方法简介Stereo-modeling方法与传统有限差分方法的主要区别在于空间离散形式，时间推进可以与任一时间格式结合。传统方法的空间高阶偏导数近似可用下面的公式表示8： （3.1）其中m是被近似的高阶导数的阶数。当应用n+1个节点近似函数f(x)的m阶导数时，是在节点j处的系数。m阶偏导数的近似，需要的节点个数依赖于m。当近似的导数阶数m从1增加到M时，使用的节点个数n从m增加到N个。将这种空间离散格式从一维推广到二维和三维的过程中，每个维度只用到沿着坐标轴方向的节点，是简单的线性推广。Stereo-

23、modeling方法对于空间高阶导数的近似，可以用下面的公式表示5-7 (3.2) 从这个公式中可以看出，Stereo-modeling方法对空间高阶导数的近似，不仅应用了函数值，还应用了函数的一阶导数。所用的节点个数与具体的导数阶数m并非直接相关，而是取决于算法的整体空间精度。具体关系是，应用2L+1个节点时，可以取得4L阶精度。在下文中我们将采用具有4阶空间精度Stereo-modeling方法进行方程的空间离散。3.2Lax-Wendroff correction方法简介 为了说明新方法STEM的有效性和高效性，我们选取经典传统方法Lax-Wendroff correction（LWC）

24、与其作对比。更确切地说，是选取LWC方法的空间离散格式进行波动方程的空间域离散，对于LWC方法，这里只列出参与计算的具有二阶精度的二阶空间偏导数的离散格式 ： (3.3) 4 粘弹性介质中的波场数值模拟 接下来应用上述两种方法求解以下的二维声波方程，来进行波场数值模拟： （3.4）在开尔芬粘弹性介质中，对应规则为 因此上述二维波动方程变为 （3.5） 其中， ，Q为品质因子将（3.5）式同除，记，则 设，=，即，则， 在该数值模拟实验中,网格点数为201×201, 空间网格步长为x = z = 50 m，震源中心频率为f0=18 Hz，介质中声波速度为c0=4 km/s，时间步长为t

25、=0.002 s。在时间时的入射角。图2.1显示了由公式（2.19）计算得到的相对误差曲线，其中红色虚线为AD-LWC格式所得，黑色实线为AD-STEM格式所得。可以看到AD-STEM格式的误差小于AD-LWC格式的误差，尽管两种方法的理论误差相同。 这是一个置于计算区域中心的爆炸源，其主频率为=20 Hz，声波波速是c=4 km/s,计算区域为,空间步长为，网格点个数为241×241，时间步长为0.002 s。图2.2a和图2.2b分别显示了由AD-STEM和AD-LWC方法得到的在粗网格下的波场快照。从图中可以看出，AD-STEM方法没有明显的可见数值频散（图2.2a），而AD-

26、LWC方法则产生了严重的数值频散（图2.2b）。为了更清楚地观察波场模拟情况，在计算区域中的（6.5 km，6.5 km）位置放置一个接收点，观察接收到的波形记录图。图2.3和图2.4分别显示了AD-STEM和AD-LWC方法在相同计算条件下，对应于图2.2a和图2.2b在相同接收点处的波形记录图。从这两个图中，可以更清晰的看到，AD-STEM方法获得的数值解与解析解吻合较好（解析解采用De hoop11等提出的方法计算得到），而AD-LWC方法计算得到的波形则与解析解相差甚远。为了消除AD-LWC方法获得的波场快照中的数值频散，我们采用加细网格的方法。在保持空间步长与时间步长的比值不变的前提

27、下，将空间网格步长加细到15 m时，数值频散才彻底消失（图2.2c），此时AD-LWC方法的计算时间为250 s，网格点数为801×801，每个网格点上存储8个数组。而AD-STEM方法在粗网格50 m的条件下已无可见的数值频散，用时为133 s，网格点数为241×241，每个网格点上存储24个数组。可见在无可见数值频散的条件下，AD-STEM方法的计算速度是AD-LWC方法的1.88倍，而存储量仅为AD-LWC的72%，具体数据参见表一。表一：Stems方法与LWC方法计算效率比较数值方法网格步长（m）CPU TIME (s)存储量计算效率AD-STEM50133241&

28、#215;241×241.88AD-LWC15250801×801×80.720 距离（km） 120 深度（km） 12(a) AD-STEM（m）0 距离（km） 120 距离（km） 120 深度（km） 120 深度（km） 12 (b) AD-LWC（m） (c) AD-LWC（m）图2.2 均匀各向同性介质中AD-STEM与AD-LWC方法所得到的波场快照图2.3 AD-STEM方法得到的数值解与解析解比较图2.4 AD-LWC方法得到的数值解与解析解比较2.3.2 波场数值模拟本节将在较复杂的均匀横向各向同性介质（两层模型）中考察比较两种方法的数值性

29、质。图2.5为两层模型的速度分布示意图。如图所示，声波在上层介质(km)和下层介质（）中的速度分别为2.4 km/s和5.0 km/s。节点数为501×501，模型规模是20 km×20 km，空间步长为40 m，时间步长为0.001 s。震源仍采用上一节中用到的震源子波，其中主频率=15 Hz，震源放置在（10 km，11.2 km）处，接收器R放置在（10 km，10 km）处（分别如图2.5的五角星和R）。图2.6为在1.8 s时由AD-STEM（图2.6（a）和Ad-LWC（图2.6（b）分别得到的波场快照。从图2.6可以明显看出，在相同的空间步长和时间步长的条件下

30、，AD-STEM可以有效的压制数值频散，而传统AD-LWC方法则有明显的数值频散。为进一步观察数值频散情况，绘制接收点R处的波形记录图，且与利用王美霞等13提出的计算方法所得到的解析解进行比较。图2.7（a）和（b）分别为AD-STEM方法和AD-LWC方法对于接收器R所接收到的波形图。这两幅波形对比图显示，AD-STEM方法在粗网格（50 m）的情形下，获得的波形记录的两个相位与解析解匹配的都较好，误差很小，而AD-LWC方法与解析解吻合不好，误差很大。该结果与图2.6的波场快照一致。 图2.5 双层模型速度分布示意图0 距离（km） 20 0 距离（km） 20 (a) AD-STEM (

31、b) Ad-LWC0 深度（km） 20 图2.6 （a）AD-STEM方法和（b）Ad-LWC方法得到的双层模型波场快照图2.7 （a）AD-STEM和（b）AD-LWC方法数值解与解析解比较3 结论本文提出了一种新的数值离散方法，AD-STEM方法。区别于传统的数值离散方法，AD-STEM方法在空间上采取Stereo-modeling方法，时间上采用3阶Adams格式，从而得到了一个具有高精度、高效率、低数值频散的新方法，对于该方法，可获得以下结论：（1）数值误差分析显示，虽然理论误差相同，但AD-STEM方法比AD-LWC方法的数值精度更高；（2）计算效率比较数值实验显示，AD-STEM

32、在粗空间网格步长（50 m）下，已经能较好的压制数值频散，在波形记录图中与解析解也匹配较好。而同等条件下的AD-LWC方法数值频散十分明显，波形记录图与解析解匹配不好。在无可见数值频散的条件下，AD-STEM方法的计算速度约为AD-LWC的1.88倍，而存储量为其0.72%。（3）较复杂的双层模型波场模拟显示，在同等计算条件下，AD-STEM方法可以有效地压制数值频散，获得清晰准确的波场快照，波形与解析解吻合，而AD-LWC方法产生较严重的数值频散，波形严重震荡，与解析解不匹配。综上，本文发展的AD-STEM方法是一种数值性质优良且十分实用的数值离散方法，尤其可以利用粗网格进行大规模波场模拟。

33、参考文献1 杜世通. 地震波动力学M. 石油大学出版社，1996.2 杨桂通，张善元. 弹性动力学M. 中国铁道出版社，1988. 3 冯英杰，杨长春，吴萍. 地震波有限差分综述J. 地球物理学进展，2007.4 佘德平. 波场数值模拟技术J. 勘探地球物理进展，2004，27(1)：1621.5 Yang D，Lu M，Wu R，et al. An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations J. Bulletin of the Seismological Society

34、 of America，2004，94：1982 1992. 6 Yang D，Peng J，Lu M，et al. Optimal nearly analytic discrete approximation to the scalar wave equation J. Bulletin of the Seismological Society of America，2006，96：11141130. 7 Yang D，Teng J，Zhang Z，et al. A nearly analytic discrete method for acoustic and elastic wave e

35、quations in anisotropic media J. Bulletin of the Seismological Society of America，2003，93：882 890. 8 Fornberg B. Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of computation J. 1988，51：699 706. 9 牛滨华，孙春岩. 半空间介质与地震波传播M. 北京：石油工业出版社，2002.10 李静爽. 基于近似解析离散化算子的逆时偏移方法及其应

36、用D. 北京：清华大学数学科学系，2014.11 Dablain M. A. The application of high-order differencing to scalar wave equation J. Geophysics，1986，51：54 66.12 De Hoop A. T. A modification of Cagniards method for solving seismic pulse problems J. Appl. Sci. Res. B8，1960：349356.13 王美霞，杨顶辉，宋国杰. 二维SH波方程的半解析解及其数值模拟J. 地球物理学报，2

37、012，55(3)：914 924.致 谢毕业论文的完成，代表着四年的大学生活划上了句号。首先，我要感谢我的论文指导老师，李静爽老师。无论从前期的准备工作，或是论文整体结构的设计，乃至论文具体数学理论的讲解和程序的调试，李老师无时无刻不在给予我最大的帮助。甚至在我周中由于实习的原因而无法同老师探讨论文问题时，李老师牺牲掉周末的休息时间来帮助我解决问题，修改论文。可以说没有李老师的帮助，就不会有这篇论文的诞生。其次，我要感谢大学四年时间里传授我知识的每一位老师和辅导员。是你们不辞辛苦、呕心沥血的谆谆教诲才让我在大学的时光里开拓了眼界，增长了学识。我还要感谢大学四年里身边的每一位同学，每一个朋友。

38、是你们无私的关怀与帮助让我在最美好的年华里收获了最纯真的友谊。最后，我要感谢我的父母和家人。在我大学四年的求学生涯中一直给予我物质和精神上的帮助，让我无忧无虑的专心求学。希望通过这篇论文，能作为我大学四年学习生活的一个总结，能当作对老师，对家人，对自己的一个交代。由于我的学术水平有限，本论文难免有不足之处，学生恳请各位老师同学不吝批评和指正。 附录1 设备清单本文所进行的一切数值实验均在该电脑上完成，具体参数为：电脑型号：惠普 HP ENVY 4 NOTEBOOK PC 笔记本电脑操作系统：Windows 7 家庭普通版 64位 sp1处理器：英特尔 第三代酷睿 i5-3317U 1.70Hz

39、 双核 超低电压处理器内存：4GB主硬盘：日立 HTS545050A7E380 (500GB / 5400转/分)附录2 程序清单算法实现的部分主要Fortran 程序：1. 空间偏导数的计算：U2x=(2.0/(h*h)*(U(I+1,J)-2.0*U(I,J)+U(I-1,J)-(Ux(I+1,J)-Ux(I-1,J)/(2.0*h) U2z=(2.0/(z*z)*(U(I,J+1)-2.0*U(I,J)+U(I,J-1)-(Uz(I,J+1)-Uz(I,J-1)/(2.0*z) U3x=(15.0/(2.0*h*h*h)*(U(I+1,J)-U(I-1,J)-(3.0/(2.0*h*h)

40、*(Ux(I+1,J)+8.0*Ux(I,J)+Ux(I-1,J) U2xz=(1.0/(2.0*h*z)*(-Ux(I+1,J+1)-Ux(I-1,J-1)+Ux(I+1,J)+Ux(I-1,J)-2.0*Ux(I,J+1)+4.0*Ux(I,J)-2.0*U x(I,J-1)+(1.0/(h*h)*(Uz(I+1,J)-2.0*Uz(I,J)+Uz(I-1,J)+(1.0/(4.0*h*h*z)*(5.0*U(I+1,J+1) -5.0*U(I-1,J-1)+U(I+1,J-1)-U(I-1,J+1)-6.0*U(I+1,J)+6.0*U(I-1,J)-4.0*U(I,J+1)+4.0*U(I,J-1) Ux2z=(1.0/(2.0*h*z)*(-Uz(I+1,J+1)-Uz(I-1,J-1)+Uz(I,J+1)+Uz(I,J-1)-2.0*Uz(I+1,J)+4.0*Uz(I,J)-2.0*U z(I-1,J)+(1.0/(z*z)*(Ux(I,J+1)-2.0*Ux(I,J)+Ux(I,J-1)+(1.0/(4.0*h*z*z)*(5.0*U(I+1,J+1) -5.0*U(I-1,J-1)-U(I+1,J-1)+U