基于CT扫描图像的碳酸盐岩油藏孔隙分类方法

2020-08-14廉培庆高文彬段太忠王付勇赵华伟李宜强

石油与天然气地质 2020年4期

廉培庆，高文彬，汤翔，段太忠，王付勇，赵华伟，李宜强

[1. 中国石化石油勘探开发研究院,北京 100083；2. 中国石油大学(北京) 石油工程学院,北京 102249； 3. 中国石油大学(北京) 非常规油气科学技术研究院,北京 102249]

全球碳酸盐岩油藏中的原油储量占原油总储量的50%，产量占原油总产量的60%，储、产分别达到半壁江山，已成为油田开发生产的主战场[1]。受岩相古地理、沉积构造演化及后期成岩作用等多种因素的共同影响，碳酸盐岩油藏具有复杂的储集空间和油气渗流特征，导致开发的风险和难度较大[2-5]。为了减少开发风险，需要定量描述油藏中孔隙、裂缝、孔洞等储集空间的大小、形状、连通性，从而准确判断碳酸盐岩储层的类型，正确选择碳酸盐岩的开发模式和开发方案[6-9]。

在碳酸盐岩储层表征研究中，常规方法主要基于岩心观察、薄片鉴定和压汞资料，着重突出沉积成因以及水动力学能量，通常以岩石矿物颗粒骨架为孔隙结构分类的主要依据，难以应用于油气运移规律的研究[10-12]。恒速压汞实验通过极低的恒定速度向岩心喉道及孔隙内注入液态汞，保证了进汞过程在准静态下进行，该方法能将喉道和孔隙分开，但主要依靠进汞饱和度判断孔径分布，不能直观地观察各种孔隙类型分布特征，并且该方法会破坏岩心[13-14]。近年来，随着核磁共振在岩心分析中的广泛应用，通过岩心核磁共振T2谱分布描述孔隙结构特征参数成为岩心物性分析的新方法。通过压汞曲线和T2谱分布相态对比，可以得到二者之间的转换系数，但该方法需要岩心饱和盐水，并且价格昂贵，难以用于批量岩心的快速扫描[15-16]。

CT层析扫描(Computed Tomography)是近些年来发展起来的一种利用X射线对岩石进行扫描成像技术，在不改变岩心外部形态和内部结构的前提下，可分析储层岩心的孔隙分布，进而了解岩石内部连通性[17-21]。目前，针对CT扫描图像多采用定性分析，没有进行定量的表征，无法在实际建模过程中得到应用。本文应用CT层析扫描技术对T油田碳酸盐岩岩心样品进行扫描成像，并基于阈值分割图像处理技术实现对岩心孔隙、裂缝和孔洞的定量分析，采用支持向量机(SVM)自动对岩心的孔、缝、洞进行识别，并提出了孔隙分类指数对岩心进行分类。该方法摆脱了现有技术中依靠研究人员经验判断碳酸盐岩孔隙类型的传统模式，为批量判定储层岩心孔隙类型提供了技术支持。

1 图像分割和特征参数提取方法

CT层析扫描技术广泛应用于油气藏的储层评价中，着重于确定岩石微观孔隙结构如孔隙大小、形状、分布情况和发育程度。本文主要使用CT扫描图像对圆柱形岩心的微观孔隙结构参数进行定量化提取。

1.1 CT图像获取

在T油田全直径岩心上钻取长度为50 mm、直径为38 mm的小岩心进行孔隙度和渗透率的测定，并对359块岩心进行了CT扫描。所用CT扫描装置为SOMATOM DR3，扫描方式可参见图1所示。对每一块岩心，在横向上等距离扫描4个横截面(横截面间距为6 mm，图1a)，在垂向上选取两个互相垂直的纵截面(图1b)。由于获得的CT图像局部存在噪声，会使图像中部分区域变的模糊，因此，本文采用自适应中值滤波的方法进行降噪处理，从而更加高效地对图像孔隙参数进行提取。自适应中值滤波的工作原理是在中值滤波基础上进行选择性优选，对滤波邻域的中心像素点进行判断：若被判断的像素点为邻域极值，则对中值滤波的统计排序进行降噪处理；否则跳过此被判断的像素点，进行下一像素点的邻域判断。

图1 岩心CT扫描方式Fig.1 The CT scanning procedures of core plugsa.横截面; b.纵截面

1.2 CT图像分割方法

利用CT机对岩心不同截面处扫描原理可知，CT图像中颜色越暗的区域其密度越小，孔隙越发育；因此，利用灰度图像不同像素的亮度值的阈值分割为孔隙判断提供了理论基础。如果阈值选取不合理，很容易将目标对象理解成为背景图像，或将背景理解为目标对象，这样在对目标图像进行特征分析时就无法有效获取目标特征。

