CN108919353B - 一种微震波形初至到时的自动分级拾取与优选方法 - Google Patents

Abstract

本发明公开了一种微震波形初至到时的自动分级拾取与优选方法，其包括四级步骤，即去噪处理‑到时快速搜索‑到时精确拾取‑到时判断与优选；利用MATLAB自带的巴特沃思Butterworth带通10‑200Hz滤波器，对原始信号进行去噪处理；利用峰值判断法快速判断信号内的峰值点，然后利用峰值点左右各推移500采样点，再利用STA/LTA法初步确立初至到时位置；初至到时精确拾取选用了两种方法，分别是AIC算法和MER算法，确立微震波形的初至到时位置；有效通道的判断与优选过程简述为三个步骤：聚类分析求取参照通道、时差判断剔除无效通道、包络曲线判断提出误差通道。本发明能够快速、准确拾取微震信号的初至到时，并能自动判断和剔除干扰或误差较大通道。

Description

一种微震波形初至到时的自动分级拾取与优选方法

技术领域

本发明涉及结构应用领域，特别涉及一种微震波形初至到时的自动分级拾取与优选方法。

背景技术

微震技术是监测、预警深部矿产开采、岩土工程诱发的冲击地压、岩爆灾害的有效手段，如何进一步提升微震技术快速、准确分析能力，如微震自动定位、无人值守快速预警，首先需要解决微震信号初至到时的准确拾取这一重要前提。快速、准确自动拾取初至到时是实现的微震数据快速处理、提升灾害应急处理的关键。

矿山微震监测或岩土领域微震监测不同于地震监测，由于监测范围小，P、S波难以区分，单分量传感器居多，因此如何进行P波拾取至关重要。而微震数据更易受到工况环境的影响，信号成分复杂、频率趋向高频，且触发的各类数据数量大，加之现场对微震数据处理时效性和定位精度要求，这给现场工程人员的处理和分析带来了很大难度。一方面，不同于常规的机械振动信号，工程现场的干扰信号种类繁多、特点各异，不仅有现场施工设备影响、人工干扰，还有微震***自身的传输干扰，同时，传播介质的不均质也会带来干扰，这些背景干扰的存在，使得微震信号内含成分复杂，并非单一的、周期性的。理论方法在理想信号或高信噪比信号中能去的较好效果，但在低信噪比信号中具有一定难度，识别精度有待提高，为此国内外专家开展了大量研究，如图1所示，为常用初至波到时拾取方法对比分析。

发明内容

为解决上述背景技术中存在的问题，本发明的目的在于提供一种微震波形初至到时的自动分级拾取与优选方法，以达到能够快速、准确拾取微震信号的初至到时，并能自动判断和剔除干扰或误差较大通道的目的。

为达到上述目的，本发明的技术方案如下：

一种微震波形初至到时的自动分级拾取与优选方法，其包括四级步骤，即去噪处理-到时快速搜索-到时精确拾取-到时判断与优选；

第一级：去噪处理，利用MATLAB自带的巴特沃思Butterworth带通10-200Hz滤波器，对原始信号进行去噪处理；

第二级：到时快速搜索，利用峰值判断法快速判断信号内的峰值点，然后利用峰值点左右各推移500采样点，再利用STA/LTA法初步确立初至到时位置；

第三级：到时精确拾取，初至到时精确拾取选用了两种方法，分别是AIC算法和MER算法，确立微震波形的初至到时位置；

第四级：到时判断与优选，有效通道的判断与优选过程简述为三个步骤：聚类分析求取参照通道、时差判断剔除无效通道、包络曲线判断剔除误差通道。

优选的，所述到时快速搜索具体为：首先利用500ms的时窗在整个5000ms信号内搜索，时窗不重叠，首次时窗的末尾采样点作为第二次搜索窗口的起点，快速在整个信号内搜索；在搜索的同时，计算时窗内振幅值绝对值的和，并取平均值，得到一系列搜索特征值；处理完整个信号后，求取特征值的最大值位置，以此为中心，分别向该次时窗两侧扩展，得到一个1500ms的大时窗。

