孔隙结构图像分析方法及其在岩石图像中的应用

张吉群, 胡长军, 和冬梅, 常军华, 李心浩, 李华

(1.北京科技大学, 北京 100083; 2.中国石油规划总院, 北京 100083; 3.中国石油勘探开发研究院, 北京 100083)

0 引 言

国内外研究孔隙结构的方法主要分为3类。①室内实验法。该方法主要包括毛细管力曲线法、铸体薄片法、扫描电镜法、CT扫描法,这4种方法如果不借助第3类方法,一般都是采用人工鉴定,分析速度慢,准确度取决于个人鉴定水平。②测井资料现场评价法。该方法主要是用电阻率测井资料或核磁共振测井研究岩石孔隙结构,但该类方法主要是反映大面积、纵向上岩石宏观孔隙性质。③图像分析法。该类方法主要针对第1类方法中的铸体薄片法、扫描电镜法、CT扫描法所生成的图像,应用数学形态学等学科分析孔隙的微观结构。

图像分析法对算法的准确性要求较高,特别是对识别喉道算法要求极高。因为喉道识别的准确性直接影响每个孔隙的属性数据(包括面积、周长、配位数等),以至影响岩石孔隙结构参数。识别喉道的数学形态学的方法主要有流域分割方法[1-2]、中轴线算法[3]和采用数学形态学的灰度膨胀和腐蚀方法[4]等3种。这3种识别喉道算法的共同特点是先由极限腐蚀算法[1]确定孔隙位置,进而寻找连接孔隙间喉道的位置。每种方法都有其优势和局限性。前2种方法对图像的噪声较为敏感、抗干扰性差,因而过分割严重,第3种方法的误分割较多。

本文提出了一种改进的识别喉道算法,并结合计算机图形学算法分析计算孔隙的各项属性数据,提高识别喉道和分析孔隙各项属性的准确性和速度,并以岩石背散射扫描电镜图像为例,说明其他图像如铸体薄片图像等也可做类似处理。

1 算法流程与数学描述

对岩石孔隙结构分析的算法流程如图1所示。

图1 岩石孔隙结构分析算法流程图

1.1 去除岩石孔隙图像的噪声

岩石背散射扫描电镜图像受离散脉冲和椒盐噪声的影响较为严重。去除脉冲干扰及椒盐噪声最常用的算法是图像处理学中的中值滤波[5]。中值滤波是一种去除噪声的非线性处理方法,在某些条件下既可去除噪声又可保护图像细节和边缘,能获得较好的图像效果。

1.2 提取岩石中的孔隙-图像分割

图像分割是指根据灰度、颜色、纹理、几何形状等信息把图像分成各具特点、互不重叠的区域并提取出感兴趣目标的技术和过程。即在一幅画图像中把目标从背景中分离出来,以便于对感兴趣目标进行更进一步分析。图2是所要处理的原始岩石孔隙结构图像,目标是提取出图像中的孔隙。

图2 原始图像(岩石背散射扫描电镜图像)

图3 原始图像的灰度直方图

图3是图2的灰度直方图,从直方图中可以看到明显的双峰,但双峰的灰度值相差极大,有宽且平的谷底。

根据直方图的特征,采用基于迭代的图像阈值分割算法[5]对原始孔隙图像进行分割,得到图像中的孔隙。迭代的图像阈值分割算法步骤如下。

第1步:求出图像最小灰度值Rmin和最大灰度值Rmax,计算初始阈值

T0=(Rmin-Rmax)/2

(1)

第2步:根据阈值将图像分割成目标和背景2部分,示出2部分的平均灰度值。

(2)

(3)

式中,R(i,j)为图像上(i,j)点的灰度值;N(i,j)为(i,j)点的权重系数,一般N(i,j)为R(i,j)的个数;TK为阈值。

第3步:重新选择阈值TK+1。

TK+1=(R0-RG)/2

(4)

第4步:循环第2步到第4步,当TK=TK+1则结束,即可获得最佳阈值来对进行分割。

图4是图2提取孔隙后的图像,其中绿色部分表示孔隙,非绿色部分(灰色部分)表示岩石。

图4 提取孔隙

1.3 去除孔隙二值图像的噪声

通常,有噪声的图像用阈值二值化所得到的边界往往是很不平滑的,物体区域具有一些错判孔洞,背景区域上则散布着一些小的噪声物体。这些噪声被称为前景噪声(如沙眼噪声)和背景噪声(如胡椒状噪声)。在识别喉道前,先要对孔隙的二值图像去噪。采用数学形态学的开运算[1][先腐蚀,见式(5),后膨胀,见式(6)、式(7)]去除胡椒状噪声

