一种冲击矿压的微震多维信息综合时序预警方法

1.一种冲击矿压的微震多维信息综合时序预警方法，其特征在于：包括如下步骤：
(1)对矿震数据进行分区分级筛选；
(2)建立微震多维信息指标体系，包括频次时序指标W1、空间时序指标W2和强度时序指标W3，所述频次时序指标W1包括频次指标W11，所述空间时序指标W2包括震源集中程度指标W21，所述强度时序指标W3包括最大应力当量指标W31和平均应力当量指标W32；
(3)根据下式计算微震多维信息综合时序预警指标值W：
W＝W1+W2+W3，W1＝ω11·W11，
W2＝ω21·W21，W3＝ω31·W31+ω32·W32
其中，ω11,ω21,ω31,ω32为相应指标的预测信度，且满足ω11+ω21+ω31+ω32＝1，λij(t)为相应指标在监测时段t内的异常隶属度，取值范围为0-1；
(4)根据微震多维信息综合时序预警指标值W确定危险等级和对应危险等级的防治对策；微震多维信息综合时序预警指标值W越大危险等级越高；
所述步骤(1)中，采用任意多边形分区方法对矿震数据进行分区筛选；
所述步骤(1)中，根据矿震能量对矿震数据进行分级筛选，即矿震能量上、下限分别以Gurtenberg-Richter幂率曲线高能量端和低能量端偏离幂率的拐点作为能量上、下限的分界线，所述拐点根据下式计算识别：
能量上限分界线根据下式计算识别：
其中，i＝2,3,…,n，Ri-High为Gurtenberg-Richter幂率曲线中横坐标lg Ei对应的相关系数值，参与计算的x序列为lg Ei-1,lg Ei,lg Ei+1,…,lg En，y序列为lg Ni-1,lg Ni,lg Ni+1,…,lg Nn， 为x序列的样本均值， 为y序列样本均值；
当 时，能量上限的分界线为
Gurtenberg-Richter幂率曲线中纵坐标Rk-High对应横坐标值lgEk中的能量Ek；
能量下限分界线根据下式计算识别：
其中，i＝2,3,…,n，Ri-Low为Gurtenberg-Richter幂率曲线中横坐标lgEi对应的相关系数值，参与计算的x序列为lg E1,…,lg Ei-1,lgEi，y序列为lg N1,…,lgNi-1,lg Ni，为x序列的样本均值， 为y序列样本均值；
当 时，能量下限的分界线为Gurtenberg-Richter幂
率曲线中纵坐标Rk-Low对应横坐标值lgEk中的能量Ek；
所述步骤(2)中，频次指标W11指统计时段内发生的矿震事件数；震源集中程度指标W21计算式为： 最大应力当量指标W31计算式为： 平均应力当
量指标W32计算式为： 其中，λ1、λ2、λ3为震源坐标参量(x，y，z)的协方差矩阵的特征根，Emax为统计时段内矿震事件的最大能量，Eavg为统计时段内矿震事件的平均能量；
所述步骤(3)中，相应指标的预测信度ωij的确定采用如下方案：当存在历史训练样本数据时，采用预测效能评分R值确定；当没有历史训练样本数据时，默认各指标的预测信度ωij相等，即
所述步骤(3)中，λij(t)采用如下归一化方法计算：
对于正向异常指标的W11、W31和W32：
λij(t)＝(Qij-Qmin)/(Qmax-Qmin)
对于负向异常指标的W21：
λij(t)＝[(Qmax-Qij)/(Qmax-Qmin)]
其中，Qij为指标序列值；Qmax为指标序列最大值；Qmin为指标序列最小值。
2.根据权利要求1所述的冲击矿压的微震多维信息综合时序预警方法，其特征在于：
所述步骤(4)中，根据微震多维信息综合时序指标值W确定危险等级和对应危险等级的防治对策，具体方案为：
当W<0.25时，危险等级为A，危险状态为无危险，防治对策为：所有的采掘工作正常进行；
当0.25≤W<0.5时，危险等级为B，危险状态为弱危险，防治对策为：在采掘过程中，加强冲击矿压危险的监测预报；
当0.5≤W<0.75时，危险等级为C，危险状态为中等危险，防治对策为：在采掘工作的同时，采取强度弱化减冲治理措施，消除冲击危险；
当W≥0.75时，危险等级为D，危险状态为强危险，防治对策为：首先停止采掘作业，让人员撤离危险地点，然后采取强度弱化减冲治理措施，消除冲击危险。

