基于分数阶微分及灰度和形状的节理裂隙提取

周 震 王卫星 王珊珊

（1.洛阳师范学院信息技术学院，河南 洛阳 471934；2.长安大学信息工程学院，陕西 西安 710064）

近年来，由于大规模的矿产资源开发，导致边坡滑坡及隧道坍塌等自然灾害的频发，严重影响了人们的生命财产安全，因此对于矿山的岩石工程研究具有重要意义。而岩石的稳定性对于矿山灾害具有很大的影响。为了对矿山的岩石边坡等安全进行防护与管理，通过检测矿山的岩石节理裂隙来判断岩石边坡的稳定性，进而可预防矿山的岩石工程灾害发生，此外，矿山中的爆破工程及铲装运输工程的优化设计和安全生产也都与节理裂隙有关［1-4］。

节理裂隙是一种很常见的构造地质现象，也是一项非常重要的岩石物理性质。它是由于岩石受力而出现的，但裂开面的两侧没有发生明显的（眼睛能看清楚的）位移。其存在于大多数岩体中，形成很多断面，这不仅影响岩体的连续性和完整性，而且可能会使岩体变形或破碎，导致岩体裂隙发生变化，从而可能导致矿山边坡坍塌提前发生［1-2］。因而在岩石工程应用领域中，岩石节理裂隙检测具有重要的现实意义。通常情况是通过测量和分析岩体节理裂隙的情况来判断其危害程度，据此预测岩层裂隙的稳定性，及时采取相应的防护措施，从而防患于未然［3-4］。随着图像处理和识别技术的日新月异以及计算机行业的高速发展，如何应用图像自动检测技术识别岩石中的各种节理裂隙，准确高效地获得清晰的节理裂隙信息，对于精确预测地震的发生、山体崩塌、矿山安全以及岩石工程顺利实施起着至关重要的作用［4］。

目前相对成熟的岩石节理裂隙提取的方法有基于支持向量机、基于统计模式识别和形态学等方法［5-7］。但由于岩石节理裂隙形状复杂无规则，且岩石表面粗糙给图像带来大量噪声，利用传统的提取方法很难达到预期效果，而且提取出的节理裂隙断断续续，连续性较差。为此，本文提出了一种基于节理裂隙走向和距离信息的岩石节理裂隙提取方法。首先对预处理的岩石图像进行初步分割，在提取出岩石节理裂隙的同时尽可能抑制噪声。其次对初步分割后的节理裂隙进行基于中轴变换的骨架提取。最后基于裂隙线段的走向和距离信息进行间隙缝合。

1 基于分数阶微分的岩石节理裂隙图像预处理

由于岩体表面粗糙，存在大量噪声干扰，所以在分割之前，需要对岩石节理裂隙图像进行预处理［8-9］。本文主要采用分数阶微分进行图像增强［10-11］，使节理裂隙图像边缘更加突出、纹理细节更加清晰，同时去除大量噪声，非线性地保留平滑区域的低频信息。不过，最优的分数阶微分阶次一般通过人工实验指定，耗时且不智能，故本文采用一种基于G-L理论、图像梯度和视觉特征的自适应分数阶微分算法［11-12］。

该方法可以自适应地调整不同灰度图像和梯度特征的算子掩模系数，从而得到相应的掩模模板，依次处理每个像素，并以灰度值、信息熵、平均梯度作为图像纹理评价参数，对其进行定量分析和实验验证。

关于分数阶微分，有多种定义。与其他定义相比，G-L定义在数值计算方面有着明显的优势：

为了得到8个对称方向上的分数阶微分，并使得其掩模具有旋转不变性，分别取0，π 4，π 2，3π 4及其各自反方向上的分数阶微分掩模。同时为了避免模糊那些灰度值变化明显的区域，降低计算复杂度，本文采用3×3模板来处理图像f。

因此，自适应分数阶微分可以使用不同的掩模算子对每个像素进行滤波操作，不仅增强了灰度变化明显区域中的像素信息，同时保留了平滑区域的信息，从而实现很好的图像增强效果。

2 基于灰度信息的节理裂隙提取—图像初步分割

图像分割是指把图像分成各具特性的区域并提取出感兴趣目标的技术和过程，其算法主要包括基于阈值的方法、基于边缘的方法和基于区域的方法。由于实现过程简单直观，基于阈值的方法在图像分割中起着重要的作用，其中基于灰度直方图的多阈值法和最大类间方差法由于其良好的分割结果［3，5］，计算简单且应用范围广，被广泛应用于图像处理过程中。

