基于改进分水岭算法的岩屑扫描图像粒径分析新方法*

戴 宗 罗东红 谢明英 唐 放 李海龙 涂志勇 王晨晨 陈 强

(1. 中海石油(中国)有限公司深圳分公司研究院 广东深圳 518000； 2. 数岩科技(厦门)股份有限公司 福建厦门 361000)

在剩余油潜力较大的老油田再开采过程中，由于部分室内实验基础数据时间久远参考意义不大；而岩心库中仅存的储层代表性岩心往往由于老化严重、剩余量少，不具备重新进行物性参数、粒径分析等室内实验测试的条件，进而无法准确体现储层岩石物理特征，这对老油田挖潜带来了较大的困难和挑战。随着计算机图像处理技术、计算方法的发展，以及CT扫描等技术在岩石物理实验中的应用，可以构建能够反映储层岩石真实孔隙空间结构特征的数字岩心，为岩石微观结构特征和渗流模拟提供研究平台[1-4]。数字岩心建模方法主要包括以CT扫描为代表的物理建模法和数值重建法[5-8]。常见的数值重建算法主要有随机法和地质过程模拟法两大类。随机重建法是一种借助少量的岩心平面图像资料，通过图像分析技术提取建模信息，然后应用某种数值计算方法建立数字岩心的方法，常用的随机算法主要有：完全随机法、高斯场法、顺序指示模拟法、模拟退火算法、多点统计法和马尔科夫链-蒙特卡洛算法等。地质过程模拟法通过对颗粒的沉积过程及后续的压实作用、成岩作用进行模拟建立数字岩心[9-10]。

对于不具备完整的三维空间结构特征和二维孔隙空间特征的岩屑样品，不能直接通过CT扫描或者随机重建方法来构建三维数字岩心，因此需要通过地质过程模拟法来构建代表性的三维数字岩心。开展地质过程模拟前，需要首先测量岩心的粒度组成，在此基础上模拟颗粒的沉积、压实和成岩过程，构建三维数字岩心模型。本文针对储层代表性的岩屑样品，通过扫描电镜成像可生成MAPS，在此基础上提出一种粒径分析新方法获取岩屑粒径分布，并通过地质过程模拟法构建代表性的三维数字岩心，分析储层微观孔隙结构特征，为老油田剩余油挖潜提供基础数据。

1 基于3D距离图的分水岭算法

MAPS图像是在选定区域内排布扫描出一系列连续且边缘重叠的大量高分辨率的小图像，扫描完成后将这些小图像进行拼接进而得到一张超高分辨率、超大面积的二维背散射电子图像。基于MAPS图像通过分水岭算法进行分割，可以研究不同尺度下的岩石颗粒及孔隙之间的关联关系。分水岭算法是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。

然而对于岩屑样品，基于传统的分水岭算法对MAPS图像分割不能有效地识别出不同的颗粒边缘。本文通过计算图像的3D距离图，提出一种基于3D距离图的分水岭算法，利用该算法可以对岩屑MAPS图像粒径进行精确分割和计算，具体步骤如下：根据照片中的特征灰度值或形态，对原始照片进行特征分割，划分出颗粒相、基质相和孔隙相；对要统计的特征进行标记，如统计矿物颗粒时，需要对颗粒相进行标记；对特征属性进行计算，如颗粒大小(面积、体积、周长等参数)；对特征参数进行分组统计，如颗粒的等效直径等[11-12]。下面分别对以上过程进行阐述。

1.1 特征分割

图1 基于改进分水岭算法的粒径分割工作流程Fig .1 Particle segmentation workflow with modified watershed algorithm

1.2 特征标记

特征分割后，需要对图中的各个特征(每个连通区域为一单独的特征)进行编号，然后再分别计算各个特征的属性。Haralick给出了特征标记算法[11]。图2为特征标记示意图，图中分别用1、2、3标记了三个特征区域;此外,0为背景区域，作为特征之间的分隔线。

1.3 特征属性计算

二维图像特征通常有以下属性，如面积、边界周长、形状因子、重心位置等。在分析颗粒相时，常用等效圆半径R作为颗粒大小的一个度量。其定义为

(1)

式(1)中：A为颗粒面积。

图2 特征标记示意图Fig .2 Characteristic mark schematic diagram

1.4 特征参数分组统计

