盛 军，阳 成，徐 立,丁晓军，杨晓菁，李 纲

(1.中国石油青海油田分公司 勘探开发研究院，甘肃 敦煌 736202； 2.中国石油青海油田分公司 采油一厂,甘肃 敦煌 736202)

引 言

储层微观渗流特征研究是储层研究中的一项重要内容，从微观的角度研究储层孔隙空间中各相流体的渗流特征，对掌握开发过程中的油水矛盾、水驱特征等方面具有重要意义[1-5]。目前实验室内研究储层微观渗流的评价方法，比如相渗实验、核磁工作实验、驱油效率实验等已经日趋成熟，但这些方法多少都存在一些问题，例如样品制备繁琐，对岩心样品损耗严重，对于低渗透岩心驱替时间过长等。近年来，随着X-CT扫描技术的发展，其在石油行业的运用也愈加广泛，利用多种尺度X-CT扫描融合建立的数字岩心，在研究储层微观渗流特征方面具有重要意义[6-13]。

本文以国内某油田陆相致密储层样品为研究对象，借助多尺度X-CT扫描技术建立样品数字岩心模型，实现对储层三维孔隙空间网络的提取，模拟两相渗流过程，并得到渗流特征参数。

1 多级融合数字岩心的建立

多级融合数字岩心的建立首先依靠不同尺度的X-CT扫描获取高质量的二维灰度原始图像，利用专业软件的图像分割技术，对重构出的三维灰度图像进行二值化分割，划分出孔隙与颗粒基质，从而获取岩石骨架及孔隙的基本信息(图1)。通过多级尺度孔隙融合技术，将纳米级高分辨率CT图像中提取的孔隙网络结构有机融合到微米级别的孔隙网络中，进一步建立数字图像模型以及孔隙网络模型，多级融合技术改进了原先的微米级别孔隙网络中无法辨别微小孔隙的问题，使得建立的数字岩心孔隙网络更加接近真实岩心。

图1 二值化处理过程Fig.1 Binaryzation process

2 孔隙网络提取

孔隙网络模型建立，是指通过某种特定的算法，从二值化的三维岩心图像中提取出结构化的孔隙和喉道模型，同时该孔隙结构模型保持了原三维岩心图像的孔隙分布特征以及连通性特征。采用Silin、Blunt 和 Dong 所提出的最大球算法(Maxima-Ball)是目前三维数字岩心对孔隙网络模型构建的主要方法，运用最大球算法进行孔隙网络结构的提取与建模，不但能够提高网络提取的速度，而且同时可以保证三维空间中孔隙分布特征与连通特征的准确性[14-17]。

最大球算法实际上是将一系列不同尺寸的球体填充到三维岩心图像的孔隙空间中，而各尺寸填充球之间根据半径从大到小存在着连接关系(图2(a))。整个岩心内部孔隙结构将通过相互交叉链接的一系列球串来表征。孔隙网络结构中的“孔隙”和“喉道”的确立是通过在球串中寻找局部最大球与2个最大球之间的最小球，从而形成“孔隙-喉道-孔隙”的配对关系来完成。最终整个球串结构简化成为以“孔隙”和“喉道”为单元的孔隙网络结构模型。“喉道”是连接2个“孔隙”的单元；每个“孔隙”所连接的“喉道”数目，则称之为配位数。(图2(b),图2(c))

图2 最大球算法提取孔隙网络模型过程Fig.2 Extraction process of pore network model using maximum ball algorithm

在用最大球法提取孔隙网络结构的过程中，形状不规则的真实孔隙和喉道被规则的球形填充，进而简化成为孔隙网络模型中形状规则的孔隙和喉道。在这一过程中，利用形状因子G来存储不规则孔隙和喉道的形状特征。形状因子的定义为G=A/P2,其中A为孔隙的横截面积，P为孔隙横截面周长。在用规则几何体来代表岩心中的真实孔隙和喉道时，要求规则几何体的形状因子与孔隙和喉道的形状因子相等。尽管规则几何体在直观上与真实孔隙空间差异较大，但他们具备了孔隙空间的几何特征(图3)。

最大球算法在识别孔隙与喉道的过程中，由于算法自身的限制，认为2个最大球之间的连接都是喉道。这就导致一部分原本应该属于孔隙的空间被划分为喉道空间，使得孔隙空间相对与实际情况偏小，针对这个问题，Øren and Bakke[18]提出了校正孔隙与喉道的方法(图3)，校正后的孔隙与喉道长度分别为：

lt=ltt-lj-li；

(1)

(2)

(3)

式中，lt为喉道长度，lj为孔隙j长度，li为孔隙i长度，ltt为总喉道长度，lj&t为孔隙j圆心到喉道圆心的距离，li&t为孔隙i圆心到喉道圆心的距离。经过上述方法校正以后提取到的孔隙网络在一定程度能够较好地反映真实岩心的孔隙网络。

图3 孔隙与喉道划分校正示意图Fig.3 Schematic diagram for division and correction of pore and throat

