CN102326098A - 用于估计地球模型参数的地球物理数据随机反演 - Google Patents

Abstract

一种用于估计地球模型的模型参数的计算机实现随机反演方法。在一个实施例中，该方法利用基于采样的随机技术来确定定义地下结构的基于边界的多维模型的模型参数的概率密度函数(PDF)。在一些实施例中，利用了称为马尔可夫链蒙特卡罗(MCMC)的采样技术。MCMC技术属于“重要性采样”技术的类别，其中与模型拟合或匹配指定获取几何条件的能力成比例地采样后验概率分布。在另一个实施例中，该反演包括多个地球物理数据组的联合反演。本发明的实施例还涉及配置成执行用于估计模型参数以便精确解释地球地下结构的方法的计算机系统。

Description

用于估计地球模型参数的地球物理数据随机反演

技术领域

本发明涉及使用基于采样的随机方法来导出模型误差和模型参数的精确估计值的地球物理数据反演方法。

背景技术

从地球物理数据中估计用于油气勘探的模型参数是具有挑战性的，并且受高度不确定性影响。像地震振幅对角度(amplitude versusangle，AVA)和振幅对偏移(amplitude versus offset，AVO)反演那样的地震成像技术可以产生潜在储层岩石的物理位置和孔隙率的高度精确的估计值，但是在许多情况下，只有辨别储层内的流体的有限能力。像电磁(EM)方法那样的其它地球物理数据可以添加有关水饱和度，乃至烃饱和度的信息，因为岩石的电导率对水饱和度非常敏感。然而，单独使用EM数据来估计流体饱和度是不可行的，因为EM数据的空间分辨率很低。地震和EM方法对储层物质的不同物理性质敏感：地震数据是地震P-和S-波速度以及储层密度的函数，并且EM数据是储层的电阻率的函数。因为岩石的弹性和电性质这两者通过岩石-物理模型与流体饱和度和孔隙率建立起物理联系，所以像地震数据和EM数据那样的多个地球物理数据组的联合反演有可能比各个数据组的反演更好地估计像流体饱和度和孔隙率那样的地球模型参数。

用于导出模型误差和模型参数的估计值的现有技术地球物理数据反演通常依赖于基于梯度的技术，该技术使包含数据失配项并可能包含附加模型正则化或平滑项的目标函数最小化。例如，方程(1)是常用的一般目标函数φ：

φ(m，d)＝[D(d^o-d^p)]^H[(D(d^obs-d^p))]+λ(Wm)^H(Wm)    (1)

其中，D是数据协方差矩阵，d^o和d^p分别是观察到的和预测的数据，W是模型正则化矩阵，m是模型参数的矢量，它可以是电导率，以及λ是相对于数据失配衡量模型平滑的重要性的权衡参数。由于数据d是复数，所以H表示转置共轭算子。线性化有关第i次迭代时的给定模型m_i的方程(1)得出二次形式：

(J^TS^TSJ+λW^TW)m_i+1＝J^TS^TSJm_i+J^TS^TSδd_i    (2)

其中，m_i+1可以使用许多技术来求解，二次规划算法是一种可能。J是数据对模型参数的偏导数的雅可比(Jacobian)矩阵，S是包含数据的标准偏差的倒数，使得S^T＝D^-1的矩阵。计算出的数据(d^p)与观察到的数据(d^obs)之间的当前差通过δd_i＝d^obs-d_i给出。随着迭代继续进行下去，从大到小调整权衡参数λ。

当方程(1)和(2)所描述的算法收敛到目标函数φ的最小值时，产生单个模型m。这种现有技术导出的模型不能以任何方式保证是“全局的”或真实的模型。模型参数误差，也称为模型参数标准偏差(方差的平方根)，从像方程(1)和(2)所描述那样的模型参数协方差计算中导出的估计值是不精确的，不能提供真实模型参数误差的充分量化。

