CN103270430A - 用于地震数据反演的系统和方法 - Google Patents

Abstract

本发明公开了用于从地震数据确定感兴趣地下区域的属性的系统和计算机实现的方法。在一种实施例中，该方法包括：获得代表地下区域的地震数据和用于该地下区域的初始地球属性模型；使用该初始地球属性模型来执行正演模拟，以创建与实际地震数据具有相似采集规范的模拟的地震数据；把模拟的和实际的地震数据变换到时间傅立叶频域，以创建频域模拟的和实际地震数据，其中频域模拟的和实际地震数据包括振幅部分和相位部分；测量频域模拟的地震数据和频域实际地震数据之间的错配，以产生频域残余地震数据；执行频域残余地震数据中的某些观察到的频率分量的相位部分的相位展开，以创建展开的残余相位部分；反演展开的残余相位部分，以确定感兴趣地下区域的期望属性，其中，反演最小化被定义成测量错配的目标函数。此外，该方法还可以包括使用第一反演步骤的结果作为初始模型的第二反演步骤。所述系统包括数据源、用户接口和配置成执行设计成执行所述方法的计算机模块的处理器。

Description

用于地震数据反演的系统和方法

技术领域

本发明总体上涉及用于反演地震数据以便计算地球地下物理属性的方法和系统，尤其涉及用于执行唯相全波形反演以便从地震数据中计算速度模型的方法和系统。

背景技术

地下勘探，尤其是对烃储层的勘探，一般使用诸如地震数据迁移的方法来产生可解释的地球地下图像。在其中地下由于断层作用、盐体等而复杂的地区中，传统的迁移方法常常不能产生适当的图像。此外，传统的迁移方法需要地下的合理准确速度模型；这种速度模型也可从地震数据中确定，但是在专家和计算成本方面都非常昂贵。

存在许多从地震数据中计算速度模型的常规方法，包括NMO速度分析、迁移速度分析、层析X射线成像以及全波形反演。有些方法，例如全波形反演，是非常计算昂贵的，并且随着计算能力的增加，最近才变得可行。常规的全波形反演是在时域中或者在诸如时间傅立叶变换域或拉普拉斯变换域的变换域中执行的。这些方法常常由于地震数据中缺乏低频，一般小于3赫兹，而失败。如本领域技术人员将认识到的，速度模型是低频模型，因此很难从缺乏低频信息的地震数据为其反演。

确定速度模型并把它们用于迁移以产生地球地下图像的传统方法是昂贵的并充满了困难，尤其是在复杂地区中。当对烃的搜索移到这些复杂地区时，需要找到更好的途径来处理地震数据并改进速度模型。

发明内容

根据本发明的一种实现，公开了一种确定地下区域的属性的计算机实现方法。该方法包括：获得代表地下区域的实际地震数据和用于该地下区域的初始地球属性模型；利用该初始地球属性模型执行正演模拟以创建与实际地震数据具有相似采集规格的模拟的地震数据；把模拟的地震数据和实际地震数据变换到时间傅立叶频域，以创建频域模拟的地震数据和实际地震数据，其中，频域模拟的地震数据和实际地震数据包括振幅部分和相位部分；测量频域模拟的地震数据和频域实际地震数据之间的错配，以产生频域残余地震数据；对频域残余地震数据中的某些观察到的频率分量的相位部分执行相位展开，以创建展开的残余相位部分；以及反演展开的残余相位部分，以确定感兴趣地下区域的期望属性，其中，所述反演最小化被定义成测量所述错配的目标函数。

该方法进一步包括相位外推。此外，该方法可以利用从第一反演步骤确定的属性作为初始模型来执行第二反演步骤。

在另一种实施例中，公开了一种用于确定地下区域的属性的系统。该系统包括数据源、用户接口和配置成执行被设计成执行所述方法的计算机模块的处理器。

在另一种实施例中，公开了一种用于确定地下区域的属性的制品。该制品包括其中并入计算机可读代码的计算机可读介质，所述代码被配置成执行所述方法。

以上概述部分的提供是为了以简化形式介绍在以下具体描述部分中进一步描述的概念的选择。本概述不是要识别所保护主题的关键特征或基本特征，也不是试图用于限定所保护主题的范围。此外，所保护的主题不限于解决本公开内容任何部分中所指出的任何或全部缺陷的实现。

附图说明

关于以下具体描述、未决的权利要求及附图，本发明的这些及其它特征将变得更好理解，附图中：

图1是说明全波形反演方法的流程图；

图2说明了各个频率处的梯度带宽；

图3说明了从好的初始地球属性模型开始的常规全波形反演过程；

图4说明了从差的初始地球属性模型开始的常规全波形反演过程；

图5是说明根据本发明实施例的方法的流程图；

图6说明了在非常低的频率有和没有预处理器的相位展开方法；

图7说明了在适中的低频率有和没有预处理器的相位展开方法；

图8说明了唯相全波形反演的实施例的结果；

图9说明了唯相全波形反演跟随有常规全波形反演的另一种实施例的结果；

