CN103329009A - 用于利用相位展开进行数据反演的***和方法 - 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)是在时间频率ω处在接收机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

梯度或搜索方向可以被写为：

公式7

其中

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

用于线性化的算子的伴随。算子L将速度扰动的向量映射为波场扰动的向量，并且伴随算子

将波场扰动的向量映射为速度扰动的向量（公式8）。

Lδv₁=δψ₁

公式8

一旦计算了搜索方向，我们需要确定在该方向采取多大的步长，其是如何在图1的步骤18中更新地球特性模型。存在至少两个替换方式：非线性的线搜索，或使用高斯-牛顿方法（举例而非限制）求解线性问题。

公布的大部分常规方法采用最陡下降或预处理的最陡下降来进行非线性优化。一旦估计了搜索方向，这些方法不考虑当前线性问题并且使用非线性的线搜索来估计在搜索方向采取的最佳的“步长”。如果我们使用δv用于搜索方向（通常目标函数相对于速度参数的梯度）并且使用α用于步长，则我们可以将非线性的线搜索表达为：

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

非线性的线搜索的一个严重的缺点是采取如此大的步长以致建模的数据相对于观察的数据变为周期跳跃的。这可能导致更小的残差并且导致收敛到局部极小值而不是真正的全局解。

使用非线性的线搜索的可替换方式是在非线性演进的每个连续的线性化中求解线性问题。求解线性问题避免了对线搜索的需要，因为在线性最优化的机器中，例如在共轭梯度法中，步长选择是隐式的。求解线性问题需要准确的线性化机器：传递伴随测试的正演（forward）并且伴随的线性化的算子。这通常需要显著的工作，但是可能导致收敛的显著的改进。使用如上所述的线性化的算子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处的传统的全波形反演的结果。在与画板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的中频是失败的。鉴于此，发明人已经发展了用于相位展开以便为反演准备频域数据的新方法。新方法使用特定类型的左预处理，其将大的相位跳变的影响去加权。观察的相位和建模的相位可以被单独地展开，或者他们的差（残差相位）可以被展开。后者是优选的，因为相邻的数据点之间的相位差将更小。

我们用于相位展开的进程是由矢量计算的基本定理（也称为亥姆霍兹分解）启发的。亥姆霍兹分解可以用于将矢量场分解为无旋度（curl-free）的分量和无散度（divergence-free）的分量。我们仅仅对无旋度的分量感兴趣，因此我们不需要精确的亥姆霍兹分解。无旋度的分量是标势的梯度，并且是保守场。保守场是任意点之间的线积分是路径独立的矢量场。我们识别具有标势的展开的残差相位，标势的梯度是亥姆霍兹分解的保守场。

我们开始于采取输入的包裹的相位的梯度，并且通过加2π或减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|^-x

W_k,r=|g_k,r|^-x  公式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，直到目标函数足够小或已经达到预定的迭代次数。

在实施例中，因为我们迭代仅相位全波形反演，所以我们可以通过使用连续逼近以调整连续的迭代并且将它们约束到低波数更新来改善我们的恢复长的空间波长的能力，诸如用于速度的那些。连续逼近是将同伦（homotopy）应用到平滑的正则化以用于非线性优化。这里同伦意思是以用于平滑的正则化的大量值开始，并且在非线性演进的过程中逐渐减小平滑的正则化的量值。

平滑的正则化可以通过向线性***增加行以惩罚被优化的模型的粗糙度来实现的。存在许多其它实现粗糙度惩罚的方式。在一个实施例中，连续逼近可以使用表示缓慢度的多项式的解析导数。对平滑函数的基的变化，例如辐射基函数，也起作用。其它可能性包括但是不局限于具有随波数缩放的右预条件算子的空间傅里叶基、以及第一或第二数值导数（或中心或非中心）。在另一实施例中，粗糙度惩罚可以通过向像素化的模型应用第一正演数值差来应用。这些示例不意指限制；本领域技术人员将理解，存在很多可能的可以用在连续逼近的背景下的正则化算子，它们属于本发明的范围。

通过利用使用第一阶数值差的导数惩罚对平滑的正则化的构思的扩展，让我们从简单的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方向的粗糙度代价）。注意，形式的正演数值导数被写为

