基于数字岩心构建碳酸盐岩多重孔裂隙结构模型的方法

技术领域

本发明涉及数字岩石物理技术及致密油气储层岩石物理模型技术领域，尤其涉及基于数字岩心构建碳酸盐岩多重孔裂隙结构模型的方法。

背景技术

碳酸盐岩的储集空间，通常分为原生孔隙、溶洞和裂缝三类。与砂岩储集层相比，碳酸盐储集层储集空间类型多、次生变化大，具有更大的复杂性和多样性。

岩石物理模型通过等效介质理论来模拟实际储层性质，通过等效方法逼近储层岩石的实际情况，对地下储层进行表征与预测。Golikov等(2012)通过对比北海测井资料中斜井与垂直井的速度，对古新统页岩油藏的各向异性进行了定量分析，建立了各向异性的岩石物理模型，解释储层流体与孔隙度的变化。Pang等(2019，2020)建立了多尺度等效岩石物理模型，并利用实验及地震数据的验证，定量预测碳酸盐岩储层中的流体饱和度及裂隙含量。常规的岩石物理模型在构建模型时，通常假设一种或两种固定形状的孔隙类型来模拟岩石内部的孔隙结构(Gupta et al.,2012；Tan et al.,2020)，然而，对于深层碳酸盐岩储层来说，储层微观结构是复杂多变的，基于等效假设的理论模型与实际地下真实结构不一致，导致所建立的模型可能与实际情况不符，因此，针对深层碳酸盐岩储层岩石复杂的孔隙结构，需要提出更精确的岩石物理建模方法。

发明内容

本发明所要解决的技术问题：提供基于数字岩心构建碳酸盐岩多重孔裂隙结构模型的方法，用以克服常规岩石物理建模基于等效假设的孔隙结构参数的限制，模拟岩石真实孔隙空间情况下孔隙结构及流体对弹性波的影响作用，实现数字岩心技术和理论岩石物理模型两种技术方法的结合。

本发明为解决以上技术问题而采用以下技术方案：

本发明提出的基于数字岩心构建碳酸盐岩多重孔裂隙结构模型的方法，包括如下步骤：

S1、选定碳酸盐岩储层岩石作为岩石样本，针对岩石样本开展CT扫描、铸体薄片鉴定、矿物分析和超声波测试，分析岩石样本通过CT扫描获得的三维灰度图像和特征，特征包括弹性波速度、孔隙度和衰减。

S2、基于CT扫描得到的岩石样本的三维灰度扫描图像，进行阈值分割，构建碳酸盐岩数字岩心，提取岩石样本的体积分数和孔隙纵横比数据。

S3、根据岩石样本，利用超声波测试得到弹性波速度，利用裂隙参数反演方法，计算岩石裂隙孔隙度。

S4、基于步骤S2中提取的结果和步骤S3中得到的裂隙孔隙度，结合微分有效介质理论及双重孔隙结构波传播方程，构建碳酸盐岩多重孔裂隙结构模型。

S5、基于碳酸盐岩多重孔裂隙结构模型，模拟岩石真实孔隙空间情况下的波响应特征，并与岩石样本及实际地层数据进行对比验证。

进一步的，步骤S1中，分析岩石样本通过CT扫描获得的三维灰度图像和特征包括以下内容：

对岩石样本开展铸体薄片鉴定和X射线计算机断层CT扫描，分别得到铸体薄片和三维灰度扫描图像，并对其进行观察分析，得到碳酸盐岩储集空间，该空间包括粒间孔隙、溶孔、微裂隙；通过孔渗测量仪测量出岩石样本的孔隙度和渗透率。

利用X射线衍射开展矿物分析，得到岩石样本的矿物组分。

使用超声波实验仪器开展变压力下的超声波测试，记录通过岩石样本的超声波波形，通过提取波形的初至时间得到弹性波速度；利用超声波谱比法估算岩石样本的品质因子，以品质因子在10万以上的标准铝块作为参考材料，计算岩石样本的衰减，具体公式为：

其中，f表示频率，A

将孔渗测量结果和超声波测试的结果进行交汇，发现弹性波速度及衰减随有效压力呈单调变化的趋势，速度随着压力的增大而增大，衰减降低，并且饱水时的衰减高于干燥状态下的岩石。

进一步的，步骤S2中，提取岩石多重孔裂隙纵横比及体积分数包括以下内容：

