基于小波散射分解变换的煤矿微震信号智能识别

樊 鑫，程建远，王云宏，栗 升，段建华，王 盼,

(1.西安科技大学 地质与环境学院，陕西 西安 710054；2.中国煤炭科工集团 西安研究院有限公司， 陕西 西安 710077)

微震监测技术是一种利用煤岩破裂产生的震动信息来研究煤岩结构稳定性的实时、动态、连续监测方法。随着我国煤矿开采深度的增加，微震监测技术在煤矿安全高效生产中得到广泛应用，并被写入《煤矿防治水细则》。通过在监测区域对微震信号进行采集、处理和分析，可以研究煤岩体变形破坏的活动规律及其内部应力分布，并对突水等煤岩动力灾害进行监测预警。但是，由于微震监测系统经常工作在高噪声环境下，各道采集的数据包含大量噪声，严重干扰到微震事件拾取和震源定位等后续处理工作。因此，将微震有效信号与噪声加以分离并进行特征提取具有重要意义。

微震信号是一种具有随机性的时变非平稳信号，且具有震源信号和干扰噪声多样的特点。近年来，相关研究人员对微震信号的分析做了很多努力。Fourier 变换是从全局时域提取特征的一种信号分析方法，无法提取局部特征，对于非平稳信号的分析有很大的局限性；短时傅里叶变换方法易于实现，但分辨率单一；加窗Fourier分析同样无法描述微震信号的局部特征。以小波变换为代表的时-频分析方法是基于时间和频率参照系的信号分析方法，被证明可用于描述非平稳信号。由于小波变换具有局部形变稳定性和多尺度性的特点，它可以有效提取信号的局部特征，但因其会随时移变化，容易遗漏信号特征，亦不适于时变非平稳信号的分析以及特征向量的构建。

2010年，ANDEN等提出了一种基于小波变换改进的时-频分析方法——小波散射变换(Wavelet Scattering Transform)。它具有平移不变性、局部形变稳定性和信息完整性等特点，解决了小波变换随时移变化的缺陷，并在相关领域得到了有效应用。2013年，HIRN和BRUNA等通过构建小波散射分解变换网络，实现对音频信号、手写文字和图像纹理的识别；2014年，ANDEN和JOAKIM等分别在经典音乐数据集GTZAN和语音通话数据集TIMIT中，通过小波散射分解变换网络提取到有效特征信息，取得良好的分类效果；同年，CHUDACEK等将该方法应用到心律失常数据的分析中。2018年，WIATOWSKI等基于小波散射的性质，通过严格的数学推导，证明了其优越性并进行了推广，在不同的小波框架下取得了良好效果。

继MALLAT之后，国内学者围绕小波分析进行了相应的研究，不仅在文字、图像等二维数据的识别中取得良好效果，还在音频、心律等时变非平稳信号的识别中得到有效验证。其中，在二维数据的识别中，吴华娟等通过小波散射变换网络提取每个子图像块的能量分布特征，对各阶散射特征进行分类，实现了图像的无监督纹理分割；伍家松等利用小波散射变换网络提取特征并形成训练数据集，分析了在不同色彩空间进行图像分类的性能；文介华等以小波散射变换网络所得特征系数的均值和方差作为新的特征，实现了大规模图像识别，并起到有效的降维作用；WANG等利用小波散射变换网络提取合成孔径雷达图像的特征，有效识别了移动和固定目标。

LI等提出了心音信号分类算法，利用小波散射变换网络获取到心音信号特征，能够有效表达信号对应的特征信息，进而得到信号的特征矩阵，用于支持向量机分类；罗雨帆改进了经典小波散射变换网络，构建了分数阶散射变换网络，将频域小波变换扩展到分数域，并在一维音频信号分析及特征提取中取得了良好的应用效果；许周乐将小波散射变换网络作为混合网络特征提取的一部分，提取机械振动故障信号特征，从而实现了轴承故障诊断。

