CN111426616B - 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质 - Google Patents

碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质 Download PDF

Info

Publication number
CN111426616B
CN111426616B CN202010104071.1A CN202010104071A CN111426616B CN 111426616 B CN111426616 B CN 111426616B CN 202010104071 A CN202010104071 A CN 202010104071A CN 111426616 B CN111426616 B CN 111426616B
Authority
CN
China
Prior art keywords
image
rock
digital core
digital
gradient
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010104071.1A
Other languages
English (en)
Other versions
CN111426616A (zh
Inventor
赵建国
潘建国
李劲松
孙朗秋
胡洋铭
刘欣泽
欧阳芳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China National Petroleum Corp Science And Technology Research Institute Co ltd
Northwest Branch Of Petrochina Co Ltd Exploration And Development Research Institute
China University of Petroleum Beijing
Original Assignee
China National Petroleum Corp Science And Technology Research Institute Co ltd
Northwest Branch Of Petrochina Co Ltd Exploration And Development Research Institute
China University of Petroleum Beijing
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China National Petroleum Corp Science And Technology Research Institute Co ltd, Northwest Branch Of Petrochina Co Ltd Exploration And Development Research Institute, China University of Petroleum Beijing filed Critical China National Petroleum Corp Science And Technology Research Institute Co ltd
Priority to CN202010104071.1A priority Critical patent/CN111426616B/zh
Publication of CN111426616A publication Critical patent/CN111426616A/zh
Application granted granted Critical
Publication of CN111426616B publication Critical patent/CN111426616B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume, or surface-area of porous materials
    • G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
    • G01N15/088—Investigating volume, surface area, size or distribution of pores; Porosimetry
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume, or surface-area of porous materials
    • G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume, or surface-area of porous materials
    • G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
    • G01N2015/0846—Investigating permeability, pore-volume, or surface area of porous materials by use of radiation, e.g. transmitted or reflected light
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/03—Investigating materials by wave or particle radiation by transmission
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/60—Specific applications or type of materials
    • G01N2223/616—Specific applications or type of materials earth materials
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/60—Specific applications or type of materials
    • G01N2223/649—Specific applications or type of materials porosity

Abstract

本说明书实施例提供了一种碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质，该方法包括：获取指定岩样的岩石物理参数；基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。本说明书实施例可以获得更准确的碳酸盐岩弹性性质与孔隙结构。

Description

碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质
技术领域
本说明书涉及油气勘探开发技术领域，尤其是涉及一种碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质。
背景技术
碳酸盐岩储层囊括了全球60％的石油储量，且开采潜力巨大。相较于常规砂岩储层，碳酸盐岩储层的特点是后期的成岩改造作用会导致其次生孔隙非常发育，使其孔隙结构极其复杂。碳酸盐岩中常见的孔隙结构类型有铸模孔、粒内溶孔、粒间溶孔、晶间孔、裂缝等，这些复杂多变的孔隙结构会对碳酸盐岩的弹性性质产生显著的影响，其对弹性性质的影响效果甚至比孔隙度还大。比如，对于孔隙度均为13％的碳酸盐岩储层，由于孔隙结构类型的不一致，会出现相同孔隙度下不同孔隙结构类型的碳酸盐岩储层其纵波速度会因孔隙结构不同而有约一千米每秒的速度差异。因此，针对于碳酸盐岩储层，我们不但需要研究碳酸盐岩的弹性性质与孔隙度之间的关系，而且需要更加精确地获取碳酸盐岩储层孔隙结构类型表征与分类下的弹性性质与孔隙度二者之间的关系。反之，如果不能相对精确地考虑孔隙结构类型对碳酸盐岩储层弹性性质的影响因素，碳酸盐岩储层反演误差将会非常大。
为达成上述目标，现有技术以及现有技术的解决方案如下所述：
1)、常规测井：利用测井仪器在储层条件下获得碳酸盐岩储层的弹性性质，包括纵波速度、横波速度、密度等，同时利用孔隙度测井技术获取碳酸盐岩储层孔隙度属性，在此基础上，进行弹性性质与孔隙度交会分析，通常获得一个比较粗略的线性拟合关系，用作后续储层反演与预测分析。该种方案通过测井工具在储层原位地层一定深度范围内进行弹性性质以及孔隙度等的测量，因此可以获得大量的数据点，特别是孔隙度分布范围比较宽，有利于后续实施弹性性质与孔隙度的交会分析，比如“纵波速度-孔隙度”规律分析。然而，一般情况下研究人员利用弹性性质与孔隙度的交会分析只获得一个比较粗略的线性关系，用作后续储层反演研究，此时，对碳酸盐岩孔隙结构的影响不做或无法进行考虑。即便在一些情形下有成像测井工具的参与，也只能获得比较定性的储层孔隙结构类型分类，此时只能实施相对定性的孔隙结构类型分类下的弹性性质与孔隙度关系研究。
2)、取芯及实验室测量。利用碳酸盐岩储层钻井取芯，在实验室利用常规技术获得碳酸盐岩样品的纵波速度、横波速度、密度、以及孔隙度等属性，从而获得弹性性质(如纵波速度)与孔隙度交会图，通常也获得一个比较粗略的线性拟合关系，用作后续储层反演与预测分析。取芯及实验室测量方案中存在的问题在于：首先，在该方案中主要的问题在于碳酸盐岩储层取芯比碎屑岩储层更加困难，因此一般来说可用于实验室测量的碳酸盐岩样品数量有限，对后续研究孔隙度及孔隙结构类型对弹性性质的影响造成制约。其次，由于取芯有限，因此，在很多情况下供测量的碳酸盐岩样品其孔隙度分布范围很窄，这也造成后续规律分析的困难。再次，实验室虽然可以对铸体薄片进行一定的孔隙结构分析，但是这种相对定性的孔隙结构分析，也只能实施相对定性的孔隙结构类型分类下的弹性性质与孔隙度关系研究。
综上所述，现有技术只能相对定性的获取碳酸盐岩储层的弹性性质及孔隙度属性。
发明内容
本说明书实施例的目的在于提供一种碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质，以获取更准确的碳酸盐岩弹性性质与孔隙结构。
为达到上述目的，一方面，本说明书实施例提供了一种碳酸盐岩弹性性质与孔隙结构获取方法，包括：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
另一方面，本说明书实施例还提供了一种碳酸盐岩弹性性质与孔隙结构获取装置，包括：
参数获取模块，用于获取指定岩样的岩石物理参数；
数字岩芯构建模块，用于基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；
岩芯网格化模块，用于网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
弹性性质模拟模块，用于对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
孔隙结构确定模块，用于根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
另一方面，本说明书实施例还提供了另一种碳酸盐岩弹性性质与孔隙结构获取装置，包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时执行如下步骤：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
另一方面，本说明书实施例还提供了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
由以上本说明书实施例提供的技术方案可见，本说明书实施例可以获得网格化的数字岩芯体中，每个子网格对应的数字岩芯体的弹性性质和孔隙结构，从而实现了对碳酸盐岩弹性性质和孔隙结构的精确获取，因而提高了获取碳酸盐岩弹性性质和孔隙结构的准确性，进而可为下一步实施更加精确的储层预测以及储层分类提供了有利基础。
附图说明
为了更清楚地说明本说明书实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本说明书中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本说明书一些实施例中碳酸盐岩弹性性质与孔隙结构获取方法的流程图；
图2为本说明书一实施例中碳酸盐岩岩样的CT扫描横向切片示意图；
图3为本说明书一实施例中碳酸盐岩岩样的纵、横波速度随围压变化的关系示意图；
图4为本说明书一实施例中碳酸盐岩岩样的弹性模量随围压变化的关系示意图；
图5为本说明书一实施例中线性拉伸法的灰阶重新分布结果；
图6a为本说明书一实施例中原始CT扫描图像的取样示意图(灰色矩形框为取样部分)；
图6b为从图6a所示的原始CT扫描图像中取出的CT扫描图像切块的示意图；
图6c为对图6b所示的CT扫描图像切块进行对比度增强处理后的结果示意图；
图7a为本说明书一实施例中CT扫描图像在进行对比度增强处理后的结果示意图；
图7b对图7a所示的CT扫描图像进行各向异性滤波处理后的结果示意图；
图8a为本说明书一实施例中CT扫描图像在进行各向异性滤波处理后的结果示意图；
图8b对图8a所示的CT扫描图像进行边缘增强处理后的结果示意图；
图9a为本说明书一实施例中预处理后的CT扫描图像的示意图；
图9b对图9a所示的CT扫描图像进行图像二值化处理后的结果示意图；
图10a为本说明书一实施例中3-1号碳酸盐岩岩样的数字岩芯示意图；
图10b为图10a所示的数字岩芯的孔隙结构示意图；
图11为本说明书一实施例中数字岩芯的网格化示意图；
图12为本说明书一实施例中获得的碳酸盐岩岩样在分类下的弹性性质与孔隙度关系示意图；
图13为本说明书一些实施例中碳酸盐岩弹性性质与孔隙结构获取装置的结构框图；
图14为本说明书另一些实施例中碳酸盐岩弹性性质与孔隙结构获取装置的结构框图。
具体实施方式
为了使本技术领域的人员更好地理解本说明书中的技术方案，下面将结合本说明书实施例中的附图，对本说明书实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明书一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本说明书保护的范围。
参考图1所示，本说明书一些实施例的碳酸盐岩弹性性质与孔隙结构获取方法，可以如下步骤包括：
S101、获取指定岩样的岩石物理参数。
在本说明书的实施例中，所述岩石物理参数例如可以包括孔隙度、成岩矿物组成、孔隙类型、密度及指定频段下随围压变化的纵、横波速度等。在一实施例中，所述的获取例如可以接收外部输入的指定岩样的岩石物理参数。在另一实施例中，所述的获取例如还以是从指定存储路径读取记录有上述指定岩样的岩石物理参数。
在本说明书一实施例中，上述指定岩样的岩石物理参数可以预先通过以下方式得到：
(1)、可以选取储层代表性碳酸盐岩样若干，一般选取至少5个岩样，岩样孔隙类型粗略地覆盖裂缝型、孔洞型、孔洞-裂缝型，并可以对岩样进行洗盐、洗油、以及烘干干燥等预处理。
(2)、可以对选取的代表性碳酸盐岩岩样进行氦气测量孔隙度，或者可以通过核磁共振测量孔隙度，获得如上代表性岩样比较精确的孔隙度，以用作后续标定。
(3)、可以对选取的代表性碳酸盐岩岩样进行X射线衍射(X-Ray Diffraction，简称XRD)矿物分析，以获取成岩矿物，用作后续弹性性质模拟的矿物骨架模量参数输入。
(4)、可以对选取的代表性碳酸盐岩岩样进行镜下薄片分析，以获取定性半定量的孔隙类型，用作后续标定。
(5)、可以对选取的代表性碳酸盐岩岩样进行密度测量及高频超声频段下随围压变化的纵横波速度测量，以用作后续标定。
其中，根据需要，以上步骤(2)～(5)可按实际情况在实验室测量时微调。
S102、基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯。
在本说明书一实施例中，所述基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，可以包括如下步骤：
(1)、获取岩样的CT扫描图像。其中，CT扫描图像可以是预先对选取的代表性碳酸盐岩岩样进行CT扫描成像获得原始图像剖面。
(2)、将所述CT扫描图像进行图像二值化处理，并根据获得的二值化图像构建第一数字岩芯。
在本说明书一实施例中，所述将CT扫描图像进行图像二值化处理可以包括如下步骤：
(21)、对比度增强
首先可以对CT扫描的原始图像进行剪裁，一般选取边长为1000体素的正方体。剪裁后的CT图像灰阶主要分布在非常低值区域，从图像的直观观测上显示为图像信息都比较暗，这样窄分布的灰阶直方图分布会给后续的图像分割造成困难，因此可以对图像进行对比度增强。对比度增强操作通过灰度变换将原本窄分布的灰阶直方图变换成宽分布的直方图。选用灰度变换公式(1)实现上述功能：
Figure BDA0002387885850000061
其中，f(x,y)表示原始二维图像上坐标(x,y)处像素点的灰度值，g(x,y)表示同一点灰度变换后的值，a、b、c、d均为常数。
(22)、滤波去噪
在本说明书一些实施例中，可以选用中值滤波、均值滤波或高斯滤波等全局滤波类算法对图像进行去噪处理。全局滤波类算法是统一对待图像上每一个像素点，容易对裂缝纹理等信息平滑掉。而碳酸盐岩相比于其他类别的岩石类型具有更加丰富的孔隙结构类型，特别是该类岩石中遍布裂缝、纹理、孔洞等复杂的结构，而这种复杂孔隙结构是如何影响碳酸盐岩弹性性质，是碳酸盐岩地球物理勘探中的重要研究内容之一。因此，获取既能合理去除图像中噪声，又能有效保留裂缝、纹理、孔洞等孔隙结构信息的滤波去噪算法，是碳酸盐岩数字岩芯图像处理中的重要任务。在本说明书一实施例中，可以使用各项异性滤波算法实现上述期望功能，具体可以如公式(2)～(4)所示：
Figure BDA0002387885850000071
Figure BDA0002387885850000072
Figure BDA0002387885850000073
上述公式(2)中，It(s)表示原始图像的灰度函数，It+1(s)表示滤波后图像的灰度函数，s则是像素点在图像上的空间位置坐标，λ为控制扩散强度的常数，λ值越大滤波后图像越平滑，gK(x)为扩散函数，其中x为图像梯度值，K为判断噪声与微观结构信号的图像梯度阈值参数，梯度低于K值判断其为噪音信号，梯度高于K值就判断其为有用的微观结构信号，
Figure BDA0002387885850000074
为像素点与周围各方向相邻像素的差值，其表达式见公式(4)，ηs＝{N,S,E,W}代表像素点s周围北、南、东、西等四个方向的符号。K设置为等于迭代公式(2)中每次迭代中的梯度直方图积分(累积和)90％值所对应的梯度。
(23)、边缘增强
对数字岩芯进行图像边缘增强可以使图像边缘更清晰。边缘增强主要困难在于锐化滤波器很容易放大图像中所有噪声。在本说明书一实施例中，可以采用边缘检测法来进行边缘增强。由于固体颗粒与孔隙间的边缘一般对应于高梯度区域，因此可以使用梯度来提取边缘信号。通过图像梯度的累加求和函数，可以确定梯度累加占总值90％的梯度点，高于这个梯度点，就可以认为是高梯度指示相边界。如此，可以很完整地将固体颗粒与孔隙间的边缘区域给识别出来了，在此基础上进行图像分割将降低分割不确定性，从而获得更好的图像分割效果。
(24)、图像分割
经过上述(21)～(23)处理步骤后，可以对碳酸盐岩数字岩芯图像进行图像二值化处理，从而可以建立初始的数字岩芯，以便于后续实施弹性性质模拟。
在本说明书实施例中，将孔隙结构从CT图像中提取出来，识别每张图像中的孔隙像素，这个识别过程在图像处理中叫做图像分割。图像分割技术一般分为几何形状分割法和聚类分割法。前者是分析灰度直方图几何形状来确定分割阈值的方法，而后者是通过图像中两组数据间的统计属性差异来确定最佳阈值的方法。聚类分割法主要包括以下三个步骤：
(241)、由于图像中每个像素点都对应一个灰度值，聚类算法则将图像中所有像素点的灰度值视作一个数据集合。
(242)、由于CT图像是灰度图像，其灰度分布一般有8bit(0～255)或16bit(0～65535)这两种，以8bit为例，阈值从0取到255共256个值，每个阈值都将数据集合分为两类。
(243)、计算每个阈值分类后两类数据的类间方差与类内方差之比，比值最大点所对应的阈值便是聚类分割法法找到的阈值。与普通阈值分割思路类似，将分出的两类数据的灰度值进行重新赋值，一般是分别赋予最大差距值，这个过程称为图像二值化处理。聚类分割算法中由于Otsu算法基于灰度直方图上单一阈值的特性，尤其适合单矿物岩样。而碳酸盐岩矿物成分比较单一，使用Otsu算法对CT图像进行分割，获得了比较好的效果。
(3)、优化所述第一数字岩芯，以使其孔隙度与所述岩石物理参数中的孔隙度相一致。
基于如上步骤获取的二值化数字岩芯数据体，可以获取孔隙相对于总体积的占比，即基于数字岩芯确定的孔隙度，可以对于所有选取的碳酸盐岩岩样，比较基于数字岩芯获取的孔隙度与此前利用氦气与NMR测量孔隙度，从而实现孔隙度标定工作。如果基于数字岩芯估计的孔隙度与实测度超出误差允许，则可以重新返回二值化图像处理流程，对其中的各个子步骤，即上述步骤(21)～(24)进行参数优选及微调，直至基于数字岩芯估计孔隙度与实测孔隙度相一致，从而可以获得优化后的数字岩芯。此外，优化后的数字岩芯还可以与上述步骤S101中基于薄片分析获得的孔隙类型进行对比，如此，还可以判断出优化后的数字岩芯的孔隙结构类型，供后续分析所用。
S103、网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯。
在本说明书实施例中，经步骤S103后可以获取相对精确的3D数字岩芯模型。此时的3D数据模型就类似于一个魔方的正方体数据，而网格化就类似于把魔方分割成了很多构成魔方的小立方体网格，而每个小立方体网格也是一个小的3D数字岩芯数据体(也可以称为子网格或数字岩芯体)。
经过网格化后，网格化的数字岩芯比较容易地获得各子网格所代表的数字岩芯体的孔隙度。由于网格化的数字岩芯中的子网格数据巨大，这样可以获得一个孔隙度跨度非常宽的子网格集合。在取芯异常困通常只能获得有限数量岩芯的情况下，如果直接依据有限数量岩芯预测孔隙度，则孔隙度分布通常非常窄。而在本说明书实施例中，由于基于网格化的数字岩芯可以获得一个孔隙度跨度非常宽的子网格集合，从而有利于获得更为全面的孔隙度分布。
在本说明书一实施例中，每个子网格可以根据CT扫描分辨率用六面体单元进行网格剖分，每个子网格(即体素)等于CT扫描分辨率。此外，还可根据步骤S101中获得的XRD矿物分析分析结果，给每个体素的矿物骨架模量赋值。
S104、对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质。
在本说明书一实施例中，所述对网格化的第二数字岩芯进行有限元弹性性质模拟，可以包括如下步骤：
(1)、确定子网格骨架参数。由于步骤S103中已为每个体素的矿物骨架模量赋值。据此可以确定每个子网格的骨架参数。
(2)、接收输入的矿物骨架矿物参数。
(3)、接收输入的流体模量。
(4)、利用有限单元法(Finite Element Method，简称FEM)弹性性质模拟器(即有限元算法)，可以模拟获得每一个子网格体的弹性性质。
在本说明书一实施例中，FEM弹性性质模拟器例如可以是Garboczi(1995)提出的线弹性有限元法：将三维数字岩芯模型中的每个体素视作一个线弹性有限元网格，这样就免去了繁琐的网格剖分过程(本质上这个过程等同于用标准的六面体网格剖分求解区域)。对模型施加周期边界条件，便可将弹性位移分布的求解问题转化为最小位能原理求解的物理问题，数学上求解就是有限元法中典型的泛函求极值问题。例如，Arns(2002)提到可以使用快速共轭梯度法来对泛函求导方程进行迭代收敛以逼近系统最小能量，最后求得数值解即为六个方向的应力σxx,σyy,σzz,σxz,σyz,σxy与六个方向的应变εxx,εyy,εzz,εxz,εyz,εxy，进而求得等效弹性模量。
S105、根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
在本说明书一实施例中，所述根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构，包括：
根据公式
Figure BDA0002387885850000101
确定所述数字岩芯体的孔隙结构；
其中，ρ为数字岩芯体的密度，ρm为数字岩芯体分割的孔隙相中固体相矿物骨架的密度，ρf为数字岩芯体分割的孔隙相中充填的流体或气体的密度，φ为数字岩芯体的孔隙度，K为数字岩芯体的体积模量，VP为数字岩芯体在指定频段下随围压变化的纵波速度，VS为数字岩芯体在指定频段下随围压变化的横波速度，μ为数字岩芯体的剪切模量，Fk为中间变量，Km为数字岩芯体分割的孔隙相中固体相矿物骨架的体积模量，Kf为数字岩芯体分割的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为数字岩芯体的孔隙结构因子，γ越大代表了该数字岩芯体代表的岩石中以容易压缩的裂缝为主，γ越小代表了该数字岩芯代表的岩石中以不容易压缩的圆形铸模孔或者晶间孔为主。
由此可见，基于本说明书实施例的碳酸盐岩弹性性质与孔隙结构获取方法，可以获得网格化的数字岩芯体中，每个子网格对应的数字岩芯体的弹性性质和孔隙结构，从而实现了对碳酸盐岩弹性性质和孔隙结构的精确获取，因而提高了获取碳酸盐岩弹性性质和孔隙结构的准确性，进而可为下一步实施更加精确的储层预测以及储层分类提供了有利基础。
为便于理解，下面对本说明书实施例的碳酸盐岩弹性性质与孔隙结构获取方法进行举例说明。
本示例性实施例中选取了7块白云岩作为代表性碳酸盐岩储层岩石，分别记为Dolo-1、Dolo-2、Dolo-3-1、Dolo-3-2、Dolo-4-1、Dolo-4-2、Dolo-5，岩样直径38mm(或者直径25mm)，CT扫描成像分辨率为每体素20.7678μm。图2中示出了这7块碳酸盐岩岩样数字岩芯数据横向切片。通过对这7块岩样进行气测孔隙度，可以得到如下表1所示的实测孔隙度，从表1中可以看出：这些岩样有非常窄的孔隙度分布(0.6％～3.71％)。此外，通过对这7块岩样进行XRD矿物分析，表明这些岩样含平均约7.4％的石英和92.6％的白云石。就矿物成分而言，这就是碳酸盐岩岩样相对特殊的地方，可以认为碳酸盐岩为单矿物岩石，这在后续利用有限元进行弹性性质模拟时，可以简单地把数字岩芯二值化图像的骨架相其弹性性质设置成此单矿物岩石的模量。
表1岩样实测孔隙度与计算孔隙度物性表
岩样编号 实测孔隙度(％) 计算孔隙度(％)
Dolo-1 3.14 1.26
Dolo-2 2.89 1.92
Dolo-3-1 3.71 2.61
Dolo-3-2 2.67 1.09
Dolo-4-1 0.97 0.81
Dolo-4-2 2.25 2.14
Dolo-5 0.60 0.46
以Dolo-3-1岩样为例，图3中示出了Dolo-3-1岩样的纵、横波速度随围压变化的关系示意图；图4中示出了Dolo-3-1岩样的弹性模量随围压变化的关系示意图。利用公式1所示的灰度变换公式，对Dolo-3-1岩样的CT扫描图像进行对比度增强处理，图像的灰阶重新分布结果可如图5所示，原来非常窄的灰阶分布(深黑色数据点)被拓宽到全灰阶分布(浅黑数据点)，且孔隙相和固体相的波峰分布特征被明显增强。孔隙与固体过渡带都得到很好的灰阶分布拓宽，从而完成了合理的图像对比度增强。
图6a示出了Dolo-3-1岩样原始扫描图像，对CT扫描的原始图像进行剪裁，即选取边长为1000体素的正方体，从而可以得到如图6b所示切块后的原始图像，由于剪裁后的CT图像灰阶主要分布在非常低值区域，从图像的直观观测上显示为图像信息都比较暗。这样窄分布的灰阶直方图分布会给后续的图像分割造成困难，因此可以对图像进行对比度增强。对比度增强可通过如公式(2)实现，也就是通过灰度变换将原本窄分布的灰阶直方图变换成宽分布的直方图，拓宽孔隙相与固体相间的过渡带的灰阶分布，对比度增强后的结果如图6c所示，这样就能更好地方便后续图像分割阈值的选取。
利用公式(2)～(4)所示的各项异性滤波方可以获得如图7a所示的结果。结合图7a和图7b所示，经各向异性滤波后，固体颗粒内部与孔隙内部的噪音信号被成功平滑掉，而固体颗粒与孔隙间的边缘信号得到了很好地保留。
利用边缘检测法处理流程来进行边缘增强。由于固体颗粒与孔隙间的边缘一般对应于高梯度区域，因此可以使用梯度来提取边缘信号。通过图像梯度的累加求和函数可以确定梯度累加占总值90％的梯度点，高于这个梯度点就可以认为是高梯度指示相边界。结合图8a和图8b所示，可以看到该方法效果非常好，很完整地将固体颗粒与孔隙间的边缘区域给识别出来了，在此基础上进行图像分割可有利于降低分割不确定性，从而获得更好的图像分割效果。
通过对Dolo3-1号岩样的一个边缘增强后的图像切片(如图9a所示，)进行二值化图像处理，可以得到如图9b所示的二值化图像。在图9b中黑色表示孔隙，灰白色表示固体颗粒。如此，通过对Dolo3-1号岩样各切片进行相同的二值化图像处理，最终可获得经二值化图像处理后的结果，也即如图10a所示的3D岩石骨架图，此时图10a展示的既有岩石固体相也有孔隙相。此时，在该步骤可以用实验室实测孔隙度与数字岩芯估算孔隙度进行相互验证，标定图像处理流程的合理性，如果超出误差则与重新返回图像处理流程进行参数调整优化等处理，重新实施图像处理步骤。从而可以得到如图10b所示的数字岩芯的孔隙结构示意图。
本示例性实施例中，实验室氦气测得孔隙度为3.711％，而分割后得到成像孔隙度略小于氦气孔隙度，其值为2.297％。这是因为碳酸盐岩具有非常强的非均质性，其孔径可以小至几百个纳米而大至几个厘米。而且，图像分辨率与成像视域总有着不可协调的矛盾，若要使得分辨率越高则成像视域就会变小。所以一旦视域固定，那么其成像分辨率也是唯一不变的。对于成像孔隙度略小于实验孔隙度这样的现象一个较合理的解释是：碳酸盐岩中总有一些孔隙，其孔径甚至比CT分辨率还要小而无法成像，因此由CT成像估算的孔隙度会小于实验室氦气测量的孔隙度。由于微孔结构不能被完全分辨，CT图像上颗粒与微孔接触区域会表现出模糊状。这样的话分割过程就会产生不可避免的误差。而相对可信的分割过程通常会使分割后的孔隙度偏低。
用完全相同的处理方法对剩余六块白云岩CT数据进行处理，并得到相应二值化后的三维数字岩芯模型。对这七块白云岩在实验室条件下利用氦气孔隙度渗透率测量仪对白云岩岩样的孔隙度进行了测试，并借此来验证所建立三维数字岩芯模型的准确性，测得数据见上表1。总体而言，可以看到这七块白云岩的数字岩芯计算孔隙度都要略低于实验室氦气孔隙度测量值。由于碳酸盐岩非均质性强，其孔径分布非常广(几百纳米至几厘米)。一般会有一部分尺寸小于成像分辨率的孔无法在CT图像上清晰显示，所以造成数字岩芯成像孔隙度会略小于实验室氦气测量孔隙度。然而，实测孔隙度比较好地证明了数字岩芯图像建立过程中图像处理流程的合理性，另外，即便略有误差但对进一步进行规律性研究本项发明所提议的研究思路与流程均已足够。
另外，本示例性实施例中，还可以讲薄片分析所获得的孔隙结构类型，与如图10b所示的孔隙网络模型进行比较，从而可以给出该岩样的孔隙结构大致类型。
图11示出了某个岩样的网格化的3D数字岩芯，通过网格化，每块岩样可以得到125个子网格数据集，这样对于本示例性实施例中的7块岩样则可以得到125*7＝875个子网格数据集。对于875个子网格中的每个子图像可以容易获取其孔隙度(即利用数字岩芯估算的孔隙度)。由于子网格数据比较大，这样可以获得一个孔隙度跨度非常宽的子网格集合。此外，还可以根据XRD矿物分析结果，给每个子网格对应的数字岩芯体的进行矿物骨架模量赋值。
对获取的875个子网格所代表的每一个数字岩芯体，可以利用FEM弹性模拟器对每一个数据体进行有限元数值模拟得到了子网格数据体代表的数字岩芯体的弹性性质。与此同时，还可以定量地计算得到每个子网格数据体代表的数字岩芯体的孔隙结构参数γ。对分割所得的大量子网格数据体，对于每一个子网格可以计算得到其孔隙度φ、弹性性质(比如纵、横波速度等)和孔隙结构因子γ。对于本示例性实施例而言，还可获得如图12所示的碳酸盐岩岩样在分类下的弹性性质与孔隙度关系。
与上述的碳酸盐岩弹性性质与孔隙结构获取装置对应，本说明书实施例还提供了碳酸盐岩弹性性质与孔隙结构获取装置。参考图13所示，本说明书一些实施例的碳酸盐岩弹性性质与孔隙结构获取装置可以包括：
参数获取模块131，可以用于获取指定岩样的岩石物理参数；
数字岩芯构建模块132，可以用于基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；
岩芯网格化模块133，可以用于网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
弹性性质模拟模块134，可以用于对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
孔隙结构确定模块135，可以用于根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
在本说明书一些实施例的碳酸盐岩弹性性质与孔隙结构获取装置中，所述岩石物理参数包括：
孔隙度、成岩矿物组成、孔隙类型、密度及指定频段下随围压变化的纵、横波速度。
在本说明书一些实施例的碳酸盐岩弹性性质与孔隙结构获取装置中，所述基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，包括：
获取岩样的CT扫描图像；
将所述CT扫描图像进行图像二值化处理，并根据获得的二值化图像构建第一数字岩芯；
优化所述第一数字岩芯，以使其孔隙度与所述岩石物理参数中的孔隙度相一致。
在本说明书一些实施例的碳酸盐岩弹性性质与孔隙结构获取装置中，所述根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构，包括：
根据公式
Figure BDA0002387885850000141
确定所述数字岩芯体的孔隙结构；
其中，ρ为数字岩芯体的密度，ρm为数字岩芯体分割的孔隙相中固体相矿物骨架的密度，ρf为数字岩芯体分割的孔隙相中充填的流体或气体的密度，φ为数字岩芯体的孔隙度，K为数字岩芯体的体积模量，VP为数字岩芯体在指定频段下随围压变化的纵波速度，VS为数字岩芯体在指定频段下随围压变化的横波速度，μ为数字岩芯体的剪切模量，Fk为中间变量，Km为数字岩芯体分割的孔隙相中固体相矿物骨架的体积模量，Kf为数字岩芯体分割的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为数字岩芯体的孔隙结构因子。
参考图14所示，在本说明书另一些实施例中，碳酸盐岩弹性性质与孔隙结构获取装置可以包括存储器、处理器、以及存储在所述存储器上的计算机程序，其所述计算机程序被所述处理器运行时执行如下步骤：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
与上述的碳酸盐岩弹性性质与孔隙结构获取装置对应，本说明书实施例还提供了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
虽然上文描述的过程流程包括以特定顺序出现的多个操作，但是，应当清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行(例如使用并行处理器或多线程环境)。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本说明书时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。
内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁盘式存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法或者设备中还存在另外的相同要素。
本领域技术人员应明白，本说明书的实施例可提供为方法、系统或计算机程序产品。因此，本说明书可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本说明书的实施例而已，并不用于限制本说明书。对于本领域技术人员来说，本说明书可以有各种更改和变化。凡在本说明书的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本说明书的权利要求范围之内。

