CN106662664A - 快速粘声波和粘弹性全波场反演 - Google Patents

Abstract

一种方法，其包括：获得初始地球物理模型；利用粘声波或粘弹性波方程对前向波场建模；利用伴随粘声波或伴随粘弹性波方程对伴随波场建模，其中分别地，伴随粘声波波方程基于是压力和记忆变量二者的函数的辅助变量，或伴随粘弹性波方程基于应力和记忆变量的组合；基于前向波场的模型和伴随波场的模型的组合获得成本函数的梯度；以及使用成本函数的梯度来更新初始地球物理模型并获得更新的地球物理模型。

Description

快速粘声波和粘弹性全波场反演

相关申请的交叉参考

本申请要求2014年6月17日提交的名称为“FAST VISCOACOUSTIC ANDVISCOELASTIC FULL-WAVEFIELD INVERSION”的美国临时专利申请62/013,455的权益，所述美国临时专利申请的全部通过引用并入本文。

技术领域

本文描述的示例性实施例涉及地球物理勘探，并且更具体地涉及地震数据处理，所述地震数据处理包括使用全波场反演(FWI)技术从地震数据恢复粘声波(viscoacoustic)和/或粘弹性参数，并且还包括使用逆时偏移(RTM)对地下界面成像。

背景技术

该部分旨在介绍可以与本发明的示例性实施例相关联的技术的各个方面。相信该讨论有助于提供促进对本发明的特定方面的更好理解的框架。因此，应理解的是，该部分应从这个角度来阅读，而不一定作为对现有技术的承认。

粘声波和粘弹性波场建模

在时域中，固有衰减(吸收和分散(dispersion))利用应变和松弛函数之间的卷积算子来公式表达。此类卷积对于使用时间推进方法的大规模波传播模拟是计算上不可实现的。将实际的衰减并入到时域计算最先由Day和Minster(Day和Minster，1984)使用帕德(Padé)逼近来实现。后来，帕德逼近演变成复杂的有理形式的粘声波和粘弹性流变学模型，例如广义麦克斯韦尔(Maxwell)模型(GMB)和标准线性固体(SLS)模型。时域中的衰减建模方法是基于由Emmerich和Korn(1987)以及Carcione等人(1988)分别给出的GMB或SLS公式，并且使用基于松弛机构(mechanism)的流变学模型。

松弛机构是表示在窄频带中的粘滞效应的时域衰减模型的单元。多个松弛机构被组合以对在期望频带上的衰减建模，将附加状态变量和偏微分方程(PDE)引入到声学和弹性前向波方程。这些附加变量和方程在前向波模拟期间消耗相当大的量的计算时间和存储空间，并且在伴随模拟期间消耗至甚至更大程度。

仅逼近SLS公式的简化版本已并入到FWI中以便以减少的计算成本推断出频率不变的衰减，例如由Robertsson等人于1994年、Charara等人于2000年、Hestholm等人于2006年以及Royle于2011年提出。此类逼近的精确度不仅随降低的品质因子值(较强衰减)而退化，而且它们不能被用于对随频率变化的品质因子(quality factor)建模。另一方法是将松弛机构的数目限制为一个或两个，以便最小化SLS模型的复杂性，这导致在感兴趣的地震频带上的目标品质因子的不佳表示。例如，Charara(Charara 2006)建议将仅两个松弛机构用于SLS模型以减少计算梯度时的计算成本。它们的公式要求计算在伴随模拟中的记忆变量的空间导数，这在可以在其逼近SLS模型中使用的松弛机构的总数目上限制它们。并且，Bai(Bai等人，2012)提出一种用于使用单个机构SLS模型来补偿FWI的衰减效应的方法，其中假定频率不变的品质因子。

尽管仅SLS衰减模型已并入到时域粘声波和粘弹性FWI技术中，但是SLS和GMB模型两者已被用于波的前向粘声波和粘弹性建模。例如，等人(等人，2007)应用GMB以用于利用不连续Galerkin法来对粘弹性波建模。

衰减的常用量度是品质因子，即一种定义声学或弹性模量的频率依赖性的无量纲的量。品质因子本身可以是频率依赖的，尤其是对于含流体(fluid-bearing)岩石诸如烃(hydrocarbon)，并且对于干燥岩石，品质因子通常被假定为频率不变的(Muller等人，2010；Quintal，2012)。因此，推断品质因子的频率依赖性以识别诸如储层的含流体岩石是有用的。

SLS衰减模型通过以下关系式表示品质因子Q的倒数(也被称为损耗因子)(Carciones等人，1988)：

其中

Q＝品质因子，

τ_∈l＝SLS模型中的机构l的应变松弛时间，

