裂缝性页岩油储层三维数字岩心电阻率模拟

赵岩龙，李星宇，方正魁，张永安，张明科，代 星

(1.中国石油大学(北京)克拉玛依校区 石油学院，新疆 克拉玛依 834000；2.中国石油天然气股份有限公司新疆油田分公司 石西油田作业区，新疆 克拉玛依 834000)

引 言

页岩油储层岩石通常具有较高脆性，在高强度的构造应力场作用下易产生大量尺度不同的裂缝[1]。准确表征页岩油储层中微裂缝的赋存规模及其相应的电学特性，对页岩油的高效开发具有重要意义。吉木萨尔凹陷二叠系芦草沟组页岩油储层中，其岩心孔隙度普遍低于10%，渗透率低于0.5×10-3μm2，微裂缝的存在连接了原本互不连通的孤立孔隙，且自身具有平均0.76%的孔隙度和30.04×10-3μm2的高渗透率决定了其对储层渗流能力的贡献不可或缺。微裂缝能够将导电介质串联起来，从而减小区域电阻率，尤其在高含水饱和度的储层中更加显著，需要在电阻率成像测井资料中加以考虑[2-4]。常规电阻率成像测井能够解释构造应力下的大型裂缝特征，但对于失水收缩、溶蚀等作用下岩性微缝表征难度较大，且常规岩石物理实验取样较困难，无法定量评价岩心中微裂缝参数对电阻率的影响。基于数字岩心的数值模拟方法能够解决上述难题。

数字岩心技术[5-6]是一种基于岩心结构信息、从微观尺度出发分析岩石宏观物理属性的新兴手段。近十年来国内外发展出了众多的数字岩心建模方法，为基于数字岩心的精细建模与数值模拟提供了条件[7-8]。四参数随机生长法(QSGS,Quarter-Sweep Gauss-Seidel)属于数字岩心完全随机法的一种，是较为成熟的一种数值重构方法，在对各向异性岩石的模拟上具有优势，能够生成一定孔隙度条件下的强各向异性数字岩心[9]。本文在岩心连通性与分选性方面对QSGS方法作出改进并进行应用。孔强夫等[10]总结了数字岩心电性数值模拟方法与发展方向，指出有限元方法具有模拟精度高，对不规则区域适应强，能够考虑各向异性且便于处理复杂边界的优势。李潮流等[11]采用有限元法对致密砂岩电阻率进行了数值模拟计算，模拟结果在高含水区间与驱替实验结果基本一致。张龙海等[12]分析了研究地区低孔低渗储集层的孔隙结构特征及其对储集层产能和电性的影响。刘学锋等[13]利用有限元方法，基于三维数字岩心研究了油层岩石的电阻率特性。本文采用有限元方法及数字岩心技术，对微尺度条件下不同的润湿性以及不同裂缝排列方式对页岩油储层电性的影响进行了探讨。

1 FIB-SEM扫描图像分析与数字岩心重构

选取吉木萨尔凹陷二叠系芦草沟组页岩油储层取心岩样，利用蔡司聚焦离子束扫描电子显微镜(FIB-SEM，Crossbeam 550)扫描得到原始岩心序列切片图像(图1(a)为典型图像)。通过校正消除了切割扫描过程中相邻图像间的错位并叠加成像，通过中值滤波去除了部分图像噪声，如图1(b)所示。叠加得到的真实岩心三维数据体(体素数1 723×782×281)为后续数值重建环节提供了孔喉参数依据。采用赵岩龙等[14]提出的基于改进BP神经网络的岩心图像分割方法，以各像素点灰度值、色值分量、局部领域灰度均值组成输入向量，隐藏层节点数设为7，输出层节点为岩石各组分，并分别基于连通领域、局部领域、边缘领域对训练样本进行了修正，分割结果如图1(c)所示。得出的实际岩样的基本构成为：方解石及石英基质(92.35%，白色)、有机质(6.09%，红色)、黄铁矿(0.033 6%，黄色)及孔隙(1.52%，黑色)。由于岩样中黄铁矿含量极低，本文将其与方解石视为同一种组分考虑。

图1 岩心图像处理与BP神经网络分割

