利用多重-大津阈值算法和扫描电镜分割CT图像

高衍武，吴伟，张虔，赵燕红，邵广辉，李国利，毛超杰

(中国石油集团测井有限公司测井应用研究院，陕西 西安 710000)(中石油玉门油田分公司勘探开发研究院，甘肃 酒泉 735000)(中国石油集团测井有限公司测井应用研究院，陕西 西安 710000)(中石油玉门油田分公司勘探开发研究院，甘肃 酒泉 735000)

页岩复杂的孔隙分布导致储层非均质性强，孔隙度、渗透率、孔隙大小分布是页岩孔隙结构评价中重要的参数[1～4]。作为典型的致密储层，有必要评价页岩孔隙结构特征，以便更好地理解页岩气体的储集、运移机理[5，6]。评价页岩孔隙结构方法主要有3种：高压压汞法、气体吸附法、场发射环境扫描电镜法(SEM)[7，8]。上述方法可以提供一些局部信息，如孔隙类型、孔隙大小和比表面积等，但都不能从三维尺度评价页岩孔隙结构，且均对样品有破坏作用。作为一种不破坏样品的测试方法，XCT(X-ray computed tomography)提供了一种高效的、定量评价页岩三维孔隙结构的方法，其主要由图像获取、图像处理、结构探测3个阶段组成。在图像获取阶段，分析图像特征来获取孔隙结构特征，其主要由图像分割和孔隙-裂缝提取2部分组成。图像分割是图像处理中最重要也最具有挑战意义的一步[9]，双扫描算法[10，11]和阈值算法[12～19]都可以用于孔隙结构分割。阈值通常作为图像分割体的标准来提取背景图像中的孔隙结构[20]。阈值技术又分为全阈值方法和局部阈值方法[21]，局部阈值方法对灰色信息局部化来选取多重阈值，该方法处理过程慢、复杂；全阈值方法对整个图像使用单一阈值，因此广泛用于CT(computed tomography)图像分割中。

多种阈值方法用来提取CT图像中的孔隙-裂缝[16,22,19]，如大津阈值算法[14，23]，但对多矿物组分CT图像方面的研究很少。利用数字岩心对多矿物组分CT图像进行研究，对于理解不同矿物组分对流体渗流的影响有重要作用。多矿物组分的CT图像和扫描电镜(SEM)图像也可以通过人工观测的方法分割[24]，但是人工观察的方法具有耗时、易受到人为主观因素影响等缺点。大津阈值算法[25]对图像直方图上的最大类间方差或最小类间方差选择阈值，是一种对图像选择阈值最常用的方法[26]。分割图像上的多重元素，则利用多重-大津阈值算法[27]。基于多重-大津阈值算法，笔者提出了一种多重阈值算法联合场发射环境SEM的方法，用以分割CT图像中的多重元素。

1 样品和方法

1.1 样品

研究样品主要取自2口井，S1样品取自一口井，S2、S3样品来自另一口井。研究区位于渤海湾盆地东南部的东营凹陷，页岩样品深度在2500～3000m，位于古近系沙河街组，是东营凹陷页岩油勘探开发的主要层位。3块样品都切削为直径25mm的圆柱状，利用氦气扩散法测量干岩样在实验室条件下的有效孔隙度，同时对样品碎屑进行了XRD(X-ray diffraction)矿物组分分析和总有机碳质量分数(w(TOC))测试分析(见表1)。该次研究中S1、S2样品利用多重-大津阈值算法和SEM进行CT图像分割，S3样品用于验证该方法。

表1 3块样品矿物组分质量分数、w(TOC)及有效孔隙度统计

1.2 试验

使用美国Xradia公司制造的Micro XCT-200仪器进行XCT测量，其对长度1mm物体的扫描精度可达0.7μm。圆柱状页岩样品垂直于样品台放置在扫描仪的中央位置。扫描样品的X射线源电压为150kV、功率10W。S1号样品长度5.2cm，直径2.512cm；S2号样品长度6.1cm，直径2.512cm；S3号样品长度4.8cm，直径2.512cm。样品的重构切片图像分辨率为2048×2048，切片长度为13.7135μm，产生的立体像素的体积为2578.96μm3。总体积为16.8275cm3的3个页岩样品在该试验条件下扫描产生1571个二维切片图像，并对二维切片图像叠加产生三维数据。

