基于数字岩心的古龙页岩油储层含油饱和度建模方法及其应用

2022-05-17闫伟林章华兵郑建东王文娟张兆谦朱建华

大庆石油地质与开发 2022年3期

闫伟林章华兵郑建东王文娟张兆谦朱建华

（1.中国石油大庆油田有限责任公司勘探开发研究院，黑龙江大庆 163712；2.黑龙江省陆相页岩油重点实验室，黑龙江大庆 163712）

0 引言

近年来，随着勘探开发技术的发展，页岩油已经成为各油田接替资源的重要勘探领域。松辽盆地青山口组、嫩江组是页岩油主要赋存层系，其中古龙页岩油是目前勘探开发的重点。目前大庆油田在青山口组获得勘探重大战略突破，多口水平井获得了工业油气流，实现了陆相页岩从“生”油到“产”油的历史性跨越。松辽盆地古龙页岩油具有巨大的资源量，是大庆油田重要的战略接替资源［1-2］。

古龙页岩是松辽盆地陆相湖盆内富含有机质的细粒沉积岩，以粒径小于3.9 μm 的泥级页岩为主，其中夹持的粉砂岩、碳酸盐岩等其他岩石累计厚度占总厚度的比例小于10%，通过技术改造可获得工业油气流的页岩［3］。前人在古龙页岩油的地质条件、富集机理、参数评价等研究上已经取得了很多成果［4-17］。

含油饱和度是储量评价和开发方案设计中非常重要的一个参数。对于传统砂岩或致密储层，可以通过常规岩石物理实验获得岩电参数，然后利用阿尔奇公式或其他饱和度公式计算含油饱和度参数［18-20］。但古龙页岩为陆相深水沉积，与常规储层相比具有高黏土含量、纹层页理和纳米孔发育、孔隙结构复杂、实验难度大等特点，在油气运移、赋存方式、储集空间类型等方面具有较大差异，应用常规储层饱和度评价方法计算饱和度精度较低。

前人在古龙页岩油饱和度方面目前仅在核磁共振技术确定含油饱和度方面做了相关研究［4］，岩石物理实验确定古龙页岩储层含油饱和度参数方面未有新的进展，对于无核磁资料的井不能满足生产需要。如何利用通过高精度岩石物理实验确定的岩电参数来精确地计算页岩油储层含油饱和度，认识富集层段空间分布，评价资源潜力，是页岩油勘探取得重大突破的关键。因此亟需开展古龙页岩油含油饱和度参数评价研究。

本文通过系统配套的高精度岩石物理实验，充分发挥物理实验、数值模拟各自的优势，将页岩电阻率物理实验及以高分辨率数字岩心为基础的电性数值模拟两种手段相结合，深入研究页岩储层电性的总体特征、影响因素及规律，建立了基于四维高分辨率数字岩心的含油饱和度模型，并用一维现场核磁饱和度数据标定，确定了古龙页岩油饱和度模型参数，在古龙页岩新井解释、老井复查及井位部署中均取得了较好的应用效果，为松辽盆地古龙页岩油储量提交提供了可靠的饱和度数据。

1 基于数字岩心的古龙页岩油储层含油饱和度计算方法

1.1 饱和度评价方法优选

从岩石物理特征的角度看，页岩储层与富含泥质的粉砂岩储层类似，其电阻率高低不仅与地层水电阻率、含水饱和度、有效孔隙度有关，而且还与泥质电阻率（或黏土电阻率）、泥质含量（或黏土含量）以及泥质（或黏土）的分布形式有关。

从岩矿组成上看，古龙页岩主要为黏土质长英页岩，相较国内外其他页岩油气田，古龙页岩的黏土矿物含量高、长石含量高、碳酸盐含量低。薄片分析表明，古龙地区目的层的黏土矿物与长英质矿物杂乱堆积，属于典型的分散状黏土沉积（图1）。

图1 井A3、S8青一段储层扫描电镜照片Fig.1 SEM photographs of Member Qing-1 in Well A3 and Well S8

近几十年来测井学科针对含分散状黏土的泥质砂岩，提出了如表1 所示的多种饱和度计算模型。

