在煤岩体中存在各种微裂隙、孔隙，在外力作用时会在这些缺陷部位产生应力集中，发生突发性破裂，使积聚在煤岩体中的能量得以释放，并以弹性波的形式向外传播. 在极短的时间内，由煤体向巷道或采场喷出大量的瓦斯流及碎煤，形成动力效应，严重地威胁着煤矿的安全生产.

微震监测技术近年来被广泛用来研究冲击地压活动. 利用煤岩体破裂所产生的振动信号来研究煤岩体破碎的过程，对煤岩体破裂的时间、地点、强度等特征实时监测，使该技术成为研究煤与瓦斯突出和冲击地压的有效手段^[1-4]. 与其他技术相比，微震监测具有距离远、动态、三维、实时监测的特点. 目前，国内开展的基于微震监测的地质灾害预警研究主要集中在案例统计、动力灾害预测和灾害防控技术等方面，结合微震监测与地震波波速层析成像，是应力异常监测和冲击地压防治研究的最新方向^[5-10]. 窦林名等^[11]在应力与纵波波速的试验关系模型基础上，建立冲击危险性动态预测评价的震动波 CT 探测的技术，构建标波速异常系数和波速梯度变化系数等指标，给出指标的判别准则. 夏永学等^[12]采用微震和地音监测系统，对千秋煤矿煤岩震动信息进行联合监测，建立冲击危险性联合评价模型. 姜福兴等^[13]采用微震监测技术和应力实时在线监测技术，研究综采工作面煤岩破坏、应力变化和瓦斯涌出三者之间的关系. 王书文等^[14]根据波速异常率和波速分区特征，研究采煤工作面内断层展布的探测和工作面超前支承压力分布的划定，提出冲击危险性评价模型. 将微震波速层析成像技术用于探测地压冲击或矿震危险分布是可行的，也是未来国内外煤矿地压灾害预测的一个重要手段和研究方向.

结合震动层析成像和微震监测技术对矿区冲击地压危险性进行评估的研究，较好地解决了冲击危险区域的探测评价问题，但受到煤矿地质结构复杂、反演耗时等因素的影响，时效性较差，在应用中受到一定的限制. 本文针对煤矿层状煤系地层的典型特征，对工作面微震活动进行连续监测. 采用子空间分段法进行走时层析成像，研究监测区域内地震波速度的三维空间动态分布及变化，为采煤工作面冲击地压的危险性判断提供一种手段.

1 煤系地层地震波传播与应力场

微震在矿山冲击地压灾害的应用上，具有独特的特征. 对矿山来说，监测的范围小，尤其是在生产矿井，单一微震的发生具有极大的随机性，一个微震不能产生冲击地压，即使模型能够准确预报下一微震发生的时间地点，也不能确定冲击地压能否发生. 地震波传播速度是一个区域的特征，符合冲击地压高危区的空间分布描述，高波速区与应力集中区对应一致性好^[15]. 利用微震成像技术，根据监测区域内的地质力学特征，建立三维动态模型；根据地震波速度的变化，分析工作面应力场的变化，监测应力异常区，对冲击地压动力灾害的预防研究具有积极的意义.

1.1 地震波传播速度与应力场

地震波在岩石中的传播速度受到岩石应力的影响，对于均质、各向同性连续介质体，矿震震动波的传播与煤岩体物理力学参数及其在煤岩体中产生的动载荷^[16]可以表示为

$\left.{{{v_{\rm p}}}}\right/{{{v_{\rm s}}}} = \sqrt {\left.{{2(1 - \mu )}}\right/{({1 - 2\mu })}} ;$

(1)

$E ={{\rho v_{\rm s}^2\left(3v_{\rm p}^2 - 4v_{\rm s}^2\right)}}\left/\left[{{2\left(v_{\rm p}^2 - v_{\rm s}^2\right)}}\right]\right.;$

(2)

