基于蚁群和Canny 边缘检测算子混合算法的二维岩石图像裂隙特征提取与修复研究

张愉玲 , 邢会林 , , 李三忠 , , 逄 硕 , , 刘骏标 , 张熔鑫

(1.深海圈层与地球系统前沿科学中心, 海底科学与探测技术教育部重点实验室, 中国海洋大学 海洋高等研究院/海洋地球科学学院, 山东 青岛 266100; 2.中国海洋大学 海底科学与工程计算国际中心, 山东 青岛266100; 3.青岛海洋科学与技术国家实验室 海洋矿产资源评价与探测技术功能实验室, 山东 青岛 266100)

0 引 言

利用CT 扫描等对岩石内部结构进行成像是数字岩心、矿物乃至底层构造识别的基础。随着岩石成像技术的快速发展, 数字图像被广泛应用。在岩石图像处理中, 裂隙或者薄构造等重要结构相对于岩体的体积要小得多, 需要高精度与高分辨率的成像数据。如果对采集的全部数据进行存储处理, 不仅占用大量的存储空间, 而且面对如此庞大的数据存储量, 可视化及分析处理也将遇到前所未有的挑战(翟明国等, 2018)。另外, 图像特征提取过程中往往会丢失部分重要的特征, 且可能在后续去噪过程中完全丢失。因此, 需要理清岩石中的主要特征, 并采取有效的图像处理技术对岩石图像特征进行精确提取,从而在保留主要特征的前提下降低数据存储量。

近年来, 图像处理技术快速发展, 蚁群算法(Dorigo et al., 1999; Dorigo and Stüzle, 2003)、神经网络技术(Zhang et al., 2018)、遗传算法(Rodrigues et al.,2017)、免疫算法(Yu et al., 2016)等被广泛运用到计算机图像处理中。其中, 蚁群算法最早由 Alberto Colorni 在1991 年提出(Colorni and Dorigo et al.,1991), 它作为一种具有正反馈性、鲁棒性和并行性的高效率仿生算法, 适用于复杂的图像数据处理中。已有学者将蚁群算法运用到图像特征提取(Aydın and Uğur, 2011)、图像分割(Ma et al., 2009;Arnay et al., 2017)等复杂图像处理问题。但蚁群算法搜索时间长, 容易停滞而陷入局部最优, 且修复裂隙特征容易产生多余特征。为了降低数据复杂度并提高效率, Ester et al. (1996)提出了DBSCAN 聚类算法, 该方法能够自动分簇, 不需要人工输入簇数就可将一组庞大的数据分为多个簇。此外, 在图像识别中,Canny 边缘检测算子(Canny, 1986)被提出, 并用来准确识别图像的边界特征。以上这些方法已经广泛的运用于人脸图像识别(孙珊珊, 2016)、医学(Yoursefi et al., 2012)、遥感(赵芳等, 2020)和航空航天(陈侠等,2019)等领域, 在地学中, 蚁群算法主要用于断层识别(叶涛等, 2018), 但未有应用于岩石学的图像处理中。岩石中的裂隙所占的体积虽小, 但其对于岩石的宏观力学性质及渗流特征有着重要影响。比如裂隙是含油气储层主要的储集空间和运移通道, 准确提取岩石裂缝特征有助于评估储层的连通性和渗透性, 进而提高油气的采收率。本文将综合利用DBSCAN 聚类算法、蚂蚁算法和Canny 算子三种算法, 对岩石精细裂隙特征进行提取、修复和连接, 确保在图像特征提取过程中保留岩石精细结构的主要几何形状特征, 为下一步相关分析计算奠定基础。

下面简要介绍蚁群算法和Canny 算子的基本理论, 提出蚁群算法、DBSCAN 聚类算法和Canny 算子的混合算法的框架; 随后基于PANDAS 软件平台(Xing, 2017)对岩石图像特征提取, 并运用混合算法对提取特征进行修复, 对比单独使用蚁群算法和混合算法的应用效果; 最终, 获得较完整的裂隙特征数据, 并对计算结果进行讨论。

1 基于蚁群和Canny 边缘检测算子混合算法的图像处理方案

1.1 蚁群算法原理

