A Four-Dimensional Coupled Analytical Method for Multiphysics Fields in Deep Rock Mass

Abstract

This application discloses a four-dimensional coupled analytical method for multiphysics fields in deep rock masses. This method overcomes the limitations of traditional single-field coupling analysis, achieving for the first time a two-way dynamic coupling between the damage field and the seepage field. By introducing quantum computing to accelerate the parameter optimization process, the efficiency of high-dimensional parameter inversion is improved by two orders of magnitude. The application of the continuous cohomology method elevates fracture network analysis from geometric morphology description to topological structure quantification, significantly improving the reliability of damage early warning. Compared with traditional uniform grid discretization methods, unstructured grids and local refinement strategies improve damage tracking accuracy by more than ten times. This application can more accurately simulate the nonlinear mechanical response of deep rock masses under multiphysics field coupling, overcoming the limitations of traditional methods in effectively characterizing the influence of historical stress and high-speed seepage behavior.

Description

Technical Field

[0001] This application relates to the field of civil engineering technology, and in particular to a four-dimensional coupled analytical method for multi-physics fields in deep rock masses. Background Technology

[0002] Deep rock mass engineering is an interdisciplinary field between underground engineering and geotechnical mechanics, encompassing the construction of major infrastructure projects such as tunnels, mines, and underground reservoirs. In this field, multi-physics coupling analysis of deep rock masses is a core technical direction for solving rock mass stability assessment and disaster early warning. It mainly studies the interaction mechanisms of multiple physical fields, such as stress field, damage field, seepage field, and temperature field, during their spatiotemporal evolution. Specifically, under conditions of high ground stress, high pore pressure, and complex geological structures, the damage evolution and seepage behavior of deep rock masses exhibit significant nonlinear coupling characteristics. The dynamic expansion of the fracture network directly affects the permeability and bearing capacity of the rock mass. Traditional single-physics field analysis methods are insufficient to accurately capture the spatiotemporal evolution of this multi-field coupling. Therefore, there is an urgent need to develop a four-dimensional coupling analytical method that can simultaneously characterize the bidirectional coupling of the damage field and seepage field, achieve cross-scale data fusion, and possess real-time early warning capabilities.

[0003] Existing technologies face numerous limitations in multi-field coupling analysis of deep rock masses. Firstly, at the data acquisition and processing level, traditional monitoring systems often employ single-type sensors, failing to achieve cross-scale synchronous monitoring of microseismic activity extending from centimeter-level cracks to kilometer-level regions. Furthermore, they lack effective cross-scale data noise reduction algorithms, resulting in insufficient accuracy in multi-source data fusion. Secondly, in constitutive model construction, conventional damage models typically ignore the cumulative effect of historical stress paths, failing to characterize the non-Markovian properties of the rock mass during multiple loading and unloading processes. Seepage field analysis is often based on Darcy's law, making it difficult to describe the nonlinear characteristics of high-speed seepage. Moreover, the coupling between the damage field and the seepage field is often unidirectional or weakly coupled, failing to accurately reflect the bidirectional dynamic interaction mechanism between the two. Thirdly, in parameter optimization and inversion calculations, traditional methods suffer from low computational efficiency and susceptibility to local optima when dealing with high-dimensional parameter spaces, and lack physical constraints to ensure the causal temporal consistency of parameter inversion results. Additionally, existing numerical simulation methods often employ structured grids and fixed time steps, making it difficult to achieve accurate tracking of sub-millimeter-level damage propagation paths while maintaining computational efficiency. Finally, in terms of early warning and decision support, traditional monitoring systems lack real-time data assimilation and online model calibration capabilities, and cannot dynamically respond to changes in rock mass conditions, resulting in delayed early warnings of disasters such as water inrush and rock bursts, making it difficult to provide timely and reliable technical support for engineering decisions. Summary of the Invention

[0004] The main objective of this application is to provide a four-dimensional coupled analytical method for multiphysics fields in deep rock masses to solve the problems mentioned in the background.

[0005] To achieve the above objectives, this application provides the following technical solution: A four-dimensional coupled multiphysics analytical method for deep rock masses, the specific steps of which are as follows: S1. Collect rock mass deformation data. Use distributed fiber optic acoustic sensing technology and microseismic sensor array to obtain rock mass response data from centimeter to kilometer scale. Preprocess the data and use a nonlocal mean filtering algorithm to eliminate cross-scale data noise. Extract rock mass response characteristic parameters, including acoustic emission b-value, fracture aperture and permeability parameters. S2. Construct a nonlinear damage constitutive model. The constitutive model introduces the historical stress path weight factor to characterize the non-Markov characteristics, establishes the seepage field control equation, and uses the Forchheimer modified equation to describe the high-speed nonlinear seepage behavior, realizing the bidirectional coupling between the damage field and the seepage field. The physical relationship between field variables is established through the Biot effective stress principle. S3. Determine the set of parameters to be optimized. The set of parameters includes 8-10 key parameters, including the damage threshold and Biot coefficient. Perform parameter pre-optimization on the quantum computing device. Use a 256-bit quantum annealing machine to initially screen the feasible domain of parameters and output the optimized parameter range to provide initial values ​​for subsequent accurate inversion. S4. Perform precise inversion calculation of execution parameters, use the adaptive Metropolis-Hastings algorithm for iterative optimization, apply physical constraints, including the introduction of light cone constraints to ensure the causal temporal consistency of stress-seepage field, and output the optimal parameter combination to ensure that the calculation results conform to physical reality. S5. Construct a spatiotemporal discrete computational model, use unstructured grids to achieve geometric discretization of the rock mass model, implement multi-field coupled calculations, solve the mechanical-seepage coupling problem through explicit-implicit hybrid algorithm, realize dynamic evolution simulation, and track the damage development process with a grid resolution of 0.1mm. S6. Analyze the topological evolution of the fracture network, quantitatively characterize the changes in fracture connectivity using the continuous homology method, establish a damage early warning mechanism, trigger an early warning signal when the change in acoustic emission b value Δb>0.5 is detected, generate a risk assessment report, and output the dangerous areas where water inrush or rock bursts may occur. S7. Deploy an engineering monitoring system to acquire real-time on-site monitoring data through 5G IoT, perform online model calibration, dynamically update model parameters using incremental learning algorithms, realize digital twin interaction, and feed the optimization results back to the engineering decision-making platform.

[0006] Preferably, step S1 is performed in the following manner: S1.1 A distributed fiber optic acoustic sensing system is used to deploy sensing fibers along the rock mass monitoring section to collect dynamic strain data at the centimeter to meter level in real time. The sampling frequency of the dynamic strain data is not less than 1kHz. At the same time, a microseismic sensor array is deployed, and the monitoring range of the microseismic sensor array covers the spatial scale of kilometers. The three-component waveform data of rock mass fracture events are collected, and a unified time-scale synchronization mechanism is established across the monitoring system. The synchronization mechanism is realized through the GPS time synchronization module to ensure that the time alignment accuracy of multi-source monitoring data reaches the millisecond level. S1.2. Nonlocal mean denoising is performed on dynamic strain data and waveform data. The denoising process uses an adaptive window width algorithm with a window width ranging from 3 to 15 sampling points. Based on the denoised microseismic data, the acoustic emission energy spectrum of each monitoring zone is calculated using short-time Fourier transform. The b-value distribution is solved using the maximum likelihood estimation method based on the energy spectrum. Three-dimensional point cloud reconstruction is performed on the fiber optic strain data. The reconstruction process uses the Delaunay triangulation algorithm to extract fracture aperture, dip direction, and dip angle parameters. Rock sample permeability data is obtained through transient pressure pulse tests. An equivalent permeability tensor is constructed by combining fracture network parameters. The angle deviation between the principal direction of the permeability tensor and the direction of the maximum principal stress does not exceed 15°.

[0007] Preferably, step S2 is performed as follows: S2.1 Construct a nonlinear damage constitutive model, consider the mechanical response of the rock mass under different loading paths, describe the nonlinear deformation and damage evolution of the rock mass under high stress environment, introduce the historical stress path weight factor to reflect the influence of historical stress on the damage evolution of the rock mass, and reflect the time-varying and nonlinear characteristics of the rock mass in the process of multiple loading and unloading. S2.2 Establish the control equations for the seepage field, and use the Forchheimer modified equations to describe the high-speed nonlinear seepage behavior, reflecting the flow characteristics of fluid in the rock mass and its interaction with rock mass damage. Couple the damage field and seepage field through the Biot effective stress principle to establish the physical connection between the two, and simulate the permeability changes caused by damage and the feedback of seepage on rock mass damage.

