基于微震群震动波能量衰减特性的冲击地压危险预测方法

刘耀琪，曹安业,2,3，王崧玮，杨 耀，郭文豪，薛成春，白贤栖，李 震

(1.中国矿业大学 矿业工程学院，江苏 徐州 221116；2.中国矿业大学 江苏省矿山地震监测工程实验室，江苏 徐州 221116；3.徐州物硕信息技术有限公司，江苏 徐州 221116；4.陕西小保当矿业有限公司，陕西 榆林 719302)

我国煤矿开采深度持续增加，地质条件与开采布局等愈加复杂，冲击地压灾害风险日益增加，成为制约矿井安全高效生产的主要灾害之一。如何建立有效的监测预警方法对降低冲击地压风险以及冲击地压防治等方面尤为重要。

煤岩体破裂向外辐射震动波，通过微震系统可进行捕捉并记录为微震事件。微震事件可定义为煤岩体一定体积内的突然非弹性变形，伴随于采掘空间周围，通过微震监测技术可实现冲击地压等动力灾害的监测预警。诸多学者依据微震特性在冲击地压预测预警方面进行了尝试，如震动波波速层析成像技术、值、震源参数、冲击变形能等。现场应用方面，丛利等建立了一种基于动静载叠加的综合预警方法，并将其应用于强矿震的预警。WANG等建立了基于微震定位误差的聚类分析方法，并将其用于冲击危险区域的预测，初步确定了预警指标。HE等建立了动载应力集中系数以定量评价矿震对煤岩体的动载扰动程度。田向辉等分析了某冲击地压频发工作面的微震能量、频次变化特征，并建立了一种定量-趋势冲击危险预警方法，有效提高了预警效能。李宏艳等基于统计损伤力学原理建立了以响应能量异常系数和无响应时间异常系数为指标的冲击危险性动态评价技术，可实现对现场冲击危险的实时连续评价。然而，由于煤矿井下条件复杂，采掘扰动、顶板破断、断层滑动等均会对煤岩体造成动载扰动，如何科学评估震动波扰动下的冲击危险程度是预测预警方面亟需开展的工作。

震动波扰动会使煤岩体内部裂隙发育扩展，造成煤岩体损伤劣化，整体强度降低；另一方面震动波携带的能量会在煤岩体内部产生集聚，使得煤岩体内部聚集大量弹性能，两者均会导致冲击危险性迅速升高。震动波对煤岩体的损伤研究方面，周朝等提出了一种微震损伤效应的围岩稳定性分析方法，可有效反应施工扰动和地质构造活化的微震信息。刘书贤等量化分析了采动区域开采变形和地震对建筑物的损害程度。李铁等提出了一种基于采动顶、底板岩层损伤的冲击地压预测方法。震动波的扰动研究方面，曹安业等研究得出矿震动载作用于煤岩体，造成煤岩体能量标量叠加和应力矢量叠加，使得煤岩系统“差能”增加，当满足冲击地压发生条件时可诱发冲击。王桂峰等构建了微震累积能量密度指标并用于冲击危险性预测。上述研究主要针对震动对煤岩体的损伤或者扰动效应单方面进行研究，未考虑两者的叠加效应，并且未考虑震动波能量衰减效应造成的影响。因而亟待提出一种基于震动波能量衰减特性，并综合考虑震动对煤岩体的损伤以及动载扰动的冲击危险预测方法及指标。

笔者以甘肃某矿冲击地压频发工作面为背景，统计分析了工作面回采期间冲击明显显现位置和冲击震源的分布特征；在统计了工作面不同震源能级的能量衰减系数的基础上，理论分析了影响震动波能量衰减的因素，特别讨论了震源能量、震源半径、受载点距震源的距离以及能量衰减系数对能量衰减的影响。在考虑震动波能量衰减的基础上，分别构建了基于震动波损伤作用的静载强度指标和动载强度指标，将两者结合并赋予不同权重，提出一种基于震动波能量衰减特性的危险性预测方法及指标，并进行了方法应用及验证。

