深埋隧道时间延迟型岩爆诱发灾害机理研究：基于微震数据驱动的方法

《TUNNELLING AND UNDERGROUND SPACE TECHNOLOGY》：Mechanism of Time-Delayed Rockburst-Inducing disasters in Deep-Buried Tunnels: A microseismic Data-Driven Methodology

　　

•微震能量反演与威布尔统计重构围岩的异质性劣化特征。
•局部软化导致拱顶非对称荷载分布，并引发衬砌外压内拉的破坏模式。
•爆破扰动触发弹性应变能的快速释放，导致时间滞后岩爆的动破坏。
•案例研究表明，隐性岩爆事件发生前，b值与D值呈同步下降趋势。

## 1 引言

随着地下工程向更深地质层延伸，岩爆事件的发生频率不断攀升（Akdag 等，2018；Li 等，2025；Luo 等，2025；Gong 等，2025）。对于采用钻爆法施工的深埋长隧道，部分岩爆并非在掌子面爆破时立即显现，而是延迟数天甚至更久后才发生，出现在已完成初期支护或二次衬砌的后方区段（Ma 等，2025；Hu 等，2026）。更为关键的是，这类潜在破坏可能一直潜伏至隧道竣工并投入长期运营阶段。西南某隧道在竣工两年后发生了延迟岩爆，导致无砟轨道严重隆起翻浆，二次衬砌出现局部破坏。由于支护表面呈现完整状态，此类事件隐蔽性极强，一旦突然爆发，极易造成支护系统的灾难性溃塌及人员伤亡。因此，深入探究被高刚度支护所掩盖的时间滞后型岩爆的致灾机理，具有重要的学术价值和工程意义。

在时间滞后型岩爆的触发机制方面，当前研究主要聚焦于高地应力作用下岩石的时变变形与应力重分布特征。已有研究广泛采用蠕变模型来理解卸载岩体的流变行为（Schoenball 等，2014；Liu 等，2026；Wang 等，2026）。对于深埋双线地下隧道，时间滞后型岩爆是内外致灾因素相互耦合的结果（Li 等，2017；Keneti & Sainsbury，2018；Ren 等，2025）。先期开挖爆破不可避免地在邻域围岩中形成开挖损伤区（EDZ），导致邻近岩体力学性质不同程度劣化（Chen 等，2025）。在该软化区采用高刚度二次衬砌，会造成系统内部显著的刚度失配（Mei 等，2026；Tian 等，2026）。在外部扰动方面，除施工期间邻近隧道的偶发爆破干扰外，隧道运营后的列车周期性振动构成有害动力荷载，导致显著的能量释放。现有静态理论尚不能充分解释这种由局部弱化区与刚性衬砌组成的隐蔽结构在爆破冲击、列车高频振动等多种外部动力荷载作用下能量积聚与突发性失稳的机理（Manouchehrian & Cai，2017；He 等，2023；Wu 等，2025；Chen 等，2026）。

微震（MS）监测技术能够实时探测岩体中的微裂隙信号，已成为岩爆早期预警的关键工具（Ma 等，2015；Zhao 等，2019；Ma 等，2020；Ma 等，2022）。目前，基于微震数据的统计参数分析已较为成熟（Dai 等，2017；Xue 等，2021）；然而，在数值计算领域，如何将现场监测的微观数据与岩体本构模型的演化进行物理表征仍然是一个重大挑战。传统连续介质力学模型通常采用预设的静态几何缺陷来处理开挖损伤（Wang & Qian，1995；Yang 等，2018；Yu 等，2022；He 等，2026）。这种简化方法无法准确再现岩体实际的非均质初始物理状态——即经历先期扰动后内部力学性质呈非均匀退化。尽管已有研究提出利用威布尔分布函数（Tang 等，2000）建模岩体力学参数的非均质性（如RFPA方法），以模拟岩体的脆性破坏特征（Tang 等，1998），但要完全将非均质统计模型与实际现场微震能量相结合，以精确重构局部损伤区，仍需进一步深入研究。

本文提出一种数据驱动的数值模拟策略，利用微震技术克服上述挑战，聚焦于深埋高原铁路隧道拱顶附近一次典型的延迟岩爆事件。首先，计算地震学参数（b值）和空间分形维度（D值），研究岩爆前兆的时空演化规律。随后，展示一种基于微震能量的初步损伤等效反演方法。通过Python二次开发，在动态模拟软件（Ansys/LS-DYNA）中引入威布尔分布，重现经历先期开挖影响后岩体的非均质初始物理状态。模拟围岩的动态响应，即由高刚度衬砌限定的局部软化区在远场冲击下发生的响应。本研究旨在探讨特殊埋深结构在动静组合应力作用下的能量演化与剥落失稳机理，为钻爆法施工条件下的岩爆灾害防治提供坚实的理论基础。

## 2 工程背景

### 2.1 工程概况与地质条件

本研究以西南某在建高原铁路深埋双线隧道为对象。所研究段（里程DK104+100至DK106+640）长度约2.54 km，最大埋深1999 m。两隧道中心间距26 m，均采用钻爆法施工，平行交替掘进。隧道常规开挖宽度9.47 m，高度9.92 m，开挖断面面积约73.51 m2。

