Intelligent Decision-Making Method and System for Tunnel Support in Cold Regions Based on Dynamic Feedback

Abstract

This invention provides an intelligent decision-making method and system for tunnel support in cold regions based on dynamic feedback, relating to the field of tunnel engineering technology. The method includes: acquiring historical monitoring data; performing feature decoupling to obtain freeze-thaw damage characteristics and deformation field components; predicting the evolution of the surrounding rock damage field based on embedded constraints; constructing a risk assessment map by performing graph topology mining on the prediction results; optimizing the support scheme by combining geological weak surface information and existing support status; and dynamically updating constraint parameters through surrounding rock response data to achieve intelligent optimization of support decisions. This invention can accurately identify the propagation law of surrounding rock damage, improving support efficiency and safety.

Description

Technical Field

[0001] This invention relates to the field of tunnel engineering technology, and in particular to an intelligent decision-making method and system for tunnel support in cold regions based on dynamic feedback. Background Technology

[0002] As my country's transportation infrastructure construction extends to high-altitude and cold regions, the number of tunnel projects in cold regions is increasing. These tunnels are exposed to seasonal freeze-thaw cycles, leading to significant deterioration of the mechanical properties of the surrounding rock. Rock stability has become a major risk factor for the safe operation of tunnel projects in cold regions. The freeze-thaw damage process of the surrounding rock in cold regions is complex, involving multi-field coupling of heat, water, and force, and exhibits significant spatiotemporal heterogeneity. Therefore, the rationality of the support design directly affects the safety and durability of the tunnel structure.

[0003] Currently, tunnel support decisions in cold regions mainly rely on experience and static design specifications. However, there are still problems such as a lack of accurate characterization and prediction of the dynamic evolution of freeze-thaw damage in the surrounding rock, difficulty in adapting to local differences and spatiotemporal evolution characteristics of the surrounding rock, a focus on apparent indicators such as deformation, failure to deeply explore the mechanism and propagation law of freeze-thaw damage, difficulty in identifying potential risk connectivity areas and critical propagation paths, and a general lack of dynamic feedback mechanisms based on implementation results, making it impossible to continuously optimize based on actual support effects. Summary of the Invention

[0004] This invention provides an intelligent decision-making method and system for tunnel support in cold regions based on dynamic feedback, which can at least solve some of the problems existing in the prior art.

[0005] A first aspect of this invention provides an intelligent decision-making method for tunnel support in cold regions based on dynamic feedback, comprising:

[0006] Historical monitoring data of the surrounding rock of tunnels in cold regions is acquired and time-aligned and spatially mapped to obtain a monitoring dataset. Feature decoupling is performed on the historical monitoring data to obtain freeze-thaw damage features and deformation field components. Based on the freeze-thaw damage features and pre-set embedded pore water phase transformation constraints and strength degradation constraints, coupled calculations are performed to predict the evolution trajectory of the surrounding rock damage field and obtain evolution prediction results.

[0007] Graph topology mining is performed on the evolution prediction results to construct an evolution graph with monitoring points as nodes and damage propagation as edges, and high-risk connected regions and critical propagation paths are identified. The stress concentration degree of the connected regions is calculated in combination with the deformation field components, and a risk assessment map is constructed.

[0008] Based on pre-acquired geological weak surface information and existing support status, and combined with the risk assessment map, a multi-objective optimization algorithm is used to obtain the support timing sequence and support type scheme corresponding to each connected area, thus obtaining the zonal support decision results;

[0009] Based on the zonal support decision results, support operations are performed on the surrounding rock and the response data of the surrounding rock after support is collected. Based on the surrounding rock response data, the suppression effect of the support operation is determined, and based on the suppression effect, Bayesian updates are performed on the embedded pore water phase change constraint and strength degradation constraint to obtain optimized constraint parameters. Real-time monitoring data is collected, and the instability characteristic parameters are solved based on the optimized constraint parameters to determine the corresponding optimal support decision scheme.

[0010] In one alternative implementation,

[0011] Historical monitoring data of the surrounding rock of tunnels in cold regions is acquired and subjected to time alignment and spatial mapping to obtain a monitoring dataset. Feature decoupling is then performed on the historical monitoring data to obtain freeze-thaw damage features and deformation field components, including:

[0012] Raw data collected by various sensors located at different cross sections and burial depths over a historical period are extracted and combined with the timestamps and spatial coordinates of the data transmission to construct historical monitoring data.

[0013] Cross-correlation analysis is performed on the timestamps of each sensor in the historical monitoring data to calculate the response delay difference of different sensors. Based on the response delay difference, the time offset of each sensor is determined and offset compensation is performed to obtain corrected data. Sensor channels with inconsistent sampling frequencies in the corrected data are resampled to obtain unified sampling data.

[0014] The unified sampling data is combined with the spatial coordinates of the corresponding sensors to divide the data into three-dimensional grid nodes. The Euclidean distance from each sensor to each grid node is calculated, and the mapping weight is determined based on the Euclidean distance. The temperature field value, displacement field value, and strain field value of each grid node are solved based on the mapping weight and the unified sampling data and combined to obtain the monitoring dataset.

[0015] The temperature field values ​​in the monitoring dataset are decomposed into time-frequency components to extract periodic fluctuation components and calculate the time derivative and spatial gradient of the temperature field values. Based on the time derivative and spatial gradient, the freeze-thaw interface is identified. The freeze-thaw damage characteristics are obtained by statistical analysis of the temperature change amplitude, time derivative peak value, and spatial gradient extreme value at the freeze-thaw interface. The deformation field components are obtained by spatial differentiation of the displacement field values ​​and strain field values ​​and extraction of the time evolution trend.

[0016] In one alternative implementation,

[0017] Based on the freeze-thaw damage characteristics and pre-set embedded pore water phase transformation constraints and strength degradation constraints, coupled calculations are performed to predict the evolution trajectory of the surrounding rock damage field. The evolution prediction results include:

[0018] The freeze-thaw interface in the freeze-thaw damage features is used as the phase change trigger boundary. The temperature change amplitude and time derivative peak value in the freeze-thaw damage features are extracted as phase change parameters. Based on the phase change trigger boundary and the phase change parameters, the pore water volume expansion rate is calculated in the pre-set embedded pore water phase change constraint and converted into an expansion stress field.

[0019] The stress concentration location is determined based on the spatial gradient extrema in the expansion stress field and freeze-thaw damage characteristics. The strain field value history sequence of the corresponding grid node in the monitoring dataset is extracted at the stress concentration location and accumulated over time to obtain the cumulative strain. Based on the cumulative strain and the pre-set strength degradation constraint, the elastic modulus attenuation coefficient and compressive strength reduction coefficient are calculated and coupled with the expansion stress field to obtain the equivalent stress field.

[0020] Using the equivalent stress field as the initial stress, and combining the freeze-thaw damage characteristics to set the evolution time step, stress redistribution is performed on the equivalent stress field within each time step, and damage points are determined. The damage points are connected and tracked in the spatial domain, and the occurrence order and position changes of each damage point are recorded and combined to obtain the evolution trajectory. The spatial distribution characteristics and temporal evolution characteristics of the damage field corresponding to the evolution trajectory are extracted and summarized to obtain the evolution prediction results.

[0021] In one alternative implementation,

[0022] Graph topology mining is performed on the evolution prediction results to construct an evolution graph with monitoring points as nodes and damage propagation as edges. High-risk connected regions and critical propagation paths are identified. The stress concentration degree of the connected regions is calculated by combining the deformation field components, and a risk assessment map is constructed, including:

[0023] Based on the spatial distribution characteristics of the damage field in the evolution prediction results, the spatial coordinates of each damage point are extracted and spatially matched with the grid nodes in the monitoring dataset. According to the matching results, the damage points are mapped to the nearest grid node and marked as graph nodes.

[0024] The position change vector of each damage point between adjacent time steps is obtained and the corresponding orientation angle and length are calculated. Damage points with the same orientation angle and continuous length are identified and combined to obtain a propagation chain. Graph nodes corresponding to adjacent damage points in the propagation chain are connected and directed edges are created. The weight of the directed edges is set based on the time interval between adjacent damage points to obtain an evolution graph. The directed edges in the evolution graph are traversed and fast edges are determined in combination with a pre-set weight threshold. The connectivity of the fast edges is detected and the set of nodes reachable through the fast edges is marked as a high-risk connected region.

[0025] Traverse all paths from the initial damage point to the boundary node of the high-risk connected region in the evolution graph and determine the velocity index based on the sum of the weights of all directed edges on the path. The path with the largest velocity index is taken as the critical propagation path.

[0026] The high-risk connected regions are matched with the deformation field components in the spatial dimension, and the strain rate and strain gradient at the corresponding positions are extracted. The spatial second derivatives of the strain rate and strain gradient are calculated and the stress concentration is determined. The risk assessment map is obtained by combining the high-risk connected regions and the evolution map.

[0027] In one alternative implementation,

[0028] Based on pre-acquired geological weak surface information and existing support status, and combined with the risk assessment map, a multi-objective optimization algorithm is used to obtain the support timing sequence and support type scheme corresponding to each connected region, resulting in the following zonal support decision results:

[0029] The fault strike and joint density in the pre-acquired geological weak surface information are superimposed and matched with the high-risk connected regions in the risk assessment map to identify the coupling zone and extract the weakening coefficient.

[0030] Extract the spatial location and strength grade of the constructed support from the pre-acquired existing support status, calculate the shortest distance between the coupling zone and the spatial location of the constructed support to obtain the coverage, and calculate the gap amount by the difference between the strength grade and the pre-acquired stress concentration degree.

[0031] Using the weakening coefficient, coverage, and gap amount as constraints, the spatial orientation of the pre-obtained critical propagation path as boundary conditions, and the stress concentration degree of each high-risk connected area as risk weight, an optimization function is constructed with the objectives of minimizing the overall risk exposure time and maximizing the support resource utilization rate. The optimal intervention time and optimal strength configuration are obtained by iteratively solving the optimization function through a multi-objective optimization algorithm.

[0032] The optimal intervention time is sorted according to the order of each high-risk connected area on the critical propagation path to obtain the support timing sequence. The optimal strength configuration is matched with the preset support type library. Support structure and support parameters are assigned to each high-risk connected area to obtain the support type scheme. The support timing sequence and the support type scheme are combined according to the spatial location of the high-risk connected areas to obtain the zonal support decision result.

