时间空间域全波形反演

1、时间-空间域全波形反演任浩然，王华忠同济大学海洋与地球科学学院波现象与反演成像研究组，上海，200092摘 要 地震波反演成像的理论研究地震波场在地球介质中的传播现象，从地震数据直接反求物性参数。而全波形反演是利用全波场进行地震反演的方法。在反演框架中，正演方法是反演成败的基础。基于二次型泛函，非线性寻优可以采用梯度导引类的解法，也可以采用牛顿类的解法。我们主要讨论梯度导引类方法求解时间-空间域全波形反演的方案。基于模型试验，讨论了全波形反演方法的优缺点，并就全波形反演的实用化提出了相应策略。关键词：地震波反演，全波形反演，特征波场，非线性反演1 引言从Tarantola(1984)建立以波动

2、方程为基础的地震反演方法以后的近三十年来，地震波反演方法已经建立了较为完善的理论体系。这些方法可以分别从不同的模型参数化、不同的正演方法和不同的反演策略来区分。其中，介质模型的参数化与正演方法相互对应。对模型有什么样的认识，才能建立什么样的正演方程,有什么样的正演方程才会有什么样的反演策略。逆散射层析技术基于波场的线性近似，将地下的成像点看作是互不相干的散射点。传统上的散射层析只针对均匀背景模型（Devaney，1982；Wu和 Toksöz，1987）。后来发展的广义散射层析基于De Wolf近似(De Wolf，1985)，这种对散射波场进行多次前向散射和单次后向散射的近似使得广

3、义散射层析在非均匀背景上进行反演(朱小三，2010)。基于旅行时的层析反演方法可以分为：利用初至波或折射波的近地表速度反演(刘玉柱，2007)，利用反射波的深层速度反演（Billette和Lambare，1998），利用井间透射波旅行时反演井旁精细速度反演（Harris等，1995）等。基于射线理论的各种反演方法本质上还是模拟波场的旅行时信息，因此可以反演影响旅行时的速度等参数。然而，勘探工业发展的现实表明，仅仅研究速度参数是不够的。后期的地震解释要求偏移成像是一个保振幅的，即按照偏移或者反演估计出的反射系数随角度的变化应当能够反应真实的地下物性。这就要求对模型的参数化要更加贴近物理现实、正问

4、题过程要更加真实的描述地震波场的传播过程。另外，近年来随着计算机技术的发展，基于全波动方程的逆时偏移技术已经得到工业界大规模应用。因此，有理由相信，本世纪第二个十年将会是基于波动方程的反演方法大规模应用的阶段。当前，基于声波全波方程的反演还没能得到很好的应用，制约的因素主要有：如何估计一个较为精确的初始地球物理模型；如何构造更好的泛函以减少各种因素造成的振幅误差敏感性；正则化的问题如何处理，重建地震数据的低频信息。当然，对实际地震数据的FWI，仍然存在诸多制约因素，比如低信噪比、低频缺失、初始模型偏差太大、各向异性、粘性衰减，这些因素都有可能造成寻优过程的不收敛。在FWI的实用化过程中，面临的

5、主要问题可以归纳为：（1）建立一个较为精确的初始地球物理模型，体现对模型的某种假设（地质的或地球物理的）；（2）构建一个更好的泛函，以减少FWI对各种因素造成的振幅误差的敏感性，提高多参数反演中的稳定性；（3）构造适当的正则化参数，突出地质和偏移成像的先验信息参与度；（4）重建地震数据中的低频信息。从计算方法上，随着过去三十年中地球物理学家对利用迭代策略的全波形反演的研究,全波形反演的实现方式上主要有梯度导引类方法（如，Tape等，2007）和牛顿类方法（如高斯-牛顿法）。本文，首先从梯度导引类的全波形反演方法出发,分析声波方程下，全波形反演的梯度和步长的求取，并讨论迭代的优化算法。基于优化的

