基于微地震和计算机断层扫描（MS-CT）的煤层Regional应力场定量表征方法研究：一个案例分析

《International Journal of Coal Science & Technology》：Research on quantitative characterization method of regional stress field in coal seams based on microseismic and computed tomography (MS-CT): A case study

　　摘要：随着埋藏深度的增加，煤炭开采中由地应力引起的煤和岩石动力灾害的潜在风险也在增加。本研究试图利用微地震和计算机断层扫描（MS-CT）技术来定量表征煤层区域的应力场。为了实现这一目标，采用了一种综合方法，包括实验室实验、数值模拟和现场调查。实验结果表明，煤的波速与应力之间存在幂函数关系。将波速比与应力相关联的模型更适合用于应力场量化研究。该模型能够捕捉波速与应力变化之间的正相关关系，并通过标准化消除了样本差异的影响，使得拟合参数b可以在现场应用。本文以金家煤矿的22号煤层作为研究对象，使用MS-CT技术生成11224工作面的波速场分布云图，并结合实验建立的关系模型，将波速比与应力相关联，从而计算出区域应力的定量表征结果。同时，建立了一个详细的三维数值模型来模拟11224工作面的应力分布。将模拟得到的应力场结果与MS-CT的定量结果进行比较，发现两者的大致趋势是一致的，埋藏深度较大的区域应力较高。受上覆采空区影响的煤层区域应力较低。两种方法的最大误差值在±2 MPa以内，相关系数为0.79，表明它们之间存在强相关性。因此，MS-CT技术与实验得到的波速-应力耦合关系模型相结合，可以确定应力场分布并实现应力的定量表征。这为分析煤层不稳定性的力学模型和计算煤岩动力灾害激发判据提供了重要的基础。

1 引言
中国浅层煤炭资源的枯竭促使人们转向更深的开采方式，导致地应力和气体压力的显著增加。因此，地下矿山中动力灾害（包括岩爆和瓦斯突出）的风险也随之上升（Zhang等人2017；Xie等人2019）。这类灾害通常伴随着能量的剧烈释放，对采矿人员和设备构成威胁，甚至可能导致严重的事故，对煤矿的安全性和效率构成重大挑战（Li等人2021；Kong等人2021）。尽管已经有许多关于预测和预防这些重大风险的研究并取得了一些成果，但由于发生机制的复杂性，偶尔仍会发生事故，造成重大人员伤亡（Ma等人2020；Yin等人2024）。特别是，由于高地应力和复杂的地质结构导致的区域应力分布不均是造成这些挑战的主要因素。不均匀的应力会在局部区域产生应力梯度，严重影响采矿过程中煤和岩体的稳定性。特别是在瓦斯突出灾害中，局部异常应力会改变煤体的固有渗透性，导致采掘困难，无法达到标准。这增加了瓦斯突出矿井中煤和气体灾害的发生可能性（Zhao等人2023b, a）。因此，迫切需要有效的方法来确定煤层区域的应力分布。利用技术手段动态监测应力场分布和演变对于提高预测和预警煤岩动力灾害的能力至关重要。同时，这种方法还可以为精确的灾害预防和控制提供基础，大大减轻动力灾害的预防工作量，并显著提高采矿效率（Hosseini等人2017；Niu等人2022）。
目前，局部应力监测采用各种直接或间接方法，包括钻孔应力计、支护阻抗等直接监测方法，或钻屑量指数、声发射和电磁辐射等地球物理技术进行间接表征（Deng等人2023；Ma等人2023；Jia等人2009；Yang等人2021）。这些方法可以有效监测局部煤和岩体的应力状态，对于监测和提前预警煤岩动力灾害起着关键作用，并已得到广泛应用。然而，它们在实现大规模区域应力状态监测方面存在局限性（Zhao等人2022, 2023b, a）。通常，区域应力场研究依赖于地质力学仪器在特定点测量煤岩体内的原位应力，并将数据与数值模型进行比较以验证（Maxwell和Young 1996）。然而，要满足政府对于主动安全预防的严格要求，还需要更多的努力。补充和整合各种维度和方法可以显著提高煤岩动力灾害的预测效率。MS数据为三维监测岩体提供了可能性。值得注意的是，地震断层扫描是一种有效的技术，可以检测煤炭开采中煤和岩体的大规模区域应力状态，在评估煤岩动力灾害的风险水平方面发挥着关键作用（Chlebowski和Burtan 2021；Dou等人2014）。该技术利用MS数据，通过P波到达时间来计算地震波穿过煤和岩体的传播速度，最终生成描述空间速度分布的图像（Zeng等人2016）。根据MS事件的来源，地震波断层扫描可以分为“主动”或“被动”类型（Luxbacher 2008）。在他们的研究中，Young等人（1992）应用了主动地震断层扫描技术来研究加拿大Sudbury地区Falconbridge矿的P波速度结构与诱发地震活动之间的相关性。他们还通过钻探监测了高速和低速异常区域的应力，确认高速区域的岩应力升高以及低速区域的显著断裂。波速是表征地下矿山内应力状态的一个重要指标。

