Method for predicting stability of surrounding rock of tunnel

Abstract

The present application relates to surrounding rock stability prediction, in order to realize the advanced prediction of surrounding rock fracture under deep buried condition, a tunnel surrounding rock stability prediction method is provided, by adopting the coupling of finite element and discrete element, the overall deformation and microscopic crack propagation of surrounding rock can be considered at the same time, the accuracy and calculation efficiency of surrounding rock stability analysis are improved; the double-channel hypergraph neural ODE fracture mode prediction model can give the fracture volume and instability probability 3-7 days in advance, and the physical consistency problem is considered at the same time, that is, the energy and momentum residual error converges in real time, which ensures the accuracy of the prediction result; real-time monitoring data is obtained through the finite element-discrete element coupling model, and stability prediction is carried out based on the double-channel hypergraph neural ODE fracture mode prediction model, without manual adjustment of parameters, the system automatically completes the data analysis-prediction-warning closed loop; the stability analysis of deep buried tunnel surrounding rock is upgraded from static and empirical analysis to real-time and dynamic advanced prediction.

Description

Technical Field

[0001] This invention relates to the prediction of surrounding rock stability, specifically a method for predicting the stability of surrounding rock in tunnels. Background Technology

[0002] Rock stability analysis has always been an important research topic in the design and construction of underground space engineering projects. With the increasing depth of transportation, hydropower, and mining projects, tunnels are generally buried at depths exceeding 1,000 meters, placing the surrounding rock in a multi-field coupled environment characterized by high ground stress (>20 MPa), high osmotic pressure, strong unloading, and complex tectonic stress fields. Rock stability analysis can clarify the safety of tunnel construction, thereby ensuring safe tunnel construction. Common methods for rock stability analysis include field tests, physical tests, theoretical analysis, and numerical simulation; however, due to the unique dynamic mechanical properties of deep rock masses, each method has its own limitations.

[0003] Sampling deep rock masses is difficult, and in-situ tests (such as geostress testing and hydraulic fracturing) are costly and time-consuming. Furthermore, the test results are highly variable due to excavation disturbances, making it difficult to reflect the true stress path. Physical tests are limited by size effects (scaled-down models cannot reproduce the rock mass joint network and stress gradient) and loading capacity (existing true triaxial equipment is difficult to simulate stress fields at kilometer-level burial depths), resulting in significant deviations between test results and actual engineering conditions.

[0004] In terms of theoretical analysis, deep rock masses exhibit obvious nonlinear, discontinuous and anisotropic characteristics. Existing analytical solutions (such as the Kirsch equation and strain softening model) are all based on the ideal elastoplastic assumption and cannot take into account the temporal evolution of crack dynamic propagation and stress redistribution.

[0005] In numerical simulation, software-based numerical simulation has become an important research method in engineering practice. While the traditional finite element method (FEM) can simulate the overall deformation and stress distribution of surrounding rock relatively well, it is difficult to capture the initiation, propagation, and penetration process of micro-cracks within the rock mass. The discrete element method (DEM), although it can reveal microscopic failure mechanisms such as rock bridge fracture and block slip, has computational efficiency that decreases exponentially with the number of particles, making it difficult to handle large-scale continuous medium problems. In addition, existing methods generally adopt the "prior parameters + offline calculation" mode, and do not effectively integrate monitoring data during construction (such as multi-point displacement gauges and microseismic signals). The model parameters are no longer updated during construction and cannot reflect the dynamic response of the surrounding rock, resulting in a 3-7 day time lag between the calculated analysis results and the actual surrounding rock deformation / failure events, which fails to provide effective prevention and control. Summary of the Invention

[0006] To enable advanced prediction of surrounding rock fracturing under deep burial conditions, this invention provides a method for predicting the stability of surrounding rock in tunnels.

[0007] The technical solution adopted by the present invention to solve the above problems is:

[0008] Methods for predicting the stability of surrounding rock in tunnels include:

[0009] Step 1: Construct the FDEM model of the tunnel. The FDEM model is a finite element-discrete element coupled model.

[0010] Step 2: Collect monitoring data at the tunnel site using the data acquisition device;

[0011] Step 3: Update the parameters of the finite element-discrete element coupled model based on the monitoring data and obtain the forward modeling results;

[0012] Step 4: Establish a dual-channel hypergraph neural ODE rupture pattern prediction model based on monitoring data and forward modeling results;

[0013] Step 5: Predict the potential fracture volume and instability probability of the surrounding rock based on the dual-channel hypergraph neural ODE fracture mode prediction model.

[0014] Further, step 1 includes:

[0015] Step 11: Establish a deep-buried tunnel model based on the tunnel centerline; divide the tunnel model into concentric circle partitions according to the distance from the tunnel wall, and divide it into a near-field continuous zone, a far-field discrete zone, and a non-reflection absorption zone.

