CN114294062B - 一种冲击地压时空动态综合预警方法 - Google Patents

Abstract

本发明公开了一种冲击地压时空动态综合预警方法，包括：基于微震现场实际监测数据计算冲击地压前兆预警指标；以各预警指标的时序变化趋势作为预警准则，计算各预警指标的预警效能，筛选出预警效能高的指标作为优选指标；将优选指标进行数据融合，得到时序冲击地压危险量化预警指标，判定矿井大区域范围内冲击危险等级；对矿井空间进行网格划分，计算各网格内优选指标值；基于各网格内优选指标值和各优选指标的预警效能计算各网格内空间冲击地压危险指标，绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区，根据局部冲击危险区和区域危险等级采取针对性防冲措施。本发明可实现冲击地压的精准预警，并有效降低其防治成本。

Description

一种冲击地压时空动态综合预警方法

技术领域

本发明涉及地下开挖工程和煤岩动力灾害预警技术领域，特别涉及一种冲击地压时空动态综合预警方法。

背景技术

冲击地压是井下生产过程中面临的一种典型动力灾害，其具有突发性、猛烈性等特征，一旦发生则易出现群死群伤事故，给井下相关工作人员的生命和财产安全带来了极大威胁。并且随着浅部矿产资源的日近枯竭，煤炭开采逐年向地球深部扩展，受地应力、地质条件等多种因素的综合影响，冲击地压事故发生频次逐年上升，其发生机理也愈发复杂，给煤矿的安全生产造成严重挑战

为保障煤矿的安全生产，近年来，针对冲击地压监测预警的技术取得了诸多优秀成果，研究表明微震、地音、电磁辐射、地应力等在线监测系统具有良好的现场应用效果，并且其中微震监测系统可以对矿井大范围区域内煤岩体破断进行实时监测，被认为是冲击地压防治效能最高的监测手段之一。

但是，井下相关人员应用微震系统进行监测预警的过程中，常需要将独立的微震事件在空间中简化为“点源”，并且受系统定位误差、台网布置等因素的影响，其预警准确率还有较大上升空间，为保障井下工作人员的生命和财产安全，常需大范围施行各项防冲措施，其效果难以检验，造成了部分资源的浪费和开采成本的上升。因此，有必要提出一种能适应井下复杂开采条件，为相关人员提供高效、准确防冲决策依据的预警方法。

发明内容

本发明提供了一种冲击地压时空动态综合预警方法，以解决现有技术所存在的预警准确率不够高，造成了部分资源浪费和开采成本上升的技术问题。

为解决上述技术问题，本发明提供了如下技术方案：

一种冲击地压时空动态综合预警方法，包括：

基于微震现场实际监测数据，计算冲击地压前兆预警指标；

判定各预警指标的时序变化趋势，以各预警指标的时序变化趋势是否符合预警前兆表征规律作为其是否触发预警的标准，并计算各预警指标的预警效能，筛选出预警效能大于预设预警效能阈值的预警指标，作为优选指标；

将所述优选指标进行数据融合，得到时序冲击地压危险量化预警指标，基于所述时序冲击地压危险量化预警指标，判定矿井大区域范围内冲击危险等级；

对矿井空间进行网格划分，计算各网格内优选指标值；基于各网格内优选指标值和各优选指标的预警效能，计算各网格内空间冲击地压危险指标，绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区；

根据局部冲击危险区和区域危险等级采取针对性防冲措施。

进一步地，所述微震现场实际监测数据指的是利用安装在矿井的微震监测系统采集到的一系列独立微震事件的相关信息，包括：各微震事件的发生时间、X坐标、Y坐标、Z坐标、能量以及事件类型；

所述预警指标包括正向预警指标和负向预警指标，所述正向预警指标指的是冲击地压事件发生前出现增长趋势的指标，所述负向预警指标指的是冲击地压事件发生前出现降低趋势的指标；其中，所述正向预警指标包括：日最大能量、日总能量、日总频次、日平均能量、能量偏差值、频次偏差值、平均总能量、微震活动度、缺震、A(b)值、断层总面积、平均总频次、微震活动标度和算法复杂度；所述负向预警指标包括：震源集中程度、b值、P(b)值和时间信息熵。

