采用机器学习分割算法和扫描电镜分析页岩微观孔隙结构

刘学锋, 张晓伟, 曾 鑫, 程道解, 尼 浩,李潮流, 俞 军, 胡法龙, 李长喜, 魏宝君

(1.中国石油大学(华东)理学院,山东青岛 266580; 2.四川省兴文县发展和改革局,四川宜宾 644400;3.中国石油集团测井有限公司, 陕西西安 710077; 4.中国石油勘探开发研究院,北京 100083)

页岩孔隙类型和结构参数影响储层流体的富集和运移,也影响页岩储层的渗透性和电阻率特性,准确评价页岩微观结构有助于提高储层评价精度和开发效率[1-2]。在传统扫描电子显微镜(简称为扫描电镜,SEM)技术基础上发展起来的大视域高分辨率SEM图像自动采集技术(MAPS)和聚焦粒子束扫描电镜(FIB-SEM)可在纳米尺度分别建立页岩的二维和三维灰度图像,为定量分析页岩有机质含量、孔隙度和孔隙结构提供了新途径[3-9],但分析结果的适用性尚未得到充分讨论。另外,灰度图像分割的精度决定了基于SEM图像定量分析页岩微观结构的准确性[10-11],但由于孔隙边缘粗糙和颗粒的充电效应造成图像局部灰度异常,基于图像灰度阈值的分割算法难以准确识别矿物组分和孔隙。目前,机器学习已成功应用于生物和医学灰度图像分割领域[12],该方法有望准确自动分割页岩SEM图像中的灰度异常区域。笔者以五峰-龙马溪组页岩为例,采用MAPS和FIB-SEM测试分别建立样品的二维和三维灰度图像,采用散射能谱扫描电镜(EDS-SEM,又称为QEMSCAN)建立样品端面二维彩色矿物嵌布图像,应用基于机器学习的算法分割灰度图像,识别孔隙类型,计算面积分数和孔径分布曲线。通过与气测孔隙度和低温氮气吸附等实验结果对比,评价采用SEM技术定量计算页岩孔隙结构参数的适用性。

1 样品来源和物性测试

四川盆地奥陶系五峰组和志留系龙马溪组沉积了一套海相页岩地层,龙马溪组自下而上可以划分为龙一段、龙二段和龙三段,其中龙一段又可划分为三个亚段,分别为龙一段一亚段、龙一段二亚段和龙一段三亚段。本文中选取3块龙一段亚一段页岩样品,样品编号分别为D1-3、J2-2和M1-2。

为评价采用SEM分析页岩微观结构的准确性,分别测量样品的孔隙度、总有机碳含量(TOC)、矿物组分和含量、低温氮气吸附脱附曲线,结果如表1所示。样品的氦气孔隙度采用Coretest公司生产的AP608型孔隙度仪测量,孔隙度整体偏低,均小于5%。采用PANalytical XRD设备测量页岩样品的矿物组分含量,样品骨架主要由石英、长石和方解石组成,含有少量的黄铁矿,黏土含量较高,主要黏土矿物为伊利石、绿泥石和伊/蒙间层。采用Elab-TOC总有机碳分析仪测量页岩样品的TOC,结果显示平均TOC为2.47%。

表1 样品物性参数测量结果

选用ASAP 2460比表面与孔隙度分析仪测量页岩样品的氮气(N2)吸附脱附曲线表征页岩的孔隙结构,结果如图1(a)所示,横坐标相对压力为无量纲参数。吸附-脱附曲线特征反映页岩孔隙形态,样品D1-3和J2-2页岩吸附与脱附曲线的滞回环较宽大,反映页岩中存在细颈广体的墨水瓶形孔隙。

图1 低温氮气吸附实验结果Fig.1 Experimental results of low-temperature nitrogen adsorption