根据表1 中对不同模型的适用性分析，结合古龙地区地质特征，优选Simandoux 提出的适合分散黏土的泥质砂岩饱和度模型，其表达式为

表1 含分散状黏土的泥质砂岩饱和度模型及适用条件Table 1 Saturation models and application conditions for shaly sand containing dispersed clay

式中：Rt——实测地层电阻率，Ω·m；

φcl——黏土体积分数；

Rcl——黏土电阻率，Ω·m；

Sw——地层含水饱和度；

φ——有效孔隙度；

Rw——地层水电阻率，Ω·m；

m——胶结指数；

n——饱和度指数；

a——系数，通常取1。

式（1）中，m、n是关键的饱和度计算参数，需要通过岩电实验测试或数值模拟方法来确定。在古龙地区，根据地质背景资料选取Rw=0.2 Ω·m；对于青一段，可以选取Rcl=2.0 Ω·m。

1.2 数字岩心模拟实验

由于页岩岩电实验驱替困难、饱和度参数难以准确获取，因此开展了基于高分辨率数字岩心物理实验和配套的数值模拟实验，从而确定古龙页岩油储层岩电参数，开展含油饱和度的测井计算。

1.2.1 数值模拟方法概述

岩石物理实验为测井信息向储层信息转换提供了桥梁，但在页岩油气等非常规储层中，这项工作面临极大挑战，主要体现在测试工艺流程、测试精度和分析周期3 个方面。随着计算机技术的发展，对储层的某些岩石物理属性进行数值模拟已成为一种经济有效的研究手段，一方面可以节省工作量和时间，另一方面还可以实现实验室无法测量的某些微观属性特征，如储层润湿性模拟、三相流相渗特征模拟等。

数值计算方法作为一种对流体的流动、传热、质量和能量传输进行预测与评价的手段之一，近年来在多孔介质流动、能量与质量的应用上也越来越广泛，例如岩石中渗透率模拟、电磁场传输特性模拟等。

随着X 射线CT 扫描技术在岩石物理实验中的广泛应用和计算机技术的发展，利用岩心CT 扫描图像构建表征岩石三维微观结构的数字岩心，并提取出只包含孔隙喉道的网络格架，为基于微观结构模拟其宏观物理性质提供了桥梁。数字岩石物理研究（即基于三维数字岩心的岩石物理数值模拟）将在岩石物理理论研究和测井解释与评价，特别是对复杂储层岩石电阻率特性研究具有极其重要的作用。

本文采用针对多孔介质储层比较常用的数值模拟方法有限元法模拟岩石电阻率的基本过程。根据电流连续性方程，对于稳态电流，表达式为

式中：j——电流密度，A/m2；

ρ——体电荷密度，C/m3；

t——时间，s。

根据电流密度与电场强度关系式

式中：σ——电导率，S/m；

E——电场强度，V/m。

若电导率为常数，根据E= -∇φ，得到

式中φ——电势，V。

式（4）即为Laplace 方程。对于一个待求解电阻率的三维数字岩心，将每一像素视为一个单元，通常假设骨架严格不导电，仅仅孔隙流体是导电介质，首先施加一个外电场，然后求解Laplace 方程，可得稳态时三维数字岩心中的电势、电流，进而计算三维数字岩心的有效电导率。

1.2.2 古龙页岩油储层数字孔隙格架构建

1.2.2.1 微米CT扫描

采用CT 扫描构建数字岩心是目前应用最为广泛的建模方法，已在石油工业得到成功应用。微米CT 扫描是利用锥形的X 射线穿透物体，根据射线衰减图像重构得到三维立体图像，由于孔隙的密度明显低于岩石骨架，在图像中孔隙对应的亮度低于矿物骨架，通过设定一个阈值分割、用数字0 和1分别代表骨架和孔隙，就得到一个二值化的图像，代表数字孔隙格架，其中每一个像素点代表着一定尺度的岩石体元，它是当前分辨率能够识别的孔隙或矿物骨架，也可能还包含当前分辨率无法识别的更小孔隙。为了提高后期的数模精度，必须获取高精度、高分辨率、尽可能囊括所有孔隙的CT 图像。

