量化评估生产预测不确定性的系统和方法

摘要

一种更新后验地质模型的系统和方法，其通过整合各种储层数据来支持动态量化数据反演、随机不确定性管理和智能储层管理。

说明书

相关领域的交叉引用

不适用

关于联邦政府资助研究的声明

不适用

技术领域

本发明大体涉及各种储层生产情况下的油气采收估计。更具体的 说，本发明涉及用于量化估计油气储层生产预测不确定性的新的储层 数据整合和反演技术。

背景技术

在大规模油气勘探和开采中，不同开采情况下的采收性能的统计 分布构成了投资决策过程的量化基础。精确的储层预测主要取决于精 确估计动态性能统计数据以及通过储层模拟评估采收性能。通过储层 模拟评估采收性能是动态数据整合和反演的一部分，其通常被称为 “历史拟合”，最终被称为储层生产预测。

历史拟合

历史拟合表示调整储层模型直到其近似地再现油田观测数据的 行为。历史的生产和压力尽可能匹配。历史拟合的精确度取决于储层 模型的质量以及压力和生产数据的质量和数量。当对模型进行历史拟 合后，尤其在根据储层内已知的地质特性而对调整进行限制时，该模 型可被用于对未来储层走向(趋势)进行模拟，并且具有高可信度。

通常，历史拟合由人工执行，且仍在广泛使用二十世纪八十年代 的方法。该方法的过程耗时且有人为偏误的倾向。例如，一次计算一 个参数这种类型的灵敏度研究通常缺少量化精确性，难以核查，且通 常必定会产生单匹配储层模型。对于地质复杂的储层，对于储层模拟 来说，地层和结构的不确定性对动态性能统计数据的影响的量化是必 需的。为了充分考虑这类不确定性，选择覆盖了足够大的范围的不同 地质实现(geological realization)来接近预测统计数据，后者成为量 化储层管理(QRM)的相关基础。

近来，计算上的发展产生了新的技术学科，其通常被称为计算机 辅助历史拟合，其可被分为三种主要方法：确定型、随机型和混合型。 历史拟合是一种可逆的反演问题，具有高度非线性和不适定性。对该 问题的解决方案不可避免地求助于递归算法，递归算法依赖于先验信 息，其可以揭示整组模型(即实现)来拟合动态数据。

过去十年中，随机型方法在油气产业中被广泛熟知，其在领域内 的主要发展在于：i)演化算法(例如，遗传算法或演化策略)；ii)使 用集合卡尔曼滤波的数据同化；以及iii)基于贝叶斯理论的序列蒙特 卡洛或马尔可夫链蒙特卡洛(MCMC)方法。在贝叶斯抽样时，使用 贝叶斯公式来构造目标函数：

$P_{m|d}\left(m\right|d)=\frac{P_{d|m}\left(d\right|m)P_m\left(m\right)}{P_d\left(d\right)}$

此公式将观测数据与先验地质信息结合，其中P_m|d(m|d)， P_d|m(d|m)和P_m(m)分别表示后验、似然度和先验概率分布。归一化因子 P_d(d)表示与数据相关的概率，它与模型参数无关且通常被视为常数。

基于贝叶斯理论的序列蒙特卡洛或马尔可夫链蒙特卡洛 (MCMC)方法为后验概率密度函数(pdf)抽样提供了在统计上最 严谨且最精确的方式，但在直接(模拟)模拟中应用该方法时，由于 高判退率以及必须对提出的每个实现模型运行完整的模拟流程，会产 生令人望而却步的高计算成本。该方法对于线性问题有效，但对于非 线性的计算机辅助历史拟合应用不佳。这一问题已经通过一种两阶方 法的提出得到解决，其具有M-H算法的高接受概率。该方法在线性 问题下运行良好，但在非线性计算机辅助历史拟合的应用下运行不 佳。为了改进上述问题，提出了一种严谨的两阶MCMC方法，该方 法使用借助于流线追踪技术计算得到的灵敏度(即：广义行程时间相 对于给定储层性质的导数)来增加可接受率，并减少计算量。

量化和储层生产不确定性预测

不确定性表示主要由于数据噪声、系统建模误差、以及反演问题 求解过程的非唯一性而存在于动态储层建模中的固有分量。通常，分 析静态数据(例如钻井记录、地核和地震道)和动态数据(例如，生 产历史，瞬变压力测试等)。近来在静态建模技术上的改进使得能够 快速构造大量的地质实现，但是与全物理储层模拟相关的计算成本仍 是惊人的。由于所使用的数据在异质性的长度、标度上都存在差异， 且通常具有不同的精密度，整合不同数据源的数据是很繁琐的工作。 历史拟合中实施的数据整合是数据同化过程的一部分，该过程涉及基 于测量数据来估计模型的参数和动态变量。所述模型指的是“油储层 模型”，并根据基于性能准则来建立生产公式。

在地表下建模中整合附加数据可能会降低不确定性，并被视为是 智能储层管理的重要部分之一。现有技术中不确定性管理和生产预测 的工艺应用组合使用了实验设计、最大熵估计、集合卡尔曼滤波或邻 域算法。邻域算法是已知的用于解决非线性地球物理反演问题的算 法，通常包括两个阶段：i)搜索阶段----在多维参数空间内进行搜索 的方法，其目标是在提供给用户的目标函数的可接受数值内寻找模 型，以及ii)评价阶段----其中搜索阶段中生成的整个模型集合被用于 导出一些贝叶斯测量形式(例如协方差或临界PDF)。

因此需要整合各种储层数据以支持动态量化数据反演、随机不确 定性管理以及智能储层管理。

