基于光反射率突变的受力岩石裂纹提取方法

高祥， 高亚楠

（1. 安徽理工大学 矿业工程学院，安徽 淮南 232001；2. 安徽理工大学 深部煤矿采动响应与灾害防控国家重点实验室，安徽 淮南 232001）

0 引言

在巷道、硐室等地下工程中，围岩裂纹发育、贯通是导致围岩结构失稳的主要原因，而裂纹识别与提取是认识围岩裂纹演变规律的基础。常见的裂纹提取方法有可见光拍照[1]、CT扫描[2]、电镜扫描[3]、声学成像[4]、红外热成像[5]、薄片显微观测[6]等。其中，CT扫描和声学成像法能够探测到岩石内部裂纹发育信息，但设备昂贵，现场应用困难；可见光照片在岩石力学实验中应用最为广泛。此外，在智能矿山建设过程中，视频监测装置已被广泛应用[7]，获得了大量视频数据，为研究基于可见光照片的受力岩石裂纹提取方法提供了数据来源。

基于可见光照片的受力岩石裂纹提取方法主要包括阈值分割二值化法[8-9]、边缘检测法[10-11]、区域生长法[12-13]等。这些方法对单张照片的分析较方便快捷，对由颜色相近的矿物组成的岩石（如砂岩）裂纹的提取效果较好[14]，但对于由多色矿物组成的复杂岩石（如花岗岩），由于不同矿物灰度存在较大差异，矿物边界会被一并搜出，对裂纹提取造成干扰。另外，对视频监测数据进行分析时，需要首先确定视频图像中是否存在裂纹。实际上，岩石裂纹之所以能被肉眼识别，根本原因在于裂纹反射率与岩石反射率不同。在岩石受力变化过程中，不仅红外反射率[15-17]、微波反射率会发生变化[18-19]，可见光反射率也会发生变化。这种可见光反射率的变化可作为裂纹提取的重要依据。而且，基于光反射率的裂纹识别方法受环境光辐射的影响较小，更适用于光照条件复杂多变的井下。因此，本文在分析试件破坏后的光反射率响应及其空间分异特征的基础上，提出一种基于光反射率突变的受力岩石裂纹提取方法。

1 基于光反射率突变的受力岩石裂纹提取方法

基于光反射率突变的受力岩石裂纹提取流程如图1 所示。关键步骤如下：① 光反射率计算，得到由光反射率表示的数字图像。② 数字图像时间序列差分计算和光反射率突变点识别。③ 数字图像光反射率突变点空间分异性分析。④ 裂纹信息提取。

图 1 受力岩石裂纹提取流程Fig. 1 Crack extraction process of stressed rock

1.1 光反射率计算

物质实测反射率与数字图像灰度之间具有显著的线性相关性[20]。根据反射率波谱野外实地测量方法，当目标物体和标准板的环境条件相同时，目标物的光反射率为[21]

式中：Kb为 标准板的光反射率；N，Nb分别为目标物、标准板的可见光数字图像灰度。

一般地，标准板由硫酸钡、聚四氟乙烯等材料制成，具有固定的反射率，是用于反射率测量的标准漫反射参考物。

根据式（1）可方便地将数字图像各像元的灰度转换为光反射率，进而得到由光反射率表示的数字图像。数字图像中位置（i，j）处像元的光反射率为

式中f为光反射率关于位置的空间函数。

某一时刻t数字图像中位置（i，j）处像元的光反射率为

式中g为光反射率的时空函数。

1.2 光反射率突变点识别

裂纹发育与光反射率的突变信息相对应。为了便于反映裂纹发育前后光反射率的突变性，对数字图像时间序列进行差分处理。

式中：Vti j为t～t+Δt时间段内数字图像中位置（i，j）处像元的光反射率变化速率；Δt为时间变化值。

一般地，Δt取值不能太大，因为Δt较小时才能反映裂纹发育时刻光反射率的突变性。