$ \begin{gathered} {\sigma _{{\rm dp}}} = \rho {v_{\rm p}}{{\rm{(}}{v_{{\rm pp}}}{\rm{)}}_{\rm p}} , {\tau _{{\rm ds}}} = \rho {v_{\rm s}}{{\rm{(}}{v_{{\rm pp}}}{\rm{)}}_{\rm p}} . \end{gathered} $

(3)

式中： ${v_{\rm p}}$ 、 ${v_{\rm s}}$ 分别为纵波和横波的速度；μ为泊松比， $0 \leqslant \mu \leqslant 0.5$ ；E为弹性模量；ρ为介质密度； ${\sigma _{{\rm dp}}}$ 、 ${\tau _{{\rm ds}}}$ 分别为纵波和横波产生的动载； ${{\rm{(}}{v_{{\rm pp}}}{\rm{)}}_{\rm p}}$ 、 ${{\rm{(}}{v_{{\rm pp}}}{\rm{)}}_{\rm s}}$ 分别为纵波、横波传播引起的质点峰值振动速度.

巩思园等^[17]研究不同受载方式下应力-应变与纵波波速的关系. 在单轴压缩和三轴压缩条件下，在弹性阶段，对煤岩进行压密时，应力使细微裂隙闭合，促使颗粒间的接触紧密，从而导致纵波波速增加. 随着压力的增加，缝隙密度相比初始阶段变化较小，波速增加的幅度降低. 在进入塑性区后，试样发生损伤，主要表现为颗粒间的位错和滑移，纵波波速的增加微弱，应力与纵波波速间具有幂函数关系；在三轴压缩条件下，平行于加载方向产生的微观裂纹对纵波波速的传播影响微弱，垂直于加载方向产生的微观裂隙使得纵波波速迅速降低. 根据岩石全应力-应变曲线可知，在岩石破坏阶段，内部形成宏观裂隙，纵波速度下降. 根据地震波纵波速度的变化，可以判断岩石所受应力，分析监测区域的应力异常.

不同的岩石中地震波速度都有一定的范围. 在淮南地区煤系地层中，煤层中的纵波速度为2 000~2 500 m/s，瓦斯突出煤层中的纵波波速通常小于900 m/s，页岩中的纵波波速为2 000~4 500 m/s,砂岩中的纵波波速为2 500~5 000 m/s，石灰石中的纵波波速为4 000~6 000 m/s.

1.2 地震波在煤系地层中的传播

微震成像技术不同于垂直地震剖面(vertical seismic profiling)或者地震勘探，它不使用人工震源，而是利用矿山开采过程中发生的微震作为震源，计算地震波传播速度. 由于煤系地层具有典型的层状特点，而且层与层之间具有明显的界面，各层之间的速度差别较大. 另外，微震的产生具有随机的特征，位置是不确定的，因此在使用初始到达时间进行地震波速度模拟时，需要考虑从震源到传感器之间的所有传播路径以及沿这些路径传播所需的时间.

如图1所示为层状煤系地层地震波传播路径示意图，每层的介质是均匀的，介质4为介质3提供反射系数作用. 各介质的参数如表1所示. 在这4条路径中，地震波沿哪条路径首先到达传感器是未知的. 当纵波沿路径1到达传感器时，表明介质1在震源与传感器之间无裂隙、构造煤发育区和流体渗入区等低速异常区. 当纵波沿路径2首先到达传感器时，表明在震源和传感器之间存在低速异常区，纵波经介质1和介质2之间的界面反射后首先被传感器接收. 当纵波沿路径3首先到达传感器时，表明在地震波反射线经过的区域中，介质2中纵波波速大于介质1中的纵波波速. 路径4是一种煤层中特有的传播方法，在地震波进入煤层后，地震波的主要能量会在煤层的顶底板之间反复多次反射，在煤层中传播较远距离，仍能被传感器接收，这是一种极端的情况，表明震源和传感器之间的岩层存在较大的垮落带，煤层存在应力异常，有煤层突出的危险. 在实际应用中，将地质勘探的结果作为模型的初始值，在每次模拟计算中，震源到传感器之间的所有路径都须进行计算和比较；根据检波器的初至时间，选择地震波传播的最优路径. 在反演计算过程中，每改变一次速度模型，重新计算震源位置，并修改传播射线的起点，重新选择路径.