[0016] Step 12: Create a joint network based on the field data and assign values, then map the joint network to the far-field discrete region;

[0017] Step 13: Divide the near-field continuous region into finite element meshes, fill the far-field discrete region with discrete particles, and set shared nodes for the transition layer;

[0018] Step 14: Insert an interface transition unit at the interface between the near-field continuous region and the far-field discrete region;

[0019] Step 15: Parameter assignment, including: setting rock physical and mechanical parameters and strain softening model for FEM region; configuring parallel bond model and initial geostress for DEM region; generating random field and interpolating using Karhunen-Loève expansion; setting time step and mapping DEM to FEM; automatically generating mesh using Python+Gmsh and setting up interface for monitoring data updates.

[0020] Furthermore, the data acquisition device includes: a micro-vibration sensor, an acoustic emission sensor, a distributed optical fiber, and a multi-point displacement meter.

[0021] Furthermore, step 3 updates parameters based on both inner and outer time scales, with the outer loop updating parameters every... Update the macroscopic parameter field of the rock mass once; internal circulation in Internal press frequency Update the cohesion and friction angle of the fracture surface in front of the face.

[0022] Furthermore, the outer loop uses an ensemble Kalman filter to update the parameters, while the inner loop uses an unscented Kalman filter to update the parameters.

[0023] Furthermore, step 3 also includes calculating the global displacement residual and the UKF covariance. If the global displacement residual is not < 1 mm and the UKF covariance is not < 0.01, then the time limit is shortened. .

[0024] Furthermore, the implementation steps of the dual-channel hypergraph neural ODE rupture mode prediction model are as follows:

[0025] Hypergraph construction: Finite element elements and discrete particles are all treated as nodes, and three types of hyperedges are used to capture the stratigraphic structure, microseismic event clusters and displacement gradient clusters respectively;

[0026] Dual-channel hypergraph neural ODE modeling: Each rock mass node is split into two channels: intact rock mass channel and damaged rock mass channel. Based on the hypergraph construction results, the right side of the ODE is driven by hypergraph aggregation messages, allowing the difference between the intact and damaged channels to directly control the damage rate. At the same time, physical constraint loss design is performed to correct the right side of the ODE in real time, and the fracture volume and the corresponding instability probability are output.

[0027] Furthermore, it also includes: retraining the dual-channel hypergraph neural ODE rupture mode prediction model based on the forward modeling results of the FDEM model.

[0028] Furthermore, during retraining, the training version number is incremented and archived.

[0029] Furthermore, it also includes step 6: visualizing and issuing early warnings based on the prediction results.

[0030] The advantages of this invention compared to existing technologies are as follows: By employing the coupling of finite element and discrete element methods, it can simultaneously consider the overall deformation and micro-fracture propagation of the surrounding rock, improving the accuracy of surrounding rock stability analysis and avoiding the low efficiency problem of the single discrete element method in large-scale calculations, thus improving computational efficiency; the dual-channel hypergraph neural ODE fracture mode prediction model can provide the fracture volume and instability probability 3-7 days in advance, while also considering the physical consistency problem, i.e., the energy and momentum residuals converge in real time, ensuring the accuracy of the prediction results; and the damage rate is pre-written into the ODE through the integrity-damage dual-channel approach. On the right side of DE, a hypergraph is used to model the potential fracture surface in advance. With one integration, the trajectory 3-7 days later can be seen, achieving advanced prediction. Real-time monitoring data is obtained through a finite element-discrete element coupled model, and stability prediction is performed based on a dual-channel hypergraph neural ODE fracture mode prediction model. No manual parameter adjustment is required, and the system automatically completes the data analysis-prediction-early warning closed loop. This realizes the upgrading of the stability analysis of the surrounding rock of deep buried tunnels from static and empirical analysis to a real-time, dynamic, and advanced digital twin system, which is applicable to the stability prediction of the surrounding rock of deep buried tunnels with different geological conditions and tunnel sizes. Attached Figure Description

[0031] Figure 1 This is a flowchart of a method for predicting the stability of surrounding rock in tunnels. Detailed Implementation

[0032] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0033] like Figure 1 As shown, the method for predicting the stability of surrounding rock in tunnels includes:

[0034] Step 1: Construct the FDEM model of the tunnel. The FDEM model is a coupled finite element-discrete element model. This includes the following steps:

[0035] Step 11: Establish a deep-buried tunnel model based on the tunnel centerline; divide the tunnel model into zones using a concentric circle partitioning method, and according to the distance from the tunnel wall, divide it into a near-field continuous zone, a far-field discrete zone, and a non-reflection absorption zone.