6、梯度类全波形反演策略进行程序实现，在简单模型上进行数值实验。针对全波形反演相关的一些问题进行了分析。最后，对全波形反演的实用化给出了我们的策略。2 全波形反演理论基础地震反演问题本质上是一个非线性问题。但人们在研究非线性问题时,常常在局部甚至全局线性化,以简化问题的求解过程。地震波场的正传播过程可以记为，(1)这里，为地震地球物理参数矢量，如速度、密度、弹性参数等，为观测到的地震数据, 描述了依赖于的地震波场正传播过程，它表达了地震波场传播的系统。在此系统下，研究特定模型中地震波的传播过程以得到数据的过程为地震正演。反之，由数据反推模型参数的过程为地震反演。反演的过程可以记为，(2)这里，即描

7、述了利用各种数学工具进行地震反演的过程。相对地，地震反演问题比正演问题要复杂，这是由于地震正传播得到的数据只把模型的部分信息带入了数据中。而且，模型和正问题的传播过程都只是实际情况的近似。是依赖于模型参数的非线性函数。定义误差泛函为，(3)这里，上标表示矩阵转置。定义Fréchet微商为，(4)在非线性问题中，Fréchet微商是依赖于模型的。相关问题在下面的讨论中会逐步涉及。非线性问题的求解主要有牛顿法和梯度法。牛顿法认为误差函数在初始模型附近满足二次型，因此可以用泰勒公式将误差函数在模型附近展开，(5)这里，是误差函数对模型的导数在初始模型处的值。是误差函数在处对模型的

8、二阶导数，有，(6)又称为Hessian算子，它有两部分构成。其中，为其非线性项目。这是因为对于非线性问题，Fréchet微商是依赖于模型的。由公式(5)，可以建立迭代公式，(7)不过这里的Hessian项本身是一个十分庞大的算子，其求解与求逆都是非常耗时的，本论文集中的另一篇文章中对此给予了详细论述。梯度法是泛函问题中最基本的一个数值解法，它把方程(7)的牛顿法反演转换成，(8)这里，为迭代步长，为迭代次数。迭代的方向即是目标函数梯度的反方向。步长的求解可以将误差函数在第次的模型上展开求解。这样就建立的梯度法求解全波形反演问题的理论基础。2.1 梯度的求取在声波方程意义下，设模型参

9、数为速度场（当然还可以定为慢度或慢度的平方，同样可以导出类似的方程）。记观测得到的炮集波场为，模拟得到的全波场为是与速度有关的量，表示正传播过程。其中，为炮点坐标，为与炮点相关的检波点坐标。由于我们不可能在整个模型空间观测数据，只能在地表（或井中等）观测，因此需要用一个算子来表述检波点位置限定，。数据残差为。测量这一残差受制于许多因素，如地震子波的影响、噪声的影响等。这里，首先考虑最基本的情况，即子波正确、噪声符合高斯分布（在统计意义下可以忽略）。误差泛函写，(9)可以用变分原理来推导梯度。由于限定算子和观测数据对模型无关，误差泛函的变分可以写为，(10)问题转化为求取。正传播波场的变分可以写

10、为，(11)这里，为速度场扰动。上式中，和均满足波动方程，记，代入波动方程，(12)(13)基于线性化假设，将方程（13）在处泰勒展开，并取其一阶项，合并式（12）可得，(14)方程式（14）就是波动方程的表达式。因此，其含义表示以右端项为源的波场传播出的波场为，即，(15)因此，(16)方程（16）的推导过程应用了两次共轭。其中，定义将限定在检波器上的数据残差空间扩展到整个模型空间（只需要将非检波器位置充零即可）。的应用基于共轭状态法，即局部扰动近似。即是将剩余波场在模型空间反传播，即逆时偏移。则泛函对速度模型的梯度为， (17)可见，梯度方向的值由正传播场在时间上的二阶导数与剩余数据的反传

