具有两个或更多个不混溶相的岩石的岩石物理性质的直接数值模拟的制作方法

本申请要求2016年10月31日提交的美国临时申请号62/415,050的优先权。该申请的整个公开内容通过引用并入本文中。

背景技术：

在烃的生产中，获得岩层的岩石物理性质的准确地下估计对于评估岩层中含有的烃体积和制定从岩层提取烃的策略来说是重要的。传统上，对岩层样品(例如来自岩心样品或钻屑)进行物理实验室测试以测量岩石物理性质，如渗透率、孔隙率、地层因子、弹性模量等。这些测量中的一些需要很长时间，在一些情况下延续几个月，这取决于岩石本身的性质。用于进行这些测量的设备也可能非常昂贵。

由于直接测量岩石物理性质所需的成本和时间，可以应用“直接数值模拟”技术来有效地估计岩石样品的物理性质，如孔隙率、绝对渗透率、相对渗透率、地层因子、弹性模量等，所述岩石样品包括来自难以处理的岩石类型如致密气砂岩或碳酸盐岩的样品。根据这种方法，例如通过计算机断层摄影(ct)扫描获得岩石样品的三维断层摄影图像。将三维图像体积中的体素“分割”(例如通过对其亮度值进行“阈值处理”或通过另一种方法)以区分岩石基质与空隙空间。接着执行流体流动或其它物理行为如弹性或电导率的直接数值模拟，从中可以导出孔隙率、渗透率(绝对和/或相对)、弹性性质、电性质等。可以应用各种数值方法来求解或逼近模拟适当行为的物理方程。这些方法包括格子玻尔兹曼(lattice-boltzmann)、有限元、有限差分、有限体积数值法等。

附图说明

为了详细描述示例性实施方式，现在将参考附图，所述附图可能未按比例绘制，其中：

图1a是说明用于根据本文公开的原理构造和操作的测试系统的岩石样品来源的实例的示意性层级图；

图1b显示了根据本文公开的原理分析岩石样品的测试系统的方框图；

图1c显示了适用于根据本文公开的原理分析岩石样品的测试系统中的计算装置的方框图；

图2显示了根据本文公开的原理分析在各种饱和度条件下的岩石样品的方法的流程图。

图3a显示了适用于本文公开的各种实施方式的岩石样品的横截面显微照片；

图3b显示了根据本文公开的各种实施方式的图3a的岩石样品的数字表示。

图4显示了根据本文公开的原理在岩石样品的数字表示上执行的欧几里德距离变换的实例；

图5显示了根据本文公开的原理在岩石样品的数字表示上执行的曼哈顿距离变换的实例；

图6显示了根据本文公开的原理在岩石样品的数字表示上执行的最大内切球变换的实例；

图7显示了根据本文公开的原理通过连通组分变换鉴定的岩石样品的数字表示中的连通的实例；

图8显示了根据本文公开的原理在岩石样品的数字表示上执行的排驱过程的流程图；

图9显示了根据本文公开的原理在流体从任何一个或多个侧面侵入的岩石样品的数字表示上执行的排驱过程的应用；

图10a显示了基于根据本文公开的原理从一个侧面侵入的排驱模拟的数字岩石中的流体分布；

图10b显示了基于根据本文公开的原理从六个侧面侵入的排驱模拟的数字岩石中的流体分布；

图11显示了根据本文公开的原理为数字岩石确定的有效气体渗透率的图。

图12显示了根据本文公开的原理通过准静态相对渗透率法产生的数据图；并且

图13显示了通过排驱模拟产生的水分布的实例，该排驱模拟防止侵入水被分离的区域。

注释和命名

在以下讨论和权利要求书中，术语“包括”和“包含”以开放方式使用，因此应该被解释为是指“包括但不限于……”。任何形式的术语“连通”、“接合”、“耦合”、“附着”或描述元件之间的相互作用的任何其它术语的任何使用不旨在将相互作用限于元件之间的直接相互作用，并且还可以包括所述元件之间的间接相互作用。术语“软件”包括能够在处理器上运行的任何可执行代码，而不管用于存储软件的介质如何。因此，存储在存储器(例如非易失性存储器)中并且有时称为“嵌入式固件”的代码被包括在软件的定义内。叙述“基于”旨在指“至少部分基于”。因此，如果x基于y，则x可以基于y和任何数量的其它因子。

详细描述

在附图和下面的描述中，相同的部件通常在整个说明书和附图中用相同的参考数字标记。附图不一定按比例绘制。实施方式的某些特征可以以按比例放大的形式或以某种示意形式显示，并且为了清楚和简明起见，可能未显示常规元件的一些细节。本公开容许不同形式的实施方式。具体实施方式被详细描述并且显示于附图中，其中应理解本公开应被认为是本公开的原理的示例，并且不旨在将本公开限于本文所说明和描述的内容。应充分认识到，下文讨论的实施方式的不同教导和组成部分可以单独使用或以任何合适的组合使用以产生希望的结果。

通常，从直接数值模拟导出的岩石物理性质假设岩石的孔隙空间被一种流体完全饱和。然而，地下岩石的孔隙空间通常由两种或更多种流体填充，接着烃回收过程如注水/气体驱动和烃提取改变了地下孔隙流体分布。岩石的孔隙空间内多种不混溶孔隙流体如水、油和/或气体的存在可显著影响岩石物理性质，如用于表征石油和天然气储层的电导率、有效渗透率和相对渗透率。遗憾的是，用于估计岩石物理性质的常规直接数值模拟没有考虑孔隙空间中多种流体的存在。

