Three-dimensional numerical simulation method for seepage field evolution of soft rock tunnel based on digital twinning

Abstract

This invention discloses a three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins, belonging to the field of numerical simulation technology for tunnel engineering. The method includes: constructing a follow-up in-situ sensing network integrated into the initial support and advanced pipe roof; combining advanced geological prediction to obtain the morphology of the preceding fault zone, generating an initial digital twin containing a three-dimensional discrete fracture network; establishing a dynamic fracture aperture evolution model controlled by four mechanisms: tunneling unloading, support reaction force, water absorption expansion, and hydraulic scouring; using graph neural networks to topologically reshape the fracture network map based on the spatial advancement coordinates of the tunnel face, dynamically updating the edge connection states; and utilizing real-time monitoring data transmitted from the steel arch fiber optic cable and deep-hole piezometer, asynchronously assimilating the twin parameters through ensemble Kalman filtering, ultimately predicting the extreme values ​​of the seepage field and the breakthrough path for the next tunneling cycle.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for tunnel engineering, specifically a three-dimensional numerical simulation method for the evolution of seepage field in soft rock tunnels based on digital twins. Background Technology

[0002] When soft rock tunnels are excavated in high-water-head sections, the surrounding rock contains numerous irregular fissures, and the seepage state is strongly disturbed by excavation. Existing numerical simulation methods for seepage fields typically treat the fissure network as a static geometric object, setting its opening as a constant value or adjusting it only according to a simple stress-permeability relationship during the simulation process. This fails to simultaneously consider the time-dependent mechanisms unique to soft rock, such as water swelling and fine particle piping loss, and the coupling effects of excavation unloading and support application. Furthermore, once the conventional fissure network model is discretized, its edge-node topology remains fixed. Even if local fissures are completely closed or scoured into through channels, the breakage and reconstruction of seepage paths cannot be reflected at the topological level, leading to distortion of the conduction network on which subsequent seepage calculations are based. In terms of utilizing monitoring data, traditional methods often involve offline inversion of field measured information, which cannot achieve rolling correction and advanced prediction of the seepage field in unexcavated areas as the tunnel face moves forward and new monitoring data continuously emerges. The aforementioned deficiencies mean that reliable technical tools are lacking for the prediction of water and mud inrush disasters during soft rock tunnel construction. There is an urgent need for a numerical simulation method that can realistically depict the dynamic competition of multiple mechanisms of fracture aperture and reconstruct the seepage topology in real time as the aperture changes, in order to solve the problem of advanced identification of seepage disaster paths. Summary of the Invention

[0003] This paper presents a three-dimensional numerical simulation method for the evolution of seepage field in soft rock tunnels based on digital twins. This method addresses the problems that existing methods cannot simultaneously reflect the dynamic evolution of the opening of soft rock fissures under multiple influences such as unloading, support, expansion and scouring, and the inability to adjust the seepage network topology in real time. It enables the advanced prediction of the extreme values ​​of the seepage field ahead and the breakthrough path.

[0004] The objective of this invention can be achieved through the following technical solutions:

[0005] This invention discloses a three-dimensional numerical simulation method for the evolution of seepage field in soft rock tunnels based on digital twins. It constructs an initial digital twin by integrating a servo-driven in-situ sensing network with a three-dimensional discrete fracture network. A dynamic fracture aperture evolution model is established, controlled by four mechanisms: tunneling unloading, support reaction force, water absorption and expansion, and hydraulic scouring. Furthermore, a graph neural network driven by the spatial advancement coordinates at the tunnel face is used to topologically reshape the fracture network map, dynamically updating the connection state of the fracture edges. Finally, monitoring data transmitted in real-time from the steel arch fiber optic cable and deep-hole piezometer is used to asynchronously assimilate the twin parameters through ensemble Kalman filtering. This allows for the advanced prediction of the seepage field extremes and breakthrough paths for the next tunneling cycle, thereby achieving intelligent prediction and precise control of the risk of water and mud inrush in soft rock tunnels.

[0006] As a technical solution of this invention, the process of constructing a follow-up in-situ sensing network integrated in the initial support and advanced pipe roof to generate an initial digital twin specifically involves: embedding fiber optic strain sensors along the longitudinal main reinforcement direction into the flanges and webs of each steel arch frame; simultaneously binding the fiber optic strain sensors to the inner wall of the steel pipes of the advanced pipe roof, forming a distributed strain sensing array that moves synchronously forward with the tunnel face; and embedding miniature multi-parameter sensors composed of permeability sensors and moisture content sensors in the sidewalls of the completed support section behind the tunnel face and in the advanced exploratory boreholes in front of the tunnel face. Sensors are evenly spaced along the depth direction within the exploratory borehole. Using seismic wave reflection and ground-penetrating radar methods from tunnel advanced geological prediction, the geometry and spatial distribution of the fault fracture zone and primary fracture group ahead of the tunnel face are jointly inverted. The coordinates of all measuring points of the fiber optic strain sensors and miniature multi-parameter sensors, along with the inverted geometric parameters of the fracture group, are uniformly mapped to the same three-dimensional spatial coordinate system, generating a three-dimensional discrete fracture network containing initial fracture aperture, strike, dip angle, and roughness attributes. This fracture network and the surrounding rock medium together constitute an initial digital twin. Preferably, the fiber optic strain sensors are demodulated using wavelength scanning optical frequency domain reflection technology. The sensors implanted on each steel arch frame form an independent strain monitoring section, thereby accurately capturing the stress state of each part of the entire steel arch frame.

[0007] As a further technical solution of the present invention, when establishing a dynamic fracture aperture evolution model in a digital twin, the hydraulic aperture of each discrete fracture is regarded as a dynamic variable that changes with time and spatial coordinates. This variable is composed of the initial hydraulic aperture and the aperture increment generated by four driving forces. The first driving force is the excavation unloading opening amount. Based on the three-dimensional spatial distance between the current working face and the fracture location, the normal tensile displacement of the fracture surface caused by the redistribution of surrounding rock stress is calculated. This displacement gradually accumulates as the working face advances. The second driving force is the support reaction closure amount. Based on the dynamic support force applied to the surrounding rock surface by the steel arch and shotcrete, the compressive displacement generated by the support force in the normal direction of the fracture surface is calculated. Among them, the radial support force of the steel arch is obtained by combining the real-time strain data of the fiber optic strain sensor with the elastic modulus and the moment of inertia of the section. The bonding force and surface constraint force of the shotcrete are obtained based on its mix ratio and the strength growth curve at age. The two are combined into a dynamic support resultant force, and then the compressive displacement is calculated by combining the spatial distance between the fracture and the steel arch using a semi-infinite plane force model of elasticity. The third driving force is the water absorption and expansion closure. Based on the local water content data transmitted in real time by a micro multi-parameter sensor, the compressive displacement generated on the two sides of the fracture by the volume expansion of soft rock upon contact with water is calculated. This displacement is in the opposite direction to the tensile displacement. The fourth driving force is the hydraulic mudification expansion capacity. When the hydraulic gradient of seepage within the fracture exceeds the critical initiation value of the soft rock particles, the piping erosion module is activated to calculate the mass of soft rock particles stripped from the fracture wall and to convert it into the increased fracture aperture value due to particle loss. The particle stripping rate is proportional to the cube of the local seepage velocity, and the expansion capacity calculation automatically stops when the stripped mass reaches the upper limit of the total erodible volume of the fracture wall. Through the synergistic effect of these four mechanisms, the model can realistically reflect the dynamic evolution of fracture aperture in soft rock tunnel surrounding rock under the coupled action of excavation, support, and seepage.

[0008] As another technical solution of the present invention, when using a graph neural network to reshape the fracture network map, each discrete fracture in the three-dimensional discrete fracture network is first mapped to an edge in the graph data structure, and the fracture intersection point or the contact surface between the fracture and the surrounding rock matrix is ​​mapped to a node, forming an initial fracture map topology. Each node is assigned a feature vector containing spatial coordinates and pressure head, and each edge is assigned a feature vector containing dynamic hydraulic aperture, current permeability coefficient, and water content on both sides. The three-dimensional spatial coordinates of the tunnel face at the end of the current tunneling cycle are used as external input and fed into the multilayer graph neural network model to iteratively aggregate and update the feature vectors of the edges. When the dynamic hydraulic aperture of a certain edge converges to zero or is lower than a preset closure threshold, the graph neural network model marks the state of the edge as open, indicating that the fracture is completely closed and water flow cannot pass through. When the hydraulic aperture of a certain edge exceeds a preset widening threshold, the edge is marked as thickened, indicating that it has been scourned into a dominant seepage channel. The graph neural network model outputs the updated fracture map topology, and the state changes of open or thickened edges trigger the redistribution calculation of pressure head on adjacent nodes. Preferably, the graph neural network model uses a graph attention layer to weighted aggregate the features of adjacent nodes and edges. The attention coefficient is dynamically determined by the relative distance between the current face advance coordinates and the spatial orientation of the fracture, so that the topology reshaping of the model can focus on the changes in fracture connectivity in the face excavation disturbance area.