图10是说明利用相位外推的本发明另一种实施例的流程图；

图11说明了利用相位外推的实施例的结果；以及

图12示意性地说明了用于执行根据本发明实施例的方法的系统。

具体实施方式

本发明可以在由计算机执行的系统与计算机方法的通用背景下描述和实现。这种计算机可执行的指令可以包括可以用于执行特定任务并处理抽象数据类型的程序、例程、对象、组件、数据结构和计算机软件技术。为了在各种计算平台和环境中应用，本发明的软件实现可以用不同语言编码。应当认识到，本发明的范围和根本原理不限于任何特定的计算机软件技术。

而且，本领域技术人员将认识到，本发明可以利用硬件和软件配置中的任何一种或者其组合来实践，包括但不限于具有单个和/或多个计算机处理器的系统、手持式设备、可编程消费者电子产品、迷你计算机、大型计算机等。本发明还可以在分布式计算环境中实践，其中任务是由通过一个或多个数据通信网络链接的服务器或其它处理设备执行的。在分布式计算环境中，程序模块可以既位于本地又位于远端计算机存储介质中，存储介质都包括存储器存储设备。

而且，供计算机处理器使用的制造品，例如CD、预先录制的盘或其它等效设备，可以包括计算机程序存储介质和记录在其上的程序装置，所述程序装置指示计算机处理器方便本发明的实现与实践。这种设备与制造品也在本发明的主旨与范围内。

现在将参考附图描述本发明的实施例。本发明可以许多种方式实现，包括例如作为系统（包括计算机处理系统）、方法（包括计算机实现的方法）、装置、计算机可读介质、计算机程序产品、图形用户界面、门户网站或者有形地固定在计算机可读存储器中的数据结构。以下讨论本发明的几种实施例。附图仅仅说明了本发明的典型实施例而且不应当因此就认为是对其范围和广度的限定。

本发明涉及计算地球地下的物理属性而且，作为例子但不是限制，可以利用唯相全波形反演计算速度模型。

为了开始本发明的解释，首先考虑在图1的流程图中说明的基本全波形反演方法100。在步骤10，我们获得地球属性的初始模型，作为例子但不是限制，这个地球属性是速度。全波形反演是一种局部优化方法，因此强烈依赖于优化开始的位置。对于常规的全波形反演，关于让非线性演进收敛到实际解需要什么，对于初始模型存在严格条件：初始模型必须在最低可用时间频率生成所观察到的数据的半波周期内的数据。重要的是注意到通过常规方法，没有容易的途径来确定初始模型是否满足这个条件，而且优化很容易由于差的初始模型而失败。

在步骤12，地震模拟引擎使用地球属性的初始模型来生成模拟的地震数据。总的来说，依赖于各种因素，像模拟域的尺寸/广度和可用存储器的量，模拟能够没有障碍地在时域或频域（时间傅立叶变换）中执行。大型3D勘测一般需要时域模拟，因为频域模拟对于大量的模型参数是极端存储器密集的。频域模拟的一个显著优点是可以直接访问振幅和相位两者，而且这允许使用“唯相”方法，这种方法可被调整成受运动学支配而非受振幅支配。

在步骤14，我们计算将测量记录的地震数据和模拟的地震数据之间的错配的目标函数。常规全波形反演最广泛使用的目标函数是简单最小二乘方：对于所有源、接收器和记录的时间样本，观察到的数据和模拟的数据之间的差的平方和。然而，这不意味着限制；其它目标函数也可被使用，包括相关、L1范数以及混合或长尾范数。可在时域中或者在诸如频域的变换域中构造所述目标函数。

在时域中，最小二乘方目标函数可以采取下列形式：

$E=\frac{1}{2}{\mathrm{\Σ}}_s{\mathrm{\Σ}}_r{\mathrm{\Σ}}_t{[{\mathrm{\ψ}}_{\mathrm{obs}}(t,r,s)-{\mathrm{\ψ}}_{\mathrm{mod}}(t,r,s)]}^2$       等式1

其中，E是目标函数，s是源，r是接收器，t是时间，ψ_obs是记录的数据，而ψ_mod是模拟的数据。这个目标函数遭受地震数据是频带有限的严重缺陷。频带有限信号的求差引入了“周期跳跃”的可能性，其中模拟的数据和观察到的数据的波形状相似，足以引起小的差别，但是从绝对意义上来说却失准（至少）一个波周期。与全波形反演的局部本质一起，这有可能导致非线性优化将失败并收敛到一个局部最小值而不是全局解的可能性。

改变这个问题的特性的一种途径是改变目标函数。如果我们变换到频域，则我们能够单独地（单一地）考虑一个或多个频率分量处的目标函数。在时域中，由于对前面时间的依赖性，我们不能考虑单个时间样本。在频域中，不同频率处的响应被去耦合：一个频率处的解不依赖于任何其它频率处的解。重要的是，我们还能够区别对待振幅和相位。采用等式1的时间傅立叶变换，目标函数变成：

    等式2

其中，A_obs(ω,r,s)是在接收器r处从源s在时间频率ω观察到的数据的振幅，

