CN103306722B - 一种冲击危险区的微震多维信息综合区域探测评价方法 - Google Patents

Abstract

本发明公开了一种冲击危险区的微震多维信息综合区域探测评价方法，包括如下步骤：（1）矿震数据的空间坐标旋转及分级筛选；（2）建立微震多维信息探测评价指标体系；（3）空间统计滑移模型创建；（4）据下式计算微震多维信息综合探测评价指标值W；（5）根据微震多维信息综合探测评价指标值W确定危险等级和对应危险等级的防治对策；微震多维信息综合探测评价指标值W越大危险等级越高。本发明综合考虑了微震时空强三要素，并定量描述了冲击危险性强弱程度，能够实时反映当前时段内的冲击危险区域及等级和指导现场采取对应的防治对策，预测效能极高；建立的探测评价指标体系及综合指标物理意义明确，可操作性强，运算速度快。

Description

一种冲击危险区的微震多维信息综合区域探测评价方法

技术领域

本发明涉及一种冲击危险区的探测评价方法，尤其涉及一种冲击危险区的微震多维信息综合区域探测评价方法。

背景技术

目前，矿井冲击危险区域评价方法主要有：实验室冲击倾向性测定、综合指数法、数值模拟法、覆岩空间结构理论、钻屑法、应力监测法、电磁辐射法、声发射法和微震法等。其中，实验室冲击倾向性测定、综合指数法、数值模拟法和覆岩空间结构理论主要是预判采掘工作面发生强矿震或冲击矿压的危险性，若危险性较高，则采取相应的监测、预警与防治措施。钻屑法和应力监测法属于局部点监测，只能在巷道附近进行监测，虽然可通过多个点监测的数据插值获得监测参量的分布云图，但其评价范围仍有限。电磁辐射法与声发射法属于局部探测法，探测范围为工作面尺度，且只能绘制局部区域内统计参量的变化趋势，例如地音事件能量、频次、脉冲数和电磁辐射强度等，并根据这些参量的趋势变化进行预警，仍存在监测范围小的弊端。微震监测方法，能够对全矿范围进行实时监测，是一种区域性、及时监测手段，能够给出震动后的各种信息，具有不损伤煤体、劳动强度小、时间和空间连续等优点，该技术目前被公认为煤岩动力灾害，特别是对于煤矿动力现象监测最有效和最有发展潜力的监测方法之一。尤其是近年来基于微震监测系统提出的CT探测评价技术，很好地解决了冲击危险区域的探测评价问题，其唯一存在的缺点就是反演计算周期较长。另外，基于微震事件分布提出的微震活动性空间演化仅考虑了微震频次因素，以及基于微震事件分布提出的微震能量空间演化和应力等值线分布仅考虑了微震强度的因素。

由于冲击矿压发生的突发性和复杂性，仅仅依靠单一影响因子解决冲击危险区探测评价问题是不可能的，必须走综合探测评价的道路。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种冲击危险区的微震多维信息综合区域探测评价方法，是一种综合考虑微震时空强三要素、定量描述冲击危险性强弱程度、运算速度快、可操作性强的冲击危险区域探测评价方法。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种冲击危险区的微震多维信息综合区域探测评价方法，包括如下步骤：

（1）矿震数据的空间坐标旋转及分级筛选；

（2）建立微震多维信息探测评价指标体系，包括频次因子W₁、空间因子W₂和强度因子W₃，所述频次因子W₁包括频次指标W₁₁，所述空间因子W₂包括震源集中程度指标W₂₁，所述强度因子W₃包括最大应力当量指标W₃₁和总应力当量指标W₃₂；

（3）空间统计滑移模型创建；

（4）根据下式计算微震多维信息综合探测评价指标值W：

W=ω₁₁·W₁₁+ω₂₁·W₂₁+ω₃₁·W₃₁+ω₃₂·W₃₂

其中，ω₁₁,ω₂₂,ω₃₁,ω₃₂为相应指标的预测信度（即权重），且满足ω₁₁+ω₂₂+ω₃₁+ω₃₂=1，λ_ij(t)为相应指标在监测时段t内统计区域的异常隶属度，取值范围为0-1；

（5）根据微震多维信息综合探测评价指标值W确定危险等级和对应危险等级的防治对策；微震多维信息综合探测评价指标值W越大危险等级越高。

