基于改进的Hessian 和Live-wire 算法的岩石节理裂隙检测

王 艳 李 晗 陈佳悦 陈卫卫 王梦菲 闫 迪 李宏霞 王卫星

(1.西安航空职业技术学院人工智能学院,陕西 西安 710089;2.辽宁省交通高等专科学校汽车工程系,辽宁 沈阳 110122;3.伦敦大学学院计算机科学系,伦敦 WC1E 6BT;4.西安航空职业技术学院自动化工程学院,陕西 西安 710089;5.长安大学信息工程学院,陕西 西安 710064)

在各种岩石工程应用中,岩石节理裂隙会显著影响岩体的物理性质[1-3]。地质科学中已经开发并采用了许多测量方法。同时,岩石裂隙对矿山工程中的岩石爆破质量的影响也很大。同时,对于公路、铁路、隧道和大坝,精确的岩石节理裂隙数据对于评估岩石边坡的稳定性至关重要。与此相关的研究有很多,如在混凝土和水泥工程中,裂缝将主导设施的破坏过程。为了解决由裂缝引起的工程问题,许多研究人员和工程师研究并开发了用于测量这些裂缝的各种设备和数学算法。

工程地质学中有一套完整的理论和方法来测量和描述岩体中各种节理裂隙的方向、密度和其他特征[1]。早些年,1976 年Priest 和Hudson 引入了岩石不连续间距[2];然后,在1981 年,他们测量了一个方向上的裂隙间距[3];之后,进行了更精确的节理裂隙测量[4]。近年来,在野外测量、岩心钻探、钻孔电视和摄像机测量中获得岩石节理裂隙的1D 和2D 的测量和分析可以用图像技术进行。一个典型的例子:岩石质量指标(RQD)被定义为由≥0.1 m 的全长组成的钻孔岩芯或扫描线的百分比,并已被广泛用作岩体岩石节理裂隙发育程度的工程测量[5-6]。

通过从岩石钻孔中获取三维RQD 信息,研究人员使用蒙特卡洛模拟方法创建了一个假设的三维裂缝网络。该方法允许对所获得的数据进行测量和分析。以验证该方法,将计算的RQD值与实际RQD值进行比较,得出了令人满意的拟合[5],这可能与岩石表面而非钻孔的3D 裂缝特征有关。建立了几个模型中的网络可信渲染,用于节理裂隙网络表征。然后,应用最佳拟合Ferret 方法自动确定节理裂隙面积,然后计算裂隙长宽、裂隙线端点及交叉点、颗粒分布等参数来综合评价节理裂隙的分布和发育状况。最终得到网络的孔隙率、密度、连接性和复杂性。然而,了解岩石裂隙的三维分布要复杂得多。因此,大多数研究人员都进行了三维建模和模拟,以避免直接在岩体中进行测量。

自20 世纪80 年代以来,图像分析技术得到了迅速发展,并研究了许多用于岩石节理裂隙测量和分析的算法[6]。在早期对岩石节理裂隙的图像处理和分析的研究中,主要工作是对单条裂隙的检测和识别(如孔径检测)[7]。近年来,研究已扩展到多裂隙处理和分析,包括岩石节理裂隙方位测量、空间分析等[8]。以前的研究是或已经是基于与不同岩石工程相关的岩石特征而进行的。通常,在岩石节理裂隙图像处理和分析之前,图像是通过不同的传感器而获取的,所以图像的类型和设备有关。然后根据图像的类型和特征,利用计算机软件对图像进行预处理、分割和裂缝测量与分析。最后,将裂隙分析数据应用于估算应力、应变、位移、渗流计算等。在过去的10 年里,对于3D 岩石节理裂隙的测量,激光扫描技术已在现场使用,主要用于整个岩石表面或单个大裂隙的测量和分析[9-10]。而在实验室中,对于小型的3D 节理裂隙岩石样本,可以更详细地对节理裂隙参数进行测量和分析,一般都应用了CT 技术[10]。无论应用何种硬件,都可以使用1D 和2D 裂缝信息来估计3D 岩石节理裂隙情况。通常,在涉及单个或多个岩石节理裂隙的图像处理中,首先将图像增强算法应用于岩石节理裂隙图像[11-12],然后根据对象(节理裂隙)的特征研究提取算法[13-14]。在图像分析中,岩石节理裂隙的标准参数包括2D 孔隙度、孔径、节理裂隙密度和方向等。