2 地震波走时层析成像反演方法 2.1 目标函数的建立

子空间分段法属于迭代的弱非线性反演方法，实现方式为在模型的若干个子空间，同时沿着多个方向去寻找目标函数的最小值^[18]. 该次反演采用子空间分段法和Frechet偏导结合的方法. 反演的目标函数可以表达为

${{S({ m})}} = {\rm{[}}{{\varPsi ({ m})}} + \varepsilon {{\varPhi ({ m})}} + \eta {{\varOmega ({ m})}}{\rm{]}}/2.$

(4)

式中：ε为衰减系数，η为平滑系数，m为速度模型矢量. Ψ(m)可以定义为

${{\varPsi ({ m})}} = {{\rm{[}}{{{ g}}} - {{{{ D}}}_{{\rm obs}}}{\rm{]}}^{\rm T}}{{{ C}}}_{\rm{d}}^{ - {\bf{1}}}{\rm{[}}{{{ g}}} - {{{{ D}}}_{{\rm obs}}}{\rm{]}}.$

(5)

式中：C_d为数据收敛矩阵，g为根据速度模型计算得到的传播时间矩阵，D_obs为传感器初始到达时间矩阵. 当Ψ(m)达到最小值时，代表了获取的速度模型最接近观察数据所反映的真实模型. 若假设数据的误差不相互关联，则 ${{{{C}}}_{{\rm{d}}i,j }}= {\delta _{ij}}{(\sigma _{\rm d}^j)^2}$ （i =1,2,···,N，j=1,2,···,K），其中N为反演迭代次数，K为数据个数，δ为数据局部变化值，σ_d为观测值与模拟值的差.

反演时，观测的数据不能决定所有的模型参数，存在模拟参数影响被低估或被其他参数掩盖的问题，为此，在目标函数加了2项限制因子. 第1项限制因子Φ (m)定义为

${{\varPhi ({ m})}} = {{{({ m}}} - {{{{ m}}}_{\bf{0}}}{\bf{)}}^{\rm T}}{{{ { C}}}}_{\rm m}^{ - 1}{{({ m}}} - {{{ m}}_{\bf{0}}}{\bf{)}}.$

(6)

式中：C_m为局部收敛矩阵. 假设初始模型的误差是随机的，则 ${{{{C}}}_{\rm m}}_i,_j = {\delta _{ij}}{(\sigma _m^j)^2}$ ，C_m中的数值取自反演过程中的阶段成果；Φ (m) 的作用是加速反演结果m更趋近输入模型m₀.

另一项限制因子Ω(m)是关于模型结构最小变化(光滑约束)的正则化约束项，定义为

${{\varOmega ({ m})}} = {{{{ m}}}^{\rm T}}{{{{ D}}}^{\rm T}}{{{{ D}}}_{\rm m}}.$

(7)

式中：D_m为中间模型的空间导数的差分近似解，D为二阶导数差分算子. 当模型中的数据比较均衡时，Ω(m)会快速减小.

对于ε和η的选择，从理论上说，应当选择数据拟合最佳时的ε和η，但ε和η的选取要考虑两者间的平衡. 当采用光滑约束时，反演算法的实现相对简单，但是数据拟合和模型扰动间的关系可能是强非线性的，此时ε和η的选择是不稳定的，相对较好的选择是根据图2所示ε和η的关于数据拟合度和模型扰动曲线进行权衡^[19]，通过一定次数的反演试算来选择最佳的ε和η.

2.2 目标函数的求解

