CN114137600A - 一种利用微震监测数据反演岩石破裂机理及失稳预测方法 - Google Patents

Abstract

本发明公开了一种利用微震监测数据反演岩石破裂机理及失稳预测方法，根据微震波形的特征，确定微震波形的初至幅值与初至到时；基于微震波形初至到时，构建新的定位目标函数，求解破裂源的坐标，通过考虑微震波形数和初至到时与反算初至到时差值的标准差，剔除不符合要求的定位点；计算微震波形的理论法向位移值，通过理论法向位移值与测量位移值的差值表达式，求解破裂源矩张量并将其正则化，通过破裂面的性质直接判断岩石破裂机理；在此基础上，建立岩石破裂失稳预测指标，根据风险评价等级对岩石破裂失稳风险性进行预测，为工作面冲击地压、岩爆、煤与瓦斯突出突水、崩塌、滑坡这些灾害提供监测预警。

Description

一种利用微震监测数据反演岩石破裂机理及失稳预测方法

技术领域

本发明涉及岩体变形破裂监测评价技术领域，具体涉及一种利用微震监测数据反演岩石破裂机理及失稳预测方法。

背景技术

随着近地表煤炭资源逐渐枯竭，深部开采已成为一种必然趋势。深部条件下，煤岩受高地应力、高瓦斯及工作面采动影响，煤岩动力灾害日趋严重且复杂，严重影响煤矿的产能，威胁煤矿安全高效生产。准确监测预警是实现煤岩动力灾害有效防控的前提和基础。目前，微震监测技术由于其安装简单便捷、不损害被监测物体、监测空间范围大且可实现连续监测而被广泛使用。通过深入挖掘分析传感器采集的煤岩破裂产生的微震波形，可以确定煤岩破裂的空间坐标、得到煤岩周围应力状态、判断其破坏程度、获得破裂瞬间能量大小等，从而提前为可能发生的煤岩动力灾害发出预警。近些年来，大多数矿井都安装了微震监测系统，但煤矿生产环境复杂且煤层较具有松软、节理裂隙发育、非均匀性强的特点，致使岩爆、工作面冲击地压、煤与瓦斯突出等煤岩动力灾害机理仍不明确。

为此，不少学者开展了小尺度煤岩样加卸载破坏实验，研究煤岩破裂产生声发射信号的时空演化规律。声发射信号与地震波在传播方式、产生机理、波形信号的特征具有相似性，因此，煤岩小尺度破裂诱发产生的矿震事件又称为：“小尺度地震”。由于两者的相似性，许多学者直接采用地震学中的知识，研究小尺寸煤岩样破裂产生的信号，揭示煤岩破裂力学机理。使用最多的是地震矩张量理论，如中国专利申请CN112964787A公布了一种基于声发射的脆性材料裂纹类型检测方法，该方法通过求解震源矩张量的三个特征值，并根据其占比来确定震源破裂类型；中国专利申请CN113218766A公布了一种基于矩张量分析的岩石起裂应力与损伤应力辨识方法，中国专利申请CN113050159A公布了一种煤岩水力压裂裂缝微震定位及扩展机理监测方法，都采用矩张量分解后纯剪切部分的占比来判断岩石破裂机理，并没有考虑到分解部分中非线性偶极子CLVD意义指代不明，致使震源破裂类型判断结果值得质疑这一问题。此外在初至到时及振幅拾取方向，都采用了局部AIC法，容易出现多个局部最小值致使拾取结果不准确。中国专利申请CN106154307A公布了一种煤岩冲击失稳模式的微震识别方法，通过震源矩张量分析来获取破裂面产状，但并没有对煤岩冲击失稳模式做出任何的论述。上述文献都提到了如何求解矩张量，然而使用的方法单一，都采用了简化矩张量反演方法且针对岩石破裂机理但仅停留在理论和实验方面，并没有将其合理的利用起来，为岩石破裂失稳提供指导。事实上，岩石破裂机理及定位事件点蕴含丰富的信息。岩石破裂失稳前期，变形量较小，产生的裂缝较少，故定位事件点偏少，一旦失稳短时间内会产生大量裂缝，出现较多的定位事件点。此外，通过大量的实验室实验发现岩石破裂失稳时发生力学机理以剪切破裂为主，剪切破裂产生的微震波形能量较大。合理的利用这些信息可以对岩石破裂失稳风险性进行评价，为工作面冲击地压、岩爆、煤与瓦斯突出这些灾害提供监测预警。

