碳酸盐岩铸体薄片面孔自动提取研究

杜相仪 ，董少群，曾联波，何娟，孙福亭

1)油气资源与探测国家重点实验室，北京，102249； 2)中国石油大学(北京)地球科学学院，北京，102249；3)中国石油大学(北京)理学院，北京，102249；4)中国海洋石油国际有限公司，北京，100028

内容提要:岩石孔隙是储层储集油气的重要场所，通过铸体薄片进行岩石孔隙结构分析对储层质量评价具有重要意义。碳酸盐岩孔隙结构较为复杂，铸体薄片图像中存在大量噪声和干扰因素，致使常规方法面孔提取效果不佳，因此本文引入一种多阈值铸体薄片面孔自动提取方法(ctsPore方法)进行孔隙区域提取和面孔率估计，方法综合利用HSV色彩空间中不同参数提取孔隙区域。本文针对碳酸盐岩储层的特点确定了一套ctsPore阈值参数，以解决溶蚀孔与噪声易混淆的问题。实验以伊拉克A油田中新统 Asmari组A段碳酸盐岩储层为例对所使用算法进行了效果检验。算法精度和效率分析实验表明，基于ctsPore方法的碳酸盐岩铸体薄片面孔提取误差小于0.24%。全井段铸体薄片分析结果表明，A段储层面孔率范围主要分布在6%～8%，储层较为致密，从上到下A1～A3小层的储层平均面孔率分别为6.9%、11.7%、5.4%，与储层物性、沉积微相和岩性等纵向分布规律较为吻合。

孔隙结构特征是进行储层质量评价的重要参数，直接决定着储层的储集和渗流性能(赖锦等,2014;王翼君等,2020)。孔隙结构分析可以为储层质量评价、油田储量计算和开发方案设计提供基础。孔隙提取是孔隙结构分析的首要任务，其提取效果对于孔隙结构复杂且非均质性强的碳酸盐岩储层质量评价有着很大影响(谢淑云等,2019)，对于油田未来开发方案制定十分重要。岩石物性实验是研究储层物性特征的一种有效方法，但其成本较高，并且过程较为复杂；而铸体薄片数据相对成本低，可以获得更多的孔隙信息，由于铸体薄片分析目标特征单一、突出，易于快速准确地对孔隙区域进行识别和提取，所以在对孔隙区域进行提取时多数选用铸体薄片(赵明等,2009)。实际应用中，常采用目测法、计数统计法、图像处理法对偏光条件的铸体薄片图片进行分析，以提取孔隙区域、孔隙结构、面孔率等信息。

目测估计法(Folk,1951)需要地质人员根据视觉百分比估算图版来估计薄片中矿物和孔隙度的占比，该方法简单易操作，但主观性强，对解释人员的经验要求较高。计数统计法(Chayes,1949)利用点计数法对薄片中的矿物和孔隙占比进行分析，这种方法经济实惠且易于操作但同样因为主观原因使得识别精度受限。图像处理法根据铸体特定的颜色进行面孔提取，相较传统方法更高效和准确。图像处理法包括阈值法(张莹,2012;范宜仁等,2014;刘颜等,2017)和神经网络方法(Wardaya et al.,2013;蔡宇恒等,2020)等。神经网络方法根据人工标定的面孔样本训练面孔预测模型，以实现对面孔率的快速估计(Wardaya et al.,2013)，该类方法识别精度高，但是对人工标定样本数量和质量要求相对较高。阈值法是通过铸体薄片图片中像素的特征阈值，把图像像素划分为若干类的一种图像分割技术，较为常用的开源程序是基于ImageJ开发的jPor方法，该方法基于一个特定的256色调调色板和单阈值法区分孔隙和非孔隙区域，适用于被蓝色树脂染色的铸体薄片，其在冰岛的玄武质火山碎屑岩和纳米比亚的风成砂岩的面孔提取应用中均取得了较好效果(Grove and Jerram,2011)。此外，孔隙提取也可借助Photoshop等商业软件实现，PSQ方法是根据面孔区域颜色特性对岩石孔隙进行快速自动提取的一种方法(Zhang Xuefeng et al.,2014)，其可以进行岩石结构方面的定量化研究(郝哲敏等,2020)，该方法在二维面积百分含量分析上具有易操作、准确性高、统计分析方便等优点(张学丰等,2009)。