本公开的实施方式通过使用数值技术在岩石的孔隙空间内填充两种或更多种流体来使用直接数值模拟确定用两种或更多种不混溶相饱和的岩石的岩石物理性质。用于填充孔隙空间的数值技术基于形态学算法和逼近流体饱和过程，如排驱和渗吸。接着在计算域上操作用于计算岩石物理性质的直接数值模拟，所述计算域明确地将流体的相关物理性质，如适合于每种流体类型的流体粘度、流体密度和/或流体传导率，和/或相关物理性质如接触角指派给与计算模型中的一个或多个流体相单元相邻的固体面。

本文公开的实施方式所利用的用不同流体填充孔隙空间的形态学算法与在地下发生的饱和过程有关。重要的是使饱和岩石的方法反映真实的流体饱和过程，因为这些过程控制孔隙空间内孔隙流体的体积分数和空间分布两者。用不产生物理上合理的流体分布模式的孔隙流体填充岩石孔隙空间的算法不会产生高保真度的模拟的岩石物理性质。例如，仅假设每个孔隙被两种或更多种流体部分占据将产生使人误解的电性质。

图1a以高层级说明了根据本文公开原理的岩石样品的获取和岩石样品的分析。本公开的实施方式可能特别有益于分析来自石油和天然气生产中重要的地下地层的岩石样品。因此，图1a说明了可以根据各种实施方案从中获得待用测试系统102分析的岩石样品104的环境100。在这些说明的实施例中，可以从地面钻探系统106或从海洋(大洋、海、湖等)钻探系统108获得岩石样品104，两种系统都用于提取如烃(石油、天然气等)、水等的资源。作为本领域的基础，石油和天然气生产操作的优化很大程度上受地面钻探系统106或海洋钻探系统108钻入的或过去已钻入的岩层的结构和物理性质影响。

获得岩石样品104的方式以及这些样品的物理形式可以广泛变化。可用于本文公开的实施方式的岩石样品104的实例包括全岩心样品、侧壁岩心样品、露头样品、钻屑和实验室产生的合成岩石样品如砂包和胶结包。

如图1a所示，环境100包括测试系统102，该测试系统102被配置成分析岩石样品104的图像128(图1b)，以便确定对应的地下岩石的物理性质，此类性质包括石油和天然气勘探和生产背景下的岩石物理性质。

图1b以通用方式说明了分析图像128的测试系统102的组成部件。在一般意义上，测试系统102包括成像装置122，该成像装置122用于获得岩石样品104的二维(2d)或三维(3d)图像以及其它表示形式，此类图像和表示形式包括岩石样品104的内部结构的细节。成像装置122的实例是x射线计算机断层摄影(ct)扫描仪，如本领域中已知的，所述x射线计算机断层摄影(ct)扫描仪发射与物体相互作用的x射线辐射124，并测量物体对该x射线辐射124的衰减，以便产生物体内部结构和成分的图像。ct扫描仪122的特定类型、构造或其它属性可以对应于能够产生代表岩石样品104的内部结构的图像的任何类型x射线装置如微ct扫描仪的特定类型、构造或其它属性。成像装置122产生岩石样品104的一个或多个图像128，并将这些图像128转送到计算装置120。

由成像装置122产生的图像128可以是由岩石样品104的多个二维(2d)截面组成或由所述截面产生的三维(3d)数字图像体积(即数字岩石)的形式。在这种情况下，将每个图像体积分成称为体积元素或更通常“体素”的3d规则元素。通常，每个体素是立方体，其在x、y和z方向上具有相等长度的边。另一方面，数字图像体积128本身可以在x、y和z方向上含有不同数量的体素。数字体积内的每个体素具有相关的数值或幅度，所述数值或幅度表示在数字体积所表示的介质的位置处成像样品的相对材料性质。这些数值的范围，通常称为灰度范围，取决于数字体积的类型、值的颗粒度(例如8位或16位值)等。例如，16位数据值使得x射线断层摄影图像体积的体素具有在0至65,536范围内的幅度，颗粒度为1。

如上文所提及，成像装置122将图像128转送到计算装置120，在图1b的实例中，计算装置120可以是任何类型的计算装置，例如台式计算机或工作站、膝上型计算机、服务器计算机、平板计算机等。因此，计算装置120将包括常规计算装置中通常存在的硬件和软件部件。如图1b所示，计算装置120的这些硬件和软件部件包括测试工具130，该测试工具130被配置成分析图像128以确定岩石样品104在一个或多个模拟的流体饱和度条件下的岩石物理性质，所述模拟的流体饱和度条件包括岩层在地下可能遇到的流体饱和度条件。在这方面，测试工具130可以作为软件、硬件或两者的组合实施，包括用于执行本文进一步详细描述的功能和过程的必要且有用的逻辑、指令、例程和算法。在一般意义上，测试工具130被配置成分析岩石样品104的图像体积128，以在代表岩层的地下条件的流体饱和度条件下执行岩石物理性质的直接数值模拟，包括用多种流体改变饱和度。

图1c一般性地说明了根据各种实施方式的测试系统102中的计算装置120的架构。在该示例性架构中，计算装置120包括一个或多个处理器902，其可以具有工业中可用的不同核心配置和时钟频率。用于存储由一个或多个处理器902执行的数据和/或程序指令的计算装置120的存储器资源包括在计算装置120操作期间用作主存储器的一个或多个存储器装置904，和一个或多个存储装置910，例如所述存储装置910作为非易失性固态存储器、磁盘或光盘驱动器或随机存取存储器中的一个或多个来实现。提供一个或多个外围接口906以与对应的外围装置如显示器、键盘、鼠标、触摸板、触摸屏、打印机等耦合。提供可以是以太网(ethernet)适配器、无线收发器、串行网络部件等的形式的网络接口908，以便于通过一个或多个网络在计算装置120之间进行通信，所述网络如以太网、无线以太网、全球移动通信系统(globalsystemformobilecommunications；gsm)、增强型数据速率gsm演进(enhanceddataratesforgsmevolution；edge)、通用移动电信系统(universalmobiletelecommunicationssystem；umts)、全球微波接入互操作性(worldwideinteroperabilityformicrowaveaccess；wimax)、长期演进(longtermevolution；lte)等。在该示例性架构中，处理器902被显示为通过单个总线与部件904、906、908和910耦合；当然，可以将不同的互连架构如多个总线、专用总线等并入计算装置120中。

