CN121328362A - 隧洞围岩稳定性预测方法 - Google Patents

Abstract

本发明涉及围岩稳定性预测，为了实现深埋条件下对于围岩破裂的超前预测，提供了隧洞围岩稳定性预测方法，通过采用有限元与离散元的耦合，能够同时考虑围岩的整体变形和微观裂隙扩展，提高了围岩稳定性分析的准确性及计算效率；双通道超图神经ODE破裂模式预测模型能够提前3‑7天给出破裂体积及失稳概率，同时考虑了物理一致性问题，即能量、动量残差实时收敛，保证与预测结果的准确性；通过有限元‑离散元耦合模型获取实时监测数据并基于双通道超图神经ODE破裂模式预测模型进行稳定性预测，无需人工调整参数，系统自动完成数据分析‑预测‑预警闭环；实现了把深埋隧洞围岩稳定性分析从静态的、经验的分析升级为实时的、动态的超前的预测。

Description

隧洞围岩稳定性预测方法

技术领域

本发明涉及围岩稳定性预测，具体是一种隧洞围岩稳定性预测方法。

背景技术

围岩稳定性分析一直是地下空间工程设计和施工过程中的重要研究内容。随着交通、水电、矿山等工程向深部发展，隧洞埋深普遍超过千米，围岩处于高地应力（>20MPa）、高渗透压、强卸荷及复杂构造应力场等多场耦合环境中。通过进行围岩稳定性分析可以明确隧洞建设的安全性，进而保障隧洞的安全施工。通常地，围岩稳定性分析方法主要有现场试验、物理试验、理论分析和数值模拟，但受制于深部岩体特殊的动态力学特性，均存在着各自的缺陷。

深部岩体取样困难，现场原位试验（如地应力测试、水压致裂）成本高、周期长，且受开挖扰动影响导致测试结果离散性大，难以反映真实应力路径。物理试验则受制于尺寸效应（缩尺模型无法还原岩体节理网络及应力梯度）及加载能力（现有真三轴设备难以模拟千米级埋深应力场），试验结果与工程实际偏差较大。

在理论解析方面，深部岩体呈现明显的非线性、非连续及各向异性特征，现有解析解（如Kirsch方程、应变软化模型）均基于理想弹塑性假设，无法考虑裂隙动态扩展及应力重分布的时序演化。

在数值模拟方面，通过软件进行数值模拟已经成为工程实践中重要的研究方法，传统有限元法（FEM）虽然能够较好地模拟围岩的整体变形和应力分布，但难以捕捉岩体内部微裂隙的萌生-扩展-贯通过程；离散元法（DEM）虽可揭示岩桥断裂、块体滑移等细观破坏机制，但计算效率随颗粒数量呈指数下降，难以处理大范围的连续介质问题。此外，现有方法普遍采用"先验参数+离线计算"模式，同时未有效融合施工期监测数据（如多点位移计、微震信号），模型参数在施工过程中不再更新，无法反映围岩动态响应，导致计算分析结果与实际围岩变形/破坏事件之间存在3~7天时间滞后，未能进行有效防控。

发明内容

为了实现深埋条件下对于围岩破裂的超前预测，本发明提供了一种隧洞围岩稳定性预测方法。

本发明解决上述问题所采用的技术方案是：

隧洞围岩稳定性预测方法，包括：

步骤1：构建隧洞FDEM模型，FDEM模型为有限元-离散元耦合模型；

步骤2：基于数据采集装置采集隧洞现场的监测数据；

步骤3：基于监测数据对有限元-离散元耦合模型进行参数更新并获取正演结果；

步骤4：基于监测数据及正演结果建立双通道超图神经ODE破裂模式预测模型；

步骤5：基于双通道超图神经ODE破裂模式预测模型对围岩潜在破裂体积及失稳概率进行预测。

进一步地，步骤1包括：

步骤11：以隧洞中心线为基准建立深埋隧洞模型；采用同心圆分区方式将隧洞模型进行分区，依据与洞壁的距离，依次划分为近场连续区、远场离散区及无反射吸收区；

步骤12：根据现场资料创建节理网络并赋值，然后将节理网络映射到远场离散区；

步骤13：对近场连续区划分有限元网格，远场离散区填充离散颗粒，并设置过渡层共享节点；

步骤14：在近场连续区与远场离散区交界面插入界面过渡单元；

