CN109061727A - 一种频率域粘声介质全波形反演方法 - Google Patents

Abstract

本发明提供一种频率域粘声介质全波形反演方法，包括：采用截断牛顿法对速度与Q值进行联合反演选用最优化混合交错差分格式对波动方程进行离散化。本发明方法改善了伪影现象，提高了联合反演的效果和效率。

Description

一种频率域粘声介质全波形反演方法

技术领域

本发明涉及地震勘探技术领域，尤其涉及一种频率域粘声介质全波形反演方法。

背景技术

在地层Q值(或吸收系数)与速度联合反演方面，国内外研究都相对较少。Tarantola(1988)最早推导了粘弹介质中吸收系数与速度联合反演的理论基础，指出模型可以通过正向传播波场与反传波场的互相关进行更新。Liao和McMechan(1996)在频率域对其进行了算法实现并做了简单的模型测试。Causse等(1999)提出了针对粘声介质的预条件方法，补偿由于频散和振幅衰减效应对梯度法反演结果的影响。Malinowski等(2011)对粘声介质反演目标函数进行了敏感性分析，指出震源子波估计的好坏对反演结果的重要性，最后通过模型试算及实际资料处理验证了方法的可行性。Kamei和Pratt(2013)针对粘声介质提出了初步的反演策略，并利用井间地震模型对方法进行了评估。而国内关于粘声介质全波形反演方面的文献更是鲜有所见。于长澎(2011)在其硕士论文中，基于反射率法对层状介质地层吸收系数进行了反演算法研究及模型试算，但基于有限差分算法的吸收系数反演仍然只停留在总结国外的经验和成果上。Ren和Liu(2014)等利用变网格时空域有限差分方法和基于第二代小波变换的多尺度FWI算法进行了粘声介质波形反演的研究，有效的提高了反演准确性和计算效率，并且通过模型试算取得了较好效果。虽然国内学者已针对粘声介质的反演方面做出了有益尝试并获得一定成果，但总的看来，国内外针对粘声介质的全波形反演方法仍然处于初期阶段。

发明内容

本发明主要目的就是针对以现有技术的不足，提供一种粘声介质全波形反演方法、以及多参数反演的优化方法和稀疏约束方法。

一种频率域粘声介质全波形反演方法，包括：采用截断牛顿法对速度与Q值进行联合反演。

进一步地，如上所述的方法，所述截断牛顿法实现的过程分为内外两个循环结构，外部为传统的模型更新形式，内部的循环为共轭梯度形式的线性问题。

进一步地，如上所述的方法，所述联合反演包括：

首先反演速度，将反演的速度结果和初始Q值模型作为第二步反演的初始模型；

其次，对速度模型与Q值模型同时更新，得到最终的反演结果。

进一步地，如上所述的方法，使用稀疏约束方法对速度与Q值进行联合反演。

有益效果：

本发明为适应日渐复杂的油气勘探目标，选择更接近实际地层性质的粘声介质方程对地层速度与品质因子(Q值)进行全波形反演，对描述油气储层性质、落实勘探目标、深化地质认识具有重要理论研究意义和应用价值。

附图说明

图1为一维散射示意图；

图2为本发明多尺度频率域全波形反演流程图；

图3为目标函数随Q值相对变化的变化曲线；

图4为目标函数随速度相对变化的变化曲线；

图5(a)为异常体模型的速度参数模型；

图5(b)为异常体模型的Q值参数模型；

图6(a)为异常体模型联合反演得到的速度结果；

图6(b)为异常体模型联合反演得到的Q值结果；

图7(a)为速度反演结果中x＝660m处的速度曲线与真实模型值对比图；

图7(b)为Q值反演结果中x＝1320m处的Q值曲线与真实模型值对比图；

图8(a)为复杂模型的真实速度模型图；

图8(b)为复杂模型的初始速度模型图；

图8(c)为复杂模型的真实Q值模型图；

图8(d)为复杂模型的初始Q值模型图；

图9(a)为复杂模型联合反演得到的速度反演结果；

图9(b)为复杂模型联合反演得到的Q值反演结果；

