CN101218517B - 预测非唯一解集合中的最佳和最坏解的方法 - Google Patents

Abstract

用于为地震波反射振幅的物理数据的矩阵求逆生成的孔隙度和页岩体积分数的模型参数的值决定最佳和最坏情况的方法(图1)。矩阵被对角化，且与不重要的对角元素相关的正交单位基向量被用于生成解的上下边界(6)。最佳和最坏情况解被确定为零基向量的线性组合，其展开系数通过对上下边界(J)作最佳拟合来确定。

Description

预测非唯一解集合中的最佳和最坏解的方法

技术领域

本申请要求2005年7月13号提交的美国临时申请第60/698,760号的权益。

本发明一般地涉及地球物理建模领域，尽管本发明具有更广泛的应用。具体地，本发明是一种方法，其用于当模型反演产生非唯一解时，预测最佳和最坏解。

背景技术

在石油工业中，常见的是面对一套数据，人们希望从其中推断出某种感兴趣信息。也相当常见的是，这种反演问题是非唯一的，即不同的解同等好地解释数据。虽然获取用户认为最有可能的单一解是简单的，但经常希望知道除最有可能的解之外，拟合数据的最佳和最坏情况的解，以充分理解所给行动过程的风险。在石油工业中，这种问题的示例是预测油藏中沙子和孔隙的分布，其中人们想知道除最有可能的碳氢化合物体积之外的可能的最大和最小碳氢化合物体积(即最佳和最坏情况模式)。准确理解给潜在油藏排水包含的潜在风险将减少总成本(合适大小的平台，最优排水策略等)。

用于确定可替换模式的普通方法是进行多种不同模型的预期或前向模拟，在这些模型中，被认为影响最后结果的变量从某些预定义分布中随机地选择。前向模型然后与所观察的数据相比较，以查看各种前向模型中的哪些匹配。将拟合数据的参数的分布从模型集中提取出来，这些模型被认为很好地拟合了数据。从这个分布中，最佳和最坏情况原则上可能被确定。这个方法很耗时间，因为其需要大量的前向模型。此外，因为仅有的所试模型是用户已经考虑或认为相关的那些模型，所以它因用户基础而蒙受损害。

另一种方法是取最有可能的模型并只是将它按比例放大和缩小一些量，并称这为最佳和最坏情况。本方法导致的结果通常在预期建模意义上并不与所观察的数据相匹配，且不必然是最佳和最坏情况答案。

所需要的是一种方法，其中最佳和最坏情况模式作为反演问题的数学解来获取。

发明内容

在一个实施例中，本发明是计算机实现的方法，用于(看图1流程图的框1)确定矩阵方程式的最大和最小可能解，该矩阵方程式可以表达成以下形式：

$\left[G\right]\left[\begin{array}{c}{m}_1\\ m_2\\ \·\\ \·\\ \·\\ m_N\end{array}\right]=\left[\mathrm{data}\right]$

其中，m₁...m_N是所要求解的物理参数，G是基于将m_i与已测数据相联系的物理系统模型的矩阵，其中该方程式可以通过产生无限数量的可能解的数值方法进行非唯一逆运算，所有这些可能解都基本同样好地拟合所述数据，并从这些可能解中确定最可能的解，所述方法包含(a)寻找(图1中的步骤2)使矩阵G对角化的正交单位基向量，并使用所述向量将G对角化；(b)选择(步骤3)阈值，在对所述矩阵方程式的最可能解的影响方面，低于阈值的被对角化的G中元素的值被认为不重要；(c)从正交单位基向量中识别(步骤4)与不重要的对角元素相关的那些向量(“零”向量)；(d)选择(步骤5)L_P数学范数，其中p∈[0，∞]；(e)通过根据所选范数求取所述零向量之对应分量的总和来决定(步骤6)可能解m₁，m₂...m_N的上界和下界，所述下界由所述总和取负值给出；(f)查找(步骤7)最接近逼近所述上界的零向量的线性组合，并为所述下界重复这种查找；和(g)依次将所述两个线性组合中的每一个加到(步骤8)最有可能的解中，产生最大和最小解，按需要进行按比例缩放以去掉任何物理上不真实的结果。

附图说明

