一种煤岩冲击失稳模式的微震识别方法

摘要

一种煤岩冲击失稳模式的微震识别方法，属于冲击矿压微震特性分析技术领域。其特征在于：包括如下步骤：步骤1001，在同一采场不同位置设置多个微震记录仪对采场进行微震监测，每个微震记录仪形成一个微震信号记录；步骤1002，在微震记录仪监测到的微震记录中，筛选多个微震波形；步骤1003，根据震源远场位移计算步骤，计算每个微震波形的震源远场位移；步骤1004，计算得到震源矩张量；步骤1005，对步骤1004中计算得到的震源矩张量进行分解；步骤1006，对步骤1005中分解后的震源矩张量进行分析得到：震源破裂模式，震源P、T轴以及震源破裂面产状。

说明书

技术领域

一种煤岩冲击失稳模式的微震识别方法，属于冲击矿压微震特性分析技术领域。

背景技术

冲击矿压灾害是一种开采诱发的矿山地震，不仅造成井巷破坏、人员伤害、地面建筑物破坏，而且会引起瓦斯、煤尘爆炸。由于这种灾害发生时间、地点、位置等的复杂多样性和突发性，对其防治，特别是预测是世界性的难题。

目前，冲击矿压的监测预警方法主要包括：微震监测、电磁辐射监测、声发射监测、钻屑量监测、工作面矿压监测、采动应力监测等。尤其是微震监测方法，能够对全矿范围进行实时监测，是一种区域性、及时监测手段，能够给出震动后的各种信息，具有不损伤煤体、劳动强度小、时间和空间连续等优点。该技术目前被公认为煤岩动力灾害，特别是对于煤矿动力现象监测最有效和最有发展潜力的监测方法之一。在微震监测的基础上，目前用矩张量方法分析煤矿采场微震震源机制主要存在以下问题：矩张量分解分量物理意义不明确、适用性有限，冲击矿压事件中存在大量非剪切破裂形式震源，需要考虑破裂面运动方向与破裂面不共面等问题。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种在综合考虑煤矿采场尺度的基础上，对微震波形数据进行筛选，定量描述冲击矿压震源破裂形式、震源应力状态、破裂面产状的煤岩冲击失稳模式的微震识别方法。

本发明解决其技术问题所采用的技术方案是：该煤岩冲击失稳模式的微震识别方法，其特征在于：包括如下步骤：

步骤1001，在同一采场不同位置设置多个微震记录仪对采场进行微震监测，每个微震记录仪形成一个微震信号记录；

步骤1002，在微震记录仪监测到的微震信号记录中，根据波形筛选原则筛选出多个微震波形；

步骤1003，根据震源远场位移计算步骤，计算每个微震波形的震源远场位移；

步骤1004，计算得到震源矩张量；

步骤1005，对步骤1004中计算得到的震源矩张量进行分解；

步骤1006，对步骤1005中分解后的震源矩张量进行分析得到：震源破裂模式，震源P、T轴以及震源破裂面产状。

优选的，步骤1003中所述的震源远场位移计算步骤，包括：

步骤a1：根据所述的步骤1002中筛选出的微震波形，拾取每一个微震波形的P波初振幅值，并记作U1；

步骤a2，在筛选的微震波形中剪切P波的时域波形；

步骤a3，根据微震记录仪固有采样频率将步骤a2中剪切的P波时域波形进行快速傅里叶变换，转化为频域波形；

步骤a4，对步骤a3中转化得到的频域波形利用如下公式进行衰减修正：

$A^{new}\left(f\right)=A\left(f\right)\exp \left(\frac{\pi fD}{vQ}\right)$

其中，A(f)为时域速度谱进行快速傅里叶变换的结果；f为步骤a3中转换得到的频域波形的频率；v为P波速度；Q为衰减因子；

步骤a5，在步骤a4中计算得到的衰减修正值得到地动速度功率谱并根据地动速度功率谱得出地动位移功率谱

$D_0^2\left(f\right)=\frac{v_0^2\left(f\right)}{{\left(2\pi f\right)}^2}$

步骤a6，计算得到地动速度功率谱积分S_V2和地动位移功率谱积分S_D2：

$S_{V2}=2{\int }_0^{\infty }{V}_0^2\left(f\right)df$

$S_{D2}=2{\int }_0^{\infty }{D}_0^2\left(f\right)df$

其中：为考虑自由面影响乘以1/4的修正速度功率谱；为对应的位移功率谱；

步骤a7，计算得到各个微震波形的低频位移U2：