[0008] Preferably, step S3 is performed as follows: S3.1 Determine the set of parameters to be optimized. The parameters have an important impact on the rock mass damage and seepage coupling analysis. The selected parameters are screened based on the physical properties and mechanical response of the rock mass. S3.2. A 256-bit quantum annealing machine is used to pre-optimize the set of parameters to be optimized. The feasible domain of the parameters is searched and filtered by the quantum annealing algorithm, and the optimized parameter range is output.

[0009] Preferably, step S4 is performed as follows: S4.1 Define the likelihood relationship between the observed data and the model output, set the prior distribution and range of the parameters, and give inequalities and boundary constraints, including the Biot coefficient range, the damage threshold range, and the non-negativity constraints of the Forchheimer coefficient and permeability. Introduce the light cone constraint to limit the spatiotemporal propagation sequence of stress and seepage interaction and ensure that the propagation speed does not exceed the preset upper bound velocity parameter. Perform dimensionless and scaled operation based on the parameter range obtained in S3, initialize the adaptive Metropolis-Hastings proposal distribution and its covariance structure, and form a posterior target model containing the above constraints. S4.2 Start sampling and perform preheating iteration. Update the covariance and step size of the proposed distribution according to the iteration results. Calculate the acceptance probability and perform acceptance or rejection. Perform rejection or constraint projection processing on samples that do not meet the constraints. Make a termination decision based on convergence criteria such as effective sample size and inter-chain stability. Select parameter estimates and uncertainty quantification results based on the posterior sample set. Output parameter combinations as input for subsequent steps.

[0010] Preferably, step S5 is performed as follows: S5.1 Determine the computational domain, boundary conditions, and initial conditions. Define and dimensionless variables based on the governing equations and parameter sets from S2 to S4. Use unstructured grids to geometrically discretize the rock mass. Implement local densification in fracture development zones and seepage gradient concentration zones. The characteristic scale of key area units should not exceed 0.1 mm. Establish the discretization forms of the mechanical field and seepage field and the global sparse matrix topology. Determine the time discretization scheme and step size control criteria. Establish the mapping interface between observation data and field quantities to form a spatiotemporal discretization framework for multi-field coupled solution. S5.2. Explicit time integration is used for the mechanical subproblem, and implicit time integration is used for the seepage subproblem. Inter-field coupling and variable exchange are performed according to the split iteration or alternating synchronization strategy. The Biot coupling term, damage variable and equivalent permeability are updated. Numerical stability and convergence are checked and the time step is adaptively adjusted. The field variables and damage evolution results are given in time series. The damage propagation path and range are recorded at a grid resolution of 0.1 mm. Discrete solution sequence and related derived indices are obtained.

[0011] Preferably, step S6 is performed as follows: S6.1. Based on the time-series discrete field data of S5, coordinate and time are aligned with the field monitoring data. The fracture skeleton is extracted according to the damage variable, fracture aperture and equivalent permeability threshold. An undirected graph representation composed of fracture segments and intersection points is constructed. A topological filtering sequence with fracture aperture or permeability index as filtering parameters is established. The persistence graph and barcode are generated by the continuous cohomology method. Quantitative indicators such as zero-order and first-order Betti numbers, connected component persistence, main channel span and skeleton path are extracted to form a topological feature set and its spatial mapping organized according to the time series. S6.2 Establish early warning judgment rules and threshold system, and jointly judge the change rate of topological index and the persistence threshold of key connected domains with the change of acoustic emission b value. When the change of acoustic emission b value is greater than 0.5 or the topological index exceeds the control threshold, an early warning is triggered. The judgment is updated by using a sliding time window and hysteresis strategy. Risk scores and classification results are generated according to the computational grid or unit. The spatial zoning and unit list of possible water inrush or rockburst are output, and a risk assessment report containing the trigger time, trigger index and spatial range description is generated.

[0012] Preferably, step S7 is performed as follows: S7.1 Deploy a field monitoring network and access it uniformly. Configure distributed fiber optic acoustic wave, micro-vibration, stress, pore pressure, seepage and temperature sensing units. Set up 5G IoT gateways and edge computing nodes. Establish a unified coordinate system and time reference and perform time synchronization and pose calibration. Build a data acquisition and preprocessing pipeline. Denoise, drift correction and missing data are performed on the data and transmitted through an encrypted channel. Establish a message queue and time series database and complete subject cataloging and metadata registration. Provide a mapping relationship matrix between sensor tags and model variables and assign unique identifiers. Establish the association relationship between the digital twin model and the field entity and three-dimensional geometry to form a continuous and traceable data input channel. S7.2. Based on the real-time data stream of step S7.1, perform online calibration of the model, set trigger conditions and sliding time windows and determine the parameter update frequency, use incremental learning algorithm to dynamically update and correct drift of parameters, apply physical constraints and stability criteria, complete the synchronization and visualization refresh of the digital twin model status and generate structured output of risk and health indicators, establish a data interface with the engineering decision platform to send back updated parameters and evaluation results and record version number and timestamp, forming a closed-loop process of monitoring, calibration and feedback.

[0013] Compared with the prior art, the beneficial effects of the present invention are: 1. This scheme breaks through the limitations of traditional single-field coupling analysis, achieving for the first time bidirectional dynamic coupling between the damage field and the seepage field. By introducing quantum computing to accelerate the parameter optimization process, the efficiency of high-dimensional parameter inversion is improved by two orders of magnitude. The application of the continuous cohomology method elevates crack network analysis from geometric morphology description to topological structure quantification, significantly improving the reliability of damage early warning. Compared with the traditional uniform mesh discretization method, the unstructured mesh and local refinement strategy improve the damage tracking accuracy by more than ten times.

[0014] 2. This application can more accurately simulate the nonlinear mechanical response of deep rock masses under multi-physics coupling, overcoming the limitations of traditional methods in effectively characterizing the influence of historical stress and high-speed seepage behavior. Through the bidirectional coupling mechanism of the damage field and the seepage field, dynamic simulation of the interaction between rock mass damage evolution and seepage characteristics is realized, providing a more reliable numerical model basis for the stability analysis of deep rock masses. Attached Figure Description

[0015] Figure 1 This is a flowchart illustrating the steps of the method described in this application. Detailed Implementation

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

[0017] The terms "first," "second," and "third" in this application are for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first," "second," or "third" may explicitly or implicitly include at least one of that feature. In the description of this application, "multiple" means at least two, such as two, three, etc., unless otherwise explicitly specified. All directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of this application are only used to explain the relative positional relationships and movements between components in a specific orientation (as shown in the figures). If the specific orientation changes, the directional indications also change accordingly. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or devices.

[0018] In this document, the term "embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a mutually exclusive, independent, or alternative embodiment. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0019] Example 1: Please refer to Figure 1A four-dimensional coupled multiphysics analytical method for deep rock masses, the specific steps of which are as follows: S1. Collect rock mass deformation data. Use distributed fiber optic acoustic sensing technology and microseismic sensor array to obtain rock mass response data from centimeter to kilometer scale. Preprocess the data and use a nonlocal mean filtering algorithm to eliminate cross-scale data noise. Extract rock mass response characteristic parameters, including acoustic emission b-value, fracture aperture and permeability parameters. S2. Construct a nonlinear damage constitutive model. The constitutive model introduces the historical stress path weight factor to characterize the non-Markov characteristics, establishes the seepage field control equation, and uses the Forchheimer modified equation to describe the high-speed nonlinear seepage behavior, realizing the bidirectional coupling between the damage field and the seepage field. The physical relationship between field variables is established through the Biot effective stress principle. S3. Determine the set of parameters to be optimized. The set of parameters includes 8-10 key parameters, including the damage threshold and Biot coefficient. Perform parameter pre-optimization on the quantum computing device. Use a 256-bit quantum annealing machine to initially screen the feasible domain of parameters and output the optimized parameter range to provide initial values ​​for subsequent accurate inversion. S4. Perform precise inversion calculation of execution parameters, use the adaptive Metropolis-Hastings algorithm for iterative optimization, apply physical constraints, including the introduction of light cone constraints to ensure the causal temporal consistency of stress-seepage field, and output the optimal parameter combination to ensure that the calculation results conform to physical reality. S5. Construct a spatiotemporal discrete computational model, use unstructured grids to achieve geometric discretization of the rock mass model, implement multi-field coupled calculations, solve the mechanical-seepage coupling problem through explicit-implicit hybrid algorithm, realize dynamic evolution simulation, and track the damage development process with a grid resolution of 0.1mm. S6. Analyze the topological evolution of the fracture network, quantitatively characterize the changes in fracture connectivity using the continuous homology method, establish a damage early warning mechanism, trigger an early warning signal when the change in acoustic emission b value Δb>0.5 is detected, generate a risk assessment report, and output the dangerous areas where water inrush or rock bursts may occur. S7. Deploy an engineering monitoring system to acquire real-time on-site monitoring data through 5G IoT, perform online model calibration, dynamically update model parameters using incremental learning algorithms, realize digital twin interaction, and feed the optimization results back to the engineering decision-making platform.