是所观察到的数据的相位，A_mod(ω,r,s)是模拟的数据的振幅，而

是模拟的数据的相位。

在频域中，我们能够独立于振幅部分来考虑相位部分。对于全波形反演的唯相情况，作为例子但不是限制，最小二乘方目标函数变成：

       等式3

等式1-3中的模拟的数据可在时域或频域中生成。等式1-3的目标函数测量观察到的数据和模拟的数据之间的失配，并且在每次迭代时减少。反演可以作为时域或频域中的唯相反演而进行，只要失配能够关于一个或多个频率分量的相位被直接或间接地测量。

一旦在图1的步骤14中计算出了目标函数，就在步骤16中计算搜索方向。为了更新地球属性模型并减小观察到的数据和模拟的数据之间的错配，目标函数的梯度被用于生成用于改进模型的搜索方向。然后，地球属性模型沿着相继的搜索方向被迭代地扰动，直到达到某个满意标准。

如果我们把模拟的数据看作非线性地震模拟算子对地球属性模型的作用，则搜索方向的计算变得更加清晰。利用作为地球属性的速度（v）的例子，算子非线性意味着速度的线性变化并不一定导致模拟的数据的线性变化。

利用符号N来表示把速度模型映射到地震数据中的非线性地震模拟算子，而且这个算子对当前速度模型的作用为N(v)，我们可以重写等式1：

$E=\frac{1}{2}{\mathrm{\Σ}}_s{\mathrm{\Σ}}_r{\mathrm{\Σ}}_t{[{\mathrm{\ψ}}_{\mathrm{obs}}(t,r,s)-N\left(v\right)]}^2$         等式4

因此，关于速度的导数变成：

$\frac{\∂}{\∂v}E=-{\mathrm{\Σ}}_s{\mathrm{\Σ}}_r{\mathrm{\Σ}}_t\left(\right[{\mathrm{\ψ}}_{\mathrm{obs}}(t,r,s)-N\left(v\right)\left]\frac{\∂}{\∂v}N\left(v\right)\right)$      等式5

等式5示出用于更新地球属性模型的导数非常重要地依赖于模拟算子、模拟算子关于速度的导数以及当前地震数据残余。

全波形反演的非线性问题是通过逐次线性化来解决的。对于为速度进行反演的例子，在迭代k，这是通过在速度v_k周围进行线性化并且寻找对速度的更新δv来进行的，使得更新后的模型是：v_k+1=v_k+δv。我们需要线性化，以便计算搜索方向。给定通用的线性最小二乘方系统：

E＝||y-Ax||²        等式6

梯度或者搜索方向可被写成：

$\frac{\∂}{\∂x}E=A^f[y-\mathrm{Ax}]$             等式7

其中，A^f是线性算子A的伴随矩阵（共轭转置）。对于我们全波形反演的非线性问题，我们具有非线性算子N，而且我们需要线性化算子的伴随矩阵以便计算梯度。我们使用L作为线性化算子，而L^f是该线性化算子的伴随矩阵。算子L把速度扰动的矢量映射到波场扰动的矢量中，而且伴随矩阵算子L^f把波场扰动的矢量映射到速度扰动的矢量中（等式8）。

Lδv₁＝δψ₁

L^fδψ₂=δv₂           等式8

一旦计算出搜索方向，我们就需要确定在那个方向上采取多大的一步，这是在图1的步骤18中如何更新地球属性模型。存在至少两个可供选择的方法：非线性线搜索，或者利用，作为例子而不是限制，高斯-牛顿方法来求解线性问题。

大部分已经发表的常规方法都采用最陡下降或预处理的最陡下降用于非线性优化。一旦估计出搜索方向，这些方法就忽略当前线性问题并且使用非线性线搜索来估计在搜索方向上采取的最佳“步长”。如果我们使用δv作为搜索方向（通常是目标函数相对于速度参数的梯度），而α作为步长，则我们可把非线性线搜索表示为：

${\min }_{\mathrm{\α}}\left\{\frac{1}{2}{\mathrm{\Σ}}_s{\mathrm{\Σ}}_r{\mathrm{\Σ}}_t{[{\mathrm{\ψ}}_{\mathrm{obs}}(t,r,s)-N(v+\mathrm{\α\δv})]}^2\right\}$         等式9

非线性线搜索的一个严重缺陷是采取这么大的步长使模拟的数据关于观察到的数据变得周期跳跃。这会导致更小的残余并且导致收敛到局部最小值而非真正的全局解。

相对于使用非线性线搜索的一种可选方法是在非线性演进的每次相继线性化求解线性问题。求解线性问题避免了对线搜索的需要，因为步长选择隐含在线性优化的方法中，就像例如在共轭梯度方法中。求解线性问题需要线性化的准确方法：通过伴随矩阵测试的正演和伴随矩阵线性化算子。这常常需要大量的工作，但是会导致收敛的显著改进。利用上述线性化算子L和L^f，作为例子而不是限制，我们可以使用关于正则方程的共轭梯度来求解线性系统。我们希望求解的线性系统是：

