中国南海蚀变火成岩导电机理研究与低阻成因分析

郭宇航, 潘保芝*, 冯进, 管耀, 王清辉, 张丽华, 阿茹罕, 王欣茹, 张乃毓, 张鹏济, 李岩, 贾瑞龙, 高云翔

1 吉林大学地球探测科学与技术学院, 长春 130026 2 中海石油(中国)有限公司深圳分公司, 深圳 518000

0 引言

近年来,随着国家经济高速发展和国际局势的变化,油气资源的需求越来越大.为了满足发展需求,近几年各大油田将勘探重点转移到火成岩气藏上,松辽盆地、渤海湾盆地、四川盆地和南海等地区相继发现大型火成岩气藏(刘军等,2022;杨晓利等,2022).岩石导电机理一直是岩石物理性质研究中的重难点,以往岩石导电机理研究对象主要为砂泥岩和页岩,并且仅考虑岩石孔隙结构、矿物组成和地层水电阻率等影响因素,这导致传统岩石导电机理在火成岩领域存在一定的局限性.同时普遍存在于火成岩的蚀变现象使其的孔隙结构等储层特性发生了复杂的变化(刘思慧,2016),令蚀变火成岩导电机理研究比以往的砂泥岩和页岩都更为复杂(Guo et al.,2022).尽管研究火成岩导电机理困难重重,但是广泛存在的火成岩储层是中国石油勘探的重要领域之一,蚀变火成岩导电机理和饱和度模型的研究是天然气勘探过程中无法忽视的问题,因此研究蚀变火成岩导电机理可为油气勘探和储层评价提供可靠的支撑.

然而,影响火成岩导电的因素较多,如:岩性、矿物组成、导电矿物分布、蚀变程度、孔隙度、润湿性和流体分布等,通过常规岩石物理实验想要单独研究某一项影响因素是比较困难的,而且对于真实岩心的需求也是大量的,无论成本和时间上都是一种巨大的消耗.随着CT扫描技术的不断发展,岩心扫描图像经过处理后可以得到数字岩心,基于数字岩心的岩石物理数值模拟被称为数字岩石物理实验(陶果等,2005;孙建孟等,2012;赵建国等,2014,2021;孔强夫等,2015;Yang et al., 2021;吴建彪等,2023),相比于传统岩石物理实验,数字岩石物理实验计算速度快且不会破坏岩心.相比于传统岩石物理数值模拟,数字岩心技术能建立与实际岩心十分接近的岩石物理模型,重现岩石内部复杂的孔隙网络结构(刘学锋,2010;聂昕,2014;王建福,2017;吴丰等,2019;范雨霏,2021).目前针对蚀变火成岩的数字岩心电性模拟的研究十分稀少.基于数字岩心模型可以模拟孔隙发育、溶蚀、蚀变等地质活动,从而进行更深入的岩石物理性质模拟研究.目前针对数字岩心电性数值模拟主要采用的方法有随机游走法、格子玻尔兹曼方法和有限元法(孔强夫等,2015).考虑到蚀变矿物中含有黏土矿物,国内目前针对黏土矿物的数字岩心电性模拟主要利用有限元法(刘学锋,2010;聂昕,2014;范雨霏,2021)和格子气自动机方法(岳文正等,2004),上述研究主要是针对泥质砂岩、页岩和碳酸盐,目前还没有针对蚀变火成岩的数字岩心电性模拟研究.

鉴于目前蚀变火成岩导电机理研究中存在的诸多问题,本文采用数字岩心模拟技术研究南海某区块蚀变火成岩储层岩石导电机理.研究成果可以为厘清蚀变程度、蚀变矿物类型、润湿性、流体分布和孔隙度对火成岩导电性的影响机制提供可靠的依据.考虑地质客观规律的矿物定向膨胀、孔隙溶蚀和混合润湿数字岩心构造法,为蚀变岩石数值模拟提供了新的思路.综上所述,本文从蚀变火成岩岩石导电机理的科学问题研究出发,旨在突破蚀变火成岩储层测井评价的技术瓶颈,在油气勘探开发技术层面具有广泛的应用前景.

1 数字岩心构建