一种冲击矿压的微震多维信息综合时序预警方法\n技术领域\n[0001] 本发明涉及一种冲击矿压的微震多维信息综合时序预警方法，是一种微震多维信息综合时序监测预警冲击危险新指标，属于一种冲击矿压时序监测预警技术。\n背景技术\n[0002] 冲击矿压灾害是一种开采诱发的矿山地震，不仅造成井巷破坏、人员伤害、地面建筑物破坏，而且会引起瓦斯、煤尘爆炸。由于这种灾害发生时间、地点、位置等的复杂多样性和突发性，对其防治，特别是预测是世界性的难题。\n[0003] 目前，冲击矿压的监测预警方法主要包括：微震监测、电磁辐射监测、声发射监测、钻屑量监测、工作面矿压监测、采动应力监测等。尤其是微震监测方法，能够对全矿范围进行实时监测，是一种区域性、及时监测手段，能够给出震动后的各种信息，具有不损伤煤体、劳动强度小、时间和空间连续等优点。该技术目前被公认为煤岩动力灾害，特别是对于煤矿动力现象监测最有效和最有发展潜力的监测方法之一。在此基础上，国内外学者提出了大量微震监测预警指标，如能量、频次、参数b值、A(b)、P(b)、η值、时空扩散性、赫斯特指数、能量指数EI、视在应力、视在体积、地震黏度、空间集中度、蠕变等。然而上述指标仅单方面反映了微震的时、空、强三要素，同时各指标对异常的分析也仅停留在定性层面，如“指标曲线上升异常”，“下降异常”、“低值异常”、“高值异常”等。\n发明内容\n[0004] 发明目的：为了克服现有技术中存在的不足，本发明提供一种冲击矿压的微震多维信息综合时序预警方法，是一种综合考虑微震时空强三要素、定量描述冲击危险性强弱程度、可操作性强的冲击危险监测预警综合新指标。\n[0005] 技术方案：为实现上述目的，本发明采用的技术方案为：\n[0006] 一种冲击矿压的微震多维信息综合时序预警方法，包括如下步骤：\n[0007] （1）对矿震数据进行分区分级筛选；\n[0008] （2）建立微震多维信息指标体系，包括频次时序指标W1、空间时序指标W2和强度时序指标W3，所述频次时序指标W1包括频次指标W11，所述空间时序指标W2包括震源集中程度指标W21，所述强度时序指标W3包括最大应力当量指标W31和平均应力当量指标W32；\n[0009] （3）根据下式计算微震多维信息综合时序预警指标值W：\n[0010] W=W1+W2+W3，W1=ω11·W11，\n[0011] W2=ω21·W21，W3=ω31·W31+ω32·W32\n[0012] 其 中，ω11,ω22,ω31,ω32为相应指标的预测信度（即权重），且满足ω11+ω22+ω31+ω32=1， λij(t)为相应指标在监测时段t内的异常隶属度，取值范围为0-1；\n[0013] （4）根据微震多维信息综合时序预警指标值W确定危险等级和对应危险等级的防治对策；微震多维信息综合时序预警指标值W越大危险等级越高。\n[0014] 具体的，所述步骤（1）中，采用任意多边形分区方法对矿震数据进行分区筛选。\n[0015] 具体的，所述步骤（1）中，根据矿震能量对矿震数据进行分级筛选，即矿震能量上、下限分别以Gurtenberg-Richter幂率曲线高能量端和低能量端偏离幂率的拐点作为能量上、下限的分界线，所述拐点根据下式计算识别：\n[0016] 能量上限分界线根据下式计算识别：\n[0017] \n[0018] 其中，i=2,3,…,n，Ri-High为Gurtenberg-Richter幂率曲线中横坐标lgEi对应的相关系数值，参与计算的x序列为lgEi-1,lgEi,lgEi+1,…,lgEn，y序列为lgNi-1,lgNi,lgNi+1,…,lgNn， 为x序列的样本均值， 为y序列样本均值；\n[0019] 当 时，能量上限的分界线为\nGurtenberg-Richter幂率曲线中纵坐标Rk-High对应横坐标值lgEk中的能量Ek；\n[0020] 能量下限分界线根据下式计算识别：\n[0021] \n[0022] 其中，i=2,3,…,n，Ri-Low为Gurtenberg-Richter幂率曲线中横坐标lgEi对应的相关系数值，参与计算的x序列为lgE1,…,lgEi-1,lgEi，y序列为lgN1,…,lgNi-1,lgNi， 为x序列的样本均值， 为y序列样本均值；\n[0023] 当 时，能量下限的分界线为\nGurtenberg-Richter幂率曲线中纵坐标Rk-Low对应横坐标值lgEk中的能量Ek。\n[0024] 具体的，所述步骤（2）中，频次指标W11指统计时段内发生的矿震事件数；震源集中程度指标W21计算式为： 最大应力当量指标W31计算式为：\n平均应力当量指标W32计算式为： 其中，λ1、λ2、λ3为协方差矩\n阵的特征根，Emax为统计时段内矿震事件的最大能量，Eavg为统计时段内矿震事件的平均能量。\n[0025] 具体的，相应指标的预测信度ωij的确定采用如下方案：当存在大量历史训练样本数据时，采用预测效能评分R值确定；当没有历史训练样本数据时，默认各指标的预测信度ωij相等，即\n[0026] 具体的，所述步骤（3）中，λij(t)采用如下归一化方法计算：\n[0027] 对于正向异常指标的W11、W31和W32：\n[0028] λij(t)=(Qij-Qmin)/(Qmax-Qmin)\n[0029] 对于负向异常指标的W21：\n[0030] λij(t)=[(Qmax-Qij)/(Qmax-Qmin)]\n[0031] 其中，Qij为指标序列值；Qmax为指标序列最大值；Qmin为指标序列最小值。\n[0032] 具体的，所述步骤（4）中，根据微震多维信息综合时序指标值W确定危险等级和对应危险等级的防治对策，具体方案为：\n[0033] 当W<0.25时，危险等级为A，危险状态为无危险，防治对策为：所有的采掘工作正常进行；\n[0034] 当0.25≤W<0.5时，危险等级为B，危险状态为弱危险，防治对策为：在采掘过程中，加强冲击矿压危险的监测预报；\n[0035] 当0.5≤W<0.75时，危险等级为C，危险状态为中等危险，防治对策为：在采掘工作的同时，采取强度弱化减冲治理措施，消除冲击危险；\n[0036] 当W≥0.75时，危险等级为D，危险状态为强危险，防治对策为：首先停止采掘作业，让人员撤离危险地点，然后采取强度弱化减冲治理措施，消除冲击危险。\n[0037] 有益效果：本发明提供的冲击矿压的微震多维信息综合时序监测预警方法，其创建的指标体系物理意义明确，综合考虑了微震时、空、强三要素，定量描述了冲击危险性强弱程度；且该方法中涉及的各项参数均可以采用自动识别算法实现，适于编程实现智能化实时预警；同时本发明能够实时反映当前统计区域内的冲击危险状态和指导现场采取对应的防治对策。\n附图说明\n[0038] 图1为微震多维信息指标体系示意图；\n[0039] 图2为矿震数据分级筛选中能量上下限求解示意图；\n[0040] 图3为实例中的矿震分布图；\n[0041] 图4为微震多维信息指标曲线图；其中(a)为频次指标W11曲线图，(b)为震源集中程度指标W21曲线图，(c)为最大应力当量指标W31曲线图，(d)为平均应力当量指标W32曲线图；\n[0042] 图5为微震多维信息综合时序指标W曲线图。\n具体实施方式\n[0043] 下面结合附图对本发明作更进一步的说明。\n[0044] 一种冲击矿压的微震多维信息综合时序预警方法，具体为一种多维信息异常综合预警新指标，综合考虑了微震时、空、强三要素，并定量描述了冲击危险性强弱程度，能够实时反映当前统计区域内的冲击危险状态和指导现场采取对应的防治对策，且该指标物理意义明确，可操作性强，适于编程实现智能化实时预警。\n[0045] 关于矿震数据分区筛选方法的说明。\n[0046] 由于微震监测系统记录的是全矿井尺度内发生的矿震信号，而每个矿区可能又不止一个生产区域，各生产区域在不同地质条件下产生震源的机制不尽相同，并且即使同一工作面在不同时期也存在规律上的变化，所以非常有必要对监测区域在监测时段内的矿震数据进行分区筛选。同时，考虑到现场工作面形状近似矩形，并且一般情况下工作面布置位置与绝对坐标轴存在一定夹角，所以采用任意多边形分区较为合理。其分区原则建议采用文献（蔡武，2011，基于分区监测的矿震时空强演化规律分析，《煤矿安全》）给出的分区原则。\n[0047] 关于矿震数据分级筛选方法的说明。\n[0048] 通常，煤岩体发生冲击矿压（危险性矿震）前一定会在应力作用下产生众多小能量级别的矿震（也可称为“微震”），两者之间具有伴生关系，而这些小矿震则是研究并预警冲击矿压或危险性矿震的重要信息源，所以监测预警分析时应剔除已发生过的危险性矿震，即认为危险矿震事件能量为上限。同时，由于矿震事件的监测受到仪器灵敏度、记录条件、台网控制能力等影响，微震观测和数据处理能力有限，即存在一个能量下限。本发明以Gurtenberg-Richter幂率曲线高能量端和低能量端偏离幂率的拐点分别作为能量上、下限的分界线。\n[0049] 关于微震多维信息指标体系的说明。\n[0050] 综合考虑微震时、空、强三要素，建立如图1所示的微震多维信息指标体系。\n[0051] 频次时序指标W1：在给定的时空范围内，矿震频次越大则矿震活动性越强，冲击危险性越大，反之矿震活动性越弱，冲击危险性越小，矿震频次为零或接近于零，则说“平静”。\n[0052] 空间时序指标W2：在一定的研究范围内，当矿震分布很密（成丛条带分布）时，矿震活动性强，或说冲击危险性大，如果是一种正常的分散分布，则安全，矿震活动性低。\n[0053] 令∑为变量X1,X2,X3的协方差矩阵，X=(X1,X2,X3)T，u=(u1,u2,u3)T。考虑到(X-u)T -1 2\n∑ (X-u)=d（d为常数），为方便，设u=0，因此有\n[0054] \n[0055] 式中，λ1、λ2、λ3为协方差矩阵∑的特征根，Y1、Y2、Y3为特征根对应的主成分。可知式（1）是一个椭球方程。设变量X1,X2,X3遵从三元正态分布，则其概率密度函数为：\n[0056] \n[0057] 式中，|∑|为协方差矩阵∑的行列式。\n[0058] 很明显，式（1）为三元正态分布的等概率密度椭球曲面，即椭球体积越大，说明椭球表面处样本出现的概率越小，分布的离散程度越高；反之，椭球表面处样本出现的概率越大，集中程度越高。\n[0059] 因此，在三维空间可采用等概率密度椭球的体积 来反映矿震事\n件分布的震源集中程度，通过消除常量及量纲影响，得出震源集中程度指标W21为：\n[0060] \n[0061] 强度时序指标W3：除矿震频次及震源集中程度以外，矿震能量的大小也是一个重要指标。按照岩石力学原理，一个矿震事件被定义为在一定体积内的突然非弹性变形，该变形引起可检测的地震波，且每次矿震所释放能量的平方根与这次矿震发生前岩体内的应变呈正比，在弹性范围内应力和应变成正比。于是，矿震所释放能量的平方根就是冲击矿压发生前岩体内应力状态的一个测度。因此，可采用如下公式表示平均应力当量指标W32：\n[0062] \n[0063] 其中，Eavg为统计时段内矿震事件的平均能量。\n[0064] 在所讨论的时空范围内，如有两组矿震事件，其频度相同，平均能量也相同，但其最大能量仍可能不同。此时，可认为最大能量大的矿震事件组活动性强。因此，强度时序指标还应包含最大应力当量指标：\n[0065] \n[0066] 其中，Emax为统计时段内矿震事件的最大能量。值得注意的是，在有了指标Q32后，Q31的预测信度ω31不会太大。\n[0067] 关于公式 的说明。\n[0068] 在可靠性分析理论当中，指数分布函数描述了一种产品失效的分布：\n[0069] \n[0070] 其中，F(t)为失效分布函数，即产品寿命的分布函数。λ>0为产品的失效率，即工作到时刻t的产品，在时刻t后单位时间内发生失效的概率。将产品的失效比喻为出现异常矿震的概率，即失效率（异常概率）越高，产品失效（异常矿震发生）的可能性越大。进一步推导可得出适用于把各指标时序数据统一转换为异常指数时序数据的表达式：\n[0071] \n[0072] 其中，λij(t)为相应指标在监测时段t内的异常隶属度，取值范围为0-1。\n[0073] 具体λij(t)采用如下归一化方法计算：\n[0074] 对于正向异常指标的W11、W31和W32：\n[0075] λij(t)=(Qij-Qmin)/(Qmax-Qmin) （8）\n[0076] 对于负向异常指标的W21：λij(t)=[(Qmax-Qij)/(Qmax-Qmin)] （9）[0077] 其中，Qij为指标序列值；Qmax为指标序列最大值；Qmin为指标序列最小值。\n[0078] 关于微震多维信息综合时序指标判别准则构建的说明。\n[0079] 在矿井监测区域内，在一定的时间内，已进行了一定的矿震观测。在这种情况下，就可以通过微震多维信息综合时序预警指标，对当前的冲击矿压危险等级进行预警。冲击矿压危险程度分为四级，根据不同的危险程度，可采用相应的防治措施，见表1。\n[0080] 表1对应危险等级的防治措施\n[0081] \n[0082] \n[0083] 下面结合实例对本发明作出进一步的说明。\n[0084] 实例背景分析。\n[0085] 实例分析选取河南义马跃进煤矿25110工作面回采过二次见方及断层危险区期间（2011-5-1～2011-10-1）的微震监测数据，总共得到25110工作面微震事件757个，其\n6.93\n中能量上限10 J以上事件11个，其分布如图2所示。另外，期间共发生有记录的冲击矿压显现4次，当中包括一次工作面透水事件，见表2。\n[0086] 表2实例的背景数据\n[0087] \n[0088] 按照本发明思想实施本发明：\n[0089] （1）采用本发明方法求出实例矿震数据的上下限，如图3所示，即下限为101.33J，上\n6.93\n限为10 J，然后结合任意多边形分区方法，进行矿震数据的分区分级筛选，得出矿震事件\n701个；\n[0090] （2）采用5天时间窗，1天滑移步长，根据式（3）-（9）绘制出微震多维信息指标曲线，如图4所示；\n[0091] （3）赋予各指标相等预测信度，即 绘制出微震多\n维信息综合时序指标曲线，如图5所示。\n[0092] 为便于说明，将需预测的10次危险性矿震依次按①-⑩编号。实例表明，频次时序指标仅报准④号矿震，但它是四项指标体系中唯一报准④号矿震的指标；震源集中程度指标报准①②⑦⑧⑨⑩号矿震；最大应力当量指标报准③⑦⑧⑨⑩号矿震；平均应力当量指标报准⑦⑧⑨⑩号矿震；微震多维信息综合时序指标报准③④⑤⑦⑧⑨⑩号矿震。分析表明，微震多维信息综合时序指标报准率最高，达70%，综合考虑所有指标时，除⑥号矿震未报出异常外，其余9次均出现异常信息，报准率达90%。因此，建议在实际现场应用中，当指标体系中四项指标和综合指标中任意一项指标出现异常时，就应引起注意并及时进行现场危险性排查。综上所述，本发明建立的微震多维信息时序指标体系及综合指标考虑信息全面，预测效能极高。\n[0093] 以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。