min||Lδv-δψ||²          等式10

其中，δψ是当前残余δψ=ψ_obs-N(v_k)。

在更新地球属性模型之后，处理循环回到步骤12，在那里更新后的模型被用来生成模拟的地震数据。执行步骤14，而且，如果模拟的地震数据和记录的地震数据之间的差大，则还执行步骤16和18并且循环回到步骤12，直到步骤14的差值充分小或者循环或迭代次数达到预定的次数。

当尝试传统的全波形反演时，图1的方法100具有严重的限制。首先，全波形反演是一种局部优化方法，这意味着它对于非线性演进开始的位置是敏感的。如果初始模型远离真正的模型，则局部方法失败。这个问题影响着所有的局部方法，包括牛顿和准牛顿方法。对于常规的全波形反演，获得好的起始模型绝对是关键的。总的来说，不存在明显的途径来量化地确定一个给定的起始模型是否将收敛到真正的全局最小值。

常规全波形反演的另一个严重限制是带宽限制。在被用于生成梯度（搜索方向）的数据的时间带宽和由等式5的评估获得的梯度的空间带宽之间存在直接的关系。数据中的低时间频率产生梯度中的长空间波长。考虑图2，这个图通过在空间X和Z坐标中绘制在四个频率计算出的梯度来证明这种关系。应当指出，在0.5Hz的最低频率（面板20），计算出的梯度在空间上平滑得多。在1Hz（面板21）、1.5Hz（面板22）和2Hz（面板23），梯度变得逐渐不太平滑。地震数据的带宽是有限的，而且，如果初始模型中不存在速度的正确长空间波长，则常规的全波形就不能恢复它们，而且总的来说将会失败并且收敛到一个局部最小值而非真正的全局解。这直接暗示我们应当在最低的可用频率反演地震数据，以便采用修改速度的长空间波长的梯度。然而，地震数据的最低可用频率常常没有低到足以恢复最长空间波长并导致全局最小值–这是本发明要解决的现有技术的关键限制因素。

用于常规全波形反演的初始地球属性模型的重要性的例子可以在图3和4中看到。在图3中，初始速度模型可以在面板30中看到。这是面板38中的真正速度模型的平滑版本。面板31-37示出了常规全波形反演在8个相继的频率：1、3、5、7、9、11和13Hz的结果。当与面板38中的真正速度模型比较时，面板37中的最终结果是相当准确的。

在图4中，面板40中的初始速度模型是恒定的并被设定成水速度。这远离面板48中的真正速度模型。面板41-47示出了常规全波形反演在8个相继频率：1、3、5、7、9、11和13Hz的结果。尽管该模型最上面的部分被准确地恢复，但是较深的部分收敛到离真正解非常远的局部最小值。我们可以从图3和4中得出结论，为了收敛到正确的解，常规的全波形反演必须具有好的初始地球属性模型。

基于图1的方法100、图2的梯度带宽以及图3和4中所证明的常规全波形反演的初始模型需求，本发明人已经确定需要一种用于全波形反演的新方法。本发明克服了常规方法的带宽和初始模型限制。

本发明的一种实施例在图5中通过方法500来描述。方法500的许多步骤都与图1中的方法100的步骤相似，但是方法500不遭受传统全波形反演的限制。开始，在步骤50，本发明设定任意的初始地球模型，例如，作为例子但不是限制，把整个初始模型设定成1500m/s的水速度。在步骤51，这个初始模型被用于生成模拟的地震数据。模拟的地震数据的正演模拟可以通过许多已知的正演模拟算法中的任意一种，例如有限差分模拟，在时域中或者在频域中进行。如果正演模拟在时域中进行，那么然后它可被变换到频域。在步骤52，获得记录的地震数据，并且在步骤53中把它变换到频域中。当模拟的地震数据和记录的地震数据都在频域中时，可以在步骤54中计算残余相位，这是模拟的地震数据和记录的地震数据的相位部分之间的差。在步骤55，所述残余相位被相位展开。还可能分别展开模拟的地震数据和记录的地震数据的相位。然后，展开的相位可被用于计算展开的残余相位。

相位展开确保2π的全部适当倍数都已被包括在数据的相位部分中，这意味着相位是连续的，而没有跳过2π。存在多种相位展开的方法，但是许多对于甚至适中的频率，例如大于2Hz的那些频率，都失败了。为此，本发明人开发出一种用于相位展开的新方法来准备用于反演的频域数据。该新方法使用特定类型的左预处理，这种处理减轻大相位跃变的影响。观察到的相位和模拟的相位可被独立地展开，或者它们的差值：残余相位可被展开。后者是优选的，因为相邻数据点之间的相位差将更小。

我们用于相位展开的过程是通过矢量微积分的基本定理，也被称为Helmholtz分解，来激发的。Helmholtz分解可被用于把矢量场分解成无旋度分量和无散度分量。我们只对无旋度分量感兴趣，因此我们不需要精确的Helmholtz分解。无旋度分量是标量势的梯度，而且是守恒场。守恒场是任意点之间的线积分独立于路径的矢量场。我们识别出具有标量势的展开残余相位，其梯度是Helmholtz分解的守恒场。