[0033] In one alternative implementation,

[0034] Based on the zonal support decision results, support operations are implemented on the surrounding rock, and the response data of the surrounding rock after support is collected. Based on the surrounding rock response data, the suppression effect of the support operation is determined, and based on the suppression effect, Bayesian updates are performed on the embedded pore water phase change constraint and strength degradation constraint to obtain optimized constraint parameters, including:

[0035] Based on the zonal support decision results, support operations are carried out in each high-risk connected area, and surrounding rock displacement, temperature, and stress are collected after support is completed to obtain surrounding rock response data. Based on the surrounding rock response data, the displacement difference and stress difference before and after support are calculated to obtain the displacement suppression amount and stress release amount. The displacement suppression amount and stress release amount are weighted and summed to obtain the suppression effect. The surrounding rock temperature is compared with the phase change threshold in the embedded pore water phase change constraint, and the proportion of positions where the surrounding rock temperature is lower than the phase change threshold is statistically analyzed.

[0036] The phase transition threshold and the intensity degradation constraint are adjusted according to the relationship between the suppression effect and the position ratio to obtain the updated phase transition threshold and the updated degradation coefficient;

[0037] Based on the updated phase transition threshold and Bayesian inference algorithm, the temperature change rate at the phase transition position in the surrounding rock temperature is fitted with the latent heat parameter in the embedded pore water phase transition constraint. The updated latent heat parameter is obtained by combining the least squares method and the embedded pore water phase transition constraint is updated. Based on the updated degradation coefficient and Bayesian inference algorithm, the maximum value of the attenuation ratio of the surrounding rock stress relative to the initial stress of the surrounding rock when no support operation is performed is calculated. When the maximum value of the attenuation ratio is greater than the preset degradation upper limit, the degradation upper limit is adjusted upward and the strength degradation constraint is updated.

[0038] The updated embedded pore water phase change constraint and the updated strength degradation constraint are combined to obtain the optimized constraint parameters.

[0039] In one alternative implementation,

[0040] The process of collecting real-time monitoring data, solving for instability characteristic parameters based on the optimization constraint parameters, and determining the corresponding optimal support decision scheme includes:

[0041] Real-time monitoring data is obtained by extracting the surrounding rock temperature, displacement, and stress from sensors in the tunnel surrounding rock. An updated phase transition threshold and an updated degradation coefficient are extracted from the optimized constraint parameters. The displacement change rate is calculated based on the real-time monitoring data, and the degradation rate is obtained by correlating the displacement change rate with the updated degradation coefficient. The unstable region is identified in the real-time monitoring data where the surrounding rock temperature is lower than the updated phase transition threshold and the degradation rate is greater than a preset rate threshold.

[0042] The stress concentration index is obtained by calculating the spatial and temporal gradients of the surrounding rock stress corresponding to the unstable region, and the instability characteristic parameters are obtained by combining the deterioration rate and the unstable region.

[0043] Based on the instability characteristic parameters, the support target is determined. Based on the support target, the support type is selected from the preset support type library and the support time is determined. The updated embedded pore water phase change constraint and the updated strength degradation constraint are used as constraints. The optimal support decision scheme is obtained by minimizing the instability characteristic parameters.

[0044] A second aspect of this invention provides an intelligent decision-making system for tunnel support in cold regions based on dynamic feedback, comprising:

[0045] The decoupled prediction unit is used to acquire historical monitoring data of the surrounding rock of tunnels in cold regions and perform time alignment and spatial mapping to obtain a monitoring dataset. It performs feature decoupling on the historical monitoring data to obtain freeze-thaw damage features and deformation field components. Based on the freeze-thaw damage features and pre-set embedded pore water phase transformation constraints and strength degradation constraints, it performs coupled calculations to predict the evolution trajectory of the surrounding rock damage field and obtain the evolution prediction result.

[0046] The graph construction unit is used to perform graph topology mining on the evolution prediction results, construct an evolution graph with monitoring points as nodes and damage propagation as edges, identify high-risk connected regions and critical propagation paths, calculate the stress concentration degree of the connected regions in combination with the deformation field components, and construct a risk assessment graph.

[0047] The support optimization unit is used to obtain the support timing sequence and support type scheme corresponding to each connected area based on the pre-acquired geological weak surface information and existing support status, combined with the risk assessment map, through a multi-objective optimization algorithm, and thus obtain the zonal support decision results;

[0048] The feedback update unit is used to implement support operations on the surrounding rock based on the partitioned support decision results and collect the surrounding rock response data after support. Based on the surrounding rock response data, it determines the suppression effect of the support operation and performs Bayesian update on the embedded pore water phase change constraint and strength degradation constraint based on the suppression effect to obtain optimized constraint parameters. It also collects real-time monitoring data and solves the instability characteristic parameters based on the optimized constraint parameters to determine the corresponding optimal support decision scheme.

[0049] A third aspect of the present invention provides an electronic device, comprising:

[0050] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0051] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0052] In this invention, by decoupling historical monitoring data, freeze-thaw damage characteristics and deformation field components are separated. Coupled calculations are then performed using embedded pore water phase transition constraints and strength degradation constraints, enabling accurate prediction of the evolution trajectory of the damage field in the surrounding rock of tunnels in cold regions. This improves the understanding and prediction accuracy of the deformation and failure mechanism of surrounding rock under special geological conditions in cold regions. An evolution map is constructed using graph topology mining to identify high-risk connected regions and critical propagation paths. Stress concentration is calculated by combining deformation field components, and a risk assessment map is constructed, achieving a refined assessment and visual representation of the risk of damage to the surrounding rock of tunnels in cold regions. Based on geological weak surface information, existing support status, and the risk assessment map, a multi-objective optimization algorithm determines the support timing sequence and support type scheme for each connected region, achieving precise allocation of support resources and targeted construction, avoiding blindness and resource waste. The suppression effect is evaluated by collecting surrounding rock response data after support, and Bayesian updates are performed on the constraints, achieving adaptive optimization of support parameters. This allows the decision-making scheme to be dynamically adjusted with construction progress and environmental changes, improving the safety and economy of tunnel construction in cold regions. Attached Figure Description

[0053] Figure 1 This is a flowchart illustrating the intelligent decision-making method for tunnel support in cold regions based on dynamic feedback, as described in an embodiment of the present invention.

[0054] Figure 2 This is a flowchart illustrating the optimization of freeze-thaw constraint parameters driven by surrounding rock response data in an embodiment of the present invention, based on a dynamic feedback-based intelligent decision-making method for tunnel support in cold regions. Detailed Implementation

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

[0056] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0057] Figure 1This is a flowchart illustrating the intelligent decision-making method for tunnel support in cold regions based on dynamic feedback, as described in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0058] Historical monitoring data of the surrounding rock of tunnels in cold regions is acquired and time-aligned and spatially mapped to obtain a monitoring dataset. Feature decoupling is performed on the historical monitoring data to obtain freeze-thaw damage features and deformation field components. Based on the freeze-thaw damage features and pre-set embedded pore water phase transformation constraints and strength degradation constraints, coupled calculations are performed to predict the evolution trajectory of the surrounding rock damage field and obtain evolution prediction results.

[0059] Graph topology mining is performed on the evolution prediction results to construct an evolution graph with monitoring points as nodes and damage propagation as edges, and high-risk connected regions and critical propagation paths are identified. The stress concentration degree of the connected regions is calculated in combination with the deformation field components, and a risk assessment map is constructed.

[0060] Based on pre-acquired geological weak surface information and existing support status, and combined with the risk assessment map, a multi-objective optimization algorithm is used to obtain the support timing sequence and support type scheme corresponding to each connected area, thus obtaining the zonal support decision results;

[0061] Based on the zonal support decision results, support operations are performed on the surrounding rock and the response data of the surrounding rock after support is collected. Based on the surrounding rock response data, the suppression effect of the support operation is determined, and based on the suppression effect, Bayesian updates are performed on the embedded pore water phase change constraint and strength degradation constraint to obtain optimized constraint parameters. Real-time monitoring data is collected, and the instability characteristic parameters are solved based on the optimized constraint parameters to determine the corresponding optimal support decision scheme.

[0062] In one alternative implementation,

[0063] Historical monitoring data of the surrounding rock of tunnels in cold regions is acquired and subjected to time alignment and spatial mapping to obtain a monitoring dataset. Feature decoupling is then performed on the historical monitoring data to obtain freeze-thaw damage features and deformation field components, including:

[0064] Raw data collected by various sensors located at different cross sections and burial depths over a historical period are extracted and combined with the timestamps and spatial coordinates of the data transmission to construct historical monitoring data.

[0065] Cross-correlation analysis is performed on the timestamps of each sensor in the historical monitoring data to calculate the response delay difference of different sensors. Based on the response delay difference, the time offset of each sensor is determined and offset compensation is performed to obtain corrected data. Sensor channels with inconsistent sampling frequencies in the corrected data are resampled to obtain unified sampling data.

[0066] The unified sampling data is combined with the spatial coordinates of the corresponding sensors to divide the data into three-dimensional grid nodes. The Euclidean distance from each sensor to each grid node is calculated, and the mapping weight is determined based on the Euclidean distance. The temperature field value, displacement field value, and strain field value of each grid node are solved based on the mapping weight and the unified sampling data and combined to obtain the monitoring dataset.

[0067] The temperature field values ​​in the monitoring dataset are decomposed into time-frequency components to extract periodic fluctuation components and calculate the time derivative and spatial gradient of the temperature field values. Based on the time derivative and spatial gradient, the freeze-thaw interface is identified. The freeze-thaw damage characteristics are obtained by statistical analysis of the temperature change amplitude, time derivative peak value, and spatial gradient extreme value at the freeze-thaw interface. The deformation field components are obtained by spatial differentiation of the displacement field values ​​and strain field values ​​and extraction of the time evolution trend.

