基于数字岩心的致密砂岩储层孔隙结构与渗流机理

汪新光，郇金来，彭小东，张 冲，袁 雯，王 寅

（1.中海石油（中国）有限公司湛江分公司，广东湛江 524057；2.科吉思石油技术咨询（北京）有限公司，北京 100176）

南海西部海域目前已发现大量高温高压致密 砂岩气藏，作为一种天然的多孔介质气藏，其储层本身具有低孔、低—特低渗透的特征，给这类气田的商业开发带来了巨大挑战。储层孔隙结构是微观岩石物理研究的核心内容［1］，对提高气藏产能和采收率有着积极的作用，而致密砂岩储层在经历构造、沉积、成岩等作用过程中形成更为复杂的微观孔隙结构。随着大量致密油气藏的发现和开发，众多学者借助不同实验方法，从微观角度对致密砂岩储层孔隙结构与渗流规律进行了研究［2-15］，但传统的压汞、铸体薄片、扫描电镜等评价手段难以实现对这类储层内部结构的精细表征，而且经过压汞实验的岩样无法二次利用，对海上高成本的取心样品损耗很大。此外，基于达西定律所建立的渗流方程，只能对流体的运动特征进行宏观表征，无法对岩心孔隙结构（微米尺度）进行精细刻画［16］，更无法实现岩心微-纳米尺度的微观渗流机理研究。近年来快速发展的数字岩心技术［17-21］，在一定程度上实现了对储层内部微观结构的定量表征，尤其在微观层面对低渗透致密砂岩储层内部不同尺度的孔喉空间分布特征有着更为细致的刻画，从而进一步揭示了这类储集体岩石物理属性的变化规律和微观渗流机理。

笔者以莺歌海盆地乐东10 区高温高压气田致密砂岩储层为研究对象，开展岩石组构定量分析，借助多尺度Micro-CT 扫描技术建立三维孔隙结构模型，实现对致密砂岩储层数字岩心的两相渗流过程模拟，明确致密砂岩储层微观渗流机理，以期为乐东10区高温高压气田的合理开发提供理论基础。

1 区域地质概况

莺歌海盆地总体形态呈北北西走向的长条形，由东南部的莺歌海凹陷和西北部的河内凹陷组成，两凹陷间为临高凸起带分隔。根据构造变形特征，进一步将莺歌海盆地划分为莺东斜坡带、临高凸起带、中央坳陷和莺西斜坡带4 个一级构造单元（图1）。乐东10区气田位于莺歌海盆地莺东斜坡，主要含气层位黄流组为大型水道中细砂岩储层，砂岩厚度大、分选好。岩心孔隙度主要为7.7%～14.5%，中值为10%；渗透率主要为0.05～10.58 mD，中值为0.51 mD，整体属于特低孔、特低渗透储层，局部发育物性相对较好的“甜点”储层。气藏埋深为3 916.0～4 322.3 m，地层压力系数高，为2.174～2.306，原始地层温度为181.25～192.45 ℃，为高温超高压弹性气驱岩性气藏。

图1 乐东10区构造位置Fig.1 Structure location at study area

2 岩石组构定量分析

2.1 矿物组分特征

选取研究区2 口井7 块致密砂岩样品进行数字岩心分析，样品的常规物性实验测试显示，孔隙度为7.7%～11.3%，渗透率为0.059～0.268 mD。采用纳米级扫描电镜能谱（EDS）和矿物标准库的岩石定量测试系统进行矿物组分及含量测试，并与具备更高分辨率且可表征密度分布的背散射电子成像（BSE）对照，可获得精度远高于EDS 分辨率的矿物分布，准确区分粒径为几十纳米的矿物。测试结果表明（表1），岩石类型为长石质石英砂岩，主要矿物成分为石英、长石和碳酸盐与黏土胶结物，构成了致密砂岩孔隙空间发育的物质基础。长石以较为稳定的钾长石、钠长石为主，基性斜长石消耗殆尽。碳酸盐胶结物含量较高，主要发育方解石、白云石和铁白云石。黏土矿物含量较低，以伊利石为主，另有少量绿泥石、高岭石和蒙皂石。碎屑颗粒之间呈线-凹凸接触。