[0009] Simultaneously, this invention utilizes real-time monitoring data transmitted via steel arch fiber optic cables and deep-hole piezometers to asynchronously assimilate twin parameters using ensemble Kalman filtering. Within each assimilation time window, distributed fiber optic sensors on the steel arch collect surrounding rock deformation pressure values ​​at each measuring point, while piezometers collect pore water pressure values ​​at each deep-hole measuring point, forming an observation vector. The unloading modulus, soft rock expansion coefficient, and piping erosion rate from the dynamic fracture aperture evolution model are used as state variables to be assimilated, and an initial set following a priori distribution is generated for each state variable. Before the observation data is input into the ensemble Kalman filter, a time window alignment operation between the observation data and model prediction is performed. This involves extracting observation values ​​at the same physical time point from the original data stream to form the original observation vector, which is then compared with the prediction vector composed of simulated pore water pressure values ​​and simulated surrounding rock deformation values ​​output by the model at that time point. If the difference between a sensor's observation value and the predicted value exceeds a preset fault threshold, the sensor data is discarded and replaced with historical observation values, resulting in a cleaned observation vector. The filter calculates the Kalman gain matrix for each ensemble member, and uses the gain matrix to correct the predicted values ​​of the state variables, obtaining updated ensemble mean and covariance. The updated unloading modulus, expansion coefficient, and erosion rate are then re-introduced into the dynamic fracture aperture evolution model to complete the rolling optimization of twin parameters within the current time window. Preferably, the assimilation time window of the ensemble Kalman filter is consistent with the duration of one tunnel excavation cycle, and the number of ensemble members is dynamically adjusted to the square root of the current fracture network node number, thus balancing computational efficiency and assimilation accuracy.

[0010] After completing the twin parameter rolling optimization, this invention implements an advanced seepage field simulation for the next tunneling cycle: the virtual excavation advance length and the final spatial coordinates of the tunnel face in the next cycle are used as boundary conditions input to the already parameter-assimilated dynamic fracture aperture evolution model; within each micro-hour step of the virtual excavation, the dynamic hydraulic aperture of each discrete fracture is updated sequentially by calling the four-fold driving force calculation function, and the topology of the fracture map is updated simultaneously by calling the graph neural network model; using the pressure head of the fracture network intersection nodes as the solution variable, the seepage field solution based on the finite volume method is performed on the updated fracture map to obtain the pore pressure cloud map and seepage velocity vector field at the moment the virtual excavation advance is completed. Subsequently, all connected paths in the seepage velocity vector field that exceed the preset critical flow velocity for water inrush are searched, and the discrete fracture sets traversed by the path are marked as potential water and mud inrush channels. The simulation results can reveal in advance the fracture combinations that may form concentrated seepage in front of the tunnel face, providing clear guidance for advanced treatment.

[0011] As a preferred technical solution of the present invention, after identifying potential water and mud inrush channels, an adaptive control and early warning parameter output step is further performed. All discrete fracture nodes within the virtual excavation area in front of the working face of the potential channel are extracted, and a three-dimensional high-brightness early warning cloud map of the area is generated and displayed as a red transparent body in the digital twin visualization interface. Based on the seepage velocity vector direction of each potential channel, the source position of the channel is traced in reverse, and its spatial coordinates are used as the target point for advanced grouting reinforcement. Centered on the target point, and combined with the empirical value of the diffusion radius of the grout in the fracture network, a minimum convex hull geometry covering all target points is automatically calculated and output as the optimal advanced grouting reinforcement range. Simultaneously, the total water pressure value within the upstream water storage structure of the potential water and mud inrush channel is calculated. When the total water pressure value exceeds a preset pressure relief threshold, a set of pressure relief and drainage hole layout coordinates and hole diameter parameters radially distributed on the working face contour line are output. Through the above steps, the present invention achieves integrated closed-loop decision-making from seepage field prediction and risk identification to automatic generation of control schemes.

[0012] To ensure the reliability of long-term predictions, this invention also introduces a closed-loop iterative process of prediction-verification-re-assimilation. The pore pressure cloud map and seepage velocity vector field at the moment the virtual excavation advance is completed are stored as the prediction results of this advanced simulation. After the actual excavation advance is completed, the measured pore water pressure values ​​from the deep-hole piezometer are collected and compared with the pore pressure values ​​at the corresponding locations in the prediction results to calculate the residual matrix. The residual matrix is ​​then used as additional observation information and input again into an ensemble Kalman filter to perform a secondary assimilation update on the unloading modulus, expansion coefficient, and erosion rate in the current twin, resulting in a corrected parameter set. This corrected parameter set will serve as the initial parameters for the advanced simulation of the next excavation cycle. Thus, through a continuous cycle of measurement-prediction-feedback-correction, the digital twin learns and adapts as the tunnel face advances, continuously approximating the actual mechanical behavior of the surrounding rock, significantly improving the accuracy and timeliness of the evolution and simulation of the seepage field in soft rock tunnels.

[0013] The beneficial effects of this invention are:

[0014] A dynamic fracture aperture evolution model was established, controlled by four mechanisms: tunneling unloading, support reaction force, water absorption and expansion, and hydraulic scour. The hydraulic aperture of each discrete fracture was defined as a spatiotemporally related dynamic variable, which is the sum of the initial hydraulic aperture and the aperture increment generated by the four driving forces. The excavation unloading opening amount was calculated based on the spatial distance between the tunnel face and the fracture, which was used to calculate the normal tensile displacement caused by the redistribution of surrounding rock stress and was gradually accumulated during the advance. The support reaction force closure amount was calculated by combining the dynamic support force applied to the surrounding rock surface by the steel arch frame and shotcrete, and using a semi-infinite plane force model to solve for the normal compressive displacement of the fracture surface. The water absorption and expansion closure amount was calculated by using the local water content transmitted in real time by a miniature multi-parameter sensor to calculate the squeezing displacement on both sides of the fracture caused by the volume expansion of soft rock. The hydraulic mudification expansion amount was calculated by using the peeling rate, which is related to the cube of the local seepage velocity, to calculate the aperture increment caused by the loss of particles from the fracture wall under piping erosion, and the calculation was stopped when the total amount of erodible material was exhausted. The four mechanisms work in parallel at each time step, achieving a dynamic approximation of the combined effects of tension, compression, expansion and squeezing, and scouring on fracture aperture. This avoids system inaccuracies caused by relying on only a single or simplified mechanism, making the aperture evolution process more consistent with the actual response of soft rock tunnels.

[0015] Based on the spatial advancement coordinates of the tunnel face, a graph neural network is used to reshape the fracture network graph, converting the three-dimensional discrete fracture network into an edge-node graph data structure. Each fracture is mapped to an edge, and fracture intersections are mapped to nodes. The edge feature vectors contain state information such as dynamic hydraulic aperture and permeability coefficient. A multi-layer graph attention model takes the tunnel face coordinates as input and performs weighted aggregation of the features of adjacent nodes and edges. Finally, the fracture state is determined based on the updated edge features: when the hydraulic aperture of an edge is lower than a preset closure threshold, it is marked as open; when the aperture exceeds a preset widening threshold, it is marked as a thickened dominant seepage channel. Opening or thickening events trigger the redistribution calculation of pressure head of adjacent nodes. This dynamic topology reshaping mechanism enables the fracture network connection relationship on which seepage calculations depend to be adjusted in real time according to the changes in fracture aperture caused by unloading, support, expansion and scouring. It accurately captures the formation of dominant channels and the exit of failed channels, ensuring that the seepage field solution grid is synchronized with the current physical seepage structure, thereby improving the ability to predict the location and scale of water and mud inrush. Attached Figure Description

[0016] The invention will now be further described with reference to the accompanying drawings.

[0017] Figure 1 This is a flowchart of a three-dimensional numerical simulation method for the evolution of seepage field in soft rock tunnels based on digital twins;

[0018] Figure 2 This is a flowchart of a method for predicting the dynamic evolution of discrete fracture networks based on graph neural networks;