选取中国南海惠州某区块具有代表性的火成岩储层岩石样品,图1为该区块蚀变火成岩岩石薄片,图中能够看到斑状变晶结构、变余斑状结构;基质为微-细粒状鳞片变晶结构,片状(流状)构造.变斑晶为角闪石,含量 20%,有单晶和聚斑两种;聚斑者多呈透镜状,具多颗成串似骨排状分布特征(图1a);分析是早期糜棱岩化变形的辉石碎斑后期被角闪石交代形成,部分碎斑中含有辉石的交代残余单晶变斑晶(图1b).部分角闪石斑晶呈半自形柱状,边缘参次不齐,是后期变成(图1a、c).基质由:细-微粒的角闪石、黑云母(50%)和少量的长英质微粒组成,角闪石和黑云母强烈定向,显示出片状构造.长英质微粒轮廓不清,显示有明显的泥化蚀变特征,这种明显的蚀变特征在研究区内较为普遍,需要在数字岩心建模时考虑蚀变矿物的类型和含量.为了后续模拟不同蚀变程度的岩石,本文选用了蚀变程度较低的安山岩进行CT扫描,样品的具体信息如表1所示,不同组分的含量是根据全岩分析和孔隙度实验室测量得到的.数字岩心图像处理,阈值分割,电性模拟,流体分布构建,模拟蚀变和孔隙溶蚀过程相关程序模块都为自主研发,本文研究流程图如图2所示.

表1 岩石样品信息Table 1 The basic information of core

图1 南海某地区蚀变火成岩岩石薄片Fig.1 Thin section of altered igneous rocks in an area of the South China Sea

图2 蚀变火成岩导电机理研究流程图Fig.2 Workflow of the study of the conductive mechanism of altered igneous rocks

1.1 图像预处理

岩石样品经过切割后被放入CT扫描仪器中进行成像扫描测量,根据之前所观察到岩心外观特征,为了最佳量化分析岩石内部的孔隙结构,样品使用微米CT扫描仪时的扫描分辨率大致范围初步定为0.76 μm/像素.制备样品时,样品的尺度不超过4 mm×4 mm×4 mm,可使整个样品可以完全呈现在显微镜视野范围内,以便观察扫描图像与三维结构.微米CT平板接收器的可接受像素范围约为2000×2000×2000,在设定了分辨率范围之后,即可大致推知整个视野范围内可以容纳何种尺寸的样品了,本次扫描样品切割尺寸为1.52 mm×1.52 mm×0.99 mm,像素范围为2000×2000×1299.由于在CT扫描图像常常存在各种噪声、尾影、光圈和坏点等问题,所以需要对所有切片图像进行统一的预处理.一般地,通过滤波处理,除去图像中的噪声,保留岩心本身的内容,对于坏点、尾影和光圈,可以分别采用插值、剪裁的方式.由于滤波、插值和剪裁算法众多,本文在此不在赘述和比较,本文选用中值滤波、三次样条插值和切边剪裁的方式,图像预处理前后对比如图3所示(刘学锋,2010;范雨霏,2021).

图3 图像预处理前(a)后(b)对比图Fig.3 Comparison of images before (a) and after (b) pre-processing

1.2 阈值分割与多组分三维数字岩心构建

经过X射线CT扫描的岩石三维灰度图像,像素的灰度值反映了该像素点对应岩石组分对X射线吸收的强弱,和密度呈现很好的对应关系.通常地,对于矿物组分单一的岩石,通过图像二值化阈值分割来实现骨架和孔隙的分割,灰度值小,亮度暗的像素点,对应密度小,表示岩石孔隙空间,相反的,灰度值大,亮度高的像素点,对应密度大,表示岩石骨架.采用阈值分割法效率高,阈值依据相关实验数据来确定,是一种可靠高效的图像分割方法.