本文首先对预处理后的图像进行Sobel边界扫描，然后进行采用基于二维直方图的Otsu阈值分割算法，能够充分利用图像的边缘信息和细节特征，使分割结果更加准确，抗噪性好且简单快速。其基本原理是：假定二维直方图的横坐标表示中心像素的灰度级c(x,y)，纵坐标表示像素邻域的灰度梯度，即像素灰度级c(x,y)与其邻域平均灰度级d(x,y)差分的绝对值 |c(x,y)−d(x,y)|。如果灰度级e和图像中每个像素的邻域灰度梯度f组成一个二元组(e,f)，则相应的联合概率密度为

被动拉伸训练对大鼠腓肠肌HGF、MGF基因表达的影响…………………………………………吴国梁，李 娜（25）

式中，M×N为图像的大小；L为灰度级；cef表示二元组(e,f)出现的次数。

根据以上二维直方图的定义，任意阈值向量(g,h)都可以将二维直方图划分成四个区域。由于梯度值反映了邻域内的灰度差，所以小邻域的梯度幅值对应的灰度级表明该灰度级较多出现在目标或背景中，少量像素在边缘处，而大邻域的梯度幅值则说明该灰度级大多出现在边缘上，极少存在于目标或背景中。

在二维直方图的基础上，假设在二维直方图中，用I0表示目标类，用I1表示背景类，如果P0和P1分别为目标类和背景类的概率，则P0、P1可表示如下：

由式（7）、式（8）可得到目标类和背景类的均值向量μ0和μ1分别为

则总体的均值向量μz可表示如下：

3 节理裂隙的跟踪

3.1 骨架提取

骨架是一个目标物的重要拓扑信息，如何进行骨架提取是图像处理研究的热点问题之一。骨架提取是对二值图像的一种细化，具有重要的应用价值。常用的骨架提取算法是基于形态学的骨架提取算法［13-16］，主要分为两种：击中击不中变换和中轴变换。相对现有的其他算法，基于形态学提取的骨架既满足了连通性又具有稳定性，相比之下，基于中轴变换提取的骨架［7］连通性更好，因此本文采用中轴变换进行骨架提取，取得了很好的效果。

中轴就是精确定义的骨架，是所有与物体在两个或更多非邻接边界点处相切的圆心的轨迹。中轴变换是寻找满足以下条件的所有点及对应参数：以该点为中心存在一个包含于物体内的且与物体边界相切于两点的圆盘，该点的参数为相应圆盘的半径。因此中轴变换是一个得到物体骨架和骨架上每点到物体边界的最短距离的过程。确定图像A的骨架S(A)应满足两个条件：在求S(A)过程中，A应有规律地缩小；在A逐步缩小的过程中应使其连通性质保持不变。对骨架的定义如下：

式中，S(A)表示图像A的骨架；B是一个结构元素；(AΘkB)表示对A连续进行k次腐蚀；第K次是A被腐蚀为空集合前进行的最后一次迭代，即

为了得到完整的目标骨架，首先必须从骨架起点开始确定重要的拓扑节点，然后从每个拓扑节点开始调用提取算法，得到一个骨架，直至到达源点或者已经提取过的骨架线为止，以防出现重叠的路径。如果一个目标物可以表示为图形，那么它的重要拓扑节点就很容易被识别出来，因此，可以通过将目标物转换成图形的方式来区分突出部分的拓扑节点。假设图形生成过程由参数k控制，采用下面的方法能够自动提取出拓扑节点。

首先，用一个中速波从边界到中心循序渐进地计算最短距离域L(x)；接着通过中速传播一个演变的峰面k-面来自动选择一个骨架点作为源点Po，表面峰的运动是由短距离方程确定的，并且解为一个新的距离域L1(x)。峰面的速度由下式给定：

根据上面的方法已经得到了拓扑节点，那如何提取单条骨架线成为关键。为了提取连接两个骨架点R和S之间的骨架线，初始化骨架点R的传播时间为0，然后选择高速模式到达骨架点S；最后在ΔT内从S到R进行回溯追踪，提取过程可简化为求下面常微分方程的解：

式中，Y(t)用来描述骨架线。

3.2 基于节理裂隙走向和间距的间隙缝合

