CN112649512B

CN112649512B - 一种岩体声发射初至自适应识别方法

Info

Publication number: CN112649512B
Application number: CN202011475811.9A
Authority: CN
Inventors: 谢学斌; 叶永飞; 王小平; 刘涛; 姜伟; 张欢; 赵伟; 李少乾; 高山
Original assignee: Central South University
Current assignee: Central South University
Priority date: 2020-12-14
Filing date: 2020-12-14
Publication date: 2022-05-13
Anticipated expiration: 2040-12-14
Also published as: CN112649512A

Abstract

本发明涉及声发射初至拾取方法，公开了一种岩体声发射初至自适应识别方法，包括以下步骤：S1、获得监测信号的包络曲线；S2、确定所述监测信号的包络周期；S3、对所述包络曲线进行去除周期变化和异常值后获得趋势曲线；S4、通过判定曲线对所述趋势曲线进行判断，以自动识别声发射初至点。本发明能自动判别监测信号有无声发射，并自动识别声发射初至，且对初至拾取的适应性强、准确性高。

Description

一种岩体声发射初至自适应识别方法

技术领域

本发明涉及声发射初至拾取方法，具体地，涉及一种岩体声发射初至自适应识别方法。

背景技术

近年来，岩体发射监测技术在地压监测中得到了广泛的运用，相应的声发射监测信息分析主要的研究内容包括波形识别、初至拾取、声源定位、参数分析，波形反演等，其中初至拾取是声源定位和地压灾害预警、预报最为基础和重要的一步。在声发射实时监测中数据量庞大复杂，人工识别难以适用，而运用较多的阀值法缺乏灵活性，准确率较低。为此，各种根据信号与噪声间差异特征而提出的初至识别方法相继产生，根据各方法的特点可归纳为两类：一是根据监测信号的波形特征变换识别初至，如通过波形振幅、能量、偏振、分形维数等特征变换判断初至，比较经典的方法有：STA/LTA法、AIC法高阶统计法、MER法；二是联合判别分析，该类方法先通过某种方法粗略拾取初至，再结合其他方法确定精确初至，如STA/LTA与AIC法联合确定初至。现今存在的各种初至拾取方法都具有一定的准确性和适用性，但它们不可避免的存在对参数设置的依赖，当参数设置不合理时，这些方法的准确率较低甚至还会出现漏选和误选。

不同于实验室，在地下矿山实际监测工程实践中背景噪音复杂多变，强度大，对于这种低信噪比噪音环境，目前多采用降噪后进行初至识别的策略，或先大致识别初至然后对该区域进行降噪再进行初至识别的策略，这两种策略均涉及到降噪，但不管是采用何种降噪处理均会造成有用信号的损失，致使拾取精度较差。特征函数是众多识别方法的核心，特征函数选择的好坏直接决定了初至拾取的准确性，尤其是在低信噪比时，特征函数起到简化信号数据波动变化的干扰，放大初至特征的作用。但现今的众多方法中，特征函数依赖预先设置，不能根据监测信号本身进行自适应改变，缺乏足够的稳健性，容易在低信噪比或者高频噪音环境中造成拾取结果的错误。为此在声发射监测中亟需一种无需预先设置参数且能自适应识别初至的岩体声发射初至拾取方法。

发明内容

本发明所要解决的技术问题是提供一种岩体声发射初至自适应识别方法，能自动判别监测信号有无声发射，并自动识别声发射初至，且对初至拾取的适应性强、准确性高。

为了解决上述技术问题，本发明提供一种岩体声发射初至自适应识别方法，包括以下步骤：

S1、获得监测信号的包络曲线；

S2、确定所述监测信号的包络周期；

S3、对所述包络曲线进行去除周期变化和异常值后获得趋势曲线；

S4、通过判定曲线对所述趋势曲线进行判断，以自动识别声发射初至点。

优选地，所述步骤S1中，对所述监测信号进行希尔伯特变换来获得所述包络曲线。

优选地，所述步骤S2中，采用周期图法和自相关系数相结合的方法来确定所述监测信号的包络周期。

优选地，所述监测信号的包络周期的确定方法为：

(1)通过周期图法计算所述监测信号的功率谱，并选择前五个功率谱对应的频率f，其中，功率谱的计算公式为：