$U_2=\sqrt{\frac{4S_{D2}^{3/2}}{S_{V2}^{1/2}}}.$

优选的，步骤1004中所述的震源矩张量为：

$\left(\begin{array}{c}{M}_{11}\\ M_{12}\\ M_{13}\\ M_{22}\\ M_{23}\\ M_{33}\end{array}\right)=4{{\mathrm{\pi \rho rv}}_p}^3{\left(\begin{array}{cccccc}{\gamma }_3^1{\gamma }_1^1{\gamma }_1^1& 2{\gamma }_3^1{\gamma }_1^1{\gamma }_2^1& 2{\gamma }_3^1{\gamma }_1^1{\gamma }_3^1& {\gamma }_3^1{\gamma }_2^1{\gamma }_2^1& 2{\gamma }_3^1{\gamma }_2^1{\gamma }_3^1& {\gamma }_3^1{\gamma }_3^1{\gamma }_3^1\\ {\gamma }_3^2{\gamma }_1^2{\gamma }_1^2& {\gamma }_3^2{\gamma }_1^2{\gamma }_2^2& 2{\gamma }_3^2{\gamma }_1^2{\gamma }_3^2& {\gamma }_3^2{\gamma }_2^2{\gamma }_2^2& 2{\gamma }_3^2{\gamma }_2^2{\gamma }_3^2& {\gamma }_3^2{\gamma }_3^2{\gamma }_3^2\\ {\gamma }_3^3{\gamma }_1^3{\gamma }_1^3& 2{\gamma }_3^3{\gamma }_1^3{\gamma }_2^3& 2{\gamma }_3^3{\gamma }_1^3{\gamma }_3^3& {\gamma }_3^3{\gamma }_2^3{\gamma }_2^3& 2{\gamma }_3^3{\gamma }_2^3{\gamma }_3^3& {\gamma }_3^3{\gamma }_3^3{\gamma }_3^3\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ {\gamma }_3^n{\gamma }_1^n{\gamma }_1^n& 2{\gamma }_3^n{\gamma }_1^n{\gamma }_2^n& 2{\gamma }_3^n{\gamma }_1^n{\gamma }_3^n& {\gamma }_3^n{\gamma }_2^n{\gamma }_2^n& 2{\gamma }_3^n{\gamma }_2^n{\gamma }_3^n& {\gamma }_3^n{\gamma }_3^n{\gamma }_3^n\end{array}\right)}^{-1}\left(\begin{array}{c}{u}_p^1\\ u_p^2\\ u_p^3\\ ·\\ ·\\ ·\\ ·\\ ·\\ ·\\ u_p^n\end{array}\right)$

其中：v_p为P波传播速度；r为震源到微震记录仪的距离；ρ为岩石密度；γ_i为震源至微震记录仪的震动波射线对应于各坐标轴的分量，即γ_i＝(x_i-x_0i)/r(x_i为微震记录仪各坐标分量，x_0i为震源各坐标分量，i＝1,2,3)；M_ij为作用于震源的矩张量。

优选的，所述的震源矩张量分解为以下各部分：一个各向同性部分、一个补偿线性矢量偶极子部分及一个双力偶部分，即：

$\left(\begin{array}{c}M=PI-{M_1}^{,}\left(\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 2\end{array}\right)+\left({M_3}^{,}+2{M_1}^{,}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& -1& 0\\ 0& 0& 1\end{array}\right)\\ =PI-{M_1}^{,}\left(2a_3{a}_3-a_2{a}_2-a_1{a}_1\right)+\left({M_3}^{,}+2{M_1}^{,}\right)\left(a_3{a}_3-a_2{a}_2\right)\\ =PI+{M_3}^{,}F\left(2a_3{a}_3-a_2{a}_2-a_1{a}_1\right)+{M_3}^{,}\left(1-2F\right)\left(a_3{a}_3-a_2{a}_2\right)\end{array}\right)$

式中PI为各项同性部分，M₃'F(2a₃a₃-a₂a₂-a₁a₁)为补偿线型偶极部分，

M₃'(1-2F)(a₃a₃-a₂a₂)为双力偶部分；

其中：F＝-M₁'/M₃'，且0≤F≤1/2，a₁，a₂，a₃为矩阵M的本征矢量，M₁，M₂和M₃为对应的特征值，P＝(M₁+M₂+M₃)/3，I为单位矩阵，M′表示矩张量的偏张量部分，本征值为M_i'＝M-P，i＝1,2,3。

优选的，在所述的步骤1005中，所述的震源破裂模式的计算公式为：