A⊖B={x:B+x⊂A}

(5)

A⊕B=[AC⊖(-B)]C

(6)

A∘B=[A⊖B]⊕B

(7)

采用闭运算[1](先膨胀,后腐蚀)[见式(8)],去除沙眼噪声

A·B=[A⊕(-B)]⊖(-B)

(8)

如果连续采用开和闭运算可以显著改善二值图像中的前、背景噪声的情况。

1.4 识别喉道-分割孔隙

岩石图像中,多个孔隙被喉道连接成连通区域,应用改进的识别喉道算法,把连通的孔隙分割开,以便能更准确计算每个孔隙的各项属性数据。改进识别喉道中轴线算法流程图见图5。

图5 改进的识别喉道中轴线算法流程图

A*B=(A⊖E)∩(AC⊖F),E∩F=∅

(9)

B∶A=∪{(B+h)∩A∶h∈H}

(10)

采用上述算法识别图4中孔隙的喉道,识别结果如图6所示。

图6 识别喉道后的图像(红色像素表示喉道)

1.5 计算孔隙各项属性-图像分析与理解

1.5.1 计算孔隙总数及孔隙面积

计算孔隙总数及孔隙面积的算法是用不同的颜色表示不同的孔隙,这为计算每个孔隙的周长及喉道提供基础。

(1) 查找未计算的孔隙像素。从上一个种子点(若刚开始扫描图像,则从图像的(0,0)像素点)开始对图像进行由左向右、由上向下扫描,如果扫描到未计算的孔隙像素(绿色像素)时,设置该像素为种子点,且孔隙总数加1,并执行(2),否则表示所有的孔隙已经计算完成。

(2) 计算孔隙颜色。计算得到一个图像中未使用过的颜色,且该颜色不能为绿色、红色和灰度色。

(3) 计算孔隙面积。根据(1)中查找到的种子点,使用(2)中计算出的孔隙颜色,采用计算机图形学的种子填充算法,记录下孔隙被填充的像素个数即为孔隙包含的像素个数。根据像素个数及比例值计算得到孔隙面积值。

1.5.2 计算每个孔隙的配位数及喉道长度

一个孔隙的配位数就是与一个孔隙连接的喉道个数。

第1步:查找喉道。从上一个红色像素(若刚开始扫描图像,则从图像(0,0)像素)开始对图像进行由左向右、由上向下扫描,如果扫描到喉道像素(红色像素)时,执行第2步,否则继续扫描图像,直至扫描完整个图像。

第2步:得到该喉道分割的孔隙。得到第1步找到的红色像素的4个近邻的颜色值,从这4个颜色值中找到该喉道分割开的孔隙颜色值,对这2个孔隙的喉道各加1,并根据喉道的2个端点值计算出喉道的长度。

第3步:删除喉道。把图像中该喉道所经过的像素点的颜色设置为白色(或设置为原始岩石孔隙结构图像中的颜色)。返回执行第1步。

1.5.3 计算孔隙周长

运用数学形态学的边缘跟踪(又称轮廓跟踪)算法计算每个孔隙的周长。边缘跟踪算法如下。

第1步:按从上到下、从左到右的顺序扫描图像,寻找没有标记跟踪结束记号的第1个边界起始点A0,A0是具有最小行和列值的边界点。定义一个扫描方向变量dir,该变量用于记录上一步中沿着前一个边界点到当前边界点的移动方向,其初始化取值:对4连通区域取dir=3;对8连通区域取dir=7。

第2步:按逆时针方向搜索当前象素的3×3邻域,其起始搜索方向设定:对4连通区域取(dir+3) mod 4;对8连通区域,若dir为奇数取(dir+7) mod 8;若dir为偶数去(dir+6) mod 8;

在3×3邻域中搜索到的第1个与当前像素值相同的像素便为新的边界点An,同时更新变量dir为新的方向值。

第3步:如果An等于第2个边界点A1且前1个边界点An-1等于第1个边界点A0,则停止搜索,结束跟踪,否则重复步骤2继续搜索。

第4步:由边界点A0、A1、A2、…、An-2构成的边界便为要跟踪的边界。

图7是以图6为基础计算孔隙各项属性后的图像。

图7 计算孔隙各项属性后的图像

2 岩石孔隙结构分析结果

根据每个孔隙的各项属性数据,可计算出如表1所示的18项孔隙结构参数表、表2所示的孔径分布测定数据。

表1 孔隙结构参数

表2 孔径分布测定数据表

3 在岩石粒度分析中的应用

该算法除了可以分析岩石孔隙图像中的孔隙外,还可分析岩石颗粒。

