1、 一种新的频率选取策略及其在频域波形反演中的应用专 业 名 称： 微电子学与固体电子学摘 要近年来，全波形反演已经成为估计地下介质参数的有力工具，其具体实现方法分为频域和时域两大类。频域波形反演方法相对于时域反演方法而言，数据的选取更加灵活，考虑到地震数据在频域表示中的冗余性，可知前者具有更高的计算效率。频域全波形反演是一种能够获得高分辨率反演结果的反演方法，但在模型较大时进行全频点反演的计算量很大。Sirgue 和 Pratt (2004)提出一种波形反演中的频率点选取的方法，在不影响反演精度的前提下，可以减少频域波形反演所需计算的频率点数目，从而节省计算量。目前，大多数的频率点选取方法都是
2、基于地下介质速度为均匀分布这一假设前提的，为了研究在更加接近实际地下介质的背景模型下的情形，本文在线性递增速度模型的基础上，利用散射波场的 Born 近似、WKBJ 近似以及地震波的一般传播规律，根据模型修正量中竖直方向上的波数覆盖，提出一种新的频率点选取策略。结果表明，在接近实际的地下介质速度分布情况下，可以采用比 Sirgue 和 Pratt 方法更少的频点来进行频域反演，而不会影响反演的精度，从而可以进一步提高反演计算的效率。本文针对一维速度模型、Marmousi 模型和 Overthrust 三个模型进行反演实验，模拟算例结果验证了该优化的频率点选取方法的有效性。关键词：频率点选取 ；
3、波数覆盖 ； 线性递增模型 ； 频域波形反演ABSTRACTIn recent years, waveform inversion has become a powerful tool for estimating the parameters of the underground medium, and its implementation is divided into two major categories of frequency domain and time domain method. Compared to time domain inversion method, the
4、frequency domain inversion method has more flexible data selection, considering the redundancy of seismic data in the frequency domain, show that the frequency domain inversion method has higher computational efficiency.The method of frequency-domain full-waveform inversion is capable of obtaining h
5、igh resolution results, but it suffered a heavy computation when inverting large-scale models with all frequencies. Sirgue and Pratt (2004) proposed a discretization strategy that reduces the input frequencies required for inversion and increases efficiency without loss accuracy of the inversion res
6、ults. At present, the most frequency selection methods are based on the assumption that the velocity of underground medium follows uniform distribution, in order to study the situation that the background is more close to the actual underground medium, in this paper we use the velocity of medium lin
7、early increases with depth, apply the Born approximation of scattered field, the WKBJ approximation, and the general law of seismic wave propagation, according to the coverage of vertical direction wave number in model correction, and proposed a new frequency selection strategy. Experiment results s
8、howed that under the more realistic assumption, our method could obtain the same level of accuracy with fewer frequencies than that required by Sirgue and Pratts strategy. In this paper we conducted experiments with 1D model, Marmousi model and Overthrust model, and the results demonstrated the vali
9、dity of our method.Key words: Frequency selection; Wave-number coverage; Linear increasing velocity; Frequency domain waveform inversion目 录摘 要ABSTRACT1 绪论1.1 研究背景和意义1.2 国内外发展现状1.3 主要研究内容2 频域波形反演算法简介2.1 波形反演的基本步骤2.2 频域正演模拟2.2.1 二维频域波动方程2.2.2 有限差分表示2.2.3 PML 层边界处理方法2.3 频域波形反演3 频率点的选取策略3.1 散射波场的 Born 近
10、似及 Green 函数求解3.2 频域波形反演中的波数覆盖原理3.3 线性递增速度介质中地震波的传播3.4 频率点的选取方法4 实验验证4.1 一维模型合成数据反演4.2 二维模型合成数据反演4.2.1 Marmousi 模型合成数据反演4.2.2 Overthrust 模型合成数据反演5 总结与展望参考文献1绪论随着人们对地下非均匀介质认识的不断深入以及计算机性能的大幅度提升，地震波反演成像技术也相应地得到了非常迅猛地发展。然而，在使用传统频率点选取策略所得频率点进行频率从低到高的迭代频域波形反演的前提下，在勘探的目标区域尺度非常大的情况下，计算效率就会成为此种反演方法所面临的挑战。此时，基
11、于一种新的梯度图像中波数覆盖的频率点选取策略可能会为地震波形反演带来新的发展。本章将回顾相关的研究背景及现状，简述传统的频率点选取策略的原理，阐述它在实际应用中的局限性及其改进的可能性，给出本文的总体研究思路和主要研究内容。1.1研究背景和意义地震勘探(seismic exploration)是地球物理勘探中的一种重要方法，其主要目标是使用人工震源激发的地震波来定位地下介质中的矿藏(如石油，天然气，煤，矿石，矿物，地热能源，水资源等)，获得地下的岩石特性、碳氢化合物属性等地质信息。在陆地地震勘探中，由地表或近地表的人工震源激发的地震波经由地下介质向地球内部传播，当遇到地下介质的物性参数快速变化
12、区域时地震波便会产生散射波(包括反射波，折射波，绕射波等)，一部分返回到设置在地表的高灵敏度检波器的散射波会被此检波器记录下来，所有的这些被不同检波器(hydrophone)记录的地震记录构成地震数据集；同理，在海洋地震勘探中，返回到水面下方沿拖曳电缆方向设置的高灵敏度检波器处的散射波会被此检波器记录下来，所有这些不同检波器的散射波记录构成了地震数据集。对于上述两种勘探情形下所得的地震数据集，其中主要记录了地震波传播的旅行时信息以及振幅响应信息。通过分析地震数据集，人们一方面想建立一种从地震数据集到地下结构概貌的映射，这一过程通常被称作是地震成像或地震偏移成像；另一方面想通过地震数据集来对目标
13、区域中未知的物理参数(如速度，弹性模量，密度，阻抗等)进行适当估算，这一过程通常被称作是地震反演；当然，人们更希望可以通过分析地震数据集一次性地获得目标区域的结构概貌和物性参数。由于地震数据集的复杂性，要想完全正确地实现上述的想法依然是目标我们所面临的困难。地震数据集中包括了所有检波器所能接收的散射波的能量信息，每一个能量信息代表了从震源到检波器的地震波所特有的散射传播方式，根据散射波在地下介质中的具体传播方式可以将其分为单次反射波和多次反射波。根据微扰理论，可以将地震波的单次散射波表示为目标区域物性参数扰动的线性映射，相应地多次反射则是一个非线性映射。考虑到非线性反问题相对于线性问题的复杂性
14、，目前大多数地震成像和地震反演算法都是基于数据中只包含地震波单次散射波的假设。由于地震数据在频域中存在信息冗余(Mulder and Plessix 2004)，因此可以使用所有频率点的一个子集来表示整体信息，考虑到频域中声波方程的有限差分的计算量与频率点数目成正比(Plessix et al. 2001)，利用此结论就可以将时间域中的波形反演转化到频率域中从而提高反演的计算效率。就从低频率点到高频率点的频率迭代反演来说，其计算用时线性正比于波形反演时的频率点数，因此如何高效地去除地震数据集中的冗余频率点，也即如何合理地选择最好(或最少)的频率点组合，就成了提高此时计算效率的关键。现在常用的方
15、法有，(1)基于均匀介质假设的前提下采用固定间隔频率点；(2)基于均匀介质假设的前提下使用间隔随频率值线性递增的频率点；(3)通过速度值随深度线性递增模型中地震波的最大传播深度以及多尺度反演中频率点对局部极点的影响来对方法(2)中的线性系数进行简单的修正。鉴于实际的地下介质往往是非均匀的复杂情形，以及方法(3)只是对方法 (2)的简单修正，因此我们有必要研究一种新的频率点选取策略。本论文旨在充分借鉴现有的频率点选取策略，发现它们的不足之处，提出一种完全基于速度值随深度递增的线性速度模型下的新的频率点选取方法，为下一步的研究打下良好的基础。1.2国内外发展现状人们对于地震数据的处理方法一般分为地
16、震偏移成像和地震反演，其中偏移成像的主要目的是通过恢复光滑背景模型上的不连续区域来构建地下结构概貌，而地震反演的目的则是定量地恢复地下介质的物理参数。如果基于最小二乘法的失配函数(error functional)的成像方法所得的黑塞 (Hessian)矩阵是一个对角(或近对角 )阵，那么此失配函数的最小值可以通过对梯度矩阵前乘黑塞矩阵(或伪(pseudo-)黑塞矩阵) 的逆矩阵来获得，此时的成像方法可以用于表示真振幅偏移成像( Mulder and Plessix 2004; Plessix and Mulder 2004)，考虑到基于最小二乘法的频域波形反演的情形，我们可以发现地震成像和地
17、震反演有时并没有明确的界线( Tarantola 1986)。由于地震数据集在频域表示时存在信息冗余(Mulder and Plessix 2004; Plessix et al. 2001)，因此频域波形反演中可以使用全频率点的一个子集来提高反演的效率。Mulder 和 Plessix (2004)指出在频率点选择时需要考虑两个主要准则： (1)由等价海森堡测不准原理(Heisenbergs principle)所确定的偏移图像中的反射体分辨率；(2)频域奈奎斯特抽样定理(Nyquists theorem)，考虑到偏移包括在给定背景模型上的地震波传播以及多个震源检波器的求和，因此尽管有时不符
18、合奈奎斯特抽样定理但依然可以获得一个没有混叠的结果。在均匀背景的假设下，利用上述的两准则，得到了一个频率点间隔为固定值的频率点选取策略。Sirgue 和 Pratt (2004)根据频域最速下降反演方法，在均匀介质假设的情况下，利用偏移成像和频域波形反演的等价性(Tarantola 1986)，借助平面波近似及波恩近似(Born approximation)原理，揭示了利用多震源检波距地震数据进行频域反演时在垂直方向上波数覆盖存在冗余的现象，据此提出了一种适用此情形下的频率点间隔随频率线性递增的选取策略。同年，Yokota 和 Matsushima (2004)在均匀介质假设基础上，根据观测数
19、据与正演数据之间的相位误差大于 时就会在结果上出现跳周现象(Cycle Skipping，即一个异常高的渡越时间记录)的情形，进而提出一种相对应的频率点选取策略。Mulder 和 Plessix (2008)利用速度值随深度线性递增模型中地震波传播的最大深度以及多尺度反演中频率点对局部极点的影响 (Bunks et al. 1995)，对 Sirgue 和 Pratt (2004)所提出的频率点选取策略中的线性系数进行了简单修正。Anagaw 和 Sacchi (2014)以及 Kim 等 (2011)对 Sirgue 和 Pratt (2004)的方法所得的频率点进行简单地分组，从而验证不同
20、的频率点组合对反演的影响。目前来说，对于频率点选取策略的研究还不是很多，前面所述的方法基本上都是基于均匀背景(或只是对均匀背景简单修正)的情形。随着反演理论的发展以及勘探的需要，人们自然要问是否存在基于相对于均匀背景而言更一般的频率点选取策略。目前，该课题的研究还处于初级阶段，对于从事反演研究的学者来说它是一个具有一定挑战性的研究方向。1.3主要研究内容在认真分析总结现有频率点选取策略的基于基础上，本文结合傅里叶变换的意义，首先建立基于一阶 Born 近似的散射波场积分表达式；然后在 WKBJ近似的情形下获得 Green 函数的解析表达式；在此基础上结合频域波形反演算法通过梯度图像中的波数覆盖
21、情况推导出线性递增速度背景模型下相邻两个频率点之间的关系式，进而得到新的频率点选取策略。本论文各章节的内容安排如下：第一章 介绍了本论文的研究目的及意义，对相关研究现状进行总结，并简要介绍了本论文的研究内容及安排。第二章 回顾了频域波形反演的基于原理和方法，并对预处理梯度法进行了简要介绍。第三章 针对线性递增速度背景模型，提出一种基于散射波场 Born 近似及WKBJ 近似的新的频率点选取策略，并在不同源检距组合情形下，推导了反演中相邻两个频率点之间的关系。第四章 对于本文提出的新的频率点选取策略，在不同的背景速度模型下与Sirgue 和 Pratt (2004)提出的方法进行模型测试。详细对
22、比两种频率点选取策略下反演结果的分辨率及计算用时。第五章 总结本论文的相关研究结果，分析研究中存在的问题，并展望下一步研究的工作重点。2频域波形反演算法简介地震波形反演是一个非线性问题，为了降低问题的求解难度，通常会使用基于梯度的线性处理方法来对其进行求解。常用了求解方法有：全牛顿法，高斯牛顿法，最速下降法，预处理梯度法等。在实际使用中，由于黑塞矩阵的尺寸及其较高的计算消耗，一般很难对其进行直接计算(Mulder and Plessix 2004; Pratt et al. 1998; Tarantola 1984)。通常，最速下降法具有较低的收敛速度(Tarantola 1986)，考虑到黑
23、塞矩阵的求解难度，我们可以通过使用伪黑塞矩阵预处理来加速梯度法的收敛过程(Mulder and Plessix 2004; Plessix and Mulder 2004)。本章以频域反演为前提，主要介绍预处理梯度反演方法，为后续的研究奠定基础。2.1 波形反演的基本步骤目前地震成像中绝大多数的偏移成像方法都是采用线性化简化，其思想为对散射波场进行线性化处理(也即 Born 近似)，然后将 Born 近似的积分形式自然地与广义的 Radon 变换结合起来，进而通过逆广义 Radon 变换来求取近似地反演解，进而恢复地下介质结构。就其中的数学部分而言，Radon 已于 1917 年首次证明了利用
24、投影函数可以恢复出满足条件的唯一的图像函数，因此利用上述的偏移成像方法可以唯一地重建地下介质结构，以达到地震勘探的最终目的。而我们这里要介绍的波形反演方法也遵从上面所述的理论，它们都是通过建立勘探所得的实际数据与用理论计算得到的模拟数据的数据残差，并将此残差与实际地下模型与先验地下模型之间的扰动之间形成映射，通过不断修正先验地下模型来达到最终获得一个在其上的模拟数据与勘探数据拟合度最好的地下模型( Anagaw and Sacchi 2014; Gauthier et al. 1986; Mora 1987; Pratt 1999; Tarantola 1984)。原理虽然相同，但是对于具体的
25、实现方法之间则有较大的区分，为了论文的需要这里仅就波形反演的方法进行讨论，其实现过程基本可分为以下几个部分：(1)在某一地理区域设置勘探的震源及相应的检波器(也即接收点)。(2)采集各个震源所对应的检波器上的波形记录，形成实测(观测)地震记录。由于实测地震记录是波形反演中的表征地下结构的最重要的数据，故实测地震记录的质量显著地影响了波形反演的效果及可信程度。(3)根据实测地震记录中震源及检波器的坐标位置以及检波器中的地震道数据来建立初始介质模型，也即先验模型。在通常的基于梯度的局部最优化反演方法中，如果先验模型与实际勘探区域的实际模型差别过大，则在反演过程中陷入局部极小点的可能性大大增加，在实
26、际的应用过程中一般可以使用走时反演的反演结果作为初始模型，而对于基于全局优化的方法则可以减弱对初始模型的先验限制。(4)在初始介质模型上按实际勘探时的震源及检波器的设置情形进行其上的正演模拟，并记录检波器处的波形信息，即预测(正演)地震记录。(5)把正演地震记录与观测地震记录按每个检波器做差得到地震记录残差，由波恩近似可以看出地震记录残差为先验模型与目标模型之间的模型扰动的映射，由 Radon 变换可知可用此地震记录残差来解得模型的扰动也即恢复地下介质结构。将地震记录残差表示成列矩阵的形式，其中的每个元素代表相应的检波器上的地震记录的残差，把地震残差矩阵按一定方法映射到一个标量目标函数(如 F
27、 范数) 。(6)在目标函数的形式上，使用泛函分析的相关理论获得目标函数极小时所对应的相对于初始介质模型的模型修正量，用此修正量修正初始模型。把修正后的初始模型当作先验介质模型重复步骤(2)-(6)，直到修正后的初始模型上的目标函数满足设定的阈值时停止对介质模型的迭代更新。概括起来波形反演的基本步骤可分为：获取观测地震记录，设置初始介质模型，在初始介质模型上正演获得正演地震记录，对正演地震记录与观测地震记录做差并将结果表示成波场残差矩阵，然后将波场残差矩阵映射成标量目标函数，再通过泛函方法获得介质模型修正量并按此对介质模型进行更新，具体的流程如图 2.1 所示。图 2.1 波形反演流程图2.2
28、 频域正演模拟地震波的速度是指地震波在岩层中的传播速度，简称地震速度，有时又叫做岩石速度。地震速度是地震勘探中最重要的一个参数，也是一个复杂的参数，其不仅与岩石的物理性质(包括：岩性、孔隙度、孔隙中的充填物、岩石的密度、埋藏深度和地质年代等)有关，而且还与介质的结构和求取速度的方法紧密相关。本文中我们选用基于地震速度的函数作为介质模型的参数。2.2.1 二维频域波动方程非均匀介质中的地震波的传播是一个十分复杂的过程，在实际的传播过程中一般可以将地震波分为纵波和横波两种形式，在传播过程中纵波的传播速度大概是横波的二倍左右。考虑到在实际野外勘探记录中的记录时长一般不太长，因此我们完全可以在数据集中
29、滤掉横波然后用只能描述纵波的声波方程来对地震波纵波的传播过程进行模拟，这样可以提高波形正反演的精度。这里我们只考虑地震波在二维恒定密度速度介质中传播的情形，其具体的声波方程可以表示如下：, (2.1)tzxstzpxvtzp,1,22 其中 x 表示水平方向的坐标，z 表示竖直方向的坐标，t 为时间，s(x, z, t)表示在坐标(x, z)处的时域震源，p(x, z, t)表示由震源 s(x, z, t)激发的地震波在坐标(x, z)处的接收点所接收到的时域波形记录。 表示拉普拉斯算子，其具体形式可表2示如下. (2.2)22zx对式(2.1)的两边同时进行时域傅里叶变换，也即把声波方程从时
30、间域变换到频率域，可得, (2.3),22 zxSzxPvzxP其中 为角频率， P(x, z, )为时域接收点波场 p(x, z, t)的傅里叶变换，表示频域接收点波场，S(x, z, )为时域震源 s(x, z, t)的傅里叶变换，表示频域震源。在对式(2.3)所表示的微分方程进行有限差分求解时，可以看到拉普拉斯算子离散时所使用到的所有坐标点构成一个十字线，也即只使用了左、右、上、下四个方向上的离散点，为了在有限差分的具体实现中包含左上、左下、右上、右下方向上的离散点，建立一个两条坐标轴分别为离散网格不同对角线的新的坐标系(仿射坐标系)，其具体情形可由图 2.2 给出。若记原来 XZ 坐标
31、系的 x 轴正方向单位矢量为 i，z 轴正方向单位矢量为 j, 在新 坐标系的 轴正方向单位矢量为 ， 轴正方向单位适量为 。并记ZX i 离散的网格的水平间距为 ，竖直间距为 ，对角线间距为 ，可看出此三者xz满足下式. (2.4)22zx从图 2.2 可以看出新旧两个坐标系轴向上的单位向量之间应满足如下关系式, (2.5)jiizx. (2.6)jz假设一个向量在原坐标系中表示为(x, z)，而在新坐标系中则相应地表示为 ，zx,由于这两个表示均代表同一个向量，故对于同一向量而言其在新旧坐标系下的坐标之间应满足如下关系：, (2.7)zx. (2.8)xz由式(2.7) 和式(2.8) 可
32、知，在原坐标系下的拉普拉斯算子在新坐标系中可以表示为如下形式. (2.9)zxzx zxxzx 2222222 14141 由式(2.9)可以看出拉普拉斯算子在新坐标系中的形式相对于式(2.2)的形式而言，出现了混合偏导形式，为了使之具有与式(2.2)一样的简洁的形式，我们在此人为地令网格的水平间距和竖直间距相同，并记, (2.10)zx则式(2.9)可以简洁地表示为如下形式. (2.11)22zx由式(2.11) 可以看出其与式(2.2) 具有非常类似的表达形式，而相应地如果用新坐标来表示声波方程在地下速度介质中的传播规律则可将声波方程表示为如下形式：, (2.12),22 zxSzxPvz
33、xP其中 为拉普拉斯算子在新坐标系下的表示，也即. (2.13)22zx由前面讨论可知，式(2.3)和式(2.12) 分别为频域声波方程在原坐标系及新坐标系下的表达式，若 和 表示同一个接收点的波形记录，也即,P,zx坐标(x, z)和坐标 满足式(2.7)和式(2.8) ，那么我们可以将式 (2.3)和式(2.12)zx通过加权的方式将其合并为下式. (2.14) ,1, 222 zxSzxPvzzxP 其中表示在合成的表达式中新旧坐标系下声波方程所占的比例，在实际应用中，可以通过最优化方法来求得此参量的使用值。在后文中可以看到，这时在声波方程的有限差分中不再只是使用十字线上的点同时使用离散
34、网格对角线上的点(而相应的不同位置的点所点的权重则由 反映 )的差分式，并且当 时则恢复到式 所表示的情形。图 2.2 两种坐标系的表示2.2.2 有限差分表示有限差分法是使用差分原理来解决微分方程的一种常用的数值求解方法，其中心思想为：将需要进行数值计算的区域离散成很多的点的集合，离散方法通常为离散为正方形或矩形网格(此时后续的求解在形式上较为简洁)，然后对方程进行微分或偏微分向差商形式的代换，在最终的差分表达形式相对于原来的微分或偏微分方程而言是相容及收敛的，并且差分形式本身是稳定的时候那么就可以使用此差分形式来对原微分或偏微分方程进行数值求解。对于频域的声波方程而言通常有多种离散方法(G
35、u et al. 2013; Jo et al. 1996; Liu 2014)，本文中我们选用一种新型的 9 点离散策略。首先把速度模型按水平及竖直方向间隔均为 的正方形网格进行离散，并把此对角线的长度记为 (明显有 )，离散后水平方向的总的网格点数2目记为 ，相应的竖直方向的总的网格点数目记为 。为了用有限差分法求xNzN解式(2.3)、式 (2.12)及式(2.14)所表示的声波方程，需要使用泰勒级数来完成微分形式到差分形式的转换，对于一元函数 由高数知识可知若函数在 的某xf 0x个领域中具有直到(n+1) 阶的导数，则在该领域内 的 n 阶泰勒公式为：f. xRfnxxfxfxfx
36、nn 0002000 !1(2.15)其中 为拉格朗日余项，在不影响计算或影响很小的情况下可以将其直Rn接计为 0。为了求解的方便，我们把 简记为 ，同理把,mPnmP,简记为 。把波场 及波场 按式(2.15)的形式表示为关,mPnm, n,1n,1于坐标点( m, n)的函数，可得：, (2.16)nRxxxxPnmnmnmn 1!4!3!24,1 . (2.17)PP nnnmnn ,! 4,3,2,1 把式(2.16)和式 (2.17)等号两端分别相加，再经过简单的化简可得：, (2.18)nmRnxPPx nmnmnmn 1,1124424,2,1,2, 同理，可得：. (2.19)
37、annzz nmnmnmn 1,1,124424,2,1,2, 把式(2.18)与式 (2.19)在等号两边分别求和可得：. (2.20) 1,1,1,12444424,42,1,1,2 nmRnmRn zPxPPnnnmnnmn同理可求得在新坐标系下拉普拉斯算子的相应的表达为：. (2.21) 1,1,1,12444424,42,1,1,1,1,2 nmRnmRn zPxPPP nmnnmnnnmn由式(2.7)及式 (2.8)，可将式 (2.21)中的在新坐标系下表示的四阶偏导转化为在原来坐标系下的情形，可得：. (2.22) 1,1,1,1222444442 2,4,4,4 ,1,1,1
38、,1,2 nmRnmRnzxPzPxPnnnmmnmn如果 和 表示同一个接收点的数据，则可利用式(2.20)和式(2.22)把在nmP,混合新旧两种坐标系的波动方程式(2.14)离散为如下所示形式：, (2.33) RSPvPnmnmnnmnmnn ,2,2,1,1,1,12,1 4其中的 R 为式(2.14) 离散时的截断误差，其具体形式可根据前面的表达表示如下：. (2.34) 1,1,1,11 , 12221 44442 2,44,4,4,4,2 nmRnmRnR zxPzPxzPx nmnmnmnnm由图 2.3 中离散点的分布情形，以及离散时的离散网格为正方形，则式(2.33)及式
39、(2.34)可化为, (2.35) RSPvPPP nmnmnnmnnmnnm ,2,2,1,1,1,1,2,1 44. (2.36) 1,1,1,1,212 4444 42,44,4, nmRnRnmRnRzxzxnnn由式(2.36) 可知，当 时，R 的值也趋于零，故可忽略式(2.35)中的 R 进0而将式(2.35) 表示为：. (2.37) nmnmnnmnnmnnm SPvPPP ,2,2,1,1,1,1,2,1 44 为了减少声波方程离散化的频散，进一步提升正演的精度，可对式(2.37)中所表示的质量加速项进行优化(Jo et al. 1996)，具体的优化可以有多种，nmPv,
40、2,这里选用 上、下、左、右及其本身五个位置上的波场值在一阶近似的情况,下来近似地表示 ，也即：nm, (2.38)1,1, 41nmnmPP其中的参量 表示在近似时中心点所占的比例，在实际应用中，可以通过最优化方法来求得此参量的使用值。把式(2.38)代入式(2.37)可得：, (2.39)nmnmnmnmn nnnnnmmnm SPPvP ,1,1,2, 2,1,1,1,1,2,1444 显然上式相对于式(2.14) 而言是相容且收敛的，当离散网格间距、模型的速度值以及相应的频率点之间满足一定关系时(简单情形下一个波长覆盖 4 个或者更多的网格点时)，式(2.39) 所表示的离散方程本身是
41、稳定的(Boonyasiriwat et al. 2009; Chen 2014; Chen and Cao 2014; Huang and Schuster 2014; Jo et al. 1996; 刘璐等 2013)，故可以将式(2.39)作为声波方程的有限差分形式来对其进行数值求解。对于离散后的每一个网格点而言，均存在一个如式(2.39)所示的线性方程，就整体而言可形成一个由 个方程组成的线性方程组。这样在速度模型上zxN的正演模拟过程可以表示为如下的线性方程组形式：, (2.40)SAP其中 为表示模型空间的压强场的列向量，速度模型中沿 X 轴正方向变动过程中沿 Z 轴正方向从上到下
42、的网格点分别与 中从上到下的元素一一对应； 为PA离散尺度、频率点和差分格式等有关的复数方阵，一般情况下，矩阵 是高度稀疏、非对称及非正定的； 为震源向量。S通过对(2.40) 进行求解可以获得此时的压强场的列向量 ，也即正演结果，P而求解方法一般分为两类：迭代法和直接求解法(Pratt 1999)，本文中则通过调用 mumps 库来进行直接求解。对于解出的结果 ，由于我们只对设置了检波器的位置上的波场值感兴趣，故我们需要对列矩阵 进行处理，即保留检波器位P置所对应的元素值，而把其余的所未设置检波器位置的网格点所对应的元素值设置为 0，经此处理后的列矩阵 记录了给定频率 点上的波场记录，对于所
43、有P频率点的波场数据进行反傅里叶变换即可获得给定先验速度模型上的正演波场记录。图 2.3 在不同坐标系中的离散网格点2.2.3 PML层边界处理方法对于真实的地下介质模型而言，一般而言是一个无限区域，但是由于使用计算机计算的有限差分法只能对有限的计算的计算量进行计算，因此在离散时只能考虑我们感兴趣的区域，这样由于人为边界所产生边界反射就不可避免有限差分所得的结果与与真实情况之间出现差异，就离散后的模型而言越靠近边界区域这种差异表现的越明显。为了尽量地消除边界反射，在这里我们采用一种常用的也比较有代表性的方法完全匹配层边界吸收法，这种方法最开始主要用于电磁波传播规律的有限差分方程的边界处理中，后
44、来逐渐被用在地震层析成像中，相对于其它的边界处理方法而言，PML 方法(Berenger 1994; Hustedt et al. 2003; Liu et al. 2011; 郭念民 和 吴国忱 2012; 刘璐等 2013; 张毅 2007)是目前来说最有效并且相对简洁的方法，这一点在相关的有很多实验中已经得到了验证。由于 PML 层是对有限差分法离散后的区域的边界之外的区域进行处理，因此可知此区域中并没有任何的震源，故此时的时域地震波方程相对于式(2.1)而言只是使后者的震源变为零值，也即：, (2.41)tzxpvtzxp22,1,为了使上式降阶，也即由两个一阶偏微分方程来表示 PML
45、 层中地震波传播规律，在这里我们引入两个中间变量 和 ，中间变量中的下标 x 和 z 只tzxq,tzx,是一个记号标记并不指代任何其它的方向信息，且这两个中间变量满足以下方程：, (2.42)(,tzxptzxq; (2.43)z对于时域的波场响应 而言，我们在这里将其拆分成 和),(tzx tzxp,两部分，相应分量的下标 x 和 z 只是一个记号标记并不指代任何其它tzxp,的方向信息，也即：, (2.44)tptzxtz,由数学关系可知一定存在上述的波场响应的拆分方式，使得下两式成立，也即：, (2.45)xtzqzxvtzxp),(,2. (2.46)zttz ,),(,2式(2.4
46、5)及式 (2.46)即为我们用来表示 PML 层中声波方程的两个一阶偏微分方程，由形式上我们可以看出上两式所描述的地震波在传播过程中并没有能量的衰减，也即能量是守恒的，这一点使得相对于没有边界处理时的情形而言只是相当于扩大离散区域的大小，当然如果我们能把区域扩展到无穷大那么在物理上来说也是对真实问题的较好的近似，但是前面我们也提过这一点对于实际计算而言是没有任何意义的，因此式(2.45)和式(2.46) 的边界处理对于边界处的地震波反射而言并没有得到根本地解决。因此为了消除边界处的波形反射，PML 区域中传播的地震波必须在能量上出现一定程度上的衰减，根据一阶常微分方程的解的一般规律，在这里我们对式(2.42)、式(2.43) 、式(2.45)和式(2.46) 所表示的PML 区域中的地震波传播的规律中加入衰减项(刘璐等 2013)，也即：, (2.47)(, tzxptzxqtzxqx , (2.48),),(, tzttzzz,
