CN110869814B - 用于地震数据的全波形反演的系统和方法 - Google Patents

Abstract

描述了一种使用基于树的贝叶斯方法进行全波形反演的方法，该方法自动选择模型复杂度，从而降低了计算成本。该方法可以由计算机系统执行。

Description

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

相关申请的交叉引用

本申请要求2017年7月6日提交的美国临时申请号62/529297的权益。

关于联邦赞助的研究或开发的声明

不适用。

技术领域

所公开的实施例总体上涉及使用用于从地震数据确定地下速度的技术进行地震成像，并且尤其涉及通过使用基于树的贝叶斯方法的全波形反演来确定地下速度的方法，该方法导致减少的参数以及用于描述地下速度(或其他地震性质)的基函数，从而降低了计算成本。

背景技术

地震勘探包括勘测地下地质介质中的烃沉积物。勘测通常涉及在预定位置部署地震源和地震传感器。这些源产生地震波，地震波传播到地质介质中，从而产生压力变化和振动。地质介质的物理性质的变化会引起地震波的某些性质(诸如其传播方向和其他性质)的变化。

地震波的一部分到达地震传感器。一些地震传感器对压力变化敏感(例如，水听器)，而其他地震传感器对粒子运动敏感(例如，地震检波器)，并且工业勘测可以部署一种类型的传感器或两者。响应于检测到的地震波，传感器生成相应的电信号(称为迹线)，并将它们记录在存储介质中作为地震数据。地震数据将包括多个“炮点”(被激活的地震源的各个实例)，每个“炮点”都与记录在多个传感器处的多条迹线相关联。

处理地震数据以创建地震图像，该地震图像可以被解释为识别包括烃沉积物的地下地质特征。该过程可以包括确定地下地层的速度以便执行成像。可以通过多种方法(诸如外观分析、断层成像或全波形反演)来确定速度。全波形反演(FWI)是计算量巨大的过程，其需要大量的模型参数化。一些常规的FWI方法假定最佳参数化，并且不对可变数量的参数进行尝试和采样。没有人使用基于树的概率方法。Hawkins等(2017)使用类似的思想进行机载电磁反演，Dettmer等(2016)使用类似的思想来量化海啸海平面位移的不确定性，Hawkins和Sambridge(2015)使用类似的思想来评估2D环境噪声和3D远震断层成像。但是，这些工作都是基于对地震数据无效的假设。

来自改善的地下速度的改善的地震图像允许更好地解释岩石的位置和流体性质的变化。定义岩石的位置和流体性质变化的能力对于我们做出最合适的选择来购买材料、安全操作并成功完成项目的能力至关重要。项目成本取决于对地球内部的物理边界的位置的准确预测。决策包括但不限于预算计划、获得矿产和租赁权、签署井下承诺、许可钻机位置、设计井路径和钻井策略、通过计划适当的套管和固井策略来防止地下完整性问题以及选择和购买合适的完井和生产装备。

存在对更精确的有成本效益的FWI方法的需求，以允许更好的地震成像，这反过来又允许对潜在的烃储层进行更好的地震解释，以进行油气勘探(hydrocarbonexploration)和生产。

发明内容

根据一些实施例，公开了一种使用基于树的贝叶斯方法的跨维度(transdimensional)地震全波形反演(FWI)的方法。在这种方法中，观测到的地震数据会告知模型似然性。还需要指定有关地下结构的适度信息先验作为输入。然后使用跨维度或可逆跳跃马尔可夫链蒙特卡洛(RJ-McMC)方法对所得的地震速度(或其他地震性质)的后验模型分布进行采样。使用基于树的结构表示地震速度模型，在感兴趣的地震性质的小波变换域中进行采样。快速收敛到后验模型的平稳分布，同时需要有限数量的小波系数来定义采样模型。因此，此方法提供了从远距离起始模型的更好收敛性以及量化模型不确定性的能力。通过FWI方法确定的地下速度可用于地震成像。

在本发明的另一方面，为了解决上述问题，一些实施例提供了一种存储一个或多个程序的非暂时性计算机可读存储介质。所述一个或多个程序包括指令，该指令在由具有一个或多个处理器和存储器的计算机系统执行时，使计算机系统执行本文提供的任何方法。

在本发明的又一方面，为了解决上述问题，一些实施例提供了一种计算机系统。该计算机系统包括一个或多个处理器、存储器以及一个或多个程序。所述一个或多个程序被存储在存储器中并且被配置为由所述一个或多个处理器执行。所述一个或多个程序包括操作系统和指令，该指令在由所述一个或多个处理器执行时使计算机系统执行本文提供的任何方法。

附图说明

图1示出了针对同一图像在五个截断级别处的离散小波逆变换(DWT)；

图2示出了DWT的小波系数；

图3示出了跨维度树结构；

图4示出了用于传播主导的研究的模型和噪声数据；

图5示出了DWT域中的树结构和对应的节点位置；

图6示出了跨维度马尔可夫链蒙特卡洛(McMC)采样进度；

图7示出了并行退火；

图8示出了采样统计以及真实模型和反演结果；

图9示出了当使DWT树的级别5对采样模型可访问时的采样统计；

图10示出了通过模型的速度的边际概率密度函数的切片；

图11是后验不确定性的比较；

图12示出了对Marmousi模型的表面反射实验的模型和噪声数据；

图13示出了反射示例的小波系数值；

图14示出了具有并行退火的跨维度McMC采样的进度；

图15示出了Marmousi实验的后验速度的边际分布；

图16还示出了Marmousi实验的后验速度的边际分布；

图17还示出了Marmousi实验的后验速度的边际分布，示出了深度分辨率；

图18是在最大允许的DWT截断级别处的真实模型与均值后验模型的比较；

图19示出了来自随机选择的后验速度模型的模型响应；

图20是图19的放大轨迹；

图21示出了“蝴蝶图”以比较数据与后验匹配；

图22示出了归一化反演残差；

图23示出了“蝴蝶图”以比较数据与后验匹配；以及

图24是示出根据一些实施例的地震成像系统的框图。

在整个附图中，相似的附图标记指代对应的部分。

具体实施方式

以下描述的是提供使用全波形反演(FWI)的地震成像方式的方法、系统和计算机可读存储介质。这些实施例被设计为特别用于地质复杂区域中的地下体积的地震成像。

现在将详细参考各种实施例，其示例在附图中示出。在以下详细描述中，阐述了许多具体细节以便提供对本公开和本文描述的实施例的透彻理解。然而，可以在没有这些具体细节的情况下实践本文描述的实施例。在其他情况下，未详细描述公知的方法、过程、组件和机械设备，以免不必要地使实施例的各方面不清楚。

地下的地震成像用于识别潜在的烃储层。地震数据是在地表(例如地球表面、海洋表面或海底)获取的，作为地震轨迹，其共同构成了地震数据集。地震数据可用于全波形反演(FWI)方法中以确定地下速度，以便可以对地震数据进行正确成像。