图10为截断牛顿法的算法流程图；

图11(a)为截断牛顿法异常体模型联合反演得到的速度结果；

图11(b)为截断牛顿法异常体模型联合反演得到的Q值结果；

图12(a)为异常体模型的速度反演曲线与真实模型值对比图；

图12(b)为异常体模型的Q值反演曲线与真实模型值对比图；

图13为异常体模型的两种收敛算法的反演误差对比图；

图14为异常体模型的两种收敛算法的归一化目标函数收敛曲线对比图；

图15(a)为复杂模型采用截断牛顿法一步法联合反演得到的速度结果；

图15(b)为复杂模型采用截断牛顿法一步法联合反演得到的Q值结果；

图16为复杂模型采用两步法反演策略时第一步的速度反演结果；

图17(a)为复杂模型采用两步法反演策略时得到的速度反演结果；

图17(b)为复杂模型采用两步法反演策略时得到的Q值反演结果；

图18为复杂模型采用一步法和两步法反演策略得到的反演模型误差对比图；

图19(a)为Marmousi模型的真实速度模型图；

图19(b)为Marmousi模型的初始速度模型图；

图20(a)为1hz频带下的多源频率域正演波场快照图；

图20(b)为10hz频带下的多源频率域正演波场快照图；

图21(a)为Marmousi模型无约束条件下的反演结果；

图21(b)为Marmousi模型Seislet稀疏约束条件下的反演结果；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的发明思想是：介绍地层吸收系数的参数化方式，引入品质因子Q，建立了粘声介质波动方程。在频率域正演的基础上，结合模型对双参数联合反演进行的试算，并提出了优化方法，在收敛算法上使用截断牛顿法，反演策略上使用两步法反演。同时为了提高反演的效率和稳定性，介绍了基于“多源”全波形反演的稀疏约束方法。

本发明主要构建了频率域粘声介质全波形反演的流程，对速度与Q值进行了联合反演，并采用多种方法针对联合反演结果进行优化。本发明主要的创新点体现在，将截断牛顿法和稀疏约束方法应用于频率域粘声介质全波形反演过程中，获得精确的Hessian算子信息的同时可以应用于“多源”反演，使得联合反演的反演结果和效率取得了一定程度上的优化。

本发明是通过以下方案实现的：

无约束/无权重下的最小二乘目标函数C(m)可以写成如下形式：

其中式子中^T表示转置符号，^*表示共轭符号。局部优化方法需要在初始模型m₀附近寻找最小值，因此第一步需要将目标函数在m₀线性化函数。需要存在一个很小的扰动量，在Born近似的框架内(Born and Wolf,1993；Beydoun and Tarantola,1988),m可以写成初始波形m₀与扰动模型δm之和，写成m＝m₀+δm。

对目标函数C(m)进行二阶泰勒展开如下所示：

其中M为向量m的维度。

当目标函数的一阶导数为零的时候，其能达到最小值。关于模型参数的一阶导数表达为：

其中为梯度，为Hessian矩阵。设定为零，给出下降方向的表达式为：

该下降方向对应于牛顿方向。接下来分别介绍梯度和Hessian的表达式，由于本发明考虑的是频率域下的表达式，所以数据是复数值。

目标函数的一阶导数表示为：

其中N代表数据向量的数目并且与采集参数有关，即炮点和检波器的数量。其中为数据的实数部分，上述式子可以改成更紧凑的形式：

其中J₀表述为Fr′echet矩阵。

依照梯度的表达式，Hessian可以写成：

其中第一项成为Hessian近似项。根据对角线零滞后自相关定义，该项是对角占优的这个对角线项减少了几何扩散的影响。因此，在地表采集的系统中，其有助于衡量深部扰动(大偏移距/小振幅)。

其中表示的波场二阶偏导数与数据残差之间的零滞后相关性。

将上述式子整合在一起，更新方向方程(牛顿方向)可以改写成：

牛顿方向具有局部二次收敛性质。通常情况下，Hessian矩阵的第二项可以被忽略(Pratt et al.,1998)，这就是高斯-牛顿方向，表示如下：