图像处理法可以在不同色彩空间中进行铸体薄片面孔提取，例如RGB空间(范宜仁等,2014;刘颜等,2017)、YUV空间(王子敏,2010)、HLS空间(Dai He et al.,2017)、HSV空间(张莹,2012;Dong Shaoqun et al.,2019)等，通过色彩空间中属性阈值的组合实现面孔信息的获取。其中HSV空间是一种较为常用的色彩空间，可以较好地反映人类视觉系统对色彩的理解，便于提取铸体薄片中被染色的树脂部分。因此本文使用HSV色彩空间进行面孔提取研究，具体采用一种基于HSV多阈值的铸体薄片孔隙自动识别方法(Pore extraction from cast thin sections by multiple thresholds，简称ctsPore方法)，该方法通过粗提取、错误识别带校正、降噪处理3个步骤，实现面孔区域的自动提取。原始ctsPore方法针对砂岩储层提出，已在准噶尔盆地南缘地区的砂岩铸体薄片中进行了实践应用，且效果良好(Dong Shaoqun et al.,2019)，但在碳酸盐岩铸体薄片孔隙提取中仍存在一些问题，例如一些溶蚀孔隙连续性差、孔隙边缘不分明等。因此，本文将以伊拉克A油田的碳酸盐岩储层为例，针对ctsPore方法在碳酸盐岩铸体薄片面孔提取的不足，对算法中阈值的选取进行调整，并确定了一套适用于碳酸盐岩储层的阈值参数。文章将通过算法精度分析、算法效率分析和储层特征分析等多组实验，检验该方法对于碳酸盐岩储层铸体薄片面孔率自动提取的效果。

1 研究区地质概况

伊拉克A油田位于伊拉克东南部与伊朗接壤处米桑油田北部(图 1)，构造处于扎格罗斯山前坳陷南部边缘的低角度褶皱带(欧成华等,2016;Wang Hongqiang et al.,2018,2019)。本文研究目的层为古近系渐新统—新近系中新统Asmari组A段，地层厚度约为70 m。受扎格罗斯造山运动产生的巨大北东—南西向水平挤压应力影响，油田表现为北西—南东向的长轴挤压背斜构造，南北共包含两个构造高点(韩耀祖等,2014;李开奇,2015;文磊等,2015;Ghanadian et al.,2017;Garzic et al.,2019)。Asmari组由上向下划分为A、B、C、D段，其中A段以白云岩、含硬石膏团块白云岩为主，硬石膏呈团块状出现，包含A1、A2、A3 三个小层 (图 1；Wang Hongqiang et al.,2019)。Asmari组为局限—半局限台地沉积，以碳酸盐岩沉积为主，陆源碎屑沉积较少，上部覆盖广泛分布的蒸发岩，A段的沉积相为局限台地，白云岩中蒸发岩(一般为硬石膏团块)丰度普遍呈上升趋势(Mossadegh et al.,2009;李开奇,2015;Ni June et al.,2017;Luo Sicong et al.,2019)。研究区储集空间复杂(邓虎成等,2014)，A段储层非均质性强，主要受原始沉积作用控制，为低孔、中—低渗储层。通过岩心物性测试得到A段平均孔隙度为8.7%，平均渗透率为10.61×10-3μm2，储层平均厚度为28 m，储地比为42%(Wang Hongqiang et al.,2019)。研究区原生孔隙较少，孔隙类型主要为成岩作用中溶蚀作用产生的次生孔隙，包含粒间溶孔、铸模孔、生物体腔孔等作为研究区碳酸盐岩的储集空间，且有微裂缝发育。

图1 伊拉克A油田位置及地层柱状图Fig.1 Location of A Oilfield in Iraq and stratigraphic column

铸体薄片来自研究区A-1井Asmari组A段储层(2969.56～3032.59 m)，共计122个薄片样本。岩心经过切割、染色、研磨等过程制成尺寸为25 mm × 25 mm × 0.5 mm的铸体薄片，在单偏光下使用多功能显微镜(Leixa DM4500 P)和附加的数码相机拍摄，铸体薄片图像分辨率为2560×1920。研究区薄片由蓝色树脂染色，孔隙区域呈蓝色，非孔隙区域为矿物和基质本身颜色。

2 铸体薄片孔隙自动识别方法

2.1 方法原理