[0019] Figure 3 This is a flowchart of the Kalman filter assimilation process for the parameter set of the dynamic fracture aperture evolution model.

[0020] Figure 4 This is a flowchart of a method for advanced simulation and early warning of water and mud inrush in soft rock tunnels based on digital twins. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] See Figure 1 This invention provides a three-dimensional numerical simulation method for the evolution of seepage field in soft rock tunnels based on digital twins, comprising: constructing a follow-up in-situ sensing network integrated in the initial support and advanced pipe roof; combining the morphology of the preceding fault zone obtained from the tunnel's advanced geological prediction to generate an initial digital twin containing a three-dimensional discrete fracture network; establishing a dynamic fracture aperture evolution model in the twin controlled by four mechanisms: tunneling unloading, support reaction force, water absorption expansion, and hydraulic scouring; using graph neural networks to reshape the fracture network map based on the spatial advancement coordinates of the tunnel face, and dynamically updating the connection state of the edges; using the monitoring data transmitted in real time from the steel arch fiber optic and deep hole piezometers, asynchronously assimilating the twin parameters through ensemble Kalman filtering, and finally predicting the extreme values ​​of the seepage field and the breakthrough path for the next tunneling cycle.

[0023] In practical implementation, the process of constructing a follow-up in-situ sensing network integrated into the initial support and advanced pipe roof is as follows: Fiber optic strain sensors are implanted along the longitudinal main reinforcement direction at the flange and web positions of each steel arch frame, forming an independent strain monitoring section on each steel arch frame. Simultaneously, the fiber optic strain sensors are bound to the inner wall of the steel pipe of the advanced pipe roof. As the tunnel face advances, the fiber optic strain sensors move synchronously with the steel arch frame and the advanced pipe roof, forming a distributed strain sensing array. The fiber optic strain sensors are demodulated using wavelength scanning optical frequency domain reflection technology. During demodulation, a scanning laser source injects a linear frequency modulated optical signal into the fiber optic strain sensor. The backscattered Rayleigh signals generated at various positions inside the fiber optic strain sensor are returned to the demodulator. Through Fourier transform, the frequency domain interference signal is converted to the spatial domain, obtaining a continuous strain distribution along the entire length of the fiber optic strain sensor, with a spatial resolution on the order of 0.1 millimeters.

[0024] Miniature multi-parameter sensors, consisting of a pressure sensor and a moisture content sensor, are installed in the sidewall of the completed support section behind the tunnel face and in the advance exploratory borehole ahead of the tunnel face. In the sidewall section behind the tunnel face, the miniature multi-parameter sensors are embedded at the contact surface between the shotcrete layer and the surrounding rock, with the pressure sensor's sensing surface facing the surrounding rock and the moisture content sensor's probe inserted into the surrounding rock. In the advance exploratory borehole ahead of the tunnel face, the miniature multi-parameter sensors are evenly spaced along the borehole depth, with a spacing of 0.5 meters between adjacent sensors. The depth of the advance exploratory borehole is determined based on the distance to the anomaly zone initially identified through tunnel advance geological prediction. The miniature multi-parameter sensors are fixed to the borehole wall within the borehole using airbag-type packers to ensure close contact between the pressure sensor and the moisture content sensor and the borehole wall rock mass. The miniature multi-parameter sensors integrate a temperature compensation module to perform real-time temperature correction on the raw readings of the pressure sensor and the moisture content sensor.

[0025] This study utilizes seismic reflection and ground-penetrating radar (GPR) methods in advanced tunnel geological prediction to jointly invert the geometry and spatial distribution of fault fracture zones and primary fracture groups ahead of the tunnel face. In practice, the seismic reflection method employs multiple seismic source points and receivers spaced evenly at the tunnel sidewalls. Elastic waves generated by the seismic source points propagate towards the tunnel face, are reflected at impedance interfaces, and the receivers record the travel time and amplitude of the reflected waves. The GPR method transmits high-frequency electromagnetic pulses towards the tunnel face via a transmitting antenna, and a receiving antenna receives the reflected signals from interfaces with different dielectric constants. The P-wave migration profile obtained from the seismic reflection method and the time profile obtained from the GPR method are jointly inverted. The joint inversion objective function includes both seismic travel time residuals and electromagnetic wave envelope differences, and a cross-gradient constraint is introduced to ensure consistency in interface morphology between the two physical field inversion results. The inversion outputs the spatial coordinates of the fault fracture zone boundary control points and the endpoint coordinates and attitude parameters of each fracture trace in the primary fracture group.

[0026] The coordinates of the measuring points of the fiber optic strain sensor and the miniature multi-parameter sensor, along with the inverted geometric parameters of the fracture group, are mapped together to the same three-dimensional spatial coordinate system. In this three-dimensional spatial coordinate system, the tunnel design centerline mileage direction is used as... The axis is defined by the center point of the tunnel face as the origin, and the horizontal and vertical midline directions are defined by... The axis, vertical direction is A right-handed rectangular coordinate system was established. The coordinates of the fiber optic strain sensor measuring points were calculated from the total station measurements taken during the installation and positioning of the steel arch frame. The spatial trajectory of the fiber optic strain sensor on each steel arch frame was fitted with a cubic spline curve and discretized into sampling points spaced 0.05 meters apart. The coordinates of the measuring points of the miniature multi-parameter sensor were obtained by accurately measuring the center position of its sensing surface using a total station during installation. The trace of each crack obtained through inversion was expressed in the form of a line segment or polysegment, and the endpoint coordinates were directly provided by the joint inversion results. Within the same three-dimensional spatial coordinate system, a three-dimensional discrete fracture network is generated based on the coordinates of the fracture endpoints. Each discrete fracture is assigned an initial value for fracture aperture, strike, dip angle, and roughness attributes. The initial value for fracture aperture is estimated based on the residual aperture of the fracture under confining pressure conditions obtained during inversion. The strike and dip angle are determined by the angle between the fracture trace direction vector and the coordinate axis. Roughness is expressed using the joint roughness coefficient JRC value, ranging from 0 to 20, with the specific value determined based on lithological analogy of the surrounding rock or borehole image analysis. The three-dimensional discrete fracture network and the surrounding rock medium together constitute the initial digital twin. The surrounding rock medium fills the space outside the fracture network with continuum units and is assigned parameters such as the rock mass elastic modulus, Poisson's ratio, and permeability tensor, thus completing the construction of the initial digital twin.

[0027] In practical implementation, the hydraulic aperture of each discrete fracture in the three-dimensional discrete fracture network is defined as a dynamic variable that varies with time and spatial coordinates. This dynamic variable is formed by the superposition of the initial hydraulic aperture and the aperture increment generated by the four driving forces. The superposition method is expressed as follows:

[0028]

[0029] in, Indicates time Time and spatial coordinates The dynamic hydraulic aperture of the discrete fracture at the location, in meters; Represents spatial coordinates The initial hydraulic aperture of the discrete fracture at the location is assigned by the joint inversion results when the initial digital twin is constructed, specifically the residual aperture value estimated under confining pressure conditions; This represents the amount of excavation unloading and opening caused by the first driving force, i.e., in time. Time and spatial coordinates The normal tensile displacement of the fracture surface caused by the redistribution of stress in the surrounding rock; This represents the closing amount of the second driving force support reaction force, i.e., in time... Time and spatial coordinates The compressive displacement generated in the normal direction of the crack surface by the dynamic support force applied by the steel arch frame and shotcrete; This represents the third driving force, the amount of water absorption and expansion closure, i.e., in time... Time and spatial coordinates The displacement caused by the compression of the two sides of the fissure due to the volume expansion of soft rock after it comes into contact with water; This represents the fourth driving force, hydraulic siltation, which expands capacity over time. Time and spatial coordinates The increased fracture aperture value due to the loss of soft rock particles caused by piping erosion. The sign of the aperture increment generated by the four driving forces is defined with the fracture opening direction as positive, where... It is a positive value. Negative value Negative value It is a positive value.