如果Hessian矩阵简单用常数的来表示，则是最速下降法：

考虑到Hessian矩阵的计算和存储成本，对其的使用大部分方法的思路都是取其的近似值。其中的一种途径就是利用Hessian矩阵的对角近似。该近似方法的Hessian矩阵可以定义为：

灵敏度矩阵元素可以写成：

其中k(s,r)表示为震源接收点耦合对，其中的s表示震源位置，r表示为接收点位置。δ_r表示激发震源在接收点r处。为了建立灵敏度矩阵，我们需要对每个震源和接收点位置建立一个正演问题。这也就意味着近似Hessian矩阵的计算严重依赖于采集系统。

Shin et al.(2001)提出了一种对角伪Hessian矩阵，这一种对已经提出的近似Hessian矩阵的进一步近似。事实上，在操作过程中只是使用到了虚震源而不是A^-1f_i，伪Hessian矩阵Η_a可以写为：

其中F_(i,j)＝(f_i)_j表示虚震源矩阵。该表达式将近似Hessian矩阵的计算成本降低到与梯度计算成本相当的水平。对于实际操作过程中，对于伪Hessian矩阵的求取可以与目标函数梯度的运算同时进行，这大大的减少了计算与储存的工作。

有些情况下，需要考虑Hessian矩阵的非对角项才能达到理想的收敛效果。要达到这一目的，就要使用到BFGS方法，Nocedal(1980)BFGS方法的高效形式，称之为L-BFGS方法，较低了计算与储存的要求。

准Hessian矩阵的逆写成BFGS格式为：

其中s_k＝m_k+1-m_k，这种关系是递归的。在给定迭代下，准Hessian矩阵的逆仅仅取决于之前迭代的模型和目标函数的梯度。

在L-BFGS算法中，模型向量和目标函数的梯度仅迭代有限的次数。另外，准Hessian矩阵的逆不会直接被存储。

对于一个震源一个频带来说，该方程关于模型参数m的偏导数表示为：

其中表示为偏导数波场。在接收点附近提取的偏导数波场可以建立雅克比矩阵J，因此，其中R是在接收点提取的限制算子，在梯度表达式中加入雅克比矩阵以及关系式，得到：

其中A^-1的每一列代表激发震源在模型空间节点处的格林函数。由于格林函数在空间上的互异性，矩阵A^-1是对称的，因此可得到：

其中R_b为反传剩余波场。共轭符号*表示反传，公式中可以看出，对梯度的计算变为两个正演的问题，第一个是计算波动方程的解，第二个是波场残差的反传波场。

这种计算梯度的策略称为伴随法，它的理论源于约束优化的拉格朗日公式，对于多源和多个频带的梯度表达式为：

在实际应用的过程中，震源子波是未知的。为了更新震源，最小化未知震源下的目标函数。地震波场和震源之间的关系是线性的，因此：

P＝Gs

其中G是格林函数，s为震源子波。这目标函数可以改写成：

其中G_R是接收点位置处的格林函数，由上式导出震源的更新表达式：

在每次迭代的过程中，震源随着模型扰动的更新而更新。

全波形反演的最佳空间分辨率是波场的一半。首先考虑一个速度为c的均匀背景模型，一入射的平面波沿着震源-散射体方向传播如图1所示，入射和散射的格林函数可以近似写成：

G₀(x,s)≈exp(iks,x)

G₀(x,r)≈exp(ikr,x)

其中k＝ω/c，x,s,r分别表示为散射点，震源点和接收点向量。则目标函数的梯度表达式为：

梯度表达式具有截断傅里叶级数的形式，其中积分变量是散射波数向量k＝k(s+r)，并且级数的系数是数据残差。在给定频率的情况下，采集系统通过震源点与激发点数量和来控制该傅里叶级数的采样点与带宽的。为了更清楚的看到空间分辨率与采集系统之间的关系，式子中的波数k可以表示为如下形式：