为保证实验结果的配套性，15 块样品直径为2.54 cm，CT 图像分辨率约为7.9 μm。考虑到页岩油主要发育纳米级微孔，其尺度远小于该分辨率，为确保数字岩心模拟结果的可靠性，对样品再次细加工，切割成直径为0.15 cm 的圆柱体样品4 块，编号分别为C2、C5、C6 和F33，采用三英精密NanoVoxel 3 000 型微米CT 仪器对直径为0.15 cm 柱塞样进行扫描，将其像素分辨率提高到1.05 μm。

针对CT 图像，通常采用固定阈值的方法来提取孔隙。以C5 样品为例，采用灰度截止值59 得到三维图像（图2），确定的孔隙度为3.10%。该样品对应的直径为2.54 cm 柱塞样气测孔隙度为9.02%，这说明由于分辨率不够，页岩油储层中仍有2/3 左右的纳米孔没有被微米CT 识别。为了确保后续数值模拟的准确计算，必须采用其他方法确定出纳米孔隙含量及其分布特征，这也是影响古龙地区页岩油导电特征的关键。

图2 井S8青一段C5样品微米CT及阈值分割后的三维图像Fig.2 Micron CT and 3D images after threshold segmentation of Sample C5 of Member Qing-1 in Well S8

1.2.2.2 利用聚焦离子束扫描分析页岩油储层纳米孔隙

聚焦离子束扫描（Focused Ion Beam，简称FIB）是将液态金属离子源产生的离子束经过加速、聚焦后照射到样品表面产生二次电子信号获得高分辨率图像的方法。

对直径为2.54 cm 样品圆柱端面进行抛光和镀碳层处理后，选择一块1 cm×1 cm 的方形区域，再精细划分成41×41 个子区域，对每一子区域采用FIB 扫描获得10 nm 分辨率的图像，将这些图像拼接成一张超高分辨率、超大面积的图像，这一过程称为MAPS 成像，据此可以对样品中的各种孔隙，包括残余粒间孔、溶孔、有机质微孔的形状、分布进行精细分析（图3）。可以看出，青一段页岩油储层发育大量黏土微孔隙，其尺寸基本小于1 μm，CT 图像不能识别。沿矿物层理面发育的微裂缝开度基本在微米级，延展长度大，CT 图像能够识别。

图3 井S8青一段C5样品MAPS图像显示粒间孔、黏土晶间孔及裂缝Fig.3 Intergranular pores,clay intercrystalline pores and fissures in Sample C5 of Member Qing-1 in Well S8

为了进一步明确古龙地区青一段不同矿物的纳米孔发育程度，对样品端面同一个区域开展MAPS成像与扫描电镜矿物定量评价（QemScan）矿物测量，通过对比分析，就可以确定页岩油储层中各类无机、有机矿物的孔隙类型与含量（面孔率），如图4 所示。

表2 给出了4 块样品的QemScan 扫描矿物含量结果。应该强调的是，虽然都是矿物含量分析，但是QemScans 仅针对圆柱形岩样端面很小的区域进行扫描，考虑到页岩油储层微观的非均质性，表2的结果与X 衍射黏土矿物定量分析（XRD）的全岩矿物并没有可比性，但揭示的本区页岩油储层主要岩矿类型是一致的。

表2 4块样品矿物含量电镜扫描分析结果Table 2 SEM analysis of minerals contents of 4 samples

根据图4 的成果，通过对比将纳米级分辨率MAPS 灰度图像划分为孔隙、有机质、伊利石、石英/钠长石、方解石和黄铁矿等高密度矿物。据此确定不同矿物的MAPS 图像灰度分布区间及孔隙发育比例，结果见表3。在表3 中，阈值1、阈值2 分别表示在MAPS 灰度图像中不同组分的灰度分布下限、上限值；体积比例是指在直径2.54 cm 柱塞样品微米CT 图像中该组分所占的比例；孔隙率是指基于MAPS 纳米分辨率图像确定的该组分孔隙发育程度，或者理解为面孔率。

表3 井S8青一段C5样品四维数字图像中矿物—孔隙映射关系Table 3 Mineral-pore mapping relationship in digital image of Sample C5 of Member Qing-1 in Well S8

图4 井S8青一段C5样品MAPS与QemScan图像对比Fig.4 Comparison between MAPS and QemScan images of Sample C5 of Member Qing-1 in Well S8