[0036] Using the tunnel centerline as a reference, a deep-buried tunnel model (diameter D, length L≥30D) was established in Abaqus, and divided into concentric circles: FEM zone (near-field continuous zone 0~5D from the tunnel wall); DEM zone (far-field potential fracture zone 5D~20D from the tunnel wall); ABS zone (artificial damping boundary zone, outermost layer 20D~25D); among which, the 0-5D region has stress concentration in the tunnel wall but no through cracks have yet appeared, and the deformation is mainly in the continuous medium, so pure finite element method (FEM) is used; 5D-2 In the 0D region, micro-cracks are most likely to appear in the surrounding rock under high ground stress. If joint propagation occurs, it will lead to block slippage. Discrete element method (DEM) is required to simulate discontinuous phenomena such as cracking and slippage. In the 20D-50D region, external stress disturbance can be ignored, but the boundary must absorb wave fluctuations to prevent reflection. Therefore, the ABS region is set. The ABS region is an unaffected finite element region that is different from the FEM region. It is only used to absorb reflected waves and does not participate in the calculation of the mechanical response of the surrounding rock. Its material parameters are set as viscoelastic medium with the elastic modulus decreasing layer by layer to avoid stress wave reflection interference.

[0037] Step 12: Create a joint network based on the field data and assign values, then map the joint network to the far-field discrete region.

[0038] Geological strength index (GSI), rock quality index (RQD), and joint occurrence statistics from the field data were imported into Python to generate a three-dimensional Voronoi joint network. Then, macroscopic parameters were assigned to the entire rock mass, and the Monte Carlo algorithm was used to assign random cohesion c and internal friction angle to each joint. The normal stiffness Kn and tangential stiffness Ks are calculated, and the joint network is finally mapped onto the DEM region as the initial particle bonding location for subsequent discrete element analysis. This mapping process uses a bidirectional hash index structure to maintain the one-to-one correspondence between Voronoi joint surface IDs and DEM particle bonds, ensuring that joint numbers and spatial coordinates can be retrieved in subsequent fracture surface evolution.

[0039] Step 13: Divide the near-field continuous region into finite element meshes, fill the far-field discrete region with discrete particles, and set shared nodes for the transition layer.

[0040] 1) Finite element mesh generation is performed in the FEM zone. The finite element mesh adopts 8-node hexahedral reduced integral (C3D8R) elements, with 0.1D near the cavity wall and increasing by 1.3 times outward. One layer of common node elements is reserved at the interface between the FEM zone and the DEM zone.

[0041] 2) In the DEM region, spherical / polyhedral particles are generated using the radius expansion method to ensure that the local porosity is <3% and the particle radius r satisfies r≤0.05D, which can distinguish the joint thickness. At the interface between the DEM region and the FEM region, common node particles are forcibly generated within a thickness of 1.5r, that is, the center of these particles falls on the finite element node coordinates.

[0042] Step 14: Insert an interface transition unit at the interface between the near-field continuous region and the far-field discrete region.

[0043] Using a shared-node mapping approach, the KD-Tree algorithm is employed to perform nearest-neighbor matching between finite element nodes and discrete particle centers, generating 0.3 m thick interface elements at the interface. Goodman joint elements combined with a damage evolution constitutive relation are used, with initial normal stiffness set. ,in The ratio of the elastic modulus to the characteristic length of the element in the FEM region reflects the axial stiffness per unit thickness of the FEM region. The ratio of the equivalent elastic modulus to the average particle size of the particle aggregate in the DEM region reflects the axial stiffness per unit thickness of the DEM region. A smaller value is chosen to prevent excessive stiffness in the transition layer from causing stress concentration or numerical reflection. The tangential stiffness Ks = 0.4Kn. When the interface displacement > 0.5 mm, it automatically degenerates into zero-thickness contact, allowing slippage or cracking.

[0044] Step 15: Parameter assignment, including: setting rock physical and mechanical parameters and strain softening model for FEM region; configuring parallel bond model and initial geostress for DEM region; generating random field and interpolating using Karhunen-Loève expansion; setting time step and mapping DEM to FEM; automatically generating mesh using Python+Gmsh and setting up interface for monitoring data updates.

[0045] 1) The physical and mechanical parameters of the rock in the FEM zone, such as the elastic modulus E and Poisson's ratio ν, are assigned using a strain softening model. The softening modulus is calibrated by the triaxial test curve in the field.

[0046] 2) Parallel bonding model (PBM) is used between particles within the DEM region, and the normal bond strength is... Tangential bond strength The initial geostress field is mapped from joint network statistics using Voronoi element volume weighting. Calculation, where This is the lateral pressure coefficient, derived from hydraulic fracturing tests. For vertical stress, by =γh is calculated (γ is the average unit weight of the rock mass, h is the burial depth);