研究区域位于活跃的地质构造带，地处欧亚板块、印度板块与太平洋板块的挤压带，节理裂隙发育。图1表明，隧道穿越地层主要为燕山期花岗岩和砂质板岩。花岗岩体因深部构造运动而具有显著各向异性。超前地质预报和室内岩石力学试验结果表明：P波速（v?）为5140–5590 m/s，S波速（v?）为2950–3440 m/s，泊松比为0.14–0.35，静态弹性模量为30–58 GPa，动态弹性模量为43–70 GPa，单轴抗压强度（UCS）在92–200 MPa之间。岩体分类为较坚硬岩。然而，受深部高地应力和断层带影响，围岩完整性较差，综合划分为Ⅲ级围岩。隧道掌子面和拱顶在开挖过程中尤为容易发生挤压剥落和块体坠落。

图1 工程地质条件（a）地质剖面（b）隧道断面尺寸

此外，在地下工程中，构造活动常使水平应力超过垂直应力（Ghorbani 等，2020；Du 等，2022）。现场水力压裂地应力测试表明，该研究段具有典型的高地应力特征。最大水平主应力为39.16 MPa（方位角N56°E），垂直主应力为35.16 MPa，最小水平主应力为27.56 MPa。初始地应力场与隧道轴线之间复杂的夹角关系加剧了深部岩爆的发生可能性。

### 2.2 隧道破坏特征

本节对隧道中发生的一次典型且剧烈的时间滞后型岩爆事故进行分析。2025年11月17日晚，按照施工顺序对左右线隧道进行了连续爆破，左线爆破时间为20:58，右线爆破时间为21:37。左线爆破结束后约2小时，即23:27，传来一声巨响，随后在DK1061+472至+496加固区段发生了严重岩爆。破坏区跨越施工缝两侧，位于左线掌子面后方259.5 m、右线掌子面后方97 m处，如图2所示。隧道采用钻爆法开挖，每个循环平均进尺约3 m。在正常施工条件下，一般每日爆破一次。因此，当掌子面通过受影响区段后，该区域在衬砌破坏发生前，长时间承受了反复的施工爆破扰动。

图2 岩爆发生概况

图3表明，此次时间滞后型岩爆局限于二次衬砌从拱顶到拱肩的区域。宏观表现为混凝土衬砌大面积剥落和块体脱落，并伴有沿拱墙环向轴线向上延伸的拉伸裂缝。2026年4月2日，在距掌子面210至220 m处发生了类似的拱顶衬砌开裂事件。这表明在高刚度支护后方发生的时间滞后型岩爆并非偶然或孤立事件。

图3 破坏区现场情况

早期开挖记录和详细的现场勘察发现，这些易发时间滞后型岩爆的区段存在显著的地质-工程前兆。在初步爆破开挖阶段，拱顶2点至5点位置（对应严重破坏区DK1061+484.4–487.4）出现了严重的超挖凹坑，如图4所示。凹坑径向深度达1至1.2 m，轴向长度约2.7 m。由于拱顶定位精度和回填施工方法（包括振实不充分等原因），这些凹坑极易在二次衬砌下方形成物理空洞。这些隐蔽空洞不仅为开挖损伤区（EDZ）的退化扩展提供了自由面，而且局部削弱了围岩对衬砌的径向支撑。在上述边界条件下，叠加邻近隧道爆破的周期性动力荷载扰动，围岩损伤逐步累积，最终导致累积应变能超过临界阈值，引发灾难性的时间滞后型岩爆。

图4 岩体初始损伤情况

## 3 隧道开挖期间的微震活动特征

### 3.1 微震监测系统布置

本项目采用加拿大ESG微震监测系统，实时捕捉深部围岩微裂隙的动态演化信息，实现岩爆实时预警。该系统由Paladin数据采集单元、灵敏度为40 V/g、频率响应范围为20 Hz至5 kHz的高灵敏度加速度传感器，以及Hyperion数字信号处理平台组成（Cheng 等，2017；Duan 等，2024）。当岩体发生微裂隙时，释放的应力波被传感器捕获并转化为电信号。信号通过光缆传输至数据采集设备，再经通信链路中继至Hyperion平台进行数字信号转换。处理后的数据最终上传至云平台，技术人员可下载波形并生成风险报告，供项目总工审阅。

图5、图6表明，传感器阵列布置在平行隧道的左侧边界，同时覆盖两条隧道的开挖区域。监测网络采用动态部署方式，随掌子面推进而同步扩展。共设置3个监测断面，包含9个传感器。第一排距隧道壁50至100 m，排间距为20至25 m。该技术采用神经网络方法，通过智能波形识别和滤波，对真实裂隙信号进行精确识别和源参数的有效反演。图7清晰展示了该技术处理的标准微震信号与干扰噪声之间的区别。

图5 微震监测流程图

图6 微震监测系统布置

图7 典型波形

### 3.2 微震事件的时空演化