基于CT扫描获取的三维灰度扫描图像，利用中值滤波算法去除CT扫描过程中的背景噪声，提高灰度图像显示质量及对比度，增强岩石骨架与孔隙空间的边界过渡；根据矿物和孔隙灰度的对比，选择阈值进行图像分割完成矿物和孔隙的划分，使三维灰度扫描图像中的全体像素仅代表孔隙空间和矿物，获得岩石样本的三维数据体，得到碳酸盐岩数字岩心。

基于构建的碳酸盐岩数字岩心，标记所有连接的孔隙，统计孔隙面积所占的比例，得到各个孔隙的体积分数，利用最大内切圆和最小外接圆算法，分别计算孔隙的长短轴半径，获得孔隙纵横比数据。

进一步的，步骤S3中，计算岩石裂隙孔隙度包括以下子步骤：

S301、基于通过超声波测试得到的岩石样本数据，即实验室所测量的干燥时高有效压力下岩石的纵横波速度，估算岩石高压下的等效体积模量和剪切模量，具体公式为：

其中，K

S302、以包含硬孔隙的介质为背景相，考虑裂隙对岩石弹性性质的影响，基于Mori-Tanaka理论建立岩石弹性模量与微观孔隙结构的关系，得到岩石介质的等效弹性模量，具体内容为：

其中，K

S303、基于各有效压力下的累积裂隙密度，建立裂隙密度随有效压力变化的定量关系，具体公式为：

其中，

S304、根据裂隙孔隙度与裂隙密度的关系，计算岩石的裂隙孔隙度，具体公式为：

其中，φ

进一步的，步骤S4中，构建碳酸盐岩多重孔裂隙结构模型包括以下内容：

基于岩石裂隙孔隙度，将数字岩心提取的孔隙空间划分为具有多种纵横比的孔隙及裂隙，设定岩石骨架包括两种不同孔隙的介质，背景相包括硬孔隙，包体包括软孔隙，硬孔隙和软孔隙具有不同的纵横比，硬孔隙纵横比大于软孔隙纵横比。

利用微分有效介质理论将数字岩心提取的孔隙和裂隙分别加入背景相和包体，将包体加入到背景相中，模拟岩石骨架的等效弹性模量，具体公式为：

其中，K

针对深层碳酸盐岩储层岩石内部混合流体，利用双重孔隙结构波传播方程，模拟局部流体流动造成的频散和衰减效应，进行流体替换，估算含流体及复杂孔隙结构的饱和岩石波响应特征，构建碳酸盐岩多重孔裂隙结构岩石物理模型。

进一步的，步骤S5中，模型的验证包括以下内容：

将碳酸盐岩多重孔裂隙结构模型的模拟结果与岩石样本的数据进行对比，并求取各个压力下的数据平均值来验证岩石的速度频散及衰减模拟结果。对比结果表明，模型可以较好的模拟出流体在孔隙结构中造成的速度频散及衰减；提取实际地层数据进一步验证碳酸盐岩多重孔裂隙结构模型，将样本数据的平均值及模型结果放置在取样深度处。对比实验及测井数据和模型拟合结果，该模型与岩石样本及地层数据整体吻合较好。研究结果表明，模拟结构与实验及地层数据吻合较好，构建的模型可以有效地模拟岩石波响应特征。

进一步的，本发明还提出了一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现前文所述的基于数字岩心构建碳酸盐岩多重孔裂隙结构模型的方法的步骤。

进一步的，本发明还提出了一种计算机可读的存储介质，所述计算机可读的存储介质存储有计算机程序，所述计算机程序被处理器运行时执行前文所述的基于数字岩心构建碳酸盐岩多重孔裂隙结构模型的方法。

本发明采用以上技术方案，与现有技术相比，其显著技术效果如下：

1.本发明基于X射线CT扫描、铸体薄片鉴定技术、矿物分析及超声波室内实验，综合分析储层岩石岩性、孔隙结构、弹性及衰减特征，揭示了深层碳酸盐岩储层岩石的孔隙结构、流体及波响应特征规律。

2.本发明基于CT扫描及铸体薄片，构建碳酸盐岩数字岩心，在微观微纳米尺度精细表征岩石内部孔喉几何尺寸、拓扑结构等特征，结合裂隙反演方法和实验数据计算样本的裂隙结构参数，精细刻画了岩石微观孔隙结构。

3.本发明利用数字岩心技术在岩石微观结构精确表征方面的优势，获得了更多的岩石微观结构信息，弥补理论岩石物理模型的不足，从而建立起更精确的储层岩石物理模型，结合数字岩心技术和理论岩石物理模型两种技术，可以较好的模拟岩石波响应特征，有效地应用到实际碳酸盐岩储层。

附图说明

图1为本发明的整体流程示意图。