[0047] 3) Use Karhunen-Loève expansion to generate random fields (correlation length = 2D); interpolate the random fields to the finite element integration points and discrete bonds to generate initial spatial random field samples of the macroscopic physical and mechanical parameters of the rock mass. This includes elastic modulus E, Poisson's ratio v, and lateral pressure coefficient. wait;

[0048] 4) Set the time step and finite element step. Discrete element step size Furthermore, a coarsening mapping is performed every 200 discrete element steps to transmit the resultant force and resultant torque of the particles to the finite element nodes.

[0049] 5) Use Python + Gmsh to automatically generate the grid and set up an interface for monitoring data updates to ensure that subsequent monitoring data can be directly interpolated to nodes or particles.

[0050] Step 2: Collect monitoring data at the tunnel site using the data acquisition device.

[0051] During the excavation and operation phases of deeply buried tunnels, sensors, multi-point displacement gauges, stress gauges, and other devices are deployed along the tunnel axis and typical cross-sections to collect and store data, establishing a multi-source monitoring data stream interface encompassing space, time, and frequency. This multi-source monitoring data stream includes:

[0052] (1) Space: 32 channels of microseismic (MS), 16 channels of acoustic emission (AE), 2km of distributed optical fiber (DAS) and 12 points of multi-point displacement meter (MPBX) form a double-ring monitoring network of "circumferential + longitudinal".

[0053] (2) Time: Real-time sampling from 1 Hz to 1 kHz, micro-vibration / AE 1 kHz, DAS 10 kHz, MPBX 1 Hz;

[0054] (3) Frequency: Extract 20-dimensional feature vectors such as microseismic b-value, displacement velocity spectral index β_disp, and dispersion curve. .

[0055] Step 3: Update the parameters of the finite element-discrete element coupled model based on the monitoring data and obtain the forward modeling results.

[0056] This embodiment implements FDEM model updates driven by monitoring data through a dual-loop assimilation algorithm across both inner and outer time scales. The outer loop employs an ensemble Kalman filter (EnKF-AE) based on adaptive covariance dilation. The macroscopic parameters of the rock mass are updated every 6 hours (e.g., every 6 hours); the internal circulation is updated every 6 hours. Internally, the cohesion c and friction angle of the fracture surface are analyzed using unscented Kalman filtering (UKF). Perform high-frequency correction (e.g., 5-minute intervals). The specific steps are as follows:

[0057] (1) First, the initial physical and mechanical parameters of the rock mass are specified based on the borehole information and the statistical strength of the joint network. and fracture surface strength parameters For the i-th member Its initial parameter set is denoted as: ,in This is a vector of macroscopic rock mass parameters. Let N() be the vector of fracture surface strength parameters, which follows a normal distribution, and the prior mean of the macroscopic rock mass parameter vector is... Prior mean of fracture surface strength parameter vector The covariance Σθ of the macroscopic rock mass parameter vector and the prior mean Σφ of the fracture surface strength parameter vector are given by field tests and prior statistics. It is the largest member of the set.

[0058] (2) Use ensemble Kalman filtering (EnKF-AE) based on adaptive covariance dilation in the outer loop. The specific steps to update the macroscopic parameter field of the rock mass are as follows:

[0059] 1) For each set member Running FDEM forward modeling for 6 hours yields prior estimates of the macroscopic parameter field of the rock mass. , ,in This is the FDEM forward modeling operator.

[0060] 2) Read 6 hours of cumulative observation data from the "space-time-frequency" interface, and then extract virtual observation vectors with the same dimensions and location as the field monitoring based on prior estimates. For example, cave wall convergence and micro-seismic energy. ,in This is the functional form of the observation operator, used to extract displacement, microseismic energy, etc.

[0061] 3) Set the covariance inflation coefficient Calculate the observation error covariance matrix among the members of the set. In this embodiment To avoid filter divergence. ,in The mean vector of the set of virtual observations. , This is the matrix transpose symbol.

[0062] 4) Calculate the Kalman gain This is used to quantify the weights of model error and observation error; the larger the matrix, the more trust is placed in the observations. Where H is the observation operator matrix and R is the instrument noise covariance, given by the manufacturer (displacement gauge ±0.1 mm, micro-vibration energy ±5%).

[0063] 5) Observation residuals By flipping back to the rock mass parameter space using the gain, a model calibration is completed, obtaining the updated macroscopic rock mass parameter vector of the i-th set member at time k. , , The values ​​are the actual measured values ​​(displacement, microseismic energy, etc.) at time k.

[0064] (3) Internal circulation Internal use of unscented Kalman filtering by frequency The cohesion and friction angle of the fracture surface are corrected. The specific steps are as follows:

[0065] 1) Using a window of 30 m in front of the tunnel face as the reference, extract all fracture surface elements and set the state vector. State vector Defined as the equivalent cohesion c and internal friction angle of all activated fracture surfaces within a 30 m window in front of the tunnel face at time k. Its initial value is obtained by weighting the statistical intensity of the joint network by the volume of Voronoi elements. Where c is the cohesive force. It is the internal friction angle.

[0066] 2) Set up Sigma points with dimension n, using (2n+1) Sigma points to capture the mean and covariance of the state distribution. Set a center point for calculating the mean weight. Covariance weights The predicted mean weight is calculated with the other 2n Sigma points. Covariance weights .

[0067] ,

[0068] ,

[0069] ,

[0070] ,

[0071] ,

[0072] ,

[0073] in, The parameter is used to adjust the distance between the Sigma point and the mean. To adjust the scaling factor for converting parameters into actual spatial step size, This is the j-th positive deviation point (along the covariance principal axis + γ direction). This is the j-th negative deviation point (along the covariance principal axis -γ direction). Decompose the covariance matrix for Cholesky Column vectors.

[0074] 3) Perform a 5-minute fast FDEM sub-loop for each Sigma point, which can update the loop status. And extract virtual observations Calculate the prior mean of the sigma point at time k+1. and prior covariance .

[0075] ,

[0076] ,

[0077] ,

[0078] ,

[0079] in, For the j-th Sigma point generated at time k, For a 5-minute mini FDEM sub-model, Displacement, energy, and dominant frequency are extracted; Q represents the process noise covariance. Where diag() is the symbol for a diagonal matrix. Let c be the noise variance. Friction angle The noise variance.

[0080] 4) Calculate the Kalman gain The predicted value is transformed into the optimal estimate for the state vector. With covariance Perform optimal updates.

[0081] ,

[0082] ,

[0083] ,

[0084] in, These are the measured values ​​(displacement, microseismic energy, etc.) at time k+1. The model predicts the observed values ​​at the same time (obtained from the prior state through the observation operator). For state-observation cross-covariance, For mutual covariance, This is the matrix transpose symbol.

[0085] (4) Write the macroscopic parameters obtained from the outer circulation and the local intensity obtained from the inner circulation back into the finite element mesh and discrete element particles of the FDEM model, and ensure energy conservation and coordinate consistency.

[0086] Furthermore, the system compares the latest FEM / DEM forward modeling results with the measured values ​​from multiple displacement gauges and optical fibers in real time to obtain the global displacement residuals. Simultaneously, the covariance matrix P output by the UKF is read, and its trace Tr(P) is calculated to quantify the current parameter uncertainty; if the global displacement residual The parameters are maintained when the distance is < 1 mm and the UKF covariance Tr(P) < 0.01; otherwise, it is marked as a "strong disturbance" and the outer loop is shortened by ΔT = 3 h in the next iteration.

[0087] Existing technologies in engineering applications often employ traditional static information processing methods, such as obtaining rock physical and mechanical parameters through core drilling and triaxial compression tests and inputting them into the model. The numerical model is not updated during construction, and manual parameter adjustment is required for data correction, resulting in significant data uncertainty. Furthermore, parameter adjustment is limited to the elastic modulus E. In contrast, the dual-loop assimilation algorithm designed in this invention utilizes multi-source information fusion processing technology to drive the updating of the numerical simulation model using monitoring data. It employs displacement, microseismic, and other monitoring data for autonomous real-time updates of the finite element and discrete element models, eliminating the need for manual parameter adjustment. The outer loop (EnKF-AE) transforms macroscopic rock mass parameters into a random field form, performing autonomous correction every 6 hours. The inner loop (UKF) performs high-frequency correction of the cohesion c and internal friction angle φ of the fracture surface every 5 minutes, achieving a cumulative deviation of <1mm over 30 days, compared to >10mm for traditional static models during the same period.

[0088] Step 4: Establish a prediction model for the ODE rupture pattern of a dual-channel supergraph neural network based on monitoring data and forward modeling results.

[0089] The implementation steps of the dual-channel hypergraph neural ODE rupture mode prediction model are as follows:

[0090] (1) Hypergraph construction. Setting time windows. , At the current external circulation time (in h), all FEM elements and DEM particles are treated as nodes, and three types of heterogeneous hyperedges are used to capture the stratigraphic structure, microseismic event clusters, and displacement gradient clusters.

[0091] 1) All nodes within the same stratum or joint group are pulled together to form a higher-order edge, naturally accommodating interlayer slip. Based on the joint occurrence and Geological Intensity Index (GSI) zoning, geological bedding equations are established: ;in This is the number of the m-th geological layer (rock stratum or joint surface), used to distinguish different strata; is the unit normal vector of this plane, pointing upwards from the plane; x is the coordinate of any point in space, used to calculate its relative position to the plane. M represents the coordinates of a known reference point on the stratum, used to pinpoint the specific location of the stratum in three-dimensional space; M is the stratum number, numbered sequentially from 1 to M, covering all major rock strata and joint groups within the model area.