阈值定义如下：假设F(i,j)表示图像像素灰度值，经分离后灰度值为1的像素表示目标对象，为0的表示图像背景，由下式获得目标对象：

(1)

式中：g(i,j)为目标对象的值;i,j分别为x,y方向像素点位置;M,N分别为图像在x,y方向的像素点数；T即为阈值，又称门限值。阈值的确定是图像分割的基础和关键。

对处理后的图像进行分析并分割，可达到区分图像中的孔隙区域与基质区域的效果。本文采用两种不同的图像分割算法来区分图像中的基质区域和孔隙区域[22-23]。

1) 基质区域为生长点的种子生长算法

使用种子生长算法主要针对CT图像中的基质区域进行生长识别，得出的二值图像进行反转便可得到图像中多数连通孔隙的分割。该算法以交互式选择的种子生长点的阈值差异为原则进行生长，预先通过人机交互的方式定义起始生长点，按照相似准则的要求进行生长。种子生长点设置在发育普遍的基质像素上，然后对灰度图像进行孔隙区域分离，获得的二值图像中的黑色区域(即为0区域)为孔隙区域，白色区域则为岩石骨架基质区域。图2显示了通过种子生长算法分离出的孔隙区域与基质区域，相似生长原则的阈值设置为40。图2a1—a3是基于裂缝型CT扫描图像分离出的孔隙区域与基质区域，可识别出具有较为明显的裂缝。图2b1—b3是基于孔洞型CT扫描图像分离出的孔隙区域与基质区域，也可识别出明显的溶蚀性孔隙区域。

图2 通过种子生长算法分离出的孔隙区域Fig.2 The pore area separated by seed growing algorithma1—a3.裂缝型；b1—b3.孔洞型

2) 图像边缘检测算法

图像中边缘点的特征主要表现为图像中周围像素灰度有阶跃变化或屋顶变化的那些像素点，即灰度值导数值最大或者极大值的位置处，图像属性中的显著变化通常反映了属性的重要意义和特征，利用图像的边缘检测算法对图像进行分割处理，对关键区域的边

缘进行识别，进而了解岩心中的重要特征。孔隙区域边缘提取算法可以用来验证基质区域为生长点的种子生长算法的分割效果。本文针对低精度CT图像和局部含噪的CT图像的边缘检测使用Canny算子，不但对图像噪声有较好的抗干扰性，而且对图像不同方向上的边缘检测具有较高的精度。

图3显示了通过孔隙区域边缘提取算法分离出的孔隙区域与基质区域。其中，图3a是基于裂缝型CT扫描图像分离出的孔隙区域与基质区域；图3b是基于孔洞型CT扫描图像分离出的孔隙区域与基质区域。可发现无论裂缝还是孔洞，Canny算子均表现出良好的效果，原图中的主要裂缝和明显孔洞的边缘均得到准确的识别。边缘检测算法是对图像分割的再次验证，也是后续图像孔隙参数提取的重要依据。

图3 通过Canny算子分离出的边缘效果Fig.3 The edges separated by canny operatora1—a3.裂缝型；b1—b3.孔洞型

1.3 孔隙特征参数提取

基于分离出的孔隙区域图像提取孔隙特征参数，并识别孔隙区域中的裂缝和孔洞，可以为孔隙类型分布、发育程度等信息建立无因次参数团进而定量表征提供基础。孔隙参数提取主要包括形态学处理和图像特征计算2个步骤。

1) 形态学处理

图像分割后产生的孔隙形状通常是不规则的，使用形态学处理不但对分割后的孔隙形状进行规范，而且可以去除图像分割中多余的噪声像素，使孔隙形状更加真实平滑。对分割出的孔隙区域进行形态学开闭运算处理，选择性地间断或连接孔隙，以对分割出的孔隙区域完成形态学上的校正和平滑处理，以及去除无效噪声点。对二值化的CT图像进行分割后，本文在识别出连通域的基础上，同时使用开闭运算进行形态学处理，对微观孔隙进行平滑处理。形态学处理并未从根本上改变岩心孔隙参数，仅对多余的像素点选择性的删除，从而保留了图像中的重要特征(图4b)。

2) 图像特征参数提取

主要是从CT扫描图像中提取特征参数。通过该步骤可以获得孔隙区域的面积和长度两大类型的参数，包括面积、凸面积、比表面、有效长度、迂曲长度、等效宽度、最大内切圆半径、等效直径、离心度和方位等10余项参数。

以最大内切圆半径和方向为例说明特征参数的提取过程。最大内切圆半径是孔隙评价中的重要参数，通过欧氏距离图获得孔隙区域中不同位置处的内切圆半径数值，取孔隙区域中的最大值作为最大内切圆半径，本文通过对二值图像的背景色区域进行欧式距离图的转化，选取其中的最大值为最大内切圆半径，如图4c所示。