步骤15：参数赋值，包括：为FEM区设置岩石物理力学参数和应变软化模型；为DEM区配置平行黏结模型和初始地应力；使用Karhunen-Loève展开生成随机场并插值；设置时间步长并进行DEM到FEM的映射；利用Python+Gmsh进行自动生成网格，并设置监测数据更新的接口。

进一步地，数据采集装置包括：微震传感器、 声发射传感器、分布式光纤及多点位移计。

进一步地，步骤3基于内、外两个时间尺度进行参数更新，外循环每更新一次岩体宏观参数场；内循环在内按频率更新掌子面前方破裂面黏聚力与摩擦角。

进一步地，外循环使用集合卡尔曼滤波更新参数，内循环使用无迹卡尔曼滤波更新参数。

进一步地，步骤3还包括计算全局位移残差及UKF协方差，若不满足全局位移残差<1mm且UKF协方差< 0.01，则缩短。

进一步地，双通道超图神经ODE破裂模式预测模型实现步骤为：

超图构建：将有限元单元与离散颗粒全部视为节点，采用三类超边分别捕捉地层结构、微震事件簇与位移梯度簇；

双通道超图神经ODE建模：将每个岩体节点分裂成两条通道：完整岩体通道和已损伤岩体通道，根据超图构建的结果，采用超图聚合消息驱动ODE右侧，让完整-损伤两通道差异直接控制损伤速率，同时进行物理约束损失设计对ODE右侧实时校正，输出破裂体积及对应的失稳概率。

进一步地，还包括：基于FDEM模型正演结果对双通道超图神经ODE破裂模式预测模型进行重训练。

进一步地，重训练时，训练版本号自增并存档。

进一步地，还包括步骤6：基于预测结果进行可视化显示及预警。

本发明相比于现有技术具有的有益效果是：通过采用有限元与离散元的耦合，能够同时考虑围岩的整体变形和微观裂隙扩展，提高了围岩稳定性分析的准确性，避免了单一离散元方法在大范围计算中的低效率问题，提高了计算效率；双通道超图神经ODE破裂模式预测模型能够提前3-7天给出破裂体积及失稳概率，同时考虑了物理一致性问题，即能量、动量残差实时收敛，保证与预测结果的准确性，通过完整-损伤双通道将损伤速率提前写进ODE右侧，再使用超图把潜在破裂面提前建模，伴随积分一次就能看到3-7天后的轨迹，实现了超前预测；通过有限元-离散元耦合模型获取实时监测数据并基于双通道超图神经ODE破裂模式预测模型进行稳定性预测，无需人工调整参数，系统自动完成数据分析-预测-预警闭环；实现了把深埋隧洞围岩稳定性分析从静态的、经验的分析升级为实时的、动态的超前的数字孪生系统，可适用于不同地质条件和隧洞规模的深埋隧洞围岩稳定性预测。

附图说明

图1为隧洞围岩稳定性预测方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图1所示，隧洞围岩稳定性预测方法，包括：

步骤1：构建隧洞FDEM模型，FDEM模型为有限元-离散元耦合模型。具体包括以下步骤：

步骤11：以隧洞中心线为基准建立深埋隧洞模型；采用同心圆分区方式将隧洞模型进行分区，依据与洞壁的距离，依次划分为近场连续区、远场离散区及无反射吸收区。

以隧洞中心线为基准，在Abaqus中建立深埋隧洞模型（直径D，长度L≥30D），采用“同心圆”分区的方式进行分区：FEM区（距洞壁0~5D的“近场”连续区）；DEM区（距洞壁5D~20D的“远场”潜在破裂区）；ABS区（为人工阻尼边界区，最外层 20D~25D）；其中0-5D区域洞壁应力集中但尚未出现贯通裂缝，变形以连续介质为主，因此采用纯有限元（FEM）手段；5D-20D区域由于高地应力下围岩最易出现微裂隙，若发生节理扩展会导致块体滑移，需采用离散元（DEM）手段才能模拟开裂、滑移等非连续现象；20D-50D区域由于外应力扰动可忽略，但边界必须吸收波动防止反射，因此设置ABS区，ABS区为区别于FEM区的无影响有限元区域，仅用于吸收反射波，不参与围岩力学响应计算，其材料参数设定为弹性模量逐层递减的黏弹性介质，以避免应力波反射干扰。

步骤12：根据现场资料创建节理网络并赋值，然后将节理网络映射到远场离散区。