[0092] For any node i, compute the surface distance .like Then i will be added to the set. ; δ represents the distance tolerance from the node to the rock layer / joint surface, reflecting the accuracy of rock layer thickness identification. In this embodiment, it is taken as 0.5m.

[0093] Each structural hyperedge Right now ,in, Let m be the set of all nodes captured by the m-th rock layer (or joint surface). As the structural hyperedge weight, the average microseismic energy density within the edge is accumulated on top of 1 to highlight the interaction intensity of high-energy rock layers. For node i in The accumulated microseismic energy density within the body.

[0094] 2) Connect nodes within the same rupture event or stress zone into higher-order edges to capture damage localization. Extract from monitoring flow. All microseismic events ,in Let be the spatial coordinates of the s-th microseismic event, representing the location of the rupture source. The time of event occurrence is used for time window filtering and sorting. 's' is the event sequence number, numbered sequentially from 1 to S, and S is the total number of valid events in the current window. The energy released for the event determines its contribution to the hyperedge weights, and a Delaunay simplex is constructed in 3D space, i.e.: For any simplex ,like Then the simplex is preserved, and its vertices are expanded to the nearest set of FEM / DEM nodes. The energy threshold for microseismic events is taken in this embodiment. Only events whose average energy within the simplex exceeds this value are retained, and this is used to filter out noisy events.

[0095] Each event exceeds the edge Right now ,in It is a simplex microseismic pattern Center, event radius The set of all FEM / DEM nodes within the system. The event super-edge weight is accumulated by adding logarithmic energy to the base of 1, so that the high-energy fracture clusters can obtain a greater propagation weight.

[0096] 3) Cluster nodes with similar displacement gradients into a hyperedge to capture continuous deformation zones. Node-level displacement vectors are obtained through interpolation of multi-point displacement gauges and fiber optic data. , , where x i Let i be the spatial coordinates of node i.

[0097] Calculate the second-order displacement gradient tensor (Spherical neighborhood Ωi least squares fitting is used), and invariant features are extracted. Then, for all node features Execute K-Means (K=64) to obtain cluster labels. ,use Indicates the cluster number, The value range is 1, ..., K.

[0098] ,in The spatial position of node i The matrix for finding the partial derivatives of the displacement at a given point.

[0099] ,

[0100] in, The gradient magnitude (Frobenius norm) is the overall size of the gradient; a larger value indicates a more severe deformation. It is a gradient determinant that reflects the trend of volume change (its value > 0 indicates expansion, its value ≈ 0 indicates shear, and its value < 0 indicates compression). It is the gradient partial norm, representing the magnitude of pure shear deformation after eliminating volume changes.

[0101] Each gradient hyperedge Right now .

[0102] 4) Hypergraph assembly. Merge all hyperedges to obtain the hyperedge set. Correspondence matrix Hyperedge weight diagonal matrix .

[0103] ,

[0104] ,

[0105] .

[0106] Here, ∪ represents the union, where the three elements are distinct and cover each other, forming the complete hypergraph edge set at the current moment. E×V is the size (E superedges, V nodes). The function is for determining whether a node i belongs to hyperedge e, and its value is either 0 or 1: =1; otherwise =0, diag() is the symbol for a diagonal matrix, w() stores each hyperedge sequentially on the diagonal. … weights, E is the dimension (E = total number of hyperedges).

[0107] (2) Dual-channel hypergraph neural ODE modeling. This is achieved through initial observation features of nodes. Set two initial states, namely channel A (intact rock mass). Channel B (damaged rock mass) .