优选的，所述聚类分析求取参照通道，为对到时差进行空间聚类分析，即k＝2，选择数据个体多的一类作为主簇，并将该簇内与聚心距离最小的点确立为最佳参照点；所述时差判断剔除无效通道，选取参照通道，分别计算其他各通道与参照通道的空间距离和实际到时时差，由空间距离和标定波速可以计算参照点与任一通道的理论最大到时时差，对比理论到时时差与实际到时时差的大小，选定判断阈值，即可剔除到时值误差大的通道，得到初步筛选的结果；所述包络曲线判断剔除误差通道，从前一部分初步筛选的结果中选择最先触发，即到时最小的通道C₀为参照通道，以这个点的到时为基准，求取其他各通道C_i(i为除0以外的各通道编号)与这个通道的实际到时差，以到时差和标定炮获得的波速求取半径值R_i，再以半径R_i画圆，对于内场的微震事件而言，其空间位置必然处于圆内或圆形之上，外场事件可通过初步定位，求取微震震源与各台站的距离U_i，当sum(U_i)最小时，即为最佳位置，然后反推R_i。

优选的，利用初步确立震源位置的方式对参与计算的通道进行筛选，所述计算采用单纯形法进行，初步计算出的定位结果为O’(x_f,y_f,z_f,t_f)，则该震源位置与各台站之间的传播时间可表述为：

时距法计算最佳震源实际上是求取到时差之和λ的最小值，即

其中N为台站数量，R_i为震源位置与各台站之间的距离，即半径。

通过上述技术方案，本发明提供的一种微震波形初至到时的自动分级拾取与优选方法，能够快速、准确拾取微震信号的初至到时，并能自动判断和剔除干扰或误差较大通道。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。

图1为现有技术中常用初至波到时拾取方法对比分析图；

图2为本发明AIC计算初至波到时图；

图3为本发明初至到时分级拾取与优选流程图；

图4为本发明台站及震源位置平面示意图；

图5为本发明带通滤波前后微震波形及其时频图对比图；

图6为本发明典型微震事件各通道波形及到时拾取结果图；

图7为本发明到时聚类及参照通道的确立图；

图8为本发明理论最大到时时差、实际时差对比图；

图9为本发明内场事件分布的包络曲线图；

图10为本发明自动拾取结果示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

本发明提供的一种微震波形初至到时的自动分级拾取与优选方法，如图3所示，其包括四级步骤，即去噪处理-到时快速搜索-到时精确拾取-到时判断与优选；

第一级：去噪处理，考虑到矿山微震信号的频率范围相对集中这一特点，利用MATLAB自带的巴特沃思Butterworth带通10-200Hz滤波器，对原始信号进行去噪处理；

第二级：到时快速搜索，利用峰值判断法快速判断信号内的峰值点，然后利用峰值点左右各推移500采样点，再利用STA/LTA法初步确立初至到时位置，具体为：首先利用500ms的时窗在整个5000ms信号内搜索，时窗不重叠，首次时窗的末尾采样点作为第二次搜索窗口的起点，快速在整个信号内搜索；在搜索的同时，计算时窗内振幅值绝对值的和，并取平均值，得到一系列搜索特征值；处理完整个信号后，求取特征值的最大值位置，以此为中心，分别向该次时窗两侧扩展，得到一个1500ms的大时窗，如图2中(a)、(b)所示，灰色背景方框为初步确立的初至时刻所处时窗，图2(b)中直线区域即为搜索特征值的曲线；

第三级：到时精确拾取，初至到时精确拾取选用了两种方法，分别是AIC算法和MER算法，进一步确立微震波形的初至到时位置，如图2(c)中所示，中部位置即为上述方法拾取得到的初至到时和终止时刻位置；

第四级：到时判断与优选，有效通道的判断与优选过程简述为三个步骤：聚类分析求取参照通道、时差判断剔除无效通道、包络曲线判断剔除误差通道，所述聚类分析求取参照通道，为对到时差进行空间聚类分析，即k＝2，选择数据个体多的一类作为主簇，并将该簇内与聚心距离最小的点确立为最佳参照点；所述时差判断剔除无效通道，选取参照通道，分别计算其他各通道与参照通道的空间距离和实际到时时差，由空间距离和标定波速可以计算参照点与任一通道的理论最大到时时差，对比理论到时时差与实际到时时差的大小，选定判断阈值，即可剔除到时值误差大的通道，得到初步筛选的结果；所述包络曲线判断剔除误差通道，从前一部分初步筛选的结果中选择最先触发，即到时最小的通道C₀为参照通道，以这个点的到时为基准，求取其他各通道C_i(i为除0以外的各通道编号)与这个通道的实际到时差，以到时差和标定炮获得的波速求取半径值R_i，再以半径R_i画圆，对于内场的微震事件而言，其空间位置必然处于圆内或圆形之上，外场事件可通过初步定位，求取微震震源与各台站的距离U_i，当sum(U_i)最小时，即为最佳位置，然后反推R_i，利用初步确立震源位置的方式对参与计算的通道进行筛选，所述计算采用单纯形法进行，初步计算出的定位结果为O’(x_f,y_f,z_f,t_f)，则该震源位置与各台站之间的传播时间可表述为：