为了获取页岩样品矿物组分的灰度成分分布，对S1、S2样品进行SEM测试分析。在SEM观察之前，2种页岩样品岩心碎屑首先进行氩离子抛光，根据GB/T 5162—1997标准对岩心碎屑进行涂碳，目的是提高岩心表面的导电性进而增加图像的质量和分辨率。利用FEI Quanta 200F场发射环境扫描电镜进行观测，图像被放大8000倍，观测分辨率可达10nm，通过搭载能谱仪的X射线光谱仪定量分析页岩矿物组分。

1.3 方法

1.3.1 CT和背散射电子(BSE)成像技术原理

X射线照射到不同矿物时，由于光电效应、康普顿效应、电子对效应等复杂物理过程的影响，造成其信号衰减幅度不同[18,28]，可以用Lamber-Beer定律来表示：

I=I0e-μΔx

(1)

(2)

式中：I为传输强度，cd；I0为离子束发射强度,cd；μ为线性衰减幅度系数,1；Δx为X射线照射的样品厚度,cm；ρb为物体的骨架密度,g/cm3；a为Kling-Nishian系数,1；b取值9.8×10-24；Z是原子数目,1；E为X射线能量,kV。

在低X射线能量(<100kV)下，μ是一个与Z有主要关系的函数；在高X射线能量(>100kV)下，μ主要是与E有关的函数[29]。该次研究中X射线电压设置在150kV，因此X射线衰减幅度主要取决于物体的密度。不同的物质由于其密度不同产生不同衰减幅度的信号。CT值是用来定量评价X射线穿透物质时的衰减幅度，在高密度矿物中，X射线衰减量大，CT值也相应变大。CT值计算公式为：

(3)

式中：μw为水的衰减幅度系数,1。

该次研究的初始CT切片图像被转换成具有256位灰色范围值的灰度图像，在该范围中0代表最黑色，255代表最亮色。CT图像像素的灰度值主要与岩石样品密度有关：对于高密度物体，其图像为亮色；对于低密度物体，图像为暗色。

场发射环境SEM的原理是利用高能量电子轰击样品表面积产生次生电子和背散射电子，次生电子对样品表面敏感，可以有效揭示其微孔图像特征，因此被广泛用于孔隙研究。背散射电子是指从样品表面反射角度小于90°的原子束，其反射强度依赖于样品密度，物体的密度越高，产生的背散射电子反射强度越高，其生成的背散射电子图像就越亮。背散射电子图像像素的灰度值与岩石样品密度有关，因此可以用来印证CT图像不同矿物组分灰度值的分布。

1.3.2 多重-大津阈值算法

大津阈值算法是对图像选择阈值的一种最常用方法[25，26]，其核心思想是发现最大类间方差的阈值。一个图像可以描述为I(x，y)，其灰度值为0～L-1，L为灰度区间最大值，灰度值为I的像素值为ni，灰色图像中像素总数为n。该灰度值出现的概率Pi为：

(4)

整个图像的平均灰度范围μT为：

(5)

当一个单一的阈值t使用时，一个图像被分成2个部分——D0和D1，D0像素的灰度区间为[0，t]，D1像素的灰度区间为[t+1，L-1]。P0(t)、P1(t)表示D0和D1的累加概率，μo(t)和μ1(t)表示D0和D1的灰度平均值:

(6)

(7)

(8)

(9)

对于阈值t，其类间方差σB2(t)定义为：

σB2(t) =P0(t)(μ0(t)-μT)2+P1(t)(μ1(t)-μT)2

(10)

最优化阈值T为：

T= argmaxσB2(t)

(11)

当一个图像被分割成k+1部分时，k个不同的阈值(t1，t2，t3，…，tk)被使用，变量k的不同阈值的类间方差为：

(12)

多重-大津阈值算法中最优化阈值为：

(13)

1.3.3 CT图像和SEM图像处理

使用ImageJ数字图像处理分析软件对单一的BSE图像和CT图像进行分析，利用BSE图像中灰度值模式识别不同的元素相(包括孔隙-裂缝、有机质、矿物)，识别结果与XRD、EDX(搭载能谱仪的X射线光谱仪)测试结果进行相互验证。根据直方图法求取不同元素相的平均灰度值，利用多重-大津阈值算法求取CT切片图像中不同元素相的阈值；CT图像被分割之后，在ImageJ软件中分析所有元素相，分割的图像在ImageJ软件中导出后叠加用于重构三维数字岩心。

1.4 矿物和有机质体积分数

