裂隙砂岩变形破裂过程中应变场灰度及纹理特征分析

张 科，张 凯

(1.昆明理工大学 电力工程学院，云南 昆明 650500; 2.昆明理工大学 建筑工程学院，云南 昆明 650500)

在经历了漫长复杂的地质构造作用后，岩石内部不可避免地发育有大量的地质缺陷，这些缺陷对岩石的力学特征产生了重要的影响[1-2]。地下采矿场的矿柱中普遍存在裂隙，在外部荷载作用下容易诱发裂隙扩展、贯通，最终导致矿柱的失稳破坏[3-5]。

数字图像相关技术(Digital Image Correlation Method，DIC)是一种非接触式全场变形观测方法[6-7]，可以连续观测试件表面的位移场和应变场，具有操作简单、测量精度高、试验消耗小、可重复测量等优点，已应用于岩石力学测试。李露露等[8]利用DIC技术得到单轴压缩下三叉裂隙类岩石试件的应变场。雷冬等[9]基于DIC技术计算出花岗岩在动力破坏过程中的应变场，与数值模拟结果一致。潘红宇等[10]结合DIC技术，分析煤岩表面应变场变化特征，深入研究了轴向应变率和裂纹扩展速率之间的关系。大久保诚介等[11]和YANG等[12]利用3D-DIC技术实时跟踪安山岩和煤岩件在单轴压缩破坏过程中的位移场和应变场演化规律。上述研究表明，通过DIC技术可以得到试件加载全过程的应变场。但是，在后续分析中，有必要结合一些具有物理意义的、能够反映应变场动态变化的定量指标，以便更直观、更准确地把握岩体变形破裂规律。

近些年，研究发现灰度和纹理等图像特征可以表征岩石破裂状态，因而受到了学者们的关注与尝试。徐金明等[13]和XU等[14]利用图像纹理特征参数对单轴压缩试验的视频图像进行分析，认为这些参数可以反映出石灰岩的变形状态和花岗岩的裂纹萌生、扩展过程。宋勇军等[15]分析了CT图像的灰度直方图，研究了单轴加载下冰冻红砂岩的裂隙扩展规律。朱红光等[16]通过对CT图像的灰度统计，分析了裂纹演化过程中图像灰度变化的统计规律。

目前的图像特征研究主要集中于视频拍摄图像和CT图像，对应变场云图却鲜有研究。而基于DIC技术得到的应变场云图能够更加真实地反映出岩体在加载过程中的变形破裂特征。因此，笔者以地下开采中矿柱的力学问题为工程背景[3-5]，结合数字图像相关技术和图像特征分析方法，通过应变场云图灰度直方图以及10个图像特征描述指标，量化研究了单轴压缩过程中含裂隙的矿柱岩体应变场灰度及纹理特征演化规律，为精细化研究裂隙岩体变形破裂机制提供一种新的思路和方法，对矿柱稳定性评价及灾变预警也有重要的参考价值。

1 研究方法

1.1 试验设备与方案

选取云南省昆明市地区的红砂岩作为试验对象，现场采集红砂岩岩块，将其加工成120 mm×60 mm×20 mm(长×宽×厚)的试件，如图1所示，端面不平整度误差<0.05 mm。挑选质地均匀、外观完整的试件，采用高速水刀切割技术加工贯穿裂隙，裂隙形心与试件形心重合，裂隙长度2b=28 mm。本次试验共设计3种裂隙倾角β，分别为30°，45°和60°，对应的试件编号分别记为S30，S45和S60，见表1。在预制裂隙中填充云母片，模拟充填闭合型裂隙。制备完成的不同裂隙倾角砂岩试件如图2所示。另外，制备完整砂岩试件，编号记为SW，用于对比分析。砂岩试件的物理力学参数见表1。

表1 砂岩试件物理力学参数Table 1 Physical and mechanical parameters of sandstone specimens

图2 待试验的裂隙砂岩试件Fig.2 Fractured sandstone specimens ready to be tested

为增强图像纹理效应[17]，通过喷漆的方法制作试件表面的人工散斑场:对试件正面喷涂白漆，待其干透后再喷涂黑漆做随机散斑处理，喷涂过程尽可能均匀，保证最后制作出优良的人工散斑，如图2所示。

