基于GMM+SVM的声发射花岗岩裂纹识别研究

黄晓红，董诗琪，李 静，卢 晔

(华北理工大学 人工智能学院，河北 唐山 063200)

在深部矿产资源开采中，岩土工程灾害时有发生，给国家经济和职工生命财产造成巨大损失。岩石失稳破坏是由岩体内部张拉与剪切破裂逐步发育导致的。声发射(AE)作为岩石破裂的伴生现象，可通过其参数特征和波形特征描述岩石内部拉剪破裂演化过程[1-5]，目前虽然有一些裂纹识别技术，但还不能实现裂纹模式识别。岩石的破坏类型存在张拉破坏和剪切破坏，张拉和剪切模式下的裂纹比是评估破坏状态的一个有用参数。分析各加载阶段中，岩石张拉型破裂和剪切型破裂之间的变化关系，对分析岩石破坏机理及工程实践有一定的指导意义。

基于声发射参数特征判别时，较为常用的就是基于AF和RA值的参数法。刘文德统计了花岗岩、灰岩、泥岩和砂岩在单轴压缩和劈裂条件下各加载阶段中张拉破坏和剪切破坏的AE计数的数量及占比，定额给出了各破坏阶段张拉裂纹和剪切裂纹的比例[6]；王林均等通过采用AF值和RA值推导岩石类材料张拉和剪切裂纹模式的合理性及其阈值[7]；何满潮等根据声发射参数RA值与平均频率分布情况，对不同卸荷速率下岩爆试验产生的裂纹类型进行了分析[8]。普遍发现张裂AE信号呈高AF低RA,而剪裂呈低AF和高RA，张裂的声发射波形频率比剪裂的频率高，而且大量的能量在很短时间内释放出来[9-13]。然而传统的基于声发射参数判别法没有明确定义张裂和剪裂的边界，需要着重考虑比值k(AF与RA的比值)。张正虎等针对这一问题，提出了一种基于主频特征统计分析的方法，用于确定声发射参数分析中的裂纹分类标准[14]。

近5年，关于基于机器学习的裂纹识别方法被研究提出。闫召富将幅度和模态分解第四模态熵值作为识别特征，利用高斯混合模型进行分类[15]；Farhidzadeh最早提出基于高斯混合模型的概率方法对钢筋混凝土裂纹进行分类，给出了曲线的分类分割线[16]；周逸飞利用高斯混合模型对加载过程声发射信号参数特征进行深入分析，探索性识别灰岩破裂失稳过程的裂纹模式及前兆特征，揭示灰岩在单轴加载过程中先后主要存在张拉和剪切两种破裂模式[17]。然而在实际情况下，通过高斯混合模型得到的每个簇的聚类是交叉的，无法确定一条将两种裂纹聚类明确分割开的分界线，如果没有适当的方法将数据分成簇，就无法得到反映破坏状态的两个簇中准确的裂纹占比。为了弥补这一问题，笔者提出一种确定声发射参数特征分析中裂纹分类最佳分界线的新方法。通过花岗岩单轴压缩试验获得声发射信号特征，首先进行应力—时间—累计振铃数及应力—时间—累计能量全过程分析，通过JCMS-IIIB5706[18]裂纹分类法得到AF和RA两个参数，再利用GMM对参数分布进行概率密度处理，最后通过基于支持向量机(SVM)算法创建的超平面分离数据中的交叉簇，最终验证了花岗岩破裂模式与GMM+SVM识别方法存在相对应关系，这有助于提高声发射源机理、损伤过程分析的准确性和效率，为定量制订岩体灾害预警方案提供一系列可靠的监测阈值。

1 花岗岩声发射试验

1.1 花岗岩试样制备

试验所采用的花岗岩来自辽宁沈阳,将花岗岩加工成标准长方体，尺寸为100 mm×50 mm×50 mm，两端的不平行度误差均小于0.02 mm,共制备7块试样,编号为GUD-1～GUD-7。

1.2 试验设备与方法

