CN103221842A - 用于利用相位外推的数据反演的系统和方法 - Google Patents

Abstract

公开了用于反演来自受关注区域数据以确定受关注区域的物理性质的系统和计算机实施的方法。所述方法包括把所述数据变换到傅里叶频率域以获得频率域数据，其中所述频率域数据包括振幅部分和相位部分，进行所述频率域数据的相位部分的相位展开以产生展开数据的展开的相位部分，外推所述展开的相位部分以产生外推的展开数据，以及反演所述外推的展开数据以确定所述受关注区域的物理性质。所述反演的数据可以是例如地震数据或合成孔径雷达数据。所述系统包括数据源、用户接口以及的被配置为执行设计为实施所述方法的计算机模块的处理器。

Description

用于利用相位外推的数据反演的系统和方法

技术领域

一般来说，本发明涉及用于反演地震数据以计算地下的物理性质的方法和系统，尤其是用于进行纯相位全波形反演以从地震数据计算速度模型的方法和系统。

背景技术

地下勘探，尤其是对油气藏的勘探，典型情况下使用若干方法，比如偏移地震数据，以产生地下的可解释图像。在断层、盐体等造成的复杂地下区域中，传统的偏移方法往往不能产生适当的图像。另外，传统的偏移方法要求相当准确的地下速度模型；这样的速度模型也可以从地震数据确定但是在专门知识和计算成本上可能都非常昂贵。

从地震数据计算速度模型有许多常规方法，包括NMO速度分析、偏移速度分析、层析和全波形反演。某些方法，比如全波形反演，在计算上非常昂贵并且因为计算能力已经上升仅仅在最近才变成了实用的方法。常规的全波形反演在时间域中或在变换域中进行，比如时间傅里叶变换域或拉普拉斯变换域。这些方法的失败往往由于地震数据中缺失低频，典型情况下低于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)是在时间频率ω，来自震源s的在接收器r观测到的数据的振幅，