矿柱是支撑地下结构稳定的关键结构，近似处于一维受压应力状态(图3)，故可通过单轴压缩试验探究其破坏特征[4-5]。本次试验系统如图4所示。采用WDW-100E万能材料试验机进行单轴压缩试验，该设备最大轴向荷载为100 kN，测量精度为0.5%。试验加载采用位移控制方式，加载速率设置为0.3 mm/min。试验时，采用图像分辨率为2 592×1 944像素的工业相机，对试件表面进行散斑图像采集，采集速率设置为1张/s。试验前，保证加载设备和监测设备的时间一致。

图3 含裂隙的矿柱[3]Fig.3 A pillar containing flaw[3]

图4 试验系统示意Fig.4 Schematic diagram of experimental system

1.2 图像灰度特征参数

灰度是描述像素敏感程度的整数量，反映了图像色度或亮度等级的分布特征[18]。灰度等级可划分为0～255，图像中像素点越白，相应的灰度级越高，反之灰度级越低。依据数字图像处理理论[19]，灰度直方图定义为不同灰度等级概率密度的分布函数，其积分函数称为累计分布函数，反映了图像中某灰度等级与其出现频率之间的关系。灰度直方图的横坐标为灰度等级，纵坐标为对应灰度级在图像中出现的频数，数学表达式为

(1)

式中，pr(k)为k灰度级出现的概率；k为灰度级;nk为灰度级为k的像素点个数;n为像素点总数。

笔者对应变场云图的灰度直方图进行计算，采用均值、标准差、三阶矩、平滑度、一致性以及熵6种灰度特征参数定量描述其灰度特征[20]。

(1)灰度均值m，表示图像的平均灰度值，数学表达式为

(2)

(2)灰度标准差σ，表示图像灰度的离散程度，数学表达式为

(3)

(3)灰度三阶矩μ，表示图像灰度分布的不对称程度，数学表达式为

(4)

(4)灰度平滑度R，表示图像灰度集中于均值附近的程度，数学表达式为

(5)

(5)灰度一致性U，表示图像灰度分布的统一程度，数学表达式为

(6)

(6)灰度熵h，表示灰度分布的不均匀、混乱程度，数学表达式为

(7)

1.3 图像纹理特征参数

在图像空间中，一定距离内的2个像素之间存在一定的灰度关系，即图像中灰度的空间相关特性。它们呈现出的某种规律变化，形成了纹理特征。目前，提取纹理特征的方法很多[21]，灰度共生矩阵(Gray Level Co-occurrence Matrix，GLCM)对描述纹理有很强的适应性和鲁棒性[22]，故本文采取灰度共生矩阵提取纹理特征。

基于GLCM，笔者采用对比度、相关性、角二阶矩以及熵4种纹理特征参数定量描述应变场云图的纹理特征[23-24]。

(1)纹理对比度Con，表示图像中相邻像素点间的灰度级差异，数学表达式为

(8)

式中，矩阵P(i,j,d,θ)为像素点(x,y)的所有θ方向，相邻间隔为d的像素对中1个取i值，另1个取j值的相邻像素对出现的次数。

(2)纹理相关性Cor，表示图像中相邻像素点间灰度级的相似程度，数学表达式为

(9)

(3)纹理角二阶矩ASM，表示像素灰度级分布的均匀程度，数学表达式为

(10)

(4)纹理熵H，表示像素灰度级的复杂程度，数学表达式为

(11)

为了抑制方向分量，使计算的纹理特征与方向无关，通常引入4个方向(θ一般取0°，45°，90°和135°)的均值和标准差[25]来表示，本文采用标准差S表征纹理特征:

(12)

2 计算结果与分析

单轴压缩作用下3种裂隙倾角的砂岩试件以及完整砂岩试件的应力-应变曲线如图5所示。表1列出了对应的峰值应力σc和弹性模量E，从表1可以看出，裂隙砂岩试件的极限承载力显著低于完整砂岩试件(σc=75.05 MPa)，这说明预制裂隙的存在改变了岩石的力学响应。以加载前的数字散斑图像为参考图像，采用数字图像相关软件Ncorr[26]进行计算，得到试件加载全过程的应变场。研究发现，各裂隙倾角试件变形破裂过程中的应变场演化规律大致相似;相比于其他类型的应变场，水平应变场的变化规律最为明显。限于篇幅，以β=45°砂岩试件的水平应变场灰度和纹理分析结果作为典型进行分析。结合WONG等[1]以及MIAO等[27]对裂纹类型的划分(图6)，本次试验主要观察到张拉裂纹和剪切裂纹，分别记为T和S。图7和8分别为5种应力水平下的典型裂隙砂岩试件(β=45°)裂纹扩展和水平应变场云图(拉应变为正，压应变为负)。