本次试验采用的加载系统是TAW-3000型常规三轴实验机，该仪器精度高且性能稳定。为了确保样品和加载表面之间完全接触，首先对样品加载到 2 kN，然后以0.15 mm/min的速率增加，直至岩样破坏。在试验过程中，声发射系统使用全信息多通道声发射监测系统DS5-16，传感器频率为120 kHz，阈值电压为0.5 V，放大器增益为 40 dB。试验设备装置如图1所示。此外，在样品的左侧和右侧分别布置了AE探头(见图2)，并用气球固定，在岩样与AE探头之间涂抹一层凡士林以便达到最佳耦合。在花岗岩样品断裂过程中收集AE信息。声发射系统的幅度阈值设置为45 dB，以消除环境噪声对声发射信号的影响。

图1 试验设备装置

图2 花岗岩试样

1.3 花岗岩声发射特征

首先分析花岗岩试样累计声发射计数及累计声发射能量的变化规律。限于篇幅，选取典型试样进行分析，如图3所示。

从图3(a)中可以看出，花岗岩试样在到达峰值应力前声发射现象较少，而在接近峰值应力时声发射大量发生。从图3(b)中同样发现，试样在峰值应力前AE能量释放很少,然而到达峰值之后大量AE能量释放。原因是花岗岩构成颗粒之间的结合相对较牢，只有当应力超过颗粒的胶结强度时才会形成裂纹。而且，花岗岩结构较均匀，不同部位的断裂时间相差不大。因此，当应力达到试样的峰值强度时，将产生大量声发射现象，同时释放大量能量。

图3 GUD-7试样应力—时间—累计声发射计数及能量曲线

2 基于声发射参数的花岗岩裂纹分类

对累计声发射计数、累计声发射能量、AF值、RA值等声发射参数进行分析和运用。JCMS-IIIB5706建议使用声发射波的平均频率(AF)和RA值，具体计算公式如下：

此前很多学者利用JCMS-IIIB5706方法进行裂纹分类，此方法是通过AF和RA值之间的关系，将裂纹分为张拉和剪切2种类型，如图4所示。

图4 JCMS-IIIB5706方法裂纹分类

由于加载前期声发射活动较少，故选取GUD-7加载应力为0.25P、0.50P、0.60P、0.70P、0.80P、0.90P和1.00P(P为峰值强度)对应的数据进行研究，得出单轴压缩条件下花岗岩试样各阶段的AF-RA散点分布图，如图5所示。

图5 单轴压缩条件下花岗岩各阶段AF-RA散点分布图

从图5中可以看出，花岗岩试样张拉裂纹区域散点在各阶段内呈现出高AF值、低RA值的特点，而剪切裂纹区域内散点呈现高RA值对应低AF值的特点。可以得出在该应力加载水平条件下，岩石破裂产生张拉型裂纹和剪切型裂纹，整个加载过程先以张拉裂纹为主，剪切裂纹后期增量较高。RA值主要集中分布在 0～5 ms/V，平均频率AF值主要集中分布在 0～600 kHz。

3 基于GMM+SVM花岗岩裂纹分类

3.1 GMM+SVM算法原理

GMM是一种有限混合模型无监督学习算法。假设所有数据点都是在有限数量的高斯分布和未知参数的混合中生成的，GMM算法可利用每个(高斯)分量的加权因子将数据分类为不同的簇[19]。式(1)为用于输入特征向量x的GMM的广义方程：

(1)

(2)

式中：k为裂纹类型的总数；ωi为材料中观察到的破裂机制的相对密度；x为所选波形参数的数据集；k*1 是均值向量μi的维数，k*k是裂纹类型协方差矩阵Σi的维数。

由于均值和协方差未知，因此模型的下一个目标是找到权重因子(ωi)、均值向量(μ)和协方差矩阵(Σi)的最佳值，使得高斯混合密度模型p(x)可以最佳地拟合所有的特征向量(x)。为了找到最优解，使用了最大似然估计方法。最大似然法最大化给定x数据集的对数似然，以获得最佳GMM参数。对于一组测试训练样本T,T={x1,x2,…,xT}，假设向量之间是独立的，GMM似然方程可以写成公式(3)的形式：

(3)

而后使用期望最大化(EM)的迭代算法[20]，以获得参数值。EM算法基于新模型参数，从之前的参数获得的结果应遵循公式(4)：

p(X|ωi,μi,Σi)m+1>p(X|ωi,μi,Σi)

(4)

第一步，称为期望步骤(或E步骤)，包括在给定模型参数的情况下，计算由每个数据点的分量i得出的训练数据点的期望值：

(5)

式中m为迭代次数。