近5 年来,语义分割在路面裂缝检测中已经成为了一个研究热点。其分割网络是基于深卷积网络和神经网络的模型,语义分割最基本的任务是对路面图像中的不同类别图像像素点进行归类,分类的目的以区分图像中不同的目标对象(如:裂缝和非裂缝)[15-16]。进一步的发展在 U-net 等模型结构中引进了扩张卷积的运算,可在不降维数的情况下增大局部感受野信息和收集多尺度信息。在 Deeplabv1 等网络中,能够更进一步地增加展开卷积之比例,并且在相邻的像素关系推理中加入条件随机场模型。当前,Deeplab 系列网络仍然在不断地发展中,其骨干网络也得到了应用[17]。但天然的岩石节理裂隙比人造的混凝土裂缝要复杂得多,所以深度学习等语义分割只能用在特殊场景中,很难有普遍的适应性。

综上所述,岩石节理裂隙的测量和分析在岩石工程中非常重要,图像测量手段有一定的优越性,但由于节理裂隙图像千变万化及噪声太多,导致二维和三维的节理裂隙测量和分析都会遇到很多问题。基于此,本文研究了一种复杂岩石节理裂隙边界线的检测和跟踪算法,研究内容主要包括三个步骤:① 根据岩石节理裂隙图像的特点,提出一种新的图像预处理方法来平滑和增强节理裂隙信息,研究内容主要包括图像缩小,Geodesic Shadow Removal 和基于Hessian 矩阵的图像增强;② 根据Live-wire Contour 的算法思想,研究了一种新的图像分割算法,用于提取节理裂隙边界线上的特征点;③ 在点距离和角度差都满足基本要求的前提下,综合获取的特征点方向、亚像素精度定位及强度进行临近特征点连接。

1 图像预处理

在通常情况下,岩石表面比较粗糙,噪声较多,直接进行目标边缘线的跟踪比较困难,所以要进行图像预处理。预处理主要包括3 个方面:图像缩小,如果图像尺寸太大,为了减少计算量,要进行图像缩小,这里要进行的是有选择地缩小,图像缩小本着保持边缘信息及减少噪声的原则进行;图像平滑,主要是尽可能多地去除图像中的各种噪声,如:黑白点、阴影、三维体等;图像增强,也就是节理裂隙边缘的增强,增强边缘的原则是尽可能同时抑制噪声的产生和增强。3个方面的具体叙述如下。

1.1 图像缩小

为了减少图像存储量及后续算法运行时间,我们需要设计一种算法,既能缩小图像尺寸又不丢失节理裂隙信息,这种算法与多尺度空间理论有关。我们将在原始图像中,以每个检测点为中心,取3×3 模板的3 个最小灰度值的平均灰度值作为缩小图像对应点的值,如图1 所示,图1(b)～(e)是图4(a)4 种不同算法的缩小结果。

1.2 图像平滑

进行Geodesic Shadow Removal 平滑[18],该算法主要包括5 个部分:数学形态学“闭合”,高斯平滑,分水岭,直方图均衡化,对比度增强。

利用形态学的“闭合”函数来消除由于裂隙具有与阴影区域相似的强度而产生的裂隙。在阴影区域估计中,它可以防止裂隙被计入阴影。如果应用高斯平滑滤波器去除噪声,可进行下面三步骤。

(1)利用分水岭法将图像分割成一系列具有不同平均光照的测地线区域,{Ri|i=1,…,N,…,M},以S={Ri|i=1,2,…,N}区域为阴影区域,以B={Ri|i=N+1,N+2,…,M}区域为背景区域。这里,N=M。

(2)选择亮度补偿的最佳参考区域。我们把一个测地区域设为Hi=max(Ii)-min(Ii) 和把一个测地区域的灰度值的跨度设为H。它是指与半影区域相对应的测地区域的H值通常大于背景区域。它通过公式将图像大致分为阴影区和背景区,其中HS和HB分别表示低值和高值的类别。我们可以通过公式{HS,HB} =OTSU({Hi|i=1,2,…,M}) 选择最佳参考区域其中,Ai=ave(Ri) 和Di=dev(Ri) 分别表示图像中像素的平均值和标准差。