时距法计算最佳震源实际上是求取到时差之和λ的最小值，即

其中N为台站数量，Ri为震源位置与各台站之间的距离，即半径。

其判断方法可简述为：如果震源位置与台站之间的到时误差越小，表明该台站参与定位计算的可靠性越高；反之，如果计算出的震源位置与台站之间的到时误差越大，表明该台站的可靠性越低。

在传播介质均质条件下，以台站位置为圆心，台站到震源处距离为半径画圆，则各台站的圆曲线在震源处聚合，如图4所示，△为监测台站，*为震源，由于传播介质不均匀、波速变化，该设想只存在于理论之中。但考虑上述因素影响，给定一定的富余ΔR，虚构半径R+ΔR，以此来判断到时值的有效性。

实例验算：

为了验证上述方法，以河北某矿微震***监测到的一次微震事件2011-10-07T13:42:17为例，对上述算法进行验证。实验场所为该矿某回采工作面，工作面内风巷、运巷共布置有12个微震传感器，参数信息参见表1，用于对该工作面回采工作进行实时动态监测。

表1现场台站布置参数

带通滤波Bandfilter

现场数据统计分析表明，矿山微震信号的频率范围受多方面因素影响，如传输距离、围岩岩性、震源自身等，这些成分的存在给到时拾取带来难度，同时也影响了拾取的精度。因此，在进行到时拾取之前，有必要首先对微震数据进行去噪处理，带通滤波在地震信号去噪中应用较为广泛，这其中尤以巴特沃思滤波器(Butterworth)最为著名。利用MATLAB编制了Bandfilter带通滤波器，考虑到上述因素，设置有效信号的频率分布范围为10～200Hz，利用滤波器对微震信号进行处理。图5所示为典型微震信号带通滤波的前、后效果，可以看出，信号的主要成分和能量分布于100～200Hz范围。在滤波前原始信号在高频区域还有稀疏分布，在去噪后高频区域成分已被消除。通过计算得到去噪后信号的信噪比SNR＝35.6126db，表明去噪后信号信噪比得到大幅提升。

初至到时拾取

图6为该事件各通道的波形图及初至拾取结果。图中标示了该次事件12个通道内的“时间-振幅”曲线，竖直线为本文方法所拾取得到的初至到时。Station1～12分别表示12个监测传感器。从图可以看出，在一定的信噪比条件下，本文方法可以快速拾取初至到时，如图中Station1、2、3、5、6、7、8、9以及11号通道，而Station4、10、12未能准确拾取到时。这其中4号通道内信噪比较低，本文算法拾取到时为0.003s，表示该次拾取无法准确获取到时，这说明信噪比高低严重影响到时的拾取；10和12号通道内无有效的波形起伏，表明上述两通道在事件发生时刻并未触发。上述通道拾取的到时及振幅如下表2所示。通道10和12中未能有效触发，由模型计算得到的初至到时为1s。

表2各通道内波形的初至到时

从图6和表2可以看出，信噪比的高低对初步拾取到时影响较大，信噪比越高，拾取到时约准确，反之，信噪比或异常信号的初至很难准确拾取或判断。如图6中4、6、10、12，这是通过肉眼能够直接观察的，但到时拾取算法如果返回值，而不进行判断和优选，这将对严重影响后续的定位计算精度。初至到时拾取完毕后，各到时值将被带入定位模型之中进行震源计算。上图中，除了Station1、2、3、5、6、7、8、9以及11号通道，通过肉眼初步判断为有效，另外三个通道为无效。

有效通道优选

有效通道的判断和选择，实际上是对各个初至时刻进行甄选，剔除无效的或误差较大的初至值，其目的是为了保证定位计算的精度。

(1)参照通道的选择

由于传感器于同一工作面内，按间隔30～50m间距布置，因此，对于台网布置内的微震事件而言，触发传感器的时间间隔不大。如果对这些触发时刻进行聚类计算，选用k-means聚类算法(按k＝2计算)，可以首先剔除误差比较大的通道。参照通道的选择是借助聚类分析方法实现的，如表3、图7所示，将12个初至时刻聚为两类，图7中“*”号为聚类质心(聚心)。

表3各站点初至到时聚类结果