将现场资料中的地质强度指标（GSI）、岩石质量指标（RQD）、节理产状统计导入Python中生成三维Voronoi节理网络，再对整个岩体进行宏观参数赋值，利用Monte-Carlo算法给每条节理赋予随机的粘聚力c、内摩擦角、法向刚度Kn、切向刚度Ks，最后将节理网络映射到DEM区，作为后续离散元的初始颗粒黏结位置。该映射过程采用双向哈希索引结构，保留Voronoi节理面ID与DEM颗粒键之间的一一对应关系，确保在后续破裂面演化中可反演节理编号与空间坐标。

步骤13：对近场连续区划分有限元网格，远场离散区填充离散颗粒，并设置过渡层共享节点。

1）在FEM区进行有限元网格划分，有限元网格采用8节点六面体减缩积分（C3D8R）单元，洞壁附近0.1D，向外按1.3倍递增，在FEM区与DEM区交界面预留1层共节点单元。

2）DEM区使用半径膨胀法生成球体/多面体颗粒，确保局部孔隙率<3%，颗粒半径r满足r≤0.05D，能分辨节理厚度；在DEM区与FEM区交界面，强制在1.5r厚度内生成共节点颗粒，即这些颗粒圆心落在有限元节点坐标上。

步骤14：在近场连续区与远场离散区交界面插入界面过渡单元。

以共节点映射的方式，用KD-Tree算法把有限元节点与离散颗粒中心做最近邻匹配，在交界面生成0.3 m厚的界面过渡单元（Interface Element），采用Goodman节理单元+损伤演化的本构关系，设初始法向刚度，其中为FEM区材料的弹性模量与单元特征长度之比，反映FEM区单位厚度轴向刚度，为DEM区颗粒集合体的等效弹性模量与平均粒径之比，反映DEM区单位厚度轴向刚度，取较小值，防止过渡层刚度过大导致应力集中或数值反射。切向刚度Ks=0.4Kn，当界面位移>0.5 mm时，自动退化为无厚度接触，允许滑移或开裂。

步骤15：参数赋值，包括：为FEM区设置岩石物理力学参数和应变软化模型；为DEM区配置平行黏结模型和初始地应力；使用Karhunen-Loève展开生成随机场并插值；设置时间步长并进行DEM到FEM的映射；利用Python+Gmsh进行自动生成网格，并设置监测数据更新的接口。

1）FEM区赋值岩石物理力学参数如弹性模量E、泊松比ν，采用应变软化模型，软化模量由现场三轴试验曲线标定；

2）DEM区内颗粒之间采用平行黏结模型（PBM），法向黏结强度、切向黏结强度通过Voronoi单元体积加权，从节理网络统计值映射，初始地应力场通过计算，其中为侧压力系数，是由水压致裂试验得出的，为垂直应力，由=γh计算得到(γ为岩体平均容重，h为埋深)；

3）使用Karhunen-Loève展开生成随机场（相关长度=2D）；将随机场插值到有限元积分点和离散黏结键上，生成岩体宏观物理力学参数的初始空间随机场样本，包括弹性模量E、泊松比v、侧压力系数等；

4）设置时间步长，有限元步长，离散元步长，并且每200步离散元进行一次粗化映射，把颗粒合力、合力矩传递给有限元节点。

5）利用Python+Gmsh进行自动生成网格，并设置监测数据更新的接口，确保后续监测数据可直接插值到节点或颗粒。

步骤2：基于数据采集装置采集隧洞现场的监测数据。

在深埋隧洞开挖与运行阶段，沿着隧洞轴线及典型断面布设传感器、多点位移计、应力计等装置进行数据的采集与保存，建立“空间-时间-频率”多源监测数据流接口。“空间-时间-频率”多源监测数据流包括：

（1）空间：微震（MS）32通道、声发射（AE）16通道、分布式光纤（DAS）2km及多点位移计（MPBX）12点，形成“环向+纵向”双环监测网；

（2）时间：1Hz~1 kHz实时采样，微震/AE 1 kHz、DAS 10 kHz、MPBX 1 Hz；

（3）频率：提取微震b值、位移速率谱指数β_disp、频散曲线等20维特征向量。

步骤3：基于监测数据对有限元-离散元耦合模型进行参数更新并获取正演结果。