将微米CT 未识别的孔隙归位到有机质和伊利石、方解石等无机矿物中，最终确定数字岩心总孔隙度。数字岩心总孔隙度就是将不同组分所占比例乘以其对应的孔隙率，然后求和得到数字岩心的总孔隙度，计算表达式为

式中：φt——总孔隙度；

φp——粒间孔隙所占比例，%；

φorg——有机质所占比例，%；

φil——伊利石所占比例，%；

φca——方解石所占比例，%。

采用上述处理方法得到的四维数字岩心，其孔隙度约为8.44%，与气测孔隙度9.02%基本接近，说明已经将微米CT 图像丢失的纳米级孔基本上找到并准确归位，据此可以开展数值模拟。

之所以称为四维数字岩心，是指在直径为2.54 cm 柱塞样品的四维数字岩心图像中（图5），每一个像素点（对应三维坐标X、Y、Z）并不是传统上的非0 即1、非骨架即孔隙，而是每一个像素点都对应某种有一定孔隙率的组分，将孔隙率作为第4 个坐标，样品C5 的处理结果见图5。

图5 井S8青一段C5样品的高分辨率四维数字岩心图像Fig.5 High resolution 4D digital core image of Sample C5 of Member Qing-1 Well S8

需要说明的是，对于每一块岩心都需要重复上述测试和分析过程，得到类似表3 的映射关系。在不同地区、不同层位，甚至于同一口井同一层位不同深度的样品，上述矿物与孔隙映射关系可能都不同。

2 饱和度模型参数的确定

2.1 各组分等效电导率确定

在获得了囊括页岩油储层中所有类型孔隙的数字岩心格架基础上，就可以利用前面介绍的有限元法计算其不同方向的电阻率，分析其导电特征。

根据以上分析，古龙页岩油多组分三维数字岩心包括5 种相：粒间孔隙、石英/钠长石、方解石、伊利石、有机质。其中，粒间孔隙部分的电导率由孔隙流体确定；石英/钠长石、方解石两种相具有一定的面孔率，其等效电阻率根据地层水电导率和Archie 公式计算，由于前者的次生溶孔发育程度低，其面孔率很小，因此只需要讨论方解石；伊利石为黏土矿物，具有附加导电性，其等效电阻率利用Waxman-Smith 模型计算。

2.1.1 方解石

该组分孔隙度为7%，则其电阻率根据地层水电阻率和Archie 公式共同确定，表达式为

式中：Rkfe+ca——方解石等效电阻率，Ω·m；

φca——方解石孔隙度，%。

电导率与电阻率互为倒数，然后可以获得该组分的等效电导率Ckfe+ca。

2.1.2 伊利石

多组分孔隙格架中的伊利石部分与泥质砂岩类似，区别在于该相骨架全部为黏土矿物，因此采用泥质砂岩的Waxman-Smith 模型（简称W-S 模型）计算伊利石部分电导率。根据W-S 模型，岩石饱和地层水时电导率表达式为

式中：C0——含水泥质砂岩的电导率，S/m；

Ccl——黏土交换阳离子产生的电导率，S/m；

Cw——自由电解液（地层水）的电导率，S/m；

x，y——系数。

自由电解液与可交换阳离子适用于同一个几何常数，即系数x=y，因此得到考虑地层因素的含水泥质砂岩电导率公式为

式中F——在Cw足够高时，泥质砂岩的地层因素。

其中

式中：Qv——阳离子交换容量，mol/L；

CCEC——阳离子交换能力，取值0.3 mmol/g；

φsh——黏土体积分数，取值0.85；

φ——孔隙度，取值0.15；

ρ——平均颗粒密度，g/cm3。

2.1.3 有机质

有机质发育微孔隙，MAPS 分析结果显示古龙页岩油有机质面孔率约为20%，其等效电阻率可根据Archie 公式确定，有机质等效电阻率的表达式为

式中Rorg——有机质等效电阻率，Ω·m。

电导率与电阻率互为倒数，然后可以获得该组分的等效电导率Corg。考虑到有机质孔隙表面以亲油为主，存在油相吸附层，所以设定胶接指数m的取值为2.3。