1 工程实例背景

1.1 250106-1工作面概况

甘肃某矿主采煤层为5煤层，煤层厚度为17.58～48.01 m，平均厚度37.51 m，煤层具有强冲击倾向性，顶、底板分别具有弱冲击倾向性。采用分层综采放顶煤采煤法开采，共分上、中、下3个分层。

250106-1工作面为5煤层上分层工作面，平均开采厚度为12.9 m，煤层平均倾角7°。工作面布置如图1所示，综合柱状图如图2所示。

图1 250106-1工作面布置Fig.1 Longwall face layout of Longwall 50106-1

图2 250106-1工作面综合柱状Fig.2 Synthesis column around Longwall 250106-1

该矿已安装SOS微震系统用于冲击危险的监测与预警，250106-1工作面回采期间(2016年7月—2018年7月)监测到的大量微震数据可为后续分析提供基础。

1.2 冲击显现概况

250106-1工作面自2016年7月—2017年11月期间发生24起冲击事件，造成工作面两回采巷道不同程度损坏，冲击震源定位如图3所示。

图3 250106-1工作面回采期间冲击震源定位Fig.3 Location of coal burst sources duringretreat of Longwall 250106-1

如图4所示，250106-1工作面冲击显现形式以底臌为主，底臌量可达0.1～1.0 m；部分区域伴随顶板冒落或两帮移近、锚杆失效等现象，顶板下沉量0.15～0.20 m，两帮移近量0.2～0.7 m。工作面开采初期破坏程度较为严重，严重时可造成巷道整体发生收缩。

图4 工作面典型冲击显现特征素描Fig.4 Sketch of typical impact characteristicsof Longwall 250106-1

统计了冲击明显显现位置距工作面的距离如图5所示。由图5可知，24起冲击事件中有20起发生于临空侧(回风巷)，此外同样有20起事件发生于工作面超前300 m范围内，说明工作面超前300 m内为冲击发生的高风险区域，并且临空侧回采巷道的冲击地压发生概率显著高于实体煤侧。

图5 24起冲击事件明显显现范围统计Fig.5 Statistics of the apparent occurrence range of24 coal burst incidents

由图6(a)可知，冲击震源能量介于10～10J，其中能量10～10J的事件有22起(91.22%)，能量10～10J的事件有2起，因此当震源能量大于10J时便存在诱发冲击的可能。由图6(b)可知，冲击震源中23个震源距明显显现区的距离小于300 m，说明当震源距离巷道小于300 m时诱发冲击的概率迅速升高。

图6 冲击震源的能量及距离显现区域的最短距离统计Fig.6 Energy of coal burst sources and the shortestdistance statistics from the occurrence area

不难发现，巷道冲击的可能性与震源能量和距震源的最短距离密切相关，当震源能量越大，距离巷道距离越短时诱发冲击的可能性也越大。

2 震动波能量衰减规律

2.1 震动波能量衰减的理论模型

采掘过程中顶板的破断失稳、断层滑移等会产生矿震并向外辐射震动波，其震源破裂尺度往往较大，例如震源能量为10J的矿震破裂尺度可达上百米，震源的破裂尺度范围内、外煤岩体变形差异巨大，甚至有学者提出震源破裂尺度范围内煤岩体产生永久非弹性变形，而破裂尺度范围外仅为弹性变形，意味着例如图5中统计的冲击事件，冲击震源对不同距离的巷道或者煤岩体的作用是显著不同的。

此外，由于震动波的几何和固有衰减，震源能量由点传播到点的残余能量可表达为

=e-

(1)

式中，为受载点的残余能量；为微震事件在点的能量；为点和点间的空间距离；为能量衰减系数。

为体现震源破裂尺度内外震动效应的差异，将震源假设为半径为的震源球，假设震动波能量在球体内不变，能量从球面开始衰减，则震源能量的衰减过程可表述为