通过XRD分析测试得到的矿物组分含量是在没有有机质和流体的存在下矿物的质量分数，但是w(TOC)作为有机质在岩石中占据一定的比例，通过CT图像叠加分析计算的矿物组分和有机质含量是其占岩石总体积的体积分数，因此矿物组分和有机质含量应转换为体积分数，然后归一化处理。

为了消除有机质含量的影响，首先对矿物组分的质量分数进行归一化处理：

wnm,a=wm,a×(100-wom)

(14)

式中：wnm,a为矿物归一化后的质量分数,%；wm,a为XRD分析测试的矿物质量分数,%；wom为有机质质量分数,%;下标a为不同种类的矿物。

矿物组分和有机质含量根据其密度归一化后的体积分数为：

(15)

(16)

式中：Vnm,a、Vnom分别为归一化后矿物组分和有机质的体积分数,%；ρa为矿物组分的密度,g/cm3；ρom为有机质的密度,g/cm3；φ为孔隙度,%。

矿物组分、有机质、孔隙-裂缝归一化后的体积分数见表2。

表2 3块样品矿物组分、有机质及孔隙-裂缝归一化体积分数统计

2 结果和讨论

2.1 矿物成分鉴定及其灰度值

在BSE图像中，黄铁矿(最亮的多边形颗粒)、有机质(最暗的无定形颗粒)、孔隙-裂缝由于其自身特别高或低的密度被识别出来，其他的矿物组分,如石英、碳酸盐矿物、长石、菱铁矿、黏土矿物等由于其具有相同的元素组分，不能仅仅根据其自身的灰度值而被区分出来，图像中特别暗或亮的区域可能是由于部分矿物颗粒导电性较弱或表面不规则性而导致局部电子电荷释放造成的。

根据BSE图像、EDX图像和XRD分析测试结果确定其他矿物，每一种矿物和有机质根据其自身形状和灰度值来识别(见图1，其中(a)～(c)为S1样品；(d)～(f)为S2样品)。最亮的多边形颗粒为黄铁矿，在BSE图像中呈草莓状(见图1(c)～(e))；有机质和孔隙-裂缝在BSE图像中以无定形态呈黑色区域(见图1(b)～(f))；菱铁矿在BSE图像中呈最亮的平行四边形，由于其存在Fe元素，利用EDX可以将其区分出来(见图1(c)、(e)、(f))；碳酸盐矿物大多为方解石，在BSE图像中呈较亮的灰色不规则矿物颗粒(见图1(a)、(b)、(d)～(f))；石英呈黑色圆形颗粒(见图1(b)、(d)、(e))；正长石呈亮的灰色圆状颗粒(见图1(a)、(e))；斜长石在EDX图像中呈暗的灰色圆状颗粒(见图1(a)、(f))；黏土矿物通常为暗的灰色不规则长条状的集合体，其内部发育微构造(见图1(a)、(d))。

图1 S1、S2样品在BSE模式下的矿物组分、有机质和孔隙-裂缝及其灰度值

根据图1中每种矿物组分的平均灰度值，将7种矿物组分、有机质和孔隙-裂缝分成6种不同元素部分——有机质和孔隙-裂缝、石英、黏土矿物和斜长石、碳酸盐矿物和正长石、菱铁矿、黄铁矿。S1样品6种不同元素部分根据其平均灰度值依降序排列为：黄铁矿(平均灰度值254)>菱铁矿(平均灰度值237)>碳酸盐矿物和正长石(平均灰度值221～233)>石英(平均灰度值197～199)>黏土矿物和斜长石(平均灰度值191～193)>有机质和孔隙-裂缝(平均灰度值<78)。S2样品6种不同元素部分根据其平均灰度值依降序排列为：黄铁矿(平均灰度值为253～255)>菱铁矿(平均灰度值230～235)>碳酸盐矿物和正长石(平均灰度值190～199)>黏土矿物和斜长石(平均灰度值165～167)>石英(平均灰度值160～161)>有机质和孔隙-裂缝(平均灰度值<56)。

2.2 CT图像分割

在原始CT切片图像中央提取像素为1300×1300的局部图像以消除岩心外部区域造成的误差，为了区分6种不同元素部分，确定出5种不同的全阈值(T1、T2、T3、T4、T5)。根据大津阈值算法求取了5种全阈值并获取了分割图像(见图2)，为了获取5个全阈值，对每50个切片图像进行多阈值分割，每一个CT图像被分割为30个切片图像，最终通过30个切片图像中的平均阈值求取CT图像中每个元素的全阈值(见表3)。每个元素的阈值可以通过BSE图像中灰度值分布图来确定。如S1样品，黄铁矿、菱铁矿、碳酸盐矿物和正长石、石英、黏土矿物和斜长石、有机质和孔隙-裂缝的阈值分别为58～255、47～57、45～46、43～44、41～42、0～40。