有利地，本领域普通技术人员将认识到，例如，本文提供的实施例可用于生成更精确的数字地震图像(即，校正后的数字地震图像)。更精确的数字地震图像可以改善油气勘探并提高油气生产。更精确的数字地震图像可以提供地下的细节，这些细节在传统的地震图像中显示得很差或根本没有显示。此外，更精确的数字地震图像可以更好地描绘出不同特征的开始、结束或其任何组合的位置。作为一个示例，更精确的数字地震图像可以更精确地示出断层。作为另一个示例，假设更精确的数字地震图像指示出烃沉积物的存在。更精确的数字地震图像可以更精确地描绘出烃沉积物的边界，从而可以对烃沉积物进行生产。

本领域普通技术人员将理解，例如，可以在烃勘探和烃生产中利用更精确的数字地震图像来进行决策。例如，可以利用更精确的数字地震图像来选择井眼的位置。本领域普通技术人员将理解，可以基于更精确的数字地震图像来做出以下决定：(a)在何处钻探一个或多个井眼以生产烃沉积物；(b)钻探多少井眼以生产烃沉积物等。甚至可以利用更精确的数字地震图像来选择要钻探的每个井眼的轨迹。此外，与描绘指示出较小的烃沉积物相比，如果描绘指示出较大的烃沉积物，则可以选择更多数量的井眼位置并且可以钻更多数量的井眼。

本领域普通技术人员将理解，例如，可以在烃勘探和烃生产中利用更精确的数字地震图像来进行控制。例如，可以利用更精确的数字地震图像来操纵工具(例如，钻探工具)钻探井眼。可以操纵钻探工具钻探一个或多个井眼以对烃沉积物进行生产。取决于所需结果，操纵工具可以包括围绕或避开某些地下特征(例如，断层、盐底劈(salt diapir)、页岩底劈(shale diapir)、页岩脊(shale ridge)、麻坑(pockmark)、埋藏河道、气烟囱(gaschimney)、浅气穴和塌陷)钻探，钻探通过某些地下特征(例如，烃沉积物)或其任意组合。作为另一示例，可以利用更精确的数字地震图像来控制注入地下、井眼或其任何组合中或从地下、井眼或其任何组合中接收的流体的流动。作为另一示例，可以利用更精确的数字地震图像来控制注入到地下的至少一个烃生产区中或从地下的至少一个烃生产区接收的流体的流动。位于地面或井下的油嘴或井控装置可用于控制流体的流进和流出。例如，更精确的数字地震图像中的某些地下特征可以促使油嘴或井控制装置的激活、去激活、修改或其任意组合，从而控制流体的流动。因此，可以利用更精确的数字地震图像来控制注入速率、生产速率或其任何组合。

本领域普通技术人员将理解，例如，可以利用更精确的数字地震图像来选择井眼的完井、组件、流体等。可以基于更精确的数字地震图像为每个要钻的井眼选择各种套管、油管、封隔器、加热器，防砂筛、砾石过滤层、用于微粒迁移的物品等。此外，可以基于更精确的数字地震图像来选择用于生产烃沉积物的一种或多种采收技术。

简而言之，本领域普通技术人员将理解，在油气行业中需要进行许多决定(例如，在(a)操纵决定、(b)着陆决定、(c)完井决定、(d)在包括但不限于拖缆、海底传感器、VSP、DASVSP以及具有初次和自由表面多次波的成像的工程控制系统和储层监控等的背景下)，并且基于更精确的数字地震图像做出正确的决定应提高安全可靠操作的可能性。为了简单起见，包括井眼位置、用于井眼的组件选择、采收技术选择、控制流体的流动等的许多可能性可以统称为管理地下储层。

本发明包括使用基于树的贝叶斯方法进行FWI的方法和系统的实施例，所述基于树的贝叶斯方法自动选择模型复杂度，从而允许适当的参数化。有限的照明、不足的偏移、噪声数据以及不良的启动模型可能会对地震全波形反演提出挑战。本发明包括基于树的贝叶斯反演方案，其试图通过在使用关于地下结构的适度信息性先验的同时考虑数据不确定性来减轻这些问题。该方法在速度的小波变换域中使用跨维度(trans-D)或可逆跳跃马尔可夫链蒙特卡罗方法对得到的压缩速度的后验模型分布进行采样。这允许快速收敛到后验模型的平稳分布，同时需要有限数量的小波系数来定义采样模型。基于Trans-D树的方法与并行退火一起用于找到正确的崎岖的似然性(rugged likelihood)(即失配)地形应对方法为解决难以优化的大规模地球物理反演问题提供了一种有希望的易于推广的方法，但其中真实模型包含多尺度特征的层次。除了对数字地震成像的改进之外，这种计算机实现的方法在计算效率上也比常规方法高得多。

有源源地震全波形反演(FWI)方法原则上是一个简单的思想。通过最少的处理或手动干预，它的目的不仅是提供地下图像，此外还提供速度模型，当经历正演算子时，该模型“紧密地”匹配观测到的地震场。这需要解决反演问题，其中通过地震波方程支配正演物理学。然而，这种具有有限的接收机覆盖范围和频率带宽的反演问题是非常非线性的，因此很难解决。此外，在不适当的频率处存在噪声会使许多优化方法混淆，并且复杂的地球模型会导致很难以计算有效方式使用的非常高维的模型空间。所暗示的非线性表现为局部失配极小值，从而导致模型在FWI看来不是最佳收敛或“周期跳过”的模型。存在多种改善收敛性的有前途的方法，诸如估计时移以最小化初始建模数据和观测数据之间的运动学差异、使用扩展模型域和/或非局部波物理学。另一种方法是解决一系列约束的局部凸起子问题。还有其他方法试图通过使用最佳传播距离、通过向数据添加人工低频、迭代使用Wiener滤波器或使用二次罚分法来改善失配函数的凸度。所有这些方法的一个共同点是努力使优化算法更易于找到正确的失配地形应对方法。所有这些方法在不同程度上都可以在不同情况下正常工作，但不能保证收敛。此外，考虑到所涉及的各个步骤，这些方法不容易适用于解决方案评估或不确定性估计。当这种解决方案本身不容易被找到时，本发明量化了我们为FWI问题提供解决方案的可信度(在贝叶斯意义上)。此外，该算法会自动选择速度模型的有限的一组离散小波变换系数并使用该有限的一组离散小波变换系数进行操作。与正演建模有限差分网格中的单元相比，这导致减少数量的未知数，从而容许在2D和潜在的3D FWI中，以最小的先验假设来进行易处理的不确定性估计。

在大多数常规的地球物理反演方案中，模型网格的几何结构是固定的，也就是说，在反演期间，单元的大小及其数量是不允许变化的。传统上，解决方案着重于最小化以下目标函数：

其中m是模型向量，d是观测数据，并且f(m)提供由于m而导致的正演建模预测。λ²是正则化参数，R是一旦应用到m就会在p范数中产生被认为感觉保持较小的长度的度量的任何算子。(1)中的第一项是数据失配(由数据精度W加权)，并且第二项是设计用于使模型(或模型与优选模型的偏差)保持较小的正则化项。两者之间的折衷由所谓的Tikhonov正则化参数λ²控制。这类似于脊回归的统计技术，也就是说，取决于λ²的值，对于线性问题和p＝2范数，(1)的解位于数据失配项的最小值与模型长度项的最小值之间的“脊”，以便同时最小化两者。显然，需要对算子R、赋予模型长度的权重以及模型范数的选择进行选择。涉及梯度下降或结合Tikhonov正则化使用雅可比矩阵(或其逆)的非线性最小二乘FWI解很容易概念化，很好理解，但是如果R选择不当，则收敛速度会很慢。选择较少数量的参数，或将p＝1范数与稀疏模型“框架”结合使用，可以消除某些超参数选择。