(2)

其中，为震源的视半径；为震源的视体积；为震源的地震矩；为震源区域的剪切模量。估算震源半径的典型模型有Brune模型、Madariaga模型、拉伸破裂模型以及视半径模型等，但是前3者多基于一定的震源破裂机制的假设，在煤矿中适用性较差，而视半径为震源发震时期内煤岩体塑性变形的标量估计，与描述震源强度的地震矩和震源能量相关，与震源破裂机制、震源模型无关，可对震源破裂半径(扰动范围)进行较好的估计，因此本文选择视半径估算震源半径。从式(2)中不难发现，影响震动波能量衰减的因素主要有能量衰减系数、震源能量、半径以及受载点距震源的距离。

2.2 不同能级震源的能量衰减系数特征

在波前几何扩散和介质的阻尼作用下，震动波传播过程中质点的质点峰值速度()与能量均会发生衰减，可描述为

(3)

式中，p为受载点的质点峰值速度；p为震源处的质点峰值速度；为振幅衰减系数；质点处的震动速度和动能满足关系：

(4)

式中，为质点的能量；为质点介质的单位质量；为质点的质点峰值速度。

联立式(2)～(4)，可得

(5)

该矿SOS微震系统使用速度型微震探头，仅可捕捉到探头布置点的速度波形，通过拟合距震源不同距离探头的可拟合出振幅衰减系数，以250106-1工作面回采期间监测到的一起微震事件为例进行说明该拟合过程(图7)。

图7 vp传播衰减特征拟合曲线Fig.7 Fitting curve of vp propagation attenuation characteristic

如图7为衰减拟合曲线，由图7可知该次事件的衰减系数=0.000 471 3，则能量衰减系数=2=0.000 942 6。

图8为250106-1工作面2016年7月—2017年11月共11 578次微震事件的衰减系数拟合结果。由图8可知能级为10～10J时衰减系数最大，为6.983 9×10；能级为10～10J时衰减系数最小，为2.785 7×10；随震源能级不断增大，衰减系数持续减小，但是减小的幅度越来越小，并趋于平缓。

图8 250106-1工作面不同能级的vp衰减系数拟合结果Fig.8 Fitting results of vp attenuation coefficient fordifferent energy levels of Longwall 250106-1

根据能量衰减系数和衰减系数的换算关系，可得不同能级微震事件的能量衰减系数见表1。

表1 250106-1工作面不同能级的能量衰减系数

2.3 震动波能量对不同因素的衰减响应规律

从冲击地压与煤岩体震动的关系来看，诱发冲击地压的最低能量为1×10J。但是如何准确理解诱发冲击地压的最低能量却是一个问题，对一个震源能量为1×10J的微震事件来说，其震源能量为1×10J，但是其传播过程中的衰减导致震动波残余能量(即真正作用于煤岩体的能量)急剧减少，如图9(a)中④号曲线所示，其震源能量为10J，当传播100 m时，残余能量约为100 J；当传播400 m时，残余能量约为10 J，此时仅凭10 J的震动能量是难以诱发冲击的，对图9(a)中其余情况的分析也可得到类似的结论；同时从图9(a)中不难发现震动波传播相同距离时，震源能量越大，则残余能量也越大。因此，涉及从能量角度解释冲击地压的相关问题时，不应单考虑震源能量，应从震动波传播的能量衰减效应入手，重点考量衰减后的残余能量对煤岩体的影响。

图9 不同因素对震动波能量衰减的影响Fig.9 Effect of different factors on the energy attenuation of shock wave

如图9(b)所示，当震源能量、能量衰减系数相同，震源半径不同时，其能量衰减差异仍较大。当震源半径为20，40 m时，以残余能量1×10J为界限，其具有冲击地压风险的最大距离分别为200，400 m左右；当震源半径为60，80或100 m时，传播距离为400 m时，残余能量仍远大于1×10J。说明震源半径越大，震动波能量衰减也越慢，易诱发煤岩体发生冲击的范围也越广。