假设在受力岩石裂纹未发育时，对于用光反射率变化速率表示的数字图像，其灰度呈正态分布，则光反射率变化速率太大或太小都是小概率事件。利用数理统计理论可以确定数字图像中光反射率突变的判别阈值。

式中： [Vt]up，[Vt]down分别为t时刻光反射率突增、突降的判别阈值；E(Vt)，D(Vt)分别为用光反射率变化速率Vt表示的数字图像灰度的均值和方差；k为常数，取3或4。

当Vti j＞[Vt]up或Vti j＜[Vt]down时，认为位置（i,j）处像元的光反射率发生了突变，称该类像元为光反射率突变点。

1.3 光反射率突变点空间分异性分析

受力岩石裂纹发育可引起光反射率突变，但光反射率突变并不意味着裂纹发育，二者不是一一对应的关系。裂纹发育引起的光反射率突变点在空间上的分异性是识别裂纹的关键特征。根据光反射率突变点的空间分异性可以判断数字图像中是否存在裂纹。为此，引入反映空间点分布集中程度的指标-分异系数C[22]。

对数字图像中的光反射率突变点进行样方分析[23]，即通过突变点分布密度的变化来评定其分布模式（包括聚集分布、均匀分布和随机分布）。将数字图像划分成方格网（方格尺寸应结合裂纹宽度确定，建议方格对应的空间分辨率为 2 mm×2 mm），然后统计落入每个方格中的突变点数量，得到样方点量数列{x(k),k=1,2,···,m}（m为样方数）。分异系数为样方点量数列的方差D(x)与 均值E(x)的比值：

利用式（6）可判断数字图像时间序列中每帧数字图像的光反射率突变点是否呈聚集分布。当C＜1时，光反射率突变点呈均匀分布；当C=1 时，光反射率突变点呈随机分布；当C＞1 时，光反射率突变点呈聚集分布，C值越大，对应的空间分异性越强。

1.4 裂纹信息提取

利用试件加载过程中光反射率突变点分异系数的变化曲线，可判断出试件裂纹发育时刻，为后续裂纹信息提取和空间特征分析提供支持。

判断出数字照片中存在裂纹后，一方面，通过光反射率变化速率反映裂纹活动的瞬时信息，以研究试件在何时、何处有裂纹产生或扩展；另一方面，将考察时间段内所有突变点的空间信息在同一张数字图像中显示，提取出裂纹活动的累积信息，以分析一段时间内岩石表面裂纹的发育、扩展情况及贯通程度。

2 受力岩石裂纹提取实验系统

受力岩石裂纹提取实验系统主要由压力机、快速照相机、声发射系统及标准板等组成，如图2所示。实验所用试件由150 mm×150 mm×50 mm的花岗岩岩块经四角切割、中心钻孔制备而成，孔直径为40 mm。压力机型号为RLW-3000，可进行双轴加载，最大轴压为3 000 kN，侧压可达1 000 kN，用于对试件施加荷载，模拟地应力环境。快速照相机型号为GX3300，最高分辨率为3 296像素×2 472像素，最大帧率为17 帧/s，用于可见光照片实时拍摄。采用8CHs PCI-2声发射系统[24]对加载过程中试件的裂纹发育信息进行实时监测。

图 2 受力岩石裂纹提取实验系统Fig. 2 Experimental system of crack extraction of stressed rock

实验前，将试件放进压力机加载仓内，将声发射探头布置在试样背面圆孔上方。照相机由三脚架支撑，架设在距离试件大于2 m的位置，对试件正面进行垂直观测，采样帧率设为10 帧/s。在照相机监测过程中，由固定光源照射。根据试件尺寸和数字图像中试件所占用像素，经计算可知，每个像素的大小为0.088 7 mm。采用聚四氟乙烯作为标准板（Kb为0.93）。

采用双轴加载、应力控制方式进行实验，其中垂直轴为主轴。实验开始后，先将水平荷载预加载至100 kN（20 MPa）保持恒定，以提供稳定的水平应力，随后垂直方向以1.2 kN/s 匀速加载，通过快速照相机、声发射系统同步监测，直至试件破坏。实验结束后，对可见光照片中试件和标准板进行截取分析。