本实施例在内、外两个时间尺度上，通过设计双循环同化算法实现监测数据驱动FDEM模型更新。外循环采用基于自适应协方差膨胀的集合卡尔曼滤波（EnKF-AE），每（如6h）更新一次岩体宏观参数场；内循环则在内，使用无迹卡尔曼滤波（UKF）对破裂面内聚力c、摩擦角进行高频（如5min）修正。具体步骤如下：

（1）首先根据钻孔信息及节理网络统计强度分别规定岩体物理力学初始参数和破裂面强度参数。对第i个成员，其初始参数集合记为：，其中为宏观岩体参数向量，为破裂面强度参数向量，N()为正态分布，且宏观岩体参数向量先验均值、破裂面强度参数向量先验均值及宏观岩体参数向量协方差Σθ、破裂面强度参数向量先验均值Σφ由现场试验和先验统计给定，为集合最大成员。

（2）在外循环中使用基于自适应协方差膨胀的集合卡尔曼滤波（EnKF-AE）每更新一次岩体宏观参数场，具体步骤如下：

1）对每个集合成员运行FDEM正演6h，获得岩体宏观参数场的先验估计值 ，，其中为 FDEM正演算子。

2）从“空间-时间-频率”接口读取6 h累积观测数据，然后基于先验估计值提取与现场监测同量纲、同位置的虚拟观测向量，例如洞壁收敛、微震能量。，其中为观测算子的函数形式，用于提取位移、微震能量等。

3）设置协方差膨胀系数计算集合成员之间的观测误差协方差矩阵，在本实施例中，避免滤波发散。 ，其中为集合虚拟观测的均值向量，，为矩阵转置符号。

4）计算卡尔曼增益，用于量化模型误差与观测误差的权重；矩阵越大，越信任观测。，其中H为观测算子矩阵形式，R为仪器噪声协方差，由厂家标定给出（位移计±0.1 mm，微震能量±5 %）。

5）把观测残差 按增益折回岩体参数空间，完成一次模型校正，获得第i个集合成员在k时刻更新后的岩体宏观参数向量，，为k时刻现场实测观测值（位移、微震能量等）。

（3）内循环在内使用无迹卡尔曼滤波按频率对破裂面内聚力及摩擦角进行修正。具体步骤如下：

1）以掌子面前30 m为窗口，提取所有破裂面单元，设置状态向量。状态向量定义为k时刻掌子面前方30 m窗口内所有已激活破裂面的等效粘聚力c与内摩擦角，其初始值由节理网络统计强度按Voronoi单元体积加权平均得到。，其中，c为粘聚力，为内摩擦角。

2）设置Sigma点，维数为n，用（2n+1）个Sigma点捕获状态分布的均值和协方差。设置中心点用于计算均值权重和协方差权重，与其他2n个Sigma点计算预测均值权重和协方差权重。

，

，

，

，

，

，

其中，为调节参数，用于控制Sigma点与均值的距离，为调节参数转成实际空间步长的缩放系数，为第j个正向偏离点（沿协方差主轴+γ方向），为第j个负向偏离点（沿协方差主轴-γ方向），为Cholesky分解协方差矩阵的列向量。

3）对每个Sigma点进行5 min快速FDEM子循环，可更新循环状态并提取虚拟观测值，计算该sigma点在k+1时刻的先验均值和先验协方差。

，

，

，

，

其中，为在第k时刻生成的第j个Sigma点，为5 min迷你FDEM子模型，提取位移、能量、主频，Q为过程噪声协方差，，其中diag()为对角矩阵符号，为黏聚力c的噪声方差，为摩擦角的噪声方差。

4）计算卡尔曼增益，将预测值转为最优估计，对状态向量与协方差进行最优更新。

，

，

，

其中，为k+1时刻的实测观测值（位移、微震能量等），为同一时刻的模型预测观测值（由先验状态经观测算子得到），为状态-观测互协方差，为互协方差，为矩阵转置符号。

（4）把外循环得到的宏观参数和内循环得到的局部强度写回到FDEM模型的有限元网格与离散元颗粒中，并保证能量守恒、坐标一致。

进一步地，系统把最新FEM/DEM正演结果与现场多点位移计、光纤实测值实时对比，得到全局位移残差，同时读取UKF输出的协方差矩阵P，计算其迹Tr(P)，用于量化当前参数不确定度；若全局位移残差< 1 mm且UKF协方差Tr(P) < 0.01时保持参数，否则标记为“强扰动”，下一次外环缩短ΔT = 3 h。