Claims (10)

1.一种碳酸盐岩弹性性质与孔隙结构获取方法，其特征在于，包括：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，包括：获取所述岩样的CT扫描图像；对所述CT扫描图像进行图像二值化处理，并根据获得的二值化图像构建第一数字岩芯，所述图像二值化处理依次包括对比度增强、滤波去噪、边缘增强和图像分割处理，所述对比度增强处理采用灰度变换公式实现，所述灰度变换公式为：
Figure FDA0003538879710000011
其中，f(x,y)表示原始二维图像上坐标(x,y)处像素点的灰度值，g(x,y)表示同一点灰度变换后的值，a、b、c、d均为常数；
所述滤波去噪采用各项异性滤波算法实现，公式为：
Figure FDA0003538879710000012
Figure FDA0003538879710000013
Figure FDA0003538879710000014
其中，It(s)表示原始图像的灰度函数，It+1(s)表示滤波后图像的灰度函数，s则是像素点在图像上的空间位置坐标；λ为控制扩散强度的常数，λ值越大滤波后图像越平滑；gK(x)为扩散函数，其中x为图像梯度值，K为判断噪声与微观结构信号的图像梯度阈值参数，梯度低于K值判断其为噪音信号，梯度高于K值就判断其为有用的微观结构信号；
Figure FDA0003538879710000015
为像素点与周围各方向相邻像素的差值；ηs＝{N,S,E,W}代表像素点s周围北、南、东、西四个方向的符号，K设置为公式(2)中每次迭代中梯度直方图积分(累计和)90％值所对应的梯度；
所述边缘增强采用边缘检测法实现，通过图像梯度的累加求和函数，确定梯度累加占总值90％的梯度点，根据高梯度指示相边界识别固体颗粒与孔隙间的边缘区域，所述高梯度指示相边界高于所述梯度点；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
2.如权利要求1所述的碳酸盐岩弹性性质与孔隙结构获取方法，其特征在于，所述岩石物理参数包括：
孔隙度、成岩矿物组成、孔隙类型、密度及指定频段下随围压变化的纵、横波速度。
3.如权利要求1所述的碳酸盐岩弹性性质与孔隙结构获取方法，其特征在于，所述基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，还包括：
优化所述第一数字岩芯，以使其孔隙度与所述岩石物理参数中的孔隙度相一致。
4.如权利要求1所述的碳酸盐岩弹性性质与孔隙结构获取方法，其特征在于，所述根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构，包括：
根据公式
Figure FDA0003538879710000021
确定所述数字岩芯体的孔隙结构；
其中，ρ为数字岩芯体的密度，ρm为数字岩芯体分割的孔隙相中固体相矿物骨架的密度，ρf为数字岩芯体分割的孔隙相中充填的流体或气体的密度，φ为数字岩芯体的孔隙度，K为数字岩芯体的体积模量，VP为数字岩芯体在指定频段下随围压变化的纵波速度，VS为数字岩芯体在指定频段下随围压变化的横波速度，μ为数字岩芯体的剪切模量，Fk为中间变量，Km为数字岩芯体分割的孔隙相中固体相矿物骨架的体积模量，Kf为数字岩芯体分割的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为数字岩芯体的孔隙结构因子。
5.一种碳酸盐岩弹性性质与孔隙结构获取装置，其特征在于，包括：
参数获取模块，用于获取指定岩样的岩石物理参数；
数字岩芯构建模块，用于基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，包括：用于获取所述岩样的CT扫描图像；用于对所述CT扫描图像进行图像二值化处理，并根据获得的二值化图像构建第一数字岩芯，所述图像二值化处理依次包括对比度增强、滤波去噪、边缘增强和图像分割处理，所述对比度增强处理采用灰度变换公式实现，所述灰度变换公式为：
Figure FDA0003538879710000031
其中，f(x,y)表示原始二维图像上坐标(x,y)处像素点的灰度值，g(x,y)表示同一点灰度变换后的值，a、b、c、d均为常数；
所述滤波去噪采用各项异性滤波算法实现，公式为：
Figure FDA0003538879710000032
Figure FDA0003538879710000033
Figure FDA0003538879710000034
其中，It(s)表示原始图像的灰度函数，It+1(s)表示滤波后图像的灰度函数，s则是像素点在图像上的空间位置坐标；λ为控制扩散强度的常数，λ值越大滤波后图像越平滑；gK(x)为扩散函数，其中x为图像梯度值，K为判断噪声与微观结构信号的图像梯度阈值参数，梯度低于K值判断其为噪音信号，梯度高于K值就判断其为有用的微观结构信号；
Figure FDA0003538879710000035
为像素点与周围各方向相邻像素的差值；ηs＝{N,S,E,W}代表像素点s周围北、南、东、西四个方向的符号，K设置为公式(2)中每次迭代中梯度直方图积分(累计和)90％值所对应的梯度；
所述边缘增强采用边缘检测法实现，通过图像梯度的累加求和函数，确定梯度累加占总值90％的梯度点，根据高梯度指示相边界识别固体颗粒与孔隙间的边缘区域，所述高梯度指示相边界高于所述梯度点；
岩芯网格化模块，用于网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
弹性性质模拟模块，用于对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
孔隙结构确定模块，用于根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
6.如权利要求5所述的碳酸盐岩弹性性质与孔隙结构获取装置，其特征在于，所述岩石物理参数包括：
孔隙度、成岩矿物组成、孔隙类型、密度及指定频段下随围压变化的纵、横波速度。
7.如权利要求5所述的碳酸盐岩弹性性质与孔隙结构获取装置，其特征在于，所述基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，还包括：
优化所述第一数字岩芯，以使其孔隙度与所述岩石物理参数中的孔隙度相一致。
8.如权利要求5所述的碳酸盐岩弹性性质与孔隙结构获取装置，其特征在于，所述根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构，包括：
根据公式
Figure FDA0003538879710000041
确定所述数字岩芯体的孔隙结构；
其中，ρ为数字岩芯体的密度，ρm为数字岩芯体分割的孔隙相中固体相矿物骨架的密度，ρf为数字岩芯体分割的孔隙相中充填的流体或气体的密度，φ为数字岩芯体的孔隙度，K为数字岩芯体的体积模量，VP为数字岩芯体在指定频段下随围压变化的纵波速度，VS为数字岩芯体在指定频段下随围压变化的横波速度，μ为数字岩芯体的剪切模量，Fk为中间变量，Km为数字岩芯体分割的孔隙相中固体相矿物骨架的体积模量，Kf为数字岩芯体分割的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为数字岩芯体的孔隙结构因子。
9.一种碳酸盐岩弹性性质与孔隙结构获取装置，包括存储器、处理器、以及存储在所述存储器上的计算机程序，其特征在于，所述计算机程序被所述处理器运行时执行如下步骤：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，包括：获取所述岩样的CT扫描图像；对所述CT扫描图像进行图像二值化处理，并根据获得的二值化图像构建第一数字岩芯，所述图像二值化处理依次包括对比度增强、滤波去噪、边缘增强和图像分割处理，所述对比度增强处理采用灰度变换公式实现，所述灰度变换公式为：
Figure FDA0003538879710000051
其中，f(x,y)表示原始二维图像上坐标(x,y)处像素点的灰度值，g(x,y)表示同一点灰度变换后的值，a、b、c、d均为常数；
所述滤波去噪采用各项异性滤波算法实现，公式为：
Figure FDA0003538879710000052
Figure FDA0003538879710000053
Figure FDA0003538879710000054
其中，It(s)表示原始图像的灰度函数，It+1(s)表示滤波后图像的灰度函数，s则是像素点在图像上的空间位置坐标；λ为控制扩散强度的常数，λ值越大滤波后图像越平滑；gK(x)为扩散函数，其中x为图像梯度值，K为判断噪声与微观结构信号的图像梯度阈值参数，梯度低于K值判断其为噪音信号，梯度高于K值就判断其为有用的微观结构信号；
Figure FDA0003538879710000055
为像素点与周围各方向相邻像素的差值；ηs＝{N,S,E,W}代表像素点s周围北、南、东、西四个方向的符号，K设置为公式(2)中每次迭代中梯度直方图积分(累计和)90％值所对应的梯度；
所述边缘增强采用边缘检测法实现，通过图像梯度的累加求和函数，确定梯度累加占总值90％的梯度点，根据高梯度指示相边界识别固体颗粒与孔隙间的边缘区域，所述高梯度指示相边界高于所述梯度点；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
10.一种计算机存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现以下步骤：
获取指定岩样的岩石物理参数；
基于所述岩石物理参数及所述指定岩样的CT扫描图像构建第一数字岩芯，包括：获取所述岩样的CT扫描图像；对所述CT扫描图像进行图像二值化处理，并根据获得的二值化图像构建第一数字岩芯，所述图像二值化处理依次包括对比度增强、滤波去噪、边缘增强和图像分割处理，所述对比度增强处理采用灰度变换公式实现，所述灰度变换公式为：
Figure FDA0003538879710000061
其中，f(x,y)表示原始二维图像上坐标(x,y)处像素点的灰度值，g(x,y)表示同一点灰度变换后的值，a、b、c、d均为常数；
所述滤波去噪采用各项异性滤波算法实现，公式为：
Figure FDA0003538879710000062
Figure FDA0003538879710000063
Figure FDA0003538879710000064
其中，It(s)表示原始图像的灰度函数，It+1(s)表示滤波后图像的灰度函数，s则是像素点在图像上的空间位置坐标；λ为控制扩散强度的常数，λ值越大滤波后图像越平滑；gK(x)为扩散函数，其中x为图像梯度值，K为判断噪声与微观结构信号的图像梯度阈值参数，梯度低于K值判断其为噪音信号，梯度高于K值就判断其为有用的微观结构信号；
Figure FDA0003538879710000065
为像素点与周围各方向相邻像素的差值；ηs＝{N,S,E,W}代表像素点s周围北、南、东、西四个方向的符号，K设置为公式(2)中每次迭代中梯度直方图积分(累计和)90％值所对应的梯度；
所述边缘增强采用边缘检测法实现，通过图像梯度的累加求和函数，确定梯度累加占总值90％的梯度点，根据高梯度指示相边界识别固体颗粒与孔隙间的边缘区域，所述高梯度指示相边界高于所述梯度点；
网格化所述第一数字岩芯，并根据所述岩石物理参数对网格化的第一数字岩芯进行矿物骨架模量赋值，获得网格化的第二数字岩芯；
对所述网格化的第二数字岩芯进行有限元弹性性质模拟，获得所述网格化的第二数字岩芯中每个子网格对应的数字岩芯体的弹性性质；
根据所述弹性性质及所述岩石物理参数确定所述数字岩芯体的孔隙结构。
CN202010104071.1A 2020-02-20 2020-02-20 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质 Active CN111426616B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010104071.1A CN111426616B (zh) 2020-02-20 2020-02-20 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010104071.1A CN111426616B (zh) 2020-02-20 2020-02-20 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质