在生产矿井中，通常优先选择被动地震源进行地震断层扫描，因为数据来源是采矿操作不可或缺的一部分，可以在不产生额外成本的情况下实现微地震的动态连续监测（Zhao等人2023b, a）。Cai（2014）和Dou（2014）应用被动断层扫描技术连续跟踪应力变化，并识别出采矿过程中的高应力区域。结果表明，被动速度断层扫描作为一种工具，适用于持续监测高应力区域或评估地下煤矿中岩爆的风险（Cao等人2015；Zhao等人2025）。Friedel等人（1995）使用简单的机械模型和在地下环境中进行三维断层扫描，建立了应力与速度分布之间的强相关性（Zhou等人2022）。然而，大多数当前的研究只是基于波速分布对煤和岩体中异常应力区域进行定性分析，关于利用波速分布量化原位应力场分布的报告较少（Hanson等人2002）。
本文对应力与波速之间的相关性进行了实验探索，建立了一个将波速比与应力关联的耦合模型，适用于实际场景。随后，使用MS-CT技术检测了金家煤矿的区域波速场。通过将这些数据与实验室获得的应力波速比模型相结合，得到了煤层区域的应力场定量表征云图。接着，在金家煤矿11224工作面附近创建了一个精细的三维模型，模拟了煤层中的应力场分布云图。与模拟结果相比，我们分析了使用MS-CT定量表征应力场的精度和有效性。最终，我们探讨了提高MS-CT在量化应力场分布方面的准确性的方法，以及相关的预防措施和该技术的应用前景。

2 实验、分析与结果
2.1 实验设置
本研究使用的实验系统包括YAW-600电液伺服岩石测试机和DS5系列声波采集与分析系统，共同构成了实验平台。该平台具备煤动态加载功能，能够在实验过程中收集多种数据。这些数据包括加载的开始和结束时间、加载应力变化、声发射信号的初始激发时间，以及接收穿过煤层的声波的声发射传感器的信息。实验原理和物理图示如图1所示。其中，机械测试系统采用YAW-600微计算机控制的电液压力测试机，最大测试力为600kN，负载分辨率为3N，位移分辨率为0.3 μm。声发射传感系统的响应频率为50～400 kHz，前置放大器的放大倍数为20、40、60 dB可调，内置滤波器频带为20 kHz～1.5 MHz，输入阻抗＞10MΩ，输出阻抗为50ω。高速数据采集仪器具有16个数据采集通道，信号同步触发功能，最大频率为10 MHz，A/D转换精度为16位，输入信号电压范围为±10 V，可以实时、同步采集声发射信号。

图1 实验系统图
在开始实验之前，对样品施加500 N的预加载力，以确保加载装置与样品之间保持一致接触。加载过程在应力控制下进行。实验前，声发射系统配置为单通道重复激发模式，重复周期为3秒。激发信号脉冲的数量设置为两个，以便于识别活跃的激发信号。脉冲宽度为20微秒，脉冲周期为2,000微秒。声发射激发探头在整个加载过程中刺激振动信号，而接收探头捕捉穿过加载煤层的声波。前置放大器放大这些信号后将其传输到数据采集仪器。随后，使用计算机进行波速分析和计算。

