基于数字岩心的致密砂岩微观孔喉结构定量表征

胡 渤，蒲 军，苟斐斐

（1.中国石化油田勘探开发事业部，北京 100728；2.中国石化石油勘探开发研究院，北京 101599）

中国致密砂岩油藏资源丰富，初步评估可采储量达35×108～40×108t［1-6］，约占世界致密油总可采储量的十分之一，是未来现实的资源接替类型。致密砂岩储层多发育微米级孔隙和纳米级喉道，孔喉半径级差大、形态不规则、结构复杂、微观非均质性强，极大地影响了储层的渗流能力［7-13］，定量描述微观孔喉结构特征可以为致密砂岩储层质量差异化表征提供依据，也可以为致密砂岩储层分类及开发对策研究奠定基础，具有重要的意义。但常规实验方法（压汞、核磁、驱替等）存在实验周期长、精度低等问题［14-16］，难以准确地确定微观孔喉结构参数，制约了致密砂岩孔喉结构特征研究。

近年来，随着扫描成像设备与计算机技术的快速发展，基于高精度扫描图像构建数字岩心，建立多孔介质模型，开展统计分析，为研究致密砂岩储层孔喉结构特征提供了新的技术手段［17-19］。目前业界内主要采用孔隙网络模型方法来分割孔喉并计算孔喉特征参数。该方法是将岩石孔喉空间简化为球棍模型，用圆球代表孔隙，孔隙之间用细棍代表喉道。常用的孔隙网络模型方法主要有最大球法［20-21］和中轴线法［22-26］。最大球法需要计算每个像素的最大球，最终生成的孔隙网络模型是对真实孔喉的一种等效假设，计算量较大，孔喉定量表征精度低。中轴线法受孔喉形态影响大，致密砂岩孔喉空间复杂多变，提取的中轴线会产生大量分支，也影响了孔喉分割与表征的精度。

直接在数字岩心图像的基础上开展孔喉结构定量表征，无需对孔喉空间进行简化处理，可以最大限度地保证孔喉表征的精度与可靠性，成为当前研究的热点与前沿，但目前仍然存在孔喉提取效率低、精细分割难度大等问题。笔者从图像形态学出发（以SEM-Maps 扫描图像为例，分辨率为10 nm），开发了自适应孔喉识别与提取方法，研发了直接基于图像形态学算法的孔喉分割方法及孔喉结构特征参数计算方法，并以红河油田长8 储层典型井致密砂岩岩样为例，分析了不同储层的微观孔喉结构及连通性差异，研究成果可以为致密砂岩储层质量表征提供依据和借鉴。

1 致密砂岩数字岩心图像处理

数字岩心图像（包括SEM-Maps 和CT 等）的灰度值为0～255，灰度值较高的像素反映的是相对密度和原子序数较大的物质，如岩石骨架；灰度值较低的像素反映的是相对密度和原子序数较低的物质，如孔喉空间［9-10］。数字岩心图像处理的目的是提高图像清晰度并提取孔喉空间，为下一步孔喉结构定量表征做好准备。

1.1 图像降噪

数字岩心图像在数字化和传输过程中常受到成像设备与外部环境噪声的影响，图像中含有大量噪点，如图1所示，这些噪点极大地影响了孔喉识别与提取的准确性［11-12］。因此，孔喉空间识别与提取的首要步骤是进行降噪处理。

图1 数字岩心图像噪点Fig. 1 Noises in digital core images

系统分析数字岩心图像噪点分布特征（以SEM-Maps 扫描图像为例），发现图像中的噪点主要为高斯噪点（图2），且不同矿物颗粒和孔喉空间具有不同的噪点标准差。图1 中1—4 位置处的噪点标准差分别为7.82，8.48，0 和6.87，选择合适的降噪算法有效去除噪点是获取高质量数字岩心图像的关键。

图2 SEM-Maps扫描图像噪点分布（图1位置1处）Fig. 2 Distribution of noises in SEM-Maps scanning image（location1 on Fig.1）

双边降噪算法是在高斯降噪算法中加入图像灰度权重项，不但可以有效去除高斯噪点，同时能够保持孔隙和骨架边界清晰［13］，成为数字岩心图像降噪的首选算法。

设f(i,j)表示原图像灰度值，g(i,j)表示降噪后图像灰度值，则双边降噪公式为：

根据图像噪点统计结果，取统计标准差的最大值σd=σr=8.48。双边降噪前后数字岩心图像如图3所示，矿物颗粒和孔喉空间的噪点大幅度减少，相同灰度阈值下提取的孔隙更加准确（灰度阈值为150，红色为孔喉）。

图3 降噪前后数字岩心图像孔喉空间提取结果对比Fig. 3 Comparison of extraction results of pore throat spaces in digital images before and after noise reduction（Threshold is150，and the red part is pore throat）