[0068] Raw monitoring data was extracted from various sensors deployed within tunnels in cold regions. These sensors included temperature sensors, displacement sensors, and strain sensors, which were positioned at different cross-sections and depths within the tunnel. For example, in a highway tunnel project in a cold region, temperature sensors were deployed at 5-meter intervals on the outer surface of the tunnel lining in contact with the surrounding rock; displacement sensors were deployed at 10-meter intervals on the tunnel arch, arch waist, and bottom; and strain sensors were deployed on the secondary lining reinforcement at 8-meter intervals. Raw data collected by the sensors over a continuous 180-day period was extracted. Each sensor's data included a timestamp of the acquisition and corresponding spatial coordinates, thus constructing a historical monitoring database.

[0069] Because the acquisition clocks of each sensor may have deviations, cross-correlation analysis is performed on the timestamps of each sensor in the historical monitoring data. A reference sensor is selected, and the cross-correlation function between the responses of other sensors and the reference sensor is calculated. The time corresponding to the maximum value of the cross-correlation function is the response delay difference between the two sensors. For example, selecting the temperature sensor at the tunnel entrance as a reference, the calculated response delay between the displacement sensor and the reference sensor is 15 minutes, and the response delay between the strain sensor and the reference sensor is 8 minutes. Based on the response delay difference, a time offset is determined for each sensor and compensated for to obtain corrected data.

[0070] The calibrated data suffers from inconsistent sampling frequencies among different sensor types. For example, temperature sensors sample every 30 minutes, displacement sensors every 10 minutes, and strain sensors every 15 minutes. Therefore, the sensor channels with inconsistent sampling frequencies are resampled. An interpolation algorithm is used to unify all sensor data to the same sampling frequency, such as every 5 minutes, resulting in uniform sampled data.

[0071] The unified sampling data is combined with the spatial coordinates of the corresponding sensors to divide the tunnel space into three-dimensional grid nodes. The grid density is determined according to the monitoring accuracy requirements. In this embodiment, the longitudinal spacing along the tunnel is 2 meters, and the horizontal and vertical spacing is 0.5 meters. For each grid node, the Euclidean distance from the grid node to each sensor is calculated. The Euclidean distance is calculated directly from the distance between the sensor coordinates and the grid node coordinates in three-dimensional space. For example, if the coordinates of a temperature sensor are (120.5, 35.2, 15.8) and the coordinates of a grid node are (121.0, 35.5, 16.0), then the Euclidean distance between the two points is 0.67 meters.

[0072] Based on the calculated Euclidean distance, a mapping weight is determined for each grid node. The mapping weight is inversely proportional to the Euclidean distance; the closer the sensor is, the greater its influence on the grid node. In the actual calculation, an inverse distance weighting method is used, where the weight of a sensor at a distance of 'd' is 1 divided by the square of 'd'. When the distance between the grid node and the sensor is less than a threshold, such as 0.1 meters, the current sensor value is directly used. By weighting the mapping weight with the uniformly sampled data, the temperature field value, displacement field value, and strain field value of each grid node are solved, and combined to obtain a complete monitoring dataset.

[0073] For the acquired monitoring dataset, the temperature field values ​​are decomposed using time-frequency analysis to extract periodic fluctuation components such as diurnal and seasonal variations. For example, after decomposition, 180 days of temperature data for a certain grid node can identify a diurnal periodic fluctuation with an amplitude of 4.5℃ and a seasonal periodic fluctuation with an amplitude of 18.2℃. Based on the decomposed temperature field values, their time derivative and spatial gradient are calculated. The time derivative represents the rate of temperature change over time, calculated by dividing the temperature difference between adjacent time points by the time interval; the spatial gradient represents the rate of temperature change in space, calculated by dividing the temperature difference between adjacent grid nodes by the spatial distance.

[0074] The location of freeze-thaw interfaces was identified using time derivatives and spatial gradients. Freeze-thaw interfaces are characterized by temperature fluctuations around 0°C, a significant peak in the time derivative, and extreme values ​​in the spatial gradient. In this embodiment, the identified freeze-thaw interfaces were mainly distributed in the area 15 to 25 cm outside the tunnel arch and 10 to 20 cm outside the sidewalls. Statistical analysis of the temperature variation amplitude, peak time derivative, and extreme values ​​of the spatial gradient at the freeze-thaw interfaces yielded freeze-thaw damage characteristics. For example, the average temperature variation amplitude at the freeze-thaw interface in the arch area was 2.8°C, the average peak time derivative was 0.15°C per hour, and the average extreme value of the spatial gradient was 4.2°C per meter.

[0075] Simultaneously, spatial differentiation operations are performed on the displacement and strain field values ​​to calculate the components of the deformation tensor. By extracting the temporal evolution trends corresponding to the components of the deformation tensor, a complete description of the deformation field components is obtained. For example, in the tunnel arch region, the radial displacement field shows a slow increasing trend over 180 days, with a growth rate of 0.05 mm per day; the tangential strain field exhibits periodic fluctuations with an amplitude of 50 microstrains.

[0076] In this embodiment, by performing cross-correlation analysis on the timestamps of each sensor and implementing response delay compensation, the time offset problem caused by differences in communication and sampling mechanisms between different acquisition channels is effectively eliminated, improving the temporal consistency and reliability of temperature, displacement, and strain changes. By uniformly resampling sensor data with inconsistent sampling frequencies, the synchronous expression of different physical quantities at the same time scale is ensured, providing a stable data foundation for subsequent multi-field coupling analysis. By introducing sensor spatial coordinates and constructing three-dimensional grid nodes, and determining the mapping weight by combining the spatial distance from the sensor to the grid node, a smooth mapping from discrete measurement points to a continuous spatial field is achieved, significantly enhancing the spatial integrity and resolution of the monitoring results, enabling the temperature field, displacement field, and strain field to truly reflect the spatial distribution characteristics inside the engineering body.

[0077] In one alternative implementation,

[0078] Based on the freeze-thaw damage characteristics and pre-set embedded pore water phase transformation constraints and strength degradation constraints, coupled calculations are performed to predict the evolution trajectory of the surrounding rock damage field. The evolution prediction results include:

[0079] The freeze-thaw interface in the freeze-thaw damage features is used as the phase change trigger boundary. The temperature change amplitude and time derivative peak value in the freeze-thaw damage features are extracted as phase change parameters. Based on the phase change trigger boundary and the phase change parameters, the pore water volume expansion rate is calculated in the pre-set embedded pore water phase change constraint and converted into an expansion stress field.

[0080] The stress concentration location is determined based on the spatial gradient extrema in the expansion stress field and freeze-thaw damage characteristics. The strain field value history sequence of the corresponding grid node in the monitoring dataset is extracted at the stress concentration location and accumulated over time to obtain the cumulative strain. Based on the cumulative strain and the pre-set strength degradation constraint, the elastic modulus attenuation coefficient and compressive strength reduction coefficient are calculated and coupled with the expansion stress field to obtain the equivalent stress field.

[0081] Using the equivalent stress field as the initial stress, and combining the freeze-thaw damage characteristics to set the evolution time step, stress redistribution is performed on the equivalent stress field within each time step, and damage points are determined. The damage points are connected and tracked in the spatial domain, and the occurrence order and position changes of each damage point are recorded and combined to obtain the evolution trajectory. The spatial distribution characteristics and temporal evolution characteristics of the damage field corresponding to the evolution trajectory are extracted and summarized to obtain the evolution prediction results.

[0082] The freeze-thaw interface in the identified freeze-thaw damage features is used as the phase transition trigger boundary. The temperature change amplitude and peak time derivative are extracted from the freeze-thaw damage features as phase transition parameters to describe the severity of the phase transition process. For example, in a tunnel project in a cold region, the temperature change amplitude at the freeze-thaw interface in the arch region is 2.8℃, and the peak time derivative is 0.15℃ per hour. The phase transition parameters directly affect the phase transition rate. Based on the phase transition trigger boundary and phase transition parameters, the pore water volume expansion rate is calculated within a pre-set embedded pore water phase transition constraint. The embedded pore water phase transition constraint considers the characteristic of approximately 9% volume expansion during the ice-water phase transition, while also taking into account factors such as the surrounding rock porosity and saturation. For example, for the case where the surrounding rock porosity in the arch region is 0.23 and the saturation is 0.85, the calculated volume expansion rate at the freeze-thaw interface is 7.2%. When converting the volumetric expansion rate into an expansion stress field, considering the surrounding rock constraints and deformation characteristics, the expansion stress is calculated through the stress-strain relationship, with typical values ​​ranging from 0.42 MPa to 1.25 MPa.

[0083] Based on the calculated expansion stress field, the stress concentration locations are determined by combining the spatial gradient extrema in the freeze-thaw damage characteristics. These spatial gradient extrema reflect the locations of the most drastic temperature changes and are typically highly correlated with stress concentration areas. In this embodiment, the region with a spatial gradient extrema of 4.2℃ per meter corresponds to stress concentration locations 20 cm outside the arch crown and 15 cm outside the sidewall. At these stress concentration locations, the historical strain field values ​​of the corresponding grid nodes in the monitoring dataset are extracted. The strain data from 180 consecutive days are accumulated to obtain the cumulative strain. For example, the cumulative strain at the arch crown stress concentration location reaches 850 microstrains, and the cumulative strain at the sidewall stress concentration location is 720 microstrains. Based on the cumulative strain and pre-set strength degradation constraints, the elastic modulus attenuation coefficient and compressive strength reduction coefficient are calculated. The strength degradation constraint considers the damage effect of freeze-thaw cycles on the rock mass strength. A degradation function is established with cumulative strain as the independent variable. Calculations show that the elastic modulus at the crown is 0.82 and the compressive strength reduction factor is 0.78; the elastic modulus at the sidewalls is 0.85 and the compressive strength reduction factor is 0.83. The degradation coefficients are coupled with the expansion stress field to obtain the equivalent stress field. This coupling calculation considers the combined effect of the decrease in bearing capacity due to strength degradation and the increase in expansion stress. The calculated equivalent stress at the crown is 1.32 times the original stress, and at the sidewalls it is 1.25 times.