进一步地，所述判定各预警指标的时序变化趋势，以各预警指标的时序变化趋势是否符合预警前兆表征规律作为其是否触发预警的标准，并计算各预警指标的预警效能，筛选出预警效能大于预设预警效能阈值的预警指标，作为优选指标，包括：

利用Mann-Kendall趋势检验法判定各预警指标的时序变化趋势；

以各预警指标的时序变化趋势是否符合预警前兆表征规律作为其是否触发预警的标准，并计算各预警指标的预警效能；其中，所述预警前兆表征规律是指预警指标随冲击地压危险程度的变化而出现的响应特征，所述响应特征为：所述正向预警指标在冲击危险程度增加时出现的异常高值、连续增加或持续在高值附近波动，所述负向预警指标在冲击危险程度增加时出现的异常低值、连续降低或持续在低值附近波动；所述触发预警的标准为：若正向预警指标存在增长趋势或负向预警指标存在降低趋势，则触发预警，否则，不触发预警；

筛选出预警效能大于预设预警效能阈值的预警指标，作为优选指标。

进一步地，所述利用Mann-Kendall趋势检验法判定各预警指标的时序变化趋势，包括：

将各预警指标前一段时间内的计算结果组成时间序列，记为X[x₁,x₂,x₃,...,x_n]，其中，n表示时间窗长度；首先计算X的检验统计量S：

其中：

式中，x_k表示第k个指标，k＝1,2,3,...,i-1，x_i表示第i个指标，i＝k+1,2,3,...,n；

再计算X的检验标准量Z：

当Z>0时，相应的预警指标具有增长趋势；当Z<0时，相应的预警指标具有降低趋势；当Z＝0时，相应的预警指标不具有明显变化趋势。

进一步地，所述预警效能的计算公式为：

其中，F表示预警指标的预警效能，TP表示预警指标触发预警，并且实际发生了冲击地压事件，FP表示预警指标未触发预警，但实际发生了冲击地压事件，FN表示预警指标触发预警，但实际未发生冲击地压事件。

进一步地，所述方法还包括：根据需求重新筛选预警指标，以适应工作环境，提高预警准确率。

进一步地，所述时序冲击地压危险量化预警指标T_Q的计算公式为：

其中，n为优选指标总数；F_k为第k个优选指标对应的预警效能；W_k表示第k个优选指标的异常隶属度，当指标触发预警时取值为1，未触发时取值为0；

基于所述时序冲击地压危险量化预警指标，判定矿井大区域范围内冲击危险等级，具体为：当0≤T_Q≤0.25时，判定矿井处于无冲击危险等级；当0.25<T_Q≤0.5时，判定矿井处于弱冲击危险等级；当0.5<T_Q≤0.75时，判定矿井处于中冲击危险等级；当0.75<T_Q≤1时，判定矿井处于强冲击危险等级。

进一步地，对矿井空间进行网格划分，计算各网格内优选指标值；基于各网格内优选指标值和各优选指标的预警效能，计算各网格内空间冲击地压危险指标，绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区，包括：

将矿井整体所在空间划分为等大小的六面体规则单元格，计算各网格内优选指标值；其中，所述单元格的边长不小于微震监测系统的定位误差；

对各网格内优选指标值进行单向归一化；其中，

正向预警指标的单向归一化方式为：

负向预警指标的单向归一化方式为：

以各优选指标的预警效能作为权重，计算各网格内空间冲击地压危险指标：

其中，

表示归一化结果，x_i表示第i个指标，x_min表示x_i的最小值，x_max表示x_i的最大值，

表示网格内第k个优选指标x_k的归一化结果，F_k表示x_k对应的预警效能，n表示优选指标的总数；

利用高斯空间光滑法绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区。

进一步地，所述利用高斯空间光滑法绘制空间冲击地压危险指标的空间分布云图，具体为：

以网格大小作为滑移步长，其值利用下式计算：

其中，

为第i个网格中S_Q对应的高斯光滑值；n_j为第i网格周围第j个网格中的S_Q值；△_ij为第i个网格与第j个网格之间的距离；c为相关距离。

进一步地，根据局部冲击危险区和区域危险等级采取相应防冲措施，包括：

