一种页岩气储层有机孔提取及测井计算新模型

欧发辉，赖富强，夏小雪，李娇，马麒，李慧

（重庆科技学院复杂油田勘探开发重庆市重点实验室，重庆 401331）

0 引言

页岩气资源丰富， 在能源发展战略中占有重要地位。页岩气通常以游离态广泛存在于纳米孔隙中，四川盆地川南地区五峰组—龙马溪组的海相页岩储层发育大量有机质孔和少量无机质孔， 是页岩气的重要存储空间［1-2］。

测井作为测量地球物理参数的手段， 在页岩储层基本特征和解释评价方面发挥了重要作用， 可对矿物含量、 地化参数及物性参数等储层关键参数进行精细评价［1-3］。 近年来，国内外学者针对储层测井评价模型进行的研究取得一定进展。Luffel 等［4］通过研究岩心总孔隙度与测井参数关系， 采用概率统计方法建立总孔隙度测井模型，此方法考虑因素单一，结果偏差较大。侯颉等［5］通过建立矿物体积模型及测井响应方程，采用最优化算法反演计算出矿物组分含量和总孔隙度，然而此模型中有许多参数需要确定， 难以准确评估且操作复杂。 Herron 等［6］利用矿物元素质量百分比计算矿物骨架参数值，随后基于矿物体积模型计算孔隙度，最终建立优化的总孔隙度测井模型。 由于建立的测井模型有机质因素未考虑在内， 导致计算结果精度比较低，普适性不高，不适用于有机质孔发育的地区。因此，有必要建立一种新的页岩气储层有机质孔隙度的评价方法。

图像处理技术是当前各个行业的重点研究手段，传统的阈值分割算法和边缘检测算法已经在多种行业中得到应用。 如李应果等［7］通过阈值分割算法有效识别了单板穿孔缺陷现象，齐瑞燕等［8］利用迭代阈值分割算法处理电成像图像，效果良好。

为更精细地对不同层段有机质孔隙度进行计算，本文以川南地区五峰组—龙马溪组页岩（来自取心井X1，X2）为例，对研究区不同深度的扫描电镜（SEM）图像进行有机质孔隙提取及面孔率计算， 进而代入面孔率和孔隙度关系式，对有机质孔隙度进行计算，最终建立测井计算模型；同时进行核磁实验，利用实验结果计算的有机质孔隙度与基于SEM 图像有机质孔隙度的计算结果进行对比，验证该方法的准确性。 结果表明，该项工作有利于页岩气储层有机质孔纵横向展布规律研究，可完善页岩气储层甜点的精细测井评价。

1 区域地质背景

研究区主体位于滇黔北坳陷的中西部区域（川南低陡褶皱带南缘带）［6-10］（见图1）， 该区五峰组—龙马溪组地层具有厚度大、总有机碳丰度高的特点，且保存良好，页岩气赋存量巨大［9-10］。

图1 研究区构造及实验取心井位置Fig.1 Structure and location of experimental coring wells in the study area

本次研究的32 个样品均采自研究区典型取心井（X1，X2），取心井位置见图1。 为反映储层物性纵向变化，达到刻度测井解释目的，基本上等间距取样，取样纵深跨度基本为30～35 m，每口井取样8 个，基本覆盖了所有层位（五峰组—龙一14层），覆盖了所有优质储层，在底部最优质储层段，加密了取样间隔。