(ω,r,s)是观测数据的相位，A_mod(ω,r,s)是模拟数据的振幅，而(ω,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

梯度或搜索方向能够写为：

    等式7

其中

是线性算子A的伴随矩阵（共轭转置）。对于我们的全波形反演的非线性问题，我们有非线性算子N，并且我们需要线性化算子的伴随矩阵以便计算梯度。我们使用L做线性化算子，而

做线性化算子的伴随矩阵。算子L将速度扰动的矢量映射为波场扰动矢量，而该伴随矩阵算子将波场扰动矢量映射为速度扰动矢量（等式8）。

Lδ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和

我们能够使用例如而非限制，标准方程上的共轭梯度求解线性系统。我们要求解的线性系统是：

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的常规全波形反演的结果。画面37中的最终结果与画面38中的真实速度模型对比时相当准确。

在图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开始，本发明设置了任意初始地下模型，例如而非限制，设置整个初始模型为1500米/秒的水速。在步骤51这个初始模型用于产生模拟地震数据。通过许多公知的正向模拟算法比如有限差分模拟都可以在时间域或在频率域进行模拟地震数据的正向模拟。如果在时间域进行正向模拟，它可以然后变换到频率域。在步骤52，获得了记录的地震数据而在步骤53它被变换到频率域。当模拟地震数据和记录的地震数据都在频率域时，在步骤54可以计算剩余相位，它是模拟与记录地震数据相位部分之间的差异。在步骤55，剩余相位被相位展开。也有可能分别地展开模拟地震数据和记录地震数据的相位。展开的相位然后可以用于计算展开的剩余相位。

相位展开确保了2π的一切适合的倍数都已经被包括在数据的相位位部分，意味着相位是连续的而不是按2π跳跃的。相位展开有若干方法但是甚至对于适中频率比如大于2Hz的频率许多都失败了。由于这个原因，本发明已经开发相位展开的新方法为反演准备频率域数据。新方法使用特定类型的左预处理减轻大相位跳跃的影响。或者可以单独地展开所观测的相位或模拟的相位，或者可以展开其差异──剩余相位。优选的是后者，因为邻近数据点之间的相位差将更小。

对于相位展开我们使用的过程受到矢量演算的基本定理启示，也称为赫尔姆霍茨分解。赫尔姆霍茨分解能够用于把矢量场分解为无旋分量和无散分量。我们仅仅关注无旋分量，因此我们不要求精确的赫尔姆霍茨分解。无旋分量是标量势的梯度，并且是守恒场。守恒场是这样的矢量场，任意点之间的线积分都是路径无关的。我们使用其梯度是赫尔姆霍茨分解的守恒场的标量势识别展开的剩余相位。

我们从获得输入的重叠相位的梯度开始，并且通过增加或减去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中0.5Hz和图7中1.5Hz的数据，能够看到有没有预处理器的相位展开的实例。图6和图7都显示了画面A中的重叠相位、画面B中没有使用预处理器的展开相位以及画面C中使用左对角线预处理器的展开相位。在图6的低频情况下，有没有预处理器的展开结果差异不大。不过在图7中，没有预处理器的结果已经错误地改变了由D和E所指示的区域中的相位，表明随着频率逐渐更高，为了得到好的结果预处理是必要的。

我们注意到，为了从重叠相位梯度的主值获得展开的相位，这种相位展开方式不要求边界条件的完整或规范。

在另一个实施例中，相位展开可以用在速度更新的搜索方向已经预先确定的非线性线性搜索中。至少有两个备选方案。在一个备选方案中，使用了常规的目标函数，但是排除了剩余相位幅度超过π的数据。这暗示线性搜索仅仅对没有周期跳跃的数据敏感。在另一个备选方案中，以展开的剩余相位的最小平方和替换了非线性线性搜索的目标函数。这意味着线性搜索将正确地处理周期跳跃的数据。这引起目标函数非常类似于等式3所示，但是具有等式14所示的展开的剩余相位

我们进一步注意到，展开的剩余相位有可能用作随机或贝叶斯反演的目标函数以便正确地处理周期跳跃的数据。

    等式14

尽管按照准备用于反演的地震数据已经讲解了具有预处理器的相位展开的本方法，但是这不意味着限制。本领域技术人员将认识到，展开的地震数据也可以用在其他处理流程中，比如层拉平、同态反褶积、折射波静校正和剩余对齐；以及其他类型的数据，比如合成孔径雷达，也可以从这种具有预处理器的相位展开方法中受益。

再次参考图5，一旦展开的剩余相位可用，步骤55就计算目标函数，度量在记录数据与模拟地震数据的相位之间的不吻合。在某实施例中，这种目标函数有可能是等式3。在这种情况下，我们执行纯相位全波形反演。为了这样做，我们在步骤57计算搜索方向、在步骤58更新地下性质模型，并且迭代步骤51、54、55、56、57和58直到目标函数足够小或者已经到达了迭代的预定次数。

在某实施例中，当我们在纯相位全波形反演中迭代时，我们能够改进我们恢复空间长波长的能力，比如对于速度的长波长，方式为使用延拓方式调整依次迭代并使其约束到低波数更新。延拓方式为对非线性最优化的平滑正则化的同伦应用。这里同伦意味着平滑正则化以大的幅度开始，并且在非线性演变进程上逐步地减小平滑正则化的幅度。

为了补偿被优化模型中的粗糙，通过向线性系统增加若干行能够实现平滑正则化。有无数种其他方式实现粗糙度补偿。在一个实施例中，延拓方式可以使用表示慢度的多项式的解析导数。改变平滑函数的基，例如径向基函数，也起作用。其他可能性包括但不限于具有以波数缩放的右预处理器的空间傅里叶基，以及一阶或二阶数值导数，或者居中或者不居中。在又一个实施例中，通过向像素化模型应用第一向前数值差分可以施加粗糙度补偿。这些实例不意味着限制；本领域技术人员将认识到，存在着许多更可能的在延拓方式的语境中可以使用的正则化算子，它们落入本发明的范围。

使用一阶数值差分，通过使用导数处罚，对平滑正则化的思路进行扩展，让我们从简单的3×3像素化速度模型开始。在二维空间中，9个速度(v_x,z)将显现为：

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

表1：3×3速度模型

[0001]将这个速度模型写为列矢量，我们得到：

$\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+δv)=0    等式17

这要求非零的右边，但是通过把导数算子D_x和D_z应用到当前速度，很容易获得该右边：

λ_xD_xδv=-λ_xD_xv_k

λ_zD_zδv=-λ_z    等式18

图8展示了本发明实施例的结果，纯相位全波形反演使用了具有左预处理器的相位展开、延拓方式以及求解依次的线性问题。画面80是初始模型，是常数1500米/秒（水速）。这与图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是初始模型，在这种情况下是1500米/秒的不变水速而画面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成像方法。此外，使用预处理器已经被展开的和相位外推的数据然后可以用于估计成本函数。一个实例是使用展开的相位计算随机或贝叶斯最优化的目标函数，优点是成本函数将正确地处理周期跳跃的数据。