前人利用小波散射分解技术在二维图像识别取得了一定的成果，同时还在一维时序非平稳信号等领域得到了有效应用，如音频信号。作为一种典型的时序非平稳信号，微震信号中携带大量地质信息，对其进行特征提取不仅可以提高震源定位的精度，而且能够实现震源的准确定位。在页岩气开发领域，TAN等提取了微震事件的时长、频域和统计特征，并基于主成分分析法构建了微震事件分类识别模型，识别精度达到90%以上，取得良好效果。矿山微震信号识别领域的相关研究目前处于发展阶段。张少泉等指出要以矿山地震学理论为基础并选择区分度较高的判据，进行矿山震动信号识别。朱权洁等利用微震信号的小波包能量特征和分形特征构建了相应的特征向量，通过SVM模型对爆破、机械振动和岩石破裂波形进行分类。

与前述研究不同，笔者针对煤矿微震波形，提取信号特征、降低原始信号维数、构成特征矩阵，进而建立基于小波散射分解和SVM分类的识别模型，为实现微震事件的智能识别提供了一种新的技术思路。

1 小波散射分解变换基本原理

小波变换是时变非平稳信号分析的有效工具之一，因其具有尺度可变性和多分辨率的特点，能同时描述信号的时域和频域特征，因此对信号的局部分析具有良好效果。对连续有限时间内的信号，小波变换定义为

(1)

实际采集的微震信号通常会受到很多干扰，即使整体并没有发生质变，局部的变化也会使提取的信号特征受到干扰，从而影响信号的分析和识别。因此，需要一种兼具平移不变性和局部形变稳定性的信号分析和特征提取方法。

对小波变换进行取模运算，去掉所有小波散射系数的复相位，则可以得到算子||，与输入信号做卷积可得到小波模变换算子：

||=(*,|*|)

(2)

其中，为低通滤波器；为高频小波；()=*称为关于信号的局部平移不变描述符，即散射系数，该描述符是对输入信号低通滤波，具有平移不变性，提取输入信号的低频信息，去除了所有的高频信息；而这些高频信息由模变换()=|*()|恢复，该变换表示在尺度上的高频信息，通过对非线性小波变换取模，具有形变稳定性。故小波散射变换0阶的低频信息(散射系数)和高频信息分别为

(3)

将0阶高频信息部分()作为一阶散射变换的输入，可得

|||*1|=(|*|*,||*|*|)

(4)

则一阶散射系数

()=|*|*

(5)

以此类推，重复以上迭代过程即可得任意阶的散射系数。

对于任意的≥1，信号的小波模变换卷积可表示为：

()=|||*|*…|*|

(6)

将()作为下一阶的输入进行低通滤波，得到阶散射系数：

()=|||*|*…|*|*=()*

(7)

将|+1|应用到()可同时计算得到()和+1()，即

|+1|()=((),+1())

(8)

通过初始化()=可定义散射分解的最高阶，当0≤≤且1≤≤时，对式(1)～(8)进行迭代运算。

最终由0≤≤上的散射系数构成特征向量：()={(),(),…,()}，即：()=()，(0≤≤)。

综上，小波散射变换的过程可描述为：在小波模算子||上进行散射变换迭代，卷积计算次的小波模变换()并输出经低通滤波后的散射系数()(图1)。

图1 小波散射分解结构Fig.1 Wavelet scattering decomposition structure

2 微震信号的小波散射分解系数分析

2.1 微震事件的小波散射分解系数

图2(a)为使用中煤科工集团西安研究院有限公司研发的KJ959微震监测系统于西部某煤矿采集的底板破坏深度微震事件信号，采用频率响应范围为10～1 000 Hz的单分量检波器，采样率为2 kHz。

根据前文所述小波散射分析方法进行处理，得到每一次散射分解变换产生的散射分解系数矩阵维数(表1)，其各阶散射系数矩阵由不同的散射行向量组成，每个行向量包含4个散射系数。图2(b)为微震事件信号波形及其对应的经4次分解所得的0～4阶散射系数数值(从各阶系数矩阵中随机抽取一行)分布。

图2 微震事件波形及其0～4阶小波散射分解系数分布Fig.2 Waveform of a microseimic event signal and its 0-4th order wavelet scattering decomposition coefficients

表1 微震事件小波散射分解系数矩阵维数Table 1 Matrix dimensions of wavelet scattering decomposition coefficients of microseismic events

