基于Canny优化算法物理模拟实验覆岩裂隙场自动识别研究

王顺阳，贺易航，韩 博

(1.河南理工大学 测绘与国土信息工程学院，河南 焦作 454000； 2.河南省岩石矿物测试中心，河南 郑州 450014)

大规模煤炭资源开采破坏了围岩的应力平衡状态，导致上覆岩层发生移动变形进而形成裂隙，影响煤矿的安全生产[1-2]。因此，研究采动覆岩裂隙场发育规律对于保障矿井安全生产，促进我国煤炭工业健康稳步发展具有现实意义。

通过阅读国内外参考文献发现，针对沥青路面裂隙、岩土体裂隙等开展识别研究的比较多[3-5]。针对相似材料模拟实验开展的研究主要可分为2个方面：①覆岩裂隙发育特征和时空演化规律研究。Wang等[6]采用相似材料模拟实验和颗粒流数值模拟实验，对孔隙变化和裂隙发育规律进行了研究，并与现场原位测试的结果进行了对比。Behrooz等[7-8]针对上组煤和下组煤存在部分重叠的2个长壁开采工作面采动覆岩形变特征进行了分析，采用包括近景摄影测量和激光扫描仪在内的监测手段对覆岩位移量进行了监测，并将监测结果与有限元法所得结果进行了对比。郭昌贵等[9]利用FLAC软件建立煤层开采数值模拟模型，研究了不同高度的导水裂隙带的发育规律。柴敬等[10]采用光纤光栅传感器对采场覆岩的内部应变进行了监测，以此来研究关键层的初次破断和周期破断特征。②覆岩导水裂隙带发育高度预测研究。Wang[11]基于关键层沉降的分析和评价提出了采动覆岩裂隙率的理论分布模型，并采用UDEC数值模拟方法证实了裂隙的“拱形”发育特征。此外，有关学者基于图像分割对采矿过程中裂缝等问题进行研究，张昊[12]针对模糊C均值图像分割方法存在大量噪声干扰问题提出了一种融合空间噪声信息的改进模糊C均值图像分割方法，结果表明新的分割方法能够实现图像的精准分割。徐平等[13]针对采用阈值分割算法分割图像时，背景存在干扰的问题，提出利用分割算法对图像进行初始分割和二次分割，研究结果表明分割出的图像边界清晰能够较好地抑制噪声。邹盛等[14]针对现阶段煤矿带式输送机机辊缺乏有效的监测，传统人工巡检方法效率低、自动化程度低等问题，提出基于Hough变换的带式输送机托辊图像分割方法，实验结果表明该算法可以实现不同位置托辊的自动分割，能够为煤矿安全开采带来一定的技术基础。邵佳园等[15]针对油基泥浆裂缝识别问题，选用具有各向异性的曲波变换对超声波图像进行了图像分割，此方法能够获取图像边缘信息，对图像分割有一定的作用。综上，前人从不同行业背景和专业视角研究了采动覆岩裂隙的时空演化规律，构建了形式各异、功能多样的数学模型，丰富了裂隙/裂隙场的识别理论。但针对相似材料模拟实验在开挖过程中形成的裂隙场，开展自动识别研究的还不多。

基于此，本文以某矿15235工作面为工程研究背景，开展物理模拟实验，反演煤层开挖过程中上覆岩层的运移过程，就采动过程中发育的裂隙场开展识别和定量分析研究。

1 Canny算法裂隙识别的不足及其改进

1.1 基于Canny算法的裂隙识别

Canny边缘检测算法是一种数字图像裂隙识别方法，它是基于变分原理导出的一种用高斯模板导数逼近的最优算子。裂隙识别包括4个步骤：高斯滤波、梯度及梯度方向计算、非极大值抑制和边缘连接。Canny算法首要步骤就是利用高斯函数对图像进行滤波，图像上任意点(x，y)二维高斯函数可表达为：

(1)

式中，σ为高斯滤波器参数，控制图像滤波的平滑程度。