[0084] The aforementioned equivalent stress field is used as the initial stress, and the evolution time step is set in conjunction with the characteristics of freeze-thaw damage. The time step is selected based on the rate of temperature change; in this embodiment, it is set to 12 hours to ensure that key changes during the freeze-thaw process can be captured. Stress redistribution calculations are performed on the equivalent stress field within each time step. Based on the principles of elasticity, when the stress at a certain point exceeds the strength threshold, the excess stress is redistributed to the surrounding area. The stress distribution is calculated iteratively until an equilibrium state is reached. During the stress redistribution process, when the stress at a certain point exceeds the material strength at that point, it is marked as a damage point. For example, in the arch region, as the freeze-thaw cycle proceeds, the first damage point appears at the 5th time step, located 20 cm outside the arch; at the 8th time step, the damage area expands to 15 to 25 cm outside the arch; at the 15th time step, the damage area further expands to 10 to 30 cm outside the arch. Damage points are connected and tracked in the spatial domain, recording the order of occurrence and positional changes of each damage point, and combined to obtain the damage evolution trajectory.

[0085] Spatial distribution characteristics of the damage field were extracted from the evolution trajectory, revealing an elliptical distribution of damaged areas with the major axis aligned with the temperature gradient direction, primarily concentrated near the freeze-thaw interface of the arch and sidewalls. Damage depth increased with the number of freeze-thaw cycles, reaching 35 cm in the arch and 25 cm in the sidewalls after 30 cycles. Temporal evolution characteristics of the damage field showed a non-linear development pattern, initially slow, with a significantly increased rate of damage area growth after 10 freeze-thaw cycles, and connectivity of damaged areas appearing after 20 cycles. Combining spatial distribution and temporal evolution characteristics yielded a complete evolution prediction result. The prediction indicates that without support measures, a through-crack will form in the arch area after 40 freeze-thaw cycles, potentially leading to localized rockfall; through-cracks will also appear in the sidewalls after 60 cycles, severely threatening the overall stability of the tunnel.

[0086] In this embodiment, the phase transition expansion behavior of pore water is constrained by using the freeze-thaw interface as the phase transition trigger boundary and introducing the temperature change amplitude and time derivative peak value as phase transition parameters. This allows the calculation of expansion stress to be directly controlled by the actual freeze-thaw process, improving the accuracy and sensitivity of the expansion stress field to the real freeze-thaw environment. By combining the expansion stress field with the spatial gradient extremum to determine the stress concentration location and performing time accumulation analysis based on the historical strain field, a quantitative characterization of the latent damage accumulation effect under repeated freeze-thaw action is achieved. This can identify potential high-risk areas and accurately reflect the gradual process of material internal performance degradation. By introducing strength degradation constraints and calculating the elastic modulus attenuation coefficient and compressive strength reduction coefficient, the degradation process of material mechanical parameters is coupled with the expansion stress induced by freeze-thaw, overcoming the problem of the separation between stress analysis and damage evolution, and significantly enhancing the engineering applicability of the damage assessment results.

[0087] In one alternative implementation,

[0088] Graph topology mining is performed on the evolution prediction results to construct an evolution graph with monitoring points as nodes and damage propagation as edges. High-risk connected regions and critical propagation paths are identified. The stress concentration degree of the connected regions is calculated by combining the deformation field components, and a risk assessment map is constructed, including:

[0089] Based on the spatial distribution characteristics of the damage field in the evolution prediction results, the spatial coordinates of each damage point are extracted and spatially matched with the grid nodes in the monitoring dataset. According to the matching results, the damage points are mapped to the nearest grid node and marked as graph nodes.

[0090] The position change vector of each damage point between adjacent time steps is obtained and the corresponding orientation angle and length are calculated. Damage points with the same orientation angle and continuous length are identified and combined to obtain a propagation chain. Graph nodes corresponding to adjacent damage points in the propagation chain are connected and directed edges are created. The weight of the directed edges is set based on the time interval between adjacent damage points to obtain an evolution graph. The directed edges in the evolution graph are traversed and fast edges are determined in combination with a pre-set weight threshold. The connectivity of the fast edges is detected and the set of nodes reachable through the fast edges is marked as a high-risk connected region.

[0091] Traverse all paths from the initial damage point to the boundary node of the high-risk connected region in the evolution graph and determine the velocity index based on the sum of the weights of all directed edges on the path. The path with the largest velocity index is taken as the critical propagation path.

[0092] The high-risk connected regions are matched with the deformation field components in the spatial dimension, and the strain rate and strain gradient at the corresponding positions are extracted. The spatial second derivatives of the strain rate and strain gradient are calculated and the stress concentration is determined. The risk assessment map is obtained by combining the high-risk connected regions and the evolution map.

[0093] Based on the spatial distribution characteristics of the damage field in the evolutionary prediction results, the spatial coordinates of each damage point are extracted. For each monitored damage point, the corresponding three-dimensional spatial coordinates are recorded. For example, in the arch area of ​​a tunnel in a cold region, the spatial coordinates of the first damage point appearing at the 5th time step are (52.3, 10.5, 8.2), and the coordinates of the damage point appearing at the 8th time step are (52.5, 10.6, 8.0). The damage point coordinates are spatially matched with the grid nodes in the monitoring dataset. The Euclidean distance from each damage point to all grid nodes is calculated, and the grid node with the smallest distance is selected as the mapping node for the current damage point. The aforementioned two damage points are mapped to grid nodes numbered N1258 and N1265, respectively. The mapped grid nodes are marked as graph nodes, serving as the basic units for constructing the evolutionary graph. In this embodiment, a total of 128 damage points are extracted and mapped to 112 graph nodes.

[0094] The position change vector of each damage point between adjacent time steps is obtained by calculating the spatial coordinate difference between damage points appearing in adjacent time steps. Taking the two aforementioned damage points as an example, the position change vector is (0.2, 0.1, -0.2). The orientation angle and length of the position change vector are calculated. The orientation angle is determined by the angle between the vector and the coordinate axis, which is (42°, 27°, -53°), and the vector length is 0.3. By comparing the orientation angle and length of all damage points, damage points with the same orientation angle and continuous length are identified. The orientation angle similarity threshold is set to ±15°, and the length continuity criterion is that the length difference between adjacent damage points does not exceed 25%. Damage points that meet the conditions are combined into a propagation chain. For example, a propagation chain containing 7 damage points is identified in the arch area, with an orientation angle of approximately (45°, 25°, -50°) and a length increasing from 0.3 to 0.8.

[0095] Directed edges are created by connecting graph nodes corresponding to adjacent damage points in the propagation chain. The direction of the directed edges points from the earlier damage point to the later damage point, representing the direction of damage propagation. Weights are assigned to the directed edges based on the time interval between adjacent damage points; shorter time intervals indicate faster damage propagation and thus a larger weight. The weight is calculated by multiplying the reciprocal of the time interval by a standardization factor of 100. For example, a time interval of 3 time steps between two damage points corresponds to an edge weight of 33.3; a time interval of 1 time step corresponds to an edge weight of 100. A complete evolutionary graph is constructed using all graph nodes and directed edges, containing 112 nodes and 156 directed edges.

[0096] All directed edges in the evolution graph are traversed, and fast edges are determined by combining them with a pre-set weight threshold. The weight threshold is set to 75, meaning that edges with a time interval of less than or equal to 1.33 time steps are marked as fast edges. In this embodiment, 42 fast edges are identified. Connectivity detection is performed on the fast edges using a depth-first search algorithm, starting from each initial damage point, to find all nodes reachable through the fast edges. The set of nodes reachable through fast edges is marked as a high-risk connected region. A high-risk connected region containing 23 nodes is identified in the vault region, and two high-risk connected regions are identified in the sidewall region, containing 15 and 12 nodes respectively.

[0097] The evolution graph is traversed, exploring all possible paths from the initial damage point to the boundary node of the high-risk connected region. The initial damage point is the earliest damage point, and the boundary node of the high-risk connected region is the node connected to nodes outside the high-risk connected region. An improved Dijkstra algorithm is used to calculate all paths, resulting in 28 valid paths. A velocity index is determined based on the sum of the weights of all directed edges on the path; a larger sum of weights indicates a faster propagation speed. In the vault region, the path with a sum of weights of 582 is identified as the critical propagation path. This critical propagation path contains 8 nodes and reaches the boundary of the high-risk connected region from the initial damage point through the key nodes. The critical propagation paths corresponding to the two high-risk connected regions in the sidewall region have sums of weights of 468 and 425, respectively.

[0098] High-risk connected regions were spatially matched with deformation field components to extract strain rate and strain gradient at corresponding locations. Strain rate represents the strain change per unit time, calculated by dividing the strain difference between adjacent time steps by the time interval. Strain gradient represents the rate of strain change in space, calculated by dividing the strain difference between adjacent grid nodes by the spatial distance. In the high-risk connected region of the arch crown, the average strain rate was 12.5 microstrains per hour, and the average strain gradient was 35 microstrains per meter. The spatial second derivatives of the strain rate and strain gradient were calculated. The second derivative reflects the acceleration of strain change and is obtained by dividing the gradient difference between adjacent points by the distance. The peak value of the second derivative in the arch crown region was 120 microstrains per square meter, located in the central part of the high-risk connected region, indicating the highest stress concentration at this location, which is the potential starting point for failure.

[0099] A risk assessment map was constructed by combining high-risk connected regions and an evolution diagram. The map is based on grid nodes, with color depth representing the risk level and lines indicating potential propagation paths. Risk levels were categorized into five levels according to stress concentration. Level 1 risk areas require immediate support and reinforcement, Level 2 risk areas require enhanced monitoring and appropriate reinforcement, and Level 3 and below risk areas maintain routine monitoring. In this embodiment, 6 nodes in the high-risk connected region of the arch were marked as Level 1 risk, and 12 nodes were marked as Level 2 risk; in the two high-risk connected regions of the sidewalls, 3 and 2 nodes were marked as Level 1 risk, and 8 and 7 nodes were marked as Level 2 risk, respectively.