图4 形态学图像、欧式距离和方位角玫瑰花图Fig.4 Rose diagrams showing morphological image,Euclidean distance,and azimuth a.原图；b.形态学图像；c.欧氏距离；d.方位角玫瑰花图

方向参数在岩心CT图像中起着辅助性的作用，方向的识别时专门针对于具有明显特征性岩心而言的，如裂缝型岩心，方向的识别对区域岩心的局部应力场分析具有指示意义。方位角的识别依赖于图像连通域中可匹配所有像素的转矩大小，由所有像素点在x-y坐标系中的净距可计算出连通域中的方位角，具体计算如下：

(2)

Sx=∑xi；Sy=∑yi

式中：xi,yi为连通域所有像素点坐标位置，mm；A为连通域面积，mm2。

基于识别出的方位角做方位角玫瑰花图(图4d)，可有效看出图像中的主要应力方位，若要弄清的该岩心截断面的主应力方位仍需要清楚知道该岩心在地下产状，结合产状可分辨出区域的应力椭圆分布情况。

2 裂缝、孔洞自动识别和分类

在前期图像特征参数提取的基础上，选取已经过图像分割、形态学处理后的CT图像中典型裂缝和孔洞参数作为特征向量，将选择的特征向量进行支持向量机训练，根据支持向量机训练模板类对新的孔隙参数进行训练并识别孔隙类型。

2.1 特征向量建立

选取经过图像分割、形态学处理后的CT扫描图像中具有典型特征的长条状裂缝和似圆状孔洞参数为特征向量。表1是24块岩心经过区域特征参数提取后对孔隙类型分类的结果，其中类型Ⅰ为裂缝型，类型Ⅱ为基质型，类型Ⅲ为孔洞型。以相同的方法对其余编号岩心进行区域特征参数提取和孔隙类型初步标定等处理。

表1 不同连通域岩心孔隙分类及特征参数Table 1 Pore classification and characteristic parameters of cores from different connected domains

1) 长宽比和纵横比

长宽比和纵横比均为评价孔隙空间条状形态的发育程度的重要参数，区别在于计算中使用的参数不同，长宽比采用的计算参数为迂曲长度，其值为迂曲长度与等效宽度之比。与之相对应的为纵横比，计算参数则变为有效长度，其值为有效长度与等效宽度之比。

2) 迂曲度

迂曲长度主要评价孔隙空间的中轴骨架的弯折程度，其涉及参数主要为迂曲长度和有效长度两大参数，其值为迂曲长度与有效长度之比。迂曲度的数值变化均在1以上，数值越大反应弯折程度越大。

3) 离心度

离心度是评价孔隙空间形状特征的重要参数，计算方法为：将孔隙空间近似成椭圆状，分别计算得到近似椭圆长轴a和短轴b的长度大小，计算公式为：

(3)

式中：ω为离心度，无因次，数值范围介于0～1，数值越接近于1，孔隙空间形状越扁；a为椭圆长轴长度，mm；b为椭圆短轴长度，mm。

4) 形状因子

形状因子主要用于描述孔隙边缘的平滑和规则性，其数值大小用来表征孔隙形状均匀规则的程度。由于形状因子可以对孔隙的形态作定量化的表征，因此它是后续判断裂缝、孔洞类型的重要参数。形状因子的计算方法如下：

(4)

式中：σ为形状因子，无因次；S为孔隙面积，mm2；C为孔隙周长，mm。

5) 实心度

依据区域面积和凸面积的大小可以进一步确定区域孔隙实心度。实心度值大小介于0～1，能够反应形状规则程度、内部孔洞发育情况，为后续的裂缝、孔洞识别提供重要参数。区域孔隙实心度的具体计算方法如下：

(5)

式中：λ为实心度，无因次；SA为孔隙凸面积，mm2。

通过以上无因次参数团的计算后，可生成由无因次参数构成的标准特征向量，为后续的支持向量机的样本学习提供基础，无因次标准特征向量具体可见表2。

表2 孔隙结构参数无因次化后的特征向量Table 2 Dimensionless eigenvector of pore structure parameters

2.2 支持向量机分类识别

支持向量机通过构造多个分界面处理有监督的数据分类问题，是一种基于统计学习的机器学习方法，具有优秀的样本学习能力，并较好地解决了非线性、高维度、局部极小值等问题。本文采用支持向量机对孔隙类型进行自动化的分类筛选，从而为后续孔隙类型的判断提供数据支撑[24-25]。优选采用两参数组合的方式对裂缝、孔洞类型分别作训练，依次评价分类效果。