尽管已经按照二维模型讲解了以上的实施例，但是本方法不难扩展到三维和多参数的地下模型。本发明公开的相位展开、相位外推和纯相位全波形反演的方法可以扩展到多维并保持在本发明的范围之内。

图12示意地展示了执行本方法的系统1200。该系统包括数据存储设备或存储器130。数据存储设备130包含记录的数据并可以包含初始模型。可以使记录的数据对处理器131可用，比如可编程通用计算机。处理器131被配置为执行初始模型模块132以便在必要时创建初始模型或者从数据存储器130接收初始模型。处理器131也被配置为执行域变换模块133把记录的和可选地模拟的数据变换到频率域、数据模拟模块134根据初始和更新的模型正演模拟数据、相位准备模块135对记录的数据以预处理器相位展开和可选地相位外推、目标函数模块136计算目标函数将模拟数据与相位展开的记录数据对比、搜索方向模块137确定搜索方向以及模型更新模块138更新模型。处理器131也被配置为重复地执行模块134、135、136、137和138到来自目标函数模块136的结果满足了用户需求或者达到了迭代的最大次数。处理器131可以包括接口组件比如用户接口139，显示器和用户输入设备它都可以包括，并且被用于根据本发明的实施例执行以上介绍的变换。用户接口可以被用于既显示数据和处理后数据产品，也允许用户在若干选项当中选择以实施本方法的若干方面。

虽然在以上说明书中已经关于本发明一定的优选实施例介绍了它，以及为了展示目的已经阐述许多细节，但是对于本领域技术人员来说显而易见，本发明能容许变化并且本文介绍的某些其他细节能够有相当大改变而不脱离本发明的基本原理。此外，应当认识到，在本文任何一个实施例中显示和介绍的结构特征或方法步骤也能够在其他实施例中使用。

Claims (13)

1.一种用于反演来自受关注区域的数据以确定所述受关注区域的物理性质的计算机实施的方法，包括：

a．把所述数据变换到傅里叶频率域以获得频率域数据，其中所述频率域数据包括振幅部分和相位部分；

b．进行所述频率域数据的相位部分的相位展开以产生展开数据的展开的相位部分；

c．外推所述展开的相位部分以产生外推的展开数据；以及

d．反演所述外推的展开数据以确定所述受关注区域的物理性质，其中变换、进行相位展开、外推以及反演步骤由计算机处理器执行。

2.根据权利要求1所述的方法，其中，所述展开的相位部分被外推以产生在比所述数据中常规可用频率更低频率的相位值。

3.根据权利要求2所述的方法，其中，所述常规可用频率高于2Hz。

4.根据权利要求1所述的方法，其中，所述外推使用傅里叶频率域中的线性关系。

5.根据权利要求1所述的方法，其中，所述反演包括全波形反演。

6.根据权利要求1所述的方法，其中，所述数据包括地震数据。

7.根据权利要求1所述的方法，其中，所述数据包括合成孔径雷达数据。

8.一种用于反演来自受关注区域的数据以确定所述受关注区域的物理性质的系统，包括：

a．数据源，包含计算机可读数据；

b．处理器，被配置为执行来自计算机模块的计算机可读代码，所述计算机模块包括：

i．用于把数据变换到傅里叶频率域的域变换模块；

ii．用于相位展开的相位准备模块；

iii．相位外推模块；以及

iv．反演模块；以及

c．用户接口。

9.根据权利要求8所述的系统，其中，所述相位外推模块使用傅里叶频率域中的线性关系。

10.根据权利要求8所述的系统，其中，所述反演模块执行全波形反演。

11.根据权利要求8所述的系统，其中，所述数据源包含地震数据。

12.根据权利要求8所述的系统，其中，所述数据源包含合成孔径雷达数据。

13.一种包括计算机可读介质的制成品，在该计算机可读介质中包含计算机可读程序代码，所述计算机可读程序代码适于执行用于反演来自受关注区域的数据以确定所述受关注区域的物理性质的方法，所述方法包括：

a．把所述数据变换到傅里叶频率域以获得频率域数据，其中所述频率域数据包括振幅部分和相位部分；

b．进行所述频率域数据的相位部分的相位展开以产生展开的相位部分；

c．外推所述展开的相位部分；以及

d．反演所述展开的相位部分以确定所述受关注区域的物理性质。

Images