(3)将像素(i,j)∈S分配为新值I'i,j=qIi,j+e。其中,Ii,j表示像素(i,j)处的灰度值,I'i,j表示亮度补偿的结果,和DB分别表示S和B的标准差,和分别表示S和B的平均灰度值。这种算法可以增强图像的对比度。

1.3 节理裂隙增强

Hessian 矩阵I在点p(x,y)定义为

相应的节理裂隙定义如下:

式中,td(p,θ;σ)=gxxcos2θ+gyysin2θ+gxysin2θ.

对于本文提出的一套方向滤波器—前向过滤器tf(p;σ,ψ1) 和后向过滤器tb(p;σ,ψ2) ,其定义如下:

在上述公式中,采用角度参数ψ1、ψ2检测相邻像素之间的裂纹迹象;常数d是偏移量参数,尽可能设置为适当的值,若其值较低,其贡献就小,将引起图像分割错误,以至于将不是节理裂隙的像素错误地归为节理裂隙像素。

图像通过定向滤波器的响应分别为Tf(p;σ,ψ1)=tf(p;σ,ψ1)×I(p)和Tb(p;σ,ψ2)=tb(p;σ,ψ2)×I(p)。增强节理裂隙的整体响应如式(5)所示:

然后在多个方向上搜索最大响应作为结果的输出,如图2 所示。

2 岩石节理裂隙边界检测基于Live-wire 算法的改进

Live wire 是在1992 年由Barrett 和Mortensen 利用交互式边界提取的一种图像分割方法[19]。其基本思想是利用动态规划来获取图像中给定的2 点间最短路径,然后构造合理的代价函数来检测目标的边界(Live-wire contour),但该算法不是全自动的。本文为了增强其自动化程度,首先对掩模(template)及中心线(skeleton)进行图像的去噪声处理及平滑处理得到了新掩模(template_pro)及新的中心线(skeleton_pro);然后把Live-wire 的代价函数写成:

式中,fZ(q) 是Laplace 的过零点;fG(q) 是梯度的大小值;fD(p,q) 是梯度的方向;ωZ、ωG、ωD是相应权重。

式中,IL(q)≠0 是图像及通过LOG 边界检测算子的卷积后在图像中的q点处的新值。

式中,IG(q)是原始图像与梯度算子卷积运算后在点q处的值。

定义D(p)= (Gy(p),-Gx(p))为点p处的梯度方向。

其中,dp(p,q)=D(p)·L(p,q),dq(p,q)=L(p,q)·D(q)。

L(p,q)定义为

因为Canny 边界检测器可以在降低噪声和精准定位目标边界点之间选择一种平衡的方案,故本文的代价函数是基于Canny 检测器而构造的,从而可以降低噪声的干扰和提高定位的精度。

节理裂隙边界种子点搜索原理如图3 所示,3 点S、P、E为1 条中心线上的种子点;而T、Q、F和L、R、G分别为节理裂隙边界种子点左侧和右侧的对应种子点。

图3 边界种子点的搜索原理Fig.3 Principle of seed point searching on object boundary

在图3 中可将边界线的权重定义为

根据上述,那么对一幅无方向的加权G(V,E)连通图,可定义其最短路径为G(V,E),设其最短路径图为Gpath(Vpath,Epath)。首先从G(V,E)中,提取出源点u0,然后放进Vpath,这样,我们有Vpath= {u0},Epath=Ø;然后从任一的u∈Vpath,v∈V-Vpath路径当中,我们选择权重weight(u,v) 为最小,把顶点记为weight(u,v)为最小的,其顶点记为u*,v*,而把边记为＜u*,v*＞,Vpath←Vpath∪{v*},Epath←Epath∪ {＜u*,v*＞},修正源点u0到V-Vpath各个顶点的权值weight(u0,vi),修正方法如下:

最终,一直重复上述过程,当Vpath=V时,停止。这样通过Epath中的那条路径就被确定为最短的路径。具体步骤及检测阶段效果可以图4 为例。

图4 节理裂隙的边界检测Fig.4 Fracture edge tracing

3 实验与分析

