一种岩爆信号阈值动态调整的多时窗简化形式识别算法

摘要

本发明公开了一种岩爆信号阈值动态调整的多时窗简化形式识别算法，包括以下步骤：获取各类典型岩爆信号和典型噪音信号，建立岩爆信号阈值动态调整数据库；确定必要参数和判断阈值；自动识别处理；定期进行判断阈值动态调整。本发明适用于微震数据实时处理，满足工程实际需要，提高波形拾取效率，减少人工识别工作量，进而提高岩爆、矿震、塌方等地质灾害预警的实时性。

说明书

技术领域

本发明涉及微震监测技术领域。具体涉及一种岩爆信号阈值动态调整的多时窗简化形式识别算法，该方法可广泛用于矿业工程、水利水电工程、石油工程、岩土工程以及地下工程。

背景技术

岩爆信号自动识别技术是微震源定位的关键。现有的岩爆信号识别方法主要有：根据在时间域能量和能量变化构建特征函数时间域的STA/LTA算法；根据岩爆信号与噪音信号波形特征的差别，如Fisher判别法、快速傅里叶变换、最大似然分类法、逻辑回归和神经网络、形态分形维数、统计方法、能量极值法等。然而这些方法大都仅对爆破信号具有较好的滤除效果，对实际工程中常见的锚杆钻机信号、电气信号、二次爆破信号等少有或没有效果，这导致了岩爆信号淹没在大量的数据中，很难通过软件自动分析处理。

通过实际工程应用，发现在自动识别岩爆信号时算法需要考虑以下三点：1)实时处理中计算量是否能够满足软硬件条件；2)阈值如何快速有效地确定满足工程中的复杂波形；3)是否能够适应工程中岩爆信号的变化。

因此，针对上述问题，兼顾算法实用条件，发明一种岩爆信号阈值动态调整的多时窗简化形式识别算法，适用于微震数据实时处理，满足工程实际需要，提高波形拾取效率，减少人工识别工作量，进而增强岩爆、矿震、塌方等地质灾害预警的实时性。

发明内容

本发明的目的在于克服现有技术存在的问题，兼顾算法实用条件，提供了一种岩爆信号阈值动态调整的多时窗简化形式识别算法，适用于微震数据实时处理，满足工程实际需要，提高波形拾取效率，减少人工识别工作量，进而增强岩爆、矿震、塌方等地质灾害预警的实时性。

一种岩爆信号阈值动态调整的多时窗简化形式识别算法，包括以下步骤：

步骤1、挑选典型岩爆信号和典型噪音信号。

岩爆信号一般指岩爆发生前岩体内部渐进破坏产生的弹性波信号，噪音信号一般包括有爆破信号、锚杆钻机信号、电气信号、现场作业信号等等。各类型岩爆信号/噪音信号的波形振幅和频率变化比较相似，波形能够代表某一类岩爆信号/噪音信号的信号为典型信号。

步骤2、求取典型岩爆信号的RK函数和典型噪音信号的RK函数。

步骤3、以时间为横坐标，RK函数值为纵坐标，

确定参考点的个数a，选取参考点r_j，参考点r_j的横坐标为延迟位置t₀+d_j，纵坐标为判断阈值R_j，判断阈值R_j位于延迟位置t₀+d_j处的典型岩爆信号的RK函数值和延迟位置t₀+d_j处的典型噪音信号的RK函数值之间，延迟位置t₀+d_j处的典型岩爆信号的RK函数值和延迟位置t₀+d_j处的典型噪音信号的RK函数值不相等；其中，a∈{1，2}，j为参考点的序号，t₀为触发时刻，d_j为延迟长度；

在本发明中，参考点r_j可以为1个或者2个，即参考点r_j可以为参考点r₁，也可以为参考点r₁和参考点r₂；判断阈值R_j对应为1个或者2个，即可以为判断阈值R₁，可以为判断阈值R₁和判断阈值R₂。

步骤4、对待识别微震数据的RK函数值进行识别，包括以下步骤：

步骤4.1、实时读取产生的待识别微震数据；

步骤4.2、实时计算待识别微震数据的RK函数值，待识别微震数据的RK函数值第一次大于触发阈值R_q对应的时刻为触发时刻t₀，触发阈值R_q为预设值；

步骤4.3、计算延迟位置t₀+d_j待识别微震数据的RK函数值，

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值大于典型噪音信号的RK函数值的情况下：

若延迟位置t₀+d_j处的待识别微震数据的RK函数值大于判断阈值R_j，则识别本次读取的待识别微震数据为待定岩爆信号；否则，本次读取的待识别微震数据为噪音信号；

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值小于典型噪音信号的RK函数值的情况下：

