本说明书实施例提供了一种碳酸盐岩储层孔隙结构预测方法、装置、设备及存储介质，该方法包括：获取基于地震叠前反演的纵波速度数据体和密度数据体；根据测井数据确定孔隙度与纵波速度、密度的第一关系模型，根据数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型，并根据这两个模型形成第三关系模型；基于测井数据获取第一孔隙结构因子，基于数字岩芯数据获取第二孔隙结构因子，并根据这个两个孔隙结构因子形成第三孔隙结构因子；将纵波速度数据体和密度数据体输入第三关系模型，获得孔隙度数据体；将纵波速度数据体与孔隙度数据体交会，获得孔隙结构数据体。本说明书实施例可以提高碳酸盐岩储层孔隙结构预测的准确性。
碳酸盐岩储层孔隙结构预测方法、装置、设备及存储介质
技术领域
[0001] 本说明书涉及油气勘探开发技术领域，尤其是涉及一种碳酸盐岩储层孔隙结构预测方法、装置、设备及存储介质。
背景技术
[0002] 目前，现有技术中的储层特征参数预测，常见于生产中采用的属性交会的方式。特别地对于碎屑岩储层，利用属性交会的方法能够相对准确地获得储层孔隙度与饱和度，从而由此出发获得储层的结构信息。
[0003] 然而，对于碳酸盐岩储层，由于其后期的成岩改造作用会导致其次生孔隙非常发育，使其孔隙结构极其复杂。碳酸盐岩中常见的孔隙结构类型有铸模孔、粒内溶孔、粒间溶孔、晶间孔、裂缝等，这些复杂多变的孔隙结构会对碳酸盐岩的弹性性质产生显著的影响，致使常用的属性交会方式无法获得准确的孔隙度与饱和度参数，更无法获得孔隙结构参数。因此，针对碳酸盐岩储层，如何提高碳酸盐岩储层孔隙结构预测的准确性已成为目前亟待解决的技术难题。
发明内容
[0004] 本说明书实施例的目的在于提供一种碳酸盐岩储层孔隙结构预测方法、装置、设备及存储介质，以提高碳酸盐岩储层孔隙结构预测的准确性。
[0005] 为达到上述目的，一方面，本说明书实施例提供了一种碳酸盐岩储层孔隙结构预测方法，包括：
[0006] 获取目标区块基于地震叠前反演的纵波速度数据体和密度数据体；
[0007] 根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，并根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型；
[0008] 将所述第一关系模型和所述第二关系模型的加权和作为第三关系模型；
[0009] 基于所述测井数据获取所述目标区块的第一孔隙结构因子，并基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子；
[0010] 将所述第一孔隙结构因子和所述第二孔隙结构因子的加权和作为第三孔隙结构因子；
[0011] 将所述纵波速度数据体和所述密度数据体输入所述第三关系模型，获得所述目标区块的孔隙度数据体；
[0012] 将所述纵波速度数据体与所述孔隙度数据体进行交会，获得所述目标区块的孔隙结构数据体。
[0013] 另一方面，本说明书实施例还提供了一种碳酸盐岩储层孔隙结构预测装置，包括：
[0014] 第一数据体获取模块，用于获取目标区块基于地震叠前反演的纵波速度数据体和密度数据体；
[0015] 第一模型获取模块，用于根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，并根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型；
[0016] 第二模型获取模块，用于将所述第一关系模型和所述第二关系模型的加权和作为第三关系模型；
[0017] 第一因子获取模块，用于基于所述测井数据获取所述目标区块的第一孔隙结构因子，并基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子；
[0018] 第二因子获取模块，用于将所述第一孔隙结构因子和所述第二孔隙结构因子的加权和作为第三孔隙结构因子；
[0019] 第二数据体获取模块，用于将所述纵波速度数据体和所述密度数据体输入所述第三关系模型，获得所述目标区块的孔隙度数据体；
[0020] 孔隙结构获取模块，用于将所述纵波速度数据体与所述孔隙度数据体进行交会，获得所述目标区块的孔隙结构数据体。
[0021] 另一方面，本说明书实施例还提供了一种计算机设备，包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时执行如下步骤：
[0022] 获取目标区块基于地震叠前反演的纵波速度数据体和密度数据体；
[0023] 根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，并根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型；
[0024] 将所述第一关系模型和所述第二关系模型的加权和作为第三关系模型；
[0025] 基于所述测井数据获取所述目标区块的第一孔隙结构因子，并基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子；
[0026] 将所述第一孔隙结构因子和所述第二孔隙结构因子的加权和作为第三孔隙结构因子；
[0027] 将所述纵波速度数据体和所述密度数据体输入所述第三关系模型，获得所述目标区块的孔隙度数据体；
[0028] 将所述纵波速度数据体与所述孔隙度数据体进行交会，获得所述目标区块的孔隙结构数据体。
[0029] 另一方面，本说明书实施例还提供了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
[0030] 获取目标区块基于地震叠前反演的纵波速度数据体和密度数据体；
[0031] 根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，并根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型；
[0032] 将所述第一关系模型和所述第二关系模型的加权和作为第三关系模型；
[0033] 基于所述测井数据获取所述目标区块的第一孔隙结构因子，并基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子；
[0034] 将所述第一孔隙结构因子和所述第二孔隙结构因子的加权和作为第三孔隙结构因子；
[0035] 将所述纵波速度数据体和所述密度数据体输入所述第三关系模型，获得所述目标区块的孔隙度数据体；
[0036] 将所述纵波速度数据体与所述孔隙度数据体进行交会，获得所述目标区块的孔隙结构数据体。
[0037] 由以上本说明书实施例提供的技术方案可见，在本说明书实施例中，结合目标区块的测井数据和数字岩芯数据，获得了目标区块的孔隙度与纵波速度、密度的关系模型；并结合目标区块的测井数据和数字岩芯数据，获得了目标区块的孔隙结构因子；在此基础上，利用目标区块的地震叠前反演的纵波速度数据体和密度数据体作为输入，输入到目标区块的表征孔隙度与纵波速度、密度的关系模型中，从而可以得到更为精确的目标区块的孔隙度数据体；而通过将目标区块的地震叠前反演的纵波速度数据体与该目标区块的孔隙度数据体进行交会，从而可以得到更为精确的目标区块的孔隙结构数据体。进而为碳酸盐岩储层的储量预测、前期勘探开发、后期生产监控等提供了更为可靠的基础数据。
附图说明
[0038] 为了更清楚地说明本说明书实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本说明书中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。在附图中：
[0039] 图1为本说明书一些实施例中碳酸盐岩储层孔隙结构预测方法的流程图；
[0040] 图2为本说明书一实施例中碳酸盐岩岩样的CT扫描横向切片示意图；
[0041] 图3为本说明书一实施例中碳酸盐岩岩样的纵、横波速度随围压变化的关系示意图；
[0042] 图4为本说明书一实施例中碳酸盐岩岩样的弹性模量随围压变化的关系示意图；
[0043] 图5为本说明书一实施例中线性拉伸法的灰阶重新分布结果；
[0044] 图6a为本说明书一实施例中原始CT扫描图像的取样示意图(灰色矩形框为取样部分)；
[0045] 图6b为从图6a所示的原始CT扫描图像中取出的CT扫描图像切块的示意图；
[0046] 图6c为对图6b所示的CT扫描图像切块进行对比度增强处理后的结果示意图；
[0047] 图7a为本说明书一实施例中CT扫描图像在进行对比度增强处理后的结果示意图；
[0048] 图7b对图7a所示的CT扫描图像进行各向异性滤波处理后的结果示意图；
[0049] 图8a为本说明书一实施例中CT扫描图像在进行各向异性滤波处理后的结果示意图；
[0050] 图8b对图8a所示的CT扫描图像进行边缘增强处理后的结果示意图；
[0051] 图9a为本说明书一实施例中预处理后的CT扫描图像的示意图；
[0052] 图9b对图9a所示的CT扫描图像进行图像二值化处理后的结果示意图；
[0053] 图10a为本说明书一实施例中3-1号碳酸盐岩岩样的数字岩芯示意图；
[0054] 图10b为图10a所示的数字岩芯的孔隙结构示意图；
[0055] 图11为本说明书一实施例中数字岩芯的网格化示意图；
[0056] 图12为本说明书一实施例中获得的碳酸盐岩岩样在分类下的弹性性质与孔隙度关系示意图；
[0057] 图13a～图13c为本说明书一实施例中获得测井数据示意图；
[0058] 图14为本说明书一实施例中的神经网络示意图；
[0059] 图15为本说明书一实施例中基于神经网络训练出的模型所预测的孔隙度与真实孔隙度的验证误差示意图；
[0060] 图16为本说明书一实施例中基于地震叠前AVO反演获得的纵波速度属性数据体的剖面示意图；
[0061] 图17为本说明书一实施例中基于地震叠前AVO反演获得的密度属性数据体的剖面示意图；
[0062] 图18为本说明书一实施例中获得的目标区块的孔隙度数据体的剖面示意图；
[0063] 图19为本说明书一实施例中目标区块的孔隙度数据体与纵波速度数据体的交会示意图；
[0064] 图20为本说明书一实施例中获得的目标区块的孔隙结构数据体的剖面示意图；
[0065] 图21为本说明书一些实施例中碳酸盐岩储层孔隙结构预测装置的结构框图；
[0066] 图22为本说明书一些实施例中计算机设备的结构框图。
具体实施方式
[0067] 为了使本技术领域的人员更好地理解本说明书中的技术方案，下面将结合本说明书实施例中的附图，对本说明书实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明书一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本说明书保护的范围。
[0068] 参考图1所示，本说明书一些实施例的碳酸盐岩储层孔隙结构预测方法可以包括以下步骤：
[0069] S101、获取目标区块基于地震叠前反演的纵波速度数据体和密度数据体。
[0070] S102、根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，并根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型。
[0071] S103、将所述第一关系模型和所述第二关系模型的加权和作为第三关系模型。
[0072] S104、基于所述测井数据获取所述目标区块的第一孔隙结构因子，并基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子。
[0073] S105、将所述第一孔隙结构因子和所述第二孔隙结构因子的加权和作为第三孔隙结构因子。
[0074] S106、将所述纵波速度数据体和所述密度数据体输入所述第三关系模型，获得所述目标区块的孔隙度数据体。
[0075] S107、将所述纵波速度数据体与所述孔隙度数据体进行交会，获得所述目标区块的孔隙结构数据体。
[0076] 由此可见，在本说明书实施例中，结合目标区块的测井数据和数字岩芯数据，获得了目标区块的孔隙度与纵波速度、密度的关系模型；并结合目标区块的测井数据和数字岩芯数据，获得了目标区块的孔隙结构因子；在此基础上，利用目标区块的地震叠前反演的纵波速度数据体和密度数据体作为输入，输入到目标区块的表征孔隙度与纵波速度、密度的关系模型中，从而可以得到更为精确的目标区块的孔隙度数据体；而通过将目标区块的地震叠前反演的纵波速度数据体与该目标区块的孔隙度数据体进行交会，从而可以得到更为精确的目标区块的孔隙结构数据体。进而为碳酸盐岩储层的储量预测、前期勘探开发、后期生产监控等提供了更为可靠的基础数据。
[0077] 在本说明书一实施例中，基于地震叠前反演的纵波速度数据体和密度数据体例如可以为，基于地震AVO(Amplitude variation with offset，振幅随偏移距的变化)反演获得的纵波速度数据体和密度数据体。可以预先对目标区块进行地震AVO反演，从而获得包括目标区块的纵波速度数据体和密度数据体等在内的反演结果。如此，通过读取反演结果就可以得到目标区块的纵波速度数据体和密度数据体。当然，本说明书不限于此，在其他实施例中，根据需要也可以采用其他地震叠前反演技术。
[0078] 在本说明书一实施例中，可以预先对目标区块进行测井，从而获得相应的测井数据。所述测井数据一般可以包括密度、纵波速度、横波速度和孔隙度等数据。通过读取测井数据就可以得到目标区块基于测井数据的孔隙度、纵波速度和密度。在此基础上，可以以测井数据中的纵波速度和密度作为输入，以测井数据中的孔隙度作为目标输出，训练预设的机器学习模型(例如神经网络等)，从而可以得到孔隙度与纵波速度、密度的关系曲线，即确定孔隙度与纵波速度、密度的第一关系模型。由于在拟合该关系式时，神经网络相比常用的非线性拟合算法，更容易加入数据集，并且能够表征的关系式也更加复杂，理论上只要孔隙度与纵波速度和密度间存在一个近似关系式，无论多么复杂，通过神经网络训练总能找到一个表征孔隙度与纵波速度、密度的关系曲线。
[0079] 例如，在本说明书一示例性实施例中，某工区内3口井(即井1、井2和井3)的测井数据(纵波速度、密度与孔隙度)可以如图13a～图13c所示。在测井数据的基础上，使用神经网络进行训练，如图14所示，从而可以得到表征孔隙度与纵波速度、密度的关系曲线。在一示例性实施例中，可以将上述测井数据对应的数据集分成训练集与验证集，验证集占总数据的20％。当利用训练后得到的关系曲线出的孔隙度后，可以将该孔隙度与真实孔隙度进行对比，通过对比可以验证关系曲线的预测准确性，例如图15所示。
[0080] 在本说明书一实施例中，所述根据目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型可以包括以下步骤：
[0081] (一)、获取所述目标区块的指定岩样的岩石物理参数。
[0082] 在本说明书的实施例中，所述岩石物理参数例如可以包括孔隙度、成岩矿物组成、孔隙类型、密度及指定频段下随围压变化的纵、横波速度等。在一实施例中，所述的获取例如可以接收外部输入的指定岩样的岩石物理参数。在另一实施例中，所述的获取例如还以是从指定存储路径读取记录有上述指定岩样的岩石物理参数。
[0083] 在本说明书一实施例中，上述指定岩样的岩石物理参数可以预先通过以下方式得到：
[0084] (1)、可以选取储层代表性碳酸盐岩样若干，一般选取至少5个岩样，岩样孔隙类型粗略地覆盖裂缝型、孔洞型、孔洞-裂缝型，并可以对岩样进行洗盐、洗油、以及烘干干燥等预处理。
[0085] (2)、可以对选取的代表性碳酸盐岩岩样进行氦气测量孔隙度，或者可以通过核磁共振测量孔隙度，获得如上代表性岩样比较精确的孔隙度，以用作后续标定。
[0086] (3)、可以对选取的代表性碳酸盐岩岩样进行X射线衍射(X-Ray Diffraction，简称XRD)矿物分析，以获取成岩矿物，用作后续弹性性质模拟的矿物骨架模量参数输入。
[0087] (4)、可以对选取的代表性碳酸盐岩岩样进行镜下薄片分析，以获取定性半定量的孔隙类型，用作后续标定。
[0088] (5)、可以对选取的代表性碳酸盐岩岩样进行密度测量及高频超声频段下随围压变化的纵横波速度测量，以用作后续标定。
[0089] 其中，根据需要，以上步骤(2)～(5)可按实际情况在实验室测量时微调。
[0090] 例如，在说明书一示例性实施例中，选取了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％的白云石。就矿物成分而言，这就是碳酸盐岩岩样相对特殊的地方，可以认为碳酸盐岩为单矿物岩石，这在后续利用有限元进行弹性性质模拟时，可以简单地把数字岩芯二值化图像的骨架相其弹性性质设置成此单矿物岩石的模量。
[0091] 表1岩样实测孔隙度与计算孔隙度物性表
[0092]岩样编号 实测孔隙度(％) 计算孔隙度(％)
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
[0093] 以Dolo-3-1岩样为例，图3中示出了Dolo-3-1岩样的纵、横波速度随围压变化的关系示意图；图4中示出了Dolo-3-1岩样的弹性模量随围压变化的关系示意图。
[0094] (二)、基于所述岩石物理参数及所述指定岩样的CT扫描图像构建数字岩芯。
[0095] 在本说明书一实施例中，所述基于所述岩石物理参数及所述指定岩样的CT扫描图像构建数字岩芯，可以包括如下步骤：
[0096] (1)、获取岩样的CT扫描图像。其中，CT扫描图像可以是预先对选取的代表性碳酸盐岩岩样进行CT扫描成像获得原始图像剖面。
[0097] (2)、将所述CT扫描图像进行图像二值化处理，并根据获得的二值化图像构建数字岩芯。
[0098] 在本说明书一实施例中，所述将CT扫描图像进行图像二值化处理可以包括如下步骤：
[0099] (21)、对比度增强
[0100] 首先可以对CT扫描的原始图像进行剪裁，一般选取边长为1000体素的正方体。剪裁后的CT图像灰阶主要分布在非常低值区域，从图像的直观观测上显示为图像信息都比较暗，这样窄分布的灰阶直方图分布会给后续的图像分割造成困难，因此可以对图像进行对比度增强。对比度增强操作通过灰度变换将原本窄分布的灰阶直方图变换成宽分布的直方图。选用灰度变换公式(1)实现上述功能：
[0101]
[0102] 其中，f(x,y)表示原始二维图像上坐标(x,y)处像素点的灰度值，g(x,y)表示同一点灰度变换后的值，a、b、c、d均为常数。
[0103] (22)、滤波去噪
[0104] 在本说明书一些实施例中，可以选用中值滤波、均值滤波或高斯滤波等全局滤波类算法对图像进行去噪处理。全局滤波类算法是统一对待图像上每一个像素点，容易对裂缝纹理等信息平滑掉。而碳酸盐岩相比于其他类别的岩石类型具有更加丰富的孔隙结构类型，特别是该类岩石中遍布裂缝、纹理、孔洞等复杂的结构，而这种复杂孔隙结构是如何影响碳酸盐岩弹性性质，是碳酸盐岩地球物理勘探中的重要研究内容之一。因此，获取既能合理去除图像中噪声，又能有效保留裂缝、纹理、孔洞等孔隙结构信息的滤波去噪算法，是碳酸盐岩数字岩芯图像处理中的重要任务。在本说明书一实施例中，可以使用各项异性滤波算法实现上述期望功能，具体可以如公式(2)～(4)所示：
[0105]
[0106]
[0107]
[0108] 上述公式(2)中，It(s)表示原始图像的灰度函数，It+1(s)表示滤波后图像的灰度函数，s则是像素点在图像上的空间位置坐标，λ为控制扩散强度的常数，λ值越大滤波后图像越平滑，gK(x)为扩散函数，其中x为图像梯度值，K为判断噪声与微观结构信号的图像梯度阈值参数，梯度低于K值判断其为噪音信号，梯度高于K值就判断其为有用的微观结构信号， 为像素点与周围各方向相邻像素的差值，其表达式见公式(4)，ηs＝{N,S,E,W}代表像素点s周围北、南、东、西等四个方向的符号。K设置为等于迭代公式(2)中每次迭代中的梯度直方图积分(累积和)90％值所对应的梯度。
[0109] (23)、边缘增强
[0110] 对数字岩芯进行图像边缘增强可以使图像边缘更清晰。边缘增强主要困难在于锐化滤波器很容易放大图像中所有噪声。在本说明书一实施例中，可以采用边缘检测法来进行边缘增强。由于固体颗粒与孔隙间的边缘一般对应于高梯度区域，因此可以使用梯度来提取边缘信号。通过图像梯度的累加求和函数，可以确定梯度累加占总值90％的梯度点，高于这个梯度点，就可以认为是高梯度指示相边界。如此，可以很完整地将固体颗粒与孔隙间的边缘区域给识别出来了，在此基础上进行图像分割将降低分割不确定性，从而获得更好的图像分割效果。
[0111] (24)、图像分割
[0112] 经过上述(21)～(23)处理步骤后，可以对碳酸盐岩数字岩芯图像进行图像二值化处理，从而可以建立初始的数字岩芯，以便于后续实施弹性性质模拟。
[0113] 在本说明书实施例中，将孔隙结构从CT图像中提取出来，识别每张图像中的孔隙像素，这个识别过程在图像处理中叫做图像分割。图像分割技术一般分为几何形状分割法和聚类分割法。前者是分析灰度直方图几何形状来确定分割阈值的方法，而后者是通过图像中两组数据间的统计属性差异来确定最佳阈值的方法。聚类分割法主要包括以下三个步骤：
[0114] (241)、由于图像中每个像素点都对应一个灰度值，聚类算法则将图像中所有像素点的灰度值视作一个数据集合。
[0115] (242)、由于CT图像是灰度图像，其灰度分布一般有8bit(0～255)或16bit(0～65535)这两种，以8bit为例，阈值从0取到255共256个值，每个阈值都将数据集合分为两类。
[0116] (243)、计算每个阈值分类后两类数据的类间方差与类内方差之比，比值最大点所对应的阈值便是聚类分割法法找到的阈值。与普通阈值分割思路类似，将分出的两类数据的灰度值进行重新赋值，一般是分别赋予最大差距值，这个过程称为图像二值化处理。聚类分割算法中由于Otsu算法基于灰度直方图上单一阈值的特性，尤其适合单矿物岩样。而碳酸盐岩矿物成分比较单一，使用Otsu算法对CT图像进行分割，获得了比较好的效果。
[0117] (3)、优化所述数字岩芯，以使其孔隙度与所述岩石物理参数中的孔隙度相一致。
[0118] 基于如上步骤获取的二值化数字岩芯数据体，可以获取孔隙相对于总体积的占比，即基于数字岩芯确定的孔隙度，可以对于所有选取的碳酸盐岩岩样，比较基于数字岩芯获取的孔隙度与此前利用氦气与NMR测量孔隙度，从而实现孔隙度标定工作。如果基于数字岩芯估计的孔隙度与实测度超出误差允许，则可以重新返回二值化图像处理流程，对其中的各个子步骤，即上述步骤(21)～(24)进行参数优选及微调，直至基于数字岩芯估计孔隙度与实测孔隙度相一致，从而可以获得优化后的数字岩芯。此外，优化后的数字岩芯还可以与上述步骤中基于薄片分析获得的孔隙类型进行对比，如此，还可以判断出优化后的数字岩芯的孔隙结构类型，供后续分析所用。
[0119] 例如，在说明书一示例性实施例中，利用公式1所示的灰度变换公式，对Dolo-3-1岩样的CT扫描图像进行对比度增强处理，图像的灰阶重新分布结果可如图5所示，原来非常窄的灰阶分布(深黑色数据点)被拓宽到全灰阶分布(浅黑数据点)，且孔隙相和固体相的波峰分布特征被明显增强。孔隙与固体过渡带都得到很好的灰阶分布拓宽，从而完成了合理的图像对比度增强。
[0120] 图6a示出了Dolo-3-1岩样原始扫描图像，对CT扫描的原始图像进行剪裁，即选取边长为1000体素的正方体，从而可以得到如图6b所示切块后的原始图像，由于剪裁后的CT图像灰阶主要分布在非常低值区域，从图像的直观观测上显示为图像信息都比较暗。这样窄分布的灰阶直方图分布会给后续的图像分割造成困难，因此可以对图像进行对比度增强。对比度增强可通过如公式(2)实现，也就是通过灰度变换将原本窄分布的灰阶直方图变换成宽分布的直方图，拓宽孔隙相与固体相间的过渡带的灰阶分布，对比度增强后的结果如图6c所示，这样就能更好地方便后续图像分割阈值的选取。
[0121] 利用公式(2)～(4)所示的各项异性滤波方可以获得如图7a所示的结果。结合图7a和图7b所示，经各向异性滤波后，固体颗粒内部与孔隙内部的噪音信号被成功平滑掉，而固体颗粒与孔隙间的边缘信号得到了很好地保留。
[0122] 利用边缘检测法处理流程来进行边缘增强。由于固体颗粒与孔隙间的边缘一般对应于高梯度区域，因此可以使用梯度来提取边缘信号。通过图像梯度的累加求和函数可以确定梯度累加占总值90％的梯度点，高于这个梯度点就可以认为是高梯度指示相边界。结合图8a和图8b所示，可以看到该方法效果非常好，很完整地将固体颗粒与孔隙间的边缘区域给识别出来了，在此基础上进行图像分割可有利于降低分割不确定性，从而获得更好的图像分割效果。
[0123] 通过对Dolo3-1号岩样的一个边缘增强后的图像切片(如图9a所示，)进行二值化图像处理，可以得到如图9b所示的二值化图像。在图9b中黑色表示孔隙，灰白色表示固体颗粒。如此，通过对Dolo3-1号岩样各切片进行相同的二值化图像处理，最终可获得经二值化图像处理后的结果，也即如图10a所示的3D岩石骨架图，此时图10a展示的既有岩石固体相也有孔隙相。此时，在该步骤可以用实验室实测孔隙度与数字岩芯估算孔隙度进行相互验证，标定图像处理流程的合理性，如果超出误差则与重新返回图像处理流程进行参数调整优化等处理，重新实施图像处理步骤。从而可以得到如图10b所示的数字岩芯的孔隙结构示意图。
[0124] 本示例性实施例中，实验室氦气测得孔隙度为3.711％，而分割后得到成像孔隙度略小于氦气孔隙度，其值为2.297％。这是因为碳酸盐岩具有非常强的非均质性，其孔径可以小至几百个纳米而大至几个厘米。而且，图像分辨率与成像视域总有着不可协调的矛盾，若要使得分辨率越高则成像视域就会变小。所以一旦视域固定，那么其成像分辨率也是唯一不变的。对于成像孔隙度略小于实验孔隙度这样的现象一个较合理的解释是：碳酸盐岩中总有一些孔隙，其孔径甚至比CT分辨率还要小而无法成像，因此由CT成像估算的孔隙度会小于实验室氦气测量的孔隙度。由于微孔结构不能被完全分辨，CT图像上颗粒与微孔接触区域会表现出模糊状。这样的话分割过程就会产生不可避免的误差。而相对可信的分割过程通常会使分割后的孔隙度偏低。
[0125] 用完全相同的处理方法对剩余六块白云岩CT数据进行处理，并得到相应二值化后的三维数字岩芯模型。对这七块白云岩在实验室条件下利用氦气孔隙度渗透率测量仪对白云岩岩样的孔隙度进行了测试，并借此来验证所建立三维数字岩芯模型的准确性，测得数据见上表1。总体而言，可以看到这七块白云岩的数字岩芯计算孔隙度都要略低于实验室氦气孔隙度测量值。由于碳酸盐岩非均质性强，其孔径分布非常广(几百纳米至几厘米)。一般会有一部分尺寸小于成像分辨率的孔径在CT图像上清晰显示，所以造成数字岩芯成像孔隙度会略小于实验室氦气测量孔隙度。然而，实测孔隙度比较好地证明了数字岩芯图像建立过程中图像处理流程的合理性，另外，即便略有误差但对进一步进行规律性研究本项发明所提议的研究思路与流程均已足够。
[0126] 另外，本示例性实施例中，还可以讲薄片分析所获得的孔隙结构类型，与如图10b所示的孔隙网络模型进行比较，从而可以给出该岩样的孔隙结构大致类型。
[0127] (三)、读取所述数字岩芯的指定数据；所述指定数据包括孔隙度、纵波速度和密度。
[0128] (四)、基于预设的机器学习模型训练所述指定数据，以获取孔隙度与纵波速度、密度的第二关系模型。
[0129] 例如，在本说明书一实施例中，可以数字岩芯的纵波速度和密度作为输入，以数字岩芯的孔隙度作为目标输出，训练预设的神经网络模型，从而可以得到孔隙度与纵波速度、密度的关系曲线，即确定孔隙度与纵波速度、密度的第二关系模型。
[0130] 由于碳酸盐岩孔隙结构的复杂性，尤其次生孔隙的影响，其孔隙度和速度以及密度之间的关系也会变的复杂。因此仅基于地震叠前反演预测的孔隙度会变得非常不可靠。而测井数据或数字岩芯数据中的孔隙度信息是已知的，并且也较为准确。在本说明书一些实施例中，为了使预测的孔隙度准确，可以从这两种数据中来寻找孔隙度与纵波速度、密度之间的关系。因此，通过将第一关系模型和第二关系模型的加权和作为第三关系模型，可以有利于进一步提高碳酸盐岩储层孔隙结构预测的准确性。其中，在将第一关系模型和第二关系模型进行加权求和时，第一关系模型和第二关系模型的权重系数可以根据需要设定。
例如，在本说明书一示例性实施例中，第一关系模型和第二关系模型的权重系数可以各为
0.5。
[0131] 在本说明书一实施例中，所述基于所述测井数据获取所述目标区块的第一孔隙结构因子可以根据公式 获得。
[0132] 其中，ρ为目标区块基于测井数据的密度，ρm为目标区块基于测井数据的孔隙相中固体相矿物骨架的密度，ρf为目标区块基于测井数据的孔隙相中充填的流体或气体的密度，φ为目标区块基于测井数据的孔隙度，K为目标区块基于测井数据的体积模量，VP为目标区块基于测井数据在指定频段下随围压变化的纵波速度，VS为目标区块基于测井数据在指定频段下随围压变化的横波速度，μ为目标区块基于测井数据的剪切模量，Fk为中间变量，Km为目标区块基于测井数据的孔隙相中固体相矿物骨架的体积模量，Kf为目标区块基于测井数据的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为目标区块基于测井数据的孔隙结构因子。
[0133] 在本说明书另一实施例中，所述基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子也可以根据上述计算公式 获得。其中，ρ为目标区块基于数字岩芯的密度，ρm为目标区块基于数字岩芯的孔隙相中固体相矿物骨架的密度，ρf为目标区块基于数字岩芯的孔隙相中充填的流体或气体的密度，φ为目标区块基于数字岩芯的孔隙度，K为目标区块基于数字岩芯的体积模量，VP为目标区块基于数字岩芯在指定频段下随围压变化的纵波速度，VS为目标区块基于数字岩芯在指定频段下随围压变化的横波速度，μ为目标区块基于数字岩芯的剪切模量，Fk为中间变量，Km为目标区块基于数字岩芯的孔隙相中固体相矿物骨架的体积模量，Kf为目标区块基于数字岩芯的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为目标区块基于数字岩芯的孔隙结构因子。
[0134] 在本说明书一些实施例中，上述基于数字岩芯数据获取目标区块的第二孔隙结构因子中，数字岩芯可以为网格化的数字岩芯，且网格化的数字岩芯体对应的弹性性质，可以通过有限单元法(Finite Element Method，简称FEM)弹性性质模拟器模拟得到，具体如下：
[0135] 网格化之前，数字岩芯作为3D数据模型就类似于一个魔方的正方体数据，而网格化就类似于把魔方分割成了很多构成魔方的小立方体网格，而每个小立方体网格也是一个小的3D数字岩芯数据体(也可以称为子网格或数字岩芯体)。
[0136] 经过网格化后，网格化的数字岩芯比较容易地获得各子网格所代表的数字岩芯体的孔隙度。由于网格化的数字岩芯中的子网格数据巨大，这样可以获得一个孔隙度跨度非常宽的子网格集合。在取芯异常困通常只能获得有限数量岩芯的情况下，如果直接依据有限数量岩芯预测孔隙度，则孔隙度分布通常非常窄。而在本说明书实施例中，由于基于网格化的数字岩芯可以获得一个孔隙度跨度非常宽的子网格集合，从而有利于获得更为全面的孔隙度分布。
[0137] 在本说明书一实施例中，每个子网格可以根据CT扫描分辨率用六面体单元进行网格剖分，每个子网格(即体素)等于CT扫描分辨率。此外，还可根据上述步骤获得的XRD矿物分析分析结果，给每个体素的矿物骨架模量赋值。
[0138] 例如，在说明书一示例性实施例中，图11示出了某个岩样的网格化的3D数字岩芯，通过网格化，每块岩样可以得到125个子网格数据集，这样对于本示例性实施例中的7块岩样则可以得到125*7＝875个子网格数据集。对于875个子网格中的每个子图像可以容易获取其孔隙度(即利用数字岩芯估算的孔隙度)。由于子网格数据比较大，这样可以获得一个孔隙度跨度非常宽的子网格集合。此外，还可以根据XRD矿物分析结果，给每个子网格对应的数字岩芯体的进行矿物骨架模量赋值。
[0139] 在本说明书一实施例中，所述对网格化的数字岩芯进行有限元弹性性质模拟，可以包括如下步骤：
[0140] (1)、确定子网格骨架参数。由于步骤S103中已为每个体素的矿物骨架模量赋值。据此可以确定每个子网格的骨架参数。
[0141] (2)、接收输入的矿物骨架矿物参数。
[0142] (3)、接收输入的流体模量。
[0143] (4)、利用有限单元法(Finite Element Method，简称FEM)弹性性质模拟器(即有限元算法)，可以模拟获得每一个子网格体的弹性性质。
[0144] 在本说明书一实施例中，FEM弹性性质模拟器例如可以是Garboczi(1995)提出的线弹性有限元法：将三维数字岩芯模型中的每个体素视作一个线弹性有限元网格，这样就免去了繁琐的网格剖分过程(本质上这个过程等同于用标准的六面体网格剖分求解区域)。
对模型施加周期边界条件，便可将弹性位移分布的求解问题转化为最小位能原理求解的物理问题，数学上求解就是有限元法中典型的泛函求极值问题。例如，Arns(2002)提到可以使用快速共轭梯度法来对泛函求导方程进行迭代收敛以逼近系统最小能量，最后求得数值解即为六个方向的应力σxx,σyy,σzz,σxz,σyz,σxy与六个方向的应变εxx,εyy,εzz,εxz,εyz,εxy，进而求得等效弹性模量。
[0145] 例如，在说明书一示例性实施例中，对获取的875个子网格所代表的每一个数字岩芯体，可以利用FEM弹性模拟器对每一个数据体进行有限元数值模拟得到了子网格数据体代表的数字岩芯体的弹性性质。与此同时，还可以定量地计算得到每个子网格数据体代表的数字岩芯体的孔隙结构参数γ。对分割所得的大量子网格数据体，对于每一个子网格可以计算得到其孔隙度φ、弹性性质(比如纵、横波速度等)和孔隙结构因子γ。对于本示例性实施例而言，可获得如图12所示的碳酸盐岩岩样在分类下的弹性性质与孔隙度关系。
[0146] 在本说明书一些实施例中，通过将第一孔隙结构因子和第二孔隙结构因子的加权和作为第三孔隙结构因子，也可以有利于进一步提高碳酸盐岩储层孔隙结构预测的准确性。其中，在将第一孔隙结构因子和第二孔隙结构因子进行加权求和时，第一孔隙结构因子和第二孔隙结构因子的权重系数可以根据需要设定。例如，在本说明书一示例性实施例中，第一孔隙结构因子和第二孔隙结构因子的权重系数可以各为0.5。
[0147] 在本说明书一示例性实施例中，将如图16所示的纵波速度数据体以及如图17所示的密度数据体输入到第三关系模型中，可以获得如图18所示目标区块的孔隙度数据体。
[0148] 在本说明书一些实施例中，所述将纵波速度数据体与孔隙度数据体进行交会，获得目标区块的孔隙结构数据体，是指利用交会图法将纵波速度数据体与孔隙度数据体在平面图上进行交会，从而根据交会点的坐标定出孔隙结构的数值或范围。
[0149] 如图19所示，在本说明书一示例性实施例中，将目标区块的孔隙度数据体与基于地震叠前反演的纵波速度数据体进行交会，可以预测出不同交会点对应的孔隙结构类型，进而可以得到对应的孔隙结构数据体(例如图20所示)。而将三维纵波速度数据体的所有测线的孔隙结构都计算出，就可以得到整个区块的孔隙结构数据体。
[0150] 参考图21所示，与上述碳酸盐岩储层孔隙结构预测方法对应，本说明书一些实施例的碳酸盐岩储层孔隙结构预测装置可以包括：
[0151] 第一数据体获取模块211，可以用于获取目标区块基于地震叠前反演的纵波速度数据体和密度数据体；
[0152] 第一模型获取模块212，可以用于根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，并根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型；
[0153] 第二模型获取模块213，可以用于将所述第一关系模型和所述第二关系模型的加权和作为第三关系模型；
[0154] 第一因子获取模块214，可以用于基于所述测井数据获取所述目标区块的第一孔隙结构因子，并基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子；
[0155] 第二因子获取模块215，可以用于将所述第一孔隙结构因子和所述第二孔隙结构因子的加权和作为第三孔隙结构因子；
[0156] 第二数据体获取模块216，可以用于将所述纵波速度数据体和所述密度数据体输入所述第三关系模型，获得所述目标区块的孔隙度数据体；
[0157] 孔隙结构获取模块217，可以用于将所述纵波速度数据体与所述孔隙度数据体进行交会，获得所述目标区块的孔隙结构数据体。
[0158] 在本说明书一些实施例中的碳酸盐岩储层孔隙结构预测装置中，所述根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，可以包括：
[0159] 获取所述目标区块的测井数据；所述测井数据包括孔隙度、纵波速度和密度；
[0160] 基于预设的机器学习模型训练所述测井数据，以获取孔隙度与纵波速度、密度的第一关系模型。
[0161] 在本说明书一些实施例中的碳酸盐岩储层孔隙结构预测装置中，所述根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型，包括：
[0162] 获取所述目标区块的指定岩样的岩石物理参数；
[0163] 基于所述岩石物理参数及所述指定岩样的CT扫描图像构建数字岩芯；
[0164] 读取所述数字岩芯的指定数据；所述指定数据包括孔隙度、纵波速度和密度；
[0165] 基于预设的机器学习模型训练所述指定数据，以获取孔隙度与纵波速度、密度的第二关系模型。
[0166] 在本说明书一些实施例中的碳酸盐岩储层孔隙结构预测装置中，所述基于所述测井数据获取所述目标区块的第一孔隙结构因子，包括：
[0167] 根据公式 确定所述目标区块的第二孔隙结构因子；
[0168] 其中，ρ为目标区块基于测井数据的密度，ρm为目标区块基于测井数据的孔隙相中固体相矿物骨架的密度，ρf为目标区块基于测井数据的孔隙相中充填的流体或气体的密度，φ为目标区块基于测井数据的孔隙度，K为目标区块基于测井数据的体积模量，VP为目标区块基于测井数据在指定频段下随围压变化的纵波速度，VS为目标区块基于测井数据在指定频段下随围压变化的横波速度，μ为目标区块基于测井数据的剪切模量，Fk为中间变量，Km为目标区块基于测井数据的孔隙相中固体相矿物骨架的体积模量，Kf为目标区块基于测井数据的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为目标区块基于测井数据的孔隙结构因子。
[0169] 在本说明书一些实施例中的碳酸盐岩储层孔隙结构预测装置中，所述基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子，包括：
[0170] 根据公式 确定所述目标区块的第二孔隙结构因子；
[0171] 其中，ρ为目标区块基于数字岩芯的密度，ρm为目标区块基于数字岩芯的孔隙相中固体相矿物骨架的密度，ρf为目标区块基于数字岩芯的孔隙相中充填的流体或气体的密度，φ为目标区块基于数字岩芯的孔隙度，K为目标区块基于数字岩芯的体积模量，VP为目标区块基于数字岩芯在指定频段下随围压变化的纵波速度，VS为目标区块基于数字岩芯在指定频段下随围压变化的横波速度，μ为目标区块基于数字岩芯的剪切模量，Fk为中间变量，Km为目标区块基于数字岩芯的孔隙相中固体相矿物骨架的体积模量，Kf为目标区块基于数字岩芯的孔隙相中充填的流体或气体的体积模量，f为中间变量，γ为目标区块基于数字岩芯的孔隙结构因子。
[0172] 为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本说明书时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
[0173] 参考图22所示，与上述碳酸盐岩储层孔隙结构预测方法对应，本说明书一些实施例的计算机设备可以包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时执行如下步骤：
[0174] 获取目标区块基于地震叠前反演的纵波速度数据体和密度数据体；
[0175] 根据所述目标区块的测井数据确定孔隙度与纵波速度、密度的第一关系模型，并根据所述目标区块的数字岩芯数据确定孔隙度与纵波速度、密度的第二关系模型；
[0176] 将所述第一关系模型和所述第二关系模型的加权和作为第三关系模型；
[0177] 基于所述测井数据获取所述目标区块的第一孔隙结构因子，并基于所述数字岩芯数据获取所述目标区块的第二孔隙结构因子；
[0178] 将所述第一孔隙结构因子和所述第二孔隙结构因子的加权和作为第三孔隙结构因子；
[0179] 将所述纵波速度数据体和所述密度数据体输入所述第三关系模型，获得所述目标区块的孔隙度数据体；
[0180] 将所述纵波速度数据体与所述孔隙度数据体进行交会，获得所述目标区块的孔隙结构数据体。
[0181] 虽然上文描述的过程流程包括以特定顺序出现的多个操作，但是，应当清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行(例如使用并行处理器或多线程环境)。
[0182] 本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0183] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0184] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0185] 在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。
[0186] 内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。
[0187] 计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁盘式存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。
[0188] 还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法或者设备中还存在另外的相同要素。
[0189] 本领域技术人员应明白，本说明书的实施例可提供为方法、系统或计算机程序产品。因此，本说明书可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
[0190] 本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
[0191] 本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
[0192] 以上所述仅为本说明书的实施例而已，并不用于限制本说明书。对于本领域技术人员来说，本说明书可以有各种更改和变化。凡在本说明书的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本说明书的权利要求范围之内。