尽管作为单个计算装置示出，计算装置120可以包括几个计算装置，所述几个计算装置一起协作以提供计算装置的功能。同样地，尽管作为实体装置示出，计算装置120也可以表示抽象计算装置如虚拟机和“云”计算装置。

如图1c的示例性实施方案中所示，计算装置120包括软件程序912，所述软件程序912包括一个或多个操作系统、一个或多个应用程序等。根据实施方式，软件程序912包括对应于测试工具130(图1b)的程序指令，所述程序指令作为独立的应用程序实施，作为程序模块实施，该程序模块是另一个应用或程序的一部分，作为适当的插件或用于访问经由网络接口908与计算装置120联网的远程计算机上的测试工具软件的其它软件部件实施，或者以其它形式和其组合形式实施。

存储与测试工具130的功能对应的软件程序912的可执行指令的程序存储器可以实体地存在于计算装置120内或计算装置120可访问的其它计算资源中，即在存储器装置904和存储装置910的本地存储器资源内，或在服务器或其它网络可访问的存储器资源内，或分布在多个位置之间。在任何情况下，该程序存储器组成存储可执行计算机程序指令的非暂时性计算机可读介质，根据所述计算机程序指令，计算装置120或服务器或与计算装置120通过网络接口908耦合的其它计算机执行本说明书中描述的操作(例如在输入从计算装置120传送的数据之后以交互式应用的形式，用于由与计算装置120耦合的外围装置显示或输出)。对应于与测试工具130相关的软件程序912的计算机可执行软件指令可以最初存储在可移动或其它非易失性计算机可读存储介质(例如dvd盘、闪存等)上，或可以软件包的形式作为关于电磁载波信号的编码信息下载，计算装置120以常规软件安装方式从该软件包安装所述计算机可执行软件指令。预期本领域技术人员能够以适合于每个特定应用的方式容易地实施可用于该实施方式的适用数据、程序指令和其它信息的存储和检索，而无需过度的实验。

组成与测试工具130相关的软件程序912的特定计算机指令可以是一个或多个可执行程序的形式，或源代码或更高级代码的形式，其中从所述代码获得、组装、解释或编译出一个或多个可执行程序。根据执行所希望的操作的方式，可以使用许多计算机语言或协议中的任一种。例如，用于产生根据实施方式的模型的这些计算机指令可以用常规高级语言如python、java、fortran或c++编写成常规的线性计算机程序，或可以被布置成以面向对象的方式执行。还可以将这些指令嵌入在更高级别的应用中。在任何情况下，预期参考本说明书的本领域技术人员能够容易地以对于所希望的安装合适的方式实现实施方式而无需过度的实验。

现在将参考图2以及图1a至图1c来描述测试工具130的特定功能，包括通过软件程序912实施的那些功能，以根据实施方式分析在各种饱和度条件下的岩石样品。

图2显示了用于根据本文公开的原理分析在各种饱和度条件下的岩石样品的方法200的流程图。尽管为方便起见依序描绘，但是所示的至少一些操作可以以不同的顺序执行和/或平行地执行。另外，一些实施方式可以仅执行所示操作中的一些。在一些实施方式中，方法200的至少一些操作以及本文所述的其它操作可以作为存储在计算机可读介质中并由一个或多个处理器902执行的指令实施。

在方框202中，测试系统102获取待分析的岩石样品104，该岩石样品104如来自通过地面钻探系统106或海洋钻探系统108获得的地下岩层或来自其它来源。特定岩石样品104可以从较大体积的地下岩层制备，以例如通过钻凿或切出较大体积的目标岩层的一部分而具有可以由成像装置122(例如ct扫描仪)成像的大小、尺寸和配置。

在方框204中，测试系统102的成像装置122与计算装置120组合产生数字图像体积128，该数字图像体积128代表岩石样品104，包括其内部结构。例如，如果成像装置122是ct扫描仪，则执行岩石样品104的x射线成像(即发射朝向岩石样品104的辐射并测量衰减)以产生2d切片图像的图像体积128或从2d切片图像产生图像体积128。用于方框204中获取和处理岩石样品104的3d数字图像体积128的特定常规技术包括但不限于x射线断层摄影术、x射线显微断层摄影术、x射线纳米断层摄影术、聚焦离子束扫描电子显微术和核磁共振。在一些实施方式中，数字图像体积128可以是计算产生的，而不是通过扫描实体样本产生的。在通过扫描岩石样本产生数字图像体积128的实施方式中，岩石样本可以是天然存在的岩石或人造多孔材料(例如合成岩石)。

图3a说明了岩石样品的3d图像的一个2d切片图像300的实例，其显示了该岩石样品的结构细节的横截面切片，包括固体材料302和孔隙或空隙空间304的特征。此时的图像数据可以是灰度值的形式，所述灰度值代表岩石样品104的成分对x射线辐射的衰减。尽管图3a说明了一个2d切片图像300，但岩石样品104的3d数字图像体积128通常由沿着岩石样品104的一个轴步进的位置处的多个2d切片图像组成，所述多个2d切片图像一起形成岩石样品104的3d图像。将2d切片图像组合成3d数字图像体积128可以根据测试系统102的特定架构，从成像装置122产生的一系列2d切片图像128由成像装置122本身内的计算资源执行，或由计算设备120执行。