DC％＝(1-2*F)*(1-norm(Miso)/norm(M))

其中：norm为求解矩阵最大奇异值函数，M为所述的步骤1004中得到的震源矩张量，Miso是震源矩张量M中的各项同性部分，F＝-M₁'/M₃'，M′表示震源矩张量M的偏张量部分。

优选的，在所述的步骤1006中，所述的震源P、T轴的求取方法为：

将所述的震源矩张量进行分解后得到震源矩张量的偏张量部分M’，在柱坐标(b₁，b₂，b₃)中，震源矩张量的偏张量部分M’对角化为：

$M^{\prime }=\left(\begin{array}{ccc}{M_{11}}^{\prime }& {M_{12}}^{\prime }& {M_{13}}^{\prime }\\ {M_{21}}^{\prime }& {M_{22}}^{\prime }& {M_{23}}^{\prime }\\ {M_{31}}^{\prime }& {M_{32}}^{\prime }& {M_{33}}^{\prime }\end{array}\right)=\left(\begin{array}{ccc}{M_1}^{\prime }& 0& 0\\ 0& {M_2}^{\prime }& 0\\ 0& 0& {M_3}^{\prime }\end{array}\right)$

其中，b₁，b₂，b₃为矩阵M′的本征矢量，M₁′，M₂′和M₃′为对应的特征值。假设有M₁′<M₂′<M₃′，则b₁对应拉张方向轴T，b₃对应压缩方向轴P。

优选的，在所述的步骤1006中，所述的震源破裂面产状的求取方法包括如下步骤：

步骤b1，震源矩张量可以由破裂面的位置与运动矢量表示为：

式中：i，j，k＝1，2，3，为空间坐标系中的3个方向；Mij为地震矩张量形式，根据i，j取值不同，分别代表矩张量不同分量；u为破裂面运动方向的位移量；S为破裂面表面积；λ和μ为拉梅常数；v为运动方向，割据i,j,k取值不同，vi,vj,vk表示破裂面运动矢量在各坐标系方向上的分量；n为破裂面法向方向，其分量取值意义通v；

步骤b2，对步骤b1的公式进行本征值化，得到：

$M=uS\left(\begin{array}{ccc}(\lambda +\mu )nv+\mu & 0& 0\\ 0& \lambda nv& 0\\ 0& 0& (\lambda +\mu )nv-\mu \end{array}\right)$

则有(假设M1>M2>M3)：

步骤b3，根据步骤b2得到的特征矢量与破裂面运动方向和方向的关系：

式中，e₁⊥e₂⊥e₃，绝对值符号表示矢量大小；×表示矢量乘法；e₁、e₂、e₃分别为矩张量的最大特征值、中间特征值和最小特征值；

步骤b4，设步骤b3中得到的矢量v和矢量n的夹角为β，v和e1的夹角以及n和e1的夹角均为β/2，得到：

$cos\frac{\beta }{2}=\sqrt{\frac{M_1-M_2}{M_1-M_3}}$

$sin\frac{\beta }{2}=\sqrt{\frac{M_2-M_3}{M_1-M_3}}$

并得到破裂面运动方向和法向方向与矩张量最大特征值对应特征矢量和最小特征值对应特征矢量的关系：

$n=cos\frac{\beta }{2}{e}_1+sin\frac{\beta }{2}{e}_3=\sqrt{\frac{M_1-M_2}{M_1-M_3}}{e}_1+\sqrt{\frac{M_2-M_3}{M_1-M_3}}{e}_3$

$v=cos\frac{\beta }{2}{e}_1-sin\frac{\beta }{2}{e}_3=\sqrt{\frac{M_1-M_2}{M_1-M_3}}{e}_1-\sqrt{\frac{M_2-M_3}{M_1-M_3}}{e}_3$

步骤b5，根据破裂面法向方向的空间矢量值到破裂面的几何方程表达式，并确定破裂面的方位角和倾角的计算公式：

方位角：

式中，n(3)为破裂面法向矢量Z轴分量，|n|表示法向向量的模

若dip>90°，则破裂面方位角为取值为(π-dip)

倾角:

式中，n(1)为破裂面法向向量X轴分量，n(2)为破裂面法向向量Y轴分量。

优选的，在所述的步骤1001中，微震记录仪的数量至少为6个；在所述的步骤1002中，至少筛选6个微震波形。

优选的，在执行完所述的步骤1006后，绘制微震的震源球。

优选的，在步骤1002中所述的波形筛选原则为：在距离同一采场2Km范围内筛选至少6个微震波形。