图6 裂纹类型示意[1,27]Fig.6 Schematic diagram of crack patterns[1,27]

图7 加载过程典型裂隙砂岩试件裂纹扩展Fig.7 Crack propagations of typical fractured sandstone specimen during loading process

图8 加载过程典型裂隙砂岩试件水平应变场Fig.8 Horizontal strain fields of typical fractured sandstone specimen during loading process

2.1 应变场和灰度直方图演化规律

将水平应变场的数据矩阵转化为灰度图像，而后进行二值化处理，最后利用Matlab自编程序，求解灰度图像对应的灰度直方图。图9为不同应力水平下的典型裂隙砂岩试件(β=45°)灰度直方图以及频率峰值变化曲线。

由图7～9可以看出：① 当加载至0.08σc，应变场分布较均匀，如图8(a)所示;对应的灰度直方图表现为“矮胖”型(峰值频率约为0.019)，此时像素点灰度大部分集中在灰度区间中部，如图9(a)所示。② 当加载至0.30σc，此时试件处于弹性变形状态，未出现新生裂纹，如图7(b)所示;但从应变场云图上可以看出，应变集中区已在预制裂隙周围形成，裂隙尖端处的应变值最大(图8(b))，远大于其他区域的应变，成为最薄弱的部位，易萌生裂纹;对应的灰度直方图基本维持不变，如图9(b)所示。③ 当加载至0.93σc，预制裂隙尖端出现张拉裂纹，其扩展方向大致与轴向压应力方向平行，如图7(c)所示，与WONG等[1]以及MIAO等[27]对张拉裂纹的分类结果吻合(图6(a));裂隙及新生裂纹周边产生应变局部化带(图8(c))，此时，灰度直方图中峰值频率增加为0.032，应变场灰度整体向高灰度级方向移动，如图9(c)所示。④ 当加载至0.95σc～1.00σc阶段，出现剪切裂纹，与WONG等[1]、MIAO等[27]对剪切裂纹的分类结果吻合(图6(b));应变局部化带也随之发生改变，沿试件对角方向快速发展如图8(d),(e)所示，此时应变场灰度直方图中峰值频率增加至0.054，从“矮胖”型转变为“高瘦”型，应变场灰度整体表现为向更高等级的灰度方向移动，如图9(d),(e)所示。

图9 加载过程典型裂隙砂岩试件应变场灰度直方图Fig.9 Gray scale histograms of strain fields of typical fractured sandstone specimen during loading process

2.2 灰度与纹理特征参数演化规律

计算加载过程中灰度图像的均值、标准差、三阶矩、平滑度、一致性以及熵6种灰度特征参数;为便于比较分析，对各参数进行最大最小值归一化处理。典型裂隙砂岩试件(β=45°)计算结果如图10所示。

将上述灰度图像转化为灰度共生矩阵，分别计算0°，45°，90°和135°方向对应的对比度、相关性、角二阶矩和熵4种纹理特征参数。为了抑制方向分量的影响，求解4个方向纹理特征参数的标准差，并进行最大最小值归一化处理。典型裂隙砂岩试件(β=45°)计算结果如图11所示。

图10 典型裂隙砂岩试件灰度特征参数-时间曲线Fig.10 Gray feature parameters-time curves of typical fractured sandstone specimen

图11 典型裂隙砂岩试件纹理特征参数-时间曲线Fig.11 Texture feature parameters-time curves of typical fractured sandstone specimen