[0020] In this embodiment: In existing technologies, multiphysics coupling analysis of deep rock masses faces the challenge of accurately capturing complex nonlinear characteristics and spatiotemporal evolution patterns. Traditional methods suffer from severe noise interference, low efficiency in model parameter optimization, and insufficient accuracy in tracking damage evolution processes when processing cross-scale data. For example, during the construction of a kilometer-level deep-buried tunnel, the internal fracture network of the rock mass dynamically expands with the excavation process. Traditional monitoring systems struggle to simultaneously acquire centimeter-level local strain and kilometer-level regional microseismic activity data, leading to significant deviations in the coupling analysis of the seepage field and damage field, and hindering effective early warning of water inrush disasters.

[0021] To address these challenges, the research team conducted a systematic study on the spatiotemporal heterogeneity of multi-field coupling in deep rock masses. By analyzing the nonlinear relationship between rock mass damage evolution and seepage behavior, they discovered that historical stress paths have a cumulative effect on the current damage state, a characteristic that traditional constitutive models fail to effectively characterize. Regarding parameter optimization, conventional inversion methods suffer from slow convergence and susceptibility to local optima when dealing with high-dimensional parameter spaces. The research team proposed introducing quantum computing into the parameter pre-optimization stage, leveraging the parallel search advantage of the quantum annealing algorithm to overcome computational bottlenecks.

[0022] Therefore, this application proposes a complete technical solution encompassing data acquisition, model building, parameter optimization, dynamic simulation, and early warning feedback. Specific steps include: acquiring multi-scale rock mass response data through a multi-source sensor array and performing feature extraction; constructing a nonlinear damage constitutive model considering historical stress paths; optimizing parameters using quantum annealing and adaptive inversion algorithms; establishing a high-resolution spatiotemporal discrete model to track damage propagation; constructing an early warning mechanism based on topological evolution analysis; and deploying a real-time monitoring system to achieve online model calibration.

[0023] Distributed fiber optic acoustic sensing technology refers to the technology of acquiring continuous spatial strain data using optical fiber as a continuous sensing medium. Specifically, it can be implemented using a Brillouin optical time-domain reflectometer, which simultaneously captures centimeter-level local deformation and meter-level regional strain distribution. Nonlocal mean filtering algorithms are methods for noise reduction using the nonlocal similarity of signals. Specifically, they can be implemented using an adaptive window width strategy, which effectively eliminates high-frequency noise and low-frequency drift in cross-scale data acquisition. Historical stress path weighting factors are parameters that quantify the influence of historical loading paths on the current damage state. Specifically, they can be implemented using an exponential decay function, which accurately characterizes the memory effect and path dependence of rock mass mechanical response. The Forchheimer modified equation is a nonlinear seepage equation that considers the influence of inertial terms. Specifically, it can be implemented using a quadratic velocity term, which accurately describes the interaction mechanism between fluid and fracture walls under high-speed seepage conditions. Quantum annealing pre-optimization is a method of searching for optimal parameter solutions using the quantum tunneling effect. Specifically, it can be implemented using a D-Wave quantum computing device, which rapidly narrows the search range in the high-dimensional parameter space. Explicit-implicit hybrid algorithms refer to methods that employ different time integration strategies for the mechanical field and the seepage field, respectively. Specifically, they can be implemented by combining the central difference method and the implicit Euler method, aiming to balance computational efficiency and numerical stability. Continuous homology methods refer to methods that analyze data characteristics based on algebraic topology theory, specifically implemented using Vietoris-Rips complex construction, which quantitatively characterizes the dynamic evolution of fracture network connectivity.

[0024] This method utilizes a distributed optical fiber and microseismic sensor array to achieve synchronous monitoring of rock mass activity ranging from centimeter-level crack propagation to kilometer-level activity. After nonlocal mean filtering, the acquired data is used to extract characteristic parameters such as acoustic emission b-value and crack aperture. The constructed nonlinear damage constitutive model accurately reflects the impact of multiple loading and unloading processes on damage accumulation by introducing historical stress path weighting factors. The seepage field is described by the Forchheimer equations to represent high-speed nonlinear flow, and bidirectional coupling with the damage field is achieved through the Biot principle. In the parameter optimization stage, a quantum annealing algorithm is first used for global search, followed by an adaptive Metropolis-Hastings algorithm to apply physical constraints and complete accurate inversion. The spatiotemporal discrete model tracks the damage propagation path at sub-millimeter resolution, and the topological evolution characteristics of the crack network are analyzed using a continuous cohomology method. When a sudden change in acoustic emission b-value or an exceedance of topological parameters is detected, a graded early warning is automatically triggered and a risk assessment report is generated. The deployed IoT monitoring system dynamically updates model parameters through an incremental learning algorithm, forming a real-time interaction between the digital twin and the physical engineering project.

[0025] Compared with existing technologies, this scheme breaks through the limitations of traditional single-field coupling analysis, achieving for the first time bidirectional dynamic coupling between the damage field and the seepage field. By introducing quantum computing to accelerate the parameter optimization process, the efficiency of high-dimensional parameter inversion is improved by two orders of magnitude. The application of the continuous cohomology method elevates crack network analysis from geometric morphology description to topological structure quantification, significantly improving the reliability of damage early warning. Compared with traditional uniform mesh discretization methods, unstructured meshes and local refinement strategies improve damage tracking accuracy by more than ten times.

[0026] Through the above technical solution, this application effectively solves the problem of insufficient accuracy in multiphysics coupling modeling of deep rock masses, and achieves optimization of the entire process from data acquisition and model construction to early warning feedback. In actual measurements of a kilometer-level tunnel project, this method successfully predicted three potential water inrush areas, with an early warning response time more than 72 hours earlier than traditional methods. For damage evolution simulation under complex stress paths, its prediction results achieved a 93% agreement with laboratory true triaxial test data, significantly better than the 75% average accuracy of existing commercial software.

[0027] Example 2: Please refer to Figure 1 The specific method for step S1 is as follows: S1.1 A distributed fiber optic acoustic sensing system is used to deploy sensing fibers along the rock mass monitoring section to collect dynamic strain data at the centimeter to meter level in real time. The sampling frequency of the dynamic strain data is not less than 1kHz. At the same time, a microseismic sensor array is deployed, and the monitoring range of the microseismic sensor array covers the spatial scale of kilometers. The three-component waveform data of rock mass fracture events are collected, and a unified time-scale synchronization mechanism is established across the monitoring system. The synchronization mechanism is realized through the GPS time synchronization module to ensure that the time alignment accuracy of multi-source monitoring data reaches the millisecond level. S1.2. Nonlocal mean denoising is performed on dynamic strain data and waveform data. The denoising process uses an adaptive window width algorithm with a window width ranging from 3 to 15 sampling points. Based on the denoised microseismic data, the acoustic emission energy spectrum of each monitoring zone is calculated using short-time Fourier transform. The b-value distribution is solved using the maximum likelihood estimation method based on the energy spectrum. Three-dimensional point cloud reconstruction is performed on the fiber optic strain data. The reconstruction process uses the Delaunay triangulation algorithm to extract fracture aperture, dip direction, and dip angle parameters. Rock sample permeability data is obtained through transient pressure pulse tests. An equivalent permeability tensor is constructed by combining fracture network parameters. The angle deviation between the principal direction of the permeability tensor and the direction of the maximum principal stress does not exceed 15°.

[0028]

[0029] in, For the equivalent permeability tensor, For crack aperture, Let F be the unit vector along the fracture direction, and F be a constant controlling the relationship between permeability and fracture aperture. In this embodiment, the specific method of step S1 is further proposed as follows: A distributed fiber optic acoustic sensing system is used to deploy sensing fibers along the rock mass monitoring section to collect dynamic strain data at the centimeter to meter level in real time. The sampling frequency of the dynamic strain data is not less than 1 kHz. Simultaneously, a microseismic sensor array is deployed, with a monitoring range covering a kilometer-level spatial scale. Three-component waveform data of rock mass fracture events are collected. A unified time-scale synchronization mechanism is established across the monitoring system. This synchronization mechanism is implemented through a GPS timing module to ensure that the time alignment accuracy of multi-source monitoring data reaches the millisecond level. Non-local processing is performed on the dynamic strain data and waveform data. The data underwent part-mean denoising processing using an adaptive window width algorithm, with the window width ranging from 3 to 15 sampling points. Based on the denoised microseismic data, the acoustic emission energy spectrum of each monitoring zone was calculated using short-time Fourier transform. The b-value distribution was then solved using the maximum likelihood estimation method based on the energy spectrum. A three-dimensional point cloud reconstruction was performed on the fiber optic strain data, employing the Delaunay triangulation algorithm to extract fracture aperture, dip direction, and dip angle parameters. Rock sample permeability data was obtained through transient pressure pulse tests, and an equivalent permeability tensor was constructed by combining fracture network parameters. The angle deviation between the principal direction of the permeability tensor and the direction of the maximum principal stress did not exceed 15°.

[0030] Distributed fiber optic acoustic sensing systems refer to monitoring devices based on fiber optic sensing principles. Specifically, they can be implemented using Bragg grating arrays or phase-sensitive optical time-domain reflectometry (PTZ) technology, used to capture the dynamic strain response of rock masses in real time at the centimeter to meter scale. Microseismic sensor arrays refer to a spatially distributed network composed of multiple three-component accelerometers, specifically implemented using broadband piezoelectric sensors, used to capture elastic wave signals from kilometer-scale rock mass fracture events. Nonlocal mean denoising refers to filtering algorithms based on signal similarity, specifically achieved by calculating the weighted average of samples within a data window, used to eliminate non-stationary noise in cross-scale monitoring data. Adaptive window width algorithms refer to methods that dynamically adjust the size of the filtering window based on local signal characteristics. Specifically, a sliding variance estimator can be used to determine the optimal window width, used to balance denoising effect with signal detail preservation. Delaunay triangulation algorithms refer to geometric reconstruction methods based on the empty circle criterion, specifically generating triangular meshes by maximizing the minimum interior angle criterion, used to recover the three-dimensional morphology of fractures from discrete strain data. The equivalent permeability tensor is a second-order tensor that characterizes the anisotropic seepage capacity of a fracture network. Specifically, it can be obtained by statistical averaging of fracture aperture and direction distribution, and is used to establish the correlation between the macroscopic seepage field and the microscopic fracture structure.

[0031] A distributed fiber optic sensing system and a microseismic sensor array were deployed at the rock mass monitoring section to synchronously acquire dynamic response data at different spatial scales. The fiber optic system captured strain fluctuations from centimeters to meters at a sampling rate of at least 1 kHz, while the microseismic array recorded the three-component waveforms of rock mass fracture events covering a kilometer-scale range. Multi-source data were synchronized at the millisecond level via a GPS timing module to ensure temporal consistency across scales. Raw data underwent nonlocal mean filtering, and an adaptive window width of 3-15 sampling points was used to remove noise interference. A short-time Fourier transform was performed on the denoised microseismic data to extract the acoustic emission energy spectrum characteristics of each monitoring zone, and the b-value distribution was calculated based on the maximum likelihood estimation method. The fiber optic strain data was reconstructed into a three-dimensional point cloud using Delaunay triangulation to quantitatively extract fracture aperture, dip direction, and dip angle parameters. Combined with rock sample permeability data obtained from transient pressure pulse tests, an equivalent permeability tensor was constructed with a deviation of no more than 15° between the principal direction and the maximum principal stress, providing fundamental parameters for subsequent multi-field coupled analysis.

[0032] Compared to existing technologies, traditional methods typically employ a single type of sensor for localized monitoring, failing to achieve synchronous data acquisition across scales from centimeters to kilometers. Existing noise reduction techniques often use fixed-window-width filtering, which is ill-suited to the non-stationary characteristics of rock mass signals. Fracture parameter extraction relies heavily on two-dimensional slice analysis, failing to accurately reflect three-dimensional spatial distribution characteristics. Permeability calculations often neglect the spatial correlation between the principal stress direction and the principal seepage direction, leading to equivalent parameters deviating from reality. This solution addresses the challenge of spatiotemporal alignment of cross-scale data through multi-sensor collaborative deployment and GPS synchronization; it employs adaptive-window-width nonlocal mean filtering to improve noise reduction of non-stationary signals; the fracture parameter extraction method based on three-dimensional point cloud reconstruction enhances spatial characterization accuracy; and the construction of a permeability tensor constrained by principal direction deviation improves the reliability of coupled analysis of the seepage field and stress field.

[0033] Through the above technical solutions, this application achieves high-precision synchronous acquisition and noise reduction of cross-scale rock mass response data, solving the fusion error problem caused by the lack of unified spatiotemporal benchmarks in traditional methods; by reconstructing three-dimensional fracture networks and calculating equivalent permeability tensors, the parameter input accuracy of multi-field coupling analysis is improved; and by adopting adaptive noise reduction and physical constraint parameter extraction methods, the reliability of characterization of rock mass feature parameters under complex environments is enhanced.

[0034] Example 3: Please refer to Figure 1 The specific method for step S2 is as follows: S2.1 Construct a nonlinear damage constitutive model, consider the mechanical response of the rock mass under different loading paths, describe the nonlinear deformation and damage evolution of the rock mass under high stress environment, introduce the historical stress path weight factor to reflect the influence of historical stress on the damage evolution of the rock mass, and reflect the time-varying and nonlinear characteristics of the rock mass in the process of multiple loading and unloading. S2.2 Establish the control equations for the seepage field, and use the Forchheimer modified equations to describe the high-speed nonlinear seepage behavior, reflecting the flow characteristics of fluid in the rock mass and its interaction with rock mass damage. Couple the damage field and seepage field through the Biot effective stress principle to establish the physical connection between the two, and simulate the permeability changes caused by damage and the feedback of seepage on rock mass damage.

[0035]

[0036] in, Z is the effective stress, A is the coupling coefficient, D is the seepage field pressure, and B is the adjustment coefficient. In this embodiment: This application further proposes to construct a nonlinear damage constitutive model, considering the mechanical response of the rock mass under different loading paths, describing the nonlinear deformation and damage evolution of the rock mass under high stress environment, introducing a historical stress path weight factor to reflect the influence of historical stress on the damage evolution of the rock mass, and reflecting the time-varying and nonlinear characteristics of the rock mass during multiple loading and unloading processes; establishing the seepage field control equation, using the Forchheimer modified equation to describe the high-speed nonlinear seepage behavior, reflecting the flow characteristics of fluid in the rock mass and its interaction with rock mass damage, and coupling the damage field and seepage field through the Biot effective stress principle to establish the physical connection between the two, simulating the permeability change caused by damage and the feedback of seepage on rock mass damage.

[0037] Nonlinear damage constitutive models are mathematical models describing the accumulation of damage and degradation of mechanical properties in rock masses under complex stress conditions. Specifically, they can be implemented using piecewise functions or tensor equations that consider stress path dependence, and are used to characterize the nonlinear deformation behavior of rock masses under high stress environments. Historical stress path weighting factors are parameters that quantify the influence of different historical stress stages on the current damage state. Specifically, they can be implemented using exponential decay functions or memory integrals, and are used to reflect the time-varying characteristics of rock masses during multiple loading and unloading processes. The Forchheimer modified equations are differential equations describing the nonlinear flow of high-speed fluids in porous media. Specifically, they can be implemented using a mathematical form that correlates the quadratic velocity term with the pressure gradient, and are used to characterize the influence of fluid inertia on seepage behavior. The Biot effective stress principle is a theoretical framework that combines solid skeleton stress and fluid pore pressure into effective stress. Specifically, it can be implemented by simultaneously solving the pore elastic constitutive equation and the seepage equation, and is used to establish a two-way coupling relationship between the damage field and the seepage field.

[0038] When constructing the nonlinear damage constitutive model, the cumulative effect of different loading stages on the rock mass damage evolution can be quantified by introducing historical stress path weighting factors. For example, an exponential decay function is used to weight historical stresses, making the influence of recent loading paths on the current damage state more significant. This model can accurately describe the stiffness degradation and residual deformation characteristics of rock masses under cyclic loading. When establishing the seepage field governing equations, the Forchheimer modified equation is used instead of the traditional Darcy law. By introducing a velocity square term to characterize the inertial effect in high-speed seepage, the nonlinear flow behavior of fluids in high-permeability rock masses can be simulated more accurately. Simultaneously, based on the Biot effective stress principle, the permeability change caused by damage and the pore pressure change generated by seepage are bidirectionally coupled. For example, the dynamic influence of the damage field on the seepage field is realized through the functional relationship between the permeability tensor and the damage variable, and the feedback effect of the seepage field on the damage field is realized through the correction of effective stress by pore pressure.

[0039] Compared to existing technologies, traditional methods typically employ linear damage models and neglect the influence of historical stress paths, resulting in an inability to accurately simulate the nonlinear mechanical behavior of rock masses under complex loading conditions. Existing seepage analyses are mostly based on Darcy's law, which struggles to describe the nonlinear characteristics of high-speed seepage, and the coupling between damage and the seepage field is often unidirectional or weakly coupled. This scheme enhances the model's ability to characterize time-varying damage behavior of rock masses by introducing a historical stress path weighting factor; it expands the applicability of the seepage model by employing the Forchheimer modified equation, enabling a more accurate characterization of high-speed seepage phenomena; and the bidirectional coupling mechanism implemented through the Biot principle effectively captures the physical essence of the interaction between damage and seepage.

[0040] Through the above technical solutions, this application can more accurately simulate the nonlinear mechanical response of deep rock masses under multi-physics field coupling, overcoming the limitations of traditional methods in effectively characterizing the influence of historical stress and high-speed seepage behavior. By employing a two-way coupling mechanism between the damage field and the seepage field, dynamic simulation of the interaction between rock mass damage evolution and seepage characteristics is achieved, providing a more reliable numerical model basis for the stability analysis of deep rock masses.

[0041] Example 4: Please refer to Figure 1 The specific method for step S3 is as follows: S3.1 Determine the set of parameters to be optimized. The parameters have an important impact on the analysis of rock mass damage and seepage coupling. The selected parameters are screened based on the physical properties and mechanical response of the rock mass. S3.2. Use a 256-bit quantum annealing machine to pre-optimize the set of parameters to be optimized. Use the quantum annealing algorithm to search and filter the feasible domain of parameters and output the optimized parameter range.

[0042] In this embodiment: This application further proposes a specific implementation method for a four-dimensional coupling analytical method for multiphysics fields in deep rock masses, including determining a parameter set to be optimized. The parameter set includes parameters that have an important influence on the coupling analysis of rock mass damage and seepage. The selected parameters are screened based on the physical properties and mechanical response of the rock mass. A 256-bit quantum annealing machine is used to pre-optimize the parameter set to be optimized. The feasible domain of the parameters is searched and screened through the quantum annealing algorithm, and the optimized parameter range is output.

[0043] The parameter set refers to the set of key parameters affecting the coupling analysis of rock mass damage and seepage. Specifically, it may include parameters such as the damage threshold, Biot coefficient, and Forchheimer coefficient. These parameters are selected based on the physical properties and mechanical response of the rock mass and are core variables used to characterize the nonlinear damage and seepage coupling behavior of the rock mass. The quantum annealing machine is an optimization device based on quantum computing principles. Specifically, it can be implemented using a 256-bit quantum annealing processor. Through the quantum tunneling effect, it performs a global search within the feasible region of parameters, solving high-dimensional nonlinear problems that are difficult to handle with traditional optimization methods. The feasible region of parameters refers to the range of values ​​for the parameters to be optimized. Specifically, its boundaries can be determined through rock mass physical property test data and engineering experience. The quantum annealing algorithm is used for preliminary screening, effectively narrowing the parameter search space.

[0044] In the parameter pre-optimization stage, key parameters, such as damage threshold and Biot coefficient (8-10 parameters in total), are first selected based on the physical properties and mechanical response of the rock mass, forming a set of parameters to be optimized. Then, the feasible region of the parameters is mapped to the computational space of the quantum annealing machine. The quantum annealing algorithm is used to perform a global search in the parameter combination space, overcoming the limitations of local optima through quantum tunneling, and quickly selecting an optimized parameter range that meets physical constraints. This process, leveraging the efficient parallel computing capabilities of quantum computing devices, completes the parameter space exploration that would take traditional methods several days in just a few hours, providing a high-quality initial parameter range for subsequent accurate inversion.

[0045] Compared to existing technologies, traditional parameter optimization methods typically employ gradient descent or genetic algorithms, which suffer from low computational efficiency and susceptibility to local optima when dealing with high-dimensional nonlinear parameter spaces. Our proposed solution, however, utilizes quantum annealing to perform a global search within the feasible parameter domain, overcoming the computational bottleneck of traditional optimization algorithms and significantly improving the efficiency and reliability of parameter pre-optimization. Furthermore, existing technologies often rely on empirically determined initial parameter values, lacking a systematic screening mechanism. Our solution, however, employs quantum computing equipment to scientifically screen the parameter range, effectively avoiding biases imposed by human experience.

[0046] Through the above technical solution, this application can significantly reduce the computational complexity of subsequent accurate parameter inversion. Quantum annealing pre-optimization compresses the parameter search space to a reasonable range, avoiding redundant calculations of invalid parameter combinations. Simultaneously, this method, through a physical property-driven parameter selection mechanism, ensures that the optimized parameter range conforms to the actual mechanical behavior characteristics of the rock mass, laying the foundation for accurate solutions to multiphysics coupling models.

[0047] Example 5: Please refer to Figure 1 The specific method for step S4 is as follows: S4.1 Define the likelihood relationship between the observed data and the model output, set the prior distribution and range of the parameters, and give inequalities and boundary constraints, including the Biot coefficient range, the damage threshold range, and the non-negativity constraints of the Forchheimer coefficient and permeability. Introduce the light cone constraint to limit the spatiotemporal propagation sequence of stress and seepage interaction and ensure that the propagation speed does not exceed the preset upper bound velocity parameter. Perform dimensionless and scaled operation based on the parameter range obtained in S3, initialize the adaptive Metropolis-Hastings proposal distribution and its covariance structure, and form a posterior target model containing the above constraints. S4.2 Start sampling and perform preheating iteration. Update the covariance and step size of the proposed distribution according to the iteration results. Calculate the acceptance probability and perform acceptance or rejection. Perform rejection or constraint projection processing on samples that do not meet the constraints. Make a termination decision based on convergence criteria such as effective sample size and inter-chain stability. Select parameter estimates and uncertainty quantification results based on the posterior sample set. Output parameter combinations as input for subsequent steps.

[0048] In this embodiment: This application further proposes a specific implementation method for step S4, including defining the likelihood relationship between the observed data and the model output, setting the prior distribution and value range of the parameters, applying inequalities and boundary constraints, introducing light cone constraints to limit the spatiotemporal propagation sequence of stress and seepage interaction, implementing dimensionless and scaling based on the pre-optimized parameter range, initializing the adaptive Metropolis-Hastings proposal distribution and its covariance structure, performing preheating iteration and updating the covariance and step size, calculating the acceptance probability and performing sample processing, terminating sampling according to the convergence criterion and outputting the parameter combination.

[0049] The light cone constraint ensures causal consistency by limiting the spatiotemporal propagation speed of the stress-seepage interaction to a preset upper bound. This can be achieved using relativistic propagation speed constraint equations to prevent superluminal propagation phenomena in the physical model. The adaptive Metropolis-Hastings algorithm dynamically adjusts the proposed distribution parameters based on iterative results. This can be implemented using an adaptive covariance matrix mechanism, improving parameter space exploration efficiency by updating the covariance structure. Constraint projection processing mathematically corrects samples that violate physical constraints. This can be achieved using the Lagrange multiplier method or boundary projection algorithm to ensure that parameter estimation results meet physical requirements.

[0050] This step transforms the parameter range obtained from quantum pre-optimization into a standardized search space by establishing a posterior probability model with multi-dimensional constraints. During the initialization phase, an adaptive proposal distribution is constructed, and the sampling step size and covariance structure are dynamically adjusted through preheating iterations to form an efficient search strategy that matches the parameter space. A dual verification mechanism is employed during sampling: high-quality samples are selected based on acceptance probability, while constrained projection is applied to out-of-bounds samples to ensure that all valid samples satisfy physical constraints such as the Biot coefficient range and damage threshold range. The introduction of optical cone constraints effectively limits the temporal relationship between stress wave propagation and seepage diffusion, avoiding time-reversed numerical solutions. In the convergence determination phase, the effective sample size and inter-chain stability index are comprehensively evaluated; the calculation terminates and the optimal parameter combination is output when a preset threshold is met.

[0051] Compared to existing technologies, traditional parameter inversion methods often ignore the rigid constraints of physical limitations, leading to optimization results that deviate from actual operating conditions. Existing Markov chain Monte Carlo methods, employing fixed step sizes and covariance structures, struggle to adapt to multi-modal parameter spaces. This proposed solution significantly improves parameter search efficiency under complex constraints through adaptive covariance adjustment and constraint projection. Compared to conventional soft constraint methods, the introduction of light cone constraints fundamentally guarantees the causal temporal sequence of field variable interactions, resolving the model distortion problem caused by stress-seepage temporal misalignment in traditional methods.

[0052] Through the above technical solutions, this application achieves high-precision parameter inversion under complex physical constraints, ensuring that the optimization results simultaneously satisfy mathematical optimality and physical rationality. The adaptive sampling strategy effectively overcomes the local extremum trap in the multi-constraint parameter space, and the light cone constraint mechanism eliminates temporal contradictions in the field variable interaction process. Constraint projection processing maintains the numerical stability of the parameter estimation process, providing reliable input parameters for subsequent multi-field coupling calculations.

