基于数字岩心模拟的深层页岩气储层可压性评价模型

赖富强，刘粤蛟，张海杰，王海涛，毛海艳，张国统

(1.重庆科技学院复杂油田勘探开发重庆市重点实验室，重庆 401331；2.重庆页岩气勘探开发有限责任公司，重庆 401147)

深层页岩气储层一般是指深度大于3 500 m的页岩气储层[1]，由于沉积作用、成岩作用、埋深、温度及压力等因素影响，深层页岩气储层相较中浅层层理、断裂及天然微裂缝更为发育，应力差及岩石力学参数变化更大，孔隙结构更加复杂，导致储层的可压性主控因素不清[2-4]。而且由于页岩储层为多尺度、多组构储集空间，具有非均质性强、易风化等特点[5-7]，受限于岩心来源、收获率等因素，难以取得具有代表性的样品，开展传统岩石物理实验的难度较大，所以需要构建深层页岩气储层典型的多尺度、多组分三维数字岩心模型来进行岩石参数模拟[8-10]。国内外学者从岩石力学性质、矿物组分、孔隙结构、层理及裂缝发育和地质因素等方面对储层进行可压性评价分析[11-12]，但是由于储层各异，每种模型考虑的可压性影响因素存在差别，而且并未给出如何确定每种影响因素对于可压性变化的敏感性大小，主控因素权重不清楚，只是定性地规定了各影响因素之间的重要程度，所以构建的可压裂性评价模型的准确性和适用性较差[12]。笔者提出综合考虑研究区深层页岩气储层特点，由于没有页岩岩心三维X射线CT扫描，利用已有的二维岩心图片，利用马尔科夫链-蒙特卡洛法[14](MCMC法)构建多尺度、多组分三维数字岩心，基于数字岩心物理模型，对深层页岩气储层可压性主控因素采用有限元计算方法进行弹性参数数值模拟[15]，然后利用摩尔斯分类筛选法定量地分析裂缝倾角、裂缝密度、层理密度、裂缝长度、脆性矿物、层理倾角、孔隙半径、含气饱和度这8种影响因素引起弹性参数变化的敏感性相对大小。最后综合考虑水平应力差异系数的影响，利用层次分析法确定这9种影响因素(上述8种影响因素和水平应力差异系数)各自的权重，建立一种深层页岩储层的可压性评价计算模型。

1 区域概况

研究区为大足区块，位于渝西地区北部，地处四川盆地丘陵地带，东南部为低山，最高海拔934 m，中部为浅丘带坝，西部为深丘，最低海拔为185 m[16]。四川盆地大地构造分区包括川西北坳陷带、川中隆起带和川东南坳褶带(图1)，大足区块在构造上属于川东南坳褶带、川中隆起带交界[17]。区内稳定发育有埋深介于3 500～4 500 m的五峰组-龙马溪组页岩层，是页岩气开发的主要目标层序。

图1 区域构造位置示意图Fig.1 Schematic diagram of regional structure location

根据研究区域X1井、X2井、X3井、X5井的资料分析，龙马溪组发育有深灰色、黑灰色、灰黑色及黑色页岩，自上而下颜色逐渐变深，主要为黑色碳质页岩，有机质丰富，岩心上还有笔石富集与少量黄铁矿充填物。岩心整体页理较发育，少量局部发育。微裂缝较发育且少量微裂缝呈高角度分布。发育有大量的微米—纳米级无机孔隙和有机质孔隙。储层以石英和黏土矿物为主，两者之和的平均体积分数为40%，其中白云石、斜长石、黄铁矿、方解石等矿物仅占少量，地层非均质性较强，石英含量由下往上逐渐减少。

2 页岩气储层可压性模拟与主控因素敏感性

通过对研究区大量的测井资料对比分析认为，研究区可压性评价的关键参数包括了矿物组分、水平地应力差、层理特征、裂缝、孔隙结构，含油气饱和度等[11-13]。在没有岩心三维CT扫描信息的情况下，本文中针对研究区岩心铸体薄片(图2(a))进行二维扫描电镜与X射线衍射实验，对薄片矿物组分(石英、黏土、方解石、白云石、斜长石、黄铁矿)及有机质进行识别得到二维图片(图2(b))，对二维图像进行遍历扫描获得每个方向上的条件概率，然后依据MCMC法对3个相互垂直的二维图片(图2(c))进行三维重构[18]单组分数字岩心(图2(d)～(g))。