与现有技术相比，本发明所具有的有益效果是：

1、本煤岩冲击失稳模式的微震识别方法，在对远场位移优化的基础上，利用矩张量表示震源受力，无需事先对震源机制做出任何假定，且远场位移与矩张量是线型关系式，通过矩张量分解计算，定量描述了震源破裂形式、震源应力状态、破裂面产状。

2、本煤岩冲击失稳模式的微震识别方法中，选择矿震波形低频位移振幅表征矿震远场位移,该计算方法无需任何源模型假设,且波形双对数频率域的振幅值在低频部分是相对平直的常数,其大小正比于波形标量矩,提高了矿震远场位移的精准度；

3、本煤岩冲击失稳模式的微震识别方法中，破裂面产状求取方法,针对工程实践中存在大量非剪切破裂形式,使用破裂面法向矢量、运动矢量和矩张量特征矢量之间的关系，定量推导一般矩张量破裂面产状，提高了矩张量求解矿震破裂面产状方法的适用性。

4、本煤岩冲击失稳模式的微震识别方法中，通过绘制震源球，将地质构造分析的赤平极射投影方法与矩张量反演结果相结合，增强了反演结果的可读性。

附图说明

图1为煤岩冲击失稳模式的微震识别方法流程图。

图2为煤岩冲击失稳模式的微震识别方法震源矩张量分析示意图。

图3为煤岩冲击失稳模式的微震识别方法实例震源筛选波形图。

图4为煤岩冲击失稳模式的微震识别方法实例震源球示意图。

具体实施方式

图1～4是本发明的最佳实施例，下面结合附图1～4对本发明做进一步说明。

如图1所示，一种煤岩冲击失稳模式的微震识别方法，包括如下步骤：

步骤1001，布置微震记录仪。

在同一采场不同位置设置至少6个微震记录仪对采场进行微震监测，每个微震记录仪形成一个微震信号记录。

步骤1002，在微震记录仪监测到的微震记录中，选择6个相对清晰地微震波形。

在步骤1002中，清晰的地震波形指波形存在明显的主振，优选从距离采区2km以内的微震记录仪中进行微震波形的筛选。

步骤1003，计算震源远场位移。

在计算震源远场位移时，包括如下步骤：

步骤a1：需要根据步骤1002中筛选出的微震波形，拾取每一个微震波形的P波初振幅值，并记作U1。

步骤a2，在筛选的微震波形中剪切P波的时域波形；

步骤a3，根据微震记录仪固有采样频率将步骤a3中剪切的P波时域波形进行快速傅里叶变换，转化为频域波形；

步骤a4，对步骤a3中转化得到的频域波形利用如下公式(1)进行衰减修正：

$A^{new}\left(f\right)=A\left(f\right)\exp \left(\frac{\pi fD}{vQ}\right)---\left(1\right)$

其中，A(f)为时域速度谱进行快速傅里叶变换的结果；f为步骤a3中转换得到的频域波形的频率；v为P波速度；Q为衰减因子。

在拾取微震波形的P波初振幅值后，计算微震记录仪记录的微震波形的远场位移。由于不同位置微震记录仪距离噪源远近不同，受到噪声的影响也不尽相同。低频位移振幅计算无需任何源模型假设，并且双对数频率域的振幅值在低频部分是相对平直的常数，其大小正比于地震标量矩。因此采用低频位移振幅来表示岩石破裂的远场位移。

步骤a5，在步骤a4中利用公式(1)进行衰减修正后得到地动速度功率谱与相应频率的关系：并根据地动速度功率谱得出地动位移功率谱

$D_0^2\left(f\right)=\frac{v_0^2\left(f\right)}{{\left(2\pi f\right)}^2}---\left(2\right)$

步骤a6，根据微震波形的地动速度功率谱和地动位移功率谱得出地动速度功率谱积分S_V2和地动位移功率谱积分S_D2：

$S_{V2}=2{\int }_0^{\infty }{V}_0^2\left(f\right)df---\left(3\right)$

$S_{D2}=2{\int }_0^{\infty }{D}_0^2\left(f\right)df---\left(4\right)$

其中：为考虑自由面影响乘以1/4的修正速度功率谱；为对应的位移功率谱。

步骤a7，由上述公式(2)～(4)得出各个微震波形的低频位移U2：

$U_2=\sqrt{\frac{4S_{D2}^{3/2}}{S_{V2}^{1/2}}}---\left(5\right)$

利用上述公式计算出的微震记录仪记录的微震波形低频位移U₂，可近似表示为该波形对应的P波远场位移。