[0100] In this embodiment, by spatially matching damage points with monitoring grid nodes and constructing graph nodes, a unified mapping relationship is formed between discrete damage information and the original three-dimensional monitoring data system. This avoids the problem that damage results are difficult to correspond to actual monitoring locations, and improves the interpretability and feasibility of damage location results in engineering applications. By introducing the position change vector of damage points between adjacent time steps, automatic identification of the directionality and continuity of damage propagation is achieved, which can accurately depict the real propagation pattern of damage from point to line and from local to global. By constructing a directed evolution graph containing time weights and identifying fast edges and their connected regions based on weight thresholds, quantitative discrimination of rapid damage expansion channels is achieved, avoiding the problem of difficulty in distinguishing between slow degradation and sudden damage. This significantly improves the sensitivity and reliability of identifying high-risk connected regions. By traversing all propagation paths in the evolution graph and introducing velocity indicators to screen critical propagation paths, risk analysis is upgraded from static spatial distribution to dynamic assessment based on propagation velocity and path characteristics. This can identify the key propagation paths most likely to cause structural instability in advance, enhancing the foresight of risk warning.

[0101] In one alternative implementation,

[0102] Based on pre-acquired geological weak surface information and existing support status, and combined with the risk assessment map, a multi-objective optimization algorithm is used to obtain the support timing sequence and support type scheme corresponding to each connected region, resulting in the following zonal support decision results:

[0103] The fault strike and joint density in the pre-acquired geological weak surface information are superimposed and matched with the high-risk connected regions in the risk assessment map to identify the coupling zone and extract the weakening coefficient.

[0104] Extract the spatial location and strength grade of the constructed support from the pre-acquired existing support status, calculate the shortest distance between the coupling zone and the spatial location of the constructed support to obtain the coverage, and calculate the gap amount by the difference between the strength grade and the pre-acquired stress concentration degree.

[0105] Using the weakening coefficient, coverage, and gap amount as constraints, the spatial orientation of the pre-obtained critical propagation path as boundary conditions, and the stress concentration degree of each high-risk connected area as risk weight, an optimization function is constructed with the objectives of minimizing the overall risk exposure time and maximizing the support resource utilization rate. The optimal intervention time and optimal strength configuration are obtained by iteratively solving the optimization function through a multi-objective optimization algorithm.

[0106] The optimal intervention time is sorted according to the order of each high-risk connected area on the critical propagation path to obtain the support timing sequence. The optimal strength configuration is matched with the preset support type library. Support structure and support parameters are assigned to each high-risk connected area to obtain the support type scheme. The support timing sequence and the support type scheme are combined according to the spatial location of the high-risk connected areas to obtain the zonal support decision result.

[0107] The pre-acquired geological weak surface information is overlaid and matched with high-risk connected areas in the risk assessment map. The geological weak surface information includes key parameters such as fault strike, joint density, and rock mass integrity. Taking a tunnel in a cold region as an example, the pre-acquired geological weak surface information shows that there is a fault in the tunnel arch area with a strike of 60 degrees east and a dip angle of 35 degrees, and a joint density of 3.2 joints per meter; there are two sets of joints on the left wall, with the main joint set striking 45 degrees west and a joint density of 2.8 joints per meter; the joint density on the right wall is relatively low, at 1.5 joints per meter. This information is spatially overlaid with the high-risk connected areas in the risk assessment map, and a distance threshold method is used to identify coupling zones. The distance threshold is set to 0.5 meters, meaning that when the spatial distance between a geological weak surface and a high-risk connected area is less than 0.5 meters, it is determined to be a coupling zone. Spatial analysis revealed that eight nodes in the high-risk connected region of the arch crown were coupled with the fault, five nodes in the high-risk connected region of the left side wall were coupled with the main joint group, and only two nodes in the high-risk connected region of the right side wall were coupled with the joint.

[0108] The weakening coefficient is extracted from the coupling zone. This coefficient reflects the degree to which the geological weak surface weakens the stability of the surrounding rock. The calculation of the weakening coefficient comprehensively considers factors such as the aperture of the fault or joint, the nature of the infill material, and its continuity. For the crown fault coupling zone, due to the large fault aperture (average 3.5 mm) and the predominantly soft clayey material, the weakening coefficient is 0.68. In the left side joint group coupling zone, with a smaller joint aperture (average 1.2 mm), the weakening coefficient is 0.82. The weakening coefficient for the right side coupling zone is 0.91. A smaller weakening coefficient indicates a greater degree of weakening of the surrounding rock stability by the geological weak surface, requiring stronger support measures.

[0109] The spatial location and strength grade of the existing support structure are extracted from the pre-obtained existing support status. The existing support status is obtained through tunnel construction records and on-site inspections, including parameters such as support structure type, thickness, and strength. In this embodiment, the existing tunnel support includes C25 shotcrete in the crown area, with a thickness of 6 cm and 2-meter long system anchor bolts spaced 1 meter by 1 meter, with a strength grade of 3; the side walls use C25 shotcrete, with a thickness of 5 cm and a strength grade of 2. The coverage is calculated by determining the shortest distance between the coupling zone and the spatial location of the existing support. Coverage represents the degree to which the existing support covers high-risk areas. For the crown coupling zone, the coverage of the existing support is 85%; the coverage of the left side coupling zone is 92%; and the coverage of the right side coupling zone is 95%. Higher coverage indicates more comprehensive protection of high-risk areas by the existing support.

[0110] The gap is calculated by subtracting the strength grade from the pre-obtained stress concentration level. The gap represents the difference between the existing support strength and the actual requirements. The stress concentration level in the crown area is 4.5, which differs from the existing support strength grade of 3 by 1.5, resulting in a gap of 1.5. The stress concentration level on the left wall is 3.8, with a gap of 1.8; the stress concentration level on the right wall is 2.5, with a gap of 0.5. A larger gap indicates insufficient existing support strength, requiring more timely reinforcement measures.

[0111] Using the weakening coefficient, coverage, and gap amount as constraints, the spatial orientation of the critical propagation path as boundary conditions, and the stress concentration degree of each high-risk connected area as risk weight, an optimization function is constructed. The optimization function aims to minimize the overall risk exposure duration and maximize the support resource utilization rate. The risk exposure duration is defined as the time interval from risk identification to support implementation, and the support resource utilization rate is defined as the degree of matching between support strength and actual needs. The optimization function adopts a weighted summation form, with a risk exposure duration weight of 0.7 and a support resource utilization rate weight of 0.3. The constraints are set as follows: areas with a weakening coefficient below 0.75 must be supported within 48 hours of identification; areas with coverage below 90% require supplementary support; and areas with a gap amount greater than 1 require increased support strength.

[0112] The optimal intervention time and intensity configuration were obtained through iterative solutions using a multi-objective optimization algorithm. An improved particle swarm optimization algorithm was employed, with a population size of 100, a maximum number of iterations of 200, and a convergence threshold of 0.001. After 136 iterations, the algorithm converged, yielding the optimal intervention time of 24 hours and the optimal intensity configuration of level 5 for the arch coupling zone; 36 hours and level 4 for the left side coupling zone; and 72 hours and level 3 for the right side coupling zone.

[0113] The optimal intervention time was ordered according to the sequence of high-risk connected areas along the critical propagation path to obtain the support timing sequence. Critical propagation path analysis showed that the damage propagation sequence was from the crown area to the left wall area to the right wall area; therefore, the support timing sequence was: 24 hours for the crown area, 36 hours for the left wall area, and 72 hours for the right wall area. The optimal strength configuration was matched with a pre-set support type library to assign support structures and parameters to each high-risk connected area. The pre-set support type library contains various support types, such as shotcrete, anchor bolts, steel arch frames, reinforced steel grids, and waterproof membranes, as well as different parameter combinations for each type.

[0114] In this embodiment, by overlaying and matching geological weak surface information such as fault strike and joint density with high-risk connected areas, the coupling zone between damage propagation and natural structural weakening is identified. This allows support decisions to explicitly consider the amplification effect of geological structures on damage expansion, avoiding the risk underestimation caused by ignoring geological heterogeneity. This improves the reliability of high-risk area identification and subsequent support design. By quantifying the coverage between the coupling zone and the existing support spatial location, and combining gap analysis of support strength level and stress concentration degree, an objective assessment of the effectiveness of existing support is achieved. This can accurately identify support blind spots and areas with insufficient strength, providing a clear basis for support adjustment. Through iterative solution of multi-objective optimization algorithm, the optimal intervention time and optimal strength configuration are obtained, realizing the transformation of support measures from ex-post reinforcement to early intervention based on risk evolution. Intervention can be implemented in the critical time window before rapid damage propagation, effectively shortening the high-risk exposure time and reducing the probability of sudden instability.

[0115] In one alternative implementation,

[0116] Based on the zonal support decision results, support operations are implemented on the surrounding rock, and the response data of the surrounding rock after support is collected. Based on the surrounding rock response data, the suppression effect of the support operation is determined, and based on the suppression effect, Bayesian updates are performed on the embedded pore water phase change constraint and strength degradation constraint to obtain optimized constraint parameters, including:

[0117] Based on the zonal support decision results, support operations are carried out in each high-risk connected area, and surrounding rock displacement, temperature, and stress are collected after support is completed to obtain surrounding rock response data. Based on the surrounding rock response data, the displacement difference and stress difference before and after support are calculated to obtain the displacement suppression amount and stress release amount. The displacement suppression amount and stress release amount are weighted and summed to obtain the suppression effect. The surrounding rock temperature is compared with the phase change threshold in the embedded pore water phase change constraint, and the proportion of positions where the surrounding rock temperature is lower than the phase change threshold is statistically analyzed.

[0118] The phase transition threshold and the intensity degradation constraint are adjusted according to the relationship between the suppression effect and the position ratio to obtain the updated phase transition threshold and the updated degradation coefficient;

[0119] Based on the updated phase transition threshold and Bayesian inference algorithm, the temperature change rate at the phase transition position in the surrounding rock temperature is fitted with the latent heat parameter in the embedded pore water phase transition constraint. The updated latent heat parameter is obtained by combining the least squares method and the embedded pore water phase transition constraint is updated. Based on the updated degradation coefficient and Bayesian inference algorithm, the maximum value of the attenuation ratio of the surrounding rock stress relative to the initial stress of the surrounding rock when no support operation is performed is calculated. When the maximum value of the attenuation ratio is greater than the preset degradation upper limit, the degradation upper limit is adjusted upward and the strength degradation constraint is updated.

[0120] The updated embedded pore water phase change constraint and the updated strength degradation constraint are combined to obtain the optimized constraint parameters.

