基于主动轮廓的岩心FIB-SEM序列图像孔隙提取方法

户瑞林，滕奇志，何小海，龚剑

（1.四川大学电子信息学院图像信息研究所，成都 610065；2.成都西图科技公司，成都 610024）

0 引言

近来，由于世界常规油气产量的持续下降，致密气、页岩气、致密油等非常规油气资源越来越受到世界各国的关注[1]。非常规油气储集层以纳米级孔喉为主，局部发育在微米-毫米级，研究储集层的孔隙结构对于非常规油气的开采具有重要意义[2]。

聚焦离子束扫描电镜(FIB-SEM)是研究非常规油气储集层结构的一种新方法。FIB-SEM使用离子束对岩石进行连续剥蚀切割，同时利用电子束进行成像，能够还原真实岩心三维孔隙结构[2]。对岩心进行物相区分，提取岩心孔隙结构是数字岩石分析的核心步骤，但在对岩石的FIB-SEM图像进行孔隙提取时，存在着如下的难点：

（1）FIB-SEM成像时，由于观察面与电子束不垂直，底部信号比上部弱，扫描图由上至下会有变暗的效果，这在利用量化方式进行物相区分时会产生很大的干扰[4]；

（2）在对岩石成像时，由于每张图片都是分别拍摄，当前帧相对上一帧图像存在着一定的偏移，这在利用序列图像层间相关性进行孔隙提取时，会产生极大的影响；

（3）岩石中除孔隙之外的其他结构，如有机质、黏土矿物等，在FIB-SEM图像中也会呈现深色区域，这提高了区分孔隙和岩石的难度；

（4）由于SEM的分辨率较高，当岩石的观察面出现较大孔隙时，SEM成像会将其内部细节表现出来，再加上荷电的作用，孔隙内部一般还伴随着高亮的特点[4]，这意味着孔隙在图像中所呈现的形式有所不同，增加了鉴别孔隙的难度。

图1 致密碳酸岩FIB-SEM图片

由于各种岩石结构差别很大，FIB-SEM自身成像的特点又对孔隙提取造成干扰，目前尚没有一种统一的提取岩石孔隙的方法[4]，在实践中，经常需要根据岩石样品的特点，使用各种方式尝试。Martin Salzer[5]等在2012年提出了一种两阶段提取多孔材料FIB-SEM图像的方法，其利用FIB-SEM图像孔隙出现的高光效果，提取出高光区域，然后进行反向传播，但这是针对多孔材料的分割方法，要求孔隙外的其他部分不会出现明显灰度的差异，无法对含有多种杂质干扰的岩心提取孔隙。王羽[6]等在2016年对页岩FIB-SEM图像分别进行了边缘检测分割法、流域分割法和手动或自动阈值分割方法的实验，但效果依赖人为调整，序列图中孔隙提取效果表现不佳。

本文利用岩石FIB-SEM图像具有孔隙存在较为明显的边缘，每张图像变化差异不大的特点，针对孔隙个数较少的岩石，提出了一种基于主动轮廓的岩心序列图像孔隙提取方法。我们先在单张图像上初始化孔隙的边缘控制点，利用Snake算法收敛到当前孔隙的边缘，然后将该区域利用图像形态学操作，将边缘映射到下一张序列图，继续采用Snake算法，进行孔隙提取。由于岩石FIB-SEM图像中含有黏土矿物等干扰，在各帧图像进行Snake算法收敛时，部分会出现错误分割，本文针对该种情况提出了一种自动寻找最佳分割效果的方法，最终提取出序列图像孔隙，提取效果良好。

1 预处理

由于FIB-SEM的成像特点，所生成的FIB-SEM图像会产生位移。我们可以发现，相邻两张岩石FIBSEM图像的真实观察面变化很小，若图像未产生位移，则两张图像的差值图整体灰度应较小，故我们可将后一张图像逐像素移动，寻找使得其与前一张图像的差值图整体灰度最小的位置，则此时应为配准后的图片。

