Water conveyance tunnel construction face risk early warning method based on digital twinning

Abstract

The application relates to a water conveying tunnel construction face risk early warning method based on digital twinning, in particular to the field of digital twinning, an agent model is established by BIM initial geological exploration parameters to generate an initial parameter set, observation vectors are constructed by combining real-time monitoring data of microseismic sensors and osmotic pressure gauges, parameter set is updated for multiple iterations by using an iterative set smoother to realize parameter calibration; the calibrated parameter set is directly used as a Monte Carlo sample to drive the agent model to deduce the surrounding rock response, the safety factor is calculated by using the strength reduction method, the collapse probability is statistically calculated by combining the shear failure zone through criterion, and the result is labeled to the BIM twin and a closed loop is formed; the scheme realizes real-time parameter calibration and risk deduction driven by multi-source data, significantly improves the accuracy and timeliness of collapse early warning, reduces the misjudgment caused by parameter mismatch, guarantees construction safety and reduces engineering risk.

Description

Technical Field

[0001] This invention relates to the field of digital twins, and more specifically, to a method for risk early warning at the working face of a water conveyance tunnel based on digital twins. Background Technology

[0002] With the large-scale construction of inter-basin water transfer projects, water conveyance tunnels often traverse areas with complex geological conditions. During construction, the stability of the tunnel face directly affects the safety and progress of the project. Currently, drill-and-blast and tunnel boring machine methods are the main construction methods. Microseismic sensors and piezometers are usually installed behind the excavated section of the tunnel face to monitor the micro-fracture activity of the surrounding rock and the groundwater seepage status in real time. However, the mechanical parameters of the unexcavated rock mass in front of the tunnel face, such as the elastic modulus, internal friction angle, and permeability coefficient, cannot be directly measured because they are not visible. These parameters are strongly affected by internal factors such as the development of joints and significant anisotropy of the surrounding rock, and are also subject to the continuous effects of external factors such as fluctuations in groundwater recharge and blasting vibration disturbances, resulting in extremely high spatial variability. The number of traditional geological exploration boreholes is limited, and the parameters obtained are difficult to fully reflect the true mechanical behavior of the rock mass in front. Against this background, digital twin technology has been introduced into construction risk management, aiming to achieve early warning by constructing a high-fidelity virtual model, integrating real-time monitoring data, and dynamically extrapolating the probability of tunnel face collapse.

[0003] The application of existing digital twin models in water conveyance tunnel construction still has significant shortcomings. Most models adopt an open-loop simulation architecture, establishing finite element or finite difference models based solely on initial geological survey parameters, and keeping these parameters fixed throughout the construction process. These models cannot utilize the massive monitoring data accumulated during the construction phase, including the spatiotemporal distribution characteristics of microseismic events and dynamic changes in seepage pressure curves, for automatic calibration and updates. As tunneling progresses, the mismatch between model parameters and the actual rock mass condition gradually accumulates, leading to a sharp decline in the accuracy of collapse probability estimation. Although some studies have attempted to use response surface methodology for parameter inversion, its computational cost is high, making it difficult to meet the real-time requirements during construction. Furthermore, manual parameter adjustment is not only inefficient but also highly dependent on operational experience, failing to establish a standardized and automated closed-loop process. Therefore, a lightweight, online parameter calibration method is urgently needed. This method should use continuous monitoring data such as microseismic event rate and seepage pressure change rate as observations, and update the mechanical parameter field of the unexcavated rock mass in the digital twin in real time through a data assimilation algorithm. This would ensure that the collapse probability projection model always remains consistent with the actual situation on site. If this problem is not effectively solved, the digital twin will continue to use incorrect parameters, leading to an overestimation or underestimation of the collapse probability, which could result in unnecessary work stoppages or missed disaster reports, causing serious project delays and safety accidents. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a method for risk early warning at the working face of a water conveyance tunnel based on digital twins, thereby solving the problems mentioned in the background art.

[0005] The technical solution of this invention to solve the above-mentioned technical problems is as follows: specifically, it includes the following steps: Step S1: Based on the initial geological parameters of BIM, the unexcavated rock mass in front of the tunnel face within the preset range is divided into multiple parameter estimation sub-regions along the axial and radial directions. Each parameter estimation sub-region is assigned a parameter vector consisting of elastic modulus, internal friction angle and permeability coefficient. At the same time, an axisymmetric simplified mechanical model based on strain softening constitutive model is established as a surrogate model to generate an initial parameter set. The initial parameter set consists of the parameter vectors of all parameter estimation sub-regions. Step S2: In each construction cycle, acquire real-time monitoring data from the microseismic sensors and piezometers deployed in the excavated section behind the tunnel face, and construct microseismic channel observation operators and piezometer observation operators respectively. Map the spatial density distribution of microseismic events to the volume increment of the shear failure zone, and map the piezometer changes to the pore pressure boundary conditions to form the observation vector of the current construction cycle. Among them, the shear failure zone is defined as the set of units in the surrogate model where the damage variable exceeds the critical damage threshold. Step S3: For the first construction cycle, the initial parameter set is used as the parameter set to be updated; for construction cycles other than the first construction cycle, the parameter set calibrated by the previous construction cycle is used as the parameter set to be updated; through an iterative set smoother, combined with the observation vector of the current construction cycle, the microseismic channel observation operator and the seepage pressure channel observation operator, the parameter set to be updated is iteratively updated multiple times to obtain the parameter set calibrated for the current construction cycle. Step S4: The parameter set calibrated for the current construction cycle is directly used as a Monte Carlo sample to drive the proxy model to infer the surrounding rock response for the next construction cycle. The strength reduction method is used to calculate the safety factor, and the collapse criterion is that the safety factor is lower than the preset safety threshold and the shear failure zone extends to the crown. The collapse probability for the next construction cycle is statistically obtained, and the collapse probability is marked in the BIM twin. The parameter set calibrated for the current construction cycle is used as the parameter set to be updated in step S3 of the next construction cycle, forming a closed loop. In a preferred embodiment, in step S1, the initial BIM geological survey parameters consist of advanced geological prediction data of the tunnel face, borehole radar detection data, surface mapping data, and regional geological profile maps. Based on the initial BIM geological survey parameters, surrounding rock classification information and joint development statistics are extracted. The surrounding rock classification information reflects the differences in rock mass quality in different sections, and the joint development statistics include joint density, occurrence, and connectivity indicators.

[0006] In a preferred embodiment, the specific operation of dividing the unexcavated rock mass within a preset range in front of the tunnel face into multiple parameter-evaluated sub-regions along the axial and radial directions is as follows: The area five times the tunnel diameter in front of the tunnel face is divided into multiple axial layers along the tunnel axis. The thickness of each axial layer is equal to the advance length of a single construction cycle. The division is based on the surrounding rock classification information and joint development statistics in the initial geological parameters of BIM. Within each axial layer, the rock mass is divided into multiple radial partitions according to the radial distance from the center of the tunnel face and the magnitude of the joint development index extracted from the joint development statistics. Finally, multiple parameter sub-regions are formed by the product of the number of axial layers and the number of radial partitions. Each parameter sub-region to be estimated is assigned a unique spatial coordinate identifier, which is composed of an axial layer number and a radial partition number. Each parameter sub-region to be estimated corresponds to a parameter vector consisting of the elastic modulus, the internal friction angle, and the permeability coefficient.

[0007] In a preferred embodiment, the specific process of generating the initial parameter set is as follows: First, based on the statistical ranges of elastic modulus, internal friction angle, and permeability coefficient, the lower and upper limits of the values ​​for each parameter vector to be estimated are determined. Then, the Latin hypercube sampling method is used to generate initial samples. The specific execution process of the Latin hypercube sampling method is as follows: For the elastic modulus, internal friction angle and permeability coefficient in the parameter vector to be estimated, the complete range of values ​​for each parameter is divided into several non-overlapping sub-intervals of equal length. For each parameter dimension, a uniform random number between zero and one is generated. The specific representative value is calculated by using the relative position of the uniform random number in the corresponding equal-length subinterval. The sampling order of each dimension is randomly permuted so that each equal-length subinterval is selected exactly once in each dimension, thereby generating a preset number of parameter vectors covering the full parameter space distribution. All parameter vectors to be estimated together constitute the initial parameter set. In the stage of establishing the surrogate model, an axisymmetric simplified mechanical model based on strain softening constitutive model is introduced. The constitutive relation of the surrogate model takes the exponential decay of the internal friction angle with the equivalent plastic strain as the core mechanism. The internal friction angle under the current equivalent plastic strain is determined as follows: first, the difference between the peak internal friction angle and the residual internal friction angle is multiplied by one and the result calculated by raising the product of the negative softening modulus and the equivalent plastic strain to the natural base. Then, the result obtained by multiplying is subtracted from the peak internal friction angle. The surrogate model defines the damage variable as the ratio of the equivalent plastic strain to the critical plastic strain. When the damage variable is greater than one, the element is deemed to have lost its bearing capacity. The state equation of the surrogate model describes the dynamic evolution relationship between the current state vector and the previous state vector and the pore pressure boundary conditions. The state vector includes stress, damage, pore pressure and displacement. The constitutive relation and state equation of the surrogate model are used to determine the damage variable.