[0121] Based on the zoning support decision-making results, support operations were implemented in each high-risk connecting area. Taking a tunnel in a cold region as an example, antifreeze shotcrete with a thickness of 12 cm was applied to the arch crown area, coupled with 2.5-meter-long high-strength anchor bolts spaced 0.8 meters x 0.8 meters, and a waterproof membrane and steel mesh were added. Antifreeze shotcrete with a thickness of 10 cm was applied to the left side wall area, coupled with 2.2-meter-long anchor bolts spaced 0.9 meters x 0.9 meters, and a waterproof membrane was added. Shotcrete with a thickness of 8 cm was applied to the right side wall area, coupled with 2-meter-long anchor bolts spaced 1 meter x 1 meter. After the support was completed, monitoring points were set up in each high-risk connecting area to collect data on surrounding rock displacement, temperature, and stress, forming a surrounding rock response dataset. Fifteen monitoring points were set up in the arch crown area, ten in the left side wall area, and eight in the right side wall area. Data was collected hourly for 72 consecutive hours.

[0122] Based on the surrounding rock response data, the displacement and stress differences before and after support were calculated to obtain the displacement suppression and stress release amounts. The displacement suppression amount represents the degree of displacement reduction after support, while the stress release amount represents the magnitude of stress reduction after support. In the crown region, the average displacement before support was 32.5 mm, and the average displacement after support decreased to 8.3 mm, with a displacement suppression amount of 24.2 mm; the average stress before support was 12.8 MPa, and the average stress after support decreased to 5.7 MPa, with a stress release amount of 7.1 MPa. In the left wall region, the displacement suppression amount was 18.5 mm, and the stress release amount was 5.3 MPa; in the right wall region, the displacement suppression amount was 12.3 mm, and the stress release amount was 3.8 MPa.

[0123] The suppression effect was obtained by weighted summation of displacement suppression and stress release. The weight of displacement suppression was set to 0.6, and the weight of stress release was set to 0.4. The comprehensive suppression effect was calculated through normalization. The suppression effect in the crown area was 0.85, in the left wall area it was 0.72, and in the right wall area it was 0.63. The higher the suppression effect, the more effective the support measures. The surrounding rock temperature was compared with the phase change threshold embedded in the pore water phase change constraint, and the percentage of locations where the surrounding rock temperature was lower than the phase change threshold was statistically analyzed. The initial phase change threshold was set to -3 degrees Celsius, representing the temperature at which pore water begins to freeze. Monitoring results showed that the percentage of locations with temperatures lower than the phase change threshold was 18% in the crown area, 25% in the left wall area, and 31% in the right wall area. The higher the percentage of locations, the greater the risk of frost damage.

[0124] The phase transition threshold and strength degradation constraint were adjusted based on the relationship between the suppression effect and the location proportion. Analysis of the correlation between the suppression effect and the location proportion revealed that the suppression effect significantly decreased when the location proportion exceeded 20%. Based on this pattern, the phase transition threshold was adjusted from -3°C to -2.5°C to improve the sensitivity of the frost damage risk warning. Simultaneously, the degradation coefficient in the strength degradation constraint was adjusted. The initial value of the degradation coefficient was set to 0.25, representing the proportion of surrounding rock strength reduction caused by freeze-thaw action. Based on feedback on the suppression effect, the degradation coefficient was updated to 0.31 to more accurately reflect the degradation characteristics of surrounding rock strength under freeze-thaw conditions.

[0125] Based on updating the phase transition threshold and using a Bayesian inference algorithm, the temperature change rate at the phase transition location in the surrounding rock temperature is fitted with the latent heat parameter embedded in the pore water phase transition constraint. The phase transition location refers to the region where the temperature is near the phase transition threshold, and the temperature change rate represents the temperature change per unit time. The average temperature change rate at the phase transition location in the arch region is 0.12 degrees Celsius per hour, in the left wall region it is 0.15 degrees Celsius per hour, and in the right wall region it is 0.18 degrees Celsius per hour. A prior probability distribution is constructed using a Bayesian inference algorithm, with the initial latent heat parameter set at 334 kJ / kg. The prior distribution of the Bayesian model is set to a normal distribution with a mean of 334 and a variance of 25. The posterior distribution is calculated using the observed data to obtain the updated latent heat parameter. Parameter optimization is performed using the least squares method to minimize the sum of squares error between the predicted and actual observed temperature change rates. After 50 iterations of optimization, the updated latent heat parameter is 328 kJ / kg, slightly lower than the initial value, which better reflects the phase transition characteristics in the actual environment of tunnels in cold regions.

[0126] Based on the updated degradation coefficient and Bayesian inference algorithm, the attenuation ratio of the surrounding rock stress relative to the initial stress of the surrounding rock before support was calculated. The initial stress was obtained through inversion calculation. The initial stress in the crown region was 15.3 MPa, and the stress after support was 5.7 MPa, with an attenuation ratio of 62.7%. The initial stress in the left wall region was 12.6 MPa, and the stress after support was 5.2 MPa, with an attenuation ratio of 58.7%. The initial stress in the right wall region was 9.8 MPa, and the stress after support was 4.5 MPa, with an attenuation ratio of 54.1%. The maximum attenuation ratio for each region was 62.7%, which is greater than the preset degradation upper limit of 60%, indicating that the actual degree of degradation exceeded expectations. Therefore, the degradation upper limit was adjusted from 60% to 65%, and the strength degradation constraint was updated so that the model could more accurately reflect the severe degradation of the surrounding rock strength under extreme freeze-thaw conditions.

[0127] The updated phase transition constraints for embedded pore water and the updated strength degradation constraints are combined to obtain optimized constraint parameters. The updated phase transition constraints for embedded pore water include a phase transition threshold of -2.5 degrees Celsius and a latent heat parameter of 328 kJ / kg; the updated strength degradation constraints include a degradation coefficient of 0.31 and a degradation limit of 65%.

[0128] In this embodiment, by comparing and analyzing the changes in displacement and stress of the surrounding rock before and after support, and constructing a comprehensive suppression effect index that reflects the effects of displacement suppression and stress release, the problem of relying solely on experience or a single indicator to evaluate the support effect is avoided. This improves the objectivity and stability of the parameter adjustment basis. By comparing the surrounding rock temperature with the phase change threshold and statistically analyzing the spatial proportion of areas below the threshold, the freeze-thaw influence range can be accurately quantified. This allows for the identification of the actual changes in the distribution of freeze-thaw conditions caused by the support action, improving the environmental adaptability of phase change determination. By dynamically adjusting the phase change threshold and strength degradation constraint based on the relationship between the suppression effect and the proportion of low-temperature locations, the synchronous update of freeze-thaw phase change conditions and material performance degradation rules is achieved. By introducing Bayesian inference and least squares fitting, the temperature change rate and latent heat parameters at the phase change location are updated through inversion. This allows the energy parameters of the pore water phase change process to be corrected by measured data, significantly improving the physical consistency and prediction accuracy of the phase change constraint parameters.

[0129] Figure 2 This is a flowchart illustrating the optimization of freeze-thaw constraint parameters driven by surrounding rock response data in an embodiment of the present invention, based on a dynamic feedback-based intelligent decision-making method for tunnel support in cold regions.

[0130] In one alternative implementation,

[0131] The process of collecting real-time monitoring data, solving for instability characteristic parameters based on the optimization constraint parameters, and determining the corresponding optimal support decision scheme includes:

[0132] Real-time monitoring data is obtained by extracting the surrounding rock temperature, displacement, and stress from sensors in the tunnel surrounding rock. An updated phase transition threshold and an updated degradation coefficient are extracted from the optimized constraint parameters. The displacement change rate is calculated based on the real-time monitoring data, and the degradation rate is obtained by correlating the displacement change rate with the updated degradation coefficient. The unstable region is identified in the real-time monitoring data where the surrounding rock temperature is lower than the updated phase transition threshold and the degradation rate is greater than a preset rate threshold.

[0133] The stress concentration index is obtained by calculating the spatial and temporal gradients of the surrounding rock stress corresponding to the unstable region, and the instability characteristic parameters are obtained by combining the deterioration rate and the unstable region.

[0134] Based on the instability characteristic parameters, the support target is determined. Based on the support target, the support type is selected from the preset support type library and the support time is determined. The updated embedded pore water phase change constraint and the updated strength degradation constraint are used as constraints. The optimal support decision scheme is obtained by minimizing the instability characteristic parameters.

[0135] Real-time monitoring datasets are obtained by extracting data on surrounding rock temperature, displacement, and stress collected in real time by sensors in the tunnel surrounding rock. For example, in a tunnel project in a cold region, a distributed fiber optic sensing monitoring system was installed. A monitoring section was set up every 10 meters along the tunnel's longitudinal direction, with 12 monitoring points deployed at the arch crown, arch shoulder, sidewalls, and bottom. The sensors collected data every 30 minutes, achieving a temperature accuracy of 0.1 degrees Celsius, a displacement accuracy of 0.01 millimeters, and a stress accuracy of 0.05 MPa. An updated phase transition threshold and an updated degradation coefficient were extracted from the optimized constraint parameters. The updated phase transition threshold was set to -2.5 degrees Celsius, and the updated degradation coefficient was set to 0.31. The displacement change rate was calculated based on the real-time monitoring data, representing the change in displacement per unit time. At the tunnel mileage marker K25+600, the displacement in the arch area increased from 5.2 mm to 8.7 mm within 24 hours, with a displacement change rate of 0.146 mm per hour; the displacement in the left wall area increased from 4.5 mm to 6.8 mm, with a displacement change rate of 0.096 mm per hour; and the displacement in the right wall area increased from 3.9 mm to 5.6 mm, with a displacement change rate of 0.071 mm per hour.

[0136] The degradation rate is calculated by correlating the displacement change rate with the updated degradation coefficient. The degradation rate calculation considers the product of the displacement change rate and the updated degradation coefficient, and then corrects for it using a temperature influence factor. The temperature influence factor is determined based on the difference between the current temperature and the phase transition threshold; the closer the temperature is to the phase transition threshold, the larger the influence factor. For example, the temperature in the crown area is -2.2 degrees Celsius, the temperature influence factor is 0.92, and the calculated degradation rate is 0.042 per hour; the temperature in the left wall area is -1.8 degrees Celsius, the temperature influence factor is 0.75, and the degradation rate is 0.022 per hour; the temperature in the right wall area is -1.5 degrees Celsius, the temperature influence factor is 0.63, and the degradation rate is 0.014 per hour. Areas in the real-time monitoring data where the surrounding rock temperature is below the updated phase transition threshold and the degradation rate is greater than the preset rate threshold are identified as unstable areas. The preset speed threshold is set to 0.03 per hour. Based on this, the arch area is determined to have reached the instability condition and is therefore an unstable area; the left and right side wall areas have not yet reached the instability condition.