如图9(c)所示，当震源能量、震源半径相同，能量衰减系数不同时，其能量衰减差异较小。结合表1中不同能级的能量衰减系数统计结果可发现，不同能级的微震事件传播过程中衰减系数差异较小，究其原因，能量衰减系数主要与震动波的固有衰减有关，即主要与震动波频率、波速与耗散品质因子等相关，而煤矿井下不同能级的微震事件传播环境类似，因此衰减系数的差异较小，进而导致其对能量衰减的影响较弱。

综上得出：① 微震事件能否诱发冲击与震源能量、震源半径以及距离震源的距离密切相关，而不同能级微震事件的衰减系数差异对其影响较小。② 震动波的残余能量随震源能量、震源半径的增加而增加，诱发冲击的有效范围相应增加，冲击危险性升高。③ 受震动波固有性质以及井下传播环境的影响，不同能级微震事件的能量衰减系数差异不明显，但利用能量衰减效应进行冲击地压的监测预警时，为进一步提高监测预警的准确率，则必须考虑该差异性。

3 基于震动波能量衰减特性的冲击危险预测指标构建

3.1 基于震动波损伤效应的静载强度指标

震动波作用导致煤岩体原有缺陷进一步劣化，造成煤岩体强度降低。根据连续介质损伤力学，损伤条件下的煤岩体全应力应变曲线如图10所示，表达式为

=(1-)

(6)

式中，为煤岩体所受应力；为弹性模量；为应变；为损伤参量。

图10 煤岩体受载过程中全应力应变及损伤参量变化Fig.10 Variation diagram of total stress-strain and damageparameters of coal and rock mass under loading

矿井尺度下运用式(6)描述煤岩体的应力状态时，如何定量表述煤岩体的应变，以及损伤参量的构建是关键。BENIOFF研究发现地震的弹性应变增量与其能量平方根(即Benioff应变)成正比，并且可线性叠加，因此震动波作用下煤岩体产生的累积应变可采用能量的平方根之和来近似估计。

如图11所示，将研究区域划分成若干网格，假设在研究时间段内共发生次微震事件，可根据式(2)计算出单个震源对所有网格节点(即受载点)的残余能量，依次对每个震源进行计算可得出受载点在次微震事件累积作用下的Benioff应变。

图11 任一受载点i的累积Benioff应变计算示意Fig.11 Calculation of cumulative Benioff strain for any load point i

据此，可通过Benioff应变来估算煤岩体产生的应变以及构建损伤参量：

(7)

式中，e为受载点的累积Benioff应变；为比例系数；为微震次数；临界应变的计算公式为

(8)

式中，为研究区域内e的最大值；为临界损伤值，对应于煤岩体完全损伤状态，=0.95。

联立式(6)～(8)，可得基于震动波损伤效应的静载应力表达式：

(9)

由于研究区域内不同位置的应力水平可能相差较大，即不同节点的应力数值相差较大，绘制云图时降低可视化效果，因而为增强计算结果的可视化程度，对s进行归一化处理，得到基于震动波损伤效应的静载强度指标：

(10)

式中，和分别为研究区域内s的最大和最小值。

3.2 描述弹性能集聚程度的动载强度指标

煤岩体应力调整和破裂伴随弹性能的聚散并向外辐射震动波，可被微震系统捕捉为微震事件。据此，前人利用微震事件的异常集聚来判断冲击危险性，但缺点是该类图像便于主观判断，无法定量分析。因而本文通过构建能量聚集指数来定量描述煤岩体的弹性能集聚程度：

(11)

其中，d为研究区域内某一受载点的能量聚集指数;为研究区域一定时间内的微震事件总数；为根据式(2)计算所得的受载点的残余能量。同样为增强可视化效果，对d归一化后得到受载点的动载强度指标d为

(12)

式中，和分别为研究区域内d的最大和最小值。

3.3 综合预测指标的构建及计算流程

