Prediction method of coal burst risk based on focal mechanism and location error calibration
-
摘要:
冲击地压智能监测技术迅速发展大背景下，如何利用海量微震监测信息挖掘震源破裂性状并将其应用于冲击地压的预测预警或防控是未来研究的重点与难点。为解决目前微震群破裂机制聚合规律挖掘不充分、震源机制反演技术难以直接应用于冲击地压的预测预警难题，采用理论分析、数值仿真等手段研究建立基于震源机制与定位误差校准的冲击地压危险预测方法。以内蒙古某冲击地压矿井大能量矿震事件频发2215工作面为研究背景，总结了煤矿微震波形数量有限、噪声影响大、反演效率要求高的特点，比较分析了不同矩张量方法反演震源机制的特点及优缺点，确立了混合矩张量反演法更适用于煤矿应用场景；2215工作面顶板破断诱发事件以拉剪以及压剪破坏为主；大能量事件发生前微震事件剪切成分占比不超过40%，顶板破断主导的微震事件走向角主要分布于3°～80°以及150°～270°，倾角主要分布于70°～90°，存在明显的集聚特征；基于煤矿常见的震源破裂机制及定位误差分析了微震事件聚合的可能性并建立了微震聚合准则，构建了基于震源破裂机制及定位误差校准的危险预测方法及指标；预测实例结果表明考虑震源破裂机制及定位误差后可显著增强微震异常聚合区的识别；微震数据较为完整情况下，基于震源破裂机制及定位误差校准的危险预测指标IP阈值设置为0.7时可实现对不同破裂类型大能量事件的准确预测。
Abstract:Under the background of the rapid development of intelligent coal burst monitoring technology, how to use massive seismic monitoring information to investigate source rupture and apply it to the prediction, early warning, or prevention of coal burst is the focus and difficulty of future research. To solve the problems of insufficient investigation of seismic clustering fracture mechanism aggregation law and overcome the difficulty of direct application of seismic source mechanism inversion technology to the prediction and early warning of coal burst, the theoretical analysis and numerical simulation methods are used to establish a prediction method of coal burst risk based on seismic source mechanism and location error calibration. Taking the 2215 working face with frequent high magnitude seismic events in a coal burst mine in Inner Mongolia as an example, this paper summarizes the characteristics of the limited number of seismic waves, large noise impact, and high inversion efficiency. Also, it compares and analyzes the characteristics, advantages, and disadvantages of different moment tensor inversion methods, and establishes that the hybrid moment tensor inversion method is more suitable for coal mine application scenarios. The statistical results of the rupture mechanism of high magnitude events in the 2215 working face show that the main inducing events of roof rupture are a tensile shear failure and compressive shear failure. The shear component statistics of seismic events before the occurrence of high magnitude events show that the proportion is not more than 40%. The occurrence statistics of source rupture surface show that the strikes of seismic events dominated by roof breaking are mainly distributed in the range of 3°−80° and 150°−270°, and the dips are mainly distributed in the range of 70°−90°, with prominent agglomeration characteristics. Based on the common source rupture mechanism and location errors in the coal mine, the possibility of seismic event aggregation is analyzed, and the criteria of seismic aggregation are established. The risk prediction method and index based on the source rupture mechanism and location error calibration are constructed. The results of prediction examples show that the identification of seismic anomaly aggregation areas can be significantly enhanced by considering the source rupture mechanism and location error. When the seismic data is relatively complete, the IP threshold of the risk prediction index based on the source rupture mechanism and location error calibration is set to 0.7, and the accurate prediction of high magnitude events of different rupture types can be achieved.
-
冲击地压作为煤矿中强烈的动力灾害现象，致灾时间短、破坏范围大，并且往往伴随强矿震的产生[1]。就冲击地压灾害的防控而言，在国家有关部门和科研院所加大管控和防治研究力度的背景下已得到有效遏制。然而近年来，随着我国煤炭开采深度整体向着深部转移，在我国新开发的陕西、内蒙古等特大矿区也涌现出一批冲击地压矿井或受冲击地压威胁的矿井。例如仅2020年底以来，榆林金鸡滩煤矿、内蒙古石拉乌素煤矿、红庆河煤矿等陆续发生大能量矿震(最大震级3.1级)[2]。由于冲击地压与大能量矿震的相互作用关系尚不明晰，而冲击地压发生时必然伴随大能量矿震事件[1,3]。因此，如何做到对大能量事件的有效监测预警，既是大能量矿震防控的关键，也是有效防治冲击地压的前提与保障。
随着微震监测系统在冲击地压矿井的广泛部署，基于微震监测的震源破裂机制解析技术已广泛应用于冲击地压机制及预警研究中。陈栋等[4]基于Brune模型和网格搜索法求解了千秋煤矿震源机制及震源参数，并尝试评估和量化震源破裂的尺度和机制。曹安业等[2]采用相对矩张量算法研究了褶皱构造区、深埋断层构造区以及高应力大巷煤柱区等3个冲击地压频发区矿震破裂机制与震源参量响应规律。李铁等[5]基于震源机制解分析了矿井尺度采动应力场的演化规律并研究了深部开采断层活化异变特征。LISLE等[6]提出了基于震源机制解中的P-T轴反演震源区域主应力方向的方法。FUJII等[7]分析了冲击地压发生前最大剪切矩张量释放率的变化特性，建立了基于矩张量的冲击地压预测方法。MA等[8]利用震源机制反演方法研究了顶板关键层的采后破断运移状态和震源破裂面产状分布特征。SONG等[9]研究了不规则煤柱区不同区域的大能量事件破裂机制与破裂参数分布特征。WU等[10]研究了厚硬顶板破断滑移诱发大能量事件的破裂机制及破裂面空间分布特征。然而，已有研究多针对大能量事件破裂参数的统计分析或大能量事件发生前某一震源参量的变化特征分析，未对大能量事件与其发生前微震事件的破裂面参数的相关性及聚集特征作深入分析。
根据冲击地压监测预警目的的差异，可分为冲击地压时序预警以及危险区域空间预测[11]。就冲击危险的空间预测而言，为解决冲击地压海量监测数据利用不足的困境，聚类分析算法已被诸多学者引入冲击前兆规律分析及危险源的探测之中。MEYER等[12]结合微震事件的定位以及震源机制描述断层可能的核位置，并采用期望最大化聚类算法识别断层活动段的概率核。LURKA[13]提出微震事件时-空分层聚类分析方法，并结合Gutenberg and Richter定律进行大能量矿震的风险评估。FALMAGNE等[14]考虑微震事件的空间距离以及震源半径并提出聚类指数函数量化事件间的聚合可能性。HUDYMA等[14]采用综合聚类方法以识别不同破裂机制的煤岩体破坏并进行冲击风险评估。COULSON[15]基于主成分方法分析了微震事件距离分布特征并用以识别断层构造。OUILLON等[16]基于震源破裂面的特征向量和特征值识别断层段的方向和大小并实现断层产状的可视化。CEACA等[17]提出了断层方向识别的震源机制聚类方法。已有研究在分析时仅从微震事件的发生时间、空间距离或震源机制方面进行分析，未充分考虑3者的内在联系。因而亟待提出一种综合微震时间、空间以及震源机制等信息的冲击地压危险分析预测方法。
笔者以内蒙古某冲击地压矿井大能量矿震事件频发工作面为背景。总结了煤矿应用场景下震源机制反演的特点，比较分析了不同矩张量方法反演震源机制的特点及优缺点，确立了适用于煤矿应用特点的震源机制反演方法；求解了顶板破断诱发大能量事件的震源破裂机制，统计分析了大能量事件发生前微震破裂机制及破裂面参数的集聚规律；分析了考虑震源破裂机制及定位误差条件下微震事件聚合的可能性并建立了聚合准则，构建了基于震源破裂机制及定位误差校准的危险预测方法及指标，并进行了方法应用及验证。
1. 工程背景
内蒙古某冲击地压矿井主采2-2煤层，为近水平煤层开采，煤层及其顶底板均具有弱冲击倾向性。直接顶以砂质泥岩为主，基本顶以中、细砂岩为主，巨厚层状，水平层理，成分以石英、长石为主。
本文研究对象2215工作面位于矿井22采区东南部，北为2213设计工作面，南为2217采空区，西为井底广场及大巷区域，东为实体煤区域，如图1所示。2215工作面及主要采掘区域附近已覆盖微震监测台网，可实时采集微震波形并进行定位及震源能量计算。
![]()
图 1 2215工作面及微震监测台网布置Figure 1. Layout of 2215 working face and seismic monitoring network2215工作面煤层平均采深约为731 m，煤厚5.64～7.57 m，平均厚度6.41 m，平均倾角为2°。由工作面综合柱状图(图2)可知，煤层上方100 m范围内赋存有多层厚硬砂岩层，工作面回采过程中厚硬砂岩层的悬顶可造成工作面前方煤体支承压力急剧升高，是无明显构造条件下冲击地压发生的主要静载应力源，并且厚硬砂岩层的破断可诱发大能量事件，对工作面安全回采造成极大威胁。
如图3所示，矿方在进行顶板预裂爆破等多种卸压措施的情况下，2215工作面2022−01−09—2022−05−17回采期间仍监测到多达17起大能量事件(震源能量 > 105 J)。大能量事件的分布区域可明显划分为3个区域，区域Ⅰ的事件偏向工作面侧向煤柱侧，区域Ⅱ的事件趋向工作面中部，而区域Ⅲ的事件趋向采空区煤柱侧。根据已有研究[18]，位于工作面中部的大能量事件多由本工作面顶板的“O-X”破断诱发，而煤柱侧的事件多由相邻采空区与本工作面顶板协同“O-X”破断造成，采空区后方的事件则由顶板“砌体梁”结构的滑移失稳导致。可见工作面不同区域大能量事件的内在机制存在差异，但不同地质与采掘条件下可能存在特殊性，亟需开展相关研究。
2. 震源机制反演方法
矩张量反演作为震源机制反演方法之一，近年来发展迅速，在计算方法和可靠性上都大大增强，但是矩张量反演法存在众多类型[19-20]，根据反演对象和目的的不同，应对反演方法进行比较和选择。
2.1 震源机制反演方法比较及选择
矩张量模型采用二阶矩张量表示震源所受的力，由9个大小和方向不同的力偶(矩张量分量)组成[2]。当震源满足点源假设和同步震源假设时，震源破裂激发的位移场、格林函数以及矩张量的关系[21]可表示为
$$ {{\boldsymbol{u}}_k} = {{\boldsymbol{M}}_{ij}} \frac{{\partial {G_{ki}}}}{{\partial {X_j}}} = {{\boldsymbol{M}}_{ij}} {G_{ki,j}}\begin{array}{*{20}{c}} {}&{(i,j,k = 1,2,3)} \end{array} $$ (1) 式中，${{\boldsymbol{u}}_k}$为震源激发的位移场矩阵；${{\boldsymbol{M}}_{ij}}$为矩张量；${G_{ki}}$为格林函数；${G_{ki,j}} $为格林函数在k方向的分量；${X_j}$为笛卡尔坐标系下某一坐标轴；i，j分别为笛卡尔坐标系下某一坐标轴方向；k为震动波震相(k=1，2，3，分别代表P波、SH波以及SV波)。
对输入位移场以及格林函数处理方法的不同，矩张量反演法可按2种原则进行分类。按照反演所用位移场不同可分为振幅法、振幅比法以及全波形法[19]；按照对格林函数处理方式不同可分为绝对、相对以及混合矩张量法[20]。不同方法特点及优缺点比较见表1。
表 1 不同震源机制反演方法及优缺点比较Table 1. Comparison of different source mechanism inversion methods and their advantages and disadvantages分类原则 反演方法 特点 优点 缺点 备注 按对格林函数处理方式不同 绝对法 根据理论或经验模型直接构建格林函数，以单个震源进行反演 无需基于额外假设、震源可单独求解 构建精确格林函数相对困难、误差不可控、对波形数量要求较高 绝对矩张量法、相对矩张量法、混合矩张量法均属于振幅法 相对法 基于震动波射线传播理论对格林函数进行转化与分解，以震源群为单位进行反演 避免直接求解格林函数、对波形数量要求相对较低 无法单独求解震源机制、误差不可控、要求参与反演的震源群具有多样的震源机制 混合法 融合绝对法的格林函数构建以及相对法的震源群思想，迭代修正绝对矩张量解 求解精确度高，误差可控，反演迭代速度快 构建精确格林函数相对困难 按所选用位移场的不同 振幅法 使用P波初动振幅作为位移场信息求解矩张量 仅使用P波初动振幅作为约束
条件，反演方程较为简单对波形信噪比要求较高，反演精度受P波初动拾取准度影响大 振幅比法 使用P波、S波初动极性及振幅比作为位移场信息求解矩张量 使用P波、S波极性及振幅比作为约束条件，一定程度上提高反演精度 对P波、S波到时及极性拾取准度要求高 全波形法 通过特定频率段实测波形与理论波形的拟合迭代确定最优的
震源机制解求解精确度高，误差可控 对波形信噪比要求高、计算量大反演周期长、精确格林函数构建困难 由表1可知，不同反演方法优缺点不同，反演的精确度与计算量上也有差异。目前各种反演方法在地震、水压致裂诱发地震研究中均有广泛应用，但其是否适用于煤矿仍需进一步研究。煤矿井下微震台网随采掘位置而进行挪移、导致不同时期微震事件波形数量有限，同时微震波形易受噪声污染、信噪比较低，煤矿微震数据量大也是其典型特点，因此煤矿震源机制反演方法的选择具有其特殊性。
煤矿震源机制反演方法的选择，须满足以下3点：① 克服单个事件波形数量有限，可支持微震事件大规模反演；② 尽量降低噪声、格林函数构建误差的影响，反演准确度高；③ 反演效率高，可快速处理大量微震数据。对比表1中所列几类反演方法，绝对法误差不可控，对波形质量要求较高；相对法以震源群为单位进行反演，反演效率较高，但误差不可控；振幅比法和全波形法虽反演精度较高，但是均对波形质量要求较高，同时反演周期较长；而混合矩张量法虽在构建格林函数方面存在困难，但是可通过迭代对波形噪声、不准确的格林函数造成的误差进行修正，反演迭代速度较快，可满足煤矿震源机制反演的需求。
混合矩张量反演流程如图4所示，主要分为2步。首先是矩张量初始解的获取，获取微震事件矩张量初始解后，将所有震源筛选为不同的震源群，每个震源群中所有震源联合进行迭代更新，当震源群中的理论振幅与预测振幅间的误差小于阈值时停止迭代，结束反演。其中，振幅更新服从式(2)[22]，震源群中某一震源的理论与预测振幅的误差计算公式为式(3)[22]。
$$ u_{pq}^{{\rm{new}}} = u_{pq}^{{\rm{old}}} + {w_p}u_{pq}^{{\rm{old}}}({\tilde r_p} - 1) $$ (2) $$ {\rm{RMS}} = \sqrt {{{\displaystyle\sum_{q = 1}^N {{{\left( {u_q^{{\rm{new}}} - u_q^{{\rm{old}}}} \right)}^2}} }}\left/{{\displaystyle\sum_{q = 1}^N {u_q^{{\rm{new}}}} }}\right.} $$ (3) 式中，$u_{pq}^{{\rm{new}}}$和$u_{pq}^{{\rm{old}}}$分别为本轮和上一轮迭代时第q个震源的第p个通道的波形振幅；${w_p}$为权重系数，取值为0～1；${\tilde r_p}$为震源群中通道p的振幅比，${\tilde r_p}{\text{ = }}u_{pq}^{{\rm{new}}}/u_{pq}^{{\rm{old}}}$；RMS为震源群中所有震源的位移误差的平方根；p和q分别为通道和震源编号。
2.2 矩张量分解及表示
获得矩张量M后，进一步对其本征值化可获得矩张量M的3个特征向量(${e_1}$，${e_2}$，${e_3}$)及其对应的特征值($ {M_1} $，$ {M_2} $，$ {M_3} $)，震源破裂的各向同性(ISO)、轴向拉伸(CLVD)以及剪切成分(DC)等3个成分与矩张量特征值的关系如式(4)所示，进而可计算震源破裂各向同性、轴向拉伸、剪切成分的占比($ {P_{{\rm{ISO}}}} $，$ {P_{{\rm{CLVD}}}} $，$ {P_{{\rm{DC}}}} $)如式(5)所示，3者满足关系$ {P_{{\rm{ISO}}}} + {P_{{\rm{CLVD}}}} +{P_{{\rm{DC}}}}{\text{ = }}1 $，$ {P_{{\rm{ISO}}}} $和$ {P_{{\rm{CLVD}}}} $的取值范围为−1～1，$ {P_{{\rm{DC}}}} $的取值范围为0～1。
$$ \left\{ {\begin{array}{*{20}{l}} {{M_{{\rm{ISO}}}} = \dfrac{1}{3}\left( {{M_1} + {M_2} + {M_3}} \right)} \\ {{M_{{\rm{CLVD}}}} = \dfrac{2}{3}\left( {{M_1} + {M_3} - 2{M_2}} \right)} \\ {{M_{{\rm{DC}}}} = \dfrac{1}{2}\left( {{M_1} - {M_3} - \left| {{M_1} + {M_3} - 2{M_2}} \right|} \right)} \end{array}} \right. $$ (4) $$ \left\{ {\begin{array}{*{20}{l}} {{P_{{\rm{ISO}}}} = \dfrac{{2\left( {{M_1} + {M_2} + {M_3}} \right)}}{{\left| {2\left( {{M_1} + {M_2} + {M_3}} \right) + \left| {4\left( {{M_1} + {M_3} - 2{M_2}} \right)} \right|} \right| + 3\left( {{M_1} + {M_3} - 2{M_2}} \right)}}} \\ {{P_{{\rm{DC}}}} = \dfrac{{3\left( {{M_1} - {M_3} - \left| {{M_1} + {M_3} - 2{M_2}} \right|} \right)}}{{\left| {2\left( {{M_1} + {M_2} + {M_3}} \right)} \right| + \left| {4\left( {{M_1} + {M_3} - 2{M_2}} \right) + 3\left( {{M_1} - {M_3} - \left| {{M_1} + {M_3} - 2{M_2}} \right|} \right)} \right|}}} \\ {{P_{{\rm{CLVD}}}} = \dfrac{{4\left( {{M_1} + {M_3} - 2{M_2}} \right)}}{{\left| {2\left( {{M_1} + {M_2} + {M_3}} \right)} \right| + \left| {4\left( {{M_1} + {M_3} - 2{M_2}} \right)} \right| + 3\left( {{M_1} - {M_3} - \left| {{M_1} + {M_3} - 2{M_2}} \right|} \right)}}} \end{array}} \right. $$ (5) 震源破裂面法向向量${\boldsymbol{n}}$与位错向量υ与矩张量特征值($ {M_1} $，$ {M_2} $，$ {M_3} $)满足关系式(6)，获得法向向量${\boldsymbol{n}}$与位错向量υ后可进一步计算出破裂面方位角、倾角与滑动角等产状信息以便后续分析。
$$ \left\{ \begin{gathered} {\boldsymbol{n}} = \sqrt {\frac{{{M_1} - {M_2}}}{{{M_1} - {M_3}}}} {{\boldsymbol{e}}_1} + \sqrt {\frac{{{M_1} - {M_2}}}{{{M_1} - {M_3}}}} {{\boldsymbol{e}}_3} \\ { {\textit{υ}}} = \sqrt {\frac{{{M_1} - {M_2}}}{{{M_1} - {M_3}}}} {{\boldsymbol{e}}_1} - \sqrt {\frac{{{M_1} - {M_2}}}{{{M_1} - {M_3}}}} {{\boldsymbol{e}}_3} \\ \end{gathered} \right. $$ (6) 式中，${\boldsymbol{e}}_1 $、${\boldsymbol{e}}_3 $分别为矩张量的最大、最小特征值对应的特征矢量。
目前常用的震源破裂类型统计分析方法有CLVD–ISO钻石图[23](图5)、Hodson倾斜钻石图[24]与Riedesel-Jordan月形图[25]。其中，CLVD–ISO钻石图表示直观、物理含义清楚被广泛使用。如图5所示，纯剪切破裂将落入图的中央，图的正上方与正下方分别对应外爆和内爆破裂类型，第一与第三象限边缘分别对应纯拉伸与压缩破裂。$ {P_{{\rm{DC}}}}\leqslant 0.4 $对应拉伸破裂(第一象限)或压缩破裂(第三象限)；$ 0.4 < {P_{{\rm{DC}}}}< 0.6 $对应拉伸−剪切混合破裂(第一象限)或压缩−剪切混合破裂(第三象限)；$ {P_{{\rm{DC}}}} \geqslant 0.6 $对应剪切破裂。本文后续将采用CLVD–ISO钻石图分析大能量事件的破裂类型统计规律。
![]()
图 5 震源破裂类型的CLVD−ISO钻石图Figure 5. Representation of source mechanism using CLVD−ISO diamond plot3. 大能量事件震源机制解统计分析
为研究2215工作面顶板破断诱发大能量事件的震源机制响应特征，采用混合矩张量反演法求解了17次大能量事件的震源机制解；并进一步求解了大能量事件发生前一周微震事件的震源机制解，分析不同破裂类型主导的大能量事件的震源机制解前兆集聚规律，为后续预测方法的提出提供依据。
3.1 震源破裂类型统计分析
图6为17次大能量事件震源机制沙滩球分布。由图6可知17次大能量事件中无纯剪切事件，说明顶板破裂运移过程较为复杂，不存在简单的剪切滑移破裂，均为拉张、压缩与剪切破坏相混合的模式；在混合破坏模式下，大能量事件的震动波辐射方向差异明显(从沙滩球的辐射方向来看)，说明顶板岩层受非均质性、原生裂隙及采动后主应力变化等因素的影响[26]，不存在明显的优势破裂方向。可见在无明显构造条件下，本工作面“O-X”破断、相邻采空区与本工作面顶板协同“O-X”破断、顶板砌体梁结构滑移失稳等顶板破断方式主导下震源破裂类型差异不明显，说明顶板破断方式与其震源破裂类型间不存在必然联系。
图7为17次大能量事件CLVD−ISO成分占比。由图7可知，17次大能量事件以拉剪和压剪破坏为主(13/17)，说明在工作面无明显构造条件下，顶板破断诱发的大能量事件以顶板的拉伸断裂和回转挤压为主[18]。此外，有2次剪切破裂事件(No.13和16事件)发生于临空侧，说明相邻采空区与本工作面顶板协同“O-X”破断方式主导下顶板更容易在在临近煤柱区域发生剪切滑移失稳。综合来看，临空侧顶板受采空区悬顶影响，破裂类型更为复杂，预测及防治的难度也越高。
17次大能量事件发生前一周微震事件剪切成分占比如图8所示。由图8可知17次大能量事件发生前微震事件剪切成分占比不超过40%，说明在无构造情况下，顶板破断产生的微震事件剪切成分占比较少，顶板破断运移主要以拉剪以及压剪破坏为主。
![]()
图 8 大能量事件发生前一周微震事件剪切成分占比Figure 8. Proportion of shear components of seismic events in the week before the occurrence of high magnitude events3.2 震源破裂面产状统计分析
图9为17次大能量事件发生前一周微震事件破裂面产状分布。由于矩张量的对称性，矩张量分解后对应的破裂面参数有2组，在无地应力测量以及现场观测的情况下很难对真实破裂面进行判别[27]，因此此处在统计时对2组破裂面参数均进行了统计。由图9可知破裂面1的走向角主要分布于3°～80°以及150°～270°，倾角主要分布于70°～90°；破裂面2的走向主要分布于215°～312°，倾角主要分布于30°～60°。由于多数微震事件由顶板破断引起，其产生的裂隙更可能趋向于垂直层面，即破裂面倾角更接近于90°，因此破裂面1的产状参数更接近于真实破裂面。这说明顶板破断诱发的微震事件相互之间并不孤立，究其原因，微震事件的产生伴随煤岩体采动裂隙的扩展和滑移，当多条裂隙扩展至一处使得裂隙发生贯通时，可诱发煤岩体瞬间失稳，煤岩体内部积聚的弹性能快速释放，进而诱发大能量事件。
![]()
图 9 大能量事件发生前一周微震事件破裂面产状分布注：左侧为破裂面1的走向分布，右侧为破裂面2的走向分布。Figure 9. Distribution of rupture surface of seismic events occurred one week before the high magnitude events综上，顶板破断时采动裂隙的扩展、滑移以及贯通等诱发一系列微震事件的产生，在震源机制的破裂类型与破裂面产状上表现出集聚性。从冲击地压危险预测的角度来看，如何量化大能量事件发生前微震事件破裂的集聚特征是实现大能量事件预测的关键。
4. 基于震源破裂机制及定位误差校准的危险预测方法
本节基于断裂力学理论从预测裂隙贯通的角度实现大能量事件的预测，分析了考虑定位误差条件下裂隙贯通的可能性并建立了相关准则，建立了定位误差校准的方法，并构建了基于震源机制与定位误差校准的冲击危险预测方法，从而实现大能量事件的准确预测。
4.1 震源破裂机制与聚类可能性分析
如图10所示，根据震源作用力方式和煤岩体破坏处与工作面的相对位置，震源破裂机制可分为6种(其中，红色圆点为微震震源)：顶板冒落型、煤柱失稳型、顶板破断型、顶板错动型、底板错动型以及顶板水平滑移型[28]。前3种震源破裂方式主要为应力主导，后3种主要为煤岩体中地质弱面的滑移错动造成。不论是应力主导型还是弱面错动型，煤岩体在大的破裂发生之前均会诱发一系列微震事件，这些微震事件往往在发生时间与空间位置上服从一定的规律，即大能量破裂事件发生前的前兆信息。
此外，研究表明微震系统的定位误差难以避免，如国内冲击地压矿井普遍装备的微震系统水平定位误差在20 m以内，垂直定位误差在50 m左右。不难发现，定位误差的存在导致采用震源定位信息分析的过程中存在较大的不确定性。例如聚类分析时，不同的标波结果使事件定位发生改变，进而导致同一震源可能属于不同聚类类别。另外，从裂隙贯通诱发煤岩失稳或者大能量事件的角度来说，聚类分析时仅考虑震源定位而忽略震源机制是不可取的。例如当震源破裂机制以内爆或外爆(即剪切成分占比较低)为主导时，当2个震源空间位置足够接近到一定距离时，直接将震源聚为一类是可行的；但是当2个震源机制为剪切破裂占主导时(即剪切成分占比较高)，此时震源存在明显的破裂面，破裂时裂隙往往沿破裂面扩展或错动，例如图10(d)～(f)对应的震源破裂情况。此时即便2个震源空间位置足够接近，但是其破裂面若是平行的话，裂隙也难以相互贯通，显然此时将2个震源聚为一类是不合理的。因此在聚类时需同时考虑定位误差及震源破裂机制。
已有研究表明，由定位误差导致的震源位置服从椭圆概率分布[29]。如图11所示，当考虑定位误差和震源破裂机制时，震源A和B存在3种相对关系。其中，位置A和B分别为震源的真实位置，A′和B′分别为由于定位误差导致的震源可能的位置，其分布在黑色椭圆范围中(误差椭圆)震源A和B的3种相对位置关系可描述为：① 震源A(或A′)与B(或B′)的距离d(或d′)大于裂隙可能贯通的最短距离D且震源A与震源B的破裂机制相同；② 震源A(或A′)与B(或B′)的距离d(或d′)不大于裂隙可能贯通的最短距离D且震源A与震源B的破裂机制相同；③ 震源A(或A′)与B(或B′)的距离不大于裂隙可能贯通的最短距离d(或d′)且震源A与震源B的破裂机制不同。
![]()
图 11 结合震源破裂机制与定位误差的聚类可能性注：红色和蓝色虚线圆分别代表震源A和B的震源半径。Figure 11. Clustering possibility combining source mechanism and source location error不难发现，当震源A和B满足关系①时，由于相距距离较远而产生的破裂机制难以相互影响；震源A和B满足关系②时，两震源虽在裂隙可相互贯通的范围之内，但是由于震源A和B的破裂机制相同，意味震源A和B产生的裂隙面相互平行，裂隙也难以相互贯通；只有当震源A和B满足关系③时，两震源才满足裂隙相互贯通的条件。因此可得出震源A和B破裂裂隙相互贯通的条件为
$$ \left\{ {\begin{array}{*{20}{l}} {d \leqslant D\;\;(D = 2\left( {{r_A} + {r_B}} \right))} \\ {{F_A} \ne {F_B}} \end{array}} \right. $$ (7) 式中，d为震源A和B之间的距离；D为震源A和B破裂时裂隙可能贯通的最短距离；rA和rB分别为震源A和B的震源半径；FA和FB分别为震源A和B的破裂机制。
4.2 危险预测指标构建及计算流程
微震事件诱发微裂隙扩展、汇聚进而增大诱发冲击的可能，因此结合微震定位误差与震源破裂机制的聚类分析的关键是量化微震事件聚集的程度，其物理含义为采掘活动区域某处受微震影响下微裂隙聚集的疏密程度。在同等应力条件下，显然微裂隙越发育，导致顶板完整性越差，发生整体性失稳的可能性也越大，对应冲击危险性越高。
根据4.1节提出的思路，聚类分析的前提是量化微震定位误差，即考虑定位误差条件下微震事件可能被定位的位置。为使研究结论更具有代表性，本节采用在煤矿微震定位算法中使用最广泛的P波初至定位算法进行定位误差数值仿真试验，其核心公式如式(8)～(10)所示。由于影响定位精度的因素主要有人工标波误差和速度模型精度，为量化2者对于定位误差的影响，参照前人研究[29]，笔者特引入人工标波误差因子$ {\delta _{\text{P}}} $和速度模型误差因子$ {\delta _{{v}}} $量化2者的影响(其中，$ {\delta _{\text{P}}} $和$ {\delta _{{v}}} $均为服从正态分布的随机数)，并在此基础上进行数值仿真，计算受标波和速度模型误差影响下微震事件的分布概率，数值仿真试验关键参数取值见表1。
$$ \varPhi ({t_0},{x_0},{y_0},{{\textit{z}}_0}) = \sum\limits_{\alpha = 1}^m {R_\alpha ^2} $$ (8) $$ {R_\alpha } = {t_\alpha } + {\delta _{{\text{P}}\alpha }} - {t_0} - {T_\alpha }({x_0},{y_0},{{\textit{z}}_0}) $$ (9) $$ {T_\alpha }\left( {{x_0},{y_0},{{\textit{z}}_0}} \right) = \frac{{\sqrt {{{\left( {{x_0} - {x_\alpha }} \right)}^2} + {{\left( {{y_0} - {y_\alpha }} \right)}^2} + {{\left( {{{\textit{z}}_0} - {{\textit{z}}_\alpha }} \right)}^2}} }}{{{v_{\text{P}}} + {\delta _{{{v}}\alpha }}}} $$ (10) 式中，$ \varPhi ({t_0},{x_0},{y_0},{z_0}) $为定位目标函数；$ {R_\alpha } $为第$\alpha $个台站理论到时和标记到时的残差；$ {t_\alpha } $为第$\alpha $个台站的P波标记到时；t0为震源的发震时刻；$ {T_\alpha }\left( {{x_0},{y_0},{{\textit{z}}_0}} \right) $为第$\alpha $个台站的理论到时；$ \left( {{x_0},{y_0},{{\textit{z}}_0}} \right) $为震源空间坐标；$ {v_{\text{P}}} $为P波平均波速；$\delta_{\mathrm{P} \alpha} $、δvα分别为第α个台站的人工标波误差因子和速度模型误差因子。
图12为2215工作面某一微震事件2 000次重定位结果及概率密度分布，试验关键参数见表2。由图12(a)可知使用数值仿真后重定位的80%震源在X轴和Y轴方向上均分布于距真实震源20 m范围之内，说明了数值仿真试验结果与微震系统定位误差一致。由图12(a)可计算得到重定位后微震事件的概率密度分布(图12(b))。由图12(b)可知，越靠近真实位置概率密度越高，离真实位置越远则概率密度越低，符合真实标波定位的情况，由概率密度分布可进一步计算重定位震源的分布概率。
表 2 数值仿真试验关键参数取值Table 2. Key parameters in numerical simulation test网格
尺寸/mP波平均波
速/(m·s−1)人工标波误差因子
$ {\delta _{\text{P}}} $的标准差速度模型误差因子
$ {\delta _{{v}}} $的标准差试验
次数20×20 4 500 0.006 100 2 000 获取震源机制和定位误差计算结果后，依据4.1节提出的基于震源破裂机制和定位误差校准的聚类依据，可对微震事件进行聚类分析，计算流程如图13所示。
![]()
图 13 基于震源破裂机制和定位误差校准的冲击危险预测指数IP计算流程Figure 13. IP calculation flowchart of rockburst risk based on source mechanism and location error calibration5. 方法验证和实践
本节首先对比分析了有无考虑震源破裂机制及定位误差的微震聚合效果；并采用冲击危险预测指数IP对2215工作面不同震源机制的大能量事件进行预测以验证提出的冲击危险预测方法的有效性和合理性。
5.1 基于震源破裂机制及定位误差校准的微震事件聚合效果分析
图14为2215工作面某次大能量事件(压缩−剪切破裂)发生前一周微震在是否考虑震源机制及定位误差的微震聚合效果。由图14(a)可知，未考虑震源破裂机制及定位误差时在工作面后方约200 m处识别出微震事件异常聚合区，微震聚合峰值为12，峰值区位于工作面侧但较不明显；考虑震源破裂机制与定位误差后(图14(b))，微震异常聚合区识别效率显著增加，微震聚合峰值为25，峰值区集中于采空区后方临空侧区域。对比2种结果不难发现，2种方法虽均可在采空区后方识别出微震异常聚合区，但在考虑震源破裂机制及定位误差条件下识别出微震异常聚合数更高，异常聚集区域也更明显；同时后者可识别出临空侧巷道(回风巷)微震高聚集区，对后续顶板破断诱发的大能量事件的预测更为有利。
![]()
图 14 是否考虑震源机制及定位误差的微震聚合效果Figure 14. Seismic clustering effect whether to consider the source mechanism and location error5.2 预测实例及效果验证
受篇幅所限，本节选取了2215工作面17次大能量事件中4次代表性的大能量事件进行预测效果验证。4次大能量事件具有不同的震源破裂机制，分别为压缩−剪切破裂、拉伸破裂、剪切破裂以及拉伸−剪切破裂机制。
如图15所示，顶板破断运移主导的大能量事件发生前在采空区后方均显现出不同程度的微震异常聚合区，且主要集中于临空侧。不同破裂类型大能量事件发生前微震异常聚合区与工作面的相对位置不同：压缩−剪切破裂的大能量事件发生前(图15(a))微震聚合区距离工作面较远(200 m左右)，而拉伸破裂(图15(b))、拉伸−剪切破裂(图15(d))的大能量事件发生前微震聚合区距离极近(工作面端头附近)，剪切破裂的大能量事件(图15(c))发生前微震聚合区距工作面较近(100 m左右)。
![]()
图 15 不同震源类型大能量事件发生前微震聚合情况Figure 15. Seismic clustering situation before high magnitude events of different source types说明距离工作面不同距离的采空区后方顶板破裂方式存在显著差异，距工作面由近及远存在拉伸、拉伸−剪切、剪切、压缩−剪切的转换过程，分别对应顶板的拉伸破断、剪切失稳以及回转挤压过程。此外，根据大能量事件与微震聚合区的对应关系来看，大能量事件虽然破裂机制不同，但是总发生于微震密集聚合区域；换言之，微震密集聚合区可预测大能量事件发生的高风险区域。然而，受不同采掘区域应力环境及微震部署的差异影响，微震事件活动强度不同，如图15所示，大能量事件发生区域微震峰值聚合数显著不同，若将微震聚合数作为冲击危险预测指标，难以确立合适的危险阈值。
因此，为检验基于震源破裂机制及定位误差校准的危险预测方法的有效性，采用4.2节中图13的计算流程计算得到冲击危险预测指数IP分布如图16所示。由图16可知，不同破裂机制的大能量事件均位于IP指数大于0.7的范围内，因此将0.7可作为利用IP指数预测大能量事件的阈值；同时说明在微震数据较为完整的情况下，IP指数可实现对大能量事件的有效预测，进而证明结合震源破裂机制及定位误差校准的危险预测指标可有效改善微震聚合分析对冲击地压危险预测的准确率，提高微震监测预警水平。
![]()
图 16 不同震源类型大能量事件的IP预测情况Figure 16. IP prediction of high magnitude events with different source types6. 结 论
(1) 求解了顶板破断诱发大能量事件的震源破裂机制，分析表明该大能量事件以拉剪以及压剪破坏为主；大能量事件发生前微震事件剪切成分占比不超过40%，顶板破断主导的微震事件走向主要分布于3°～80°以及150°～270°，倾角主要分布于70°～90°，集聚特征明显。
(2) 基于煤矿常见震源破裂机制及定位误差分析了微震事件聚合的可能性，建立了基于震源裂隙贯通原则下的微震聚合准则，构建了基于震源破裂机制及定位误差校准的危险预测方法及指标。
(3) 预测实例表明考虑震源破裂机制及定位误差后可显著增强微震异常聚合区的识别；微震数据较为完整的情况下，基于震源破裂机制及定位误差校准的危险预测指标IP阈值设置为0.7时可实现对大能量事件的准确预测。
震源破裂机制的深度解析是冲击地压危险源头判识及全时空智能预警重要组成部分。本文是将震源机制反演技术用于冲击地压危险预测成像的初步尝试，针对冲击危险源的智能预警和实时成像，未来如何进一步提高反演效率，研究适用于煤矿的震源机制实时反演技术是前提，提出综合岩层结构、地质构造、采掘布置、震源破裂等信息于一体的危险源实时成像技术是关键。
-
![]()
图 1 2215工作面及微震监测台网布置
Figure 1. Layout of 2215 working face and seismic monitoring network
![]()
图 5 震源破裂类型的CLVD−ISO钻石图
Figure 5. Representation of source mechanism using CLVD−ISO diamond plot
![]()
图 8 大能量事件发生前一周微震事件剪切成分占比
Figure 8. Proportion of shear components of seismic events in the week before the occurrence of high magnitude events
![]()
图 9 大能量事件发生前一周微震事件破裂面产状分布
注：左侧为破裂面1的走向分布，右侧为破裂面2的走向分布。
Figure 9. Distribution of rupture surface of seismic events occurred one week before the high magnitude events
![]()
图 11 结合震源破裂机制与定位误差的聚类可能性
注：红色和蓝色虚线圆分别代表震源A和B的震源半径。
Figure 11. Clustering possibility combining source mechanism and source location error
![]()
图 13 基于震源破裂机制和定位误差校准的冲击危险预测指数IP计算流程
Figure 13. IP calculation flowchart of rockburst risk based on source mechanism and location error calibration
![]()
图 14 是否考虑震源机制及定位误差的微震聚合效果
Figure 14. Seismic clustering effect whether to consider the source mechanism and location error
![]()
图 15 不同震源类型大能量事件发生前微震聚合情况
Figure 15. Seismic clustering situation before high magnitude events of different source types
![]()
图 16 不同震源类型大能量事件的IP预测情况
Figure 16. IP prediction of high magnitude events with different source types
表 1 不同震源机制反演方法及优缺点比较
Table 1 Comparison of different source mechanism inversion methods and their advantages and disadvantages
分类原则 反演方法 特点 优点 缺点 备注 按对格林函数处理方式不同 绝对法 根据理论或经验模型直接构建格林函数，以单个震源进行反演 无需基于额外假设、震源可单独求解 构建精确格林函数相对困难、误差不可控、对波形数量要求较高 绝对矩张量法、相对矩张量法、混合矩张量法均属于振幅法 相对法 基于震动波射线传播理论对格林函数进行转化与分解，以震源群为单位进行反演 避免直接求解格林函数、对波形数量要求相对较低 无法单独求解震源机制、误差不可控、要求参与反演的震源群具有多样的震源机制 混合法 融合绝对法的格林函数构建以及相对法的震源群思想，迭代修正绝对矩张量解 求解精确度高，误差可控，反演迭代速度快 构建精确格林函数相对困难 按所选用位移场的不同 振幅法 使用P波初动振幅作为位移场信息求解矩张量 仅使用P波初动振幅作为约束
条件，反演方程较为简单对波形信噪比要求较高，反演精度受P波初动拾取准度影响大 振幅比法 使用P波、S波初动极性及振幅比作为位移场信息求解矩张量 使用P波、S波极性及振幅比作为约束条件，一定程度上提高反演精度 对P波、S波到时及极性拾取准度要求高 全波形法 通过特定频率段实测波形与理论波形的拟合迭代确定最优的
震源机制解求解精确度高，误差可控 对波形信噪比要求高、计算量大反演周期长、精确格林函数构建困难 
下载: 导出CSV
表 2 数值仿真试验关键参数取值
Table 2 Key parameters in numerical simulation test
网格
尺寸/mP波平均波
速/(m·s−1)人工标波误差因子
$ {\delta _{\text{P}}} $的标准差速度模型误差因子
$ {\delta _{{v}}} $的标准差试验
次数20×20 4 500 0.006 100 2 000
下载: 导出CSV
-
[1] 朱建波,马斌文,谢和平,等. 煤矿矿震与冲击地压的区别与联系及矿震扰动诱冲初探[J]. 煤炭学报,2022,47(9):3396−3409. ZHU Jianbo, MA Binwen, XIE Heping, et al. Differences and connections between mining seismicity and coal bursts in coal mines and preliminary study on coal bursts induced by mining seismicity[J]. Journal of China Coal Society,2022,47(9):3396−3409.
[2] 曹安业,陈凡,刘耀琪,等. 冲击地压频发区矿震破裂机制与震源参量响应规律[J]. 煤炭学报,2022,47(2):722−733. CAO Anye, CHEN Fan, LIU Yaoqi, et al. Response characteristics of rupture mechanism and source parameters of mining tremors in frequent coal burst area[J]. Journal of China Coal Society,2022,47(2):722−733.
[3] 窦林名,田鑫元,曹安业,等. 我国煤矿冲击地压防治现状与难题[J]. 煤炭学报,2022,47(1):152−171. DOU Linming, TIAN Xinyuan, CAO Anye, et al. Present situation and problems of coal mines rock burst prevention and control in China[J]. Journal of China Coal Society,2022,47(1):152−171.
[4] 陈栋,王恩元,李楠. 千秋煤矿微震震源参数特征以及震源机制分析[J]. 煤炭学报,2019,44(7):2011−2019. CHEN Dong, WANG Enyuan, LI Nan. Analysis of microseismic source parameters and focal mechanism in Qianqiu Coal Mine[J]. Journal of China Coal Society,2019,44(7):2011−2019.
[5] 李铁,蔡美峰,孙丽娟,等. 基于震源机制解的矿井采动应力场反演与应用[J]. 岩石力学与工程学报,2016,35(9):1747−1753. LI Tie, CAI Meifeng, SUN Lijuan, et al. Inversion of mining-induced stress field and its application based on focal mechanism solution[J]. Chinese Journal of Rock Mechanics and Engineering,2016,35(9):1747−1753.
[6] LISLE R. New method of estimating regional stress orientations: application to focal mechanism data of recent british earthquakes[J]. Geophysical Journal International,2007,110:276−282.
[7] FUJII Y, ISHIJIMA Y, DEGUCHI G. Prediction of coal face rockbursts and microseismicity in deep longwall coal mining[J]. International Journal of Rock Mechanics and Mining Sciences,1997,34:85−96. doi: 10.1016/S1365-1609(97)80035-4
[8] MA K, YUAN F, WANG H, et al. Fracture mechanism of roof key strata in Dongjiahe Coal Mine using microseismic moment tensor[J]. Geomatics, Natural Hazards and Risk,2021,12(1):1467−1487. doi: 10.1080/19475705.2021.1933615
[9] SONG C, LU C, ZHANG X, et al. Moment tensor inversion and stress evolution of coal pillar failure mechanism[J]. Rock Mechanics and Rock Engineering,2022,55(4):2371−2383. doi: 10.1007/s00603-022-02783-1
[10] WU K, ZOU J, JIAO Y, et al. Focal mechanism of strong ground seismicity induced by deep coal mining[J]. Rock Mechanics and Rock Engineering,2023,56(1):779−795. doi: 10.1007/s00603-022-03079-0
[11] 曹安业, 刘耀琪, 杨旭, 等. 物理指标与数据特征融合驱动的冲击地压时序预测方法[J/OL]. 煤炭学报: 1−16[2023−05−11]. https://doi.org/10.13225/j.cnki.jccs.2022.0680. CAO Anye, LIU Yaoqi, YANG Xu, et al. Physical index and data fusion-driven method for coal burst prediction in time sequence[J/OL]. Journal of China Coal Society: 1−16[2023−05−11]. https://doi.org/10.13225/j.cnki.jccs.2022.0680.
[12] MEYER S G, BASSOM A P, READING A M. Delineation of fault segments in mines using seismic source mechanisms and location uncertainty[J]. Journal of Applied Geophysics,2019,170:103828. doi: 10.1016/j.jappgeo.2019.103828
[13] LURKA A. Spatio-temporal hierarchical cluster analysis of mining-induced seismicity in coal mines using ward's minimum variance method[J]. Journal of Applied Geophysics,2021,184:104249. doi: 10.1016/j.jappgeo.2020.104249
[14] FALMAGNE V. Quantification of rock mass degradation using microseismic monitoring and applications for mine design[D]. Kingston: Queen’s University, 2001.
[15] COULSON A L. Investigation of the pre to post peak strength state and behaviour of confined rock masses using mine induced microseismicity[D]. Toronto: University of Toronto, 2010.
[16] OUILLON G, DUCORBIER C, SORNETTE D. Automatic reconstruction of fault networks from seismicity catalogs: three-dimensional optimal anisotropic dynamic clustering[J]. Journal of Geophysical Research:Solid Earth,2008,113:B01306. doi: 10.1029/2007JB005032
[17] CESCA S, ŞEN AT, DAHM T. Seismicity monitoring by cluster analysis of moment tensors[J]. Geophysical Journal International,2014,196(3):1813−1826. doi: 10.1093/gji/ggt492
[18] 张全,邹俊鹏,吴坤波,等. 深部采煤上覆关键层破断诱发矿震特征研究[J]. 岩石力学与工程学报,2023,42(5):1150−1161. ZHANG Quan, ZOU Junpeng, WU Kunbo, et al. On the characteristics of mine earthquakes induced by key strata breaking during deep mining[J]. Chinese Journal of Rock Mechanics and Engineering,2023,42(5):1150−1161.
[19] VAVRYCUK V, KÜHN D. Moment tensor inversion of waveforms: A two-step time-frequency approach[J]. Geophysical Journal International,2012,190:1761−1776. doi: 10.1111/j.1365-246X.2012.05592.x
[20] 吴顺川,黄小庆,陈钒,等. 岩体破裂矩张量反演方法及其应用[J]. 岩土力学,2016,37(S1):1−18. WU Shunchuan, HUANG Xiaoqing, CHENG Fan, et al. Moment tensor inversion of rock failure and its application[J]. Rock and Soil Mechanics,2016,37(S1):1−18.
[21] AKI K, RICHARDS P. Quantitative seismology theory[M]. Mill Valley: University Science Books, 2002.
[22] LINZER L. A relative moment tensor inversion technique applied to seismicity induced by mining[J]. Rock Mechanics and Rock Engineering,2005,38:81−104. doi: 10.1007/s00603-004-0041-4
[23] VAVRYCUK V. Moment tensor decompositions revisited[J]. Journal of Seismology,2014,19:231−252.
[24] HUDSON JA, PEARCE RG, ROGERS RM. Source type plot for inversion of the moment tensor[J]. Journal of Geophysical Research:Solid Earth,1989,94(B1):765−774. doi: 10.1029/JB094iB01p00765
[25] RIEDESEL M, JORDAN T. Display and assessment of seismic moment tensors[J]. Bulletin of the Seismological Society of America,1989,79(1):85−100.
[26] 庞义辉,王国法,李冰冰. 深部采场覆岩应力路径效应与失稳过程分析[J]. 岩石力学与工程学报,2020,39(4):682−694. PANG Yihui, WANG Guofa, LI Bingbing. Stress path effect and instability process analysis of overlying strata in deep stopes[J]. Chinese Journal of Rock Mechanics and Engineering,2020,39(4):682−694.
[27] VAVRYCUK V. Tensile earthquakes: Theory, modeling, and inversion[J]. Journal of Geophysical Research (Solid Earth),2011,116:12320. doi: 10.1029/2011JB008770
[28] 曹安业. 采动煤岩冲击破裂的震动效应及其应用研究[D]. 徐州: 中国矿业大学, 2009. CAO Anye. Research on seismic effort of burst and failure of coal-rock mass associated with mining and its application[D]. Xuzhou: China University of Mining & Technology, 2009: 101-104.
[29] WANG C, SI G, ZHANG C, et al. Location error based seismic cluster analysis and its application to burst damage assessment in underground coal mines[J]. International Journal of Rock Mechanics and Mining Sciences,2021,143:104784. doi: 10.1016/j.ijrmms.2021.104784
-
期刊类型引用(2)
1. 窦林名，曹安业，杨耀，贺虎，杨垚鑫，白贤栖，顾倩悦，李松徽，付相超，顾颖诗，吴震，张帝. 巨厚覆岩矿震孕育破裂特征与应力触发机制. 煤田地质与勘探. 2024(10): 1-13 . 
百度学术
2. 邓志刚，李云鹏，赵善坤，王寅，李少刚，孔令海，蒋军军，苏振国，李一哲，秦凯. 冲击地压煤层智能安全高效开采系统架构与关键技术. 煤矿安全. 2024(12): 1-9 . 百度学术
其他类型引用(1)
计量
- 文章访问数: 112
- HTML全文浏览量: 17
- PDF下载量: 103
- 被引次数: 3