扫描电镜图像显示方解石的溶蚀孔连通较差，而伊利石中晶间孔发育，连通性好。在确定各矿物组分等效电阻率过程中，还涉及式（6）、式（8）、中胶结指数m的选取。这是根据样品在实验室加压饱和水时测量其电阻率Ro，通过调整m计算等效电导率，利用有限元法多次模拟多矿物组分三维数字岩心的电导率，并与该岩样的实验结果对比，从而确定m的取值。本文模拟中，通过多次反复测试，最终确定方解石组分的m取值为2.0；伊利石组分m取值为1.6。

2.2 地层因素模拟

根据前面确定的各组分等效电导率，利用有限元法模拟确定在z方向上电压为Uz、电流为Iz时对应的电导率公式为

式中：Uz——电压，V；

Iz——电流，A。

已知多组分三维数字岩心的孔隙度φ，设a=1，可得到胶结指数m。利用该方法对4 块样品模拟结果与实测结果对比（图6）。

图6 古龙地区4块样品数值模拟与实测的F-Φ图版Fig.6 Numerical simulated and measured F vs.φ of 4 samples of Gulong area

从图6 可以看出，利用增水法实验测量确定的古龙页岩油m实验=1.116，利用数字岩心模拟得到的m数模=1.50，两者相差明显。对于本区页岩油储层而言，其黏土体积分数一般在35%以上，广泛发育纳米级微孔，考虑到其中纳米级孔的饱和水过程很难彻底完成，因此实验测量与数值模拟结果的差异是必然存在的。

图6 中数值模拟的结果显示4 块样品的F与φ关系较为发散，这主要是由于取样井段的岩性岩相差异所致。从表2 可以看出，样品C6 主要矿物为石英、方解石和钠长石，特别是方解石含量远高于另外3 块样品，因此样品C6 应该代表本区钙质硅质混合页岩相，单独确定其胶结指数为1.73。将样品C2、C5 和F33 的数据单独分析，得到如图7 所示的图版。

图7 C2、C5和F33三块样品数值模拟F -Φ图版Fig.7 Numerical simulated F vs.φ of C2,C5 and F33 samples

综合上述分析，根据4 块样品的数字岩心模拟结果，确定本区岩电参数为：黏土质硅质页岩相储层，a=1，m=1.375；钙质硅质混合页岩相储层，a=1，m=1.73。

2.3 电阻率增大率模拟

对于页岩这类非常规储层，试图通过驱替测量来解决电阻率增大率的测试是不可能的，一方面岩心无法洗油，另一方面，由于其极低的渗透率，现有实验测试装置都无法实现在页岩储层中用油驱替水的过程，因此只有通过数字岩心，利用图像运算来模拟油驱水的过程，这可以通过数字图像的运算来实现。

图像的腐蚀和膨胀也是依据数学形态学方法发展起来的一种图像处理算法，其基本思想就是用具有一定形态的结构元素去度量和提取图像中对应的形状，实现图像特征分析与目标识别。

2.3.1 结构元素

结构元素是图像腐蚀和膨胀运算的基本组成部分，通常比待处理图像小得多，在二维平面上就是一个数值为0 或1 的矩阵。其原点指向图像中需要处理的像素位置，大小决定需要处理的范围，数值为1 的点决定结构元素的邻域像素在运算时是否需要参与计算。

2.3.2 腐蚀运算

腐蚀运算就是用一个结构元素，如3×3 的矩阵，扫描数字图像的每一个像素，结构元素与其覆盖的图像进行“与”操作，如果都为1 输出图像的像素为1，否则为0。显然，这种算法使得处理目标减小一圈。

2.3.3 膨胀运算

膨胀运算也称为扩张运算，过程同上，如果都为0 输出图像的像素为0，否则为1。算法处理结果是使得处理目标扩大一圈。

2.3.4 开运算

开运算将腐蚀和膨胀结合起来，对图像作先腐蚀后膨胀处理就是开运算，它可以平滑图像的轮廓，消弱图像狭窄的部分。开运算的主要作用与腐蚀运算相似，但具有基本保持目标原有大小不变的优点。

2.3.5 闭运算