本发明和其优势将通过涉及如下详细描述和附图更好地得到理解，其中：

图1是流程图，示出本发明方法的基本步骤；和

图2A-C是由本发明方法得到的3-D图像的照片，示出倒转的页岩体积分数(vshale)体积中沙体的“最坏”情况(2A)和“最佳”情况(2C)模式，图2B示出最有可能的情况。

本发明将以其有关的优选实施例进行描述。然而，在一定程度上，以下描述针对特殊实施例或本发明的特殊使用，其目的仅在于示意说明，并不被解释为限于本发明的范围。相反，旨在覆盖所有可替换方案、修改方案和等价方案，这些方案都可以包含在本发明的精神和范围内，正如所附加的权利要求所定义的。

具体实施方式

本发明是一种方法，其通过求解将观察(数据)与感兴趣参数相联系的方程式系统，获得“最佳”和“最坏”情况的解。对于以上提到的岩性预测问题，这些方程式可以包括例如卷积方程式、Aki&Richards(1980)反射方程式和线性化的岩石物理方程式(见2005年6月24日由Saltzer、Finn和Lu提交的美国专利申请)。在矩阵符号中，它们采用以下形式：

$\left[G\right]\left[\begin{array}{c}\mathrm{\φ}\\ \mathrm{vsh}\end{array}\right]=\left[\mathrm{data}\right]---\left(1\right)$

其中，φ和vsh是作为时间的函数的孔隙度和vshale(shale volume fraction，页岩体积分数)值，data是与不同源接收器孔相关的地震轨迹，G是将模型参数(在此示例应用中是指页岩体积分数vshale和孔隙度)和数据参数相联系的矩阵，典型地，数据参数指地震波反射数据。然而，本发明可以应用于任何物理系统，其物理模型存在以提供将模型参数m与所测或者另外获取的数据相联系的G矩阵，其中在一般情况下方程式(1)可以写为Gm＝data。矩阵G可以分成两块：第一区域块的特征在于数据对模型参数的敏感性，而第二区域块几乎不具有敏感性。这两块区域通过定义将矩阵G对角化的标准正交基或正交单位基建立。一旦这些基建立，分离点值(cut-off value)被选取，小于分离点值的对角化矩阵G的元素不重要。

与对角化矩阵G中这些不重要的分量相关的正交单位向量是“零”向量。因此，如果基向量是u_k和v_k，它们就可以用于构造矩阵U和V以满足G＝USV′，其中S在其对角线上只具有非零元素。熟悉线性代数的人都知道，可以称作对角化G矩阵的矩阵S能够被确立。典型地，G进而S将不是方矩阵，但当i＝j时，元素G_i，j和元素S_i，j被认为是对角元素。如果S的对角元素称为λ₁，λ₂，…λ_N，且如果λ_k低于为重要性选择的阈值，那么v_k是零向量。数学上，零向量对应于Gm＝0       (2)的解(如上所述，m是列矩阵或向量，其分量是在方程式(1)所表示的实施例中的

值和vsh值。)结果，它们可以被加到最可能的解上而不改变该模型对所测数据的拟合，因为它们并不投影到数据空间。本发明的基础理论是，由于解的非唯一性，几乎同等好地拟合数据的无限数量的解可以被当作最可能解的扰动，且该扰动受到可以被构建的零向量的不同可能的线性组合的驱动。这根据这样的事实得出：零向量是生成相关空间的一部分的向量的基集合。因此，事实上并不影响最可能解的方程式部分导致了任何给定解和最可能解之间的差异。

那么问题演变成寻找零向量组合的问题，其将产生最大的孔隙度和沙子体积

这可以在数学上表达成

${\stackrel{\→}{m}}_{\mathrm{biggest}}=\underset{k}{\mathrm{\Σ}}{\mathrm{\α}}_k{\stackrel{\→}{v}}_k^{\mathrm{null}}---\left(3\right)$

其中，(α_k)是系数向量，用以将每个零向量的相对重要性加权。如果地表下的最佳可能模型

是先验地已知的，获得合适的权重因子(α_k)的最小平方解就是一件简单的事情。最困难的部分是确定孔隙度和vshale扰动的上界(即，已知求解

的什么模型)。

可能的扰动