2.2 样品制备
根据国际岩石力学委员会关于实验室和现场实验标准化的规定（ISRM标准），在确定煤芯中的声速时，样品的轴向和径向尺寸应至少是波长的十倍（Culshaw和Ulusay 2015）。在本实验中，所使用的声发射传感器的基频为0.55 MHz。考虑到煤样中的波速通常低于2500米/秒，所产生波的波长小于4.5毫米。因此，煤样的尺寸应当超过45毫米。对于这项研究，我们选择制造轴径为100.0毫米、径向尺寸为50.0毫米的煤样。这些煤样的原材料来源于晋家煤矿的22号煤层。由于煤的质地较软，获取完整的原煤样本存在挑战。先前的研究（Wang等人，2015年）表明，使用水合腐殖酸钠溶液作为粘合剂，并与粉碎的煤混合制成的煤球与原煤非常相似。因此，本研究采用这种方法来制备煤球样本。具体的生产 process 如下：1. 将原煤样本粉碎后，选择粒径在0到1毫米和1到3毫米之间的煤粉作为骨料，确保以76:24的比例均匀混合，然后彻底搅拌混合物，加入含水量为9%、粘合剂质量分数为20%的水合腐殖酸钠溶液。2. 将混合材料转移到样品模具中进行压实。对压机施加30千牛顿的压力并保持20分钟，从而製得直径为50毫米、长度为100毫米的圆柱形样品，端面不平行度不超过0.05毫米。3. 将制备好的样品放入红外干燥箱中，在105摄氏度下干燥24小时，然后自然干燥七天。图2展示了样品制备所需的原材料和成品样本。

图2 煤球样本示意图
全尺寸图像

2.3 纵波速度与应力之间关系的分析
纵波速度Vp通过公式（1）确定，其中Tm-Tn表示样品在加载过程中沿长度L（1-ε）的传播时间（Cao等人，2015年）。
$$V_{\text{p}} = \frac{L(1 - \varepsilon)}{T_{m} - T_{n}}$$ （1）
其中Vp是纵波速度，单位为米/秒；Tm是m探头激发纵波的时间；Tn是n探头接收纵波的时间；L是接收传感器与激发传感器之间的距离；ε是轴向应变率。
实验前测量的样品基本参数和从实验数据获得的机械参数如表1所示。纵波速度（Vp）随应力变化的规律如图3所示。显然，应力与纵波速度呈正相关趋势。速度的增加率是非线性的。通过拟合，确定使用幂函数可以准确捕捉到样品在初始应力下波速度的急剧增加，随后增长速率会减缓。

表1 煤球的基本参数
全尺寸表格

图3 样品中应力与P波速度（Vp）和波速比（r）之间的关系
全尺寸图像
拟合参数如表2所示，平均相关拟合系数达到0.964，这表明耦合关系模型可以很好地指示波速与应力之间的关系。P波速度在应力作用下的增加归因于裂缝和孔隙空间的闭合。随着应力的持续增加，由于新的裂缝（屈服点）的形成，速度增加率减小，这一点在之前的研究中已有说明（Friedel等人，1997年）。这些发现意味着高速度区域可能表明应力集中的区域，而低速度区域可能表示应力得到显著缓解的区域。

表2 波速和应力拟合参数
全尺寸表格

由于煤的异质性和声学时间标记错误的影响，从不同样本获得的波速存在显著波动。为了减少这些波动对拟合结果的影响，并确保这些结果在实际应用中的有效性，引入了波速比r作为替代原始波速的指标，波速比r计算为样品在不同压力条件下的波速除以其平均波速。波速比r与施加在样品上的应力之间的关系如图3中的黑点和线条所示。值得注意的是，波速比r随应力变化的趋势与波速随应力变化的趋势一致。使用这种方法有助于结果的标准化。基于以上考虑，我们研究了波速比与应力之间的关系。因此，在构建波速与应力之间的耦合模型时，使用了波速比的参数来建立耦合模型。