火成岩矿物成分复杂,不仅有密度较高的暗色矿物,还有蚀变产生的低密度黏土矿物,由于黏土矿物中含有微孔隙,而受限于扫描分辨率,该部分纳米级微孔无法被识别出来,所以在CT成像上黏土矿物呈现较低的灰度值.针对多矿物阈值分割,目前在实验方面可以采用QEMSCAN(Quantitative Evaluation of Minerals by Scanning Electron Microscopy)和MAPS(Modular Automated Processing System)技术作为配套,以此为基准构建多组分矿物的三维数字岩心,这些方法具有较高的分割精度,但是成本较高.本研究拟采用全岩矿物分析与CT扫描图像联合确定多组分分割阈值的方法.对于全岩矿物分析的结果,需要对物理属性相近的矿物进行合并处理,例如石英和长石合并为长英类,伊利石、高岭石、蒙脱石和绿泥石合并为黏土类矿物,辉石和角闪石合并为暗色矿物等.值得注意的是,这并不是一个可以泛化的合并方案,在火成岩中例如花岗岩,存在大量的钾长石,而钾长石物理属性比较特殊,需要单独划分为一种组分.本文研究以安山岩为例,具体矿物合并和分割示意图如图4所示,表2为分割结果与全岩分析对比结果.图5是三维多组分蚀变火成岩数字岩心分割结果.

表2 数字岩心多矿物组分分割结果Table 2 Multi-component segmentation results of digital core

图4 多组分阈值分割后各组分三维空间分布图Fig.4 3D spatial distribution of each component after multi-component threshold segmentation

图5 三维多组分蚀变火成岩数字岩心Fig.5 3D multi-component altered igneous digital rock

2 有限元电性模拟与流体分布构建方法

岩石的电性一般包括电阻率特性、复电阻率特性和介电特性等,其中岩石的复电阻率和介电在不同频率下具有一定的频散特性.从狭义上讲,本文讨论的岩石的电性指在直流电下岩石的电阻率特性.

2.1 有限元电性模拟原理

(1)

取第k个像素体为立方体单元k,该单元内电势函数取为单元顶点处电势的线性插值,即:

(2)

(3)

再令:

(4)

则整个岩心的功率为

(5)

单元总顶点数大于单元总数,设岩心的尺寸为nx×ny×nz,则x=nx、y=ny、z=nz等处的位移量是需要外加条件来确定.

对整个岩心施加给定的稳恒电场e0,则岩心的电势可转化为周期边界,总功率改写为

(6)

根据能量极小原理,在稳恒电场e0作用下,岩心中各处电势最终将松弛到使整个岩心电势能即总功率最小的值,满足方程:

(7)

(8)

2.2 不同蚀变黏土矿物类型影响

在本文多组分数字岩心模型中,导电相为孔隙中的水相和蚀变黏土矿物,其他矿物由于导电性与上述两相相差过大,模拟过程中设置电导率为0.0001 S·m-1.地层水电导率和黏土矿物附加导电性参考刘学锋(刘学锋,2010),聂昕(聂昕,2014)和范雨霏(范雨霏,2021)研究中的参数布置.采用这套参数模拟了绿泥石、高岭石、伊利石和蒙脱石4种蚀变黏土矿物,如图6所示,在饱和地层水的情况下,蚀变火成岩整体电导率随着地层水电导率的增大而增大,其中蚀变产生蒙脱石、伊利石、绿泥石和高岭石的火成岩电导率依次减小,黏土矿物的附加导电性明显.

图6 蚀变火成岩数字岩心不同黏土矿物Cw-C0 电性模拟结果Fig.6 Simulation results of Cw-C0 electrical properties of different clay minerals in digital core of altered igneous

2.3 流体分布构建方法

上述模拟是对单向流体饱和的岩石展开的,对于双相流体饱和岩石(本文讨论油水两相),导电相流体与非导电相流体在岩石孔隙空间的分布会影响岩石整体的导电性.所以研究不同含水饱和度下的流体分布是多相流体饱和数字岩心电性模拟的重中之重.本文采用的方法是基于数学形态学来构造不同含水饱和度下的数字岩心.

2.3.1 数学形态法构建