式中，x(t)为监测信号，N为监测点数，j为傅里叶变换系数，ω为傅里叶变换频率；

(2)对频率f求导取整得到周期t_p，

式中，Round为取整函数；

(3)计算t_p中对应的自相关系数，通过所述自相关系数来确定包络周期S，其中自相关系数最大的且满足大于0.5的所对应的周期即为包络周期S，如果不存在满足条件的周期则包络周期S取两倍的采样间隔，自相关系数计算公式为：

式中，

为x(t)的平均数，k为滞后阶数。

优选地，所述步骤S3中，采用STL分解方法将所述包络曲线分解为趋势分量T_t、周期分量S_t和余项R_t，其表达式为：

X_t＝T_t+S_t+R_t t＝1，2,3...n

所述STL分解方法的分解过程为迭代过程，其包括外循环和嵌套在所述外循环中的内循环，每完成一次所述外循环需要进行多次所述内循环，在每次所述内循环过程中，更新所述周期分量和趋势分量，每次所述外循环过程中会计算鲁棒权重，所述鲁棒权重会用在内循环中，以减少异常值的影响；所述内循环和所述外循环均设定有循环次数，当所述内循环完成其设定的循环次数时，进行下一次所述外循环；当所述外循环完成其设定的循环次数时，会判断是否收敛，若收敛则输出所述趋势分量、周期分量和余项，并停止迭代，若不收敛，再次进行循环直到收敛。

优选地，在第v次迭代过程中，所述内循环包括以下步骤：

(1)去趋势：用原始序列减去在第v-1次迭代过程中获得的趋势分量，其表达式为：X_t-T_t(v)；

(2)周期子序列平滑：对每个子序列做LOESS(q＝n(s)，d＝1)回归，并向前后各延展一个周期，得到结果序列C_t(v+1)，t＝-n(p)+1,...,-n+n(p)；

(3)周期子序列的低通过滤：对结果序列C_t(v+1)依次做长度为n(p)、n(p)、3的滑动平均，然后做LOESS(q＝n(l)，d＝1)回归，得到结果序列L_t(v+1)，t＝1,2,...,n；

(4)平滑周期子序列的去趋势，其表达式为：S_t(v+1)＝C_t(v+1)-L_t(v+1)；

(5)去周期：原始序列减去周期分量，其表达式为：X_t-S_t(v+1)；

(6)趋势平滑，对于去周期之后的序列做LOESS(q＝n(t)，d＝1)回归，得到趋势分量T_t(v+1)。

优选地，在第一次迭代时，所述趋势分量的初始值T_t(0)＝0。

优选地，在第k次迭代过程中，所述外循环包括以下步骤：

(1)计算所述鲁棒权重，所述鲁棒权重的计算式为：

式中，R_t为上一次外循环输出的余项，median为中位值函数，B为Bisquare函数，Bisquare函数的计算公式为：

(2)进行多次所述内循环，并将所述鲁棒权重ρ_t用于所述内循环的LOESS回归中；当所述内循环完成其设定的循环次数后，输出所述所述趋势分量、周期分量和余项，然后进行下一次所述外循环；

(3)当所述外循环完成其设定的循环次数时，计算收敛值并判断是否收敛，其中收敛值为：

U_t＝T_torS_t

判断是否收敛需要满足的收敛条件为：

若满足所述收敛条件则输出所述趋势分量、周期分量和余项，并停止迭代，若不满足所述收敛条件，所述内循环和所述外循环设定的循环次数均进行加一并再次进行循环直到满足收敛条件。

优选地，在所述内循环的所述步骤(2)和步骤(6)中，在每次做LOESS回归时产生的邻域权重需要乘以所述鲁棒权重。

优选地，所述步骤S4中，所述判定曲线的计算式为：

Pan(t)＝Pan1+Pan2

式中，Pan1为对趋势曲线T(t)前h-1个值的平均估计；Pan2为趋势曲线T(t)前t个值的最大波动估计；

为T(t)前h-1个值的平均值；var为方差函数；diff为一阶差分；max、min分别为最大值、最小值函数；T(1:t)为趋势曲线T(t)中1到t的值；h为滞后阶数；abs为绝对值函数；所述判定曲线与所述趋势曲线T(t)的第一个交点t_g为所述声发射初至点。