在数字岩心分析过程中，所有的模拟都是在不同尺度的数字岩心或孔隙网络模型下完成的。获得结果后，如何将小尺度的岩心属性数据合理地加载到较大尺度，从而获得最准确的全尺度下岩心属性，这就是所谓的粗化技术。粗化技术是现代油藏模拟技术的重要研究内容,这也是当前国际上研究的一个难点与热点问题[19-20]。

本次实验选取的样品孔隙度6.76%，渗透率为0.02×10-3μm2，属于典型的致密储层，对该岩心进行切片扫描，完成样品的全岩心扫描到局部精细扫描，扫描精度实现微米级—纳米级扫描，总计获取不同尺度切片数量3 863张(图4)。

利用微米CT、纳米CT扫描结果，对此次分析的样品建立数字岩心孔隙网络，纳米级CT扫描结果得到的三维数字岩心模型在微细孔喉的刻画方面明显优于微米级CT，所提取出的三维孔隙空间尺度相对更小，数量相对更多。利用最大球算法、孔喉校正方法以及多尺度CT扫描融合技术，得到最终的孔隙网络模型，其中蓝色球体代表了孔隙，连接球体之间的白色管状代表喉道(图5、图6)。

图4 选取样品局部精细CT扫描示意图Fig.4 Partial fine CT scanning sketch map of sample

图5 微米CT扫描分割后三维数字岩石模型与提取出的三维孔隙空间Fig.5 3D digital rock model after segmenting by micron CT scanning and extracted 3D pore space

图6 微米/纳米尺度融合后的三维数字岩石模型、三维孔隙空间以及三维孔隙网络Fig.6 3D digital rock model,3D pore space and 3D pore network after Micro / Nano scale fusion

3 数字岩心两相渗流模拟

孔隙网络油水两相渗流模拟的基本过程如下：首先将孔隙网络模型充满水，此时网络模型由于饱和水而具有强亲水性；然后进行油驱水吸吮过程至束缚水饱和度，油驱过程会导致网络模型的润湿性发生变化；最后进行水驱油驱替过程模拟油田开发过程中的水驱采油过程。在油驱水和水驱油过程中可计算油水驱替相对渗透率曲线和吸吮相对渗透率曲线。

孔隙网络模型起初充满水，网络具有强亲水性。油驱水过程开始之前，孔隙网络模型入口端与油相相连。在进行油驱水的过程中，保持孔隙网络模型中的水相压力不变，逐步增大油相压力，油相以类似活塞推进方式进行驱替。油驱水过程结束后，油相作为非润湿相占据网络单元体的中央部位，水相作为润湿相占据角落位置。与非润湿相油相直接接触的单元体壁面，可通过改变油水接触角改变其润湿性。润湿性的改变使得水驱油过程中发生以下3种驱替方式:

(1)活塞式驱替，润湿滞后的存在使得前进接触角θa大于后退接触角θr，油水界面压力差降低导致曲率半径增大，在接触角达到θa之前，油水界面不能移动。

(2)孔隙体充填，驱替相水相通过喉道进入孔隙单元体时会发生孔隙体充填。对于配位数为Z的孔隙，如果该孔隙的相邻喉道中只有一个充满油，则其驱替过程与活塞式驱替类似。

(3)卡断，由于角落的水在膨胀之后与相邻其他角落的水相遇，直接导致整个单元体被水充满。通常存在2种情形:一是两油水界面同时向孔隙中央运动，在界面相碰时发生卡断，二是最小角处的油水界面移动到最大角处时发生卡断。

3.1 渗流模拟参数计算

在孔隙级流动模拟模型的基础上定义流体流动模拟规则，可模拟流体在孔隙网络的渗流特征。目前，拟静态流动模拟模型被广泛采用进行孔隙网络的渗流模拟。拟静态模型假设流动完全由毛细管力控制，模型中由黏滞力所造成的压降同毛管压力相比，可以忽略；流体是不可压缩的牛顿流体，且多相流体间不混相；根据侵入-逾渗理论，流体从一个孔隙流动到另一个孔隙是瞬时的，不考虑孔喉中的流动过程。

利用孔隙级网络模型可以模拟渗流的驱替和吸吮过程。将模型饱和一种流体，给模型施加一个驱动压力(pI-pO)，统计流体流量，则可由达西公式求解模型的绝对渗透率

(4)

式中：K为绝对渗透率，μm2；μi为i相流体的黏度，mPa·s；Qi为模型完全饱和i相流体时在模型所加压差下的流量，cm3/s；A为模型截面积，cm2。

油进入充满水的单元时所需要克服的毛细管入口压力，可通过Young-Laplace方程进行计算。当孔隙的形状和油水界面的接触角确定后，就可以计算出相应的毛细管压力。模型是用规则的几何形状来表征孔隙空间，所以具体到每一个孔隙孔喉中的油水分布可以定量地利用初等几何进行求解。计算出每一孔隙孔喉中的油水量后，就可以计算整个模型的含水饱和度：

(5)