数学形态学中的两种基本运算膨胀和腐蚀,如图7和图8所示,在得到目标图像和结构元素后,需要确立结构元素的原点,然后依据顺序用结构元素扫描目标图像,最后会对原点所占位置的对应图像进行判断,在膨胀运算时结构元素中为1的元素位置上对应的图像值至少有一个为1,则原点所占位置图像元素被赋值为1,否则赋值为0;在腐蚀运算时结构元素中为1的元素位置上对应的图像值全部为1,则原点所占位置图像元素被赋值为1,否则赋值为0.所以从效果上看膨胀运算扩大图像,腐蚀运算收缩图像.

图8 腐蚀运算示意图Fig.8 Diagram of erosion operation

如图9所示,将这两种运算方式组合得到了开运算和闭运算.开运算能平滑图像的边界,削弱图像狭窄的部分.以A为待处理的图像,用结构元素B作开运算,记为A○B,用集合表示为:

图9 开闭运算示意图Fig.9 Diagram of opening and closing operation

A○B=(A⊖B)⊕B,

(9)

其中⊕为膨胀算子,⊖为腐蚀算子,孔隙空间开运算结果代表油驱水过程中的油相,如果孔隙的最大半径是Rmax,结构元素半径为Rmax的时候,孔隙空间经过开运算后填充物质完全由水变为油.开运算可以模拟排驱过程,进而确定在不同含水饱和度状态下孔隙空间中油和水的分布状态.为准确模拟岩石孔隙空间的流体分布,必须保证先被油侵入的孔隙空间在后续排驱过程中均被油占据.因此,结构元素半径小的开运算结果必须包含半径大的结构元素开运算结果,即:

ϑS(R)⊂ϑS(R′),

(10)

其中,ϑ代表开运算的结果,S(R)代表结构元素.R和R′为不同结构元素的半径,且R>R′.由于三维数字岩心为数字化的二值图像,因此三维结构元素—球体在数值模拟过程中也只能通过离散方式表示.球形结构元素S′(R)的表达式为

S′(1)=S(1),

(11)

S′(R)=S(R)∪S(R-1).

(12)

开运算计算过程与毛管饱和度模型和蚀变校正-毛管饱和度模型流体饱和过程一致,油优先占据大孔隙中间位置.获得不同含水饱和度下孔隙空间的流体分布后,采用有限元方法计算其电阻率,进而结合饱和度模型计算岩石的电阻率相关参数,最后与岩心电阻率实验进行比对,验证模型效果.

2.3.2 考虑混合润湿性的流体分布构建

基于多组分蚀变火成岩数字岩心模型,矿物和孔隙的空间分布已知,岩石整体的润湿性也会对流体分布形式产生影响,从而影响岩石电性.本文分别模拟了完全油湿、完全水湿和混合润湿三种情况,图10左为岩石流体分布切片图,图中绿色为骨架,红色为油相,蓝色为水相,当前含水饱和度为52%,对应图10右油湿电阻率模拟结果中黑色的点.当岩石整体表现为油湿时,在高含水饱和度段,油湿岩石孔隙内部表面被油相附着,水相则占据大孔中心被油相包裹,导致导电连通性不好,电阻率随着含水饱和度的降低而迅速升高,当含水饱和度下降到一定程度时,由于孔隙整体连通性较差,孔隙度不高,流体影响逐渐减弱,黏土矿物附加导电逐渐占据主导地位,这导致油湿和水湿岩石在低含水饱和度段趋于一致.无论油湿还是水湿都呈现了明显的非阿尔齐现象.

图10 不同润湿性蚀变火成岩I-Sw交会图Fig.10 The relationship between I and Sw of different wettability in alteration igneous cores

岩石润湿性是各种矿物表面润湿性的综合反应,常规岩石物理实验无法获得孔隙周围矿物组分信息,因此通常采用岩石整体润湿性这种表述方式.然而,实际火成岩内部矿物分布是非均质的,各种矿物的润湿性是不同的,例如白云石疏水,石英亲水(王贤晨,2020).这就造成了非均质润湿的情况,基于多组分数字岩心能够获得岩石孔隙空间附近矿物的种类信息,依据矿物润湿性认识可以定义混合润湿流体分布状态,由于本地区蚀变黏土矿物以绿泥石为主,所以蚀变黏土矿物按包膜绿泥石润湿性设定为疏水(闫顶点,2021),长英类矿物为亲水,碳酸盐矿物以方解石和白云石为主设定为疏水,暗色矿物主要以辉石为主设定为亲水.图10中黑色曲线为混合润湿的电阻率指数随含水饱和度变化的模拟结果,整体接近油润湿曲线,这是由于孔隙周围主要以油润湿矿物为主,基于数字岩心的混合润湿模拟结果更贴近岩石的真实情况.