[0008] In a preferred embodiment, the specific operation of constructing the microseismic channel observation operator in step S2 is as follows: The raw waveform data in the real-time monitoring data recorded by the microseismic sensor in each construction cycle is automatically picked up and located to obtain the set of spatiotemporal coordinates of microseismic events occurring in the excavated section behind the working face in the current construction cycle. An adaptive threshold screening based on the energy attenuation rate of microseismic events was adopted to retain only microseismic events whose energy was greater than twice the moving median of the energy of microseismic events in the first three construction cycles, thus forming a set of screened microseismic events. The weighted covariance matrix is ​​calculated for the three-dimensional coordinates of all microseismic events in the selected microseismic event set. The construction process of the weighted covariance matrix is ​​as follows: First, the weighted centroid of all monitored microseismic events is determined. The spatial location coordinates of the weighted centroid are obtained by multiplying the occurrence coordinates of each microseismic event by its corresponding weight, summing them up, and then dividing by the sum of all weights. After determining the weighted centroid, each element in the weighted covariance matrix is ​​obtained by calculating the outer product of the difference between the coordinates of each microseismic event and the coordinates of the weighted centroid, multiplying the calculated outer product value by the corresponding weight of the microseismic event, and finally summing the calculated weighted outer product values ​​of all microseismic events and dividing by the sum of all weights. Subsequently, the three eigenvalues ​​of the weighted covariance matrix are extracted and arranged in descending order. The lengths of the three semi-axis are equal to the square roots of the three eigenvalues. An ellipsoid is constructed based on the lengths of the three semi-axis. The volume of the ellipsoid is equal to four-thirds of pi multiplied by the product of the three semi-axis lengths. The volume of the shear failure zone at the end of the previous construction cycle is then subtracted. The difference is the volume increment of the shear failure zone in the current construction cycle. The volume increment of the shear failure zone corresponds to the volume change of the shear failure zone formed by the set of units in the surrogate model whose damage variables exceed the critical damage threshold. The volume increment of the shear failure zone is directly used as the first component of the observation vector.

[0009] In a preferred embodiment, the construction process of the osmotic pressure channel observation operator is as follows: Real-time monitoring data of multiple piezometers deployed in the excavated section behind the working face during the current construction cycle are obtained. Each piezometer has its own axial distance from the working face. A control equation describing the pore pressure attenuation along the axial direction is introduced. The fractional derivative of the control equation with respect to time is equal to the diffusion coefficient multiplied by the second partial derivative of the axial distance. The frequency domain relationship of the governing equations is obtained through Laplace transform. The frequency domain relationship is as follows: the frequency domain component of the real-time monitoring data of each piezometer after Laplace transform is multiplied by the product of the negative axial distance and the first half of the frequency domain component, and then divided by the first half of the diffusion coefficient raised to the natural base. The real-time monitoring data of each piezometer is then subjected to Laplace transform, substituted into the frequency domain relationship, and then the inverse Laplace transform is performed to obtain the equivalent pore pressure at the working face. During the inversion calculation, spatial weights are assigned based on the axial distance of each piezometer from the working face. The inversion results of each piezometer are then weighted and summed according to the assigned spatial weights to obtain the equivalent pore pressure at the working face in the current construction cycle. The equivalent pore pressure at the working face is directly used as the second component of the observation vector. The observation vector consists of two components: the volume increment of the shear failure zone and the equivalent pore pressure at the working face. These components are input into the iterative set smoother in step S3 to update the set of parameters to be updated.

[0010] In a preferred embodiment, in step S3, the specific process of iteratively updating the parameter set to be updated multiple times using an iterative set smoother, combining the observation vector of the current construction cycle, the microseismic channel observation operator, and the seepage pressure channel observation operator, is as follows: First, each parameter vector to be estimated in the parameter set to be updated is input into the surrogate model established in step S1, and the predicted observation vector corresponding to each parameter vector to be estimated is calculated through the microseismic channel observation operator and the pressure channel observation operator constructed in step S2. All the predicted observation vectors constitute the predicted observation set. Then, the mean vector of the parameter set to be updated and the mean vector of the predicted observation set are calculated, and parameter perturbation matrix and observation perturbation matrix are constructed respectively. Each column of the parameter perturbation matrix is ​​the square root of the difference between each estimated parameter vector and the parameter mean vector in the parameter set to be updated, divided by the preset sample number minus one. Each column of the observation perturbation matrix is ​​the square root of the difference between each predicted observation vector and the observation mean vector in the predicted observation set, divided by the preset sample number minus one. Next, a matrix decomposition operation is performed on the observed perturbation matrix, which is decomposed into three components: a left orthogonal matrix, a diagonal matrix, and a right orthogonal matrix. Singular values ​​are extracted from the diagonal matrix, and the cumulative energy ratio of these singular values ​​is calculated. The number of singular values ​​to be retained is determined based on the cumulative energy ratio, so that the proportion of energy contained in the retained singular values ​​to the total energy of all singular values ​​reaches a preset cumulative energy ratio threshold. The columns and rows corresponding to the retained singular values ​​in the left orthogonal matrix, diagonal matrix, and right orthogonal matrix are extracted together to form the truncated left orthogonal matrix, the truncated diagonal matrix, and the truncated right orthogonal matrix. Subsequently, the gain matrix of the iterative ensemble smoother is calculated using the truncated left orthogonal matrix, the truncated diagonal matrix, the truncated right orthogonal matrix, and the parameter perturbation matrix. Finally, each parameter vector in the parameter set to be updated is added to the product of the gain matrix and the difference between the observed vector of the current construction cycle and the predicted observed vector corresponding to the parameter vector to be updated, to obtain the updated parameter vector. All the updated parameter vectors to be updated constitute the parameter set to be updated for the next iteration. All operations in step S3 from calculating the predicted observed vector to obtaining the updated parameter vector are repeated. The parameter set to be updated obtained in the last iteration is the parameter set after calibration for the current construction cycle.

[0011] In a preferred embodiment, in step S3, after each iteration, boundary clipping is performed on each parameter vector to be estimated in the updated parameter set. The specific operation of boundary clipping is as follows: The elastic modulus, internal friction angle, and permeability coefficient in the parameter vector to be estimated are each compared with the lower limit and upper limit determined in step S1. If the value of the elastic modulus is less than the corresponding lower limit, the value of the elastic modulus is forcibly assigned to the lower limit; if the value of the internal friction angle is less than the corresponding lower limit, the value of the internal friction angle is forcibly assigned to the lower limit; if the value of the permeability coefficient is less than the corresponding lower limit, the value of the permeability coefficient is forcibly assigned to the lower limit; if the value of the elastic modulus is greater than the corresponding upper limit, the value of the elastic modulus is forcibly assigned to the upper limit; if the value of the internal friction angle is greater than the corresponding upper limit, the value of the internal friction angle is forcibly assigned to the upper limit; if the value of the permeability coefficient is greater than the corresponding upper limit, the value of the permeability coefficient is forcibly assigned to the upper limit; if the value of the elastic modulus, internal friction angle, or permeability coefficient is between the corresponding lower limit and upper limit, the value of that parameter remains unchanged, and the parameter vector to be estimated after boundary trimming constitutes the parameter set to be updated in the next iteration.