五峰组—龙马溪组， 岩性整体由硅质页岩逐渐过渡为灰质页岩到灰岩（见图2），表明五峰组—龙马溪组沉积期，水体逐渐变浅，盐度逐渐变高，氧化程度逐渐增强，总有机碳质量分数逐渐减小［11-17］。 研究区五峰组—龙马溪组一段局部地区发育深水砂泥质陆棚浊积岩砂体。龙马溪组二段沉积期，水体随深度减少而变咸且富含氧气，水动力条件也随之增强，总体上不利于有机质的形成与保存。 滇黔北地区中西部以浅水泥质陆棚沉积为主，五峰组—龙马溪组，岩性颜色逐渐变浅，灰质含量增加， 岩性整体上由硅质页岩逐渐过渡为灰质页岩到灰岩。龙一段主要发育深水泥质陆棚、灰泥质陆棚及深水硅泥质陆棚微相。 龙马溪组二段灰质含量增加，笔石化石含量相对减少，黄铁矿含量减少，主要发育为浅水泥质陆棚及浅水灰质陆棚微相。因此，微相控制了页岩的脆性矿物含量， 由龙一11亚段—龙一13亚段硅质页岩向上逐渐过渡为泥质页岩、灰质页岩。优质页岩气储层主要发育在五峰组—龙一11亚段， 该层段水平层理发育［11］。

2 有机质面孔率计算及核磁实验标定

有机质孔隙是四川盆地川南地区五峰组—龙马溪组海相页岩储层非常重要的储集空间［12］，所以对该地区的有机质孔隙进行识别尤为重要。 在SEM 图像中，有机质孔与无机质孔在亮度上有差异， 导致其成像后的阈值也呈现出明显差异。基于这种情况，本研究以川南地区X1，X2 井不同井段的扫描电镜图像阈值差异为出发点，利用Matlab 软件，构建改进Canny-Otsu 算法，对有机质孔进行提取并计算出其面孔率，同时搭建 核磁实验，计算有机质孔隙占比。

2.1 有机质孔隙的识别及相关参数计算

2.1.1 阈值分割法原理

阈值分割主要采用Otsu 阈值分割法［13-14］，基本思路为：根据岩石扫描电镜的实验结果，利用Matlab 分析软件编制相应的命令流， 通过计算得到扫描电镜图像不同区域的类间方差σ2。当图像分割阈值在灰度范围内顺序取值时，计算SEM 图像的σ2；当计算的σ2值最大时，此时的分割阈值即为最佳分割阈值。 σ2计算方法见式（1）。

式中：ω0，ω1分别为不同灰度值的概率之和；uT为扫描电镜图像灰度均值；u0，u1分别为扫描电镜图像不同区域灰度均值。

利用Otsu 阈值分割算法对有机质孔隙进行识别及参数表征的步骤为：1）对扫描电镜图像进行获取（见图3a）；2）对获取的扫描电镜图像，利用Otsu 阈值分割法对有机质孔与有机质进行阈值提取， 并利用灰度值绘制灰度直方图［15-16］（见图3b）；3）利用扫描电镜图像，通过特定的阈值，对有机质进行提取（见图3c）；4）利用扫描电镜图像， 通过特定的阈值对有机质孔隙进行提取（见图3d）；5）通过有机质孔隙的像素点和有机质的像素点数量的比值进行面孔率求取。

图3 有机质孔隙提取及参数获取Fig.3 Organic pore extraction and parameter acquisition

统计样品的图像数据， 可以进一步获得样品所处页岩段的微观孔隙结构特征， 但在有机质孔隙提取和面孔率计算的过程中，会有误差产生。产生误差的原因主要有： 扫描电镜图像中有机质孔隙与有机质、 无机质的阈值区分不明显， 使得阈值分割无法达到理想效果； 在有机质孔的提取中， 图像中无机质孔与有机质孔的阈值过于近似， 使用阈值分割算法无法剔除裂缝和无机质孔的影响， 对有机质面积和有机质孔面积计算造成误差； 在计算面孔率时， 由于有机质孔边缘存在高亮区域，导致有机质提取时存在空白区域，进而计算有机质的面积偏小。

当前软件还不能解决以上问题， 为了尽可能精确表征孔隙微细结构，需要将灰度变换增强［17］、边缘检测算法、阈值分割算法进行结合改进，才能使算法更适用于有机质孔面孔率的计算，有效解决存在的问题。

2.1.2 Canny 边缘检测算法