样品M1-2 吸附曲线与脱附曲线几近平行,反映页岩中存在平行板状的狭缝型孔隙[13]。基于低温氮气吸附与脱附数据,采用Barrett-Joyner-Halenda(BJH)方程计算孔径和孔体积分布[14],采用Brunner-Emmett-Teller(BET)方程计算比表面积[15],结果如表1所示。样品J2-2的TOC和黏土含量最高,N2吸附结果表明比表面积和孔体积最大,微孔隙含量多。孔隙尺寸分布曲线如图1(b)所示。以孔隙体积随孔径的变化趋势(dv/dr)表征孔隙分布,可以反映孔隙空间在不同孔径范围内的集中程度,结果表明页岩样品发育大量半径小于10 nm的微孔隙,其中样品J2-2微孔隙含量最高,样品D1-3次之,样品M1-2微孔隙含量最少。

2 样品SEM扫描测试

为定量计算页岩孔隙结构特征参数,分别完成了3块页岩样品的SEM扫描测试,包括MAPS、QEMSCAN和FIB-SEM测试,测试参数如表2所示。MAPS测试获取样品端面的高分辨率大视域二维灰度图像,QEMSCAN测试建立样品相同端面的二维彩色矿物嵌布图像,两者结合用于划分页岩孔隙类型和统计孔径分布。FIB-SEM建立有机质发育部位的高精度三维灰度图像,用于分析局部孔隙连通性。

表2 页岩样品SEM测试参数

2.1 MAPS测试

选用蔡司(ZEISS)公司的Crossbeam 540仪器开展页岩样品MAPS测试。页岩样品端面采用氩离子抛光,将选定的扫描区域划分为一系列规整的网格,选择二次电子成像模式,自动获取每一网格的SEM图像,经拼接形成大视域高精度的MAPS图像,既具备了纳米孔隙识别能力,也满足了样品代表性的要求。图2为样品D1-3、J2-2和M1-2的MAPS扫描整体视域,灰度图像中不同灰度值(亮度)代表不同的矿物组分,高亮的区域表示矿物的原子序数高,亮度较低的区域代表原子序数小,但仅通过亮度差异无法判断不同亮度区域所代表的矿物类型。

图2 MAPS测试的整体视域Fig.2 Overall field of view for MAPS scanning

在整体视域下的MAPS图像基本无法识别孔隙,将MAPS图像放大,可以观察到页岩样品的显微特征,如图3所示。图中“OM”代表有机质,“OP”代表有机孔,“IP”代表矿物骨架孔隙,“Frac”代表微裂缝。图3(a)显示样品D1-3中有机质发育气泡状有机孔。由于孔隙边界平整度和充电效应的综合作用,孔隙边缘处亮度异常。图3(b)显示样品J2-2中发育蜂窝状有机孔,与气泡状有机孔相比,蜂窝状有机孔尺寸小。颗粒界面发育微裂缝,孔隙周边也存在亮度异常的现象。图3(c)显示样品M1-2发育草莓状黄铁矿,黄铁矿单晶间还存在部分有机质。

图3 MAPS测试表征的样品精细结构Fig.3 Fine structure of shales imaged by MAPS

2.2 QEMSCAN测试

为确定MAPS图像中不同灰度区域对应的矿物类型,选用QEMSCAN测试方法扫描页岩样品端面。首先获取扫描像素点的散射能谱(EDS),然后对照岩石矿物的元素组成识别该像素点代表的矿物类型,最后自动逐点扫描样品端面建立二维矿物彩色嵌布图像。选用蔡司(ZIESS)公司的MERLINSEM扫描页岩样品,扫描视域与MAPS测试视域相同。图4为D1-3、J2-2和M1-2样品的测试结果,3块样品骨架主要由石英、长石、方解石和白云石组成,黏土矿物以伊利石和绿泥石为主,与样品XRD矿物分析结果一致。统计QEMSCAN测试图像中各矿物组分的像素点个数可得面积分数Farea,i,即

(1)

图4 QEMSCAN获取的二维矿物嵌布图像Fig.4 2-D images of mineralogy maps measured by QEMSCAN

式中,设定有机质密度为1.3 g/cm3,石英的密度为2.62 g/cm3。设定各矿物的密度ρi,根据式(1)可计算质量分数Fmass,i,结果如表3所示。