图2为本发明实施例中岩石样本的铸体薄片图。

图3为本发明实施例中岩石样本的三维灰度扫描图像。

图4为本发明实施例中超声波实验装置图。

图5为本发明实施例中岩石样本的纵、横波速度分别在干燥和饱和水情况下随有效压力的变化情况图。

图6为本发明实施例中岩石样本的衰减分别在干燥和饱和水情况下随有效压力的变化情况图。

图7为本发明实施例中岩石样本的三维灰度扫描图像和三维空间孔隙结构图与分布图。

图8为本发明实施例中岩石样本的孔裂隙纵横比计算流程图。

图9为本发明实施例中岩石样本的孔裂隙纵横比提取结果图。

图10为本发明实施例中岩石样本的裂隙孔隙度计算结果图。

图11为本发明实施例中岩石样本的纵横波速度和衰减随裂隙孔隙度变化情况图。

图12为本发明的碳酸盐岩多重孔裂隙结构模型构建流程图。

图13为本发明实施例中岩石样本的纵波速度频散及衰减响应特征图。

图14为本发明实施例中模拟结果与岩石样本数据的对比图。

图15为本发明实施例中模拟结果与岩石样本数据及测井数据对比结果图。

具体实施方式

以下结合附图对本发明公开的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本公开，并不用于限制本公开。

本发明提出了基于数字岩心构建碳酸盐岩多重孔裂隙结构模型的方法，如图1所示，包括如下步骤：

S1、选定四川盆地G区龙王庙组深层碳酸盐岩储层岩石作为岩石样本，发育不同类型的孔隙与裂隙，孔隙内常含有不相混溶流体，储层埋深超过4.5千米，地层的温度约为120℃，上覆地层压力约为110MPa，气藏的孔隙压力约为65～70MPa。储层岩石孔隙度较低，岩性为较纯净的白云岩，含极少量的方解石及泥质，孔隙空间主要为粒间孔隙和裂隙，孔隙流体主要为气体和水。

从岩石样本中提取2块碳酸盐岩岩石样本(DS3及DS17)，样本采集深度分别为4627.15m及4625.6m，将样本加工为直径38mm和长度50mm左右的圆柱体，孔隙度分别为1.96％和3.46％，渗透率分别为0.16mD和0.001mD，密度分别为2.802g/cm

针对岩石样本开展CT扫描、铸体薄片鉴定、矿物分析和超声波测试，分析岩石样本通过CT扫描获得的三维灰度图像和特征，特征包括弹性波速度、孔隙度和衰减，具体内容为：

如图2所示，对岩石样本开展铸体薄片鉴定，得到铸体薄片，图2的a)为岩石样本DS3的铸体薄片图，图2的b)为岩石样本DS17的铸体薄片图，可以看到样本DS3发育有裂隙，样本DS17发育孔隙和溶孔相对较多。利用X射线衍射开展矿物分析，得到岩石样本的矿物组分，岩石为纯净的白云岩，白云石矿物含量在98％以上，含有少量的方解石与泥质。三维灰度图像和铸体薄片鉴定结果与孔隙度测量结果表现一致。

如图3所示，对岩石样本开展X射线计算机断层CT扫描，得到三维灰度扫描图像，图3的a)为岩石样本DS3的三维灰度扫描图像，图3的b)为岩石样本DS17的三维灰度扫描图像，分别包含1980和1954张二维灰度图像，图像分别为1398×1404、1408×1402像素，扫描分辨率(每个像素点对应的物理尺寸)为27.6μm。图像是由一定数目像素组成的灰阶图像，灰度值为0～65535，灰度值较高的像素反映的是岩石骨架、矿物；灰度值较低的像素反映的是孔喉空间。可以看到，岩石样本DS3相较于岩石样本DS17更为致密，孔隙含量较低，但发育有裂隙，扫描图像与样本铸体薄片及孔隙度测量结果一致。

如图4所示，使用超声波实验仪器开展变压力下的超声波测试，记录通过岩石样本的超声波波形，通过提取波形的初至时间得到弹性波速度。进行干燥状态下的实验，将岩石样本烘干后置入实验装置中，在围压5-60MPa、温度为26.5℃及大气压下，进行超声波实验测量(1MHz)，记录通过岩石样本的超声波波形；对于饱水实验观测，用加压饱和的方法使样本饱水，设置围压为25-60MPa、孔隙压力为20MPa、温度为27℃，进行超声波实验测量，记录波形；制备与岩石样本相同尺寸的铝块作为参考材料，通过超声波脉冲法测量铝块中的波形数据；利用超声波谱比法估算岩石样本的品质因子，以品质因子在10万以上的标准铝块作为参考材料，计算岩石样本的衰减，具体公式为：