因此，目前迫切需要一种利用微震监测数据反演岩石破裂机理及失稳预测方法。

发明内容

本发明的目的是提供一种利用微震监测数据反演岩石破裂机理及失稳预测方法。

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

一种利用微震监测数据反演岩石破裂机理及失稳预测方法，具体包括以下步骤：

a.根据实际的地质赋存条件及监测需要，确定微震传感器安装个数G并将其编号为第i＝1,2,3…G，同时利用已知的破裂源坐标标定微震传感器的坐标；

b.采集岩石破裂产生的微震波形并进行滤波；所述微震波形滤波完成后，基于波形特征，获取第i个微震传感器采集的微震波形的初至幅值Tⁱ _P及初至到时Bⁱ _p；

c.根据微震波形的初至到时、P波传播速度、微震传感器的坐标信息，构建破裂源定位的目标函数，并选择最先接收到微震波形的微震传感器作为初始体中心，通过迭代求解所述目标函数最小值来确定破裂源位置坐标，得到初始的破裂源定位点个数A；

d.将破裂源位置坐标作为输入数据，通过考虑同一微震事件的微震波形数、所述微震波形的初至到时与反算初至到时差值的标准差，对不符合条件的破裂源定位点进行剔除，得到最终破裂源定位点的个数H；

e.在微震传感器位置处进行锤击测试，采集锤击测试时产生的微震波形并得到其起跳幅值F_p，然后计算第i个微震传感器位置处微震波形引起的垂直法向位移值

f.计算所述垂直法向位移值

与锤击测试时产生微震波形的起跳幅值F_p的比值W；将比值W与岩石破裂产生微震波形的初至幅值Tⁱ _P相乘积，即得第i个微震传感器位置处岩石破裂产生微震波形的理论法向位移值Dⁱ _p；

g.计算第i个微震传感器位置处岩石破裂产生微震波形引起的测量位移值

，引入L2范数法，写出所述垂向理论位置值与测量位移值差值表达式，计算出破裂源矩张量M_jk，j＝1,2,3；k＝1,2,3；

h.在矩张量M_jk特征向量构成的直角坐标系下，将破裂源矩张量正则化，求取矩张量M_jk的特征值x_i，i＝1,2,3以及特征向量y_i，i＝1,2,3，然后利用特征值进一步获取破裂源形成的破裂面方向向量与法向向量的夹角；根据破裂面方向向量与法向向量的关系，确定破裂源发生的破裂机制；

i.通过考虑破裂源定位事件点的个数与所述破裂源发生的破裂机制，结合岩石的单轴抗压强度Rc、岩石的弹性能指数W_ET这四个预测指标，建立岩石破裂失稳预测指标W_i，i＝1，2，3，4，并利用风险评价等级W_ti对岩石破裂失稳危险性进行评价。

优选的，所述步骤b中，微震波形的初至幅值拾取方法如下：

首先根据微震波形特征，找出一个微震波形的最大幅值并确定其位于波峰还是波谷，然后选取前后相邻的波峰或波谷之间的距离作为时窗长度并记为L；确定时窗长度便可采用能量比法获取第j个时间窗口内其他微震传感器拾取的同一个微震事件对应的微震波形；

计算第j个时间窗口内同一个微震事件对应的所有微震波形的峰值度E及偏离度R；峰值度E及偏离度R的计算公式如下；

其中L为时窗长度(单位：ms×10⁴)，σ_x为微震波形的标准差，x_i为微震波形数据点，

为微震波形的平均值；

然后对比偏离度和峰值度的极值，同时计算极值点前峰度值与偏离度曲线斜率的最大值，其最大值位置处所对应的时间点即为初至时间；

最后在已知微震波形初至时间的基础上，进一步确定所述波形的初至幅值。

优选的，所述步骤c中，破裂源定位的目标函数的建立过程如下：

首先利用理论上每个微震传感器接收到微震波形的开始时间与实际微震波形从破裂源传播至微震传感器时间之间的关系，建立如下表达式：

其中，

为理论上第i微震传感器接收到微震波形的初至到时，ms；(x_i,y_i,z_i)表示微震传感器的坐标，mm；(x,y,z)表示破裂源的坐标，mm；v表示材料的纵波速度，mm/ms；