表3 QEMSCAN测试获取的主要矿物组分含量

2.3 FIB-SEM测试

选用蔡司(ZEISS)Crossbeam 540开展FIB-SEM测试,选择二次电子成像模式,定量分析孔隙连通性。由于扫描视域小,选取有机质较为发育的位置开展FIB-SEM测试。经过图像预处理后的三维灰度图像如图5所示,扫描视域的物理尺寸约为18 μm×18 μm×7 μm。

图5 FIB-SEM扫描获取的页岩三维灰度图像Fig.5 3-D volume rendering of shale samples acquired by FIB-SEM

图5(a)和(c)显示样品D1-3和M1-2中发育草莓状黄铁矿,而J2-2样品黄铁矿含量少,视域中不包含黄铁矿,如图5(b)所示。FIB-SEM图像显示M1-2样品中有机质的尺寸最小。由于视域和选区的限制,无机矿物骨架颗粒尺寸较大,图像中包含的无机矿物颗粒少,FIB-SEM图像难以包括无机孔隙和裂缝信息。

在综合研究的基础上，根据古地理条件、古气候条件、古生物群特征、地壳活动等众多因素，对地层区域进行划分。一个时代的地层划分主要基于地层沉积特征的分布，类似于沉积划分。对于整个地质时期的地层，主要是基于每个地区地层发育的总面貌，称为综合划分。一般地层划分是指综合划分。

3 基于机器学习的SEM灰度图像分割

在页岩SEM测试阶段,分别采用MAPS测试和FIB-SEM测试建立页岩二维和三维灰度图像。若定量计算页岩孔隙结构参数,需要分割灰度图像,将相同组分的像素点标记为一种相态。因此灰度图像分割的精度决定了孔隙结构分析结果的准确性。在图3所示的MAPS图像中已观察到在孔隙边缘处会出现亮度突变,图6为选取的典型亮度突变区域,图6(b)为图6(a)中红色方框范围内放大后的图像。图像显示孔隙边缘处高亮区域与黄铁矿等高密度矿物灰度值重叠,常规分割算法难以准确区分。另外,灰度图像中任一像素点的灰度值不仅取决于该位置处的矿物组分,也受到该像素点周围像素点的影响。图6(a)中绿色方框中发育微裂缝或狭长的孔隙,由于裂缝边缘处灰度异常高亮,且裂缝宽度小,导致裂缝的灰度与有机质灰度重叠,常规分割方法也不能准确识别。

图6 SEM图像中灰度异常区域Fig.6 Anomalous contrast around grains in a SEM image

本文中采用基于机器学习算法分割灰度图像,分两步实现:①将MAPS图像和QEMSCAN图像配准,明确MAPS图像中不同灰度区域代表的矿物组分,为人工选取先验训练数据集奠定基础;②输入人工标记的训练数据集,采用基于机器学习的图像分割算法分割MAPS图像和FIB-SEM图像。

3.1 灰度图像与矿物嵌布图像配准

QEMSCAN根据散射能谱确定矿物组分,因此可准确识别大多数矿物组分;而MAPS灰度图像受对比度的限制,无法通过灰度区分所有矿物组分。综合考虑主要矿物类型和后续孔隙类型划分需要,本文中将页岩划分为孔隙、有机质、无机矿物骨架和黄铁矿4种组分。配准相同端面的MAPS图像和彩色矿物嵌布图,可确定上述4种组分所对应的灰度值范围。以样品J2-2为例,图2(b)为该样品的MAPS二维图像,图4(b)为该样品的相同视域二维彩色矿物嵌布图像。通过图像配准发现黄铁矿亮度最高,石英、长石和黏土矿物骨架亮度次之,有机质亮度较低,孔隙空间亮度最低。上述图像配准过程基于图像的灰度,所以不能选取灰度异常区域。

3.2 基于机器学习的分割算法