与现有技术的反演方法不同，使用基于采样的随机模型的地球物理数据组反演可以提供所有模型参数值的概率密度函数(PDF)的精确估计值。并且，基于采样的随机方法可以用于像地震和EM数据那样的多个地球物理数据组的联合反演，以便比各个数据组的反演更好地估计地球模型参数。术语随机反演广泛用于涵盖用于确定模型参数变量的PDF的许多不同手段。模型参数PDF提供了每个模型参数的方差的精确估计值以及各个模型参数的均值、最频值和中值。当比较多个模型以便为精确解释地球的地下结构(subsurface)确定最可能模型时可以使用精确的模型参数方差。

发明内容

本发明的一个方面涉及一种用于估计地球模型的模型参数的随机反演方法，其具有如下操作：获取对感兴趣地下地质体的一部分进行采样的至少一个地球物理数据组，每个地球物理数据组定义感兴趣地下地质体的获取几何条件(acquisition geometry)；生成感兴趣地下地质体的指定数目的基于边界的多维模型，所述模型通过模型参数来定义；生成每种指定获取几何条件的模型的正演模型响应；生成正演模型响应与每种指定获取几何条件的地球物理数据组匹配的似然值；将模型参数保存成每个模型的马尔可夫链(Markov Chain)的一个元素；对马尔可夫链的收敛性进行测试；更新每个模型的模型参数的值并且顺序地或并行地重复上面的操作，直到达到收敛；导出形成收敛马尔可夫链的模型的每个模型参数的概率密度函数；以及从每个模型的每个模型参数的概率密度函数中计算方差、均值、最频值和中值，以便生成感兴趣地下地质体的地球模型的模型参数方差和模型参数的估计值，所述估计值用于确定感兴趣地下地质体的特性。

本发明的另一个方面涉及一种被配置成生成感兴趣地质体的多维模型的系统。在一个实施例中，该系统被配置成运行包含程序的计算机可读媒体，所述程序当被执行时完成包含如下步骤的操作：获取对感兴趣地下地质体的一部分进行采样的至少一个地球物理数据组，每个地球物理数据组定义感兴趣地下地质体的获取几何条件；生成感兴趣地下地质体的指定数目的基于边界的多维模型，所述模型通过模型参数来定义；生成每种指定获取几何条件的模型的正演模型响应；生成正演模型响应与每种指定获取几何条件的地球物理数据组匹配的似然值；将模型参数保存成每个模型的马尔可夫链的一个元素；对马尔可夫链的收敛性进行测试；更新每个模型的模型参数的值并且顺序地或并行地重复上面的操作，直到达到收敛；导出形成收敛马尔可夫链的模型的每个模型参数的概率密度函数；从每个模型的每个模型参数的概率密度函数中计算方差、均值、最频值和中值，以便生成感兴趣地下地质体的地球模型的模型参数方差和模型参数的估计值，所述估计值用于确定感兴趣地下地质体的特性。

通过参照形成本说明书的一部分、相同标号在各个图形中表示相应部件的附图考虑如下描述和所附权利要求书，本发明的这些和其它目的、特征、和特性，以及操作方法、相关结构元件的功能、部件的组合、和制造经济性将变得更加显而易见。但是，不言而喻，这些附图只用于例示和描述的目的，而无意作为本发明的限制的定义。正如用在说明书和权利要求书中的那样，除非上下文另有明确指明，单数形式“一个”、“一种”、和“该”也包括复数指示物。

附图说明

图1例示了根据本发明一个或多个实施例的估计模型参数的方法的流程图。

图2A和2B例示了根据本发明一个或多个实施例的用于模型参数化的方法。

图3A至3E例示了按照本发明一个或多个实施例的电导率模型以及所计算模型参数方差、均值、最频值和中值。

图4例示了按照本发明一个或多个实施例的地下结构的基于边界的多维模型。

图5例示了根据本发明一个或多个实施例的用于执行随机反演方法的系统。

具体实施方式