其中，f表示频率，A

图5的a)为岩石样本的纵、横波速度在干燥情况下(Dry)随有效压力的变化情况，图5的b)为岩石样本的纵、横波速度在饱和水情况下(Wet)随有效压力的变化情况，结合图6中岩石样本的衰减在变压力下的实验结果图，图中散点形状代表不同的样本，V

干燥情况下，在低压时，样本DS17的速度低于样本DS3，而随着压力的增大，岩石内部的微裂隙闭合，样本DS17的速度超过样本DS3。可知，样本DS17随着压力的增加，内部闭合的裂隙较多，导致速度迅速增大。饱和状态下，样本DS17的横波速度低于样本DS3，而对于纵波速度，样本DS17高于样本DS3。这是由于流体在孔隙中在波的激励下会产生局部流体流动，产生高频频散，高频时速度增大，这种现象在孔隙度较大时表现更为明显，这与图6中样本DS17的衰减更高以及饱水与干燥间差值更高表现一致。

S2、基于通过CT扫描得到的岩石样本的三维灰度扫描图像，进行阈值分割，构建碳酸盐岩数字岩心，提取岩石样本的体积分数和孔隙纵横比数据，具体内容为：

为方便处理，从岩心三维灰度图像中选取一块长方体区域，进行图像分割得到数字岩心，如图7所示，给出了岩石样本和对应的三维灰度扫描图像以及三维空间孔隙结构与分布。基于CT扫描获取的三维灰度扫描图像，利用中值滤波算法去除CT扫描过程中背景噪声在灰度图像中形成的干扰，提高灰度图像显示质量及对比度，增强岩石骨架与孔隙空间的边界过渡；根据矿物和孔隙灰度的对比，手动选择阈值进行图像分割完成矿物和孔隙的划分，使三维灰度扫描图像中的全体像素只代表孔隙空间和矿物，获得岩石样本的三维数据体，得到碳酸盐岩数字岩心。

基于构建的碳酸盐岩数字岩心，定量表征岩心微观孔隙结构，计算岩石内部孔裂隙纵横比，如图8所示，具体内容为：

标记所有连接的孔隙，统计孔隙面积所占的比例，得到各个孔隙的体积分数，利用最大内切圆和最小外接圆算法，分别计算孔隙的长短轴半径，获得孔隙纵横比数据，最终结果如图9所示。图9的a)为岩石样本DS3的孔裂隙纵横比提取结果图，图9的b)为岩石样本DS17的孔裂隙纵横比提取结果图，可以看到，样本的孔裂隙纵横的分布大致符合高斯分布，并且，样本DS3小纵横比下的孔裂隙数量较多，密度大。

S3、根据岩石样本，利用超声波测试得到弹性波速度，利用裂隙参数反演方法，计算岩石裂隙孔隙度，具体内容为：

S301、基于通过超声波测试得到的岩石样本数据，即实验室所测量的干燥时高有效压力下岩石的纵横波速度，高压下的等效体积模量和剪切模量，具体公式为：

其中，K

利用最小二乘法求取最吻合的硬孔隙的纵横比α

S302、以包含硬孔隙的介质为背景相，考虑裂隙对岩石弹性性质的影响，基于Mori-Tanaka理论建立岩石弹性模量与微观孔隙结构的关系，得到岩石介质的等效弹性模量，具体内容为：

其中，K

基于各个有效压力下测量的纵横波速度，利用最小二乘法拟合得到各个压力下的累积裂隙密度Γ

S303、基于各有效压力下的累积裂隙密度，建立裂隙密度随有效压力变化的定量关系，具体公式为：

其中，

S304、根据裂隙孔隙度与裂隙密度的关系，计算岩石的裂隙孔隙度，具体公式为：

其中，φ

图10为岩石样本的裂隙孔隙度结果图，图11的a)为纵波速度随裂隙孔隙度的变化情况，图11的b)为横波速度随裂隙孔隙度的变化情况，图11的c)为衰减随裂隙孔隙度的变化情况，从图10和图11可以看到，样本的裂隙孔隙度随着压力的增加而降低，在低压时样本DS17裂隙含量高于样本DS3，而随着压力的增加，样本DS3的裂隙含量逐渐高于样本DS17，显示出随着压力的增大，样本DS17闭合的裂隙更多，这与图5中样本的速度结果表现一致。另外，纵横波速度随着裂隙的增加而降低，衰减增大，并且，干燥时的纵波速度低于饱水时的速度，而干燥状态下的横波速度高于饱水状态。