在方框206中，测试系统102对岩石样品104的数字图像体积128执行分割或其它图像增强技术，以从图像的灰度值中区分和标记图像体积128的不同组分或相。更具体地，计算装置120执行该分割以便鉴定组分如孔隙空间和矿物组分(例如粘土和石英)。在一些实施方式中，测试工具130被配置成将图像体积128分割成超过两个主要相，所述相代表如孔隙空间、粘土级分、石英级分和其它各种矿物类型的材料成分。

计算装置120可以使用多种类型的分割算法中的任一种。一种分割方法是对图像体积128应用“阈值处理”过程，其中计算装置120在体素幅度范围内选择阈值。具有低于阈值的幅度的那些体素被指派一个表示孔隙空间的特定数值，而具有高于阈值的幅度的那些体素被指派另一个表示基质空间(即固体材料)的数值。在该方法中，阈值处理将灰度图像体积转化为具有两个可能数值之一的体素的分割体积，所述两个可能数值通常选择为0和1。图3b说明了通过阈值处理对3d数字图像体积执行的分割的实例。如所说明的，分割使得可区分岩石样品的结构细节，在该实例中，以浅灰色显示固体材料302，并且以黑色显示孔隙或空隙空间304。还可以应用分割一次或多次以区分灰度图像内的各种特征。如果使用简单的阈值处理，则多个阈值可以区分表现出不同x射线衰减特征的不同材料，如粘土、石英、长石等。