依据时序冲击地压危险量化预警指标得出的矿井大范围危险等级，判定需要采取的措施类型；再结合空间冲击地压危险指标高斯光滑云图判定的局部危险位置，根据确定的措施类型实施防冲措施，以最终实现冲击地压危险的防治。

本发明提供的技术方案带来的有益效果至少包括：

本发明提供的冲击地压时空动态综合预警方法，首先利用Mann-Kendall趋势检验法判定冲击地压各前兆预警指标的时序变化趋势，以其趋势是否符合预警前兆表征规律作为是否触发预警的标准，之后评价各指标预警效能并进行指标优选，优选后的指标进行数据融合得时序冲击地压危险量化预警指标T_Q，以确定矿井大范围内冲击危险等级，再计算网格划分后空间内优选指标的值并单向归一化，利用高斯空间光滑法绘制空间冲击危险指数S_Q的云图，以结合矿图确定局部危险区，最终辅助井下相关人员进行准确、高效的防冲决策。本发明方法依据冲击地压前兆预警指标的时序变化信息和空间分布信息，首先从矿井区域范围内确定了整体冲击危险等级，之后结合矿图确定局部冲击危险区，针对局部范围内采取适应其危险等级的措施，极大限度得提高了井下防冲措施的实施效能并一定程度上降低了开采成本，同时该方法在现场实际监测数据的基础上定期更新，不断优选适应现场实际的冲击预警指标，适应于矿山井下复杂的开采条件，为冲击地压的精准防治提供了准确的决策依据以实现其精准预警。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明第一实施例提供的冲击地压时空动态综合预警方法的流程图；

图2是本发明第二实施例提供的冲击地压时空动态综合预警方法的流程图；

图3是本发明第二实施例提供的日总频次F_sum预警指标的时序变化曲线；

图4是本发明第二实施例提供的频次偏差值D_F预警指标的时序变化曲线；

图5是本发明第二实施例提供的平均总频次

预警指标的时序变化曲线；

图6是本发明第二实施例提供的缺震M_m预警指标的时序变化曲线；

图7是本发明第二实施例提供的A(b)值预警指标的时序变化曲线；

图8是本发明第二实施例提供的微震活动标度F预警指标的时序变化曲线；

图9是本发明第二实施例提供的算法复杂度AC预警指标的时序变化曲线；

图10是本发明第二实施例提供的P(b)值预警指标的时序变化曲线；

图11是本发明第二实施例提供的时间信息熵Q(t)预警指标的时序变化曲线；

图12是本发明第二实施例提供的微震活动度S_D预警指标的时序变化曲线；

图13是本发明第二实施例提供的震源集中程度λ预警指标的时序变化曲线；

图14是本发明第二实施例提供的日最大能量E_max预警指标的时序变化曲线；

图15是本发明第二实施例提供的日总能量E_sum预警指标的时序变化曲线；

图16是本发明第二实施例提供的日平均能量E_avg预警指标的时序变化曲线；

图17是本发明第二实施例提供的能量偏差值D_E预警指标的时序变化曲线；

图18是本发明第二实施例提供的平均总能量

预警指标的时序变化曲线；

图19是本发明第二实施例提供的断层总面积A(t)预警指标的时序变化曲线；

图20是本发明第二实施例提供的b值预警指标的时序变化曲线；

图21是本发明第二实施例提供的时序冲击危险指数T_Q的时序变化曲线；

图22是本发明第二实施例提供的冲击地压危险指数S_Q在冲击地压前后的云图演化情况示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

第一实施例

本实施例提供了一种冲击地压时空动态综合预警方法，该方法的执行流程如图1所示，该冲击地压时空动态综合预警方法包括以下步骤：

S1，基于微震现场实际监测数据，计算冲击地压前兆预警指标；

S2，判定各预警指标的时序变化趋势，以各预警指标的时序变化趋势是否符合预警前兆表征规律作为其是否触发预警的标准，并计算各预警指标的预警效能，筛选出预警效能大于预设预警效能阈值的预警指标，作为优选指标；

S3，将各优选指标进行数据融合，得到时序冲击地压危险量化预警指标，基于时序冲击地压危险量化预警指标，判定矿井大区域范围内冲击危险等级；