具体的，所述步骤（1）中，矿震数据的空间坐标旋转，采用下式计算：

$\left[\begin{array}{ccc}{x^{\′}}_1& {y^{\′}}_1& z_1^{\′}\\ {x^{\′}}_2& {y^{\′}}_2& z_2^{\′}\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ {x^{\′}}_i& {y^{\′}}_i& z_i^{\′}\end{array}\right]=\left[\begin{array}{ccc}{x}_1-x_0& y_1-y_0& z_1-z_0\\ x_2-x_0& y_2-y_0& z_2-z_0\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ x_i-x_0& y_i-y_0& z_i-z_0\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\α}& -\sin \mathrm{\α}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\\ \sin \mathrm{\α}& \cos \mathrm{\α}\cos \mathrm{\β}& -\cos \mathrm{\α}\cos \mathrm{\β}\\ 0& \sin \mathrm{\β}& \cos \mathrm{\β}\end{array}\right]$

其中，（x′_i，y′_i，z′_i)为转化后的相对坐标值，（x_i,y_i,z_i）为原始绝对坐标值，（x₀,y₀,z₀）为工作面下巷切眼绝对坐标值，α为绝对坐标系与相对坐标系的夹角，以绝对x坐标轴向相对x坐标轴逆时钟旋转为负，顺时针旋转为正，β为煤层夹角，与绝对坐标系正向夹角为正，反之为负。

具体的，所述步骤（1）中，根据矿震能量对矿震数据进行分级筛选，即矿震能量上、下限分别以Gurtenberg-Richter幂率曲线高能量端和低能量端偏离幂率的拐点作为能量上、下限的分界线。

具体的，所述步骤（2）中，频次指标W₁₁指统计区域内单位时间发生的矿震事件数；震源集中程度指标W₂₁计算式为：最大应力当量指标W₃₁计算式为：总应力当量指标W₃₂计算式为：其中，λ₁、λ₂、λ₃为协方差矩阵的特征根，E_max为统计区域内矿震事件的最大能量，E_i为统计区域内第i个矿震事件的能量，T为统计时间天数。

具体的，所述步骤（3）中，空间统计滑移模型创建时，网格划分间距S和统计滑移半径R满足如下关系：其中统计滑移半径R根据震源定位误差确定。

具体的，所述步骤（4）中，相应指标预测信度ω_ij的确定采用如下方案：当存在大量历史训练样本数据时，采用许绍燮（1989）提出的预测效能评分方法确定；当没有历史训练样本数据时，默认各指标预测信度ω_ij相等。

具体的，所述步骤（4）中的λ_ij(t)采用如下归一化方法计算：

对于正向异常指标的W₁₁、W₃₁和W₃₂：

λ_ij(t)=(Q_ij-Q_min)/(Q_max-Q_min)

对于负向异常指标的W₂₁：

λ_ij(t)=[(Q_max-Q_ij)/(Q_max-Q_min)]

其中，Q_ij为指标序列值；Q_max为指标序列最大值；Q_min为指标序列最小值。

具体的，所述步骤（5）中，根据微震多维信息综合探测评价指标值W确定危险等级和对应危险等级的防治对策，具体方案为：

当W<0.25时，危险等级为A，危险状态为无危险，防治对策为：所有的采掘工作可正常进行；

当0.25≤W<0.5时，危险等级为B，危险状态为弱危险，防治对策为：在采掘过程中，加强冲击矿压危险的监测预报；

当0.5≤W<0.75时，危险等级为C，危险状态为中等危险，防治对策为：在采掘工作的同时，采取强度弱化减冲治理措施，消除冲击危险；

当W≥0.75时，危险等级为D，危险状态为强危险，防治对策为：首先停止采掘作业，让人员撤离危险地点，然后采取强度弱化减冲治理措施，消除冲击危险。

有益效果：本发明提供的微震多维信息综合区域探测评价冲击危险方法，其创建的指标体系物理意义明确，综合考虑了微震时、空、强三要素，定量描述了冲击危险性强弱程度；且该方法运算速度快，适于编程实现智能化实时探测评价；同时本发明能够实时反映当前统计时段内的冲击危险区域及等级和指导现场采取对应的防治对策。

附图说明