本发明的实施例提供了一种用于估计地球模型的模型参数的计算机实现的随机反演方法。在一个实施例中，该方法利用基于采样的随机技术来确定定义地下结构的基于边界的多维模型的模型参数的概率密度函数(PDF)。在一些实施例中，利用了称为马尔可夫链蒙特卡罗(MCMC)的采样技术。MCMC技术属于“重要性采样”技术的类别，其中与模型拟合或匹配指定数据组的能力成比例地采样后验概率分布。重要性采样导致模型空间的不均匀采样，该模型空间的不均匀采样与更传统的采样技术相比以减少的正演函数调用次数表征高概率区域。在另一个实施例中，该反演包括多个地球物理数据组的联合反演。本发明的实施例还涉及一种被配置成执行用于估计模型参数以便精确解释地球地下结构的方法的计算机系统。

现在参照图1，这个图形示出了用于估计地球模型的模型参数的方法10。下面给出的方法10的操作是示例性的。在一些实施例中，方法10可以借助未描述的一个或多个附加操作完成，和/或可以在没有所讨论的一个或多个操作的情况下完成。另外，在图1中例示并且在下面描述的方法10的操作的次序并非是限制性的。

方法10从获取至少一个地球物理数据组的操作12开始。每个地球物理数据组采样感兴趣地下地质体的某个部分，并用于定义感兴趣地质体的地球物理获取几何条件。该地球物理数据组可以包括，例如，受控源电磁数据(CSEM)、大地电磁数据、重力数据、磁性数据、地震数据、井生产数据或上述数据的任何组合。

该地球物理获取几何条件是源和接收机的空间位置(X，Y，Z)，以及像源拖曳速度那样的任何系统操作参数和将源波形完整描述为时间函数的参数。

在操作14中，生成感兴趣地下地质体的指定数目的基于边界的多维模型。在一个实施例中，每个模型包括在感兴趣地质体中的称作节点的多个位置上定义的模型参数。指定数目的基于边界的多维模型是这样生成的：与用户定义的模型参数方差一起将初始模型参数取作平均值，以便从用户定义的分布(通常是高斯分布，但也可以使用其它分布)中随机选择新的一组模型参数值。每组唯一模型参数定义感兴趣地质体的二维(2D)，3D或4D模型。

在另一个实施例中，每个模型参数具有包含用于定义节点空间位置的X，Y和Z值和至少一个地球物理性质值的值。在一个实施例中，该地球物理性质值可以包括电导率或电阻率(σ)、压缩速度(V_p)、剪切速度(V_s)和密度(ρ)。另外，该地球物理性质值可以包括像流体饱和度、压力、温度和/或孔隙率那样的储层参数。这些储层参数可以通过包含将地球物理参数与储层参数联系在一起的岩石物理模型(从测井数据中导出)而被包括进来。举非限制性的例子来说，该模型可以用如上所述的地球物理性质和/或像流体饱和度和孔隙率那样的储层性质来参数化。在一个实施例中，在唯一一组节点之间的空间中内插地球物理性质值，以定义模型边界和生成计算正演模型响应所需的模型参数。然后，可以将基于边界的多维模型投影在有限差或有限元网格上，以便计算正演模型响应。投影过程是这样完成的：将从边界上面和下面的性质的空间加权平均中导出的性质值指定给跨过该边界的单元(cell)。

根据本发明一个或多个实施例的模型参数化的概括表示显示在图2中。图2A示出了基于标准单元的模型参数化。每个单元32上的地球物理性质是对所有相邻单元34应用正则化或平滑的模型参数。该单元模型通常用在有限差或有限元计算中，以便生成模型响应。图2B示出了基于边界的模型参数化。控制边界38的节点36可以是模型参数，每个节点具有一个位置和一个或多个地球物理性质值。可以在节点之间横向地和垂直地内插模型参数。在一个实施例中，可以将所得到的基于边界的模型投影在有限差或有限元网格上，以便计算正演模型响应。