本发明无需预先设置窗口宽度、阀值大小、特征函数等参数，因此摆脱了由窗口长度、特征函数、阀值等参数设置不合理所引起的误差；本发明通过对包络周期的确定、包络曲线的去周期变化和异常值后，去除了监测信号周期变化和异常值对初至拾取的影响，且平滑处理后的趋势曲线，与其他方法所选取的特征函数曲线相比更为平滑，并放大了初至的特征，使本发明在不同强弱的噪音环境下对初至拾取均有较强的适应性和准确性；本发明中的判定曲线是对趋势曲线波动的数学统计，其能对声发射前的噪音包络趋势变化进行自适应反馈，能够自动判别监测信号有无声发射，当岩体声发射出现时能自动捕捉这种趋势突变，进而自动识别到初至。

本发明的其他特征和优点将在随后的具体实施方式部分予以详细说明。

附图说明

图1是含有声发射的监测信号示意图；

图2是本发明的流程框图；

图3是本发明的具体流程图；

图4是本发明中STL分解的流程图；

图5是采用本发明识别方法对ω取0.2时对应五种不同噪音环境下的6组仿真信号拾取声发射初至的过程图；

图6至图11是分别采用人工法和采用本发明识别方法对26组实测信号的声发射初至的拾取结果对比图。

具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。其中，图1、图5、图6至图11中的横坐标Time为时间，纵坐标Amplitude为振幅。

如图2和图3所示，本发明提供一种岩体声发射初至自适应识别方法，包括以下步骤：

S1、获得监测信号的包络曲线

波形幅值变化是人工法、STA/LTA法等识别初至的重要依据，如图1所示，波形的幅值变化可以通过信号的波形包络来反映，声发射(简称AE)信号包络和噪音包络存在明显区别，根据众多实验结果和现场监测经验可知，岩石AE信号是呈哑铃状衰减波形，其波形包络先呈上升趋势后呈下降趋势，根据AIC准则可知在声发射到达前，噪音信号具有局部稳定性，其波形包络虽有一定的周期起伏，但变化趋势相对稳定，噪音与声发射包络趋势间的这种区别是本发明初至识别的重要依据。包络曲线反映了振幅的变化情况，因此，本发明首先对监测信号进行希尔伯特变换(Hilbert)来获得包络曲线。

具体地，对于一个监测信号x(t)，其希尔伯特变换为：

式中，*表示卷积运算，而监测信号x(t)可表示为

其中，f_s是载波频率，

是x(t)的相位调制信号，a(t)是x(t)的包络曲线。

对x(t)进行希尔伯特变换，可得解析信号：

其与x(t)对比后可求得包络曲线：

S2、确定所述监测信号的包络周期

噪音信号的幅值变化存在一定周期特征，这对于包络线变化趋势的判断存在一定的干扰，根据时间序列的加性模型可知，对于包络线趋势的准确判断需事先去周期变化的影响。包络线的周期特征本质上由监测信号本身决定，所以包络线的周期可以通过求解监测信号本身的周期来确定。

具体地，本发明采用周期图法和自相关系数相结合的方法来确定所述监测信号的包络周期，所述监测信号的包络周期的确定方法为：

(1)通过周期图法计算所述监测信号的功率谱，并选择前五个功率谱对应的频率f，其中，周期图法的基本原理是对观测到的数据直接进行傅立叶变换，然后取模的平方就是功率谱，功率谱的计算公式为：

(2)对频率f求导取整得到周期t_p，

式中，Round为取整函数；

式中，

为x(t)的平均数，k为滞后阶数。

S3、对所述包络曲线进行去除周期变化和异常值后获得趋势曲线

具体地，本发明采用STL分解方法将所述包络曲线a(t)分解为趋势分量T_t、周期分量S_t和余项R_t，其表达式为：

X_t＝T_t+S_t+R_t t＝1，2,3...n