发明内容

本发明提供了更新后验地质模型的系统和方法，该系统和方法通过整 合各种储层数据以支持动态量化数据反演、随机不确定性管理以及智能储 层管理，满足了上述需求并克服了现有技术中一个或多个缺陷。

一个实施例中，本发明包括一种更新后验地质模型的方法，其包括：i) 为先验地质模型计算新的地质实现；ii)使用先验地质模型的新的地质实现 或者使用先验地质模型的可接受数量的新的地质实现来计算目标函数的精 确似然性；iii)基于目标函数的精确似然性来定义连续的蒙特卡洛链的初始 状态；iv)基于序列蒙特卡洛的初始状态和先验地质模型的随机样本，定义 新样本；v)使用新样本来计算目标函数的近似似然性；vi)如果目标函数的 近似似然性不满足接受准则，则仅基于先验地质模型的另一随机样本重复 上述定义新样本的步骤；vii)如果新样本满足所述接受准则，则使用新样本 来计算目标函数的另一精确似然性；viii)如果目标函数的另一个精确似然 性不满足另一接受准则，则仅基于先验地质模型的另一随机样本重复上述 定义新样本的步骤；ix)仅基于先验地质模型的另一随机样本重复上述定义 新样本的步骤，直到满足收敛准则；以及，x)存储满足接受准则和另一接 受准则的每个新样本，每个新样本表示各个更新的后验地质模型。

另一个实施例中，本发明包括一种程序携载设备，其携载有用于更新 后验地质模型的计算机可执行指令，所述指令可被执行以实现：i)为先验地 质模型计算新的地质实现；ii)使用先验地质模型的新的地质实现或者使用 先验地质模型的可接受数量的新的地质实现来计算目标函数的精确似然 性；iii)基于目标函数的精确似然性来定义序列蒙特卡洛链的初始状态；iv) 基于序列蒙特卡洛的初始状态和先验地质模型的随机样本，定义新样本； v)使用新样本来计算目标函数的近似似然性；vi)如果目标函数的近似似然 性不满足接受准则，则仅基于先验地质模型的另一随机样本重复上述定义 新样本的步骤；vii)如果新样本满足接受准则，则使用新样本来计算目标函 数的另一精确似然性；viii)如果目标函数的另一精确似然性不满足另一接 受准则，则仅基于先验地质模型的另一随机样本重复上述定义新样本的步 骤；ix)仅基于先验地质模型的另一随机样本重复上述定义新样本的步骤， 直到满足收敛准则；x)存储满足接受准则和另一接受准则的每个新样本， 每个新样本表示各个更新的后验地质模型。

根据下文中对各个实施例和相关附图的说明，本发明的其他方面、优 点和实施例对本领域技术人员来说将显而易见。

附图说明

下面参照附图来具体描述本发明，其中相同的元件标示有相同的参考 标记，其中：

图1是执行本发明的方法的一个实施例的流程图。

图2是执行本发明的方法的另一个实施例的流程图。

图3是执行图1和2中步骤103和203的方法的一个实施例的流程图。

图4是执行图3中步骤306的方法的一个实施例的流程图。

图5是分别执行图1和2中步骤114和204，217的方法的一个实施例的 流程图。

图6A是图3中步骤305的3D岩石物性实现的生成的多个2D图像。

图6B是图4中步骤404的3D结果的多个2D图像。

图6C是图4中步骤407中存储的3D结果的多个2D图像。

图7是执行本发明的系统的一个实施例的方框图。

具体实施方式

以具体说明的方式描述了本发明的主题，然而，该说明本身不旨 在限制本发明的范围。因此，能够以其他方式实施本发明的主题，从 而包括类似于在此描述的、与其他技术结合的不同步骤或者步骤的组 合。此外，尽管在此使用术语“步骤”来描述应用的方法的不同要素， 该术语不应该被解释为暗示在此描述的不同步骤之中或之间的任何 特定顺序，除非另外通过说明明确限定了特定顺序。虽然下述说明涉 及石油和天然气产业，但是本发明的系统和方法并不限于此，还可以 被应用到其他产业而获得相同的结果。

方法说明

下述方法，根据流线模拟，提出了一种在不确定性的动态量化中 的积分闭环方式，作为地质模型的再筛选技术。所述方法将快速数据 反演算法与高效模型参数化，以及用于在储层管理过程中快速作出决 策的排序技术进行结合。

现在参看图1，其为例示了执行本发明的方法100的实施例的流程图。

步骤101，对于方法100，使用参照图7所描述的客户界面和/或视频界 面来输入先验地质模型的参数。所述参数可包括，例如渗透率，其可以3D 阵列[Nx，Ny，Nz]的形式被输入。

步骤102，使用参照图7所描述的客户界面和/或视频界面，为先验地 质模型定义多个(N个)新的地质实现。

步骤103，以参照图3所进一步描述的方式为先验地质模型计算所述(N) 个新的地质实现。