机器学习方法处理大型数据集,以导出统计规则和功能映射,从而自动执行特定任务[16],机器学习技术已成功应用于图像处理领域[17]。基于机器学习图像分割算法包括两个阶段:训练随机森林分类器和应用分类器分割SEM图像。训练随机森林分类器包含4个主要步骤:①手动选择训练像素并分配其矿物标签;②提取特征;③建立训练和测试数据集;④在训练数据集上训练分类器;⑤在测试数据集上测试所训练分类器的性能。在应用分类器分割SEM图像阶段,将训练好的分类器应用于页岩的SEM图像进行图像分割。

基于机器学习图像分割算法中,合理选取图像特征至关重要,决定了图像分割精度和效率。通常选取的图像特征越多,图像分割精度越高,但计算效率越低。综合考虑精度和效率,本文中提取6个类别,共计12个特征以构造数据集。

双边滤波考虑了空间信息和灰度相似度,以达到保留边缘和消除噪声的目的。选用空间半径5与色彩半径20、30和50结合起来以生成3个特征。

高斯滤波为非均匀的低通滤波,可以保持较低的空间频率,并减少图像噪声和图像中不必要的细节。通过将图像与高斯核卷积来实现,高斯核的二维形式表示为

(2)

式中,σ为分布的标准偏差;x和y为位置索引。σ值确定像素周围的模糊程度。随着σ值的增加,高频信息的内容会减少。根据最佳分割结果,将σ的值设置为1、2和3,生成3个特征。

最小值滤波是将当前像素的值替换为距该像素半径3以内的所有像素的最小值,生成1个特征。

均值迁移模糊是保留图像边缘的滤波算法之一,通常用于在分割图像之前消除噪声,可提高分割效果并保留图像边界。本文中坐标空间的半径设置为5,色彩空间的半径分别为30、50和60,生成3个特征。

上述11个特征,再叠加该像素点的灰度值(1个特征),共计生成12个特征,用于训练机器学习分类器。本文中使用随机森林模型来训练分类器,基于Scikit-Learn软件包实现。在随机森林分类器的训练过程中,树的数量为200,用于构建决策树的最大特征为3。

3.3 分割效果对比

由于SEM图像在孔隙边缘处存在灰度异常,因此在孔的边缘和骨架中的裂纹边缘选取灰度异常的像素作为训练数据,如图7所示。将所选择的训练像素分配组分标签,图7中不同颜色填充的区域代表不同的组分,提取每个训练像素的12个上述特征建立训练和测试数据集。

图7 划分典型的组分训练分类器Fig.7 Classifier trainer for typical components in shale samples

应用训练后分类器分割MAPS图像,结果如图8(a)所示。采用基于标记的分水岭算法分割相同的图像,结果如图8(b)所示。通过对比发现:① 基于机器学习算法识别的有机质边界清晰,如图8(a)黄色区域所示,较好地解决了孔隙边缘处的灰度异常问题,而分水岭算法分割结果中,有机质孔隙边缘处亮度异常区域被划分为黄铁矿,与实际不符,引起后续孔隙类型自动划分的误差,如图8(b)黄色区域所示;②基于机器学习的分割算法准确识别了微裂缝,如图8(a)中右下角黑色狭长区域所示,而分水岭算法识别的裂缝区域明显减少,在微裂缝或狭长孔隙识别方面存在误差。采用相同的流程,选取适当的训练数据集和特征参数,建立分类器,分割FIB-SEM图像,样品D1-3分割结果如图9所示。通过对比发现:基于机器学习算法分割结果中孔隙、有机质、无机矿物和黄铁矿的边界清晰光滑,而在分水岭算法分割结果中无机矿物骨架区域存在大量的孤立小孔隙,与实际情况不符。

图8 MAPS图像分割结果对比Fig.8 Segmented MAPS images comparison

图9 FIB-SEM图像分割结果对比Fig.9 Segmented FIB-SEM images comparison

4 孔隙结构分析

在准确分割灰度图像的基础上,可定量分析页岩储层的孔隙类型与尺寸分布。结合样品的区域地质特征,将页岩孔隙划分为有机孔、无机孔和裂缝,统计孔径分布。孔隙结构分析流程包括:采用簇标记算法识别孔隙;根据孔隙周围矿物组分识别孔隙类型;统计孔隙特征参数。