利用上述公式计算出P波初振幅值U1和微震波形的低频位移U2后，若U1>0，且微震记录仪在震源Z轴上方，则U2值取正；若U1>0，且微震记录仪在震源Z轴下方，则U2值取负；若U1<0，且微震记录仪在震源Z轴上方，则U2值取负；若U1<0，且微震记录仪在震源Z轴下方，则U2值取正。

步骤1004，计算震源矩张量。

对于单向微震记录仪，由于S波受噪声和P波尾波影响较为严重，一般不用S波作矩张量计算，而P波远场位移相关的激励矩阵可以表示为：

F＝γ₃γ_iγ_j(6)

式中γ_i为震源至微震记录仪的震动波射线对应于各坐标轴的分量，即γ_i＝(x_i-x_0i)/r(x_i为微震记录仪各坐标分量，x_0i为震源各坐标分量，i＝1,2,3)。

根据弹性波理论推导可以得出P波的远场位移：

$u_{p,k}=\frac{{\gamma }_k{\gamma }_i{\gamma }_j}{4{{\mathrm{\pi \rho v}}_p}^{2}r}{\stackrel{·}{M}}_{ij}(t-\frac{r}{v_p})---\left(7\right)$

v_p为P波传播速度；r为震源到微震记录仪的距离；ρ为岩石密度；k为微震记录仪的第k(k＝1,2,3)分量；γ_i为震源至微震记录仪的震动波射线对应于各坐标轴的分量，即γ_i＝(x_i-x_0i)/r(x_i为微震记录仪各坐标分量，x_0i为震源各坐标分量，i＝1,2,3)；M_ij为作用于震源的矩张量。

对于单分量微震记录仪无需进行偏振处理，则上式可以表示为：

$u_p=\frac{{\gamma }_3{\gamma }_i{\gamma }_j}{4{{\mathrm{\pi \rho v}}_p}^{2}r}{\stackrel{·}{M}}_{ij}(t-\frac{r}{v_p})---\left(8\right)$

其矩阵表示为：

$\left(\begin{array}{c}{u}_p^1\\ u_p^2\\ u_p^3\\ ·\\ ·\\ ·\\ ·\\ ·\\ ·\\ u_p^n\end{array}\right)=\frac{1}{4{{\mathrm{\pi \rho rv}}_p}^3}\left(\begin{array}{cccccc}{\gamma }_3^1{\gamma }_1^1{\gamma }_1^1& 2{\gamma }_3^1{\gamma }_1^1{\gamma }_2^1& 2{\gamma }_3^1{\gamma }_1^1{\gamma }_3^1& {\gamma }_3^1{\gamma }_2^1{\gamma }_2^1& 2{\gamma }_3^1{\gamma }_2^1{\gamma }_3^1& {\gamma }_3^1{\gamma }_3^1{\gamma }_3^1\\ {\gamma }_3^2{\gamma }_1^2{\gamma }_1^2& {\gamma }_3^2{\gamma }_1^2{\gamma }_2^2& 2{\gamma }_3^2{\gamma }_1^2{\gamma }_3^2& {\gamma }_3^2{\gamma }_2^2{\gamma }_2^2& 2{\gamma }_3^2{\gamma }_2^2{\gamma }_3^2& {\gamma }_3^2{\gamma }_3^2{\gamma }_3^2\\ {\gamma }_3^3{\gamma }_1^3{\gamma }_1^3& 2{\gamma }_3^3{\gamma }_1^3{\gamma }_2^3& 2{\gamma }_3^3{\gamma }_1^3{\gamma }_3^3& {\gamma }_3^3{\gamma }_2^3{\gamma }_2^3& 2{\gamma }_3^3{\gamma }_2^3{\gamma }_3^3& {\gamma }_3^3{\gamma }_3^3{\gamma }_3^3\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ {\gamma }_3^n{\gamma }_1^n{\gamma }_1^n& 2{\gamma }_3^n{\gamma }_1^n{\gamma }_2^n& 2{\gamma }_3^n{\gamma }_1^n{\gamma }_3^n& {\gamma }_3^n{\gamma }_2^n{\gamma }_2^n& 2{\gamma }_3^n{\gamma }_2^n{\gamma }_3^n& {\gamma }_3^n{\gamma }_3^n{\gamma }_3^n\end{array}\right)\left(\begin{array}{c}{M}_{11}\\ M_{12}\\ M_{13}\\ M_{22}\\ M_{23}\\ M_{33}\end{array}\right)---\left(9\right)$