[0012] In a preferred embodiment, the specific process of calculating the safety factor using the strength reduction method in step S4 is as follows: For each parameter vector to be estimated in the parameter set after the current construction cycle calibration, the internal friction angle and elastic modulus in the parameter vector to be estimated are taken as reduction objects; an initial reduction factor is set, with an initial value of 1.0, and the reduction step size is set to 0.01; in each reduction iteration, firstly, according to the coupling strength reduction principle, the internal friction angle in the current parameter vector to be estimated is reduced. The specific processing method is as follows: obtain the tangent value of the current internal friction angle, divide the tangent value by the current reduction factor, and then take the arctangent operation of the result of the division to calculate the reduced internal friction angle; at the same time, the elastic modulus in the current parameter vector to be estimated is reduced synchronously. The specific processing method is as follows: directly divide the current elastic modulus by the reduction factor to calculate the reduced elastic modulus; Substitute the reduced internal friction angle and the reduced elastic modulus into the surrogate model established in step S1, recalculate the surrounding rock response for the next construction cycle, extract the damage variable field output by the surrogate model, and mark the elements whose damage variables exceed the critical damage threshold as shear failure elements. All shear failure elements constitute the current shear failure zone. The following steps determine the continuity status of the shear failure zone: First, all units within a vertical distance of 0.5 times the tunnel diameter near the tunnel arch are designated as the target area. Next, starting from the unit at the tunnel face, a breadth-first search algorithm is used to traverse all units marked as shear failure units to detect whether there exists a continuous spatial path consisting entirely of adjacent shear failure units. This continuous spatial path starts from the tunnel face and extends to the target area. If such a continuous path is detected, the shear failure zone is determined to have reached the arch. At this point, the reduction factor is stopped from increasing, and the current reduction factor value is recorded as the safety factor corresponding to the parameter vector to be estimated. If no continuous continuity path is detected, the reduction factor is increased by one step, and all operations of reducing the internal friction angle and elastic modulus, extrapolating the surrogate model response, and determining the continuity status are repeated. A complete reduction and judgment process is performed on all the parameter vectors to be estimated, and finally the safety factor corresponding to each parameter vector is obtained.

[0013] In a preferred embodiment, the process of statistically obtaining the collapse probability for the next construction cycle specifically involves: First, for each parameter vector to be estimated in the parameter set after calibration in the current construction cycle, failure sample determination is performed based on the safety factor calculated in step S4 and the state of whether the shear failure zone extends to the crown as determined in step S4. If the safety factor corresponding to the parameter vector to be estimated is lower than the preset safety threshold, and the shear failure zone corresponding to the parameter vector to be estimated extends to the crown, then the parameter vector to be estimated is marked as a failure sample. The number of all samples marked as failure samples is counted and recorded as the number of failure samples. The total number of parameter vectors to be estimated in the parameter set after calibration in the current construction cycle is obtained and recorded as the total number of samples. The number of failure samples is divided by the total number of samples to calculate the collapse probability of the next construction cycle.

[0014] The beneficial effects of this invention are as follows: A proxy model is established using initial BIM geological parameters to generate an initial parameter set. An observation vector is constructed by combining real-time monitoring data from microseismic sensors and piezometers. An iterative ensemble smoother is used to iteratively update the parameter set multiple times to achieve parameter calibration. The calibrated parameter set is directly used as a Monte Carlo sample to drive the proxy model to extrapolate the surrounding rock response. The safety factor is calculated using the strength reduction method, and the collapse probability is statistically analyzed using the shear failure zone penetration criterion. Simultaneously, the results are labeled in the BIM twin to form a closed loop. This scheme achieves real-time parameter calibration and risk extrapolation driven by multi-source data, significantly improving the accuracy and timeliness of collapse early warning, reducing misjudgments caused by parameter mismatch, ensuring construction safety, and lowering engineering risks. Attached Figure Description

[0015] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

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

[0017] In the description of this application, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more features. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0018] In the description of this application, the term "for example" is used to mean "used as an example, illustration, or description." Any embodiment described as "for example" in this application is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is provided to enable any person skilled in the art to make and use the invention. Details are set forth in the following description for purposes of explanation. It should be understood that those skilled in the art will recognize that the invention can be made without using these specific details. In other instances, well-known structures and processes will not be described in detail to avoid obscuring the description of the invention with unnecessary detail. Therefore, the invention is not intended to be limited to the embodiments shown, but is consistent with the broadest scope of the principles and features disclosed in this application.

[0019] Example 1: This example provides the following... Figure 1 The method for risk early warning at the working face of a water conveyance tunnel based on digital twins, as shown, specifically includes the following steps: Step S1: Based on the initial geological parameters of BIM, the unexcavated rock mass in front of the tunnel face within the preset range is divided into multiple parameter estimation sub-regions along the axial and radial directions. Each parameter estimation sub-region is assigned a parameter vector consisting of elastic modulus, internal friction angle and permeability coefficient. At the same time, an axisymmetric simplified mechanical model based on strain softening constitutive model is established as a surrogate model to generate an initial parameter set. The initial parameter set consists of the parameter vectors of all parameter estimation sub-regions. Step S2: In each construction cycle, acquire real-time monitoring data from the microseismic sensors and piezometers deployed in the excavated section behind the tunnel face, and construct microseismic channel observation operators and piezometer observation operators respectively. Map the spatial density distribution of microseismic events to the volume increment of the shear failure zone, and map the piezometer changes to the pore pressure boundary conditions to form the observation vector of the current construction cycle. Among them, the shear failure zone is defined as the set of units in the surrogate model where the damage variable exceeds the critical damage threshold. Step S3: For the first construction cycle, the initial parameter set is used as the parameter set to be updated; for construction cycles other than the first construction cycle, the parameter set calibrated by the previous construction cycle is used as the parameter set to be updated; through an iterative set smoother, combined with the observation vector of the current construction cycle, the microseismic channel observation operator and the seepage pressure channel observation operator, the parameter set to be updated is iteratively updated multiple times to obtain the parameter set calibrated for the current construction cycle. Step S4: The parameter set calibrated for the current construction cycle is directly used as a Monte Carlo sample to drive the proxy model to infer the surrounding rock response for the next construction cycle. The strength reduction method is used to calculate the safety factor, and the collapse criterion is that the safety factor is lower than the preset safety threshold and the shear failure zone extends to the crown. The collapse probability for the next construction cycle is statistically obtained, and the collapse probability is marked in the BIM twin. The parameter set calibrated for the current construction cycle is used as the parameter set to be updated in step S3 of the next construction cycle, forming a closed loop.