从图7(a)中可以看出，除了10、12通道，其余通道为一类。由于A簇中元素远远多于B簇的，因此，选择初至时刻较多的A类继续研究，并以这一簇中离该类聚心(该簇内所有到时的平均值)最近的点(台站)作为参照物，剔除B类。当两簇元素数量持平时，可选择将两簇的中心作为聚心，求取与各台站的距离。图7(b)中显示的是各元素与聚心(A簇)的距离，可以看出，8#对应的距离最小，距离为7.57×10-4，因此，选择8号作为参照台站。

(2)无效通道的剔除

在进行第一步处理后，还需要解决初至到时值有误或误差较大的问题。图8为Station 8与其他台站的实际到时时差ΔTA(Actual Value)、理论最大到时时差ΔTT(Theoretical Max Value)对比图。

从图中可以看出，4、10、12三个通道的理论最大到时时差ΔTT远远小于实际到时时差ΔTA。为了描述理论最大到时时差与实际到时时差的关系，设定阈值λ＝ΔTA/ΔTT，当λ>10时则判断该通道为到时值误差较大通道。到时时差的相关参数计算结果如表4所示。

表4通道8与其他通道的到时时差

从上表计算结果可以看出，通道4、10和12中的λ值异常，远远大于理论最大值，这是不符合实际的。因此，这三个通道的到时值为废道，其到时拾取结果不宜参与后期的定位计算，应首先剔除这三个通道的到时值。除了剔除到时异常通道，还应考虑到时值误差较大通道。

(3)误差大通道的剔除

在剔除第4、10、12通道后，剩余1～3、5～9以及11共9个通道，再利用包络曲线判断法对剩余通道进行优化选取。在绘制包络曲线前，需要计算包络曲线的半径。半径值Ri可由真实到时差ΔTT乘以微震波传播波速求得，计算结果如表5中所示。包络曲线图如图9所示，从图中可以看出，各包络曲线密切分布于震时震源(“*”号标记)，多为两两相交状。

表5包络曲线最大半径值

对比图4中所示情况，结合理论分析可知，包络线相交的公共区域即为震源分布的区域。图9中同样印证了这一结论，除了通道6，其他包络线均将震源包裹在内。由此可见，通道6的到时值不适合参与定位计算。从波形角度来看，通道6中波形疑似电脉冲波形，这也直接证实了该通道到时值的不可取。

上述分析结果表明，4、10、12号通道的初至到时无法准确检测、6号通道的到时误差较大，均不宜参与该事件的定位计算，作为无效通道处理。1～3、5、7～9以及11号通道为有效通道，因此，最终选取上述通道参与该事件的定位计算。

对比分析

利用手动拾取方法对上述事件的初至到时进行拾取，并对比分析手动拾取与自动拾取的结果。经过对比，人工与自动方法的拾取精度和耗时对比如下表表6所示。通过实测，对上述12通道波形的初至到时进行拾取耗时3min，而利用计算求取仅需0.039s。

表6人工拾取与本文方法拾取结果对比

现场数据分析与应用

为了对上述分级检测与拾取算法的实际应用效果进行评判，选取矿山现场微震事件(265组波形数据)进行到时拾取计算，使用的计算环境为：操作***Windows 10专业版，处理器Interl(R)Core(TM)i7-7700K，CPU 4.20GHz，内存16GB。

到时拾取算法检验

为了定量对比几种方法的准确性，对上述265组数据进行了人工手动拾取(总耗时0.828s)，并以手动拾取到时作为参照，将该方法与其进行对比分析，结果如图10所示。其中，图10(a)为初至自动拾取结果与手动拾取结果的分布关系，图10(b)为到时误差分布直方图(仅显示50ms误差范围内结果)。

从两图中可以看出：利用上述所提出的拾取方法，拾取结果与手工拾取结果相近，大量自动拾取结果(*号)落入手动拾取圆圈之内；通过最终的误差分析也可以看出这一结论，如误差在-5～5ms区间的分布率达到71.7％(190组)，在-10～10ms之间的分布率达到83.4％(221组)，-20～20ms之间的分布率达到90.57％(240组)。这一结果再次证明上述所提出的到时拾取方法具有较高的拾取精度，基本可以代替人工拾取方法。

因此，总体上讲，本申请所提出的分级拾取方法一定程度上提高了自动拾取的精确度，可大大改善人工手动拾取的这一现象。

定位结果前后对比

为了验证通道优化对最终定位计算结果的影响，以前述标定炮为参照，对比真实值与优化前、后计算值的差异，以此评判文中所提出的分级检测和拾取方法的优劣性。所选择的定位方法为四四组合定位法，最终计算的结果如表7中所示。优化前，选择所有通道到时进行计算；优化后，选择优化选取的通道到时参与定位计算。