其中，γ上标表示通道号，下表表示坐标分量。

震源矩张量可以用公式(10)表示为：

$\left(\begin{array}{c}{M}_{11}\\ M_{12}\\ M_{13}\\ M_{22}\\ M_{23}\\ M_{33}\end{array}\right)=4{{\mathrm{\pi \rho rv}}_p}^3{\left(\begin{array}{cccccc}{\gamma }_3^1{\gamma }_1^1{\gamma }_1^1& 2{\gamma }_3^1{\gamma }_1^1{\gamma }_2^1& 2{\gamma }_3^1{\gamma }_1^1{\gamma }_3^1& {\gamma }_3^1{\gamma }_2^1{\gamma }_2^1& 2{\gamma }_3^1{\gamma }_2^1{\gamma }_3^1& {\gamma }_3^1{\gamma }_3^1{\gamma }_3^1\\ {\gamma }_3^2{\gamma }_1^2{\gamma }_1^2& {\gamma }_3^2{\gamma }_1^2{\gamma }_2^2& 2{\gamma }_3^2{\gamma }_1^2{\gamma }_3^2& {\gamma }_3^2{\gamma }_2^2{\gamma }_2^2& 2{\gamma }_3^2{\gamma }_2^2{\gamma }_3^2& {\gamma }_3^2{\gamma }_3^2{\gamma }_3^2\\ {\gamma }_3^3{\gamma }_1^3{\gamma }_1^3& 2{\gamma }_3^3{\gamma }_1^3{\gamma }_2^3& 2{\gamma }_3^3{\gamma }_1^3{\gamma }_3^3& {\gamma }_3^3{\gamma }_2^3{\gamma }_2^3& 2{\gamma }_3^3{\gamma }_2^3{\gamma }_3^3& {\gamma }_3^3{\gamma }_3^3{\gamma }_3^3\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ {\gamma }_3^n{\gamma }_1^n{\gamma }_1^n& 2{\gamma }_3^n{\gamma }_1^n{\gamma }_2^n& 2{\gamma }_3^n{\gamma }_1^n{\gamma }_3^n& {\gamma }_3^n{\gamma }_2^n{\gamma }_2^n& 2{\gamma }_3^n{\gamma }_2^n{\gamma }_3^n& {\gamma }_3^n{\gamma }_3^n{\gamma }_3^n\end{array}\right)}^{-1}\left(\begin{array}{c}{u}_p^1\\ u_p^2\\ u_p^3\\ ·\\ ·\\ ·\\ ·\\ ·\\ ·\\ u_p^n\end{array}\right)---\left(10\right)$

步骤1005，震源矩张量的分解。

在主坐标轴(a₁，a₂，a₃)中，公式(10)中的矩张量矩阵可对角化为：

$M=\left(\begin{array}{ccc}{M}_{11}& M_{12}& M_{13}\\ M_{21}& M_{22}& M_{23}\\ M_{31}& M_{32}& M_{33}\end{array}\right)=\left(\begin{array}{ccc}{M}_1& 0& 0\\ 0& M_2& 0\\ 0& 0& M_3\end{array}\right)---\left(11\right)$

a₁，a₂，a₃为矩阵M的本征矢量，M₁，M₂和M₃为对应的特征值，上式可分解为：

$M=\left(\begin{array}{ccc}P& 0& 0\\ 0& P& 0\\ 0& 0& P\end{array}\right)+\left(\begin{array}{ccc}{M_1}^{\prime }& 0& 0\\ 0& {M_2}^{\prime }& 0\\ 0& 0& {M_3}^{\prime }\end{array}\right)=PI+M^{\prime }---\left(12\right)$

P＝(M₁+M₂+M₃)/3，I为单位矩阵，PI表示矩张量的各向同性部分，M′表示矩张量的偏张量部分，本征值为M_i'＝M-P，i＝1,2,3。

就矩张量分解其实质而言，矩张量分解不是一个数学问题，而是一个物理问题。对于矿井采动诱发微震的震源破裂机理研究来讲，纯双力偶(DC)与补偿线性矢量偶极(CLVD)的分解方式是最为合理。即：矩张量可分解为以下三个部分：一个各向同性部分、一个补偿线性矢量偶极子部分及一个双力偶部分。