3 基于地质客观规律的蚀变火成岩数字岩心构建

3.1 基于定向膨胀法构建变蚀变程度的火成岩三维数字岩心

在本文中为了更好的研究蚀变程度对火成岩岩石导电性的影响,需要在保证其他主要影响因素如孔隙结构不变的情况下进行控制变量模拟研究.如图11所示,由于研究区主要蚀变产物为绿泥石、高岭石、蒙脱石、绢云母和方解石,主要是由长石、辉石和角闪石蚀变而成,其中黏土矿物具有附加导电性,其他矿物不导电(Wood and Hazra, 2017).所以,基于三维数字岩心矿物分割结果,采用定向膨胀运算,将蚀变矿物所在像素点连接的边界向长石、辉石和角闪石所在像素点膨胀,并替换像素点的组分属性,根据地质规律,如图12所示,辉石和角闪石会蚀变产生方解石和绿泥石,而长石不会蚀变产生方解石,由于蚀变作用是由于热液流入孔隙空间,所以可以设定发生蚀变的矿物围绕在孔隙周围.设长石蚀变出现高岭石和蒙脱石的概率各为50%,基性斜长石和橄榄石蚀变产生两种矿物的概率分别是50%.辉石和角闪石蚀变产物相同,产生方解石和绿泥石的概率也各为50%.定义蚀变程度为蚀变矿物所占像素点与总矿物所占像素点的比值.在采用定向膨胀法时,从孔隙周围矿物出发,根据随机概率模型,对原位矿物类型进行判断生成蚀变矿物,定向膨胀后可以统计出蚀变程度、黏土矿物和方解石等蚀变产物含量.

图11 蚀变矿物与一次蚀变产物示意图 (Wood and Hazra, 2017)Fig.11 Schematic diagram of alteration minerals and primary alteration products (Wood and Hazra, 2017)

图12 定向膨胀示意图0:孔隙;1:不发生蚀变的矿物;2:蚀变矿物;3:长石类矿物;4:辉石;5:角闪石.Fig.12 Schematic diagram of directional swelling0:Pores; 1:Non-altered minerals; 2:Altered minerals; 3:Feldspar-like minerals; 4:Pyroxene; 5:Hornblende.

图13是通过定向膨胀法构造的不同蚀变程度的三维数字岩心电性模拟结果,值得注意的是,图中蚀变黏土矿物是按伊利石电导率模拟的.如图所示,蚀变程度越高,电阻率指数随含水饱和度变化越小,并且在高含水饱和度处(约为95%)弯曲现象越不明显.当蚀变黏土含量高于一定值时,I-Sw曲线近乎直线,这说明在构建蚀变火成岩导电模型时需要考虑高低蚀变程度带来的差异.此外,模拟结果符合低阻储层的特征,这与许多学者的研究规律相符(刘学锋,2010;聂昕,2014;范雨霏,2021).

图13 变蚀变程度数字岩心I-Sw交会图Fig.13 The relationship between I and Sw of variable alteration level of igneous cores

3.2 基于定向溶蚀法构建变孔隙度的火成岩三维数字岩心

按照常规研究方法一般是在一个地区通过压汞、核磁等孔隙结构相关表征实验对岩石孔隙结构进行分类评价,然后结合岩电实验进行孔隙结构对岩石导电性的影响分析.这种方法的弊端是无法完全控制变量,不同岩石孔隙结构改变后无法保证蚀变程度、蚀变矿物类型等参数不变.所以本项目拟对多矿物分割后的数字岩心数据体进行处理,对孔隙所占的像素点属性进行替换,替换成非蚀变矿物属性.扩大孔隙度时,采用定向膨胀法只将孔隙附近非蚀变矿物属性变为孔隙.讨论孔隙度影响时要保证其他蚀变程度和蚀变矿物类型等因素不变,原孔隙所占像素点属性替换为非蚀变矿物,原有蚀变程度和蚀变矿物类型不变,如图14二维切片示意图所示,保持孔隙基本形态不变.图15为相同蚀变程度不同孔隙度的三维数字岩心电性模拟结果.随着地层水电阻率电导率的降低,F-φ的非阿尔奇现象越明显.