4.1 孔隙簇标记和半径计算

为统计孔隙特征参数,需要将相函数中相邻的孔隙像素点标记为同一孔隙。本文中采用Hoshen-Kopelman(HK)标记算法实现孔隙识别[18],将具有相同孔隙簇标记的所有像素视为一个孔隙,如图10所示,图中不同色彩表征不同的孔隙簇。之后,利用数字图像处理技术计算各个孔隙簇的最大内切圆半径R1、最小外接圆半径R2和等效圆半径R。

图10 HK标记算法孔隙簇标记结果Fig.10 Pore clusters labelled by HK algorithm

经过分割之后的FIB-SEM图像是三维数据体,包含了孔隙连通信息。采用三维HK算法标记孔隙体,若经过簇标记的三维数据体中两个相对的端面存在相同孔隙簇标记的像素点,则说明该孔隙簇贯穿整个样品,形成了连续的流体运移路径,否则FIB-SEM识别的孔隙不连通。该方法可分别判断X、Y、Z三个方向上的孔隙连通性。

4.2 孔隙类型识别

首先,根据孔隙簇周围矿物组分区分有机孔和无机孔。对于标记后的每一孔隙簇,通过数学形态学的开运算寻找孔隙周围的像素点,若孔隙周围像素点位于有机质,则划分为有机孔;若孔隙周围像素点属于无机矿物,则划分为无机孔。孔隙簇类型划分结果如图11所示:图11(a)显示有机孔,孔隙圆度高;图11(b)显示无机孔,孔隙形状扁平,主要由界面微裂缝和黄铁矿晶间孔构成。

图11 有机孔和无机孔划分结果Fig.11 Organic and inorganic pores distinguished results

然后,利用孔隙簇的纵横比区分无机孔中的裂缝和孔隙。孔隙的纵横比(α)定义为长轴和短轴之比,即

(3)

裂缝形状扁平,纵横比α大;孔隙形状接近圆形,纵横比α小。经过与扫描灰度图像的对比分析,设定孔隙簇纵横比阈值为5,将α>5的孔隙簇划分为裂缝。

5 结果讨论

5.1 图像分割结果

基于分割的SEM图像统计图像中孔隙、有机质和黄铁矿的面积分数,结果如表4所示。MAPS图像识别的孔隙度均低于氦气孔隙度,这是由于成像检测方法无法识别孔径小于分辨率的孔隙。3块样品中,J2-2样品的图像孔隙度和氦气孔隙度相差最大。根据低温氮气吸附结果,该样品中半径小于10 nm的孔隙含量最多,所以导致两种孔隙度差异最大。综合分析图像孔隙度、氦气孔隙度和低温氮气吸附结果发现:MAPS测试不能识别页岩所有的孔隙空间,若要分析页岩中所有尺度的孔隙,需结合其他测试方法。

表4 基于SEM图像计算的主要组分面积分数Table 4 Area content of main components calculated by SEM images %

与MAPS图像结果相比,FIB-SEM图像识别的孔隙度更低,这是由于FIB-SEM图像整体视域小,选取有机质发育的区域,仅反映了有机孔的局部特征,不足以识别样品中的无机孔隙。

MAPS图像统计的有机质面积分数如表4所示,QEMSCAN图像识别3块样品的有机质面积分数分别为0.91%、4.11%和0.36%(表3),QEMSCAN识别有机质含量远小于MAPS图像的计算结果。两种SEM测试扫描视域基本相同,但扫描分辨率不同,QEMSCAN图像的分辨率为98 nm,而3块页岩样品中均发育大量尺寸小于分辨率的有机质,无法被识别,如图3和图5所示。上述结果表明:受分辨率的限制,采用QEMSCAN图像识别的页岩有机质面积分数偏低。