但是我们可以清除分母。这导致下面显示的水平导数代价的矩阵：

以及类似地构造的垂直导数代价的矩阵：

注意，存在比列更少的行，因为导数仅仅涉及横向或纵向相邻像素。

这些水平和垂直导数矩阵也可以被写为：

λ_xD_xv＝O

λ_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=O

λ_zD_zδv＝O  公式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＝-λ_zD_zv_k  公式18

图8显示本发明的实施例的结果，使用具有左预条件算子的相位展开、连续逼近和求解顺序的线性问题的仅相位全波形反演。画板80是初始模型，其是常数1500m/s（水的速度）。这是与图4所示的画板40相同的初始模型。图8中的画板88显示真正的速度模型。画板81-87显示开始于初始模型的以1Hz的连续的非线性迭代。画板81显示了，在一次迭代之后，在反演的模型中存在准确的长的空间波长，并且随着迭代进行通过画板82-87，它们被细化。七个非线性迭代允许速度的错过的长的空间波长的恢复，这是使用常规方法不可能的，如图4所示。

在本发明的另一个实施例中，由仅相位全波形反演生成的模型可以被用作传统的全波形反演的初始模型。这显示在图9中，其中画板90中的用于传统的全波形反演的初始模型是由图8中的画板87的仅相位全波形反演的7个迭代生成的模型。以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、135136、137和138，直到由目标函数模块136的结果满足用户要求或达到最大迭代数目。处理器131可以包括诸如用户接口139之类的接口组件，其可以包括显示器和用户输入设备二者，并且用于实现根据本发明的实施例的上述变换。用户接口可以被使用来显示数据和处理后的数据产物二者并且允许用户在用于实现该方法的各方面的选项当中选择。

虽然在上述说明书中，本发明已被相对于某些优选实施例进行了描述，并且许多细节已被阐述用于例示，但是本领域技术人员将显然知道，本发明可以被改变并且这里描述的某些其它细节可以显著地变化，而不背离本发明的基本原则。此外，应该理解，在这里的任何一个实施例中显示或描述的结构特征或方法步骤也可被用于其它实施例中。

Claims (15)

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

a.将数据变换到傅里叶频域以获得频域数据，其中该频域数据包括幅度部分和相位部分；

b.执行频域数据的相位部分的相位展开以生成展开的相位部分，其中该相位展开包括：

获得相位部分的梯度，

调整梯度以位于主[–π,+π]范围中以产生调整后的梯度，

将调整后的梯度设置为等于应用于展开的相位部分的梯度的离散化，以及

通过向线性方程的集合应用预条件算子来求解展开的相位部分；以及

c.反演展开的相位部分以确定关注区域的物理特性，其中变换、执行相位展开和反演步骤通过计算机处理器执行。

2.如权利要求1所述的方法，其中该预条件算子与提升到用户定义的正的幂的调整后的梯度的倒数成比例。

3.如权利要求1所述的方法，进一步包括：外插展开的相位部分。

4.如权利要求3所述的方法，其中该外插使用傅里叶频域中的线性关系。

5.如权利要求1所述的方法，其中该相位展开用在用于反演步骤的非线性的线搜索中。

6.如权利要求1所述的方法，其中该反演包括全波形反演。

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

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

9.一种用于反演来自于关注区域的数据以确定关注区域的物理特性的***，包括：

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

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

i.域变换模块，用于将数据变换到傅里叶频域；

ii.相位准备模块，用于相位展开；和

iii.反演模块；和

c.用户接口。

10.如权利要求9所述的***，其中该相位准备模块还执行相位外插。

11.如权利要求10所述的***，其中该外插使用傅里叶频域中的线性关系。

12.如权利要求9所述的***，其中该反演模块执行全波形反演。

13.如权利要求9所述的***，其中该数据源包含地震数据。

14.如权利要求9所述的***，其中该数据源包含合成孔径雷达数据。

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

a.将数据变换到傅里叶频域以获得频域数据，其中该频域数据包括幅度部分和相位部分；

b.执行频域数据的相位部分的相位展开以生成展开的相位部分，；和

c.反演展开的相位部分以确定关注区域的物理特性。

Images