第二步，称为最大化步骤(或称M步)，包括最大化在E步骤中计算出的模型参数期望值，也包括更新后的值：

(6)

(7)

(8)

(9)

(10)

基于后验概率p(xt)，根据公式(11)完成聚类(C)分配。其中，Findmax找出与概率向量元素的最大值相对应的聚类：

C=Findmaxp(xt)

(11)

在使用GMM将数据分类到不同的聚类之后，每个数据点获得与该聚类对应的唯一标签。为了在张拉和剪切簇中定量分布声发射数据，标记数据是前提条件。尽管GMM可以分配数据标签，但是由于物理过程的原因，始终伴随着群集相交的情况。因此，在超平面的实际方程中总是出现混淆，在张拉和剪切模式的比率上也是如此，这对于识别岩石破裂表征状态判别是十分重要的。使用聚类标记数据的线性SVM来计算最佳的超平面。根据聚类的数量，线性SVM重复k次，在每组计算中，都会计算超平面以将每个聚类分开。

给定数据点(x1，C1)，…，(xn，Cn)，其中Ci的值为1或 -1 (如果数据点属于聚类I则值为1，否则值为-1)，每个数据点表示一个向量xi。目标是定义一个具有最大边缘值的分离超平面，该超平面将数据分成两部分，Ci的值为1或-1。公式(12)给出了满足xi的超平面的一般方程：

w·x-b=0

(12)

式中：w为超平面的法线向量;b为一个实数，决定了超平面与原点的距离。

为了得到最佳的超平面，需要使聚类的两个(线性)边界与其余数据(即不属于第i个群集的数据)之间的距离最大化。这些线性边界被称为支持超平面，定义支持超平面的向量x1被称为支持向量xSV。

这两个支持超平面方程可以写成：

w·xSV(Ci=1)-b=1

(13)

w·xSV(Ci=-1)-b=-1

(14)

上述超平面的定义表明，所有x必须位于支持超平面的正确侧。从数学上，该条件可以写为等式(14)。w被估计为最小值，等式(15)成立：

Ci(w·xi-b)≥1,(∀,1≤i≤n)

(15)

所提出的GMM+SVM算法流程图如图6所示。

图6 GMM+SVM算法流程图

3.2 基于GMM的参数分布聚类

不同应力时的花岗岩试样GMM特征向量概率分布如图7所示，其中红色区域代表较高的概率密度，蓝色区域代表较小的概率密度。

图7 花岗岩在不同应力时的GMM特征向量概率分布

由图7可见，花岗岩在整个加载过程中张拉裂纹占比较高，张拉聚类的椭圆颜色较深呈高概率分布，剪切聚类的形状颜色较浅呈低概率分布，在前期加载步骤中(0～0.50)P，声发射信号主要集中分布在平均频率纵轴的一侧，大部分声发射信号的RA值小于2 ms/V，平均频率主要分布在300～600 kHz。此时产生了较高的AF值和较小的RA值，这是因为微裂纹在加载初期载荷较小的情况下一般是成核的。加载到中间步骤(0.50～0.70)P时为张拉至剪切的过渡阶段，此时由于加载载荷不断增加，AF值通常较高，RA值在较小范围内变化，张拉聚类高概率区域开始向剪切聚类区域转移。张拉和剪切两个分离的簇逐渐开始形成交叉。最后阶段(0.80～1.00)P声发射AF值开始快速下降，声发射信号逐渐集中分布在RA值横轴的一侧，呈条带状，其中声发射信号RA值主要分布在0～8 ms/V，平均频率集中在0～400 kHz。

由于GMM自身无法完成数据转换和更改数据分布，且因其对角矩阵和组件独立性假设，模型通常受到限制，并且模式匹配的准确性有时不能满足预期的要求。基于GMM的聚类算法可以预测整个过程的类似密度特征，如图8(a)所示。

图8 花岗岩的GMM特征概率密度聚类及其分离

但由于从GMM获得的密度函数显示聚类是相交的，张拉和剪切聚类的可能视觉选择如图8(b)所示。对于同一组数据，可选择3个(或多个)可能的主观超平面，其在分离数据的效果方面似乎同样好。对于不同的超平面，对裂纹模式分类至关重要的两个聚类中的事件比率将是不同的。按照所提出 GMM + SVM算法，SVM将计算唯一的最佳超平面来分离相交的聚类。