Publications (2)

Publication Number Publication Date
CN111426616A CN111426616A (zh) 2020-07-17
CN111426616B true CN111426616B (zh) 2022-05-24

Family

ID=71547208

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010104071.1A Active CN111426616B (zh) 2020-02-20 2020-02-20 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质

Country Status (1)

Country Link
CN (1) CN111426616B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111239821B (zh) * 2020-02-20 2021-02-12 中国石油大学(北京) 碳酸盐岩储层孔隙结构预测方法、装置、设备及存储介质
CN112051200B (zh) * 2020-08-26 2022-02-15 中国石油大学(北京) 一种致密砂岩储层孔隙结构的定量评价方法和装置
CN115359348B (zh) * 2022-07-20 2023-07-25 中南大学 岩芯特征识别统计方法及系统、设备、存储介质

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102288986A (zh) * 2011-05-16 2011-12-21 中国石油大学(北京) 一种地震尺度下碳酸盐岩储层弹性模量获取方法
CN103822865A (zh) * 2014-03-20 2014-05-28 中国石油大学(华东) 一种高分辨率三维数字岩心建模方法
CN104516017A (zh) * 2013-09-29 2015-04-15 中国石油化工股份有限公司 一种碳酸盐岩岩石物理参数地震反演方法
CN106290105A (zh) * 2016-07-20 2017-01-04 中国石油大学(华东) 一种碳酸盐岩储层溶蚀孔隙体积含量预测方法
CN106368691A (zh) * 2015-07-24 2017-02-01 中国石油化工股份有限公司 基于岩石物理地震信息三维异常孔隙压力预测方法
CN107179562A (zh) * 2017-04-27 2017-09-19 恒泰艾普集团股份有限公司 相控岩石物理模型指导下的储层预测方法
CN109063348A (zh) * 2018-08-09 2018-12-21 中国石油天然气股份有限公司 基于孔隙喉道网络模型的驱替模拟方法及装置
CN109425913A (zh) * 2017-08-22 2019-03-05 中国石油化工股份有限公司 碳酸盐岩储层含气敏感弹性参数优选方法及系统
CN109887073A (zh) * 2018-08-16 2019-06-14 清能艾科（深圳）能源技术有限公司 岩心的三维数字模型构建方法及装置

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10830713B2 (en) * 2017-11-20 2020-11-10 DigiM Solution LLC System and methods for computing physical properties of materials using imaging data

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102288986A (zh) * 2011-05-16 2011-12-21 中国石油大学(北京) 一种地震尺度下碳酸盐岩储层弹性模量获取方法
CN104516017A (zh) * 2013-09-29 2015-04-15 中国石油化工股份有限公司 一种碳酸盐岩岩石物理参数地震反演方法
CN103822865A (zh) * 2014-03-20 2014-05-28 中国石油大学(华东) 一种高分辨率三维数字岩心建模方法
CN106368691A (zh) * 2015-07-24 2017-02-01 中国石油化工股份有限公司 基于岩石物理地震信息三维异常孔隙压力预测方法
CN106290105A (zh) * 2016-07-20 2017-01-04 中国石油大学(华东) 一种碳酸盐岩储层溶蚀孔隙体积含量预测方法
CN107179562A (zh) * 2017-04-27 2017-09-19 恒泰艾普集团股份有限公司 相控岩石物理模型指导下的储层预测方法
CN109425913A (zh) * 2017-08-22 2019-03-05 中国石油化工股份有限公司 碳酸盐岩储层含气敏感弹性参数优选方法及系统
CN109063348A (zh) * 2018-08-09 2018-12-21 中国石油天然气股份有限公司 基于孔隙喉道网络模型的驱替模拟方法及装置
CN109887073A (zh) * 2018-08-16 2019-06-14 清能艾科（深圳）能源技术有限公司 岩心的三维数字模型构建方法及装置

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
A rock-physical modeling method for carbonate reservoirs at seismic scale;Li Jing-Ye et al.;《APPLIED GEOPHYSICS》;29130331;第10卷(第1期);第1-14页 *
使用微米CT成像的孔隙结构对碳酸盐岩声学性质的影响研究;胡洋铭等;《CPS/SEG北京2018国际地球物理会议暨展览电子论文集》;20180424;第1104-1105页 *
基于微CT技术的砂岩数字岩石物理实验;刘向君等;《地球物理学报》;20140430;第57卷(第4期);第1133-1139页 *