[0030] For the first driving force excavation unloading opening, the normal tensile displacement of the fracture surface caused by the redistribution of surrounding rock stress is calculated based on the three-dimensional spatial distance between the current tunnel face and the location of the discrete fractures. In practice, at the end of each tunneling cycle, the spatial coordinates of the tunnel face after advancement are read, and the three-dimensional spatial distance between the tunnel face and the center point of each discrete fracture in the three-dimensional discrete fracture network is calculated. The three-dimensional spatial distance is substituted into a pre-calibrated unloading displacement attenuation function. The form of the unloading displacement attenuation function is constructed based on the analytical solution of the excavation disturbance stress field in the elastic half-space, and the function value decreases negatively exponentially with the increase of the three-dimensional spatial distance. When the three-dimensional spatial distance of the discrete fractures is less than the unloading influence radius, the normal tensile displacement of the fracture surface is accumulated according to the unloading displacement attenuation function. When the three-dimensional spatial distance is greater than the unloading influence radius, the accumulated amount of normal tensile displacement of the fracture surface is zero. The unloading influence radius is taken as three times the tunnel excavation diameter, and the attenuation coefficient in the unloading displacement attenuation function is determined by back-calculation using the surrounding rock convergence monitoring data in the early stage of tunnel excavation. The excavation unloading opening amount gradually increases as the tunnel face advances. The excavation unloading opening amount obtained by each discrete crack at the end of the current tunneling cycle is equal to the newly added normal tensile displacement value since the end of the previous tunneling cycle.

[0031] For the second driving force support reaction closure amount, the compressive displacement generated by the dynamic support force in the normal direction of the fracture surface is calculated based on the dynamic support force applied to the surrounding rock surface by the steel arch frame and shotcrete. In specific implementation, the real-time strain data of the steel arch frame after installation is read from the fiber optic strain sensor embedded in each steel arch frame. Based on the elastic modulus and moment of inertia of the steel arch frame, the real-time strain data is converted into the radial support force of the steel arch frame on the surrounding rock. The strength growth curve of the shotcrete is obtained from the mix proportion and age information of the shotcrete. The strength growth curve is established based on the uniaxial compressive strength test data of the reserved test blocks at the same age. Based on the strength growth curve, the bond force and surface constraint force of the shotcrete on the surrounding rock at the current age are calculated. The radial support force and surface constraint force are combined into the dynamic support resultant force applied to the surface of the surrounding rock. To obtain the spatial distance between the location of the discrete crack and the adjacent steel arch, a semi-infinite plane force model in elasticity is adopted. The resultant force of the dynamic support is used as the surface force boundary condition acting on the surrounding rock surface. The normal displacement solution of the semi-infinite plane force model at the location of the discrete crack is solved to obtain the compressive displacement. The compressive displacement is opposite in direction to the tensile displacement, and the compressive displacement takes a negative value when superimposed.

[0032] For the third driving force, water absorption and expansion closure, the compressive displacement on both sides of the fracture is calculated based on the local water content data transmitted in real time by the micro multi-parameter sensor. In specific implementation, the local water content time series of the surrounding rock is obtained from the water content sensor in the micro multi-parameter sensor, and the local water content time series is compared with the initial water content of the soft rock to obtain the water content increment. The initial water content of the soft rock is obtained after the initial installation of the micro multi-parameter sensor. According to the expansion rate-water content relationship curve of the soft rock, the water content increment is converted into the volume expansion rate of the soft rock. The expansion rate-water content relationship curve of the soft rock is pre-determined through indoor expansion rate test of soft rock samples in the tunnel. In the test, multiple groups of samples with different initial water contents were set up, and the change of volume expansion rate with water content during the immersion saturation process was recorded. When converting the volume expansion rate into the compressive displacement in the normal direction of the fracture surface, half of the fracture spacing is taken as the equivalent thickness of the expansion deformation, and the compressive displacement is obtained by multiplying the volume expansion rate by the equivalent thickness. The extrusion displacement and the tensile displacement are in opposite directions, and the extrusion displacement takes a negative value when superimposed.

[0033] For the fourth driving force, hydraulic mudification, when the hydraulic gradient of seepage within the fracture exceeds the critical initiation value of soft rock particles, the piping erosion module is activated. This module calculates the mass of soft rock particles stripped from the fracture wall and converts this stripped mass into the increased fracture aperture due to particle loss. The critical initiation value of soft rock particles is determined based on the particle size, density, and roughness coefficient of the fracture wall, using a modified Shields critical shear stress formula. In practice, the piping erosion module calculates the local hydraulic gradient within each discrete fracture at each time step. When the seepage shear stress corresponding to the local hydraulic gradient exceeds the critical initiation value of soft rock particles, particle stripping calculation is activated. The stripping rate of soft rock particles is proportional to the cube of the local seepage velocity. The stripping rate constant is calibrated through an indoor small-scale piping model test. This test uses soft rock remolded samples from the tunnel, and the particle stripping mass per unit time is measured under set seepage velocity conditions. The stripping rate constant is obtained by fitting the cube of the seepage velocity as the independent variable. The capacity expansion calculation stops when the stripping quality reaches the upper limit of the total erodible area of ​​the fracture wall. The total erodible area of ​​the fracture wall is estimated based on the thickness of the loose particle layer and the surface area of ​​the fracture wall. The thickness of the loose particle layer is determined by the depth of the loosened zone of the fracture wall in the borehole image, and the surface area of ​​the fracture wall is calculated by multiplying the fracture trace length by the element thickness. After obtaining the stripping quality, the stripping quality is divided by the skeleton density of the soft rock particles and the surface area of ​​the fracture wall to calculate the increase in fracture aperture due to particle loss.

[0034] In specific implementation, please refer to Figure 2 The initial fracture graph topology is formed by mapping each discrete fracture in the 3D discrete fracture network to an edge in a graph data structure, and mapping the intersection points between discrete fractures or the contact surfaces between discrete fractures and the surrounding rock matrix to nodes in the graph data structure. The node definition rules for the graph data structure are as follows: when two or more discrete fractures intersect at the same spatial location, a node is generated at the intersection point; when a discrete fracture terminates inside the surrounding rock matrix and does not intersect with any other discrete fractures, an end node is generated at the terminating endpoint of the discrete fracture. The edge definition rules for the graph data structure are that each discrete fracture is treated as an edge, connecting the nodes at both ends of the discrete fracture. After mapping all discrete fractures in the 3D discrete fracture network according to the above rules, an initial fracture graph topology consisting of a set of nodes and a set of edges is obtained. The adjacency matrix of the graph topology is constructed based on the connection relationships between nodes and edges. An element value of 1 in the adjacency matrix indicates that there is an edge between two nodes, and an element value of 0 indicates that there is no connection between two nodes.

[0035] Assign an initial feature vector to each node and each edge in the graph data structure. The feature vector of a node includes the node's spatial coordinates and pressure head; the spatial coordinates are taken from the node's position in a three-dimensional spatial coordinate system. The coordinate values ​​and pressure head are taken from the initial pore water pressure head value at the node location. The initial pore water pressure head value is determined based on the groundwater level and hydrostatic pressure distribution in the tunnel area. If there is no groundwater, the pressure head is set to zero. The feature vector of the edge includes the dynamic hydraulic aperture of the discrete fracture, the current permeability coefficient, and the water content of the soft rock on both sides. The dynamic hydraulic aperture is calculated and output by the dynamic fracture aperture evolution model at the current moment. The current permeability coefficient is obtained by converting the dynamic hydraulic aperture through the cubic law. The water content of the soft rock on both sides is taken as the average of the water content readings of the two micro multi-parameter sensors at the two nodes of the discrete fracture. If the discrete fracture exceeds the deployment range of the micro multi-parameter sensors, the water content is taken as the extrapolated value of the reading of the nearest adjacent sensor. The three components of the feature vector of the edge—dynamic hydraulic aperture, current permeability coefficient, and water content—are continuously updated during the model operation as the tunneling cycle progresses.

[0036] The 3D spatial coordinates of the tunnel face after the current tunneling cycle ends are used as external input parameters and fed into a multilayer graph neural network model. The multilayer graph neural network model consists of three graph attention layers and one edge state prediction layer stacked sequentially. The first graph attention layer receives the node feature vectors and edge feature vectors of the initial fracture graph topology as input, aggregates the feature vectors of all neighboring nodes and connecting edges in the first-order neighborhood of each node, and outputs an updated node feature vector. The second graph attention layer receives the updated node feature vectors from the first graph attention layer and the current edge feature vectors, performs neighborhood aggregation again, and outputs a further updated node feature vector. The third graph attention layer receives the updated node feature vectors from the second graph attention layer and the current edge feature vectors, completes neighborhood aggregation, and outputs the final updated node feature vector. The edge state prediction layer receives the updated node feature vectors corresponding to the two ends of each edge output by the third graph attention layer, as well as the edge feature vector of that edge, and outputs a state determination result for each edge through a multilayer perceptron. The multilayer perceptron with edge state prediction layer contains two hidden layers and one output layer. Each of the two hidden layers has 128 neurons and uses the ReLU activation function. The output layer has 1 neuron and uses the Sigmoid activation function. The output value represents the probability that the edge remains connected, and the value ranges from 0 to 1.