岩体中微裂隙的形成与存在始终伴随着能量的耗散与释放（Zhang 等，2021；Liu & Gu，2026）。

本节通过微震事件的时空聚类度和能量释放特征，重建了延时岩爆的完整过程。图8(b)展示了研究区域微震事件的空间分布及能量密度等高线图。从八月至十二月，微震事件主要集中在隧道拱顶和拱肩区域，与观测到的收敛变形及破坏位置精确吻合。图8(a)统计展示了微震频次与累计能量释放的时间演化规律。

在初步开挖阶段（即阶段I），由于掌子面爆破扰动剧烈及围岩卸荷效应显著，微震频次持续处于较高水平。最大单次微震事件能量达到约2×10? J，导致累计能量曲线呈阶梯式快速上升。随后，初始支护体系（包括锚杆和喷射混凝土）有效限制了围岩径向位移。系统进入稳定阶段（即阶段II），其特点是微震频次大幅下降，能量释放趋于平稳并维持在平台状态。

然而，在十月底二次衬砌完成后，远处爆破产生的应力波虽然不足以直接形成大规模新裂隙，却作为低频循环荷载，导致围岩内部持续发生疲劳损伤累积。在二次衬砌后方岩体中，长时间积聚的应变能逐渐逼近其临界阈值。这一过程引发了微震活动的剧烈突增，呈现群集式暴发特征，标志着阶段III的开始。在此期间，单次微震事件强度上升至4×10? J，构成延时破坏的重要前兆。随后，宏观裂隙发育与岩体剥落发生，储存的应变能被强烈释放。系统随后进入残余应力调整阶段（即阶段IV），其特点是事件数量稀少，单次事件能量降至1×10? J以下。

### 3.3 微震参数演化特征

本工作在直观的震源密度和能量指标之外，引入地震学b值来指示裂隙规模，以及非线性动力分形维数D值来描述空间分布特征。这些判据被用于对岩体内部损伤演化机理进行系统研究。

#### （1）裂隙规模参数

微震活动的震级-频度关系遵循古登堡-里克特定律（G-R律），该定律由Gutenberg和Richter（1944）在全球地震活动研究中提出：

$$\lg N = a - bM \quad (1)$$

式中，M表示震级，N表示震级大于或等于M的微震事件累计频度，a和b为拟合常数。其中，b值是反映微震事件震级相对分布斜率的关键参数。b值降低通常表明岩体中小尺度微裂纹逐渐减少，同时大尺度宏观裂隙占比显著增加。

#### （2）空间分形维数

分形理论研究具有自相似结构的无限复杂不规则几何形态，分形维数是定量表征分形复杂性的关键参数（Xie等，2011）。相应地，岩体破坏过程中微震事件的空间分布表现出明显的分形特征（Wang等，2023）。基于盒子计数法的计算原理，微震空间分形维数D可通过以下公式拟合得到：

$$\ln N(\varepsilon) = -D \ln(\varepsilon) + C \quad (2)$$

式中，$\varepsilon$表示用于划分微震空间的网格尺度，$N(\varepsilon)$表示包含微震事件的非空网格总数，C为常数项。拟合直线斜率的绝对值即为D值，其大小反映微震事件在三维空间中的分散或聚集程度。

为兼顾前兆特征的时效性与统计有效性，本研究采用滑动窗口法，窗口大小为30个事件，滑动步长为5个事件（Zhang等，2015）。采用该方法对关键监测断面的b值和D值进行动态计算，结果如图9所示。在开挖扰动期和支护稳定期，两参数均维持高位波动，b值在1.0至1.1之间，D值在2.4至2.6之间。此阶段岩体损伤以小尺度微裂隙为主，裂隙在三维空间中呈随机无序弥散分布。

进入以远场动力荷载为主导的延时灾害孕育阶段后，参数表现出显著的同步劣化趋势。一方面，b值发生不可逆的单调下降，降至0.7以下，表明微裂纹正在加速扩展并发生贯通。另一方面，D值打破了此前的高维散射状态，阶梯式骤降至约1.7。这一空间维数降低现象表明，原本孤立弥散的微裂纹群经历了强烈的空间定位化过程，向特定的二维剪切带或断裂面汇聚与聚集。

上述以b值骤降为特征的裂隙规模贯通与以D值骤降为特征的空间维数降低和聚集现象之间的强时间同步性，构成了识别高刚度支护背后隐蔽延时岩爆的定量前兆预警。

### 4. 数值模拟分析

#### 4.1 几何模型与边界条件

为再现延时岩爆的动力触发过程，本研究基于实际隧道数据和开挖条件，使用LS-DYNA建立了一系列双线隧道三维实体模型，如图10(a)所示。模型整体尺寸为70 m × 40 m × 17 m，分别代表长度、宽度和高度，衬砌厚度设为2 m。模型采用六面体实体单元，共计748,578个单元。单元质量评估结果（图10(b)和图10(c)）表明，整体网格质量优于0.8，长宽比控制在1.0至2.0之间，满足显式动力计算的收敛要求。围岩与衬砌之间的界面采用LS-DYNA中的面-面接触算法处理。界面在受拉荷载作用下或局部破坏后允许张开，分离后的接触面在后续碰撞中允许重新接触。同时，在模型边界施加无反射边界条件，消除应力波在人工边界上的反射与叠加效应，有效模拟深部岩体的半无限空间。初始地应力场根据现场水力压裂测试结果施加。