所述STL分解方法的分解过程为迭代过程，其包括外循环和嵌套在所述外循环中的内循环，每完成一次外循环需要进行多次所述内循环，在每次所述内循环过程中，更新所述周期分量和趋势分量，每次所述外循环过程中会计算一次鲁棒权重，鲁棒权重为ρ_t(初始时ρ_t＝1)，鲁棒权重会用在内循环的LOESS回归中，用于减少周期分量和趋势分量中短暂的异常行为，以减少异常值对回归的影响；所述内循环和所述外循环均设定有循环次数，设定所述内循环的循环次数为n(i)，当所述内循环完成其设定的循环次数时，进行下一次所述外循环；设定所述外循环的循环次数为n(o)，当所述外循环完成其设定的循环次数时，会判断是否收敛，若收敛则输出所述趋势分量、周期分量和余项，并停止迭代，若不收敛，n(i)、n(o)均进行加一并再次进行循环直到收敛。

所述STL分解方法的具体分解过程如图4所示，在第v次迭代过程中，所述内循环包括以下步骤：

其中，所述内循环在第一次迭代时，所述趋势分量的初始值T_t(0)＝0。

在第k次迭代过程中，所述外循环包括以下步骤：

(1)计算鲁棒权重ρ_t，所述鲁棒权重的计算式为：

(2)进行多次所述内循环，并将所述鲁棒权重ρ_t用于所述内循环的LOESS回归中，以减小异常值的影响；当所述内循环完成其设定的循环次数后，输出所述所述趋势分量、周期分量和余项，然后进行下一次所述外循环；

U_t＝T_torS_t

判断是否收敛需要满足的收敛条件为：

其中，T_t和S_t分别为趋势分量和周期分量，判断收敛时可以将两者中的任一个代入上述左边的式子中，以判断是否满足收敛条件。

若满足所述收敛条件则输出所述趋势分量、周期分量和余项，并停止迭代，若不满足所述收敛条件，n(i)、n(o)均进行加一并再次进行循环直到满足收敛条件。

在所述内循环的所述步骤(2)和步骤(6)中，在每次做LOESS回归时会产生邻域权重，邻域权重需要乘以所述鲁棒权重ρ_t来得到相应步骤的最终权重，以减少异常值对回归的影响。

其中，n(i)为内循环次数，一般取2；n(o)为外循环次数，即鲁棒性迭代的次数，一般取10可以保证足够收敛；n(p)为一个周期的样本数，n(p)＝S；n(s)为季节项的平滑参数，n(s)＝7；n(l)为低通滤波平滑参数，取大于或等于n(p)的最小奇数；n(t)为趋势项的平滑参数，取大于(7～10)S的奇数。

S4、通过判定曲线对所述趋势曲线进行判断，以自动识别声发射初至点

通过上述步骤S3对包络曲线进行去除周期变化和异常值后获得趋势曲线T(t)，趋势曲线T(t)将呈现出明显的波峰，该波峰即是对应的声发射区域，峰脚即为声发射初至，为此，所述步骤S4中，本发明将定义以下判定函数对声发射初至进行自适应识别：

Pan(t)＝Pan1+Pan2

式中，Pan1为对趋势曲线T(t)前h-1个值的平均估计，Pan2为趋势曲线T(t)前t个值的最大波动估计，

为T(t)前h-1个值的平均值，var为方差函数；diff为一阶差分；max、min分别为最大值、最小值函数；T(1:t)为趋势曲线T(t)中1到t的值；h为滞后阶数，一般为10-15个采样间隔；abs为绝对值函数。

通过上式即可求得判定曲线，通过判定曲线对趋势曲线的波动情况进行衡量，其中所述判定曲线与所述趋势曲线T(t)的第一个交点t_g为所述声发射初至点。

综上所述，本发明是一种新的岩体声发射初至自适应拾取方法，其不同于现有的初至识别方法，由于其无需预先设置窗口宽度、阀值大小、特征函数等参数，因此摆脱了由窗口长度、特征函数、阀值等参数设置不合理所引起的误差。本发明通过对包络周期的确定、包络曲线的STL分解，去除了监测信号周期变化和异常值对初至拾取的影响，且平滑处理后的趋势曲线，与其他方法所选取的特征函数曲线相比更为平滑，并放大了初至的特征，使本发明在不同强弱的噪音环境下对初至拾取均有较强的适应性和准确性；本发明中的判定曲线是对趋势曲线波动的数学统计，其能对声发射前的噪音包络趋势变化进行自适应反馈，能够自动判别监测信号有无声发射，当岩体声发射出现时能自动捕捉这种趋势突变，进而自动识别到初至。