[0108] ,in, The values ​​of the macroscopic rock mass parameters random field at node i (including elastic modulus E, Poisson's ratio) Lateral pressure coefficient (equal subvectors) The value of the random field of fracture surface strength at node i (sub-vectors such as joint cohesion c and friction angle φ). The current stress tensor components at node i (usually 3 normal stresses + 3 shear stresses, or invariant form). The current strain tensor components of node i (and) Correspondingly, the elastic energy density can be calculated. This is a scalar representation of the microseismic energy density interpolated to node i, reflecting the local rupture activity. is the Frobenius norm of the displacement gradient tensor, used to quantify the severity of deformation in the neighborhood of a node.

[0109] ,

[0110] in , It consists of two layers of MLP (including LayerNorm & ReLU).

[0111] Collect aggregate message vectors for any hyperedge e. :

[0112] ,

[0113] in Let e ​​be the aggregated message vector of the hyperedge, which is composed of the hidden states of all nodes within the edge. The result, after averaging and gating, is used for subsequent transmission to nodes. For small-scale MLPs, the input is the relative coordinate offset of the nodes. With the average microseismic energy density within the boundary Output gate vector To implement gated weighting for the Hadamard product, .

[0114] Node i receives messages from all superedges containing i. : Where e is the hyperedge number containing node i, which determines the range of message sources. The weight of the hyperedge e reflects the microseismic energy or the degree of deformation within that edge; the higher the energy, the greater the weight.

[0115] The aggregated message is split back into channels A and B, and then channels A and B are updated independently respectively. .

[0116] Define the latent space difference norm Node-level damage variables Satisfying the initial value It follows an evolutionary pattern. When the difference between the two channels exceeds the threshold β, damage accumulates rapidly and is integrated in a timely manner.

[0117] ,

[0118] ,

[0119] in, σ is the damage rate gain coefficient; σ() is the ReLU activation function, i.e., σ(x) = max(0,x), ensuring that damage accumulates only when the difference between the two channels exceeds the threshold β. Let i represent the physical volume (m³). For FEM elements, the volume is the element volume, and for DEM particles, the volume is the particle volume. The cumulative fracture volume of the entire model at time t is used to measure the degree of damage to the surrounding rock and to calculate the instability probability P(f).

[0120] (3) Physical constraint loss (PIL) design. Physical laws such as energy conservation, momentum balance, and constitutive relations are embedded in the loss function, so that the prediction results of the constraint model are consistent with the physical reality, thereby improving the credibility and generalization ability.

[0121] To avoid additional losses in weighted parameter tuning, the energy residual is... With momentum residual Write it as the right-hand correction.

[0122] ,

[0123] in, The total power integral (input energy) of the external force on the rock mass at time t is calculated from the boundary displacement and the external force integral. The elastic energy increment within the rock mass is calculated using stress-strain analysis. The kinetic energy increment of the rock mass is estimated from the nodal accelerations. The fracture energy (unit: J / m²) is estimated from the predicted fracture region volume. The material fracture energy (constant) is estimated by weighting the nodal-level damage variable Di(t) with the element surface area.

[0124] ,

[0125] Right-side correction (applies only to channel A, ensuring the complete side is closer to physics).

[0126] ,

[0127] in, The stress divergence at node i is obtained by the neighborhood stress difference approximation and reflects the internal force distribution. For volumetric density, gravity is usually taken as the unit. The direction is vertically downward; The lumped mass of node i is obtained by converting the element or particle volume and density. Let be the acceleration vector of node i, which is calculated in real time from the second-order difference of displacement or the derivative of velocity; Provide the internal force gradient for the stress tensor difference between neighbor j and center node i; The distance between the two nodes is used for differential scale normalization; The unit direction vector from i to j determines the difference sign and projection direction; This is the intensity coefficient for energy residual correction (dimensionless), with an initial value of 0.01 that is adjustable. Energy residual The sign function is +1, which indicates that the total energy of the system is excessive (work done by external forces > internal dissipation + elastic energy), -1 indicates that the total energy of the system is deficient, and 0 indicates that the system is in balance; Let be the hidden state gradient of the A-channel of the energy residual at node i.

[0128] Step 5: Predict the potential fracture volume and instability probability of the surrounding rock based on the dual-channel hypergraph neural ODE fracture mode prediction model.

[0129] Furthermore, it also includes step 6: visualizing and issuing early warnings based on the prediction results.

[0130] Based on the prediction results, a rupture volume-instability probability curve is plotted. The corresponding instability probability is found by comparing the rupture volume V(f) with the curve, yielding the instability probability P(f) at each time step. A red alert is triggered if the instability probability P(f) ≥ 0.7, and a reinforcement plan is automatically recommended; a yellow alert is triggered if the instability probability P(f) ≤ 0.5 < 0.7, and monitoring frequency is increased; a green alert is granted if the instability probability P(f) < 0.5. The system also provides a 95% confidence band for decision-making reference.

[0131] To improve the accuracy of model predictions, it is also necessary to update the model; this includes the following methods:

[0132] (1) The system automatically retrieves the latest data from the previous 24 hours every day at midnight and performs incremental fine-tuning of the dual-channel ODE mapping network over 10 epochs, using cosine annealing for the learning rate. Only the prediction head and gating parameters are updated, while most convolutional and attention layers are frozen to prevent catastrophic forgetting. The version number is auto-incremented and archived, supporting one-click rollback to any historical version. The version number is the number corresponding to each prediction model. Through the version number, it can be ensured that under the same model skeleton, each training result has a unique, incremental, and non-repeatable label, avoiding the overwriting of old weights and facilitating rollback or comparison.

[0133] (2) If the global displacement residuals for three consecutive days If the UKF covariance Tr(P) < 0.01, then maintain the current model; otherwise, trigger double cyclic assimilation again, recalibrate the parameters, and then train again.

[0134] (3) After the outer loop is completed, a consistency check is performed. That is, the FDEM forward modeling and the dual-channel hypergraph neural ODE rupture mode prediction model are run under the same initial conditions. If the relative error of the rupture volume of the two is >10%, the hypergraph retraining is triggered (incremental fine-tuning for 20 epochs) to ensure that the deviation between the prediction model and the physical model is <5%.

[0135] By updating model parameters driven by monitoring data, the model parameters are made consistent with the actual rock mass, improving the accuracy of the forward modeling results of the finite element-discrete element coupled model and the prediction results of the dual-channel hypergraph neural ODE fracture mode prediction model. The prediction model is then corrected using the forward modeling results, further enhancing the prediction accuracy. Using the method of this invention, fracture prediction can be made 3-7 days earlier, with a location error of <1.5 m and a time error of <12 h, achieving truly real-time, dynamic, and advanced surrounding rock stability analysis.

Claims

1. A method for predicting the stability of surrounding rock of a tunnel, characterized by, The application comprises the following steps: Step 1: constructing a tunnel FDEM model, wherein the FDEM model is a finite element-discrete element coupling model; Step 2: collecting monitoring data of the tunnel site based on a data acquisition device; Step 3: updating parameters of the finite element-discrete element coupling model based on the monitoring data and obtaining a forward result; Step 4: establishing a double-channel hypergraph neural ODE fracture mode prediction model based on the monitoring data and the forward result; Step 5: predicting a potential fracture volume and a stability probability of surrounding rock based on the double-channel hypergraph neural ODE fracture mode prediction model. Specifically, step 1 comprises the following steps: Step 11: establishing a deep-buried tunnel model based on a tunnel center line; adopting a concentric circle partitioning mode to partition the tunnel model; and sequentially dividing the tunnel model into a near-field continuous zone, a far-field discrete zone and a non-reflection absorption zone according to distances from the tunnel wall; Step 12: creating a joint network and assigning values according to field data, and then mapping the joint network to the far-field discrete zone; Step 13: dividing the near-field continuous zone into a finite element grid, filling the far-field discrete zone with discrete particles, and setting a transition layer shared node; Step 14: inserting an interface transition unit at an interface between the near-field continuous zone and the far-field discrete zone; Step 15: parameter assignment, including: setting rock physical and mechanical parameters and a strain softening model for the FEM zone; configuring a parallel bonding model and an initial ground stress for the DEM zone; generating a random field using Karhunen-Loève expansion and interpolation; setting a time step and mapping DEM to FEM; automatically generating a grid using Python+Gmsh, and setting an interface for monitoring data updating; The implementation steps of the double-channel hypergraph neural ODE fracture mode prediction model are as follows: Hypergraph construction: regarding all finite element units and discrete particles as nodes, and adopting three types of hyperedges to capture stratum structure, microseismic event clusters and displacement gradient clusters; Double-channel hypergraph neural ODE modeling: splitting each rock mass node into two channels: an intact rock mass channel and a damaged rock mass channel; according to the results of the hypergraph construction, using hypergraph aggregation message to drive the right side of the ODE, letting the intact-damaged two channels directly control the damage rate, and simultaneously performing physical constraint loss design to real-time correct the right side of the ODE, and outputting a fracture volume and a corresponding stability probability.

2. The tunnel surrounding rock stability prediction method according to claim 1, characterized by, The data acquisition device comprises: a microseismic sensor, an acoustic emission sensor, a distributed optical fiber and a multi-point displacement meter.

3. The tunnel surrounding rock stability prediction method according to claim 1, characterized by, Step 3 updates the parameters based on the inner and outer time scales, the outer loop is updated every macroscopic parameter field of rock mass is updated once; the inner loop is updated in the inner loop is updated according to the frequency cohesion and friction angle of the fracture surface in front of the working face.

4. The tunnel surrounding rock stability prediction method according to claim 3, characterized by, The outer loop uses a set Kalman filter to update parameters, and the inner loop uses an unscented Kalman filter to update parameters.

5. The tunnel surrounding rock stability prediction method according to claim 4, characterized by, Step 3 further comprises calculating global displacement residual and UKF covariance, and if global displacement residual < 1mm and UKF covariance < 0.01 are not satisfied, then shortening .

6. The tunnel surrounding rock stability prediction method according to claim 1, characterized by, The application further comprises the following steps: Based on the forward result of the FDEM model, the double-channel hypergraph neural ODE fracture mode prediction model is retrained.

7. The tunnel surrounding rock stability prediction method according to claim 6, characterized by, During retraining, the training version number is incremented and archived.

8. The tunnel surrounding rock stability prediction method according to claim 1, characterized by, The application further comprises step 6: visualizing and warning based on the prediction result.

Images