当然，使用具有较小长度度量的稀疏模型表示不仅有助于FWI收敛到(1)的合适解，还可以进行关于简约模型参数化的另一观测-更简单的理论(我们模型案例)提供更清晰的见解。这是Occam的反演使用的方法，其目的是产生与数据噪声兼容的最平滑模型或最稀疏模型。但是，这些模型是极值模型，并且不应视为真正代表地球。明智地，我们应该考虑适当简单，但也适合数据的模型。从统计上讲，这称为模型选择问题。目的是避免产生具有低方差但具有高偏差的简单模型，或避免产生具有高方差但具有低偏差的复杂模型。对于地球物理反演理想地是，我们应该对不是一个，而是与我们的数据以及我们先验的地球及其复杂度概念兼容的一系列模型进行采样。

在到目前为止概述的方法中，目标是找到(1)的最小值，希望它是全局最小值。如前所述，不存在这种收敛保证。此外，即使要找到全局最小值，也不会排除存在具有符合数据噪声的类似失配的其他模型。这些模型可能会表现出非常不同的速度结构，这是非线性问题的典型特征。继续前面提到的地球物理思想，希望使用一系列超参数(例如正则化参数、单元数等)、一系列模型进行采样，以使得模型本身具有适当的复杂度，具有符合测井数据、露头和实验室实验的地震速度，同时与噪声地震观测兼容。由于大量的可能性，通过反复试验手动执行此操作将被证明是不可能的。但是，即使使用系统的方法也很复杂，因为由于超参数和反演模型的每种组合，仍然需要对结果进行定量地加权。

本发明通过在贝叶斯意义上重新检查(1)来完成该任务。概括地说，对于每个采样模型m，失配项提供了模型的似然性的度量，而模型向量的长度概述了我们对模型的先验知识，包括模型的复杂度。

更严格地说，贝叶斯公式是

p(m|d)∝p(d|m)p(m), (2)

为了我们的目的，最好从右到左如下阅读：p(m)是m的先验概率，我们独立于观测值d知道它。我们通过进行向我们展示m与观测值相符的可能性的地震实验重新评估我们对m的先验概念。该权重由似然函数p(d|m)给出。通过似然性对我们的先验概念进行加权或更新的结果提供了观测模型m的后验概率。后验概率用项p(m|d)表示。然后，我们对我们的先验概念可以接受的各种模型重复此过程，直到我们获得由概率密度函数或PDF p(m|d)表示的模型的集合。因此，我们可以将具有许多可能解的优化问题(1)转换为采样问题(2)。本领域技术人员将注意到，(2)缺少确保其积分为整数的归一化常数，因此并不是真正的概率密度函数。实际上，(2)可以更好地表示多维直方图，直到我们通过对右侧的所有模型进行积分来对其进行归一化：

p(d)＝∫_mP(d|m)p(m)dm， (3)

其中p(d)被称为证据。但是，对于模型评估，我们仅对各种模型的相对概率感兴趣。因此，出于我们的目的，我们可以使用(2)来采样到比例常数。重要的是要注意，我们的(2)中的先验包括关于各种复杂度的规范(包括具有不同数量变量的参数化)，因此p(d)是“总”证据。

对于优化问题(1)，如对任何地球物理问题都适用，从许多不同的角度出发，模型正则化都是必须的，无论是为了提高矩阵逆的稳定性，还是保持模型波动(或更新)较小，或使模型(或更新)接近优选模型。但是，用于描述模型的参数数量(模型复杂度的度量)也可以视为未知进行采样，而无需明确进行正则化。通过这种方法，我们不仅考虑用于描述数据的最简单或最平滑的模型，还考虑具有与观测到的数据兼容的不同数量参数的模型的集合。基于出生/死亡蒙特卡罗的跨维度反演方法和更通用的可逆跳跃马尔可夫链蒙特卡罗(RJ-McMC)方法可以完成此任务。对于一维模型，这将意味着在可变数量的层上进行采样。对于二维模型，具有不同数量单元的Voronoi单元已被广泛使用。实际上，基于贝叶斯定理的trans-D算法就模型的复杂度执行模型选择的任务。模型既不是过拟合也不是不足拟合的事实是基于贝叶斯简约性的思想。当适当地制定公式时，贝叶斯定理的框架会内置“Occam因子”，该因子会对过于复杂的模型进行惩罚。为了检验该论点，我们注意到，将trans-D模型矢量定义为m＝[m_k，k]，其中m_k是具有描述压缩速度的k个参数的模型(对于本发明的FWI应用而言)。可以从数据和模型的联合概率得出：

将总证据p(d)视为常数，我们得到

(5)左侧的项是推断参数k的数量的后验概率(进行实验后)。右侧的第一项是k个参数充分拟合数据的似然性。为了检查右侧括号中的第二项，我们首先从联合概率和条件概率的定义中注意到p(m_k，k)＝p(m_k|k)p(k)。因此，右侧括号中的项是k参数模型的先验模型概率与后验概率之比。参数k的数量越多，先验概率的分布就越薄(即更低)，因为先验PDF需要在更大的体积上积分为1。由于可接受的k参数模型占用后验的少量先验空间，因此k参数的后验概率通常比先验概率更高(即更尖)。因此，使用的参数k越多，括号中的分数就越少。但是，我们使用的参数k越多，k参数拟合的似然性就越大。在trans-D公式中，括号中的系数和数据似然性相互折衷，从而在很大程度上取决于数据，自动提供类似于正则化的解决方案。由于p(k)的概率相同，并且在(Ray等2016)中讨论了一些简化的假设，括号中的分数可以被解释为后验可及体积与先验可及体积的比率，有时称为“Occam因子”。此公式允许反演自参数化以达到良好效果，从而在具有更好数据覆盖率和低噪声的区域中提供更高的模型分辨率。

既然我们已经解释了trans-D公式，以免出现(5)的右侧取决于m_k，而左侧则不取决于m_k，我们可以再次简单地使用条件概率的定义来验证(5)的右侧等于p(k，d)。这与(4)完全一致，因为根据定义，p(k|d)＝p(k，d)/p(d)。对于地震面波断层成像应用，Trans-D优于基于使用B样条的子空间变换的反演。已知基于评估经由边际似然性p(d|k)的不同参数化的证据或给定假设的证据(在我们的情况下为k参数模型)的trans-D公式的替代方法。但是，这涉及为每个k参数查找寻找证据的计算上禁止的任务，并且仅对某些种类的地球物理反演是可行的。