本文所使用的铸体薄片孔隙自动识别方法ctsPore是由Dong Shaoqun等(2019)提出的一种修正的多阈值方法，结合了误差校正和图像细化。由于色调—饱和度—明度(Hue—saturation—value，HSV)色彩空间与人类视觉感知的颜色更为相似，所以方法选用HSV空间进行面孔提取(Dong Shaoqun et al.,2019)。色调(H)范围是0～360°(图 2a)，反映铸体薄片照片的颜色；饱和度(S)由圆的半径表示(图 2a)，范围是0～1，低值对应铸体薄片照片中颜色浅的部分，高值对应颜色饱和度高的部分；明度(V)代表颜色的亮度，范围 0～1，低值反映铸体薄片照片中偏暗的部分。该方法主要针对砂岩孔隙提取，而本研究区为碳酸盐岩储层，如果直接把原始ctsPore方法用于碳酸盐岩孔隙提取则会存在一系列问题，例如生物铸模孔连续性偏低、孔隙边缘区域识别效果欠佳等。因此本文将根据碳酸盐岩中孔隙的几何形态、矿物成分等影响因素，选定适合的HSV阈值范围，以准确地提取碳酸盐岩孔隙。

图2 ctsPore孔隙提取原理图(据Dong Shaoqun et al.,2019):(a)色调原理图;(b)阈值原理示意图Fig.2 Schematic diagram of pore extraction by ctsPore method (from Dong Shaoqun et al.,2019)：(a)Schematic diagram of hue;(b)schematic diagram of threshold

2.2 孔隙自动提取流程

基于ctsPore方法的铸体薄片孔隙自动提取主要分为3步：通过色调(H)进行孔隙粗提取、利用饱和度(S)和明度(V)的阈值将错误识别带校正、用滤波法进行降噪处理(图 3)，最终实现对孔隙区域的准确提取。本节以研究区的一张铸体薄片图片为例解释ctsPore孔隙提取的流程。

图3 孔隙自动提取流程图Fig.3 Flow chart of automatic pore extraction①—孔隙区域；②③—矿物内被错误识别成孔隙的区域；④—矿物与孔隙区域过渡带；⑤—矿物颗粒；⑥—孔隙区域内小噪点① —the pore area;②③ —the area that is wrongly identified as pore in the mineral;④—a transition zone between mineral and pore area;⑤—mineral particles;⑥—small noise in the pore area

(1)孔隙粗提取：铸体薄片中孔隙区域被染色树脂所充填，而这些区域具有相近的色调H范围，因此通过设定H阈值区间即可提取铸体薄片图像中孔隙区域(图 2a)。研究区碳酸盐岩原始铸体薄片单偏光下的照片(图4a)，来自A-1井Asmari组A段2984.08 m，薄片图像中深色为白云石，白色的为硬石膏，蓝色的为被染色的孔隙区域，多发育生物铸模孔，并有一条微裂缝发育。根据本研究区铸体薄片中孔隙被染色为蓝色的特点，并结合碳酸盐岩染色的特征，改进ctsPore方法中的阈值H，通过多次试验H选择180°～240°为宜。

对图4a的原始铸体薄片图片进行粗提取，得到如图4b的孔隙提取图片，图中蓝色部分为识别的孔隙，红色部分为孔隙以外部分。由图4b可见蓝色的孔隙区域被大致识别出来，但孔隙边缘十分粗糙，基质中噪点很多。此时识别的为图 3a中①～④的面孔，而矿物中的噪点(图 3a②③)、孔隙与岩石和基质之间的过渡带(图 3a④)被错误识别。

(2)错误识别带校正：在第一步粗提取的基础上，需要对图片移除低明度(V)或低饱和度(S)的面孔，即为去除图4b中错误识别的过渡带和噪点。根据研究区碳酸盐岩铸体薄片的特点，改进ctsPore方法中的参数，移除阈值dt=S×V≤0.12的部分得到图4c的孔隙提取图片，使得面孔率减少24.4%，去除了错误识别的过渡带以及非孔隙区域的噪点。图4c中孔隙区域的边界已经较为清晰，图4b中因过渡边界不清晰而错误识别的过渡带(图 3a④)被正确识别，但是仍存在很多不真实的小面积孔隙或岩石基质，在矿物颗粒表面或已识别的孔隙区域仍有一些不真实的小噪点(图 3b②③)。

(3)降噪处理：在去除错误识别带后，针对错误识别的小噪点问题，采用一种滤波法——洪水填充(Flood filling)算法，对噪点进行滤波。其中包含T1、T2两个阈值：T1为从基质中去除孤立的噪点；T2为从孔隙中去除孤立的噪点。滤波后得到如图4d的比较纯净的孔隙提取图片，面孔率为7.13%，相比图4c减少了0.5%。此时矿物颗粒内部和孔隙区域中的噪点被去除，孔隙区域和非孔隙区域边界识别清晰平整(图 3c)，孔隙区域被准确的识别出来。

3 伊拉克A油田碳酸盐岩铸体薄片孔隙自动识别结果