对图像中所有颗粒特征的等效圆半径进行计算，得到等效圆半径列表。然后对颗粒的等效圆半径按照大小分成若干组，利用直方图工具开展统计，获得颗粒相尺寸的频率分布和累积分布图。

2 实例验证与分析

图3(左)为南海东部海域某老油田A油田目标区块岩屑样品的背散射MAPS图像，图像中颜色较深的部分为孔隙，灰色的部分为矿物颗粒。背散射图像的灰度与矿物的平均原子序数相关，平均原子序数较小的其灰度值也较小，可以据此来观察/粗估矿物/孔隙在岩心表面的分布情况。图3(右上)是选取的一小块区域进行观察，并进行方法及参数测试。测试可行后应用到整幅图像中，获得整张MAPS图像的统计分析结果。对这一小块区域进行特征分割，获得了颗粒相(图3右下白色区域)和孔隙相(黑色区域)。需要注意的是，由于部分颗粒之间有接触，在划分颗粒相时采用常规的阈值分割无法分离这些有接触的矿物碎屑颗粒，这时如果不进行处理，势必会高估矿物颗粒的大小，并低估颗粒的数量，从而影响计算结果。因此，在进行颗粒相的统计时，还需要进行分离接触的颗粒相这一操作。

图3 南海东部A油田岩心样品的背散射MAPS图像及其孔隙、颗粒分割图(MAPS图像大小16 mm×15 mm)Fig .3 Back-scattered MAPS and pore-grain image segmentation of A oilfield in eastern South China Sea (MAPS image size 16 mm×15 mm)

图4为示例背散射图像中矿物颗粒的分离效果图，左图为原始背散射图像，中图为颗粒和孔隙相分割图(其中白色为颗粒相、黑色为孔隙相)，右图为采用基于3D距离图的分水岭方法得到的颗粒分离图，互相接触的颗粒之间已经用分水岭线分隔开。

图4 南海东部A油田MAPS图像中矿物颗粒分离Fig .4 Grain separation of MAPS image of A oilfield in eastern South China Sea

基于分离的颗粒图像，进行特征分割、特征标记、特征属性计算、特征参数分组统计，最后应用到整张MAPS图像中可以获取岩屑样品的颗粒粒径频率分布图(图5)。同时，利用传统室内实验粒径分析结果对其进行了验证，两者较为吻合(图5)。

图5 南海东部A油田颗粒粒径的频率分布图Fig .5 Frequency distribution of particle size of A oilfield in eastern South China Sea

基于岩屑样品的颗粒粒径分布，通过地质过程模拟法模拟颗粒的沉积过程、压实作用和成岩作用，通过已有的地层孔隙度数据建立储层代表性的三维数字岩心(图6a)，并提取相应的孔隙网络模型[13](图6b)。

图6 基于地质过程模拟法构建的南海东部A油田三维 数字岩心(模型大小：10 mm×10 mm×10 mm)Fig .6 Reconstruction of 3D digital rock of A oilfield in eastern South China Sea with geological process-based method(model size:10 mm×10 mm×10 mm)

以数字岩心和孔隙网络模型为研究平台，可以计算出储层的孔隙度为24.36%，渗透率为8 105 mD，并分析储层微观孔隙结构特征。储层孔喉半径分布11～200 μm(图7a)，主平均值为40 μm；形状因子分布在0.022～0.079(图7b)，平均值为0.055；配位数分布在1～12(图7c)，平均值为4.48。其中，形状因子是为了定量表征孔隙空间中孔隙、喉道单元体的形状特征，其定义如下：

G=Ap/P2

(2)

式(2)中：Ap为孔隙、喉道的横截面面积，m2；P为孔隙、喉道单元体横截面的周长，m。

图7 南海东部A油田孔隙结构特征分布Fig .7 Pore structure characteristic distribution of A oilfield in eastern South China Sea

3 结论

通过改进分水岭算法对颗粒图像进行计算，包括特征分割、特征标记、特征属性计算、特征参数分组统计等，最后应用到整张MAPS图像中可以获取岩屑样品的颗粒粒径频率分布图。该方法有效解决了岩屑粒径测量困难的现状，并通过地质过程模拟法建立储层代表性的三维数字岩心，提取相应的孔隙网络模型，计算孔渗物性参数和孔隙-喉道结构特征，可为缺少岩心数据的老油田剩余油挖潜提供基础数据。