在上面的基础上，可以得到岩石节理裂隙的基本信息和大致走向，但同时存在大量断点间隙，导致节理裂隙断续不完整，因此需要进行间隙缝合。常用的方法为基于形态学膨胀腐蚀的操作方法，该方法只是不断缩小断点间的距离从而将缝隙缝合［16-18］，对复杂的图像不能正确的连接。对此，本文提出了一种基于裂隙间隙距离和走向的缝合方法，具体步骤如下。

（1）首先去除骨架提取后的毛刺即小区域，依据提取骨架线段的宽度或线段所在的行列坐标位置去除小的噪声区域。

（2）去除噪声后的断点连接：①从图像的最左边开始进行扫描［16-17］，找到断裂线段的端点，即Ai(x,y)(i=1,2,3...)，在Ai的8邻域像素中寻找端点所在的边缘线，沿着边缘线回退n（2≤n≤6）个像素，找到像素Bi(x,y)，在该断点Ai与像素Bi之间做一条直线记为ALi；②按距离在断点的8邻域进行搜索，按给定的预定距离长度范围找到邻近的几个断点Ci(x,y)；③在各个断点Ci中沿着断点所在的边缘线回退n个像素，找到像素Di(x,y)，把该断点Ci与像素Bi之间做一条直线CLi；④计算ALi和CLi之间的夹角，按实际图像给定角度θ，判断要连接的断点，用直线进行连接。

（3）如果断点没连接好且有多余的噪声，则返回（1），否则返回（2）继续连接。当所有断点连接好，则结束。

具体间隙缝合程序如下：

4 实验结果与分析

实验选取了有代表性的岩石节理裂隙图像如图1（a）和图2（a）所示。图1（a）为较简单的节理裂隙图像，本文先用分数阶微分对岩石图像进行增强去噪处理，保留图像的高频信息，同时去除低频噪声；然后用Sobel边界检测算子进行边界检测，再对梯度图像进行阈值分割，检测出大致的节理裂隙区域，在灰度图像中，对该区域进行阈值处理，将岩石节理裂隙大概提取出来如图1（b）所示。为了精确地获取岩石的节理裂隙，利用基于形态学的中轴变换对其进行骨架提取，结果如图1（c）所示。可见骨架提取后的图像仍存在诸多干扰噪声，直接影响了裂隙的连接，为此本文提出一种新的裂隙连接方法，能很好地进行裂隙连接。首先依据线段的宽度去除裂隙周边的噪声，再依据噪声所在的行列坐标进行二次去除，最后依据裂隙线段的方向和角度对岩石节理裂隙进行连接，本图选取连接距离范围为1～100像素，结果如图1（d）所示。由图可以看出该方法能很好地提取出岩石节理裂隙。

图2（a）为较复杂的岩石节理裂隙图像，图2（a）～图2（d）为其裂隙提取过程，本图选取连接距离范围为1～100，很明显可以看出，对于比较复杂的节理裂隙，本文方法也能达到很好的提取效果。

为了进一步比较分析，将本文的算法与Canny边缘检测、图论最小生成树、聚类分析及模糊聚类算法进行比较，结果如图1和2所示。由图1（e）～图1（h）可以看出，对于较简单的岩石节理裂隙，这4种传统算法也能较好地提取出节理裂隙的边缘信息，但由于岩石表面存在着大量的噪声，导致结果图像也存在许多噪声。但本文提出的算法不仅能很好地去除噪声，而且能提取出完整的裂隙信息，有效地进行间隙缝合，提取效果较好。从图2（e）～ 图2（h）可以看出，对于较复杂的岩石节理裂隙图像，对比的4种传统算法就很难提取出节理裂隙，而本文的算法依旧能够很好地提取出节理裂隙，充分显示出了本文算法的优越性。

5 结 论

由于岩石节理裂隙形状复杂，且岩石表面粗糙给图像带来大量噪声，利用传统的提取方法很难达到预期效果，而且提取出的裂隙断断续续，连续性较差。本文提出了一种针对岩石节理裂隙的提取方法。为了去除噪声和增强节理裂隙边界信息，首先对其进行分数阶微分处理；为了快速图像分割，然后结合简单的梯度图像和灰度图像信息初步提取节理裂隙；再对分割后的图像进行骨架提取，根据裂隙间的距离和走向信息进行裂隙间隙的缝合，本方法在连接的同时能很好地去除噪声，连接效果较好。该算法简单、快速，对岩石矿山的稳定性检测具有实际的指导意义和应用价值。