表示第i个微震传感器采集的微震波形在空间中的传播时间，ms；t₀表示破裂源的发生时刻；

然后由于现场环境的复杂性及煤岩的结构会存在断层及节理，同时获取的微震波形的初至到时Tⁱ _P会存在误差，使得求解各微震传感器间的破裂源的发生时刻存在误差，将β_i定义各微震传感器的误差值，其表达式如下：

其中β_i为各微震传感器的误差值，Tⁱ _P表示第i个微震传感器采集微震波形的初至到时，ms；

为理论上第i个微震传感器接收到微震波形的初至到时，ms。

最后将微震传感器的总误差定义为各微震传感器误差值的平方求和，同时对其进行微分运算，消除未知变量破裂源的发生时刻t₀即得破裂源定位的目标函数，其表达式如下：

其中，

为各微震传感器误差值的平方求和，G表示微震传感器的个数。

优选的，所述步骤d中，条件指的是微震事件的微震波形数不小于6个，所述微震波形的初至到时与反算初至到时差值的标准差在3μm之内。

优选的，所述步骤e中，所述垂直法向位移值

的计算公式如下：

其中，G(τ)为垂直法向位移的时间相关函数，T为锤击测试时产生的力，N；μ表示剪切模量，Pa；r为微震传感器与锤击点的距离，mm；c_p，c_s分别表示岩石的P波波速、S波波速，m/s。

优选的，所述步骤(g)中，所述测量位移值

的计算公式如下：

其中r为破裂源到微震传感器的空间距离m，ρ为岩石材料的密度kg/m³，α为P波波速m/s，R_P表示反射系数，

表示矩张量对时间的导数。

进一步的，所述步骤g中，差值表达式如下：

其中，G表示微震传感器的个数，Dⁱ _p为第i个微震传感器位置处岩石破裂产生微震波形的理论法向位移值，

为第i个微震传感器位置处岩石破裂产生微震波形的测量位移值；

根据差值表达式的平方和最小即可求解破裂源矩张量。

进一步的，所述步骤h中，破裂源发生的破裂机制是根据破裂面方向向量与法向向量夹角关系来确定，其中，破裂面方向向量与法向向量相互垂直时，为张拉破裂；破裂面方向向量与法向向量相互平行时，为剪切破裂；破裂面方向向量与法向向量既不平行也不垂直时，为混合破裂。

进一步的，所述步骤i中，所述岩石破裂失稳预测指标及风险评价等级具体如下：

将破裂源定位事件点增加的速率定义为事件增长率GR，并记为岩石破裂失稳的第一个预测指标W₁，其中GR的表达式如下：

其中，RP₂表示在t₂时刻破裂源定位事件点的个数，RP₁表示在t₁时刻破裂源定位事件点的个数；当0＜GR≤500个时，定义其评价指数为0；当500＜GR≤1500个时，定义其评价指数为1；当1500＜GR≤2000个时，定义其评价指数为2；当GR＞2000个时，定义其评价指数为3；

将定位事件点剪切破裂类型的占比SP，并记为岩石破裂失稳的第二个预测指标W₂；当0＜SP≤10％时，定义其评价指数为0；当10＜SP≤20％时，定义其评价指数为1；当20＜SP≤30％时，定义其评价指数为2；当SP＞30％，定义其评价指数为3；

将岩石的单轴抗压强度Rc，并记为岩石破裂失稳的第三个预测指标W₃，当0＜Rc≤7MP时，定义其评价指数为0；当7＜Rc≤10MP，定义其评价指数为1；当10＜Rc≤14MP时，定义其评价指数为2；当Rc＞14MP时，定义其评价指数为3；

将岩石的弹性能指数W_ET，并记为岩石破裂失稳的第四个预测指标W₄，当0＜W_ET≤2MP时，定义其评价指数为0；当2＜W_ET≤3.5MP时，定义其评价指数为1；3.5＜W_ET≤5MP时，定义其评价指数为2；当W_ET＞5时，定义其评价指数为3；

风险评价等级的计算公式如下：

其中W_i表岩石破裂失稳各项预测指标，n表示岩石破裂失稳预测指标的总个数，W_imax表示岩石破裂失稳的第i个预测指标对应的评价指数最大值；