由图10和11可以看出，典型裂隙砂岩试件(β=45°)的灰度/纹理特征参数-时间曲线可大致划分为3个阶段:① 波动段(阶段I:1～89 s)，试件处于微缺陷压密阶段，此时形成早期的非线性变形，反映在灰度/纹理特征参数上表现为大幅波动。② 平稳段(阶段II:90～198 s)，试件进入弹性变形状态，试件表面各像素点的灰度/纹理特征参数虽然仍在变化，但整体分布趋于稳定，相应的纹理特征参数-时间曲线呈现出平稳化趋势。③ 突变段(阶段III:199～235 s)，该阶段开始时，试件内部萌生张拉裂纹(图10,11中标识点A)，相应的灰度/纹理特征参数产生突变，标识点A1,A2,A3;伴随着新生裂纹的继续产生和扩展，特征参数再次突变;试件最终破坏时(图10,11中标识点B)，曲线表现为近乎陡立，标识点B1,B2,B3。

结合上述分析可知，灰度特征参数-时间曲线和纹理特征参数-时间曲线变化趋势有一定的相似性，故计算两者的皮尔逊相关系数确定其相关程度。考虑到波动段变化规律不明显，所以只对平稳段和突变段数据进行计算。3种裂隙倾角的砂岩试件计算结果见表2～4,可以看出，90.28%的皮尔逊相关系数绝对值大于0.5，说明绝大部分灰度特征参数与纹理特征参数变化存在显著或高度的相关关系。

表2 β=30°时应变场灰度与纹理特征参数的皮尔逊相关系数Table 2 Pearson’s correlation coefficients of gray and texture feature parameters of strain fields for β=30°

表3 β=45°时应变场灰度与纹理特征参数的皮尔逊相关系数Table 3 Pearson’s correlation coefficients of gray and texture feature parameters of strain fields for β=45°

表4 β=60°时应变场灰度与纹理特征参数的皮尔逊相关系数Table 4 Pearson’s correlation coefficients of gray and texture feature parameters of strain fields for β=60°

3 讨 论

3.1 灰度与纹理特征参数突变机制及前兆点识别

裂隙砂岩试件的灰度/纹理特征参数-时间曲线在突变段表现出明显的规律性和趋势性。根据张拉裂纹出现时(对应于图10和11中的标识点A)，灰度/纹理特征参数-时间曲线不同的趋势变化，可将特征参数划分为:增长型、降低型以及变异不明显型3种，分别对应于标识点A1,A2,A3。其中，灰度标准差、灰度平滑度、灰度三阶矩、灰度熵、纹理熵、纹理对比度以及纹理角二阶距属于增长型特征参数，灰度一致性和纹理相关度属于降低型特征参数，而灰度均值属于变异不明显型特征参数。产生上述突变现象的原因是由于伴随着宏观裂纹的突然出现，试件抗力结构发生改变，在应变场上表现为新生裂纹处产生应变局部化带，导致应变场云图的灰度及纹理特征也随之发生突变。① 灰度标准差、灰度平滑度、灰度三阶矩以及灰度熵与灰度分布的不均匀程度呈正相关关系;纹理角二阶矩、纹理对比度以及纹理熵与纹理的混乱程度呈正相关关系;所以，在裂纹产生和扩展过程中，这些特征参数表现为增长型。② 灰度一致性则表示灰度分布的均匀程度，图像的像素点灰度值越趋于一致，对应的数值最大，故灰度一致性与灰度分布的不均匀程度呈负相关关系;纹理相关性表示图像中相邻像素点间灰度级的相似程度，其与纹理局部变化差异呈负相关关系;所以，在试件变形破裂过程中，这些特征参数表现为降低型。③ 灰度均值只能反映图像的整体灰度等级，对于局部灰度变化并不敏感，研究结果与毛灵涛等[28]对混凝土CT图像的灰度均值变化趋势类似。

从上述分析可以看出，灰度标准差、灰度平滑度、灰度三阶矩、灰度熵、灰度一致性、纹理熵、纹理对比度、纹理角二阶矩以及纹理相关度等特征参数在突变段起始处会出现陡变现象，对应于宏观裂纹起裂，可将其视为前兆点。分别统计3种裂隙砂岩试件前兆点对应的前兆应力、前兆时间、前兆应力/峰值应力比和前兆/最终破坏时间比，见表5。对于同一试件，这些灰度/纹理特征参数前兆点所对应的应力和时间均一致。

表5 基于应变场灰度与纹理特征参数演化规律的前兆点Table 5 Precursory points based on the evolution law of gray and texture feature parameters of strain fields