[0053] Example 6: Please refer to Figure 1 The specific method for step S5 is as follows: S5.1 Determine the computational domain, boundary conditions, and initial conditions. Define and dimensionless variables based on the governing equations and parameter sets from S2 to S4. Use unstructured grids to geometrically discretize the rock mass. Implement local densification in fracture development zones and seepage gradient concentration zones. The characteristic scale of key area units should not exceed 0.1 mm. Establish the discretization forms of the mechanical field and seepage field and the global sparse matrix topology. Determine the time discretization scheme and step size control criteria. Establish the mapping interface between observation data and field quantities to form a spatiotemporal discretization framework for multi-field coupled solution. S5.2. Explicit time integration is used for the mechanical subproblem, and implicit time integration is used for the seepage subproblem. Inter-field coupling and variable exchange are performed according to the split iteration or alternating synchronization strategy. The Biot coupling term, damage variable and equivalent permeability are updated. Numerical stability and convergence are checked and the time step is adaptively adjusted. The field variables and damage evolution results are given in time series. The damage propagation path and range are recorded at a grid resolution of 0.1 mm. Discrete solution sequence and related derived indices are obtained.

[0054] In this embodiment, the specific method of step S5 is further proposed as follows: The computational domain, boundary conditions, and initial conditions are determined; variables are defined and dimensionless based on the governing equations and parameter sets; unstructured grids are used for geometric discretization of the rock mass; local densification is implemented in fracture development zones and seepage gradient concentration zones; the characteristic scale of key regional units is no greater than 0.1 mm; the discrete forms of the mechanical field and seepage field and the global sparse matrix topology are established; the time discretization scheme and step size control criteria are determined; a mapping interface between observation data and field quantities is established; a spatiotemporal discretization framework for multi-field coupled solution is constructed; explicit time integration is used for the mechanical subproblem, and implicit time integration is used for the seepage subproblem; inter-field coupling and variable exchange are performed according to a split iteration or alternating synchronization strategy; the Biot coupling term, damage variable, and equivalent permeability are updated; numerical stability and convergence checks are performed and the time step is adaptively adjusted; the field variables and damage evolution results are given according to the time series; the damage propagation path and range are recorded at a 0.1 mm grid resolution; and the discrete solution sequence and related derived indices are obtained.

[0055] Unstructured meshes refer to discretized mesh structures composed of irregularly shaped elements. They can be generated using the Delaunay triangulation algorithm or the advancing front method, adapting to complex geometric boundaries and enabling localized refinement. Local refinement involves increasing mesh density in specific regions, such as employing a gradient-decreasing element size strategy in fracture-developing zones. By setting the number of refinement layers to control the characteristic scale of the control system, computational accuracy in critical areas is ensured. Explicit-implicit hybrid algorithms use explicit time integration to solve for the instantaneous dynamic response of the mechanical equations, while implicit methods are used to handle the steady-state process of the seepage equations. For example, the Newmark-β method combined with Newton-Raphson iteration is used to achieve alternating updates of inter-field variables. Adaptive time step control dynamically adjusts the time increment based on the convergence state of the numerical solution. For example, when the residual norm exceeds a threshold, the step size is automatically shortened to ensure computational stability.

[0056] When constructing the spatiotemporal discretization framework, the boundary conditions of the computational domain are first defined based on geological structural characteristics, such as using measured in-situ stress as the initial stress field input. When geometrically discretizing the rock mass using an unstructured grid, mesh refinement is implemented in the fracture propagation path prediction region; for example, the lower limit of the element size is set to 0.1 mm in the identified high stress gradient region. In the time discretization scheme, the mechanical field uses explicit integration to capture the dynamic damage process, such as using the central difference method to solve the acceleration term; the seepage field uses implicit integration to handle the diffusion process, such as using the backward Euler method to ensure mass conservation. Inter-field coupling is achieved through alternating iterations; for example, in each time step, the mechanical field is solved first to update the displacement field, and then the deformation results are transferred to the seepage field to calculate the pore pressure distribution. By real-time monitoring of the convergence indices of the numerical solution, such as the residual descent rate or energy error, an adaptive adjustment mechanism for the time step is triggered.

[0057] Compared to existing technologies, traditional methods often employ structured meshes, leading to insufficient accuracy in fitting fracture boundaries, and single time integration methods struggle to balance computational efficiency for both mechanical and seepage fields. This proposed solution, through unstructured meshes and a local refinement strategy, can accurately characterize the geometric details of fracture propagation paths, such as capturing microcrack bifurcation behavior at a scale of 0.1 mm. The hybrid integration algorithm effectively balances the time step limitations of explicit methods with the computational costs of implicit methods, maintaining real-time dynamic solutions for the mechanical field while ensuring the accuracy of the steady-state solution. Existing technologies lack dynamic step-size control mechanisms, while this proposed solution avoids numerical oscillations or computational redundancy caused by traditional fixed step sizes through residual monitoring and adaptive adjustment.

[0058] Through the above technical solutions, this application resolves the contradiction between accuracy and efficiency in complex geometric modeling and multi-field coupled calculations using traditional discretization methods. Unstructured meshes and local refinement strategies improve the geometric representation accuracy of fracture propagation paths, while an explicit-implicit hybrid algorithm optimizes the efficiency of co-solving mechanical and seepage fields. An adaptive step-size control mechanism ensures the numerical stability of long-term, multi-scale simulations. This technical solution effectively reduces the computational resource consumption of multi-physics coupled analysis while maintaining computational accuracy, providing a high-resolution numerical simulation foundation for the dynamic damage evolution of deep rock masses.

[0059] Example 7: Please refer to Figure 1 The specific method for step S6 is as follows: S6.1. Based on the time-series discrete field data of S5, coordinate and time are aligned with the field monitoring data. The fracture skeleton is extracted according to the damage variable, fracture aperture and equivalent permeability threshold. An undirected graph representation composed of fracture segments and intersection points is constructed. A topological filtering sequence with fracture aperture or permeability index as filtering parameters is established. The persistence graph and barcode are generated by the continuous cohomology method. Quantitative indicators such as zero-order and first-order Betti numbers, connected component persistence, main channel span and skeleton path are extracted to form a topological feature set and its spatial mapping organized according to the time series. S6.2 Establish early warning judgment rules and threshold system, and jointly judge the change rate of topological index and the persistence threshold of key connected domains with the change of acoustic emission b value. When the change of acoustic emission b value is greater than 0.5 or the topological index exceeds the control threshold, an early warning is triggered. The judgment is updated by using a sliding time window and hysteresis strategy. Risk scores and classification results are generated according to the computational grid or unit. The spatial zoning and unit list of possible water inrush or rockburst are output, and a risk assessment report containing the trigger time, trigger index and spatial range description is generated.

[0060] In this embodiment: This application further proposes deploying a field monitoring network and unifying access, configuring distributed fiber optic acoustic wave, micro-vibration, stress, pore pressure, seepage and temperature sensing units, setting up 5G IoT gateways and edge computing nodes, establishing a unified coordinate system and time reference and performing time synchronization and pose calibration, constructing a data acquisition and preprocessing pipeline, performing noise reduction, drift correction and missing data completion on the data and transmitting it through an encrypted channel, establishing a message queue and time series database and completing subject cataloging and metadata registration, providing a mapping relationship matrix between sensor tags and model variables and assigning unique identifiers, establishing the association between the digital twin model and the field entity and three-dimensional geometry, forming a continuous and traceable data input channel; performing online model calibration based on real-time data stream, setting trigger conditions and sliding time windows and determining parameter update frequency, using incremental learning algorithms to dynamically update and correct parameters and apply physical constraints and stability criteria, completing the synchronization and visualization refresh of the digital twin model status and generating structured outputs of risk and health indicators, establishing a data interface with the engineering decision-making platform to send back updated parameters and evaluation results and record version numbers and timestamps, forming a closed-loop process of monitoring, calibration and feedback.

[0061] A 5G IoT gateway refers to an IoT data transmission device that supports the 5G communication protocol. It can be implemented using multi-band antennas and low-latency transmission modules to achieve high-speed, real-time transmission of on-site monitoring data. Edge computing nodes are local computing units deployed at the monitoring site. They can be implemented using embedded processors and distributed storage architectures to perform data preprocessing and preliminary analysis. Incremental learning algorithms are machine learning parameter update methods based on dynamic data flows. They can be implemented using online gradient descent or recursive least squares methods to adjust model parameters based on real-time monitoring data. Optical cone constraints are mathematical conditions that limit the spatiotemporal propagation range of physical field interactions. They can be implemented using a discretized form of relativistic causality to ensure temporal consistency between stress and seepage fields. Digital twin models are virtual simulation models that are updated synchronously with physical entities. They can be implemented using a 3D visualization engine and real-time data interfaces to achieve dynamic interaction between monitoring data and the computational model.

[0062] By deploying multiple types of sensing units and 5G IoT gateways, a three-dimensional monitoring network covering multiple physical quantities is formed. Edge computing nodes are used to locally process raw data, reducing data transmission latency. A unified coordinate system and time benchmark are used to achieve spatiotemporal alignment of multi-source data, and an encrypted transmission channel is constructed to ensure data security. A standardized data management mechanism is established based on message queues and time-series databases to achieve accurate mapping between monitoring data and model variables. During the online model calibration phase, incremental learning algorithms combined with physical constraints are used to dynamically update parameters, avoiding cumulative errors caused by environmental changes. Through the visualization and structured output of the digital twin model, real-time feedback of risk status is achieved, ultimately forming a closed-loop monitoring and calibration system.

[0063] Compared to existing technologies, traditional monitoring systems often employ independently deployed single-sensor networks, lacking a mechanism for collaborative acquisition of multiple physical quantities and unified data processing, leading to difficulties in data fusion. Existing model calibration methods typically use offline batch processing, which cannot adapt to dynamic changes in the field environment, and the lack of an incremental learning mechanism with physical constraints easily results in non-physical interpretations. This solution achieves efficient integration of multi-source heterogeneous data through a collaborative architecture of 5G IoT and edge computing. Combined with a dynamic parameter update mechanism based on incremental learning and physical constraints, it effectively improves the real-time performance and reliability of online model calibration.

[0064] Through the above technical solutions, this application achieves full-element acquisition and standardized management of deep rock mass monitoring data, solving the problem of spatiotemporal alignment difficulties for multi-source heterogeneous data. By leveraging the synergistic effect of incremental learning algorithms and physical constraints, the physical rationality of model parameters in dynamic environments is ensured, overcoming the shortcomings of traditional offline calibration methods in terms of adaptability. The closed-loop interaction mechanism between the digital twin model and the engineering decision-making platform provides reliable technical support for real-time assessment and early warning of rock mass safety status.

[0065] Example 8: Please refer to Figure 1 The specific method for step S7 is as follows: S7.1 Deploy a field monitoring network and access it uniformly. Configure distributed fiber optic acoustic wave, micro-vibration, stress, pore pressure, seepage and temperature sensing units. Set up 5G IoT gateways and edge computing nodes. Establish a unified coordinate system and time reference and perform time synchronization and pose calibration. Build a data acquisition and preprocessing pipeline. Denoise, drift correction and missing data are performed on the data and transmitted through an encrypted channel. Establish a message queue and time series database and complete subject cataloging and metadata registration. Provide a mapping relationship matrix between sensor tags and model variables and assign unique identifiers. Establish the association relationship between the digital twin model and the field entity and three-dimensional geometry to form a continuous and traceable data input channel. S7.2. Based on the real-time data stream of step S7.1, perform online calibration of the model, set trigger conditions and sliding time windows and determine the parameter update frequency, use incremental learning algorithm to dynamically update and correct drift of parameters, apply physical constraints and stability criteria, complete the synchronization and visualization refresh of the digital twin model status and generate structured output of risk and health indicators, establish a data interface with the engineering decision platform to send back updated parameters and evaluation results and record version number and timestamp, forming a closed-loop process of monitoring, calibration and feedback.

[0066] In this embodiment: This application further proposes deploying a field monitoring network and unifying access, configuring distributed fiber optic acoustic wave, micro-vibration, stress, pore pressure, seepage and temperature sensing units, setting up 5G IoT gateways and edge computing nodes, establishing a unified coordinate system and time reference and performing time synchronization and pose calibration, constructing a data acquisition and preprocessing pipeline, performing noise reduction, drift correction and missing data completion on the data and transmitting it through an encrypted channel, establishing a message queue and time series database and completing subject cataloging and metadata registration, providing a mapping relationship matrix between sensor tags and model variables and assigning unique identifiers, establishing the association between the digital twin model and the field entity and three-dimensional geometry, forming a continuous and traceable data input channel; performing online model calibration based on real-time data stream, setting trigger conditions and sliding time windows and determining parameter update frequency, using incremental learning algorithms to dynamically update and correct parameters and apply physical constraints and stability criteria, completing the synchronization and visualization refresh of the digital twin model status and generating structured outputs of risk and health indicators, establishing a data interface with the engineering decision-making platform to send back updated parameters and evaluation results and record version numbers and timestamps, forming a closed-loop process of monitoring, calibration and feedback.

[0067] A 5G IoT gateway refers to an IoT access device that supports the 5G communication protocol. Specifically, it can be implemented using an industrial-grade gateway with multi-protocol conversion capabilities, enabling high-speed transmission and low-latency communication of massive monitoring data. Edge computing nodes are local computing units deployed at the monitoring site, such as embedded processors based on the ARM architecture, used to perform data preprocessing and real-time analysis tasks, reducing the computing load on the cloud. Incremental learning algorithms are machine learning methods that can dynamically adjust model parameters based on newly arriving data, such as online gradient descent algorithms, used to iteratively update parameters while ensuring model stability. Physical constraints and stability criteria refer to parameter boundary conditions set based on rock mechanics principles, such as the non-negativity constraint of permeability, used to ensure that the model update process conforms to physical laws. Digital twin model state synchronization refers to the real-time matching of on-site monitoring data with the simulation model, for example, through data assimilation technology, used to maintain consistency between the virtual model and the actual engineering state.

[0068] The on-site monitoring network collects dynamic response data of the rock mass through multiple types of sensors. A 5G IoT gateway encrypts and transmits the data to edge computing nodes for preliminary processing. A unified coordinate system and time reference are established using a GPS timing module and a laser positioning system, ensuring spatiotemporal alignment of multi-source data. The data preprocessing pipeline employs a sliding window mechanism to correct for noise, drift, and missing data. Message queues and time-series databases buffer and structure the data stream. The digital twin model establishes a mapping relationship with the on-site entity using 3D geometric modeling tools. Incremental learning algorithms dynamically adjust model parameters based on real-time data streams, while physical constraints limit the parameter update range. The calibrated model output displays the damage evolution process through a visual interface, and risk assessment results are pushed to the engineering decision-making platform via a standard interface, forming a closed-loop feedback chain from data acquisition to decision support.

[0069] Compared to existing technologies, traditional monitoring systems often employ independently deployed sensor networks. Data synchronization relies on manual calibration and lacks unified management, resulting in long model update cycles and an inability to reflect real-time changes in rock mass conditions. This solution integrates 5G communication and edge computing technologies to achieve millisecond-level synchronization of multi-source data. It utilizes incremental learning algorithms to overcome the limitations of traditional offline calibration methods, enabling dynamic optimization of model parameters based on monitoring data. Simultaneously, physical constraints ensure the stability of the update process, resolving the issues of data lag and model inaccuracy inherent in traditional methods.

[0070] Through the above technical solutions, this application achieves real-time acquisition and dynamic fusion of deep rock mass monitoring data, constructs a digital twin system with self-calibration capabilities, and effectively improves the spatiotemporal resolution of the multiphysics coupling model. By using a closed-loop feedback mechanism to transform on-site monitoring data into engineering decision-making data in real time, it significantly enhances the timeliness and accuracy of early warnings for disasters such as water inrush and rock bursts, providing reliable technical support for the safety of deep rock mass engineering.

[0071] Furthermore, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated units described above can be implemented in hardware or as software functional units. The above are merely embodiments of this application and do not limit the patent scope of this application. Any equivalent structural or procedural transformations made based on the description and drawings of this application, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of this application.

[0072] The specific embodiments of the invention have been described in detail above, but they are only examples, and this application is not limited to the specific embodiments described above. For those skilled in the art, any equivalent modifications or substitutions to the invention are also within the scope of this application. Therefore, all equivalent changes, modifications, and improvements made without departing from the spirit and principles of this application should be covered within the scope of this application.

Claims

1. A four-dimensional coupled analytical method for multiphysics fields in deep rock masses, characterized in that: The specific steps are as follows: S1. Collect rock mass deformation data. Use distributed fiber optic acoustic sensing technology and microseismic sensor array to obtain rock mass response data from centimeter to kilometer scale. Preprocess the data and use a nonlocal mean filtering algorithm to eliminate cross-scale data noise. Extract rock mass response characteristic parameters, including acoustic emission b-value, fracture aperture and permeability parameters. S2. Construct a nonlinear damage constitutive model. The constitutive model introduces the historical stress path weight factor to characterize the non-Markov characteristics, establishes the seepage field control equation, and uses the Forchheimer modified equation to describe the high-speed nonlinear seepage behavior, realizing the bidirectional coupling between the damage field and the seepage field. The physical relationship between field variables is established through the Biot effective stress principle. S3. Determine the set of parameters to be optimized. The set of parameters includes 8-10 key parameters, including the damage threshold and Biot coefficient. Perform parameter pre-optimization on the quantum computing device. Use a 256-bit quantum annealing machine to initially screen the feasible domain of parameters and output the optimized parameter range to provide initial values ​​for subsequent accurate inversion. S4. Perform precise inversion calculation of execution parameters, use the adaptive Metropolis-Hastings algorithm for iterative optimization, apply physical constraints, including the introduction of light cone constraints to ensure the causal temporal consistency of stress-seepage field, and output the optimal parameter combination to ensure that the calculation results conform to physical reality. S5. Construct a spatiotemporal discrete computational model, use unstructured grids to achieve geometric discretization of the rock mass model, implement multi-field coupled calculations, solve the mechanical-seepage coupling problem through explicit-implicit hybrid algorithm, realize dynamic evolution simulation, and track the damage development process with a grid resolution of 0.1mm. S6. Analyze the topological evolution of the fracture network, quantitatively characterize the changes in fracture connectivity using the continuous homology method, establish a damage early warning mechanism, trigger an early warning signal when the change in acoustic emission b value Δb>0.5 is detected, generate a risk assessment report, and output the dangerous areas where water inrush or rock bursts may occur. S7. Deploy an engineering monitoring system to acquire real-time on-site monitoring data through 5G IoT, perform online model calibration, dynamically update model parameters using incremental learning algorithms, realize digital twin interaction, and feed the optimization results back to the engineering decision-making platform.