上面详细地介绍了基于改进的Live-wire Contour算法的边缘检测算法及图像预处理算法[20]。为了表述清楚,这里选择一幅典型的有明显裂隙宽度的图像进行测试,为了体现新算法的优越性,我们选择了常用的11 种算法进行对比试验。如图5 所示,原始的岩石节理裂隙图像是一幅光照分布不均的图像,根据试验数据统计,裂隙边缘提取过程中用到的参数设置为:根据同类100 幅图像的训练,得出Frangi 滤波器对图像增强的过程中估计的σmax=σmin= 3,二值化过程中选取的T=16,节理裂隙筛选过程中所选取的参数S=10,r=3。

图5 节理裂隙检测算法的比较Fig.5 Comparison of joint crack detection algorithm

图5(b) ～图5(d)为3 种不同的阈值检测结果[7],由于光照不均,大津阈值使得近50%的面积为背景(黑色,包括裂隙),但实际上,裂缝的面积很少,估计约有10%;自动4 个阈值分割结果要好一些,如果考虑黑色为裂隙,大部分裂隙被检测到,但有些裂隙没有完全被检测出来,而非裂隙区域也有黑色的像素;动态阈值分割得要比前面的阈值分割稍微好一些,但也存在着类似的问题,并且产生了许多小目标,为后续的图像处理造成困难。图5(e)～图5(h)是4种边界扫描结果[21]。一阶导数的Sobel 和二阶导数的Laplacian 的边界检测结果是梯度图像,图像中的假边缘点太多,而真的边缘点又没有连成线,这会给后续的二值阈值处理造成很大的难度;而Canny-1 是低阈值的边界检测,其结果给出了太多的噪声边界;当设置高阈值时,Canny-2 产生了欠分割的问题。图5(i)～图5(l)是4 种较为复杂和精确的图像分割算法结果。基于图论的MST (最小生成树)[22]显然有着欠分割的问题,而聚类分析给出的结果要好一些,但仍有欠分割的问题,在FCM(模糊聚类)[7]的结果图中,棕黄色的区域包括了大部分裂隙,但有过分割和欠分割问题,合并与分离(Merging & Splitting)基本上能够检测出大部分裂隙的区域,但仍有少量的欠分割现象,还有就是设置参数时,要做多次试验。

总之对于这种噪声较多且光照不均的岩石节理裂隙图像,无论是传统的图像分割算法(如基于相似性的阈值法及基于不连续性的边界扫描法) 还是更复杂的图像分割算法 (如基于图论及基于聚类分析的算法) 都很难直接凑效。

如果用本文研究的预处理和边缘检测算法对同样的图像进行处理,则可以得到比较好的结果。当对图5(a)的图像进行预处理后,可以得到图6(a)中的图像,光照度变得均匀,噪声会减少。在此基础上进行裂隙边缘的特征点检测及连接,得到初步的裂隙边缘检测结果,如图6(b)所示,最后进行裂隙扩展和过滤,得到裂隙充填的图像,如图6(c)所示。

图6 新算法对图5(a)的检测过程和结果Fig.6 Fracture image (Fig.5(a)) detection procedure and results by new algorithm

4 结 论

本文研究了一种节理裂隙边缘线跟踪的新算法。在跟踪之前,若图像太大,进行有选择的缩小,图像缩小本着保持边缘信息及减少噪声的原则进行。然后基于改进Geodesic Shadow Removal 和Hessian 矩阵等算法将具有粗糙面的岩石节理裂隙图像进行平滑和增强。根据图像预处理的结果进行节理裂隙边缘跟踪。跟踪分成2 个步骤:① 基于Live-wire 的思想检测出节理裂隙边缘上的特征点;② 根据特征点之间的距离及夹角来连接成长短不一的节理裂隙线段,然后再连接线段以形成完整的节理裂隙边缘线。为了验证该算法的优缺点,选择了200 幅图像进行试验,均得到了较为满意的结果。为了分析该算法的有效性,选择了11 种传统的算法和更精细和复杂的算法对200 幅图像进行处理分析。结果表明:对于有明显节理裂隙宽度且具有较多噪声的图像,新算法有着较强的适应性,能够准确地提取节理裂隙边缘线。今后进一步的研究是选取更多特性不同的图像进行试验,以改进算法来扩展其适应范围。