当W_RE小于等于0.25时，定义为无失稳危险性；当0.25＜W_RE≤5时，定义为弱失稳危险性；当0.5＜W_RE≤0.75时，定义为中失稳危险性；当0.75＜W_RE≤1时，定义为强失稳危险性。

与现有技术相比，本发明具有如下有益效果：

(1)在拾取微震波形的初至幅值及初至到时，本发明利用微震波形的特征，提出自适应长度的时间窗口来扫描微震波形，可以更好的获取有效微震波形信号，避免遗失。此外，通过求解峰值度与偏离度的极小值，再根据极小值点前峰度值与偏离度曲线斜率的最大值，可进一步确定初至时间并获取初至幅值；与传统的方法相比，本发明有效地防止了低信噪比造成微震波形到时拾取不准确的问题，同时也解决了拾取过程中存在多个极值的问题，减少了人工拾取工作量，极大的提升拾取效率。

(2)在破裂源空间坐标确定的过程中，本发明提出全新的破裂源定位目标函数，消除了未知变量破裂源的发生时刻t₀，使得求解的过程从原来的四个变量降为三个变量，求解时更加稳健。此外，本发明给出了初始体中心为目标函数的迭代输入初始数据，极大的提高了运算的效率并且使得迭代的过程容易收敛，避免方程出现无解或病态方程的问题。

(3)本发明利用微震波形数及微震波形的初至到时与反算初至到时差值的标准差，进一步精准确定破裂源的坐标，真实的反映煤岩破裂的空间位置，刻画了破裂源的空间分布。这为后续破裂源力学机理的确定提供高精度输入坐标，使得破裂源力学机理的反演更加精确同时反映真实的破裂情况。

(4)与传统的简化矩张量反演方法相比，本发明提出一种新的求解破裂源矩张量分量方法。该方法将矩张量分量的求解过程转化为简单方程的求解问题，简化了计算过程，避免迭代时方程出现无解或病态方程的问题。

(5)本发明提出的反演岩石破裂力学机理的方法，克服了简化矩张量反演得到破裂源的走向、倾向、矩震级、滑移角等震源参数仅适用于发生剪切破裂断层的地震，其反演的结果同样适用于发生张拉破裂的断层，使得对裂缝破裂类型的判断更加全面具体，更具有说服力。

(6)相比于常规的确定破裂源力学机理的方法，本发明直接利用破裂面的性质，直观判断震源破裂类型，避免了使用矩张量分解后剪切破裂成分的占比，来定性判断破裂类型，解决了破裂源破裂机理解释不准确的问题。

(7)相比于传统的岩石破裂失稳预测方法，本发明充分考虑到岩石破裂的定位事件点及力学机理，结合岩石单轴抗压强度及弹性能指数，构建了多个预测指标综合评价岩石破裂失稳危险性。极大的改善了单一指标片面性强，不具有说服力的问题，同时使得评价的结果更加可靠。为边坡、桥隧的失稳、工作面冲击地压、岩爆提供监测预警。

附图说明

图1为本发明方法流程图；

图2为某煤矿工作面推进过程中第j(j＝18)个时间窗口内第i(i＝3)个微震传感器采集到的岩石破裂失稳诱发的一个微震波形自适应时间窗口长度示意图；

图3(a)为某煤矿工作面推进过程中第j(j＝18)个时间窗口内第i(i＝3)个微震传感器采集到的岩石破裂失稳诱发的一个微震波形；

图3(b)为该微震波形对应的峰度值曲线图；

图3(c)为该微震波形对应的偏离度曲线图；

图4某矿井岩石破裂失稳示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示为一种利用微震监测数据反演岩石破裂机理及失稳预测方法流程图。

本实施例以某矿井为例，根据实际的地质赋存条件及工作面回采情况，确定安装12个微震传感器并将其标号为第i＝1,2,3…12(微震传感器的个数不能少于6，这是由于破裂源矩张量中有6个未知数，只有6个微震传感器采集的微震波形方可形成方程组)，同时利用人工爆破的方法，确实12个微震传感器的坐标依序为：(3666.5,4347.4,-998.8)、(3661.7,4893.2,-1094.2)、(4021.1,4572.2,-1034.8)、(4021.3,4873.5,-1023.2)、(6234.7,7200.6,-893.2)、(4423.5,5123.4,-1023.7)、(6218.3,7323.2,-677.4)、(6212.3,7813.2,-921)、(6327.8,6878.9,-721.1)、(6637.7,8778.3,-871.4)、(6127.6,5328.9,-679.9)、(7327.8,7893.5,-811.1)。