[0020] In this embodiment, it is specifically necessary to explain that in step S1, the initial BIM geological survey parameters consist of face-to-face advanced geological prediction data, borehole radar detection data, surface mapping data, and regional geological profile maps. Based on the initial BIM geological survey parameters, surrounding rock classification information and joint development statistics are extracted. The surrounding rock classification information reflects the differences in rock mass quality in different sections, and the joint development statistics include joint density, occurrence, and connectivity indicators. The face-to-face advanced geological prediction data is obtained through ground-penetrating radar or seismic wave reflection methods, with a detection depth of 30 to 50 meters in front of the face. The face-to-face advanced geological prediction data records areas of rock mass wave impedance anomalies in the form of two-dimensional profiles. The borehole radar detection data is collected using single-hole or multi-hole radar systems with a borehole spacing of 5 to 10 meters. The borehole radar detection data includes electromagnetic wave attenuation coefficients and dielectric constants, used to identify areas with dense fractures. The surface mapping data, obtained through UAV photogrammetry or laser scanning with a resolution better than 0.1 meters, is used to establish the distribution of surface topography and structural traces. The regional geological profile map is derived from the geological mapping in the early exploration stage, with a scale of 1:5000 to 1:10000, providing information such as stratigraphic lithology, fault location, and fold axis. The surrounding rock classification information is based on the surrounding rock classification standard in the "Code for Geological Investigation of Water Conservancy and Hydropower Engineering" GB50487, classifying the rock mass into categories I to V. The statistical ranges of elastic modulus, internal friction angle, and permeability coefficient corresponding to each category are pre-stored in the BIM database. The joint density in the joint development statistics is measured by the number of joints per cubic meter of rock mass, the attitude is expressed by dip and dip angle, and the connectivity is defined as the ratio of joint trace length to rock bridge length. All three indicators are statistically derived from borehole core logging and surface outcrop measurement. The specific operation of dividing the unexcavated rock mass within a preset range in front of the tunnel face into multiple parameter-evaluated sub-regions along the axial and radial directions is as follows: The area five times the tunnel diameter in front of the tunnel face is divided into multiple axial layers along the tunnel axis. The thickness of each axial layer is equal to the advance length of a single construction cycle, which is two to three meters. The division is based on the surrounding rock classification information and joint development statistics in the initial geological parameters of BIM. Within each axial layer, the rock mass is divided into multiple radial zones according to the radial distance from the center of the tunnel face and the magnitude of the joint development index extracted from the joint development statistics. The larger the joint development index, the more radial zones are divided. Finally, multiple parameter-to-be-estimated sub-areas are formed by the product of the number of axial layers and the number of radial zones. Each parameter sub-region to be estimated is assigned a unique spatial coordinate identifier, which is composed of an axial layer number and a radial partition number. Each parameter sub-region to be estimated also corresponds to a parameter vector to be estimated consisting of the elastic modulus, internal friction angle, and permeability coefficient, so as to ensure that the parameter vector to be estimated for each parameter sub-region can be uniquely identified and called. The specific calculation method for the five-times-diameter range is as follows: If the tunnel excavation diameter is a specific value, the calculation range in front of the tunnel face is five times that specific value. For example, when the excavation diameter is 6 meters, the calculation range is 30 meters. The number of axial layers is determined by dividing the calculation range by the advance length. When the advance length is 2.5 meters, the number of layers is 12. The rule for determining the number of radial zones is as follows: Divide the distance from the center of the tunnel face to the tunnel wall into three equal parts as the base number. Then, densify the zone according to the joint development index: When the joint development index is less than 5 joints per cubic meter, the number of radial zones is 3; when the joint development index is between 5 and 15 joints per cubic meter, the number of radial zones is 4; when the joint development index is greater than 15 joints per cubic meter, the number of radial zones is 5. The joint development index is determined by the joint density. The joint development index is calculated by weighting the degree of joint development, the dispersion of attitude, and the connectivity rate. The calculation formula is: the joint development index equals the first weight coefficient multiplied by the joint density, plus the second weight coefficient multiplied by one minus the concentration of attitude, plus the third weight coefficient multiplied by the connectivity rate. The first weight coefficient is 0.4, the second weight coefficient is 0.3, and the third weight coefficient is 0.3. The concentration of attitude is a dimensionless value between 0 and 1. The spatial coordinate identifier uses a string format, for example, "Lq03_Rq02" represents the 2nd radial partition of the 3rd axial layer. The parameter vector of each parameter to be estimated sub-region is stored as a three-dimensional array in the BIM model. The indexing method is that the first dimension is the axial layer number, and the second dimension is the radial partition number. The axial layer number ranges from 1 to the layer number, and the radial partition number ranges from 1 to the radial partition number. The specific process of generating the initial parameter set is as follows: First, based on the statistical ranges of elastic modulus, internal friction angle, and permeability coefficient, the lower and upper limits of the values ​​for each parameter vector to be estimated are determined. Then, the Latin hypercube sampling method is used to generate initial samples. The specific execution process of the Latin hypercube sampling method is as follows: For the elastic modulus, internal friction angle and permeability coefficient in the parameter vector to be estimated, the complete value range of each parameter is divided into several non-overlapping equal-length sub-intervals, and the total number of equal-length sub-intervals is equal to the preset total number of samples. For each parameter dimension, a uniformly random number between zero and one is generated. The specific representative value is calculated using the relative position of the uniformly random number within the corresponding equal-length subinterval. The sampling order of each dimension is randomly permuted so that each equal-length subinterval is selected exactly once in each dimension, thereby generating a preset number of parameter vectors covering the full parameter space distribution. The preset number of samples ranges from eighty to one hundred. All parameter vectors to be estimated together constitute the initial parameter set. In the stage of establishing the surrogate model, an axisymmetric simplified mechanical model based on strain softening constitutive model is introduced. The constitutive relation of the surrogate model takes the exponential decay of the internal friction angle with the equivalent plastic strain as the core mechanism. The internal friction angle under the current equivalent plastic strain is determined as follows: first, the difference between the peak internal friction angle and the residual internal friction angle is multiplied by one and the result calculated by raising the product of the negative softening modulus and the equivalent plastic strain to the natural base. Then, the result obtained by multiplying is subtracted from the peak internal friction angle. The surrogate model defines the damage variable as the ratio of the equivalent plastic strain to the critical plastic strain. When the damage variable is greater than one, the element is deemed to have lost its bearing capacity. The state equation of the surrogate model describes the dynamic evolution relationship between the current state vector and the previous state vector and the pore pressure boundary conditions. The state vector includes stress, damage, pore pressure and displacement. Through the constitutive relation and state equation of the surrogate model, the initial parameter set has the ability to drive the surrogate model to deduce the surrounding rock response, providing a uniform and sufficient prior distribution for the iterative update in step S3. The statistical range of the elastic modulus is taken from the corresponding value of the surrounding rock category in the BIM database. For example, the elastic modulus range of Class II surrounding rock is 10 GPa to 20 GPa, and that of Class III surrounding rock is 5 GPa to 10 GPa. The statistical range of the internal friction angle is determined based on the shear strength test results of the rock mass structural surface, with a typical range of 25 degrees to 45 degrees. The statistical range of the permeability coefficient is converted from the Lurong value of the water pressure test, ranging from 1 x 10⁻⁸ m / s to 1 x 10⁻⁶ m / s. In the Latin hypercube sampling, the preset sample size is 90. The length of each equal-length subinterval is the upper limit of the parameter minus the lower limit of the parameter, divided by the preset sample size. The uniform random number is generated by a pseudo-random number generator in the interval of zero to one. The calculation method for the representative value is: the lower limit of the parameter plus the random permutation index minus the uniform random number, divided by the preset sample size, and then multiplied by the difference between the upper limit and the lower limit of the parameter. The surrogate model is an axisymmetric plane strain finite element model. The element type is a four-node quadrilateral element, with radial mesh refinement and a minimum element size of 0.1 meters. The peak internal friction angle in the strain softening constitutive model is taken from the internal friction angle value in the parameters to be estimated, and the residual internal friction angle is taken as 0.6 to 0.8 times the peak internal friction angle. The softening modulus is taken as 0.5 to 2.0. The critical plastic strain is taken as 0.02 to 0.05, specifically set according to the surrounding rock type: 0.05 for Class II surrounding rock, 0.03 for Class III surrounding rock, and 0.02 for Class IV surrounding rock. When the damage variable exceeds 1.0, the elastic modulus of the element is reduced to one percent of the original value. The state equation is solved by explicit time integration with a time step of 0.001 seconds. The total duration of each construction cycle simulation corresponds to the actual tunneling time of 4 to 8 hours. The actual time is mapped to the number of calculation time steps through a scaling factor. The pore pressure boundary condition is output by the seepage channel observation operator and applied to the tunnel face and tunnel wall surface in the form of nodal loads.