为了验证和对比薄片孔隙提取的效果，本文将ctsPore方法识别的结果和jPor方法、PSQ方法的面孔率结果从准确率以及效率方面进行比对分析。目前实际应用中普遍使用的是人工提取面孔率的方法，根据地质学家的经验进行提取。本节将从准确率和自动提取速度方面对比分析面孔率的提取效果，并且展示研究区铸体薄片面孔率提取结果。

3.1 面孔率提取准确率分析

通过对于同一张的研究区典型碳酸盐岩铸体薄片图像(图4a)利用不同的方法进行面孔提取，从而对识别效果进行对比分析。jPor孔隙提取结果如图5a所示，面孔率为8.73%，白色部分为非面孔区域，黑色部分为提取的面孔区域，从结果图可以看出大量噪点被识别为孔隙区域，因此计算的面孔率偏大。ctsPore方法孔隙提取结果如图5b所示，面孔率为7.13%，红色部分为非面孔区域，蓝色部分为提取的面孔区域，通过对比可以发现其结果对于溶蚀孔隙的识别更为准确，且噪点去除效果明显提升。

图5 不同自动提取方法效果对比：(a)jPor方法提取效果；(b)ctsPore方法提取效果Fig.5 Comparison of different automatic extraction methods：(a)The extraction effect of jPor method;(b)the extraction effect of ctsPore method)

选取实验区的37个不同深度点共70张铸体薄片图像进行孔隙提取实验，因为同一薄片不同视域下面孔率也会有一定差异，所以同一薄片取不同视域提取面孔率均值作为该薄片的面孔率。利用PSQ方法对研究区的铸体薄片进行了面孔率提取，得到如图 6a分布直方图，面孔率为0～15%，主要分布在4%～8%，平均数为5.4%；利用ctsPore法提取的面孔率如图 6b直方图，面孔率为0～14%，主要分布在3%～7%，平均数为5.14%。图 6中的交会图展示了试验区铸体薄片图片PSQ方法和ctsPore法的对比效果，散点靠近对角线表明解释结果较为接近真实情况。此处相关系数R2=0.95，平均绝对误差为0.24%，平均相对误差为5%，与实际情况较为吻合。综合而言，ctsPore法根据孔隙染色区域与背景区域的差异自动确定阈值，相对客观，可以避免解释人员主观因素造成的误差。

图6 面孔率提取效果对比Fig.6 Comparison of porosity extraction results

3.2 面孔率自动提取速度分析

选取研究区5张分辨率为2560×1920的铸体薄片照片作为测试数据，使用ctsPore方法重复对每个图像进行面孔提取实验100次，以分析面孔率自动提取速度。本节实验使用的电脑配置为Intel(R)Core(TM)i5-10400F CPU处理器、主频2.90 GHz、16.0 GB内存，实验结果如图 7所示，每张照片的提取速度在1.8～3.6 s，平均为2.5 s。总体而言，ctsPore方法可以快速高效地自动拾取铸体薄片面孔，为基于铸体薄片大样本数据分析储层质量评价提供了技术支撑。

图7 ctsPore 法提取速度直方图Fig.7 Histogram of running time used by ctsPore method

3.3 研究区铸体薄片面孔率提取结果

通过ctsPore法对研究区A-1井的目的层段Asmari组A段的铸体薄片进行面孔率的提取，可得到研究区A段总体以及各小层的面孔率分布。A段面孔率分布如图 6b直方图所示，面孔率平均值为5.14%，但其中包含了非储层以及砂质含量较高导致储层物性极差的铸体薄片，在仅保留储层物性相对较好的储层时，得到如图8的面孔率分布直方图，面孔率主要分布在6%～8%，和研究区通过岩心物性测试得到的平均孔隙度为8%比较吻合，由于测试方法的不同会有略微的差异。研究区的储层是较为致密的中—低孔储层，孔隙发育程度比较低。

图8 ctsPore面孔率提取结果直方图Fig.8 Histogram of porosity extracted by ctsPore method

研究区铸体薄片面孔率提取结果如图 9所示，A段储层的ctsPore平均面孔率为6.7%，为局限台地沉积环境，岩性主要为白云岩、含硬石膏团块的白云岩等，且白云岩中硬石膏团块丰度普遍呈上升趋势，上部覆盖蒸发岩。各小层的平均面孔率的分布情况如图 9所示，其中A1小层沉积微相为潮坪，储层的ctsPore平均面孔率为6.9%，受上部蒸发岩的影响，此时白云岩中多发育硬石膏团块(如图9中A1小层岩心图片所示)，部分孔隙被硬石膏充填(如图9中A1小层铸体薄片图片所示)；A2小层的储层平均面孔率为11.7%，此时随着深度增加，水深逐渐变深，硬石膏含量呈逐渐减少的趋势，发育大量未被充填的生物铸模孔(如图9中A2小层铸体薄片图片所示)；A3小层沉积微相为潮坪、蒸发坪、局限泻湖，储层平均面孔率为5.4%，硬石膏含量上升，孔隙度受影响降低(如图9中A3小层铸体薄片和岩心图片所示)。