的上界可以通过求与被确定为不重要的对角化G矩阵的部分(零向量v_k)相关的正交单位基向量的各元素或分量之和确定。该求和运算根据所选择应用的特别范数来完成。L₁范数用相关正交单位向量各元素的绝对值计算出：

${\stackrel{\→}{m}}_{\mathrm{biggest}}=\mathrm{\Σ}\left|{\stackrel{\→}{v}}_k\right|---\left(4\right)$

(即每个零向量第一个元素的绝对值被相加以得到扰动向量的第一元素，而每个零向量的第二元素被相加以得到(扰动向量的)第二元素等等)。然而，L₂范数使用相同元素的平方计算得到。L₃范数使用各分量绝对值的立方计算得到等等，其中p∈[0，∞]。L_∞范数将使用相同元素的最大的绝对值(即每个零向量第一个元素的最大值取作扰动向量的第一个元素，每个向量的第二个元素的最大值取作第二元素等。)。下界是为上界所计算出的总和的负值。下一步，求解方程式(3)得到零向量的组合，这个组合最接近逼近上(下)界，以及通过将结果扰动向量按比例缩放一常数和将该结果加到最可能的解上来完成(求解)。本领域的普通技术人员将知道求解方程式(1)形式的矩阵方程最可能或最佳估测解的方法。例如，Menke在(1984年)由美国学术出版社(Academic Press)出版的《Geophysical Data Analysis：Discrete Inverse Theory》中描述了标准反演方法。按比例缩放用于防止非物理的结果，并优选地以迭代方式递增执行，直到刚好足够的按比例缩放已经用于防止参数m_i物理上并不真实的值。先验信息在特例中或许支持更强的按比例缩放。

这个方法适用于任何问题，其中合适的物理模型可以用于描述所观察的内容和推断内容之间的关系。例如，如果随时间所观察的AVO行为上的差别(时延地震的)可能与油藏中的压力和水饱和度的变化相关，那么，如果给定地震数据之间的所观察的差别，零空间方法就可以用于求解最佳和最坏情况可能的模式。另一个可能的应用是产生历史数据，从历史数据可以推测最好和最坏情况油藏渗透性。通常，测井曲线被处理且单一最佳答案(例如，vshale曲线)被产生。然而，如果无论什么数据被真实地记录下来，该零空间方法则可以用于感兴趣的某个属性(例如，渗透性、水饱和度等)的不同测井曲线数据的联合反演，以产生最佳和最坏情况的可能曲线。构造G矩阵时，当确定油藏中孔隙度和沙子的分布时，各向异性项可以包括在反射和岩石物理方程中。可替换地，方程式可以根据感兴趣的其他属性进行参数化(例如，如阻抗或速度的弹性属性)。

在本发明的一些实施例中，奇异值分解(SVD)被用于将G矩阵分解成两个感兴趣的子空间。这种分解在不同的教科书中都有描述，包括G.Strang，Introduction to Linear Algebra，Wellesley-Cambridge Press，321-333(1993)；C.Lanczos，Linear Differential Operators，Van Nostrand，(1961)；前面所引用的Menke的书；和Press，et al.，Numerical Recipes in C：the art of scientific computing，Cambridge University Press，1992。该运算产生包含数据和模型特征向量的两个正交单位矩阵U和V以及包含奇异值(特征值)的第三个矩阵S。然后执行以下步骤：

1.计算零空间特征值(选用L₁范数)各分量的绝对值总和。该总和通过每个特征向量(沿着特征向量矩阵的行)的相应元素求得，且等效于确定从数据中所提取的全部扰动的上界。下界是该和的负值。

2.选择页岩分离点值。最可能解中对应于大于这个分离点的值的扰动向量元素排除在步骤3之外。

3.求解最接近逼近该上(下)界的特征向量组合。

4.将结果的扰动向量按比例缩放某一常数并将该结果加到最可能的解上。在岩性问题中，比例常数被选择以使得所获得的最终孔隙度和vshale值物理上真实存在(孔隙度或vshale都不小于0％，并且vshale不大于100％，孔隙度不大于40％)。