图11(a)展示了模型内部视图，显示了内部衬砌和前方掌子面结构。对于左线延时区的围岩，采用齐次性指数m=8，弹性模量范围为6至35 GPa，单轴抗压强度分布在0至90 MPa之间，如图11(c)和图11(d)所示。

现场爆破荷载通过作用在开挖边界上的衰减应力波表示。峰值压力由现场使用的炸药参数结合经验爆破压力公式计算得出：

$$P_0 = \frac{n \cdot P_e \cdot 2\rho C_p}{\rho C_p + \rho_e D} \quad (3)$$

$$P = \frac{P_0}{R^\gamma} \quad (4)$$

式中，P为所考虑位置的冲击波压力，P?为初始峰值压力，R = r/r?为相对传播距离，r为距装药中心的距离，r?为装药半径，γ为传播衰减系数。此外，ρC?和ρ?D分别表示岩体和炸药的波阻抗，P?为爆轰压力，ρ?为炸药密度，D为爆轰速度，C?为岩石纵波速度。

现场使用的炸药为Φ38 mm 2号岩石乳化炸药，如图11(b)所示。这些波的时程函数可表示如下（Ma和An，2008）：

$$P = P_0 \cdot \frac{e^{-\alpha t} - e^{-\beta t}}{e^{-\alpha t_0} - e^{-\beta t_0}} \quad (5)$$

$$t_0 = \frac{\ln(\beta/\alpha)}{\beta - \alpha} \quad (6)$$

在这些方程中，数值计算中峰值应力P?设为80 MPa，加载率参数α和β分别为0.00811和0.01217。每次爆破荷载持续时间为100 μs。现场左线和右线隧道分别在20:58和21:37进行爆破，对应间隔为39 min。在数值模型中，右线爆破荷载在左线爆破引起的动力响应稳定后施加，两个加载阶段之间的间隔为10 min。

#### 4.2 基于微震能量的损伤反演与非均质重构

传统数值模型通常将开挖损伤区简化为静态几何缺陷，无法准确捕捉其力学性能的渐进劣化。为克服这一局限性，本研究提出了一种结合微震能量反演与Weibull统计分布的初始状态重构方法，具体技术流程如图12所示。数值计算采用逐步隐式-显式求解程序。首先，在隐式分析阶段施加初始地应力场以建立静力平衡。当隐式解满足残差力收敛准则，且代表性节点的位移和应力在连续计算步中趋于稳定时，认为初始应力场已达到平衡状态。

随后，将平衡后的隐式解所得节点位移写入文件，用作显式动力分析的初始变形构型。同时，将隐式单元转换为相应的显式单元，移除静力平衡阶段使用的冗余约束，并根据微震损伤反演结果更新延时区内的单元材料参数。最后，在此基础上施加爆破动力荷载，利用显式动力分析模拟后续响应。

##### （1）微震损伤反演

假设岩体单元为闭热力学系统，根据热力学第一定律，作用在其上的功满足方程U = U? + U_d，其中U?表示弹性应变能，U_d表示耗散能。考虑损伤变量D，单元可释放的弹性应变能可简化为：

$$U_{eD_i} = \frac{1}{2} E_{D_i} \left[\sigma_1^2 + \sigma_2^2 + \sigma_3^2 - 2\nu(\sigma_1\sigma_2 + \sigma_2\sigma_3 + \sigma_1\sigma_3)\right] \quad (7)$$

式中，E_{D_i}为损伤状态下的弹性模量，ν为泊松比。引入地震效率η，定义为微震系统捕获的能量ΔU'与实际释放应变能ΔU之比。现场每次爆破平均装药量约为450 kg，对应估算总爆破能为1.575×10? J。ESG微震监测系统对相应爆破事件记录的平均地震能约为2.19×10? J。据此，通过定点爆破标定计算出地震效率η = 0.001389%。对于延时爆破爆腔周围的软化区，其累积损伤变量D可定义为（Xie等，2009；Wang等，2025）：