本文实验图为致密碳酸岩的FIB-SEM图像，我们在图像的中间设定一个大小为500×500像素的模板窗口，取当前帧图像在模板窗口位置的部分为I，令为将下一帧沿x方向移动m个像素，沿y方向移动n个像素后模板内的图像，令：

Gm,n表示两帧的差值图像的所有像素灰度绝对值之和，当Gm,n最小时，我们认为此时的I'm,n为与I配准的图像，一般而言FIB-SEM图像的偏移不会超过50个像素，故我们可以取m∈[- 50,50],n∈[- 50,50]。

我们取配准前后序列图的一个x-z切面观察，可以发现，原始图像的切面图像连续性很差，这就是由于图像间产生了位移导致的结果，而配准后的图像明显比配准前的图像更加平滑，连续性更好，这说明了配准的图像更加接近真实切面的结构。

图2 配准前后的x-z面图像

由于后续操作是需要利用边缘来提取孔隙，我们需要对配准后的图像进行边缘增强操作。FIB-SEM图像存在一些小噪点，我们先进行中值滤波，然后再求图像的梯度图像，简单起见，我们取梯度图像：

其中gx，gy为图像分别在x，y方向上的微分，将梯度图像叠加到原图，这样就得到了增强后的图像。

图3 增强前后的图像

2 FIB-SEM孔隙提取

2.1 Snake模型

Snake模型是Kass[7]在1987年提出的一种主动轮廓线的图像分割算法，其定义了一条受到内、外部约束力的轮廓曲线，通过轮廓曲线的变形，让其能量极小化，最终使得曲线逼近目标的边缘。Snake模型由一组控制点：

组成，其中x(s)和y(s)是控制点的坐标位置，s是以傅里叶变换形式描述边界的自变量。轮廓曲线的能量函数定义为：

一个作家创作观的形成，是一个长期积淀的过程，他的生活经历对他的个性、气质、思维方式的形成和发展起着决定性的作用。“特别是那些印象深刻的经验往往给艺术家的一生涂上一种特殊的基调和底色，并在相当程度上决定着艺术家对于创作题材的选择和作品的情绪情感或情节。”[8]

在该定义中，Eint为内部约束能量，Eimage为图像力，Econ为外部约束力，而内部约束能量可以表示为：

在极限的情况下，弹性能量会把轮廓线压缩成一个点，弯曲能量会把轮廓线变为一个光滑的曲线或直线[8]，图像力会驱使各控制点到其附近梯度最大的位置。通过移动各控制点的位置使总能量最小化的过程中，轮廓线就会逐渐逼近目标的边缘。Kass提出通过求解欧拉方程来使得总能量最小化，但该方法需要计算高阶导数，计算量较大，对于FIB-SEM序列图像的处理速度较慢。故本文采用了贪婪算法，每次移动一个控制点，将移动控制点到使得能量函数最小的位置，然后再移动下一个控制点，如此循环移动，最终可使得轮廓线逼近孔隙边缘。

2.2 Snake模型提取单张图像孔隙

图4 主动轮廓提取孔隙

由于各类岩心的图像相差很大，孔隙结构也大有不同，甚至部分孔隙人眼也无法识别，所以第一张种子图像的轮廓我们由人为给定，如图4，在孔隙周围设定控制点，然后自动插值，增加控制点个数，以尽可能逼近孔隙轮廓。然后Snake模型进行收缩，最终逼近孔隙真实轮廓。

2.3 提取序列图像的孔隙

由于FIB-SEM图像每帧间差异变化不大，而Snake算法会自动收缩至边缘，故我们可以考虑将当前提取完成的区域进行膨胀，然后将其轮廓映射到下一帧图像作为其初始轮廓。设我们当前帧所提取的孔隙区域为M，膨胀结构元为B，下一帧图像的初始轮廓为G'，G'所围成的区域为M'，定义⊕为膨胀操作，则M'=M⊕B,G'=boundary(M')，其中boundary为取区域边界轮廓。