1.2 孔喉提取

消除噪点后，可以对孔喉空间进行识别和提取。数字岩心图像孔喉提取的关键在于确定合适的灰度分割阈值，常规方法是通过人工观察确定分割阈值，但准确性较低。笔者提出了自适应孔喉识别与提取方法（图4），通过对降噪后数字岩心图像的灰度值进行统计，结合大津法，自动计算得到合适的灰度分割阈值，避免了人为因素的干扰。

图4 自适应孔喉识别与提取方法流程Fig. 4 Process of adaptive pore throat recognition and extraction

对降噪后数字岩心图像中的所有像素点的灰度值进行统计，得到灰度累积分布曲线（图5）。灰度累积分布计算方法与面孔率（孔隙度）的计算方法一致，因此，图5的纵坐标可以理解为对应灰度值的面孔率（孔隙度）。

图5 数字岩心图像灰度与孔隙度的对应关系Fig. 5 Relationship between gray value and porosity in digital core image

致密砂岩岩心孔隙度通常大于4%，通过数字岩心图像观察，结合岩心孔隙度分析，确定灰度值小于50 的像素点为孔喉空间（图5）；当灰度值大于165 时，孔隙度急剧增大，这是因为灰度值大于165时，大量岩石矿物颗粒被误识别为孔喉。因此，合适的孔隙/岩石骨架灰度分割阈值为50～165。

为准确地确定灰度分割阈值，可以对数字岩心图像进行灰度拉伸：

对灰度拉伸后的图像（图6a）采用Ostu 算法［14］，以像素灰度值级差最大为分割依据，得到孔喉空间与岩石骨架灰度分割阈值界限为141，分割得到的结果如图6b 所示，白色为岩石矿物颗粒，黑色为孔喉空间。

图6 灰度阈值分割后的数字岩心图像Fig. 6 Digital core images after grey level threshold segmentation

2 基于数字岩心图像的孔喉分割

常规砂岩孔喉级差小、形态规则，采用孔隙网络模型分割的孔隙和喉道与实际情况基本相符（图7）。但致密砂岩岩心孔喉级差大、形态复杂多变，采用孔隙网络模型易将大孔隙识别为多个小孔隙的组合，且易将喉道中略宽的部分识别成小孔隙，这将极大地影响孔喉定量表征结果的准确性。

图7 不同类型砂岩岩心孔隙网络模型提取结果对比Fig. 7 Comparison of extraction results of pore network models in different sandstone cores

从致密砂岩孔喉形态差异性出发，基于图像形态学原理，采用腐蚀、膨胀及逻辑算法建立了一套适用于致密砂岩的孔喉分割方法。以致密砂岩数字岩心典型孔喉图像（图8）为例，将孔喉分割的具体流程说明如下。

步骤1：采用腐蚀算子消除喉道。在图像形态学中，设U为一个平面区域，设A1为区域U内部的1个目标区域，S1为指定大小和形状的结构元素1，定义位于坐标(x,y)的S1所表示的区域为S1(x,y)，则对目标区域A1的腐蚀操作可以表示为：

S1的大小取决于数字岩心图像中喉道的大小，通过观察可以看到：由于喉道可以沿着任意方向延伸，结构元素的形状通常采用圆盘状，这样可以保证所有喉道均可以被腐蚀（图9）。

图9 直径为5个像素的圆盘状结构元素示意Fig. 9 Disk-shaped structural elements with a diameter of 5 pixels

在图8 中，A1为原始孔喉，采用S1腐蚀后，不仅所有喉道被腐蚀（图10a），孔隙周边像素也会被腐蚀。因此，图10a中显示的并非真实孔隙，这里称之为残缺孔隙（A2）。

图8 致密砂岩数字岩心典型孔喉图像（A1）Fig. 8 Typical pore throat image of tight sandstone based on digital core technology（A1）

步骤2：采用膨胀算子修复残缺孔隙。在图像形态学中，设U为一个平面区域，S2为指定大小和形状的结构元素2，定义位于坐标()x,y的S2所表示的区域为S2()x,y，则对目标区域A2的膨胀运算可以表示为：

对残缺孔隙（图10a）进行膨胀运算修复后，被腐蚀的孔隙边界像素得到了恢复（图10b）；由于膨胀运算并不是腐蚀运算的逆运算，因此这种恢复并不准确。为了确保膨胀运算后的孔隙能够恢复到原始孔隙，膨胀运算的S2的直径要略大于步骤1 中腐蚀运算的S1，但S2的结构与S1保持一致即可，膨胀运算修复后的孔隙为A3。

图10 数字岩心孔喉图像腐蚀膨胀处理结果Fig. 10 Processing results of corrosion and expansion in pore throat images

步骤3：采用交集运算提取原始孔隙。采用膨胀运算修复后的孔隙A3与原始孔喉A1进行交集运算，即可得到准确的孔隙A4（图11a）。交集运算可以表示为：