图14 孔隙填充示意图Fig.14 Schematic diagram of pore filling

图15 构造变孔隙度数字岩心不同Cw下的F-φ交会图Fig.15 The relationship between F and φ of variable porosity digital cores with different Cw

4 蚀变火成岩导电机理和低阻成因分析

一般地,高地层水矿化度、高泥质含量、岩石强亲水、骨架导电和良好的物性等因素是造成储层低阻的主要原因.在目标研究区的储层主要位于古潜山风化带和内幕带顶部,如图16所示,通过最优化处理后的结果显示B井蚀变程度较高,物性也要好于A井.黑色线框内B井平均电阻率为0.95 Ωm,而A井平均电阻率为27.51 Ωm.根据图17中A井三维孔隙空间和蚀变黏土矿物空间分布可以看出,本区块安山岩孔隙具有一定的连通性,但孤立孔隙大量存在,蚀变黏土矿物包裹在其他矿物(暗色矿物和长英类矿物)颗粒表面并且位于孔隙旁,其附加导电性对于孤立孔隙起到了一定的连通作用,利于提高岩石导电性.

图17 孔隙网络模型和黏土矿物空间分布Fig.17 Pore network model and spatial distribution of clay minerals

图15变孔隙度数字岩心模拟结果显示,随着孔隙度增大,饱和地层水岩石电阻率降低,随着地层水电导率的升高,孔隙度带来的影响增大.图13模拟变蚀变程度模拟结果显示,当孔隙度一定时,随着蚀变黏土含量的升高,岩石电阻率受流体饱和度的影响变弱,这也是造成油层低阻的原因.值得注意的是,虽然孔隙具有一定的连通性,但是孔隙分布主要以次微米级微孔为主(受扫描分辨率限制).结合图10混合润湿的模拟结果可以看出,蚀变安山岩孔隙附近矿物以亲油矿物为主,整体呈现油湿的现象,这在一定程度上降低岩石导电性,油分布在微孔和大孔孔壁处,不利于油的产出.

此外,数字岩心模拟结果也显示蚀变安山岩具有较为明显的非阿尔奇现象,阿尔奇公式在蚀变火成岩储层并不适用,需要结合导电机理研究结果,考虑润湿性、蚀变程度以及孔隙结构等诸多因素,建立相应的导电模型和饱和度方程.

5 结论

本文采用有限元和数学形态法对南海某区块蚀变火成岩储层进行三维数字岩心的电阻率模拟,其中三维多组分数字岩心通过CT扫描成像、全岩分析和柱塞样物性测量等多尺度岩石物理实验共同建立,多组分分割后的结果与实际岩心实验结果相符.基于定向膨胀法和溶蚀法构造了变蚀变程度和变孔隙度的蚀变火成岩数字岩心,该方法能够遵循客观地质规律,为控制变量研究导电机理提供了可靠的支撑.模拟结果显示,蚀变黏土矿物的附加导电性按照蒙脱石,伊利石,绿泥石和高岭石的顺序递减;其他因素不变时,蚀变程度越高,黏土矿物含量越大,电阻率指数受含水饱和度的影响越小;随着孔隙度增大,饱和地层水岩石电阻率降低,随着地层水电导率的升高,孔隙度带来的影响增大;采用混合润湿表征岩石润湿性更加贴近客观规律,润湿性对于岩石整体电导率的影响是不可忽略的,也呈现了非常明显的非阿尔齐现象.模拟结果揭示了蚀变火成岩导电机理的复杂性,蚀变黏土附加导电性和对岩石导电路径的桥梁作用是造成南海某区块储层低阻的主要因素之一,为后续蚀变火成岩导电模型和饱和度方程的建模提供了可靠的依据.