对于勘探地震FWI问题，由于正演评估的巨大计算成本，最近才提出表征完全非线性不确定性的解决方案。一些方法使用基于随机源二次采样的贝叶斯解决方案，但在假设关于最大后验(MAP)模型的高斯分布的同时，利用固定参数化。其他人则将遗传算法(GA)与使用邻域算法的模型重采样结合起来，然后再进行Gibbs采样。他们使用两个网格方法，对反演模型使用粗略方法，并且对正演模型使用精细方法。但是，数据不能确定反演模型网格的粗糙度，并且估计的不确定性的质量还取决于从GA到NA+GS算法的输入集合。其他方法提出了一种两网格方法，尽管反演模型网格是固定的，但该方法涉及算子放大。所有这些方法都有望量化地震FWI的不确定性，但并未解决模型选择问题。我们了解到的唯一尝试使用tran-D反演的尝试是针对垂直地震剖面(VSP)反演问题和弹性FWI问题，但两者都假定横向不变的地球模型，而这当然不代表为油气勘探和生产目的而必须获得的真实地球模型。理论上，针对一维和二维地球模型演示的贝叶斯模型选择原则同样适用于三维反演。然而，正如Hawkins和Sambridge(2015)所指出的那样，对于三维中的trans-D问题而言，计算上高效的参数化并不容易构建，并且很难包含有关几何结构的先验知识。

Hawkins和Sambridge(2015)的最新工作表明，可以由基于树的结构表示的任何基函数集都可以用作tran-D反演的有效模型表示。使用此公式的主要优点是，从理论和实践效率的角度来看，它都与地球模型的空间维数无关，无论是一维、二维还是三维。在实施例中，我们具体地使用小波基函数和离散小波变换(DWT)，其易于服从基于分层树的表示。具有适当基函数集(例如CDF 9/7)的小波变换通常用于压缩图像信息(例如JPEG 2000)。这将使变换域对简约的地球物理模型表示具有吸引力，这将利用综合示例得到证明。如前所述，曲线波或小波基函数集已用于勘探地震FWI，但采用了优化设置。正如Hawkins和Sambridge(2015)所讨论的，未完全填充的有效小波树可以表示从低空间波数到高空间波数的特征层次。结合trans-D算法，这提供了一种多分辨率方法，其根据观测到的数据进行自适应参数化。已经进行了自适应反演网格划分，但是这些方法使用固定的准则进行自适应，而不是对一系列参数化进行采样，在参数化中，模型复杂度由数据决定。这种基于trans-D树的方法的最近成功应用可以在Hawkins等(2017)进行机载电磁反演、Dettmer等(2016)量化海啸海面位移的不确定性、Hawkins和Sambridge(2015)进行二维和三维地震断层成像中发现，而本发明是该方法首次用于地震FWI。

对于一维、二维和三维模型，树表示需要分别使用修改的二叉树、四叉树和八叉树结构。对于小波变换域中的所有这些表示，第一节点系数(在树的最高层)表示模型中速度的平均值(将提供给有限差分算子)。该节点在第二级分支为1、3和7个节点(再次，分别对于一维、二维和三维模型)，该级的系数代表波长约为整个模型的长度比例的一半的基函数的强度。从此级向下，每个节点分支成纯二叉树、四叉树或八叉树，其中每个子级精确地具有2、4和8个子级。树的深度受正演模型网格的大小限制。每个连续的深度级别(在逆小波变换域中)代表速度模型中更精细的缩放特征。在这里介绍的所有工作中，我们在二维中工作时都使用了经过修改的受限四叉树，但是本领域技术人员将认识到，通过使用适当的树结构，这些方法同样适用于一维或三维中。

使用基于树的小波变换表示的另一个优点是，取决于当前的问题，可以使用不同的小波基。对于传播主导的问题，平稳基函数(诸如CDF 9/7)可能是合适的。对于反射主导的问题，可以使用诸如Haar小波的尖基函数。图1显示了128×128像素图像，在变换域中的五个截断级别(level)处，使用这两种基函数逆变换回图像域。对于128×128的正方形图像，级别8对应于无截断，如2^level-1＝128。使用离散小波变换(DWT)的局限在于所有维数必须是2的幂。虽然我们在这项工作中使用平方速度模型，但Hawkins和Sambridge(2015)已经展示了如何通过将不止一个根节点用于小波树模型来将DWT用于矩形模型域。图2显示了在小波变换域中使用CDF 9/7在全小波变换和级别5的截短版本之间进行的比较。级别5表示需要最多16×16的系数中少量为非零，同时在图1的第5行第1列中提供近似值。

通过trans-D McMC算法对后验模型PDF(2)进行采样，其细节在下面提供。特别地，我们根据先验规范和似然函数添加节点、删除节点或修改节点系数来对不同的小波树进行采样。

我们从一个非常简单的模型开始该算法，通常是只有一个根节点的树。然后，我们允许算法将活动节点迭代地添加到树(“出生”)，修剪它们(“死亡”)，或者简单地修改现有活动节点上的系数值(“更新”)。由于数据可能通过接受概率α要求而完成所有这些操作。重复此过程，直到McMC链收敛到固定的样本链为止。Ray等(2016)详细介绍了trans-D反演的收敛监测和用于避开局部似然最大值(失配最小值)的并行退火算法的细节。

遵循Hawkins和Sambridge(2015)的表示法，我们需要保持跟踪活动节点的集合S_v，要生出的节点的集合S_b以及死亡的没有子级(“树的叶子”)的活动节点的集合S_d。图3显示了示例树模型，其中有k＝2个活动节点，并说明了活动、出生和死亡集合。

在McMC的每个步骤中，以相等的概率随机选择三个可能的动作之一：更新节点系数、出生或死亡。对于节点更新，从节点的集合S_v中随机选择节点，并使用高斯建议对系数值进行扰动。通常，我们将更新的标准偏差设置为特定节点深度处的统一边界宽度的5％。此举不会改变模型维度。

出生动作涉及以下步骤：

如果k<k_max

(1)复制初始树模型m(即系数值和三个节点集合)的副本m’；

(2)从初始模型的出生集合S_b中随机选择要激活的节点；

(3)从建议的模型的出生集合S’_b中删除所选节点；

(4)针对所选节点的深度级别，从统一的先验系数范围中统一提出系数值v’；

(5)将所选节点添加到建议的模型的活动集合S’_v中；

(6)将所选节点添加到建议的模型的死亡集合S’_d中；

(7)除非父级是根节点，否则从死亡集合S’_d中删除选定节点的父级(如果父级在死亡集合中)，因为父级不再是叶节点。

(8)找到所选节点的子级，然后将它们添加到建议的模型的出生集合S’_b中(如果子级在最大树深度限制内)。

该动作建议增加维度k’＝k+1。

死亡动作涉及以下步骤，并且与出生步骤相反：

如果k>k_min，

(1)复制初始树模型m(即系数值和三个节点集)的副本m'；

(2)随机选择树节点从初始模型的死亡集合Sd中删除；

(3)为选定的节点系数分配零(仅出于完整性考虑)；

(4)从建议的模型的死亡集合S’d中删除所选节点；

(5)从建议的模型的活动集S’_v中查找并删除所选节点；

(6)从出生集合S’_b中查找并删除所选节点的子级(如果子级在深度限制内)

(7)将所选节点添加到建议的模型的出生集合S’_b中；

(8)除非所选节点的父级是根节点，否则将其添加到死亡集合S’d中(如果它现在是叶节点)。

该动作建议减少维度k’＝k–1。