标准板灰度分析结果如图3所示。图3（a）中，X，Y为像素坐标。可看出标准板的灰度在一定范围内波动，中位数在230附近，经计算，灰度均值为231。这说明在实验光照条件下，可取 231作为标准板的灰度，代入式（1）进行光反射率计算。

图 3 标准板灰度分析结果Fig. 3 Grayscale analysis results of standard plates

3 实验结果分析

3.1 光反射率突变响应

裂纹活动包括裂纹发育、扩展和贯通等。裂纹活动过程中光反射率突变与否，一般从2个角度来衡量：① 变化是否在较短时间内进行；② 变化量是否较大。为了验证裂纹发育、扩展可引起光反射率突变，选择代表区域AR0，拍摄加载过程中第653，654 s时AR0区域的裂纹发育照片，如图4所示。图4中，AR0区域大小为17 mm×8 mm，面积较小。将照片放大后，由于分辨率较小，照片较模糊，但仍可以看出：在第653 s，AR0区域裂纹不够明显；到第654 s，裂纹依稀可见。这说明在653～654 s，AR0区域有裂纹产生或扩展。

图 4 AR0区域裂纹活动前后照片Fig. 4 Photographs of AR0 region before and after crack activity

裂纹的光反射率变化速率空间分布如图5所示。图5 用灰度反映光反射率变化速率，图中像素越亮，代表其光反射率变化速率越大。对比图4和图5可看出，在数值上，裂纹发育的区域呈白色，对应的光反射率变化速率较大，最大超过0.2/s；在空间上，裂纹活动区域的光反射率变化具有明显的空间分异性。

图 5 AR0区域裂纹活动前后光反射率变化速率空间分布Fig. 5 Spatial distribution of change speed of light reflectivity before and after crack activity in AR0 region

AR0区域无裂纹活动时刻（第300 s）光反射率变化速率直方图如图6所示。可看出无裂纹活动时刻，试件表面光反射率变化速率近似呈正态分布，3倍标准差为0.03/s。结合图5 和图6 可看出，裂纹活动引起的光反射率变化速率（接近 0.2/s）远大于无裂纹活动时刻的光反射率变化速率（小于0.03/s）。由以上分析可看出，试件受力后的裂纹发育、扩展可在较短时间内引起光反射率的显著变化。

图 6 AR0 区域无裂纹活动时刻光反射率变化速率直方图Fig. 6 Histogram of light reflectivity variation at the time without crack activity in AR0 region Histogram of light reflectivity variation

3.2 光反射率变化的空间分异性

3.2.1 光反射率突变点空间分布

光反射率突变点空间分布如图7所示。其中白色小点为突变点，黑色为非突变点。由图7（a） 可看出，在无裂纹发育时刻，试件表面光反射率突变点分布较分散，相对均匀，单位面积内突变点的数量变化较小。由图7（b）可看出，在裂纹发育时刻，试件表面光反射率突变点多呈簇状或条带状分布，单位面积内突变点的数量变化较大。

图 7 光反射率突变点空间分布Fig. 7 Spatial distribution of mutation points of light reflectivity

3.2.2 光反射率突变点分异系数变化特征

试件加载过程中光反射率突变点的分异系数变化曲线如图8 所示。可看出在0～597 s（阶段I和II），光反射率突变点的分异系数在小范围内变化，标准差为1.10，平均值为3.01。分异系数与临界值1接近，说明该时间段内试件表面光反射率突变点接近随机分布，无裂纹活动或活动不剧烈。在597～657 s（阶段III），光反射率突变点的分异系数变化幅度较大，最大值可达189，远远大于临界值1，说明该时间段内试件表面光反射率突变点多呈聚集分布，裂纹活动剧烈。

图 8 试件加载过程中光反射突变点分异系数变化曲线Fig. 8 The variation curve of differentiation coefficient of mutation points of light reflectivity during specimen loading process