商业软件产品“Matlab”(www.mathworks.com)使用奇异值分解，容易地将所给G矩阵分解成两个正交单位矩阵U和V以及一个对角阵S，即G＝USV′。S包含奇异值(G′G的特征值)，U包含与数据空间相关的正交单位基向量，V包含与模型空间相关的正交单位基向量。(回忆G将模型参数与数据参数相联系)。Press等人的书包含类似计算机程序的源码清单。典型地，SVD算法将产生S矩阵，该矩阵使其对角元素从最大到最小的顺序排列。然后，选择分离点值并从与低于分离点的值相关的V矩阵中选择基向量是简单的事。这些向量将作为连续列出现在V中，并且人们使用所选范数只是计算这些列向量的相应分量之和。本发明将可使用任意范数。

示例

本发明方法被应用到在潜在油田中所获取的一些地震数据。图2A-C分别示出黑色背景下曲折流过倒vshale体积的，“最坏”、“最可能”和“最佳”(1-vshale)推测的沙道的三维图像(页岩部分已不可见)。包含在“最好”和“最坏”情况模式中的总储油量的估算变化了几乎2倍。此外，在“最坏”和“最好”情况下的沙体之间连通性的差别显著不同，并暗含对于这个矿区可能最佳的不同排水策略。

为解释本发明，前述申请针对本发明的特殊实施例。然而，对本领域普通技术人员明显的是，此处所描述的实施例的许多修改和变化都是可能的。所有的这些修改和变化都在本发明的范围内，正如所附权利要求所定义的。

Claims (6)

1.一种从地下区域生产碳氢化合物的方法，包含：

(a)从所述地下区域的地震勘测或测井勘测中获得所测数据；

(b)获得所述地下区域的模型，该模型将地表区域的物理参数m_i与所测地震数据或测井曲线数据相联系；

(c)获得将物理参数m_i与所测数据相联系的如下形式的矩阵方程：

$\left[G\right]\left[\begin{array}{c}{m}_1\\ m_3\\ .\\ .\\ .\\ m_N\end{array}\right]=\left[\mathrm{data}\right]$

其中G是基于所述模型的矩阵；

(d)从将所述矩阵方程求逆获得物理参数m_i的可能的非唯一解中的最大和最小解，所述最大和最小解通过如下步骤获得：

(i)寻找将所述G矩阵对角化的正交单位基向量，并运用所述正交单位基向量将G矩阵对角化；

(ii)选择一个阈值，低于该阈值的被对角化的G矩阵中元素的值在它们对所述矩阵方程的最可能的解的影响方面被认为不重要；

(iii)从所述正交单位基向量中识别与不重要的对角元素相联系的那些向量，相联系的那些向量即零向量；

(iv)选择L_p数学范数，其中p∈[0，∞]；

(v)根据所选范数，通过计算所述零向量的相应分量之和，确定可能解m₁…m_N的上界和下界，所述下界由所述和的负值给出；

(vi)寻找所述零向量的线性组合，该线性组合最接近逼近所述上界，并为所述下界重复此寻找步骤；和

(vii)依次将所述两个线性组合中的每一个加到所述最可能的解上，以产生最大和最小解，根据需要进行按比例缩放以去掉任何物理上不真实的结果；和

(e)使用所述最大和最小解从所述地下区域进行碳氢化合物的生产。

2.根据权利要求1所述的方法，其中m_i是孔隙度和页岩体积分数的值，data是地震波振幅值。

3.根据权利要求1所述的方法，其中奇异值分解被用于寻找所述正交单位基向量并将所述G矩阵对角化。

4.根据权利要求1所述的方法，其中所述零向量的相应分量的绝对值被求和，对应于选择一个范数L₁。

5.根据权利要求1所述的方法，其中寻找所述零向量的线性组合是求解方程式

${\stackrel{\→}{m}}_{\mathrm{biggest}}=\underset{k}{\mathrm{\Σ}}{\mathrm{\α}}_k{\stackrel{\→}{v}}_k^{\mathrm{null}}$

对于加权因子α_k，使用最小平方方法，其中

是第k个零向量且

是可能的m个解的边界。

6.根据权利要求3所述的方法，其中G是一个非方矩阵，所述正交单位基向量是基向量的两个不同集合，一个集合与data向量的空间相关联，另一集合与m向量的空间相关联，其中所述零向量来自与m向量的空间相关联的基向量。

Images