现有技术在工程应用上往往只采用传统的静态信息处理手段，即通过钻孔取芯、三轴压缩试验获得岩石的物理力学参数输入至模型中，在施工期间不再更新数值模型，而且如果需要数据校正，需要手动调参，数据的不确定性较强，调参对象只针对弹性模量E。而本发明所设计的双循环同化算法利用多源信息融合处理技术进行监测数据驱动更新数值模拟模型，采用位移、微震等监测数据进行有限元与离散元模型的自主实时更新，无须手动调参，外循环(EnKF-AE）将宏观岩体参数转变为随机场的形式，每6h进行自主校正，内环（UKF）将破裂面的粘聚力c、内摩擦角φ进行5min高频修正，能够30天累计偏差<1mm，而传统静态模型同期＞10mm。

步骤4：基于监测数据及正演结果建立双通道超图神经ODE破裂模式预测模型。

双通道超图神经ODE破裂模式预测模型实现步骤为：

（1）超图构建。设置时间窗，为当前外循环时刻（单位 h），将FEM单元与DEM颗粒全部视为节点，分别采用三类异构超边捕捉地层结构、微震事件簇与位移梯度簇。

1）将同一岩层或同一节理组内的所有节点拉成一条高阶边，天然容纳层间滑移。根据勘察节理产状与地质强度指标GSI 分区，建立地质层面方程：；其中为第m个地质层面（岩层或节理面）的编号，用于区分不同层位；为该层面的单位法向量，指向层面上方；x为空间任意一点的坐标，用来计算它到层面的相对位置，为层面上已知参考点的坐标，用于锁定层面在三维空间中的具体位置；M为层面序号，从1到M逐个编号，覆盖模型区域内所有主要岩层与节理组。

对任意节点i，计算到面的距离。若，则将i归入集合；，δ为节点到岩层/节理面的距离容差，反映岩层厚度识别精度，本实施例取0.5m。

每条结构超边即，其中，为被第m个岩层（或节理面）捕获的所有节点集合，为结构超边权重，在1基础上累加边内平均微震能量密度，用于突出高能岩层的相互作用强度，为节点i在内累积的微震能量密度。

2）把同一破裂事件或应力区内的节点连成高阶边，捕捉损伤局部化。从监测流提取内所有微震事件，其中为第s个微震事件的空间坐标，表示破裂源位置，为事件发生时刻，用于时间窗筛选与排序，s为事件序号，从1到S逐个编号，S为当前窗内有效事件总数，为事件释放能量，决定该事件在超边权重中的贡献大小，并在3D空间构建Delaunay单纯形，即：。对任意单纯形，若则保留该单纯形，同时将单纯形顶点扩展到最近的一组FEM/DEM节点，为微震事件能量阈值，本实施例取，单纯形内平均能量超过该值才保留，用于筛掉噪声事件。

每条事件超边即，其中是以微震单纯形为中心、事件半径内所有FEM/DEM节点的集合，事件超边权重，在1基础上累加对数能量，使高能破裂簇获得更大传播权重。

3）把位移梯度相似的节点聚成一条超边，捕捉连续变形带。由多点位移计与光纤数据插值得到节点级位移向量，，其中x_i为节点i的空间坐标。

计算二阶位移梯度张量（采用球邻域Ωi最小二乘拟合），并提取不变量特征，然后对所有节点特征执行 K-Means（K=64），得到簇标签，用表示簇编号，取值范围为1，…，K。

，其中是对节点i的空间位置处的位移求偏导数的矩阵。

，

其中，为梯度整体大小（Frobenius 范数），越大表示变形越剧烈，为梯度行列式，反映体积变化趋势（其值>0，表示膨胀，其值≈0表示剪切，其值<0表示压缩），为梯度偏量范数，剔除体积变化后纯剪切变形的大小。

每条梯度超边即。

4）超图组装。合并全部超边，得到超边集、对应关联矩阵、超边权重对角矩阵。

，

，

。

其中，∪表示并集，三者无重复、全覆盖，组成当前时刻的完整超图边集，E×V为尺寸（E条超边，V个节点），为判断函数，取值只有0或1：若节点i属于超边e，则=1；否则=0，diag()为对角矩阵符号，w()为对角线上依次存放每条超边…的权重, 为维度（E=超边总数）。

（2）双通道超图神经ODE建模。通过节点初始观测特征，设置双初始态，即通道A（完整岩体），通道B（已损伤岩体）。

，其中，为宏观岩体参数随机场在节点i的值（包括弹性模量E、泊松比、侧压力系数等子向量），为破裂面强度随机场在节点 i 的值（节理黏聚力c、摩擦角φ等子向量），为节点i当前应力张量分量（通常取3个正应力+ 3个剪应力，或不变量形式），为节点i当前应变张量分量（与对应，可计算弹性能密度），为微震能量密度插值到节点i的标量，反映局部破裂活跃度，为位移梯度张量的Frobenius范数，用于量化节点邻域变形剧烈程度。

，

其中，为两层 MLP（含 LayerNorm & ReLU）。

对任意超边e收集聚合消息向量：

，

其中为超边e的聚合消息向量，由边内所有节点隐状态平均并门控后得到，用于后续向节点传递，为小型MLP，输入为节点相对坐标偏移与边内平均微震能量密度，输出门控向量，为Hadamard积，实现门控式加权，。

节点i收到所有含i的超边消息：，其中，e为包含节点i的超边编号，决定消息来源范围，为超边e的权重，反映该边内微震能量或变形剧烈程度，能量越高权重越大。

将聚合消息拆回A、B通道，随后分别独立更新通道A、通道B，。

定义隐空间差异范数，节点级损伤变量，满足初值，并服从演化。当双通道差异超过阈值β时损伤快速累积，并进行及时积分。

，

，

其中，为损伤速率增益系数；σ()为ReLU激活函数，即σ(x)=max(0,x)，确保损伤仅在双通道差异超过阈值β时累积，为节点i所代表的物理体积（m³），FEM单元即单元体积，DEM颗粒即颗粒体积，为整个模型在时刻t的累计破裂体积，用于衡量围岩损伤程度并计算失稳概率P(f)。

（3）物理约束损失（PIL）设计。将能量守恒、动量平衡、本构关系等物理规律嵌入损失函数，约束模型预测结果符合物理实际，提升可信度与泛化能力。

为避免额外损失权重调参，把能量残差与动量残差写成右侧修正项。

，

其中，为外力在 t 时刻对岩体做的总功率积分（输入能量），由边界位移与外力积分计算；为岩体内部弹性能增量，由应力-应变计算弹性能，为岩体整体动能增量，由节点加速度估计动能；为断裂能（单位：J/m²），由预测破裂区域体积估算；为材料断裂能（常数），由节点级损伤变量Di(t) 与单元表面积加权估算。

，

右侧修正（仅作用于A通道，保证完整侧更贴近物理），

，

其中，为节点i处应力散度，由邻域应力差分近似得到，反映内力分布；为体力密度，通常取重力，方向竖直向下；为节点i的集中质量，由单元或颗粒体积与密度换算得到；为节点i的加速度向量，由位移二阶差分或速度导数实时计算；为邻居j与中心节点i的应力张量差，提供内力梯度；为两节点间距，用于差分尺度归一化；从i指向j的单位方向向量，决定差分符号与投影方向；为能量残差修正强度系数（无量纲），初始值0.01可调，为能量残差的符号函数，+1表示系统总能量过剩（外力做功 >内部耗散+弹性能），系统总能量-1表示亏缺，0表示平衡；为能量残差对节点i的A通道隐状态梯度。

步骤5：基于双通道超图神经ODE破裂模式预测模型对围岩潜在破裂体积及失稳概率进行预测。

进一步地，还包括步骤6：基于预测结果进行可视化显示及预警。

基于预测结果绘制破裂体积-失稳概率曲线，通过破裂体积V(f)去对应曲线找对应的失稳概率，得到每时刻的失稳概率P(f)。失稳概率P(f)≥0.7，触发红色预警，自动推送加固方案；0.5≤失稳概率P(f)<0.7，黄色提示，增加监测频次；失稳概率P(f)< 0.5，绿色通行。系统同时给出95 %置信带，供决策层参考。

为了提高模型预测准确性，还包括对模型进行更新；包括以下方式：

（1）每日凌晨自动拉取前24h新数据，对双通道ODE的映射网络进行10个epoch的增量微调，学习率用cosine退火。只更新预测头与门控参数，冻结大部分卷积与注意力层，防止灾难性遗忘；版本号自增并存档，支持一键回滚至任意历史版本。版本号为每次预测模型所对应的编号，通过版本号，可以保证同一模型骨架下，每次训练结果都有唯一、递增、不可重复的标签，避免覆盖旧权重，方便回滚或对比。

（2）若连续三天全局位移残差且UKF协方差 Tr(P) < 0.01，则维持当前模型；否则再次触发双循环同化，重新校正参数后再训练。

（3）外循环结束后，进行一致性校验，即以相同初始条件分别运行FDEM正演与双通道超图神经ODE破裂模式预测模型，若二者破裂体积相对误差>10%，则触发超图重训练（增量微调20 epoch），确保预测模型与物理模型偏差<5%。

通过监测数据驱动模型参数更新，使得模型参数与真实岩体保持一致性，提高了有限元-离散元耦合模型正演结果及双通道超图神经ODE破裂模式预测模型的预测结果的准确性，采用正演结果修正预测模型，进一步地提高了预测准确性。采用本发明的方法可以使破裂预测提前3-7天、位置误差< 1.5 m、时间误差< 12 h，实现真正意义上的实时、动态、超前围岩稳定性分析。

Claims (10)

1.隧洞围岩稳定性预测方法，其特征在于，包括：

步骤1：构建隧洞FDEM模型，FDEM模型为有限元-离散元耦合模型；

步骤2：基于数据采集装置采集隧洞现场的监测数据；

步骤3：基于监测数据对有限元-离散元耦合模型进行参数更新并获取正演结果；

步骤4：基于监测数据及正演结果建立双通道超图神经ODE破裂模式预测模型；

步骤5：基于双通道超图神经ODE破裂模式预测模型对围岩潜在破裂体积及失稳概率进行预测。

2.根据权利要求1所述的隧洞围岩稳定性预测方法，其特征在于，步骤1包括：

步骤11：以隧洞中心线为基准建立深埋隧洞模型；采用同心圆分区方式将隧洞模型进行分区，依据与洞壁的距离，依次划分为近场连续区、远场离散区及无反射吸收区；

步骤12：根据现场资料创建节理网络并赋值，然后将节理网络映射到远场离散区；

步骤13：对近场连续区划分有限元网格，远场离散区填充离散颗粒，并设置过渡层共享节点；

步骤14：在近场连续区与远场离散区交界面插入界面过渡单元；

步骤15：参数赋值，包括：为FEM区设置岩石物理力学参数和应变软化模型；为DEM区配置平行黏结模型和初始地应力；使用Karhunen-Loève展开生成随机场并插值；设置时间步长并进行DEM到FEM的映射；利用Python+Gmsh进行自动生成网格，并设置监测数据更新的接口。

3.根据权利要求1所述的隧洞围岩稳定性预测方法，其特征在于，数据采集装置包括：微震传感器、 声发射传感器、分布式光纤及多点位移计。

4.根据权利要求1所述的隧洞围岩稳定性预测方法，其特征在于，步骤3基于内、外两个时间尺度进行参数更新，外循环每更新一次岩体宏观参数场；内循环在内按频率更新掌子面前方破裂面黏聚力与摩擦角。

5.根据权利要求4所述的隧洞围岩稳定性预测方法，其特征在于，外循环使用集合卡尔曼滤波更新参数，内循环使用无迹卡尔曼滤波更新参数。

6.根据权利要求5所述的隧洞围岩稳定性预测方法，其特征在于，步骤3还包括计算全局位移残差及UKF协方差，若不满足全局位移残差<1mm且UKF协方差< 0.01，则缩短。

7.根据权利要求1所述的隧洞围岩稳定性预测方法，其特征在于，双通道超图神经ODE破裂模式预测模型实现步骤为：

超图构建：将有限元单元与离散颗粒全部视为节点，采用三类超边分别捕捉地层结构、微震事件簇与位移梯度簇；

双通道超图神经ODE建模：将每个岩体节点分裂成两条通道：完整岩体通道和已损伤岩体通道，根据超图构建的结果，采用超图聚合消息驱动ODE右侧，让完整-损伤两通道差异直接控制损伤速率，同时进行物理约束损失设计对ODE右侧实时校正，输出破裂体积及对应的失稳概率。

8.根据权利要求1所述的隧洞围岩稳定性预测方法，其特征在于，还包括：基于FDEM模型正演结果对双通道超图神经ODE破裂模式预测模型进行重训练。

9.根据权利要求8所述的隧洞围岩稳定性预测方法，其特征在于，重训练时，训练版本号自增并存档。

10.根据权利要求1所述的隧洞围岩稳定性预测方法，其特征在于，还包括步骤6：基于预测结果进行可视化显示及预警。