图2中，微震信号具有瞬时突变的特点，其能量在各阶散射分解变换后的分布不一致，由各阶散射系数的数值可以看出：微震信号的能量主要由1～3阶散射系数体现，且与微震事件具有相同的突变特点，而2阶散射系数具有最大的能量分布，1阶和3阶系数的能量分布较小。

图3为微震事件信号经4次小波散射分解之后得到的4张时频谱，对比后可以发现：信号能量在微震事件发生的时刻集中出现，并随着散射分解次数的增加，低频部分的能量逐渐凸显(表2)。

图3 微震事件信号小波散射变换各阶时频谱Fig.3 Scattergram-scalogram coefficients filter banks of wavelet scattering decomposition for a microseismic event signal

表2 微震事件4次小波散射分解对应的频率范围Table 2 Frequency range of 4 times wavelet scattering decomposition of microseismic events

综上，小波散射分解变换可以有效分析微震事件信号的能量分布和时频特征，并由计算得到的各阶散射分解系数表征，进而形成散射系数矩阵。

2.2 噪声信号的小波散射分解系数

利用小波散射分析方法，可得到微震噪声信号波形及其对应的、经4次分解所得的0～4阶散射系数数值(从各阶系数矩阵中随机抽取一行)(图4)。表3为每一次散射分解变换产生的散射分解系数矩阵维数，其各阶散射系数矩阵由不同的散射行向量组成，每个行向量包含4个散射系数。由于截取的噪声信号长度与事件信号长度一致，故对应得到相同维数的散射系数矩阵。

图4 微震噪声波形及其0～4阶小波散射分解系数分布Fig.4 Waveform of a noise signal and its 0-4th order wavelet scattering decomposition coefficients

从图4(b)噪声信号的各阶小波散射分解系数可以看出：能量在各阶散射分解变换后的分布大致平均，主要由1～2阶散射系数体现，并无明显的突变特征及能量聚集的体现，符合噪声事件能量分散、成分杂乱的特点，与具有瞬时突变特点的微震事件信号的小波散射分解系数有明显区别(图2)。

图5为微震噪声信号经4次小波散射分解之后得到的4张时频谱，对比后可以发现：信号能量在噪声信号中分散分布，并未在某一时刻集中出现，且包含较宽频带的信息。

图5 微震噪声信号小波散射变换各阶时频谱Fig.5 Scattergram-scalogram coefficients filter banks of wavelet scattering decomposition for a noise signal

由此可见，通过小波散射分解变换技术可以获得微震事件信号和噪声信号的不同特征，进而对2类信号加以描述，为微震信号的识别提供了一种有效的特征提取方法。

2.3 微震事件和噪声分解系数对比

2.1节和2.2节利用小波散射分解技术分别对微震事件和噪声信号进行了分析，本节将2类信号的分析结果进行对比，进一步验证相关结论。

表3 微震噪声小波散射分解系数矩阵维数Table 3 Matrix dimensions of wavelet scattering decomposition coefficients of noise events

图6分别从微震事件(Signal of microseismic events，简称“sm”)和噪声信号(Signal of noise events，简称“sn”)中任意选取3条信号进行小波散射分解变换，得到的0～3阶散射系数分布(从每条信号的各阶系数矩阵中随机抽取一行)，其中微震事件信号的散射系数由“实心圆”表示、点划线连接；噪声信号的散射系数由“菱形”表示、虚线连接。由图6可知，2类信号的散射系数自一阶散射分解变换开始，逐渐表现出各自信号的差异特征，其中微震事件信号的散射系数的变化趋势表现为瞬时突变，而噪声信号的散射系数变化趋势是杂乱无序的且数值远小于微震事件信号的散射系数，符合前述信号分析结果，在数值大小和能量表示方面都有明显差异，可以进一步形成特征矩阵，用于机器学习分类模型的训练。

图6 微震事件与噪声信号小波散射分解系数对比Fig.6 Comparison of microseismic event and noise signals wavelet scattering decomposition coefficients

3 小波散射分解系数特征矩阵