求方程（4）的极小值时，采用迭代的解法：

$\Delta {{{ m}}} = {{{{ m}}}^n} - {{{{ m}}}^{n - 1}}.$

(8)

在每次迭代中，子空间分段法目标函数 S(m)的最小值搜索均限定在p维的子空间开展，模型扰动 $\Delta {{{ m}}}$ 可以用一组长度为M的p维基矢量 $ {{{{{ a}}}^j}} $ 表示：

$\Delta {{{ m}}} =\sum {_{j = 1}^p }\;{{\mu _j}{{{a}}^j} = {{{{ A}{{\mu}} }}}} .$

(9)

${{{ A}}} = {{{{{ a}}}^j}} $ ，分量 ${\mu _j}$ 决定了S(m)最小化时搜寻基矢量的长度. 对目标函数求取关于 μ 的偏导，并令偏导数为零，可得模型扰动的解^[20]：

$\Delta {{{ m}}} = - { {A}}{\left[ {{{ {A}}^{\rm T}}\left( {{{ {G}}^{\rm T}}{ {C}}_{\rm d}^{ - 1}{ {G}} + \varepsilon { {C}}_{\rm m}^{ - 1} + \eta {{ {D}}^{\rm T}}{ {D}}} \right){ {A}}} \right]^{ - 1}}{{{{ A}}}^{\rm T}}{\bf{\hat \gamma }}.$

(10)

式中： ${ {G}} = \partial { {g}}/\partial { {m}}$ ，为Frechet偏导数矩阵向量，

${{\hat {\gamma} }} = {{ {G}}^{\rm T}}{ {C}}_{\rm d}^{ - 1}({{{ g}}} - {{ {D}}_{{\rm obs}}}) + \varepsilon { {C}}_{\rm m}^{ - 1}({ {m}} - {{ {m}}_{\bf{0}}}) + \eta {{ {D}}^{\rm T}}{{{ D}_{\rm m}}}.$

(11)

利用子空间分段法时，计算模型更新量时求解的是一个较小的p×P维的线性方程组，因此计算的速度更快. 当选取基向量A时，可以通过分割模型空间的梯度向量获得. 为了避免基向量的相关性，需要对基向量进行Gram-Schmidt正交化.

Frechet偏导数表示走时关于模型参数的变化率，走时 ${t_j}$ 可以用射线 ${R_i}$ 传播路径的积分表示：

${t_j} = \sum\limits_{k = 1}^N {\int_{{l_{j,k}}} {\frac{{{\rm d}s}}{{{V_{{\rm c},k}}(x,y,z)}}} } .$

(12)

式中： ${l_{j,k}}$ 为 ${R_i}$ 穿过第k个单元时的长度， ${V_{c,k}} $ $(x,y,z)$ 为地震波在 ${l_{j,k}}$ 上的速度，N为射线穿过的单元个数. Frechet偏导数计算时，用射线段两端点导数的平均值来近似计算走时对速度的一阶导数^[21]：

$\begin{split} \frac{{\partial {t_j}}}{{\partial {v_k}}} \approx & - \sum\limits_k {\frac{{{l_{i,k}}}}{2}} \left\{ {{{\left[ {\frac{1}{{V_{{\rm c},k}^2(x,y,z)}}\left( {\frac{{\partial {V_{{\rm c},k}}(x,y,z)}}{{\partial {v_k}}}} \right)} \right]}_{{k_i}}}}+ \right. \\& \left. { {{\left[ {\frac{1}{{V_{{\rm c},k}^2(x,y,z)}}\left( {\frac{{\partial V_{{\rm c},k}^{}(x,y,z)}}{{\partial {v_k}}}} \right)} \right]}_{{k_j}}}} \right\} . \\ \end{split}$

(13)

式中： ${k_i}$ 、 ${k_j}$ 为射线的两个端点.

3 试验地质条件与传感器台网的布置 3.1 试验地质条件