3.3 基于GMM+SVM的花岗岩裂纹演化分析

在达到峰值应力之前，花岗岩单轴压缩的破坏机制可分为3种不同的状态:初始压密、弹性和塑性破坏3个阶段。不同状态下的计算超平面(separating hyperplane)如图9所示。

图9 花岗岩不同应力时间下GMM+SVM超平面分界线

从图9中可以看出，花岗岩由初始压密阶段到线弹性阶段(0～0.50)P只有少量声发射事件发生且张拉裂纹(Tensile)占有主导地位，此时，岩石样品中没有明显的破裂，声发射信号是由样品本身中的微孔和裂缝的压实引起的，此阶段是岩石破裂的平静期；弹性阶段(0.50～0.80)P声发射事件数量逐渐增加，但仍以张拉裂纹为主，AF值大部分集中在350～600 kHz，RA值为0～4.5 ms/V，此时是经过一段时间的加载，轴向应力增大，花岗岩的内部应力集中区域已发生轻微断裂，但在实际施工领域，这一现象不易被观测到；塑性破坏阶段(0.80～1.00)P剪切裂纹(Shear)产生大量声发射信号，AF值集中在 150～400 kHz，相比之前阶段，AF值大幅度下降，然后RA值开始上升，RA值集中在2～6 ms/V，大量声发射事件产生，出现大量的剪切裂纹，这是因为随着轴向应力达到最大值，花岗岩内部破坏严重，此时伴有大能量释放。

随着应力加载过程各阶段张拉和剪切裂纹占比及SVM超平面的数学表达式结果如表1所示，可以得到剪切裂纹所占比例在80%～90%峰值应力阶段出现突变。

表1 各应力区间花岗岩张拉、剪切裂纹占比及SVM超平面分界线

岩石试样GUD-1～GUD-7在(0.80～0.90)P应力区间剪切裂纹占比分别为43.1%、42.9%、40.9%、43.7%、44.2%、45.3%、44.7%，为确保安全性，可将最低值作为可靠的监测阈值，由此可以得到当剪切裂纹占比超过40.9%时可触发预警警报。最终验证了所提出的基于GMM+SVM算法获得的裂纹分析结果与JCMS-IIIB5706 提出的参数分析方法结果基本一致。

4 结语

1)通过花岗岩单轴加载试验得到花岗岩声发射AF-RA散点分布图，得出在应力加载时间内花岗岩试样的AF-RA散点呈现低AF值与高RA值对应、低 RA值与高AF值对应的特征，这是张拉型裂纹和剪切型裂纹破裂机制区分的重要表征。AF-RA的数据点在加载初期靠拢AF轴，后期逐渐向RA轴靠拢，即数据点由张拉裂纹区域向剪切裂纹区域转移。

2)通过高斯混合模型(GMM)对单轴压缩条件下花岗岩破裂裂纹的转变规律进行研究分析，相比于花岗岩试样的AF-RA散点图，GMM特征向量概率分布能够更直观地看出，随着应力加载的不断进行，AF值逐渐减小，由前期加载步骤中的300～600 kHz 下降至后期的0～400 kHz，而RA值逐步增大，由前期的0～2 ms/V上升至0～8 ms/V。此现象表明了花岗岩破裂产生裂纹的过程遵循以张拉裂纹为主，剪切裂纹不断增加的破裂机制。

3)通过使用GMM+SVM算法获得最佳分离聚类的超平面，简化了裂纹的分类，因为GMM聚类没有提供拉伸和剪切区域之间的精确区分，消除了GMM密度曲线的棘手使用。得出了花岗岩张剪裂纹在各应力加载阶段的占比。得出剪切裂纹所占比例的最大值出现在(0.80～0.90)P的阶段，该阶段为岩样加载进程的非稳定扩展阶段的中后期，该占比值可用作岩石破坏的有效监测预警阈值，为岩石失稳破坏预测提供依据。

4)需要指出的是，本文是基于室内花岗岩单轴加载试验进行的研究，岩样尺寸较小，与实际工程现场的各现象和条件会存在一定的差别，同时针对工程现场监测也需要更高灵敏的仪器。未来还需要结合工程现场进一步研究，以期获得普适性的规律。本文所提出的研究方法仍需在大数据条件下进行进一步验证。