采集工作面生产过程中半个月内产生的微震波形并进行陷波处理，滤除掉所述微震波形中的环境噪音。

微震波形陷波完成后，根据微震波形特征，找出一个微震波形的最大幅值，然后选取前后相邻的波峰或波谷之间的距离作为时窗长度并记为L。如图2所示，为某煤矿工作面推进过程中第j(j＝18)个时间窗口内第i(i＝3)个微震传感器采集到的岩石破裂失稳诱发的一个微震波形，所示微震波形的最大峰值点相邻两个波谷的长度即为该微震波形的时窗长度L＝0.2s。

如图3(a)所示，针对某煤矿工作面推进过程中第j(j＝18)个时间窗口内第i(i＝3)个微震传感器采集到的岩石破裂失稳诱发的一个微震波形，该微震波形的长度为2.5s。通过以下公式计算该微震波形的峰值度E和偏离度R曲线，峰值度E及偏离度R的计算公式如下；

从图3(b)可知该微震波形峰值度曲线的极大值为11.1mv，从图3(c)可知该微震波形偏离度曲线的极小值为0.9mv。通过对比偏离度和峰值度的极值，同时计算极值点前峰度值与偏离度曲线斜率的最大值。峰值度曲线0.92s位置处，曲线的斜率最大。因此，该微震波形的初至时间T³ _P为0.96s，初至幅值B³ _p为0.089mv。

重复上述步骤，获得第j(j＝18)个时间窗口内岩石破裂失稳诱发同一个微震事件对应所有微震波形的初至时间和初至幅值分别为T¹ _P＝0.79、T² _P＝0.83T⁴ _P＝0.98、T⁵ _P＝0.89、T⁶ _P＝0.88、T⁷ _P＝0.99、T⁸ _P＝0.86、T⁹ _P＝0.86、T¹⁰ _P＝0.88、T¹¹ _P＝0.93、T¹² _P＝0.99s；B¹ _p＝0.088、B² _p＝0.083、B⁴ _p＝0.079、B⁵ _p＝0.085、B⁶ _p＝0.088、B⁷ _p＝0.089、B⁸ _p＝0.087、B⁹ _p＝0.090、B¹⁰ _p＝0.083、B¹¹ _p＝0.082、B¹² _p＝0.087。

依据上述方法计算某煤矿工作面推进过程中，半个月内岩石破裂诱发所有微震波形的初至时间与初至幅值。

利用理论上微震传感器接收到微震波形的初至到时与实际微震波形从破裂源传播至微震传感器时间之间的关系，建立如下表达式：

由于现场环境的复杂性及煤岩的结构会存在断层及节理，同时获取的微震波形数据的初至到时Tⁱ _P会存在误差，使得求解各微震传感器间的破裂源的发生时刻存在误差，将β_i定义各微震传感器的误差值，其表达式如下：

将微震传感器的总误差定义各微震传感器误差值的平方求和，同时对其进行微分运算，消除未知变量破裂源的发生时刻t₀即得破裂源定位的目标函数，其表达式如下：

选择最先接收到微震波形的微震传感器作为初始体中心，通过迭代求解所述目标函数最小值来确定破裂源位置坐标，然后得到破裂源定位点的个数A＝4653个。选取最先接收到微震波形的微震传感器作为初始体中心其目的在于，基于初始体中心，利用随机函数可生成空间的四面体边长，进而为目标函数的初始迭代提供输入数据。

将破裂源位置坐标作为输入数据，通过考虑同一微震事件的微震波形数不小于6个，所述微震波形的初至到时与反算初至到时差值的标准差在3μm之内，对不符合条件的破裂源定位点进行剔除，最终破裂源定位点的个数H＝2376个。

接下来在第i(i＝3)个微震传感器位置处进行锤击测试，采集锤击测试时产生的一个微震波形并得到其起跳幅值F_p＝62.7mv，然后计算第i(i＝3)微震传感器位置处该微震波形引起的垂直法向位移值

计算公式如下：