[0021] In this embodiment, the specific operation of constructing the microseismic channel observation operator in step S2 is as follows: The raw waveform data in the real-time monitoring data recorded by the microseismic sensor in each construction cycle is automatically picked up and located to obtain the set of spatiotemporal coordinates of microseismic events occurring in the excavated section behind the working face in the current construction cycle. An adaptive threshold screening based on the energy attenuation rate of microseismic events was adopted to retain only microseismic events whose energy was greater than twice the moving median of the energy of microseismic events in the first three construction cycles, thus forming a set of screened microseismic events. The weighted covariance matrix is ​​calculated for the three-dimensional coordinates of all microseismic events in the selected microseismic event set. The construction process of the weighted covariance matrix is ​​as follows: First, the weighted centroid of all monitored microseismic events is determined. The spatial position coordinates of the weighted centroid are obtained by multiplying the occurrence coordinates of each microseismic event by its corresponding weight, summing them, and then dividing by the sum of all weights. The weight of each microseismic event is the natural logarithm of the energy of the microseismic event plus one. After determining the weighted centroid, each element in the weighted covariance matrix is ​​obtained by calculating the outer product of the difference between the coordinates of each microseismic event and the coordinates of the weighted centroid, multiplying the calculated outer product value by the corresponding weight of the microseismic event, and finally summing the calculated weighted outer product values ​​of all microseismic events and dividing by the sum of all weights. Subsequently, the three eigenvalues ​​of the weighted covariance matrix are extracted and arranged in descending order. The lengths of the three semi-axis are equal to the square roots of the three eigenvalues. An ellipsoid is constructed based on the lengths of the three semi-axis. The volume of the ellipsoid is equal to four-thirds of the volume of pi multiplied by the product of the three semi-axis lengths. This volume is then subtracted from the volume of the shear failure zone at the end of the previous construction cycle. The difference is the volume increment of the shear failure zone in the current construction cycle. The volume increment of the shear failure zone corresponds to the volume change of the shear failure zone formed by the set of units in the surrogate model whose damage variables exceed the critical damage threshold. The volume increment of the shear failure zone is directly used as the first component of the observation vector to characterize the extent of the expansion of the surrounding rock shear failure zone caused by microseismic activity. The energy attenuation rate of microseismic events is calculated as follows: for each microseismic event, the square of the peak amplitude of its waveform envelope is taken as the event energy. The energy attenuation rate is defined as the ratio of the current event energy to the previous event energy. This ratio is used to help determine whether the event is rock mass rupture rather than mechanical noise. The calculation window for the moving median is fixed at all event energy values ​​of the first three construction cycles, and the median is taken as the middle value after sorting (if the number of events is even, the average of the two middle values ​​is taken). The threshold for adaptive threshold filtering is set to twice the moving median of event energy of the first three construction cycles, and this threshold is updated once after each construction cycle. The eigenvalues ​​of the weighted covariance matrix were extracted using the Jacobi iteration method, with an iteration precision of 1 x 10^-6 and a maximum number of iterations of 50. Pi was taken as 3.1415926535. The lengths of the three semi-axis in the ellipsoid volume calculation formula corresponded to the square roots of the three eigenvalues. This square root calculation was implemented using the Newton iteration method, with the initial guess value being half of the eigenvalue, and convergence occurring after three iterations. The volume of the shear failure zone at the end of the previous construction cycle was stored in a temporary variable in the BIM database, overwriting the old value after each update. The unit for the volume increment of the shear failure zone was cubic meters, with three decimal places. The construction process of the osmotic pressure channel observation operator is as follows: Real-time monitoring data of multiple piezometers deployed in the excavated section behind the working face during the current construction cycle are obtained. Each piezometer has its own axial distance from the working face. A control equation describing the pore pressure attenuation along the axial direction is introduced. The fractional derivative of the control equation with respect to time is equal to the diffusion coefficient multiplied by the second partial derivative of the axial distance. The frequency domain relationship of the governing equations is obtained through Laplace transform. The frequency domain relationship is as follows: the frequency domain component of the real-time monitoring data of each piezometer after Laplace transform is multiplied by the product of the negative axial distance and the first half of the frequency domain component, and then divided by the first half of the diffusion coefficient raised to the natural base. The real-time monitoring data of each piezometer is then subjected to Laplace transform, substituted into the frequency domain relationship, and then the inverse Laplace transform is performed to obtain the equivalent pore pressure at the working face. During the inversion calculation, spatial weights are assigned based on the axial distance of each piezometer from the working face. The greater the axial distance, the smaller the spatial weight, and the closer the axial distance, the larger the spatial weight. The inversion results of each piezometer are weighted and summed according to the assigned spatial weights to obtain the equivalent pore pressure at the working face in the current construction cycle. The equivalent pore pressure at the working face is directly used as the second component of the observation vector to characterize the impact of groundwater seepage on the stability of the working face. The observation vector consists of two components: the volume increment of the shear failure zone and the equivalent pore pressure at the working face. These components are input into the iterative ensemble smoother in step S3 to update the set of parameters to be updated. In the fractional derivative diffusion model, the fractional order is fixed at 0.85, determined by curve fitting of field pumping test data, with the fitting error controlled within 5%. The diffusion coefficient Dα ranges from 1 x 10⁻⁴ square meters per second to 1 x 10⁻² square meters per second, with specific values ​​determined by previous hydrogeological tests. If no measured value is available, the default value of 5 x 10⁻³ square meters per second is used. The Laplace transform employs a numerical Laplace transform algorithm with 256 sampling points, a frequency range from 1 x 10⁻⁵ Hz to 1 Hz, and a logarithmically uniform distribution. The inverse Laplace transform uses the Talbot algorithm, which is implemented through deformed coaxial integrals, with the accuracy level set to high precision mode, and the relative error controlled within 1 x 10⁻⁶. Within the power; the calculation of the natural base power in the frequency domain relationship adopts the Taylor series expansion of the exponential function, with 20 expansion terms to ensure accuracy; the specific formula for spatial weight allocation is: the spatial weight γ_i of the i-th piezometer is equal to the reciprocal of the axial distance of the piezometer from the working face divided by the sum of the reciprocals of the axial distances of all piezometers from the working face, i.e., γ_i=(1 / z_i) / Σ(1 / z_j), where z_i is the axial distance of the i-th piezometer in meters; the unit of the equivalent pore pressure at the working face is Pascal, rounded to the nearest integer; the two components of the observation vector are normalized before being fed into the iterative ensemble smoother, with the normalization coefficients being the volume increment of the shear failure zone in the previous construction cycle and the equivalent pore pressure at the working face, respectively, to prevent assimilation imbalance caused by differences in magnitude.

[0022] In this embodiment, it is specifically necessary to explain the process in step S3, where the iterative set smoother, combined with the observation vector of the current construction cycle, the microseismic channel observation operator, and the seepage pressure channel observation operator, performs multiple iterative updates to the parameter set to be updated as follows: First, each parameter vector to be estimated in the parameter set to be updated is input into the surrogate model established in step S1, and the predicted observation vector corresponding to each parameter vector to be estimated is calculated through the microseismic channel observation operator and the pressure channel observation operator constructed in step S2. All the predicted observation vectors constitute the predicted observation set. Then, the mean vector of the parameter set to be updated and the mean vector of the predicted observation set are calculated, and parameter perturbation matrix and observation perturbation matrix are constructed respectively. Each column of the parameter perturbation matrix is ​​the square root of the difference between each estimated parameter vector and the parameter mean vector in the parameter set to be updated, divided by the preset sample number minus one. Each column of the observation perturbation matrix is ​​the square root of the difference between each predicted observation vector and the observation mean vector in the predicted observation set, divided by the preset sample number minus one. Next, a matrix decomposition operation is performed on the observed perturbation matrix, which is decomposed into three components: a left orthogonal matrix, a diagonal matrix, and a right orthogonal matrix. Singular values ​​are extracted from the diagonal matrix, and the cumulative energy ratio of these singular values ​​is calculated. The number of singular values ​​to be retained is determined based on the cumulative energy ratio, so that the proportion of energy contained in the retained singular values ​​to the total energy of all singular values ​​reaches a preset cumulative energy ratio threshold. The cumulative energy ratio threshold is set to 0.95. The columns and rows corresponding to the retained singular values ​​in the left orthogonal matrix, diagonal matrix, and right orthogonal matrix are extracted together to form the truncated left orthogonal matrix, the truncated diagonal matrix, and the truncated right orthogonal matrix. Subsequently, the gain matrix of the iterative ensemble smoother is calculated using the truncated left orthogonal matrix, the truncated diagonal matrix, the truncated right orthogonal matrix, and the parameter perturbation matrix. The specific calculation process of the gain matrix is ​​as follows: First, multiply the truncated right orthogonal matrix with the inverse of the truncated diagonal matrix. Then, multiply the result of the multiplication with the transpose of the truncated left orthogonal matrix. Finally, multiply the parameter perturbation matrix with the result of the multiplication. The final product is the gain matrix. Finally, add the gain matrix to each parameter vector in the parameter set to be updated, multiplied by the difference between the observation vector of the current construction cycle and the predicted observation vector corresponding to the parameter vector to be updated, to obtain the updated parameter vector. All the updated parameter vectors to be updated constitute the parameter set to be updated for the next iteration. Repeat all the operations in step S3 from calculating the predicted observation vector to obtaining the updated parameter vector. The number of iterations is three to five. The parameter set to be updated obtained in the last iteration is the parameter set calibrated for the current construction cycle. The preset sample size is 90, therefore the square root of the preset sample size minus one equals the square root of 89, approximately 9.43. During the construction of the parameter perturbation matrix and the observation perturbation matrix, if the preset sample size minus one becomes zero (i.e., the preset sample size equals one), no iterative update is performed; the current set of parameters to be updated is directly used as the calibrated parameter set. However, in practical applications, the preset sample size is always greater than or equal to 80, so this situation will not occur. The matrix decomposition operation uses the singular value decomposition algorithm. Specifically, the observation perturbation matrix is ​​transformed into a bidiagonal matrix through Haushold transformation, and then the singular values ​​are extracted using the implicit symmetric QR iteration method. The iteration precision is set to 1 x 10^-10. The cumulative energy percentage of the singular values ​​in the diagonal matrix is ​​calculated using... The formula is as follows: Singular values ​​are arranged in descending order, and the square of each singular value is accumulated sequentially until the sum reaches 95% of the sum of the squares of all singular values. The number of singular values ​​retained at this point is the truncation number. If there are singular values ​​less than 1 x 10 - 8 in the truncated diagonal matrix, they are replaced with 1 x 10 - 8 before inverting to prevent numerical overflow. The result of the gain matrix is ​​a matrix whose dimension is the total dimension of the parameter vector to be estimated multiplied by the dimension of the observation vector. This matrix maps the residual in the observation space to the correction amount in the parameter space. During repeated iterations, if the change in the Euclidean norm of the mean vector of the parameter set between two consecutive iterations is less than 1 x 10 - 4, the iteration can be terminated early to reduce computational overhead. In step S3, after each iteration, boundary clipping is performed on each parameter vector to be estimated in the updated parameter set. The specific operation of boundary clipping is as follows: The elastic modulus, internal friction angle, and permeability coefficient in the parameter vector to be estimated are each compared with the lower and upper limits determined in step S1. If the value of the elastic modulus is less than the corresponding lower limit, the value of the elastic modulus is forcibly assigned to the lower limit; if the value of the internal friction angle is less than the corresponding lower limit, the value of the internal friction angle is forcibly assigned to the lower limit; if the value of the permeability coefficient is less than the corresponding lower limit, the value of the permeability coefficient is forcibly assigned to the lower limit; if the value of the elastic modulus is greater than the corresponding upper limit, the value of the elastic modulus is forcibly assigned to the upper limit; if the internal friction angle is greater than the corresponding upper limit, the value of the elastic modulus is forcibly assigned to the upper limit. If the value of the friction angle is greater than the corresponding upper limit, the value of the internal friction angle is forcibly assigned to the upper limit; if the value of the permeability coefficient is greater than the corresponding upper limit, the value of the permeability coefficient is forcibly assigned to the upper limit; if the values ​​of the elastic modulus, internal friction angle, or permeability coefficient are between the corresponding lower limit and upper limit, the value of the parameter remains unchanged, and the parameter vector to be estimated after boundary trimming constitutes the set of parameters to be updated in the next iteration, so as to ensure that all parameter vectors to be estimated are always within the physical feasible domain, and avoid the abnormal calculation of the surrogate model caused by the parameter values ​​deviating from the reasonable range due to iterative updates. The specific values ​​for the lower and upper limits are pre-set according to the surrounding rock category: For Class II surrounding rock, the lower limit for the elastic modulus is 10 GPa and the upper limit is 20 GPa; the lower limit for the internal friction angle is 30 degrees and the upper limit is 45 degrees; the lower limit for the permeability coefficient is 1 x 10^-8 m / s and the upper limit is 1 x 10^-6 m / s; for Class III surrounding rock, the lower limit for the elastic modulus is 5 GPa and the upper limit is 10 GPa; the lower limit for the internal friction angle is 25 degrees and the upper limit is 35 degrees; the lower limit for the permeability coefficient is 1 x 10^-7 m / s and the upper limit is 1 x 10^-5 m / s; for Class IV surrounding rock, the lower limit for the elastic modulus is 2 GPa and the upper limit is 5 GPa; the lower limit for the internal friction angle is 20 degrees. The upper limit is 30 degrees, and the lower limit of the permeability coefficient is 1 x 10^-6 meters per second, and the upper limit is 1 x 10^-4 meters per second. The boundary pruning operation is performed in the last step of each iteration, that is, after the updated parameter vector to be estimated is calculated, but before it is used as the parameter set to be updated in the next iteration. If the pruned parameter value is forcibly assigned to the boundary value, it will not participate in the subsequent calculation in this iteration, but it is still allowed to move into the feasible region in the next iteration. If all parameter vectors to be estimated are pruned to the boundary value and there is no change for three consecutive iterations, it is determined that the parameter set has converged to the boundary, the iteration stops, and the current parameter set is output as the calibrated parameter set.