McMC链从模型m到m'的概率由接受概率α给出。对于基于树的trans-D McMC，三种不同动作类型均采用不同的形式，并且下面给出的表达式由Hawkins和Sambridge(2015)详细推导。

对于更新动作，维度没有变化，并且当如我们所做的那样从统一的先验系数范围提出时，它仅仅是似然比：

对于出生动作，接受概率为

其中|S_x|是集合S_x中的元素的数量，h是最大深度级别限制。对于死亡动作，接受概率为

如果关于节点数的先验概率是均匀的，则

但是，如果像实施例中那样使用了Jeffrey的先验，则

以及

如果建议的模型以概率α接受，则将其存储为下一个样本。如果建议被拒绝，则McMC链中的先前模型将被保留为下一个样本。

从概念上讲，该算法最困难的部分是在给定为出生和死亡建议计算α所需的活动节点数k的情况下，对树的可能排列的数目进行计数。对于二叉树，如果有n个节点，则对于节点i，例如我们可以在其前的节点具有C_i-1排列。这使得剩余节点可能有C_n-i排列。由于排列是独立的，因此节点i的排列的总数为C_i-1·C_n-i。但是由于有n个节点，因此我们必须对所有i求和，因此n个节点的排列的总数为

对于n＝1，我们设置C₀＝1，因为只有一种方法可以制作只有1个节点的树。这通过递归关系定义了加泰罗尼亚数列(Catalan number sequence)，基本情况定义为C₀＝1。人们可以使用此逻辑来构造更高阶和更通用树的排列数(Hawkins和Sambridge 2015)。C_n可以很容易地通过递归来求解，但是仔细研究我们可以发现，要获得C₃，我们需要计算C₂和C₁。但是，如果我们已经计算了C₂，则可以存储此值并重新使用它，而无需进行其他递归调用。这被称为备忘录，是一种在动态编程中广泛使用的技术。这在进行许多递归调用时很有用，如在纯四叉树的情况下，其中因此可以将排列数Y_n写为：

除了对Yn进行记忆外，我们还可以对j和k上的每个部分和进行记忆，因为部分和是求和上限的函数。笛卡尔DWT所需的修改后的四叉树具有一个根节点和三个子级，这三个子级中的每一个都遵循纯四叉树结构。因此，利用我们可以再次记忆部分和的事实，我们可以将排列数写为：

最后，我们可以使用表示深度级别限制的另一个索引来处理限制的树深度。对于二叉树的情况，对深度h的限制通过以下表达式给出：

与

我们可以将完全相同的受限二叉树排列逻辑应用于修改后的受限四叉树排列计数。我们需要做的就是简单地使计算依赖于先前的级别h-1来修改任意级别h的排列数。

对于加性噪声，通过中心限制为渐近高斯(尤其是在频域中，如Ray等2016所示)，我们将模型似然函数L(m)＝p(d|m)定义为

其中C_d是数据错误的协方差矩阵。由于DWT是线性变换，因此我们可以写

f(m)＝F(Hm), (7)

其中F是地震正演算子，H是DWT逆算子，m是由小波树上的系数值表示的模型矢量。换句话说，Hm是馈入可变密度、声学和各向同性有限差分引擎的二维速度模型。假定源签名是已知的，或者可以将其作为模型的函数作为最大似然估计来推导，如Ray等(2016)所示。

在这个实施例中，我们仅关注地球中的速度的变化，假设密度变化是已知的或密度没有变化。如本领域技术人员将认识到的，这不是方法的限制，该方法容易概括为更多变量。现有模型需要指定树上的节点的概率。因此我们可以写

p(m)＝p(v,T,k), (8)

其中v是小波变换域中的速度的矢量(在此实施例中)，需要注意的一点是，这使基于树的公式与基于层或基于单元的trans-D有所不同。T是一种特殊类型的小波树(对于我们的2D应用，是经过修改的受限四叉树)，而k是表示有效树结构的活动节点数。使用概率的链式规则，我们可以编写：

p(v,T,k)＝p(v|T,k)p(T|k)p(k),

该乘积下的最后一项假设在给定特定树类型T的k个活动节点的情况下，每个节点处的小波系数彼此独立。Hawkins和Sambridge(2015)和Dettmer等(2016)简单地在每个节点位置使用宽的统一边界。但是，如图3所示，这些系数值跨越多个数量级，但是在特定的深度级别，这些值都在有限的范围内。详细地说，对于大多数自然发生的图像，树的顶层(代表较粗的特征)的值通常比离原点较远的深度级别(代表较细的特征)的值更极端。这完全类似于大多数自然图像的傅立叶光谱，这些傅里叶光谱包含更强的低波数内容而不是高波数内容。p(k)只是节点数的先验。它可以是恒定的(即均匀的)，也可以是与活动节点的数量成反比的Jeffrey的先验(Hawkins和Sambridge，2015)。Hawkins和Sambridge(2015)的关键发展是p(T|k)的定义。如果我们假设给定k个活动节点，则具有从顶部的根节点到指定的最大深度的活动节点的所有有效树结构都是同等的概率-这是一个合理的假设，因为这将暗示任何速度模型都将具有粗略以及精细尺度的特征-然后要定义此概率，我们需要对此类有效树的排列数N_k进行计数。概率简单地由以下表达式给出

常规方法无法对修改后的受限树的排列数进行计数。对于一般的受限树，Hawkins和Sambridge(2015)中提出了一种有效的递归方法来计算N_k。在本发明中，我们提供了一种用于二维小波树结构的不太通用的可能更容易实现的有效递归伪代码。对于DWT的一维或三维小波树，可以轻松对其进行修改。

在另一个实施例中，获得后验模型PDF需要使用Metropolis-Hastings-Green算法进行采样(2)。接受或拒绝模型建议的标准由以下概率给出

其中q(m’|m)是从模型m步进到m’的提议概率，并且|J|是改变维度时变量变换的雅可比行列式。对于这种情况下使用的出生-死亡算法，它计算为一(unity)。为了避开局部失配最小值(似然最大值)，使用并行退火算法在不同的“温度”τ并行运行各种交互的McMC链。使用无偏τ＝1链进行后验推断。先前提供了采样方法和模型建议的详细信息。

在用于传播主导用例的实施例中，模型和噪声合成数据如图4所示。将62个接收器以20m的间隔放置在地面上，将两个源放置在模型的边缘处的1260的深度处。该模型为网格间距为10m的128×128单元。源是以5Hz为中心的Ricker小波。在所有迹线上都添加了0.5％的最大炮点振幅的不相关的高斯噪声。在该实施例中未示出，但是在本发明的范围内，对于真实世界带通时域数据的相关噪声的存在，将需要使用修改的(6)中的似然性，其中数据协方差中具有对角项。