根据冲击地压发生的动静载叠加原理，冲击地压发生时满足条件：

+≥

(13)

式中，和分别为静载和动载应力；为冲击地压发生临界应力。

据此，应从静载和动载2方面对冲击危险进行评价，因而本文采用静载强度指标衡量煤岩体的静载集中程度，采用动载强度指标衡量煤岩体的动载扰动强度，进一步建立动静载危险指标叠加的综合预测指标，受载点的综合预测指标sd可表示为

(14)

其中，，分别为静载和动载强度指标的权重因子，取值范围为0～1。为使预测的结果更贴合实际，此处采用标准离差法确定权重，计算公式为

(15)

其中，()，()分别为研究区域内s和d的标准差。表2为，与冲击危险性的对应关系。本方法的计算流程如图12所示。

表2 φs，φd,φsd与冲击危险性的对应关系

图12 基于震动波能量衰减特性的冲击危险预测方法计算流程Fig.12 Calculation flow chart of coal burst prediction methodbased on energy attenuation characteristics of shock waves

4 方法验证和实践

受篇幅所限，选取2016-08-30发生的一起较为严重的冲击事件(巷道破坏长度可达80 m，后文记为“8·30”冲击事件，“8·30”冲击明显显现区域记为显现区域)进行方法验证。如图13所示，冲击发生前微震日频次和日总能量均有明显上升趋势，并且频繁发生大能量事件，表明工作面冲击危险性急剧升高。

图13 “8·30”冲击事件发生前微震变化情况Fig.13 Changes of seismicity before the“8·30” coal burst

以10 d为时间窗，统计了“8·30”冲击事件发生前250106-1工作面的微震分布情况。如图14(a)所示，受开切眼外错影响，工作面外错区域采空区覆岩活动剧烈，大能量事件(震源能量≥10J)频发。同时受250104-1采空区侧向支承压力和本工作面超前支承压力的影响，小能量事件(震源能量<10J)多向工作面前方以及临空侧回采巷道集聚，表明临空侧回采巷道冲击危险性较高。

图14 250106-1工作面微震分布Fig.14 Microseismic distribution of Longwall 250106-1

如图14(b)所示，当工作面继续向前推进，受开切眼外错以及采空区悬顶范围增加影响，采空区顶板活动更为剧烈，大能量事件有向工作面中部以及临空侧回采巷道转移的趋势，同时小能量事件仍在临空侧回采巷道集聚较为明显，表明工作面临空侧回采巷道冲击危险性进一步升高。

图14(c)为“8·30”冲击发生前微震定位图，由图14(c)可知冲击发生前微震活动范围缩小，并明显向工作面周围集聚，尤其在显现区域附近大能量事件集中明显，同时震源能量0～10J的微震事件也向该区域集聚，说明该区域冲击危险性较高，但是存在多个微震异常聚集区域，凭借微震聚集情况仅能主观判断冲击危险区域，无法准确预测冲击危险区域及危险程度。

为评价，以及指标预测大能量事件与显现区域的效果，后文在计算了，以及分布云图的基础上，采用未来10 d内的大能量事件位置以及显现区域与预测的冲击危险区域的重合情况进行验证。

需要注意的是，在计算，以及时，需选用合理的时间窗以便有足够的数据保证计算效果。“8·30”冲击发生前250106-1工作面微震日频次平均为25个(图13)，经笔者多次尝试并充分考虑，的物理含义，分别采用10，3 d为时间窗可满足对，的反演效果。例如图15(a)和图16(a)中“8月11日”的云图，分别是以8月1日—8月10日和8月8日—8月10日期间的微震数据计算所得；其中大能量事件定位为8月11日—8月20日的数据定位所得。

图15 250106-1工作面φs分布云图与微震定位Fig.15 φs and seismic events distribution of LW250106-1

