谢学斌，叶永飞，高 山，刘 涛，王小平

（中南大学 资源与安全工程学院，湖南 长沙 410083）

到时拾取是进行声发射参数分析、声源定位、矩张量反演的重要前提［1－3］，在实际监测中到时拾取多采用人工识别和阈值法［4］，但极易受个人经验的影响，已不适用于大量监测数据中的到时识别。为此，学者们将地震领域相关算法引进到了AE信号到时拾取中，根据各算法的特点可归纳为两类：一是根据监测信号的变化特征识别到时，比较经典的有STA／LTA法［5］、AIC法［6］、高阶统计法［7］、MER法［8］、机器学习法［9］；二是综合判别分析，该类算法先通过某种算法粗略拾取到时，再结合其他算法确定精确到时，如STA／LTA与AIC法联合确定到时［10］。

前述众多算法存在2个共同的缺点，一是各种算法过度依赖参数选择，当参数选取不合适时往往会造成漏拾、误拾；二是为保证精度，需事先进行一定降噪处理，这不可避免地造成有用信息损失，进而影响到时拾取的可靠性。总之，上述算法虽能在一定程度上实现到时自动识别，但在矿山庞大的监测数据量和复杂的噪音环境下存在明显局限性。

为了避免上述方法的局限性，本文利用时间序列处理的相关算法，通过对监测信号的幅值包络进行去周期去残差，获得监测信号的幅值趋势变化特征，再通过趋势突变检测确定到时，避免了繁杂的参数设定和降噪处理，具有较强的适应性。

1 基于波形包络趋势变化的到时拾取方法

如图1所示，AE信号幅值变化能明显区别于纯噪音部分，这是众多到时识别算法的依据，而信号幅值变化程度可通过其包络线来反映［11］，但包络线存在周期变化特征和异常残差，极易对到时判断产生干扰。从图1可以发现AE信号出现前包络线的趋势是稳定的，AE信号出现后包络线的趋势明显突变，因此可通过趋势突变点来确定AE信号到时。

图1 监测信号示意

1.1 波形包络

如图1所示，波形包络可由其外部拐点来构成［11］，拐点可分为上下拐点，对应的包络线也可分为上下包络线，上下包络线可通过监测信号均值来区分，即大于监测信号均值的拐点构成波形上包络线，小于监测信号均值的拐点构成波形下包络线，本文只选取上包络线进行研究。

1.2 包络线隐周期估计

去除上包络线X b的周期变化，需事先确定包络线的隐周期，隐周期是包络线最有可能存在的周期，包络线的隐周期S与监测信号的隐周期相同，监测信号的隐周期可通过功率谱密度图进行大致估计［12］，即功率谱密度中最大极值点对应的频率f0求导后即为所求包络线的隐周期。

1.3 基于STL分解确定包络线的趋势曲线

确定上包络线Xb的隐周期后，包络线的变化趋势可通过STL分解获得。STL分解是通过LOESS回归将包络线X b分解为趋势分量T t、周期分量S t和残差R t：

STL分为内循环与外循环，内循环嵌套在外循环中，假定内循环中第v－1次结束时各分量为Tt（v）、S t（v），内循环可分为以下6个步骤：

1）去趋势，减去上一轮结果的趋势分量Xt－T t（v）；

2）LOESS回归进行周期子序列平滑；

3）平滑周期子序列的低通滤波；

4）去周期子序列趋势，S t（v＋1）＝Ct（v＋1）－L t（v＋1）；

5）去除周期项X t－S t（v＋1）；

6）LOESS回归进行趋势平滑得到T t（v＋1）。

外层循环主要通过产生鲁棒权重来抑制包络线中的异常值，对于时间为t的数据点，其鲁棒权重为：

式中B函数为Bisquare函数：

1.4 基于MK突变点检测确定到时

MK检测是一种稳健的非参数突变检验方法，可对包络线的趋势曲线突变进行检测，具体内容如下：

对于长度为l的趋势曲线序列Tt，构造秩序列：

式中s k为i时刻数值大于j时刻数值个数的累计。在时间序列随机独立的假设下，定义统计量：

式中UF1＝0；E（s k），Var（s k）为累计数s k的均值和方差。

按趋势序列T t逆序x l，x l－1，…，x1计算逆序秩序列s k′，再重复上述过程，同时使UB k＝－UF k（k＝l，l－1，…，1），UB1＝0。MK突变点检测步骤如下：

1）计算顺序秩序列s k，并按式（5）计算UF k。

2）计算逆序秩序列s k′，并按式（5）计算UB k。

3）给定显著水平α，寻找临界值±Uα，如果在两临界值±Uα之间UF k与UB k存在交点，该点对应的时刻t g即为岩石AE信号到时；如果不存在交点即无声发射出现。

2 仿真实验

为验证本文到时拾取算法在不同噪音环境下的实际性能，通过MATLAB生成如下动态仿真信号来模拟矿山监测信号：