S4，对矿井空间进行网格划分，计算各网格内优选指标值；基于各网格内优选指标值和各优选指标的预警效能计算各网格内空间冲击地压危险指标，绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区；

S5，根据局部冲击危险区和区域危险等级采取针对性防冲措施。

综上，本实施例的冲击地压时空动态综合预警方法，首选从时间维度描述矿井区域范围内冲击地压危险前兆演化信息，然后结合高斯空间光滑法从空间维度确定矿井局部冲击危险区，辅助相关人员针对局部冲击危险区采取对应危险等级的防治措施，实现了冲击地压的精准预警并有效降低其防治成本。

第二实施例

本实施例提供了一种冲击地压时空动态综合预警方法，该方法的执行流程如图2所示，该冲击地压时空动态综合预警方法包括以下步骤：

S1，根据现场微震实时监测数据计算冲击地压前兆预警指标；

需要说明的是，在本实施例中，现场微震实时监测数据是指利用安装在矿井的微震监测系统(如ARAMIS M/E、SOS微震监测系统等)采集到的一系列独立微震事件的相关信息，包括各微震事件的发生时间、X坐标、Y坐标、Z坐标、能量、事件类型等；预警指标包括正向预警指标和负向预警指标，正向预警指标即冲击地压事件发生前出现增长趋势的指标，负向预警指标即冲击地压事件发生前出现降低趋势的指标；正向预警指标包括日最大能量、日总能量、日总频次、日平均能量、能量偏差值、频次偏差值、平均总能量、微震活动度、缺震、A(b)值、断层总面积A(t)、平均总频次、微震活动标度F、算法复杂度AC等；负向预警指标包括震源集中程度、b值、P(b)值、时间信息熵Q(t)等。

具体地，上述各预警指标的计算方式如下：

日总频次F_sum的取值由下式计算：

其中，n表示前24h内微震事件总频次；

频次偏差值D_F的取值由下式计算：

其中，F_t表示前24h内微震事件频次，

前30d内微震事件平均日频次；

平均总频次

的取值由下式计算：

其中，T表示时间窗长度，F_(t)表示t时刻的微震频次；

缺震M_m的取值由下式计算：

其中，m为能级分档总数；lgE_i为第i档能级；N_i为第i档能级实际微震数；

A(b)值的取值由下式计算：

其中，N为微震总数，M_i为微震事件能级；

微震活动标度F的取值由下式计算：

F₀＝10^6.11+1.09M，

其中，T表示时间窗长度，M_i为微震事件能级；

算法复杂度AC的取值由下式计算：

其中，n表示时间窗内微震能级的变化次数，M_max表示微震最大能级，M_min表示微震最小能级；

P(b)值的取值由下式计算：

其中，N为微震总数，M_i为微震事件能级；

时间信息熵Q(t)的取值由下式计算：

其中，n为时间窗内微震事件频次，t_i为第i个矿震发生时间；

微震活动度S_D的取值由下式计算：

其中，N为微震总数，E_i为微震事件能量，M为微震最大能级；

震源集中程度λ的取值由下式计算：

其中，λ₁、λ₂、λ₃为时间窗内微震事件坐标(x,y,z)组成协方差矩阵的特征根；日最大能量E_max的取值由下式计算：

E_max＝max(E₁,E₂,...,E_n)，

其中，E_i为前24h内第i个微震事件；

日总能量E_sum的取值由下式计算：

其中，n表示前24h内微震事件总频次，E_i表示各微震事件能量；

日平均能量E_avg的取值由下式计算：

其中，n表示前24h内微震事件总频次，E_i表示各微震事件能量；

能量偏差值D_E的取值由下式计算：

其中，E_t表示前24h内微震事件总能量，

前30d内微震事件平均日能量；平均总能量

的取值由下式计算：

其中，T表示时间窗长度，E_(t)表示t时刻的微震能量；

断层总面积A(t)的取值由下式计算：

其中，k₀为时间窗内微震的能级下限，k为各微震的能级，N(k)为时间窗内能级为k的微震事件数量；

b值取值由下式计算：