τ_σl＝SLS模型中的机构l的应力松弛时间，

x＝空间坐标，

ω＝频率，

L＝在SLS模型中使用的松弛机构的数目。

概念地，品质因子Q表示介质中的储存能量与耗散能量的比率。应变和应力松弛时间被确定以在频带上最佳地配合期望的品质因子分布。

用于模拟衰减声学介质中的波的常规的一阶形式的线性粘声波波方程(JafarGandomi等人2007)是：

其具有用于压力p、速度v以及记忆变量m_l的适合的初始和边界条件。注意

＝散度算子，

κ＝非松弛体积模量(lim_ω→∞κ(ω)→κ)，

ρ＝质量密度，

v＝速度(在三维空间中v＝{v_x v_y v_z}^T)，

p＝压力，

m_l＝机构l的记忆变量，

s_p＝压力源，

s_v＝速度源，

且其中松弛参数τ_∈l和τ_σl可以通过针对给定品质因子分布的方程(1)来确定。

注意，贯穿该文档，连续的标量变量由斜体字符表示，并且向量和矩阵由粗体非斜体字符表示。

基于计算目标函数关于参数的梯度的FWI方法通常通过使用伴随方法来有效地实施，所述伴随方法已被证明胜过其它相关方法，如直接灵敏度分析、有限差分或复变量方法。我们可以在Thevenin等人，2008中找到关于伴随方法的详细讨论。

常规粘声波系统(方程(2))的连续伴随是

其中

＝伴随压力，

＝伴随速度，

＝机构l的伴随记忆变量，并且

和分别是目标函数关于压力和速度的导数。尽管目标函数及其关于压力和速度的导数将在粘声波和粘弹性全波场反演部分中明确定义，但是和可以简单地被认为是伴随方程的源。

注意，伴随方程的求导(derivation)涉及使用分部积分法，所述分部积分法将空间和时间边界项引入到伴随方程(其分别在空间和时间边界处被求值(evaluated))。这些伴随空间边界项并不包含在方程(3)中，但是需要被包含(即，以代码形式)以便正确地计算伴随变量的唯一解，并且时间边界项是0并因此丢弃。

伴随方程(方程(3))中的记忆变量的空间导数导致大量的计算成本。图5示出用于二至十二阶(空间)准确有限差分(FD)时域方法的伴随计算相对于前向计算的成本。相对成本基于比较每秒浮点运算次数(flops)的数量。对于三个松弛机构和八阶有限差分方法，对伴随方程求解的成本是对前向方程(方程(2))求解的成本的2.25倍。还应注意，当有限差分方法的空间阶数和松弛机构的数目增加时，对伴随方程求解的相对成本增加。此外，期望的是，对于大规模分布式存储器并行处理计算，求解伴随的成本将由于记忆变量的附加通信而增加。

线性粘声波方程可以扩展到以下线性粘弹性方程：

其中在三维笛卡尔坐标系中，

σ＝应力，σ＝{σ_xx σ_yy σ_zz σ_xz σ_yz σ_xy}^T，

m_l＝机构l的记忆变量，

m_l＝{m_l，xx m_l，yy m_l，zz m_l，xz m_l，yz m_l，xy}^T，

这是各向同性非松弛系统关于Lamé常量λ和μ的弹性本构关系，