关于梯度及梯度方向计算，Canny算法利用一阶微分来计算图像f各点处的梯度幅值及梯度方向，进而得到相应的梯度幅值图像C和梯度方向。任意点(i，j)处沿x和y两个方向的偏导数Cx(i，j)和Cy(i，j)可用式(2)和式(3)来表达：

Cx(i，j)=[f(i，j+1)-f(i，j)+f(i+1，j+1)-f(i+1，j)]/2

(2)

Cy(i，j)=[f(i，j)-f(i+1，j)+f(i+1，j+1)]

(3)

根据式(2)和式(3)可解算出任意点(i，j)处的梯度及其方向：

(4)

(5)

将当前像素的梯度强度与沿正负梯度方向上的2个像素的梯度强度进行比较，如果前者与后者相比较大，则将当前像素保留为边缘点，否则该像素被抑制；然后通过阈值来识别边缘像素点，最后对识别出的边缘像素点进行连接，从而达到对裂隙边缘的自动提取。

1.2 Canny算法裂隙识别的不足及优化

Canny算法在裂隙识别方面具有一定的优势，但需要人为设定阈值，一方面费时费力，导致计算资源浪费，另一方面不可避免地影响裂隙识别的准确度。为了克服这种缺陷，论文提出一种优化算法来自动选取最优阈值，实现图像前景信息(裂隙)与背景的分割。

现假设一幅图像的像素大小为M×N，其灰度级为L，图像上某点(x，y)处的灰度值为f(x，y)，由于f(x，y)∈[0，L-1]，如果图像中灰度级 包含的像素个数为fi，则该级灰度出现的概率P(i)可表达为式(6)：

(6)