多组分三维数字岩心是在构建有机质、黏土、方解石、石英、白云石、斜长石、黄铁矿单组分数字岩心的基础上通过选取其中4种单组分数字岩心嵌套得来(图2(h))。多尺度多组分三维数字岩心(图2(j))是将高分辨率孔隙信息(图2(i))融合到低分辨率多组分三维数字岩心得来[19]。本文中所构建的数字岩心的分辨率为3 μm/像素，模型尺寸为240 μm×240 μm×240 μm。

图2 重构的多尺度多组分三维数字岩心Fig.2 Reconstructed multi-scale and multi-component 3D digital core

基于构建的多个不同4种成分的多组分数字岩心，即可研究不同矿物组分对于岩石弹性参数的影响；另外，为了考虑层理与裂缝的影响，本文中建立平板状的层理(图2(k))、平板状裂缝(图2(m))岩心模型，将其嵌入多组分数字岩心模型，从而建立层理发育(图2(l))和裂缝发育(图2(n)的三维数字岩心模型；为了考虑含油气饱和度的影响，对于多尺度多组分岩心模型(图2(j))中的孔隙空间(图2(i))进行气水两相流体驱替模拟获取不同气水分布岩心模型(图2(p))，研究气水饱和度对于岩石弹性参数的影响。

用MCMC法构建数字岩心的优势在于其根据原始图像的转移概率对新图像进行赋值重构，间接地考虑了岩心内部物质的空间分布，可以较好地表征实际岩心。

通过计算本文中所构建的数字岩心与二维图像的自相关函数，分析数字岩心模拟的平均孔径与岩心CT实验平均孔径的相对误差，验证了本文MCMC法构建的数字岩心可以较好地反映研究区实际页岩的孔隙结构和矿物组成(图3)。在建立可靠的数字岩心模型的基础上，利用有限元计算方法计算岩心弹性参数(图2(q))，采用摩尔斯分类筛选法计算各种可压性影响因素变化引起的弹性参数变化的相对大小(图2(r))，从而计算敏感性因子，确定各影响因素敏感性相对大小，为层次分析法确定权重及可压性评价模型的构建做准备(图2(s))。

图3 自相关分析与平均孔径对比Fig.3 Comparison between autocorrelation analysis and average aperture size

2.1 有限元法计算岩石弹性参数

给定数字岩心模拟时所考虑物质的弹性参数——体积模量和剪切模量[20-23]，依据有限元法即可计算出岩石体积模量和剪切模量[24]。有限元计算弹性参数的原理是：岩心离散化，对每个网格进行弹性参数赋值；在岩心X、Y、Z方向上分别施加一个宏观应力，依据应力应变关系确定每个单元格的机械能，然后依据变分原理即能量取极小值，依据共轭梯度法求解获取每个单元格上的位移分布，从而可以确定岩心总的位移，最终确定岩心应变。已知应力和应变，便可确定岩心的弹性体积模量和剪切模量(图4)。

图4 有限元法计算岩石弹性模量原理Fig.4 Principle of calculating rock elastic modulus by finite element method

2.2 有限元法计算岩石弹性参数验证

首先利用理论值进行验证，将有限元计算的弹性模量，与V-R-H(Voight-Reuss-Hill)理论计算的等效弹性模量[25](理论值)进行对比。

Voight上限值：

(1)

Reuss下限值：

(2)

V-R-H值：

MVRH=(MV+MR)/2 .

(3)

式中,N为岩石中的组分数；fi为第i种组分的体积分数；Mi为第i种组分的体积模量或剪切模量，GPa。

从图5中可以看出，模拟值均落在理论值上下限值之间，说明通过数字岩心有限元法模拟出来的弹性模量值可靠，且理论计算弹性模量与有限元计算弹性模量的平均相对误差为：体积模量5.55%，剪切模量2.60%。

图5 理论计算弹性模量与有限元计算弹性模量Fig.5 Elastic modulus by theoretical calculation and finite element calculation

然后用试验参数值进行验证，用研究区储层4 037～4 038 m深度的12块岩心的多组分岩心模拟结果与对应深度岩心的三轴岩石力学参数试验测的弹性模量、泊松比以及计算的脆性指数进行对比(图6)。结果显示三轴岩石力学参数实验测的弹性参数与有限元计算结果的平均相对误差为:弹性模量8.9%，泊松比5.6%，脆性指数7.8%。通过理论值与试验值的验证，说明基于构建三维数字岩心物理模型，利用有限元法计算岩石弹性参数的结果可靠，可开展下一步弹性参数影响因素敏感性分析工作。

图6 试验测试值与有限元计算值的平均相对误差Fig.6 Average relative error between experimental test value and finite element calculation value

2.3 页岩弹性参数影响因素敏感性

由数字岩心岩石物理模拟出体积模量K与剪切模量G。由下式计算出弹性模量E和泊松比μ:

E=9KG/(3K+G),

(4)

μ=(3K-2G)/(6K+2G).

(5)

采用摩尔斯分类筛选法[26]，定量模拟影响因素在不同数值情况下弹性参数的变化大小定量分析影响因素的敏感性，弹性模量的模拟如图7所示。

采用自变量以固定步长变化，灵敏度判别因子取多个平均值，公式如下：

(6)

式中,S为灵敏度判别因子；n为模型运行次数；Yi、Yi+1分别为模型第i次、第i+1次运行输出值；Y0为计算结果初始值；pi、pi+1分别为第i次、第i+1次模型运算参数值相对初始参数值的变化百分率。

通过层理、脆性矿物、裂缝以及孔隙含气饱和度对于岩心弹性参数的模拟及敏感性分析结果(图7、8)可知：8种影响因素对于弹性参数的敏感性由大到小依次是裂缝倾角、裂缝密度、层理密度、裂缝长度、脆性矿物、层理倾角、孔隙半径、含气饱和度。

图7 影响因素变化引起弹性模量变化的模拟Fig.7 Simulation of Youngs modulus change caused by change of influencing factors

另外，在页岩储层实际压裂施工中裂缝在纵向上受最小水平主应力的影响而变化并沿着最大水平主应力方向进行延伸，即水平主应力差较小时，人工裂缝会更多的与天然裂缝沟通，表现出较好的可压裂性，因此本文中还把应力因素同上面8种因素进行考虑，但是由于应力很难用数字岩心来模拟，本文中对其敏感性的相对大小进行了9次排列，然后计算水平应力差异系数在不同排列情况下的可压性指数。经微地震监测和压裂施工曲线验证，水平应力差异系数排在脆性矿物与层理倾角之间时与微地震监测和压裂施工曲线的对应效果最好，故本文中定性地将水平应力差异系数对可压性的敏感性定在脆性矿物与层理倾角之间。水平应力差异系数计算公式为

Kh=(σH-σh)/σh.

(7)

式中，Kh为水平差应力系数；σH为水平最大主应力，MPa；σh为水平最小主应力，MPa。

图8 X3井岩心弹性参数敏感性Fig.8 Sensitivity of core elastic parameters of well X3

3 可压性评价模型建立

采用层次分析法综合裂缝倾角、裂缝密度、层理密度、裂缝长度、脆性矿物、水平地应力差、层理倾角、孔隙半径、油气饱和度建立一种能够对页岩气储层可压性进行全面科学评价的计算模型。为得到综合评价储层可压性的无量纲常数，还要对各参数采用经验赋值或极差变换来进行标准化、归一化处理。

(1)可压性评价参数标准化。可压性指数的影响指标包括正向指标、逆向指标，正向指标的值越大，逆向指标的值越小，可压裂性越好。正向指标Si和逆向指标Sj分别为

Si=(Xi-Xmin)/(Xmax-Xmin),

(8)

Sj=(Xmax-Xj)/(Xmax-Xmin).

(9)

式中，Xmax为参数最大值;Xmin为参数最小值；Xi、Xj为考虑参数的参数值。

(2)因素权重的确定及计算。利用层次分析法，根据参数敏感性的大小，确定它们之间重要性(表1)，并根据其重要性给每一层元素给出定量表征，构造出判断矩阵A，Aij表示元素i相对于元素j的重要程度值，Aji表示元素j相对于元素i的重要程度值。然后运用数学的方法求出判断矩阵的最大特征值所对应的特征向量,从而得到权重系数。

表1 判断矩阵元素的标度极其含义Table 1 Scale and meaning of judgment matrix elements

根据前文的影响因素敏感性分析，确定本文中所研究的9种影响因素对于弹性参数变化敏感性的大小，因此根据其对可压性的相对重要性进行取值可得到可压性指标的判断矩阵，如表2所示。

表2 判断矩阵元素的标度极其含义Table 2 Scale and meaning of judgment matrix elements

求得判断矩阵的特征向量为W=(0.710 9，0.506 4，0.354 1，0.244 9，0.168 2，0.115 4，0.079 7，0.056 3，0.041 8)，最大特征值λmax=9.401 4，再对特征向量作归一化处理，便可得到各参数的权重系数为0.312、0.222、0.156、0.108、0.074、0.051、0.035、0.025、0.018。故最终可建立可压性评价模型。裂缝倾角、层理倾角、裂缝密度、裂缝长度是通过对成像测井中的裂缝和层理进行人机交互识别得到的，裂缝密度为单位长度内裂缝发育的条数，层理密度通过成像测井纹层识别后用纹层密度来表征，纹层密度为单位长度内纹层的个数；孔隙半径从核磁共振测井中获得，主要利用其横向弛豫时间(T2谱)来反映岩石的孔径分布，即岩石在饱和流体状态下T2谱越长则岩石孔隙度越大；脆性矿物含量通过岩心试验分析得到。可压性评价模型为

FI=0.312Lf,dip+0.222Lfd+0.156Cld+0.108LfL+

0.074Kw+0.051Kh+0.035Cl,dip+0.025Kj+0.018Bhd.

(10)

式中，Lf,dip为归一化裂缝倾角；Lfd为归一化裂缝密度；Cld为归一化层理密度；LfL为归一化裂缝长度；Kw为归一化脆性矿物；Kh为归一化水平应力差异系数；Cl,dip为归一化层理倾角；Kj为归一化孔隙半径；Bhd为归一化含气饱和度。

4 可压性测井评价应用

根据以上对可压性的研究内容，利用本文中所建立的可压性计算模型对X1、X3井进行可压性系数计算，结果如图9、10所示。

图9 X1井可压性指数计算结果Fig.9 Calculation diagram of compressibility index of well X1

由于成像测井的分辨率很大，所以对可压性特征明显的部分做了截取展示(图9(b)、图10(b))。从计算的可压性曲线FI可知，X1井龙-14小层到五峰组整体的可压性指数都较高，范围为0.42～0.69，平均为0.58，龙-11小层到五峰组顶部最高(图9)；而X3井在龙-11小层和五峰组的可压性指数相对较高(X3井没有核磁测井数据，没有考虑孔径参数的影响)，范围为0.42～0.69，平均为0.58(图10)。两口井在龙-11小层和五峰组的平均可压性指数为0.56。根据Rickman等[27]的可压性判断标准，综合考虑X1、X3井的可压性指数，认为研究区五峰组-龙-11储层可压性特征相对较好。

图10 X3井可压性指数计算结果Fig.10 Calculation diagram of compressibility index of well X3

利用本文模型计算X3井直改平井段的可压性指数，首先发现7、11、14段的可压性指数较高，5、21段的可压性指数较低。从可压性指数较高的第7段的加沙压裂施工曲线来看(图11(a)、(c))，第7段180 min前有多处裂缝破裂显示(有破裂显示：排量不变，油压上升又下降；油压迅速下降，排量上升；泵压不变，排量上升。180 min后砂堵进行冲砂，不做分析)；而可压性指数较低的第21段并没有什么裂缝破裂显示(无破裂显示：泵压随排量的上升而上升)[28]，油压曲线和排量曲线都没有变化(图11(b)、(c))。各段所对应的微地震监测SRV(页岩储层改造体积)值也表现出同样的高低特征(图12)。说明本文模型所计算的可压性指数可以与实际生产所测参数对应，模型具有一定的可靠性，可以较好地为页岩气储层射孔提供参考并促进页岩气开发。

图11 X3井第7段及第21段压裂施工曲线Fig.11 Fracturing operation curve of section 7th and 21th of well X3

图12 X3井(直改平)储层综合分析Fig.12 Comprehensive analysis of reservoir in well X3 (straight to flat)

5 结 论

(1)通过数字岩心模拟，可以定量地分析弹性参数各影响因素的敏感性相对大小，可以为层次分析法建立可压性评价模型时更准确地确定各参数之间的权重。

(2)经微地震监测和压裂施工曲线验证，基于数字岩心模拟的可压性评价模型可为深层页岩气储层靶点层位选取和压裂工艺优化提供指导。但是本次研究在数字岩心模拟层理与裂缝的时候是通过构建一个同分辨率的平板来实现，较为简单，不能准确地反应真实的岩心层理、裂缝发育情况，需要后期研究进一步改进。