[0037] In the multi-layer graph neural network model, during the weighted aggregation of features of adjacent nodes and edges by each graph attention layer, the attention coefficient is dynamically determined by the relative distance between the current face advancement coordinates and the spatial orientation of the fracture. Nodes in the attention layer of a graph Its neighboring nodes Attention coefficient between The calculation method is expressed as follows:

[0038]

[0039] in, This is the sequence number of the attention layer in the graph, with a value of 1, 2, or 3. Represents a node In the The node feature vector output by the layer, when hour For nodes The initial feature vector; Indicates the connection node With nodes The edge at the 1st The edge feature vector output by the layer, when hour Let be the initial feature vector of the edge; For the first The learnable weight matrix that acts on the feature vector of a node in the attention layer of the graph has a dimension of 128×128. For the first The learnable weight matrix that acts on the edge feature vector in the graph attention layer has a dimension of 128×128. For the first The learnable weight vector used to calculate the attention score in the graph attention layer has a dimension of 128. Represents a node With nodes The three-dimensional spatial distance between the spatial center point of the connected discrete fracture and the three-dimensional spatial coordinates of the tunnel face after the current tunneling cycle ends, in meters; For the first The scaling factor of spatial distance in the attention layer of each graph is determined through backpropagation optimization during model training, with an initial value set to 0.1. The spatial distance attenuation index is set to 2, which is determined based on the negative square law of the attenuation of surrounding rock stress disturbance with distance. This is the Softmax activation function, used to normalize the attention coefficients; It is a linear rectified activation function with leakage, and the negative slope is set to 0.2.

[0040] During the iterative aggregation and updating of the feature vectors of each edge in the initial fracture graph topology by the graph neural network model, the three-layer graph attention layer updates the node feature vectors sequentially, and the edge state prediction layer determines the state of each edge. When the dynamic hydraulic aperture value of an edge converges to zero or below the preset closure threshold after model calculation, the graph neural network model marks the state of that edge in the graph as open, indicating that the discrete fracture has been completely closed and water flow cannot pass through. The preset closure threshold is one-hundredth of the minimum initial hydraulic aperture value, that is, when the dynamic hydraulic aperture decreases to less than one-hundredth of the initial value, the discrete fracture is determined to be closed. When the dynamic hydraulic aperture value of an edge exceeds the preset widening threshold after model calculation, the graph neural network model marks the state of that edge in the graph as bold, indicating that the discrete fracture has been scourned into a dominant seepage channel. The preset widening threshold is five times the maximum initial hydraulic aperture value, that is, when the dynamic hydraulic aperture increases to more than five times the initial value, the discrete fracture is determined to be scoured into a dominant seepage channel. The minimum and maximum initial hydraulic opening values ​​were obtained statistically from the initial hydraulic opening values ​​of all discrete fractures in the three-dimensional discrete fracture network.

[0041] The graph neural network model outputs an updated fracture graph topology, which includes state markers for all edges. These markers are categorized into three types: normally connected, disconnected, and bolded. Changes in the state of disconnected or bolded edges trigger a redistribution of pressure head features at adjacent nodes. In practice, when an edge's state changes from normally connected to disconnected, the adjacency matrix element corresponding to that edge in the graph data structure is changed from 1 to 0. The seepage flow on all other edges connected to the nodes at both ends of that edge is recalculated, and the pressure head values ​​at the two ends are corrected based on the new connectivity. When an edge's state changes from normally connected to bolded, the current permeability coefficient in the edge feature vector of that edge in the graph data structure is recalculated using the cubic law based on the bolded dynamic hydraulic aperture. The pressure head values ​​at adjacent nodes are redistributed based on the updated current permeability coefficient. The pressure head redistribution calculation employs a graph-based steady-state seepage solver. After the graph topology is updated, a constant head boundary condition is used as the external driver to solve for the pressure head values ​​at all nodes, ensuring that the net flow at each node is zero according to the continuity equation.

[0042] During the training phase of the multilayer graph neural network model, training data is constructed from several completed tunneling cycles. The input to the training samples includes the 3D spatial coordinates of the tunnel face at the end of each historical tunneling cycle, as well as the fracture topology and feature vector at the end of that cycle. The labels of the training samples are the actual states of each edge confirmed by measured data after the next tunneling cycle. The loss function is cross-entropy loss, the optimizer is the Adam optimizer, the initial learning rate is set to 0.001, and the training epochs are set to 500. In each epoch, a complete forward and backward propagation is performed on all training samples. After training, the weight parameters of the multilayer graph neural network model are fixed and used for inference and prediction in subsequent tunneling cycles.

[0043] In specific implementation, please refer to Figure 3 Within each assimilation time window, the surrounding rock deformation pressure values ​​at each measuring point are collected from distributed fiber optic sensors on the steel arch frame, and the pore water pressure values ​​at each deep-hole measuring point are collected from the permeability sensor in the miniature multi-parameter sensor. These two types of measured data form an observation vector. The dimension of the observation vector is the total number of all valid measuring points within the current assimilation time window, with each dimension corresponding to the measured value of one measuring point. The surrounding rock deformation pressure value is collected as follows: from the strain distribution demodulated by the fiber optic strain sensor on each steel arch frame, the strain values ​​at each measuring point at the flange and web positions of the steel arch frame are extracted, and the strain value is multiplied by the elastic modulus of the steel arch frame to obtain the surrounding rock deformation pressure value at that measuring point. The pore water pressure value is directly output by the permeability sensor and corrected by the temperature compensation module built into the miniature multi-parameter sensor. At the end of the assimilation time window, the surrounding rock deformation pressure values ​​and pore water pressure values ​​of all measuring points within that time window are summarized, arranged in order of measuring point number, and form an observation vector.

[0044] The three uncertain parameters in the dynamic fracture aperture evolution model—unloading modulus, soft rock expansion coefficient, and piping erosion rate—are treated as state variables to be assimilated. The unloading modulus is the attenuation coefficient in the unloading displacement attenuation function used in the calculation of the unloading opening amount during the first driving force excavation. Its prior distribution follows a log-normal distribution, and the logarithmic mean is the attenuation coefficient value obtained from the back-calculation at the initial stage of tunnel excavation, with a logarithmic standard deviation of 0.3. The soft rock expansion coefficient is the slope parameter of the soft rock expansion rate-water content relationship curve in the calculation of the water absorption expansion closure amount during the third driving force. Its prior distribution follows a normal distribution, and the mean is the slope value obtained from the indoor expansion rate test fitting, with a standard deviation of 10% of the mean. The piping erosion rate is the proportionality coefficient of the relationship between the stripping rate constant and the cubic relationship of the local seepage velocity in the calculation of the hydraulic mudification expansion capacity during the fourth driving force. Its prior distribution follows a log-normal distribution, and the logarithmic mean is the logarithm of the indoor piping model test calibration value, with a logarithmic standard deviation of 0.5. For each state variable, an initial set that follows its prior distribution is generated. Specifically, based on the prior distribution type of the state variable, the Latin hypercube sampling method is used to perform equal probability interval sampling on the cumulative distribution function to generate a predetermined number of initial values ​​for the set members. The set members of the three state variables are combined one-to-one to form the state variable set.

[0045] The number of set members is dynamically adjusted based on the current number of fracture network nodes, rounded down to the square root of the node count. The current number of fracture network nodes is obtained from the total number of nodes in the initial fracture graph topology. At the beginning of each assimilation time window, the total number of nodes in the 3D discrete fracture network graph is counted, the square root of the total number of nodes is calculated, and the square root value is rounded down to obtain the number of set members. If the number after rounding down is less than 10, the number of set members is set to 10. Each set member is represented by a vector containing three components: unloading modulus, soft rock expansion coefficient, and piping erosion rate.

[0046] Before inputting the observation vector into the ensemble Kalman filter, a time window alignment operation is performed between the observation data and the model prediction. From the raw data stream acquired by the steel arch fiber optic sensor and deep-hole piezometer, a set of observations at the same physical time point is extracted to form the original observation vector. The same physical time point refers to the end of the current assimilation time window. The simulated pore water pressure value and simulated surrounding rock deformation value output by the dynamic fracture aperture evolution model at this physical time point constitute the prediction vector for that time point. When the dynamic fracture aperture evolution model is running, the unloading modulus, soft rock expansion coefficient, and piping erosion rate updated after the end of the previous assimilation time window are used as model parameters to simulate the surrounding rock deformation and seepage process from the end of the previous assimilation time window to the end of the current assimilation time window. The simulated pore water pressure value and simulated surrounding rock deformation value corresponding to each measuring point position in the original observation vector are extracted from the model and arranged in the same order to form the prediction vector.