其中，m为能级分档总数；lgE_i为第i档能级；N_i为第i档能级实际微震数。

S2，利用Mann-Kendall趋势检验法判定预警指标的时序变化趋势；

需要说明的是，在本实施例中，上述Mann-Kendall趋势检验法，具体为结合微震实时监测数据，将各预警指标前一段时间内的计算结果组成时间序列，记为X[x₁,x₂,x₃...x_n](n为时间窗长度)，计算X的检验统计量S：

其中：

式中，n表示X[x₁,x₂,x₃...x_n]数据集的长度，x_k表示第k个指标，k＝1,2,3,...,i-1，x_i表示第i个指标，i＝k+1,2,3,...,n；

再计算X的检验标准量Z：

当Z>0时，相应的预警指标具有增长趋势；当Z<0时，相应的预警指标具有降低趋势；当Z＝0时，相应的预警指标不具有明显变化趋势。

S3，结合冲击地压预警指标前兆表征规律确定是否触发预警并评价各指标预警效能；

需要说明的是，在本实施例中，所述冲击地压预警指标前兆表征规律是指预警指标随冲击地压危险程度的变化而出现的响应特征，具体为日最大能量、日总能量、日总频次、日平均能量、能量偏差值、频次偏差值、平均总能量、微震活动度、缺震、A(b)值、断层总面积A(t)、平均总频次、微震活动标度F、算法复杂度AC等预警指标在冲击危险程度增加时出现的异常高值、连续增加、持续在高值附近波动等响应特征，震源集中程度、b值、P(b)值、时间信息熵Q(t)等预警指标在冲击危险程度增加时出现的异常低值、连续降低、持续在低值附近波动等响应特征；所述预警触发标准为，若正向前兆预警指标存在增长趋势、负向预警指标存在降低趋势，则触发预警，其它情况不触发预警；

所述预警效能的计算方法具体为：

其中，F表示预警指标的预警效能，TP表示预警指标触发预警，并且实际发生了冲击地压事件，FP表示预警指标未触发预警，但实际发生了冲击地压事件，FN表示预警指标触发预警，但实际未发生冲击地压事件。

S4，筛选预警效能高的冲击地压前兆指标；

需要说明的是，在本实施例中，筛选预警效能高的冲击地压前兆指标，具体为结合矿井实际设定预警效能临界值，将预警效能大于该临界值的指标作为模型的优选指标，并且该筛选过程在发生冲击地压事件、地质条件有重大改变、更换工作面等情况下需重新进行以适应井下复杂工作环境，提高预警准确率。

S5，对优选指标进行数据融合，得时序冲击地压危险量化预警指标T_Q，以此判定矿井大区域范围内冲击危险等级；

需要说明的是，在本实施例中，T_Q的计算公式为：

其中，n表示优选指标总数；F_k表示第k个优选指标对应的预警效能；W_k表示第k个优选指标的异常隶属度，其取值为0或1，当该指标触发预警时，取值为1，未触发时，取值为0；

基于所述时序冲击地压危险量化预警指标，判定矿井大区域范围内冲击危险等级，具体为：当0≤T_Q≤0.25时，判定矿井处于无冲击危险等级；当0.25<T_Q≤0.5时，判定矿井处于弱冲击危险等级；当0.5<T_Q≤0.75时，判定矿井处于中冲击危险等级；当0.75<T_Q≤1时，判定矿井处于强冲击危险等级。

S6，对空间进行网格划分，计算各网格内微震事件对应的前兆预警指标值；

需要说明的是，在本实施例中，对空间进行网格划分，具体为：将矿井整体所在空间划分为等大小的六面体规则单元格，其边长不小于微震监测系统的定位误差，一般可设置为10～30m；计算各网格内冲击前兆预警指标值是指将网格空间范围内发生的微震事件作为数据源，计算S1中优选出来的前兆预警指标。

S7，对各网格内预警指标值进行单向归一化；

需要说明的是，在本实施例中，对各网格内预警指标值进行单向归一化具体为：对于正向预警指标，采用如下单向归一化方式：

对于负向预警指标，采用如下单向归一化方式：

其中，

表示归一化结果，x_i表示第i个指标，x_min表示x_i的最小值，x_max表示x_i的最大值。