式中：n为孔隙和孔喉的总数；Vi为孔隙或孔喉的体积，cm3；Viw为对应的孔隙孔喉中含水的体积，cm3。由于模型入口端面的孔隙通过孔喉与油藏相连，某一油藏压力下也就对应着网络模型的一个毛细管压力，由式(5)求出饱和度后就可得到整个模型的毛细管压力曲线。

对于模型中的每一个孔隙，由入口孔喉进入的流体量应该等于由出口孔喉流出的量，即由流量守恒可得

(6)

式中，Zi为与i孔隙相连的孔喉数，即配位数。设pi、pj分别为由喉道相连的两孔隙中的压力，则此两孔隙之间的流量

(7)

式中：Lij为两孔隙间的距离，cm；gij为两孔隙间总的传导率，它是两孔隙和孔隙间孔喉的传导系数的调和平均。如图7所示，传导率可通过式

(8)

计算。式中：gij为i、j两孔隙间的总传导率，cm4/(MPa·s)；Lij为孔隙i和孔隙j间的距离，cm；Li、Lj、Lt分别为孔隙i、孔隙j和喉道的长度，cm；gi、gj、gt分别为孔隙i、孔隙j和喉道的传导率，cm4/(MPa·s)。

图7 传导率计算示意图Fig.7 Schematic diagram for conductivity calculation

3.2 渗流模拟参数设置

利用提取出的孔隙网络模型，模拟油、水两相流体在孔隙网络结构中的渗流特征，渗流模拟参数设置的准确与否很大程度上决定最终模拟结果的精度，参考此次选取的样品所在区块的各项实际测量参数，设置本次两相渗流模拟参数，见表1。

表1 两相渗流模拟参数设置Tab.1 Parameters for Simulation of two phase seepage

3.3 两相模拟结果

数字岩心模拟结果表明：样品孔隙以晶间孔为主，总孔隙度为6.07 %，连通孔隙度为2.71 %，绝对渗透率为0.013×10-3μm2，略小于实际测量结果(孔隙度6.76%，渗透率为0.02×10-3μm2)。主要通过半径约70 nm的喉道连通，主要孔隙半径170～240 nm，模拟压汞结果与实测压汞结果无论在进汞曲线形态还是孔喉半径分布曲线分布范围上都具有较好的相似性。从进汞曲线可以看出，无论是数字岩心模拟结果还是真实岩心，曲线的平台段都对应在微细孔喉，平台段对应的孔喉半径为7.5～75 nm，在7.5～37.5 nm区间，模拟结果与实测结果接近重合(图8)。在孔喉半径分布区间方面，两者曲线结果形态存在一定差异，但总体形态趋势接近，在分布区间上的差异也并不明显(图9)。总体来说，数字岩心模拟出的孔隙结构特征与常规压汞实验测得的结果具有较好的一直性。渗流特征方面，模拟得到的两相渗流曲线，等渗点较低，同时两相共渗区较小，等渗点与端点饱和度之间的包络面积小，油、水两相渗透率干扰严重，与实际压汞实验所反映出的样品孔隙结构所对应的渗流特征较为符合， 具有典型的致密储层的渗流特征(图10、图11)。在中性偏油湿，驱替压力3 MPa的条件下，样品的等渗点饱和度约为48%，残余油饱和度约为27.5%。样品由于颗粒分散，导致边角增多，在一次驱替中容易存水，一定程度的降低了油相渗流通道的半径，因此在一次驱替完成时油相渗透率约为0.8。同样道理，水驱过程中，由于润湿性为中性偏油湿，所以边角处的油被锁住，导致虽然润湿性为中性偏油湿，但残余油饱和度却有27.5%(表2)。

图8 模拟压汞曲线与实测压汞曲线对比Fig.8 Comparison between simulation mercury injection curve and measured mercury injection curve

图9 模拟孔喉半径分布曲线与实测孔喉半径分布曲线对比Fig.9 Comparison between simulation pore throat radius distri-bution curve and measured pore throat radius distribution curve

图10 驱替阶段模拟相渗曲线Fig.10 Simulation relative permeability curves of displacement stage

图11 吮吸阶段模拟相渗曲线Fig.11 Simulation relative permeability curves of sucking stage

名称 值束缚水饱和度/%10残余油饱和度/%27.5残余油处的水相渗透率0.494电阻率指数(驱替阶段) 1.98Amott水湿指数0.44Amott油湿指数0.49Amott-Harvey润湿指数-0.05润湿性中性偏油

4 结 论

(1)数字岩心技术相对于其他常规手段，测量过程具有岩心无损的特点，在分析储层渗流特征的同时，可以提取出储层孔隙结构特征以及储层岩石物理参数，多尺度X-CT融合技术较好地解决扫描精度与三维数字岩心尺寸之间的矛盾，针对致密储层，能够保证较高精度微细孔喉刻画，使得提取出的孔隙网络模型更加接近真实岩心，从而保证渗流模拟结果的可靠性。

(2)通过多尺度融合技术所建立的数字岩心进行储层渗流特征模拟，模拟结果具有直观性的同时能够保证较高的模拟精度，相对于常规渗流实验，参数设置灵活，可以模拟不同驱替条件下两相分布特征，为储层渗流研究提供了一种新的分析方法。