步骤4：采用差集运算提取喉道。分割得到准确的孔隙后，从原始孔喉图像中用差集运算去除孔隙部分，即可得到准确的喉道A5（图11b）。差集运算可以定义为：

图11 基于数字岩心图像的孔喉分割结果Fig. 11 Segmentation results of pore throats based on digital core images

3 孔喉定量表征

孔隙和喉道均为几何图形，具有相同的几何特征，例如面积、周长等。此外，孔隙和喉道在形态上又具有较大的差异，存在特征参数，例如形状因子、配位数等。

3.1 孔喉共有参数计算

采用GRAY提出的模板匹配方法来计算孔喉的面积和周长［26］。在结构元素模板中，1 表示目标像素，0表示空白像素（图12）。

图12 结构元素模板Fig. 12 Template of structural element

采用n{ }Q表示孔喉图像中匹配的模板图像Q的个数。孔喉面积和周长的计算公式分别为：

3.2 孔隙特征参数计算

3.2.1 形状因子

形状因子用来描述孔隙接近圆的程度，即圆的形状因子为1。由于前文已经计算得到孔隙的面积和周长，所以这里用面积和周长来描述孔隙的形状因子。对于周长为Pp的圆，其面积为：

则可以定义孔隙的形状因子为：

3.2.2 面孔率

面孔率是指数字岩心孔隙的面积与总面积的比值。假设数字岩心的总像素个数为N，孔隙的像素个数为M，则面孔率的计算公式为：

3.3 喉道特征参数计算

3.3.1 长度和平均宽度

已知喉道的面积At和周长Pt，设喉道长度为L，平均宽度为d，则有：

求解（15）和（16）式，可得喉道的长度和平均宽度分别为：

3.3.2 配位数

在储层孔隙结构表征中，配位数是指每个孔隙连通喉道的数量，通过对孔喉图像进行运算，可以得到孔喉连接关系矩阵cIJ，配位数计算公式为：

3.4 孔隙连通性计算

岩心中并非所有孔隙都是相互连通的。根据孔隙的连通关系，可以将数字岩心图像划分为不同的连通体，而最大连通体的大小直接决定岩心的渗透率，因此可以最大连通体面积占比来表征孔隙连通性：

4 实例应用

基于红河油田长8 储层物性差异选取3 口典型井，分别采集岩样进行SEM-Maps扫描，构建数字岩心，定量计算孔喉结构参数，并对比分析不同储层微观孔喉结构的差异性。

4.1 岩样选取与数字岩心构建

基于沉积与成岩作用研究，结合岩心物性分析测试结果，选择红河油田长8储层典型井岩样，基本信息见表1所示。

表1 红河油田长8储层典型井岩样基本信息Table1 Information of rock samples from Chang8 reservoir in Honghe Oilfield

对3 块岩样进行表面处理，选择特征区域进行SEM-Maps 扫描成像，得到分辨率为10 nm、大小为1 mm×1 mm 的8 位灰度图（图13）。其中，46 号岩样的粒间孔较发育，且存在大孔隙；而73 号岩样主要发育晶间孔，HH1 号岩样同时发育粒间孔和晶间孔，但均未见到大孔隙。

图13 红河油田长8储层3块典型井岩样SEM-Maps扫描图像Fig. 13 SEM-Maps scanning images of three typical rock samples from Chang8 reservoir in Honghe Oilfield

按照前文所述方法，对研究区3 块典型井岩样的SEM-Maps扫描图像进行处理，构建数字岩心，并提取孔隙和喉道。结果如图14所示，白色部分为孔喉空间，黑色为矿物颗粒。

图14 红河油田长8储层3块典型井岩样数字岩心孔喉图像Fig. 14 Pore throat images of three typical rock samples from Chang8 reservoir in Honghe Oilfield

4.2 孔喉结构定量表征

在红河油田长8储层典型井岩样数字岩心图像孔喉提取的基础上，开展定量表征与分析，明确不同油井岩样的孔隙、喉道发育特征及孔喉连通性特征，为综合评价储层质量提供重要依据。

4.2.1 孔隙发育特征

基于SEM-Maps 扫描图像对岩样进行孔隙分割（图15），并对孔隙进行定量表征，对比分析3类储层的平均孔隙半径、平均形状因子及面孔率，结果（图16）表明，46 号岩样的平均孔隙半径最大，为1.92 μm；面孔率最高，为0.068；形状最规则，其平均形状因子为6.56。HH1 号岩样的平均孔隙半径最小，为0.91 μm；面孔率最低，为0.038；孔隙不规则，其平均形状因子为12.02。73 号岩样的孔隙定量表征参数介于46 号和HH1 号之间。因此，46 号岩样的储集性能明显优于73号和HH1号岩样。