表1 测试样品基本信息Table1 Basic information of test samples

2.2 孔隙分布特征

基于样品背散射电子成像下的岩石矿物与孔隙发育特征定量分析，乐东10区致密砂岩储层孔隙空间发育程度主要受压实作用、胶结作用和溶蚀作用等成岩作用影响，高温高压环境对孔隙发育存在一定的改善。整体上储层孔隙空间以次生孔隙为主，原生孔隙为辅，其中次生孔隙以长石粒内溶孔和铸膜孔为主，原生孔隙受强烈胶结作用的影响所占比例较低。粒间溶孔与铸膜孔孔隙尺寸最大，这部分孔隙整体分布较为分散，无法形成有效的孔隙网络（图2a）；长石颗粒内部发育较多溶解作用形成的粒内溶孔，保存较完整的长石颗粒中可见解理缝发育，部分解理缝内部有其他矿物填充（图2b，2c）；胶结物主要为方解石和白云石，其内部发育较多晶间孔（图2d）；黏土矿物及云母主要作为粒内/粒间充填矿物存在（图2e，2f），由于黏土矿物间的转换作用以及孔隙压力差异，黏土矿物及云母内晶间孔尺寸差异较大。

图2 背散射电子成像下样品孔隙空间发育特征Fig.2 Characteristics of pore space in BSE imaging

基于对扫描电镜图像的二值化处理，提取出样品中的孔隙并进行定量分析。从孔隙半径分布（图3）可以看出，孔隙半径分布曲线呈现双峰或三峰的特征，说明孔隙发育呈现多尺度特征。其中粒间溶孔及铸膜孔孔隙半径主要集中在40～300 μm 左右，这部分孔隙尺寸较大，是孔隙发育的主要构成部分。部分残余粒间孔、粒内溶孔及胶结物溶孔的孔隙半径主要集中在3～40 μm，为次级孔隙发育的主要范围。胶结物内的晶间孔孔隙半径最小，这部分孔隙分布较广，且数量较多，其孔隙半径主要集中在0.070～3 μm。孔隙结构整体上表现为大孔细喉的特征，孔隙半径中等、喉道细，孔隙连通性较差。

图3 不同类型孔隙对应的半径分布范围Fig.3 Distribution of different types of pore radius

3 基于数字岩心的孔隙结构模型建立

数字岩心技术是以真实岩样为基础，通过一系列的图像处理技术和数值算法将岩心数字化［22］，从而在孔隙尺度上较真实地反映岩石微观孔隙结构。目前，业界常用的数字岩心构建方法主要有物理实验方法和数值重建方法［18，23-26］，前者是借助高精度光学仪器直接获取岩心的平面图像，结果较为准确、直观，但成本高、耗时长，常用的有CT扫描法、聚焦离子束成像法和序列成像组合法等；后者是借助岩心二维切片，利用各种数学算法间接实现基于数字岩心的孔隙结构三维构建，主要包含随机法和过程模拟法，这种方法简便易行、建模速度快，但得到的孔隙结构参数与实际岩心往往存在一定的差异。

本次研究采用CT 扫描法，利用FEI Heliscan 微米CT采集预处理后柱塞样品的真实三维投影数据，选取合适的图像重建方法将投影数据重建岩心灰度图像，最后采用图像二值分割方法分离灰度图像中的孔隙空间和岩石骨架后得到岩心内部孔隙结构的真实数字化模型。考虑到样品中孔隙发育的多尺度特征，采用2 种不同扫描精度对样品内部不同尺寸范围的孔隙进行捕捉。分别以4.7 μm 扫描精度对直径为8 mm 的岩心进行图像扫描，以0.65 μm扫描精度对直径为2 mm的岩心进行图像扫描。