11、波场在时间上的内积得到。每一炮的正反传波场都是独立计算的。最终求取的梯度值为每一炮求取结果的叠加。式（17）涉及一炮的两次波场传播。可以整理梯度求取的实现方案：（1） 对每一炮，a) 计算并记录每一时刻的正传波场；b) 计算检波点位置的波场残差；c) 将波场残差反传播到模型空间，得到反传播数据；d) 求取反传波场与正传波场在时间上的二阶导数在时间上的内积，得到单炮的梯度。（2） 叠加所有炮求取的梯度值，得到与模型空间的全局梯度。2.2 迭代步长梯度导引类的反演思想除了求取迭代方向，还需要求得一个迭代步长。这里我们使用抛物插值法。抛物插值假设误差泛函在模型附近呈抛物线型，通过给定一个试探步长求取

12、一个优化的步长值。首先，给定一个全局的试探步长（如将速度更新量控制在背景速度的1%），以此更新速度得到一个伪更新速度模型。按照这一速度可以做正演并求取误差泛函，(18)将误差泛函在周围作二阶泰勒展开，有，(19)其中，为一个跟误差泛函对模型的二阶导数有关的量。已知当前误差函数、基于伪速度场的误差函数、梯度和试探步长的情况下，由方程（19）可以导出的解，(20)如果存在一个最优的步长，则在该步长周围，应当满足。基于最优化步长的误差泛函在周围展开，同样应到满足方程（19）的形式。由于误差泛函在模型周围是二次型的，则按照和得到的展开式中不变。因此，可以由和方程（19）可以得到一个最优化的步长，(21

13、)可见，利用这种抛物插值的思想求取迭代步长的过程需要进行一次波场正传播。即基于一个伪步长更新后的速度，求取一个伪的正传数据，计算出误差函数值，用以求取最优化步长。在正向波场求取的过程中，可以实现炮的并行计算。最后将每炮计算的误差函数归约即可。2.3 FWI流程得到梯度和步长以后，速度更新过程按照下式执行，(22)这样,就可以实现迭代的速度更新，梯度导引的反演方法得以实现。当然，还需要设定一个迭代终止的条件，可以为，(23)上式规定了一个迭代收敛的规则，即第n次更新的量小于n-1次迭代的一定比例，如，则认为迭代收敛，反演可以终止；如果这个条件不满足，则将更新结果作为输入，进行下一次迭代。可以给出

14、时空域声波方程全波形反演的过程如图1所示。图1表示了梯度法全波形反演的流程。在每次迭代过程中，求取梯度的部分需要计算两次全波场模拟，而基于单炮的策略容易实现炮的计算机并行。基于步长抛物寻优的策略需要计算一次波场正演，此过程中，同样可以实现炮并行以求取最优步长。每一次迭代过程中的计算耗时主要是用于三次波场模拟的过程。三次波场模拟都可以实现炮并行。基于多节点多核CPU计算机群的并行化能够为全波形反演的实用化提供硬件支持。给出一个试探步长计算并记录正传波场求取单炮数据残差给定一个初始模型计算并记录正传波场求取波场数据残差反传残差波场得到单炮梯度炮并行得到全局梯度炮并行得到全局数据残差得到迭代步长更新

15、速度，结果满足否？输出结果，反演结束是否图1 梯度法全波形反演流程图当然，由于梯度法本身的“折叠效应”，需要在具体求解过程中对求解的梯度做进一步优化。这就引出了共轭梯度法与拟牛顿法（袁亚湘，2008）。任浩然（2011）对此进行了分析，本文不再详述。 全波形反演模型试验首先，设计一个简单模型。这一模型要求能够考察方法的有效性，能够检测算法对初始模型的依赖程度，还要能够反应复杂波场传播规律。首先，我们设计一个凹陷模型。为了考察其界面信息，在凹陷的下面再增加一个反射界面，以利用反射数据更新凹陷周围速度场。如图(a)所示。凹陷模型参数为：速度网格大小，纵横向间距均为8。水平坐标范围。地表放炮，共放4