S4、基于步骤S2中提取的结果和步骤S3中得到的裂隙孔隙度，结合微分有效介质理论及双重孔隙结构波传播方程，构建碳酸盐岩多重孔裂隙结构模型，如图12所示，具体内容为：

基于岩石裂隙孔隙度，将数字岩心提取的孔隙空间划分为具有多种纵横比的孔隙及裂隙，设定岩石骨架包括两种不同孔隙的介质，背景相包括硬孔隙，包体包括软孔隙，硬孔隙和软孔隙具有不同的纵横比，硬孔隙纵横比大于软孔隙纵横比；利用微分有效介质理论将数字岩心提取的孔隙和裂隙分别加入背景相和包体，将包体加入到背景相中，模拟岩石骨架的等效弹性模量，具体公式为：

其中，K

针对深层碳酸盐岩储层岩石内部混合流体，利用双重孔隙结构波传播方程，模拟局部流体流动造成的频散和衰减效应，进行流体替换，估算含流体及复杂孔隙结构的饱和岩石波响应特征，构建碳酸盐岩多重孔裂隙结构岩石物理模型。

通过平面波分析，给出岩石的纵波速度及衰减，具体公式为：

V

其中，v表示复波速，v＝ω/k，ω表示角频率，k表示复波数；V

S5、基于碳酸盐岩多重孔裂隙结构模型，模拟岩石真实孔隙空间情况下的波响应特征，并与岩石样本及实际地层数据进行对比验证，具体内容为：

岩石样本参数为：基质(白云石)的体积和剪切模量及密度分别为94.9GPa、45GPa及2.87g/cm

如图13所示，基于碳酸盐岩多重孔裂隙结构模型，模拟了在不同流体类型情况下，纵波速度频散和衰减随频率的变化规律，其中图13的a)为纵波速度频散随频率的变化规律，图13的b)为衰减随频率的变化规律，可以看到，含水状态下的速度、频散及衰减明显高于干燥状态，相较于含气状态，含水时主频向低频移动。

如图14所示，基于碳酸盐岩多重孔裂隙结构模型，模拟岩石饱和气及饱和水下的纵波速度和衰减，并且求取各个压力下的数据平均值来验证岩石的速度频散及衰减模拟结果，其中图14的a)为纵波速度的模拟结果与岩石样本对比情况，图14的b)为衰减的模拟结果与岩石样本对比情况，可以看到，饱和水状态下的纵波速度和衰减明显高于干燥状态，尤其是孔隙度较大的样本DS17两种饱和状态下的频散及衰减更大。结果表明，模型可以较好的模拟出流体在孔隙结构中造成的速度频散及衰减。

在工区中选取一口产气井A，该井为两块样本的采集井，如图15所示，曲线为测井数据，将样本数据的平均值及模型结果放置在取样深度处，其中，黑虚线框为储层高产气段。图15的a)为储层的孔隙度及含气饱和度，图15的b)及图15的c)为别为样本DS3及DS17实验及模型与测井数据对比图。对比显示，模型与实验及地层数据整体吻合较好，另外，除饱和水时的样本DS17之外，岩石样本的实验及模型弹性波速度小于测井数据，这是由于地层埋深较深，储层原位条件的压力较高(上覆压力约为110MPa)，导致速度较大。而饱水下的样本DS17实验和模型的纵波速度高于测井数据，这是由于高孔隙度下流体造成的速度频散较大，导致速度较高。

因此，本发明提出的碳酸盐岩多重孔裂隙结构模型可以较好的模拟岩石波响应特征，模拟结果与实验及地层数据吻合较好。结合数字岩心技术和理论岩石物理模型，利用二者的技术优势，可以较好的模拟岩石波响应特征，有效地应用到实际碳酸盐岩储层。

本发明实施例还提出一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序。需要说明的是，所述处理器执行所述计算机程序时对应本发明实施例所提供的方法的具体步骤，具备执行方法相应的功能模块和有益效果。未在本实施例中详尽描述的技术细节，可参见本发明实施例所提供的方法。

本发明实施例还提出一种计算机可读的存储介质，所述计算机可读的存储介质存储有计算机程序。需要说明的是，所述计算机程序被处理器运行时对应本发明实施例所提供的方法的具体步骤，具备执行方法相应的功能模块和有益效果。未在本实施例中详尽描述的技术细节，可参见本发明实施例所提供的方法。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。一种考虑孔隙连通性的致密砂岩岩石弹性参数预测方法。