图1为矿震数据空间坐标旋转示意图；其中(a)为平面旋转关系示意图，(b)为剖面旋转关系示意图；

图2为微震多维信息区域探测评价指标体系；

图3为空间统计滑移模型示意图；

图4为实例的矿震数据空间坐标旋转效果图；

图5为实例的矿震分布区的误差分布图；

图6为微震多维信息区域探测评价指标云图；其中(a)为频次指标W₁₁云图，(b)为最大应力当量指标W₃₁云图，(c)为总应力当量指标W₃₂云图；图中实心五角星为统计时间段内发生的11次危险性矿震；

图7为微震多维信息综合探测评价指标W云图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

一种冲击危险区的微震多维信息综合区域探测评价方法，综合考虑了微震时、空、强三要素，并定量描述了冲击危险性强弱程度，能够实时反映当前统计时段内的冲击危险区域及等级和指导现场采取对应的防治对策，且该指标物理意义明确，可操作性强，运算速度快，适于编程实现智能化实时探测评价。

关于矿震数据空间坐标旋转的说明。

由于微震监测系统监测到的位置坐标是大地绝对坐标，当将矿震事件向各绝对坐标轴上投影时，都会造成矿震事件的不利叠加，不利于分析沿工作面方向、推进方向和深度方向剖面上的震源分布规律，同时为了更直观地反映矿震事件与固定工作面的相对位置，因此，在将震源沿工作面方向、推进方向及深度方向剖面投影之前，必须将所有矿震数据的x、y、z值进行线性变换，即采用相对坐标系，如图1所示，具体采用下式计算：

$\left[\begin{array}{ccc}{x^{\&prime;}}_1& {y^{\&prime;}}_1& z_1^{\&prime;}\\ {x^{\&prime;}}_2& {y^{\&prime;}}_2& z_2^{\&prime;}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ {x^{\&prime;}}_i& {y^{\&prime;}}_i& z_i^{\&prime;}\end{array}\right]=\left[\begin{array}{ccc}{x}_1-x_0& y_1-y_0& z_1-z_0\\ x_2-x_0& y_2-y_0& z_2-z_0\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x_i-x_0& y_i-y_0& z_i-z_0\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\&alpha;}& -\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& \sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}\\ \sin \mathrm{\&alpha;}& \cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}& -\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ 0& \sin \mathrm{\&beta;}& \cos \mathrm{\&beta;}\end{array}\right]---\left(1\right)$

其中，(x′_i，y′_i，z′_i)为转化后的相对坐标值，（x_i,y_i,z_i）为原始绝对坐标值，（x₀,y₀,z₀）为工作面下巷切眼绝对坐标值，α为绝对坐标系与相对坐标系的夹角，以绝对x坐标轴向相对x坐标轴逆时钟旋转为负，顺时针旋转为正，β为煤层夹角，与绝对坐标系正向夹角为正，反之为负。

关于矿震数据分级筛选方法的说明。

通常，煤岩体发生冲击矿压（危险性矿震）前一定会在应力作用下产生众多小能量级别的矿震（也可称为“微震”），两者之间具有伴生关系，而这些小矿震则是研究并预警冲击矿压或危险性矿震的重要信息源，所以监测预警分析时应剔除已发生过的危险性矿震，即认为危险矿震事件能量为上限。同时，由于矿震事件的监测受到仪器灵敏度、记录条件、台网控制能力等影响，微震观测和数据处理能力有限，即存在一个能量下限。本发明以Gurtenberg-Richter幂率曲线高能量端和低能量端偏离幂率的拐点分别作为能量上、下限的分界线。

关于微震多维信息区域探测评价指标体系的说明。

综合考虑微震时、空、强三要素，建立如图2所示的微震多维信息区域探测评价指标体系。当中频次指标W₁₁指统计区域内单位时间发生的矿震事件数；震源集中程度指标W₂₁计算式为：最大应力当量指标W₃₁计算式为：总应力当量指标W₃₂计算式为：其中，λ₁、λ₂、λ₃为协方差矩阵的特征根，E_max为统计区域内矿震事件的最大能量，E_i为统计区域内第i个矿震事件的能量，T为统计时间天数。

值得注意的是，由于本发明建立的空间统计滑移模型中的统计区域较小，为一固定形状，且统计面积相等，此时，统计区域的频次基本反映了震源分布的集中程度，即频次越大，震源集中程度越高。因此，作为近似，实际操作时，空间因子可不计入。