为了验证本发明初至拾取性能的准确性和稳定性，以下将通过仿真试验以及实测信号处理来具体说明。

一、仿真试验

通过MATLAB生成如下动态仿真信号：

y＝Ay₁+y₂

y₁＝sin(ωπt)+0.5sin(7.5ωπt)+rand

其中，Ay₁为模拟背景噪音，背景噪音的强度、主频率可以通过A、ω来控制，y₂为模拟AE信号，相应的t_g为AE初至，仿真中取700s，rand为随机函数，仿真中采样间隔为1s，采样长度为1024s。

在仿真中将通过控制信噪比强度，生成了五种不同噪音环境，其信噪比分别为0dB、3dB、5dB、10dB、20dB，每种噪音环境中又产生了3类不同主频率的噪音，共获得了15类噪音环境，在每一类噪音环境中又通过rand函数产生5组不同的随机模拟信号，共计获得75组信号，具体见表1。为反映本发明的准确性，分别与运用较多的STA/LTA法、峰度拾取法(高阶统计法)、AR-AIC法的拾取结果进行了比较，具体结果如表1所示。其中STA/LTA法与峰度拾取法的长短时窗分别为16s、4s，阀值分别为2、3.2，AR-AIC法的滑动窗为16s。

由表1可以看出，采用本发明的识别方法，在不同噪音环境中本发明均未出现漏拾，相较于其他方法本发明具有更高的稳定性和准确性，STA/LTA法与峰度拾取法对参数依赖较大，所以其漏拾率较大，且识别的初至偏差极大，基本上可以认为是误拾，尤其是在高频噪音环境下更是出现了100％的漏拾率。相比之下，AR-AIC法更为稳健，但其在纯噪音环境下(SNR＝0dB)却出现了100％的误判。相较而言，本发明在不同的噪音环境下均具有较高的准确性和稳健性，尤其是在低信噪比(SNR＝3、5dB)和高频(ω＝0.4)噪音环境下不会出现误拾和漏拾得现象，表现较好的拾取性能。

表1不同噪音环境下得到的初至拾取结果

图5是采用本发明识别方法对ω取0.2时对应不同SNR下的6组仿真信号拾取声发射初至的过程图，其中各信号对应的包络周期分别为11s、10s、10s、12s、12s，滞后阶数均取15。从图5中可以发现，噪音强度的增强(SNR缩小)，有效信号逐渐被背景噪音淹没，从幅值变化上也越来越难以准确识别初至，但在包络曲线上噪音与有效信号间依然可以得到辨别，经过STL分解后噪音与有效信号间的区别更加明显，有效信号出现前噪音包络趋势基本趋于平坦，而有效信号出现后包络趋势峰急剧上升，对应的判别曲线因对趋势变化存在滞后，进而与趋势曲线产生交点，继而拾取到了有效信号的初至。从图5上可以发现，纯噪音信号(SNR＝0)包络趋势虽存在波动变化，但因为不存在剧烈上升的趋势峰，与判定曲线未能产生相交，所以对于纯噪音信号本发明仍可以做出准确判断。

二、实测信号处理

表2为本发明与其他三种方法对于广西中金岭南矿业有限责任公司盘龙铅锌矿声发射地压监测系统日常监测所获取的26组实测信号的声发射初至拾取结果，其中每组监测信号采样数1024个，采样间隔为200μs，具体波形可见图6至图11。

从表2和图6至图11中可以看出本发明的拾取结果与人工法拾取的结果较为相近，在16、20号数据中本发明拾取结果与人工法一致，相较于其他方法本发明与人工法的拾取平均偏差明显较小，准确性明显更高。从表2可以发现，本发明和其他方法均未出现漏选现象，但其他三种方法最大偏差明显过大，分别为679×200μs、905×200μs、213×200μs，而本发明最大偏差为9×200μs，这进一步反映了本发明的准确性。一般认为偏差大于100×200μs即可认为产生了误选，相应的计算出本发明与其他三种方法的误选率分别为0％、11.5％、38.5％、57.7％，从误选率上可以看出本发明稳定性明显更强。根据实测信号拾取结果表明，相较于其他方法本发明准确性明显更高，稳定性更好。

表2实测信号初至拾取结果

以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合，为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。