具体的，G(0.62)＝0.02，T＝16.83N，c_p/c_s＝2.51，r＝27.25mm，μ＝30.89Pa，D³ _L＝12.81×10^-12m；

计算所述垂直法向位移值D³ _L＝12.81×10^-12m与锤击测试时产生微震波形的起跳幅值值F_p＝62.7mv的比值W＝0.204×10^-12m/mV；将比值W＝0.204×10^-12m/mV与岩石破裂产生微震波形数据的初至幅值相乘积，即得第i个微震传感器位置处岩石破裂诱发微震波形的理论法向位移值Dⁱ _p；

本实施例中，第j(j＝18)个时间窗口内岩石破裂失稳诱发同一个微震事件对应所有微震波形的理论法向位移值D¹ _L＝11.83×10^-12m、D² _L＝17.7×10^-12m、D³ _L＝12.81×10^-12m、D⁴ _L＝2.89×10^-12m、D⁵ _L＝11.3×10^-12m、D⁶ _L＝21.3×10^-12mD⁷ _L＝19.8×10^-12m、D⁸ _L＝6.6×10^-12m、D⁹ _L＝7.8×10^-12m、D¹⁰ _L＝7.9×10^-12m、D¹¹ _L＝6.54×10^-12m、D¹² _L＝25.53×10^-12m。

计算第j(j＝18)个时间窗口内岩石破裂失稳诱发同一个微震事件对应所有微震波形的测量位移值

计算公式如下：

本实施例中，岩石的密度为ρ＝11.29kg/m³，P波波速为4km/s，其中第j(j＝18)个时间窗口内岩石破裂失稳诱发同一个微震事件的坐标为X＝4234.6m，Y＝6895.1m，Z＝-985.3m。

引入L2范数法，写出所述垂向理论位置值与测量位移值差值表达式进行，差值表达式如下：

根据差值表达式的平方和最小即可求解破裂源矩张量。求得一个破裂源矩张量为

在矩张量M_jk特征向量构成的直角坐标系下，将破裂源矩张量正则化，求得该矩张量的特征值，x₁＝-2.2022，x₂＝4.3467，x₃＝16.9356，特征向量y₁＝[-0.4360,0.0684,0.8873]^T，y₂＝[0.2815,-0.9367,0.2082]^T，y₃＝[0.8548,0.3433,0.3892]^T，根据特征值之间的关系，可计算出破裂面方向向量与法向向量的夹角为0°，判断该破裂点的力学机理为剪切破裂。

重复上述所有步骤，分别计算剩余的2375个破裂源的理论法向位移值与测量位置值，并进一步求解出2375个破裂源的矩张量，然后根据破裂面方向向量与法向向量的夹角，确定破裂源发生的破裂机制。破裂面方向向量与法向向量相互垂直时，为张拉破裂；破裂面方向向量与法向向量相互平行时，为剪切破裂；破裂面方向向量与法向向量既不平行也不垂直时，为混合破裂。

通过考虑定位事件点的个数与所述破裂源发生的破裂机制，将破裂源定位事件点增加的速率定义为事件增长率GR，并记为岩石破裂失稳的第一个预测指标W₁，GR的表达式如下：

当0＜GR≤500个时，定义其评价指数为0；当500＜GR≤1500个时，定义其评价指数为1；当1500＜GR≤2000个时，定义其评价指数为2；当GR＞2000个时，定义其评价指数为3；

将定位事件点剪切破裂类型的占比SP记为岩石破裂失稳的第二个预测指标W₂；当0＜SP≤10％时，定义其评价指数为0；当10＜SP≤20％时，定义其评价指数为1；当20＜SP≤30％时，定义其评价指数为2；当SP＞30％，定义其评价指数为3；

将岩石的单轴抗压强度Rc，并记为岩石破裂失稳的第三个预测指标W₃；当0＜Rc≤7MP时，定义其评价指数为0；当7＜Rc≤10MP，定义其评价指数为1；当10＜Rc≤14MP时，定义其评价指数为2；当Rc＞14MP时，定义其评价指数为3；

将岩石的弹性能指数W_ET，并记为岩石破裂失稳的第四个预测指标W₄，当0＜W_ET≤2MP时，定义其评价指数为0；当2＜W_ET≤3.5MP时，定义其评价指数为1；3.5＜W_ET≤5MP时，定义其评价指数为2；当W_ET＞5时，定义其评价指数为3。