试件加载过程中声发射突变率变化曲线如图9所示。声发射突变率反映了声发射率的突变性[24]。可看出在阶段I内，声发射不活跃，说明该阶段无裂纹产生。以该阶段分异系数C的均值加3倍标准差为阈值（5.12），可对分异系数异常时刻（裂纹活动时刻）进行自动筛选。

图 9 试件加载过程中声发射突变率变化曲线Fig. 9 The variation curve of acoustic emission mutation rate during specimen loading process

另外，可见光数字图像反映了岩石表面信息，声发射反映了岩石表面和内部信息，将二者进行对比分析，可以研究岩石受力破坏过程中内部和表面损伤的关系问题。

3.3 裂纹信息提取

3.3.1 瞬时信息提取

裂纹发育、扩展瞬时信息如图10所示。可看出在第604 s 试件表面有裂纹发育、扩展，裂纹活动区域集中于孔洞左侧；在第623 s 试件表面也有裂纹发育、扩展，裂纹活动区域转移到孔洞右侧。这说明在裂纹发育过程中，伴随着应力集中的转移，裂纹活动区域也会转移，并非总沿着已有裂纹持续扩展。

3.3.2 累积信息提取

604～623 s裂纹发育、扩展累积信息如图11所示。可看出在该时间段内试件孔洞两侧都有裂纹活动，累积信息反映了考察时间段内试件表面裂纹发育、扩展的整体情况。

图 10 裂纹发育、扩展瞬时信息Fig. 10 Instantaneous information of crack development and propagation

图 11 裂纹发育、扩展累积信息Fig. 11 Accumulation information of crack development and propagation

需要注意的是，裂纹活动的累积信息不宜通过2个时刻的数字图像（像元灰度为光反射率）相减的方法得到。因为，当时间间隔较长时，该时间段内岩石应力变化可能较大，引起的岩石表面位移和变形会大于数字图像的空间分辨率，从而造成光反射率整体大幅度变化。

4 讨论

4.1 岩石破裂过程中光反射率变化原因

岩石光反射率的大小除了与岩石内部矿物的电磁性质、环境光辐射特征有关外，也与岩石表面的粗糙度紧密相关。前2个因素分别与矿物种类和实验环境条件有关，在实验中可视为常量。在岩石受力变形破裂过程中，岩石表面粗糙度不断发生变化，是导致光反射率变化的主要因素。一般地，岩石裂纹的产生具有突变性，所以，岩石表面光反射率的变化也具有一定突变性。

4.2 所提裂纹提取方法的岩石力学意义

基于光反射率突变的受力岩石裂纹提取方法可为研究受力岩石破坏演变特征、认识岩石表面与内部损伤之间的关系提供有效方法。

4.2.1 受力岩石破坏演变研究

试件破坏过程中光反射率变化空间分布如图12所示。该图反映了试件破坏演变过程：① 第566 s，破裂先在孔洞左侧发生，集中于孔洞周边边帮的中点附近，破裂范围较小，产生的岩屑会发生弹射。② 第632 s，裂纹发育区域发生转移，由孔洞左侧转移到了右侧。③ 第641 s，孔洞两侧裂纹继续扩展，扩展到一定程度之后，在孔洞左右两侧上方远离孔洞的区域产生新的裂纹，但新裂纹与孔洞周边的破裂区域未贯通。④ 第652 s，裂纹继续扩展，并趋于贯通，意味着岩石即将破坏。实际上，试件在第656 s出现了垮塌，圆形孔洞发生严重破坏，但应力仍未达到峰值，这与现场巷道围岩破坏后仍具有承载能力的现象一致。

基于光反射率突变的花岗岩试件裂纹提取结果如图13所示。在数字照片中，裂纹很难看清；而在同时刻的光反射率变化空间分布图中，清楚地显现出了裂纹（白色条纹），为下一步判断裂纹性质及其发生机制打下基础。

图 12 试件破坏过程中光反射率变化空间分布Fig. 12 Spatial distribution of light reflectivity changes during specimen failure