3.1 CT扫描数据处理

图像扫描及滤波处理 结合孔隙半径分布特征，首先对直径为8 mm 的样品进行4.7 μm 尺度的扫描得到岩心整体三维图像，用于捕捉孔隙半径较大的孔隙信息。基于对图像的非均质性分析，优选次级孔隙较发育的区域钻取直径为2 mm 的小柱塞并进行0.65 μm 尺度的高分辨率精细扫描，用于捕捉孔隙半径中等的残余粒间孔、粒内溶孔及胶结物晶间孔。经过扫描获取岩心的三维灰度图像后，应用高斯及中值滤波对图像进行优化处理，便于减少扫描系统背景噪声的干扰，提高灰度图像的质量与对比度，方便下一步的图像分割与孔隙骨架提取（图4）。

图4 滤波处理前后灰度图像对比Fig.4 Comparison of grey images before and after filtering

图像二值化处理 二值化处理是一个极为重要的图像前期处理过程，核心是通过阈值的确定来提取孔隙空间与岩石骨架信息，从而分割出图像中的孔隙部分［27-28］。以岩心实验测量孔隙度为依据，通过选择合适的灰度阈值，将岩心CT灰度图像中的孔隙与岩石骨架进行分割，使分割的岩心图像中孔隙像素数占总像素数的比例等于岩心孔隙度，确保了分割结果的可靠性（图5）。分割阈值计算式为：

图5 岩心CT图像二值化处理前后对比Fig.5 Comparison of core CT images before and after binarization

表征单元体提取 数字岩心表征性体积单元（简称表征单元体，记为REV）是岩石相关特性趋于基本稳定时的最小体积。当岩石体积小于REV 时，其特性随体积变化而变化，表现出随机波动特性。当岩石体积大于或等于REV 时，其特性趋于稳定，不会发生明显变化。同时由于岩心CT图像体积较大，为提高计算效率，通常从图像中截取代表性较强的一块区域进行研究。因此需要通过REV 分析确定对不同孔隙结构具有代表性的最小样品或模型尺寸。

由于不同扫描方法的视域范围和扫描精度不同，每种图像所对应的孔隙尺度也有所差异。SEM图像视域范围大，扫描精度高，可以表征所有尺度的孔隙；由于分辨率的限制，直径为8 mm 柱塞样品主要捕捉粒间溶孔和铸膜孔等较大孔隙；2 mm柱塞样品扫描精度较高，而视域范围较小，主要捕捉长石溶孔、岩屑粒间孔等较小孔隙。通过对3 种扫描方式所获取的图像不同体元/面元大小下的孔隙度误差计算，可以确定微观研究的表征单元体边长。本次研究设定误差值小于5%时样品孔隙度趋于稳定，误差分析结果表明，扫描电镜粗扫柱塞样品的直径大于3 mm 时，所有尺度的孔隙度基本保持稳定，基于8 与2 mm 柱塞样品的REV 分析，其岩心尺寸分别为1.5 与1 mm 左右，对应的孔隙度基本保持稳定，据此确定8 mm柱塞样品扫描建立的六面体数字岩心模型每条边长像素数大于320，2 mm 柱塞样品扫描建立的六面体数字岩心模型每条边长像素数大于1 539。

3.2 三维孔隙结构模型建立

针对研究区致密砂岩样品孔隙呈现多尺度发育特征的情况，选取2 级扫描精度对样品孔隙结构进行获取，并利用分水岭分割方法，构建2种尺度的三维孔隙结构模型（图6），其体素分别为4.7 和0.65 μm。应用连通域分析方法对三维孔隙结构模型的连通性进行分析，0.65 μm 尺度模型可以捕捉到主要的渗流通道和相应的连通孔隙空间，而4.7 μm 尺度模型捕捉的主要为相互非连通的较大一级别孔隙，该级别的孔隙发育空间对孔隙度贡献较大，对渗透率影响较小。