将图像中的像素按灰度级用阈值t(灰度值)划分为2类，即背景C0和目标C1，背景C0的灰度级为0～t-1，目标C1的灰度级为t～L-1，则背景C0和目标C1对应像素的灰度值满足条件：{f(x，y)

背景C0和目标C1出现的概率ω0和ω1可分别表达为式(7)和式(8)：

(7)

(8)

背景C0和目标C1的平均灰度值分别为：

(9)

(10)

图像的总平均灰度值为：

(11)

为了确定图像分割的最优阈值，平衡背景方差和目标方差值，构建了顾及权重系数的类间方差函数δ2(t)：

δ2(t)=ω0[μ-μ0(t)]2+aω1[μ-μ1(t)]2

(12)

式中，a为权重系数，表征图像目标像素个数与整幅图像像素个数之比。

将式(9)、式(10)和式(11)代入式(12)即可计算不同t值下的类间方差，当δ2(t)最大时的t值即为所求最优阈值，即:

(13)

2 工程验证

2.1 相似材料模型铺设

15235工作面开采2号煤层，平均开采深度约780 m，工作面走向长度870 m，倾向长度140 m，缓倾斜煤层，开采厚度5.5 m，走向长壁后退式放顶煤开采。实验采用CM400/22平面应力试验台，模型架尺寸4 000 mm×220 mm×1 500 mm，模型相似比为1/100，模型重度之比为0.6，时间相似准数为24。相似模型中骨料为河砂、石膏、腻子粉，缓凝剂为硼砂，模拟煤层上方一部分岩层，其余部分加等效载荷。模型铺设高度133.5 cm，上覆岩层厚度120 cm，下伏地层厚度8 cm。根据工作面2-2钻孔柱状图，相似材料模型共铺设41层。为了能采用数字散斑法获取模型表面任意点的形变量，需要制作人工纹理(腻子衬底)，人工纹理用毛笔蘸墨汁点涂完成，以增强形变的量测精度。

2.2 裂隙识别过程

下面以峰峰矿区15235工作面为工程背景，分别对相似材料实验中的第16次和第29次开挖(工作面推进距离分别为118、255 cm)覆岩应力稳定后模型表面出现的裂隙进行识别。首先采用Agisoft Mateshape图像处理平台，对现场采集的影像进行拼接、校正，生成正射影像(图1)。由图1可知，拟处理区域环境较为复杂，模型表面光线偏暗，而且模型架、护板、网格线、编码点和非编码点等都可能影响裂隙的识别效果。因此，首先需要对图2进行预处理(剔除模型架、护板等)，消除潜在的外部影响因素，再进行后续的图像处理。

图1 工作面不同推进度时影像对比度增强处理Fig.1 Image contrast enhancement processing under different advancing degrees of working face

为了对比Canny算法优化前后的阈值分割效果，作者基于Matlab 研发平台，编制了计算程序，对比分析算法优化前后两幅影像数据的阈值分割效果以及模型表观裂隙识别结果，验证算法的有效性。由于实验获取的模型图像整体较暗，识别目标与背景的灰度值接近，因此需要对选取的各幅影像进行对比度增强预处理，处理后的结果如图2所示。基于式(13)，对图2(a)、图2(b)进行阈值分割，分割结果如图2所示。根据式(13)，确定的两幅影像权重系数分别为5.73%、7.68%。对比优化前后各幅影像可知，算法优化前，各影像均不同程度地存在裂隙边缘噪声、裂隙轮廓模糊不清，严重影像裂隙识别的精准度。算法优化后，各幅影像表观裂隙纹理清晰，对比度明显增强，阈值分割效果显著提升。

图2 算法优化前后的阈值分割结果Fig.2 Threshold segmentation results before and after algorithm optimization

完成图像阈值分割之后，即可进行相似材料表面裂隙的边缘识别与连接，工作面推进至118 cm时算法优化前后的裂隙识别结果如图3所示，工作面推进至255 cm时算法优化前后的裂隙识别如图4所示。

图3 工作面推进至118 cm时算法优化前后的裂隙识别结果Fig.3 Fracture identification results before and after algorithm optimization when the working face advances to 118 cm

对比图4(a)和图4(b)可以看出，整个裂隙区划分为Ⅰ、Ⅱ、Ⅲ三个部分，使用Canny常规算法对模型全景图像进行处理后，第Ⅰ部分识别的裂隙模糊不清，部分裂隙未能有效识别，第Ⅱ部分主要裂隙未能识别，第Ⅲ部分裂隙模糊不清。算法优化后，裂隙的整体识别率较高，裂隙边缘轮廓较为清晰。此外，从图4(b)可以看出，覆岩裂隙主要集中在工作面两端(开采边界的内侧)，根据矿山开采沉陷学“上三带”划分方法，当工作面推进度为255 cm时，上覆岩层中垮落带高度为13.31 cm，断裂带高度为63.36 cm(由于只模拟部分岩层，弯曲带高度无法确定)，左右垮落角分别为52°和55°，裂隙场整体呈“梯形”。Canny优化算法实现了相似材料模型表面主要裂隙和细微裂隙的准确识别，裂隙损失小，裂隙场轮廓较为清晰，证实了算法的有效性和实用性。

图4 工作面推进至255 cm时算法优化前后的裂隙识别结果Fig.4 Fracture identification results before and after algorithm optimization when the working face advances to 255 cm

2.3 裂隙定量分析

采动覆岩裂隙能够响应裂纹的萌生、扩展和动态运移过程，反馈上覆岩层的受损状态。开展采动覆岩裂隙定量提取和分析研究对于井下煤层气抽放、水体(含水层)下开采具有一定的实用价值和意义。在裂隙识别的基础上，作者拟采用像素法和计盒维数法进一步从微观角度对裂隙进行定量提取、分析和相互验证。

(1)基于像素法裂隙定量分析。像素法是对采动覆岩裂隙进行定量分析的常用方法之一，评价指标采用裂隙率。裂隙率是拟识别图像裂隙区域面积与识别区域总面积之比，面积用区域内像素个数之和表示。采用像素法可分别对工作面不同推进度的采动覆岩表观裂隙进行定量分析，结果见表1。由表1可知，算法优化后较优化前上覆岩层的裂隙率和裂隙增长率均有不同程度的提高，且裂隙率和裂隙增长率与工作面开采尺寸正相关，开采尺寸越大裂隙率越大。当工作面推进至255 cm时算法优化前后的裂隙率分别为5.32%、7.41%，裂隙率增大了21.53%。这是因为算法优化前部分裂隙模糊不清，类间方差函数的引入使原来未能识别的裂隙得以识别，裂隙识别率得到了不同程度的提升。

表1 基于像素法相似材料模型采动覆岩表观裂隙定量分析结果Tab.1 Quantitative analysis results of apparent fractures of mining overburden based on pixel method similar material mod

(2)基于分形维数的裂隙定量分析。分形维数是一种分形的定量表征，能定量化描述自然界不规则、杂乱无章的复杂现象。原理是利用计算机存储技术基于点覆盖法，计算图像各个状态的分形维数，从而反映图像裂隙的微观状态。本部分拟采用计盒维数法，基于Matlab平台编制计算程序，实现工作面不同推进度时裂隙分形维数的自动计算和分析。计盒维数法首先将图像进行二值化处理，再把图像分割成边长为s的小正方形，然后将包含裂隙的小正方形记作1，其余的记作0。若包含裂隙的小正方形数量用N(s)表示，则改变s的大小，N(s)的大小也随之变化。实际上分形维数计算就是绘制lg(s)-lgN(s)之间的关系曲线，而拟合直线斜率的绝对值即为分形维数。计算公式见式(14)：

(14)

式中，D为裂隙的分形维数；s为小正方形的边长；

N(s)为包含裂隙的小正方形数量。

根据分形维数理论，对工作面不同推进度三幅覆岩裂隙场影像进行了图像分形维数计算，结果见表2和图5所示。

表2 工作面不同推进进度裂隙场分形维数计算Tab.2 Calculation of fractal dimension fracture field in different advancing progress of working face

由表2可知，两幅影像分形维数计算相关系数均较高，拟合误差较小，具有良好的线性相关性，表明相似材料模型在开挖过程中，上覆岩层裂纹生长具有自相似的分形特性，反映了不同推进状态下上覆岩层裂隙发育的表观结构特征。当工作面推进度分别为118 cm和255 cm时，采动覆岩裂隙场的分形维数分别为1.811和1.885，分形维数为1～2。分形理论认为，二维条件下若分形维数位于0～2，则分形维数值符合客观规律。因此，采用优化Canny算法提取的3个裂隙场均具有较好的分形特征。其中，工作面推进度为255 cm时的裂隙场分形维数最大，说明此状态下图像表观裂隙被等分后，所占小正方形数量最多，即裂隙率最高，这与像素法定量分析的结果一致。

图5 工作面不同推进进度时计盒数和小正方形边长之间的拟合关系Fig.5 Fitting between the number of meter boxes and the side length of small square at different progress of working face

3 结论

(1)以15235工作面为工程背景，构建了几何相似比为1/100的相似材料实验模型，通过对定焦相机批量获取的多基线近景摄影测量影像进行拼接和校正，获得了目标区的正射影像。

(2)针对复杂环境下Canny算法裂隙识别的不足，构建了顾及权重系数的类间方差函数，确定了影像信息分割最优阈值表达式，实现目标影像前景和背景的完整分割以及裂隙信息的完整识别。

(3)基于Matlab研发平台，采用Canny优化算法，得到两组图像类间方差权重系数分别为5.73%、7.68%，影像信息分割最优阈值分别为115、128，裂隙场信息提取完整，左右垮落角分别为52°和55°，覆岩垮落带和断裂带高度分别为13.31 cm和63.36 cm，算法优化后裂隙率相对增量明显提升。

(4)采用像素法和分形理论分别对工作面推进度为118 cm和255 cm时的覆岩裂隙场进行了定量分析，上覆岩层裂纹生长具有自相似的分形特性，反映了不同推进状态下上覆岩层裂隙发育的表观结构特征。