选择CDF 9/7基进行反演，因为它在级别5截断时比Haar基提供了更低的χ²失配(见图2)。根据Hawkins和Sambridge(2015)，将p(vi|T，k)的先验边界设置为统一有界。如本节所述，我们小心不要过度指定界限。参照图5的二维图像，级别1对应于树的根节点，具有编号为1的一个系数。级别2具有编号为2-4的三个(根的)子节点。从级别2开始，树遵循四叉结构，每个节点都有4个子级。因此，级别3包含编号为5-16的节点。最后，级别4包含级别3中所有节点的4个子级的每个子节点，编号为17-64。在每个级别上都找到了真实模型的最小和最大小波系数，并且该级别上所有系数的界限都设置为小于和大于极值的2％。与所有贝叶斯方法一样，可以将先验规范的必要性看作是祝福和诅咒。如果对地球结构和地球中可能的速度变化一无所知，则此方法将无用，但所有地球物理反演问题都需要做出一些约束性假设，并且这并不是我们方法的唯一局限性。但是，如果我们对可能的结构有所了解，那么实际上可以通过以这种方式设置先验界限来量化此解释性方面。对于第二个综合示例，示出了使用我们的方法的示例先验模型实现。变换域为逆变换域中的概念性地质模型(图像)提供了一种非常简练的方法来指定压缩采样界限。逆变换域是我们习惯于思考的域。节点可以方便地用以所谓的“Z”或“Morton”顺序编号的线性树阵列表示，这同样适用于三维压缩。阵列索引值遵循自相似的Z模式。二进制交织将线性索引转换为变换域图像中其适当的行和列位置。这里需要提请注意-来自小波树模型的逆变换图像可能包含非物理特征(诸如负速度)，因此将这些模型的先验概率设置为零非常重要。使用此方法必须使用具有找到困难地形的正确应对方法的机制的随机方法，因为迭代优化方法可能会陷入由无限高脊(即零后验概率)分隔的两个区域之间的目标函数格局中。我们的解决方案是为此目的使用并行退火。

该算法非常快地减少了失配，直到在仅仅400次迭代中达到接近2的RMS(均方根)失配(图6)。随着失配的减少，模型的复杂度也会增加。由于我们使用并行退火，因此在不同温度下的每个马尔可夫链均以不同的颜色表示。仅从τ＝1处的链进行后验推断。通过构造，并行退火可确保较低温度的链始终包含最低的失配，而较高温度的链避开失配的可能性较小(即较高)，从而可以避开局部失配最小值(Ray等2016)，如图7所示。将43条平行相互作用的McMC链以1到2.5之间的对数间隔温度运行，以对(11)中的似然函数进行退火。尽管从失配减少和采样模型中可以很好地快速获得拟合良好的模型(图6和7)，但是要获得模型不确定性的估计值，我们需要更长的采样，并且RMS失配保持在1.18左右(大约20000的χ²/2)，在大多数度量下，都很好地表明了收敛性(图9)。从该图可以看出，平均采样模型非常接近真实速度模型。但是，活动节点的数量经常达到64，这是对于将树深度限制为级别4所允许的最大数量。这意味着数据需要更复杂的参数化。

当我们允许小波树模型占据级别5时，对于总共256个可能的活动节点，我们采样的时间更长，并达到图9中描述的情况。RMS从1.18下降到1.004，并且尽管从未超过150个，但现在采样的系数数量增加到140个。我们可以找到与这些树模型中的每个相对应的速度模型，并且不是像图7那样计算均值，我们现在可以计算在每个地下位置的速度的边际PDF，并在图11中将其显示为“概率立方体”。真实速度模型以黑色显示，与通过概率立方体的切片一致，其中任何位置的深色表示立方体中该点处的位置和速度的较高概率。异常的边缘似乎已经很好地解决了，速度整齐地聚集在一起，但是异常区域的中心却没有，这可能是因为该点缺少照明。还要注意其中可能有一个以上的速度的某些空间位置处的多模态。速度折衷也是可见的，其中沿传播路径的较低速度导致沿传播路径的不同位置处的较高速度。如果我们决定在级别4结束采样，尽管数据拟合度稍差，但我们将获得更加乐观的不确定性图景。图11提供了这两个级别的不确定性的比较。该图再次示出了反演期间所做的选择如何影响我们对地下的结论。

在另一个实施例中，该方法可以应用于表面反射问题。该示例基于Marmousi模型的缩放版本。它是128×128像素，网格间距为20m。假设已知，源小波是峰值频率为3.75Hz的Ricker。对两个炮点进行了建模，将0.2％的最大振幅的不相关的高斯噪声添加到所有迹线中。模型和噪声数据(减去直达波)如图12所示。

类似于先前的实施例，通过找到每个级别上的最小和最大DWT系数，并以在这些界限的上方和下方2％的值进行操作，来获得先验界限。图13示出了先验的一些5节点实现。我们在本例中使用Haar基集合，因为在我们的试验中光滑的CDF 9/7基不能令人满意地工作-我们推测这是因为反射所需的边缘比较低级别的CDF小波系数能够提供的边缘更锐利。如果通过增加更多活动节点不会显着降低失配，则贝叶斯简约将不鼓励对更复杂的树进行采样。利用Haar基，我们从200到10000个模型(RMS 2.3到1.44)获得快速收敛到类似于背景速度的模型，这取决于“可接受的失配”的概念。我们应该在这里提到，将水速固定为真实值且背景速度恒定为2.6km s-1的高斯-牛顿法的简单实现无法提供更新。在选定的迭代中，在目标链(τ＝1)中的采样过程和采样的模型如图14所示。最初使用80条并行退火McMC链，对数间隔温度在1和5之间。经过200000次迭代，我们可以合理地相信已经避开了局部最小值(似然最大值)，并且仅使用前40个链(温度从1到2.213)来对后验模型PDF进行采样。1000次迭代后，失配级别渐进到RMS 2.0，允许的树深度最大值设置为级别4(最大64个节点)。在允许算法访问级别5(最大256个活动节点)之后，失配在接近预期值1的约1.37的RMS再次渐进。但是，采样的节点数接近256，并且很显然，如果要达到RMS 1.0，则至少必须使下一个深度级别对算法可用。当模型可以访问具有1024个最大节点的级别6时，在大约200000次迭代达到了非常接近1的RMS。然后允许采样再进行一百万次迭代，并且没有模型需要超过468个活动节点。

对于后验推断，我们仅使用最后的700000次迭代来从不受拟合差的模型的影响的稳定的马尔可夫链获得样本。仅目标链用于后验推断。与前面的示例相似，我们可以创建具有每个地下位置处的速度的边际PDF的概率立方体，其结果如图15-17所示。同样，在左列中，较深的颜色表示速度在立方体中的特定点处的较高的可能性。真实速度曲线用红线显示，并且每个深度的90％可信区间在两条黑线之间。最佳的速度分辨率似乎在照明源附近(图15)，向中心变差(图16)。正如预期的，分辨率会更加浅(图17)。深度超过1.5公里时，速度的PDF过于分散，无法提供有意义的信息。令人鼓舞的是，在大多数情况下，真实速度在5％和95％的可信区间内，并且当速度变化的PDF与距离和深度成批时，可以推断出速度变化。出现后验的分辨率图片与在边缘处有两处炮点的采集设置、表面记录和高水平的噪音一致。我们应该注意，采样的后验模型会自适应地参数化以提供这种分辨率的图片-从而仅在数据能够提供的地方提供详细的细节。我们还可以提供后验统计信息(诸如每个像素处的均值(图18)和分位数速度)，但在我们看来，显示带有深度的速度的边际后验PDF(图15-17)更为有用。所有这些结果都可以通过包括图形用户界面(GUI)的用户界面呈现给用户。