其中n为慢度向量(s+r)的单位方向向量，这更突出了全波形反演的空间分辨率与频率-偏移距之间的关系。数据空间中的一个频率和一个孔径对应于模型空间中一个波数。所以，其具有对波数覆盖冗余控制性，这种冗余性会随着孔径带宽而增加。有学者已经提出，在频率域全波形反演中将反演过程限制为几个频带的反演，就可以抽取反演中的波数覆盖冗余。实现这样的做法，相对于时间域全波形反演，频率域全波形反演可以减少计算量以及能够处理更为复杂的数据。对于全波形反演过程中频率的选择要求是：该频率成像的最大波数与下一个频率最小垂向波数对应，其具体的表达形式为：

Δf_n+1＝f_n+1-f_n＝(1-α_min)f_n+1

其中α_min是取决于偏移距与反射体深度的常数。根据这个原则，频率间隔随着频率增加而增加。

低频和宽孔径的数据有助于分辨介质的中、长波长部分。在波数谱的另一边，如果观测数据中有的话，垂直入射波会带来半波长的最大分辨率。

基于以上空间分辨率理论与频率的选择方程，设计多尺度反演策略。通过对频带的分组，由低频、中频、高频的顺序依次反演，实现对目标介质的有效反演。

根据图2所示的多尺度反演流程，主要分为三个循环，首先是频率组的循环，以及迭代次数的循环和频率的循环构成。当反演流程只有一个频组时，多尺度频率域全波形反演退化为普通的多频全波形反演。

在对速度与品质因子进行联合反演之前，首先针对速度与品质因子对于目标函数的敏感性进行分析。

对于敏感性的分析采用网格点法，首先设定一个初始均匀介质模型，模型大小为101*101m，初始的速度设为3000m/s，初始的Q值设为50，密度参数固定为1000kg/m³。对其中一项参数进行固定，改变另一项参数的数值，计算出改变模型后目标函数的值，因此可以得到目标函数随着反演参数变化的曲线图如图3和4所示，该曲线图可以分析该参数对目标函数的敏感性。参数模型的相对变化量可以定义为：

通过上述定义，算出速度参数与Q参数的目标函数相对变换量曲线图(图3和4)。首先我们从目标函数变化量的数据量级上可以看出，速度参数和Q参数的对应目标函数变化量数量级分别为8次方和6次方，也就是说Q参数比速度参数少了2个数量级，因此可以看出，速度在对目标函数的敏感性上要远远大于品质因子Q的。仔细观察Q值的目标函数变化量曲线可以看出，在小于初始Q值变换量下，目标函数的变化量不大，大于初始Q值的情况下，其目标函数的相对变化量较大。而速度参数的目标变化量分布类似于二次曲线，这种情况比较易于收敛。通过分析，相对于速度的反演，Q值的反演难度明显较大。针对这一情况，我们需要多级优化方法，一步步改善反演效果。

双参数的同时反演，由于参数之间的耦合效应，有可能会导致一个参数的反演结果中出现另一个反演参数的伪影现象。因此为了更好的展示反演的效果，设计了异常体模型针对联合反演的试算。

两个模型分别如图5(a)、(b)所示，大小都为101*101，网格点大小为20m*20m,其中速度模型中异常体速度为1800m/s，背景速度为1500m/s。Q值模型中异常体速度为10，背景Q值为100。反演过程中的初始模型分别为各自参数模型中的背景速度模型，密度为1000kg/m³。反演的一组频率为5hz，10hz，15hz，反演的迭代次数为30次，反演的收敛算法为L-BFGS法，得到反演的速度和Q值结果如图6(a)、(b)所示，为了更直观的对反演结果进行分析，我们选取速度反演结果x＝660m处速度曲线以及Q值反演模型的x＝1320m处的Q值曲线，并与相同位置处的真实模型值相对比，如图7(a)、(b)所示。

可以从反演结果和曲线图中可以看出，在双参数联合反演的过程中，速度与Q值都能得到较好的恢复。为了更加精细的对比各参数反演的效果，我们引入模型误差的计算对各参数反演结果进行定量的分析，定义模型误差公式如下所示：

其中nx，nz分别为横向采样点和纵向采样点，m_inv为反演得到的参数模型，m_ture真实参数模型。