显示在图2B中的基于边界的模型参数化解决了传统的基于单元的模型参数化带来的几个固有问题。举例来说，典型的2D有限差网格需要1万到10万个数量级的单元来精确计算对感兴趣模型的正演模型响应。如果在反演中将每个单元的地球物理性质用作模型参数，则对于随机反演，会引起两个主要问题。首先，如图2A所示，相邻单元的电导率高度相关，使收敛变慢。其次，基于单元的模型参数化不能良好地对应于具有其固有层状结构、边界在相对均匀性质的层之间的沉积地质结构。基于边界的模型参数化极大地减少了模型参数的数量，只有节点上的模型参数才用在反演中，同时，模型的边界对应于地质边界。

返回参照图1，在操作16中，生成每种指定获取几何条件的模型的正演模型响应。包括边界上的节点上的所有模型参数的模型通过将像电导率那样的地球物理性质投影在有限差或有限元网格上并计算正演模型响应而被正演模拟。在操作18中，使用数值计算来生成正演模型响应与每种指定获取几何条件的地球物理数据组匹配的似然值。将反演问题定义成需要定义模型似然函数以及模型参数的先验分布的贝叶斯(Bayesian)推理问题。在其最简单形式中，这得到：

f(p|d)∝f(d|p)f(p)    (3)

其中，f(p|d)是给定数据d，所有未知模型参数p的联合后验PDF。该方程右侧第一项f(d|p)是给定模型参数p，数据d的似然函数，以及右侧最后一项f(p)是所有未知模型参数p的先验分布。似然函数f(d|p)可以取决于用于描述数据中的噪声的模型而采用各种形式。举非限制性的例子来说，对于电磁数据E，当将误差指定成测量值的百分比时，它采用如下形式：

$f\left(E\right|\mathrm{\σ})=\underset{i=1}{\overset{\mathrm{nf}}{\mathrm{\Π}}}\underset{j=1}{\overset{\mathrm{no}}{\mathrm{\Π}}}\underset{k=1}{\overset{2}{\mathrm{\Π}}}\frac{1}{\sqrt{2{\mathrm{\π\β}}_j^2}}\exp \{-\frac{1}{{2\mathrm{\β}}_j^2}{\left(\frac{e_{\mathrm{ijk}}-M_{\mathrm{ijk}}^e\left(\mathrm{\σ}\right)}{e_{\mathrm{ijk}}}\right)}^2\}---\left(4\right)$

其中，模型电导率是σ。观察到的电场e是频率(频率总数等于nf)、偏移(偏移总数等于no)的函数，并且是导致2个分量(实分量和正交分量)的复数。正演模型响应是

误差β是偏移(j)的函数，并且是e的百分比。

如果在存在N_t个时间样本和N_a个反射角的发射域中模拟反射地震数据，则可以将正演模拟地震数据表示成：

$r_{\mathrm{ij}}=S_{\mathrm{ij}}^a(V_p,V_s,\mathrm{\ρ})+{\mathrm{\ϵ}}_{\mathrm{ij}}^a---\left(5\right)$

其中，i是时间索引(1到N_t)并且j是角度索引(1到N_a)，S是地震正演问题，以及ε是测量误差。如果假设误差是高斯型，那么，地震似然函数f将是：

其中，ε是数据误差的矢量，而∑是数据协方差矩阵。

在一些实施例中，可以使用MCMC采样技术来估计f(d|p)。由于模型参数的数量很大并且要求解的问题的维数很高，传统蒙特卡罗方法是不可行的。结合随机反演方程(4)和(6)的MCMC手段，例示在图2B中的基于边界的模型参数化也可以用于减少方程中的模型参数的数量。

在一些实施例中，可以使用多于一种数据类型以及它的相关联获取几何条件。例如，一个非限制性例子是CSEM和地震反射数据的组合。在这种情况下，每个边界节点上的地球物理参数是电导率、压缩速度、剪切速度和密度。上面方程(5)和(6)中的S所表示的地震正演模型可以是弹性波方程的有限差或有限元实现。当假设误差是高斯型时，地震数据的相关联似然函数将具有上面方程(6)的形式。

如果添加地震数据，那么，必须将贝叶斯反演问题，即方程(3)修改成：