QSGS方法在三维情况下的核心重构思路为：先在初始化空矩阵中随机布置一定数量的种子点，再根据原岩样统计得到的孔喉长宽比、迂曲度等形态学参数设定8个相邻生长方向的生长概率，最终以孔隙度作为终止条件。其中，长宽比及迂曲度参数通过种子生长算法计算得出，具体概括为以孔隙内种子点为中心，遍历搜索并储存邻域内为孔隙点的体素坐标，计算三坐标全距之比即得到长宽比，同时得到各孔隙区域的体积。因技术原因本文没有进行迂曲度的计算，转为采用蒙特卡洛法对参数进行充分随机，选取最符合原样本长宽比的重构结果下的参数进行后续模拟。在QSGS算法的应用过程中，本文针对页岩的低孔渗特点，在原有方法基础上进行了两处改进，以保证低孔隙度条件下孔喉结构的非均质性以及连通性：①每次生长结束时以1体素为步长转移种子点坐标后进行下一次生长；②以实时孔隙度的对数曲线为约束规定不同生长程度下参与生长的种子点数目。

天然裂缝无论大小均存在具有自相似性的形状和构造，任意尺度上局部能够包含整体，即分形特征，可采用分形维数D量化[15]。分数布朗运动(FBM,Fractional Brownian Motion)是一种时间序列下具有长记忆性以及多尺度自相似特征的随机行走过程，其二维形式具有空间上的分形特性，可用于文中粗糙裂缝面的重构。其壁面Z坐标满足下式特征[16]：

E[ZH(X+h)-ZH(X)2]=|h|2Hσ2。

(1)