16、9炮，炮间距32，第一炮的位置在第32处。利用主频为25hz的零相位雷克子波激发，地表全孔径接收。时间采样率，采样点数为2000。三层结构速度值分别为3500，4000，4500。(a) (b)图凹陷模型速度场。(a)真实速度模型；(b)光滑速度模型通过这一简单模型，我们希望验证以下几个问题：（1） 方法的有效性。可以给定一个较好的初始模型，测试梯度导引类的全波形反演能否给出一个较为精确的反演结果。（2） 算法的收敛性和不同算法的收敛性对比。在同一模型下，运用不同算法进行计算，考察误差泛函减小的速度。（3） 全波形反演对初始模型的依赖程度。给出几种不同的初始模型，测试全波形反演在什么样的速度背

17、景下能够收敛到正确值。实验一：方法有效性检验首先，检验方法的有效性。给定一个真实速度模型光滑之后的模型（Seismic Unix软件中，取加权因子为20），如图(b)所示。利用光滑速度模型进行BFGS全波形反演，设定最大151次迭代。图给出了第1次迭代、第11次迭代、第51次迭代和第151次迭代求得的梯度方向。可见，在初始阶段，梯度方向近似于逆时偏移。但第一次迭代的梯度明显较逆时偏移的像更为尖锐，我们已经在公式中进行了讨论，全波形反演天然的比逆时偏移对高频来的敏感。经过多次迭代以后，速度已经逐步收敛到真实解附近，数据残差逐渐减小，在梯度中已经观察不到像的信息。将第1次迭代、第11次迭代、第51

18、次迭代和第151次迭代更新后的速度场列在图中，可以看到，速度场能够逐渐收敛到真实解，方法是有效的。分别用初始光滑速度模型和迭代第151次的速度模型进行逆时偏移（图）可见，全波形反演后的速度场为逆时偏移提供了一个很好的初始解。(a) (b)(c) (d)图 凹陷模型BFGS法梯度方向。(a)第1次迭代；(b)第11次迭代；（c）第51次迭代；（d）第151次迭代。(a) (b)(c) (d)图 凹陷模型BFGS法更新速度。(a)第1次迭代；(b)第11次迭代；（c）第51次迭代；（d）第151次迭代。(a) (b) (c)图 逆时偏移结果对比。(a)使用真实速度；（b）适用初始光滑速度；(c)使

19、用BFGS法第151次迭代速度。 实验二：不同方法收敛效率检验为了验证方法的收敛性，将误差泛函在每一个迭代步中记录下来，进行分析。利用已经验证有效的光滑速度模型，分别用梯度法(SpD)、共轭梯度CG法、拟牛顿BFGS方法进行151次迭代，其中BFGS法中保存前十次迭代的梯度和模型。图6中分别给出了三种方法的在第151次迭代的速度值，以及误差泛函的收敛过程。将三种方法计算的泛函做对数变换后放在一张图中（图），可以看到：（1）梯度法（最速下降SpD），CG法和BFGS法在此模型上均能收敛到真实值，误差泛函的减小说明的方法的有效性；（2）SpD方法误差泛函的收敛过程存在“折叠”现象，这与理论分析是吻

20、合的；（3）CG法和BFGS法的计算效率明显优于SpD；（4）本模型试验中，BFGS法收敛最为稳定。图 三种迭代策略的收敛效率比较。自左到右分别为SpD,CG和BFGS；上为三种方法第151次迭代结果，下为收敛效率。图7 三种方法误差泛函对数变换交汇图。实验三：全波形反演对初始模型依赖性测试前面两个实验，用一个真实解附近的光滑模型作为反演的起点。然而，对于实际地震数据，如果没有一个好的速度分析预处理，很难得到这么好的初始速度。因此，需要对全波形反演对初始模型的依赖性进行测试。给出另外两种初始模型：一个是常速模型，速度取上层速度值3500；一个是常梯度模型，即速度从浅到深从3500到4500按照