原始微震监测数据的高维空间中包含冗余信息，会在模型训练中引入干扰和误差，影响训练精度。而通过对数据进行小波散射分解变换，既可以实现对监测大数据的降维，提高计算效率；又可以提取数据本质特征，减少冗余信息。小波散射分解系数特征矩阵是维数为×的特征矩阵，其中行数为总散射路径数，列数为散射系数分量。计算小波散射分解系数并形成最优维数(×)的特征矩阵，需要构建合适的小波散射分解结构。由小波散射分解的原理和性质可知，时不变尺度(即进行小波散射变换时，所取信号片段的采样时长)、小波散射分解变换的次数以及质量因子(Quality factor，即每个小波滤波器组中的小波数目，它随着小波变换次数的增加而递减，且满足1≤<32)是重要参数，以下实验将说明上述3个不同参数对小波散射分解结构性能、算法效率及最终输出特征矩阵的影响。

考虑到微震信号特征和算法的复杂性，通过计算机程序对已有数据进行分析。已知实验所使用微震信号总长度为7 000个采样点，共包含54条微震事件信号和54条噪声信号。

3.1 时间不变尺度的影响分析

时间不变尺度用以表征小波散射变换的平移不变性，其大小以采样时间计算。构建一个小波散射分解结构，并设定小波散射分解采样频率为1 024 Hz，小波变换次数为3，质量因子为3，2，1。表4为不同时间不变尺度时，得到的相应散射特征训练矩阵、测试矩阵及总特征矩阵的构成与对应的特征矩阵构成所需运行时间。

表4 不同时间不变尺度对应的散射特征训练矩阵、 测试矩阵及总特征矩阵维数Table 4 Dimensions of different time invariance scales for the train feature matrix,test feature matrix and whole feature matrix of a microseismic signal

由以上分析可知，当小波散射采样频率为确定值时，时不变尺度越小，小波散射分解系数特征矩阵维数越高，模型运行时间越长；时不变尺度越大，特征矩阵维数越低，模型运行时间越低。而当散射分解采样频率与时间不变尺度之积近似等于采样点总数时，可得到最佳维数的特征矩阵，此时所需的计算时间较短，运算效率较高。

3.2 小波散射分解变换次数的影响

小波散射分解变换次数，即小波滤波器组的个数，对于形成特征矩阵产生的不同影响，将在以下实例中说明。构建一个小波散射分解结构，设定时间不变尺度为6 s，采样频率为1 024 Hz，小波变换次数为1时，质量因子为3；小波变换次数为2时，质量因子为3和2；小波变换次数为3时，质量因子为3，2，1。表5为取不同小波散射变换次数时，得到的相应散射特征训练矩阵、测试矩阵及总特征矩阵，表6为对应的散射变换及特征矩阵构成所需运行时间。

由表5，6可知，小波散射分解变换次数越多，计算得到的散射系数越多，构成的特征矩阵维数越大，所需运算时间越长。结合文中第2部分对微震事件和噪声信号的分析，微震事件的主要能量分布在小波散射分解变换的1～3阶散射系数，且集中在2阶散射系数；3次小波散射分解变换已得到微震事件信号的主要能量和频率分布，同时有效表征与噪声信号的特征差异，并且程序运算效率高。因此，宜采用3次小波散射分解变换进行分析，并构成特征矩阵。

表5 不同变换次数对应的散射特征训练矩阵、 测试矩阵及总特征矩阵维数Table 5 Dimensions of different transform times for the train feature matrix,test feature matrix and whole feature matrix of a microseismic signal

表6 不同变换次数对应的散射变换及特征矩阵构成所需运行时间Table 6 Function time of different transform times of wavelet scattering transform and forming its feature matrix

3.3 质量因子的影响分析

构建一个小波散射分解结构，设定时间不变尺度为6，采样频率为1 024 Hz，小波变换次数为1时，质量因子分别取8，6，4，2，得到的相应散射特征训练矩阵、测试矩阵及总特征矩阵，表7，8为对应的散射变换及特征矩阵构成所需运行时间。

表7 不同质量因子对应的散射特征训练矩阵、 测试矩阵及总特征矩阵维数Table 7 Dimensions of different quality factors for the train feature matrix,test feature matrix and whole feature matrix of a microseismic signal

表8 不同质量因子对应的散射变换及特征 矩阵构成所需运行时间Table 8 Function time of different quality factors of wavelet scattering transform and forming its feature matrix