式中Ay1为模拟背景噪音，其强度、主频率可通过A、ω来控制；y2为模拟AE信号；相应的t g为其到时，仿真中取700 s；rand为随机函数，仿真中采样间隔为1 s，采样长度为1 024 s。

如表1所示，在仿真实验中通过调节幅值参数A、频率参数ω，生成了极强、强、一般、弱共4种噪音环境（信噪比分别为3 dB、5 dB、10 dB、20 dB）下的12类共60组模拟信号。为评判本文方法的拾取性能，与运用较多的STA／LTA法进行了对比，结果如表1所示，表中误判率为拾取误差大于100 s事件所占比例，漏判率为未识别到时事件所占比例。为体现STA／LTA对相关参数依赖性，设置了2 s、4 s、8 s共3种不同短时窗对有效信号到时进行拾取，此外长时窗设为32 s，自动提取阈值设为2.5。从表1可以发现，在各类噪音环境下本文方法均能拾取到时，而STA／LTA法对短时窗的选择具有明显的依赖性，且准确率较低。当短时窗取2 s时，STA／LTA法均出现了误判；当短时窗取4 s时，STA／LTA法拾取性能明显得到改善，但随着噪音增强其误判、漏判明显增加；当短时窗取8 s时，STA／LTA法拾取偏差明显过大，且ω＝0.4时STA／LTA法已无法拾取到时。从表1可以看出，在仿真中本文方法准确度明显优于STA／LTA法，且对于低信噪比（SNR＝3）和高频率（ω＝0.4）噪音均有相对较好的拾取性能。

图2为ω＝0.2时，不同噪音环境下的4组模拟信号的波形。将图2（a）、（b）与图2（c）、（d）对比后可以发现，随着背景噪音增强，部分有效信号很难从幅值上进行区分辨别，从波形幅值上判断到时容易出现较大误差，相应的在强噪音环境下人工法和阈值法拾取到时很难保证有足够的准确率。对图2中的包络线进行STL分解后，得到图3所示的变化趋势曲线。从图3可以直观地发现，到时前后包络线变化趋势明显不同，而趋势的分界点可以通过MK突变检测来判断。对图3中的包络线进行MK检测，得到如图4所示的变化趋势曲线。从图4可以发现，在不同强度的噪音环境下，在0.05显著水平内UF曲线与UB曲线均有交点，即均能拾取到时，从到时拾取结果（图4中交点横坐标）可以看出，随着噪音强度增强，本文方法拾取偏差也会增大，但相比于STA／LTA法（见图5）却有更好的性能，SNR为10 dB和20 dB时本文方法偏差分别为11 s和6 s，而STA／LTA法偏差分别为80 s和39 s，明显大于本文方法，当SNR为3 dB和6 dB时STA／LTA法出现了误判。从图5可以发现，STA／LTA法对阈值的设置也存在明显的依赖性，不同噪音环境需要设置不同的阈值，STA／LTA法才能发挥作用。

图3 STL分解后包络线趋势变化曲线

图4 MK检测后包络线趋势变化曲线

图5 短时窗4 s时STA／LTA法到时拾取结果

3 矿山地压监测系统实测信号到时拾取实例

为验证本文方法在实际应用中的有效性，抽取广西盘龙铅锌矿地压监测系统实测的10组信号作为实例进行验证，拾取结果如表2所示，表中“无”代表不存在声发射。其中每组监测信号采样数1 024个，采样间隔为200μs。

表2 实测信号拾取结果表

对实测信号的拾取结果与人工法、STA／LTA法（短窗口取8×200μs、长窗口取32×200μs、阈值取2.5）拾取结果进行了对比，可以发现本文方法可以准确识别实测信号中有无声发射信号，而STA／LTA法在对4、5、8号信号进行识别时出现了误判。在准确度上，STA／LTA法对2、3、9号信号的拾取结果优于本文方法，但STA／LTA法对7、10号信号的拾取结果却出现了较大偏差，在平均偏差上本文方法明显比STA／LTA法小。由以上综合分析可以发现，本文方法在实测信号上的拾取表现明显优于STA／LTA法，且表现出较强的稳健性和适用性。

4 结 论

1）依据波形的起伏特征提出了一种新的到时拾取方法，该方法首先通过功率谱密度图获得信号包络线的隐周期，然后通过STL分解法对信号包络线进行去周期处理得到包络线的趋势变化曲线，最后通过MK突变检测确定AE信号到达时间。

2）仿真实验结果表明，在不同强度和不同频率的噪音环境下，本文方法均表现出较强的稳健性，且准确率也明显高于STA／LTA法。本文方法不需要事先设置参数，摆脱了传统方法对参数设置的依赖性，具有更强的适应性。

3）在实测信号中应用后发现，本文方法能很好地识别监测信号有无AE信号，且与人工法拾取的结果相差较小。在实测信号应用中，本文方法表现性能明显优于STA／LTA法，在矿山实测信号处理上具有较强的稳健性。