对图像作先膨胀后腐蚀的处理就是闭运算。为了直观显示图像运算结果，以二维孔隙介质图像为例介绍上述过程（图8）。图8（a）为原始的二维数字图像，其中黑色代表颗粒，用0 表示，白色代表孔隙，用1 表示，图像尺寸为200×200 个像素。选取半径r=5 个像素的圆作为结构元素，对图8（a）中的孔隙空间（白色区域）分别进行腐蚀、膨胀和开运算，结果如图8（b）—（d）所示，图8中灰色部分代表孔隙空间经过相应运算后的结果，腐蚀运算收缩目标图像，膨胀运算扩大目标图像，开运算可以理解为结构元素B在X内滚动所能达到的最远处的B的边界所构成的空间。因此，如图8（d）中灰色区域所示，开运算结果显示所有大于半径r的孔隙空间。

图8 二维数字图像腐蚀、膨胀和开运算结果Fig.8 Diagram of corrosion,expansion and open operation results for 2D digital images

设岩石孔隙空间中最大孔隙半径为rmax，对岩石孔隙空间进行开运算，其结果表征的孔隙空间按照孔隙半径的大小依次增加。若假设孔隙空间的开运算结果表征油驱水过程中的油，其余孔隙空间表征驱替过程中的水，则该过程与水湿岩石的排驱过程相似，即非润湿相油首先占据大孔隙，随着驱替压力增大，油按照孔隙半径由大到小的顺序依次侵入小孔隙。

因此，利用岩石孔隙空间的开运算可以用来模拟水湿岩石的排驱过程，进而确定在不同含水饱和度下孔隙空间中油和水的分布特征，并计算得到不同时期相应的电阻率。图9 为这一模拟过程的实例。

图9 利用开运算模拟砂岩储层数字图像不同含水饱和度时孔隙空间油水分布特征Fig.9 Simulation of oil and water distribution in pore space of digital image in sandstone reservoir with different water saturation by open operation

仍以C5 样品为例，针对图5 中得到的高分辨率四维数字岩心图像，分别采用结构元素D=10、7和4 进行开运算处理，模拟3 种不同含水饱和度下的数字图像（图10），图中灰色为岩石骨架，孔隙中流体分别为油（红色）、水（绿色）和水膜（黄色）。

图10 井S8青一段C5样品采用不同结构元素模拟油驱水过程的图像Fig.10 Images of simulating oil flooding water process with different structural elements for Sample C5 of Member Qing-1 in Well S8

同样地，利用上述方法完成了4 块数字岩心的电阻率增大率I与含水饱和度Sw关系模拟，共完成了29 个含水饱和度节点的模拟计算（图11）。综合分析确定本区岩电参数为b=1，n=1.42。

图11 4块样品数值模拟I-Sw关系图版Fig.11 Numerical simulated I vs.Sw of 4 samples

利用常规增水法岩电测试和数值模拟方法确定了页岩油储层孔隙指数m和饱和度指数n，见表4。

表4 15块柱塞样增水法岩电测试及4块岩心数值模拟确定的岩电参数Table 4 Statistics of rock electrical parameters determined by rock-electricity test by water increment method of 15 core plugs and numerical simulation of 4 cores

3 应用实例

对有配套岩心测试数据的单井进行了饱和度计算，图12 是井A7 的饱和度解释成果，第7 道的棒状线是井场一维核磁提供的岩心含油饱和度数据。经60 块岩心含油饱和度标定，电阻率测井计算含水饱和度平均绝对误差3.9%。结果表明，利用上述模型可以准确计算古龙地区页岩油储层含水饱和度。

图12 古龙页岩油含油饱和度处理解释成果(井A7)Fig.12 Processing and interpretation of oil saturation in Gulong shale oil(Well A7)

应用上述含油饱和度评价方法开展页岩油单井评价，共解释新老井200 多口，其中完成测试直井24 口，22 口井试油见油流，其中工业油流井16口。在单井解释基础上，为富集层段落实提供了可靠的饱和度数据，为井位部署、水平井设计及实验区开发动用提供了技术支持。同时上述含油饱和度精度达到了储量规范需求，已经应用到了2021 年古龙页岩油提交预测储量评价中，并为未来古龙页岩油探明储量的顺利提交提供可靠技术保障。