我们通过取得输入包裹相位的梯度并且通过添加或减去2π而进行调整使得结果位于范围[-π,+π]中来开始。这个“调整后的相位”也称为相位的“主要值”。在这里，“梯度”意味着分别沿源和接收器的方向的导数值。我们可以如下写调整后的相位梯度到守恒场上的投影：

           等式11

其中，

是展开的残余相位，而g是调整后的包裹相位的梯度，如以上所解释的。

为了计算展开相位，我们使梯度算子相对于源和接收器坐标离散化并且通过最小二乘方求解等式12中所示的超定系统。在一种实施例中，我们发现稀疏QR因式分解是用于求解这个等式系统的一个特别有效的方法。

            等式12

这种对于相位展开投影到守恒场上的方法在比1Hz大得多的中等频率有困难。对于n_s个源和n_r个接收器，等式12的系统相对于源坐标将具有n_s＊n_r行调整后的梯度，而相对于接收器坐标将具有n_s＊n_r行调整后的梯度。因此，两倍超定。

我们发现系统的失败与调整后的梯度的条目的大量值有关，而且，通过向下给这些大量值条目加权，我们能够显著地改进健壮性，其中向下加权具有减小它们在等式系统中的重要性的效果。在一种实施例中，其条目与调整后的梯度的量值成反比的对角线左预处理器的应用大大提高了较高频率处的相位展开的性能。其它类型的预处理器也可被使用并且在本发明的范围之内。

新系统在等式13中示出，其中，左预处理器W的第k个元素与升高到幂α的调整后的梯度的第k个元素的分量的量值成反比。

w_k，s=|g_k，s|^-α

W_k，r＝|g_k，r|^-α            等式13

在一种实施例中，这种用户定义的正幂α可以设定成2.5。利用这种实施例，有和没有预处理器的相位展开的例子可以在图6中处于.5Hz的数据并且在图7中处于1.5Hz的数据看到。图6和图7都在面板A中示出包裹相位，在面板B中示出没有使用预处理器的展开相位，并且在面板C中示出利用左对角线预处理器的展开相位。在图6的低频情况下，有和没有预处理器的展开结果差别很小。然而，在图7中，没有预处理器的结果错误地改变了由D和E指示的区域中的相位，这表示随着频率变高，预处理对于获得好的结果是必要的。

我们指出，这种相位展开方法不需要边界条件的积分或规范，以便从包裹相位的梯度的主要值中获得展开相位。

在另一种实施例中，相位展开可被用在非线性线搜索中，其中，用于速度更新的搜索方向已经预先确定。至少有两种可供选择的方法。在一种可供选择的方法中，使用常规的目标函数，但是排除了其残余相位量值超过π的数据。这暗示线搜索只对未被周期跳跃的数据敏感。在另一种可供选择的方法中，用于非线性线搜索的目标函数用展开残余相位的最小二乘方之和代替。这意味着线搜索将正确地处理周期跳跃的数据。这导致与等式3中所示非常相似的目标函数，但是具有如等式14中所示的展开残余相位我们进一步指出，展开残余相位可被用作用于随机或贝叶斯反演的目标函数，以便准确地处理周期跳跃的数据。

        等式14

尽管已经关于准备用于反演的地震数据解释了利用预处理器的本相位展开方法，但是这不意味着限制。本领域技术人员将认识到，展开的地震数据在其它处理流中可能是有用的，例如层位拉平、同态反褶积、折射静校正和残余对准；而且，其它类型的数据，例如合成孔径雷达，可以受益于利用预处理器的这种相位展开方法。

再次参考图5，一旦展开的残余相位可用，步骤55就计算测量记录的数据和模拟的地震数据的相位之间的错配的目标函数。在一种实施例中，这个目标函数可以是等式3。在这种情况下，我们执行唯相全波形反演。为此，我们在步骤56中计算搜索方向、在步骤58中更新地球属性模型，并且对步骤51、54、55、56、57和58进行迭代，直到目标函数足够小或者达到预定次数的迭代。

在一种实施例中，当我们通过唯相全波形反演进行迭代时，我们可以通过使用规则化逐次迭代并且把它们约束到低波数更新的延拓方法，来提高我们恢复长空间波长（例如用于速度的那些空间波长）的能力。该延拓方法是对平滑规则化应用同伦用于非线性优化。在这里，同伦意味着对于平滑规则化以大的量值开始，并且随着非线性演进过程逐步减小平滑规则化的量值。

平滑规则化能够通过向线性系统添加行来实现，以惩罚被优化的模型中的粗糙度。有许多其它途径来实现粗糙度惩罚。在一种实施例中，延拓方法可以使用代表慢度的多项式的解析导数。对平滑函数的基变化，例如径向基函数，也起作用。其它可能性包括但不限于具有右预处理器的空间傅立叶基，该右预处理器通过居中或者不居中的波数以及1次或2次导数值进行缩放。在还有另一种实施例中，粗糙度惩罚可以通过对像素化的模型应用1次正演数值差来应用。这些例子不意味着限制；本领域技术人员将认识到，存在多得多的落在本发明范围内的可以在延拓方法的背景下使用的可能规则化算子。