图8是图2提取岩石颗粒的结果。图9为分割岩石颗粒的结果。计算岩石颗粒的各项属性包括颗粒面积、周长、粒径,图10中一种颜色表示1个颗粒。

图8 提取岩石颗粒(绿色表示颗粒)

根据每个岩石颗粒的属性值可计算出如表3所示的9项粒度结构参数以及图11所示的粒度面积分布频率图。

表3 粒度结构参数

图9 岩石分割后的图像(红色像素为分割线)

图10 计算岩石颗粒各项属性后的图像

图11 粒度面积分布频率直方图

4 结束语

提出了一套完整的分析岩石孔隙结构图像的算法流程,特别是改进了识别喉道的算法,大大加快了计算速度,如在联想W500(处理器为Intel(R) Core(TM)2 Duo CPU P8600 @ 2.40 GHz,内存2 GB)的笔记本上,对图4(1 024像素×768像素)识别喉道,该算法用时0.853 s,流域分割方法用时2.156 s,中轴线算法用时2.785 s,灰度膨胀和腐蚀方法用时4.267 s。

应用本文算法开发的孔隙结构图像分析软件已在中国石油勘探开发研究院使用多年,分析各类样品数千个。实例应用表明,这些算法抗干扰能力强,可快速处理各种复杂类型的岩石图像,并准确获取孔隙结构、粒度参数等信息,有效支撑了相关科学研究与生产任务。

符号说明。T0:初始阈值;TK:迭代第K步阈值;Rmin:图像最小灰度值;Rmax:图像最大灰度值;R0:第K步图像目标的平均灰度值;RG:第K步图像背景的平均灰度值;R(i,j):图像上(i,j)点的灰度值;N(i,j):图像上(i,j)点的权重系数;A:输入图像;AC:为A的补集;B:结构元素或结构元素对(E,F);-B:将B相对原点旋转180°;E:探测图像内部的结构元素;F:探测图像外部的结构元素;H:为A的子集,第1次膨胀的种子点。

参考文献:

[1] 崔屹. 图象处理与分析-数学形态学方法及应用 [M]. 北京: 科学出版社, 2000.

[2] 印勇, 李阿琼. 一种粘连血细胞图像分割新方法 [J]. 计算机工程与应用, 2009, 45(35): 173-175.

[3] HU Dong, Mustafa Touati, Martin J Blunt. Pore Network Modeling: Analysis of Pore Size Distribution of Arabian Core Samples [J]. SPE 105156.

[4] 刘莉莉, 王铮. 一种适合血细胞图像分割的改进流域分割算法 [J]. 微电子学与计算机, 2010, 27(11).

[5] 王志明. 数字图像处理与分析 [M]. 北京: 清华大学出版社, 2012.

[6] 王雪莉, 张吉群. 高压地层微观模型系统中的背景校正 [J]. 科学技术与工程, 2012, 12(25): 6498-6502.

[7] Dmitry B Silin, Guodong Jin, Tad W Patzek. Robust Determination of the Pore Space Morphology in Sedimentary Rocks [J]. SPE 84296.

[8] AI Ibrahim M A, Hurley N F, Zhao W, et al. An Automated Petrographic Image Analysis System: Capillary Pressure Curves Using Confocal Microscopy [J]. SPE 159180.

[9] 桑卡, 赫拉瓦卡, 博伊尔. 图像处理、分析与机器视觉 [M]. 艾海舟, 苏延超, 译. 3版. 北京: 清华大学出版社, 2011.

[10] Suicmez V S, Touati M. Pore Network Modeling: A New Technology for SCAL Predictions and Interpretations [J]. SPE 110961.

[11] Liu Jianjun, Lin Lijun, Ji Youjun. Using Rock SEM Image to Create Pore-scale Finite Element Calculation Mesh [C]∥2011 International Conference on Physics Science and Technology (ICPST 2011).

[12] 熊承仁, 唐辉明, 刘宝琛, 等. 利用SEM照片获取土的孔隙结构参数 [J]. 中国地质大学学报: 地球科学, 2007, 32(3): 415-419.

[13] 尼克松, 阿瓜多. 特征提取与图像处理 [M]. 李实英, 杨高波, 译. 2版. 北京: 电子工业出版社, 2010.

[14] 帕科尔. 景丽译. 图像处理与计算机视觉算法及应用 [M]. 2版. 北京: 清华大学出版社, 2011.

[15] 赫恩, 巴克. 计算机图形学 [M]. 宋继强, 蔡敏, 译. 3版. 北京: 电子工业出版社, 2011.

[16] 汪灿, 刘艳敏, 祝艳波. SEM照片孔隙参数提取技术研究 [J]. 安全与环境工程, 2011, 18(3): 117-120.