$\begin{array}{c}M=PI-{M_1}^{,}\left(\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 2\end{array}\right)+\left({M_3}^{,}+2{M_1}^{,}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& -1& 0\\ 0& 0& 1\end{array}\right)\\ =PI-{M_1}^{,}\left(2a_3{a}_3-a_2{a}_2-a_1{a}_1\right)+\left({M_3}^{,}+2{M_1}^{,}\right)\left(a_3{a}_3-a_2{a}_2\right)\\ =PI+{M_3}^{,}F\left(2a_3{a}_3-a_2{a}_2-a_1{a}_1\right)+{M_3}^{,}\left(1-2F\right)\left(a_3{a}_3-a_2{a}_2\right)\end{array}---\left(13\right)$

式中：F＝-M₁'/M₃'，且0≤F≤1/2。

其中：PI为各项同性部分，M₃'F(2a₃a₃-a₂a₂-a₁a₁)为补偿线型偶极部分，M₃'(1-2F)(a₃a₃-a₂a₂)为双力偶部分。

步骤1006，震源矩张量的分析。

在本煤岩冲击失稳模式的微震识别方法中，包括对震源矩张量以下几个方面的分析：震源破裂模式的分析，震源P、T轴的分析以及震源破裂面产状的分析。

在震源破裂模式的分析中，纯双力偶模型的解的偏差可以F＝-M₁'/M₃'来计算，其中分子分母分别为最小和最大的偏本征值。在本煤岩冲击失稳模式的微震识别方法中，使用DC％＝Mdc/(Mdc+Mclvd+Miso)来计算得到矩张量中剪切破裂的比重，并根据剪切破裂部分所占矩张量的比重来进行破裂类型判别，具体计算方法如下：

DC％＝(1-2*F)*(1-norm(Miso)/norm(M))

式中norm为求解矩阵最大奇异值函数。

震源的破裂类型可以根据下式进行判别：

DC％≥60％，剪切破裂

DC％≤40％，拉张破裂

40％＜DC％＜60％，混合破裂

在震源P、T轴的分析中，将震源矩张量M分解后得到M′(偏张量部分)，在主坐标轴(b₁，b₂，b₃)中，震源矩张量矩阵可对角化为：

$M^{\prime }=\left(\begin{array}{ccc}{M_{11}}^{\prime }& {M_{12}}^{\prime }& {M_{13}}^{\prime }\\ {M_{21}}^{\prime }& {M_{22}}^{\prime }& {M_{23}}^{\prime }\\ {M_{31}}^{\prime }& {M_{32}}^{\prime }& {M_{33}}^{\prime }\end{array}\right)=\left(\begin{array}{ccc}{M_1}^{\prime }& 0& 0\\ 0& {M_2}^{\prime }& 0\\ 0& 0& {M_3}^{\prime }\end{array}\right)---\left(14\right)$

b₁，b₂，b₃为矩阵M′的本征矢量，M₁′，M₂′和M₃′为对应的特征值。假设有M₁′<M₂′<M₃′，则b₁对应拉张方向轴T，b₃对应压缩方向轴P。

在震源破裂面产状的分析中，各项同性介质下震源矩张量可以由破裂面的位置与运动矢量表示为：

式中：i，j，k＝1，2，3或者x，y，z，为空间坐标系中的3个方向；Mij为地震矩张量形式，根据i，j取值不同，分别代表矩张量不同分量；u为破裂面运动方向的位移量；S为破裂面表面积；λ和μ为拉梅常数；v为运动方向，割据i,j,k取值不同，vi,vj,vk表示破裂面运动矢量在各坐标系方向上的分量；n为破裂面法向方向，其分量取值意义通v。

将上公式进行本征值化，可以得到：

$M=uS\left(\begin{array}{ccc}(\lambda +\mu )nv+\mu & 0& 0\\ 0& \lambda nv& 0\\ 0& 0& (\lambda +\mu )nv-\mu \end{array}\right)---\left(16\right)$

则有(假设M1>M2>M3)：

由于矩张量的对称性，破裂面滑动矢量v和破裂面法向矢量n互易带来矩张量的结果一致的，且根据矩张量特征值的大小关系，可以得到如下特征矢量与破裂面运动方向的关系。

式中，e₁⊥e₂⊥e₃，绝对值符号表示矢量大小；×表示矢量乘法；e₁、e₂、e₃分别为矩张量的最大特征值、中间特征值和最小特征值。

假设矢量v和矢量n的夹角为β，v和e1的夹角以及n和e1的夹角均为β/2，则有：

$cos\frac{\beta }{2}=\sqrt{\frac{M_1-M_2}{M_1-M_3}}---\left(21\right)$

$sin\frac{\beta }{2}=\sqrt{\frac{M_2-M_3}{M_1-M_3}}---\left(22\right)$