步骤104，使用计算出的新的地质实现来计算先验地质模型的目标函数 的精确似然性。例如，一种技术通过采用本领域已知技术以及为先验地质 模型计算出的多(N)个新的地质实现来执行正演模拟、流线追踪和导数计算 (calculation of derivatives)，以获得模型响应g(m)、数据错位(d_p-g(m))以及 生产数据灵敏度系数阵列Sp∈R^[Nx*Ny*Nz]*p，从而计算目标函数的精确似然性， 其中d_p＝d(p)是观测数据矢量，p是开采油井的数量。例如，另一种技术可通 过采用本领域已知技术以及为先验地质模型计算出的多(N)个新的地质实现 来执行正演模拟、流线追踪和导数计算，以获得模型响应g(m)、数据错位 (d_p-g(m))以及灵敏度系数阵列Sp∈R^[Nx*Ny*Nz]*p和灵敏度系数阵列 Ss∈R^[Nx*Ny*Nz]*S，从而计算目标函数的精确似然性，其中d_p＝d(p)是观测数据 矢量，p是开采油井的数量，其中d_s＝d(s)是观测(时移)地震数据矢量，s是(时 移)地震数据控制点的数量。在流体饱和度表格和压力示意图中使用的地震 数据可以通过例如地震反演导出。这些图提供了除生产数据之外的独立约 束集合。就生产数据P以及(时移)地震数据S而论，通过结合灵敏度导数，目 标函数的精确似然性的引伸项具有下列贝叶斯公式：

P_m|d(m|d)∝exp[-O(m)]

其中$O\left(m\right)=\frac{1}{2}{(d_p-g\left(m\right))}^T{C}_P^{-1}(d_p-g\left(m\right))+\frac{1}{2}{(d_S-g\left(m\right))}^T{C}_S^{-1}(d_s-g\left(m\right))$

Cp和Cs是分别与生产数据P和(时移)地震数据S相对应的协方差矩 阵。

步骤105，定义序列蒙特卡洛(MC)链的初始状态。所提出的该序列 MC链的初始状态的一种定义是例如mⁱ；i＝1：I，其中I等于序列MC链的 重复次数。

步骤106，使用随机行走取样器mⁱ⁺¹＝mⁱ+σε对条件概率密度函数(pdf) q(m^*|mⁱ)取样，以定义新样本，其中σ是链步长(chain step size)，ε是先 验地质模型的随机样本。可选择地，也可使用其他已知取样器，例如独立 取样器或朗之万(Langevin)取样器。

步骤107，使用步骤106中定义的新样本来计算目标函数的近似似然性。 例如，一种技术通过提出跃迂δm＝m^*-mⁱ来计算目标函数的近似似然性， 该跃迂δm＝m^*-mⁱ使用生产数据灵敏度矩阵S_p定义了在正演模型响应期 间的变化；δd_p＝S_P·δm。例如，另一种技术通过提出跃迂δm＝m^*-mⁱ来计 算目标函数的近似似然性，该跃迂δm＝m^*-mⁱ使用生产数据灵敏度矩阵 S_p；δd_p＝S_P·δm以及地震数据灵敏度矩阵S_S；δd_S＝S_S·δm定义了在正演模 型响应期间的变化。

步骤108，方法100判定是否满足接受准则α₁(m^*，mⁱ)。如果满足接受 准则，方法100则继续执行步骤110。如果不满足接受准则，方法100则继 续执行步骤108a。例如，使用本领域已知技术Matropolis-Hastmgs准则来定 义所述接受准则。从开区间U(0，1)上的伪随机数的标准均匀分布中抽样数 值，比较该接受准则α₁(m^*，mⁱ)与该抽样数值。如果该接受准则α₁(m^*，mⁱ)大 于该抽样数值，则将新样本(mⁱ)提升到序列MC链的推荐状态(即， mⁱ⁺¹＝m^*)，且方法100继续执行步骤110。如果该接受准则α₁(m^*，mⁱ)小于 该抽样数值，则方法100最终返回到步骤106，在该步骤106中定义另一个 新样本。

步骤108a，从先验地质模型中选择另一个随机样本，方法100返回到 步骤106以定义另一个新样本。因此，方法100通过步骤106，107，108和108a 反复地进行，直到满足步骤108中的接受准则。

步骤110，当满足步骤108的接受准则后，使用步骤106中定义的最后 一个新样本来计算目标函数的精确似然性。例如，一种技术通过本领域已 知技术以及步骤106中定义的最后一个新样本来执行正演模拟、流线追踪 和导数计算，以获得模型响应g(m^*)、数据错位(d_p-g(m^*))以及生产数据灵敏 度系数阵列Sp∈R^[Nx*Ny*Nz]*p，从而计算目标函数的精确似然性，其中d_p＝d(p) 是观测数据矢量，p是开采油井的数量。例如，另一种技术通过本领域已知 技术以及步骤106中定义的最后一个新样本来执行正演模拟、流线追踪和 导数计算，以获得模型响应g(m^*)、数据错位(d_p-g(m^*))以及灵敏度系数阵列 Sp∈R^[Nx*Ny*Nz]*p和灵敏度系数阵列Ss∈R^[Nx*Ny*Nz]*S，从而计算目标函数的精确 似然性，其中d_p＝d(p)是观测数据矢量，p是开采油井的数量，，其中d_s＝d(s) 是观测(时移)地震数据矢量，s是(时移)地震数据控制点的数量。在流体 饱和度表格和压力示意图中使用的地震数据可以通过例如地震反演导出。 这些图提供了除开采数据之外的独立约束集合。就生产数据P以及(时移)地 震数据S而论，通过结合灵敏度导数，目标函数的精确似然性的引伸项具有 下列贝叶斯公式：

P_m|d(m|d)∝exp[-O(m)]

其中$O\left(m\right)=\frac{1}{2}{(d_p-g\left(m\right))}^T{C}_P^{-1}(d_p-g\left(m\right))+\frac{1}{2}{(d_S-g\left(m\right))}^T{C}_S^{-1}(d_s-g\left(m\right))$

Cp和Cs是分别与生产数据P和(时移)地震数据S相对应的协方差矩 阵。