[0023] In this embodiment, the specific process of calculating the safety factor using the strength reduction method in step S4 is as follows: For each parameter vector to be estimated in the parameter set after the current construction cycle calibration, the internal friction angle and elastic modulus in the parameter vector to be estimated are taken as reduction objects; an initial reduction factor is set, with an initial value of 1.0, and the reduction step size is set to 0.01; in each reduction iteration, firstly, according to the coupling strength reduction principle, the internal friction angle in the current parameter vector to be estimated is reduced. The specific processing method is as follows: obtain the tangent value of the current internal friction angle, divide the tangent value by the current reduction factor, and then take the arctangent operation of the result of the division to calculate the reduced internal friction angle; at the same time, the elastic modulus in the current parameter vector to be estimated is reduced synchronously. The specific processing method is as follows: directly divide the current elastic modulus by the reduction factor to calculate the reduced elastic modulus; Substitute the reduced internal friction angle and the reduced elastic modulus into the surrogate model established in step S1, recalculate the surrounding rock response for the next construction cycle, extract the damage variable field output by the surrogate model, and mark the elements whose damage variables exceed the critical damage threshold as shear failure elements. All shear failure elements constitute the current shear failure zone. The following steps determine the continuity status of the shear failure zone: First, all units within a vertical distance of 0.5 times the tunnel diameter near the tunnel arch are designated as the target area. Next, starting from the unit at the tunnel face, a breadth-first search algorithm is used to traverse all units marked as shear failure units to detect whether there exists a continuous spatial path consisting entirely of adjacent shear failure units. This continuous spatial path starts from the tunnel face and extends to the target area. If such a continuous path is detected, the shear failure zone is determined to have reached the arch. At this point, the reduction factor is stopped from increasing, and the current reduction factor value is recorded as the safety factor corresponding to the parameter vector to be estimated. If no continuous continuity path is detected, the reduction factor is increased by one step, and all operations of reducing the internal friction angle and elastic modulus, extrapolating the surrogate model response, and determining the continuity status are repeated. A complete reduction and judgment process is performed on all the parameter vectors to be estimated, and finally the safety factor corresponding to each parameter vector is obtained. The critical damage threshold is set to 0.8, which is calibrated based on the post-peak strain softening curve of indoor triaxial compression tests on rocks. When the damage variable exceeds 0.8, the shear strength of the element has dropped to less than 20% of the residual strength, and it can be considered as having completely lost its bearing capacity. The reduction step size is set to 0.01 to ensure the accuracy of the safety factor calculation and to control the computational load of a single reduction iteration to no more than 100 runs of the proxy model, avoiding excessive computation time due to an excessively small step size. The physical basis of the coupled strength reduction principle is that during the instability of the surrounding rock, the decrease in the internal friction angle reflects the degradation of shear strength, and the decrease in the elastic modulus reflects the deterioration of the deformation modulus. Simultaneous reduction of both can more realistically simulate the progressive failure process. The model's inputs include the reduced internal friction angle, the reduced elastic modulus, the dynamic load of the current construction cycle, and the pore pressure boundary conditions for the next construction cycle. The outputs are the damage variable field, stress field, and displacement field at the end of the next construction cycle. The surrogate model uses an axisymmetric plane strain finite element model with a total of two to three thousand elements, radially densified, and a minimum element size of 0.1 meters. The breadth-first search algorithm has a queue capacity of five thousand nodes and a maximum search depth of one hundred layers. If no path is found after exceeding the maximum depth, it is considered not to be a breakthrough. The selection of 0.5 times the tunnel diameter is based on the fact that this range covers the key load-bearing area of ​​the arch. If the shear failure zone extends to this area, it indicates that the arch has lost its support and is very prone to collapse. The specific process for statistically determining the collapse probability of the next construction cycle is as follows: First, for each parameter vector to be estimated in the parameter set after calibration in the current construction cycle, failure sample determination is performed based on the safety factor calculated in step S4 and the state of whether the shear failure zone extends to the crown as determined in step S4. If the safety factor corresponding to the parameter vector to be estimated is lower than the preset safety threshold and the shear failure zone corresponding to the parameter vector to be estimated extends to the crown, then the parameter vector to be estimated is marked as a failure sample. The number of all samples marked as failure samples is counted and recorded as the number of failure samples. The total number of parameter vectors to be estimated in the parameter set after calibration in the current construction cycle is obtained and recorded as the total number of samples. The number of failure samples is divided by the total number of samples to calculate the collapse probability of the next construction cycle. The preset safety threshold is 1.05. This value is determined based on the minimum requirement for the safety factor of permanent slopes in the "Code for Design of Slopes in Water Conservancy and Hydropower Projects" and takes into account the allowable reduction range of temporary construction conditions. A safety factor below 1.05 indicates that the surrounding rock is in a critical instability state. The result of dividing the number of failure samples by the total number of samples is a dimensionless decimal between zero and one. For example, 0.35 means that 35% of the estimated parameter vectors correspond to conditions where collapse will occur. If the total number of samples is zero (theoretically, this will not happen because the preset sample number is always greater than or equal to 80), then the probability of collapse is high. The probability is recorded as zero; the calculated collapse probability is marked as a percentage in the BIM twin. The specific marking method is as follows: on the surface of the 3D model corresponding to each parameter sub-area in front of the tunnel face, the collapse probability value corresponding to the sub-area is displayed with a color gradient bar. Areas with a probability higher than 0.5 are marked in red, areas between 0.3 and 0.5 are marked in yellow, and areas lower than 0.3 are marked in green. At the same time, the parameter set calibrated in the current construction cycle is stored in a specified field of the BIM database as the parameter set to be updated in the next construction cycle step S3, overwriting the old value stored in the previous cycle.