Canny 边缘检测算法是当前较为出色的多边缘检测算法，在多个领域已经得到广泛应用，基本思想是通过图像的幅度值与方向的极大值提取图片的边缘信息。 主要步骤为：1）将图像灰度化后进行高斯滤波，以此达到去除高频噪声的目的；2）利用算子计算图像梯度和强度，表征边缘强度与方向；3）根据梯度方向对图像梯度强度进行非极大值抑制，得到单像素边缘点；4）设置双阈值保留强边界和潜在边界，并连接边缘。

引入边缘检测算法的目的是为了有效避免无机质中的孔隙对于有机质孔隙提取工作的影响。 但是传统的Canny 算法对SEM 图像进行处理会存在有机质提取不准确、边缘不闭合的问题，需要根据实验需求进行拓展与改进。

2.1.3 有机质孔隙识别及参数表征

因SEM 扫描电镜图像中地质信息过多，存在孔隙边缘高亮问题和无机质孔隙的影响， 现对传统Canny算法进行相应的改进，并结合Otsu 阈值分割算法和灰度变换增强，形成一种改进Canny-Otsu 算法。 利用改进Canny-Otsu 算法，对SEM 扫描电镜图像（X2 井龙一12亚段样品（深度为2 248.88 m）有机质孔扫描电镜图）进行处理（见图4），具体步骤为7 步。

1）利用Otsu 阈值分割算法对SEM 图像进行处理，绘制灰度直方图， 获取图像有机质与有机质孔隙的阈值。 对原图观察分析看出，扫描电镜图像有机质、有机质孔隙区分不明显。对直方图进行分析得出，阈值基本为70～150，阈值集中，难以区分有机质孔隙与有机质。

2）对SEM 图像的阈值进行灰度变换增强处理。已经灰度化的SEM 图像按一定变换关系逐点改变源图像中每一个像素灰度值，使图像中的有机质、无机质和孔隙的阈值特征更加明显， 这有利于后期的有机质和有机质孔隙的提取。灰度变换增强处理得到的图像（见图4c）与原图（见图4a）对比发现，有机质与有机质孔隙变得更为突出。 由增强处理的灰度直方图（见图4d）可以看出，增强处理后的阈值更加趋近于两端，处理后增强的灰度直方图显示图像量化恰当， 使有机质与非有机质的阈值差异更大，在后续工作中，根据所求阈值范围自适应提取有机质面积更加精准。

3）采用自适应中值滤波代替高斯滤波平滑图像，去除噪点。中值滤波是一种非线性平滑技术，采用滤波窗口对待处理的图像矩阵进行处理， 将窗口中的图像像素点进行排序，并采用中值替代窗口中的像素值，以达到平滑图像的效果［18］。 中值滤波对比高斯滤波可以在兼顾保留像素信息的基础上不对原图像产生较大影响，但是中值滤波受滤波窗口大小的影响较大，易造成图像不连续。 自适应中值滤波是在中值滤波的基础上引入了自适应滤波算法的思想， 可以根据不同情况设好不同条件，动态改变中值滤波器的窗口尺寸。

4）对图像进行非极大值抑制。

5）利用Otsu 算法确定双阈值，进行阈值筛选。 利用第1）步提取到的有机质与有机质孔隙的阈值范围［0，L］确定有机质与无机质阈值分割值t：有机质与有机质孔隙的阈值分割值为t0，［t+1，L］为无机质阈值范围；［0，t］为有机质阈值范围［19］。 令有机质像素点的全幅比值为a0，平均灰度级为u0；无机质像素点比值为a1，平均灰度级为u1，图像总平均灰度为u，类间方差为σ。

根据Otsu 算法，当σ 取最大值时，此时的灰度级为最佳阈值，即为高阈值Th。 低阈值T1为

6）边缘膨胀使有机质的边缘连接成封闭区域。SEM图像中少数孔隙的阈值范围接近于有机质阈值范围，影响有机质提取（见图5）。 所以，需要进行边缘膨胀操作，经过膨胀后的有机质边缘会形成闭合区域［20］。