图6 2种尺度三维孔隙结构模型Fig.6 3D graphs of pore structures for two scales

基于2 种尺度的孔隙结构模型，利用双尺度孔隙耦合的方法计算2种尺度孔隙结构模型中孔隙同时参与的情况下岩心整体的孔渗参数。双尺度孔隙耦合方法首先对大柱塞图像进行三值化分割，分别为孔隙、包含次级孔隙的区域、不含孔隙的骨架区域，而后提取孔隙网络时将包含次级孔隙的区域作为特殊的等效喉道处理，其孔隙度及渗透率可由精细模型（0.65 μm 尺度）计算得到，而等效喉道的进出口端流量变化可由达西公式计算得到，最后根据侵入-逾渗理论对孔隙结构模型进行求解，得到模型的整体孔隙度及渗透率［29］。

基于单尺度孔隙结构模型进行计算后得到4.7 μm 尺度孔隙结构模型孔隙度为5.1%，渗透率为0，即模型不连通；0.65 μm 尺度孔隙结构模型孔隙度为4.79%，渗透率为0.051 mD。利用双尺度孔隙耦合的方法计算得到岩心整体孔隙度为9.65%，渗透率为0.097 mD。该结果与实验室测试结果（孔隙度为11.33%，渗透率为0.10 mD）吻合程度较高。

为进一步验证粒间孔、铸膜孔等微米-亚毫米级孔隙对岩石孔渗参数的影响，通过减少4.7 μm 尺度孔隙结构真实模型中粒间孔、铸膜孔等微米-亚毫米级的孔隙数量，得到孔隙度分别为1.25%，2.89%，3.5%的3 个模型（表2，图7）。在0.65 μm 尺度孔隙结构模型不变的情况下，当双尺度融合孔隙结构模型中大孔隙数量上升时，整体孔隙度上升明显，但整体渗透率无明显变化。进一步说明粒间孔、铸膜孔等微米-亚毫米级孔隙对孔隙度影响明显，但由于自身无法形成有效渗流网络，对整体渗透率无明显改善效果。

图7 4套双尺度融合孔隙结构模型中的大尺度孔隙空间Fig.7 Large scale pore space for four two-scale fusion pore structure models

表2 4套双尺度融合孔隙结构模型计算的孔隙度和渗透率参数Table2 Porosity and permeability parameters calculated by four two-scale fusion pore structure models

4 渗流机理

4.1 基于孔隙结构模型的气水两相流模拟

基于数字岩心的孔隙结构模型是利用简单几何体所组成的空间网络来代表真实的复杂多孔介质的孔隙空间，利用这些简单几何体等效为真实多孔介质的孔隙和喉道，并通过拟静态流动模型模拟流体在这些简单几何体中的流动过程，以近似代替流体在真实多孔介质中的流动［30-33］。

利用已建立的孔隙结构模型进行两相流过程渗流模拟，从模拟结果可以看出，当只存在单尺度的微孔隙时（图8a），由于气相作为非润湿相，在毛细管力的限制下优先进入尺寸较大的喉道，占据了有利的渗流通道，从而导致水相渗透率迅速下降；而当模型中同时存在大孔隙及微孔隙时（图8b），气相克服最粗喉道进入孔隙，首先开始充填附近连接的大孔隙空间。由于大孔隙体积较大，延缓了气相侵入喉道的时间，因此等渗点的含水饱和度向左偏移，水相渗透率下降速度明显缓于单尺度模型的。

图8 不同尺度孔隙结构模型气水两相流模拟结果Fig.8 Simulation results of gas-water two-phase flow for pore structure models at different scales

4.2 微观孔喉发育模式与渗流机理

基于所构建的2种尺度三维孔隙结构模型连通性分析，大孔隙尺度模型中的粒间溶孔及铸膜孔互不连通，无法仅通过这部分孔隙形成有效渗流通道，而微孔隙尺度模型中的长石溶孔、岩屑粒间孔为主的微孔隙是样品的主要渗流通道。