或者，计算装置120可以使用其它分割算法。这种替代算法的实例在本领域中称为大津法(otsu'smethod)，其中基于直方图的阈值处理技术选择阈值以使灰度值的双峰分布的波瓣的组合方差(即“类内方差”)最小化。大津法可以很容易地自动化，并且还可以扩展到多次重复地对图像进行阈值处理以区分其它材料组分，如石英、粘土和长石。可选地或另外地，计算装置120可以使用具有不同复杂度的自动分割算法的其它实例来区分图像体积的不同特征，此类算法包括指示克里金法(indicatorkriging)、收敛主动轮廓法(convergingactivecontours)、分水岭法(watershedding)等。

计算装置120还可以利用其它图像增强技术来增强或改善图像体积128中定义的结构，以进一步区分结构，减少噪声影响等。同样地，虽然计算装置120可以执行分割或其它图像增强技术，但是预期测试系统102的其它部件例如成像装置122本身可以可选地整体或部分地执行图像增强。

因此，分割将数字图像体积中的体素与岩石样品104内对应位置处的特定材料(或孔隙空间，视情况而定)相关联。将一些或所有体素各自用一种或多种材料性质标记，这些材料性质对应于指派给该体素的特定材料成分。此类成分包括孔隙空间、基质材料、粘土级分、单个颗粒、颗粒接触区、矿物类型等。

在方框208中，计算装置120对代表岩石样品104的数字图像体积128的鉴定相应用形态学操作。计算装置120执行的形态学操作可以包括距离变换、内切球变换和连通组分变换中的一个或多个。由计算装置120执行的形态学操作的其它细节如图4-7和本文的相关文本中所提供。

在方框210中，计算装置120应用方框208的形态学操作的结果，以在代表岩石样品104的数字图像体积128内分布多个流体相。多个流体相的分布可以包括用不同比率的润湿流体与非润湿流体来使数字图像体积饱和。使用用两种或更多种不混溶流体相饱和的数字图像体积作为数字直接数值模拟的输入使得可预测多相岩石性质，其中多相岩石性质取决于孔隙空间中两种(或更多种)流体的相对位置。流体分布操作可以包括排驱和渗吸。由计算装置120执行的流体分布操作的其它细节如图8和本文的相关文本中所提供。

在方框212中，计算装置120对饱和的数字图像体积128进行数值分析以确定岩石104的物理性质。数字图像体积128的数值分析可以包括执行如方框210中饱和的数字图像体积128的数字直接数值模拟以分析在模拟的流体饱和度条件下岩石样品104的一种或多种物理性质。在石油和天然气勘探和生产的背景下，可以在方框212中确定目标岩石物理性质，如孔隙度、地层因子、绝对和相对渗透率、电性质(如地层因子、胶结指数、饱和度指数、曲折因子)、毛细管压力性质(如汞毛细管注入)、弹性模量和性质(如体积模量、剪切模量、杨氏模量(young'smodulus)、泊松比(poisson'sratio)、拉梅常数(lameconstant))、核磁共振性质等。这些岩石物理性质可以使用数字图像体积128的适当离散化以及适当的直接数值模拟来估计，该直接数值模拟例如直接数值模拟单相流体流动以计算绝对渗透率。确定这些岩石物理性质中的一些也可能需要使用有限元法、有限差分法、有限体积法、格子玻尔兹曼法或任何种类的其它数值方法进行直接数值模拟。还可以估计由数字图像体积128代表的材料的不同岩石物理性质与孔隙率的关系或这些性质的其它对的关系。

现在返回图2的方框208，计算装置120可以对代表岩石样品104的数字图像体积128的鉴定相应用各种形态学操作中的任一个。应用的形态学操作可以包括距离变换、内切球变换和/或连通组分变换。当应用于目标相(例如孔隙或固体)时，距离变换将用表示到背景相(即非目标相的体素)的距离的值标记目标相中的所有体素。因此，如果目标相是孔隙空间，则距离变换产生每个孔隙体素到最近的固体体素的距离图。可以使用距离图鉴定与数字岩石中的孔隙-颗粒、孔隙-粘土界面相邻的孔隙相体素。润湿膜通常低于ct扫描的分辨率，因此距离图可以用于鉴定孔隙空间中与流体膜所存在的固体表面相邻的区域。使用距离图可以调节膜的厚度。膜变厚的过程被比作渗吸过程，其中润湿相的饱和度增加。可以使用各种方法来实施距离变换。本公开描述了欧几里德距离变换和曼哈顿距离变换(也称为城市-街区(city-block)距离变换)。一些实施方式可以应用各种其它距离变换中的任一种来产生距离图。

图4显示了根据本文公开的原理对岩石样品的数字表示执行的欧几里德距离变换的实例。在图4中，为了清楚起见，欧几里德距离变换以二维进行说明，但在实践中可以在方法200中以三维进行实施。图4显示了待变换的输入图像402、输入图像402的平方值欧几里德距离变换404以及对应于平方值欧几里德距离变换404的实值欧几里德距离变换406。在输入图像402中，孔隙空间像素由值“0”表示，并且固体空间像素由值“1”表示。对于每个孔隙空间像素，欧几里德距离变换鉴定出最近的固体空间像素，并通过毕达哥拉斯定理(pythagoreantheorem)计算出两个像素之间的距离。

在输入图像402的平方值欧几里德距离变换404中，将每个孔隙空间像素值“0”用从孔隙空间像素到最近的固体空间像素的距离的平方值代替。在404中将每个固体空间像素表示为“0”。欧几里德距离变换可以假设正交方向上的相邻体素之间的距离等于一个体素尺寸，这确保了通过取欧几里德距离变换值(表示为实值欧几里德距离变换406)的平方获得的平方值欧几里德距离变换404始终是整数。使用平方值欧几里德距离变换404是有利的，因为整数值允许编程逻辑比处理实数的编程逻辑更精确。虽然欧几里德距离变换已经被描述为从孔隙空间体素到固体空间体素的距离，但是一些实施方式可以应用欧几里德距离变换来确定从每个固体空间体素到最近的孔隙空间体素的距离。

图5显示了根据本文公开的原理对岩石样品的数字表示执行的曼哈顿距离变换的实例。在图5中，为了清楚起见，曼哈顿距离变换以二维进行说明，但在实践中可以在方法200中以三维进行实施。图5显示了待变换的输入图像402、输入图像402的实值曼哈顿距离变换504以及对应于实值曼哈顿距离变换504的平方值曼哈顿距离变换506。在输入图像402中，孔隙空间像素由值“0”表示，并且固体空间像素由值“1”表示。在实值曼哈顿距离变换504和平方值曼哈顿距离变换506中，固体空间像素由值“0”表示。曼哈顿距离变换也可以称为城市-街区距离。曼哈顿距离变换假设正交方向上的相邻体素之间的距离等于一个体素尺寸。指派给像素的值等于一个笛卡尔方向(cartesiandirection)上的最小距离。对于n维图像，指派给孔隙空间像素的距离等于在n个方向中的任一个上从孔隙空间像素到最近的固体空间像素的最小距离。因此，在实值曼哈顿距离变换504中，将每个孔隙空间像素值用一个方向(x或y)上与最近的固体空间像素的距离代替。在平方值曼哈顿距离变换506中，取实值曼哈顿距离变换504的距离值的平方。

通过曼哈顿距离变换而指派给像素的距离值可以大于或等于对应的欧几里德距离变换值，因为在曼哈顿距离变换中，最近邻域搜索被限于n维图像的n个方向上。然而，因为曼哈顿距离变换将最近邻域搜索限于n维图像的n个方向上，所以曼哈顿距离变换可以比欧几里德距离变换在计算上更有效。虽然曼哈顿距离变换已经被描述为从孔隙空间体素到固体空间体素的距离，但是一些实施方式可以应用曼哈顿距离变换来确定从每个固体空间体素到最近的孔隙空间体素的距离。

最大内切球变换是另一种形态学操作，其可以应用于代表岩石样品104的数字图像体积128的鉴定相。最大内切球变换提供了关于多孔材料结构的信息，可以用于在数字岩石中分布两种或更多种不混溶流体。应用于目标相(例如孔隙或固体)的最大内切球变换向目标相的每个体素指派等于覆盖该点并且被包含在目标相内的最大球体的值。可以从距离变换如本文所述的欧几里德距离变换或曼哈顿距离变换估计给定体素的球体大小。当使用欧几里德距离变换时，球体将被完全包含在孔隙空间内。使用替代的距离变换，例如曼哈顿距离变换或估计的欧几里德距离变换，可以产生mis变换，其中最大覆盖球体与背景相(即非目标相)重叠。

图6显示了根据本文公开的原理对岩石样品的数字表示执行的最大内切球变换的实例。在图6中，为了清楚起见，最大内切球变换以二维进行说明，但在实践中可以在方法200中以三维进行实施。图6显示了输入图像402的平方值欧几里德距离变换404和对应于平方值欧几里德距离变换404的平方值最大内切球变换606。在平方值欧几里德距离变换404和平方值最大内切球变换606中，固体空间像素由值“0”表示。

在图6中，对于每个孔隙空间像素，半径等于像素的平方欧几里德距离变换值的平方根的圆以像素中心为中心。像素的最大内切球变换值被认为是产生覆盖像素的最大圆的平方欧几里德距离变换值。在三维情况下，采用球体而不是圆。可以基于用于计算不同相的体素之间的距离的距离图来改变最大内切球变换。

当应用于数字岩石的孔隙空间时，允许最大内切球体变换的覆盖球体与背景相重叠(如图6所示)，当假设背景相为固体或微孔相时，在背景相附近产生流体-流体界面，该界面可以通过角度来表征。通过改变与背景相的重叠量，可以改变表征背景相附近的流体-流体界面的角度。

最大内切球变换可以用于通过称为阈值处理的图像处理技术产生非润湿流体相的分布。使用给定阈值，基于最大内切球变换中含有的值，可以将大于或等于该阈值的所有像素/体素值指派给非润湿相。类似地，可以将所有小于阈值并且存在于最初被认为是最大内切球变换的目标相的相中的值指派给润湿相。通过选择不同的阈值，可以产生数字岩石内两个或更多个不混溶流体相的不同分布。通常将相的饱和度定义为百分比。通常将该百分比定义为相的体积除以可供该相占据的体积(例如流体占据的孔隙空间的体积除以孔隙空间的总体积)。例如，水饱和度通过占据孔隙空间的水的体积除以孔隙空间的总体积来确定。

连通组分变换是另一种形态学操作，其可以应用于代表岩石样品104的数字图像体积128的鉴定相。连通组分变换鉴定连通的像素/体素组。例如，可以应用连通组分变换来鉴定连通的孔隙空间体素。连通组分变换的结果可以用于将两种或更多种不混溶流体的分布限于连通的孔隙空间的区域或限于不连通的孔隙空间的区域。

连通组分变换(也称为群集鉴定)根据三种可能的连通性规则之一鉴定目标相中连通的体素组，所述规则可以由连通数6、18或26定义。连通数鉴定包括在每个点处的邻域搜索中的体素数，并且可以通过将每个体素视为具有8个节点的立方元素来可视化。图7显示了根据本文公开的原理对岩石样品的数字表示执行的连通组分变换中应用的三个不同连通性规则的实例。702中说明基于面的连通性(选项6)。在基于面的连通性中，两个体素被鉴定为共享一个共用面或最少四个节点(在702中以虚线表示)。704中说明基于边的连通性(选项18)。在基于边的连通性中，两个体素被鉴定为共享一个共用边或最少两个节点(在704中以虚线表示)。706中说明基于点的连通性(选项26)。在基于点的连通性中，两个体素被鉴定为共享一个共用点或单个节点(在706中以放大点表示)。

连通组分分析与侵入模拟结果的组合可以用于预测初始润湿相饱和度和残余非润湿相饱和度。侵入是流体进入可渗透岩石的过程。在侵入过程中，两种或更多种不混溶流体将存在于多孔材料的孔隙空间内。对于数字岩石应用，侵入过程对两种不混溶流体的位置的变化和孔隙空间内这些流体的量的变化进行数值模拟。孔隙空间内每种流体量的变化被认为是饱和度变化。

现在转到图2的方框210，应用方框208的形态学操作的结果以在数字岩石的孔隙空间中分布多个流体相(例如润湿相和非润湿相)。流体进入孔隙空间通常被称为“侵入”。在石油和天然气工业中，侵入过程通常分为排驱过程或渗吸过程。侵入过程是排驱还是渗吸通常由系统的宏观反应以及系统表现出水湿特征或油湿特征的趋势决定。然而，将响应归为排驱或渗吸是系统对岩石的侵入流体和固相之间的接触角差异的响应。侵入流体是饱和度增加的流体。侵入流体的接触角由与固体相邻的流体-流体界面的形状限定。如果接触角小于90°，则侵入流体被认为是润湿流体。如果接触角大于90°，则侵入流体被认为是非润湿流体。在排驱过程中，侵入流体被认为是非润湿流体。在渗吸过程中，侵入流体被认为是润湿流体。

在排驱过程中，岩石的孔隙空间中的润湿流体被非润湿流体从孔隙空间中驱替。与石油和天然气工业相关的排驱过程的实例包括：油迁移到储层中，在钻探过程期间非润湿钻探液从钻孔侵入储层中，以及在生产期间将气体注入储层中。可以使用三种操作的组合来执行排驱过程：侵蚀、连通组分分析和扩张。在驱替模拟开始时，用润湿流体如水使数字岩石的孔隙空间饱和。非润湿流体在零毛细管压力下与岩石的一个或多个侧面(即入口)连通。在侵蚀操作中，增加毛细管压力以使非润湿流体侵入孔隙空间，并且根据杨-拉普拉斯方程(young-laplaceequation)确定与毛细管压力对应的球体的半径。在给定的毛细管压力下，确定孔隙空间中半径为r的球形结构元素r的中心的位置。之后，在连通组分分析中，鉴定与入口连通的侵蚀的孔隙空间。最后，使用在侵蚀中鉴定的半径为r的球体，对连通组分分析中鉴定的孔隙空间应用扩张。

在常规的排驱过程中，在每个施加的毛细管压力下迭代地执行侵蚀、连通组分分析和扩张。本文公开的系统和方法的实施方式应用了一种改善的排驱过程，其减少模拟数字岩石中的排驱过程所需的时间和计算资源。通过减少模拟排驱所需的时间和计算资源，本公开的实施方式使得可比常规模拟过程更快速地且更低成本地确定岩石的岩石物理性质。本文公开的排驱过程以单程方式执行侵蚀和扩张以显著减少模拟排驱所需的计算资源。

图8显示了根据本文公开的原理对岩石样品的数字表示执行的排驱过程的流程图。尽管为方便起见依序描绘，但是所示的至少一些操作可以以不同的顺序执行和/或平行地执行。另外，一些实施方式可以仅执行所示操作中的一些。在一些实施方式中，方法800的至少一些操作以及本文所述的其它操作可以作为存储在计算机可读介质中并由一个或多个处理器902执行的指令实施。

在方框802中，对数字岩石孔隙空间的每个体素计算欧几里德距离变换。欧几里德距离变换用与最近的固体空间体素的欧几里德距离值标记每个孔隙体素。

在方框804中，对方框802中产生的欧几里德距离变换值应用阈值处理。保留所有大于或等于阈值的距离值。可以改变阈值以改变通过排驱产生的非润湿流体的饱和度。在排驱过程800中，对欧几里德距离变换应用值r的阈值产生的结果与用半径为r的球形结构元素常规侵蚀孔隙空间产生的结果相同，但具有低得多的计算复杂度。

在方框806中，分析孔隙体素与入口的连通性。连通性分析可以包括迭代地分析对应于方框804的阈值处理中保留的欧几里德距离值的体素位置的连通性。将每次迭代的连通性分析的结果添加到跟踪与入口连通的体素集的演变的数据集中。入口可以是数字岩石的6个侧面的任何组合，其中非润湿流体通过该入口侵入数字岩石。

在方框808中，对通过方框806的连通性分析产生的数据集应用扩张。对所述数据集的元素应用扩张以便将元素添加到所述数据集。也就是说，扩张在入口处开始并且在孔隙空间中继续进行。与常规侵蚀不同，方法800的扩张过程以单程方式而不是迭代地执行。

方法800的排驱过程通过用应用于欧几里德距离变换的阈值代替常规排驱过程中耗时的侵蚀步骤而提供了改善的计算性能。使用欧几里德距离变换改善了模拟的毛细管压力曲线的保真度。这是因为对于任何给定的孔隙体素，该体素处的欧几里德距离变换的值提供以该体素为中心的最大球体的值。因此，欧几里德距离变换中离散值的数量依次限定了给定图像的孔隙空间可以支持的不同球体大小和毛细管压力的最大数量。因为方法800的排驱过程自动探测欧几里德距离变换中含有的每个值，所以方法800确保能优化数字岩石应用的模拟毛细管压力曲线的分辨率。

图9显示了数字岩石图像900的六个侧面的示意图。数字岩石图像900可以对应于图像体积128。方法800的排驱模拟允许非润湿相从数字岩石900的面的任意组合侵入。在电导率模拟中，模拟了通过数字岩石900的电流。例如，可以通过在z方向上施加电压降来感应出通过数字岩石900的电流。可以通过相对于出口平面(称为zmax)，在入口平面(称为zmin)处施加较高的电压来建立电压降。通过从数字岩石900的不同侧面侵入来执行排驱模拟可以产生沿着电流方向显著不同的饱和度分布，并且这又产生不同的电导率。

基于排驱毛细管压力(pc)的分布使用排驱毛细管压力模拟的结果来确定两个不混溶相的分布，其中一个优先是润湿相。基于pc结果的分布模式需要侵入域的非润湿相保持与入口连通。使用方法800在排驱模拟器中设计的灵活性使得非润湿相可从6个侧面的任意组合侵入孔隙空间。

图10a和10b显示了通过砂岩体积的3d图像的2d切片，以及基于使用方法800进行离散排驱模拟的两种流体分布图。两个图像中的饱和度分布是在相同的毛细管压力下。图10a显示了从zmin面(参见图9)的单侧侵入，并且图10b显示了6侧侵入。在图10a和10b中，白色表示固体空间，黑色表示被润湿流体占据的孔隙空间，并且灰色表示被非润湿流体占据的孔隙空间。

基于排驱模拟的结果，计算装置120可以将润湿性分布指派给数字岩石。润湿性分布为与孔隙空间体素或每个固体空间体素相邻的每个固体空间体素设定流体接触角值。在油-水储层系统中，储层可以被称为水湿、混湿或油湿的。确定给定表面是水湿还是油湿通常是基于油-水界面与固体表面的接触角。在水湿系统的情况下，接触角是小于90°的值(水铺展在表面上)。对于油湿系统，水倾向于坐立在表面上，使得接触角大于90°。指派给体素的接触角可以基于排驱过程800的结果。例如，在根据方法800完成排驱模拟时，可以确定与每个固体空间体素接触的流体。如果流体是润湿流体，则可以将小于90°的接触角指派给体素。如果流体是非润湿流体，则可以将大于90°的接触角指派给体素。在一些实施方式中，例如在分割期间区分颗粒材料类型的情况下，还可以基于固体空间的每个颗粒的确定材料来指派接触角。通常假设，由于油中存在的有机物质的沉积，岩石表面变成油湿的，而与水接触的表面仍为水湿表面。当假设油向储层的迁移是排驱过程时，建立良好的油-水相分布为在数字岩石中分布润湿性提供了物理方法。

渗吸是可以在图2的方框210中模拟的另一种侵入过程。在渗吸过程期间，用润湿流体从数字岩石的孔隙空间中驱替非润湿流体。渗吸过程的模拟容易从用于模拟排驱的框架修改而来。与侵入相优先占据较大孔隙的排驱模拟不同，在渗吸过程期间，侵入相优先占据数字岩石的临近壁区和较小孔隙。与石油和天然气工业相关的渗吸过程的实例是：来自储层下方的含水层的水可以进入储层的含烃区的情形；水注入水湿储层中；水注入烃区与储层下方含水层之间的接触区附近的区域内；以及在钻探过程期间钻探液侵入储层中。

现在返回图2的方框212，给定通过应用本文公开的排驱和/或渗吸过程产生的流体饱和度，计算装置120对数字图像体积128进行数值分析以确定岩石104的物理性质。在石油和天然气勘探和生产的背景下，各种岩石物理性质可能为人们所专注，并且可以使用数字图像体积128和适当的直接数值模拟来估计此类性质。模拟的性质可以包括有效气体渗透率、准静态相对渗透率、核磁共振和电阻率指数。

有效的气体渗透率计算使用排驱模拟的结果(以及在数字岩石图像的孔隙相中分布两种流体以产生离散饱和态(独特的气体饱和度sg)的方法。接着将每个饱和态用作单相流动直接数值模拟的输入。将给定气体饱和度下的有效气体渗透率定义为：

其中：

kabs表示岩石的绝对渗透率(流体流过岩石的容易性的度量)；

kg(sg)是给定气体饱和度下岩石对气体的渗透率；并且

kg,eff是两个值的比率。

通过在多个气体饱和度下执行直接数值模拟来预测数字岩石的有效气体渗透率随水饱和度的变化。有效气体渗透率应用的实例包括：评价气体产生速率随水饱和度的变化；和鉴定临界水饱和度，即气体产生可忽略不计(无益)时的水饱和度。有效的气体渗透率可以基于排驱模拟的输出或在数字岩石的孔隙空间中分布两种流体的任何其它方法来确定。例如，可以基于如本文所述的毛细管排驱法或最大内接球法来确定有效的气体渗透率。

图11显示了有效气体渗透率计算的结果。在图11中，相对气体渗透率数据显示了基于几个输入图像的排驱模拟对饱和度模式的响应。在产生气体的储层的情况下，需要了解气体的有效渗透率随饱和度是如何变化的。该图突出显示了不同的岩石可以如何在约70-25％之间的气体饱和度下停止流动的气体。

特定流体的相对渗透率是在给定饱和度下该流体的有效渗透率与渗透率的某个基值(绝对值)的比率。给定岩石的相对渗透率通常通过稳态方法或非稳态方法通过实验确定。在稳态和非稳态方法两者中，将流体引入某个初始饱和度条件下的岩石的孔隙空间中。可以如下计算准静态相对渗透率：

其中

kabs表示如上所述的岩石的绝对渗透率；

k(o.w)s(o.w)是在给定的油/水饱和度下对油或水的渗透率；并且

kr，即相对渗透率，是油或水的渗透率与岩石的绝对渗透率的比率。

通过分别对油和水执行单相直接数值模拟，实施方式可以确定在给定饱和度下各自的相对渗透率。相对渗透率可以用于例如评估油和水的生产速率随饱和度的变化，或鉴定初始生产速率(在初始水饱和度下的有效油渗透率)。

图12显示了基于排驱模拟(pc_1f)和用于定位数字岩石内的润湿和非润湿相的最大内切球(mis)的准静态相对渗透率计算的结果。通过假设在流体-流体(油-水)界面处没有动量传递，可以将饱和度模式用于单相流动模拟器中。在这种情况下，将流体-流体界面视为具有与流体-固体界面相同的边界条件，即无滑移(零速度)边界条件。通过在单相流动模拟器中分别使用润湿和非润湿相分布，实施方式可以产生准静态相对渗透率曲线。准静态渗透率计算可以被认为是逼近或筛选相对渗透率模拟器。虽然直接数值模拟技术可用于模拟两相流动，但此类技术与本文公开的准静态方法相比在计算方面非常昂贵。

系统102提供从不同的初始饱和度开始多相流动模拟的灵活性。此种灵活性是重要的，因为任何给定的储层将在矿田的生产寿命期间表现出一系列饱和态。在相对渗透率模拟期间，通常通过排驱过程建立初始饱和度，其中非润湿相驱替最初占据孔隙空间的水相。一旦初始饱和度建立之后，通常遵循非稳态或稳态程序进行模拟以确定相对渗透率。图13显示了具有通过排驱模拟设定的初始水分布的数字岩石。在模拟期间，方法800的实施方式可以添加标准来防止侵入孔隙空间中水变得分离或与数字岩石的出口不相连的区域。这些区域由图13的亮区鉴定，并且倾向于存在于颗粒-颗粒接触区附近和孔隙空间的较小区域中。该过程可以用于代表建立初始水饱和度的实验室程序，该初始水饱和度有时称为残余或原生水饱和度。

岩石的电阻率指数(ri)是岩石的真实电阻率(例如在给定饱和度下)与填充有水的相同岩石的电阻率的比率。因此，可以将ri定义为：

其中：

rt表示给定饱和度下的电阻率；并且

ro表示100％水饱和度下的电阻率。

使用幂律拟合进行ri数据图相对于水饱和度的最佳拟合给出了给定样品的饱和度指数(n)，其中：

ri＝sw^-n。

将饱和度指数n与地层水电阻率(rw)、真实电阻率(rt)和孔隙率(φ)的对数推导估计值相结合，以通过以下等式估计储层的水饱和度，这是地层评价的关键函数：

通过提供用两种或更多种不混溶流体饱和的数字岩石作为电导率模拟器的输入来模拟电阻率指数。系统102可以包括可基于不同数值方法使用的多种不同电导率模拟器中的任一种。两种常见电导率模拟器利用有限差分或有限元法来求出电流。系统102可以通过提供用两种或更多种不混溶流体饱和的数字岩石作为电导率模拟器的输入来模拟电阻率指数。相对于被视为导电的盐水或微孔相，可以将油(或气体)和固相有效地视为是不导电的。

孔隙空间中流体的分布对电阻率指数有显著影响，又对饱和度指数(n)有显著影响。系统102的实施方式提供了使数字岩石饱和的灵活性，这是重要的，因为了解特定储层的特征性分布对于地层评价是重要的。确定饱和度指数的适当值特别有用，因为电阻率测井通常使用饱和度指数将电测井工具采集的数据转化为水饱和度。了解含烃岩层中水的位置和量对于优化给定储层的开发是重要的。

以上讨论旨在说明本公开的各种原理和实施方式。虽然已经显示并描述了某些实施方式，但是本领域技术人员可以在不背离本公开的精神和教导的情况下对其进行修改。本文所述的实施方式仅是示例性的，而并非是限制性的。因此，保护范围不受上述描述的限制，而是仅受随后的权利要求书限制，其范围包括权利要求书主题的所有等效物。