利用风险评等级W_RE对岩石破裂失稳危险性进行评价，风险评价等级的计算公式如下：

当W_RE小于等于0.25时，定义为无失稳危险性；当0.25＜W_RE≤5时，定义为弱失稳危险性；当0.5＜W_RE≤0.75时，定义为中失稳危险性；当0.75＜W_RE≤1时，定义为强失稳危险性。

如图4所示，图中的三角形表示发生张拉破裂的定位点，实心原点表示发生剪切破裂的定位点，正方形表示发生混合破裂的定位点。该矿井半个月内，破裂源定位事件点从300增加至3520个，其中增加的3220个破裂源定位点中剪切破裂占比为59％；此外该矿井岩石的单轴抗压强度Rc为13.32MPa，岩石的弹性能指数W_ET＝3.13。

经判定，该矿井在某一时刻，岩石破裂失稳的第一个预测指标W₁对应的评价指数为3，岩石破裂失稳的第二个预测指标W₂对应的评价指数为3，岩石破裂失稳的第三个预测指标W₃对应的评价指数为2，岩石破裂失稳的第四个预测指标W₄对应的评价指数为1，因此该时刻矿井的风险评价等级W_RE＝0.75，为中失稳风险，该矿井必须采取相应的防范措施，预防岩石破裂失稳导致的事故。

Claims (9)

1.一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于，该方法包括以下步骤：

a.根据实际的地质赋存条件及监测需要，确定微震传感器安装个数G并将其编号为第i＝1,2,3…G，G≥6，同时利用已知的破裂源坐标标定微震传感器的坐标；

b.采集岩石破裂产生的微震波形并进行滤波；所述微震波形滤波完成后，基于波形特征，获取第i个微震传感器采集的微震波形的初至幅值Tⁱ _P及初至到时Bⁱ _p；

c.根据微震波形的初至到时、P波传播速度、微震传感器的坐标信息，构建破裂源定位的目标函数，并选择最先接收到微震波形的微震传感器作为初始体中心，通过迭代求解所述目标函数最小值来确定破裂源位置坐标，得到初始的破裂源定位点个数A；

d.将破裂源位置坐标作为输入数据，通过考虑同一微震事件的微震波形数、微震波形的初至到时与反算初至到时差值的标准差，对不符合条件的破裂源定位点进行剔除，得到最终破裂源定位点的个数H；

e.在微震传感器位置处进行锤击测试，采集锤击测试时产生的微震波形并得到其起跳幅值F_p，然后计算第i个微震传感器位置处微震波形引起的垂直法向位移值

f.计算所述垂直法向位移值

与锤击测试时产生微震波形的起跳幅值F_p的比值W；将比值W与岩石破裂产生微震波形的初至幅值Tⁱ _P相乘积，即得第i个微震传感器位置处岩石破裂产生微震波形的理论法向位移值Dⁱ _p；

g.计算第i个微震传感器位置处岩石破裂产生微震波形引起的测量位移值

引入L2范数法，写出所述理论法向位移值Dⁱ _p与测量位移值

差值表达式，并计算出破裂源矩张量M_jk，j＝1,2,3；k＝1,2,3；

h.在矩张量M_jk特征向量构成的直角坐标系下，将破裂源矩张量正则化，求取矩张量M_jk的特征值x_i，i＝1,2,3以及特征向量y_i，i＝1,2,3，然后利用特征值进一步获取破裂源形成的破裂面方向向量与法向向量的夹角；根据破裂面方向向量与法向向量的关系，确定破裂源发生的破裂机制；

i.通过考虑破裂源定位事件点的个数与剪切破裂类型占比，结合岩石的单轴抗压强度、岩石的弹性能指数，建立岩石破裂失稳预测指标，并利用风险评价等级对岩石破裂失稳危险性进行预测。

2.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：所述步骤b中，微震波形的初至幅值拾取方法如下：

首先根据微震波形特征，找出一个微震波形的最大幅值并确定其位于波峰还是波谷，然后选取前后相邻的波峰或波谷之间的距离作为时窗长度并记为L；采用能量比法获取第j个时间窗口内其他微震传感器拾取的同一个微震事件对应的微震波形；

计算第j个时间窗口内同一个微震事件对应的所有微震波形的峰值度E及偏离度R；峰值度E及偏离度R的计算公式如下；