根据上述阈值，利用ImageJ软件计算出6种不同元素部分所占图像的比例，对单一切片进行处理产生平面含量，最后对所有切片的平均平面含量进行计算，产生6种不同元素部分的体积分数。在原始图像中提取局部图像后，总共有1571张像素为1300×1300的灰色图像被用来作为样品。在灰色图像中，每个像素的灰度值范围为0～255，每个元素部分有一定的灰度区间值，利用Boolean方法计算每个元素部分的灰色值[19]：

(17)

式中：g(i,j)为图像的二进制值；f(i,j)为坐标(i,j)的像素灰度值；Ti-1、Ti为每一元素部分的阈值区间范围。

图2 CT切片图像切割过程

表3 CT图像中每个元素部分的全阈值

根据Boolean 方法，CT图像被转换为二进制的黑白图像，黑色图像(g(i,j)=0)代表需要研究的部分(region of interest，ROI)，通过计算黑色区域的像素值和图像中的像素总和进一步求取黑色区域平面含量：

(18)

式中：Ci为切片中每一元素部分的平面含量，%；NROI为黑色区域像素值，pt；NT为切片中所有像素总和,pt。

利用ImageJ离子分析工具分别计算出S1、S2样品的6个不同元素部分的体积分数(见图3(a)、(b))，与XRD、气测孔隙度等实际测量结果相比可以看出,CT计算值与CT实际测量值相关性较好，说明多重-大津阈值算法求取的阈值可以准确进行CT图像分割。

图3 页岩样品CT计算值与实际测量值关系图

3 S3样品验证

S3样品用于验证联合多重-大津阈值算法、多元素分割方法、BSE图像等方法对CT图像进行处理的结果。S2、S3样品取自同一口井、相同层位、有相同的物源和矿物组分，因此S3样品与S2样品的不同元素灰度值分布图类似。利用上述方法求取了S3样品的黄铁矿、菱铁矿、碳酸盐矿物和正长石、黏土矿物和斜长石、石英、有机质和孔隙-裂缝的阈值和体积分数(见表3，图3(c))。由图3(c)可以看出，S3样品的CT计算值与CT实际测量值有较好的正相关性(R2=0.90)，说明相同物源和矿物组分的样品在CT图像不同元素部分中有相同的灰度值分布，同时说明该分割方法可以用来准确、定量评价CT图像。

4 三维模型的重构

为了获取3块样品中6种不同元素部分在三维状态下的分布特征，对CT图像进行叠加重构三维模型。考虑到计算机运行处理速度，在叠加的CT图像中提取像素500×500×500的立体模型，根据6个不同元素部分的阈值把CT切片图像转换为分割图像，在ImageJ软件中处理分割图像产生不同元素部分的三维模型(见图4)。分割图像重构的三维模型可以直观地展现不同元素的分布特征和非均质性特征。由图4可以看出，S1样品有较小的非均质性，孔隙、有机质、黏土矿物在空间分布上有明显的非均质性弱的特点(见图4(a))；S2、S3样品有较强的非均质性，孔隙-裂缝、有机质、黏土矿物分布范围差异较大(见图4(b)、(c))。

图4 6种不同元素部分的三维重构模型

CT图像的多元素分割方法可以准确地确定6种不同元素部分的灰度值分布和阈值，该方法比人工观测更高效，能够消除人为因素的影响，但是也存在部分缺陷：不能有效区分CT图像中每一种矿物组分，特别是孔隙-裂缝与有机质，因为两者的灰度值分布相重叠[13]；每一个CT图像中，像素的灰度值代表了表面物体体积的衰减特征，当物体体积由多种不同物质组成时，其灰度值代表了平均衰减特征[19]，因此，部分孔隙体积效应对页岩矿物组分的鉴别带来困难。

5 结语

笔者根据BSE图像求取每个矿物组分灰度值分布范围，确定出每个矿物的阈值，然后根据多重-大津阈值算法分割CT图像。利用多组分分割方法叠加图像可有效切割CT图像，S1、S2样品计算的CT值与实际测量值相符。利用S3样品的CT图像验证该方法，结果表明：有相同物源和矿物组分的样品在BSE和CT图像中有相同的灰度值分布，利用该方法可以准确、高效地切割CT图像。