21、常梯度变化。这两种速度都较严重地偏移了真实速度模型，前者在第二个速度层偏离了原速度值的12.5%；后者在浅层的速度值错误，垂直入射至凹陷底部的旅行时差约为7.1%。图8和图9给出了这两种方法的速度更新过程。从图中可以看到，当初始速度偏离真实速度较多时，FWI难以收敛到真实解。在本模型中，尤其是常背景速度在更新最下层位置时无能为力；尽管经过151次迭代，下层界面的位置仍不能收敛到正确位置，而是在一个假的位置上“徘徊”，即陷入局部极小值。常梯度模型在凹陷部位同样出现类似状况。可见，FWI对初始模型的准确度是要求较高的。(a) (b)(c) (d)图8 常速模型BFGS法更新速度。(a)第1次迭代；

22、(b)第11次迭代；（c）第51次迭代；（d）第151次迭代。(a) (b)(c) (d)图9 常梯度模型BFGS法更新速度。(a)第1次迭代；(b)第11次迭代；（c）第51次迭代；（d）第151次迭代。以上在凹陷模型上的试验证明了FWI的有效性。作为一种十分“娇贵”的方法，其在模型试验上能够得到不错的效果。国内外现有许多研究成果都是在Marmousi模型上进行的。作者的前期研究也在该模型上取得了效果（任浩然，2011）。可以对于实际地震数据，其实用化却面临着巨大的挑战。 全波形反演实用化的途径全波形反演在实际数据上应用面临着巨大的挑战，这些挑战主要包括：子波未知；噪声的干扰；其他波现象的干

23、扰；正传播算子与实际波现象之间的差距；以及局部寻优条件的不满足。图10给出了某实际数据的初始模型以及模拟数据和接收数据之间的巨大差距。这一模型很典型的说明了实际数据中的噪声和复杂波现象的干扰。在这种情况下，单靠数据匹配的传统FWI是完全不能得到实用效果的。那么FWI实用化的出路又在何方呢？既然基于地震反演框架的FWI是一个完整的体系，加入其将来能够实用化，那么我们就需要把实际数据的处理过程纳入FWI的框架中。首先，在数据处理方面，需要重视和加强数据的预处理，目的就是提供消除了不能模拟成份的数据给FWI。一方面，我们需要通过对正演算子的设计尽可能多的模拟出多种波现象；另一方面，消除正算子不能模拟

24、的波场成分。当然，在这个过程中，需要特别注意保护有效信号，尤其是数据的低频成分。与正演过程相关的，还需要提供一个能够与实际传播过程匹配的子波。对子波的估计有几种策略。一是利用海上数据的震源记录；二是利用直达波来估计震源子波；三是还可以把子波纳入反演的框架中。这几种方法中以海上数据记录最为直接。在有效子波估计的基础上，通过考虑各向异性和有效边界条件等因素确保正传播算子与实际传播过程的一致性。图10 模拟数据与实际数据之间的巨大差距。(a)初始模型；(b)该区某单炮；(c)基于初始模型的模拟数据。在反演策略方面，还可以进行一系列的优化。如通过采用敏感度更低的泛函降低其对初始模型的依赖性，通过将数据

25、拟合类泛函改造成数据相关泛函来消除子波估计的误差。在泛函过程中，可能需要提取波场的特征成分来进行特征波场反演（CWI）来代替FWI。这种特征化表现在几个方面：1. 特定（低的）频率的地震数据；2. P波数据；3. 单向波数据；4. 时间-空间局部化数据。需要指出的是，多尺度反演的思想可以涵盖在这种特征化中。通过合适的数据预处理和巧妙设计的反演流程，特征波形反演是波形反演思想能够在陆上数据实用化的途径。 结论与讨论本文重点对梯度导引类的全波形反演实现方案进行分析。正过程采用目前实用的时间-空间域有限差分正演方法。梯度的求取需要进行两次正演，一次沿时间方向的正向传播，一次反向传播。步长的求取我们采