2.4 波速比与应力之间的耦合模型
实验得到的六组波速比与应力之间的关系以及拟合结果如图4所示，相应的拟合参数详细列在表3中。值得注意的是，所有六组样本的拟合系数均超过0.9，表明幂函数能够有效描述波速比与应力之间的关系。关于拟合参数的变化，考虑到幂函数的特性，我们认为参数a与应力大小无关，是实验系统的固有属性。不同的实验系统对应不同的参数a。另一方面，参数b表示波速比对应力变化的敏感性。随着波速的更快增加，其值也随之增加，这表明了煤的固有特性。先前对紧密砂岩的实验表明，在岩石样本峰值前的轴向应力增长段，波速与轴向应力之间的关系可以用公式V = cσ + d来拟合（Yang等人，2018年；Luxbacher等人，2008年）。通过比较不同应力状态下煤体的固有波速与固有波速的平均值，可以得到固有波速比rin。

图4 波速比与应力之间的关系及拟合曲线
全尺寸图像
表3 波速比-应力拟合参数
全尺寸表格

基于这一结论，图4中波速在加载后期几乎呈线性变化的段被认为是切线，可以用来确定不同应力下的煤的固有波速比。图4中测量波速比与固有波速比之间的绿色区域有效地阐明了波速随应力变化的机制。在初期阶段，微裂缝经历了压实，波速随着应力的增加而迅速增加。然而，随着微裂缝压实的进行，波速的增加逐渐减缓，测量的波速比与自然波速比趋于一致，遵循应力的变化趋势直到断裂。

使用Trimmean函数（在去除数据组中的最大值和最小值后计算剩余数据的平均值）对拟合参数结果进行处理，得到参数‘b’的平均值为0.266。公式（2）描述了样品应力与波速比之间的耦合关系模型。利用公式（2）可以推导出样品在不同波速下的应力值。

$$ r = a\sigma^{0.266} $$ （2）

3 煤层区域应力场的定量表征方法
在第二章中，获得了纵波速度与应力之间的响应规律，并建立了波速比与应力之间的耦合模型。在实验室尺度上证明了可以通过波速来表征应力。在工程现场尺度上，可以利用MS-CT技术检测煤层区域的波速场。然而，关于基于区域波速场对煤层区域应力场进行定量表征的研究很少。本章使用晋家煤矿作为工程试验场，研究基于MS-CT技术对煤层区域应力场的定量表征方法。

3.1 地震波层析成像原理
在地下巷道安装传感器接收MS信号时，地震波（v）的速度计算为射线路径距离除以信号从源到接收器所需的传播时间（t）。因此，时间（t）是源到接收器的慢度（s）的积分（Luxbacher 2008a）。假设第i个地震波的传播路径长度表示为li（米），传播时间为ti（秒），则关系式表示如下（Luxbacher 2008b）：
$$ v = \frac{l}{t} \to vt = l $$ （3）
$$ t_{i} = \int_{l_{i}} \frac{1}{v} dl = \int_{l_{i}} s dl $$ （4）
$$ t_{i} = \sum\limits_{j = 1}^{M} {s_{j} l_{ij}} \,\,\,\,\,\,\,\,\,\,\,i = 1, \cdots ,N $$ （5）
其中v表示速度，l表示源与传感器之间的距离，t表示传播时间，s表示慢度，N表示射线总数，M表示体素数量，ti表示第i条射线的传播时间，sj表示第j个体素的慢度，lij表示第i条射线在第j个体素中的距离。

在被动地震速度层析成像中，需要计算源位置和射线路径。初始速度模型可以计算射线路径上的距离和时间。方程（5）中的每个体素的时间、距离和慢度可以组织成矩阵。基于反演理论的速度计算公式如下（Luxbacher等人，2008）：
$$ T = LS \to S = L^{-1} T $$ （6）
其中T是每条射线的传播时间矩阵（1×N），L是每条射线每个体素的距离矩阵（N×M），P是每个网格单元的慢度矩阵（1×M）。

存在多种方法来解决反问题。其中，同步迭代重建（SIRT）是一种适合处理与速度模型相关的反问题的技术，因为它倾向于获得缓慢且收敛的解，从而产生相对稳定的结果。SIRT算法根据慢度分布计算射线的时间残差（观测时间减去计算时间），并使用公式（7）调整每个单元的慢度。模型参数不断迭代更新，直到射线的时间残差达到满意水平或预定的迭代次数（Hosseini 2017, 2012）：
$$ s_{j}^{k + 1} = s_{j}^{k} + \sum\limits_{i}^{N} {\frac{e_{i}^{k} l_{ij}}{\sum\limits_{j}^{M} {l_{ij}^{2}}} \frac{l_{ij}}{\sum\limits_{i}^{N} {l_{ij}}} $$ （7）
其中sk+1和sk分别表示第k+1次和第k次迭代后的慢度向量（M×1）和（N×1），N是射线总数（方程数），M是体素数量，sjk是第k次迭代后第j个体素的慢度，eik是第i条射线在k次迭代后的残差。

3.2 矿场和MS部署概述
本文在贵州省晋家煤矿的工程现场进行了振动波层析成像测试。晋家煤矿场位于盘关向斜西翼的南部；地层走向几乎为南北方向，属于向东倾斜的单斜构造，北部地层倾角较缓，矿区主要为山区，地形有明显的起伏。本文选择在22号煤层的11224工作面进行工程测试。22号煤为无烟煤，损伤类型属于II型，煤体的凝聚力系数为0.68。11224工作面的地面高度在+1900至+2155米之间，对应于地表垂直高度170至470米。该工作面是第22煤层第2段南部翼的第一条开采面，回风巷从1721大巷道开口，运输巷从211回风石门开口。工作面的总长度为450米，煤层的倾角在20°至24°之间，平均为22°，煤层厚度为0.93米至1.0米。工作面内没有显著的和中型的地质构造。图5显示了11224工作面附近井的上下控制图，该区域对应的地面上没有村庄或河流，是一个植被茂密的山坡地形。图5展示了11224工作面中MS传感器的布置和上下孔图对比。

在本文中，Jinjia煤矿的11224工作面部署了一个MS监测系统。MS传感器的布置方式如图5中的1#至7#所示，它们环绕着11224工作面。在正式建立MS监测系统之前，我们进行了现场爆破实验来校准MS系统的准确性。实验结果表明，当定位波速为4400米/秒时，MS系统的定位误差小于5米。因此，MS事件的平均传播速度为4400米/秒。

3.3 区域应力场的定量表征结果
在本节中，我们对采集到的MS数据进行了层析成像分析，以揭示地下煤岩体结构中的速度和应力场的详细信息。本文用于反演的数据包括2017年7月18日至2017年8月2日期间11224运输巷道的非驱动记录。反演过程包括以下步骤：（1）对地震数据进行分类，以分离出定位误差低、到达时间残差和最大到达时间小的MS事件，确保MS射线在反演区域内有最佳的覆盖度；本研究选择了65个合格的MS监测数据点。MS事件和射线覆盖分布分别如图6和图7所示。（2）提取MS到达时间和位置数据。（3）划分反演区域的三维网格，并建立三维反演模型。本文的反演模型在水平x方向、垂直y方向和垂直z方向上被分割为31×23×11个网格，网格间距分别为20米×20米×20米。（4）使用SIRT算法进行解析反演的迭代计算。

图8a显示了从MS数据中获得的11224工作面的波速场分布。与图8b中煤层的埋深相比，可以看出波速分布受到埋深的影响很大。埋深越深，地应力越大，导致煤和岩体的密度增加，构造不连续性闭合，波速增加。因此，图8a中的高速区域代表高应力区域。这与Li等人（2019年）通过类似模拟实验得到的结果一致。尽管波的传播受到岩石物理和机械性质以及岩体结构特征的影响，但在相同的埋深下，这些特性保持不变。因此，图8a可以反映应力分布的状态。

第2.4节建立了波速比与应力变化之间的耦合关系，其中参数a与测试系统相关。Jinjia煤矿MS系统中用于定位MS事件的平均波速为4400米/秒。地震波层析反演计算显示，11224工作面范围内的平均波速为4410米/秒。因此，两者的平均波速比对应于1.002。波速-应力测试结果的分析表明，当波速比为1时，应力基本上等于平均应力。11224工作面煤层的平均埋深约为240米，因此煤层的平均应力为6兆帕。将（1.002, 6）代入波-应力耦合方程（2），我们可以确定在野外条件下的检测系统参数‘a’为0.621。在实际应用中，应力与波速比的转换公式表示为方程（8）：$$ \sigma = 1.61r^{ - 0.266} $$ （8）基于方程（8）和11224工作面上检测到的速度场，图9定量表征了11224工作面上的应力场分布。

3.4 煤层区域应力场的定量表征方法
本文第2章和第3章的研究表明，基于实验室样品中波速和应力的实验结果，结合现场振动波层析检测，可以实现对突出煤层应力的定量表征。基于上述研究结果，图10展示了突出煤层应力场定量表征的流程图。
具体步骤如下：
步骤1：获取应力定量表征的参数。选择要进行应力定量表征的区域，制作该区域的标准煤样，对样品进行单轴加载损伤试验，并通过本文建立的波速比-应力耦合关系模型获得现场波速确定系统的参数a和波速比对应于应力的灵敏度参数b。
步骤2：现场波速场检测。在待测试区域内部署MS监测系统，监测煤岩体内部发生的MS事件，并筛选出符合振动波层析成像要求的100至150个MS事件。根据成像精度要求划分模型网格，并为地震波成像设置合理的参数，以获取待测试区域的波速分布数据。计算波速场的平均波速，并进一步计算待测量区域的波速比r。
步骤3：根据检测区域的平均埋深计算该区域的平均应力σ。基于MS定位的平均波速与检测区域的平均波速之间的速度比计算，结合对应于平均埋深的σ和实验得到的b'值，然后将其代入波速-应力耦合关系模型中，计算出与现场波速检测系统相对应的参数a'。
步骤4：通过将参数a'、参数b以及从现场检测得到的波速比r代入耦合的波速-应力方程中，可以得到煤层的区域应力场。

4 应力定量表征的可靠性分析
为了评估通过现场波速检测表征煤层应力场分布的准确性，本文首先通过数值模拟分析了11224工作面的煤层区域应力场分布。随后，将其与定量表征的应力场分布进行比较，以验证煤层区域应力场定量表征方法的可靠性。

4.1 三维数值模型和边界条件
第3.2节已经总结了11224工作面的情况。工作面内没有显著的断层或地质构造，只有在工作面开口右侧的滑动不良区域内有异常地质体，其边缘距离11224工作面的切割眼大约30米。根据早期采矿区域的勘探数据，建模区域主要涵盖了Jin 1采矿区内的11224工作面，包括11224工作面及上方的已开采工作面。建模范围为：X=4000至5200，Y=1100至1600，Z=+1550至地表。
建模范围涵盖了复杂的条件，包括表层厚度、煤层厚度、倾角和顶板厚度的变化。从11224工作面附近的勘探线中提取了不同煤层的高程数据和剖面数据，并从勘探阶段的地质钻孔记录中进行了处理。这些数据被组织成点或线格式，并导入Rhino建模软件中，构建了一个三维数值模型，如图11所示。该模型包括了11224工作面及其上方的开采工作面，每个工作面的分布都在图11中展示。模型包含789,852个区域和624,816个网格。

图11 三维模型。a 11224工作面附近的数值模型；b 11224工作面附近的工作面分布
根据Liupanshui采矿区煤层深度范围在135至632米内的现场应力测试结果，水平应力介于2.14兆帕至12.53兆帕之间，而垂直应力在3兆帕至14.76兆帕之间。模型的边界条件根据这些现场应力测试结果进行配置。具体来说，模型在X和Y方向的两侧都限制了正常位移，并施加了水平应力。模型底部的平均位移为零，意味着平均速度和位移为零，切向力也为零。模型的顶部代表地表，并建立了自由边界。重力加速度（g）设定为10米/秒2。

4.2 物理和机械参数
本构模型必须能够有效地捕捉煤和岩石的机械特性，并反映它们在拉伸剪切应力下的破坏行为。在煤岩体的弹性范围内，应力和弹性应变表现出线性关系。然而，当超过屈服点时，煤和岩体会发生破坏，强度减弱。该数值模拟采用了基于经典Mohr-Coulomb准则的应变软化模型来更准确地模拟这种行为（Tan等人，2021年；Shen等人，2020年）。该模型特别适合模拟煤和岩石的机械行为，其屈服函数可以按照方程（9）表示：当fs<0时，发生剪切破坏；当ft<0时，发生拉伸破坏；当fs>0且ft>0时，材料不会破坏（Tan等人，2021年）。
$$\left\{ \begin{gathered} f_{\text{s}} = \sigma_{1} - \sigma_{3} \xi + 2c\sqrt \xi \hfill \\ f_{\text{t}} = \sigma_{\text{t}} - \sigma_{3} \hfill \\ \end{gathered} \right. $$ （9）
其中fs是剪切屈服函数，σ1和σ3分别是最大和最小主应力（单位：兆帕）；σt是抗拉强度（单位：兆帕）；c是凝聚力（单位：兆帕）；\(\xi = \frac{1 + \text{sin}\varphi }{{1 - \text{sin}\varphi }}\)；φ是以度为单位的摩擦角。
为了确保应变软化模型输入参数的合理性，需要进行彻底的参数检查。在FLAC 3D软件中，使用高度为100毫米、直径为50毫米的圆柱形试样进行了系列的单轴压缩数值模拟测试。通过反复试验调整了试样的凝聚力和软化比例，以使其与实验室实验结果一致。校准过程和结果如图12所示。模型中的煤和岩石力学参数详细列在表4中。

表4 数值模型中煤和岩层的物理和机械性质

4.3 应力场的数值模拟结果
在第4.3节建立的3D模型被导入FLAC 3D 7.0软件中，并考虑了矿中煤和岩体的物理和机械参数以及地面应力的分布特性，计算了初始地应力。图13显示了达到平衡后的三向应力解的结果。在初始地应力分布解析之后，模型的应力在地面上传递均匀，没有应力集中或突然变化，模型的垂直和水平力是平衡的。模拟结果与煤岩体随埋深变化的应力模式相符，证实了该模型和机械参数适用于本研究。这还表明了在解决模型时边界条件和参数的合理性。图13展示了模型三向应力求解的结果。

由于11224工作面覆盖了许多已完成的采矿面，因此根据11224工作面附近的3D数值模型提取了已开采的煤层，遵循上层煤层的采矿条件。当上层煤层完整时，如图14a所示，应力分布与工作面的埋藏深度非常吻合。由于地表起伏，11224工作面中煤层的埋藏深度各不相同，埋藏深度较大的区域应力也更高。在11224工作面填充了止壁后，所得到的煤层应力分布如图14b所示。应力场分布的整体趋势与未开采时相似，埋藏深度较大的区域应力更高。然而，上层煤层的开采对11224工作面产生了减压效果，导致煤层应力整体下降，并且在局部区域有显著降低。

4.4 定量表征的可靠性
在本节中，将对数值模拟和MS-CT的数据进行比较分析，以得出有意义的结论。为了更准确地评估MS-CT应力场与模拟应力场分布之间的差异，我们对从数值模拟得到的11224工作面的应力云图值与从MS-CT检测得到的值进行了图像减法处理。两者之间的差异如图15所示。可以发现，在4400米至4500米的范围内，应力偏差较为显著，这主要是由于上层煤层开采后11224工作面上的压力卸载不均导致的定量表征应力与模拟应力之间的显著差异。同时，在4600米至4800米的范围内，应力偏差不仅受到上层煤层开采后11224工作面上压力卸载的影响，11224工作面切割眼右侧的无煤滑动带也对该区域的应力有一定影响。但大多数区域应力误差在±2MPa范围内，仅有极小部分超过2MPa。皮尔逊相关系数（P）为0.79，表明变量之间存在强相关性。通常，相关系数分为以下几类：0.8–1.0表示强相关性，0.6–0.8表示较强相关性，0.4–0.6表示中等相关性，0.2–0.4表示弱相关性，0.0–0.2表示非常弱或无相关性。因此，从MS-CT检测得到的应力与数值模拟得到的应力之间存在稳健的相关性。

5 讨论
应力是导致煤岩动力灾害的主要原因之一。数值模拟是分析工作面应力的重要方法。在收集了采矿区域的地质数据和采矿布置后，可以进行模拟并分析区域应力场。然而，由于地质数据无法非常详尽，且精细建模的工程工作量很大，数值模拟通常是对分析区域内地质体的简化建模和分析。先前的研究表明，地震层析成像的速度分布特性可以定性解释应力异常。然而，仍需要更多研究来定量预测区域应力场。此外，目前缺乏基于此建立煤岩动力灾害预测指标值的研究。在本文的第二部分，我们建立了应力与样品波速比之间的实验关系。随后，我们在样品的采样地点进行了地震层析成像。现场工程实验的结果表明，P波层析成像可以反映反演区域内的垂直应力，并且与理论预测的应力有很好的相关性，相关系数约为0.79。本文的结果表明，基于MS-CT的应力场定量表征方法是可行的。其使用的先决条件是现场有MS设备，并且分析区域需要被MS设备覆盖。该技术基于现场的MS监测数据进行分析和计算，从而获得的应力将更符合地下复杂地质条件下的应力分布。随着MS设备在矿井中的普及，本文提出的煤层区域应力表征方法具有更广泛的应用前景。

尽管基于MS-CT的应力定量表征具有可靠性，我们也认识到定量应力表征方法仍有改进的空间。这主要是由于实验室测量的参数a'和耦合波速-应力模型中的参数b与现场实际参数之间存在差异。因此，确保现场MS-CT技术测量出的应力绝对准确是一个挑战。为了解决这个问题，我们建议在现场安装至少两个在线应力监测传感器。每次反演后，将传感器位置的波速比和应力值代入实验关系模型中，调整参数'a'和'b'以提高准确性。结合这两种信息类型，少量的离散测量可以获得更精确的应力分布。此外，在现场煤岩体发生大规模破坏时，断裂带可以通过低速区域来显示，这与Friedel等人（1997年）的发现一致。然而，由于我们使用的应力波速比模型是基于破坏前的样品，因此在量化应力识别时可能存在较大误差。在实际应用中需要仔细考虑这一方面。

研究表明，忽略异常应力场分布的影响往往会增加煤岩动力事故的风险（Zhou等人，2022年；Li和Pi，2019年）。然而，并非所有异常区域都会导致此类灾害。因此，在现场生产过程中，精确识别危险区域至关重要。以往的研究仅限于对应力异常的定性分析，导致现场应用中危险区域范围过广且存在许多不准确性，限制了实际效果。从检测到的波速分布、应力波速实验关系模型和实时应力监测数据得出的区域应力场的定量结果为进一步开发考虑应力异常的煤层不稳定力学模型提供了重要支持。此外，这些结果有助于计算考虑应力异常的分析解，构建判断指数和风险预测方法，对于有针对性的灾害预防和控制具有重要意义。

6 结论
我们以金家煤矿的11224工作面为重点，实验研究了煤岩体中的纵波速度与应力之间的关系。我们建立了应力与波速比之间的耦合模型。此外，我们进行了数值模拟和地震层析成像检测，以分析11224煤层中的应力场分布。比较并分析了这两种方法得到的应力场分布特征。得出以下结论：
（1）煤的纵波速度随着应力的增加而稳定增加，在达到煤的峰值应力之前，两者之间的关系遵循幂函数。波速比与应力之间的关系反映了波速随应力的变化特性。这种方法有助于消除由于不同样品之间的差异而导致的拟合结果的不稳定性，提高了其在现场的应用性。
（2）对11224工作面的区域应力场分布进行数值模拟，并与地震波层析成像检测得到的应力场分布特征进行比较，发现两者之间具有一致性。区域应力误差值在±2MPa范围内，皮尔逊相关系数（P）为0.79，表明存在强相关性。波速比与应力之间的耦合关系以及波速场的检测结果，允许量化应力分布和大小。
（3）区域应力场的定量检测可以支持进一步开发考虑煤层不稳定性中应力异常的力学模型。此外，它有助于计算判断指数和考虑应力异常的煤岩动力灾害发生风险预测方法。这对于实施有针对性的灾害预防和控制以及制定合理的采矿计划具有重要意义。