f_SE(p|d)∝f(S|V_p，V_s，ρ)f(E|σ)f(V_p，V_s，ρ)f(σ)，(7)

其中，f_SE是组合PDF，f(S|V_p，V_s，ρ)是地震似然函数(方程6)，

f(E|σ)是CSEM似然函数(方程4)，f(V_p，V_s，ρ)是地震参数的先验分布，以及f(σ)是电导率参数的先验分布。

在图1中例示的操作20中，将模型参数保存成每个模型的马尔可夫链的一个元素。马尔可夫链的数量等于反演开始时指定模型的数量。本领域技术人员应该懂得，存在许多MCMC采样算法，例如，可以使用Metropolis-Hastings(Hasting，1970)采样和切片采样(Neal，2003)的组合来生成形成马尔可夫链的模型参数的序列。

在操作22中确定马尔可夫链是否在规定容限内达到收敛。当在给定噪声水平下，模型参数的PDF是真实分布的精确表示，使得分布的统计动差(均值、中值、最频值和方差)精确时，马尔可夫链是收敛的。可以多于一次，要么顺序地要么并行地运行马尔可夫链，以便确定收敛性。像Gelman和Rubin(1992)开发的方法和Raftery和Lewis(1992)开发的方法那样，可以使用许多方法来确定马尔可夫链的收敛性。在一个实施例中，使用Gelman和Rubin的方法，其中收敛性测试按如下计算序列内方差W和序列间方差B/n：

$B/n=\frac{1}{m-1}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{(p_j-\mathrm{\μ})}^2---\left(8\right)$

$W=\frac{1}{m(n-1)}{\mathrm{\Σ}}_{j=1}^m{\mathrm{\Σ}}_{t=1}^n{(p_{\mathrm{jt}}-{\mathrm{\μ}}_j)}^2---\left(9\right)$

其中，m是马尔可夫链的数量，n是链中样本的数量，以及p_jt是链j中的p的n次迭代的第t次迭代。

在计算了(8)和(9)之后，像如下那样通过B和W的加权平均计算模型参数方差σ²的估计值：

${\mathrm{\σ}}_{+}^2=\frac{n-1}{n}W+\frac{B}{n}---\left(10\right)$

像如下那样计算可能缩小因子：

$R=\frac{m+1}{m}\frac{{\mathrm{\σ}}_{+}^2}{W}-\frac{n-1}{\mathrm{mn}}---\left(11\right)$

如果方程(11)中的R接近1，则认为马尔可夫链已经收敛。关于如何接近才算“足够接近”的确定通过在真实数据的反演之前对合成模型进行测试来完成。例如，在本发明的一些实施例中，对于R来说，值1.1可能就足够了。

在操作24中，更新每个模型的模型参数的值。如果在操作22中未达到收敛(R未接近1)，则通过任意数量的采样算法更新模型参数。一个非限制性例子是使用Metropolis-Hastings采样算法，它将现有参数值用作平均值以及使用假设方差的高斯分布来抽取每个模型参数的新值。另一种可能是使用切片采样算法生成新模型参数。在一个实施例中，每次迭代时采样算法的选择通过随机抽取反演开始时指定用在任何迭代上的每种采样技术的概率的统一变量来确定。例如，为多方差Metropolis-Hasting采样、多方差切片采样、单方差Metropolis-Hasting采样、和单方差切片采样分别指定0.4、0.4、0.1、0.1的概率，然后，在每次迭代时，生成统一随机数，它的值决定使用哪一种采样技术。一旦生成新模型，工作流10就转到操作26，重复操作14到22，直到在规定容限内达到收敛，即，当模型参数的PDF是模型参数的真实分布的精确表示时。如果迭代不充分，或如果马尔可夫链的数量不足，则马尔可夫链可能不收敛。当方程(11)中的R在值1的规定容限内时，马尔可夫链已经收敛，可以停止采样。该算法针对串行(单个处理器)和并行(集群)计算两者来实现。串行实现用于简单模型测试和算法细化，而并行实现最适用于大规模生产反演。