Also Published As

Publication number Publication date
CN111426616A (zh) 2020-07-17

Similar Documents

Publication Publication Date Title
EP2972303B1 (en) Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties
CN111426616B (zh) 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质
Andrä et al. Digital rock physics benchmarks—Part I: Imaging and segmentation
Garing et al. Pore-scale capillary pressure analysis using multi-scale X-ray micromotography
CN111239821B (zh) 碳酸盐岩储层孔隙结构预测方法、装置、设备及存储介质
Ehrlich et al. Petrography and reservoir physics I: objective classification of reservoir porosity (1)
US20200183031A1 (en) Automated seismic interpretation-guided inversion
AU2011258594B2 (en) Method for obtaining consistent and integrated physical properties of porous media
Berg et al. Generation of ground truth images to validate micro-CT image-processing pipelines
US20140270393A1 (en) Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties
Ketcham et al. Accurate measurement of small features in X‐ray CT data volumes, demonstrated using gold grains
Adeleye et al. Pore-scale analyses of heterogeneity and representative elementary volume for unconventional shale rocks using statistical tools
WO2014142976A1 (en) Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties
Van Hateren et al. Identifying sediment transport mechanisms from grain size–shape distributions, applied to aeolian sediments
Silva et al. Permeability modeling of a basin-bounding fault damage zone in the Rio do Peixe Basin, Brazil
Cappuccio et al. Three-dimensional separation and characterization of fractures in X-ray computed tomographic images of rocks
Wang et al. Joint modeling based on singularity mapping and U-statistical methods for geo-anomaly characterization
Aitken et al. Semiautomated quantification of the influence of data richness on confidence in the geologic interpretation of aeromagnetic maps
Belila et al. Petrophysical characterization of coquinas from Morro do Chaves Formation (Sergipe-Alagoas Basin) by X-ray computed tomography
Qi* et al. Segmentation of salt domes, mass transport complexes on 3D seismic data volumes using Kuwahara windows and multiattribute cluster analysis
Gupta et al. Algebraic Reconstruction Techniques
Rahimov et al. Use of local binary pattern in texture classification of carbonate rock micro-CT images
CN112231881A (zh) 基质渗透率计算方法及系统
Jing et al. DigiCoal: a numerical toolbox for fractured coal characterisation
CN113325474B (zh) 生物礁判别方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant