A hard rock cavern fracture state dynamic evaluation method with tunneling optimization

Abstract

The application discloses a hard rock cavern fracture state dynamic evaluation method with tunneling optimization, and belongs to the technical field of underground engineering and intelligent evaluation. The application generates a working condition-fracture response database to form an offline digital base map through a meso-macro coupling numerical model; a multi-branch deep neural network is used to construct a physical information guided hybrid proxy model, the digital base map is trained and embedded with a physical prior guided mechanism; microseismic, acoustic wave and displacement meter signals are collected in real time during construction, and are processed into three types of feature data, which are input into the model to output fracture depth probability distribution and fracture degree index in parallel; after each excavation cycle, the measured fracture depth is extracted as hard data by using a drilling camera or a while-drilling parameter, the posterior probability is updated through Bayesian inference, and the next cycle is iterated, so that the fracture state evaluation with tunneling dynamic optimization is realized. The application realizes real-time, probabilistic and adaptive evaluation of the surrounding rock fracture depth and degree, and effectively improves the scientificity of support design and the safety of construction.

Description

Technical Field

[0001] This invention belongs to the field of underground engineering rock mechanics and intelligent evaluation technology, specifically involving a dynamic evaluation method for the fracture state of hard rock caverns that is optimized as tunneling progresses. Background Technology

[0002] During the excavation and unloading process of deep, high-stress hard rock caverns, the original stress equilibrium state of the surrounding rock is severely disturbed, leading to stress redistribution. When the localized concentrated stress exceeds the rock mass's own strength, the surrounding rock will exhibit fracture evolution behaviors such as micro-crack initiation, propagation, and penetration, thus forming an excavation damage zone. Accurately obtaining the fracture depth and degree of the surrounding rock is the core basis for assessing the stability of the surrounding rock, scientifically designing anchor bolt support, and preventing dynamic disasters such as rockbursts. Currently, numerical simulation, in-situ geophysical testing, or data-driven models are mainly used in engineering sites to quantitatively estimate or classify the fracture state of the surrounding rock.

[0003] However, the aforementioned methods have significant limitations in dynamic evaluation during tunneling cycles. Due to the complex underground excavation environment and the highly concealed characteristics of the surrounding rock, traditional numerical simulation methods rely on complex constitutive parameter settings, resulting in lengthy large-scale element iterative calculations that are difficult to meet the real-time requirements of on-site construction. In-situ testing methods such as geophysical exploration or borehole photography are limited by the testing cycle and can often only be implemented after the current tunneling cycle is completed. This delayed detection cannot provide predictive basis for the pre-support design of the next tunneling cycle. Furthermore, existing data-driven models, when faced with the scarcity of deep hard rock fracture samples, mostly output deterministic state classifications, failing to quantify the confidence interval of the prediction results, and lacking a mechanism to effectively absorb new observation data during tunneling for adaptive iteration. Therefore, it is urgent to improve existing technologies to solve the technical problem that the surrounding rock fracture state cannot be dynamically, probabilistically, and adaptively optimized during the excavation of deep high-stress hard rock caverns. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a dynamic evaluation method for the fracture state of hard rock caverns that is optimized as tunneling progresses. This method addresses the technical problem that the fracture state of the surrounding rock cannot be dynamically, probabilistically, and adaptively optimized during the excavation of deep, high-stress hard rock caverns. It enables real-time probabilistic dynamic evaluation of the depth and extent of the surrounding rock fracture and allows the evaluation results to evolve adaptively as construction progresses.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: A dynamic evaluation method for the fracture state of hard rock caverns, optimized during tunneling, includes the following steps: Step S100: Based on the geological exploration parameters of the target cavern, establish a working condition-fracture response database generated in batches by the micro-macro coupled numerical model, extract the working condition parameters, true values ​​of surrounding rock fracture depth, three-dimensional microseismic damage density field, isosurface of acoustic wave velocity reduction amplitude and multi-point displacement gradient characteristic data corresponding to each working condition in the database, and form an offline digital base map. Step S200: Based on a multi-branch deep neural network architecture, construct a hybrid proxy model guided by physical information, and embed a physical prior guidance mechanism into the model architecture; use the three-dimensional microseismic damage density field, the acoustic wave velocity reduction isosurface, the multi-point displacement gradient feature data and the corresponding working condition parameters in the offline digital base map as training inputs, and use the true value of the surrounding rock fracture depth and the predefined fracture degree index as training outputs to train the hybrid proxy model; Step S300: Acquire the microseismic monitoring signal, acoustic test signal and multi-point displacement meter signal collected in real time during the excavation of the target cavern; process the microseismic monitoring signal into a three-dimensional microseismic damage density field; process the acoustic test signal into an isosurface of acoustic wave velocity reduction relative to the reference wave velocity; and process the multi-point displacement meter signal into multi-point displacement gradient feature data. Step S400: In each excavation cycle of the target cavern, the three-dimensional microseismic damage density field, the acoustic wave velocity reduction isosurface, and the multi-point displacement gradient feature data obtained in step S300, as well as the corresponding working condition parameters, are input into the trained hybrid proxy model. The hybrid proxy model outputs in parallel the probability distribution of the rupture depth and the rupture degree index of the current state of the target cavern. Step S500: After each excavation cycle of the target cavern is completed, the borehole camera image data or the anchor trolley drilling parameter data corresponding to the cycle are obtained, and the measured fracture depth is extracted as the hard data observation value. Step S600: The probability distribution of fracture depth output in step S400 is used as the prior probability distribution, and the hard data observations obtained in step S500 are used as the likelihood basis. Bayesian inference is used to update the probability distribution of fracture depth and fracture degree index of the target cavern in the current state. The posterior probability distribution is then used as the prior probability distribution of step S400 in the next excavation cycle, thereby realizing the evaluation of fracture state through dynamic iteration with the tunneling cycle.

[0006] Further, the micro-macro coupled numerical model in step S100 is a particle discrete element and finite difference coupled model. The offline digital base map generation process includes: extracting three parameters—surrounding rock quality index, strength-stress ratio, and excavation-support ratio—based on the geological survey parameters of the target cavern; performing hierarchical orthogonal combination of the three parameters to form a standard working condition set; for each standard working condition in the standard working condition set, performing offline simulation of excavation unloading and fracture evolution using the particle discrete element and finite difference coupled model; extracting the true value of the surrounding rock fracture depth corresponding to the working condition from the simulation results; performing three-dimensional spatial meshing processing; and generating an offline digital base map composed of the true value of the surrounding rock fracture depth, the three-dimensional microseismic damage density field, the isosurface of the acoustic wave velocity reduction amplitude, and the multi-point displacement gradient feature data.

[0007] Furthermore, in step S200, the hybrid proxy model internally constructs a depth prediction branch and a degree prediction branch; the physical prior guidance mechanism is configured to: cause the acoustic wave velocity reduction isosurface and the multi-point displacement gradient feature data related to depth prediction to flow to the depth prediction branch, and cause the three-dimensional microseismic damage density field related to degree prediction to flow to the degree prediction branch.

[0008] Furthermore, the physical prior guidance mechanism includes setting a diversion layer and an attention guidance layer in the input mapping path of the hybrid proxy model; the diversion layer is used to allocate the acoustic wave velocity reduction isosurface and the multi-point displacement gradient feature data to the depth prediction branch, and to allocate the three-dimensional microseismic damage density field to the degree prediction branch; the attention guidance layer is used to learn and strengthen the correlation weights between different feature dimensions within each prediction branch and the training output.

[0009] Furthermore, during the training process, the loss function of the hybrid proxy model includes a depth prediction loss term, a severity prediction loss term, and a physical consistency constraint term. The physical consistency constraint term is used to penalize the deviation between the mean predicted rupture depth of the hybrid proxy model and the isosurface of the acoustic wave velocity reduction amplitude, as well as the deviation between the predicted rupture severity index of the hybrid proxy model and the three-dimensional microseismic damage density field. During training, the network parameters of the hybrid proxy model are updated by optimizing the weighted sum of the depth prediction loss term, the severity prediction loss term, and the physical consistency constraint term.

[0010] Further, in step S400, the probability distribution of fracture depth and the fracture severity index are determined in the following way: the probability distribution of fracture depth is obtained based on the probability density function constructed from the predicted mean and predicted variance output by the hybrid surrogate model; the fracture severity index is a dimensionless index established based on the weighted comprehensive analysis of multiple physical quantities such as fracture density, crack slip energy release rate and P-wave velocity reduction, which is used to characterize the damage and fracture evolution state of the surrounding rock at different stages.

[0011] Furthermore, the extraction method of the hard data observations in step S500 is as follows: for the borehole camera image data, the measured fracture depth is determined by identifying the abrupt change point of fracture density or the location of the fracture surface along the borehole depth direction; for the anchor trolley drilling parameter data, the measured fracture depth is determined by identifying the abrupt change point of drilling speed at a certain borehole depth position.

[0012] Furthermore, in step S600, the Bayesian inference employs a sequential Bayesian filtering or particle filtering algorithm. After acquiring new hard data observations each time, the posterior probability distribution of the previous cycle is used as the prior probability distribution of the current cycle. The likelihood function is calculated using the hard data observations of the current cycle, and the posterior probability distribution of the current cycle is calculated using the Bayesian formula. This allows the evaluation results of the fracture depth and fracture degree index of the current state of the target cavern to gradually converge to a narrower confidence interval as the tunneling cycle increases.

[0013] Furthermore, the method includes step S700: pushing the three-dimensional microseismic damage density field obtained in step S300, the probability distribution of rupture depth and the rupture degree index output in step S400, and the posterior probability distribution updated in step S600 to a three-dimensional visualization platform. The platform displays the spatial distribution of the probability contour lines of rupture depth exceeding the preset anchor length and the rupture degree index in the form of a dynamic wrapping surface on the three-dimensional geological model of the target cavern, and synchronously displays the evolution curves of rupture depth and the rupture degree index and the changes of the confidence interval updated by Bayes on the time axis.

[0014] This invention discloses a dynamic evaluation method for the fracture state of hard rock caverns optimized during tunneling. Its core lies in generating a digital base map offline in batches using a coupled meso-macro numerical model to create a database of working conditions and fracture responses. This base map is then used to train a hybrid surrogate model embedded with a physical prior guidance mechanism. Combined with hard data extracted from borehole cameras or drilling parameters during each excavation cycle, Bayesian inference is used to update the posterior probability distribution of the prior, thus constructing an adaptive evaluation closed-loop system of "offline digital base map pre-training - online real-time feature prediction - dynamic updating of hard data during tunneling." Compared to existing technologies, this invention has significant advantages: First, a working condition-fracture response database is established using a particle discrete element and finite difference coupled model, and an offline digital base map is formed. This allows the hybrid proxy model to internalize the physical laws of surrounding rock fracture evolution before construction. During online inference, there is no need to call the time-consuming numerical simulation again, which greatly improves the timeliness of fracture state evaluation during dynamic excavation and meets the real-time requirements of on-site tunneling.

[0015] Secondly, a physical prior guidance mechanism, including a diversion layer and an attention guidance layer, is embedded in the hybrid agent model architecture. The diversion layer distributes the isosurface of acoustic wave velocity reduction and multi-point displacement gradient feature data to the depth prediction branch, and the three-dimensional microseismic damage density field to the degree prediction branch. This avoids the drawback of pure data-driven models that rely solely on statistical correlation. It also outputs the probability distribution of fracture depth and the dimensionless fracture degree index that characterizes the damage evolution state in parallel, so that the evaluation results can effectively quantify cognitive uncertainty while possessing rock mechanics interpretability.

[0016] Furthermore, a high-frequency dynamic update mechanism driven by the tunneling cycle is adopted. After each excavation cycle, the Bayesian inference update is performed using physically measurable hard data observations such as borehole cameras or drilling parameters as the likelihood basis. This makes the posterior probability distribution of the previous cycle the prior probability distribution of the current cycle. This enables the understanding of the fracture depth and degree to adaptively converge to the real state as the tunneling construction progresses, solving the defect of traditional static evaluation models that cannot self-adjust as the surrounding rock conditions migrate along the tunnel line.

[0017] Finally, by pushing real-time multi-source features, probability distributions, and updated posterior probability distributions to a 3D visualization platform for dynamic display of the wrapping surface and evolution curves, the abstract probability evaluation results are transformed into an intuitive engineering risk picture, significantly improving the scientificity and intuitiveness of support adjustment decisions in high-risk areas. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is an overall flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the offline digital base map construction result of the present invention; Figure 3 This is a schematic diagram illustrating the training and verification of the physical information-guided hybrid agent model of the present invention. Figure 4This is a schematic diagram of the dynamic iteration of Bayesian inference in this invention; Figure 5 This is a schematic diagram of the output of the 3D visualization platform of the present invention; Figure 6 This is a schematic diagram showing the convergence of the time-series dynamic evolution curve and the Bayesian confidence interval of the present invention; Figure 7 This is a schematic diagram of the hard data observations for borehole imaging and drilling parameter extraction according to the present invention. Detailed Implementation

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

[0021] This application provides a method for dynamically evaluating the fracture state of hard rock caverns that is optimized during tunneling, such as... Figure 1 As shown, the method includes the following steps: Step S100: Based on the geological exploration parameters of the target cavern, establish a working condition-fracture response database generated in batches by the micro-macro coupled numerical model, extract the working condition parameters, true values ​​of surrounding rock fracture depth, three-dimensional microseismic damage density field, isosurface of acoustic wave velocity reduction amplitude and multi-point displacement gradient characteristic data corresponding to each working condition in the database, and form an offline digital base map. Step S200: Based on a multi-branch deep neural network architecture, construct a hybrid proxy model guided by physical information, and embed a physical prior guidance mechanism into the model architecture; use the three-dimensional microseismic damage density field, the acoustic wave velocity reduction isosurface, the multi-point displacement gradient feature data and the corresponding working condition parameters in the offline digital base map as training inputs, and use the true value of the surrounding rock fracture depth and the predefined fracture degree index as training outputs to train the hybrid proxy model; Step S300: Acquire the microseismic monitoring signal, acoustic test signal and multi-point displacement meter signal collected in real time during the excavation of the target cavern; process the microseismic monitoring signal into a three-dimensional microseismic damage density field; process the acoustic test signal into an isosurface of acoustic wave velocity reduction relative to the reference wave velocity; and process the multi-point displacement meter signal into multi-point displacement gradient feature data. Step S400: In each excavation cycle of the target cavern, the three-dimensional microseismic damage density field, the acoustic wave velocity reduction isosurface, and the multi-point displacement gradient feature data obtained in step S300, as well as the corresponding working condition parameters, are input into the trained hybrid proxy model. The hybrid proxy model outputs in parallel the probability distribution of the rupture depth and the rupture degree index of the current state of the target cavern. Step S500: After each excavation cycle of the target cavern is completed, the borehole camera image data or the anchor trolley drilling parameter data corresponding to the cycle are obtained, and the measured fracture depth is extracted as the hard data observation value. Step S600: The probability distribution of fracture depth output in step S400 is used as the prior probability distribution, and the hard data observations obtained in step S500 are used as the likelihood basis. Bayesian inference is used to update the probability distribution of fracture depth and fracture degree index of the target cavern in the current state. The posterior probability distribution is then used as the prior probability distribution of step S400 in the next excavation cycle, thereby realizing the evaluation of fracture state through dynamic iteration with the tunneling cycle.

[0022] It should be noted that the hybrid proxy model is based on a multi-branch deep neural network architecture. "Multi-branch" refers to the network's simultaneous inclusion of a depth prediction branch and a severity prediction branch, used to output the probability distribution of fracture depth and the dimensionless exponent of fracture severity, respectively. The "physical prior guidance mechanism" uses a diversion layer to forcibly assign monitoring features with different physical meanings to the corresponding prediction branches (wave velocity reduction and displacement gradient features flow to the depth branch, microseismic damage density features flow to the severity branch), thus avoiding the drawbacks of purely data-driven models that rely solely on statistical correlation. Furthermore, the "working condition parameters" extracted in step S100 are a collective term for three parameters: surrounding rock quality index, strength-stress ratio, and excavation-support ratio. The surrounding rock quality index reflects the baseline properties of the rock mass's resistance to fracture; the strength-stress ratio characterizes the driving intensity of excavation unloading; and the excavation-support ratio quantifies the constraint effect of human intervention on the evolution of surrounding rock fracture. These three parameters are independent yet work together, constituting the minimum complete parameter set that determines the fracture state of the surrounding rock. In practical applications, the surrounding rock quality indicators can be directly obtained from geological survey data (with minor adjustments as needed based on the exposure of the tunnel face), the strength-stress ratio is calculated based on the geostress test report and rock mechanics test results, and the excavation-support ratio is dynamically determined by the burr span and support lag distance in the construction organization design. These three types of parameters are basically known before construction or can be determined in advance in each cycle. Therefore, they serve as fixed inputs for the hybrid proxy model to distinguish the differences in fracture response under different geological sections and construction conditions.

[0023] Specifically, before the formal excavation begins, steps S100 and S200 are completed offline. Their core purpose is to perform high-density mechanical simulations using a high-performance computer, thereby creating an offline digital base map capable of covering various complex geological conditions. This digital base map is then used to train an intelligent hybrid surrogate model with prior knowledge of rock mechanics. Once construction begins, steps S300 to S600 are executed in real-time on-site. In each excavation cycle, real-time geophysical and monitoring data collected on-site is first inferred within seconds by the hybrid surrogate model. This data, along with working condition parameters, is used as input. The model outputs the probability distribution of the fracture depth and the fracture degree index of the current surrounding rock through its internal depth prediction and degree prediction branches. Subsequently, the measured fracture depth extracted from borehole photography or drilling parameters after the current excavation cycle is used as hard data observations, and the prediction results of the surrogate model are adaptively corrected using Bayesian formulas. Since the posterior probability distribution calculated in the previous cycle is directly used as the prior condition for the reasoning in the next cycle, the entire system forms a dynamic evolutionary closed loop that can gradually narrow the evaluation confidence interval as construction progresses and data accumulates, thus solving the practical contradiction of "not being able to calculate quickly" and "not being able to calculate accurately" in the evaluation of the excavation status of hard rock caverns.

[0024] The present invention further proposes that the micro-macro coupled numerical model in step S100 is a particle discrete element and finite difference coupled model. The offline digital base map generation process includes: extracting three parameters—surrounding rock quality index, strength-stress ratio, and excavation-support ratio—based on the geological survey parameters of the target cavern; performing hierarchical orthogonal combination of the three parameters to form a standard working condition set; for each standard working condition in the standard working condition set, performing offline simulation of excavation unloading and fracture evolution using the particle discrete element and finite difference coupled model; extracting the true value of the surrounding rock fracture depth corresponding to the working condition from the simulation results; performing three-dimensional spatial meshing processing; and generating an offline digital base map composed of the true value of the surrounding rock fracture depth, the three-dimensional microseismic damage density field, the isosurface of the acoustic wave velocity reduction amplitude, and the multi-point displacement gradient feature data.

[0025] The micro-macro coupled numerical model refers to the seamless integration of particle discrete element methods (such as PFC), which can simulate microscopic constitutive behaviors such as rock grain fracture and slip, with finite difference methods (such as FLAC), which can solve large-scale continuous medium stress fields. Since the excavation and fracturing of deep, high-stress hard rock is essentially a multi-scale mechanical process involving the initiation, convergence, and macroscopic spalling of microcracks, a single continuous medium model cannot realistically reproduce the crack propagation mechanism, while a full-scale discrete element model is computationally too intensive. This scheme achieves a balance between computational accuracy and efficiency by arranging discrete element particle meshes in the strongly unloaded fracture zone around the cavern and finite difference meshes in the distant stable zone. Based on this, the extracted surrounding rock quality index, strength-stress ratio, and excavation-support ratio—three core parameters—together constitute the controlling mechanical boundary that determines the surrounding rock fracture behavior. The standard working condition set designed through hierarchical orthogonal combination can cover, to the greatest extent possible, various extreme or conventional geological scenarios that the target cavern may encounter during future excavation with the most economical combined sample size. After the simulation is completed, all physical quantities are mapped to a unified three-dimensional spatial grid, thereby transforming the discrete mechanical response into structured spatial tensor data. The resulting offline digital base map actually serves as a "virtual physical textbook" for the intelligent proxy model, providing solid and detailed physical annotation data support for eliminating the problem of insufficient samples in the neural network training process.

[0026] The present invention further proposes that, in step S200, the hybrid proxy model internally constructs a depth prediction branch and a degree prediction branch; the physical prior guidance mechanism is configured to: cause the acoustic wave velocity reduction isosurface and the multi-point displacement gradient feature data related to depth prediction to flow to the depth prediction branch, and cause the three-dimensional microseismic damage density field related to degree prediction to flow to the degree prediction branch.

[0027] In the internal architecture design of the hybrid proxy model, the division between the depth prediction branch and the degree prediction branch is not an arbitrary choice of network structure, but rather based on a fundamental understanding in rock mechanics: the evolution of fracture depth is mainly controlled by the redistribution of surrounding rock stress and deformation gradient, while the degree of fracture depends more on the cumulative damage and energy release of microcracks. Based on this physical prior, this invention configures the physical prior guidance mechanism as a specific data flow constraint: the isosurface of acoustic wave velocity reduction directly reflects the radial damage gradient change of the surrounding rock. Areas with larger wave velocity reduction usually indicate more severe structural deterioration. Therefore, it, together with multi-point displacement gradient feature data (reflecting the relative deformation of the rock mass along the depth), constitutes the core basis for judging fracture depth. These two types of features are guided to the depth prediction branch; while the three-dimensional microseismic damage density field is obtained by real-time location of microseismic events and inversion of their energy release. It directly characterizes the density of microcrack activity and energy release intensity within the surrounding rock, and has a more direct physical correlation with the degree of fracture. Therefore, it is guided to the degree prediction branch. This configuration allows the surrogate model to consciously focus on the physical causal relationship between different input features and the output target during the learning process, rather than simply mining statistical correlations, thereby significantly improving the model's generalization ability in unseen conditions.

[0028] The present invention further proposes that the physical prior guidance mechanism includes setting a diversion layer and an attention guidance layer in the input mapping path of the hybrid proxy model; the diversion layer is used to allocate the acoustic wave velocity reduction isosurface and the multi-point displacement gradient feature data to the depth prediction branch, and to allocate the three-dimensional microseismic damage density field to the degree prediction branch; the attention guidance layer is used to learn and strengthen the correlation weights between different feature dimensions within each prediction branch and the training output.

[0029] To enable the aforementioned physical prior guidance mechanism to be concretely implemented in the neural network architecture, this invention sets up a split layer and an attention guidance layer in the input mapping path. The split layer acts like an intelligent "data router." After receiving the multidimensional feature tensors extracted by the convolutional coding layer, it allocates paths according to the physical meaning of different features: the encoded feature vectors from the acoustic wave velocity reduction isosurface and displacement gradient feature curves are directed to the input port of the depth prediction branch; the encoded feature vectors from the three-dimensional microseismic damage density field are directed to the input port of the degree prediction branch. This rigid path allocation mechanism enforces the injection of physical priors from the network structure. The attention guidance layer is set inside each prediction branch, and its function is to dynamically weight the multidimensional features allocated to that branch. Specifically, the attention layer learns a set of trainable feature weight coefficients to automatically identify which feature dimensions are more critical to the prediction result under the current input conditions and assign them higher attention weights. For example, in the depth prediction branch, when the wave velocity reduction isosurface shows a significant reduction gradient within a certain radial depth range, the attention mechanism automatically strengthens the influence of the feature dimension corresponding to that depth on the output. Through the collaborative work of the diversion layer and the attention guidance layer, the proxy model both follows the guidance of physical laws and retains the adaptive capabilities of data-driven approaches.

[0030] The present invention further proposes that, during the training process, the loss function of the hybrid proxy model includes a depth prediction loss term, a severity prediction loss term, and a physical consistency constraint term. The physical consistency constraint term is used to penalize the deviation between the mean predicted rupture depth of the hybrid proxy model and the isosurface of the acoustic wave velocity reduction amplitude, and to penalize the deviation between the predicted rupture severity index of the hybrid proxy model and the three-dimensional microseismic damage density field. During training, the network parameters of the hybrid proxy model are updated by optimizing the weighted sum of the depth prediction loss term, the severity prediction loss term, and the physical consistency constraint term.

[0031] It is understandable that the training process of a hybrid surrogate model is essentially a parameter optimization process, and the design of the loss function directly determines the model's learning objective and final behavior. This invention innovatively introduces a physical consistency constraint term in addition to the conventional depth and severity prediction losses. Its design aims to penalize predictions that, while statistically fitting the training data, clearly violate fundamental laws of rock mechanics. Specifically, the physical consistency constraint term includes two sub-terms: the first sub-term focuses on the consistency between the predicted fracture depth and the spatial distribution of the acoustic wave velocity reduction. According to the basic principles of rock physics, the boundary of the surrounding rock fracture zone should roughly coincide with the location of the abrupt change in wave velocity reduction. If the model predicts a fracture depth far from the location of the maximum gradient of wave velocity reduction, it indicates that the prediction violates physical laws, and the physical consistency constraint term will penalize this. The second sub-term focuses on the degree of deviation between the predicted fracture severity and the microseismic damage density field. The cumulative energy release in high-density regions of the microseismic damage density field reflects the severity of surrounding rock fracturing, which should physically be positively correlated with the fracturing degree index. If the model's predicted fracturing degree index is inconsistent with the spatial concentration of microseismic energy release, it will also be penalized. By introducing these two physical constraints into the loss function, the model, during training, must not only minimize the error between the predicted and labeled values ​​but also ensure that its predictions conform to basic physical laws, thereby significantly improving the model's reliability and engineering applicability.

[0032] The present invention further proposes that, in step S400, the probability distribution of fracture depth and the fracture degree index are determined in the following manner: the probability distribution of fracture depth is obtained based on the probability density function constructed from the predicted mean and predicted variance output by the hybrid surrogate model; the fracture degree index is a dimensionless index established based on the weighted comprehensive analysis of multiple physical quantities such as fracture density, crack slip energy release rate and P-wave velocity reduction, which is used to characterize the damage and fracture evolution state of the surrounding rock at different stages.

[0033] The fracture depth output by the surrogate model is not a fixed numerical value, but a complete probability distribution, which is a key feature distinguishing this invention from existing deterministic evaluation methods. Specifically, two output neurons are set at the end of the fully connected layer of the depth prediction branch, corresponding to the prediction mean and prediction variance, respectively. A Gaussian probability density function constructed based on these two neurons fully describes the uncertainty of the fracture depth, allowing decision-makers to calculate the probability that the fracture depth falls within any interval. The fracture severity index integrates three physical quantities: fracture density, fracture slip energy release rate, and P-wave velocity reduction. A dimensionless index between 0 and 1 is obtained through weighted aggregation, mapping the continuous evolution of the surrounding rock from an elastic state to a rockburst state onto a unified scale. Regarding the determination of the weighting coefficients, objective weighting can be performed using entropy weighting based on sample data from offline digital base maps. For example, the three physical quantities are first normalized, and the entropy value of each index is calculated. The smaller the entropy value, the greater the dispersion of the index and the higher its weight. Alternatively, expert surveys or analytic hierarchy process (AHP) can be used to adjust the weighting coefficients based on specific engineering experience. Once determined, the weighting coefficients remain unchanged within the same project.

[0034] The present invention further proposes that the extraction method of the hard data observation value in step S500 is as follows: for the borehole camera image data, the measured fracture depth is determined by identifying the fracture density abrupt change point or fracture surface position along the borehole depth direction; for the anchor trolley drilling parameter data, the measured fracture depth is determined by identifying the abrupt change point of the drilling speed at a certain borehole depth position.

[0035] Hard data observations refer to real in-situ measurement results obtained through direct physical contact or optical visualization at the excavation site, with almost no cognitive bias. These results serve as the basis for correcting prediction biases in surrogate models. This invention provides two mutually redundant methods for extracting measured fracture depth. The first method is based on borehole camera image data: In each completed anchor bolt hole or dedicated observation hole of each excavation cycle, a panoramic camera is lowered to the bottom of the hole and raised at a constant speed to acquire a 360° sequence of borehole wall images continuously unfolding along the hole depth. By identifying and statistically analyzing fractures frame by frame in the image sequence, a fracture linear density curve along the hole depth direction can be plotted. Because fractures are densely developed in the fracture zone, while fractures are sparse in the intact rock mass area, there is a clear density abrupt change point between the two. The depth corresponding to this abrupt change point is determined as the measured fracture depth. The second method is based on the drilling parameters of the anchor bolt trolley: During the drilling process, the anchor bolt trolley records parameters such as drilling speed, propulsion pressure, and rotational torque in real time. When the drill bit is in the fractured zone, the surrounding rock is broken and the resistance is low, resulting in a faster drilling speed. However, as the drill bit passes through the fractured zone and enters the intact rock mass, the resistance increases sharply, and the drilling speed exhibits a step-like decrease (or is accompanied by a step-like increase in propulsion pressure). The fracture depth can also be accurately located by identifying the abrupt change points on the drilling speed curve. These two methods can be mutually verified on-site. When one method fails to obtain effective data due to borehole collapse, rock dust blockage, or other reasons, the other method can serve as a supplement, ensuring reliable hard data observations are obtained in each drilling cycle.

[0036] The present invention further proposes that the Bayesian inference in step S600 adopts a sequential Bayesian filtering or particle filtering algorithm. After each acquisition of new hard data observations, the posterior probability distribution of the previous cycle is used as the prior probability distribution of the current cycle. The likelihood function is calculated using the hard data observations of the current cycle, and the posterior probability distribution of the current cycle is calculated using the Bayesian formula. This allows the evaluation results of the fracture depth and fracture degree index of the current state of the target cavern to gradually converge to a narrower confidence interval as the tunneling cycle increases.

[0037] In tunnel and underground engineering, as the tunnel face advances, the spatial variability of geological conditions often leads to increasingly larger errors in traditional static model predictions as the excavation progresses. Considering that the evaluation process of this invention advances sequentially with each excavation cycle, this invention employs a sequential Bayesian filtering or particle filtering algorithm to achieve this dynamic update. Its core operating principle is to treat the excavation process as a dynamic chain of temporal information flow. In any given excavation cycle, the model no longer isolates and blindly guesses the outcome. Instead, it first inherits and extracts the surrounding rock condition knowledge refined from the previous cycle (i.e., the posterior probability distribution), using it as the basis (prior probability distribution) for the current cycle. Specifically, for each excavation cycle, step S500 is first called to obtain hard data observations, and a likelihood function (usually a Gaussian likelihood function) is constructed based on the pre-calibrated observation error. Then, the fracture depth probability distribution output in step S400 is used as the prior, multiplied by the likelihood function, and normalized to obtain the posterior probability distribution for the current cycle. Since the prior and likelihood in this application can be reasonably assumed to be Gaussian distributions, and the posterior distribution also has an analytical expression, the calculation can be completed quickly without complex numerical sampling.

[0038] When the error distribution of the observed data exhibits non-Gaussian characteristics or the model exhibits nonlinearity, a particle filtering algorithm can be used to implement Bayesian updates. The specific steps include: randomly selecting N particles (usually N is 500-1000) from the prior probability distribution, each particle representing a possible fracture depth value, and assigning each particle an equal initial weight of 1 / N. Predicting the state of each particle based on the state transition model (since the change in the surrounding rock fracture state between adjacent cycles is relatively slow, a random walk model is typically used, i.e., new state = old state + Gaussian white noise), resulting in a predicted particle set. After obtaining the hard data observations for the current cycle, calculating the likelihood function value of each particle (i.e., the probability density of observing hard data at the fracture depth represented by the particle), updating the weight of each particle to the original weight multiplied by the likelihood function value, and then normalizing the weights of all particles so that the sum of the weights is 1. Since after several updates, the weights of a few particles may become extremely large, while the weights of most particles may approach zero (i.e., particle degeneration), resampling is required at this point. The core idea of ​​resampling is to resample N particles based on their weight distribution, with higher-weighted particles being replicated more frequently and lower-weighted particles being discarded. After resampling, the weights of all particles are reset to 1 / N. The weighted mean of the posterior particle set is used as the posterior estimate of the fracture depth, and the weighted variance of the particle set is used as a measure of posterior uncertainty. Simultaneously, empirical quantiles of the posterior probability distribution (such as the 2.5% and 97.5% quantiles as the 95% confidence interval) can be output. Through these steps, particle filtering can still approximate the true posterior probability distribution even when the prior distribution is non-Gaussian or the observation model is non-linear. Furthermore, as the tunneling cycle progresses, the particle set gradually clusters near the true fracture depth, and the confidence interval gradually narrows.

[0039] The present invention further proposes to include step S700: pushing the three-dimensional microseismic damage density field obtained in step S300, the probability distribution of rupture depth and the rupture degree index output in step S400, and the updated posterior probability distribution in step S600 to a three-dimensional visualization platform, displaying the probability contour lines of rupture depth exceeding the preset anchor length and the spatial distribution of the rupture degree index in the form of a dynamic wrapping surface on the three-dimensional geological model of the target cavern, and synchronously displaying the evolution curves of the rupture depth and the rupture degree index and the changes of the Bayesian updated confidence interval on the time axis.

[0040] To enable engineers to intuitively understand and effectively use the abstract probability assessment results, this invention pushes core data to a 3D visualization platform for dynamic display. Based on the design excavation outline of the target cavern and geological survey data, the 3D visualization platform constructs a 3D geological model including stratigraphic interfaces, main structural surfaces, and high-stress area boundaries. On this model, the system presents a probability isosurface for fracture depth exceeding the preset anchor bolt length in the form of a dynamic wrapping surface: by integrating the posterior probability distribution of fracture depth in various directions around the cavern, the probability value of fracture depth exceeding the designed anchor bolt length in each direction is obtained. These probability values ​​are spatially interpolated to generate an isosurface, forming a "risk wrapping surface." Areas with higher probabilities are highlighted with warm colors (such as red and orange) to visually indicate that the anchor bolt length in that area may be insufficient. Simultaneously, the fracture degree index output in step S400 is spatially interpolated and colored on the cavern wall surface to form a heatmap of the fracture degree. The timeline linkage function allows users to trace the evolution of the entire tunneling process: as the user drags the time slider, the system synchronously updates the distribution of microseismic damage density field, risk envelope surface, and fracture degree for the corresponding construction date, and displays the posterior mean of fracture depth and its 95% confidence interval convergence curve with the number of cycles. When the evaluation result triggers the radio frequency index threshold, the system automatically overlays a flashing warning mark and pushes an early warning message, realizing a complete closed loop from data acquisition and intelligent evaluation to risk warning and auxiliary decision-making.

[0041] Example 1 To more clearly understand the specific execution logic and operational details of the present invention, this embodiment provides a detailed breakdown of the complete workflow from offline base map construction, online model training, real-time inference, to Bayesian dynamic updates. The specific execution steps of this method are as follows: Step S100: Construct an offline digital base map The core of this step is to use the computationally time-consuming micro-macro coupled numerical model in the existing technology as an offline data generator to build a massive labeled database containing rock mechanics mechanisms in batches before construction, laying the physical foundation for subsequent training of surrogate models.

[0042] S101. Determine key operating condition parameters and their classification criteria. Based on the geological survey report of the target cavern, extract three parameters that have a dominant controlling effect on the surrounding rock fracture response: The quality indicators of the surrounding rock are assessed using the Basic Quality Index (BQ) grading method or the Rock Mass Quality Rating (RMR) method, and are divided into several grades according to the preset grading standards. The strength-stress ratio is defined as the ratio of the uniaxial compressive strength of the rock block to the maximum initial principal stress of the cavern, and is graded according to a preset threshold range. The excavation-support ratio is defined as a dimensionless comprehensive parameter, which is equal to the ratio of the current tunnel span to the distance from the immediate support application location to the tunnel face, and is classified according to a preset threshold range.

[0043] S102. Generate a standard working condition set through orthogonal combination. The graded values ​​of the three parameters are orthogonally combined using three factors at multiple levels to form a standard working condition set covering the working conditions that the target cavern may encounter. Each standard working condition is uniquely determined by a set of ternaries: {surrounding rock quality grade, strength-stress ratio range, excavation-support ratio range}.

[0044] S103. Conduct meso-macro coupled numerical simulations for each working condition. For each standard working condition, establish a corresponding particle discrete element and finite difference coupled model. The coupled model is performed before construction and is not limited by construction time. The specific simulation process is as follows: 1. Based on the surrounding rock quality parameters, such as the fractal characteristics of granite grain embedding and mineral composition content, a three-dimensional microstructure model of granite reflecting the microstructure of the rock mass is constructed using the PFC particle flow platform.

[0045] 2. Using the strength-stress ratio and excavation-support ratio as boundary conditions and excavation procedure input parameters, a layered excavation and unloading simulation of the tunnel was conducted on the FLAC platform.

[0046] 3. Embed the PFC model into the fracture region of the FLAC model, and obtain the full picture of the surrounding rock fracture state under the final stable state through iterative calculation.

[0047] For each standard operating condition, the data volume output by the coupled numerical simulation includes: the crack density distribution on the three-dimensional grid nodes; the energy release value of each microseismic source point in three-dimensional space; the calculated value of the P-wave velocity on the three-dimensional grid nodes; and the depth-displacement gradient data extracted along typical borehole locations.

[0048] S104. Extract the physical fingerprint field and form labeled samples. Post-process the simulation results for each of the standard operating conditions: 1. Extract the true value of the surrounding rock fracture depth corresponding to this working condition. The fracture depth is defined as the radial distance from the tunnel wall to the point where the fracture density suddenly drops to the level of the baseline fracture density of the intact rock, and is extracted along different directions around the tunnel.

[0049] 2. Constructing the microseismic damage density field. The energy release values ​​at the microseismic source points are interpolated using kernel density estimation on a three-dimensional spatial grid to generate a three-dimensional density field characterizing the continuous spatial distribution of microseismic damage. The kernel density estimation expression is as follows: In the formula, For spatial points The estimated microseismic damage density at the location, This represents the total number of microseismic events. For bandwidth parameters, For Gaussian kernel function, For the first Spatial location of microseismic events For the first Normalized release energy of a microseismic event.

[0050] 3. Construct an isosurface for the reduction in acoustic wave velocity. The calculated longitudinal wave velocity values... The baseline wave velocity of the intact rock block in the target cavern Compare and calculate the percentage decrease in wave velocity for each grid node: Each node The value is used to generate a three-dimensional isosurface field of wave velocity reduction amplitude through spatial interpolation.

[0051] 4. Extract the depth-displacement gradient feature curve. For the displacement data extracted along the borehole, calculate the displacement gradient between adjacent measuring points: In the formula, For the first The depth of each measuring point from the orifice. This represents the radial displacement value at the corresponding measuring point. For the first The depth of each displacement measuring point from the borehole opening. For depth The radial displacement measurement value at the location is used. The displacement gradient values ​​of the entire hole depth are used to form a depth-displacement gradient characteristic curve, and the depth at the point where the displacement gradient increases sharply is extracted as an auxiliary criterion for the boundary of the fracture zone.

[0052] 5. Define and calculate the fracture severity index. The fracture severity index... This is a dimensionless parameter ranging from 0 to 1, which integrates three physical quantities: fracture density, fracture slip energy release rate, and wave velocity reduction. In the formula, The average fracture density within the fracture zone. This represents the maximum crack density under all standard operating conditions. This represents the cumulative cracking and slip energy release rate within the fracture zone. This represents the maximum energy release rate across all standard operating conditions. This represents the maximum wave velocity reduction within the rupture zone. The preset wave velocity reduction limit value is taken as the wave velocity reduction value corresponding to the complete fracturing of the rock mass. , , Preset weighting coefficients and satisfying . Corresponding to the undamaged elastic state, This corresponds to a state of severe fragmentation or rockburst.

[0053] 6. Using the microseismic damage density field, wave velocity reduction isosurface, and displacement gradient characteristic curve extracted in the above steps as input features, and combining the true value of rupture depth and rupture degree index... As a label, it forms a training sample.

[0054] S105. Summarize and form an offline digital base map database. Summarize the labeled samples of all standard working conditions to form the offline digital base map, and store it in the database for use by the proxy model for offline pre-training.

[0055] Step S200: Construct and train a physical information-guided hybrid agent model This step represents a core innovation at the algorithm level, using architectural design and physical constraint embedding to guide model training and inference with prior knowledge of rock mechanics.

[0056] S201. Construct the model architecture. The physical information-guided hybrid proxy model includes an input layer, a convolutional coding layer, a spatial alignment layer, a splitting layer, an attention-guided layer, a depth prediction branch, a degree prediction branch, a parallel output layer, and a physical consistency verification layer.

[0057] The convolutional coding layer performs 3D and 2D convolution and global average pooling operations on the input 3D microseismic damage density field and wave velocity reduction amplitude isosurface field, respectively; and performs 1D convolutional coding on the depth-displacement gradient feature curve. Subsequently, through the spatial alignment layer, the extracted multi-source heterogeneous spatiotemporal tensors of varying scales are uniformly flattened and mapped into 1D latent space feature vectors of uniform length, thereby completing the dimensional alignment of multi-source data before entering the routing allocation.

[0058] The demultiplexing layer follows the spatial alignment layer and performs feature allocation based on a physical prior guidance strategy. To avoid the uninterpretability caused by purely data-driven approaches, the demultiplexing layer employs a hard-coded tensor slicing mechanism based on physical rules, rather than a learnable routing gating. Specifically, the demultiplexing layer slices the concatenated global feature tensor: feature slices extracted from the wave velocity reduction isosurface field and displacement gradient are passed to the depth prediction branch, while feature slices extracted from the microseismic damage density field are passed to the degree prediction branch.

[0059] Before the feature slices formally enter the fully connected networks of each prediction branch, an attention guidance layer is set up. This embodiment preferably employs a channel attention mechanism to adaptively learn and strengthen the physical feature dimensions most critical to the current mechanical fracture state within each prediction branch. The initial feature vector received by any branch... For example, the formula for calculating the attention weights of this attention guidance layer is: In the formula, Here is the channel attention weight matrix. It is a non-linear activation function. The Sigmoid activation function is used to map the weights to the 0-1 range. and These are the learnable weight matrices used for dimensionality reduction and dimensionality increase, respectively. Then, the attention weights... With the initial feature vector Element-wise multiplication is performed to obtain the enhanced fused feature vector, which is then input into the subsequent prediction network.

[0060] The depth prediction branch consists of a fully connected layer and a probability output layer. The probability output layer outputs the mean of the predicted rupture depth through two neurons. and standard deviation The probability distribution constituting the fracture depth The degree prediction branch consists of a fully connected layer and a sigmoid-activated output layer, which outputs a predicted fracture degree index. .

[0061] S202. Design the loss function. The training loss function for the hybrid agent model. Loss predicted by depth Degree of loss prediction And physical consistency constraint loss Three-weighted composition: In the formula: , , The preset loss weight coefficient.

[0062] Deep prediction loss Using negative log-likelihood loss: In the formula, This represents the number of samples in the current training batch. No. i The true fracture depth values ​​for each training sample are derived from the offline simulation results of this working condition in the digital base map. For the model to the first i The mean of the Gaussian distribution of the predicted rupture depth for each sample is output by the first neuron of the depth prediction branch probability output layer. For the model to the first i The standard deviation of the Gaussian distribution of the predicted rupture depth for each sample is output by the second neuron of the depth prediction branch probability output layer.

[0063] Degree of prediction loss Mean squared error loss is used: In the formula, This represents the number of samples in the current training batch. For the model to the first i The fracture severity index predicted for each sample is output by the Sigmoid activation layer of the severity prediction branch, and its value ranges from [value range missing]. , No. i The true value of the fracture severity index for each training sample is calculated using the formula described in step S104, and its value range is [range missing]. .

[0064] Physical consistency constraint loss It contains two sub-items. The first sub-item is used to penalize the deviation between the predicted mean rupture depth and the maximum drop depth indicated by the wave velocity drop isosurface: In the formula, This represents the number of samples in the current training batch. For the model to the first i The mean predicted fracture depth for each sample. For the first The depth location where the wave velocity decreases to its maximum value in each sample.

[0065] The second sub-item is used to penalize the deviation between the predicted rupture degree and the estimated spatial concentration of microseismic energy release: In the formula, This represents the number of samples in the current training batch. For the model to the first i The fracture severity index predicted for each sample is output by the Sigmoid activation layer of the severity prediction branch, and its value ranges from [value range missing]. , This represents the accumulated energy in the region with the highest density within the microseismic damage density field. To accumulate energy for the entire region. A sliding window search method based on spatial grids can be used to divide the three-dimensional space where the three-dimensional microseismic damage density field is located into several search voxels of preset size (e.g., a cubic detection unit with a side length of 1m×1m×1m). For each search voxel, the normalized release energy of all microseismic events contained within the geometric boundary of the voxel is calculated by integration.

[0066] During training, optimization is performed using backpropagation and mini-batch gradient descent. Update all network parameters of the model.

[0067] S203. Pre-training and Online Fine-tuning. The hybrid proxy model is fully pre-trained using the offline digital base map database constructed in step S100. During the initial stage of construction of the target cavern, approximately twenty to thirty sets of high-quality monitoring data and corresponding hard data observations are collected for the first excavation cycle of the project. The parameters of the pre-trained model are then fine-tuned to eliminate systematic biases that exist when migrating from the numerical simulation environment to the actual geological environment.

[0068] Step S300: Acquire and process real-time monitoring data, and construct a field monitoring data acquisition and processing link that fully corresponds to the physical properties of the "physical fingerprint field" in the digital base map, ensuring the consistency of the proxy model input.

[0069] S301. Real-time generation of the microseismic damage density field. During the excavation of the target chamber, microseismic waveforms are acquired in real time using an array of microseismic detectors spaced longitudinally along the tunnel. The spatiotemporal coordinates and moment magnitude of each microseismic event are retrieved using a source localization algorithm. A 24-hour rolling time window is set, and for all microseismic events within this window, a three-dimensional microseismic damage density field is generated in real time using the same method as the kernel density estimation method described in step S104.

[0070] S302. Real-time generation of wave velocity reduction isosurface. In anchor bolt boreholes or dedicated acoustic boreholes, ultrasonic waves are emitted and received point-by-point along the borehole depth direction using an acoustic transducer at preset step intervals. The longitudinal wave travel time at each depth is measured and converted into wave velocity values. The measured wave velocity is compared with the reference wave velocity calibrated on the intact rock block before construction. The wave velocity reduction percentage at each measuring point is calculated, and spatial interpolation is used to generate the wave velocity reduction isosurface field at the current moment.

[0071] S303. Real-time generation of the depth-displacement gradient characteristic curve. Displacement data at different depths along the borehole are acquired in real time using multiple displacement gauges deployed along a typical cross-section of the tunnel. The displacement gradient values ​​between adjacent measuring points are calculated using the gradient calculation formula described in step S104, forming the depth-displacement gradient characteristic curve for the current moment.

[0072] Step S400: Model Prediction Output In each excavation cycle, the three-dimensional microseismic damage density field, acoustic wave velocity reduction isosurface, and multi-point displacement gradient feature data obtained in step S300, along with the corresponding working parameters (determined by the geological conditions and construction parameters of the current excavation cycle), are input into the hybrid surrogate model trained in step S200. After receiving the input data, the model performs convolutional encoding, flow guidance, and branch prediction, and outputs two results in parallel: a rupture depth probability distribution (given in Gaussian distribution form as the predicted mean and variance) and a rupture degree index (a dimensionless value between 0 and 1).

[0073] Step S500: Obtain hard data observations S501. Borehole Camera Data Extraction. In each excavation cycle, a series of images of the surrounding rock wall, continuously unfolding along the borehole depth, are acquired using a borehole panoramic camera in the anchor holes or dedicated observation holes. Fractures are identified frame by frame in the image sequence, and the fracture linear density distribution along the borehole depth is statistically analyzed. The measured fracture depth is defined as the depth at which the fracture linear density drops sharply from a high value area at the borehole opening to below a preset background threshold. .

[0074] S502. Data Extraction from Anchor Bolt Drilling Rig. During anchor bolt drilling in each excavation cycle, the drilling speed, feed pressure, and rotational torque of the anchor bolt rig are continuously recorded. When the drilling speed experiences a step increase at a certain depth accompanied by a synchronous decrease in feed pressure, it is determined that the drill bit has passed through the fracture zone and entered the intact rock mass. The depth of this abrupt change point is taken as the measured fracture depth. .

[0075] S503, Hard Data Quality Assessment. For each acquired... Label observation error , which serves as the standard deviation parameter of the likelihood function in subsequent Bayesian updates.

[0076] Step S600: Bayesian Dynamic Update S601. Construct the prior distribution. This is based on the probability distribution of fracture depth output by the hybrid surrogate model in the current tunneling cycle step S400. As a prior distribution: In the formula, Let the rupture depth be a random variable. The mean of the prior distribution is the mean of the fracture depth predicted by the hybrid surrogate model in step S300 for the current excavation cycle. The prior distribution standard deviation is the predicted standard deviation of the fracture depth output by the hybrid surrogate model in step S300 for the current excavation cycle.

[0077] S602. Construct the likelihood function. Using the hard data observations obtained in step S500... and its observation error Constructing a Gaussian likelihood function: In the formula, These are hard data observations, specifically the measured fracture depth extracted from the drilling camera or anchor bolt trolley parameters in step S400. The standard deviation of the hard data observation error is preset based on the resolution of the borehole camera image or the accuracy of the drilling parameter measurement, and is used to characterize the reliability of the hard data observation values.

[0078] S603. Calculate the posterior distribution. According to Bayes' theorem, the posterior probability distribution is proportional to the product of the prior distribution and the likelihood function: Since both the prior and likelihood are Gaussian distributed, the posterior is also Gaussian distributed, and the updated posterior mean can be analytically solved. and posterior standard deviation : S604. Synchronously update the fracture severity index. The updated posterior distribution information of fracture depth is passed to the fracture severity prediction branch. Through the physical consistency verification layer within the hybrid proxy model, the fracture severity index is synchronously adjusted using the physical correlation between depth and severity. The estimated value was made consistent with the updated understanding of the fracture depth.

[0079] S605, Recursive Iteration. The posterior probability distribution parameters obtained in the current loop are used as the prior distribution parameters for the next construction loop, i.e., let... and As the number of tunneling cycles increases, the posterior standard deviation... It shows a gradual convergence trend.

[0080] Step S700: 3D Visualization Dynamic Display The microseismic damage density field obtained in step S300, the probability distribution of fracture depth and fracture degree index output in step S400, and the posterior probability distribution updated in step S600 are pushed to a three-dimensional visualization platform to intuitively present the evolution of the surrounding rock fracture state as the tunneling cycle progresses in the form of dynamic wrapping surfaces and time-series curves.

[0081] S701. Construct a three-dimensional geological model of the target cavern. Based on the designed excavation outline and geological survey data of the target cavern, establish a three-dimensional geological model that includes stratigraphic interfaces, major structural planes, and high-stress zone boundaries. The three-dimensional geological model is discretized using triangular facet meshes or voxel meshes.

[0082] S702. Dynamic rendering of the microseismic damage density field. The real-time three-dimensional microseismic damage density field generated in step S301 is mapped onto the three-dimensional geological model. Using the density values ​​of each volume element or node as input for color mapping, a gradient color band from blue to red is used for rendering, where blue corresponds to low damage density and red corresponds to high damage density. The microseismic damage density field is refreshed as the 24-hour rolling time window progresses, creating a dynamic evolution effect.

[0083] S703. Generation and display of the fracture depth exceedance probability isosurface. Based on the updated posterior probability distribution of fracture depth in step S600. For the preset anchor bolt design length Calculate the probability that the fracture depth exceeds the anchor bolt length: In the formula, This is the cumulative distribution function of the standard normal distribution.

[0084] The exceedance probability value is spatially interpolated along the perimeter of the cavern to generate a probability isosurface, or "risk envelope surface," where the fracture depth exceeds the anchor bolt length. In the 3D visualization interface, the envelope surface area where the exceedance probability exceeds a preset probability threshold is highlighted in a prominent color, indicating a risk of inadequate support design. The preset probability threshold is set to 0.2.

[0085] S704. Spatial distribution coloring of the fracture severity index. The fracture severity index output in step S400 is then colored. Or after step S604 synchronous update Spatial interpolation is performed along the perimeter of the cavern to color and render the cavern wall surface using a preset color grading scheme. For example, Green indicates minor damage. Yellow indicates moderate damage. Severe damage is indicated in red.

[0086] S705, Timeline Linked Display. The dynamic evolution of the following timeline curves is displayed synchronously on the timeline panel of the 3D visualization interface: The temporal evolution curves of multidimensional characteristic parameters, including the changes in microseismic event rate, mean wave velocity reduction, and displacement change rate with the construction calendar; The temporal changes in the posterior mean and confidence interval of the fracture depth visually demonstrate the convergence process of cognitive uncertainty as the tunneling cycle progresses. The current position line of each grade threshold is superimposed and compared with the feature parameter curve. When the parameter curve crosses the threshold line, the curve segment is automatically highlighted.

[0087] S706. Warning Information Overlay and Report Push. When the probability of the fracture depth exceeding the anchor bolt length exceeds a preset probability threshold, or the fracture severity index exceeds a preset safety range, the system automatically overlays a flashing warning marker on the corresponding area of ​​the 3D model and pushes the warning information to the on-site broadcast system, the monitoring center's large screen, and mobile terminals. The warning information includes: the location of the risk area, the specific value of the probability of the fracture depth exceeding the threshold, the current level of the fracture severity index, and the recommended anchor bolt lengthening value or support upgrade plan.

[0088] Example 2 As a specific implementation method, to verify the technical effect and practical value of the above-mentioned technical solution in actual deep high-stress engineering scenarios, this embodiment applies the method proposed in this invention to a typical deep-buried hard rock tunnel project. The tunnel is approximately 800m deep, with a horseshoe-shaped cross-section and a span of approximately 12m. It mainly traverses granite strata, with high ground stress levels and a moderate to strong risk of rockburst. Construction employs a drill-and-blast method for layered excavation, with each excavation cycle advancing approximately 2.5m. Following steps S100 to S700 detailed in Embodiment 1, this embodiment fully executes the entire process of offline digital base map construction, surrogate model training, real-time on-site monitoring and prediction, hard data extraction, Bayesian dynamic updating, and 3D visualization in this project.

[0089] In the offline preparation phase before construction, following step S100, this embodiment extracts three main control parameters—surrounding rock quality index, strength-stress ratio, and excavation-support ratio—based on the geological survey report of the project. The surrounding rock quality index is assessed using the Basic Quality Index (BQ) grading method, the strength-stress ratio is calculated based on the ground stress test results, and the excavation-support ratio is determined according to the construction organization design. By performing a three-factor, multi-level orthogonal combination of these three parameters, a total of 27 standard working conditions are generated. For each working condition, an excavation unloading and fracture evolution simulation is performed using a particle discrete element and finite difference coupled model. After the simulation, the true value of the surrounding rock fracture depth, the three-dimensional microseismic damage density field, the isosurface of the acoustic wave velocity reduction amplitude, and multi-point displacement gradient characteristic data corresponding to each working condition are extracted from the results, forming an offline digital base map. During the construction process… Figure 2 The image shows some key results of the offline digital basemap: Figure 2 (a) shows the three-dimensional response surface of fracture depth under different combinations of surrounding rock quality grades and strength-stress ratios. It can be seen that the fracture depth increases monotonically with the deterioration of surrounding rock quality and the increase of strength-stress ratio, which is in line with the basic laws of rock mechanics. Figure 2 (b) is a slice of the three-dimensional microseismic damage density field on the horizontal plane of the cavern. The high damage density area is concentrated within 2-3m around the cavern wall. Figure 2 In the middle (c), the wave velocity reduction isosurface is shown. The wave velocity reduction gradually decreases along the radial direction and is asymmetrically distributed on both sides of the cavity, reflecting the influence of non-uniform stress redistribution. Figure 2 In the middle (d), the depth-displacement gradient characteristic curve shows obvious gradient abrupt change points near the fracture depth, which can be used to help determine the fracture boundary.

[0090] Based on the aforementioned offline digital base map, and following step S200, this embodiment constructs a physically-guided hybrid proxy model using a multi-branch deep neural network as its basic architecture. This model incorporates a split layer in the input mapping path, assigning the acoustic wave velocity reduction isosurface and depth-displacement gradient features to the depth prediction branch, and the three-dimensional microseismic damage density field to the severity prediction branch. During model training, in addition to the depth prediction loss and severity prediction loss, the loss function also introduces a physical consistency constraint term to penalize the deviation between the mean predicted rupture depth and the location of the maximum gradient of wave velocity reduction, as well as the deviation between the predicted rupture severity index and the spatial concentration of microseismic energy release. Figure 3 As shown, after sufficient pre-training using massive labeled samples from offline digital base maps, the convergence curve of the loss function during model training is as follows: Figure 3 As shown in (a), the model converged after approximately 150 steps, without significant overfitting. To verify the model's prediction accuracy, approximately 25 sets of monitoring data and corresponding hard data observations were collected during the first excavation cycle at the beginning of construction to fine-tune the model. The performance of the fine-tuned model on the validation set is shown below. Figure 3 (b) and Figure 3 As shown in (c), the coefficient of determination R² between the predicted fracture depth and the actual value in the numerical simulation reaches 0.96, and the prediction error is mainly concentrated within ±0.3m; the R² between the predicted fracture degree index and the actual value reaches 0.94, indicating that the hybrid surrogate model can predict fracture depth and fracture degree simultaneously with high accuracy, and the calculation time during online inference is reduced from several days in traditional numerical simulation to milliseconds.

[0091] After entering the construction phase, according to step S300, this embodiment implemented real-time monitoring data acquisition and processing during the excavation process. A microseismic detector array was deployed at intervals along the tunnel's longitudinal direction to acquire microseismic waveforms in real time. The spatiotemporal coordinates and moment magnitude of each microseismic event were retrieved using a source location algorithm, and a three-dimensional microseismic damage density field was generated using a 24-hour rolling time window. Simultaneously, longitudinal wave travel time was tested point-by-point using acoustic transducers within the anchor bolt boreholes, calculating the wave velocity reduction percentage at each measuring point and generating a wave velocity reduction isosurface field. Furthermore, depth-displacement gradient characteristic curves were acquired in real time using multi-point displacement gauges deployed at typical sections of the tunnel. In each excavation cycle, according to step S400, the above three types of real-time monitoring characteristic data and corresponding working condition parameters were input into the trained hybrid surrogate model. The model outputs in parallel the probability distribution of rupture depth and the rupture severity index for the current state of the cycle. For example, in a typical cycle, the model output a rupture depth prediction mean of 2.8m, a prediction standard deviation of 0.6m, and a rupture severity index of 0.65.

[0092] After each excavation cycle is completed, according to step S500, this embodiment extracts hard data observations in two ways. One is the borehole camera method: the panoramic camera is lowered to the bottom of the anchor bolt hole and raised at a constant speed to acquire a sequence of borehole wall images that unfold continuously along the hole depth, and the fracture line density curve is statistically analyzed, such as... Figure 7 As shown in (a), a sudden drop in density occurs at 2.1m, and this depth is the measured fracture depth. Another method is the anchor bolt trolley-based drilling parameter method: continuously recording the drilling speed curve, such as... Figure 7 As shown in (b), the drilling speed increases dramatically at 2.05m, accurately identifying the fracture depth. The two methods serve as backups for each other, ensuring the reliability of the hard data. In this embodiment, the measured fracture depth extracted by this loop is 3.2m, with an observation error of 0.2m.

[0093] Subsequently, following step S600, this embodiment uses the probability distribution of fracture depth output in step S400 (mean 2.8m, standard deviation 0.6m) as the prior, and uses the Gaussian likelihood function constructed from the measured fracture depth (3.2m) and its observation error (0.2m) as the basis to calculate the posterior distribution using Bayes' theorem. Figure 4 As shown in Figure (a), the updated posterior distribution has a mean of 3.06m and a standard deviation that narrows to 0.19m, significantly reducing cognitive uncertainty. Simultaneously, through the physical consistency verification layer within the model, the updated perception of rupture depth is synchronously transmitted to the severity prediction branch, adjusting the estimated value of the rupture severity index. The posterior probability distribution parameters of the current cycle are used as the prior for the next cycle, forming a recursive iteration. Figure 4 As shown in Figure (b), with the continuous accumulation of hard data in 20 tunneling cycles, the posterior standard deviation gradually decreased from the initial 0.8m to about 0.15m, and the posterior mean gradually approached the actual fracture depth (3.0m), verifying the effectiveness of the dynamic iterative evaluation of the present invention.

[0094] Finally, following step S700, this embodiment pushes the microseismic damage density field, rupture depth probability distribution, rupture degree index, and updated posterior probability distribution to the three-dimensional visualization platform. For example... Figure 5 (a) and Figure 5 As shown in (b), a probability isosurface, or "risk envelope surface," is generated on the 3D model of the cavern where the fracture depth exceeds the designed anchor bolt length (2.5m). The red area (P>0.2) is located at the left arch waist and right wall foot of the cavern, indicating that the anchor bolt length in this area is insufficient and requires reinforcement. (b) is a spatial distribution map of the fracture degree index along the cavern wall. The index in a local area near the arch reaches 0.7~0.8 (severe fragmentation), which is consistent with the location of the high-density microseismic damage area, providing a visual basis for spatially differentiated adjustment of support parameters. Figure 6The system displays a dynamic evolution curve over time, including the microseismic event rate, the mean wave velocity decrease, and the posterior mean fracture depth along with its 95% confidence interval. As the excavation face advances, the number of microseismic events initially increases and then decreases, the wave velocity decrease gradually increases, and the posterior mean fracture depth gradually increases from 2.4m to 2.9m, while the confidence interval narrows, reflecting that the understanding of the surrounding rock fracture state is gradually approaching reality. When the probability of the fracture depth exceeding the anchor bolt length in a certain cycle exceeds 0.2, the system automatically issues an early warning.

[0095] The aforementioned engineering implementation results fully reproduce the entire process of this invention, from offline digital base map construction, real-time monitoring feature extraction, hybrid surrogate model prediction, Bayesian dynamic updating to 3D visualization. This demonstrates that the method can achieve real-time, probabilistic, and adaptive dynamic evaluation of the fracture depth and degree of hard rock caverns, and that the evaluation results gradually converge to the true values ​​with each tunneling cycle. Compared with traditional numerical simulation and post-detection methods, this invention significantly improves the timeliness of the evaluation and the richness of decision-making information.

[0096] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A dynamic evaluation method for the fracture state of hard rock caverns optimized during tunneling, characterized in that, Includes the following steps: Step S100: Based on the geological exploration parameters of the target cavern, establish a working condition-fracture response database generated in batches by the micro-macro coupled numerical model, extract the working condition parameters, true values ​​of surrounding rock fracture depth, three-dimensional microseismic damage density field, isosurface of acoustic wave velocity reduction amplitude and multi-point displacement gradient characteristic data corresponding to each working condition in the database, and form an offline digital base map. Step S200: Based on a multi-branch deep neural network architecture, construct a hybrid proxy model guided by physical information, and embed a physical prior guidance mechanism into the model architecture; use the three-dimensional microseismic damage density field, the acoustic wave velocity reduction isosurface, the multi-point displacement gradient feature data and the corresponding working condition parameters in the offline digital base map as training inputs, and use the true value of the surrounding rock fracture depth and the predefined fracture degree index as training outputs to train the hybrid proxy model; Step S300: Acquire the microseismic monitoring signal, acoustic test signal and multi-point displacement meter signal collected in real time during the excavation of the target cavern; process the microseismic monitoring signal into a three-dimensional microseismic damage density field; process the acoustic test signal into an isosurface of acoustic wave velocity reduction relative to the reference wave velocity; and process the multi-point displacement meter signal into multi-point displacement gradient feature data. Step S400: In each excavation cycle of the target cavern, the three-dimensional microseismic damage density field, the acoustic wave velocity reduction isosurface, and the multi-point displacement gradient feature data obtained in step S300, as well as the corresponding working condition parameters, are input into the trained hybrid proxy model. The hybrid proxy model outputs in parallel the probability distribution of the rupture depth and the rupture degree index of the current state of the target cavern. Step S500: After each excavation cycle of the target cavern is completed, the borehole camera image data or the anchor trolley drilling parameter data corresponding to the cycle are obtained, and the measured fracture depth is extracted as the hard data observation value. Step S600: The probability distribution of fracture depth output in step S400 is used as the prior probability distribution, and the hard data observations obtained in step S500 are used as the likelihood basis. Bayesian inference is used to update the probability distribution of fracture depth and fracture degree index of the target cavern in the current state. The posterior probability distribution is then used as the prior probability distribution of step S400 in the next excavation cycle, thereby realizing the evaluation of fracture state through dynamic iteration with the tunneling cycle.

2. The method according to claim 1, characterized in that, The micro-macro coupled numerical model in step S100 is a particle discrete element and finite difference coupled model. The offline digital base map generation process includes: extracting three parameters—surrounding rock quality index, strength-stress ratio, and excavation-support ratio—based on the geological survey parameters of the target cavern; performing hierarchical orthogonal combination of the three parameters to form a standard working condition set; for each standard working condition in the standard working condition set, performing offline simulation of excavation unloading and fracture evolution using the particle discrete element and finite difference coupled model; extracting the true value of the surrounding rock fracture depth corresponding to the working condition from the simulation results; performing three-dimensional spatial meshing processing; and generating an offline digital base map composed of the true value of the surrounding rock fracture depth, the three-dimensional microseismic damage density field, the isosurface of the acoustic wave velocity reduction amplitude, and the multi-point displacement gradient feature data.

3. The method according to claim 1, characterized in that, In step S200, the hybrid proxy model internally constructs a depth prediction branch and a degree prediction branch; the physical prior guidance mechanism is configured to: cause the acoustic wave velocity reduction isosurface and the multi-point displacement gradient feature data related to depth prediction to flow to the depth prediction branch, and cause the three-dimensional microseismic damage density field related to degree prediction to flow to the degree prediction branch.

4. The method according to claim 3, characterized in that, The physical prior guidance mechanism includes setting a diversion layer and an attention guidance layer in the input mapping path of the hybrid proxy model; the diversion layer is used to allocate the acoustic wave velocity reduction isosurface and the multi-point displacement gradient feature data to the depth prediction branch, and to allocate the three-dimensional microseismic damage density field to the degree prediction branch; the attention guidance layer is used to learn and strengthen the correlation weights between different feature dimensions within each prediction branch and the training output.

5. The method according to claim 1, characterized in that, During training, the hybrid proxy model employs a loss function comprising a depth prediction loss term, a severity prediction loss term, and a physical consistency constraint term. The physical consistency constraint term penalizes the deviation between the mean predicted rupture depth and the isosurface of the acoustic wave velocity reduction, as well as the deviation between the predicted rupture severity index and the three-dimensional microseismic damage density field. During training, the network parameters of the hybrid proxy model are updated by optimizing the weighted sum of the depth prediction loss term, the severity prediction loss term, and the physical consistency constraint term.

6. The method according to claim 1, characterized in that, In step S400, the fracture depth probability distribution and the fracture degree index are determined in the following manner: The probability distribution of the rupture depth is obtained based on the probability density function constructed from the predicted mean and predicted variance output by the hybrid surrogate model; The fracture degree index is a dimensionless index established by weighting multiple physical quantities such as fracture density, crack slip energy release rate, and longitudinal wave velocity reduction, and is used to characterize the damage and fracture evolution state of the surrounding rock at different stages.

7. The method according to claim 1, characterized in that, The extraction method of the hard data observations in step S500 is as follows: for the borehole camera image data, the measured fracture depth is determined by identifying the abrupt change point of fracture density or the location of the fracture surface along the borehole depth direction; for the anchor trolley drilling parameter data, the measured fracture depth is determined by identifying the abrupt change point of drilling speed at a certain borehole depth.

8. The method according to claim 1, characterized in that, In step S600, the Bayesian inference employs a sequential Bayesian filtering or particle filtering algorithm. After acquiring new hard data observations each time, the posterior probability distribution of the previous cycle is used as the prior probability distribution of the current cycle. The likelihood function is calculated using the hard data observations of the current cycle, and the posterior probability distribution of the current cycle is calculated using the Bayesian formula. This allows the evaluation results of the fracture depth and fracture degree index of the target cavern in its current state to gradually converge to a narrower confidence interval as the tunneling cycle increases.

9. The method according to claim 1, characterized in that, The process also includes step S700: pushing the three-dimensional microseismic damage density field obtained in step S300, the probability distribution of rupture depth and the rupture degree index output in step S400, and the posterior probability distribution updated in step S600 to a three-dimensional visualization platform. The platform displays the spatial distribution of the probability contour lines of rupture depth exceeding the preset anchor length and the rupture degree index in the form of a dynamic wrapping surface on the three-dimensional geological model of the target cavern, and synchronously displays the evolution curves of rupture depth and the rupture degree index and the changes of the confidence interval updated by Bayes on the time axis.

Images