S8，以各优选指标的预警效能作为权重计算各网格内空间冲击地压危险指标S_Q；其中，空间冲击地压危险指标S_Q计算方式为：

其中，

表示网格内第k个优选指标x_k的归一化结果，F_k表示x_k对应的预警效能，n表示优选指标的总数。

S9，利用高斯空间光滑法绘制S_Q云图并结合底图判定矿井局部冲击危险区；

需要说明的是，在本实施例中，所述高斯空间光滑法，具体为以网格大小作为滑移步长，其值利用下式计算：

其中，

为第i个网格中S_Q对应的高斯光滑值；n_j为第i网格周围第j个网格中的S_Q值；△_ij为第i个网格与第j个网格之间的距离；c为相关距离。一般取i网格周围3c距离内网格的数据来计算

在执行完上述步骤后，即可根据局部冲击危险区和区域危险等级采取相应防冲措施，具体为：依据T_Q得出的矿井大范围危险等级，判定需要采取的措施类型，再结合S_Q高斯光滑云图判定的局部危险位置，实施防冲措施，以最终实现冲击地压危险的高效、精准防治。

下面，结合具体的应用场景，对本实施例的方法作进一步说明：

本实施例以某煤矿的某一工作面为例，首先搜集该工作面开采期间(2016年10月15日至2017年10月15日)的微震原始监测数据及相关信息：在工作面开采过程中共发生过21次冲击地压及大能量矿震事件(微震能量大于1E+5J)，对其进行预警，具体步骤如下：

(1)采集现场微震监测系统实时监测数据，以15天为时间窗口、1天为滑移步长对原始数据进行预处理，计算得冲击地压前兆预警指标的时间序列，其中正向预警指标包括：日最大能量、日总能量、日总频次、日平均能量、能量偏差值、频次偏差值、平均总能量、微震活动度、缺震、A(b)值、断层总面积A(t)、平均总频次、微震活动标度F、算法复杂度AC；负向预警指标包括：震源集中程度、b值、P(b)值、时间信息熵Q(t)；

(2)利用Mann-Kendall趋势检验法判断各指标变化趋势，结合其冲击地压前兆表征规律判断是否触发预警，结果如图3～图20所示，其中灰色背景表示指标触发预警；

(3)计算各指标预警效能，结果如下(按预警效能高低排列)：

1.A(b)值:0.42；

2.日总频次:0.39；

3.缺震:0.38；

4.平均总频次:0.38；

5.微震活动度:0.37；

6.时间信息熵Q(t):0.36；

7.b值:0.36；

8.断层总面积A(t):0.35；

9.微震活动标度F:0.32；

10.平均总能量:0.31；

11.日总能量:0.31；

12.日最大能量:0.28；

13.频次偏差值:0.26；

14.P(b)值:0.25；

15.算法复杂度AC:0.25；

16.震源集中程度:0.23；

17.日平均能量:0.22；

18.能量偏差值:0.17。

结合该矿井现场实际，将预警效能临界值设置为0.35，优选出以下7个高预警效能指标，分别为：A(b)值、日总频次、缺震、平均总频次、微震活动度、时间信息熵Q(t)、b值。

对优选的指标进行数据融合得时序冲击危险指数T_Q，T_Q的时序变化曲线如图21所示。T_Q的不同取值对应不同冲击地压危险等级，具体如表1所示。

表1冲击地压危险性分级

以“2016-11-24”冲击地压事件为例，在其发生前5d内T_Q达到了中等冲击危险，在事件发生当天更达到了高冲击危险，预警准确，井下相关人员需及时采取对应危险等级的防治措施。

(1)对矿井空间进行网格划分，令网格边长为20m，计算各网格内优选的预警指标值，之后将其单向归一化；

(2)结合(3)中各指标的预警效能，计算各网格内空间冲击危险指数S_Q；

(3)利用高斯空间光滑法计算空间各网格内S_Q的高斯光滑值，令c＝3；之后利用matplotlib第三方库中的contourf命令绘制云图，以“2016-11-24”冲击地压事件前后8天的云图为例，结果如图22所示；可以看出从11.19开始工作面+2100m位置周围100米左右颜色明显较其他区域深，存在应力集中区，并且随着时间推移，应力集中区范围增大，为冲击地压的发生累积了能量，因此应对此局部区域随时间采取中或高冲击危险等级对应的防治措施。