关于空间统计滑移模型建立的说明。

在地震研究中，Frankel（1995）在美国中东部地震区划中首次使用了空间光滑地震活动性的方法，其最显著的特点是不根据地震构造来划分潜在震源区，而直接使用根据地震目录进行空间光滑后的点源来进行地震危险性计算。由于该法简单易行，在世界各国的地震危险性分析中得到了广泛应用。在实际应用时，上述方法存在一重要参数需要确定，即统计滑移半径R。对于统计滑移半径的确定，Lapajne（2003）根据定位误差确定，徐伟进（2012）采用核密度估计确定，Pelaez Montilla（2003）、Hagos（2006）等采用随意选取一个值作为统计滑移半径。本发明采用定位误差确定，其中定位误差采用巩思园（2010）提出的台网定位误差分析方法确定。为避免统计滑移过程中遗漏个别矿震事件而导致最终结果的失真，网格划分间距S与统计滑移半径R需满足如下关系：最终，空间统计滑移模型示意图如图3所示。

关于微震多维信息综合区域探测评价指标判别准则构建的说明。

在矿井监测区域内，在一定的时间内，已进行了一定的矿震观测。在这种情况下，就可以通过微震多维信息综合区域探测评价指标，对当前的冲击危险等级进行预警。冲击矿压危险程度分为四级，根据不同的危险程度，可采用相应的防治措施，见表1。

表1 对应危险等级的防治措施

下面结合实例对本发明做出进一步的说明。

实例背景分析。

实例分析选取河南义马跃进煤矿25110工作面回采过二次见方及断层危险区期间（2011-5-1～2011-10-1）的微震监测数据，总共得到25110工作面微震事件757个。

按照本发明思想实施本发明：

（1）选取实例参数：α=-25°，β=12°，x₀=37580326.8850，y₀=3840377.4966，z₀=-439.372，采用式（1）对实例矿震数据进行空间旋转，效果如图4所示；

（2）采用本发明方法求出实例矿震数据的上下限：下限为10^1.33J，上限为10^6.93J，进行矿震数据的分级筛选，得出矿震事件701个；

（3）采用台网误差分析，得出实例矿震区域的误差分布云图，如图5所示，结合图4的矿震数据分布图，可知震源分布区的最大误差为30m左右，因此，取统计滑移半径R=30m，网格划分间距S=42m；

（4）不计空间因子，绘制出微震多维信息区域探测评价指标云图，如图6所示；

（5）赋予各指标相等预测信度，即{ω₁₁,ω₃₁,ω₃₂}={0.5,0.25,0.25}，绘制出微震多维信息综合探测评价指标云图，如图7所示。

实例表明，通过分析发生过的11次危险性矿震，得出各指标的预测情况如下：频次指标：强危险区域2次，中等危险区域0次，弱危险区域1次，无危险区域8次；最大应力当量指标：强危险区域2次，中等危险区域3次，弱危险区域3次，无危险区域3次；总应力当量指标：强危险区域3次，中等危险区域2次，弱危险区域1次，无危险区域5次；综合指标：强危险区域3次，中等危险区域3次，弱危险区域3次，无危险区域2次。综上所述，本发明建立的微震多维信息区域探测评价指标分别从不同程度上反映了不同区域的冲击危险程度，并与发生的危险性矿震对应效果较好，最后建立的综合探测评价指标综合考虑了时、空、强三要素，预测效果最好。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种冲击危险区的微震多维信息综合区域探测评价方法，包括如下步骤：

(1)矿震数据的空间坐标旋转及分级筛选；

(2)建立微震多维信息探测评价指标体系，包括频次因子W₁、空间因子W₂和强度因子W₃，所述频次因子W₁包括频次指标W₁₁，所述空间因子W₂包括震源集中程度指标W₂₁，所述强度因子W₃包括最大应力当量指标W₃₁和总应力当量指标W₃₂；

(3)空间统计滑移模型创建；

(4)根据下式计算微震多维信息综合探测评价指标值W：

W＝ω₁₁·W₁₁+ω₂₁·W₂₁+ω₃₁·W₃₁+ω₃₂·W₃₂