如果在操作22中达到收敛，该方法转到操作28，其中将来自马尔可夫链的所有模型参数值分箱(bin)，以便产生每个参数的PDF。PDF给出模型参数与地球物理数据组一致的概率。PDF包含模型参数的最有可能的值，并且量化估计值的不确定性。

在操作30中，使用从操作28所得到的每个模型参数的PDF来计算感兴趣地下地质体的每个模型的每个模型参数的方差、均值、最频值和中值。各自对储层物质的不同物理性质敏感的多种地球物理数据类型的联合反演比各个数据组的反演更好地估计地球模型参数。

在一个实施例中，将所计算的每个模型参数PDF的均值、最频值和中值用作模型参数值，并且进行正演模拟，以计算数据均方根(RMS)数据失配。这是计算的数据与地球物理数据组匹配得怎么样的量度。将来自模型参数PDF的均值、最频值和中值的3种模型以及相关联的RMS数据失配提供给用户。另外，可以与PDF的图形例示以及均值、最频值和中值一起计算每个模型参数的方差。提供可以用于解释感兴趣地下地质体的模型参数的更好估计值。

在另一个实施例中，用户可以将任何或所有这三种模型用作开始模型用于确定性反演(最小二乘)，即方程(1)的最小化，以找出具有最小RMS数据失配的模型。这种操作向确定性反演提供了极大地提高其达到全局最小的机会的改进开始模型。

应该懂得，图1的工作流旨在包含用于估计模型参数的几种情形。在一些实施例中，使用与随机反演耦合的基于边界的模型参数化的方法10被扩展到包括多种数据类型的联合反演，以估计单个自洽模型。在其它实施例中，可以用地球物理和/或储层性质将联合反演模型参数化。该手段可以在2D，3D和4D中实现。在一个非限制性例子中，使用基于边界的模型参数化的随机MCMC采样实现显示在方程(4)中的基于有限元的反演算法。模拟地球物理性质是电导率(σ)。在2D中，使用CSEM拖缆(水平轴)上的节点之间的线性内插函数和垂直轴上的深度来描述边界。在3D中，通过以X，Y和Z表示的节点来描述边界，其中将双线性或更高阶内插函数用于生成表面Z和3D模型内的任意位置上的模型电导率值。对于各向同性、横向各向异性或完全各向异性模型，电导率可以分别是标量，矢量或张量。对于各向同性情况，σ的值与坐标方向无关。对于横向各向异性，存在3个电导率分量(σ_x，σ_y，σ_z)，每一个分量都可以独立估计。最一般情况，即完全各向异性，导致通过下式给出的对称电导率张量：

$\left[\begin{array}{ccc}{\mathrm{\σ}}_{\mathrm{xx}}& {\mathrm{\σ}}_{\mathrm{xy}}& {\mathrm{\σ}}_{\mathrm{xz}}\\ {\mathrm{\σ}}_{\mathrm{yx}}& {\mathrm{\σ}}_{\mathrm{yy}}& {\mathrm{\σ}}_{\mathrm{yz}}\\ {\mathrm{\σ}}_{\mathrm{zx}}& {\mathrm{\σ}}_{\mathrm{zy}}& {\mathrm{\σ}}_{\mathrm{zz}}\end{array}\right]---\left(9\right)$

其中，非对角项是对称的(σ_xy＝σ_yx等)。因此，各向同性、横向各向异性和完全各向异性结果的反演在反演时每个表面每个节点分别需要估计1，3或6个电导率参数。在使用其它地球物理数据(例如地震数据)的其它实施例中，将向针对电导率所说明的那样处理相关联的地球物理参数。

传统上，将生产储层上的连续时间的3D地震数据用于储层中的图像变化称为4D。当从石油储层中生产时，通过反演在生产期间顺序时间获取的观察到的数据，可以将新的反演技术用作监视地球中随时间变化的4D技术。反演每个时步的观察到的数据，并且使用时间之间的变化确定储层在何处变化。