4.2.2 岩石表面损伤与内部损伤关系研究

声发射是研究岩石破坏演变的重要工具，通过实时监测受力岩石裂纹发育扩展中所释放的高频振动信号，可获取岩石破裂、损伤发育的剧烈程度。需要注意的是，通过声发射监测获得的是岩石内部和表面的整体信息，而通过可见光照片仅能获得表面裂纹的活动信息。

光反射率突变点的分异系数是个综合指标，如果分异系数出现异常，说明光反射率发生了突变。突变点的空间分异性是反映岩石表面有无裂纹发育、扩展的重要指标。

为了研究试件内部损伤与表面损伤之间的关系，对图8和图9 进行对比分析。在阶段I（0～420 s）、阶段III（597～657 s）、阶段IV（657～703 s），声发射突变率与光反射率突变点分异系数变化特征接近；在阶段II（420～597 s），前者活跃，而后者表现较平稳，二者出现了显著差别，声发射参数比光反射率参数活跃得早。这说明岩石内部损伤要早于表面损伤发育，表面破坏可以作为岩石整体失稳的后期前兆。

由于可见光穿透性较差，通过可见光数字图像只能获取围岩表面裂纹信息，后期可以结合微波激光雷达、地质雷达等具有穿透性的设备，定期巡检围岩内部损伤及破坏情况。

图 13 基于光反射率突变的花岗岩试件裂纹提取结果Fig. 13 Crack extraction results of granite specimens based on mutation of light reflectivity

4.3 智慧矿山巷道安全巡检中的潜在应用

当前，智慧矿山建设已经成为矿山发展的趋势[25]，视频监测系统被称为智慧矿山的“眼睛”。杨春雨等[26]提出将智能视频技术用于巷道巡检机器人，对巷道变形情况进行监测，以及时发现巷道变形与支护隐患，保障煤矿生产安全。岩石表面裂纹剧烈发育、贯通是破坏的重要前兆信息，对岩体破坏具有重要的指示意义。如果在巷道巡检机器人平台上搭载高清摄像机，对巷道进行动态巡检，可得到大量数字图像，结合本文提出的基于光反射率变化的受力岩石裂纹提取方法，可实时提取巷道围岩裂纹发育信息，进行巷道围岩稳定性评估，对保障矿山安全具有重要意义。然而，本文仅以实验为基础开展岩石裂纹提取方法研究，由于室内实验环境与现场环境存在较大区别，要真正搭载到巷道巡检机器人平台，用于实时提取巷道围岩裂纹发育信息，还存在很多基础问题需要研究。比如，巡检机器人要在移动中监测，不同时间在相同位置巡检得到的系列数字图像如何精确匹配；在巷道黑暗、多光源条件下，光源稳定性对光反射率反演精度会造成影响。

5 结论

（1） 岩石受力裂纹发育、扩展可引起光反射率突变，最大突变速率可达0.2/s，远大于由其他随机因素导致的变化速率（0.03/s）；在裂纹活动时刻，光反射率突变点呈显著的空间分异现象（C值可达189），分异系数远大于随机分布时的值（C值为1）。这种空间分异性是受力岩石裂纹提取的重要依据。

（2） 基于光反射率突变的受力岩石裂纹提取方法利用分异系数判断岩石裂纹活动的时刻，大大提高了搜索效率和自动化程度。该方法能够提取出裂纹活动的瞬时信息和累积信息，分别用于研究裂纹活动的时空特征、判断裂纹贯通程度。

（3） 讨论了岩石破裂过程中光反射率变化的原因、所提方法的岩石力学意义及其在巷道安全巡检中的应用前景：在岩石受力破坏过程中，表面粗糙度是导致光反射率变化的主要因素；所提方法可为研究受力岩石破坏演变特征、认识岩石表面与内部损伤之间的关系提供了有效方法；结合巷道巡检机器人、高清摄像机和本文方法，可实时提取巷道围岩裂纹发育信息，进行巷道围岩稳定性评估，对保障矿山安全具有重要意义。