步骤111，方法100判定是否满足接受准则α₂(m^*，mⁱ)。如果满足接受 准则，方法100则继续执行步骤113。如果不满足接受准则，方法100则继 续执行步骤111a。例如，使用本领域已知技术Matropolis-Hastings准则来定 义接受准则α₂(m^*，mⁱ)。在开区间U(0，1)从伪随机数的标准均匀分布中对数 值抽样，比较该接受准则α₂(m^*，mⁱ)与该抽样数值。如果该接受准则 α₂(m^*，mⁱ)大于该抽样数值，则将新样本(mⁱ)接受为更新后的后验样本(即， mⁱ⁺¹＝m^*)，且方法100继续执行步骤113。如果该接受准则α₂(m^*，mⁱ)小于 该抽样数值，方法100则最终返回到步骤106，在该步骤106中定义另一个 新样本。

步骤111a，从先验地质模型中选择另一个随机样本，方法100返回到 步骤106以定义另一个新样本。因此，方法100通过步骤106，107，108，110，111 和111a反复地进行，直到满足步骤111中的接受准则。

步骤113，方法100判定是否满足收敛准则。如果满足收敛准则，方法 100则继续执行步骤114。如果不满足收敛准则，方法100则继续执行步骤 113a。例如，一种收敛准则利用了最大熵准则或方法。例如，可将熵S定义 为S＝-＜p_m|d(m|d)log p_m|d(m|d))＞，其中p_m|d(m|d)是更新后的后验地质模 型。

步骤113a，从先验地质模型中选择另一个随机样本，且方法100返回 到步骤106以定义另一个新样本。因此，方法100通过步骤 106，107，108，110，111，113和113a反复地进行，直到满足步骤113中的收敛准 则。

步骤114，将满足步骤108，111中的接受准则和步骤113中的收敛准则 的每个样本存储为更新后的后验地质模型，并以参照图5进一步描述的方 式对其进行动态排序。

步骤115，使用参照图7所描述的客户界面和/或视频界面来评估排序后 的后验地质模型，这使得能够基于评估结果执行一或多个商业决策。例如， 可基于选择多个排序后的后验地质模型进行评估。可以根据相应的最终恢 复因子(URF)选择特定数量的分级排序后的后验地质模型，并将其转发给操 作单元，从而例如进一步将其考虑入闭环储层管理流程中或者执行油井设 置优化过程。

现在参看图2，其例示了执行本发明的方法200的另一个实施例的流程 图。

步骤201，对于方法200，使用参照图7所描述的客户界面和/或视频界 面来输入先验地质模型的参数。所述参数可包括，例如渗透率，其可以3D 阵列[Nx，Ny，Nz]的形式输入。

步骤202，使用参照图7所描述的客户界面和/或视频界面为先验地质 模型定义多个(N个)新的地质实现。

步骤203，以参照图3所描述的方式为先验地质模型计算多个(N个)新 的地质实现。

步骤204，以参照图5所进一步描述方式对步骤203中计算得到的新的 地质实现进行动态排序。

步骤205，使用参照图7所描述的客户界面和/或视频界面，基于可接 受的新的地质实现的排序来对其进行选择。尽管必须至少选择一个新的地 质实现，优选为不选择所有的新的地质实现。可接受的新的地质实现可为 例如在一些用户设定阈值之上的新的地质实现。

步骤206，使用所选择的新的地质实现为先验地质模型计算目标函数的 精确似然性。例如，一种技术通过本领域已知技术以及所选择的新的地质 实现来执行正演模拟，流线追踪和导数计算，以获得模型响应g(m)、数据错 位(d_p-g(m))以及生产数据灵敏度系数阵列Sp∈R^[Nx*Ny*Nz]*p，从而计算目标函 数的精确似然性，其中d_p＝d(p)是观测数据矢量，p是开采油井的数量。例如， 另一种技术通过本领域已知技术以及所选择的新的地质实现来执行正演模 拟、流线追踪和导数计算，以获得模型响应g(m)、数据错位(d_p-g(m))以及灵 敏度系数阵列Sp∈R^[Nx*Ny*Nz]*p和灵敏度系数阵列Ss∈R^[Nx*Ny*Nz]*S，从而计算目 标函数的精确似然性，其中d_p＝d(p)是观测数据矢量，p是开采油井的数量，， 其中d_s＝d(s)是观测(时移)地震数据矢量，s是(时移)地震数据控制点的数量。 在流体饱和度表格和压力示意图中使用的地震数据可以通过地震反演导 出。这些图提供了除生产数据之外的独立约束集合。就生产数据P以及(时移) 地震数据S而论，通过结合灵敏度导数，目标函数的精确似然性的引伸项具 有下列贝叶斯公式：

P_m|d(m|d)∝exp[-O(m)]

其中$O\left(m\right)=\frac{1}{2}{(d_p-g\left(m\right))}^T{C}_P^{-1}(d_p-g\left(m\right))+\frac{1}{2}{(d_S-g\left(m\right))}^T{C}_S^{-1}(d_s-g\left(m\right))$

Cp和Cs是分别与生产数据P和(时移)地震数据S相对应的协方差矩 阵。

步骤207，定义序列蒙特卡洛(MC)链的初始状态。该序列MC链的 初始状态的一种定义是例如mⁱ；i＝1：I，其中I等于序列MC链的重复次数。

步骤208，使用随机行走取样器mⁱ⁺¹＝mⁱ+σε对条件概率密度函数(pdf) q(m^*|mⁱ)取样以定义新样本，其中σ是链步长(chain step size)，ε是先验 地质模型的随机样本。可选择的，也可使用其他已知取样器，例如独立取 样器或朗之万Langevin取样器。