综上，本实施例的冲击地压时空动态综合预警方法，首先根据现场微震实时监测数据计算冲击地压前兆预警指标，利用Mann-Kendall趋势检验法判定预警指标的时序变化趋势，结合冲击地压预警指标前兆表征规律确定是否触发预警并评价各指标预警效能，筛选出预警效能高的冲击地压前兆指标；其次对优选出的预警指标进行数据融合得到时序冲击地压危险量化预警指标T_Q以判定矿井大区域范围内冲击危险等级；最后，对空间进行网格划分，计算各网格内前兆预警指标值并对其进行单向归一化，以各优选指标的预警效能作为权重计算各网格内空间冲击地压危险指标S_Q，再利用高斯空间光滑法绘制S_Q云图并结合底图判定矿井局部冲击危险区。本实施例的方法从时间和空间两个维度综合考虑冲击地压孕育演化过程的前兆信息，实现了“区域-局部”多尺度递进式监测预警，现场可根据预警的局部冲击危险区及危险等级采取对应的防治措施，这为冲击地压的精准防治提供了准确的决策依据并能有效降低其防治成本。

此外，需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。

最后需要说明的是，以上所述是本发明优选实施方式，应当指出，尽管已描述了本发明优选实施例，但对于本技术领域的技术人员来说，一旦得知了本发明的基本创造性概念，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明实施例范围的所有变更和修改。

Claims (10)

1.一种冲击地压时空动态综合预警方法，其特征在于，包括：

基于微震现场实际监测数据，计算冲击地压前兆预警指标；

判定各预警指标的时序变化趋势，以各预警指标的时序变化趋势是否符合预警前兆表征规律作为其是否触发预警的标准，并计算各预警指标的预警效能，筛选出预警效能大于预设预警效能阈值的预警指标，作为优选指标；

将所述优选指标进行数据融合，得到时序冲击地压危险量化预警指标，基于所述时序冲击地压危险量化预警指标，判定矿井大区域范围内冲击危险等级；

对矿井空间进行网格划分，计算各网格内优选指标值；基于各网格内优选指标值和各优选指标的预警效能，计算各网格内空间冲击地压危险指标，绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区；

根据局部冲击危险区和区域危险等级采取针对性防冲措施。

2.如权利要求1所述的冲击地压时空动态综合预警方法，其特征在于，所述微震现场实际监测数据指的是利用安装在矿井的微震监测系统采集到的一系列独立微震事件的相关信息，包括：各微震事件的发生时间、X坐标、Y坐标、Z坐标、能量以及事件类型；

所述预警指标包括正向预警指标和负向预警指标，所述正向预警指标指的是冲击地压事件发生前出现增长趋势的指标，所述负向预警指标指的是冲击地压事件发生前出现降低趋势的指标；其中，所述正向预警指标包括：日最大能量、日总能量、日总频次、日平均能量、能量偏差值、频次偏差值、平均总能量、微震活动度、缺震、A(b)值、断层总面积、平均总频次、微震活动标度和算法复杂度；所述负向预警指标包括：震源集中程度、b值、P(b)值和时间信息熵。

3.如权利要求2所述的冲击地压时空动态综合预警方法，其特征在于，所述判定各预警指标的时序变化趋势，以各预警指标的时序变化趋势是否符合预警前兆表征规律作为其是否触发预警的标准，并计算各预警指标的预警效能，筛选出预警效能大于预设预警效能阈值的预警指标，作为优选指标，包括：

利用Mann-Kendall趋势检验法判定各预警指标的时序变化趋势；

以各预警指标的时序变化趋势是否符合预警前兆表征规律作为其是否触发预警的标准，并计算各预警指标的预警效能；其中，所述预警前兆表征规律是指预警指标随冲击地压危险程度的变化而出现的响应特征，所述响应特征为：所述正向预警指标在冲击危险程度增加时出现的异常高值、连续增加或持续在高值附近波动，所述负向预警指标在冲击危险程度增加时出现的异常低值、连续降低或持续在低值附近波动；所述触发预警的标准为：若正向预警指标存在增长趋势或负向预警指标存在降低趋势，则触发预警，否则，不触发预警；