通过利用一阶数值差使用导数惩罚来扩展平滑规则化的思想，让我们以简单的3x3像素化速度模型开始。在两维空间中，9个速度（v_x,z）看起来像：

v1,1	v2,1	v3,1
			v1,2	v2,2	v3,2
v1,3	v2,3	v3,3

表1：3x3速度模型

把这个速度模型写成列矢量，我们得到：

$\left(\begin{array}{c}{v}_{\mathrm{1,1}}\\ v_{\mathrm{1,2}}\\ v_{\mathrm{1,3}}\\ v_{\mathrm{2,1}}\\ v_{\mathrm{2,2}}\\ v_{\mathrm{2,3}}\\ v_{\mathrm{3,1}}\\ v_{\mathrm{3,2}}\\ v_{\mathrm{3,3}}\end{array}\right)$

我们可以通过惩罚相邻速度的差值，例如（v_1,1-v_1,2），来应用水平导数惩罚（X方向上的粗糙度惩罚）。应当指出，正式的正演数值导数被写成

但是我们可以清除分母。这导致所示水平导数惩罚矩阵：

$\left(\begin{array}{ccccccccc}+1& 0& 0& -1& 0& 0& 0& 0& 0\\ 0& +1& 0& 0& -1& 0& 0& 0& 0\\ 0& 0& +1& 0& 0& -1& 0& 0& 0\\ 0& 0& 0& +1& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& +1& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& +1& 0& 0& -1\end{array}\right)\left(\begin{array}{c}{v}_{\mathrm{1,1}}\\ v_{\mathrm{1,2}}\\ v_{\mathrm{1,3}}\\ v_{\mathrm{2,1}}\\ v_{\mathrm{2,2}}\\ v_{\mathrm{2,3}}\\ v_{\mathrm{3,1}}\\ v_{\mathrm{3,2}}\\ v_{\mathrm{3,3}}\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right)$

和一个类似构造的垂直导数惩罚矩阵：

$\left(\begin{array}{ccccccccc}+1& -1& 0& 0& 0& 0& 0& 0& 0\\ 0& +1& -1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& +1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& +1& -1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& +1& -1& 0\\ 0& 0& 0& 0& 0& 0& 0& +1& -1\end{array}\right)\left(\begin{array}{c}{v}_{\mathrm{1,1}}\\ v_{\mathrm{1,2}}\\ v_{\mathrm{1,3}}\\ v_{\mathrm{2,1}}\\ v_{\mathrm{2,2}}\\ v_{\mathrm{2,3}}\\ v_{\mathrm{3,1}}\\ v_{\mathrm{3,2}}\\ v_{\mathrm{3,3}}\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right)$

注意，因为导数只涉及水平或垂直相邻的像素，所以行比列少。

这些水平和垂直导数矩阵还可被写成：

λ_xD_xv＝0

λ_zD_zv＝0            等式15

其中v是速度的列矢量，D_x是水平导数的矩阵，D_z是垂直导数的矩阵，而λ_x和λ_z是拉格朗日乘子。

延拓方法以大的拉格朗日乘子λ_x和λ_z开始，而且因此第一“延拓步骤”中的初始解非常平滑。很清楚，这有助于恢复速度的长空间波长。随着非线性演进前进，我们采取附加的延拓步骤，而且λ_x和λ_z的量值被减小。随着惩罚的量值被减小，在速度模型中允许相继缩短的空间波长。

对于设定初始的λ_x和λ_z值，有许多种可能的选项。如果选得充分大，那么模型中只允许非常长的空间波长，而且非线性演进有效地变得独立于初始模型。如果选得太小，问题将不能被足够规则化，而且丧失了与初始模型的独立性。用于这些参数的初始值的一种实施例是通过在每个相继线性化的线性化算子的算子范数来规格化它们。在第一次线性化中非线性问题的开始，如果我们具有线性系统Ax=y，则我们就把λ_x和λ_z设定成被算子范数||A||缩放。||A||可以使用例如乘幂方法获得。

在本发明中执行的唯相全波形反演还可以包括在每次迭代更准确地求解线性问题。如果我们在每次相继的线性化求解高斯-牛顿问题以获得模型更新，而不是采用最陡下降和线搜索的组合，那么我们得到改进的结果。

对于全波形反演的非线性问题，我们围绕在迭代k(v_k)的速度进行线性化，并且设法获得对速度的更新δv，使得更新后的模型是：v_(k+1)=v_(k)+δv。这是相继线性化。对线性问题应用导数惩罚暗示着我们希望对模型的更新是平滑的，如下所示：

λ_xD_xδv＝0

λ_zD_zδv＝0            等式16

更期望的方法是规则化非线性问题。这暗示着我们希望更新后的模型是平滑的：