分别优选出基于形状因子与离心度、实心度与离心度组合的分类方法。图5示例性地显示了在两种无因次参数团组合下的SVM训练的分类测试结果。图5a为形状因子与离心度的参数组合下SVM针对孔洞型孔隙和裂缝型孔隙的分类训练效果图；图5b为实心度与离心度的参数组合下SVM对孔洞型孔隙和裂缝型孔隙的分类训练效果图。其中，“+”代表孔洞型孔隙的标准特征向量，“*”代表裂缝型孔隙的标准特征向量，在形状因子与离心度或实心度与离心度这两种参数组合下具有良好的分类效果。此外，“☆”代表着待测试的孔洞型孔隙，“◇”代表着待测试的裂缝型孔隙，通过SVM分类训练可发现，“☆”和“◇”均落在了与之对应的孔隙类型的范围内，其反映出通过形状因子与离心度或实心度与离心度的参数组合下的SVM分类训练是可靠的。

图5 两种无因次参数团组合下的支持向量机分类测试结果Fig.5 SVM classification results based on two kinds of dimensionless parameter combination a.形状因子与离心度组合下分类测试；b.实心度与离心度组合下分类测试

3 岩心孔隙参数定量表征及应用

通过对岩心不同位置截面处连通孔隙进行识别和类型判断，进而对岩心整体的孔隙类型进行再判断。通过定义孔隙分类指数将单一孔隙类型的判断上升到岩心孔隙类型的判别。

3.1 孔隙分类指数

孔隙分类指数(PCI)是用来评价孔隙区域中裂缝型孔隙与孔洞型孔隙相对大小的参数，通过PCI值可将岩心孔隙类型划分为裂缝型、孔洞型、基质型3大类[26]。1)若PCI<1，碳酸盐岩孔隙类型为基质型；2)若1≤PCI<10，碳酸盐岩孔隙类型为孔洞型；3)若PCI≥10，碳酸盐岩孔隙类型为裂缝型。PCI计算公式可以表示为：

(6)

式中：Ф为岩心孔隙度，%；FL为横截面中最大裂缝长度，mm；CL为岩心长度，mm；CD为岩心直径，mm；EFL为纵截面中有效长度大于指定长度的裂缝长度总和，mm；FP为裂缝面积占孔隙区域百分数，无因次；HP为孔洞面积占孔隙区域百分数，无因次。

PCI将纵截面中的裂缝EFL和横截面中的裂缝FL作同等权重加和处理置于分子位置上，以此表示裂缝信息全面性的程度，若裂缝信息越丰富，则PCI值越大。分母上反映的是岩心储存空间中孔洞的相对大小，孔洞越发育，PCI值越小。孔隙度是同一岩心中不同位置处CT图像分割的重要依据，以孔隙度为基础的图像分割算法所获得横、纵截面裂缝信息和裂缝与孔洞之间相对大小等参数均是孔隙度的函数，因此在分母上引入孔隙度参数。另外，由于岩心中的裂缝是高渗区域，对渗透率的贡献明显，但对岩心的储集能力贡献极小，以孔隙度作为分母可有效拉大裂缝型岩心和孔洞型岩心之间的距离。基质型岩心的孔隙度和横、纵截面中裂缝丰度均较小，会导致具有较大孔隙度的孔洞型岩心与较小孔隙度的基质型岩心PCI值大小接近。为了解决这一问题，在分子上引入平方项，可使截面中的裂缝信息较小的发育程度与较低的孔隙度抵消，保证孔洞型与基质型岩心之间的数值距离。图6显示了不同大小的孔隙分类指数所对应的不同孔隙类型的岩心的图像，可明显区分出基质型、孔洞型和裂缝型岩心。

图6 基于孔隙分类指数的不同孔隙类型岩心CT扫描图像Fig.6 Core images of different pore types corresponding to different PCI values

3.2 分类方法在T油田的应用

孔隙大小是影响渗透率的主要因素，但碳酸盐岩储层往往结构复杂，非均质性强，孔隙度-渗透率之间的关系往往并不明显。在实际工作中，人们常常希望用岩心数据建立起孔隙度-渗透率换算关系来估算渗透率，这就要求我们在建立碳酸盐岩储层孔隙度-渗透率关系时对不同类型的岩心数据进行分类。

T油田主力开发层为白垩系的M组和F组灰岩。早期研究没有对微裂缝和孔洞的发育没有进行详细的分析。在孔隙结构研究方面，主要基于常规孔、渗数据，以及毛管压力资料，开展了储层岩石物理相划分。本文在基于CT图像的碳酸盐岩岩心参数提取的基础上，根据PCI评价指标对T油田M油藏和F油藏359块岩心进行了分类。根据孔渗数据测量结果，绘制了不同类型岩心孔渗数据交会图(图7)。M油藏191块岩心中基质型占比67.5%，裂缝型占比23.6%，孔洞型占比8.9%，整体上M油藏属于孔隙型油藏。F油藏168块岩心中基质型占比50.6%，裂缝型占比25.1%，孔洞型占比24.3%。F油藏岩心孔洞发育比M油藏多，但整体上仍为孔隙型油藏。