26、取局部二次寻优的方案。基于优化的梯度法，实现了全波形反演，在一个简单凹陷模型的数值实验，让我们对FWI具体实现的一些细节进行了分析。基于本文的分析，从应用的角度来看，应该形成如下的技术体系：数据处理的初期利用CMP道集中的常规叠加速度分析得到粗糙的速度模型；在PSTM道集中进行偏移速度分析，如剩余曲率分析方法。如果数据质量不能保证，譬如山地地震勘探时，PSTM扫描速度分析也是有效的方法。对于低信噪比数据，考虑Fresnel带同相叠加的成像方法要引入速度估计过程中。PSDM偏移速度分析是速度估计和建模的核心环节，扩展成像道集的深度偏移速度分析是比较精确的方法。当然，在成像道集中进行层析速度反演可

27、以进一步提高速度估计的精度，尤其是利用Beam-ray的思想建立投影与反投影关系可以更好地估计背景速度。在成像道集上进行如此的速度估计与偏移成像方法匹配得很好，这也是它成为目前最常用的深度偏移速度估计方法的基本原因。在上述各种方法提供的较为可靠的长波长速度模型的基础上，进行FWI速度估计，应该得到更为细致的、包含更多短波长速度分量的速度估计结果，最好是得到能满足储层描述需要的速度模型。目前的FWI还是与RTM一起构成一套精确成像的技术体系，还没有推进到储层描述的阶段。全波形反演的研究工作推进到实际数据的应用仍需几年的时间，但一个明显的导向是，需要用特征波场反演CWI来代替全波形反演FWI，才有

28、可能将基于全波的地震反演应用到复杂的实际地震数据处理的过程之中。参考文献1 刘玉柱,董良国. 近地表速度结构初至波层析影响因素分析.石油地球物理勘探，2007，42(5): 544-5532 任浩然声介质地震波反演成像方法研究: 博士学位论文. 上海：同济大学海洋与地球科学学院. 20113 袁亚湘.非线性优化计算方法.北京: 科学出版社,20084 朱小三. 地震散射点成像和非均匀介质中广义散射层析成像反演: 博士学位论文. 北京：北京大学地球与空间科学学院. 20105 Billette F. and Lambare G., Velocity macro-model estimation

29、from seismic reflection data by stereo-tomography. Geophysical Journal International, 1998, 135(2): 671-6806 De Wolf D. A. Renormalization of EM fields in application to large-angle scattering from randomly continuous media and sparse particle distributions. IEEE trans. Ant. And Propag, 1985, 33: 60

30、8-6157 Devaney A. J. A filtered backpropagation algorithm for diffraction tomography. Ultrasonic Imaging, 1982, 4: 336-3508 Harris J. M. Nolen-Hoeksema R. C., Langan R. T., et al., High-resolution crosswell imaging of a west Texas carbonate reservoir: Part 1-Project summary and interpretation. Geoph

31、ysics, 1995, 60(3): 667-6819 Tarantola A., Inversion of seismic reflection data in the acoustic approximation. Geophysics, 1984, 49: 1259126610 Tape C. Liu Q. and Tromp J. Finite-frequency tomography using adjoint methods: methodology and examples using membrane surface waves. Geophys. J. Int., 2007

32、, 168: 1105112911 Wu R. S. and Toksöz M. N. Diffraction tomography and multisource holography applied to seismic imaging. Geophysics, 1987, 52: 11-25胜利油田埕东二元复合驱现场实施过程中，出现了因污水中H2S的存在而导致所配聚合物溶液粘度损失严重的问题，所配聚合物溶液（质量浓度为17000mg/L）粘度仅有5mPas，远低于该二元复合驱对聚合物溶液粘度的设计要求（20 mPas），这严重制约了二元复合驱作业的顺利实施。本论文通过运用实验定量的方法，研究了不同浓度的硫化