[0137] The spatial and temporal gradients of the surrounding rock stress corresponding to the instability zone were calculated to obtain the stress concentration index. The spatial gradient represents the rate of change of stress in space, calculated by dividing the stress difference between adjacent monitoring points by the distance between the points. The spatial distance between the instability zone at the crown and adjacent monitoring points is 0.8 meters, the stress difference is 1.6 MPa, and the spatial gradient is 2.0 MPa per meter. The temporal gradient represents the rate of change of stress over time, calculated by dividing the stress difference between adjacent time points by the time interval. The stress in the instability zone at the crown increased from 7.2 MPa to 8.8 MPa within 24 hours, with a temporal gradient of 0.067 MPa per hour. The stress concentration index was calculated by combining the spatial and temporal gradients. The stress concentration index is obtained by multiplying the spatial gradient by the temporal gradient and then normalizing. The calculated stress concentration index for the instability zone at the crown is 0.74. A higher index value indicates a more severe degree of stress concentration and a greater risk of instability.

[0138] Instability characteristic parameters are calculated by combining information such as the deterioration rate and the spatial location and extent of the instability zone. These parameters comprehensively consider factors such as the deterioration rate, stress concentration index, area, and depth of the instability zone. In the instability zone at the crown, the deterioration rate is 0.042 per hour, the stress concentration index is 0.74, the area of ​​the instability zone is approximately 5.3 square meters, and the depth is approximately 1.2 meters. The instability characteristic parameters are calculated using a weighted summation method, with the following weights for each factor: deterioration rate 0.35, stress concentration index 0.3, instability zone area 0.2, and depth 0.15. The calculated instability characteristic parameter for the instability zone at the crown is 0.68. A higher instability characteristic parameter indicates a greater risk of instability and a greater need for timely support measures.

[0139] Support objectives are determined based on instability characteristic parameters, including displacement control, temperature regulation, and stress release objectives. Graded objectives are established according to the numerical range of the instability characteristic parameters. When the instability characteristic parameters are between 0.6 and 0.8, the displacement control objective is to limit the displacement growth rate to no more than 0.05 mm / h; the temperature regulation objective is to maintain the surrounding rock temperature at no less than -1.5 degrees Celsius; and the stress release objective is to reduce the stress concentration index to below 0.4. Based on these objectives, suitable support types are selected from a pre-set support type library. This library contains various support combination schemes, such as ordinary shotcrete, antifreeze shotcrete, system anchors, prestressed anchors, steel arch frames, and grouting reinforcement, with each type containing multiple sets of parameter configurations. For the unstable area at the arch crown, the antifreeze shotcrete plus prestressed anchor support type is selected, and the support time is determined to be immediate and completed no later than 24 hours.

[0140] The optimal support decision scheme is obtained by minimizing the instability characteristic parameters using updated embedded pore water phase transition constraints and updated strength degradation constraints as constraints. The embedded pore water phase transition constraints include a phase transition threshold of -2.5 degrees Celsius and a latent heat parameter of 328 kJ / kg. The strength degradation constraints include a degradation coefficient of 0.31 and a degradation upper limit of 65%. A particle swarm optimization algorithm is used for the solution, with a population size of 120, a maximum number of iterations of 300, and a convergence condition of the optimal solution's change rate being less than 0.001 over 30 consecutive iterations. Optimization variables include shotcrete thickness, strength grade, mix proportion, and parameters such as the length, spacing, and prestress value of prestressed anchor cables. After 185 iterations, the algorithm converged and obtained the optimal support decision scheme: the thickness of the antifreeze shotcrete is 15 cm, the amount of antifreeze agent in the mix is ​​3.5%, the fiber content is 0.8%; the length of the prestressed anchor cable is 4.5 m, the spacing is 1.2 m × 1.2 m, and the prestress value is 180 kN.

[0141] Theoretical analysis of the optimal support decision scheme shows that after implementation, the expected displacement growth rate can be reduced to 0.028 mm / h, meeting the displacement control target; the surrounding rock temperature can be stabilized above -1.2 degrees Celsius, meeting the temperature control target; and the stress concentration index can be reduced to 0.35, meeting the stress release target. Simultaneously, the comprehensive instability characteristic parameter can be reduced from 0.68 to 0.32, significantly reducing the risk of instability. In addition to emergency support for the instability area of ​​the arch crown, the scheme also recommends preventative support for the left side wall area, using 10 cm thick antifreeze shotcrete and 3.5 m long ordinary anchor bolts spaced 1.5 m × 1.5 m apart; routine monitoring will be conducted on the right side wall area, with no additional support required at this time.

[0142] In this embodiment, the degradation rate is calculated based on the real-time surrounding rock displacement change rate and combined with the updated degradation coefficient, enabling the accelerated characteristics of surrounding rock deformation to be captured in a timely manner. This avoids the problem of delayed early warning caused by focusing only on the absolute displacement and ignoring the evolution trend, and improves the sensitivity of instability identification to early anomalies. By jointly determining whether the surrounding rock temperature is below the updated phase change threshold and whether the degradation rate exceeds the preset rate threshold, the collaborative identification of freeze-thaw environment triggering conditions and mechanical degradation process is realized, which can effectively reduce the probability of false alarms and false alarms, and make the identification of unstable areas more accurate and reliable. By comprehensively analyzing the spatial and temporal gradients of the stress in the unstable area and constructing a stress concentration index, the instability characteristics not only reflect the deformation acceleration phenomenon, but also simultaneously depict the internal stress redistribution and concentration, enhancing the physical interpretability of the risk identification results. By integrating the degradation rate, stress concentration index and unstable area to form an instability characteristic parameter, a quantitative comprehensive characterization of the surrounding rock instability state is realized, which can provide a more refined and continuous evaluation basis for support decision-making.

[0143] A second aspect of this invention provides an intelligent decision-making system for tunnel support in cold regions based on dynamic feedback, comprising:

[0144] The decoupled prediction unit is used to acquire historical monitoring data of the surrounding rock of tunnels in cold regions and perform time alignment and spatial mapping to obtain a monitoring dataset. It performs feature decoupling on the historical monitoring data to obtain freeze-thaw damage features and deformation field components. Based on the freeze-thaw damage features and pre-set embedded pore water phase transformation constraints and strength degradation constraints, it performs coupled calculations to predict the evolution trajectory of the surrounding rock damage field and obtain the evolution prediction result.

[0145] The graph construction unit is used to perform graph topology mining on the evolution prediction results, construct an evolution graph with monitoring points as nodes and damage propagation as edges, identify high-risk connected regions and critical propagation paths, calculate the stress concentration degree of the connected regions in combination with the deformation field components, and construct a risk assessment graph.

[0146] The support optimization unit is used to obtain the support timing sequence and support type scheme corresponding to each connected area based on the pre-acquired geological weak surface information and existing support status, combined with the risk assessment map, through a multi-objective optimization algorithm, and thus obtain the zonal support decision results;

[0147] The feedback update unit is used to implement support operations on the surrounding rock based on the partitioned support decision results and collect the surrounding rock response data after support. Based on the surrounding rock response data, it determines the suppression effect of the support operation and performs Bayesian update on the embedded pore water phase change constraint and strength degradation constraint based on the suppression effect to obtain optimized constraint parameters. It also collects real-time monitoring data and solves the instability characteristic parameters based on the optimized constraint parameters to determine the corresponding optimal support decision scheme.

[0148] A third aspect of the present invention provides an electronic device, comprising:

[0149] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0150] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0151] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0152] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A smart decision-making method for tunnel support in cold regions based on dynamic feedback, characterized in that, include: Historical monitoring data of the surrounding rock of tunnels in cold regions is acquired and time-aligned and spatially mapped to obtain a monitoring dataset. Feature decoupling is performed on the historical monitoring data to obtain freeze-thaw damage features and deformation field components. Based on the freeze-thaw damage features and pre-set embedded pore water phase transformation constraints and strength degradation constraints, coupled calculations are performed to predict the evolution trajectory of the surrounding rock damage field and obtain evolution prediction results. Graph topology mining is performed on the evolution prediction results to construct an evolution graph with monitoring points as nodes and damage propagation as edges, and high-risk connected regions and critical propagation paths are identified. The stress concentration degree of the connected regions is calculated in combination with the deformation field components, and a risk assessment map is constructed. Based on pre-acquired geological weak surface information and existing support status, and combined with the risk assessment map, a multi-objective optimization algorithm is used to obtain the support timing sequence and support type scheme corresponding to each connected area, thus obtaining the zonal support decision results; Based on the zonal support decision results, support operations are performed on the surrounding rock and the response data of the surrounding rock after support is collected. Based on the surrounding rock response data, the suppression effect of the support operation is determined, and based on the suppression effect, Bayesian updates are performed on the embedded pore water phase change constraint and strength degradation constraint to obtain optimized constraint parameters. Real-time monitoring data is collected, and the instability characteristic parameters are solved based on the optimized constraint parameters to determine the corresponding optimal support decision scheme.