在任何反演之后的重要检查是对数据拟合和残差的检查。利用真实数据，相关残差表示理论误差、不正确的似然函数、相干噪声或以上各项的某种组合。不能总是避免这些问题，但是残差可以告诉我们可以信任多少反演，例如，在Ray等(2016年)中，预期残差将是相关(由于加工/采集伪像)但高斯的，并且确是这样。对于本文的合成示例，我们添加了不相关的高斯随机噪声，并期望因此我们的残差也应为不相关和高斯的。对于反射实验，我们从后验PDF中选择了100个随机模型，并正演计算了炮点道集。我们已绘制了在选择的迹线上的所有100个建模响应，如图19所示。放大到其中一条迹线(如图20所示)，我们可以看到添加的高斯噪声如何使允许的模型响应散布开来，从而对反演速度的不确定性做出贡献。

我们可以检查两个炮点的100个随机后验模型的数据拟合，如图21所示。在左侧是通过对每个炮点的后验模型响应求和而计算出的均值地震响应。右边是噪声合成数据。可以看到，模型响应的均值与观测到的数据非常相似。所有主要事件及其幅度与偏移(AVO)特征、多次波和折射大多数部分都已得到很好的再现。对于两个炮点，来自后验集合中的相同的100个随机模型的所有时间样本的归一化反演残差如图22所示。这进一步证明了采样/反演按预期工作。我们假设了高斯似然，并且采样的模型并没有过拟合数据，从而产生如下残差，这些残差通过其标准偏差归一化后，近似于分析标准正态PDF。我们还可以将模型响应的均值与真实的、无噪声的合成数据进行比较，如图23所示。

我们已经用两个合成的示例证明了在基于树的框架中进行具有自适应模型复杂度的完全非线性的二维贝叶斯反演的可行性。这样做有很多优点，其中最主要的是易于使用的参数化，其在一维、2维和3维地球模型中同样有效。使用基于树的参数化，我们可以容易地获得高达25％的出生和死亡接受率，从而确保McMC的良好混合，这对于Voronoi单元参数化而言非常困难(Hawkins和Sambridge 2015)。如此处所做的那样，指定先验系数界限将先验模型限制为仅在一定范围内的可行模型之内，而不是过于严格的约束条件。并行退火的使用使我们能够避开作为基于反射的FWI的主要障碍的局部失配最小值。最后，DWT提供了一种容易切换到最适合解决当前问题的模型基(model basis)的方式。当然，使用贝叶斯先验和不同的基函数中存在固有的主观性(Hawkins和Sambridge 2015)。但是，出于实际目的，几乎所有通过优化进行的地球物理反演都利用了合理的约束条件。如此处所示，贝叶斯反演方法自然可以将多个结构约束条件并入为先验信息。尽管毫无疑问，贝叶斯评估要比优化花费更多时间，但在地球物理学尤其是统计界中，都在积极研究将采样速度提高一个数量级的快速方法，同时越来越容易从商业供应商获得并行计算的可用性。在这种情况下，我们的分析可以扩展到更高的频率和更多的炮点。地球物理数据的贝叶斯反演提供不确定性分析这一事实非常宝贵，因为它可能是许多由地球物理数据决定的决策的风险缓解因素。

图24是示出根据一些实施例的地震成像系统500的框图。尽管示出了某些特定特征，但是本领域技术人员将从本公开内容中意识到，为了简洁起见并且以免混淆本文所公开的实施例的更多相关方面，没有示出各种其他特征。

为此，地震成像系统500包括一个或多个处理单元(CPU)502、一个或多个网络接口508和/或其他通信接口503、存储器506以及用于将这些以及其他各种组件互连的一个或多个通信总线504。地震成像系统500还包括用户界面505(例如，显示器505-1和输入装置505-2)。通信总线504可以包括互连并控制系统组件之间的通信的电路(有时称为芯片组)。存储器506包括高速随机存取存储器(诸如DRAM、SRAM，DDR RAM)或其他随机存取固态存储器装置；并且可以包括非易失性存储器，诸如一个或多个磁盘存储装置、光盘存储装置、闪存装置或其他非易失性固态存储装置。存储器506可以可选地包括远离CPU 502的一个或多个存储装置。在存储器506内包括非易失性和易失性存储装置的存储器506包括非暂时性计算机可读存储介质，并且可以存储地震数据、速度模型、地震图像和/或地质结构信息。

在一些实施例中，存储器506或存储器506的非易失性计算机可读存储介质存储以下程序、模块和数据结构或其子集，包括操作系统516、网络通信模块518和地震成像模块520。

操作系统516包括用于处理各种基本系统服务并执行与硬件有关的任务的过程。

网络通信模块518促进经由通信网络接口508(有线或无线)和一个或多个通信网络(诸如，因特网、其他广域网、局域网、城域网等)与其他装置的通信。

在一些实施例中，地震成像模块520使用基于树的贝叶斯方法执行包括FWI的地震成像方法的操作。地震成像模块520可以包括数据子模块525，其处理包括地震道集525-1至525-N的地震数据集。该地震数据由数据子模块525提供给其他子模块。

FWI子模块522包含一组指令522-1并且接受元数据和参数522-2，其将使其能够执行用于全波形反演的操作。贝叶斯子模块523包含一组指令523-1并且接受元数据和参数523-2，这将使其能够执行用于FWI方法的基于树的贝叶斯方法。成像子模块524包含一组指令524-1并接受元数据和参数524-2，这将使其能够使用由FWI子模块522和贝叶斯子模块523确定的速度执行地震成像。尽管对于本文讨论的子模块已经标识了特定操作，但这并不意味着是限制性的。每个子模块可以被配置为执行被标识为其他子模块的一部分的操作，并且可以包含其他指令、元数据和参数，以使其能够执行用于处理地震数据并生成地震图像的其他操作。例如，任何子模块可以可选地能够生成将被发送到用户界面显示器505-1并在其上显示的显示。另外，地震数据或经处理的地震数据产物中的任何一个都可以经由通信接口503(一个或多个)或网络接口508来发送，并且可以存储在存储器506中。

方法100可选地由存储在计算机存储器或非暂时性计算机可读存储介质(例如，图24中的存储器506)中的指令支配，并且由一个或多个计算机系统的一个或多个处理器(例如，处理器502)执行。该计算机可读存储介质可以包括磁盘或光盘存储装置、诸如闪存的固态存储装置或其他一个或多个非易失性存储装置。存储在计算机可读存储介质上的计算机可读指令可以包括以下一种或多种：源代码、汇编语言代码、目标代码或由一个或多个处理器解释的另一种指令格式。在各种实施例中，可以组合每种方法中的一些操作和/或可以从图中所示的顺序改变一些操作的顺序。为了便于解释，方法100被描述为由计算机系统执行，尽管在一些实施例中，方法100的各种操作分布在分立的计算机系统上。