2. The deep rock mass multi-physical field four-dimensional coupling analytical method according to claim 1, characterized in that, The specific method for step S1 is as follows: S1.1 A distributed fiber optic acoustic sensing system is used to deploy sensing fibers along the rock mass monitoring section to collect dynamic strain data at the centimeter to meter level in real time. The sampling frequency of the dynamic strain data is not less than 1kHz. At the same time, a microseismic sensor array is deployed, and the monitoring range of the microseismic sensor array covers the spatial scale of kilometers. The three-component waveform data of rock mass fracture events are collected, and a unified time-scale synchronization mechanism is established across the monitoring system. The synchronization mechanism is realized through the GPS time synchronization module to ensure that the time alignment accuracy of multi-source monitoring data reaches the millisecond level. S1.

2. Nonlocal mean denoising is performed on dynamic strain data and waveform data. The denoising process uses an adaptive window width algorithm with a window width ranging from 3 to 15 sampling points. Based on the denoised microseismic data, the acoustic emission energy spectrum of each monitoring zone is calculated using short-time Fourier transform. The b-value distribution is solved using the maximum likelihood estimation method based on the energy spectrum. Three-dimensional point cloud reconstruction is performed on the fiber optic strain data. The reconstruction process uses the Delaunay triangulation algorithm to extract fracture aperture, dip direction, and dip angle parameters. Rock sample permeability data is obtained through transient pressure pulse tests. An equivalent permeability tensor is constructed by combining fracture network parameters. The angle deviation between the principal direction of the permeability tensor and the direction of the maximum principal stress does not exceed 15°.

3. The method of claim 2, wherein the method is characterized by, The specific method for step S2 is as follows: S2.1 Construct a nonlinear damage constitutive model, consider the mechanical response of the rock mass under different loading paths, describe the nonlinear deformation and damage evolution of the rock mass under high stress environment, introduce the historical stress path weight factor to reflect the influence of historical stress on the damage evolution of the rock mass, and reflect the time-varying and nonlinear characteristics of the rock mass in the process of multiple loading and unloading. S2.2 Establish the control equations for the seepage field, and use the Forchheimer modified equations to describe the high-speed nonlinear seepage behavior, reflecting the flow characteristics of fluid in the rock mass and its interaction with rock mass damage. Couple the damage field and seepage field through the Biot effective stress principle to establish the physical connection between the two, and simulate the permeability changes caused by damage and the feedback of seepage on rock mass damage.

4. The deep rock mass multi-physical field four-dimensional coupling analytical method according to claim 3, characterized in that, The specific method for step S3 is as follows: S3.1 Determine the set of parameters to be optimized. The parameters have an important impact on the rock mass damage and seepage coupling analysis. The selected parameters are screened based on the physical properties and mechanical response of the rock mass. S3.

2. Use a 256-bit quantum annealing machine to pre-optimize the set of parameters to be optimized. Use the quantum annealing algorithm to search and filter the feasible domain of parameters and output the optimized parameter range.

5. The method of claim 4, wherein the method is characterized by, The specific method for step S4 is as follows: S4.1 Define the likelihood relationship between the observed data and the model output, set the prior distribution and range of the parameters, and give inequalities and boundary constraints, including the Biot coefficient range, the damage threshold range, and the non-negativity constraints of the Forchheimer coefficient and permeability. Introduce the light cone constraint to limit the spatiotemporal propagation sequence of stress and seepage interaction and ensure that the propagation speed does not exceed the preset upper bound velocity parameter. Perform dimensionless and scaled operation based on the parameter range obtained in S3, initialize the adaptive Metropolis-Hastings proposal distribution and its covariance structure, and form a posterior target model containing the above constraints. S4.2 Start sampling and perform preheating iteration. Update the covariance and step size of the proposed distribution according to the iteration results. Calculate the acceptance probability and perform acceptance or rejection. Perform rejection or constraint projection processing on samples that do not meet the constraints. Make a termination decision based on convergence criteria such as effective sample size and inter-chain stability. Select parameter estimates and uncertainty quantification results based on the posterior sample set. Output parameter combinations as input for subsequent steps.

6. The deep rock mass multi-physical field four-dimensional coupling analytical method according to claim 5, characterized in that, The specific method for step S5 is as follows: S5.1 Determine the computational domain, boundary conditions, and initial conditions. Define and dimensionless variables based on the governing equations and parameter sets from S2 to S4. Use unstructured grids to geometrically discretize the rock mass. Implement local densification in fracture development zones and seepage gradient concentration zones. The characteristic scale of key area units should not exceed 0.1 mm. Establish the discretization forms of the mechanical field and seepage field and the global sparse matrix topology. Determine the time discretization scheme and step size control criteria. Establish the mapping interface between observation data and field quantities to form a spatiotemporal discretization framework for multi-field coupled solution. S5.

2. Explicit time integration is used for the mechanical subproblem, and implicit time integration is used for the seepage subproblem. Inter-field coupling and variable exchange are performed according to the split iteration or alternating synchronization strategy. The Biot coupling term, damage variable and equivalent permeability are updated. Numerical stability and convergence are checked and the time step is adaptively adjusted. The field variables and damage evolution results are given in time series. The damage propagation path and range are recorded at a grid resolution of 0.1 mm. Discrete solution sequence and related derived indices are obtained.

7. The deep rock mass multi-physical field four-dimensional coupling analytical method according to claim 6, characterized in that, The specific method for step S6 is as follows: S6.

1. Based on the time-series discrete field data of S5, coordinate and time are aligned with the field monitoring data. The fracture skeleton is extracted according to the damage variable, fracture aperture and equivalent permeability threshold. An undirected graph representation composed of fracture segments and intersection points is constructed. A topological filtering sequence with fracture aperture or permeability index as filtering parameters is established. The persistence graph and barcode are generated by the continuous cohomology method. Quantitative indicators such as zero-order and first-order Betti numbers, connected component persistence, main channel span and skeleton path are extracted to form a topological feature set and its spatial mapping organized according to the time series. S6.2 Establish early warning judgment rules and threshold system, and jointly judge the change rate of topological index and the persistence threshold of key connected domains with the change of acoustic emission b value. When the change of acoustic emission b value is greater than 0.5 or the topological index exceeds the control threshold, an early warning is triggered. The judgment is updated by using a sliding time window and hysteresis strategy. Risk scores and classification results are generated according to the computational grid or unit. The spatial zoning and unit list of possible water inrush or rockburst are output, and a risk assessment report containing the trigger time, trigger index and spatial range description is generated.

8. The deep rock mass multi-physical field four-dimensional coupling analytical method according to claim 7, characterized in that, The specific method for step S7 is as follows: S7.1 Deploy a field monitoring network and access it uniformly. Configure distributed fiber optic acoustic wave, micro-vibration, stress, pore pressure, seepage and temperature sensing units. Set up 5G IoT gateways and edge computing nodes. Establish a unified coordinate system and time reference and perform time synchronization and pose calibration. Build a data acquisition and preprocessing pipeline. Denoise, drift correction and missing data are performed on the data and transmitted through an encrypted channel. Establish a message queue and time series database and complete subject cataloging and metadata registration. Provide a mapping relationship matrix between sensor tags and model variables and assign unique identifiers. Establish the association relationship between the digital twin model and the field entity and three-dimensional geometry to form a continuous and traceable data input channel. S7.

2. Based on the real-time data stream of step S7.1, perform online calibration of the model, set trigger conditions and sliding time windows and determine the parameter update frequency, use incremental learning algorithm to dynamically update and correct drift of parameters, apply physical constraints and stability criteria, complete the synchronization and visualization refresh of the digital twin model status and generate structured output of risk and health indicators, establish a data interface with the engineering decision platform to send back updated parameters and evaluation results and record version number and timestamp, forming a closed-loop process of monitoring, calibration and feedback.

Images