图15 红河油田长8储层3块典型井岩样孔隙分割图像Fig. 15 Pore segmentation images of three typical rock samples from Chang8 reservoir in Honghe Oilfield

图16 红河油田长8储层3块典型井岩样孔隙特征参数对比Fig. 16 Comparison of pore characteristic parameters of three typical rock samples from Chang8 reservoir in Honghe Oilfield

4.2.2 喉道发育特征

基于SEM-Maps 扫描图像对研究区岩样进行喉道分割（图17），量化对比不同岩样的平均喉道长度、平均喉道宽度及平均配位数，结果（图18）表明，46 号岩样的平均喉道长度最短，为14.92 μm；平均喉道宽度最大，为0.35 μm；平均配位数最大，为3.95。HH1 号岩样的平均喉道长度最长，为25.05 μm；平均喉道宽度最小，为0.22 μm；平均配位数最小，为1.72。因此，研究区46 号岩样的流体渗流能力远优于73号和HH1号岩样。

图18 红河油田长8储层3块典型井岩样的喉道特征参数对比Fig. 18 Comparison of throat characteristic parameters of three typical rock samples from Chang8 reservoir in Honghe Oilfield

4.2.3 孔喉连通性特征

基于SEM-Maps 扫描图像对研究区3 块典型井岩样提取连通体，并用不同颜色标记，结果（图19）显示，各连通体之间相互独立。46 号岩样最大连通体面积占比为33.04%，基本贯穿整个岩样，表明流体能够通过整个样品；73 号岩样最大连通体面积占比为27.69%，连通性相对较差；HH1 号岩样最大连通体面积占比为8.36%，表现为局部连通，整体不连通。因此，46 号岩样连通性明显优于73 号和HH1号岩样。

图19 红河油田长8储层3块典型井岩样的孔喉连通性分析Fig. 19 Pore throat connectivity analysis of three typical rock samples from Chang8 reservoir in Honghe Oilfield

5 结论

基于孔喉形态学差异，综合应用图像腐蚀、膨胀和逻辑运算，提出一种基于图像形态学特征的致密砂岩数字岩心图像孔喉分割方法，可以准确地分割孔隙和喉道；建立直接基于数字岩心图像的孔隙、喉道以及孔喉连通性定量表征方法，并提出各表征参数的具体计算方法，实现了致密砂岩微观孔喉结构定量表征，为分析致密砂岩储层质量提供了依据。对比分析红河油田长8 储层3 口典型井岩样的孔喉结构特征，结果表明：物性最好的46 号岩样的平均孔隙半径最大、面孔率最大，平均喉道长度最小、配位数最高，最大连通体面积占比最大；物性最差的HH1 号岩样的平均孔隙半径最小、面孔率最小，平均喉道长度最大、配位数最低，最大连通体面积占比最小。

符号解释

A——孔喉的面积；

A1——区域U内部的1个目标区域；

A2——腐蚀运算后的残缺孔隙；

A3——膨胀运算后的残缺孔隙；

A4——准确的孔隙；

A5——准确的喉道；

Ab——第b个连通体像素个数，个；

Ap——孔隙的面积；

As——图像中孔喉的总像素个数，个；

At——喉道的面积；

b——连通体的个数，个；

cIJ——孔喉连接关系矩阵；

CNm——配位数；

d——喉道的平均宽度；

f（i，j），f（k，l）——原图像灰度值；

F——孔隙的形状因子；

g（i，j）——降噪后图像灰度值；

h（i，j）——灰度拉伸后图像的灰度值；

i，k——图像x方向坐标；

I——孔喉连接关系矩阵的行数；

j，l——图像y方向坐标；

J——孔喉连接关系矩阵的列数；

L——喉道的长度；

m——孔隙的个数，个；

M——孔隙的像素个数，个；

n{ }Q——孔喉图像中匹配的模板图像Q的个数，个；

N——数字岩心的总像素个数，个；

P——孔喉的周长；

Pp——孔隙的周长；

Pt——喉道的周长；

Q1，Q2，Q3，Q4，Q5——序号为1，2，3，4，5的模板图像；

S1——结构元素1；

S2——结构元素2；

T——与1个喉道连接的孔隙个数，个；

U——平面区域；

w——降噪卷积核；

wd——空间域核；

wr——灰度值域核；

w（i，j，k，l）——降噪卷积核w在i，j，k，l处的值；

wd（i，j，k，l）——空间域核wd在i，j，k，l处的值；

wr（i，j，k，l）——灰度值域核wr在i，j，k，l处的值；

x——笛卡尔坐标x轴位置；

y——笛卡尔坐标y轴位置；

α——面孔率；

β——最大连通体面积占比，%；

σd——图像噪音空间域标准差；

σr——图像灰度值域标准差；

Φ——处理后的剩余元素。