[0024] It should be noted that the descriptions of each embodiment in the above embodiments have different focuses. For parts that are not described in detail in a certain embodiment, please refer to the relevant descriptions in other embodiments.

[0025] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0026] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0027] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0028] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0029] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0030] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method for risk early warning at the working face of a water conveyance tunnel based on digital twins, characterized in that, Specifically, the steps include the following: Step S1: Based on the initial geological parameters of BIM, the unexcavated rock mass in front of the tunnel face within the preset range is divided into multiple parameter estimation sub-regions along the axial and radial directions. Each parameter estimation sub-region is assigned a parameter vector consisting of elastic modulus, internal friction angle and permeability coefficient. At the same time, an axisymmetric simplified mechanical model based on strain softening constitutive model is established as a surrogate model to generate an initial parameter set. The initial parameter set consists of the parameter vectors of all parameter estimation sub-regions. Step S2: In each construction cycle, acquire real-time monitoring data from the microseismic sensors and piezometers deployed in the excavated section behind the tunnel face, and construct microseismic channel observation operators and piezometer observation operators respectively. Map the spatial density distribution of microseismic events to the volume increment of the shear failure zone, and map the piezometer changes to the pore pressure boundary conditions to form the observation vector of the current construction cycle. Among them, the shear failure zone is defined as the set of units in the surrogate model where the damage variable exceeds the critical damage threshold. Step S3: For the first construction cycle, the initial parameter set is used as the parameter set to be updated; for construction cycles other than the first construction cycle, the parameter set calibrated by the previous construction cycle is used as the parameter set to be updated; through an iterative set smoother, combined with the observation vector of the current construction cycle, the microseismic channel observation operator and the seepage pressure channel observation operator, the parameter set to be updated is iteratively updated multiple times to obtain the parameter set calibrated for the current construction cycle. Step S4: The parameter set calibrated for the current construction cycle is directly used as a Monte Carlo sample to drive the proxy model to infer the surrounding rock response for the next construction cycle. The strength reduction method is used to calculate the safety factor, and the collapse criterion is that the safety factor is lower than the preset safety threshold and the shear failure zone extends to the crown. The collapse probability for the next construction cycle is statistically obtained, and the collapse probability is marked in the BIM twin. The parameter set calibrated for the current construction cycle is used as the parameter set to be updated in step S3 of the next construction cycle, forming a closed loop.

2. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 1, characterized in that: In step S1, the initial BIM geological survey parameters consist of advanced geological prediction data of the tunnel face, borehole radar detection data, surface mapping data, and regional geological profile maps. Based on the initial BIM geological survey parameters, surrounding rock classification information and joint development statistics are extracted. The surrounding rock classification information reflects the differences in rock mass quality in different sections, and the joint development statistics include joint density, occurrence, and connectivity indicators.

3. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 2, characterized in that: The specific operation of dividing the unexcavated rock mass within a preset range in front of the tunnel face into multiple parameter-evaluated sub-regions along the axial and radial directions is as follows: The area five times the tunnel diameter in front of the tunnel face is divided into multiple axial layers along the tunnel axis. The thickness of each axial layer is equal to the advance length of a single construction cycle. The division is based on the surrounding rock classification information and joint development statistics in the initial geological parameters of BIM. Within each axial layer, the rock mass is divided into multiple radial partitions according to the radial distance from the center of the tunnel face and the magnitude of the joint development index extracted from the joint development statistics. Finally, multiple parameter sub-regions are formed by the product of the number of axial layers and the number of radial partitions. Each parameter sub-region to be estimated is assigned a unique spatial coordinate identifier, which is composed of an axial layer number and a radial partition number. Each parameter sub-region to be estimated corresponds to a parameter vector consisting of the elastic modulus, the internal friction angle, and the permeability coefficient.

4. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 3, characterized in that: The specific process for generating the initial parameter set is as follows: First, based on the statistical range of elastic modulus, internal friction angle and permeability coefficient, determine the lower limit and upper limit of the values ​​of elastic modulus, internal friction angle and permeability coefficient in each parameter vector to be estimated; Then, the initial samples are generated using the Latin hypercube sampling method. The specific execution process of the Latin hypercube sampling method is as follows: For the elastic modulus, internal friction angle and permeability coefficient in the parameter vector to be estimated, the complete range of values ​​for each parameter is divided into several non-overlapping sub-intervals of equal length. For each parameter dimension, a uniform random number between zero and one is generated. The specific representative value is calculated by using the relative position of the uniform random number in the corresponding equal-length subinterval. The sampling order of each dimension is randomly permuted so that each equal-length subinterval is selected exactly once in each dimension, thereby generating a preset number of parameter vectors covering the full parameter space distribution. All parameter vectors to be estimated together constitute the initial parameter set. In the stage of establishing the surrogate model, an axisymmetric simplified mechanical model based on strain softening constitutive model is introduced. The constitutive relation of the surrogate model takes the exponential decay of the internal friction angle with the equivalent plastic strain as the core mechanism. The internal friction angle under the current equivalent plastic strain is determined as follows: first, the difference between the peak internal friction angle and the residual internal friction angle is multiplied by one and the result calculated by raising the product of the negative softening modulus and the equivalent plastic strain to the natural base. Then, the result obtained by multiplying is subtracted from the peak internal friction angle. The surrogate model defines the damage variable as the ratio of the equivalent plastic strain to the critical plastic strain. When the damage variable is greater than one, the element is deemed to have lost its bearing capacity. The state equation of the surrogate model describes the dynamic evolution relationship between the current state vector and the previous state vector and the pore pressure boundary conditions. The state vector includes stress, damage, pore pressure and displacement. The constitutive relation and state equation of the surrogate model are used to determine the damage variable.

5. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins according to claim 4, characterized in that: In step S2, the specific operation for constructing the microseismic channel observation operator is as follows: The raw waveform data in the real-time monitoring data recorded by the microseismic sensor in each construction cycle is automatically picked up and located to obtain the set of spatiotemporal coordinates of microseismic events occurring in the excavated section behind the working face in the current construction cycle. An adaptive threshold screening based on the energy attenuation rate of microseismic events was adopted to retain only microseismic events whose energy was greater than twice the moving median of the energy of microseismic events in the first three construction cycles, thus forming a set of screened microseismic events. The weighted covariance matrix is ​​calculated for the three-dimensional coordinates of all microseismic events in the selected microseismic event set. The construction process of the weighted covariance matrix is ​​as follows: First, the weighted centroid of all monitored microseismic events is determined. The spatial location coordinates of the weighted centroid are obtained by multiplying the occurrence coordinates of each microseismic event by its corresponding weight, summing them up, and then dividing by the sum of all weights. After determining the weighted centroid, each element in the weighted covariance matrix is ​​obtained by calculating the outer product of the difference between the coordinates of each microseismic event and the coordinates of the weighted centroid, multiplying the calculated outer product value by the corresponding weight of the microseismic event, and finally summing the calculated weighted outer product values ​​of all microseismic events and dividing by the sum of all weights. Subsequently, the three eigenvalues ​​of the weighted covariance matrix are extracted and arranged in descending order. The lengths of the three semi-axis are equal to the square roots of the three eigenvalues. An ellipsoid is constructed based on the lengths of the three semi-axis. The volume of the ellipsoid is equal to four-thirds of pi multiplied by the product of the three semi-axis lengths. The volume of the shear failure zone at the end of the previous construction cycle is then subtracted. The difference is the volume increment of the shear failure zone in the current construction cycle. The volume increment of the shear failure zone corresponds to the volume change of the shear failure zone formed by the set of units in the surrogate model whose damage variables exceed the critical damage threshold. The volume increment of the shear failure zone is directly used as the first component of the observation vector.

6. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 5, characterized in that: The construction process of the osmotic pressure channel observation operator is as follows: Real-time monitoring data of multiple piezometers deployed in the excavated section behind the working face during the current construction cycle are obtained. Each piezometer has its own axial distance from the working face. A control equation describing the pore pressure attenuation along the axial direction is introduced. The fractional derivative of the control equation with respect to time is equal to the diffusion coefficient multiplied by the second partial derivative of the axial distance. The frequency domain relationship of the governing equations is obtained through Laplace transform. The frequency domain relationship is as follows: the frequency domain component of the real-time monitoring data of each piezometer after Laplace transform is multiplied by the product of the negative axial distance and the first half of the frequency domain component, and then divided by the first half of the diffusion coefficient raised to the natural base. The real-time monitoring data of each piezometer is then subjected to Laplace transform, substituted into the frequency domain relationship, and then the inverse Laplace transform is performed to obtain the equivalent pore pressure at the working face. During the inversion calculation, spatial weights are assigned based on the axial distance of each piezometer from the working face. The inversion results of each piezometer are then weighted and summed according to the assigned spatial weights to obtain the equivalent pore pressure at the working face during the current construction cycle. The equivalent pore pressure at the working face is directly used as the second component of the observation vector. The observation vector consists of two components: the volume increment of the shear failure zone and the equivalent pore pressure at the working face. These components are input into the iterative set smoother in step S3 to update the set of parameters to be updated.

7. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 6, characterized in that: In step S3, the specific process of iteratively updating the parameter set to be updated multiple times using an iterative set smoother, combined with the observation vector of the current construction cycle, the microseismic channel observation operator, and the seepage pressure channel observation operator, is as follows: First, each parameter vector to be estimated in the parameter set to be updated is input into the surrogate model established in step S1, and the predicted observation vector corresponding to each parameter vector to be estimated is calculated through the microseismic channel observation operator and the pressure channel observation operator constructed in step S2. All the predicted observation vectors constitute the predicted observation set. Then, the mean vector of the parameter set to be updated and the mean vector of the predicted observation set are calculated, and parameter perturbation matrix and observation perturbation matrix are constructed respectively. Each column of the parameter perturbation matrix is ​​the square root of the difference between each estimated parameter vector and the parameter mean vector in the parameter set to be updated, divided by the preset sample number minus one. Each column of the observation perturbation matrix is ​​the square root of the difference between each predicted observation vector and the observation mean vector in the predicted observation set, divided by the preset sample number minus one. Next, a matrix decomposition operation is performed on the observed perturbation matrix, which is decomposed into three components: a left orthogonal matrix, a diagonal matrix, and a right orthogonal matrix. Singular values ​​are extracted from the diagonal matrix, and the cumulative energy ratio of these singular values ​​is calculated. The number of singular values ​​to be retained is determined based on the cumulative energy ratio, so that the proportion of energy contained in the retained singular values ​​to the total energy of all singular values ​​reaches a preset cumulative energy ratio threshold. The columns and rows corresponding to the retained singular values ​​in the left orthogonal matrix, diagonal matrix, and right orthogonal matrix are extracted together to form the truncated left orthogonal matrix, the truncated diagonal matrix, and the truncated right orthogonal matrix. Subsequently, the gain matrix of the iterative ensemble smoother is calculated using the truncated left orthogonal matrix, the truncated diagonal matrix, the truncated right orthogonal matrix, and the parameter perturbation matrix. Finally, each parameter vector in the parameter set to be updated is added to the product of the gain matrix and the difference between the observed vector of the current construction cycle and the predicted observed vector corresponding to the parameter vector to be updated, to obtain the updated parameter vector. All the updated parameter vectors to be updated constitute the parameter set to be updated for the next iteration. All operations in step S3 from calculating the predicted observed vector to obtaining the updated parameter vector are repeated. The parameter set to be updated obtained in the last iteration is the parameter set after calibration for the current construction cycle.

8. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 7, characterized in that: In step S3, after each iteration, boundary clipping is performed on each parameter vector to be estimated in the updated parameter set. The specific operation of boundary clipping is as follows: The elastic modulus, internal friction angle, and permeability coefficient in the parameter vector to be estimated are each compared with the lower limit and upper limit determined in step S1. If the value of the elastic modulus is less than the corresponding lower limit, the value of the elastic modulus is forcibly assigned to the lower limit; if the value of the internal friction angle is less than the corresponding lower limit, the value of the internal friction angle is forcibly assigned to the lower limit; if the value of the permeability coefficient is less than the corresponding lower limit, the value of the permeability coefficient is forcibly assigned to the lower limit; if the value of the elastic modulus is greater than the corresponding upper limit, the value of the elastic modulus is forcibly assigned to the upper limit. If the value of the internal friction angle is greater than the corresponding upper limit, then the value of the internal friction angle will be forcibly assigned to the upper limit. If the value of the permeability coefficient is greater than the corresponding upper limit, then the value of the permeability coefficient will be forcibly assigned to the upper limit. If the values ​​of the elastic modulus, internal friction angle, or permeability coefficient are between the corresponding lower and upper limits, then the values ​​of these parameters remain unchanged, and the parameter vector to be estimated after boundary trimming constitutes the set of parameters to be updated in the next iteration.

9. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 8, characterized in that: In step S4, the specific process of calculating the safety factor using the strength reduction method is as follows: For each parameter vector to be estimated in the parameter set after the current construction cycle calibration, the internal friction angle and elastic modulus in the parameter vector to be estimated are taken as reduction objects; an initial reduction factor is set, with an initial value of 1.0, and the reduction step size is set to 0.01; in each reduction iteration, firstly, according to the coupling strength reduction principle, the internal friction angle in the current parameter vector to be estimated is reduced. The specific processing method is as follows: obtain the tangent value of the current internal friction angle, divide the tangent value by the current reduction factor, and then take the arctangent operation of the result of the division to calculate the reduced internal friction angle; at the same time, the elastic modulus in the current parameter vector to be estimated is reduced synchronously. The specific processing method is as follows: directly divide the current elastic modulus by the reduction factor to calculate the reduced elastic modulus; Substitute the reduced internal friction angle and the reduced elastic modulus into the surrogate model established in step S1, recalculate the surrounding rock response for the next construction cycle, extract the damage variable field output by the surrogate model, and mark the elements whose damage variables exceed the critical damage threshold as shear failure elements. All shear failure elements constitute the current shear failure zone. The following steps determine the continuity status of the shear failure zone: First, all units within a vertical distance of 0.5 times the tunnel diameter near the tunnel arch are designated as the target area. Next, starting from the unit at the tunnel face, a breadth-first search algorithm is used to traverse all units marked as shear failure units to detect whether there exists a continuous spatial path consisting entirely of adjacent shear failure units. This continuous spatial path starts from the tunnel face and extends to the target area. If such a continuous path is detected, the shear failure zone is determined to have reached the arch. At this point, the reduction factor is stopped from increasing, and the current reduction factor value is recorded as the safety factor corresponding to the parameter vector to be estimated. If no continuous continuity path is detected, the reduction factor is increased by one step, and all operations of reducing the internal friction angle and elastic modulus, extrapolating the surrogate model response, and determining the continuity status are repeated. A complete reduction and judgment process is performed on all the parameter vectors to be estimated, and finally the safety factor corresponding to each parameter vector is obtained.

10. The method for risk early warning at the working face of a water conveyance tunnel based on digital twins as described in claim 9, characterized in that: The process of obtaining the collapse probability for the next construction cycle through statistics is as follows: First, for each parameter vector to be estimated in the parameter set after calibration in the current construction cycle, failure sample determination is performed based on the safety factor calculated in step S4 and the state of whether the shear failure zone extends to the crown as determined in step S4. If the safety factor corresponding to the parameter vector to be estimated is lower than the preset safety threshold, and the shear failure zone corresponding to the parameter vector to be estimated extends to the crown, then the parameter vector to be estimated is marked as a failure sample. The number of all samples marked as failure samples is counted and recorded as the number of failure samples. The total number of parameter vectors to be estimated in the parameter set after calibration in the current construction cycle is obtained and recorded as the total number of samples. The number of failure samples is divided by the total number of samples to calculate the collapse probability of the next construction cycle.

Images