煤矿微震数据分析学习笔记 · Class 1

以天然地震研究为起点，学习矿震定位与能量分析

目录

一、震源定位：三种方法的本质区别

1.0 核心问题与理论基础

地震（微震）发生后，需要通过连续波形记录反演以下四个未知量：

原始数据流：

连续波形采样 → 滤波去噪（带通/Bandpass） → STA/LTA 拾取 P 波到时 → 偏振分析拾取 S 波到时(深度学习直接拾取P/S) → 结合速度模型 → 定位反演

核心方程（走时方程）： 对于第 $i$ 个台站，P 波理论到时为：

$$ T_i^{\text{calc}} = t_0 + \int_{\text{source}}^{\text{station}_i} \frac{ds}{v(x,y,z)} $$

在均匀速度模型中简化为：

$$ T_i^{\text{calc}} = t_0 + \frac{D_i}{V_p} $$

其中 $D_i$ 为震源到第 $i$ 个台站的直线距离， $V_p$ 为 P 波速度。

定位的本质是最小化观测到时与理论到时之间的残差：

$$ \min_{x,y,z,t_0} \sum_{i=1}^{n} \left[ T_i^{\text{obs}} - T_i^{\text{calc}}(x,y,z,t_0) \right]^2 $$

1.1 P-S 时差法（Wadati 法）

原理推导

P 波和 S 波同时从震源出发。由于 $V_p > V_s$ ，到达同一台站的时间差为：

$$ \Delta t_{SP} = t_S - t_P = \frac{D}{V_s} - \frac{D}{V_p} = D \cdot \left( \frac{1}{V_s} - \frac{1}{V_p} \right) $$

从中解出震源距 $D$ ：

$$ D = \Delta t_{SP} \cdot \frac{V_p \cdot V_s}{V_p - V_s} $$

核心优势： 发震时刻 $t_0$ 在上述推导中完全消掉——两个到时相减时 $t_0$ 被抵消，这是此法不需要知道发震时刻的根本原因，而非巧合。

数值示例

$$\Delta t_{SP} = 2\ \text{s},\quad V_p = 5\ \text{km/s},\quad V_s = 3\ \text{km/s}$$

$$ D = 2 \times \frac{5 \times 3}{5 - 3} = 2 \times \frac{15}{2} = 15\ \text{km} $$

WADATI 图验证

Wadati（1933）提出利用 S-P 与 P 波到时的关系检验定位可靠性。由走时方程：

$$ t_S - t_P = \left( \frac{V_p}{V_s} - 1 \right) \cdot (t_P - t_0) $$

因此，对所有台站的 $(t_p, t_s - t_p)$ 数据点做直线拟合，斜率 $\kappa = V_p/V_s - 1$ ，由此验证速度模型是否与观测一致。如果数据点偏离直线，说明速度模型不当或 S 波拾取有误。

定位流程

  1. 每个台站：拾取 $\Delta t_{SP}$ → 计算震源距 $D$
  2. 以台站为球心、 $D$ 为半径画球面（三维）或圆周（二维）
  3. 最少 3 个台站 → 3 个球面 → 交点即震源位置

局限分析：

1.2 Geiger 法（最小二乘迭代定位）

核心思路

Geiger（1912）提出的定位方法至今仍是地震定位的标准算法。与 P-S 时差法不同，它将 $t_0$ 作为第四个未知数联合求解。

未知数向量： $\mathbf{m} = (x, y, z, t_0)^T$

最少台站数： 4 个（4 个未知数 × 4 个独立方程）

注： 4 个台站只是必要条件。当台站近似共面时，方程组条件数大，实际需要更多台站。更严格地说，需要台站从至少 4 个不同方向约束震源，即方位角覆盖 > 180° 且具有垂向分量差异。

数学框架

Geiger 法的核心是对非线性走时方程做一阶 Taylor 展开，转化为线性方程组迭代求解。

Step 1：Taylor 展开

将走时方程在初始猜测 $\mathbf{m}_0$ 处展开：

$$ T_i^{\text{obs}} - T_i^{\text{calc}}(\mathbf{m}_0) = \frac{\partial T_i}{\partial x}\Delta x + \frac{\partial T_i}{\partial y}\Delta y + \frac{\partial T_i}{\partial z}\Delta z + \frac{\partial T_i}{\partial t_0}\Delta t_0 $$

即：

$$ \mathbf{r} = \mathbf{G} \cdot \Delta\mathbf{m} $$

Step 2：偏导数计算

对于均匀模型 $T_i = t_0 + D_i / V_p$ ：

$$ \frac{\partial T_i}{\partial x} = \frac{x - x_i}{V_p \cdot D_i},\quad \frac{\partial T_i}{\partial t_0} = 1 $$

$y, z$ 分量同理。注意当 $D_i \to 0$ （台站无限接近震源）时偏导数趋于奇异，初始猜测不能离真值太远。

Step 3：最小二乘求解

$$ \Delta\mathbf{m} = (\mathbf{G}^T\mathbf{G})^{-1} \mathbf{G}^T \mathbf{r} $$

Step 4：更新并迭代

$$ \mathbf{m}_{k+1} = \mathbf{m}_k + \Delta\mathbf{m} $$

当 $|\Delta\mathbf{m}| < \varepsilon$ （典型值 1 m）或迭代次数超过上限时停止。

收敛条件与发散处理

$$ \Delta\mathbf{m} = (\mathbf{G}^T\mathbf{G} + \lambda \mathbf{I})^{-1} \mathbf{G}^T \mathbf{r} $$

对矿震场景的优势

不依赖 S 波——在以下情况至关重要：

1.3 Powell 法（直接搜索优化）

本质：无需梯度的非线性优化

与 Geiger 法（依赖梯度 $\partial T/\partial x$ ）不同，Powell 法属于直接搜索法（direct search），不需要计算导数。这一差异在地震定位中具有实际意义。

核心算法：共轭方向搜索

Powell 法的精髓在于沿一组共轭方向依次做一维搜索（线搜索）。对于二次函数，沿 $n$ 个共轭方向各搜索一次即可精确到达极小点。

算法流程（简化）：

$$ \min_{\alpha} \ \sum_i \left[ T_i^{\text{obs}} - T_i^{\text{calc}}(\mathbf{m} + \alpha \mathbf{d}_k) \right]^2 $$

与 Geiger 法的比较：

特性 Geiger 法 Powell 法
需要梯度 是（ $\partial T / \partial \mathbf{m}$ ） 否
收敛速度 二次收敛（近极小点时快） 超线性收敛（近极小点时慢于 Geiger）
对初值敏感度 高（Taylor 展开假设） 低（全局方向）
局部极小值 易陷入（受台站分布影响） 相对鲁棒（全局搜索特性）
适用情况 台站覆盖好的区域 台站稀疏或初值不确定时

矿震场景下的使用建议： 两种方法可结合使用——先以 Powell 法获得粗定位（鲁棒、不依赖梯度），再将结果作为 Geiger 法的初始值做精细定位（快速、可提供误差估计）。

1.4 方法对比总结

方法 需要波型 发震时刻 最少台站 定位方式 约束
P-S 时差法 P 波 + S 波 消掉，不需知 3 个 球面交会（解析） 需要可靠 S 波
Geiger 法 仅 P 波 作为未知数联合求解 4 个（理论下限） 迭代优化（梯度） 需要初值
Powell 法 任意（取决于目标函数） 作为未知数联合求解 4 个（理论下限） 直接搜索（无梯度） 速度较慢

二、矿震定位的核心挑战

2.1 台站几何约束

天然地震台网通常覆盖震源区上方及四周数百公里范围的区域，台站呈二维网状分布。矿震监测的根本制约在于：台站只能沿巷道布置。

巷道中的台站排列：

工作面方向 →
┌──────┬──────┬──────┬──────┐
│ 台1   │ 台2   │ 台3   │ 台4   │   ← 线性排列（近似一维）
└──────┴──────┴──────┴──────┘

三个球面的交点问题：当球心共线时，三个球面交于一个圆而非一个点，即定位方程在垂向几乎无约束，深度方向的定位误差可达水平方向的 3-10 倍。

2.2 数学本质：欠定系统

定位矩阵 $\mathbf{G}$ 的条件数 $\kappa(\mathbf{G})$ 反映了台站几何的质量：

$$ \kappa(\mathbf{G}) = \frac{\sigma_{\max}}{\sigma_{\min}} $$

其中 $\sigma_{\max}, \sigma_{\min}$ 分别是 $\mathbf{G}$ 的最大和最小奇异值。当台站线性排列时，最小奇异值对应的方向（垂向）几乎没有约束， $\kappa$ 极大，解极不稳定。

方位角覆盖（Azimuthal Gap）：

在水平面内，台站对震源的张角范围决定了水平定位能力。定义方位角间隙角 $\theta_{\text{gap}}$ 为相邻台站间的最大方位角差：

2.3 解决方向

台站布局受巷道固定，以下方向是实际可行的：

  1. 优化定位算法：引入正则化（阻尼最小二乘）、双差定位法（利用事件对之间的相对走时差提高相对精度）、波形互相关精化相对到时
  2. 精细速度模型：各向异性校正（煤层 vs 围岩速度差异可超过 50%）、层状模型取代均匀模型、利用掘进揭露的岩性信息标定参数
  3. 多参数联合约束：结合矿压监测的宏观应力信息辅助收敛、利用事件丛集性（相近事件应有相近定位）作为先验约束
  4. 波形反演：不依赖到时拾取，直接拟合全波形，在射线覆盖不足时利用散射波和尾波信息

三、微震能量计算

3.1 为什么用能量（焦耳）而不用震级

震级（Magnitude）和能量（Energy）是两个不同层面的指标：

这在实际监测中极为关键：当我们需要评价"过去 8 小时内累计释放了多少能量"或"该区域能量释放是否加速"时，必须用能量而非震级。

3.2 震级-能量关系

经典的 Gutenberg-Richter 能量-震级关系：

$$ \log_{10} E = 1.5 M + 4.8 $$

其中 $E$ 以焦耳计。该公式的物理基础：震源谱的角频率 $f_c$ 与应力降 $\Delta\sigma$ 和地震矩 $M_0$ 的关系（Brune 模型，1970）。

常用换算表：

震级 $M$ 能量 $E$ (J) 类比（TNT当量）
-2 0.06 一掌击桌
0 63 一拳冲墙
1 2,000 手榴弹
2 63,000 1 kg TNT
3 2,000,000 30 kg TNT

煤矿微震事件震级多在 -3 到 2 之间，能量在 $10^{-3}$ 到 $10^5$ J 量级，远小于天然地震。

3.3 能量如何从波形计算

理论推导

震源辐射地震波的能流密度（Poynting 矢量）为：

$$ S(t) = \rho v \cdot \dot{u}^2(t) $$

其中 $\rho$ 为介质密度， $v$ 为波速， $\dot{u}$ 为质点振动速度（即波形振幅的一阶时间导数）。

对台站记录的波场在时间上积分，并对球面几何做校正，得到震源辐射能量：

$$ E = 4\pi \rho v \cdot R^2 \cdot \int_{t_1}^{t_2} \dot{u}^2(t) , dt $$

其中 $R$ 为震源距。简化表达：

$$ E \propto \rho v R^2 \int A^2(t) , dt \quad \text{（振幅平方对时间积分）} $$

注： 上式隐含了各向同性辐射假设。实际辐射具有方向性（双力偶震源机制），需做辐射花样校正。在矿震场景中，由于事件多、台站少，通常不单独做辐射图案校正，而是假设各向同性，通过统计平均逼近。

实际计算步骤

$$ E_{\text{obs}} = \int_{t_{\text{start}}}^{t_{\text{end}}} v^2(t) , dt $$

矿震能量计算的误差来源

四、G-R 关系与 $b$ 值

4.1 公式与参数释义

Gutenberg-Richter（1944）关系描述地震的频度-震级分布：

$$ \log_{10} N(M) = a - b \cdot M $$

其中：

等效幂律形式：

$$ N(M) = 10^a \cdot 10^{-bM} = A \cdot 10^{-bM} $$

事件按震级的概率密度函数为指数分布：

$$ f(M) = b \ln(10) \cdot 10^{-b(M - M_{\min})}, \quad M \ge M_{\min} $$

4.2 $a$ 值分析

$a$ 值控制整条 G-R 线的高度。由 $\log_{10} N = a - bM$ 可知：

$$ N(M=0) = 10^a $$

4.3 $b$ 值的物理含义

$b$ 值控制大小地震的比例。设 $M_1 < M_2$ ：

$$ \frac{N(M_1)}{N(M_2)} = \frac{10^{a - bM_1}}{10^{a - bM_2}} = 10^{b(M_2 - M_1)} $$

特殊情况 $b=1$ 时，震级每增加 1 级，事件数减少 10 倍。

$b$ 值 特征 物理含义
$b < 0.7$ 低 大事件占比异常高 → 应力高度集中，接近破裂临界状态
$b \approx 1.0$ 正常 全球平均、区域常态，应力分布相对均匀
$b > 1.3$ 高 小事件占绝对主导 → 应力分散，裂隙发育或卸压状态

4.4 $b$ 值的估计方法

最小二乘估计：

$$ \hat{b} = \frac{ \sum (M_i - \bar{M})(\log_{10} N_i - \overline{\log_{10} N}) }{ \sum (M_i - \bar{M})^2 } $$

缺陷：大震级区间仅含少量事件但等权重，受尾部波动影响大。

最大似然估计（Aki, 1965）——推荐方法：

$$ \hat{b} = \frac{ \log_{10} e }{ \bar{M} - M_{\min} } \approx \frac{0.4343}{\bar{M} - M_{\min}} $$

其中 $\bar{M}$ 是样本震级均值， $M_{\min}$ 是完备震级阈值（所有 $\ge M_{\min}$ 的事件都被检测到）。

$M_{\min}$ （完备震级）的确定：

G-R 曲线在低震级端偏离线性——小事件漏检。确定 $M_{\min}$ 的方法包括：

  1. 最大曲率法：取频度-震级曲线曲率最大点
  2. 拟合优度法（GOF）：取使幂律拟合吻合度达到 90% 的最低震级
  3. $b$ 值稳定性法：增加 $M_{\min}$ ，取 $b$ 值开始稳定时的震级

4.5 $b$ 值变化的统计显著性

$b$ 值的变化需要与不确定性比较才能判断是否显著。Aki（1965）给出 $b$ 值的标准差估计：

$$ \sigma_b \approx \frac{b}{\sqrt{N}} $$

其中 $N$ 为用于估计的事件总数。当 $b$ 值变化 $|\Delta b| > 2\sigma_b$ 时才具有统计显著性。例如： $N=100$ 、 $b=1.0$ → $\sigma_b \approx 0.1$ ，需要 $|\Delta b| > 0.2$ 才显著。煤矿单日事件通常 $N<50$ ，统计显著性是一个重要制约。

五、 $b$ 值与矿震预警的连接

5.1 理论链条（Lab → Field）

实验室岩石力学实验（Scholz, 1968；Lockner, 1992）表明：

$$ \text{差应力} \uparrow \quad \rightarrow \quad b \text{ 值} \downarrow $$

物理机制：

5.2 矿震中的对应关系

将实验室结论映射到采矿场景：

采掘活动（推进/爆破/放顶）
  → 局部应力升高（掘进头前方、老顶悬露区）
  → 微破裂从随机转为有序（沿主应力方向排列）
  → b 值系统性下降
  → 同时伴随：Emax 升高、事件空间聚集、能量释放加速
  → 冲击地压/煤与瓦斯突出风险显著增加

5.3 沈玉恒 PPT 预测指标与 $b$ 值的对应

四项指标监测的是同一物理过程的不同侧面：

PPT 指标 定义 与 $b$ 值下降的等价关系
$E_{\max}$ 时间窗内最大单次能量 $b$ ↓ → 大事件变多 → $E_{\max}$ ↑
$E$ （累计能量） 时间窗内能量之和 $b$ ↓ → 每事件平均能量增大 → $\sum E$ ↑
$k$ （频次变化率） $\Delta N / \Delta t$ $b$ ↓ 时总频次不一定变，但需结合 $E_{\max}$ 判断
$f$ （分布集中度） 事件空间分布离散度 应力集中 → 破裂定向 → 空间聚集 → $f$ ↓（更集中）

核心结论： $b$ 值与四项 PPT 指标不是互相独立的，而是同一应力状态的不同表达方式。 $b$ 值提供了从频度-震级分布角度观测应力积累的视角，PPT 指标提供了从能量、空间、时间维度监测危险的视角。多维交叉验证可降低误报率。

5.4 预警中的注意事项

六、CT 层析成像

6.1 核心目标

以微震事件为被动震源（无需人工激发），反演工作面周围的速度分布 $v(x,y,z)$ ，进而识别应力集中区。

物理基础（岩石物理关系）：

Birch 定律（1961）给出速度-密度的经验关系：
$$ V_p = a(M) + b \cdot \rho $$
其中 $a(M)$ 与平均原子量有关， $b \approx 3.05\ \text{km} \cdot \text{s}^{-1} \cdot \text{g}^{-1} \cdot \text{cm}^3$ 。

判据总结：

速度特征 岩石状态 应力状态 冲击地压风险
高速（ $\Delta V_p / V_p > +3%$ ） 压实 应力集中 高
正常速度（ $ \Delta V_p / V_p < 3%$ ） 原岩
低速（ $\Delta V_p / V_p &lt; -3%$ ） 破碎 卸压/采空 低（但需关注顶板稳定性）

6.2 正问题：走时计算

给定速度模型 $v(\mathbf{r})$ ，震源 $\mathbf{s}$ 到台站 $\mathbf{r}$ 的走时由程函方程（Eikonal Equation）控制：

$$ |\nabla T(\mathbf{r})| = \frac{1}{v(\mathbf{r})} = s(\mathbf{r}) $$

其中 $s(\mathbf{r})$ 为慢度（slowness），是速度的倒数： $s = 1/v$ 。

沿射线路径积分：

$$ T = \int_{\mathbf{s}}^{\mathbf{r}} \frac{dl}{v(\mathbf{r})} = \int_{\mathbf{s}}^{\mathbf{r}} s(\mathbf{r}) , dl $$

非线性来源： 射线路径 $\mathbf{r}(l)$ 本身满足射线方程（Snell 定律的连续形式），而射线方程依赖于 $v(\mathbf{r})$ 。速度变了 → 射线路径改变 → 走时改变，这是典型的非线性正问题。

射线追踪方法概览：

方法 原理 优势 劣势
Shooting（打靶法） 从震源发射射线，调整出射角命中台站 直观 收敛慢，存在影区
Bending（弯曲法） 初始直线路径，逐步弯曲满足 Fermat 原理 对初值要求低 可能收敛到局部极小
Eikonal FD（程函方程有限差分） 直接求解程函方程，全场走时 鲁棒、无影区 仅能处理初至波
Fast Marching（快速行进法） 波前扩展法求解程函方程 高效 仅能处理初至波

矿震 CT 常用 Bending 法（配合伪弯曲法 pseudo-bending，Um & Thurber, 1987），兼顾效率和适应性。

6.3 反问题：矩阵方程的建立

网格离散化

将监测区域离散为 $N$ 个三维网格块（像素/体素 voxel），每个网格块具有待反演的慢度值 $s_j$ 。网格大小需同时满足：

走时方程离散化

第 $i$ 条射线（从第 $i$ 个微震事件到某台站）在离散模型中：

$$ T_i^{\text{obs}} = T_i^{\text{calc}} + \delta t_i = \sum_{j=1}^{N} L_{ij} \cdot s_j + \delta t_i $$

其中 $L_{ij}$ 是第 $i$ 条射线在第 $j$ 个网格中的路径长度。

将所有 $M$ 条走时方程写成矩阵形式：

$$ \delta\mathbf{t} = \mathbf{L} \cdot \delta\mathbf{s} $$

矩阵维度：

量 维度 含义
$\delta\mathbf{t}$ $M \times 1$ 观测走时与理论走时的残差（已知）
$\mathbf{L}$ $M \times N$ 射线路径矩阵（正向计算得到）
$\delta\mathbf{s}$ $N \times 1$ 各网格慢度扰动（待求）

迭代求解算法

由于 $\mathbf{L}$ 通常大且稀疏（ $M$ 几千到几万， $N$ 几百到几千），直接求伪逆 $(\mathbf{L}^T\mathbf{L})^{-1}\mathbf{L}^T$ 不可行。实际采用迭代类算法：

SIRT（Simultaneous Iterative Reconstruction Technique）：

$$ \delta s_j^{(k+1)} = \delta s_j^{(k)} + \frac{\eta}{\sum_i L_{ij}} \sum_i \frac{L_{ij} \cdot r_i^{(k)}}{\sum_j L_{ij}} $$

其中 $r_i^{(k)} = \delta t_i - \sum_j L_{ij} \delta s_j^{(k)}$ 为当前残差， $\eta$ 为松弛因子（0 < $\eta$ < 1，常用 0.1-0.5）。

SART（Simultaneous Algebraic Reconstruction Technique）： SIRT 的改进版，收敛更快，对欠定区域表现更好。矿震 CT 中 SIRT 最为常用。

6.4 鸡和蛋问题：定位与 CT 的联合反演

震源定位需要速度模型，CT 反演需要震源位置——两者相互依赖，构成经典反演难题。

解决方法——迭代联合反演：

初始化：假定均匀速度模型 v₀
    ↓
① 用 v₀ 定位所有微震事件 → 获得 (x, y, z, t₀) 目录
    ↓
② 对每个震源-台站对追踪射线 → 计算 L 矩阵和理论走时 T_calc
    ↓
③ 计算走时残差 δt = T_obs - T_calc
    ↓
④ 求解 δt = L · δs → 获得慢度扰动 → 更新速度模型 v_{k+1}
    ↓
⑤ 用新速度模型 v_{k+1} 重新定位事件
    ↓
⑥ 检查收敛：|v_{k+1} - v_k| < ε 或 RMS 走时残差不再下降
    ↓ No → 回到②
    ↓ Yes
输出：最终速度模型 + 精确定位目录

收敛判据：

6.5 欠定问题与正则化

欠定的根源

矿震台站沿巷道线性排列 → 射线主要沿巷道方向（水平）→ 垂直方向无台站接收 → 垂直方向速度分辨率为零。

数学上： $\mathbf{L}$ 矩阵列满秩但行数 $M \ll N$ （有效独立方程数远小于未知数），解空间为 $N - \text{rank}(\mathbf{L})$ 维的仿射空间，存在无穷多组解。

正则化的引入

通过添加先验信息约束将欠定问题转化为适定问题：

光滑性约束（Tikhonov 一阶/二阶导数）：

$$ \Phi_{\text{smooth}} = \iint |\nabla s|^2 , dV \quad \text{或} \quad \iint |\nabla^2 s|^2 , dV $$

离散形式：

$$ \Phi_{\text{smooth}} = \sum_{\text{相邻网格 } j,k} (s_j - s_k)^2 $$

总目标函数：

$$ \Phi = \underbrace{ |\delta\mathbf{t} - \mathbf{L} \delta\mathbf{s} |^2 }_{\text{数据拟合}} + \lambda^2 \underbrace{ |\mathbf{D} \delta\mathbf{s} |^2 }_{\text{光滑性约束}} $$

其中 $\mathbf{D}$ 为离散导数算子， $\lambda$ 为正则化参数。

正则化参数的选取（L-曲线法）：

以 $\log|\delta\mathbf{t} - \mathbf{L} \delta\mathbf{s}|$ vs $\log|\mathbf{D} \delta\mathbf{s}|$ 作图，曲线呈 L 形。最佳 $\lambda$ 取在 L 形拐角处：既不过拟合数据（ $\lambda$ 太小），也不过依赖先验约束（ $\lambda$ 太大）。

约束类型 正则项 过强时的副作用 深度学习类比
光滑性 $|\mathbf{D} \delta\mathbf{s}|^2$ 异常边缘被抹平，小异常消失 图拉普拉斯正则化
阻尼 $|\delta\mathbf{s}|^2$ 结果向 $\mathbf{s}_0$ 塌陷，异常被缩小 L2 权重衰减
稀疏性 $|\delta\mathbf{s}|_1$ 可能仅恢复最大异常 L1 (LASSO)

矿震 CT 的推荐策略： 先施加较强阻尼（稳定全局结构），再逐步减小阻尼、引入光滑约束（恢复细节），即从粗到细的联合迭代策略。

6.6 棋盘测试（Checkerboard Test）

目的： 量化反演结果在不同空间位置的可信度，而非验证模型是否正确。

流程：

  1. 构造人工棋盘模型：相邻网格块交替设定 ±5% 的速度扰动
  2. 用真实射线几何（与反演相同的 $\mathbf{L}$ 矩阵）做正演，计算合成走时 $T^{\text{syn}}$
  3. 加入符合实际噪声水平（如 0.5-1 ms）的随机噪声
  4. 用与真实数据完全相同的反演流程求解
  5. 比较恢复的棋盘图案与原始棋盘：恢复好 = 射线在该处覆盖充分 → 该处结果可信

矿震 CT 的典型分辨率分布：

方向 射线覆盖 分辨率 典型可分辨尺度
沿巷道（水平X） 密集 好 10-20 m
垂直巷道（水平Y） 稀疏 差 30-50 m
深度方向（Z） 几乎无 最差 50-100 m

实际解读策略： CT 结果对工作面正前方（射线主要穿过的区域）可信；对侧面围岩和深部（射线边缘或缺失方向的区域）需谨慎解读，宜结合其它监测手段（如钻孔应力计、矿压）交叉验证。

6.7 CT 在矿震监测中的完整逻辑链

微震事件自动拾取 & 初始定位（均匀模型）
    │
    ▼
┌──────────────────────────────────────┐
│  迭代联合反演（通常 3-5 轮）         │
│  ┌─ 定位事件                         │
│  ├─ 射线追踪（Bending / Eikonal）     │
│  ├─ 求解 δt = L · δs（SIRT + 正则化） │
│  ├─ 更新速度模型                      │
│  └─ 回到定位 → 重复                   │
└──────────────────────────────────────┘
    │
    ▼
最终速度模型 v(x,y,z) + 棋盘测试 → 确定可信区域
    │
    ▼
速度异常解读：
  高速区 = 岩石压实 = 应力集中 = ⚠ 冲击地压危险区
  低速区 = 岩石破碎 = 采空区/卸压区 = 相对安全
    │
    ▼
与其他监测手段联合 → 综合预警 & 针对性治理措施

七、微震与矿压的联合监测

7.1 两套系统的本质差异

维度 微震监测 矿压监测
物理原理 弹性波传播（被动声学） 力学量直接测量
测量对象 质点振动速度 → 反演震源参数 液压支架压力、锚杆轴力、围岩位移
空间覆盖 三维体（定位能力强时覆盖整个监测区域） 一维线（沿工作面一线）
提供信息 破裂位置、深度、能量、震源机制 支架载荷、顶板下沉量、帮部变形量
不能提供 绝对应力数值（仅相对速度变化） 破裂深度、破裂源具体位置
采样特性 连续触发（100-500 Hz） 定时采集（1-5 分钟一次）
数据量 极大（波形存储） 中等（时序数据）

7.2 为什么必须联合：四种状态的物理含义

仅靠单套系统存在盲区。联合后的四种状态对应不同的物理过程和风险：

① 静默期（微震静默 + 压力升高）→ 高危：能量积累

这一状态对应实验室岩石加载至峰值应力前的扩容硬化阶段：裂隙被压密，摩擦锁定抑制了小破裂的发生，但应变能仍在持续累积。此时的物理量：

一旦某点应力达到局部破裂强度，能量瞬间释放即为大能量冲击地压事件。这个阶段矿压监测是唯一能感知危险的手段。

② 微震静默 + 压力平稳 → 待观察

两种可能：

③ 微震活跃 + 压力下降 → 安全：稳态释放

典型的安全状态。老顶周期性断裂（周期来压）后：断裂面成形 → 支架压力下降；断裂过程中微震事件活跃 → 能量以众多小事件的形式平稳释放。类比：用许多小裂缝释放应力，比突然断裂安全。

④ 微震活跃 + 压力升高 → 高危：破裂加剧且支撑失效

最危险状态——破裂正在扩展（微震频繁），且支撑系统（支架+锚杆+锚索）无法控制围岩变形（压力持续升高）。可能对应：

联合判断矩阵：

微震状态 矿压状态 风险等级 物理含义 建议措施
静默 压力升高 🔴 高危 能量积累，静默期 立即预警，加强支护
静默 压力平稳 🟡 待观察 停采或蠕变初期 延长时间窗观察
活跃 压力下降 🟢 安全 小破裂稳态释放 正常推进
活跃 压力升高 🔴 高危 破裂加剧+支撑失效 停止作业，紧急预警

7.3 联合的核心价值：补充深度信息

微震与矿压联合的最大价值在于将矿压监测的一维信息扩展为三维。

实例： 工作面第 45 号支架（距运输巷 45 m）处压力传感器显示载荷异常偏高。

此时仅靠矿压数据，只知道"这个位置的顶板压力大"，但不知道压力量来自哪个层位的哪层岩石。可能的来源包括：

微震事件的深度定位补充信息： 45 号支架附近，过去 4 小时内微震事件集中在煤层上方约 28 m 处，对应老顶厚砂岩层位。

→ 联合判读： 老顶厚砂岩层正在断裂，导致了支架压力升高。

→ 精确措施： 对老顶层位实施水力致裂（定向切顶），而非盲目加固支架（若问题来自深部硬岩，浅部加固无效）。

联合信息流程：

支架压力异常（水平位置: 45号架附近）
  │  知道：顶板载荷异常大
  │  不知道：哪层岩石在活动
  │
  ├──→ 微震事件空间分布查询
  │    事件集中在 (X: 45m, Z: -28m)
  │    对应老顶厚砂岩层
  │    → 是这层岩石在断裂
  │
  ▼
联合判读结论：
  Z = -28 m 老顶厚砂岩断裂 → 压力来源明确
  → 治理决策：对该层位实施高压水力致裂
  → 目的：切断应力传递路径，卸压

7.4 完整联合监测体系

┌───────────────────────┐     ┌──────────────────────┐
│     微震监测系统      │     │     矿压监测系统      │
├───────────────────────┤     ├──────────────────────┤
│ · 事件目录（x,y,z,t） │     │ · 支架压力时序       │
│ · 能量时序 E(t)       │     │ · 锚杆/锚索轴力       │
│ · b 值时空演化        │     │ · 围岩位移量          │
│ · 震源机制（可选）    │     │ · 顶板离层量          │
│ · CT 速度分布         │     │ · 钻孔应力            │
└───────────┬───────────┘     └──────────┬───────────┘
            │                            │
            └──────────┬─────────────────┘
                       ▼
            ┌──────────────────┐
            │    联合分析引擎  │
            ├──────────────────┤
            │ · 空间关联：微震 │
            │   位置 ↔ 矿压    │
            │   传感器位置      │
            │ · 时序关联：能量 │
            │   加速 ↔ 压力变化 │
            │ · 层位判识：微震 │
            │   深度 ↔ 岩层柱  │
            │   状图            │
            │ · CT 速度 ↔       │
            │   应力集中区      │
            └────────┬─────────┘
                     ▼
            ┌──────────────────┐
            │  三维应力状态图像│
            │  (位置+深度+强度)│
            └────────┬─────────┘
                     ▼
            ┌──────────────────┐
            │  精确预警 +       │
            │  针对性加固方案   │
            └──────────────────┘

八、附录：关键公式汇总

震源定位

名称 公式 说明
P-S 距震源距 $D = \Delta t_{SP} \cdot \frac{V_p V_s}{V_p - V_s}$ 均匀模型，消掉 $t_0$
Geiger 线性化 $\mathbf{r} = \mathbf{G} \cdot \Delta\mathbf{m}$ Taylor 一阶展开
Geiger 最小二乘解 $\Delta\mathbf{m} = (\mathbf{G}^T\mathbf{G})^{-1} \mathbf{G}^T \mathbf{r}$ 法方程
阻尼最小二乘 $\Delta\mathbf{m} = (\mathbf{G}^T\mathbf{G} + \lambda \mathbf{I})^{-1} \mathbf{G}^T \mathbf{r}$ Levenberg-Marquardt

能量与震级

名称 公式 说明
G-R 能量-震级 $\log_{10} E = 1.5M + 4.8$ $E$ 单位：J
辐射能量 $E = 4\pi \rho v R^2 \int \dot{u}^2(t) dt$ 球面扩散
G-R 关系 $\log_{10} N = a - bM$ 频度-震级分布
$b$ 值 MLE $\hat{b} = \frac{0.4343}{\bar{M} - M_{\min}}$ Aki (1965)
$b$ 值标准误 $\sigma_b \approx b / \sqrt{N}$ 显著性参考

层析成像

名称 公式 说明
程函方程 $|\nabla T| = 1/v$ 走时正问题
走时离散化 $\delta\mathbf{t} = \mathbf{L} \cdot \delta\mathbf{s}$ 反问题矩阵方程
正则化目标函数 $\Phi = |\delta\mathbf{t} - $
转载请注明出处