图15为“8·30”冲击发生前一个月内的分布云图与微震定位。由图15可知随工作面推进位置不同，工作面静载集中水平变化较大，其中工作面初采阶段(8月11日)应力集中范围主要位于采空区后方以及实体煤回采巷道一侧；当工作面继续向前推进，工作面静载应力有向临空侧回采巷道转移的趋势(图15(b))，尤其后续显现区域附近也出现了一定程度的应力集中；当工作面继续向前推进直至发生冲击前一天，工作面应力集中范围进一步发生转移，应力集中程度最高的区域位于临空侧工作面端头附近；同时显现区域附近仍处于应力集中区域内，具有中等冲击危险性。

从图15中大能量事件、显现区域与预测的中等以上冲击危险区域的重合程度来看：① 大能量事件与中等以上冲击危险区重合度较差，说明对大能量事件的预测效果较差；② 显现区域与中等冲击危险区重合度较高，说明对显现区域预测效果较好。

图16为“8·30”冲击发生前一个月内的分布云图与微震定位。由图16发现：①仅可对部分大能量事件进行预测，说明对大能量事件的预测效率较低；② 工作面前方尤其是临空侧超前区域较大范围预测为强冲击危险，侧面说明显现区域与强冲击以上危险区域重合程度较差，说明对显现区域的预测效果较差。

图16 250106-1工作面φd分布云图与微震定位Fig.16 φd and seismic events distribution of Longwall 250106-1

图17 250106-1工作面φsd分布云图与微震定位Fig.17 φsd and seismic events distribution of Longwall 250106-1

图17为对，叠加后的分布云图与微震定位情况。由图17(a)可知，工作面前方约300 m范围以及采空区临空侧>0.5，具有中等冲击危险性，工作面前方中部局部区域大于0.75，具有强冲击危险性。8月11日—8月20日共发生17起大能量事件，其中有12起位于中等冲击危险区域(>0.5)内，对大能量事件预测准确率为71%。

由图17(b)可知，此时工作面前方的分布与8月11日相比发生较大变化。较高的区域由工作面中部向临空侧区域转移，其中显现区域附近进一步增大，局部超过0.75，具有强冲击危险性。此外，8月21日—8月29日期间19起大能量事件中17起位于中等冲击危险区域(>0.5)内，对大能量事件预测准确率为89%。

由图17(c)可知“8·30”冲击震源和显现位置均位于中等冲击危险区域，其分别高达0.64和0.67，说明通过有效预测了该次冲击震源和显现区域。

综上可知：①可对冲击明显显现区域进行预测，而对大能量事件的预测效果较差；②仅可对部分大能量事件进行预测，但预测的危险区域较广，侧面说明对冲击明显显现区域预测效能较低；③可有效预测大能量事件和冲击明显显现区域，预测效果较好。

5 结 论

(1)250106-1工作面24起冲击事件统计结果表明工作面超前300 m为冲击发生的高风险区域，临空侧回采巷道的冲击显现概率显著高于实体煤侧。

(2)煤矿井下震源破裂尺度对微震能量的衰减影响较大，分析得出影响能量衰减的主要因素有震源的能量、半径、衰减系数以及距受载点的距离；震源的能量、半径以及距受载点的距离对能量衰减的影响较为明显，而衰减系数的影响相对较小。

(3)充分考虑震动波能量的衰减效应，通过基于震动波损伤作用的静载强度指标来衡量静载应力集中程度，采用描述弹性能集聚程度的指标来衡量动载扰动强度；综合两者构建了基于震动波能量衰减特性的冲击危险预测指标，并采用标准离差法量化指标权重。

(4)预测实例表明对冲击显现区域的预测效果较好，而对大能量事件的预测效能较低；仅可对部分大能量事件进行预测，同时预测的冲击危险区域较广，侧面说明对冲击显现区域的预测效能较低；对于综合预测指标，大能量事件位于中等冲击危险区域的概率超过70%，并且冲击明显显现区域全部位于中等冲击危险区内，说明可有效预测大能量事件位置和冲击显现区域，预测效果较好。