图9 伊拉克A油田A段各小层面孔率分布直方图Fig.9 Porosity distribution histogram of each layers in the Section A in the A oilfield,Iraq

4 讨论

本文针对碳酸盐岩储层的特点确定了一套ctsPore阈值参数用于铸体薄片面孔提取，其中参数包含H(色调)、dt(S饱和度与V明度的乘积)、滤波法中的阈值T1和T2。以铸体薄片图片(图 10d)为例进行孔隙提取，首先H主要依靠孔隙区域颜色范围选取(图 10a)，研究区碳酸盐岩铸体薄片中孔隙被染色为蓝色，H可选择180°～240°范围，如果受到拍摄光线、白平衡等条件影响导致颜色有所偏差，可在此范围上进行一定调整；阈值dt的确定依靠S和V的共同影响，在铸体薄片图片中S和V值较低的部分(图 10b、c)，对应于颜色饱和度和明度较低的区域。去除dt小于0.12的区域，即可去除一些错误识别的过渡带，得到如图 10e的HSV空间中孔隙区域选择效果；在此基础上，需要通过滤波法去除其中的噪点，本文针对碳酸盐岩储层孔隙溶蚀孔多、噪点多的特点，选择去除像素点小于320 mm2的点，最终可得到研究区碳酸盐岩铸体薄片相对较好的提取效果。

图10 HSV空间中参数选取示意图：(a)H色调选择原理；(b)S饱和度选择原理；(c)V明度选择原理；(d)铸体薄片图片；(e)面孔区域选取效果图Fig.10 Schematic diagram of parameter selection in HSV space：(a)H selection principle;(b)S selection principle;(c)V selection principle;(d)Cast thin section;(e)Effect picture of pore selection

研究区原生孔隙较少，孔隙类型主要为成岩作用中溶蚀作用产生的次生孔隙，包含粒间溶孔、铸模孔、生物体腔孔等作为研究区碳酸盐岩的储集空间，且吼道较窄，多为硬石膏所充填，孔隙连通性较差。面孔区域的提取有利于后续判断孔隙和吼道类型(图 11)，并且为储层质量评价和储量计算提供重要参数。同时，面孔区域提取的为孔隙及吼道区域，通过对提取后孤立区域的长宽比等进行分析，可为后续研究孔隙连通性、确定孔喉配位数以及孔喉平均直径比奠定基础。

图11 孔喉区域及未充填缝示意图：(a)铸体薄片图像；(b)提取结果Fig.11 Schematic diagram of pore area and unfilled fractures：(a)Cast thin section;(b)final pore feature extraction)

研究区裂缝类型按充填程度可以划分为全充填、半充填、部分充填和未充填裂缝，且充填物多为硬石膏。其中，未充填的裂缝以及半充填和部分充填裂缝中的未被充填部分，会在铸体薄片中呈现蓝色(如图 11中未充填裂缝所示)。在进行面孔区域提取时这些未被充填的裂缝会被识别为孤立的小区域。在后续研究中，通过对孤立小区域的长宽进行比较，可以对裂缝的开度定量分析，为裂缝有效性研究以及裂缝建模(董少群等,2018)提供数据基础。

5 结论

通过对研究区的碳酸盐岩铸体薄片利用ctsPore方法进行孔隙自动提取，得到以下4条结论：

(1)本文利用ctsPore方法优选出一套适合碳酸盐岩储层进行孔隙自动提取的方法参数：H为180～240、dt大于等于0.12、T1和T2为320 mm2，该方法与PSQ方法提取结果平均绝对误差为0.24%，平均相对误差为5%，面孔提取结果较为可靠；

(2)该方法提取一张分辨率为2560×1920的图片均只需2.5 s左右，可以用于快速提取大量铸体薄片图片的面孔率，为海量铸体薄片图像孔隙特征分析提供可行的技术支撑；

(3)通过ctsPore方法对伊拉克A油田的铸体薄片进行孔隙自动提取，得到研究区面孔率主要分布在6%～8%，储层是较为致密的中—低孔储层，孔隙发育程度比较低；

(4)A1～A3各小层的面孔率提取结果分别为6.9%、11.7%、5.4%，与研究区碳酸盐岩储层物性、沉积微相和岩性等纵向分布规律情况较为符合。