λ_xD_x(v_k+δv)＝0

λ_zD_z(v_k+dv)＝0           等式17

这需要非零的右侧，但是右侧很容易通过对当前速度应用导数算子D_x和D_z而获得：

λ_xD_xδv＝-λ_xD_xv_k

λ_zD_zδv＝-λ_zD_zv_k            等式18

图8示出了本发明实施例，利用具有左预处理器的相位展开、延拓方法并求解相继线性问题的唯相全波形反演的结果。面板80是初始模型，这是恒定的1500m/s（水的速度）。这是与图4的面板40中所示相同的初始模型。图8中的面板88示出了真正的速度模型。面板81-87示出了从该初始模型开始以1Hz的相继非线性迭代。面板81示出，在一次迭代之后，准确的长空间波长出现在反演后的模型中，而且它们随着迭代前进通过面板82-87被精炼。七次非线性迭代允许恢复缺少的速度的长空间波长，这是利用常规方法不可能实现的，如图4中所看到的。

在本发明的另一种实施例中，由唯相全波形反演生成的模型可被用作常规全波形反演的初始模型。这在图9中举例说明，其中面板90中的用于常规全波形反演的初始模型是通过图8中的唯相全波形反演的7次迭代（面板87）生成的模型。在2.5Hz执行常规全波形反演的5次迭代（面板91-95）导致与面板96中的真正速度模型非常相当的反演模型（面板95）。

图10说明了本发明的再一种实施例。在这种实施例中，唯相全波形反演流作为方法1000示出。其步骤与图5中的方法500的那些相同，但是在相位展开步骤1006之后添加了步骤1007-相位外推。步骤1001、1002、1003、1004、1005、1006、1008、1009和1010分别以与步骤50、51、52、53、54、55、56、57和58相同的方式执行。步骤1007是相位外推步骤，它可被用于把展开的相位外推到比记录的地震数据中存在的更低的频率。然后，这种非常低的频率相位信息可在步骤1008、1009和1010中被用于帮助恢复构成速度模型的非常长的空间波长。

该相位外推方法使用线性相移和行进时间之间的关系：

            等式19

其中，

是频率f₁处的相位，而t是行进时间。为了把相位外推到另一个频率f₂并且假定行进时间不变，我们求解t并且替换它：

            等式20

             等式21

在这种实施例中，相位外推到比观察到并且常规可使用的那些更低的频率。常规可使用的频率一般大于2Hz。这是通过作为频率的函数的展开相位的线性化而进行的，而且可被应用于观察到的相位、模拟的相位或者残余相位。然后，利用被定义成测量相位失配的某种目标函数来反演外推的数据。当相位是频率线性时，该方法可应用于任何情况。

图11说明了相位外推方法的一种实施例的结果。面板110是初始模型，在这种情况下是1500m/s的恒定水速度，而面板121是真正的速度模型。面板111-115分别是从2.5Hz到0.1、0.2、0.3、0.4和0.5Hz的相位外推反演。面板116-120是从面板115中的相位外推结果继续的频率2.5、4.5、6.5、8.5和10.5Hz处的常规反演。

本领域技术人员将认识到，相位外推数据存在许多其它可能的使用。作为例子但不是限制，在合成孔径雷达（SAR）成像方法之前，可以获得SAR数据、利用预处理器进行相位展开并且相位外推。此外，已经利用预处理器被相位展开和相位外推的数据可被用于评估成本函数。一个例子是使用展开的相位来计算用于随机或贝叶斯优化的目标函数，其优点是该成本函数将正确地处理周期跳跃的数据。

尽管以上已经关于两维模型解释了实施例，但是所述方法很容易扩展到三维和多参数地球模型。本发明中公开的用于相位展开、相位外推和唯相全波形反演的方法可被扩展到多维并且仍然在本发明的范围之内。

用于执行方法的系统1200在图12中示意性地说明。该系统包括数据存储设备或存储器130。数据存储设备130包含记录的数据而且可以包含初始模型。记录的数据可以让处理器131（例如可编程通用计算机）获得。处理器131配置成执行初始模型模块132，以在必要的时候创建初始模型或者从数据储存器130接收初始模型。处理器131还配置成执行用于把记录的数据及可选地模拟的数据变换到频域的域变换模块133、用于基于初始模型和更新后的模型正演模拟数据的数据模拟模块134、用于利用预处理器进行相位展开并可选地相位外推记录的数据的相位准备模块135、用于计算比较模拟的数据和相位展开的记录数据的目标函数的目标函数模块136、用于确定搜索方向的搜索方向模块137、以及用于更新模型的模型更新模块138。处理器131还配置成重复地执行模块134、135、136、137和138，直到来自目标函数模块137的结果满足用户需求或者达到最大迭代次数。处理器131可以包括诸如用户接口139的接口组件，其中用户接口可以包括显示器和用户输入设备两者，并被用于实现上述根据本发明实施例的变换。用户接口可被用于显示数据和处理后的数据产品，并允许用户在用于实现所述方法的各方面的选项之间进行选择。