其中L为时窗长度，s；σ_x为微震波形的标准差；x_i为微震波形数据点；

为微震波形的平均值；

然后对比偏离度和峰值度的极值，同时计算极值点前峰度值与偏离度曲线斜率的最大值，其最大值位置处所对应的时间点即为初至时间；

最后在已知微震波形初至时间的基础上，进一步确定所述微震波形的初至幅值。

3.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：所述步骤c中，破裂源定位的目标函数的建立过程如下：

首先利用理论上每个微震传感器接收到微震波形的初至到时与实际微震波形从破裂源传播至微震传感器时间之间的关系，建立如下表达式：

其中，

为理论上第i个微震传感器接收到微震波形的初至到时，ms；x_i,y_i,z_i表示微震传感器的坐标，mm；x,y,z表示破裂源的坐标，mm；v表示材料的纵波速度，mm/ms；

表示第i个微震传感器采集的微震波形在空间中的传播时间，ms；t₀表示破裂源的发生时刻；

然后计算各微震传感器的误差值β_i，其表达式如下：

其中β_i为各微震传感器的误差值，Tⁱ _P表示第i个微震传感器采集微震波形的初至到时，ms；

最后将微震传感器的总误差定义为各微震传感器误差值的平方求和，同时对其进行微分运算，消除未知变量破裂源的发生时刻t₀即得破裂源定位的目标函数，其表达式如下：

其中，∑β_i ²为各微震传感器误差值的平方求和，G表示微震传感器的个数。

4.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：所述步骤d中，条件指的是微震事件的微震波形数不小于6个，所述微震波形的初至到时与反算初至到时差值的标准差在3μm之内。

5.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：所述步骤e中，锤击测试产生微震波形的垂直法向位移值

的计算公式如下：

其中，G(τ)为垂直法向位移的时间相关函数，T为锤击测试时产生的力，N；μ表示剪切模量，Pa；r为微震传感器与锤击点的距离，mm；c_p，c_s分别表示岩石的P波波速、S波波速，m/s。

6.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：所述步骤g中，所述测量位移值

的计算公式如下：

其中r为破裂源到微震传感器的空间距离，m；ρ为岩石材料的密度，kg/m³；α为P波波速，m/s；R_P表示反射系数，

表示矩张量对时间的导数。

7.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：所述步骤g中，差值表达式如下：

其中，G表示微震传感器的个数，Dⁱ _p为第i个微震传感器位置处岩石破裂产生微震波形的理论法向位移值，

为第i个微震传感器位置处岩石破裂产生微震波形的测量位移值；

根据差值表达式的平方和最小求解破裂源矩张量。

8.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：所述步骤h中，破裂源发生的破裂机制是根据破裂面方向向量与法向向量夹角关系来确定，其中，破裂面方向向量与法向向量相互垂直时，为张拉破裂；破裂面方向向量与法向向量相互平行时，为剪切破裂；破裂面方向向量与法向向量既不平行也不垂直时，为混合破裂。

9.根据权利要求1所述一种利用微震监测数据反演岩石破裂机理及失稳预测方法，其特征在于：步骤i中，所述岩石破裂失稳预测指标及风险评价等级具体如下：

将破裂源定位事件点增加的速率定义为事件增长率GR，并记为岩石破裂失稳的第一个预测指标W₁，GR的表达式如下：

其中，RP₂表示在t₂时刻破裂源定位事件点的个数，RP₁表示在t₁时刻破裂源定位事件点的个数；

将定位事件点剪切破裂类型的占比SP，并记为岩石破裂失稳的第二个预测指标W₂；将岩石的单轴抗压强度Rc，并记为岩石破裂失稳的第三个预测指标W₃；将岩石的弹性能指数W_ET记为岩石破裂失稳的第四个预测指标W₄；

风险评价等级的计算公式如下：

其中W_i表岩石破裂失稳各个预测指标，n表示岩石破裂失稳预测指标的总个数，W_imax表示岩石破裂失稳的第i个预测指标对应的评价指数最大值；

当W_RE小于等于0.25时，定义为无失稳危险性；当0.25＜W_RE≤5时，定义为弱失稳危险性；当0.5＜W_RE≤0.75时，定义为中失稳危险性；当0.75＜W_RE≤1时，定义为强失稳危险性。

Images