谢一矿望峰岗井在-660 m水平以下，井田的含矿井主要可采煤层瓦斯含量最大值为8.40~17.85 m³/t，甲烷成分约为80%，同一煤层的瓦斯含量除了具有随深度增加而增高的趋势以外，还可能在局部形成瓦斯富集带. 在进入深部开采后，煤层瓦斯等级将进一步升级，属于“煤与瓦斯突出”矿井.

C5125工作面区域为2 000 m×200 m. 其中Ⅳ-Ⅴ线以北-660 m以上C15煤层未回采；下伏C13煤层-720 m 以上已回采，B10煤层正在回采；C13底巷正在掘进. 工作面内煤系地层总体为单斜构造，走向为NW方向，倾向为NE方向，平均倾角为20°，工作面内发育多条中、小型断裂构造.

煤地层为华北型石炭、二叠系，其中二叠系的山西组与上、下石盒子组为主要含煤层段，共含5组37层煤，可采煤层为17层，平均总厚为34.55 m；煤层结构以简单为主，夹矸多为泥岩、砂质泥岩和炭质泥岩，稳定性为稳定~较稳定型.

3.2 系统安装的台网布置

对于地震波走时层析成像来说，微震数量关系着地震波速度计算的精准性，强度较小的微震事件与强度较大的微震事件同样重要，因此系统选择较大的传感器矩阵及较高的传感器自然频率. 系统选取自然频率为15 Hz的动圈式检波器. 根据微震定位精度要求以及岩石地层的特点，选择30个传感器组成的传感器矩阵. 传感器安装时，在未破裂的岩石中钻孔，传感器与钻孔之间用树胶充填，树胶的声阻抗与岩石声阻抗应尽量相近，以防止地震波能量的损失.

传感器的布置设计是影响监测系统监测效果的重要因素. 为了有效地发挥监测系统的作用，传感器矩阵的布置基于以下3个原则：1）在三维空间上覆盖整个监测区；2）传感器尽可能均匀分布，使整个监测区的事件记录具有相同的效果；3）避免传感器布置在一个两维平面上，从而避免定位误差过大. 根据上述基本原则及C5125工作面的巷道分布情况，将传感器布置在C5125工作面的上、下煤巷、C15顶板运输巷3条运输巷中，传感器的分布处于立体空间中，受巷道走向和现场施工情况的影响，每个巷道中传感器的直线距离为100~200 m，如图3所示.

4 试验与分析 4.1 监测区域地震波速度异常分布与微震分布

选取4个具有代表性的三维反演截图进行分析，每个图片的时间间隔为10 d，监测区域地震波速度异常分布与微震事件位置分布的关系如图4所示.

在监测期间，中部区域接近开采工作面的岩层中，地震波速度较高，达到5 000 m/s；在工作面后方，存在岩石破碎和垮落，地震波速度较低，仅为1 000 m/s；在其他区域的岩石中，地震波速度处于正常状态，约为4 000 m/s. 在开采工作面向前推进的过程中，图4（b）中，在未开采煤层上方出现一个大面积的地震波高速区域；随着开采的开展，图4（c）、（d）中，该区域逐渐消失，仅在接近开采工作面的煤岩中存在地震波高速区.

在监测期间，微震活动主要集中在以下3个区域：C15下方的活动带（左侧）、工作面上方的活动带（中部）、C15上方的活动带（右侧），基本范围覆盖了整个工作区域. 根据工作现场的施工情况可知，中部及右部应有更多的微震活动，左段远离工作面及采空区，采动影响甚微，微震活动应较少；监测结果说明，除了开采工作区域（中部），左段巷道（上回风巷）周围的微震活动同样活跃，表明微震活动不总是与明显的波速异常（应力异常）有关，其他因素可以在没有明显的应力变化的情况下，产生微震活动. 微震活动的范围覆盖了整个区域，地震波速度异常带只限于C5125工作面周围.