[0047] The difference vector between the original observation vector and the predicted vector is calculated. Each element in the difference vector is the absolute value of the measured value minus the simulated value at the corresponding measuring point. When the difference vector for a measuring point exceeds a preset sensor fault threshold, the data of the corresponding sensor in the original observation vector is removed. The preset sensor fault threshold is set to 50% of the full scale of the steel arch fiber optic strain sensor for rock deformation pressure value sensors and 50% of the full scale of the seepage pressure sensor for pore water pressure value sensors. The removed sensor data is replaced with the historical observation value at the corresponding position in the observation vector cleaned by the assimilation time window on that sensor. If a sensor is judged to be faulty for three consecutive assimilation time windows, the observation dimension corresponding to that sensor is permanently removed in subsequent assimilation time windows. After removal and replacement, the cleaned observation vector is obtained.

[0048] The cleaned observation vectors are input into an ensemble Kalman filter. The filter calculates the Kalman gain matrix for each ensemble member, and uses the gain matrix to correct the predicted values ​​of the state variables. In practice, the predicted values ​​of the ensemble Kalman filter are obtained by integrating the set of state variables updated in the previous assimilation time window forward through the dynamic fracture aperture evolution model to the end of the current assimilation time window. For each ensemble member, the model is run using the state variable values ​​corresponding to that member to obtain the simulated values ​​at the corresponding measurement point locations. The simulated values ​​of all ensemble members form a simulation set. Based on the simulation set, the covariance matrix of the simulated value set, the cross-covariance matrix between the simulated values ​​and the state variables, and the observation error covariance matrix are calculated. The observation error covariance matrix is ​​in diagonal form, with the diagonal elements being the measurement noise variance calibrated by the sensors at each measurement point. The Kalman gain matrix is ​​calculated using the standard perturbation observation method of the ensemble Kalman filter. The gain matrix is ​​obtained by inverting the outlier covariance of the simulation set and the state variable set combined with the observation error covariance matrix. The predicted state variables of each set member are corrected using the Kalman gain matrix to obtain the updated set of state variables. The correction process is described as follows:

[0049]

[0050] in, Indicates the first The state variable vectors of each set member after assimilation correction The value of is an integer ranging from 1 to the total number of members in the set; Indicates the first The predicted state variable vector of each set member before correction, that is, the value obtained by integrating the state variable after the previous assimilation window update forward through the dynamic fracture aperture evolution model to the end time of the current assimilation time window; This represents the Kalman gain matrix calculated from the ensemble Kalman filter, whose dimension is the number of state variables multiplied by the dimension of the observation vector; This represents the cleaned observation vector; The observation operator maps the state variable vector to the observation space and outputs a vector of simulated values ​​corresponding to the observation vector. In practical implementation... That is, with the first The simulated output of the dynamic fracture aperture evolution model at the corresponding measuring point when the state variables of each set member are used as model parameters.

[0051] After obtaining the updated set of state variables, the updated set mean and set covariance are calculated. The set mean is the vector average of all set member state variable vectors, and the set covariance is the sample covariance of the deviations of all set member state variable vectors from the set mean. The updated unloading modulus, expansion coefficient, and erosion rate are then re-substituted into the dynamic fracture aperture evolution model, replacing the original parameters, to complete the rolling optimization of the twin parameters within the current assimilation time window. The assimilation time window of the ensemble Kalman filter is consistent with the duration of one tunnel excavation cycle. The duration of one tunneling cycle is from the start of the current excavation to the start of the next cycle, including the time occupied by all procedures such as excavation, muck removal, and support.

[0052] In specific implementation, please refer to Figure 4 After completing the rolling optimization of the twin parameters within the current assimilation time window, the virtual excavation advance length of the next tunneling cycle and the final spatial coordinates of the tunnel face corresponding to the virtual excavation advance are input as boundary conditions into the dynamic fracture aperture evolution model that has completed parameter assimilation. The virtual excavation advance length is equal to the advance length of one tunneling cycle in the tunnel construction design, and the final spatial coordinates of the tunnel face are calculated by adding the virtual excavation advance length to the current three-dimensional spatial coordinates of the tunnel face along the tunnel excavation direction.

[0053] Within each micro-hour step of the virtual excavation, the dynamic hydraulic opening of each discrete fracture in the three-dimensional discrete fracture network is updated sequentially using the four driving force calculation functions, and the topology of the fracture map is updated simultaneously using the graph neural network model. The length of each micro-hour step is one-hundredth of the tunneling cycle time. The four driving force calculation functions include the excavation unloading opening function, the support reaction closure function, the water absorption expansion closure function, and the hydraulic mudification expansion function. At the beginning of each micro-hour step, the first driving force, excavation unloading opening, is calculated based on the virtual face advance distance corresponding to the current micro-hour step; the second driving force, support reaction closure, is calculated based on the support force state of the steel arch and shotcrete at the current virtual moment; the third driving force, water absorption expansion closure, is calculated based on the local moisture content data transmitted from the micro multi-parameter sensor; and the fourth driving force, hydraulic mudification expansion, is calculated based on the current seepage velocity within the fracture. The results of the four driving force calculations are superimposed on the dynamic hydraulic opening of each discrete fracture. After updating the dynamic hydraulic opening, the three-dimensional spatial coordinates of the tunnel face at the end of the current micro-hour step and the updated edge feature vectors are used as input. A graph neural network model is then invoked to determine the state of all edges in the fracture graph topology and update the connectivity status of the edges. This process is repeated in each micro-hour step of the virtual excavation until the entire virtual excavation length is completed.

[0054] Using the pressure head at each intersection node in a 3D discrete fracture network as the solution variable, a seepage field solution based on the finite volume method is performed on the updated fracture map. The updated fracture map topology is transformed into a finite volume mesh, and a control volume element is established at each node, with the interface of the control volume element corresponding to the edge connecting the node. A flow conservation equation is established for each control volume element, and the interface flow is calculated by Darcy's law based on the pressure head difference between the nodes on both sides of the interface and the current permeability coefficient of the edge. Boundary conditions are applied to the finite volume mesh, with the boundary condition at the working face set to zero pressure head, and the outer boundary of the model far from the working face set to a fixed pressure head, the value of which is determined according to the groundwater level. The discretized linear equations are solved iteratively using the preprocessed conjugate gradient method, with a convergence tolerance of 1.0 × 10⁻⁻⁻⁶. 6 The pressure head values ​​of all nodes at the moment of completion of the virtual excavation are obtained. The seepage velocity on each edge is calculated from the node pressure head values. The seepage velocity is determined by the node pressure head difference, the edge length, and the current permeability coefficient of the edge. Finally, a pore pressure cloud map and a seepage velocity vector field are formed. The pore pressure cloud map is rendered with the node coordinate positions as discrete points and the pressure head value as a color level variable. The seepage velocity vector field uses vector arrows at the midpoint of each edge to indicate the direction and magnitude of the seepage velocity.

[0055] The search engine searches for all connected paths in the seepage velocity vector field that exceed a preset critical velocity for water inrush. The preset critical velocity is set at 0.01 m / s, determined by a safety margin, based on the critical starting velocity for water to carry mud particles in soft rock tunnels, and combined with on-site hydrogeological tests and laboratory erosion tests. The algorithm for searching connected paths employs a graph-based depth-first search. Starting from any edge in the seepage velocity vector field where the seepage velocity exceeds the preset critical velocity, the algorithm traces adjacent edges that also exceed the preset critical velocity along the connection direction of the edge until it reaches the tunnel face boundary. The discrete set of fractures traversed during this tracing process is marked as a potential water and mud inrush channel. The depth-first search is repeated until all edges satisfying the velocity condition are classified into a potential water and mud inrush channel.

[0056] Extract all discrete fracture nodes located within the virtual excavation area in front of the tunnel face within potential water and mud inrush channels. The virtual excavation area refers to the spatial region extending from the current tunnel face in the tunnel excavation direction to the virtual excavation advance endpoint. This spatial range is formed by stretching the tunnel face contour along the tunnel axis to the virtual excavation advance length. After extracting the discrete fracture nodes falling within the virtual excavation area, a 3D highlighted early warning cloud map of this area is generated. The 3D highlighted early warning cloud map is constructed by assigning a red color level marker to the extracted discrete fracture nodes as the center, and highlighting the edges connected to the nodes that exceed the preset critical flow velocity for water inrush with red lines. In the visualization interface of the digital twin, the red transparent body is rendered and displayed. The transparency coefficient of the transparent body is set to 0.4, allowing the red transparent body to be superimposed with the gray semi-transparent rendering of the surrounding rock medium, highlighting the spatial relationship between the potential water and mud inrush channels and the tunnel face and surrounding rock.