由表7，8可知：当时不变尺度和分解次数不变时，质量因子越大，形成的特征矩阵列数越多(图7)，矩阵维数越高，运算时间越长。因此，在满足分类算法需求的基础上，为提高运算效率，质量因子的选择不宜过大。

图7 质量因子对特征矩阵维数的影响Fig.7 Influence of the quality factor to the dimension of the feature matrix

因此，以满足运算时效性为前提，在尽可能保留不同信号主要特征信息的情况下，特征矩阵的维数不宜过大；在训练过程中可以根据时间不变尺度、散射分解次数和质量因子等主要参数的影响，调控特征矩阵的维数，便于进一步的训练和测试，以得到最优的智能分类识别模型。

4 基于SVM分类的事件智能识别

近年来，监督机器学习中的支持向量机SVM(Support Vector Machine)分类算法被广泛应用于图像和文本识别、信号处理等领域，解决非线性、高维、小样本的模式识别问题。在本文中，笔者将微震信号的智能识别作为一种非线性二分类问题处理。利用前文由小波散射系数构成的特征矩阵，基于SVM分类算法，构建微震事件自动识别模型。

4.1 SVM分类模型的建立

基于SVM算法的模式识别原理是通过对已标记标签的信号样本进行训练，得到分类模型，进而对未知信号样本进行分类处理。其中，构建特征矩阵和分类器的训练是该算法的核心。图8为SVM分类模型的训练及优化流程。

图8 SVM分类模型训练优化流程Fig.8 Workflow of training and optimizing SVM classification model

(1)样本选择。从已有数据中选取360组微震信号样本，微震事件和噪声信号各180组。随机选择70%的样本构成训练数据集，其余30%为测试数据集。分别对两类信号进行标签化，事件信号标记为“M”，噪声信号标记为“N”。

(2)特征提取。通过设定参数时不变尺度、变换次数和质量因子，计算样本信号的小波散射系数构建特征矩阵。

(3)模型训练。SVM分类器选用多项式核函数：

(,)=(1+′)

(9)

式中,′为的转置；为正整数，本文选择=2。

4.2 识别结果及讨论

利用前述特征提取方法得到训练数据集特征矩阵，训练得到微震事件SVM分类模型，然后将训练得到的分类模型对测试数据集进行分类预测。利用小波散射系数特征矩阵识别的结果见表9。从表9可以看出，基于SVM分类的微震事件智能识别算法准确率较高，能够有效、准确地对微震信号分类。

表9 SVM分类模型识别结果Table 9 Recognition results of SVM classification model

随着微震监测技术被广泛应用于煤矿地质灾害监测预警等领域，急剧增加的海量监测数据面临着传输效率低、区分难度大以及现场干扰因素多等问题。因此，如何准确划分信号类型是微震监测大数据处理效果好坏的关键。

综上，该方法可以通过智能识别代替手动识别，大大提高识别效率、减小人工工作量，通过有效识别微震事件，进一步提高震源定位精度，同时也为实现震源自动定位奠定了基础。

5 结 论

(1)基于微震信号低宽频特性、不可预知性和突发瞬态性的特点，采用小波散射分解对微震事件和噪声信号进行特征提取，构成相应的特征矩阵，该方法可以提取出具有明显区分度的微震事件和噪声的信号特征，且效果良好。

(2)通过对比分析小波散射分解结构的3个主要参数：时不变尺度、散射分解次数和质量因子对特征矩阵维数的影响，以及比较选取不同参数时的运算效率，可以发现：在特征矩阵维数相同时，为有效表征微震事件信号与噪声信号的特征差异，并且提高程序运算效率，宜采用3次小波散射分解变换进行分析，并构成特征矩阵，且各阶小波散射变换的质量因子选择不宜过大。

(3)在SVM分类模型的微震事件智能识别过程中，可以选择二次多项式核函数。该方法适用于长时间、大范围的微震监测产生的大数据信号特征提取，有效避免了冗余数据的干扰，提高了数据处理的效率。针对不同监测区域，尤其是低信噪比环境下采集的数据，如何选取合适的主要参数，使得信号特征计算效率最高，需要后续研究进一步验证、完善。