小孔小喉渗流模式 在以长石溶孔、岩屑粒间孔为主的微孔隙尺度下，孔隙与喉道半径差异不大，渗流空间由微孔隙和喉道构成，渗流模式呈现小孔小喉特征（图9a）。该模式中，孔隙连通率对残余水饱和度影响较大。在孔隙结构相近的情况下，随着孔隙连通率升高，残余水饱和度逐渐降低。孔隙连通率低，表明死孔隙或孔隙盲端较多，而此类孔隙结构无法形成渗流通道，根据孔隙形态不同，残余水以水柱或水珠形式封闭于死孔隙或盲端中，使得残余水饱和度变大。在气水两相流残余水状态下，强表面张力和高毛细管力共同作用下以厚水膜形式赋存于微孔隙及喉道壁上，残余水饱和度较高。而极细喉道在高毛细管力作用下形成卡断，残余水以水柱形式堵塞喉道，残余水饱和度极高。

图9 不同孔隙结构模型下气水两相流模式Fig.9 Gas-water two-phase flow modes for different pore structure models

大孔小喉渗流模式 在双尺度融合后的孔隙结构模型中，存在粒间溶孔、铸膜孔等大孔隙，渗流模式呈大孔小喉特征（图9b）。在大孔隙参与下，气水渗流模式较小孔小喉模式更为复杂。大孔隙作为占比较高的储集空间，在驱替阶段对含水饱和度影响较大。在驱替开始阶段，部分气相侵入孔隙空间，大孔隙中部分含水，残余水饱和度较高；大孔隙提供的储集空间减缓了气体侵入的进程，微孔隙及喉道部分仍然保持了水相的连通，极细喉道在高毛细管力作用下基本未被气体侵入。在残余水状态下，大孔隙中的残余水在表面张力作用下以薄水膜形式赋存于孔隙表面，残余水饱和度小；微孔隙及喉道中的残余水在强表面张力和高毛细管力共同作用下以厚水膜形式赋存于微孔隙及喉道壁上，残余水饱和度较高；极细喉道中的残余水在高毛细管力作用下形成卡断，残余水以水柱形式堵塞喉道，残余水饱和度极高。

5 结论

乐东10 区致密砂岩储层主要发育3 种级别的孔隙，粒间溶孔及铸膜孔的孔隙半径最大，主要集中在40～300 μm左右；部分残余粒间孔、粒内溶孔及胶结物溶孔孔隙半径中等，主要集中在3～40 μm；胶结物内的晶间孔孔隙半径最小，主要集中在0.070～3 μm。

构建4.7 和0.65 μm 尺度的三维孔隙结构模型，4.7 μm 尺度模型反映相互非连通的较大一级孔隙，对孔隙度贡献较大，对渗透率影响较小；0.65 μm 尺度模型反映主要的渗流通道和相应的连通孔隙空间。利用双尺度孔隙耦合方法进行三维孔隙结构模型的渗透率计算，模型计算结果与实验测试渗透率结果基本吻合。

以长石溶孔、岩屑粒间孔为主的微孔隙发育，岩心表现为小孔小喉渗流模式；以粒间溶孔、铸膜孔等大孔隙及次级孔隙发育，岩心表现为大孔小喉渗流模式。大孔隙尺度中的粒间溶孔及铸膜孔互不连通，无法通过这部分孔隙形成有效渗流通道，而微孔隙尺度中的长石溶孔、岩屑粒间孔为主的微孔隙是样品的主要渗流通道。

符号解释

abs——k取不同值时对图像的孔隙度与岩心孔隙度的差取绝对值，%；

i——灰度值，无量纲；

Imax，Imin——图像的最大、最小灰度值，无量纲；

k——从Imin—Imax依次取得的孔隙分割阈值，无量纲；

k*——最优阈值，无量纲；

p（i）——灰度值为i的像素数，无量纲；

ϕ——岩心孔隙度，%。