其中，Z为二维分数布朗运动过程Z轴的坐标值；E为数学期望；H为Hurst指数(0

将重构得到的三维裂缝截选50个切面，计算其二维粗糙度方差为0.284，均值为2.821，距实际截取的真实值3.027较为接近。将重构页岩数字岩心与裂缝按照等尺度叠加，经Tecplot可视化得到的结果如图2(a)所示。自相关函数能够表征二维至三维图像的内在结构关联，而局部孔隙度分布函数描述了岩心中的孔隙分选性[17]。为了验证该岩心的准确性，本文分别计算了原岩样与重构岩心的自相关函数与局部孔隙度分布函数(L为统计尺寸/体素)。在本文假设条件中，岩石基质骨架、油气组分仅有微弱的电导率，孔隙中电解质溶液是主要的导电介质，电导率为方解石、石英、黏土矿物的104～109倍，对孔隙结构的重构准确性至关重要[1]。图2(b)以及图2(c)展示的对比结果验证了该改进方法的可行性。

图2 重构得到的三维数字岩心及自相关函数验证

2 基于体素单元的有限元模型与验证

目前，国内针对页岩油储层不同微裂缝密度下的岩石电学性质的研究较少，本文使用的有限元法十分便于处理六面体体素组成的数字岩心下的微观结构。在现有网格条件下，根据变分原理，给定流体组分和固相各组分的电导率张量后，用变分法将电场问题转换为求解数字岩心整体能量泛函数极小值的状态，该能量状态对应的各节点电压即为待求解。三维数字岩心电能(能量)方程表示为：

(2)

其中，ep、eq表示电场强度e在x、y方向上的分量，σpq表示xoy面法向上的电导率张量，取决于该体素表示的组分。r表示节点编号，每个六面体单元拥有编号为1～8的节点。

电压函数U表示为：

U(x,y,z)=Nrur。

(3)

其中，ur表示处于第r个节点位置的节点电压，N为形状矩阵。对于单个体素点，其电场e在空间中具有3个分量，其中方向为p的电场分量表示为：

(4)

将(3)式代入(4)式中，电场分量变为：

(5)

表1 矩阵npr的计算公式

将(5)式代入能量计算公式(2)中得：

(6)

(7)

其中，三维情况下，Drs为有限元方法中规格为8行8列的劲度矩阵，整体劲度矩阵由单元劲度矩阵组装而成。在单元结点上对能量函数进行插值，将拉普拉斯方程中各变量用相应的插值函数组成线性表达式，最终将所有体素点的能量En求和可得整个系统所耗的电能。岩心的视电阻率由岩心电流方向进出口端电压之差除以电流密度得到。本文使用Python语言feon模块完成上述过程。

为了验证该方法的准确性，本文使用含水饱和度为100%的页岩数字岩心，进行了电阻率模拟测试。结果显示，在纯孔隙介质情况下，充满电阻率为20 Ω·m电解质溶液的不同孔隙度数字岩心电阻率R0与电解质溶液本身的电阻率Rw的比值为一常数，符合阿尔奇特征。经计算得到的地层因素F=68，岩石胶结系数m=1.012，与某页岩油区块实验室实测数据F=74，m=1.114的数量级基本相符[15]。

为了方便直观了解电压在数字岩心中的传播规律及特点，下面给出了一次模拟中某二维截面的节点电压分布云图，如图3所示，其中电压单位为毫伏。该图可以表明，电压沿裂缝传导时的衰减量较少，与理论相符，进而验证了该方法的可行性。

图3 截面电压分布云图

3 微裂缝参数对页岩电阻率的影响规律

本文基于上述模型，对相同矿物组成下含有不同开度、角度、排列方式微裂缝以及不同含水饱和度、润湿性的三维数字岩心分别进行了电阻率模拟。

设定裂缝开度为250 nm(比例尺为1体素∶50 nm)，讨论主视轴为X轴时裂缝角度分别为0°、15°、30°、45°、60°、75°、90°下的三轴轴向视电阻率变化情况。模拟结果如图4(a)所示，当Y、Z方向裂缝角度由45°过渡至60°时，电阻率数值出现了突变。此现象可以通过电学基本原理解释，即作为主要导电通路的裂缝流体空间贯穿了电流流向前后截面时，电阻率最小的元件在电路中的参与形式由串联变为并联，一定程度上提升了对岩石整体的导电贡献能力。相同裂缝条件下Z方向电阻率数值略高于Y方向，是岩心自身存在的各向异性导致的。当裂缝角度一定时，电导率随裂缝开度的增大而接近线性减小，如图4(b)。上述变化规律的具体数值见表2。

图4 不同角度、开度裂缝的存在对岩心电阻率的影响

表2 不同夹角和开度下的电阻率变化率

不同含水饱和度下的油水分布由数学形态学方法[18]确定，其原理是通过膨胀、腐蚀运算完成三维图像的开运算与闭运算，为孔隙空间铺设水膜或油滴。图5展示了该方法模拟的相同岩心不同润导率随含水饱和度增大而递减得更加显著。而当岩湿性下的孔隙内油水分布。其中图5(a)含水饱和度为76.3%，图5(b)含水饱和度为78.7%。(图中黄色代表岩石基质，红色代表油，蓝色代表水)。

图5 两种润湿性岩心的内部结构示意图

本文按照三种裂缝排列方式(平行排布、正交排布、斜交排布)建立了三种类型的孔隙-裂缝双重介质，其各自截面形态如图6所示。模拟结果表明，含水饱和度越低，岩石电阻率随岩石含水饱和度的变化越迅速；相同含水饱和度下亲水岩石具有比亲油岩石更优异的导电性。当含水饱和度Sw≥20%时，采用指数函数式

图6 模拟的三种裂缝排列方式示意图

R0=a·exp(-bSw)

(8)

能够较好地拟合亲水或亲油状态三种裂缝排列方式下岩石电阻率与含水饱和度存在的对应关系，如图7所示。

图7 岩石电阻率与含水饱和度函数拟合曲线

式(8)中，a、b的数值与裂缝密度以及排列方式有关，最终拟合结果见表3。

表3 岩石电阻率与含水饱和度函数关系拟合结果

各拟合曲线相关系数R2均在0.98以上，能够说明拟合函数的可信度。比较三种裂缝排列方式下的参数a可以明显得到aA>aB>aC的规律。说明等量的岩石体积内沿电流方向发育的平行裂缝数量越多或沿电流方向发育的裂缝交错夹角越小，电石亲水时，参数a的骤减现象则明显强于亲油时的情况(亲油时递减率为0.119而亲水时为0.210)，这一定程度上反应了孔隙空间内部水膜的连续性。

电阻率测井或自然电位测井中，电流方向通常为垂直于井筒方向，亦是生产时的渗流方向。局部高电导率对应较好的连通程度和渗流能力，在极低孔隙度条件下意味着该方向微裂缝的发育。以本文结论为基础的升尺度缝网研究可以为电成像测井、储层评价与生产提供进一步的理论依据。

4 结 论

(1)改进的QSGS方法对低孔渗页岩岩心数值重构具有比传统方法更好的适用性；岩石自身的各向异性以及裂缝的开度、角度、发育方向的变化能够分别引起岩石某方向视电阻率的不同程度变化。

(2)当含水饱和度Sw≥20%时，亲水或亲油状态三种裂缝排列方式下岩石电阻率与含水饱和度存在一定的指数关系，其关键参数a与裂缝密度以及排列方式有关，具体表现为随同电流方向发育裂缝数量的增多而增大。

(3)同含水饱和度情况下亲水岩石因具有更连续的水膜而电阻率更小。当岩石亲水时，电阻率与含水饱和度函数关系拟合式中参数a的骤减比亲油状态下更显著。