步骤209，使用步骤206中定义的新样本来计算目标函数的近似似然性。 例如，一种技术通过提出跃迂δm＝m^*-mⁱ来计算目标函数的近似似然性， 该跃迂δm＝m^*-mⁱ使用生产数据灵敏度矩阵S_p；δd_p＝S_P·δm定义了在正演 模型响应期间的变化；δd_p＝S_P·δm。例如，另一种技术通过提出跃迂 δm＝m^*-mⁱ来计算目标函数的近似似然性，该跃迂δm＝m^*-mⁱ使用生产数 据灵敏度矩阵S_p；δd_p＝S_P·δm以及灵敏度矩阵S_S；δd_S＝S_S·δm定义了在正 演模型响应期间的变化。

步骤210，方法200判定是否满足接受准则α₁(m^*，mⁱ)。如果满足接受 准则，方法200则继续执行步骤212。如果不满足接受准则，方法200则执 行步骤210a。可使用本领域已知技术Matropolis-Hastings准则来定义所述接 受准则。在开区间U(0，1)从伪随机数的标准均匀分布中对数值抽样，比较该 接受准则α₁(m^*，mⁱ)与该抽样数值。如果该接受准则α₁(m^*，mⁱ)大于该抽样数 值，则将新样本(mⁱ)提升到序列MC链的推荐状态(即，mⁱ⁺¹＝m^*)，且方 法200继续执行步骤212。如果该接受准则α₁(m^*，mⁱ)小于该抽样数值，方 法200最终返回到步骤208，在该步骤208中定义另一个新样本。

步骤210a，从先验地质模型中选择另一个随机样本，方法200返回到 步骤208以定义另一个新样本。因此，方法200通过步骤208，209，210和 210a反复进行，直到满足步骤210中的接受准则。

步骤212，当满足步骤210的接受准则后，使用步骤208中定义的最后 一个新样本来计算目标函数的精确似然性。例如，一种技术通过本领域已 知技术以及步骤208中定义的最后一个新样本来执行正演模拟、流线追踪 和导数计算，以获得模型响应g(m^*)、数据错位(d_p-g(m^*))以及生产数据灵敏 度系数阵列Sp∈R[^Nx*Ny*Nz]*p，从而计算目标函数的精确似然性，其中d_p＝d(p) 是观测数据矢量，p是开采油井的数量。例如，另一种技术通过本领域已知 技术以及步骤208中定义的最后一个新样本来执行正演模拟、流线追踪和 导数计算，以获得模型响应g(m^*)、数据错位(d_p-g(m^*))以及灵敏度系数阵列 Sp∈R^[Nx*Ny*Nz]*p和灵敏度系数阵列Ss∈R^[Nx*Ny*Nz]*S，从而计算目标函数的精确 似然性，其中d_p＝d(p)是观测数据矢量，p是开采油井的数量，其中d_s＝d(s)是 观测(时移)地震数据矢量，s是(时移)地震数据控制点的数量。在流体饱 和度表格和压力示意图中使用的地震数据可以通过地震反演导出。这些图 提供了除生产数据之外的独立约束集合。就生产数据P以及(时移)地震数据S 而论，通过结合灵敏度导数，目标函数的精确似然性的引伸项具有下列贝 叶斯公式：

P_m|d(m|d)∝exp[-O(m)]

其中$O\left(m\right)=\frac{1}{2}{(d_p-g\left(m\right))}^T{C}_P^{-1}(d_p-g\left(m\right))+\frac{1}{2}{(d_S-g\left(m\right))}^T{C}_S^{-1}(d_s-g\left(m\right))$

Cp和Cs是分别与生产数据P和(时移)地震数据S相对应的协方差矩 阵。

步骤213，方法200判定是否满足接受准则α₂(m^*，mⁱ)。如果满足接受 准则，方法200则继续执行步骤215。如果不满足接受准则，方法200则继 续执行步骤213a。可使用本领域已知技术Matropolis-Hastings准则来定义接 受准则α₂(m^*，mⁱ)。在开区间U(0，1)从伪随机数的标准均匀分布中对数值抽 样，比较该接受准则α₂(m^*，mⁱ)与该抽样数值。如果该接受准则α₂(m^*，mⁱ)大 于该抽样数值，则将新样本(mⁱ)接受为更新后的后验样本(即，mⁱ⁺¹＝m^*)， 且方法200继续执行步骤215。如果该接受准则α₂(m^*，mⁱ)小于该抽样数值， 方法200则最终返回到步骤208，在该步骤208中定义另一个新样本。

步骤213a，从先验地质模型中选择另一个随机样本，方法200返回到 步骤208以定义另一个新样本。因此，方法200通过步骤208，209，210，212，213 和213a反复地进行，直到满足步骤213中的接受准则。

步骤215，方法200判定是否满足收敛准则。如果满足收敛准则，方法 200则继续执行步骤216。如果不满足收敛准则，方法200则继续执行步骤 215a。例如，一种收敛准则利用了最大熵准则或方法。例如，可将熵S定 义为S＝-＜p_m|d(m|d)logp_m|d(m|d))＞，其中p_m|d(m|d)是更新后的后验地质 模型。

步骤215a，从先验地质模型中选择另一个随机样本，且方法200返回 到步骤208以定义另一个新样本。因此，方法200通过步骤 208，209，210，212，213，215和215a反复地进行，直到满足步骤215中的收敛准 则。