图3例示了按照本发明一个或多个实施例的电导率模型。在图3A中示出了P1-P4的简单分层1D电导率模型。对CSEM勘查获取几何条件的模拟响应显示在图3B-3E中。图3B对应于P1，图3C对应于P2，图3D对应于P3，图3E对应于P4。电磁源位于海底上方50m处，而电磁接收机被放置在海底上。源是工作在0.25，0.75和1.25Hz频率上的电偶极子。接收机偏离源500到1500m。模拟数据具有反演之前加入的10％高斯随机噪声。

虽然方程(1)和(2)给出的基于梯度的算法可以求解，在本例中，与PDF的最频值(PDF Mode)一样接近真实值(True Value)，所得到的模型参数标准误差(Cov SD)在2到5的倍数之间。如果将所得到的模型误差估计值用在任何定量风险评估过程中，则这是传统反演手段的严重缺点。在方程(1)极小时从模型协方差矩阵中计算模型参数误差与模型参数PDF标准误差(PDF SD)相比是不充分的。

图4例示了按照本发明一个或多个实施例的随机反演的方法。图形例示40示出了具有10个层的2D基于边界的模型。控制边界42深度(Z)的节点是黑色实心菱形44，控制地球物理性质的节点是空心方形46。白线是从地震数据中解释的层位(horizon)48。灰度是电导率，明亮是导电性的，灰暗是电阻性的。为2个节点显示了示范性电阻率PDF。所示的模型将每个PDF的最频值用作节点电导率。

在一些实施例中，方法10可以在一个或多个处理设备(例如，数字处理器、模拟处理器、设计成处理信息的数字电路、设计成处理信息的模拟电路、状态机、和/或用电子方式处理信息的其它机制)中实现。该一个或多个处理设备可以包括响应以电子方式存储在电子存储媒体上的指令，执行方法10的一些或所有操作的一个或多个设备。该一个或多个处理设备可以包括通过硬件、固件、和/或软件配置成为执行方法10的一个或多个操作而专门设计的一个或多个设备。配置成运行包含当被执行时完成方法10的操作的程序的计算机可读媒体的系统示意性地例示在图5中。该系统50包括数据存储设备或存储器52。可以使所存储数据可用于像可编程通用计算机那样的处理器54。处理器54可以包括像显示器56和图形用户界面58那样的接口组件。图形用户界面(GUI)可以用于显示数据和经处理数据产物，和允许用户在实现该方法的各个方面的选项当中作出选择。数据可以经由总线60直接从数据获取设备，或从中间存储设备或处理设施(未示出)传送给系统50。

尽管为了例示的目的，根据当前认为最可行的优选实施例对本发明作了详细描述，但应该明白，这样的细节仅仅为了该目的，本发明不局限于所公开的实施例，而是相反，打算涵盖在所附权利要求书的精神和范围之内的所有修改和等同安排。例如，应该理解，本发明设想，可以尽可能地将任何实施例的一个或多个特征与任何其它实施例的一个或多个特征组合。

Claims (19)

1.一种计算机实现的随机反演方法，用于估计感兴趣地下地质体的地球模型的模型参数，所述方法包含：

a)获取对感兴趣地下地质体的一部分进行采样的至少一个地球物理数据组，每个地球物理数据组定义感兴趣地下地质体的获取几何条件；

b)生成感兴趣地下地质体的指定数目的基于边界的多维模型，所述模型通过模型参数来定义；

c)生成每种指定获取几何条件的模型的正演模型响应；

d)生成正演模型响应与每种指定获取几何条件的地球物理数据组匹配的似然值；

e)将模型参数保存成每个模型的马尔可夫链的一个元素；

f)对马尔可夫链的收敛性进行测试；

g)更新每个模型的模型参数的值并且顺序地或并行地重复b)到f)，直到达到收敛；

h)导出形成收敛马尔可夫链的模型的每个模型参数的概率密度函数；以及

i)从每个模型的每个模型参数的概率密度函数中计算方差、均值、最频值和中值，以便生成感兴趣地下地质体的地球模型的模型参数方差和模型参数的估计值，所述估计值用于确定感兴趣地下地质体的特性。