筛选出预警效能大于预设预警效能阈值的预警指标，作为优选指标。

4.如权利要求3所述的冲击地压时空动态综合预警方法，其特征在于，所述利用Mann-Kendall趋势检验法判定各预警指标的时序变化趋势，包括：

将各预警指标前一段时间内的计算结果组成时间序列，记为X[x₁,x₂,x₃,...,x_n]，其中，n表示时间窗长度；首先计算X的检验统计量S：

其中：

式中，x_k表示第k个指标，x_i表示第i个指标；

再计算X的检验标准量Z：

当Z>0时，相应的预警指标具有增长趋势；当Z<0时，相应的预警指标具有降低趋势；当Z＝0时，相应的预警指标不具有明显变化趋势。

5.如权利要求3所述的冲击地压时空动态综合预警方法，其特征在于，所述预警效能的计算公式为：

其中，F表示预警指标的预警效能，TP表示预警指标触发预警，并且实际发生了冲击地压事件，FP表示预警指标未触发预警，但实际发生了冲击地压事件，FN表示预警指标触发预警，但实际未发生冲击地压事件。

6.如权利要求1所述的冲击地压时空动态综合预警方法，其特征在于，所述方法还包括：根据需求重新筛选预警指标，以适应工作环境，提高预警准确率。

7.如权利要求1所述的冲击地压时空动态综合预警方法，其特征在于，所述时序冲击地压危险量化预警指标T_Q的计算公式为：

其中，n为优选指标总数；F_k为第k个优选指标对应的预警效能；W_k表示第k个优选指标的异常隶属度，当指标触发预警时取值为1，未触发时取值为0；

基于所述时序冲击地压危险量化预警指标，判定矿井大区域范围内冲击危险等级，具体为：当0≤T_Q≤0.25时，判定矿井处于无冲击危险等级；当0.25<T_Q≤0.5时，判定矿井处于弱冲击危险等级；当0.5<T_Q≤0.75时，判定矿井处于中冲击危险等级；当0.75<T_Q≤1时，判定矿井处于强冲击危险等级。

8.如权利要求2所述的冲击地压时空动态综合预警方法，其特征在于，对矿井空间进行网格划分，计算各网格内优选指标值；基于各网格内优选指标值和各优选指标的预警效能，计算各网格内空间冲击地压危险指标，绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区，包括：

将矿井整体所在空间划分为等大小的六面体规则单元格，计算各网格内优选指标值；其中，所述单元格的边长不小于微震监测系统的定位误差；

对各网格内优选指标值进行单向归一化；其中，

正向预警指标的单向归一化方式为：

负向预警指标的单向归一化方式为：

以各优选指标的预警效能作为权重，计算各网格内空间冲击地压危险指标：

其中，

表示归一化结果，x_i表示第i个指标，x_min表示x_i的最小值，x_max表示x_i的最大值，

表示网格内第k个优选指标x_k的归一化结果，F_k表示x_k对应的预警效能，n表示优选指标的总数；

利用高斯空间光滑法绘制空间冲击地压危险指标的空间分布云图，结合矿井底图判定局部冲击危险区。

9.如权利要求8所述的冲击地压时空动态综合预警方法，其特征在于，所述利用高斯空间光滑法绘制空间冲击地压危险指标的空间分布云图，具体为：

以网格大小作为滑移步长，其值利用下式计算：

其中，

为第i个网格中S_Q对应的高斯光滑值；n_j为第i网格周围第j个网格中的S_Q值；Δ_ij为第i个网格与第j个网格之间的距离；c为相关距离。

10.如权利要求1所述的冲击地压时空动态综合预警方法，其特征在于，所述根据局部冲击危险区和区域危险等级采取相应防冲措施，包括：

依据时序冲击地压危险量化预警指标得出的矿井大范围危险等级，判定需要采取的措施类型；再结合空间冲击地压危险指标高斯光滑云图判定的局部危险位置，根据确定的措施类型实施防冲措施，以最终实现冲击地压危险的防治。

Images