步骤216，方法200判定是否执行可接受样本(更新后验地质模型)的 动态排序。如果需要对更新后的后验地质模型执行动态排序，方法200则 继续执行步骤217。如果不需要对更新后的后验地质模型执行动态排序，方 法200则继续执行步骤220。基于主观准则决定动态排序的执行，所述主观 准则为例如使用参照图7所描述的客户界面和/或视频界面来确定更新后的 后验地质模型的数量。

步骤217，将满足步骤210、213中的接受准则和步骤215中的收敛准 则的每个样本存储为更新后的后验地质模型，并以参照图5进一步描述的 方式对其进行动态排序。

步骤218，使用参照图7所描述的客户界面和/或视频界面来选择最佳 的排序后的后验地质模型。优选地，排位最高的后验地质模型是最佳的。

步骤219，用最佳的排序后的后验地质模型来替换先验地质模型，并将 其重命名为“先验地质模型”。方法200返回到步骤202以为先验地质模型 定义多个(N个)其他新的地质实现。

步骤220，存储更新后的后验地质模型，并使用参照图7所描述的客户 界面和/或视频界面、基于更新后的后验地质模型来执行商业决策。例如， 将更新后的后验地质模型传送给操作单元后，可进一步考虑将其并入闭环 储层管理流程中或者执行油井设置优化任务。

现在参看图3，其是例示了分别执行图1和2中步骤103和203的方法 的一个实施例的流程图。

步骤301，使用参照图7所描述的客户界面和/或视频界面为DESM_TM定 义框架。例如，可通过将步骤101获得的参数输入DESM_TM中来定义该框架。

步骤302，使用参照图7所描述的客户界面和/或视频界面来定义岩相 和岩石物性计算的变差函数参数。在DESM_TM中使用用于变差函数建模和计 算的技术和算法来初始计算变差函数参数。DESM_TM使用已知的克里格法的 变量来计算变差函数参数。

步骤303，使用参照图7所描述的客户界面和/或视频界面来定义岩相 和岩石物性计算的随机数/种子。例如，已知的标准随机数生成器可用于将 随机数/因子输入DESM_TM。

步骤304，使用步骤301-301中定义的参数来计算新的岩相实现。例 如，可通过已知技术-多高斯(Pluri Gaussian)模拟来计算该岩相实现。因 此，计算得到单个新的3D岩相实现，其形式为单个列矢量，具有维度[Nx^*，Ny^*，Nz^*]。

步骤305，使用上述新的岩相实现对3D阵列[Nx，Ny，Nz]形式的岩石物 性实现的生成进行约束。例如，可以通过已知技术-回转带(Turning bands) 模拟来生成岩石物性实现。因此，当生成岩石物性实现时，将岩相实现作 为约束。由此计算得到单个新3D岩相实现，其形式为单个列矢量，具有维 度[Nx^*，Ny^*，Nz^*]。图6A中用多幅表示原始渗透油田的不同层的2D图像例 示了生成3D岩石物性实现的示例。用于构造图6A中图像的数据是由被称 为Brugge综合数据集的已知数据源获取的，该数据源是对公众开放的。

步骤306，以参照图4进一步描述的方式对多(N)个新的地质实现中 的一个地质实现执行模型参数化。

步骤307，方法300判定是否每个新的地质实现都被参数化。如果每个 新的地质实现都已被参数化，方法300则继续执行步骤308。如果不是每个 新的地质实现都已被参数化，方法300则继续执行步骤307a。

步骤307a，(n)加一(1)，直到(n)等于所定义的新的地质实现的数量(N)， 方法200返回到步骤303以定义另一个随机数/因子用于模拟。因此，方法 300通过步骤303，304，305，306，307和307a反复地进行，直到每个新的地质 实现都已根据步骤306被参数化。

步骤308，将每个参数化后的新的地质实现存储为新的地质实现，然后 返回到步骤104或204。

现在参看图4，其是例示了执行图3中步骤306的方法400的一个实施 例的流程图。该方法利用了离散余弦变换(DCT)。DCT是目前用于图像压 缩和模式识别的最成功的变换方法，其被认为是使用历史拟合和集合卡尔 曼滤波(EnKF)方法对渗透域的新实现进行参数化的有效方法。

步骤401，使用参照图7所描述的客户界面和/或视频界面来定义被保 留的3D DCT模式的数量(NMODES)。基于主观判断最佳数量来确定3D DCT模式的数量(NMODES)。

步骤402，使用步骤305中生成的岩石物性实现、通过以下方式计算 3D DCT的一组基函数(α(u)，α(v)，α(w))：

$\alpha \left(u\right)=\left(\begin{array}{cc}\sqrt{\frac{1}{N_x}}& u=0\\ \sqrt{\frac{2}{N_x}}& u\ne 0\end{array}\right)$$\alpha \left(v\right)=\left(\begin{array}{cc}\sqrt{\frac{1}{N_y}}& v=0\\ \sqrt{\frac{2}{N_y}}& v\ne 0\end{array}\right)$$\alpha \left(w\right)=\left(\begin{array}{cc}\sqrt{\frac{1}{N_z}}& w=0\\ \sqrt{\frac{2}{N_z}}& w\ne 0\end{array}\right)$

因此，Nx，Ny和Nz对应于步骤305中由3D阵列表示的岩石物性实现。

步骤403，使用3D DCT的该组基函数α(u)，α(v)，α(w)通过以下方式来 计算3D DCT的系数(C(u，v，w))：