若延迟位置t₀+d_j处的待识别微震数据的RK函数值小于判断阈值R_j，则识别本次读取的待识别微震数据为待定岩爆信号；否则，本次读取的待识别微震数据为噪音信号；

若所有延迟位置t₀+d_j处的识别结果均为待定岩爆信号，则本次读取的待识别微震数据为岩爆信号。

在本发明中，参考点r_j可以为1个或者2个，对应的延迟位置t₀+d_j可以为延迟位置t₀+d₁，也可以为延迟位置t₀+d₁和延迟位置t₀+d₂。对于有两个参考点的情况下，必须两个参考点对应的延迟位置处的识别均为待定岩爆信号，才能确定本次读取的待识别微震数据为岩爆信号。

如上所述的判断阈值R_j的动态调整由以下步骤实现：

步骤5.1、建立岩爆信号阈值动态调整数据库，岩爆信号阈值动态调整数据库内按照发生时间的早晚顺序存储岩爆信号，岩爆信号阈值动态调整数据库内各个岩爆信号记为X₁，X₂，…，X_p，其中p为岩爆信号阈值动态调整数据库内岩爆信号的最大个数，发生时间最晚的岩爆信号是X_p；

一般情况下，岩爆信号阈值动态调整数据库内岩爆信号的最大个数p取400；

若现场地质条件多变，如深埋隧洞掘进或大型厂房开挖等工程可减少岩爆信号阈值动态调整数据库内岩爆信号的最大个数p，使阈值动态调整加快；

若现场地质条件变化较小，如固定点长期监测等工程可增加岩爆信号阈值动态调整数据库内岩爆信号的最大个数p，使阈值动态调整减慢。

步骤5.2、将岩爆信号阈值动态调整数据库内各个岩爆信号在延迟位置t₀+d_j的RK函数值进行排列构成的数组记为RK_j-(k)，其中k为岩爆信号的序号，j为参考点的序号；

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值大于典型噪音信号的RK函数值的情况下：岩爆信号阈值动态调整数据库中各个岩爆信号在延迟位置t₀+d_j的RK函数值从大到小进行排列；

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值小于典型噪音信号的RK函数值的情况下：岩爆信号阈值动态调整数据库中各个岩爆信号在延迟位置t₀+d_j的RK函数值从小到大进行排列；

步骤5.3、确定选取率B_j，选取满足下式的最小值为选取率B_j

其中，B_z为总选取率，B_z的范围值为50％～100％；

步骤5.4、获取k_j值，k_j＝int(p×B_j)，int为取整运算；

若k_j﹤p，则判断阈值R_j为RK_j-(k_j)和RK_j-(k_j+1)的平均值；

若k_j＝p，则判断阈值R_j为RK_j-(k_j)；

步骤5.5、获取新微震数据，新微震数据中的岩爆事件数a₁和噪音事件数b₁均为已知参数，

根据判断阈值R_j利用步骤1～4，对新微震数据进行自动识别处理获得岩爆信号和噪音信号，获得对应的岩爆事件数a₂与噪音事件数b₂；将新获得的岩爆信号按照发生时间进行排列，并记为Y₁，Y₂，…，Y_q，q为新获得的岩爆信号的最大个数，q≤p，发生时间最晚的岩爆信号是Y_q；

一般认为，设定时间内的发生的岩爆信号/噪音信号的个数大于或等于4可确定为一个岩爆事件/噪音信号，不能成为事件的所有信号都将滤除。其中，上述设定时间一般规定为0.5s。

步骤5.6、根据公式(5)计算新微震数据岩爆事件自动识别正确率E

若新微震数据岩爆事件自动识别正确率E大于预设的岩爆事件自动识别正确率阈值R_b，则进行步骤5.7，否则进行步骤5.8；其中，岩爆事件自动识别正确率阈值R_b的取值范围为50％～100％；

步骤5.7、根据公式(6)计算新微震数据噪音事件自动识别正确率G

若新微震数据噪音事件自动识别正确率G大于噪音事件识别正确率阈值R_g，则确认判断阈值R_j，否则进行步骤5.8；其中，噪音事件识别正确率阈值R_g的取值范围为0～50％；

步骤5.8、以新获得的岩爆信号Y_s替换岩爆信号阈值动态调整数据库中的岩爆信号X_s，s值加1，其中，s初值为1，s∈{1，2…，(q+1)}，判断s是否大于q，若s大于q，则表示阈值动态调整失败，将s设定为初值1，并在设定时间后返回步骤5.1；否则返回步骤5.2。