通过以上公式的计算得出，其中速度的模型误差为0.0031，Q值的模型误差为0.0553。可以得出速度反演的结果较为准确，Q值反演结果中除了Q值异常体部分得到了成像，在速度位置处，Q值的反演结果出现了伪影的现象，这是因为速度参数与Q值相互耦合造成的结果。这种耦合现象，对速度的影响也是存在的，但是较小。但是对于Q值的影响较为明显，这也是之前分析Q值和速度参数对于目标函数的敏感性的得出的结论。如何解决这个问题，是本申请的主要目的。

接下来针对复杂模型进行联合反演的试算，分别为真实速度模型、初始速度模型、真实Q值模型和初始Q值模型，如图8(a)-8(d)所示。速度模型参数和Q值模型参数分别跟之前对应的单参数反演模型一致，反演的频率为四组频率，分别为1～4hz,4～8hz,8～12hz和12～16hz的四组频率。每组频率迭代20次，反演的收敛算法为L-BFGS法，获得速度参数与Q值联合反演结果如图9(a)、9(b)所示。

从反演结果可以看出，在联合反演的过程中，速度的反演效果较为理想，但Q值的效果不好，基本上更新的程度很小，在浅层部分其结果相对于深层部分较好。这可能跟波的能量有着直接关系，在深层由于吸收衰减，其能量较弱，反演的效果也不理想。因此在联合反演中Q值的反演结果需要在一定程度上得到优化。

针对反演存在的一些问题，尤其是对双参数联合反演，有必要进行一定程度的优化处理，旨在获得更好的反演结果。本发明提供了三种主要的优化方法。截断牛顿法(Truncated Newton，TN)针对联合反演能够使参数耦合性降低，获得更好的Q反演结果；针对联合反演问题，提出更优化的反演策略，提高联合反演的效果；最后为了提高反演的效率，针对频率域“多源”全波形反演，在模型更新的过程中，使用Seislet域稀疏约束，能够在提高反演效率的同时降低串扰噪音对反演结果的影响。

截断牛顿法能够利用更为精确的Hessian算子。其通过计算Hessian矩阵和已知向量乘积的形式，巧妙的保证充分精确利用到Hessian算子的信息，又能够减轻计算的内存消耗。对于多参数反演，截断牛顿法能够减轻不同参数之间的耦合效应。

截断牛顿法实现的过程如图10所示共分为内外两个循环结构，外部为传统的模型更新形式，内部的循环为共轭梯度形式的线性问题，首先计算Hessian矩阵与提前设定的已知向量的乘积

式子中v是已知向量，h为一接近于0的正实数。在对H(m_k)v计算过程中，需要求取两次梯度，两次梯度对应为四次正演的计算，相对来说其计算量是可以接受的。

为了保证精度的前提下提高计算效率，根据Eisenstat和Walker的描述的方法设立终止条件，其终止条件与目标函数的二阶近似精度相关，设定为如下形式

式子中η_k为约束因子，η_k值描述了目标函数的近似精度，其值越小精度越高，给出η_k的定义为：

这里利用一个简单模型进行试算，使用与前面相同的采集参数，频率组，迭代次数等反演策略，仅仅改变收敛算法这一个条件，得到截断牛顿法对异常体模型的联合反演结果如图11(a)、11(b)所示。

对于反演的结果分析来看，截断牛顿法的反演效果更好，其伪影现象也得到了一定程度的改善。

通过图12(a)、(b)所示的更新值与真实值的对比曲线图和图13所示的采用两种收敛算法所得的模型误差对比分析图还可以得出，相对于L-BFGS算法，截断牛顿法在其他反演条件不变的情况下，其速度和Q值模型误差都相应的减小，尤其对于Q值模型误差，减小的幅度更大。

另外，我们针对该反演频率组内的归一化目标函数曲线进行对比分析，如图14所示。从归一化的目标函数收敛曲线对比分析出，在该异常体速度与Q值的联合反演中，截断牛顿法相对于L-BFGS法收敛的速度的更快，收敛的效果也更好。体现了截断牛顿法在联合反演中的优势。但也存在一定问题，虽然截断牛顿法对Q值有一定程度的改善，但是效果仍然不明显。值得进一步优化改善Q值的反演结果。