由表5可知，裂隙砂岩试件的前兆/峰值应力比为0.83～0.97，前兆/最终破坏时间比为0.85～0.95，所以前兆点均在岩体临近最终破坏时出现，对应于影响岩体稳定性的关键破坏事件，可将其作为矿柱失稳的预警信号;而应变局部化带的产生位置则预示着矿柱出现破坏的相应部位。需要指出的是，本文着重于提出一种识别裂隙岩体最终破坏前的异常前兆的方法，表5的结论是基于不同裂隙倾角的砂岩试件室内试验结果，而影响前兆应力和前兆时间的因素有很多，如岩性、矿物成分、裂隙特征等。但有一点可以明确，岩体一旦出现前兆点，表明其物理力学特性已经产生了重大变化。相关试验结论还需要大量不同岩性的裂隙岩石试件来进一步佐证。另外，若将本文的小尺度试验结果应用于岩体工程灾害预警，还需进一步研究前兆点的尺度效应。

3.2 裂隙岩体与完整岩石计算结果的区别

本文还计算加载过程中完整砂岩试件的6种灰度特征参数和4种纹理特征参数，对上述参数同样进行最大最小值归一化处理，与裂隙砂岩试件计算结果进行对比。完整砂岩试件的灰度/纹理特征参数-时间曲线可大致分为:波动段(阶段I)，平稳段(阶段II)以及持续变化段(阶段III)。完整砂岩试件阶段I和II的变化特征与裂隙砂岩试件类似;而后，灰度/纹理特征参数并未出现突变现象，如图12所示，灰度标准差、灰度平滑度、灰度熵、纹理熵、纹理对比度表现为持续增长，灰度一致性则表现为持续降低，其余参数属于变异不明显型，这是因为相比于完整岩石，预制裂隙的存在导致受荷裂隙岩体的变形不连续性更为显著。阶段II和阶段III之间的分界处可视为前兆点，对应的前兆应力、前兆时间、前兆应力/峰值应力比和前兆/最终破坏时间比见表5。

图12 完整砂岩试件灰度与纹理特征参数-时间曲线Fig.12 Gray and texture feature parameters-time curves of intact sandstone specimen

本文提出的6种灰度特征参数和4种纹理特征参数对加载过程中裂隙岩体和完整岩石的应变场变化特征有着不同的敏感性。建议监测的灰度/纹理特征参数见表6。相比于完整岩石，灰度/纹理特征参数对裂隙岩体破裂前兆的识别能力更强。

表6 建议监测的灰度/纹理特征参数Table 6 Recommended monitoring gray/texture feature parameters

4 结 论

(1)裂隙砂岩试件在裂纹起裂、扩展以及贯通过程中表现出不同的应变场灰度特征。从灰度直方图上可以看出，整体从低灰度级向高灰度级方向移动，形态方面则从“矮胖”型转化为“高瘦”型，峰值频率也逐渐增大，最终破坏时达到最大。

(2)根据裂隙砂岩试件应变场灰度/纹理特征参数-时间曲线的变化特征，可将其划分为波动段、平稳段以及突变段。裂纹的萌生与扩展，在应变场上表现为应变局部化带产生与发展，使得应变场灰度及纹理特征发生突变，在应变场灰度/纹理特征参数变化曲线上形成突变段。根据不同的突变段变化趋势，将特征参数划分为增长型、降低型以及变异不明显型。

(3)通过计算裂隙砂岩试件应变场灰度特征参数与纹理特征参数之间的皮尔逊相关性系数，发现90.28%的皮尔逊相关系数绝对值>0.5，说明它们之间存在显著或高度的相关关系。

(4)裂隙砂岩试件应变场灰度/纹理特征参数-时间曲线的第1次突增，此时萌生张拉裂纹，可将其视为前兆点，用于矿柱失稳预警;相应的前兆/最终破坏时间比为0.83～0.97，前兆/峰值应力比为0.85～0.95。相比于完整岩石试件，应变场灰度/纹理特征参数更易于捕捉裂隙岩体的破裂前兆。

通过以上试验研究，应变场灰度与纹理特征参数能够很好地反应裂隙岩体在单轴加载过程中的破裂演化特征，为矿柱失稳灾害预警提供了一种新思路。考虑到地下工程裂隙围岩多处于真三轴应力环境，仍需进一步开展裂隙岩体真三轴应力条件下应变场灰度及纹理特征研究。