如上所述的典型岩爆信号的RK函数和典型噪音信号的RK函数基于公式(1)，

RK(t)＝R(t)公式(1)

公式(1)中的R函数基于公式(2)，

其中，STA(t)为短时窗STA函数，LTA(t)为长时窗LTA值，t为时刻，n为STA短时窗长度，m为LTA长时窗长度，

公式(2)中的CF函数基于公式(3)，

CF(t)＝Y(t)²-Y(t-1)·Y(t+1)公式(3)

其中，Y(t)为幅值函数。

本发明相对于现有技术，具有以下有益效果：针对岩爆信号特点，兼顾算法实用条件，提供了一种岩爆信号阈值动态调整的多时窗简化形式识别算法，适用于微震数据实时处理，满足工程实际需要，提高波形拾取效率，减少人工识别工作量，进而提高岩爆、矿震、塌方等地质灾害预警的实时性。

附图说明

图1为以RK函数值为纵坐标，时间(用采样点数表达)为横坐标，RK函数的曲线图；

图2为实施例中参考点r₁和参考点r₂的示意图；其中，(a)为参考点r₁的选取；(b)为参考点r₂的选取；

图3为识别算法验证图；其中，(a)为RK(2040)值密度函数分布规律；(b)为RK(2130)值密度函数分布规律。

具体实施方式

以下结合附图对本发明的技术方案进一步详细说明：

为了使本发明的技术手段、创作特征、工作流程、使用方法达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。本发明的保护范围不受以下实例的限制。

红透山铜矿位于辽宁省抚顺市清原满族自治县境内，是我国典型的深部金属矿山之一。该铜矿深部采场多处发生的岩爆、应力型塌方等动力型破坏，已经成为制约该矿山安全开采的首要问题。现场微震监测过程中发现，由于采场及巷道多采用电气设备，微震数据中存在大量的噪音信号，包括电气信号、二次爆破信号、机车撞击铁轨信号、锚杆钻机信号等等。本实例以红透山铜矿深部采场微震数据为例加以说明。

一种岩爆信号阈值动态调整的多时窗简化形式识别算法，包括以下步骤：

步骤1、挑选典型岩爆信号和典型噪音信号。

岩爆信号一般指岩爆发生前岩体内部渐进破坏产生的弹性波信号，噪音信号一般包括有爆破信号、锚杆钻机信号、电气信号、现场作业信号等等。各类型岩爆信号/噪音信号的波形振幅和频率变化比较相似，波形能够代表某一类岩爆信号/噪音信号的信号为典型信号。

本实例中，从红透山铜矿深部采场实时监测的数据选取时间段为2016年6月1日微震数据。从该数据中挑选出典型岩爆信号、典型低振幅电气信号、典型高振幅电气信号、典型短持续电气信号、典型长持续电气信号、典型二次爆破信号、典型机车撞击铁轨信号、典型锚杆钻机信号。

步骤2、求取典型岩爆信号的RK函数和典型噪音信号的RK函数。

典型岩爆信号的RK函数和典型噪音信号的RK函数基于公式(1)

RK(t)＝R(t)公式(1)

其中，R函数为微震监测系统自带触发算法中的STA/LTA函数，常用的计算方式基于公式(2)，但不限于此公式

其中，STA(t)为短时窗STA函数，LTA(t)为长时窗LTA值，t为时刻，n为STA短时窗长度，m为LTA长时窗长度，

公式(2)中的CF函数为微震监测系统自带触发算法中的特征函数，常用的计算方式基于公式(3)，但不限于此公式

CF(t)＝Y(t)²-Y(t-1)·Y(t+1)公式(3)

其中，Y(t)为幅值函数。

步骤3、以时间为横坐标，RK函数值为纵坐标，

确定参考点的个数a，选取参考点r_j，参考点r_j的横坐标为延迟位置t₀+d_j，纵坐标为判断阈值R_j，判断阈值R_j位于延迟位置t₀+d_j处的典型岩爆信号的RK函数值和延迟位置t₀+d_j处的典型噪音信号的RK函数值之间，延迟位置t₀+d_j处的典型岩爆信号的RK函数值和延迟位置t₀+d_j处的典型噪音信号的RK函数值不相等；其中，a∈{1，2}，j为参考点的序号，t₀为触发时刻，d_j为延迟长度。