对于速度与品质因子Q的联合反演的策略，基本上可以分为一下几种情况：1.顺序反演，即先反演速度再反演Q值，或者先反演Q值再反演速度。2.同时反演，即速度与Q值同时反演。3.两步法反演策略，即第一步，固定Q值参数，只更新速度模型，得到速度模型的反演结果后，再同时反演速度和Q值。Askan(2006,2007)在研究中指出，由于速度参数和Q值参数对于目标函数的敏感程度不同，速度反演结果较少受到Q值结果的影响，但是Q值的结果受到不精确速度的影响较大，因此首先在Q值不更新的基础的上，对速度进行反演，直到速度精度达到一定程度后，将速度反演得到的结果作为初始速度场，与初始Q值同时进行反演，由于此时速度更新量相对较小，主要是对Q值的更新。

首先输入初始Q值模型，以及速度模型，但是只对速度模型进行更新，Q值模型一直保持为初始值。经过四组频率的反演之后，得到第一步速度反演的结果。

将第一步反演的速度结果和初始Q值模型，作为第二步反演的初始模型，本步骤对速度模型与Q值模型同时更新，得到最终的反演结果。

为了说明两步法反演的优化效果，我们将一步法反演策略得到的反演结果(图15(a)、(b))与两步法反演策略得到的结果(第一步更新结果如图16所示和第二步更新结果速度和Q值分别如图17(a)、(b)所示)进行对比，同时从模型误差的角度，精确的对比了两种反演策略模型误差。

从以上反演结果和反演模型误差对比(图18)中可以明显看出，相对于同时反演策略，两步法反演策略能够得到更好的联合反演结果，无论是速度参数还是Q值参数的反演曲线更加接近真实值，其模型误差也相较于一步法有所减小。其中对Q值结果的改善也是有效的。基本上起到了对联合反演的优化效果。

频率域全波形反演的过程中，一般都需要通过多次迭代过程才能最终获得最优解，尤其是数据量大的时造成全波形反演的计算负担。为了减少计算量，一般通过组合若干个震源建立超级炮的方法(Krebs et al.,2009)。但是，这种方法存在缺陷，由于同时激发的震源产生的波场在传播的过程中会引起串扰噪音，产生的串扰噪音，会使全波形反演问题性态变差，严重会导致目标函数无法收敛。

为了能使用“多源”的方法提高全波形反演的效率，需要寻求消除串扰噪音的方法。本专利使用的方法是通过在稀疏Seislet域中对参数模型施加L1范数正则化来提高全波形反演在噪音干扰下的稳定性。与其他的变换域相比，Seislet域展现更好的数据压缩能力和稀疏性。Seislet变换运用提升算法，有预测和更新两个步骤，对地震数据进行分解，可以定义如下

1)将输入数据当作一个记录序列，将数据分解成偶数分量和奇数分量e和o；

2)根据偶数分量对奇数分量进行预测，把预测的奇数分量和对应的奇数分量进行相减，可以得到残差序列r

r＝o-P[e]

残差序列，能反映原信号的高频信息，通常被称为小波系数，P为预测算子。

3)对数据的偶数分量进行更新，找到原始输入数据的近似序列c

c＝e+U[r]

其中c为尺度系数，反映原信号的低频信息，U为更新算子。

4)把近似序列c作为新数据，重复以上的步骤，再进行下个尺度运算，最终完成变换的多尺度分解过程。

由上述步骤可知，设计有效的变换关键是要得到预测算子P，使得残差足够小，而更新算子U保留原来数据的基本特征，同时进入下一个尺度水平。

Seislet变换重构是运用最高尺度系数c和全部尺度系数下的小波系数r，Seislet变换分解的逆过程步骤可以表示为如下

1)寻找最高尺度的尺度系数c以及所对应的残差r；

2)根据更新算子U，通过上式的逆运算来重构偶数分量，可以表示为如下式子

e＝c-U[r]

3)根据预测算子P，通过残差序列公式的逆运算来重构奇数分量，可以表示为如下式子

o＝r+P[e]