2. The method according to claim 1, characterized in that, Historical monitoring data of the surrounding rock of tunnels in cold regions is acquired and subjected to time alignment and spatial mapping to obtain a monitoring dataset. Feature decoupling is then performed on the historical monitoring data to obtain freeze-thaw damage features and deformation field components, including: Raw data collected by various sensors located at different cross sections and burial depths over a historical period are extracted and combined with the timestamps and spatial coordinates of the data transmission to construct historical monitoring data. Cross-correlation analysis is performed on the timestamps of each sensor in the historical monitoring data to calculate the response delay difference of different sensors. Based on the response delay difference, the time offset of each sensor is determined and offset compensation is performed to obtain corrected data. Sensor channels with inconsistent sampling frequencies in the corrected data are resampled to obtain unified sampling data. The unified sampling data is combined with the spatial coordinates of the corresponding sensors to divide the data into three-dimensional grid nodes. The Euclidean distance from each sensor to each grid node is calculated, and the mapping weight is determined based on the Euclidean distance. The temperature field value, displacement field value, and strain field value of each grid node are solved based on the mapping weight and the unified sampling data and combined to obtain the monitoring dataset. The temperature field values ​​in the monitoring dataset are decomposed into time-frequency components to extract periodic fluctuation components and calculate the time derivative and spatial gradient of the temperature field values. Based on the time derivative and spatial gradient, the freeze-thaw interface is identified. The freeze-thaw damage characteristics are obtained by statistical analysis of the temperature change amplitude, time derivative peak value, and spatial gradient extreme value at the freeze-thaw interface. The deformation field components are obtained by spatial differentiation of the displacement field values ​​and strain field values ​​and extraction of the time evolution trend.

3. The method according to claim 1, characterized in that, Based on the freeze-thaw damage characteristics and pre-set embedded pore water phase transformation constraints and strength degradation constraints, coupled calculations are performed to predict the evolution trajectory of the surrounding rock damage field. The evolution prediction results include: The freeze-thaw interface in the freeze-thaw damage features is used as the phase change trigger boundary. The temperature change amplitude and time derivative peak value in the freeze-thaw damage features are extracted as phase change parameters. Based on the phase change trigger boundary and the phase change parameters, the pore water volume expansion rate is calculated in the pre-set embedded pore water phase change constraint and converted into an expansion stress field. The stress concentration location is determined based on the spatial gradient extrema in the expansion stress field and freeze-thaw damage characteristics. The strain field value history sequence of the corresponding grid node in the monitoring dataset is extracted at the stress concentration location and accumulated over time to obtain the cumulative strain. Based on the cumulative strain and the pre-set strength degradation constraint, the elastic modulus attenuation coefficient and compressive strength reduction coefficient are calculated and coupled with the expansion stress field to obtain the equivalent stress field. Using the equivalent stress field as the initial stress, and combining the freeze-thaw damage characteristics to set the evolution time step, stress redistribution is performed on the equivalent stress field within each time step, and damage points are determined. The damage points are connected and tracked in the spatial domain, and the occurrence order and position changes of each damage point are recorded and combined to obtain the evolution trajectory. The spatial distribution characteristics and temporal evolution characteristics of the damage field corresponding to the evolution trajectory are extracted and summarized to obtain the evolution prediction results.

4. The method according to claim 1, characterized in that, Graph topology mining is performed on the evolution prediction results to construct an evolution graph with monitoring points as nodes and damage propagation as edges. High-risk connected regions and critical propagation paths are identified. The stress concentration degree of the connected regions is calculated by combining the deformation field components, and a risk assessment map is constructed, including: Based on the spatial distribution characteristics of the damage field in the evolution prediction results, the spatial coordinates of each damage point are extracted and spatially matched with the grid nodes in the monitoring dataset. According to the matching results, the damage points are mapped to the nearest grid node and marked as graph nodes. The position change vector of each damage point between adjacent time steps is obtained and the corresponding orientation angle and length are calculated. Damage points with the same orientation angle and continuous length are identified and combined to obtain a propagation chain. Graph nodes corresponding to adjacent damage points in the propagation chain are connected and directed edges are created. The weight of the directed edges is set based on the time interval between adjacent damage points to obtain an evolution graph. The directed edges in the evolution graph are traversed and fast edges are determined in combination with a pre-set weight threshold. The connectivity of the fast edges is detected and the set of nodes reachable through the fast edges is marked as a high-risk connected region. Traverse all paths from the initial damage point to the boundary node of the high-risk connected region in the evolution graph and determine the velocity index based on the sum of the weights of all directed edges on the path. The path with the largest velocity index is taken as the critical propagation path. The high-risk connected regions are matched with the deformation field components in the spatial dimension, and the strain rate and strain gradient at the corresponding positions are extracted. The spatial second derivatives of the strain rate and strain gradient are calculated and the stress concentration is determined. The risk assessment map is obtained by combining the high-risk connected regions and the evolution map.

5. The method according to claim 1, characterized in that, Based on pre-acquired geological weak surface information and existing support status, and combined with the risk assessment map, a multi-objective optimization algorithm is used to obtain the support timing sequence and support type scheme corresponding to each connected region, resulting in the following zonal support decision results: The fault strike and joint density in the pre-acquired geological weak surface information are superimposed and matched with the high-risk connected regions in the risk assessment map to identify the coupling zone and extract the weakening coefficient. Extract the spatial location and strength grade of the constructed support from the pre-acquired existing support status, calculate the shortest distance between the coupling zone and the spatial location of the constructed support to obtain the coverage, and calculate the gap amount by the difference between the strength grade and the pre-acquired stress concentration degree. Using the weakening coefficient, coverage, and gap amount as constraints, the spatial orientation of the pre-obtained critical propagation path as boundary conditions, and the stress concentration degree of each high-risk connected area as risk weight, an optimization function is constructed with the objectives of minimizing the overall risk exposure time and maximizing the support resource utilization rate. The optimal intervention time and optimal strength configuration are obtained by iteratively solving the optimization function through a multi-objective optimization algorithm. The optimal intervention time is sorted according to the order of each high-risk connected area on the critical propagation path to obtain the support timing sequence. The optimal strength configuration is matched with the preset support type library. Support structure and support parameters are assigned to each high-risk connected area to obtain the support type scheme. The support timing sequence and the support type scheme are combined according to the spatial location of the high-risk connected areas to obtain the zonal support decision result.

6. The method according to claim 1, characterized in that, Based on the zonal support decision results, support operations are implemented on the surrounding rock, and the response data of the surrounding rock after support is collected. Based on the surrounding rock response data, the suppression effect of the support operation is determined, and based on the suppression effect, Bayesian updates are performed on the embedded pore water phase change constraint and strength degradation constraint to obtain optimized constraint parameters, including: Based on the zonal support decision results, support operations are carried out in each high-risk connected area, and surrounding rock displacement, temperature, and stress are collected after support is completed to obtain surrounding rock response data. Based on the surrounding rock response data, the displacement difference and stress difference before and after support are calculated to obtain the displacement suppression amount and stress release amount. The displacement suppression amount and stress release amount are weighted and summed to obtain the suppression effect. The surrounding rock temperature is compared with the phase change threshold in the embedded pore water phase change constraint, and the proportion of positions where the surrounding rock temperature is lower than the phase change threshold is statistically analyzed. The phase transition threshold and the intensity degradation constraint are adjusted according to the relationship between the suppression effect and the position ratio to obtain the updated phase transition threshold and the updated degradation coefficient; Based on the updated phase transition threshold and Bayesian inference algorithm, the temperature change rate at the phase transition position in the surrounding rock temperature is fitted with the latent heat parameter in the embedded pore water phase transition constraint. The updated latent heat parameter is obtained by combining the least squares method and the embedded pore water phase transition constraint is updated. Based on the updated degradation coefficient and Bayesian inference algorithm, the maximum value of the attenuation ratio of the surrounding rock stress relative to the initial stress of the surrounding rock when no support operation is performed is calculated. When the maximum value of the attenuation ratio is greater than the preset degradation upper limit, the degradation upper limit is adjusted upward and the strength degradation constraint is updated. The updated embedded pore water phase change constraint and the updated strength degradation constraint are combined to obtain the optimized constraint parameters.

7. The method according to claim 1, characterized in that, The process of collecting real-time monitoring data, solving for instability characteristic parameters based on the optimization constraint parameters, and determining the corresponding optimal support decision scheme includes: Real-time monitoring data is obtained by extracting the surrounding rock temperature, displacement, and stress from sensors in the tunnel surrounding rock. An updated phase transition threshold and an updated degradation coefficient are extracted from the optimized constraint parameters. The displacement change rate is calculated based on the real-time monitoring data, and the degradation rate is obtained by correlating the displacement change rate with the updated degradation coefficient. The unstable region is identified in the real-time monitoring data where the surrounding rock temperature is lower than the updated phase transition threshold and the degradation rate is greater than a preset rate threshold. The stress concentration index is obtained by calculating the spatial and temporal gradients of the surrounding rock stress corresponding to the unstable region, and the instability characteristic parameters are obtained by combining the deterioration rate and the unstable region. Based on the instability characteristic parameters, the support target is determined. Based on the support target, the support type is selected from the preset support type library and the support time is determined. The updated embedded pore water phase change constraint and the updated strength degradation constraint are used as constraints. The optimal support decision scheme is obtained by minimizing the instability characteristic parameters.

8. A dynamic feedback-based intelligent decision-making system for tunnel support in cold regions, used to implement the method described in any one of claims 1-7, characterized in that, include: The decoupled prediction unit is used to acquire historical monitoring data of the surrounding rock of tunnels in cold regions and perform time alignment and spatial mapping to obtain a monitoring dataset. It performs feature decoupling on the historical monitoring data to obtain freeze-thaw damage features and deformation field components. Based on the freeze-thaw damage features and pre-set embedded pore water phase transformation constraints and strength degradation constraints, it performs coupled calculations to predict the evolution trajectory of the surrounding rock damage field and obtain the evolution prediction result. The graph construction unit is used to perform graph topology mining on the evolution prediction results, construct an evolution graph with monitoring points as nodes and damage propagation as edges, identify high-risk connected regions and critical propagation paths, calculate the stress concentration degree of the connected regions in combination with the deformation field components, and construct a risk assessment graph. The support optimization unit is used to obtain the support timing sequence and support type scheme corresponding to each connected area based on the pre-acquired geological weak surface information and existing support status, combined with the risk assessment map, through a multi-objective optimization algorithm, and thus obtain the zonal support decision results; The feedback update unit is used to implement support operations on the surrounding rock based on the partitioned support decision results and collect the surrounding rock response data after support. Based on the surrounding rock response data, it determines the suppression effect of the support operation and performs Bayesian update on the embedded pore water phase change constraint and strength degradation constraint based on the suppression effect to obtain optimized constraint parameters. It also collects real-time monitoring data and solves the instability characteristic parameters based on the optimized constraint parameters to determine the corresponding optimal support decision scheme.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images