[0057] Based on the seepage velocity vector direction of each potential water and mud inrush channel, the source location of the potential water and mud inrush channel is traced backward. The backward tracing process proceeds in the opposite direction of the seepage velocity vector field, starting from the exposed node at the working face, and tracing back segment by segment along the discrete fracture edges that make up the potential water and mud inrush channels. At each node, the seepage velocity direction of all edges flowing into that node is compared, and the inflow edge with the smallest angle to the mainstream direction is selected as the backward tracing path, until no edge in the seepage velocity vector field has a seepage velocity exceeding the preset critical flow velocity for water inrush. The spatial coordinates of the tracing endpoint are determined as the target point for pre-grouting reinforcement. After tracing the source of all potential water and mud inrush channels, a set of spatial coordinates of the target points is obtained.

[0058] Centered on the target points, and combining empirical values ​​of the grout diffusion radius in the fracture network, an automatic calculation and output of a minimum convex hull geometry covering all target points is used as the optimal pre-grouting reinforcement range. The empirical value of the grout diffusion radius in the fracture network is determined comprehensively based on the grouting material type, fracture aperture range, and grouting pressure. For ordinary cement grout in soft rock fractures with a fracture aperture of 0.1 mm to 1 mm, the empirical value of the diffusion radius is taken as 1.5 meters. The calculation of the minimum convex hull geometry uses the QuickHull 3D convex hull algorithm. The input is the spatial coordinates of all target points, and the output is a convex polyhedron composed of triangular facets. This convex polyhedron encloses all target points and has the smallest volume. The optimal pre-grouting reinforcement range is the spatial region enclosed by the convex hull geometry, and the output includes a list of vertex coordinates and a list of facet indices of the convex hull geometry.

[0059] Calculate the total water pressure within the reservoir structure upstream of the potential water and mud inrush channel. The total water pressure is calculated using the following formula:

[0060]

[0061] in, This represents the total water pressure within the water storage structure, expressed in kilonewtons. This indicates the specific gravity of water, and its value is 9.81 kN per cubic meter. This represents the total number of discrete triangular facets on the boundary surface of the water storage structure. The boundary of the water storage structure is defined by the closed area enclosed by all the edges connected to the upstream end of the potential water and mud inrush channels that do not exceed the preset critical flow velocity for water inrush. Indicates the first The pressure head value at the center point of each triangular facet, in meters, is read from the pore pressure cloud map obtained by solving the seepage field. Indicates the first The area of ​​each triangular facet, in square meters, is calculated from the cross product of the three vertices of the facet. When the total water pressure exceeds a preset pressure relief threshold, a set of coordinates and aperture parameters of pressure relief drainage holes radially distributed along the tunnel face contour line are output. The preset pressure relief threshold is 0.6 times the product of the minimum principal stress of the surrounding rock and the area of ​​the tunnel face. The minimum principal stress of the surrounding rock is obtained from in-situ stress testing. The number of pressure relief drainage holes is determined based on the ratio of the excess total water pressure to the preset pressure relief threshold, with a minimum of 3 holes. They are radially distributed with evenly spaced angular intervals along the tunnel face contour line. The aperture of the pressure relief drainage holes is uniformly set to 50 mm, and the length of each hole is the straight-line distance from the tunnel face to the boundary of the water storage structure. The coordinates of the pressure relief drainage holes are output as the coordinates of evenly distributed points along the tunnel face contour line.

[0062] In the calculation of the support reaction closure, real-time strain data of each steel arch frame after installation is read from fiber optic strain sensors embedded in the frame. The real-time strain data is demodulated using wavelength scanning optical frequency domain reflection technology to obtain the continuous strain distribution along the fiber optic strain sensor. Based on the elastic modulus and moment of inertia of the steel arch frame, the real-time strain data is converted into the radial support force of the steel arch frame on the surrounding rock. The conversion process is as follows: Strain values ​​are taken at measuring points on the flange and web of the steel arch frame; the axial stress at these points is calculated; the axial stress is multiplied by the cross-sectional area of ​​the steel arch frame to obtain the internal force at that location; the internal forces at each location are integrated along the circumference of the steel arch frame and then projected onto the radial direction to obtain the radial support force. The elastic modulus of the steel arch frame is taken as 206 GPa, and the moment of inertia is obtained from a table based on the cross-sectional geometric parameters of the steel arch frame model. The strength growth curve of shotcrete was obtained from the mix proportion and age information. The mix proportion information included cement type, water-cement ratio, and accelerator dosage. The age information was the time from the completion of shotcrete construction to the current virtual time. The strength growth curve was established using uniaxial compressive strength test data of concrete specimens with the same mix proportion as those on site at different ages. Based on the strength growth curve, the uniaxial compressive strength value of shotcrete at the current age was read. 10% of the uniaxial compressive strength value was taken as the bond force between the shotcrete and the surrounding rock. The bond force was integrated along the area of ​​the contact surface between the surrounding rock and the shotcrete to obtain the surface restraint force. The radial support force and the normal component of the surface restraint force on the surrounding rock surface were added together to form the dynamic support resultant force applied to the surrounding rock surface.

[0063] The spatial distance between the location of the discrete fracture and the nearest steel arch is obtained. The nearest steel arch is the one closest to the discrete fracture along the longitudinal direction of the tunnel. The spatial distance is the three-dimensional spatial distance between the three-dimensional spatial coordinates of the center point of the discrete fracture and the nearest point on the contour line of the adjacent steel arch at the same mileage. A semi-infinite plane force model from elasticity is adopted, treating the resultant force of the dynamic support as a concentrated force or a line-distributed force acting on the surface of the semi-infinite plane. Based on the magnitude of the resultant force of the dynamic support, the location of its point of application, and the spatial distance between the discrete fracture and the adjacent steel arch, the analytical solution of the normal displacement at any point within the semi-infinite plane is solved. For the case where the concentrated force acts on the surface of the semi-infinite plane, the normal displacement decays along the depth direction, and the decay function is obtained using the Boussinesq solution. Substituting the spatial coordinates of the point of application of the resultant force of the dynamic support and the discrete fracture into the Boussinesq solution, the compressive displacement value in the normal direction of the discrete fracture surface is calculated. The direction of the compressive displacement is opposite to the direction of the tensile displacement and serves as the second driving force closure quantity in the superposition of the dynamic hydraulic opening.

[0064] After obtaining the pore pressure cloud map and seepage velocity vector field at the moment of completion of the virtual excavation, these data are stored as the prediction results of this advanced simulation. The stored content includes the spatial coordinates of all nodes, the pressure head value, and the seepage velocity vector of all edges. Once the actual excavation length corresponding to the virtual excavation is completed during the actual tunneling process, the measured pore water pressure value from the deep-hole piezometer at the moment of actual tunneling completion is collected. The deep-hole piezometer is a pore pressure sensor among miniature multi-parameter sensors; the measured pore water pressure value is obtained after temperature compensation and data cleaning. The measured pore water pressure value is compared with the pore pressure value at the corresponding location in the prediction results. Corresponding locations are matched according to the nearest neighbor principle of spatial coordinates, and a residual matrix is ​​calculated. The elements of the residual matrix are the difference between the measured pore water pressure value and the predicted pore pressure value at each measuring point.

[0065] The residual matrix is ​​then input again into the ensemble Kalman filter as additional observation information to perform a secondary assimilation update on the unloading modulus, expansion coefficient, and erosion rate in the current twin. The secondary assimilation update is performed by adding a dimension corresponding to the residual matrix to the observation vector of the ensemble Kalman filter. The newly added part of the observation error covariance matrix uses the sample variance of the residual as its diagonal element. A complete ensemble Kalman filter assimilation process is then executed to obtain the corrected set of parameters for the unloading modulus, expansion coefficient, and erosion rate. This corrected set of parameters replaces the rolling optimization set and serves as the initial set for the next tunneling cycle's advance simulation. The advance simulation for the next tunneling cycle is executed again at the beginning of the next assimilation time window.