在多数情况下该方法能够获得较好的结果，但膨胀使用固定的结构元，存在提取结果不理想的情况：当黏土矿物等杂质靠近孔隙边缘时，其也具有强边缘，若初始化区域膨胀过大，Snake模型在收缩时会受到杂质干扰，容易收敛到杂质的边缘上；当孔隙变化较为剧烈时，若初始化区域膨胀过小，其轮廓覆盖不够大，则Snake模型容易收缩到孔隙内部，而出现错误分割。

图5 孔隙提取较坏情况

本文提出一种自动寻找孔隙提取最优效果的方法，其可以动态调整膨胀结构元，使得孔隙提取区域达到最佳。思路如下，取相邻两帧的同一孔隙的分割结果，其差异区域为T，T的面积为Ta，由于两帧之间的大部分孔隙变化差异不大，若Ta越小，我们越倾向于该次分割效果较好。同时，我们取相邻两帧的差值图像为Tsub，其包含两帧间的变化信息，由于两帧之间的变化集中在边缘部分，取Tsub内的T区域的灰度值之和为S，若S越大，则本次分割变化部分越接近边缘，我们倾向于本次分割效果越好。由于S和Ta的量纲不同，我们可以将其归一化后作差得R，若R值越大，则我们倾向于该次分割效果越好。算法具体步骤如下：

（1）设当前帧图像为I，下一帧为I'，两帧差值为为圆盘形结构元，n为圆盘半径，其范围可根据具体情况指定，一般可取n∈[ ]1,10。使用Bn对当前提取的孔隙区域M分别膨胀，得

（4）取Isub内的Tn区域内的所有灰度值之和为STn，分别取 ATn，STn在 n∈[ ]1,10的最大值和最小值Amax，Amin，Smax，Smin，归一化得:

图6展示了利用上述算法的过程，图6中，图(e)为当前帧的提取区域，图(f)为下一帧采用结构元半径为5的提取区域，其R值为0.18，图(g)为下一帧采用结构元半径为7的提取区域，其R值为0.05，由于R值越大，效果越好，故图(f)的提取孔隙效果更优。视觉上看，图(f)确实较图(g)更优，其原因是，图(g)与图(f)的差异部分对图(g)而言，其所增加的A远大于S，所以图(g)所提取的区域的R值较小。

3 实验结果与分析

为验证本文所提出的算法的正确性，将本文的算法与阈值分割，Martin Salzer[5]文中所提出的反向传播算法作比较。阈值分割是常用的图像分割方法，其原理是设定分割区间，通过灰度值的差异来进行图像分割，Martin Salzer[5]文中所提出的反向传播算法是针对FIBSEM图像提出的一种特定算法，在分割多孔材料中表现良好。本文以致密碳酸岩作为样本，分别使用三种算法进行试验，可以发现，由于FIB-SEM存在由上至下的灰度变暗的现象，而且孔隙内部灰度不均匀，使用阈值分割，容易受外部黏土矿物的灰度影响，表现不好。而反向传播算法，由于黏土矿物的边缘也会出现类似孔隙内部灰度突变，而误将其也认为孔隙，提取效果也比较差，而利用本文所提出的算法，观察孔隙提取结果，能够明显看出是三种算法中表现最好的。我们也可以取其中一个孔隙，人工识别标注该孔隙区域，与各算法所提取的区域相异或，这样可以看出不同的部分。观察图7(s)、(t)、(v)，可以发现本算法所提取的区域与人工识别区域差别较小，提取效果良好。

图6 自动调整孔隙提取效果

图7 不同算法的效果对比

4 结语

本文提出了一种基于主动轮廓的岩心FIB-SEM序列图像孔隙提取方法，在初始化首张图像的孔隙轮廓后，利用FIB-SEM图像的层间相关性，通过主动轮廓算法自动提取孔隙，同时也给出了一种自动寻找序列图像中孔隙最优提取效果的算法。本文方法优于其他对比算法，能够准确地提取出孔隙区域，在目前尚没有一种统一的提取岩石孔隙的方法的情况下，本文所提出的方法对岩心FIB-SEM图像具有一定的通用性。但本文方法需要人为给出孔隙的初始轮廓，还不能够完全自动化提取孔隙，对于孔隙较多的岩石而言，工作量还比较大，因此可在这些方面进一步完善。