尽管在前面的说明书中已经关于其某些优选实施例描述了本发明，而且，为了说明，已经阐述了许多细节，但是，对本领域技术人员来说很显然，在不背离本发明基本原理的情况下，本发明很容易有改变，而且在此所述的某些其它细节可以有显著变化。此外，应当认识到，在本文任何一种实施例中示出或描述过的结构特征或方法也可以在其它实施例中使用。

Claims (16)

1.一种用于从实际地震数据确定感兴趣地下区域的属性的计算机实现的方法，包括：

a．获得代表地下区域的实际地震数据和用于地下区域的初始地球属性模型；

b．使用所述初始地球属性模型来执行正演模拟，以创建与实际地震数据具有相似采集规范的模拟的地震数据；

c．把模拟的地震数据和实际地震数据变换到时间傅立叶频域，以创建频域模拟的地震数据和实际地震数据，其中，频域模拟的地震数据和实际地震数据包括振幅部分和相位部分；

d．测量频域模拟的地震数据和频域实际地震数据之间的错配，以产生频域残余地震数据；

e．执行所述频域残余地震数据中的某些观察到的频率分量的相位部分的相位展开，以创建展开的残余相位部分；以及

f．反演展开的残余相位部分，以确定感兴趣地下区域的期望属性，其中，所述反演最小化被定义成测量所述错配的目标函数，而且其中，执行正演模拟、变换、测量、执行相位展开和反演步骤是由计算机处理器执行的。

2.如权利要求1所述的方法，进一步包括在所述反演步骤之前外推展开的相位。

3.如权利要求1所述的方法，进一步包括第二反演步骤，所述第二反演步骤反演所述频域残余地震数据中的某些观察到的频率分量的振幅和相位部分，其中，所述第二反演步骤使用通过反演展开的残余相位部分而确定的期望属性作为新的初始模型。

4.如权利要求1所述的方法，其中，所述反演使用延拓方法来规则化相继迭代。

5.如权利要求1所述的方法，其中，所述反演为每次相继线性化求解线性反演问题。

6.如权利要求5所述的方法，其中，所述线性反演问题是利用高斯-牛顿方法来求解的。

7.一种用于从地震数据确定感兴趣地下区域的属性的系统，包括：

a．数据源，包含计算机可读数据，所述数据包括代表感兴趣地下区域的实际地震数据；

b．处理器，配置成从计算机模块执行计算机可读代码，所述计算机模块包括：

i．初始模型模块，配置成从数据源接收初始属性模型或者生成所述初始属性模型；

ii．域变换模块，把实际地震数据变换到频域中，以生成频域地震数据；

iii．数据模拟模块，从初始属性模型或者更新后的属性模型中生成模拟的地震数据；

iv．残余计算模块，计算残余波场；

v．相位准备模块，相位展开所述频域地震数据；

vi．反演模块，计算目标函数、搜索方向并更新初始属性模型或者之前更新的属性模型；

c．用户接口。

8.如权利要求7所述的系统，其中，所述数据源还包含初始属性模型。

9.如权利要求7所述的系统，其中，所述处理器重复执行所述数据模拟模块、残余计算模块、相位准备模块和反演模块，以最小化所述目标函数。

10.如权利要求7所述的系统，其中，所述计算机模块进一步包括相位外推模块。

11.如权利要求7所述的系统，其中，所述反演模块执行唯相全波形反演，直到唯相目标函数已被最小化或者达到用户规定的迭代次数，以生成更新后的唯相反演属性模型。

12.如权利要求11所述的系统，其中，在生成更新后的唯相反演属性模型之后，所述反演模型利用所述更新后的唯相反演属性模型作为初始属性模型来执行常规的全波形反演。

13.如权利要求7所述的系统，其中，所述反演模块使用延拓方法来规则化相继迭代。

14.如权利要求7所述的系统，其中，所述反演模块为每次相继线性化求解线性反演问题。

15.如权利要求14所述的系统，其中，所述线性反演问题是利用高斯-牛顿方法来求解的。

16.一种包括其中包含计算机可读代码的计算机可读介质的制品，所述计算机可读程序代码适于被执行以实现从实际地震数据确定感兴趣地下区域的属性的方法，该方法包括：

a．经由计算机，使用初始地球属性模型来执行正演模拟，以创建与实际地震数据具有相同或相似采集规范的模拟的地震数据；

b．经由计算机，把模拟的地震数据和记录的地震数据变换到时间傅立叶频域，以创建频域模拟的地震数据和频域记录的地震数据，其中，频域模拟的地震数据和频域记录的地震数据包括振幅部分和相位部分；

c．经由计算机，测量频域模拟的地震数据和频域记录的地震数据之间的错配，以产生频域残余地震数据；

d．经由计算机，执行所述频域残余地震数据中的某些观察到的频率分量的相位部分的相位展开，以创建展开的残余相位部分；以及

e．经由计算机，反演展开的残余相位部分，以确定感兴趣地下区域的期望属性，其中，所述反演最小化被定义成测量所述错配的目标函数。

Images