2.如权利要求1所述的方法，其中，所述多维模型是2D、3D或4D模型。

3.如权利要求1所述的方法，其中，所述模型参数被定义在节点上，并且具有包含用于定义节点空间位置的X，Y和Z值和至少一个地球物理性质值的值。

4.如权利要求3所述的方法，其中，所述模型边界通过一个唯一节点集合来定义。

5.如权利要求3所述的方法，其中，所述至少一个地球物理性质值包括电导率或电阻率、压缩速度、剪切速度、密度、流体饱和度、压力、温度和/或孔隙率。

6.如权利要求3所述的方法，其中，所述节点通过在节点之间的空间中内插地球物理性质值来定义边界，以便生成计算正演模型响应所需的模型参数。

7.如权利要求1所述的方法，其中，将所述模型投影在有限差或有限元网格上，以便计算正演模型响应。

8.如权利要求1所述的方法，其中，所述地球物理数据组包括受控源电磁数据、大地电磁数据、重力数据、磁性数据、地震数据、井生产数据或上述数据的任何组合。

9.如权利要求1所述的方法，其中，计算多于一种类型的地球物理数据组的正演模型响应，以便得出组合地球物理数据组的似然值，从而产生联合反演。

10.如权利要求1所述的方法，其中，使用马尔可夫链蒙特卡罗(MCMC)采样来生成形成马尔可夫链的模型参数的序列。

11.如权利要求10所述的方法，其中，MCMC采样算法包括Metropolis-Hastings采样和切片采样算法。

12.如权利要求1所述的方法，其中，将采样算法应用于为每次模型运算更新模型参数的值的每次迭代，并且通过随机抽取统一变量来确定所述采样算法，其中反演开始时指定用在任何迭代上的每种采样技术的概率。

13.如权利要求1所述的方法，其中，计算正演模型响应与指定获取几何条件的每个地球物理数据组匹配的似然值根据似然函数来进行。

14.如权利要求1所述的方法，其中，确定马尔可夫链的收敛性根据下述可能缩小因子来进行：

$R=\frac{m+1}{m}\frac{{\mathrm{\σ}}_{+}^2}{W}-\frac{n-1}{\mathrm{mn}}$

15.如权利要求1所述的方法，其中，确定马尔可夫链的收敛性计算序列内方差和序列间方差。

16.如权利要求1所述的方法，其中，使用通过模型参数均值、最频值和/或中值定义的模型来定义用于确定性(例如，最小二乘)反演的开始模型，以便生成具有最小均方根数据失配的模型。

17.如权利要求1所述的方法，其中，将通过模型参数均值、最频值和/或中值、和/或模型参数方差定义的模型用图形显示出来。

18.一种被配置成运行包含程序的计算机可读媒体的系统，所述程序当被执行时完成包含如下步骤的操作：

a)获取对感兴趣地下地质体的一部分进行采样的至少一个地球物理数据组，每个地球物理数据组定义感兴趣地下地质体的获取几何条件；

b)生成感兴趣地下地质体的指定数目的基于边界的多维模型，所述模型通过模型参数来定义；

c)生成每种指定获取几何条件的模型的正演模型响应；

d)生成正演模型响应与每种指定获取几何条件的地球物理数据组匹配的似然值；

e)将模型参数保存成每个模型的马尔可夫链的一个元素；

f)对马尔可夫链的收敛性进行测试；

g)更新每个模型的模型参数的值并且顺序地或并行地重复b)到f)，直到达到收敛；

h)导出形成收敛马尔可夫链的模型的每个模型参数的概率密度函数；以及

i)从每个模型的每个模型参数的概率密度函数中计算方差、均值、最频值和中值，以便生成感兴趣地下地质体的地球模型的模型参数方差和模型参数的估计值，所述估计值用于确定感兴趣地下地质体的特性。

19.如权利要求18所述的方法，其中，将通过模型参数均值、最频值和/或中值、和/或模型参数方差定义的模型用图形显示出来。

Images