$C(u,v,w)=(a\left(u\right)a\left(v\right)a\left(w\right)\underset{x=0}{\overset{\mathrm{Nx}-1}{\Sigma }}\underset{y=0}{\overset{\mathrm{Ny}-1}{\Sigma }}\underset{z=0}{\overset{\mathrm{Nz}-1}{\Sigma }}f(x,y,z)\cos \left(\frac{\pi (2x+1)u}{2*\mathrm{Nx}}\right)\cos \left(\frac{\pi (2y+1)v}{2*\mathrm{Ny}}\right)\cos \left(\frac{\pi (2z+1)w}{2*\mathrm{Nz}}\right)$

步骤404，通过选择阈值来执行3D DCT的系数截取，从而为3D DCT 保持NMODES个重要系数。用于保持NMODES个重要系数的阈值的选择 是通过以下两个步骤实现的：i)以升序对3D DCT系数的绝对值阵列进行 排序；ii)通过将(Nx*Ny*Nz-NMODES)的分量下舍取整为最接近的整 数来计算该阈值。在图6B所示的3D DCT的系数截取示例中，多幅2D图 像表示截取后的3％的log 3D DCT的不同层。用于构造图6B中的图像的数 据是由与构造图6A中图像的数据源相同的数据源获得的。

步骤405，使用本领域已知技术对不重要的3DDCT系数进行归零化， 从而丢弃不重要的3D DCT系数。可以例如通过以下方式丢弃不重要的3D DCT系数：i)保留步骤404获得的NMODES个重要系数；ii)将步骤403 中计算出的3DDCT系数映射到该NMODES个重要系数。

步骤406，对其余的截取后的3D DCT系数执行逆向3D DCT：

$f(x,y,z)=\underset{u=0}{\overset{\mathrm{Nx}-1}{\Sigma }}\underset{v=0}{\overset{\mathrm{Ny}-1}{\Sigma }}\underset{w=0}{\overset{\mathrm{Nz}-1}{\Sigma }}a\left(u\right)a\left(v\right)a\left(w\right)C(u,v,w)\cos \left(\frac{\pi (2x+1)u}{2*\mathrm{Nx}}\right)\cos \left(\frac{\pi (2y+1)v}{2*\mathrm{Ny}}\right)\cos \left(\frac{\pi (2z+1)w}{2*\mathrm{Nz}}\right)$

步骤407，将逆向3DDCT存储为参数化地质实现，并返回到图3中的 步骤307。图6C中，通过多幅表示具有3％的3D DCT系数的、恢复的渗 透域的不同层的2D图像，示出了步骤406中的逆向3D DCT的结果。用于 构造图6C中图像的数据是由与构造图6A，6B中图像的数据源相同的数据 源获得的。比较图6A和6C可看出，图6A和6C中同一层的每幅图像看起 来几乎相同。图6A和6C看起来几乎相同的事实突出显示了3DDCT模型 参数化方法400的精确性和效率。其还表明图6A中的原始渗透域可以仅通 过作为方法400的结果的、被保留的3％的3DDCT系数来恢复。

现在参看图5，其为例示了执行图1和2中步骤114和204，217的方法 的一个实施例的流程图。

步骤501，使用本领域已知技术对更新后的后验地质模型(步骤114， 217)或新的地质实现(步骤204)执行流线模拟，其产生动态系统响应， 例如恢复因子。使用本领域已知技术执行所述流线模拟。

步骤502，使用步骤501中计算出的动态系统响应数据来生成模式不一 致性距离矩阵(D)。例如，在一个实施例中，通过以下方式生成模式不一 致性距离矩阵：

(D)＝欧几里得空间E动态响应的{δ_ij}，其中i，j＝1，...N。

步骤503，执行多维度量化(MDS)，使用例如核主成分分析或K均值 聚类法的本领域已知技术将上述模式不一致性距离矩阵映射到线性高维度 特征空间(F)。

步骤505，使用作为现有技术中已知技术的核方法(Kernel method)来 划分特征空间(F)的数据点(x)。例如，在一个实施例中，可通过以下方 式划分特征空间(F)的数据点(x)：i)定义聚类(k)的先验数量；ii)将Nc 个数据点(x)随机划分到聚类(k)中；iii)计算聚类中心(Θ_i)。

步骤506，使用本领域已知技术来计算数据点(x)的每个部分的聚类 中心(Θ_i)。

步骤507，通过优化过程、基于邻近度将数据点(x)分配给聚类中心 (Θ_i)。例如，在一个实施例中，所述优化过程可以最小化目标函数O(x，Θ_i)的 “均方距离”，其中：$O=\underset{i=1}{\overset{k}{\Sigma }}\underset{j=1}{\overset{\mathrm{Nc}}{\Sigma }}{|x_j-{\Theta }_i|}^2$

步骤508，方法500判定聚类中心(Θ_i)是否保持恒定。如果聚类中心(Θ_i) 保持恒定，方法500则继续执行步骤509。如果聚类中心(Θ_i)不保持恒定， 方法500则继续执行步骤508a。

步骤508a，对于步骤507中的优化步骤，i、j加1，然后方法500返回 到步骤507。因此，方法500通过步骤507，508和508a反复地进行，直到聚 类中心(Θ_i)保持恒定。

步骤509，使用本领域已知技术对最接近每个类别中心(Θ_i)的每个更新 后的后验地质模型或每个新的地质实现执行全物理模拟，所述已知技术为 例如有限导数正演储层模拟器。以这种方式，根据全物理模拟的结果，对 更新后的后验地质模型或新的地质实现进行动态排序，并将排序结果返回 到步骤115，步骤205或步骤218。

系统说明