尽管以上描述了特定实施例，但是将理解，其并不旨在将本发明限制于这些特定实施例。相反，本发明包括在所附权利要求书的精神和范围内的替代、修改和等同物。阐述了许多具体细节以便提供对本文提出的主题的透彻理解。但是对于本领域的普通技术人员将显而易见的是，可以在没有这些具体细节的情况下实践本主题。在其他情况下，未详细描述公知的方法、过程、组件和电路，以免不必要地使实施例的各方面不清楚。

在本文的本发明的描述中使用的术语仅出于描述特定实施例的目的，并且不旨在限制本发明。如在本发明的说明书和所附权利要求书中所使用的，单数形式“一”、“一个”和“该”也意图包括复数形式，除非上下文另外明确指出。还应理解，本文所用的术语“和/或”是指并涵盖一个或多个相关联所列项目的任何和所有可能的组合。将进一步理解的是，当在本说明书中使用时，术语“包括”、“包含”、“囊括”和/或“含有”执行了所述特征、操作、元件和/或组件的存在，但是不排除存在或增加一个或多个其他特征、操作、元件、组件和/或其组。

取决于上下文，如本文中所使用的，术语“如果”可被解释为意指“当……时”或“在……时”或“响应于确定”或“根据确定”或“响应于检测到”所述先决条件为真。类似地，取决于上下文，短语“如果确定[所述先决条件为真]”或“如果[所述先决条件为真]”或“当[所述先决条件为真]时”可解释为“在确定……时”或“响应于确定”或“根据确定”或“在检测到……时”或“响应于检测到”所述条件先决条件为真。

尽管各种附图中的一些以特定顺序示出了多个逻辑阶段，但是不依赖于顺序的阶段可以被重新排序，并且其他阶段可以被组合或分解。尽管具体提到了一些重新排序或其他分组，但是其他分组对于本领域普通技术人员而言将是显而易见的，因此并未提供详尽的替代列表。此外，应当认识到，这些阶段可以以硬件、固件、软件或其任何组合来实现。

为了说明的目的，已经参考特定实施例描述了前述描述。然而，以上说明性讨论并非旨在穷举或将本发明限制为所公开的精确形式。鉴于以上教导，许多修改和变化是可能的。选择和描述实施例是为了最好地解释本发明的原理及其实际应用，从而使本领域的其他技术人员能够最好地利用本发明以及具有适于预期的特定用途的各种修改的各种实施例。

Claims (9)

1.一种油气勘探方法，包括：

a.在计算机处理器处接收表示感兴趣的地下体积的地震数据集和先验模型规范；

b.通过计算机处理器，使用基于树的贝叶斯采样方法执行所述地震数据集的全波形反演，所述基于树的贝叶斯采样方法使用维度减少的先验模型规范，自动适应模型复杂度，并生成适应所述地震数据集的空间分辨率的经采样的地球模型的集合；

c.通过计算机处理器，从地震数据和经采样的地球模型的集合生成地震图像的集合；

d.通过计算机处理器，使用经采样的地球模型的集合和地震数据执行地震图像的集合的不确定性分析；

e.使用地震图像的集合、经采样的地球模型的集合以及不确定性分析来评估感兴趣的地下体积，并确定钻井的位置以及从感兴趣的地下体积开采油气；

其中：

所述全波形反演使用关于树结构的简化模型基来进行模型压缩；

所述简化模型基和所述模型压缩使用离散小波变换；以及

通过找到树的每个级别上的最小和最大DWT系数来指定先验系数界限，并且将每个级别上的系数的界限设置成高于和低于最小值和最大值，以将先验模型限制为仅在给定范围的可行模型之内。

2.根据权利要求1所述的方法，其中，所述全波形反演使用可逆跳跃马尔可夫链蒙特卡罗(RJ-McMC)方法。

3.根据权利要求1所述的方法，其中，所述全波形反演使用并行退火来避开局部最小值。

4.根据权利要求1所述的方法，其中，所述方法自动选择速度模型的有限的一组DWT系数并使用所述有限的一组DWT系数进行操作。

5.根据权利要求1所述的方法，其中，所述经采样的地球模型的集合中的每个地球模型包括P波速度(V_p)、横波速度(V_s)和密度(ρ)中的一个或多个。

6.根据权利要求1所述的方法，其中：所述树结构包括其中第一级别对应于树的根节点的节点级别。

7.根据权利要求1或权利要求6所述的方法，其中：

每个级别上的界限是小于和大于所述级别的最小和最大DWT系数的2％。

8.一种计算机系统，包括：

一个或多个处理器；

存储器；以及

一个或多个程序，其中所述一个或多个程序被存储在存储器中并且被配置为由所述一个或多个处理器执行，所述一个或多个程序包括指令，所述指令在由所述一个或多个处理器执行时使计算机系统：

在所述一个或多个处理器处接收表示感兴趣的地下体积的地震数据集和先验模型规范；

通过所述一个或多个处理器，使用基于树的贝叶斯采样方法执行所述地震数据集的全波形反演，所述基于树的贝叶斯采样方法使用维度减少的先验模型规范，自动适应模型复杂度，并生成适应所述地震数据集的空间分辨率的经采样的地球模型的集合；

通过所述一个或多个处理器，从地震数据和经采样的地球模型的集合生成地震图像的集合；

通过所述一个或多个处理器，使用经采样的地球模型的集合和地震数据执行地震图像的集合的不确定性分析；以及

在用户界面上显示地震图像的集合、经采样的地球模型的集合以及不确定性分析中的至少一个，以允许评估感兴趣的地下体积，并确定钻井的位置以及从感兴趣的地下体积开采油气，

其中：

所述全波形反演使用关于树结构的简化模型基来进行模型压缩；

所述简化模型基和所述模型压缩使用离散小波变换；以及

通过找到树的每个级别上的最小和最大DWT系数来指定先验系数界限，并且将每个级别上的系数的界限设置成高于和低于最小值和最大值，以将先验模型限制为仅在给定范围的可行模型之内。

9.一种存储一个或多个程序的非暂时性计算机可读存储介质，所述一个或多个程序包括指令，所述指令当由具有一个或多个处理器和存储器的电子装置执行时，使装置：

在所述一个或多个处理器处接收表示感兴趣的地下体积的地震数据集和先验模型规范；

通过所述一个或多个处理器，使用基于树的贝叶斯采样方法执行所述地震数据集的全波形反演，所述基于树的贝叶斯采样方法使用维度减少的先验模型规范，自动适应模型复杂度，并生成适应所述地震数据集的空间分辨率的经采样的地球模型的集合；

通过所述一个或多个处理器，从地震数据和经采样的地球模型的集合生成地震图像的集合；

通过所述一个或多个处理器，使用经采样的地球模型的集合和地震数据执行地震图像的集合的不确定性分析；以及

在用户界面上显示地震图像的集合、经采样的地球模型的集合以及不确定性分析中的至少一个，以允许评估感兴趣的地下体积，并确定钻井的位置以及从感兴趣的地下体积开采油气，

其中：

所述全波形反演使用关于树结构的简化模型基来进行模型压缩；

所述简化模型基和所述模型压缩使用离散小波变换；以及

通过找到树的每个级别上的最小和最大DWT系数来指定先验系数界限，并且将每个级别上的系数的界限设置成高于和低于最小值和最大值，以将先验模型限制为仅在给定范围的可行模型之内。

Images