3块样品的TOC分别为2.43%、2.79%和2.21%,MAPS图像计算的有机质面积分数与TOC的比值分别为2.84,3.93和2.50。J2-2样品的比值偏大,这是因为该样品中发育蜂窝状的有机孔,孔隙尺寸小,MAPS图像中未被识别的有机孔被划分为有机质,所以MAPS图像中统计的有机质含量偏高。

5.2 孔隙结构特征参数统计

基于MAPS图像计算了3块页岩样品的孔隙结构特征参数,如表5所示。样品D1-3无机孔隙占比高,样品中的有机质发育气泡状有机孔,但总体有机质含量不高,有机孔占比较低。样品中发育大量粒缘缝,因此样品裂缝孔隙度最高。样品J2-2有机质含量最高,有机孔占比也最高,但有机孔面孔率低,这是由于受图像分辨率的限制将未识别的蜂窝状有机孔划分为有机质。样品M1-2孔隙度最低,有机质中发育气泡状有机孔,但有机孔发育不均匀,导致有机质面孔率最低。

表5 基于MAPS图像的孔隙特征参数

统计MAPS图像识别所有孔隙的孔径分布曲线,如图12(a)所示。图中纵坐标为孔体积分数,即某一孔径的孔隙所占总孔隙体积的比例,孔体积分数采用孔隙度归一化,结果表明:在10 nm分辨率下,MAPS图像可准确识别孔径分布范围为10～200 nm的孔隙,3块样品中样品D1-3的孔隙尺寸最大,样品J2-2最小,样品M1-2居中,与MAPS图像直观观察到的现象一致。

图12 SEM图像识别孔径分布Fig.12 Pore size distributions calculated from images acquired by SEM

图12(b)为3块样品MAPS图像识别有机孔孔径分布曲线,平均半径如表5所示。结果表明:样品J2-2有机孔孔隙尺寸最小,样品D1-3最大,样品M1-2居中。这是由于样品中发育的有机孔类型不同,样品J2-2中发育海绵状有机孔,而样品D1-3和M1-2中发育气泡状有机孔,前者尺寸明显小于后者。

图12(c)为3块样品FIB-SEM图像识别所有孔隙的孔径分布曲线,孔隙半径与MAPS图像识别有机孔结果基本一致,说明FIB-SEM建立的图像主要包含了有机孔隙信息。通过对比图12(b)和(c)发现,样品J2-2和M1-2的FIB-SEM图像中存在孔径大于100 nm的孔隙,这是由于两块样品中存在部分圆度低的孔隙,基于三维图像和二维图像计算的孔径存在差异,如图13(b)和(c)中显示的形状不规则且尺寸较大孔隙。

图13 有机孔连通性分析结果Fig.13 Connectivity analysis results of organic pore

基于FIB-SEM图像分析孔隙连通性的结果如图13所示,图中不同的孔隙簇用不同色彩填充,填充颜色只为区别孔隙簇,与孔隙簇尺寸无关。连通性分析结果表明:3块页岩的有机孔连通性较差,FIB-SEM图像识别的有机孔均不能在X、Y、Z三个方向上贯穿样品,但形成了局部连通的孔隙簇。低温氮气吸附曲线证实页岩样品中发育大量孔径小于10 nm微孔隙,FIB-SEM图像识别的孤立孔隙簇可通过更小尺寸的喉道连通起来。

6 结 论

(1)与常规图像分割算法相比,基于机器学习的图像分割方法解决了孔隙边缘灰度异常导致的图像分割难题,为定量统计各组分含量和区分孔隙类型奠定了基础。

(2)由于分辨率的限制,MAPS图像识别的总孔隙小于氦气孔隙度,且微孔隙含量越高,图像识别孔隙度与氦气孔隙度的差异越大;MAPS图像识别的有机质含量偏高,而QEMSCAN图像识别的有机质含量偏低。

(3)基于SEM图像实现了页岩有机孔、无机孔和裂缝的精细划分,孔径分布结果表明气泡状有机孔半径明显大于海绵状有机孔,孔隙连通性结果表明有机孔连通性差,无法构成贯穿岩心的渗流通道。