4)重复上述步骤，结合奇数和偶数分量以生成一个与前一个数据同一个粗度

Seislet变换是地震局部倾角与小波变换的提升算法相结合，而在上述的Seislet变换的分解以及重构中预测算子以及更新算子是Seislet的核心，因此，这种结合也就体现在预测算子P以及更新算子U上了，P，U算子可以定义为如下

其中，和分别表示沿着同相轴局部倾角方向对其左边以及右边地震道的预测算子。预测需要在不同的尺度上进行，在这种情况下，这意味着道之间具有不同的间隔距离。预测算子P和更新U算子沿着地震同相轴局部倾角的方向进行预测和更新。

离散的小波变换方向是沿着水平进行的，而Seislet变换是沿着地震同相轴局部倾角方向进行的，因此，在运用Seislet变换是就需要知道局部倾角信息，局部倾角信息可通过平面波分解得到。所以Seislet变换不像是小波变换只沿着水平方向进行，而是沿着同相轴的局部倾角方向进行，更有利于对地质构造特征的分析。如果同相轴的局的斜率都为零，Seislet变换就相当于小波变换提升算法。

在频率域全波形反演过程中，每一次迭代的过程中都涉及到模型更新，不管是速度模型还是Q值模型都可以写成如下形式：

m_k+1＝m_k+α_kp_k

为了对全波形反演进行稀疏约束，我们在模型更新的过程中，采用非线性整形正则化的方式，将更新迭代的过程更改为如下形式：

m_k+1＝S^-1TS(m_k+α_kp_k)

其中T表示软阈值算子，S和S^-1分别表示Seislet正、逆变换。这里面包含了整形算子，整形算子是由Seislet正、逆变化与软阈值算子的组合。整形算子的目的是利用Seislet域中参数模型的稀疏性，在每次迭代的过程中去除更新模型中的伪影。具体过程是首先将更新模型进行Seislet正变换，变换到Seislet域，这时将包含噪音的模型施加软阈值约束。再将Seislet域中的模型数据逆变换回空间域。根据以上理论，我们针对“多源”全波形反演，使用Seislet域稀疏约束条件，达到既能提高全波形反演的效率，又能消除“多源”全波形反演带来的串扰噪音。

选取Marmousi模型对以上理论进行试算，模型大小为512*188，网格点大小为16m*16m，初始模型为真实模型平滑结果,真实模型和初始模型如图19(a)、(b)所示，其中炮点位于深度90m处的地表，每次同时激发8炮，每炮相距1280m，反演的频率为一组，1～10hz，频域间隔为1hz，每个频带迭代10次，图20(a)、(b)所示分别为正演1hz和10hz的波场快照图，以及使用无约束方法的全波形反演结果和使用Seislet稀疏约束下的反演结果图(图21(a)、(b))。

通过之前建立的模型误差公式计算得出，普通全波形反演的模型误差结果为0.1135，而Seislet稀疏约束下的全波形反演模型误差为0.0941，Seislet稀疏约束下的全波形反演结果相对于普通全波形反演能得到更好的反演结果，能得到相对更准确的结果。多源同时激发带来的串扰噪音，在无Seislet稀疏约束的情况下，反演的结果出现很多噪点，同时容易使得全波形反演过程病态化。而Seislet稀疏约束下的全波形反演得到的反演结果，速度相对平滑，减轻了串扰噪音的影响，同时也增加了反演的稳定性。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (4)

1.一种频率域粘声介质全波形反演方法，其特征在于，包括：采用截断牛顿法对速度与Q值进行联合反演。

2.根据权利要求1所述的方法，其特征在于，所述截断牛顿法实现的过程分为内外两个循环结构，外部为传统的模型更新形式，内部的循环为共轭梯度形式的线性问题。

3.根据权利要求2所述的方法，其特征在于，所述联合反演包括：

首先反演速度，将反演的速度结果和初始Q值模型作为第二步反演的初始模型；

其次，对速度模型与Q值模型同时更新，得到最终的反演结果。

4.根据权利要求2所述的方法，其特征在于，使用稀疏约束方法对速度与Q值进行联合反演。