在本发明中，参考点r_j可以为1个或者2个，即参考点r_j可以为参考点r₁，也可以是参考点r₁和参考点r₂；判断阈值R_j对应为1个或者2个，即可以为判断阈值R₁，可以为判断阈值R₁和判断阈值R₂。

本实例中，微震监测系统自带的触发算法中触发时刻t₀为2000，STA短时窗长度为20，LTA长时窗长度为200，CF函数基于公式(3)，R函数基于公式(2)。如图2所示，根据典型岩爆信号和典型噪音信号差异情况，确定参考点的个数a为2，参考点r₁的坐标为(2040，2.5)，参考点r₂的坐标为(2130，1)，即确定延迟位置t₀+d₁为2040，判断阈值R₁为2.5，延迟位置t₀+d₂为2130，判断阈值R₂为1。

步骤4、对待识别微震数据的RK函数值进行识别，包括以下步骤：

步骤4.1、实时读取产生的待识别微震数据。

步骤4.2、实时计算待识别微震数据的RK函数值，待识别微震数据的RK函数值第一次大于触发阈值R_q对应的时刻为触发时刻t₀，触发阈值R_q为设立值，可以采用微震监测系统自带触发算法的参数。

步骤4.3、计算延迟位置t₀+d_j待识别微震数据的RK函数值，

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值大于典型噪音信号的RK函数值的情况下：

若延迟位置t₀+d_j处的待识别微震数据的RK函数值大于判断阈值R_j，则识别本次读取的待识别微震数据为待定岩爆信号；否则，本次读取的待识别微震数据为噪音信号；

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值小于典型噪音信号的RK函数值的情况下：

若延迟位置t₀+d_j处的待识别微震数据的RK函数值小于判断阈值R_j，则识别本次读取的待识别微震数据为待定岩爆信号；否则，本次读取的待识别微震数据为噪音信号；

若所有延迟位置t₀+d_j处的识别结果均为待定岩爆信号，则本次读取的待识别微震数据为岩爆信号。

在本发明中，参考点r_j可以为1个或者2个，对应的延迟位置t₀+d_j可以为延迟位置t₀+d₁，也可以为延迟位置t₀+d₁和延迟位置t₀+d₂。对于有两个参考点的情况下，必须两个参考点对应的延迟位置处的识别均为待定岩爆信号，才能确定本次读取的待识别微震数据为岩爆信号。

本实例中，判断条件为RK(2040)>2.5和RK(2130)<1，触发阈值R_q为6。

如上所述的判断阈值R_j由以下步骤实现动态调整：

步骤5.1、建立岩爆信号阈值动态调整数据库，岩爆信号阈值动态调整数据库内按照发生时间的早晚顺序存储岩爆信号，岩爆信号阈值动态调整数据库内各个岩爆信号记为X₁，X₂，…，X_p，其中p为岩爆信号阈值动态调整数据库内岩爆信号的最大个数，发生时间最晚的岩爆信号是X_p。

一般情况下，岩爆信号阈值动态调整数据库内岩爆信号的最大个数p取400；

若现场地质条件多变，如深埋隧洞掘进或大型厂房开挖等工程可减少岩爆信号阈值动态调整数据库内岩爆信号的最大个数p，使阈值动态调整加快；

若现场地质条件变化较小，如固定点长期监测等工程可增加岩爆信号阈值动态调整数据库内岩爆信号的最大个数p，使阈值动态调整减慢。

本实例中，岩爆信号阈值动态调整数据库内岩爆信号的最大个数p取400。

步骤5.2、将岩爆信号阈值动态调整数据库内各个岩爆信号在延迟位置t₀+d_j的RK函数值进行排列构成的数组记为RK_j-(k)，其中k为岩爆信号的序号，j为参考点的序号；

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值大于典型噪音信号的RK函数值的情况下：岩爆信号阈值动态调整数据库中各个岩爆信号在延迟位置t₀+d_j的RK函数值从大到小进行排列；

在延迟位置t₀+d_j处的典型岩爆信号的RK函数值小于典型噪音信号的RK函数值的情况下：岩爆信号阈值动态调整数据库中各个岩爆信号在延迟位置t₀+d_j的RK函数值从小到大进行排列。

本实例中，岩爆信号阈值动态调整数据库内各个岩爆信号在延迟位置t₀+d₁的RK函数值RK(2040)从大到小排列(即数组RK₁-(k)从大到小排列)，岩爆信号阈值动态调整数据库内各个岩爆信号在延迟位置t₀+d₂的RK函数值RK(2130)从小到大排列(即数组RK₂-(k)从小到大排列)。

步骤5.3、确定选取率B_j，选取满足下式的最小值为选取率B_j