于是，可以得到破裂面运动方向和法向方向与矩张量最大特征值对应特征矢量和最小特征值对应特征矢量的关系：

$n=cos\frac{\beta }{2}{e}_1+sin\frac{\beta }{2}{e}_3=\sqrt{\frac{M_1-M_2}{M_1-M_3}}{e}_1+\sqrt{\frac{M_2-M_3}{M_1-M_3}}{e}_3---\left(23\right)$

$v=cos\frac{\beta }{2}{e}_1-sin\frac{\beta }{2}{e}_3=\sqrt{\frac{M_1-M_2}{M_1-M_3}}{e}_1-\sqrt{\frac{M_2-M_3}{M_1-M_3}}{e}_3---\left(24\right)$

根据破裂面法向方向的空间矢量值就可以得到破裂面的几何方程表达式，可以确定破裂面的方位角和倾角，具体计算公式如下：其中方位角公式为：

$dip=acos\left(\frac{n\left(3\right)}{\left|n\right|}\right)---\left(25\right)$

式中，n(3)为破裂面法向矢量Z轴分量，|n|表示法向向量的模

若dip>90°，则破裂面方位角为(180°-dip)

倾角公式为：

$strike=\left(\begin{array}{c}pi-atan\left(\frac{n\left(2\right)}{n\left(1\right)}\right),n\left(1\right)>0\\ -atan\left(\frac{n\left(2\right)}{n\left(1\right)}\right),n\left(1\right)<0\end{array}\right)---\left(26\right)$

式中，n(1)为破裂面法向向量X轴分量，n(2)为破裂面法向向量Y轴分量。

步骤1007，绘制震源球。

根据P、T轴方向及破裂面产状绘制震源球。

下面结合某煤矿某工作面的一起微震事件实例，对本煤岩冲击失稳模式的微震识别方法进行进一步说明，在该次微震事件中，震源定位于X＝4614m，Y＝1726m，Z＝-528m，微震能量为2.1E+06J。

利用上述步骤1001中布置的微震记录仪进行微震的记录，微震记录仪的布置坐标如表1所示：

表1微震记录仪坐标

从上述12处微震记录仪中，筛选出其中3、5、8、9、10、11六个通道波形，其波形如图3所示。

根据步骤1003得出，P波初振幅值U1为：

U1＝[-2.57e-5 4.56e32 5.50e-5 1.72e-4 -6.11e-5 7.67e-6]^T；

根据公式(2)～(5)计算得到微震波形的低频位移U2为：

U2＝[4.70e-5 6.64e-5 -5.02e-5 6.61e-5 1.185e-04 5.16e-5]^T。

根据上述步骤1004，可得出震源矩张量为：

$M=\left(\begin{array}{ccc}1.997e15& -7.260e15& 4.566e16\\ -7.260e15& -1.355e16& -9.367e16\\ 4.566e16& -9.367e16& 3.067e17\end{array}\right)$

根据上述步骤1005，对震源矩张量进行分解得到：

矩张量各项同性部分：

$Miso=1.0e16·\left(\begin{array}{ccc}9.838& 0& 0\\ 0& 9.838& 0\\ 0& 0& 9.838\end{array}\right)$

矩张量纯剪切部分：

$Mdc=\left(\begin{array}{ccc}-6.914e12& 3.409e15& 6.015e15\\ 3.409e15& -3.022e16& -1.830e16\\ 6.015e15& -1.830e16& 3.023e16\end{array}\right)$

矩张量偏张量部分：

$Mclvd=1.0e16·\left(\begin{array}{ccc}-9.637& -1.067& 3.965\\ -1.067& -8.171& -7.536\\ 3.965& -7.536& 17.808\end{array}\right)$

其中Mdc部分，即DC％为10.64％，小于40％，因此该震源以拉张破裂为主。

根据上述步骤1006对震源矩张量进行分析：

将震源矩张量偏张量部分(dev)对角化可以求得矩张量本征矢量和特征值，从而得到拉张方向轴T和压缩方向轴P：

T＝[-0.136 0.952 -0.275]^T

P＝[-0.135 0.257 -0.957]^T

震源破裂面参数：

破裂面法向向量：n＝(0.0868 -0.0488 0.9950)^T

破裂面滑动向量：v＝[0.1709 -0.5386 0.8251]^T

基于破裂面法向量和破裂滑动向量可求得破裂面走向为258.91°，破裂面倾角为5.51°。

(8)绘制震源球：

基于以上矩张量分解信息绘制震源球，如图4所示。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。