20 世纪50 年代, 许多学者开始相继研究各种仿生算法以寻求最优化问题, 其中, 蚁群算法以独特的优势得到了许多学者的关注。在蚁群算法中, 最经典且与现实生活紧密相关的是旅行商问题(traveling salesmen problem, TSP)。它主要是描述如何找到最短路径、快速遍历所有城市。将常见的TSP 数学问题进行模型的建立, 可以更好理解蚁群算法。

一只蚂蚁随机从n座城市中的某一座城市出发,对所有的城市只访问一遍绝不重复, 直到遍历完所有城市蚂蚁将停止搜索。设n座城市的集合式中为C={C1,C2, …,Cn}, 两两城市之间的距离集合为L={li j|ci,cj⊂C}, 则G=(C,L)构成有向图可以直观地表现出蚂蚁遍历城市的场景, 将现实生活中的问题转化为对最短Hamilton 圈的求解。城市坐标为(xi,yi)和两个城市之间的距离满足下列公式:

设t时刻位于第i座城市上蚂蚁的数量为bi(t),则系统中蚂蚁的总数量为在t时刻从城市i到j的信息素大小为τij, 信息素增量为Δτij(t)。当一只蚂蚁从一个城市转移到另一个城市, 信息素挥发系数ε将随时间不断的局部更新,Δτij(t)=1/(n×Lnn),Lnn为最邻近算法求得蚂蚁遍历所有城市的距离。更新规则满足(Dorigo, 1997):

当所有蚂蚁都完成转移, 信息素挥发系数ρ将随时间不断的全局更新, Δτij(t)=1/L,L为蚂蚁经过所有城市的最短长度, 更新的规则满足:

在t时刻从城市i到j的启发函数为ηij(t)。并且启发函数与距离dij的关系为:

蚂蚁转移规则为:q为随机数,q0是0～1 之间的值, 当q≤q0时, 则从城市i转移到城市j:

当q>q0时, 则选择下一个城市的j的概率:

allowedk是蚂蚁可以选择下一个没有经过的城市的集合, 它与蚂蚁所经过的城市形成的禁忌表tabuk互为补集。一般q>q0时随机访问下一个城市,为了方便应用, 此处更改为访问allowedk的一个城市。期望启发式因子β决定了信息素在路径中的影响程度, 局部残留的信息素挥发系数ε决定蚂蚁选择其他路径的概率, 全局残留的信息素挥发系数ρ决定了全局搜索能力及收敛速度。

1.2 Canny 边缘检测算子

Canny 双阈值算子是一种图像处理技术中常用的利用双阈值和梯度方向的边缘检测算法。它先用高斯函数对初始图像进行平滑, 将平滑后的图像计算梯度大小和方向, 根据梯度差值来寻找图像边缘。

设二维高斯函数为:

将F(x,y)利用高斯函数进行卷积平滑处理后得到图像:

采用2×2 邻域内一阶偏导数有限差分计算I(x,y)的梯度幅值和梯度方向, 将点在x方向和y方向的一阶偏导数分别记为Kx(x,y),Ky(x,y):

之后对梯度幅值进行非极大值抑制, 若M(x,y)比θ(x,y)方向的两个相邻像素点的幅值小, 则归零。若像素点梯度的幅值M(x,y)大于或等于θ(x,y)方向的两个相邻像素点的梯度幅值, 则将该点判定为可能为边缘点。再对该点进行判定, 设置高、低2 个阈值, 大于高阈值的认为是强边缘点进行保留, 介于两者之间的认为是弱边缘点, 若采用8 连邻域像素检查周围有强边缘点, 则保留该点。

1.3 混合算法处理流程

为了将岩石图片特征提取缺失的细线状裂隙进行修复, 本文将会进行以下几个操作: 首先将岩石片状图片转化为含有色彩信息的像素数据, 并提取出主要的裂隙特征。利用DBSCAN 聚类算法将裂隙特征点划分为多个簇, 每个簇利用蚁群算法进行裂隙连接。然后, 将多余的裂隙利用Canny 算子得到图像求交集去除多余特征点。详细的步骤如下:

步骤1: 根据岩石片状图像像素点灰度均分为多个区域, 每个区域赋予一种颜色信息并做标记,提取出裂隙颜色信息。