图5 无机质孔放大图像Fig.5 Enlarged image of inorganic pores

7）连通域面积过滤。 经过膨胀后的有机质边缘会形成闭合，将闭合区域进行面积筛选，保留大面积，去除小面积， 可以有效规避无机质孔隙与无机质闭合区域的影响。 改进Canny-Otsu 算法的流程见图6。

图6 改进Canny-Otsu 算法计算有机质提取流程Fig.6 Flow chart of organic matter extraction calculated by improved Canny-Otsu algorithm

利用改进Canny-Otsu 算法对实例图进行处理，可以获得有机质的区域， 通过计算区域面积可以获得有机质的面积， 然后通过对有机质区域进行二值化处理提取出有机质孔隙并计算出有机质孔隙面积。 传统Otsu 阈值分割算法处理后结果与改进Canny-Otsu 算法的处理后结果对比见图7、图8。

图7 有机质区域提取区域Fig.7 Extraction results of organic matter regions

图8 有机质孔区域提取区域Fig.8 Extraction results of organic pore regions

由图7a，7b 看出，有机质中去除了无机孔的影响，并且有机质孔边缘高亮区域的影响也得到了去除。 由图8a，8b 可以清晰看出，图7 中的无机质孔已经去除。

2.1.4 表征结果

在X1 井和X2 井的电镜图像中，选取代表性的8张（具有统一比例尺）图像进行有机质和有机质孔提取，计算有机质孔面积与有机质面积的比值，获取面孔率。 最终计算结果见表1。

表1 X1 井和X2 井面孔率计算Table 1 Surface porosity calculation of Well X1 and X2

2.2 基于核磁实验有机质孔和无机质孔比例计算

饱和水、 饱和油核磁对比实验的一个目标是判断有机质孔、无机质孔的比例［21］。 如图9 所示，蓝色数据点为压力在20.68 MPa、 温度在90 ℃状态下饱和盐水的核磁T2谱减去来样状态的核磁T2谱， 目的是消除黏土束缚水的影响［20-28］。 红色数据点为压力在20.68 MPa、温度在90 ℃状态下饱和油的核磁T2谱减去来样状态的核磁T2谱。

对比2 套T2谱的包络面积可以分别得到有机质孔（亲油孔）和无机质孔（亲水孔）的孔隙度。X1 井2 个样品有机质孔与无机质孔的比例分别为56.7∶43.3 和54.2∶45.8。 X2 井2 个样品有机质孔与无机质孔的比例分别为54.6∶45.4 和51.9∶48.1。

综合分析得出，研究区有机质孔、无机质孔的比例约为55∶45。

3 测井计算模型构建

利用研究区矿物元素类型及质量分数构建的储层矿物组分体积模型见图10。 利用最优化反演计算储层的有机质骨架干酪根、 无机质骨架矿物质量分数及孔隙度，为构建有机质孔孔隙度计算模型做准备。

图10 页岩气储层矿物组分体积模型Fig.10 Volume model of mineral components in shale gas reservoirs

3.1 页岩储层矿物组分的确定

X1 井和X2 井均进行了QEMSCAN 矿物成分分析，与LithoScanner 岩性扫描测井解释结果（岩心刻度后）对比结果见图11。 由图11 可以看出，整体上岩性扫描测井的矿物质量分数与岩心测试结果一致性较强，为计算模型的建立提供了基础。

图11 X2 井QEMSCAN 矿物成分结果与标定后岩性扫描测井处理结果对比Fig.11 Comparison of QEMSCAN mineral composition results with lithological scanning logging results after calibration of Well X2

3.2 总有机质骨架干酪根体积分数计算方法

1）总有机碳质量分数计算方法。 将密度测井响应特征纳入总有机碳质量分数计算影响因素，可表示为

式中：TOC为总有机碳质量分数，g/g；R为深侧向电阻率，Ω·m；驻t为声波时差，μs/m；ρ 为测井密度，g/cm3；A，B，C为无量纲拟合系数。