其中，B_z为总选取率，B_z的范围值为50％～100％。

本实例中，总选取率B_z取95％，可得选取率B₁和B₂均为97.5％。

步骤5.4、获取k_j值，k_j＝int(p×B_j)，int为取整运算；

若k_j﹤p，则判断阈值R_j为RK_j-(k_j)和RK_j-(k_j+1)的平均值；

若k_j＝p，则判断阈值R_j为RK_j-(k_j)。

本实例中，k₁值和k₂值均为390。

步骤5.5、获取新微震数据，新微震数据中的岩爆事件数a₁和噪音事件数b₁均为已知参数，

根据判断阈值R_j利用步骤1～4，对新微震数据进行自动识别处理获得岩爆信号和噪音信号，获得对应的岩爆事件数a₂与噪音事件数b₂；将新获得的岩爆信号按照发生时间进行排列，并记为Y₁，Y₂，…，Y_q，q为新获得的岩爆信号的最大个数，q≤p，发生时间最晚的岩爆信号是Y_q。

一般认为，设定时间内的发生的岩爆信号/噪音信号的个数大于或等于4可确定为一个岩爆事件/噪音信号，不能成为事件的所有信号都将滤除。其中，上述设定时间一般规定为0.5s。

步骤5.6、根据公式(5)计算新微震数据岩爆事件自动识别正确率E

若新微震数据岩爆事件自动识别正确率E大于预设的岩爆事件自动识别正确率阈值R_b，则进行步骤5.7，否则进行步骤5.8；其中，岩爆事件自动识别正确率阈值R_b的取值范围为50％～100％。

本实例中，岩爆事件自动识别正确率阈值R_b取80％。

步骤5.7、根据公式(6)计算新微震数据噪音事件自动识别正确率G

若新微震数据噪音事件自动识别正确率G大于噪音事件识别正确率阈值R_g，则确认判断阈值R_j，否则进行步骤5.8；其中，噪音事件识别正确率阈值R_g的取值范围为0～50％。

本实例中，噪音事件识别正确率阈值R_g取30％。

步骤5.8、以新获得的岩爆信号Y_s替换岩爆信号阈值动态调整数据库中的岩爆信号X_s，s值加1，其中，s初值为1，s∈{1，2…，(q+1)}，判断s是否大于q，若s大于q，则表示阈值动态调整失败，将s设定为初值1，并在设定时间后返回步骤5.1；否则返回步骤5.2。

本实例中，阈值动态调整后，判断阈值R₁为1.306，判断阈值R₂为1.367。

本实例中，选取时间段为2016年6月2日至2016年6月12日的5109个微震数据，其中包括1602个岩爆信号，其他均为噪音信号，包括440个低振幅电气信号、419个高振幅电气信号、560个短持续电气信号、440个长持续电气信号、822个机车撞击铁轨信号、487个二次爆破信号、339个锚杆钻机信号，分别做出各类型中所有信号的RK(2040)值和RK(2130)值，如图3所示。由图3可以看出，所有岩爆信号与所有低振幅电气信号、所有高振幅电气信号、所有短持续电气信号、所有长持续电气信号的RK(2040)值密度分布存在明显差异，分别集中在2～5、1～2、1～2、1～3、1～3；所有岩爆信号与所有长持续电气信号、所有机车撞击铁轨信号、所有二次爆破信号、所有锚杆钻机信号的RK(2130)值密度分布存在明显差异，分别集中在0～1、0～3、1～3、0～2、1～3。RK(2040)值与RK(2130)值与判断条件RK(2040)>1.306、RK(2130)<1.367进行判断，可从微震数据中自动识别出绝大部分岩爆信号，并且滤除大部分噪音信号。

由图3可以分析，由于初始的判断阈值R₁和初始的判断阈值R₂十分依赖信号的典型情况，并且不能对识别效果进行评估，因此判断阈值的动态调整在判断阈值的确定方面具有较大优势。动态的判断阈值与定值的判断阈值相比有更低的、可以评估和控制的漏检率，因此能够保证更高的识别效率。

通过RK(2040)值对应判断关系RK(2040)>1.306和RK(2040)值对应判断关系RK(2130)<1.367，对时间段为2016年6月13日～19日实测数据进行自动识别处理。实测数据中已确定包含212个岩爆事件和5191个噪音事件，识别算法自动识别结果为202个岩爆事件和4222个噪音事件。岩爆事件自动识别正确率E达到95.28％，噪音事件自动识别正确率G达到81.33％，减少了78.32％人工识别的工作量。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。