从表中可以看出，不同通道的选取对最终计算结果有较大影响，优化前与真实震源误差约为21.12m，而优化后误差减小至9.41m。优化前，以手动方式拾取到时和有效通道有两大缺陷，一是耗时长，二是无法对到时值差异较大的通道进行筛选，这将严重影响微震数据的处理速率；而通过优化，可以快速拾取到时，并对参与计算的通道进行优选，保证了参与计算通道的合理性，减小了计算误差。

表7优化前后定位结果对比

利用上述所提出方法，一定程度上解决了矿山微震波初至到时拾取的难题，并为微震定位计算提供了前提。

通过理论分析和工程实例验证，得出如下结论：

(1)带通滤波可以有效去除矿山微震信号中的干扰成分，提高微震信号的信噪比，为初至到时的拾取奠定基础。

(2)利用最大峰值和STA/LTA快速判断可以快速搜寻微震信号的初至时刻，这为后续的精确拾取提供了帮助，并且大大减少了数据处理的计算量和时间。

(3)改进AIC法拾取初至到时，解决了特征值出现多个波峰或波谷难以确立初至的难题，同时，结合STA/LTA算法，提高了到时拾取计算速率。两种方法相辅相成，有效提高了到时拾取的精度和速度。

(4)利用265组现场波形数据，对人工拾取和自动拾取两类方法的对比分析，结果表明，利用上述方法拾取结果平均误差在-5～5ms区间的分布率达到71.7％(190组)，在-10～10ms之间的分布率达到83.4％(221组)，-20～20ms之间的分布率达到90.57％(240组)，此外，本申请所提出方法算法耗时为0.828s，远远快于人工手动拾取。

本发明公开的一种微震波形初至到时的自动分级拾取与优选方法，能够快速、准确拾取微震信号的初至到时，并能自动判断和剔除干扰或误差较大通道。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (1)

1.一种微震波形初至到时的自动分级拾取与优选方法，其特征在于，其包括四级步骤，即去噪处理-到时快速搜索-到时精确拾取-到时判断与优选；

第一级：去噪处理，利用MATLAB自带的巴特沃思Butterworth带通10-200Hz滤波器，对原始信号进行去噪处理；

第二级：到时快速搜索，利用峰值判断法快速判断信号内的峰值点，然后利用峰值点左右各推移500采样点，再利用STA/LTA法初步确立初至到时位置，具体为：首先利用500ms的时窗在整个5000ms信号内搜索，时窗不重叠，首次时窗的末尾采样点作为第二次搜索窗口的起点，快速在整个信号内搜索；在搜索的同时，计算时窗内振幅值绝对值的和，并取平均值，得到一系列搜索特征值；处理完整个信号后，求取特征值的最大值位置，以此为中心，分别向该次时窗两侧扩展，得到一个1500ms的大时窗；

第三级：到时精确拾取，初至到时精确拾取选用了两种方法，分别是AIC算法和MER算法，确立微震波形的初至到时位置；

第四级：到时判断与优选，有效通道的判断与优选过程简述为三个步骤：聚类分析求取参照通道、时差判断剔除无效通道、包络曲线判断剔除误差通道，所述聚类分析求取参照通道，为对到时差进行空间聚类分析，即k＝2，选择数据个体多的一类作为主簇，并将该簇内与聚心距离最小的点确立为最佳参照点；所述时差判断剔除无效通道，选取参照通道，分别计算其他各通道与参照通道的空间距离和实际到时时差，由空间距离和标定波速可以计算参照点与任一通道的理论最大到时时差，对比理论到时时差与实际到时时差的大小，选定判断阈值，即可剔除到时值误差大的通道，得到初步筛选的结果；所述包络曲线判断剔除误差通道，从前一部分初步筛选的结果中选择最先触发，即到时最小的通道C₀为参照通道，以这个点的到时为基准，求取其他各通道C_i与这个通道的实际到时差，i为除0以外的各通道编号，以到时差和标定炮获得的波速求取半径值R_i，i为除0以外的各通道编号，再以半径R_i画圆，对于内场的微震事件而言，其空间位置必然处于圆内或圆形之上，外场事件可通过初步定位，求取微震震源与各台站的距离U_i，当sum(U_i)最小时，即为最佳位置，然后反推R_i；利用初步确立震源位置的方式对参与计算的通道进行筛选，所述计算采用单纯形法进行，初步计算出的定位结果为O’(x_f,y_f,z_f,t_f)，则该震源位置与各台站之间的传播时间可表述为：

时距法计算最佳震源实际上是求取到时差之和λ的最小值，即

其中N为台站数量，R_i为震源位置与各台站之间的距离，即半径。

Images