2）总有机质骨架干酪根体积计算模型［22］。 依据密度测井，假设ρ 为主要测量无机矿物骨架密度值，TOC也可为

式中：ρtoc为有机碳密度，g/cm3；Vtoc为有机质骨架干酪根体积分数。

根据式（5），可得：

井眼环境及含气性的一些影响导致测井值偏低，或黄铁矿的影响导致测井值偏高，因此，需要对其进行适当的校正，则校正后干酪根体积分数为［23］

式中：Vkero为校正后干酪根体积分数；λ 为校正因子。

式（7）中，对于Vtoc偏低情形时，λ 一般为1.0～2.0；对于Vtoc偏高情况下，λ 为0.5～1.0。 根据前人的研究成果，ρtoc一般为1.3 g/cm3，范围为［1.2，1.5］。

3.3 有机孔孔隙度计算方法

本文构建的模型参数， 主要包括泥质体积分数（Vsh）、石英长石砂岩体积分数（VQFM）、碳酸盐岩体积分数（Vcar）、黄铁矿体积分数（Vpyr）、校正后干酪根体积分数（Vkero），还有孔隙度（ϕ）。 通过最优化反演计算出这些矿物体积分数之后，则总孔隙度ϕt为

基于本次研究所涉及的基于阈值分割算法所求得的有机质孔面孔率和利用核磁实验所求得的有机质孔隙度占比，可分别得到不同的计算模型（见式（9）、式（10）、式（11））。

式中：β 为图像处理获得的不同小层的有机质孔面孔率均值；ε 为核磁实验确定的有机质孔占比；ϕorganic为有机质孔孔隙度；ϕinorganic为无机质孔孔隙度。

3.4 有机质孔隙度计算模型实例与分析

利用图像分割算法和基于最优化反演与核磁实验结合的有机质孔孔隙度计算式（9）、式（10），分别得到了X1 井（见图12）、X2 井（见图13）的计算有机质孔孔隙度和核磁实验有机质孔孔隙度。

图12 X1 井岩心实验测得的有机质孔孔隙度与测井计算模型结果对比Fig.12 Comparison of organic porosity measured by core experiment with results of logging calculation model of Well X1

图13 X2 井岩心实验测得的有机质孔孔隙度与测井计算模型结果对比Fig.13 Comparison of organic porosity measured by core experiment with results of logging calculation model of Well X2

2 条数据曲线对比可以看出，数据曲线基本吻合，只有局部有所偏差。 造成此现象的原因可能是改进Canny-Otsu 算法所处理的图像样本不足。 计算的实验有机质孔隙度和核磁实验分析的有机质孔隙度的绝对误差均在10%以内，X1 井样品的平均绝对误差为6.30%，X2 井样品的平均绝对误差为5.02%。

4 结论

1）通过对扫描电镜图像进行图像处理与分析，构建出适用于SEM 图像的改进Canny-Otsu 算法。 分别对X1 和X2 井的样品图像进行图像处理，计算出面孔率为11.68%～21.72%，平均面孔率为15.55%。 饱和盐水、 饱和油的核磁对比实验得出，X1 井和X2 井样品的有机质孔、无机质孔比例约55∶45。

2）依据最优化反演计算得到了储层的有机质骨架干酪根及总孔隙度， 结合图像分割算法得到的研究区面孔率与核磁实验得到的有机质孔和无机质孔比例，构建了基于研究区面孔率的有机质孔孔隙度计算模型和基于有机质孔占比的有机质孔孔隙度计算模型。

3）基于图像分割算法的有机质孔孔隙度计算方法和基于核磁实验的有机质孔孔隙度计算方法所得到的有机质孔孔隙度绝对误差均在10%以内。 这表明，本文的研究成果为页岩气储层有机质孔隙度测井计算提供了一种可借鉴的方法和思路， 为海相页岩气储层甜点评价和预测奠定了基础。