[0066] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins, characterized in that, include: A follow-up in-situ sensing network integrated into the initial support and advanced pipe roof is constructed. Combined with the morphology of the forward fault zone obtained from the tunnel's advanced geological prediction, an initial digital twin containing a three-dimensional discrete fracture network is generated. A dynamic fracture aperture evolution model was established in the twin, controlled by four mechanisms: tunneling unloading, support reaction force, water absorption expansion, and hydraulic scouring. Based on the spatial advancement coordinates of the tunnel face, a graph neural network is used to reshape the topology of the fracture network graph and dynamically update the connection status of the edges. The system utilizes real-time monitoring data transmitted via fiber optic cables and deep-hole piezometers from the steel arch frame. By employing ensemble Kalman filtering to asynchronously assimilate the twin parameters, it ultimately anticipates the extreme values ​​of the seepage field and the breakthrough path for the next tunneling cycle.

2. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 1, characterized in that, The aforementioned construction of a follow-up in-situ sensing network integrated into the initial support and advanced pipe roof, combined with the morphology of the preceding fault zone obtained from tunnel advanced geological prediction, generates an initial digital twin containing a three-dimensional discrete fracture network, specifically including: Fiber optic strain sensors are implanted along the longitudinal main reinforcement direction at the flange and web positions of each steel arch frame. At the same time, the fiber optic strain sensors are bound to the inner wall of the steel pipe of the advanced pipe shed, forming a distributed strain sensing array that moves forward with the tunnel face. Miniature multi-parameter sensors consisting of a pressure sensor and a moisture content sensor are installed in the sidewall of the completed support section behind the working face and in the advance test hole in front of the working face. The sensors in the advance test hole are arranged at equal intervals along the depth direction of the hole. By using the seismic wave reflection method and ground-penetrating radar method in tunnel advanced geological prediction technology, the geometry and spatial distribution of the fault fracture zone and primary fracture group in front of the tunnel face are jointly inverted. The coordinates of the measuring points of the fiber optic strain sensor and the micro multi-parameter sensor, as well as the inverted geometric parameters of the fracture group, are mapped together to the same three-dimensional spatial coordinate system to generate a three-dimensional discrete fracture network containing the initial value of fracture aperture, strike, dip angle and roughness attributes. This three-dimensional discrete fracture network and the surrounding rock medium it contains together constitute the initial digital twin.

3. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 2, characterized in that, The fiber optic strain sensor is demodulated using wavelength scanning optical frequency domain reflection technology, and the sensor implanted on each steel arch forms an independent strain monitoring section.

4. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 1, characterized in that, The establishment of a dynamic fracture aperture evolution model in the twin, controlled by four mechanisms—tunneling unloading, support reaction force, water absorption expansion, and hydraulic scour—specifically includes: The hydraulic opening of each discrete fracture in the three-dimensional discrete fracture network is defined as a dynamic variable that varies with time and spatial coordinates. This dynamic variable is formed by the superposition of the initial hydraulic opening and the opening increment generated by the four driving forces. The first driving force is the excavation unloading opening amount. Based on the three-dimensional spatial distance between the current working face and the location of the discrete fracture, the normal tensile displacement of the fracture surface caused by the redistribution of surrounding rock stress is calculated. This tensile displacement gradually accumulates as the working face advances forward. The second driving force is the closing amount of the support reaction force. Based on the dynamic support force applied to the surrounding rock surface by the steel arch frame and shotcrete, the compressive displacement generated by the support force in the normal direction of the fracture surface is calculated. The direction of the compressive displacement is opposite to that of the tensile displacement. The third driving force is the water absorption and expansion closure amount. Based on the local water content data transmitted back in real time by the micro multi-parameter sensor, the extrusion displacement generated on both sides of the crack by the volume expansion of soft rock after it comes into contact with water is calculated. This extrusion displacement is opposite to the direction of the tensile displacement. The fourth driving force is hydraulic mudification capacity expansion. When the hydraulic gradient of seepage in the fissure exceeds the critical starting value of soft rock particles, the piping erosion module is activated to calculate the mass of soft rock particles peeled off from the fissure wall and to calculate the increase in fissure aperture due to particle loss based on the peeling mass.

5. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 4, characterized in that, In the piping erosion module, the peeling rate of soft rock particles is proportional to the cube of the local seepage velocity, and the capacity expansion calculation stops when the peeling quality reaches the upper limit of the total amount that can be eroded on the fracture wall.

6. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 1, characterized in that, Based on the spatial advancement coordinates of the tunnel face, a graph neural network is used to topologically reshape the fracture network graph and dynamically update the edge connection states. Specifically, this includes: Each discrete fracture in the three-dimensional discrete fracture network is mapped to an edge in the graph data structure, and the intersection points between discrete fractures or the contact surfaces between fractures and the surrounding rock matrix are mapped to nodes in the graph data structure, thereby forming the initial fracture graph topology. Assign an initial feature vector to each node and each edge in the graph data structure. The feature vector of a node includes the spatial coordinates and pressure head of the point, and the feature vector of an edge includes the dynamic hydraulic aperture of the discrete fracture, the current permeability coefficient, and the water content of the soft rock on both sides. The three-dimensional spatial coordinates of the tunnel face after the current tunneling cycle ends are used as external input parameters and passed into a multi-layer graph neural network model. The graph neural network model iteratively aggregates and updates the feature vectors of each edge in the initial fracture diagram topology. When the dynamic hydraulic opening value of a certain edge converges to zero or below the preset closure threshold after being calculated by the model, the graph neural network model marks the state of the edge in the graph as open, which means that the discrete fracture has been completely closed and water flow cannot pass through. When the dynamic hydraulic opening value of a certain edge exceeds the preset widening threshold after model calculation, the graph neural network model marks the state of the edge in the graph as bold, indicating that the discrete fracture has been scoured into a dominant seepage channel. The graph neural network model outputs an updated crack graph topology, where changes in the state of broken or thickened edges trigger a redistribution of pressure head features on adjacent nodes.

7. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 6, characterized in that, The graph neural network model employs a graph attention layer to weighted aggregate the features of adjacent nodes and edges, where the attention coefficient is dynamically determined by the relative distance between the current face advance coordinates and the spatial orientation of the crack.

8. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 1, characterized in that, The monitoring data transmitted in real time via the steel arch frame fiber optic cable and deep hole piezometer is asynchronously assimilated into twin parameters using ensemble Kalman filtering, specifically including: Within each assimilation time window, the surrounding rock deformation pressure values ​​at each measuring point are collected from the distributed optical fiber sensor on the steel arch frame, and the pore water pressure values ​​at each deep hole measuring point are collected from the seepage sensor in the miniature multi-parameter sensor. These two types of measured data are combined to form an observation vector. The three uncertain parameters in the dynamic fracture aperture evolution model—unloading modulus, soft rock expansion coefficient, and piping erosion rate—are taken as state variables to be assimilated, and an initial set that follows its prior distribution is generated for each state variable. The observation vector is input into the ensemble Kalman filter, which calculates the Kalman gain matrix for each ensemble member. The predicted values ​​of the state variables are corrected using the gain matrix to obtain the updated ensemble mean and ensemble covariance. The updated unloading modulus, expansion coefficient, and erosion rate are re-substituted into the dynamic fracture aperture evolution model to replace the original parameters, thus completing the rolling optimization of twin parameters within the current time window.

9. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 5, characterized in that, The assimilation time window of the ensemble Kalman filter is consistent with the duration of one tunnel excavation cycle, and the number of each ensemble member is dynamically adjusted to the square root of the current number of nodes in the fracture network.

10. The three-dimensional numerical simulation method for seepage field evolution in soft rock tunnels based on digital twins according to claim 8, characterized in that, After completing the rolling optimization of the twin parameters within the current time window, the process further includes the following steps for advance simulation of the seepage field in the next tunneling cycle: The virtual excavation advance length of the next tunneling cycle and the final spatial coordinates of the face corresponding to that advance are input as boundary conditions into the dynamic fracture aperture evolution model that has completed parameter assimilation. Within each micro-hour step of the virtual excavation, the four-fold driving force calculation function is sequentially called to update the dynamic hydraulic opening of each discrete fracture in the three-dimensional discrete fracture network, and the graph neural network model is simultaneously called to update the topology of the fracture map. Using the pressure head at each intersection node in the three-dimensional discrete fracture network as the solution variable, the seepage field solution based on the finite volume method is performed on the updated fracture map to obtain the pore pressure cloud map and seepage velocity vector field at the moment when the virtual excavation progress is completed. Search for all connected paths in the seepage velocity vector field that exceed the preset critical flow velocity for water inrush, and mark the set of discrete fractures that the connected path passes through as potential channels for water and mud inrush.

Images