4.2 应力异常区在煤层倾向剖面上的变化

图5给出4月5日至5月6日期间，煤层倾向剖面上的地震波传播速度变化，在此期间内，C5125工作面通过该切面. 如图5（a）所示为工作面通过之前的剖面，如图5（b）所示为工作面通过时的剖面，如图5（c）、（d）所示为工作面通过后的剖面. 在剖面图中，黑色为C15煤层.

在剖面方向上，C15煤的上方有一低速区域，速度小于1 000 m/s，起始于工作面的中部偏下，在整个监测期间始终存在，区域面积经历了由小变大再变小的过程，但相对位置固定. 结合开采过程中的施工情况可知，在工作面前方，煤层顶板可能开始破碎，裂隙增加，构成破坏带，造成地震波的传播速度大大降低，随着工作面的逐步接近，区域面积增加；随着工作面的远离，垮落区的破碎带受到顶板的下沉压实，地震波低速带缩小.

煤层顶板是应力异常的主要监测对象. 图5（b）中，在工作面通过时，顶板出现一个面积较大的高应力区域，位于顶板破碎带上方50 m左右，该区域的高度达到300 m. 此时的顶板岩层受到的应力高度集中，处于冲击地压灾害发生的危险时刻. 当工作面推进过后，地震波速度恢复到初始值附近，应力异常消失.

煤层底板的高应力区域与顶板不同，高应力区域延伸方向与煤层倾向方向相反，处于持续存在状态，但是区域面积动态变化. 当工作面推进时，高应力区域面积增加，随后面积降低；当工作面远离时，出现高应力异常区面积增加的状况，如图5（d）所示，可能的原因是由于采空区残留煤柱失稳，覆岩层受力不均，导致应力集中区域增大. 从应力异常区在煤层倾向剖面上的变化可以看出，利用微震成像技术可以有效地判断工作面的影响区域及工作面顶板与底板的异常应力分布状态随时间的变化.

地震速度剖面所显示的另一特点是煤层的存在对地震波速度的影响. 煤层与围岩的物理参数差异较大，在煤层的上、下，地震波速度可能出现较大的差异，甚至可以为相反的高速带和低速带.

4.3 地震波传播速度在空间上的变化

为了更全面地了解工作面应力异常在空间上的分布变化及影响范围，图6给出–800 m水平面上的微震成像监测结果. 与经线剖面相比，水平面地震波的速度变化相对较小.

从图6可以看出，在工作面前方，高应力异常不稳定；图6（a）中，它出现在工作面前沿，但在图6（b）中，应力异常区消失. 图6（c）中，高应力异常出现在工作面前方，同时伴随出现2个小的地震波低速带，显示该区域有煤岩破裂现象. 图6（d）中，高应力异常区呈现扩大趋势，表明应力集中区域向周围扩展，2个小的地震波低速带消失. 工作面后方的地震波低速带比较稳定且范围较大，它的区域远大于高速异常带，表明工作面后的垮落带一直存在，这与实际的工作现场情况一致.

5 结　论

(1)采用层状地质结构分析地震波传播路径，结合子空间分段法进行反演，是一种有效的、适合煤系地层的地震波走时层析成像的方法. 成像结果表明，其他人为因素与应力异常均可以诱发微震，地震波走时层析成像可以将微震诱因区别分类.

(2)C5125采煤工作面开采过程监测表明，工作面后方垮落区存在较稳定的低应力异常，工作面的推进对低应力异常区域面积的影响较小，工作面前方的高应力异常区动态变化较大；在工作面通过时，煤层顶板的高应力异常达到最大，煤层底板的高应力区持续存在.

(3)利用微震走时层析成像可以有效地判断工作面的影响区域以及工作面顶板与底板的应力分布状态随时间的变化. 在将来的研究中，可以与地震波全波反演技术联合，获得更精确的地质结构信息，为煤矿开采提供安全、高效、优质的生产环境.