$$D_i = \frac{\Delta U}{U_e} = \frac{\Delta U'}{\eta \cdot U_e} \quad (8)$$

D_i = 1表示岩体单元已完全破坏，D_i = 0表示完整无损状态。

基于微震数据，建立以掌子面中心为参考点的局部坐标系，将微震事件的三维坐标转换至数值模型坐标系。然后根据每个事件的坐标确定对应的有限单元。

当多个事件映射到同一单元时，材料属性按照时间顺序逐步降低，以考虑累积损伤。在第 i 次事件后，单元的弹性模量按以下公式更新：

(9) ??????=(1?????)·??????

其中 ?????? 表示第 i 次事件前该单元的弹性模量。因此，多次微震事件的累积效应通过材料参数的逐步更新来表示，相应结果见表1。

(2) 非均匀物理初始状态重建

表1 微震损伤数据库

| 序号 | 计算坐标/m | 微震能量/J | 损伤/% |
|------|-----------|-----------|-------|
| | X | Y | Z | | |
| 1 | ?2.94 | 2.62 | 5.68 | 266.98 | 5.12 |
| 2 | 5.72 | 4.79 | ?3.12 | 89.61 | 6.61 |
| 3 | ?3.98 | 2.30 | 4.77 | 1250.13 | 11.53 |
| 4 | ?5.33 | 1.32 | 136 | 1.18 | 10.49 |
| 5 | 3.04 | 1.43 | 4.39 | 1244.29 | 10.05 |
| 6 | ?3.13 | 1.59 | 4.74 | 1188 | 3.38 |
| 7 | 6.81 | 5.98 | 2.04 | ?1.36 | 9102.35 |
| 8 | 5.62 | ?0.79 | 3.92 | 5.84 | 19788.26 |
| 9 | 83.94 | ?2.69 | 3.11 | 4.52 | 14846.82 |
| 10 | 73.62 | 5.24 | 3.74 | 4.32 | 3439.07 |
| 11 | 25.65 | 5.69 | 4.07 | 2.55 | 12619.25 |
| 12 | 69.44 | 5.78 | 2.97 | 1.85 | 1962.47 |
| 13 | 16.28 | ?5.23 | 3.53 | ?2.49 | 14322.48 |
| 14 | 73.66 | 2.54 | 3.41 | 4.39 | 13521.64 |
| 15 | 70.05 | ?5.14 | 3.87 | 5.81 | 4097.55 |
| 16 | 72.15 | ?3.01 | 3.85 | 5.68 | 691.74 |
| 17 | 6.03 | ?2.83 | 4.07 | 4.94 | 1051.51 |
| 18 | 7.84 | 5.45 | 1.96 | ?2.65 | 1938.96 |
| 19 | 16.62 | ?3.20 | 0.16 | ?5.23 | 651.06 |
| 20 | 5.43 | ?5.44 | 4.43 | 4.61 | 145.06 |
| 21 | 1.29 | ?6.00 | 0.95 | ?2.93 | 2163.14 |
| 22 | 18.72 | ?4.75 | 3.25 | ?1.87 | 649.81 |
| 23 | 5.57 | 6.26 | 2.09 | 3.51 | 1011.65 |
| 7.69 | | | | | |
| … | | | | | |
| 131 | 5.21 | ?4.80 | 1.49 | 5.64 | 3421.89 |
| 132 | 25.52 | 5.44 | 1.42 | 1.51 | 5236.92 |
| 35.11 | | | | | |

考虑到现场微震数据主要反映特定局部区域的损伤特征，为进一步表征深部岩体的微观结构非均匀性并真实再现其宏观渐进破坏过程，本研究引入Weibull分布来描述网格单元力学参数的统计特征。该方法建立在基于微震损伤变量计算的力学参数均值降低的基础之上（Tang 1997）：

(10) ??(??)=?????·(??/???)^(???1)·exp[?(??/???)^??]

该方程中，?? 表示岩体单元的力学参数，如弹性模量 E、单轴抗压强度 ???? 和抗拉强度 ????；??? 表示宏观平均值，m 为均匀性指标。基于Python的二次开发，本研究根据计算的损伤状况，将均匀性指标赋给左线时间滞后区的围岩。

如图11(c)(d)和表2所示，在该统计分布下，该软化区内围岩的力学性能表现出明显的非均匀程度，平均弹性模量为22 GPa，平均单轴抗压强度为70 MPa。该方法在避免常规模型中固有的网格畸变问题的同时，客观地恢复了岩体在经历前期开挖扰动后的真实非均匀初始物理状态。需要注意的是，基于Weibull分布的随机赋值仅应用于由微震事件空间分布及相应损伤反演预定义的弱化区域内。因此，随机种子并不决定重建弱化区的宏观位置或范围，而是控制该预定义区域内单元间的局部力学性能变化。不同的随机种子可能会改变局部弱单元的位置以及详细的裂纹起裂和扩展路径，而宏观弱化区域仍受现场微震信息的约束。

表2 基于Weibull分布的力学参数

| 参数符号 | 参数规格 | 值 |
|---------|---------|-----|
| ?? | 均值 | 8 |
| ˉ?? | 平均弹性模量 | 22/GPa |
| ˉ???? | 平均抗压强度 | 70/MPa |
| ˉ???? | 平均抗拉强度 | 7/MPa |

4.3 本构模型与参数选取

采用RHT（Riedel-Hiermaier-Thoma）损伤塑性模型来描述围岩。该强度模型由德国恩斯特·马赫研究所的Riedel、Hiermaier和Thoma开发，可有效捕捉高应变率动力冲击下脆性岩石的三个力学阶段：弹性阶段、线性硬化阶段和损伤软化阶段，如图13所示。围岩的宏观物理力学参数基于实验室试验和TSP超前地质预报确定，如表3所列。二次衬砌混凝土采用HJC（Johnson-Holmquist-Cook）本构模型表示。该模型是Johnson-Cook模型的增强迭代版本（Mazars, 1986, Mazars and Pijaudier-Cabot, 1989, Papa and Taliercio, 1996），专门用于表征混凝土在动力大变形下的屈服和压碎行为，其参数详见表4。为再现宏观剥落和破裂现象，添加了\*MAT\_ADD\_EROSION关键字，采用最大主应力和极限剪应变作为破坏和剥落准则。首先采用隐式静力求解器达到地应力平衡并建立初始损伤状态。随后修改材料参数，并将模拟切换为显式求解器以模拟双线隧道交替爆破冲击过程。

图13 RHT模型的应力-应变曲线

表3 RHT模型参数

| 参数符号 | 参数规格 | 值 | 参数符号 | 参数规格 | 值 |
|---------|---------|-----|---------|---------|-----|
| ??? | 质量密度 | 2665 kg/m3 | ??? | 拉应变率依赖指数 | 0.045 |
| ???? | 压碎压力 | 1.61E-2/GPa | ?????? | 参考压缩应变率 | 3.0E-8/ms?1 |
| ??????? | 压实压力 | 6/GPa | ??? | 多项式状态方程参数 | 1.2 |
| ?? | 孔隙率指数 | 3.0 | ??? | 多项式状态方程参数 | 1.2 |
| ??? | 初始孔隙率 | 1.0 | ??? | Hugoniot多项式系数 | 13.37/GPa |
| ??? | 多项式状态方程参数 | 13.37/GPa | ??? | Hugoniot多项式系数 | 19.68/GPa |
| ??? | 多项式状态方程参数 | 0/GPa | ?????? | 参考拉伸应变率 | 3.0E-9/ms?1 |
| ??? | Hugoniot多项式系数 | 8.55/GPa | ????? | 断裂压缩应变率 | 3.0E22/ms?1 |
| ???? | 抗压强度 | 106.16/MPa | ???? | 断裂拉伸应变率 | 3.0E22/ms?1 |
| ??\*? | 相对抗拉强度 | 0.169 | ??\*?? | 压缩屈服面参数 | 0.53 |
| ??\*?? | 相对剪切强度 | 0.18 | ??\*? | 拉伸屈服面参数 | 0.7 |
| ?? | 弹性剪切模量 | 23.4/GPa | ?? | 剪切模量衰减因子 | 0.5 |
| ?? | 破坏面参数A | 1.6 | ??? | 损伤参数 | 0.04 |
| ?? | 破坏面参数N | 0.61 | ??? | 损伤参数 | 1.0 |
| ??? | Lode角依赖因子 | 0.6805 | ???? | 最小损伤残余应变 | 0.01 |
| ?? | Lode角依赖因子 | 0.0105 | ???? | 残余面参数 | 1.6 |
| ???? | 压应变率依赖指数 | 0.043 | ???? | 残余面参数 | 0.61 |

表4 HJC模型参数

| 参数符号 | 参数规格 | 值 | 参数符号 | 参数规格 | 值 |
|---------|---------|-----|---------|---------|-----|
| ??? | 质量密度 | 2440 kg/m3 | ?? | 应变率系数 | 0.007 |
| ?? | 剪切模量 | 14.86/GPa | ??? | 压力常数 | 12/GPa |
| ?? | 高黏性强度系数 | 0.79 | ??? | 压力常数 | 16/GPa |
| ?? | 压缩硬化系数 | 1.6 | ??? | 压力常数 | 25/GPa |
| ?? | 压缩硬化指数 | 0.61 | ???? | 抗拉强度 | 6.7/MPa |
| ???? | 抗压强度 | 46.6/MPa | ???? | 抗压强度 | 46.6/MPa |
| ???? | 参考应变率 | 3.0E-9/ms?1 | ??? | 阻尼常数 | 0.04 |
| ????? | 最大等效应力 | 7.0 | ??? | 阻尼常数 | 1.0 |
| ????????? | 静水压力 | 22/MPa | ???? | 失效参数 | 0.8 |
| ????????? | 体积变形 | 0.001 | ??????? | 最小塑性应变 | 0.01 |
| ????c? | 体积变形 | 0.1 | ????c? | 静水压力 | 12/GPa |

4.4 动静耦合下的失稳响应与破坏结果

数值模拟准确描绘了远场动力荷载引起弱化围岩的渐进破坏过程。根据图14(a)中的损伤云图，在初始静力加载阶段（即步骤1），时间滞后区内的拱顶呈现初始软化状态，损伤变量D约为0.45，由微震反演确定。

图14 滞后区损伤演化云图 (a) 损伤变量演化 (b) 滞后区断裂破坏

在动力荷载响应阶段，左右隧道爆破应力波依次到达（分别为步骤2和步骤3），导致弱拱顶区内损伤范围的迅速扩展。因此，核心区内的构件遭受完全损伤，D值达到1。动力荷载的激发扰乱了初始结构平衡。在围岩持续的应力调整和变形之后，系统最终发生宏观坍塌，即步骤4。

图14(b)中最终模拟的破坏模式与现场观察到的主要空间损伤特征基本一致。在受损围岩深部形成了一个明显的裂纹网络，而刚性衬砌的破坏集中在自由表面，呈现出特征性的大规模剥落形态。

恢复模型爆破后的位移、振速和有效应力场以阐明动力触发过程，如图15所示。由于局部开挖引起的损伤和软化区域的存在，系统表现出明显的非对称力学性能。右隧道的应力分布（代表未软化区域）相当均匀。左隧道中，高刚度衬砌与受损围岩之间产生了显著的刚度差异。高地应力绕过低模量软化区域，引发应力重分布，导致围岩处于高度松弛的卸荷状态，从而增强了其变形动能。因此，左线拱顶的位移场和最大振速明显大于右线。爆破应力波充当触发机制，激活储存在承载拱中的残余弹性能，并以动能的形式直接作用于衬砌结构，无任何缓冲反力。

图15 动力响应云图

此外，图16所示的静水压力分布清晰地对比了两种不同边界条件下衬砌的力学性能。在完好的隧道中，无软化损伤，围岩与衬砌表现出优越的协同承载性能，使整体衬砌结构维持相当均匀的压缩状态。相比之下，由于爆破坑洞和软化造成的空隙，左隧道的协同承载机制受到破坏。因此，拱顶线缺乏来自围岩的一致性径向支撑。刚度差异给左线衬砌施加了极度危险的非对称加载条件，即"外压内拉"。这种显著的力学差异不仅印证了无空洞环境下的标准压缩力学行为，也与He等人（2021）提出的高地应力下空洞衬砌的破坏机制高度吻合。

图16 压力分布图 (a) 左线爆破后演化 (b) 右线爆破后演化

4.5 能量耗散与时间滞后致灾机制

基于数值模拟揭示的损伤演化规律，本次时间滞后型岩爆事件的能量驱动机制可概括如图17所示。在静力加载期间（即图17(a)中的OB段），输入系统的外部功主要转化为弹性应变能???，储存在深部完好岩体中。随着应力缓慢重分布，微裂纹的大量萌生导致耗散能????呈非线性增加，从而持续降低储能上限。与纯静力流变过程相比，交替爆破产生的循环动力荷载显著加速了耗散能的比例，如图17(b)所示。

图17 能量转化与耗散机制 (a) 沿应力-应变曲线演化 (b) 应变能与耗散能演化

由于拱顶处存在超挖空洞或软弱夹层，围岩具有充足的变形空间以供长期流变行为。当承载拱内积累的弹性应变能接近临界损伤状态（即峰值点C）时，邻近隧道的一次远场爆破充当最终触发驱动，使系统走向破坏。被激活的高势能瞬时绕过峰值后软化阶段，转化为剥落岩块的动能。这种瞬时能量释放最终导致衬砌结构在长时间滞后后发生灾难性剥落。

5. 讨论

5.1 时间滞后型岩爆致灾机制讨论

传统上，深埋隧道中的时间滞后型岩爆大多归因于高地应力环境下的静力岩体流变或卸荷后的缓慢应力重分布。然而，本研究通过整合现场微震测量和数值模拟，证实了前期爆破损伤引起的局部围岩非均匀软化是触发此类特殊岩爆的关键隐含边界条件。在工程实践中，虽然爆破坑洞用水泥混凝土回填，但其周围围岩的力学性能已发生不可逆退化。施加高刚度二次衬砌后，围岩-回填-衬砌系统内部出现了显著的刚度不匹配。二次衬砌的强约束阻碍了深部岩体的径向应力释放路径，导致应变能在软化区外围异常积累。

在此背景下，邻近隧道的交替爆破已不再仅仅是简单的振动扰动，而是演变为一种低频疲劳触发布载。它不仅加剧了弱区内的损伤演化，还在积累能量突破接触面临界强度时瞬时激活系统的残余弹性能，随后触发严重动力失稳。这一发现有效弥补了传统静力蠕变理论在解释动静耦合触发岩爆方面的局限性。

5.2 累积损伤的等效处理与数值模型优势

针对深部长距离双线隧道中多循环爆破模拟的计算瓶颈问题，本文提出了一种将微震损伤反演重建与动力触发相结合的计算方法。在显式动力框架下，对数百次历史爆破进行全尺度模拟极易导致严重的网格畸变和计算发散。因此，本研究摒弃了预设静力几何缺陷的传统方法，而是通过将Weibull分布与微震能量等效原理相结合，将前期累积开挖损伤转化为空间力学参数的非均匀降低。

这种数据驱动的参数弱化方法为破坏前宏观非均匀弱化状态提供了等效表示。通过仅显式加载引发岩爆的关键强扰动，该模型避免了数值奇异问题，同时有效解耦了两个主要致灾因素：内部疲劳损伤和外部动力触发。这验证了所提方法在复杂岩爆反演中的高效性和可靠性。

5.3.

工程影响与灾害防治策略基于上述时延型破坏机理，本研究针对高应力围岩条件下深埋钻爆施工提出以下针对性的岩爆防治与防控策略：（1）弱区刚度补偿机制：数值分析结果表明，根据地质调查、超挖记录以及微震监测所识别出的明显弱化区域，在二次衬砌施工前实施必要的加固措施具有合理的力学依据。可考虑采用深孔灌浆作为可选加固措施，以提高弱化围岩的整体刚度，恢复其宏观等效模量，从而减少岩体-支护体系内的局部刚度不连续性，改善其协同承载性能。（2）多参数联合预警：在常规微震监测中，应对地震波b值与空间分形维数D值的同步演化保持高度警惕。当这两个分别表征裂隙规模汇合与空间维度降低及聚集的指标在短时间内同步突破经验阈值，例如b值降至0.7以下、D值降至1.7以下时，即使当前支护表面看似稳定，也应立即启动防灾应急响应，以防范衬砌背后可能发生的能量释放。（3）爆破时空动态优化：双线隧道平行施工期间，应动态调整隧道掌子面之间的轴向掘进间距与起爆间隔。建议依据微震能量衰减曲线划定扰动敏感区。在已识别的弱区附近，严禁左右掌子面频繁交替起爆，以降低岩体的疲劳加载速率。此外，双线隧道开挖的合理掌子面错距与起爆间隔仍需结合现场爆破振动、微震活动及围岩条件进行进一步的工程标定。6. 研究局限与未来展望尽管本研究尝试将现场微震信息引入数值模型，并通过微震能量反演与Weibull分布重建了围岩非均匀预破坏损伤状态，但仍存在若干不确定性。首先，损伤变量的反演依赖于地震效率参数，该参数系基于现场定点爆破统计结果确定，该参数的变化将直接影响实际释放能量的估计及相应损伤程度的评估。其次，微震监测中的震源定位误差不可避免，当离散微震事件映射至数值单元并用于折减材料参数时，这些空间误差可能进一步影响对局部弱化区域位置、范围及损伤程度的预测。因此，本研究重构的损伤场更应被理解为对围岩破坏前宏观空间弱化特征的表征。此外，本研究采用连续介质模型表征围岩，无法完全捕捉天然节理、裂隙及其扩展和贯通对应力波传播、散射和耗散的影响。对于隧道的长期运营阶段，列车振动等长时程、低幅值循环动荷载对既有软化区域的持续损伤效应仍需进一步研究。本方法主要在本工程的岩性、原岩应力及微震监测条件下进行评估，其在不同岩性、地质构造和应力环境下的适用性仍有待通过更多工程案例加以验证。未来工作将结合不同地质条件下的现场监测数据，并进一步引入连续-非连续耦合建模方法，以提高复杂裂隙岩体损伤重构与动力失稳分析的可靠性。7. 结论本研究将现场微震监测与非均匀动力数值模拟充分结合，深入研究了超叠加动力与静力荷载作用下局部软化围岩的致灾机理，并以深埋双线高原隧道拱顶典型时延型岩爆为研究基础，得出以下主要结论：（1）揭示了时延型岩爆前兆演化特征：根据长期微震监测，时延型岩爆的发生与裂隙规模及空间分布的极端演化密切相关。表征裂隙面空间维度降低与聚集的D值和表征微裂隙大规模汇合的b值在灾害前夕均出现了极为同步的骤降。当b值降至0.7以下且D值降至1.7以下时，岩体系统逼近动力失稳临界状态。这些前兆特征的发现为高刚度支护背后隐蔽性岩爆的量化预警准则的建立提供了可能。（2）构建了微震驱动的初始状态重构数值技术：通过将Weibull统计分布与基于微震能量的损伤变量等效反演相结合，实现了对局部软化区域围岩力学参数的动力非均匀折减。该方法克服了传统几何缺陷模型无法还原岩体原始物理损伤状态的不足，忠实再现了拱顶岩爆的完整动力破坏过程，同时避免了网格畸变。（3）阐明了刚度失配主导的非对称加载致灾机理：由于拱顶软化及前期爆破形成的空隙，二次衬砌无法获得围岩均匀的径向反力，导致下部岩体与衬砌之间产生显著刚度失配。邻近隧道的爆破动力荷载充当触发开关，激活了深部残余弹性能，造成拱顶部位出现"外压内拉"的严重局部应力集中，最终导致高刚度衬砌的宏观剥落破坏。另一方面，在均匀且无损伤的情况下，衬砌结构则保持整体受压状态。（4）提出了聚焦现场的防治体系：对于超挖或大面积损伤区域中出现的变形稳定假象，必须保持警惕。工程实践中应充分融入多参数联合预警机制，在弱区实施主动灌浆加固以避免承载刚度故障。此外，必须对双线平行掘进的掘进间距进行动态优化，以从源头上切断动力荷载的触发路径。作者贡献声明Liu Yumin：原稿撰写、形式化分析、数据整理。Ma Tianhui：方法、调查研究。Dai Le：审阅与编辑、概念设计。Zhao Jiaxu：审阅与编辑。资金支持本研究受国家自然科学基金资助（基金编号：42272314）。Liu Yumin｜Ma Tianhui｜Dai Le｜Zhao Jiaxu