步骤2: 根据细线状结构断裂处与周围点像素特征, 提取了此次需要的重要像素点, 降低数据量,提高算法的运行速率。

步骤3: 初始化蚁群参数β、ε、ρ, 并且根据城市坐标计算两座城市之间的距离, 对信息素进行初始化。设置迭代次数Tmax。

步骤4: 设置DBSCAN 聚类算法的邻域阈值Eps, 点数阈值Minpts。

步骤5:X1k,X2k表示两个数据点k维坐标, 表示两个点之间距离。当一个点与周围的点满足且周围点数大于等于Minpts时, 该点为核心点, 以核心点逐渐向外访问,分为簇C(i)并标记该点到访问集合v。直到未访问集合访问完为止。

步骤6: 将蚂蚁随机分布到n座城市上来进行搜寻, 并进入迭代次数循环,N=N+1。

步骤7: 蚂蚁数目也开始循环,k=k+1(0

步骤8: 蚂蚁根据概率选择满足的下一个目的地, 转移到该目标城市后将城市加入到禁忌表, 直到禁忌表满结束搜索。

步骤9: 如果k

步骤10: 根据公式(2)对信息素进行满足蒸发机制的全局更新, 并再一次进入到步骤6 开始下一次迭代, 直到达到最大迭代次数。最终找到全局最优路径。

步骤11: 设置灰度阈值提取裂隙图像, 设置Canny算子双阈值提取裂隙图像, 将两幅图与裂隙修复图像求解交集去除多余特征。

2 基于混合算法的岩石图像分析

2.1 岩石图像数据的选取

应用提出的混合算法, 以图1a 岩石图像为例进行分析。利用PANDAS 软件将图1a 所示的岩石不同特征提取出来, 包含矿物颗粒边界及其裂隙(为描述方便, 本文以下统称为裂隙)、岩石基质、岩石矿物, 分别如图1b～e 所示。由图1b 细线状为主的裂隙特征图像可以看出, 在特征提取的过程中, 虽然保留了主要部分及其特征, 但岩石的裂隙图像在提取过程中部分缺失, 如图1b 中红色框线标注所示。

2.2 蚁群参数的选取

较好的蚁群算法的参数可以使算法获得全局最优路径和较短的运行时间。为验证本次实验的有效性,应选取合适的蚁群参数。根据前人对bayg29, att48等常用数据库的研究可以判断出合适的蚁群参数(甘屹和李胜, 2011), 对于不同裂隙数, 为获得全局最优路径, 蚁群算法参数范围相同, 但具体的取值不同。根据全局最优路径大小, 选取最合适蚁群参数需要大量的实验, 没有规律和理论作为参考。以岩石裂隙图像1b 为例, 可通过大量的实验分析参数对算法的影响, 来选取适合的参数。

该数值实验算法主要涉及参数β、ε和ρ, 对其不同组合进行测试研究。随着蚂蚁数增加, 搜索能力增强, 但搜索速度呈指数增加。本文蚂蚁数设置为150, 迭代次数为200, 图像均匀分割为三份,Eps设置为15, 缺省参数设置为q0=0.9。为了防止信息素积累或挥发过快, 在公式(2)和(3)中,ε和ρ是介于0～1的参数。在公式(5)和(6)中,β越大计算量越大。本文选取0<β<10, 0<ε<1, 0<ρ<1 进行实验。β以1 开始步长2 递增、ε和ρ以0.1 开始步长0.2 递增, 调整蚁群参数, 判断对全局最优路径和运行时间的影响。

数值实验运行结果显示, 参数ρ固定时,β和ε两个值的变化影响全局最优路径和运行时间(图2)。ρ=0.1时, 7≤β≤9, 0.1≤ε≤0.3 时全局最优路径小(图2a)。ρ=0.1 时, 随着ε的增大, 算法运行时间不变, 运行时间随着β的增大逐渐增加, 当β≥3 时, 运行时间趋于稳定(图2b)。

而参数ε固定时,β和ρ两个的值变化影响全局最优路径和运行时间(图3)。ε=0.1 时, 7≤β≤9, 0.1≤ρ≤0.7 时全局最优路径小(图3a)。ε=0.1 时, 随着ρ的增大, 算法运行时间不变, 运行时间随着β的增大逐渐增加, 当β≥3 时, 运行时间趋于稳定(图3b)。

参数β固定时,ρ和ε两个的值变化也影响全局最优路径和运行时间(图4)。β=9 时, 0.1≤ρ≤0.7, 0.1≤ε≤0.5 时全局最优路径小(图4a)。β=9 时, 随着ρ和ε的增大, 算法运行时间不变(图4b)。

综合以上研究, 参数β、ε和ρ对全局最优路径都有一定的影响, 总的来说, 7≤β≤9, 0.1≤ε≤0.3, 0.1≤ρ≤0.7 范围内全局最优路径较小。参数ε和ρ对运行时间影响非常小,β对运行时间影响大, 但β≥3后运行时间趋于稳定。并且可以看出, 当全局最优路径较小时, 不可避免要消耗一定的时间。

为了防止不同裂隙之间相互干扰, 应选取较短的全局最优路径连接岩石裂隙, 并基于合适的蚁群算法参数进行实验。针对图1b 部分特征缺失和产生多余特征的情况进行研究, 根据大量实验得到图2～4, 本文基于β=9,ε=0.1,ρ=0.3 蚁群算法参数, 从岩石图像裂隙连接和修复两个方面, 对实验结果进行分析和讨论。

图1 图像特征提取Fig.1 Image feature extraction

2.3 岩石图像裂隙连接和修复

图2 ρ=0.1 时, β 和ε 的变化对全局最优路径的影响(a)和β 和ε 的变化对运行时间的影响(b)Fig.2 Th e influence of β and ε on the global optimal path (a) and the influence of β and ε on the running time (b)(ρ=0.1)

图3 ε=0.1 时, β 和ρ 的变化对全局最优路径的影响(a)及β 和ρ 的变化对运行时间的影响(b)Fig.3 Th e influence of β and ρ on the global optimal path (a) and the influence of β and ρ on the running time (b)(ε=0.1)

图4 β=9 时, ρ 和ε 的变化对全局最优路径的影响(a)及ρ 和ε 的变化对运行时间的影响(b)Fig.4 The influence of ρ and ε on the global optimal path (a) and the influence of ρ and ε on the running time (b)(β=9)

为了合理地判断实验效果, 本文通过定性和定量两种方式进行对比。定性即通过对比蚁群算法、混合算法连接和修复的图像与人工连接和修复的图像。定量则通过利用误差来验证蚁群算法和混合算法图像修复效果, 误差公式为:

其中, 期望图像是人工手动修复图像, 将图像修复像素数与理想修复像素数建立误差公式, 误差绝对值越小则修复效果越好。当误差为负且绝对值较大时, 裂隙大部分未修复; 当误差为正且绝对值较大时, 产生许多多余的特征。

本实验选取图1b 进行修复, 由于在实验过程中DBSCAN 算法的邻域阈值Eps对修复效果影响最大,将比较不同Eps对图像的影响(点数阈值(Minpts)设置为固定参数2 不变)。对比蚁群算法和混合算法实验效果差异, 实验结果如图5、6 所示。

图5 不同邻域阈值Eps 对岩石图形修复的影响(红色裂隙为修复产生)Fig.5 The influence of different neighborhood threshold Eps on rock image restoration

图6 基于岩石图像的蚁群算法和混合算法修复误差对比Fig.6 Comparison of r ock image repair errors of ant colony algorithm and hybrid algorithm methods

蚁群算法结果显示, 随着Eps的增加, 误差绝对值先变小后变大(图6)。Eps=5 时, 误差为负, 许多断裂特征未连接(图5b、6), 说明Eps过小会导致大部分断裂特征无法连接。Eps=8 时, 误差为0, 部分断裂无法连接, 且产生多余特征(图5c、6), 说明误差为0 只是因为期望图像和修复图像点数相同,但不是理想情况, 并且Eps较小会导致部分断裂无法连接, 且产生部分多余特征。Eps≥16 时, 误差为正且大于100%, 断裂特征基本连接但也有多余特征产生(图5d、e、6), 说明Eps适中或偏大时虽然断裂特征基本连接但产生许多多余特征。结果证明,蚁群算法虽然可以修复裂隙, 但产生的修复图像干扰过多, 与理想裂隙修复图像相似度低。

混合算法结果显示, 随着Eps的增加, 误差绝对值也先变小后变大。Eps=5 和Eps=8 时, 误差为负且绝对值较大, 许多断裂特征未连接(图5f、g, 6);说明使用混合算法时,Eps小会导致部分断裂特征无法连接, 且Eps越小连接的特征越少。Eps=16时, 误差为0, 断裂特征基本连接但也有多余特征产生(图5h、6), 说明误差为0 只是因为期望图像和修复图像点数相同, 但不是理想情况,Eps适中时连接效果好,产生极少的多余特征, 且混合算法中误差已经能够较准的反映断裂连接状况。Eps=30 时, 误差为正且小于50%, 断裂特征大部分连接但也有多余特征产生(图5i、6), 说明混合算法产生的多余特征少, 但Eps过大也会导致部分特征不能修复。结果证明, 混合算法修复裂隙且产生干扰少, 与理想裂隙图像相似度高。

结合图5、6 对比蚁群算法和混合算法可知, 混合算法降低了蚁群算法的误差, 且蚁群算法产生的误差越大, 混合算法能够降低误差的差值越大。这说明Canny 算子去除多余特征, 优化了蚁群算法连接的图像。同时, 相比于蚁群算法, 混合算法修复图像更接近理想裂隙图像。

为了验证混合算法的适用性, 本文选取了裂隙更复杂的岩石图片进行实验。由于蚁群参数的范围具有适用性, 在范围内, 这里依然选取β=9,ε=0.1,ρ=0.3 蚁群算法参数, 复杂岩石图像实验结果详见图7、8。

在复杂岩石裂隙图像中, 单独使用蚁群算法结果显示,Eps=6 时, 误差为0, 部分断裂特征无法连接且产生部分多余特征(图7c、8);Eps=10 时, 误差很大, 大部分特征修复, 同时产生部分多余特征,这说明Eps适中时, 蚁群算法效果也不理想(图7d、8)。Eps>10 时或Eps过大时, 误差很大, 断裂特征基本连接, 但产生很多多余特征(图7e、f, 8)。

在复杂岩石裂隙图像中, 使用混合算法时结果显示,Eps过小时, 误差为负, 大部分断裂特征无法连接(图7g、8,)。Eps=10 时, 误差为0, 断裂特征基本连接, 且只有少部分多余特征产生(图7h、8)。Eps过大, 误差较大, 特征大部分连接且混合算法产生的多余特征少, 但Eps过大也会导致部分特征不能修复(图7i、j, 8)。

由图7、8 可知, 混合算法误差小, 使用混合算法可以有效避免间断和去除多余的修复特征, 效果比单独使用蚁群算法更好。由于混合算法使用Canny 算子改进, 裂隙修复的结果得到了提升, 基本去除多余特征, 更接近理想的结果。总之, 实验结果证明, 混合算法同样适用于复杂的裂隙图像中,具有普适性。

3 结 论

在岩石图像修复过程中, 为了识别出裂隙等重要结构并降低数据量, 本文提出基于蚁群算法、DBSCAN 聚类算法和Canny 算子的混合图像特征修复算法, 来改善蚁群算法在岩石图像特征修复中存在的效率低、易产生多余特征的缺点, 较准确恢复岩石图像中缺失的裂隙特征。文中首先通过数值实验发现在7≤β≤9, 0.1≤ε≤0.3, 0.1≤ρ≤0.7 范围内全局最优路径较小;β对蚂蚁搜索时间有显著影响,大于等于3 时蚂蚁搜索时间增加。采用蚁群算法与混合算法得到的修复图像和期望图像进行对比发现,混合算法修复裂隙和去除多余特征的能力优于蚁群算法。实际应用结果证明, 混合算法也可以用于复杂岩石裂隙图像, 提高了图像修复的准确度, 验证了该算法的有效性。

图7 不同邻域阈值Eps 对复杂岩石图形修复的影响(红色裂隙为修复产生)Fig.7 The influence of different neighborhood threshold Eps on complex rock image restoration

图8 基于复杂岩石图像的蚁群算法和混合算法修复误差对比Fig.8 Comparison of complex rock image repair errors of ant colony algorithm and hybrid algorithm methods