本发明可通过计算机可执行指令程序执行，例如程序模块，通常被称 为由计算机执行的软件应用或应用程序。所述软件可包括例如执行特定任 务或实施特定抽象数据类型的例程、程序、对象、组件和数据结构。软件 形成一允许计算机根据输入源作出响应的界面。Landmark Graphics  Corporation公司销售的商业软件应用-决策空间地球建模(DESM_TM)和 可以和其他第三方应用(例如流线模拟和DESTINY)一起用作执行 本发明的界面应用程序。该软件还和其他代码段协作，与接收数据的源一 起响应接收到的数据，启动多个任务。软件可存储和/或携载在任何种类的 存储介质上，例如CD-ROM、磁盘、磁泡存储器和半导体存储器(例如各 种类型的RAM或ROM)。另外，软件和其结果可以通过多种载体介质传输， 例如光纤、金属线，和/或通过多种网络传输，例如因特网。

而且，本领域技术人员应理解，本发明可以用多种计算机系统设置实 施，包括手持设备、多处理器系统、基于微处理器或可编程消耗性电子设 备、微机、大型计算机以及类似设备。任何数量的计算机系统和计算机网 络都可用于本发明。本发明还可应用于借助通信网络连接的远程处理设备 来执行任务的分布式计算环境。在分布式计算环境中，程序模块可位于包 括存储器设备的本地和远程计算机存储介质上。因此可以结合各种硬件、 软件或其组合，在计算机系统或其他处理系统中执行本发明。

现在参看图7，其为例示了一种用于在计算机上执行本发明的系统的流 程图。该系统包括一计算单元，有时称为计算系统，该单元包括存储器、 应用程序、客户界面、视频界面和处理单元。计算单元仅是适合的计算环 境的一个示例，不应被认为是对本发明的应用或功能性范围产生任何限制。

存储器主要存储应用程序，后者还被描述为包含计算机可执行指令的 程序模块，由计算单元执行该应用程序来运行在此处描述的且在图1-5中例 示出的本发明。因此该存储器包括能够执行所例示的、且参照图1和2进 行过描述的方法的反演和排序平台，并整合了图7所示的其余应用程序的 功能性。例如，所述反演和排序平台集成了DESM_TM来分别执行参照图1和 2中步骤103和203所述的功能。尤其，DESM_TM可用于执行参照图3中步 骤301-305所描述的功能，而反演和排序平台用于执行参照图3中步骤306 -308和图4中步骤401-407所描述的功能。例如，所述反演和排序平台 还整合了第三方流线模拟应用程序(例如，由StreamSim Technologies，Inc提 供的3DSL)来执行参照图5中步骤501所描述的功能，而反演和排序平 台用于执行参照图5中步骤502-508所描述的功能。例如，所述反演和排 序平台还整合了由Texas A&M大学和的联合工业项目开发的应用程 序-DESTINY来执行参照图1中步骤104，110以及图2中步骤206，212 所描述的功能。还可用于执行参照图5中步骤509所描述的功能。

尽管所示的计算单元有一个通用存储器，该计算单元通常包括多种计 算机可读介质。举例来说，但不作为限制，计算机可读介质可包括计算机 存储介质。计算机系统存储器可包括：易失性或非易失性存储器形式的计 算机存储介质，例如只读存储器(ROM)和随机存取存储器(RAM)。基 本输入/输出系统(BIOS)通常存储在ROM中，其包含例如在启动过程中 有助于计算单元部件之间的信息传送的基本例程。RAM通常包含可立即由 处理单元获取和/或目前正由处理单元操作的数据和/或程序模块。举例来 说，但不作为限制，计算单元包括操作系统、应用程序、其他程序模块和 程序数据。

存储器中示出的组件还可包括在其他可移动/不可移动、易失性/非易失 性计算机存储介质上，或者这些组件可以通过应用程序界面(API)在计算 单元上执行，该应用程序界面位于独立的计算单元中，连接到计算系统或 网络上。仅作为示例，硬盘驱动器可以读取不可移动的、非易失性磁性介 质或向其写入，磁盘驱动器可以读取可移动的、易失性磁盘或向其写入， 光盘驱动器可以读取可移动的、非易失性光盘或向其写入，例如CD ROM 或其他光学介质。其他可用于该示例性操作环境的可移动/不可移动、易失 性/非易失性计算机存储介质包括但不限于：磁带、闪存卡、数字通用盘、 数字视频磁带、固态RAM、固态ROM以及类似介质。所述驱动器和上述 与其相关联的计算机存储介质用于存储计算机可读指令、数据结构、程序 模块以及用于计算单元的其他数据。

客户可通过客户界面把指令和信息输入计算单元，该界面可以是输入 设备，例如键盘和定点设备，通常指的是鼠标、跟踪球或触摸板。输入设 备可以包括麦克风、操纵杆、碟形卫星天线、扫描仪或类似设备。这些和 其他输入设备通常通过系统总线连接到处理单元，但也可通过其他接口和 总线结构连接，例如并行端口或通用串行总线(USB)。

监视器或其他类型的显示设备可通过接口，例如视频接口，连接到系 统总线。图形用户界面(GUI)还可以和视频界面一起接收客户界面的指令， 并将指令传输到处理单元。除了监视器之外，计算机还包括通过外围输出 接口连接的其他外围输出设备，例如扬声器和打印机。

尽管未示出计算单元的许多其他内部组成部分，本领域技术人员应理 解这些组成部分和其相互连接是公知的。

尽管已结合目前的优选实施例阐述了本发明，本领域技术人员应当能 够理解本发明不局限于这些实施例。因此，在不偏离由所述权利要求书或 其等同物所限定的本发明的主旨和范围的情况下，可以预见到各种替换实 施例并可对所公开实施例进行修改。