其中，ω₁₁,ω₂₁,ω₃₁,ω₃₂为相应指标的预测信度，且满足ω₁₁+ω₂₁+ω₃₁+ω₃₂＝1，λ_ij(t)为相应指标在监测时段t内统计区域的异常隶属度，取值范围为0-1；

(5)根据微震多维信息综合探测评价指标值W确定危险等级和对应危险等级的防治对策；微震多维信息综合探测评价指标值W越大危险等级越高；

所述步骤(1)中，矿震数据的空间坐标旋转，采用下式计算：

$\left[\begin{array}{ccc}{x^{\&prime;}}_1& {y^{\&prime;}}_1& z_1^{\&prime;}\\ {x^{\&prime;}}_2& {y^{\&prime;}}_2& z_2^{\&prime;}\\ .& .& .\\ .& .& .\\ .& .& .\\ {x^{\&prime;}}_i& {y^{\&prime;}}_i& z_i^{\&prime;}\end{array}\right]=\left[\begin{array}{ccc}{x}_1-x_0& y_1-y_0& z_1-z_0\\ x_2-x_0& y_2-y_0& z_2-z_0\\ .& .& .\\ .& .& .\\ .& .& .\\ x_i-x_0& y_i-y_0& z_i-z_0\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\&alpha;}& -\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}& \sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}\\ \sin \mathrm{\&alpha;}& \cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}& -\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ 0& \sin \mathrm{\&beta;}& \cos \mathrm{\&beta;}\end{array}\right]$

其中，(x'_i,y'_i,z_′i)为转化后的相对坐标值，(x_i,y_i,z_i)为原始绝对坐标值，(x₀,y₀,z₀)为工作面下巷切眼绝对坐标值，α为绝对坐标系与相对坐标系的夹角，以绝对x坐标轴向相对x坐标轴逆时针旋转为负，顺时针旋转为正，β为煤层夹角，与绝对坐标系正向夹角为正，反之为负；

所述步骤(1)中，根据矿震能量对矿震数据进行分级筛选，即矿震能量上、下限分别以Gurtenberg-Richter幂率曲线高能量端和低能量端偏离幂率的拐点作为能量上、下限的分界线；

所述步骤(2)中，频次指标W₁₁指统计区域内单位时间发生的矿震事件数；震源集中程度指标W₂₁计算式为：最大应力当量指标W₃₁计算式为：总应力当量指标W₃₂计算式为：其中，Q_ij为指标序列值，是用于计算W_ij的原始物理意义表达式，λ₁、λ₂、λ₃为震源坐标参量(x，y，z)的协方差矩阵的特征根，E_max为统计区域内矿震事件的最大能量，E_i为统计区域内第i个矿震事件的能量，T为统计时间天数；

所述步骤(3)中，空间统计滑移模型创建时，网格划分间距S和统计滑移半径R满足如下关系：其中统计滑移半径R根据震源定位误差确定；

所述步骤(4)中，相应指标预测信度ω_ij的确定采用如下方案：当存在历史训练样本数据时，采用预测效能评分方法确定；当没有历史训练样本数据时，默认各指标预测信度ω_ij相等；

所述步骤(4)中，λ_ij(t)采用如下归一化方法计算：

对于正向异常指标的W₁₁、W₃₁和W₃₂：

λ_ij(t)＝(Q_ij-Q_min)/(Q_max-Q_min)

对于负向异常指标的W₂₁：

λ_ij(t)＝[(Q_max-Q_ij)/(Q_max-Q_min)]

其中，Q_ij为指标序列值；Q_max为指标序列最大值；Q_min为指标序列最小值；

所述步骤(5)中，根据微震多维信息综合探测评价指标值W确定危险等级和对应危险等级的防治对策，具体方案为：

当W<0.25时，危险等级为A，危险状态为无危险，防治对策为：所有的采掘工作可正常进行；

当0.25≤W<0.5时，危险等级为B，危险状态为弱危险，防治对策为：在采掘过程中，加强冲击矿压危险的监测预报；

当0.5≤W<0.75时，危险等级为C，危险状态为中等危险，防治对策为：在采掘工作的同时，采取强度弱化减冲治理措施，消除冲击危险；

当W≥0.75时，危险等级为D，危险状态为强危险，防治对策为：首先停止采掘作业，让人员撤离危险地点，然后采取强度弱化减冲治理措施，消除冲击危险。