一种新的频率选取策略及其在频域波形反演中的应用

一种新的频率选取策略及其在频域波形反演中的应用 专 业 名 称： 微电子学与固体电子学 摘 要近年来，全波形反演已经成为估计地下介质参数的有力工具，其具体实现方法分为频域和时域两大类。频域波形反演方法相对于时域反演方法而言，数据的选取更加灵活，考虑到地震数据在频域表示中的冗余性，可知前者具有更高的计算效率。频域全波形反演是一种能够获得高分辨率反演结果的反演方法，但在模型较大时进行全频点反演的计算量很大。Sirgue和Pratt (2004)提出一种波形反演中的频率点选取的方法，在不影响反演精度的前提下，可以减少频域波形反演所需计算的频率点数目，从而节省计算量。目前，大多数的频率点选取方法都是基于地下介质速度为均匀分布这一假设前提的，为了研究在更加接近实际地下介质的背景模型下的情形，本文在线性递增速度模型的基础上，利用散射波场的Born近似、WKBJ近似以及地震波的一般传播规律，根据模型修正量中竖直方向上的波数覆盖，提出一种新的频率点选取策略。结果表明，在接近实际的地下介质速度分布情况下，可以采用比Sirgue和Pratt方法更少的频点来进行频域反演，而不会影响反演的精度，从而可以进一步提高反演计算的效率。本文针对一维速度模型、Marmousi模型和Overthrust三个模型进行反演实验，模拟算例结果验证了该优化的频率点选取方法的有效性。关键词：频率点选取 ；波数覆盖 ； 线性递增模型 ； 频域波形反演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 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 high 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 results. 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 linearly 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 showed 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 validity 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近似及Green函数求解3.2频域波形反演中的波数覆盖原理3.3线性递增速度介质中地震波的传播3.4频率点的选取方法4实验验证4.1一维模型合成数据反演4.2二维模型合成数据反演4.2.1 Marmousi模型合成数据反演4.2.2 Overthrust模型合成数据反演5总结与展望参考文献1绪论随着人们对地下非均匀介质认识的不断深入以及计算机性能的大幅度提升，地震波反演成像技术也相应地得到了非常迅猛地发展。然而，在使用传统频率点选取策略所得频率点进行频率从低到高的迭代频域波形反演的前提下，在勘探的目标区域尺度非常大的情况下，计算效率就会成为此种反演方法所面临的挑战。此时，基于一种新的梯度图像中波数覆盖的频率点选取策略可能会为地震波形反演带来新的发展。本章将回顾相关的研究背景及现状，简述传统的频率点选取策略的原理，阐述它在实际应用中的局限性及其改进的可能性，给出本文的总体研究思路和主要研究内容。1.1研究背景和意义地震勘探(seismic exploration)是地球物理勘探中的一种重要方法，其主要目标是使用人工震源激发的地震波来定位地下介质中的矿藏(如石油，天然气，煤，矿石，矿物，地热能源，水资源等)，获得地下的岩石特性、碳氢化合物属性等地质信息。在陆地地震勘探中，由地表或近地表的人工震源激发的地震波经由地下介质向地球内部传播，当遇到地下介质的物性参数快速变化区域时地震波便会产生散射波(包括反射波，折射波，绕射波等)，一部分返回到设置在地表的高灵敏度检波器的散射波会被此检波器记录下来，所有的这些被不同检波器(hydrophone)记录的地震记录构成地震数据集；同理，在海洋地震勘探中，返回到水面下方沿拖曳电缆方向设置的高灵敏度检波器处的散射波会被此检波器记录下来，所有这些不同检波器的散射波记录构成了地震数据集。对于上述两种勘探情形下所得的地震数据集，其中主要记录了地震波传播的旅行时信息以及振幅响应信息。通过分析地震数据集，人们一方面想建立一种从地震数据集到地下结构概貌的映射，这一过程通常被称作是地震成像或地震偏移成像；另一方面想通过地震数据集来对目标区域中未知的物理参数(如速度，弹性模量，密度，阻抗等)进行适当估算，这一过程通常被称作是地震反演；当然，人们更希望可以通过分析地震数据集一次性地获得目标区域的结构概貌和物性参数。由于地震数据集的复杂性，要想完全正确地实现上述的想法依然是目标我们所面临的困难。地震数据集中包括了所有检波器所能接收的散射波的能量信息，每一个能量信息代表了从震源到检波器的地震波所特有的散射传播方式，根据散射波在地下介质中的具体传播方式可以将其分为单次反射波和多次反射波。根据微扰理论，可以将地震波的单次散射波表示为目标区域物性参数扰动的线性映射，相应地多次反射则是一个非线性映射。考虑到非线性反问题相对于线性问题的复杂性，目前大多数地震成像和地震反演算法都是基于数据中只包含地震波单次散射波的假设。由于地震数据在频域中存在信息冗余(Mulder and Plessix 2004)，因此可以使用所有频率点的一个子集来表示整体信息，考虑到频域中声波方程的有限差分的计算量与频率点数目成正比(Plessix et al. 2001)，利用此结论就可以将时间域中的波形反演转化到频率域中从而提高反演的计算效率。就从低频率点到高频率点的频率迭代反演来说，其计算用时线性正比于波形反演时的频率点数，因此如何高效地去除地震数据集中的冗余频率点，也即如何合理地选择最好(或最少)的频率点组合，就成了提高此时计算效率的关键。现在常用的方法有，(1)基于均匀介质假设的前提下采用固定间隔频率点；(2)基于均匀介质假设的前提下使用间隔随频率值线性递增的频率点；(3)通过速度值随深度线性递增模型中地震波的最大传播深度以及多尺度反演中频率点对局部极点的影响来对方法(2)中的线性系数进行简单的修正。鉴于实际的地下介质往往是非均匀的复杂情形，以及方法(3)只是对方法(2)的简单修正，因此我们有必要研究一种新的频率点选取策略。本论文旨在充分借鉴现有的频率点选取策略，发现它们的不足之处，提出一种完全基于速度值随深度递增的线性速度模型下的新的频率点选取方法，为下一步的研究打下良好的基础。1.2国内外发展现状人们对于地震数据的处理方法一般分为地震偏移成像和地震反演，其中偏移成像的主要目的是通过恢复光滑背景模型上的不连续区域来构建地下结构概貌，而地震反演的目的则是定量地恢复地下介质的物理参数。如果基于最小二乘法的失配函数(error functional)的成像方法所得的黑塞(Hessian)矩阵是一个对角(或近对角)阵，那么此失配函数的最小值可以通过对梯度矩阵前乘黑塞矩阵(或伪(pseudo-)黑塞矩阵)的逆矩阵来获得，此时的成像方法可以用于表示真振幅偏移成像(Mulder and Plessix 2004; Plessix and Mulder 2004)，考虑到基于最小二乘法的频域波形反演的情形，我们可以发现地震成像和地震反演有时并没有明确的界线(Tarantola 1986)。由于地震数据集在频域表示时存在信息冗余(Mulder and Plessix 2004; Plessix et al. 2001)，因此频域波形反演中可以使用全频率点的一个子集来提高反演的效率。Mulder和Plessix (2004)指出在频率点选择时需要考虑两个主要准则：(1)由等价海森堡测不准原理(Heisenbergs principle)所确定的偏移图像中的反射体分辨率；(2)频域奈奎斯特抽样定理(Nyquists theorem)，考虑到偏移包括在给定背景模型上的地震波传播以及多个震源检波器的求和，因此尽管有时不符合奈奎斯特抽样定理但依然可以获得一个没有混叠的结果。在均匀背景的假设下，利用上述的两准则，得到了一个频率点间隔为固定值的频率点选取策略。Sirgue和Pratt (2004)根据频域最速下降反演方法，在均匀介质假设的情况下，利用偏移成像和频域波形反演的等价性(Tarantola 1986)，借助平面波近似及波恩近似(Born approximation)原理，揭示了利用多震源检波距地震数据进行频域反演时在垂直方向上波数覆盖存在冗余的现象，据此提出了一种适用此情形下的频率点间隔随频率线性递增的选取策略。同年，Yokota和Matsushima (2004)在均匀介质假设基础上，根据观测数据与正演数据之间的相位误差大于时就会在结果上出现跳周现象(Cycle Skipping，即一个异常高的渡越时间记录)的情形，进而提出一种相对应的频率点选取策略。Mulder和Plessix (2008)利用速度值随深度线性递增模型中地震波传播的最大深度以及多尺度反演中频率点对局部极点的影响(Bunks et al. 1995)，对Sirgue和Pratt (2004)所提出的频率点选取策略中的线性系数进行了简单修正。Anagaw和Sacchi (2014)以及Kim等 (2011)对Sirgue和Pratt (2004)的方法所得的频率点进行简单地分组，从而验证不同的频率点组合对反演的影响。目前来说，对于频率点选取策略的研究还不是很多，前面所述的方法基本上都是基于均匀背景(或只是对均匀背景简单修正)的情形。随着反演理论的发展以及勘探的需要，人们自然要问是否存在基于相对于均匀背景而言更一般的频率点选取策略。目前，该课题的研究还处于初级阶段，对于从事反演研究的学者来说它是一个具有一定挑战性的研究方向。1.3主要研究内容在认真分析总结现有频率点选取策略的基于基础上，本文结合傅里叶变换的意义，首先建立基于一阶Born近似的散射波场积分表达式；然后在WKBJ近似的情形下获得Green函数的解析表达式；在此基础上结合频域波形反演算法通过梯度图像中的波数覆盖情况推导出线性递增速度背景模型下相邻两个频率点之间的关系式，进而得到新的频率点选取策略。本论文各章节的内容安排如下：第一章 介绍了本论文的研究目的及意义，对相关研究现状进行总结，并简要介绍了本论文的研究内容及安排。第二章 回顾了频域波形反演的基于原理和方法，并对预处理梯度法进行了简要介绍。第三章 针对线性递增速度背景模型，提出一种基于散射波场Born近似及WKBJ近似的新的频率点选取策略，并在不同源检距组合情形下，推导了反演中相邻两个频率点之间的关系。第四章 对于本文提出的新的频率点选取策略，在不同的背景速度模型下与Sirgue和Pratt (2004)提出的方法进行模型测试。详细对比两种频率点选取策略下反演结果的分辨率及计算用时。第五章 总结本论文的相关研究结果，分析研究中存在的问题，并展望下一步研究的工作重点。2频域波形反演算法简介地震波形反演是一个非线性问题，为了降低问题的求解难度，通常会使用基于梯度的线性处理方法来对其进行求解。常用了求解方法有：全牛顿法，高斯牛顿法，最速下降法，预处理梯度法等。在实际使用中，由于黑塞矩阵的尺寸及其较高的计算消耗，一般很难对其进行直接计算(Mulder and Plessix 2004; Pratt et al. 1998; Tarantola 1984)。通常，最速下降法具有较低的收敛速度(Tarantola 1986)，考虑到黑塞矩阵的求解难度，我们可以通过使用伪黑塞矩阵预处理来加速梯度法的收敛过程(Mulder and Plessix 2004; Plessix and Mulder 2004)。本章以频域反演为前提，主要介绍预处理梯度反演方法，为后续的研究奠定基础。2.1 波形反演的基本步骤目前地震成像中绝大多数的偏移成像方法都是采用线性化简化，其思想为对散射波场进行线性化处理(也即Born近似)，然后将Born近似的积分形式自然地与广义的Radon变换结合起来，进而通过逆广义Radon变换来求取近似地反演解，进而恢复地下介质结构。就其中的数学部分而言，Radon已于1917年首次证明了利用投影函数可以恢复出满足条件的唯一的图像函数，因此利用上述的偏移成像方法可以唯一地重建地下介质结构，以达到地震勘探的最终目的。而我们这里要介绍的波形反演方法也遵从上面所述的理论，它们都是通过建立勘探所得的实际数据与用理论计算得到的模拟数据的数据残差，并将此残差与实际地下模型与先验地下模型之间的扰动之间形成映射，通过不断修正先验地下模型来达到最终获得一个在其上的模拟数据与勘探数据拟合度最好的地下模型(Anagaw and Sacchi 2014; Gauthier et al. 1986; Mora 1987; Pratt 1999; Tarantola 1984)。原理虽然相同，但是对于具体的实现方法之间则有较大的区分，为了论文的需要这里仅就波形反演的方法进行讨论，其实现过程基本可分为以下几个部分：(1)在某一地理区域设置勘探的震源及相应的检波器(也即接收点)。(2)采集各个震源所对应的检波器上的波形记录，形成实测(观测)地震记录。由于实测地震记录是波形反演中的表征地下结构的最重要的数据，故实测地震记录的质量显著地影响了波形反演的效果及可信程度。(3)根据实测地震记录中震源及检波器的坐标位置以及检波器中的地震道数据来建立初始介质模型，也即先验模型。在通常的基于梯度的局部最优化反演方法中，如果先验模型与实际勘探区域的实际模型差别过大，则在反演过程中陷入局部极小点的可能性大大增加，在实际的应用过程中一般可以使用走时反演的反演结果作为初始模型，而对于基于全局优化的方法则可以减弱对初始模型的先验限制。(4)在初始介质模型上按实际勘探时的震源及检波器的设置情形进行其上的正演模拟，并记录检波器处的波形信息，即预测(正演)地震记录。(5)把正演地震记录与观测地震记录按每个检波器做差得到地震记录残差，由波恩近似可以看出地震记录残差为先验模型与目标模型之间的模型扰动的映射，由Radon变换可知可用此地震记录残差来解得模型的扰动也即恢复地下介质结构。将地震记录残差表示成列矩阵的形式，其中的每个元素代表相应的检波器上的地震记录的残差，把地震残差矩阵按一定方法映射到一个标量目标函数(如F范数)。(6)在目标函数的形式上，使用泛函分析的相关理论获得目标函数极小时所对应的相对于初始介质模型的模型修正量，用此修正量修正初始模型。把修正后的初始模型当作先验介质模型重复步骤(2)-(6)，直到修正后的初始模型上的目标函数满足设定的阈值时停止对介质模型的迭代更新。概括起来波形反演的基本步骤可分为：获取观测地震记录，设置初始介质模型，在初始介质模型上正演获得正演地震记录，对正演地震记录与观测地震记录做差并将结果表示成波场残差矩阵，然后将波场残差矩阵映射成标量目标函数，再通过泛函方法获得介质模型修正量并按此对介质模型进行更新，具体的流程如图2.1所示。图2.1 波形反演流程图2.2 频域正演模拟地震波的速度是指地震波在岩层中的传播速度，简称地震速度，有时又叫做岩石速度。地震速度是地震勘探中最重要的一个参数，也是一个复杂的参数，其不仅与岩石的物理性质(包括：岩性、孔隙度、孔隙中的充填物、岩石的密度、埋藏深度和地质年代等)有关，而且还与介质的结构和求取速度的方法紧密相关。本文中我们选用基于地震速度的函数作为介质模型的参数。2.2.1 二维频域波动方程非均匀介质中的地震波的传播是一个十分复杂的过程，在实际的传播过程中一般可以将地震波分为纵波和横波两种形式，在传播过程中纵波的传播速度大概是横波的二倍左右。考虑到在实际野外勘探记录中的记录时长一般不太长，因此我们完全可以在数据集中滤掉横波然后用只能描述纵波的声波方程来对地震波纵波的传播过程进行模拟，这样可以提高波形正反演的精度。这里我们只考虑地震波在二维恒定密度速度介质中传播的情形，其具体的声波方程可以表示如下：, (2.1)其中x表示水平方向的坐标，z表示竖直方向的坐标，t为时间，s(x, z, t)表示在坐标(x, z)处的时域震源，p(x, z, t)表示由震源s(x, z, t)激发的地震波在坐标(x, z)处的接收点所接收到的时域波形记录。表示拉普拉斯算子，其具体形式可表示如下. (2.2)对式(2.1)的两边同时进行时域傅里叶变换，也即把声波方程从时间域变换到频率域，可得, (2.3)其中w为角频率，P(x, z, w)为时域接收点波场p(x, z, t)的傅里叶变换，表示频域接收点波场，S(x, z, w)为时域震源s(x, z, t)的傅里叶变换，表示频域震源。在对式(2.3)所表示的微分方程进行有限差分求解时，可以看到拉普拉斯算子离散时所使用到的所有坐标点构成一个十字线，也即只使用了左、右、上、下四个方向上的离散点，为了在有限差分的具体实现中包含左上、左下、右上、右下方向上的离散点，建立一个两条坐标轴分别为离散网格不同对角线的新的坐标系(仿射坐标系)，其具体情形可由图2.2给出。若记原来XZ坐标系的x轴正方向单位矢量为i，z轴正方向单位矢量为j, 在新坐标系的轴正方向单位矢量为，轴正方向单位适量为。并记离散的网格的水平间距为，竖直间距为，对角线间距为，可看出此三者满足下式. (2.4)从图2.2可以看出新旧两个坐标系轴向上的单位向量之间应满足如下关系式, (2.5). (2.6)假设一个向量在原坐标系中表示为(x, z)，而在新坐标系中则相应地表示为，由于这两个表示均代表同一个向量，故对于同一向量而言其在新旧坐标系下的坐标之间应满足如下关系：, (2.7). (2.8)由式(2.7)和式(2.8)可知，在原坐标系下的拉普拉斯算子在新坐标系中可以表示为如下形式. (2.9)由式(2.9)可以看出拉普拉斯算子在新坐标系中的形式相对于式(2.2)的形式而言，出现了混合偏导形式，为了使之具有与式(2.2)一样的简洁的形式，我们在此人为地令网格的水平间距和竖直间距相同，并记, (2.10)则式(2.9)可以简洁地表示为如下形式. (2.11)由式(2.11)可以看出其与式(2.2)具有非常类似的表达形式，而相应地如果用新坐标来表示声波方程在地下速度介质中的传播规律则可将声波方程表示为如下形式：, (2.12)其中为拉普拉斯算子在新坐标系下的表示，也即. (2.13)由前面讨论可知，式(2.3)和式(2.12)分别为频域声波方程在原坐标系及新坐标系下的表达式，若和表示同一个接收点的波形记录，也即坐标(x, z)和坐标满足式(2.7)和式(2.8)，那么我们可以将式(2.3)和式(2.12)通过加权的方式将其合并为下式. (2.14)其中a表示在合成的表达式中新旧坐标系下声波方程所占的比例，在实际应用中，可以通过最优化方法来求得此参量的使用值。在后文中可以看到，这时在声波方程的有限差分中不再只是使用十字线上的点同时使用离散网格对角线上的点(而相应的不同位置的点所点的权重则由a反映)的差分式，并且当a=1时则恢复到式(2.3)所表示的情形。图2.2 两种坐标系的表示2.2.2 有限差分表示有限差分法是使用差分原理来解决微分方程的一种常用的数值求解方法，其中心思想为：将需要进行数值计算的区域离散成很多的点的集合，离散方法通常为离散为正方形或矩形网格(此时后续的求解在形式上较为简洁)，然后对方程进行微分或偏微分向差商形式的代换，在最终的差分表达形式相对于原来的微分或偏微分方程而言是相容及收敛的，并且差分形式本身是稳定的时候那么就可以使用此差分形式来对原微分或偏微分方程进行数值求解。对于频域的声波方程而言通常有多种离散方法(Gu et al. 2013; Jo et al. 1996; Liu 2014)，本文中我们选用一种新型的9点离散策略。首先把速度模型按水平及竖直方向间隔均为的正方形网格进行离散，并把此对角线的长度记为(明显有)，离散后水平方向的总的网格点数目记为，相应的竖直方向的总的网格点数目记为。为了用有限差分法求解式(2.3)、式(2.12)及式(2.14)所表示的声波方程，需要使用泰勒级数来完成微分形式到差分形式的转换，对于一元函数由高数知识可知若函数在的某个领域中具有直到(n+1)阶的导数，则在该领域内的n阶泰勒公式为：. (2.15)其中为拉格朗日余项，在不影响计算或影响很小的情况下可以将其直接计为0。为了求解的方便，我们把简记为，同理把简记为。把波场及波场按式(2.15)的形式表示为关于坐标点(m, n)的函数，可得：, (2.16). (2.17)把式(2.16)和式(2.17)等号两端分别相加，再经过简单的化简可得：, (2.18)同理，可得：. (2.19)把式(2.18)与式(2.19)在等号两边分别求和可得：. (2.20)同理可求得在新坐标系下拉普拉斯算子的相应的表达为：. (2.21)由式(2.7)及式(2.8)，可将式(2.21)中的在新坐标系下表示的四阶偏导转化为在原来坐标系下的情形，可得：. (2.22)如果和表示同一个接收点的数据，则可利用式(2.20)和式(2.22)把在混合新旧两种坐标系的波动方程式(2.14)离散为如下所示形式：, (2.33)其中的R为式(2.14)离散时的截断误差，其具体形式可根据前面的表达表示如下：. (2.34)由图2.3中离散点的分布情形，以及离散时的离散网格为正方形，则式(2.33)及式(2.34)可化为, (2.35). (2.36)由式(2.36)可知，当时，R的值也趋于零，故可忽略式(2.35)中的R进而将式(2.35)表示为：. (2.37)为了减少声波方程离散化的频散，进一步提升正演的精度，可对式(2.37)中所表示的质量加速项进行优化(Jo et al. 1996)，具体的优化可以有多种，这里选用上、下、左、右及其本身五个位置上的波场值在一阶近似的情况下来近似地表示，也即：, (2.38)其中的参量b表示在近似时中心点所占的比例，在实际应用中，可以通过最优化方法来求得此参量的使用值。把式(2.38)代入式(2.37)可得：, (2.39)显然上式相对于式(2.14)而言是相容且收敛的，当离散网格间距、模型的速度值以及相应的频率点之间满足一定关系时(简单情形下一个波长覆盖4个或者更多的网格点时)，式(2.39)所表示的离散方程本身是稳定的(Boonyasiriwat et al. 2009; Chen 2014; Chen and Cao 2014; Huang and Schuster 2014; Jo et al. 1996; 刘璐等 2013)，故可以将式(2.39)作为声波方程的有限差分形式来对其进行数值求解。对于离散后的每一个网格点而言，均存在一个如式(2.39)所示的线性方程，就整体而言可形成一个由个方程组成的线性方程组。这样在速度模型上的正演模拟过程可以表示为如下的线性方程组形式：, (2.40)其中为表示模型空间的压强场的列向量，速度模型中沿X轴正方向变动过程中沿Z轴正方向从上到下的网格点分别与中从上到下的元素一一对应；为离散尺度、频率点和差分格式等有关的复数方阵，一般情况下，矩阵是高度稀疏、非对称及非正定的；为震源向量。通过对(2.40)进行求解可以获得此时的压强场的列向量，也即正演结果，而求解方法一般分为两类：迭代法和直接求解法(Pratt 1999)，本文中则通过调用mumps库来进行直接求解。对于解出的结果，由于我们只对设置了检波器的位置上的波场值感兴趣，故我们需要对列矩阵进行处理，即保留检波器位置所对应的元素值，而把其余的所未设置检波器位置的网格点所对应的元素值设置为0，经此处理后的列矩阵记录了给定频率w点上的波场记录，对于所有频率点的波场数据进行反傅里叶变换即可获得给定先验速度模型上的正演波场记录。图2.3 在不同坐标系中的离散网格点2.2.3 PML层边界处理方法对于真实的地下介质模型而言，一般而言是一个无限区域，但是由于使用计算机计算的有限差分法只能对有限的计算的计算量进行计算，因此在离散时只能考虑我们感兴趣的区域，这样由于人为边界所产生边界反射就不可避免有限差分所得的结果与与真实情况之间出现差异，就离散后的模型而言越靠近边界区域这种差异表现的越明显。为了尽量地消除边界反射，在这里我们采用一种常用的也比较有代表性的方法完全匹配层边界吸收法，这种方法最开始主要用于电磁波传播规律的有限差分方程的边界处理中，后来逐渐被用在地震层析成像中，相对于其它的边界处理方法而言，PML方法(Berenger 1994; Hustedt et al. 2003; Liu et al. 2011; 郭念民 和 吴国忱 2012; 刘璐等 2013; 张毅 2007)是目前来说最有效并且相对简洁的方法，这一点在相关的有很多实验中已经得到了验证。由于PML层是对有限差分法离散后的区域的边界之外的区域进行处理，因此可知此区域中并没有任何的震源，故此时的时域地震波方程相对于式(2.1)而言只是使后者的震源变为零值，也即：, (2.41)为了使上式降阶，也即由两个一阶偏微分方程来表示PML层中地震波传播规律，在这里我们引入两个中间变量和，中间变量中的下标x和z只是一个记号标记并不指代任何其它的方向信息，且这两个中间变量满足以下方程：, (2.42); (2.43)对于时域的波场响应而言，我们在这里将其拆分成和两部分，相应分量的下标x和z只是一个记号标记并不指代任何其它的方向信息，也即：, (2.44)由数学关系可知一定存在上述的波场响应的拆分方式，使得下两式成立，也即：, (2.45). (2.46)式(2.45)及式(2.46)即为我们用来表示PML层中声波方程的两个一阶偏微分方程，由形式上我们可以看出上两式所描述的地震波在传播过程中并没有能量的衰减，也即能量是守恒的，这一点使得相对于没有边界处理时的情形而言只是相当于扩大离散区域的大小，当然如果我们能把区域扩展到无穷大那么在物理上来说也是对真实问题的较好的近似，但是前面我们也提过这一点对于实际计算而言是没有任何意义的，因此式(2.45)和式(2.46)的边界处理对于边界处的地震波反射而言并没有得到根本地解决。因此为了消除边界处的波形反射，PML区域中传播的地震波必须在能量上出现一定程度上的衰减，根据一阶常微分方程的解的一般规律，在这里我们对式(2.42)、式(2.43)、式(2.45)和式(2.46)所表示的PML区域中的地震波传播的规律中加入衰减项(刘璐等 2013)，也即：, (2.47), (2.48), (2.49), (2.50)其中和分别表示在坐标(x, z)处下标标记为x和z的波场的衰减系数，在实际应用中有多种方法来对不同空间坐标点的衰减系数进行选取，本文中所使用的衰减系数规律为，其中为最大的衰减系数， L为所添加的PML层的厚度，d为坐标点离已经离散的区域的边界的距离，对于不同的下标标记的衰减系数其相应的和L可以选用不同的值，本文中的值我们取为90。对式(2.47)、式(2.48)、式(2.49)及式(2.50)进行傅里叶变换将它们从时间域变换到频域可得：, (2.51), (2.52), (2.53), (2.54)其中, (2.55). (2.56)对于新坐标系下的情况而言可以得到类似于式(2.51)、式(2.52)、式(2.53)及式(2.54)所表示的地震波方程，对所得到的新旧坐标系下的方程进行有限差分处理，并进行相应的简化可得类似于式(2.37)的形式，也即：. (2.57)式(2.38)代回式(2.57)可得类似于式(2.39)的形式，也即：. (2.58)如果上式中的所有衰减系数就限定为0，那么式(2.58)可化为式(2.39)中的所有震源都为零的情形。如果我们把计算分为普通区域和PML层区域，并在不同的区域中分别使用式(2.39)和式(2.58)所表示的地震波传播方程来进行相应的正演计算，那么在形式上而在无形中就给计算增加了麻烦，这里为了计算的简单及形式上的简洁，我们把普通区域和PML层区域整合成一个区域计算区域，结合式(2.39)和式(2.58)的形式，我们可以得到描述此的方程可以表示为：. (2.59)对于普通区域及PML层区域中的每个离散网格点都有一个与式(2.59)相对应的式子，它们的总体可以组成一个方程的形式，可以用式(2.40)的形式来表示这个矩阵方程。对于添加的PML层而言，一般除了表层上层区域外都会添加不同厚度的PML层，而对于表层上层的区域而言，如果是对陆地进行勘探则不需要添加任何PML层，但如果是对海底进行勘探则需要添加相应厚度的PML层。而对于实际的介质地表而言，一般都是不规则或有较多起伏的，这时就需要有相应的方法对起伏的地表进行相应地处理来使之用于正反演过程，本文由于使用的是平地表模型，故而在此不对此问题加以论述。2.3 频域波形反演上文中我们已经介绍了给定区域中地震波的正演方法，本节中则主要对本文中将要使用的频域波形反演进行简要介绍。由正演的差分形式可知，由于正演数据相对于介质参数(速度值)而言并不是一个简单的线性关系，因此波形反演问题是一个比较复杂的非线性问题。反演问题的最终目的是使我们所得的最终模型与真实的地下介质模型相同或非常接近，但由于我们并不知道真实地下结构(否则就不需要反演了)，而仅仅知道在相应地理区域中由给定位置震源在相应检波器点所激发的波形响应，也即地震数据集，因此我们通过实际数据集与正演数据集的差异目标函数来衡量反演结果的优劣，目标函数值越小则认为反演结果越好。首先对正演模拟波场及勘探所得波场之间做差，并把结果记为波场残差，也即：, (2.59)其中的P和均为列矩阵，其中从上到下的每个元素都与离散后并添加了PML层后的网格点上的波场值一一对应，但由于勘探波场只是设置了检波器的网格点上的波场响应，故对于没有设置检波器的网格点所对应的元素值我们将其设置为0。我们可以使用波场残差的F范数来表示相应的目标函数(Brossier et al. 2010; Shin et al. 2007)，也即：. (2.60)其中M表示模型参数，在声波方程中一般表示声波网格点声波速度值或声波慢度或者为前两者的函数。由前面的讨论可知，反演的过程可以等价于式(2.60)所表示的目标函数收敛到全局极点值的过程。假设先验模型位于全局极小点附近，在此我们可以使用基于梯度的方法来完成对原模型的修正，首先应先求得目标函数对模型的整体参数的梯度，也即：, (2.61)其中Re表示取其后得数的实部，上标T表示矩阵的转置，J表示雅可比矩阵，其具体形式如下：, (2.62)可以看出J是二维矩阵，其第i行第j列的元素值可以表示为. (2.62)如果使用式(2.61)式来直接求取梯度矩阵，由于当模型规模增大时求取雅可比矩阵的时耗显著变大，因此不能有效地用于大规模模型情形，为了解决这个问题，这里首先对式(2.40)所表示的正演方程两边同时对模型参数M求偏导，考虑到其中震源函数本身与模型参数无关，可得, (2.63)其中A-1表示矩阵A的逆矩阵，可以看出上式实际上是雅可比矩阵的另一种表示形式，将其代回式(2.61)并考虑到矩阵A的近似对称性(Pratt et al. 1998)，可得. (2.64)其中的上标*表示矩阵的共轭。可以看出目标函数关于模型参数的梯度可以表示为正演波场与回传波场的零偏移相关函数表示(Gauthier et al. 1986; Mora 1987; Pratt et al. 1998; Ravaut et al. 2004; Tarantola 1984)。目标函数可以在模型参数M处按泰勒级数进行展开，我们忽略级数中二阶以上的高阶小项可得下式：, (2.65)其中H为黑塞矩阵，其具体形式为, (2.66)如果模型参数M表示真实地下介质模型，那么若对式(2.65)两端同时对求一阶导数，也即, (2.67)由于全局极小点存在且为M，式(2.67)的左端值为0，故可得. (2.68)上式即为所得的实际修正量，用此式所表示的波形反演方法通常被称作是牛顿法(或全牛顿法)(Pratt et al. 1998)。相关的大量实验表明式(2.66)所示的矩阵其值主要集中在第一项上，我们可以把此部分记为Ha，也即, (2.67)其中E表示单位矩阵，上标H表示矩阵的共轭转置，此时可称Ha为近似黑塞矩阵(Pratt et al. 1998)，将其代回式(2.68)时所得到的即为高斯牛顿法的表达式。若我们用式(2.67)所代表矩阵的对角阵来近似地表示完整的黑塞矩阵(此时的矩阵称为伪黑塞矩阵)(Plessix and Mulder 2004)，则模型参数的修正量可以表示为(Operto et al. 2006; Ravaut et al. 2004; Sourbier et al. 2009), (2.68)其中表示修正步长，为正则化因子，diag表示矩阵对角线上元素保持不变而其它位置元素置零。可以使用黄金分割搜索法、抛物线法、二分搜索法等多种进行求解，本文中采用抛物线法。可以看作是对梯度的预处理。式(2.68)表示的是对于单震源多检波器点的情形，对于通常使用的多震源多检波器的情形可以通过对单震源的情形的加权求和来确定其使用的模型参数的修正量。式(2.68)所表示的预处理梯度法，虽然仍然需要计算黑塞矩阵，但由于只是计算对角元上的元素，因此相对于式(2.68)所表示的牛顿法而言具体较小的计算费用，由于对梯度矩阵进行了预处理，故相对于通常的最速下降法而言有较好的收敛效率。由式(2.64)所表示的梯度形式及式(2.67)所表示的近似黑塞矩阵的形式可以看出，预处理梯度法的反演过程可以表示成多次的正演过程，考虑到反演时设置的震源与实际情形下的震源通常是不同的，故应对震源进行源估计(Koo et al. 2011; Pratt 1999; Operto et al. 2006; Ravaut et al. 2004)。3频率点的选取策略基于均匀背景模型，根据不同的原理或假设已经提出了多种频率点选取策略(或只是使用其它背景模型下的一些结果对所得策略做简单修正)，但就真实地下非均匀介质结构而言，此模型假设存在较大的区别。为了更好地用于真实地下介质反演，本章基于线性递增速度背景模型，根据散射波场Born近似及WKBJ近似及反演时梯度图像上竖直方向上的波数覆盖情形，提出一种新