这是记忆系统的本构关系(对于SLS公式

并且和使用方程(1)从压缩和剪切波品质因子Q^P和Q^S计算出，

＝应变算子，

常规粘弹性系统(方程(4))的连续伴随是

其中在三维笛卡尔坐标系中，

＝伴随应力，

＝机构l的伴随记忆场，

并且和是目标函数分别关于应力和速度的导数。

注意，空间边界项未被包含在方程(4)中。时间边界项是零并且因此丢弃。

关于粘声波，用于粘弹性情况的伴随方程还包含记忆变量的空间导数，这将给粘弹性FWI增加大的计算成本。

粘声波和粘弹性全波场反演

FWI是迭代地最小化测量波场与计算波场之间的失配(misfit)的范数的偏微分方程约束优化方法。地震FWI涉及多次迭代，并且单个迭代可以涉及以下计算：(1)前向方程的求解，(2)伴随方程的求解，以及(3)这些前向解和伴随解的卷积，以产生成本函数的梯度。注意，对于二阶优化方法，例如高斯-牛顿法(Gauss-Newton)，还需要(4)扰动前向方程的求解。针对该情况的更稳健的数学论证可以例如在美国专利公开20130238246中找到，所述公开的全部内容通过引用并入本文。

对于粘声波和粘弹性地震FWI，每个模拟的成本由于以下事实而加重：对品质因子的频率依赖分布准确地建模所需的松弛机构引入更多的变量和方程。此外，必须计算的模拟的数目与在反演中的迭代的总数目成比例，所述迭代的总数目通常约为数百到数千。然而，使用该方法推断地下的衰减特性的益处被期望超过成本，并且引起较快周转时间的算法和工作流的发展是朝向使得该技术对于场规模数据可行的关键步骤，从而允许使用者更快地解决更大规模问题。

在地球物理学中使用的常见迭代反演方法是成本函数优化。成本函数优化包括成本函数关于模型θ的值的迭代最小化或最大化。成本函数(也被称为目标函数)是模拟数据和观察数据之间的失配的测量。模拟(模拟数据)通过以下操作来进行：首先用适合的数值方法(例如有限差分或有限元方法)来离散化介质中的源信号的物理管理传播，并且使用当前的地球物理特性模型来在计算机上计算数值解。

下文概述了FWI的局部成本函数优化过程：(1)选择开始模型；(2)计算搜索方向S(θ)；以及(3)搜索更新模型，其为在搜索方向上的模型的扰动。

成本函数优化过程通过将新的更新模型用作用于寻找另一搜索方向的开始模型来迭代，然后将其用于扰动模型以便更好地解释观察数据。所述过程继续，直到找到令人满意地解释观察数据的更新模型。常用的局部成本函数优化方法包含梯度搜索、共轭梯度、准牛顿、高斯-牛顿以及牛顿的方法。

声学逼近中的地震数据的局部成本函数优化是常见的地球物理反演任务，并且通常说明其它类型的地球物理反演。当反演声学逼近中的地震数据时，成本函数可以写作：

其中

＝成本函数，

θ＝描述地下模型的N个参数(θ₁，θ₂，...，θ_N)的向量，

g＝道集索引，

w_g＝道集g的源函数，其为空间坐标和时间的函数，对于点源而言，这是空间坐标的δ函数，

Ng＝道集的数目，

r＝道集内的接收器索引，

N_r＝道集中的接收器的数目，

t＝轨迹内的时间样本索引，

N_t＝时间样本的数目，

W＝范数函数(最小化函数，例如，用于最小二乘函数(x)＝x²)，

ψ_calc＝来自模型θ的计算的地震数据，

ψ_obs＝测量的地震数据(压力、应力、速度和/或加速度)。

道集(来自共享公共几何结构的多个传感器的数据)可以是可以在地震前向建模程序的一次运行中被模拟的任何类型的道集(共中心点、共源、共偏移、共接收器等)。通常，道集对应于地震激发(shot)，但是激发可以比点源更普遍。对于点源，道集索引g对应于独立点源的位置。该广义源数据ψ_obs可以在场中获得或可以从使用点源获得数据合成。在另一方面，计算的数据ψcalc通常可以通过在前向建模时使用广义源函数来直接地计算。

FWI试图更新离散化模型θ，使得是最小值。这可以通过局部成本函数优化来实现，所述局部成本函数优化如下更新给定模型θ^(k)：

θ⁽ⁱ⁺¹⁾＝θ⁽ⁱ⁾+γ⁽ⁱ⁾S(θ⁽ⁱ⁾)， (7)

其中i是迭代次数，γ是模型更新的标量步长，并且S(θ)是搜索方向。对于最陡下降，这是关于模型参数所采用的失配函数的梯度的负数。在此情况下，模型扰动或通过其模型被更新的值通过将目标函数的梯度与步长γ相乘来计算，这必须反复计算。对于二阶优化技术，梯度通过海赛(Hessian)函数(目标函数关于模型参数的二阶导数)来缩放。的计算需要关于N个模型参数中的每一个的导数的计算。N在地球物理学问题中通常非常大(超过一百万)，并且如果必须针对每个单独模型参数执行该计算，则其可能是非常耗时的。幸运地是，可以使用伴随方法来有效地针对所有模型参数一次执行该计算(Tarantola，1984)。尽管使用伴随方法的梯度的计算相对于其它方法是有效的，但该方法对于粘声波和粘弹性FWI成本仍然非常高。

时域中的常规的粘声波和粘弹性全波场反演方法计算记忆变量的梯度以时间推进伴随方程。已经提出在频域和时域两者中的考虑衰减的地球模型(Ursin和Toverud，2002)。频域方法的主要焦点已经是在因果原理下使用复数特性在衰减和介质速度之间建立关系，这使得它们是简单的以直接地应用于FWI(Hak和Mulder，2010)。

发明内容

一种方法，其包含：获得初始地球物理模型；利用粘声波或粘弹性波方程对前向波场建模；利用伴随粘声波或伴随粘弹性波方程对伴随波场建模，其中分别地，伴随粘声波波方程基于为压力和记忆变量的函数的辅助变量，或伴随粘弹性波方程基于应力和记忆变量的组合；基于前向波场的模型和伴随波场的模型的组合来获得成本函数的梯度；以及使用成本函数的梯度来更新初始地球物理模型并获得更新的地球物理模型。

附图说明

尽管本发明容许各种修改和替代形式，但是本发明的特定示例实施例已经在附图中示出并且在本文中详细描述。然而，应理解，本文中对特定示例实施例的描述不旨在将本公开限制到本文所公开的特定形式，本公开将覆盖如由随附权利要求书限定的所有修改和等同物。还应理解，附图未必是按比例的，而是将重点放在清楚地图示说明本发明的示例性实施例的原理上。此外，某些尺寸可以被放大以帮助在视觉上传达此类原理。

图1示出用于针对频率不变Q＝30(频带从2.5到40Hz)而优化的广义麦克斯韦尔体模型的基于1、2、3和5个松弛机构的示例性品质因子模型。

图2示出针对图1的广义麦克斯韦尔体模型的利用1、2、3和5个松弛机构计算的示例性速度分散曲线。

图3示出用于针对在从2.5到40Hz的频带上线性变化的品质因子而优化的广义麦克斯韦尔体模型的基于1、2、3和5个松弛机构的示例性品质因子模型。

图4示出用于图3的广义麦克斯韦尔体模型的利用1、2、3和5个松弛机构计算的示例性速度分散曲线。

图5示出用于针对L＝3和L＝5个松弛机构的常规三维粘声波的伴随相对于前向的成本(按照浮点运算-每秒浮点运算次数)的非限制性示例。

图6示出用于粘声波全波场反演的示例性方法。

图7示出用于粘弹性全波场反演的示例性方法。

图8是根据本技术进步对三维粘声波前向和伴随方程求解的成本相对于利用三个松弛机构衰减模型对常规的声学前向方程求解的成本的非限制性示例。

图9是根据本技术进步的三维粘声波伴随相对于针对L＝2、L＝3以及L＝5个松弛机构衰减模型的常规伴随的加速的非限制性示例。

图10是Marmousi速度模型的示例(标度单位为km/s)。

图11是真实的品质因子模型的示例(标度单位是针对品质因子的倒数)。

图12是反演速度模型的示例(标度单位为km/s)。

图13是反演品质因子模型的示例(标度单位是针对品质因子的倒数)。

图14是目标函数随迭代数目的变化(数据失配的最小二乘范数)的示例。

图15是可以执行实施本技术进步的方法的计算机的示例。

具体实施方式

本文描述示例性实施例。然而，就以下描述特定于特定实施例来说，这旨在于仅用于示例性目的并且简单地提供对示例性实施例的描述。因此，本发明不限于下文描述的特定实施例，相反，本发明包含落入随附权利要求书的真正精神和范围内的所有替代、修改以及等同物。

本技术进步的示例性实施例是一种用于执行需要粘声波或粘弹性波传播的全波场反演(FWI)或逆时偏移(RTM)的方法。所述方法包含使用多松弛流变学模型在时域中的衰减，所述多松弛流变学模型最优地表示针对给定数目的松弛机构的衰减。该衰减模型灵活用于表示在给定地震频率带宽上的品质因子的常量或变化的频率依赖性。该方法使用新辅助状态变量来用公式表达粘声波和/或粘弹性波方程，所述新辅助状态变量是针对粘声波情况的压力和记忆变量的组合或针对粘弹性情况的应力和记忆变量的组合。这些公式导致具有降低的计算复杂度的新的波方程系统。

粘声波和粘弹性FWI框架的示例性实施例依赖于GMB流变学模型来解释衰减，所述GMB流变学模型在频域ω中将损耗因子定义为

其中

φ_l＝机构l的松弛频率(φ＝{φ₁，φ₂，...，φ_L})，

α_l＝机构的松弛模量(moduli)比率((α＝{α₁，α₂，...，α_L}))。

松弛频率通常选择为在频带上指数地变化并且松弛频率是空间不变的(Robertsson等人，1994)。松弛模量通过以下方法确定。方程的超定系统通过在源的频带ω＝{ω₁，ω₂，...，ω_n}上离散化方程(8)来形成，其中n是在频带上离散化的频率的数目并且大于记忆变量的数目L，并且针对几乎频率不变的品质因子n通常为2L+1并且针对频率变化的品质因子更大。频率通常在频带上呈指数性离散化以实现计算有效性。方程的该超定系统是针对在坐标x处的离散化品质因子Q＝{Q(ω₁，x)，Q(ω₂，x)，...，Q(ω_n，x)}：

A(ω，φ，Q^-1(ω))α＝Q^-1(ω)， (9)

其中

离散化品质因子关系(7)并且A具有L乘n的维度。

α可以是反演参数。然而，如果品质因子Q在频带上是恒定的，那么所

述品质因子可以替代α直接被认为是反演参数。在此类情况下，目标函数

关于品质因子Q的导数可以使用以下链式法则来获得

其中松弛模量关于Q的导数通过对以下方程求解来获得

图1和2示出通过GMB模型的1、2、3和5个松弛机构获得的品质因子和速度分散曲线的示例，所述GMB模型被设计成表示在从2.5到40Hz的频带上的频率不变的品质因子30(目标Q)。这些衰减模型的参数使用方程(9)来计算。这些图示出建模的衰减和分散的精确度随着松弛机构的数目而增加。图3和图4示出这些模型针对线性变化的品质因子(在频带上)的品质因子和分散响应的示例，并且再次示出建模的衰减和分散的精确度随着松弛机构的数目而增加。

我们导出(derive)以连续形式的粘声波和粘弹性伴随方程(通常被称为连续伴随方法)，以便证明本技术进步的提高的效率。然而，我们可以首先离散化前向方程，并且然后导出离散化方程的伴随(Betts等人，2005)，这被称为离散伴随方法。在该申请中描述的本技术进步适用于FWI方程的两种类型的求导。

本技术进步的一个实施例提供改进的粘声波FWI方法，其基于新形式的粘声波波方程。新的辅助变量定义为

且方程(12)用于执行方程(2)的线性变换。所得的粘声波方程是

其中r是出于清楚而定义的新变量。记忆方程中使用的松弛模量比率从方程(9)获得。

方程(13)的伴随为

注意，空间边界项不包含在方程(14)中。时间边界项是零且因此丢弃。

(方程(12))引入的辅助变量(所述辅助变量为压力和记忆变量的组合)取代记忆变量，以便相比于不采用辅助变量的替代方程而降低对伴随波场方程求解的成本。

与常规粘声波伴随(方程(3))相比，记忆变量的空间导数不存在于伴随方程(14)中，这降低对伴随求解的计算成本。在实施本技术进步的示例性方法中，对伴随(方程(14))求解的成本与对前向方程(方程(13))求解的成本近似相同。图8示出对常规的和新的前向和伴随粘声波方程求解的成本相对对前向声学方程求解的成本(基于比较每秒浮点运算次数的数量)。在图8中，实线描绘对常规方程求解的成本，并且虚线描绘对实施本技术进步的方程求解的成本。即使新前向求解程序的成本实际上增加了较小的边际(margin)(相对于常规的前向求解程序)，但是新伴随求解程序的成本相对于现有技术伴随求解程序大大降低，并且因此总成本降低。例如，对于八阶FD时域方法和使用三个松弛机构的流变学模型，用于现有技术方法的相对伴随成本是2.68，并且针对新方法降低至1.24。此外，在图9中示出针对2、3以及5个松弛机构的当前进步相对于现有技术的伴随加速。

目标函数关于反演参数的连续梯度是

如果Q是反演参数的话，

因此，

二阶优化方法所需的扰动前向方程是

其中

＝扰动压力，

＝扰动速度，

＝机构l的扰动记忆变量。

同一方法可以扩展到粘弹性波方程。对于给定数目的松弛机构，计算节省将较大，因为三维粘弹性方程的记忆变量的总数将是粘声波情况的记忆变量的总数的六倍。我们定义新的辅助变量

并且利用新的辅助变量执行方程(4)的线性变换。该变换消去伴随中的记忆变量的空间导数，并且降低用于求解伴随系统的计算成本。所得的粘弹性波方程是

其中r是为了清楚而定义的新变量。记忆方程中使用的松弛模量比率从方程(9)获得。

方程(18)的伴随为

其中

＝机构l的伴随记忆场，

＝为了清楚为机构l定义的变量，

注意，空间边界项未被包含在方程(19)中。时间边界项是零并且因此丢弃。

与该申请中示出的一阶时间导数公式相反，直接导出声波或弹性波方程的替代的计算上有效的二阶时间导数公式。对于声学情况，可以通过取(13)中的第一(压力)和第三(记忆)方程两者的时间导数并且使用第二(速度)方程消去项来导出二阶形式的声学方程。对于弹性情况，可以通过取(18)中的第一(应力)和第三(记忆)方程两者的时间导数并且使用第二(速度)方程消去项来导出二阶形式的弹性方程。记忆变量的空间导数将不存在于如针对一阶系统(方程14和19)所证明的此类二阶系统的伴随方程中。

示例

图6示出用于使用方程(13)到(16)的粘声波FWI方法的示例性方法。在步骤601中，假定地球物理模型。地球物理模型依据区域中的位置给出一个或多个地下特性。对于该示例，使用Marmousi速度模型(图10)和具有若干低Q异常、中间Q值层以及平滑背景的品质因子模型(图11)。速度和品质因子模型在结构上是不相关的。

在步骤603中，使用方程(9)来计算衰减参数α。GMB模型可以被用于假定的Q频率依赖性。

在步骤605a和605b中，对前向波场模型方程(13)和伴随模型方程(14)求解。步骤605b可以将测量的数据用于建模伴随波场。

在步骤607中，成本函数的梯度从方程(13)和方程(14)的卷积获得，以便获得目标函数关于反演参数(14)的梯度。步骤607可以将测量的数据用于获得成本函数的梯度。

在步骤609中，成本函数的梯度(其提供成本函数在给定方向上的变化率)随后被用于更新地球物理模型，以便最小化成本函数。步骤609可以包括搜索更新的地球物理特性模型，该更新的地球物理特性模型是在更好地解释观察数据的梯度方向上的初始地球物理特性的扰动。图6的迭代过程可以重复，直到获得在测量的数据与更新的模型之间的预定收敛。图14示出目标函数的数据失配(数据失配的最小二乘范数)如何随迭代的数目而减少。

假定16km最大偏移和10秒的记录(合成)数据、40个源、100个接收器以及10HzRicker子波源，图12和图13中分别示出速度和Q模型的反演结果或更新结果，证明本技术进步能够以大大降低的计算成本来针对两个参数(速度和Q)进行反演。

当更新的假定模型收敛时，过程前进到步骤611。在步骤611中，使用更新的地下模型来管理烃。如本文中所使用，烃管理包含烃提取，烃产生，烃勘探，识别潜在的烃资源，识别井位置，确定井注入和/或提取率，识别储层连通性，获取、处置和/或放弃烃资源，回顾先前的烃管理决策，以及任何其它烃相关的动作或活动。

图7示出用于使用方程(18)到(19)的粘弹性FWI方法的示例性方法。除了方程已经改变之外，图7的方法(步骤701到711)类似于图6的方法(步骤601到611)。

图15是可以用于执行当前的地球物理反演技术的计算机系统2400的框图。中央处理单元(CPU)2402耦合到系统总线2404。CPU 2402可以是任何通用CPU，但是只要CPU 2402(以及系统2400的其它组件)支持如本文中描述的操作，就也可以使用其它类型架构的CPU2402(或示例性系统2400的其它组件)。本领域技术人员应认识到，尽管在图15中仅示出单个CPU2402，但是可以存在附加的CPU。此外，计算机系统2400可以包括网络化多处理器计算机系统，所述网络化多处理器计算机系统可以包括混合并行CPU/GPU系统。CPU 402可以根据本文中公开的各种教导来执行各种逻辑指令。例如，CPU 2402可以执行机器级指令以用于执行根据描述的操作流的处理。

计算机系统2400还可以包括计算机组件，诸如非暂时性计算机可读介质。计算机可读介质的示例包括随机存取存储器(RAM)2406，其可以为SRAM、DRAM、SDRAM等。计算机系统2400还可以包括附加的非暂时性计算机可读介质诸如只读存储器(ROM)2408，其可以为PROM，EPROM，EEPROM等。如本领域中已知，RAM 2406和ROM 2408保持用户和系统数据和程序。计算机系统2400还可以包括输入/输出(I/O)适配器2410、通信适配器2422、用户接口适配器2424以及显示器适配器2418。

I/O适配器2410可以将诸如(多个)存储装置2412的附加的非暂时性计算机可读介质连接到计算机系统2400，所述(多个)存储装置2412包括例如硬盘驱动器、光盘(CD)驱动器、软盘驱动器、磁带驱动器以及类似物。当RAM 2406不足以用于与存储用于本技术的操作的数据相关联的存储器需求时，可以使用(多个)存储装置。计算机系统2400的数据存储装置可以用于存储如本文中公开的所使用或产生的信息和/或其它数据。例如，(多个)存储装置2412可以用于存储根据本技术的配置信息或附加插件。此外，用户接口适配器2424将用户输入装置诸如键盘2428、定位装置2426和/或输出装置耦合到计算机系统400。显示器适配器2418由CPU 2402驱动以控制在显示器装置2420上的显示，例如，向用户呈现关于可用插件的信息。

系统2400的架构可以根据需要变化。例如，可以使用任何合适的基于处理器的装置，包含但不限于个人计算机、膝上型计算机、计算机工作站以及多处理器服务器。此外，本技术进步可以被实施在专用集成电路(ASIC)或超大规模集成(VLSI)电路上。事实上，本领域技术人员可以使用能够执行根据本技术进步的逻辑操作的任何数目的合适硬件结构。术语“处理电路”涵盖硬件处理器(例如在上述硬件装置中发现的那些)、ASIC以及VLSI电路。计算机系统2400的输入数据可以包括各种插件和库文件。输入数据可以附加地包括配置信息。

本技术可以容许各种修改和替代形式，并且以上讨论的示例已仅通过示例的方式示出。然而，本技术不旨在被限制到本文公开的特定示例。事实上，本技术包含落入随附权利要求书的精神和范围内的所有替代、修改以及等同物。

参考文献

以下文件通过引起以其全部并入本文。

1.Tarantola，A.，″Inversion of seismic reflection data in the acousticapproximation，″Geophysics 49，1259-1266(1984).

2.Hak，B.，and Mulder W.A.″Seismic attenuation imaging with causality，″Geophysical Journal International 184，439-451(2011).

3.Day，S.M.and Minster J.B.，“Numerical simulation of attenuatedwavefields using a Padé approximant method，”Geophysical Journal ofInternational，78，105-118(1984).

4.Ursin，B.and Toverud，T.，“Comparison of seismic dispersion andattenuation models，”Studia Geophysica et Geodaetica，46，293-320(2002).

5.Carcione，J.M.，Kosloff D.and Kosloff R.，“Viscoacoustic wavepropagation simulation in the Earth，”Geophysics，53，769-777(1988).

6.Carcione J.M.“Wave fields in real media：Wave propagation inanisotropic anelastic，porous and electromagnetic media”Handbook ofGeophysical Exploration：Seismic Exploration，38，Elsevier，2^nd Edition(2007)

7.Muller，T.M.，Gurevich，B.and Lebedev，M.，“Seismic wave attenuation anddispersion resulting from wave induced flow in porous rocks-A review，”Geophysics，75，147-164(2010).

8.Bai，J.，Yingst D.，Bloor R.and Leveille J.“Waveform inversion withattenuation，”Society of Exploration Geophysicists Extended TechnicalAbstract，2012

9.Royle，G.T.，“Viscoelastic orthorhombic full wavefield inversion：development of multiparameter inversion methods，”Society of GeophysicistExtended Abstact，2011.

10.Robertsson，J.A.O.，Blanch J.O.and Symes W.W.，“Viscoelastic finite-difference modeling，”Geophysics，59，1444-1456(1994).

11.Hestholm，S.，Ketcham，S.，Greenfield，R.，Moran，M.and McMechan G.，“Quick and accurate Q parameterization in viscoelastic wave modeling，”Geophysics，71，147-150(2006).

12.Charara，M.，Barnes，C.and Tarantola，A.，“Full waveform inversion ofseismic data for a viscoelastic medium”Methods and Applications of Inversion：Lecture Notes in Earth Sciences，92，68-81(2000).

13.M.，Dumbser，M.，Puente，J.and Igel，H.“An arbitrary high-orderDiscontinuous Galerkin method for elastic waves on unstructured meshes-III.Viscoelastic attenuation，”Geophysics Journal International，168，224-242(2007).

14.JafarGandomi，A.and Takenaka，H.，“Efficient FDTD algorithm for planewave simulation for vertically heterogeneous attenuative media，”Geophysics，72，43-53(2007).

15.Quintal，B.，“Frequency-dependent attenuation as a potentialindicator of oil saturation，”Journal of Applied Geophysics，82，119-128，2012.

16.Betts，J.T.and Campbell，S.L.，“Discretize then optimize，”Mathematicsin Industry：Challenges and Frontiers A Process View：Practice and Theory，Ferguson，D.R.and Peters，T.J.，eds.，SIAM Publications(2005).

17.Thevenin D.and Janiga G.，“Optimization and computational fluiddynamics，”Springer-Verlag，(2008).

Claims (14)

1.一种方法，其包括：

利用计算机获得初始地球物理模型；

利用所述计算机使用粘声波或粘弹性波方程对前向波场建模；

利用所述计算机使用伴随粘声波或伴随粘弹性波方程对伴随波场建模；

其中分别地，所述伴随粘声波波方程基于是压力和记忆变量的函数的辅助变量，或所述伴随粘弹性波方程基于应力和记忆变量的组合；

利用所述计算机基于所述前向波场的模型和所述伴随波场的模型的组合来获得成本函数的梯度；以及

利用所述计算机使用所述成本函数的所述梯度来更新所述初始地球物理模型并获得更新的地球物理模型。

2.根据权利要求1所述的方法，其进一步包括使用所述辅助变量以使记忆变量的空间导数从所述伴随粘声波或粘弹性波方程中退出。

3.根据权利要求1或2所述的方法，其中获得所述地球物理模型包括使用初始流变学模型和声学或弹性介质参数来获得针对给定品质因子的松弛参数。

4.根据权利要求1至3中的任一项所述的方法，其中使用所述伴随粘声波方程建模所述伴随波场。

5.根据权利要求1至4中的任一项所述的方法，其中所述辅助变量是其中p是压力，m_l是机构l的所述记忆变量，其中l是大于或等于1的整数，并且α_l是机构l的反演参数。

6.根据权利要求1至5中的任一项所述的方法，其中所述伴随粘声波方程是

$\frac{\∂{\stackrel{\‾}{\mathrm{\ζ}}}_l}{\∂t}+{\mathrm{\φ}}_1{\stackrel{\‾}{\mathrm{\ζ}}}_l-{\mathrm{\φ}}_1{\stackrel{\‾}{r}}_l=0,$

其中

是散度算子，

t是时间，

κ是非松弛体积模量(lim_ω→∞κ(ω)→κ)，

ρ是质量密度，

v是速度(在三维空间中(v＝{v_x v_y v_z}^T)，

是伴随压力，

是成本函数，其中θ是模型参数(即，Q或α_l、κ以及ρ的某组合)

是伴随速度，以及

φ_l是机构l的松弛频率，

和分别是目标函数关于所述压力和速度的导数。

7.根据权利要求1至6中的任一项所述的方法，其中使用所述伴随粘弹性方程建模所述伴随波场。

8.根据权利要求7所述的方法，其中所述辅助变量是其中σ是应力，m_l是机构l的记忆变量，其中l是大于或等于1的整数，C是各向同性的非松弛系统的弹性本构关系并且D_l是记忆系统的本构关系。

9.根据权利要求8所述的方法，其中所述伴随粘弹性方程是

$\frac{\∂{\stackrel{\‾}{\mathrm{\ζ}}}_l}{\∂t}+{\mathrm{\φ}}_1{\stackrel{\‾}{\mathrm{\ζ}}}_l-{\mathrm{\φ}}_1{\stackrel{\‾}{r}}_l=0,$

其中

t是时间，

是伴随应力，

是伴随速度，

是机构l的伴随记忆场，

${\stackrel{\‾}{\mathrm{\ζ}}}_l={\left\{\begin{array}{cccccc}{\stackrel{\‾}{\mathrm{\ζ}}}_{l,\mathrm{xx}}& {\stackrel{\‾}{\mathrm{\ζ}}}_{l,\mathrm{yy}}& {\stackrel{\‾}{\mathrm{\ζ}}}_{l,\mathrm{zz}}& {\stackrel{\‾}{\mathrm{\ζ}}}_{l,\mathrm{xz}}& {\stackrel{\‾}{\mathrm{\ζ}}}_{l,\mathrm{yz}}& {\stackrel{\‾}{\mathrm{\ζ}}}_{l,\mathrm{xy}}\end{array}\right\}}^T,$

是为了清楚为机构l定义的变量，

${\stackrel{\‾}{r}}_l={\left\{\begin{array}{cccccc}{\stackrel{\‾}{r}}_{l,xx}& {\stackrel{\‾}{r}}_{l,yy}& {\stackrel{\‾}{r}}_{l,zz}& {\stackrel{\‾}{r}}_{l,xz}& {\stackrel{\‾}{r}}_{l,yz}& {\stackrel{\‾}{r}}_{l,xy}\end{array}\right\}}^T,$

ρ是质量密度，

和分别是所述目标函数关于所述应力和速度的导数，

这是所述各向同性非松弛系统关于Lamé常量λ和μ的所述弹性本构关系，

这是所述记忆系统的所述本构关系(使用方程(9)从压缩和剪切波品质因子来计算)，且应变算子，

是成本函数，其中θ是模型参数(即，Q^P、Q^S或λ、μ以及ρ的某组合)。

10.根据权利要求1至9中的任一项所述的方法，其进一步包括使用所述更新的地下模型管理烃。

11.根据权利要求1至10中的任一项所述的方法，其中所述地球物理模型包括速度模型。

12.根据权利要求1至11中的任一项所述的方法，其中所述地球物理模型包括品质因子模型。

13.根据权利要求1至12中的任一项所述的方法，其中所述辅助变量基于多个松弛机构。

14.根据权利要求1至13中的任一项所述的方法，其中所述粘声波或粘弹性波方程基于所述辅助变量。