Cold region jointed rock mass fracture evolution analysis system based on visual inspection

Abstract

The invention relates to the technical field of geotechnical engineering safety monitoring, and discloses a cold region jointed rock mass fracture evolution analysis system based on visual inspection. The system comprises a multi-source space-time registration module, a freeze-thaw damage feature extraction module, a fracture network growth prediction module, a coupling effect evolution correction module and a stability control list generation module which are connected in sequence. According to the method, space-time registration is realized on a time sequence by fusing a rock mass surface image and internal point cloud data, and a fracture characteristic spectrum is extracted. And predicting a potential path by using a growth kinetic model, performing physical mechanism correction on the predicted path by applying freeze thawing-loading coupling constraint, and finally outputting a key dangerous path and a stability control list. According to the system, dynamic and accurate analysis and risk early warning of the cold region fracture evolution process are realized.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering safety monitoring technology, specifically a visual detection-based system for analyzing the evolution of jointed rock mass fractures in cold regions. Background Technology

[0002] In the safety assessment of engineering geology and geotechnical engineering in cold regions, accurately analyzing the fracture evolution process of jointed rock masses under freeze-thaw cycles is crucial to ensuring engineering safety. Existing technologies typically separate surface observation from internal detection. Surface monitoring relies on periodic manual surveys, fixed-point photography, or the deployment of traditional sensors; these methods struggle to comprehensively capture the spatial expansion morphology of fractures. Internal detection employs single-point-time 3D laser scanning, drilling, or acoustic testing, failing to achieve continuous observation with high spatiotemporal resolution. The two types of data are often independent in terms of acquisition time and coordinate system, creating information silos. Furthermore, existing prediction models are mostly based on assumptions of ambient temperature conditions or uniformly distributed loads, failing to dynamically couple the unique cyclic water heave force and thaw contraction effect of cold regions with the actual tectonic stress and engineering loads borne by the rock mass, leading to deviations between predicted paths and actual physical processes.

[0003] These shortcomings make it difficult for existing technologies to fully reproduce the three-dimensional dynamic process of crack initiation from the surface to internal penetration, and also make it impossible to reliably identify the key dangerous paths most prone to instability under the coupling of freeze-thaw cycles and stress fields. This invention aims to solve these problems by achieving dynamic fusion of visual information about the rock mass surface and internal structural data over the time dimension of freeze-thaw cycles, and establishing a physical mechanism model of the coupling between freeze-thaw action and in-situ stress fields to constrain and correct crack evolution paths, thereby enabling accurate analysis, early warning, and targeted control of the stability of fractured rock masses in cold regions. Summary of the Invention

[0004] The purpose of this invention is to provide a visual detection-based system for analyzing the evolution of fractures in cold-region jointed rock masses, in order to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, this invention provides a visual detection-based system for analyzing the fracture evolution of jointed rock masses in cold regions, the system comprising: The multi-source spatiotemporal registration module is used to acquire a set of multi-angle images and a set of internal structure point clouds generated by freeze-thaw cycles on the surface of a target jointed rock mass in a cold environment within a set time period, and to perform registration and fusion processing on the multi-angle image set and the internal structure point cloud set in a time series to generate a spatiotemporally registered multi-source data array. The freeze-thaw damage feature extraction module is used to extract crack growth features from the spatiotemporally registered multi-source data array to obtain a freeze-thaw damage feature spectrum, which includes crack surface roughness distribution, crack branch node coordinate set, and crack boundary displacement vector field. The fracture network growth prediction module is used to perform fracture path prediction processing based on the freeze-thaw damage characteristic spectrum and the fracture network growth dynamics model to generate a set of potential penetration paths and rock mass block isolation status identifiers. The coupling effect evolution correction module is used to apply the freeze-thaw-loading coupling effect constraint unique to cold regions to the set of potential through paths, perform evolution path correction calculations, and output the corrected set of critical dangerous paths and their evolution stage determination results. The stability control list generation module is used to perform cross-validation analysis between the set of critical hazardous paths and their evolution stage determination results and the rock mass block isolation status identifiers to generate a stability control list for cold-region fractured rock masses that includes priority reinforcement levels and dynamic monitoring frequencies.

[0006] Preferably, the step of extracting crack growth features from the spatiotemporally registered multi-source data array to obtain a freeze-thaw damage feature spectrum includes the following process: From the spatiotemporally registered multi-source data array, image subsets and point cloud subsets corresponding to different freeze-thaw cycles are separated; Joint feature extraction is performed on each of the image subsets and point cloud subsets respectively: Perform crack contour boundary recognition on the subset of images to obtain the set of pixel coordinates of the crack contour boundary; Perform crack space surface fitting on the point cloud subset to obtain a three-dimensional crack space surface model; The set of pixel coordinates of the crack contour boundary is mapped and matched with the three-dimensional spatial surface model of the crack. The surface roughness distribution of the crack is generated by calculating the projection deviation of the pixel coordinates of the crack contour boundary on the three-dimensional spatial surface model of the crack. On the three-dimensional spatial surface model of the fracture, the abrupt change points of the surface curvature are detected, and the spatial coordinates of the curvature abrupt change points are recorded as a set of fracture branch node coordinates; Compare the deformation differences of the three-dimensional spatial surface model of the crack during adjacent freeze-thaw cycles, calculate the displacement of each point on the three-dimensional spatial surface model of the crack in the normal direction, and form a crack boundary displacement vector field. By integrating the surface roughness distribution of the crack, the coordinate set of the crack branch nodes, and the displacement vector field of the crack boundary, a complete freeze-thaw damage characteristic spectrum is generated.

[0007] Preferably, the step of using a fracture network growth dynamics model to predict fracture paths based on the freeze-thaw damage characteristic spectrum, and generating a set of potential penetration paths and rock mass block isolation status identifiers, includes the following processing: The main fracture intersection point is identified as the growth starting point from the set of fracture branch node coordinates; Based on the directional trend provided by the fracture boundary displacement vector field, a fracture growth direction guiding field is set in the fracture network growth dynamics model; Taking the growth starting point as the origin, and under the constraint of the crack growth direction guiding field, multi-directional random growth simulation is performed according to the preset growth probability rules. In the multi-directional random growth simulation process, the surface roughness distribution of the crack is used as a resistance factor to calculate the energy consumption of each random growth path; When the simulated growth path of the fracture encounters other nodes in the set of fracture branch node coordinates, or when the energy consumption of the growth path exceeds a preset threshold, the growth simulation of the fracture path is stopped. Record all random growth paths that reach the stopping condition and integrate them into a set of potential through paths; Based on the distribution of the potential through paths in three-dimensional space, it is determined whether they divide the rock mass into independent blocks, and the divided independent blocks are spatially marked to generate rock mass block isolation status identifiers.

[0008] Preferably, the step of applying freeze-thaw-loading coupling constraints specific to cold regions to the set of potential through paths, performing evolution path correction calculations, and outputting the corrected set of critical hazardous paths and their evolution stage determination results is achieved through the following method: Historical temperature cycle data and geostress field data of the target area are obtained from external data sources. The historical temperature cycle data and geostress field data are fused together to establish a stress intensity factor field that couples freeze-thaw cycles and geostress. From the set of potential connecting paths, extract the spatial orientation and geometric shape of each path one by one; Substitute the spatial orientation and geometry of each path into the stress intensity factor field to calculate the stress concentration degree of each extracted path under the current freeze-thaw-loading coupling effect. A cumulative frost heave damage model for cold-region rock masses is introduced, and the frost heave fatigue damage value along each path is calculated based on the historical temperature cycle data. The stress concentration level and the frost heave fatigue damage value are weighted and superimposed to obtain the coupling hazard index of each potential through path; Based on the coupling hazard index, paths exceeding a preset hazard threshold are selected to form a set of critical hazard paths; For each path in the set of critical hazardous paths, an evolution stage determination result is assigned based on the growth rate of its stress concentration and the accumulation rate of its frost heave fatigue damage value. The evolution stage determination result includes the germination stage, the stable expansion stage, and the accelerated expansion stage.

[0009] Preferably, the specific steps of acquiring historical temperature cycle data and geostress field data of the target area from an external data source, fusing the historical temperature cycle data and geostress field data to establish a stress intensity factor field coupled with freeze-thaw cycles and geostress are as follows: Temperature time series data of the target area within the historical monitoring period is obtained from an external temperature sensor network. The daily maximum and minimum temperatures are extracted from the temperature time series data. The daily temperature difference and the frequency of temperature fluctuation around zero degrees are calculated, and a temperature cycle feature matrix is ​​synthesized. Data on the triaxial geostress components of the target area are obtained from external geostress monitoring devices, and combined with the regional geological structure model to generate a spatially continuous geostress field distribution model. A stress transformation function is established, which takes the temperature cycle characteristic matrix and the geostress field distribution model as inputs to calculate the combined stress value generated at different locations of the rock mass due to the superposition of frost heave force and geostress. The synthesized stress values ​​are mapped onto a three-dimensional rock mass space model according to a preset grid to form a stress intensity factor field. The value of each grid point in the stress intensity factor field represents the stress intensity factor of the spatial position of the grid point in the three-dimensional rock mass space model under freeze-thaw-loading coupling.

[0010] Preferably, the process of cross-validating the set of critical hazard paths and their evolution stage determination results with the rock mass block isolation status identifiers to generate a stability control list for cold-region fractured rock masses, including priority reinforcement levels and dynamic monitoring frequencies, is completed through the following steps: Obtain the evolution stage determination result of each path in the set of critical dangerous paths, as well as the volume and spatial location information of each isolated block in the rock mass block isolation status identifier; Each path in the set of critical hazardous paths is spatially superimposed with the rock mass block isolation status markers to identify paths that pass through critical blocks or may cause instability of critical blocks, and these paths are defined as high-risk paths. An initial priority reinforcement level is assigned to each of the high-risk paths. The assignment rule is based on the determination result of its evolution stage, where the accelerated expansion stage corresponds to the highest level, the stable expansion stage corresponds to the second highest level, and the nascent stage corresponds to the basic level. Based on the volume of the critical blocks affected by the high-risk paths, the initial priority reinforcement level is adjusted. The larger the volume of the affected blocks, the greater the increase in the priority reinforcement level of the corresponding path. Based on the priority reinforcement level of the high-risk path, a dynamic monitoring frequency is set for it. The higher the priority reinforcement level, the higher the dynamic monitoring frequency is allocated. By compiling the route information of all high-risk routes, the adjusted priority reinforcement levels, and the corresponding dynamic monitoring frequencies, a final list of stability control measures for fractured rock masses in cold regions is formed.

[0011] Preferably, the initial priority reinforcement level is adjusted based on the volume of the critical blocks affected by the high-risk path, specifically implemented as follows: In the rock mass block isolation status identification, the rock mass blocks that intersect with or are cut by each high-risk path are located; Calculate the three-dimensional spatial volume of each located rock mass block; A volume-level adjustment mapping table is established, which defines the amount of priority reinforcement level adjustment corresponding to rock mass blocks in different volume ranges. The volume-level adjustment mapping table is consulted based on the volume of each rock mass block to obtain the corresponding upward adjustment amount; The adjusted amount is added to the initial priority reinforcement level of the high-risk path to obtain the adjusted priority reinforcement level, and it is ensured that the adjusted level does not exceed the highest level threshold set by the system.

[0012] Preferably, the establishment of the stress transformation function, which takes the temperature cycle characteristic matrix and the geostress field distribution model as inputs, calculates the combined stress values ​​generated at different locations in the rock mass due to the superposition of frost heave force and geostress, specifically including: Extract the daily temperature difference and the cumulative duration of sub-zero temperature for each grid point from the temperature cycle feature matrix; Extract the three principal stress components at each grid point from the aforementioned geostress field distribution model; A frost heave force calculation sub-function is constructed. The frost heave force calculation sub-function uses the daily temperature difference and the cumulative duration of sub-zero temperature as variables, and outputs the maximum frost heave force vector generated by the current processing grid point in the temperature cycle feature matrix within a freeze-thaw cycle. A sub-function for calculating geostress contribution is constructed. The sub-function takes three principal stress components as input and combines them with the anisotropic elastic parameters of the rock mass to calculate the stress tensor caused by the geostress field at the grid point. The maximum frost heave force vector is superimposed with the stress tensor caused by the geostress field to calculate the composite stress tensor at the grid point. The maximum principal stress value is extracted from the synthetic stress tensor and used as the synthetic stress value representing the stress intensity at the grid point.

[0013] Preferably, the step of substituting the spatial orientation and geometry of each path into the stress intensity factor field to calculate the stress concentration degree of each extracted path under the current freeze-thaw-loading coupling effect is specifically implemented as follows: Along the extension direction of each potential through path, a series of path feature points are set at preset intervals; For each path feature point, the stress intensity factor value corresponding to its spatial location is queried in the stress intensity factor field. If the path feature point is not located on the precise grid node of the stress intensity factor field, its stress intensity factor interpolation is calculated by the cubic spline interpolation algorithm. Extract the stress intensity factor values ​​at all path feature points to form a local stress distribution sequence for the potential through path; Calculate the statistical characteristics of the local stress distribution sequence, wherein the statistical characteristics include at least: the maximum value of the sequence, the average value of the sequence, and the stress change gradient along the path; Based on the spatial orientation of the potential through path, calculate the average angle between its main extension direction and the direction of the maximum principal stress in the stress intensity factor field. The maximum value of the sequence, the average value of the sequence, the stress change gradient, and the average included angle are input into a preset stress concentration coefficient calculation model. The stress concentration coefficient calculation model outputs a dimensionless coefficient characterizing the overall stress concentration degree of the potential through path through a weighted fusion algorithm, which is used as the stress concentration degree.

[0014] Preferably, the method further includes a process of continuously optimizing the fracture network growth dynamics model, the process including: After a set time interval, the latest image and point cloud data of the target jointed rock mass are reacquired, and the latest freeze-thaw damage feature spectrum is obtained through the crack growth feature extraction process. Extract the observed crack propagation path information from the latest freeze-thaw damage feature spectrum; The observed crack propagation path information is matched and compared with the set of potential through paths previously predicted by the crack network growth dynamics model. Calculate the difference between the actual observed length and the model predicted length for each matching path, as well as the deviation angle between the actual observed direction and the model predicted direction; Based on the average value of the difference and deviation angle, the model growth probability rule correction parameters are generated; The growth probability rule correction parameters of the fracture network growth dynamics model are iteratively updated using the model growth probability rule correction parameters to obtain an optimized fracture network growth dynamics model, which is then used for subsequent fracture path prediction processing.

[0015] Compared with the prior art, the beneficial effects of the present invention are: By precisely registering and fusing multi-angle image sequences of the rock mass surface acquired within a set time period with point cloud data of the internal structure, a spatiotemporally aligned multi-source data array is generated. This allows the visible fracture textures and morphological changes on the surface to be correlated with the spatial location and geometric topology of the internal fractures within a unified time axis and three-dimensional coordinate system. It enables the complete reconstruction of the three-dimensional spatiotemporal evolution trajectory of fractures, from their initial initiation on the rock mass surface to their gradual extension and branching interactions, achieving three-dimensional dynamic visualization and tracing of the entire fracture system development process, overcoming the blind spots in spatiotemporal coverage of a single data source.

[0016] Based on the potential fracture paths predicted by the preliminary dynamic model, the evolution path was corrected by specifically applying freeze-thaw-loading coupling constraints under cold-region environments. The interaction and superposition effects between cyclic frost heave forces and the in-situ stress field of the rock mass were quantified. This ensures that the corrected fracture path prediction is no longer simply the result of geometric extension or a single force field, but rather truly reflects how the dynamic superposition of frost heave stress concentration areas and the original stress field under cold-region environments alters the stress intensity factor at the fracture tip, thereby driving fractures to preferentially extend along specific directions. This makes the ultimately identified "critical hazardous paths" more consistent with the actual failure mechanism of cold-region rock masses, improving the engineering physical realism and reliability of the prediction results. Attached Figure Description

[0017] Figure 1 This is a time-series diagram of the visual detection-based jointed rock mass fracture evolution analysis system for cold regions described in this invention. Figure 2 This is a flowchart for extracting freeze-thaw damage features; Figure 3 A flowchart for establishing the stress intensity factor field; Figure 4 A scatter plot showing the relationship between the block volume and the increase in reinforcement level; Figure 5 This is a graph showing the relationship between the volume of the rock mass block and the reinforcement level and monitoring frequency. Detailed Implementation

[0018] 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.

[0019] Please see Figure 1This invention provides a visual inspection-based system for analyzing the evolution of fractures in cold-region jointed rock masses. The system includes: First, a multi-source spatiotemporal registration module operates within a set monitoring period. This module uses image acquisition and 3D scanning equipment deployed around the target jointed rock mass in the cold region to simultaneously acquire a set of multi-angle images of the rock mass surface caused by freeze-thaw cycles, as well as a set of point clouds of the internal structure of the rock mass. These data are uniformly timestamped. The module uses feature matching and coordinate transformation algorithms to spatially align the image information and point cloud information at the same time point, and arranges and correlates the registration data from different time points according to a time sequence, ultimately generating a multi-source data array containing spatiotemporal dimensional information and integrating surface texture and 3D geometric structure. Subsequently, a freeze-thaw damage feature extraction module processes this multi-source data array. Its core function is to identify and quantify the growth characteristics of fractures in different freeze-thaw cycles. This module extracts the surface roughness changes of fractures, the generation and coordinates of branch nodes, and the displacement field of fracture boundaries by comparing and analyzing data from different time points. These features are structured and integrated into a freeze-thaw damage feature spectrum.

[0020] The fracture network growth prediction module receives the feature spectrum as input. Its pre-built fracture network growth dynamics model, based on the fracture initiation point, extension direction trend, and growth resistance represented by surface roughness provided in the feature spectrum, simulates multiple paths for random fracture growth in the future. It then determines whether these potential paths will cut the rock mass into isolated blocks, thus outputting a set of potential through-paths and an identifier of the rock mass's block isolation status. The coupling effect evolution correction module further introduces external constraints unique to cold-region environments. This module obtains temperature history and geostress field data of the monitored area from an external system, constructs a stress intensity factor field reflecting the coupling effect of freeze-thaw cycles and geostress, and uses this field to assess the stress concentration and frost heave damage of each path in the potential through-path set. It then selects the most dangerous path with the highest probability of expansion under coupling effects, determines its expansion stage, and outputs a corrected set of critical dangerous paths and their evolution stage determination results. Finally, the stability control list generation module comprehensively analyzes the information on critical hazardous paths and rock mass segments. It assesses which hazardous paths will endanger the stability of key rock blocks and assigns appropriate priority reinforcement levels and dynamic monitoring frequencies to each high-risk path based on the path's degree of hazard, evolution stage, and the scale of the affected rock blocks. Ultimately, it generates a stability control list for cold-region fractured rock masses that can be directly used for engineering decisions, thus completing the closed loop from data perception to control decision-making.

[0021] Example 1: See Figure 2The freeze-thaw damage feature extraction module separates image subsets and point cloud subsets corresponding to different freeze-thaw cycles from the spatiotemporally registered multi-source data array according to time index, and performs joint feature extraction processing on the image subsets and point cloud subsets of each cycle. An edge detection and contour tracking-based crack contour boundary recognition algorithm is applied to the image subsets to obtain a set of crack contour boundary pixel coordinates consisting of a series of pixel coordinates. For the point cloud subsets of the same cycle, crack spatial surface fitting based on least squares fitting or triangulated surface reconstruction is performed to generate a spatial surface model describing the three-dimensional morphology of the crack. The set of crack contour boundary pixel coordinates identified from the image subsets is mapped onto the three-dimensional spatial surface model of the crack in the same cycle through the established camera imaging model and spatial transformation relationship. By calculating the statistical distribution of the projection deviation between the mapped points and the model surface, the crack surface roughness distribution is generated. On the crack three-dimensional spatial surface model, a surface curvature analysis algorithm based on principal curvature calculation is used to detect points where the surface curvature changes abruptly, and the three-dimensional spatial coordinates of these abrupt points are recorded as a set of crack branch node coordinates. By comparing the three-dimensional spatial surface models of the cracks generated in two adjacent freeze-thaw cycles, and employing nearest-point iteration or deformation field calculation algorithms, the displacement of corresponding points on the model of the later cycle relative to the model of the previous cycle is obtained. Specifically, the displacement components of each point in the surface normal direction are calculated to form a crack boundary displacement vector field. The crack surface roughness distribution, the set of crack branch node coordinates, and the crack boundary displacement vector field are integrated and encapsulated into a complete freeze-thaw damage feature spectrum containing spatial geometry, morphology, and kinematic information.

[0022] The fracture network growth prediction module operates based on the received freeze-thaw damage characteristic spectrum. It identifies the main fracture intersection points from the fracture branch node coordinate set based on the number and length of connecting fractures, and sets these points as the starting point for fracture network growth. Based on the fracture boundary movement trend revealed by the fracture boundary displacement vector field, a fracture growth direction guidance field is constructed in the fracture network growth dynamic model to indicate the future fracture propagation direction. Using the selected growth starting point as the origin, under the probabilistic constraints of the fracture growth direction guidance field, multi-directional random growth simulation is performed according to preset growth step size and random turning angle rules. During the simulated growth process, the fracture surface roughness distribution data is used as a quantification factor for local growth resistance, and the accumulated energy consumption of each random growth path is calculated in conjunction with the path length. When a simulated growth fracture path encounters other nodes recorded in the fracture branch node coordinate set in space, or when the accumulated energy consumption value of the path exceeds a preset fracture energy threshold, the current growth simulation of that fracture path is stopped. All random growth paths that stop growing due to encounters or energy exceeding the limit are recorded, and their spatial coordinate sequences are integrated into a set of potential through-paths. Based on the distribution of the potential through paths in the 3D rock mass model, spatial topology analysis is performed to determine whether these paths are interconnected to form closed loops, thereby dividing the rock mass into independent blocks. Each identified independent block is then uniquely spatially labeled to generate a rock mass block isolation status identifier.

[0023] In practical implementation, the freeze-thaw damage feature extraction module separates image subsets and point cloud subsets corresponding to different freeze-thaw cycles from the spatiotemporally registered multi-source data array. The module performs joint feature extraction processing on each image subset and point cloud subset. For the image subsets, it performs crack contour boundary recognition based on the Canny or Sobel operator to obtain a set of crack contour boundary pixel coordinates. For the point cloud subsets, it performs crack spatial surface fitting based on moving least squares or Poisson surface reconstruction to obtain a three-dimensional crack spatial surface model. The crack contour boundary pixel coordinate set is then mapped and matched with the three-dimensional crack spatial surface model. This mapping and matching transforms the two-dimensional image through the relationship between the camera imaging model and coordinate transformation. The primitive coordinates are mapped to a 3D surface, and the standard deviation sequence of the projection deviation is calculated to generate the surface roughness distribution of the crack. A surface curvature abrupt change point detection algorithm based on principal curvature analysis is used on the 3D spatial surface model of the crack to detect abrupt changes in surface curvature and record the spatial coordinates of these abrupt changes as a set of crack branch node coordinates. The deformation differences of the 3D spatial surface model of the crack are compared within adjacent freeze-thaw cycles. An iterative nearest-point algorithm is used to calculate the displacement of each point on the 3D spatial surface model of the crack in the normal direction to form a crack boundary displacement vector field. The freeze-thaw damage feature extraction module integrates the surface roughness distribution of the crack, the set of crack branch node coordinates, and the crack boundary displacement vector field to generate a complete freeze-thaw damage feature spectrum. In some embodiments, the set of pixel coordinates of the crack contour boundary is stored in chain code or ordered point sequence format, the 3D spatial surface model of the crack is represented in the form of a triangular mesh, and the surface roughness distribution of the crack is represented by the root mean square deviation value on each triangular facet.

[0024] In practical implementation, the fracture network growth prediction module operates based on the freeze-thaw damage characteristic spectrum. It identifies the main fracture intersection point as the growth starting point from the fracture branch node coordinate set based on the number and geometric length of connected fractures. According to the directional trend provided by the fracture boundary displacement vector field, a fracture growth direction guiding field is set in the fracture network growth dynamics model. This guiding field is a three-dimensional spatial vector field. With the growth starting point as the origin, multi-directional random growth simulation is performed under the constraint of the guiding field, following a preset growth probability rule. The preset growth probability rule defines that the probability of a fracture choosing a certain direction during the next expansion step is positively correlated with the alignment degree between that direction and the guiding field vector. During the multi-directional random growth simulation, the fracture surface roughness distribution is used as a resistance factor to calculate the path consumption of each random growth path. Path consumption is a dimensionless cumulative cost indicator. When the simulated fracture growth path encounters other nodes in the fracture branch node coordinate set at a spatial Euclidean distance, or when the path consumption of the growth path exceeds a preset threshold, the fracture path growth simulation stops. All random growth paths that meet the stopping condition are recorded, and their spatial node sequences are integrated into a potential through-path set. The fracture network growth prediction module determines whether the potential through-path set will divide the rock mass into independent blocks based on its distribution in three-dimensional space, and spatially marks the divided independent blocks to generate rock mass block isolation status identifiers. In some embodiments, the multi-directional random growth simulation uses a random walk algorithm, where the random turning angle of each step follows a normal distribution with the fracture growth direction guiding the field direction as the mean. The formula controlling the probability of the growth direction in the fracture network growth dynamics model is:

[0025] in: This indicates that the crack extends from the current position along the direction... Extended growth probability, This represents a normalized partition function that ensures the sum of probabilities in all directions is 1. This represents the weighting coefficient of the orientation alignment term. This represents the field vector guiding the crack growth direction at the current location. Represents direction vector With the guiding field vector dot product, and These represent the magnitudes of the vectors, This represents the weighting coefficient of the roughness resistance term. Indicates the location The surface roughness distribution value of the crack at that location This represents the natural exponential function.

[0026] In specific implementation, the freeze-thaw damage feature extraction module and the fracture network growth prediction module are executed sequentially. The freeze-thaw damage feature spectrum serves as the sole input to the fracture network growth dynamics model, and the potential penetration path set contains hundreds to thousands of candidate paths generated by random simulation. Optionally, the identification of main fracture intersection points is achieved by calculating the average length and number of fractures connected to each node and setting a weighted score threshold. The identification of the rock mass block isolation state is achieved by calculating the connected components of the three-dimensional network graph formed by the potential penetration path set and marking different regions. The fracture growth direction guiding field is obtained by normalizing the fracture boundary displacement vector field after Gaussian smoothing and spatial interpolation. In some embodiments, the growth starting point can be set as multiple main fracture intersection points to start multiple growth simulations in parallel, and the encounter distance threshold in the stopping condition is set to one-tenth of the average grain size of the rock mass.

[0027] Example 2: The coupling effect evolution correction module obtains historical temperature cycle data and geostress field distribution data of the target area from external data sources. It then fuses the historical temperature cycle data with the geostress field data to establish a stress intensity factor field coupling freeze-thaw cycles and geostress. From the set of potential penetration paths output by the fracture network growth prediction module, the spatial orientation vector and geometric parameters of each path are extracted one by one. The spatial orientation and geometric parameters of each path are substituted into the established stress intensity factor field to calculate the stress concentration degree of each extracted path under the current freeze-thaw-loading coupling effect. A frost heave damage accumulation model for cold-region rock masses is introduced. This model uses the duration of negative temperatures and the number of cycles in the historical temperature cycle data as the main inputs to calculate the frost heave fatigue damage value of the material along each path. The calculated stress concentration degree and frost heave fatigue damage value are weighted and superimposed according to preset weighting coefficients to calculate the coupling hazard index of each potential penetration path. Based on the coupling hazard index, paths with index values ​​exceeding a preset hazard threshold are selected to form a set of critical hazard paths. For each path in the set of critical hazard paths, an evolution stage determination result is assigned based on the growth rate of its stress concentration and the accumulation rate of its frost heave fatigue damage value. The evolution stage determination result includes the germination stage, the stable expansion stage, and the accelerated expansion stage.

[0028] When substituting the spatial orientation and geometry of each path into the stress intensity factor field to calculate the stress concentration, the specific implementation is as follows: A series of path feature points are set at preset fixed intervals along the extension trajectory of each potential through path. For each path feature point, the stress intensity factor value corresponding to its three-dimensional spatial position is queried in the stress intensity factor field. If the coordinates of the feature point are not exactly located on a discrete grid node of the stress intensity factor field, a cubic spline interpolation algorithm is used to calculate its stress intensity factor interpolation based on the values ​​of surrounding grid nodes. The calculated stress intensity factor values ​​at all path feature points on the path are extracted and arranged sequentially to form a local stress distribution sequence for the potential through path. The statistical characteristics of the local stress distribution sequence are calculated, including the sequence maximum value, sequence average value, and stress change gradient along the path orientation. Based on the spatial orientation vector of the potential through path, the angle between its principal extension direction and the direction of the maximum principal stress in the stress intensity factor field at each feature point is calculated, and the average angle is obtained. The calculated maximum value, average value, stress change gradient, and average included angle are input into a preset stress concentration coefficient calculation model. This model outputs a dimensionless coefficient that characterizes the overall stress concentration of the potential through path through a predefined weighted fusion algorithm. This dimensionless coefficient is the stress concentration of the path.

[0029] In practical implementation, the coupling effect evolution correction module obtains historical temperature cycle data and geostress field data of the target area from external data sources. It integrates the historical temperature cycle data and geostress field data to establish a stress intensity factor field for the coupling of freeze-thaw cycles and geostress. It extracts the spatial orientation and geometry of each path from the set of potential connecting paths one by one. It substitutes the spatial orientation and geometry of each path into the stress intensity factor field to calculate the stress concentration degree of each path under the current freeze-thaw-loading coupling effect. It introduces the frost heave damage accumulation model of cold-region rock mass and calculates the frost heave fatigue damage value along each path based on the historical temperature cycle data. It weights and superimposes the stress concentration degree and the frost heave fatigue damage value to obtain the coupling hazard index of each potential connecting path. Based on the coupling hazard index, the paths that exceed the preset hazard threshold are selected to form a set of critical hazard paths. For each path in the set of critical hazard paths, an evolution stage judgment result is assigned according to the growth rate of its stress concentration degree and the accumulation rate of its frost heave fatigue damage value. The evolution stage judgment result includes the nascent stage, the stable expansion stage, and the accelerated expansion stage. In some embodiments, historical temperature cycling data comes from a wireless temperature sensor network deployed in the target area and is stored in a time-series format. The geostress field data comes from a deep borehole geostress monitoring device and is processed by Kriging space interpolation to generate a continuous field. The frost heave damage accumulation model adopts the Miner linear accumulation rule based on fatigue accumulation damage theory. The allocation of evolution stage determination results is based on the comparison results of the stress concentration degree growth rate and the frost heave fatigue damage value accumulation rate with the preset grading threshold. The specific implementation of the frost heave damage accumulation model is as follows: Based on the Miner linear accumulation rule of fatigue cumulative damage theory, the model calculates damage by analyzing freeze-thaw cycle events in historical temperature cycling data. First, each complete freeze-thaw cycle is identified from the temperature time series data, i.e., the process of temperature dropping from positive to negative and then rising back to positive, and the duration of negative temperature, minimum temperature value, and temperature change amplitude of each cycle are extracted as damage calculation parameters. Then, based on the frost heave fatigue test data of rock mass materials, a basic damage value is assigned to each freeze-thaw cycle, which is positively correlated with the parameters of this cycle. Finally, the basic damage values ​​of all freeze-thaw cycles within the target time period are linearly accumulated to obtain the frost heave fatigue damage value along each path, which is used for the subsequent calculation of the coupled hazard index.

[0030] In practical implementation, the spatial orientation and geometry of each path are substituted into the stress intensity factor field to calculate the stress concentration degree of each path under the current freeze-thaw-loading coupling effect. Specifically, a series of path feature points are set at preset intervals along the extension direction of each potential through path. For each path feature point, the stress intensity factor value corresponding to its spatial location is queried in the stress intensity factor field. If the path feature point is not located on the precise grid node of the stress intensity factor field, its stress intensity factor interpolation is calculated using a cubic spline interpolation algorithm. The stress intensity factor values ​​at all path feature points are extracted to form the local stress distribution of the potential through path. The sequence calculates the statistical characteristics of the local stress distribution sequence, including the maximum value, average value, and standard deviation of the local stress distribution sequence. Based on the spatial orientation of the potential breakthrough path, the average angle between its principal extension direction and the direction of the maximum principal stress in the stress intensity factor field is calculated. The maximum value, average value, standard deviation, and average angle of the local stress distribution sequence are input into a preset stress concentration factor calculation model. The stress concentration factor calculation model outputs a dimensionless coefficient representing the overall stress concentration degree of the potential breakthrough path through a weighted fusion algorithm. Optionally, the preset interval distance is set to an integer multiple of the average grain size of the rock mass. The cubic spline interpolation algorithm uses the values ​​of eight adjacent nodes in the three-dimensional grid of the stress intensity factor field for interpolation calculation. The average angle is obtained by calculating the arithmetic mean of the inverse cosine of the dot product of the principal extension direction vector of the calculated path and the direction of the maximum principal stress at each characteristic point of the path in the stress intensity factor field. In some embodiments, the preset stress concentration factor calculation model adopts a linear weighted form, with the formula:

[0031] in: The dimensionless coefficient representing the overall stress concentration degree of a potential through path is the stress concentration degree. This represents the maximum value of the local stress distribution sequence. This represents the average value of the local stress distribution sequence. The standard deviation of the local stress distribution sequence. Reference value representing the uniaxial compressive strength of rock mass. This represents the average angle between the principal extension direction of the path and the direction of the maximum principal stress in the field. Represents the cosine function. , , , This represents the preset weighting coefficients that satisfy the normalization conditions. .

[0032] In practice, the coupling hazard index is calculated by linearly weighting the dimensionless coefficient of stress concentration with the frost heave fatigue damage value. The weighting coefficient reflects the relative contribution weights of freeze-thaw cycles and mechanical loading in crack propagation. The generation of the critical hazard path set is accomplished by comparing the coupling hazard index of each potential through-path with a preset hazard threshold. The allocation of the evolution stage determination results is based on the time series slope of the stress concentration growth rate and the frost heave fatigue damage value accumulation rate within the monitoring period. Optionally, the stress concentration growth rate is obtained by taking the difference between the dimensionless coefficient of stress concentration calculated in the current monitoring period and the previous monitoring period and dividing it by the monitoring time interval. The frost heave fatigue damage value accumulation rate is obtained by linearly fitting the damage value sequence output by the frost heave damage accumulation model to obtain the slope. The germination stage corresponds to both the growth rate and accumulation rate being lower than the first threshold. The stable propagation stage corresponds to either the growth rate or the accumulation rate exceeding the first threshold but lower than the second threshold. The accelerated propagation stage corresponds to either the growth rate or the accumulation rate exceeding the second threshold.

[0033] Example 3: See Figure 3 Raw temperature data of the target area, recorded in a time series over a historical monitoring period, is acquired from an external temperature sensor network. The daily maximum and minimum temperatures are extracted from this time-series data, and the daily temperature difference and the frequency of temperature fluctuations around zero degrees Celsius are calculated. These characteristic parameters are organized according to the spatial monitoring point locations to synthesize a temperature cycle characteristic matrix. Triaxial geostress component measurement data of the target area at multiple measuring points are acquired from an external geostress monitoring device. Combined with a regional geological structure model and spatial interpolation algorithms, a spatially continuous geostress field distribution model covering the target area is generated. A stress transformation function is established, taking the temperature cycle characteristic matrix and the geostress field distribution model as input, to calculate the composite stress value generated at different locations in the rock mass due to the superposition of frost heave force and geostress. The calculated composite stress values ​​are mapped according to a grid system consistent with the three-dimensional rock mass model to form a stress intensity factor field. The value of each grid point in the stress intensity factor field represents the stress intensity factor of its corresponding spatial location under freeze-thaw-loading coupling.

[0034] The establishment of the stress transformation function, using the temperature cycle characteristic matrix and the geostress field distribution model as inputs, calculates the composite stress value. Specifically, this includes: extracting the daily temperature difference and cumulative sub-zero temperature duration corresponding to the current grid point from the temperature cycle characteristic matrix; extracting the values ​​of the three principal stress components at the same grid point from the geostress field distribution model; constructing a frost heave force calculation subfunction, which uses the extracted daily temperature difference and cumulative sub-zero temperature duration as independent variables, and outputs the maximum frost heave force vector generated by the current grid point within a complete freeze-thaw cycle based on the theory of ice-water phase change volume expansion and rock pore structure parameters; constructing a geostress contribution calculation subfunction, which uses the extracted three principal stress components as inputs, and calculates the stress tensor caused by the regional geostress field at the grid point using elasticity formulas, combined with the anisotropic elastic parameters of the rock mass; and performing tensor superposition operation between the maximum frost heave force vector output by the frost heave force calculation subfunction and the stress tensor output by the geostress contribution calculation subfunction, according to their direction of action, to calculate the composite stress tensor at the grid point. Extract the maximum principal stress value from the synthetic stress tensor and output the maximum principal stress value as the synthetic stress value representing the stress intensity of the grid point.

[0035] In practical implementation, the process of establishing a stress intensity factor field coupling freeze-thaw cycles and geostress includes specific steps of acquiring data from external data sources and performing fusion calculations. The coupling effect evolution correction module acquires temperature time-series data of the target area within historical monitoring periods from an external temperature sensor network. It extracts the daily maximum and minimum temperatures from the temperature time-series data and calculates the daily temperature difference and the frequency of temperature fluctuations around zero degrees Celsius. These calculated characteristic parameters are organized and gridded according to the spatial coordinates of the sensor network to synthesize a spatially discrete temperature cycle characteristic matrix. Three-dimensional geostress component data from multiple borehole measuring points in the target area are acquired from an external geostress monitoring device. Combined with fault and bedding attitude information provided by the regional geological structure model, a spatial interpolation algorithm is used to expand the discrete measuring point data into a spatially continuous geostress field distribution model covering the entire target rock mass area. A stress transformation function is established, taking the temperature cycle characteristic matrix and the geostress field distribution model as inputs, to calculate the composite stress values ​​generated at different locations in the rock mass due to the superposition of frost heave force and geostress. The calculated synthetic stress values ​​are mapped one-to-one according to a grid system completely consistent with the three-dimensional digital model of the rock mass, forming a stress intensity factor field. The value of each grid point in the stress intensity factor field represents the stress intensity factor of the corresponding spatial position of the grid point in the three-dimensional rock mass spatial model under freeze-thaw-loading coupling. In some embodiments, the temperature sensor network transmits data wirelessly, the temperature time series data is stored in a database, the daily temperature difference is calculated as the daily maximum temperature minus the daily minimum temperature, the zero-degree fluctuation frequency is statistically represented by the number of times the temperature time series data crosses the zero-degree line, the geostress monitoring device uses the hydraulic fracturing method or the stress relief method for measurement, and the spatial interpolation algorithm uses the inverse distance weighting method or the kriging method. It can be understood that there is a mapping relationship between the row and column indices of the temperature cycle characteristic matrix and the grid coordinates of the three-dimensional rock mass spatial model, the geostress field distribution model stores the magnitude and direction of the three principal stress components for each grid point, and the stress intensity factor field is finally stored in a three-dimensional array data structure, the dimension of which is consistent with the grid division of the three-dimensional digital model of the rock mass.

[0036] In practical implementation, the stress transformation function calculates the synthetic stress value using the temperature cycle characteristic matrix and the geostress field distribution model as inputs. Specifically, it extracts the diurnal temperature range characteristic value and the cumulative duration of sub-zero temperature characteristic value corresponding to each grid point to be calculated from the temperature cycle characteristic matrix, and extracts the magnitude and direction tensors of the three principal stress components at the same grid point from the geostress field distribution model. A frost heave force calculation subfunction is constructed, using the extracted diurnal temperature range characteristic value and the cumulative duration of sub-zero temperature characteristic value as independent variables. Based on the theory of phase change volume expansion of pore water in rock mass and the porosity, saturation, and elastic modulus parameters of the rock, it outputs the maximum frost heave force vector generated by the currently processed grid point in a complete freeze-thaw cycle from the temperature cycle characteristic matrix. A geostress contribution calculation subfunction is also constructed, using the extracted three principal stress components as inputs, combined with the anisotropic elastic parameter matrix of the rock mass obtained through laboratory experiments, and calculating the total stress tensor caused by the regional geostress field at the grid point according to the generalized Hooke's law. The maximum frost heave force vector output by the frost heave force calculation subfunction is converted into an equivalent frost heave stress tensor. This frost heave stress tensor is then superimposed with the full stress tensor output by the geostress contribution calculation subfunction to calculate the composite stress tensor at each grid point. The maximum principal stress value is extracted from the composite stress tensor through eigenvalue decomposition and output as the composite stress value representing the stress intensity at each grid point. In some embodiments, the frost heave force calculation subfunction considers the gradual freezing process of pore water, and its output exhibits a non-linear positive correlation with the cumulative duration of sub-zero temperatures. The geostress contribution calculation subfunction transforms the principal stress components from the measuring point coordinate system to the global model coordinate system for calculation. The transformation of the frost heave stress tensor is achieved by dividing the frost heave force vector by the characteristic area it acts upon. Optionally, the formula for calculating the composite stress tensor is:

[0037] in: This represents the resultant stress tensor at the grid point. This represents the total stress tensor output by the subfunction for calculating the geostress contribution. Represents the feature area associated with the current grid cell. The maximum frost heave force vector output by the sub-function for calculating frost heave force is at the [missing information - likely a typo]. Components in each direction, The unit normal vector related to the direction of frost heave force is in the th... Components in each direction, subscript and Indicates the spatial coordinate direction index.

[0038] In practice, the calculation of the synthetic stress value is performed iteratively across all grid points of the three-dimensional rock mass spatial model. The temperature cycle characteristic matrix and the geostress field distribution model are used as global input data by the stress transformation function. The frost heave force calculation sub-function and the geostress contribution calculation sub-function are executed sequentially as built-in modules of the stress transformation function. It can be understood that the frost heave force calculation sub-function relies on the physical and mechanical parameters of the rock, which are pre-obtained through laboratory tests on rock core samples. The anisotropic elastic parameter matrix of the rock mass in the geostress contribution calculation sub-function is determined through sonic logging or triaxial compression tests, and the characteristic area... The stress intensity factor field is determined based on the element size of the 3D mesh model. Optionally, the final generation of the stress intensity factor field is accomplished by filling the calculated synthetic stress value at each mesh point into a corresponding 2D array. The data format of the stress intensity factor field supports rapid spatial querying and interpolation calculations in subsequent modules. In some embodiments, the stress transformation function is executed by a dedicated computational unit in the coupling evolution correction module. The calculation process employs parallel computing technology to improve the processing efficiency of large-scale 3D mesh models, with unit normal vectors... The direction is determined by the direction of the frost heave force vector and the anisotropic characteristics of the rock mass structure.

[0039] Example 4: The stability control list generation module obtains the evolution stage determination results of each path in the critical hazard path set, as well as the volume and spatial location attribute information of each isolated block in the rock mass block isolation status identifier. The spatial trajectory of each path in the critical hazard path set is spatially superimposed with the block boundaries defined in the rock mass block isolation status identifier to identify paths that directly penetrate the interior of the marked block or constitute block boundaries, thus causing instability of the critical block. These paths are defined as high-risk paths. An initial priority reinforcement level is assigned to each identified high-risk path. The assignment rule is based on the evolution stage determination results, where paths determined to be in the accelerated expansion stage correspond to the highest level, the stable expansion stage corresponds to the second highest level, and the nascent stage corresponds to the basic level. The initial priority reinforcement level is adjusted based on the volume attributes of the critical blocks affected by the high-risk path. A dynamic monitoring frequency is set for each high-risk path according to the adjusted final priority reinforcement level; the higher the priority reinforcement level, the higher the assigned dynamic monitoring frequency. The path numbers, spatial descriptions, adjusted priority reinforcement levels, and corresponding dynamic monitoring frequencies of all high-risk paths are compiled to form a structured list of stability control measures for fractured rock masses in cold regions.

[0040] The adjustment of the initial priority reinforcement level based on the volume of key blocks affected by high-risk paths is implemented as follows: In the three-dimensional block model defined by the rock mass block isolation status identifier, the rock mass blocks that intersect or are cut by each high-risk path in space are located by spatial position determination. The three-dimensional spatial volume of each located rock mass block is calculated. A volume-level adjustment mapping table is established, which predefines the priority reinforcement level increase amount corresponding to different volume ranges. For example, if the volume is greater than a certain threshold V1, the level is increased by 2 levels, and if the volume is between V2 and V1, the level is increased by 1 level. Based on the volume calculated for each rock mass block, the volume-level adjustment mapping table is queried to obtain the corresponding increase amount. The obtained increase amount is added to the initial priority reinforcement level of the high-risk path to obtain the adjusted priority reinforcement level. During the calculation process, it must be ensured that the adjusted level does not exceed the system's preset highest level threshold.

[0041] In practical implementation, the stability control list generation module cross-validates the set of critical hazardous paths and their evolution stage determination results with the rock mass block isolation status identifiers to generate a stability control list for cold-region fractured rock masses, including priority reinforcement levels and dynamic monitoring frequencies. The stability control list generation module obtains the evolution stage determination results for each path in the critical hazardous path set, as well as the volume and spatial location information of each isolated block in the rock mass block isolation status identifiers. It spatially overlays each path in the critical hazardous path set with the rock mass block isolation status identifiers to identify paths that traverse critical blocks or cause their instability, defining these paths as high-risk paths. An initial priority reinforcement level is assigned to each high-risk path, based on the evolution stage determination results: the accelerated expansion stage corresponds to the highest level, the stable expansion stage to the second-highest level, and the nascent stage to the basic level. The initial priority reinforcement level is adjusted based on the volume of the critical blocks affected by the high-risk paths. Dynamic monitoring frequencies are assigned to high-risk paths based on their priority reinforcement levels, with higher priority reinforcement levels receiving higher dynamic monitoring frequencies. The path information of all high-risk paths, the adjusted priority reinforcement levels, and the corresponding dynamic monitoring frequencies are then compiled to form the final stability control list for fractured rock masses in cold regions. In some embodiments, spatial overlay analysis is achieved through a geometric algorithm that calculates whether path segments intersect with the three-dimensional boundary polyhedra of rock mass blocks. The instability risk of key blocks is comprehensively assessed based on block volume and spatial location, and the dynamic monitoring frequency is quantified in units of "times / day" or "times / month".

[0042] In specific implementation, the method for adjusting the initial priority reinforcement level based on the volume of key blocks affected by high-risk paths involves locating rock mass blocks that intersect with or are cut by each high-risk path in the rock mass block isolation status identifier, and calculating the three-dimensional spatial volume of each located rock mass block. A volume-level adjustment mapping table is established, defining the increase in priority reinforcement level corresponding to rock mass blocks in different volume ranges. The corresponding increase is obtained by querying the volume-level adjustment mapping table based on the volume of each rock mass block. The increase is added to the initial priority reinforcement level of the high-risk path to obtain the adjusted priority reinforcement level, ensuring that the adjusted level does not exceed the system's set maximum level threshold. In some embodiments, the three-dimensional spatial volume of the rock mass block is calculated by summing the results after tetrahedral meshing of the block boundary polyhedrons. The volume-level adjustment mapping table is pre-installed in the system configuration file in the form of a lookup table. Optionally, the initial priority reinforcement level, the adjusted priority reinforcement level, and the maximum level threshold are all represented by integer levels, such as level 1 to level 5. See Table 1.

[0043] Table 1: Configuration table of a volume-level adjustment mapping table ; The adjusted priority reinforcement level is understandable. The calculation formula is:

[0044] in: This indicates the adjusted priority level for reinforcement. This indicates the initial priority reinforcement level assigned based on the results of the evolution stage determination. This indicates the increase in the priority hardening level obtained by querying the volume-level adjustment mapping table. This indicates the highest level threshold set by the system. This represents the minimum value function. The process of querying the volume-level adjustment mapping table is achieved by comparing the calculated rock mass block volume with the volume range defined in the mapping table. When the volume falls within a certain range, the corresponding upward adjustment amount for that range is applied.

[0045] In practical implementation, the dynamic monitoring frequency is set according to the principle that the higher the level, the higher the frequency. Different priority reinforcement level intervals correspond to different monitoring frequency baseline values. Optionally, the monitoring frequency baseline value is set through a system configuration file. For example, when the adjusted priority reinforcement level is level 1, the corresponding monitoring frequency is "1 time / month", level 2 corresponds to "1 time / week", and level 3 or above corresponds to "1 time / day". In some embodiments, the process of spatial overlay identification of high-risk paths uses a three-dimensional ray intersection detection algorithm to determine whether the path segment passes through the interior of the rock mass block. The volume information of each block in the rock mass block isolation status identifier is calculated and stored as an attribute when the identifier is generated. It can be understood that the allocation rule of the initial priority reinforcement level is a predefined mapping relationship: the accelerated expansion stage is mapped to the initial level 3, the stable expansion stage is mapped to the initial level 2, and the germination stage is mapped to the initial level 1.

[0046] See Figure 4 This is a scatter plot showing the relationship between block volume and the increase in reinforcement level, a data visualization chart used in the analysis of fracture evolution in jointed rock masses in cold regions. The scatter plot colors correspond to the increase in reinforcement level; for volumes ≤ 20 cubic meters, the increase is concentrated around 2.0. The red dashed line is a trend line, showing that the increase in reinforcement level rises rapidly after the volume exceeds 20 cubic meters, reflecting the rule that "the larger the block volume, the higher the increase in reinforcement priority." This chart is used in the stability analysis of fractured rock masses in cold regions to quantify the impact of block volume on reinforcement level, assisting in generating a decision-making basis for "priority reinforcement level." Larger blocks correspond to higher reinforcement priority for high-risk paths. This type of chart is commonly used in rock mass stability assessments in geological engineering, helping engineers quickly identify reinforcement priorities for high-risk areas by visualizing the correlation between block parameters and reinforcement strategies.

[0047] Example 5: After the system runs for a set time interval, such as after several freeze-thaw cycles, the data acquisition and processing flow is restarted. The latest multi-angle images and point cloud data of the target jointed rock mass are acquired through the multi-source spatiotemporal registration module and the freeze-thaw damage feature extraction module, and the latest freeze-thaw damage feature spectrum is obtained. From the latest freeze-thaw damage feature spectrum, the information on the actual observed fracture propagation paths between the previous prediction cycle and the current time is extracted, including newly added fracture segments and their spatial coordinates. The information on the actual observed fracture propagation paths is matched and compared with the set of potential through-paths previously predicted by the fracture network growth dynamics model to find predicted-observed path pairs with similar spatial locations and orientations. The difference between the actual observed propagation length and the growth length originally predicted by the model is calculated for each matched path pair, and the deviation angle between the average orientation of the actual observed path and the average orientation of the model-predicted path is also calculated. Based on the average length difference and the average orientation deviation angle of all matched path pairs, a set of model growth probability rule correction parameters is generated through predefined parameter update rules. These correction parameters aim to reduce future prediction bias. By using the generated model growth probability rules to correct parameters, the internal parameters of the rules controlling the growth direction probability, growth step size, or energy consumption in the fracture network growth dynamics model are iteratively updated to obtain an optimized fracture network growth dynamics model. This optimized model will be used for fracture path prediction processing in subsequent time periods.

[0048] In specific implementation, the system continuously optimizes the fracture network growth dynamics model. This process includes reacquiring the latest images and point cloud data of the target jointed rock mass at set time intervals; obtaining the latest freeze-thaw damage feature spectrum through fracture growth feature extraction; extracting the observed fracture propagation path information from the latest freeze-thaw damage feature spectrum; matching and comparing the observed fracture propagation path information with the set of potential penetration paths previously predicted by the fracture network growth dynamics model; calculating the difference between the actual observed length and the model-predicted length of each matched path, as well as the deviation angle between the actual observed direction and the model-predicted direction; generating model growth probability rule correction parameters based on the average of the difference and the deviation angle; and iteratively updating the growth probability rule parameters in the fracture network growth dynamics model using these correction parameters to obtain the optimized fracture network growth dynamics model for subsequent fracture path prediction processing. In some embodiments, the set time interval is determined according to the cycle of freeze-thaw cycles in cold regions, for example, triggering an optimization process after every five complete freeze-thaw cycles. The method of reacquiring data is consistent with the initial data acquisition process of the system, and the data structure of the latest freeze-thaw damage feature spectrum is completely consistent with the previously generated freeze-thaw damage feature spectrum. It can be understood that the information on the actual observed fracture propagation path refers to the spatial coordinate sequence of newly generated or extended fracture segments between the previous prediction period and the current time. The matching and comparison process is achieved by calculating the spatial hausdorff distance or fréchet distance between the actual observed path and each potential through path.

[0049] In practice, the difference between the actual observed length and the model predicted length of each matched path, as well as the deviation angle between the actual observed direction and the model predicted direction, are calculated. The model growth probability rule correction parameters are generated based on the average of the difference and the deviation angle. The actual observed length is obtained by accumulating the Euclidean distances between adjacent nodes in the actual observed fracture propagation path information. The model predicted length is obtained by accumulating the Euclidean distances between adjacent nodes of the corresponding path in the set of matched potential through-paths. The actual observed direction is obtained by calculating the global direction vector of the actual observed path, and the model predicted direction is obtained by calculating the global direction vector of the matched potential through-paths. The deviation angle is obtained by calculating the spatial angle between the actual observed direction vector and the model predicted direction vector. The generation of the model growth probability rule correction parameters relies on statistical analysis of the differences and deviation angles of all successfully matched path pairs. Optionally, the average difference between the actual observed length and the model predicted length of all matched path pairs is calculated. And the average angle of deviation between the actual observed trajectories and the model-predicted trajectories for all matched path pairs. Model growth probability rule correction parameters The calculation formula is:

[0050] in: This represents the correction parameters for the model growth probability rule; This represents the average difference between the actual observed length and the model predicted length for all matched path pairs; This represents the length normalization factor, which is a characteristic dimension of the rock mass model; This represents the average angle of deviation between the actual observed direction and the model-predicted direction for all matched path pairs; The angle normalization factor is usually taken as... radian; and These are preset weighting coefficients for the length difference term and the angle deviation term, respectively. In some embodiments, the length normalization factor... Take the length of the longest side of the outer envelope cuboid of the rock mass model, and the weighting coefficient. and Determined and satisfied based on historical data. Model growth probability rule correction parameters Used to quantitatively guide the adjustment of core parameters in the fracture network growth kinetics model.

[0051] In practical implementation, the growth probability rule parameters in the fracture network growth dynamics model are iteratively updated using the model growth probability rule correction parameters. These parameters include at least the direction weight parameters controlling the preference for random growth directions and the basic step size parameters controlling the growth step size. The iterative update employs an incremental adjustment method, adjusting the model growth probability rule correction parameters... Multiply by a learning rate coefficient Afterwards, according to The sign of the learning rate coefficient determines whether the direction weights and the base step size are adjusted in the same or opposite direction. This can be understood as the learning rate coefficient... It is a positive number less than 1, used to control the magnitude of each parameter update and avoid model oscillation. The updated direction weight parameters and basic step size parameters will be applied to the next fracture path prediction process. In some embodiments, if the model growth probability rule corrects the parameters... If positive, then proportionally reduce the component in the direction weight parameter that promotes crack growth towards the average direction of the prediction deviation angle, and slightly reduce the base step size parameter; if the model growth probability rule correction parameter If the value is negative, the component in the directional weight parameter that promotes the development of the crack in the direction of the average prediction deviation angle is increased proportionally, and the basic step size parameter is slightly increased.

[0052] See Figure 5This is a chart showing the relationship between rock mass volume, reinforcement level, and monitoring frequency, belonging to the category of data visualization charts for engineering stability analysis. The initial reinforcement level gradually increases with the increase in rock mass volume; the increase in reinforcement level after adjustment is more significant, reflecting the engineering logic of "the larger the volume, the higher the reinforcement priority." The monitoring frequency increases linearly with the increase in rock mass volume and is positively correlated with the reinforcement level; areas with larger volumes and higher reinforcement levels also have higher monitoring frequencies. This type of chart is commonly used in the stability control of jointed rock masses in cold regions, assisting engineers in determining the reinforcement priority and monitoring intensity for rock masses of different volumes, and is one of the visualization forms of a "stability control checklist."

[0053] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0054] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A visual detection-based analysis system for crack evolution of a jointed rock mass in cold regions, characterized in that, The system comprises: A multi-source space-time registration module is configured to acquire a multi-angle image set and an internal structure point cloud set of a target jointed rock mass surface in a cold region environment generated by freeze-thaw alternation within a set time, and perform registration and fusion processing on the multi-angle image set and the internal structure point cloud set in a time sequence to generate a multi-source data array after space-time registration; A freeze-thaw damage feature extraction module is configured to perform crack growth feature extraction processing on the multi-source data array after space-time registration to obtain a freeze-thaw damage feature spectrum, wherein the freeze-thaw damage feature spectrum comprises a crack surface roughness distribution, a crack branch node coordinate set and a crack boundary displacement vector field; A crack network growth prediction module is configured to perform crack path prediction processing on the freeze-thaw damage feature spectrum by using a crack network growth dynamics model to generate a set of potential through paths and a rock mass block isolation state identifier; An evolution correction module is configured to apply a freeze-thaw-loading coupling constraint specific to a cold region to the set of potential through paths to perform evolution path correction calculation, and output a set of corrected critical dangerous paths and evolution stage determination results; A stability control list generation module is configured to cross-verify and analyze the set of critical dangerous paths and the evolution stage determination results, and the rock mass block isolation state identifier to generate a cold region jointed rock mass stability control list comprising a priority reinforcement level and a dynamic monitoring frequency.

2. The visual detection-based analysis system for fracture evolution of cold region jointed rock mass according to claim 1, characterized in that, The crack growth feature extraction processing on the multi-source data array after space-time registration to obtain a freeze-thaw damage feature spectrum comprises the following processes: Separate image subsets and point cloud subsets corresponding to different freeze-thaw cycles from the multi-source data array after space-time registration; Perform joint feature extraction processing on each of the image subsets and point cloud subsets: Perform crack contour boundary identification on the image subsets to obtain a crack contour boundary pixel coordinate set; Perform crack spatial surface fitting on the point cloud subsets to obtain a crack three-dimensional spatial surface model; Map and match the crack contour boundary pixel coordinate set and the crack three-dimensional spatial surface model, generate a crack surface roughness distribution by calculating the projection deviation of the crack contour boundary pixel coordinates on the crack three-dimensional spatial surface model; Detect the mutation points of the surface curvature on the crack three-dimensional spatial surface model, and record the spatial coordinates of the curvature mutation points as a crack branch node coordinate set; Compare the deformation differences of the crack three-dimensional spatial surface model in adjacent freeze-thaw cycles, calculate the displacement of each point on the crack three-dimensional spatial surface model in the normal direction to form a crack boundary displacement vector field; Integrate the crack surface roughness distribution, the crack branch node coordinate set and the crack boundary displacement vector field to generate a complete freeze-thaw damage feature spectrum. 3.The visual detection based analysis system for fracture evolution of cold region jointed rock mass according to claim 2, characterized in that, The crack path prediction processing on the freeze-thaw damage feature spectrum by using the crack network growth dynamics model to generate a set of potential through paths and a rock mass block isolation state identifier comprises the following processing: Identify the main crack intersection points from the crack branch node coordinate set as growth starting points; According to the direction trend provided by the fracture boundary displacement vector field, a fracture growth direction guide field is set in the fracture network growth dynamics model; With the growth starting point as the origin, under the constraint of the fracture growth direction guide field, a multi-directional random growth simulation is performed according to a preset growth probability rule; During the multi-directional random growth simulation, the fracture surface roughness distribution is taken as a resistance factor to calculate the energy consumption of each random growth path; When the simulated growth fracture path meets other nodes in the fracture branch node coordinate set, or the energy consumption of the growth path exceeds a preset threshold, the growth simulation of the fracture path is stopped; All random growth paths that meet the stopping condition are recorded and integrated into a potential through-path set; According to the distribution of the potential through-path set in the three-dimensional space, it is judged whether the potential through-path set will divide the rock mass into independent blocks, and the independent blocks divided are spatially marked to generate a rock mass block isolation state identifier.

4. The visual detection-based analysis system for fracture evolution of cold region jointed rock mass according to claim 3, characterized in that, The potential through-path set is subjected to the constraint of the freeze-thaw-loading coupling effect unique to cold regions, and an evolution path correction calculation is performed to output a corrected key dangerous path set and an evolution stage determination result, which is achieved by the following means: Temperature cycle historical data and ground stress field data of the target area are obtained from an external data source, the temperature cycle historical data and the ground stress field data are fused, and a stress intensity factor field coupled with freeze-thaw cycles and ground stress is established; The spatial trend and geometric morphology of each path are extracted from the potential through-path set one by one; The spatial trend and geometric morphology of each path are substituted into the stress intensity factor field to calculate the stress concentration degree of each extracted path under the current freeze-thaw-loading coupling effect; A frost heaving damage accumulation model of the rock mass in the cold region is introduced, and the frost heaving fatigue damage value along each path is calculated based on the temperature cycle historical data; The stress concentration degree and the frost heaving fatigue damage value are weighted and superimposed to obtain a coupling danger index of each potential through-path; According to the coupling danger index, paths exceeding a preset danger threshold are selected to constitute a key dangerous path set; For each path in the key dangerous path set, an evolution stage determination result is assigned according to the growth rate of the stress concentration degree and the accumulation speed of the frost heaving fatigue damage value, and the evolution stage determination result includes the germination stage, the stable expansion stage and the accelerated expansion stage. 5.The visual detection based analysis system for fracture evolution of cold region jointed rock mass according to claim 4, characterized in that, The specific steps of obtaining temperature cycle historical data and ground stress field data of the target area from an external data source, fusing the temperature cycle historical data and the ground stress field data, and establishing a stress intensity factor field coupled with freeze-thaw cycles and ground stress include: Temperature time series data of the target area in a historical monitoring period are obtained from an external temperature sensor network, daily maximum and minimum temperatures are extracted from the temperature time series data, daily temperature difference and frequency of temperature fluctuation around zero are calculated, and a temperature cycle feature matrix is synthesized; Three-directional ground stress component data of the target area are obtained from an external ground stress monitoring device, and a spatially continuous ground stress field distribution model is generated in combination with a regional geological structure model; establish a stress conversion function, the stress conversion function takes the temperature cycle characteristic matrix and the ground stress field distribution model as input, calculates the synthetic stress value generated at different positions of the rock mass due to the superposition of frost heaving force and ground stress; map the synthetic stress value to the three-dimensional rock mass space model according to the preset grid, form a stress intensity factor field, and the value of each grid point in the stress intensity factor field represents the stress intensity factor of the corresponding spatial position of the grid point in the three-dimensional rock mass space model under the action of freeze-thaw and loading coupling. 6.The visual detection based cold region jointed rock mass fracture evolution analysis system according to claim 1, characterized in that, The cross-validation analysis of the key dangerous path set and its evolution stage determination result and the rock mass block isolation state identification generates a cold region fractured rock mass stability control list containing priority reinforcement levels and dynamic monitoring frequencies, which is completed through the following process: obtain the evolution stage determination result of each path in the key dangerous path set, and the volume and spatial position information of each isolated block in the rock mass block isolation state identification; superimpose each path of the key dangerous path set on the rock mass block isolation state identification, identify the paths that pass through the key blocks or may cause the key blocks to lose stability, and define them as high-risk paths; assign an initial priority reinforcement level to each high-risk path, and the assignment rule is based on the evolution stage determination result, wherein the accelerated expansion stage corresponds to the highest level, the stable expansion stage corresponds to the second highest level, and the germination stage corresponds to the basic level; adjust the initial priority reinforcement level in combination with the volume of the key blocks affected by the high-risk path. The larger the volume of the affected blocks, the higher the priority reinforcement level of the corresponding path. According to the priority reinforcement level of the high-risk path, set the dynamic monitoring frequency, the higher the priority reinforcement level, the higher the dynamic monitoring frequency. Path information, adjusted priority reinforcement level and corresponding dynamic monitoring frequency of all high-risk paths are summarized to form a final cold region fractured rock mass stability control list.

7. The visual detection-based analysis system for fracture evolution in cold region jointed rock mass according to claim 6, characterized in that, The specific implementation of adjusting the initial priority reinforcement level in combination with the volume of the key blocks affected by the high-risk path is as follows: locate the rock mass blocks intersected or cut by each high-risk path in the rock mass block isolation state identification; calculate the three-dimensional spatial volume of each located rock mass block; establish a volume-level adjustment mapping table, which defines the priority reinforcement level adjustment amount corresponding to rock mass blocks in different volume intervals; query the volume-level adjustment mapping table according to the volume of each rock mass block to obtain the corresponding adjustment amount; perform addition operation on the adjustment amount and the initial priority reinforcement level of the high-risk path to obtain the adjusted priority reinforcement level, and ensure that the adjusted level does not exceed the highest level threshold set by the system. 8.The visual detection based cold region jointed rock mass fracture evolution analysis system according to claim 5, characterized in that, The specific implementation of establishing a stress conversion function, the stress conversion function takes the temperature cycle characteristic matrix and the ground stress field distribution model as input, calculates the synthetic stress value generated at different positions of the rock mass due to the superposition of frost heaving force and ground stress, is as follows: extracting a daily temperature difference value and a cumulative duration of sub-zero temperature corresponding to each grid point from the temperature cycle feature matrix; extracting three principal stress components at each grid point from the geo-stress field distribution model; constructing a frost heaving force calculation sub-function, which takes the daily temperature difference value and the cumulative duration of sub-zero temperature as variables, and outputs a maximum frost heaving force vector generated by the current processing grid point in the temperature cycle feature matrix within one freeze-thaw cycle; constructing a geo-stress contribution calculation sub-function, which takes the three principal stress components as input, and calculates a stress tensor caused by the geo-stress field at the grid point in combination with the anisotropic elastic parameters of the rock mass; performing tensor superposition operation on the maximum frost heaving force vector and the stress tensor caused by the geo-stress field to calculate a combined stress tensor at the grid point; extracting a maximum principal stress value from the combined stress tensor as a combined stress value representing the stress intensity at the grid point. 9.The visual detection based cold region jointed rock mass fracture evolution analysis system according to claim 4, characterized in that, The implementation of the method is as follows: a series of path feature points are set at a preset interval along the extension direction of each potential through path; for each path feature point, the stress intensity factor value corresponding to its spatial position in the stress intensity factor field is queried, and if the path feature point is not located on the exact grid node of the stress intensity factor field, the stress intensity factor interpolation is calculated by a cubic spline interpolation algorithm; the stress intensity factor values at all path feature points are extracted to form a local stress distribution sequence of the potential through path; statistical characteristics of the local stress distribution sequence are calculated, including at least the sequence maximum value, the sequence average value, and the stress variation gradient along the path direction; the average included angle between the main extension direction of the potential through path and the maximum principal stress direction in the stress intensity factor field is calculated; the sequence maximum value, the sequence average value, the stress variation gradient, and the average included angle are input into a preset stress concentration coefficient calculation model, which outputs a dimensionless coefficient representing the overall stress concentration degree of the potential through path as the stress concentration degree through a weighted fusion algorithm. 10.The visual detection based cold region jointed rock mass fracture evolution analysis system according to claim 1, characterized in that, The method further includes a continuous optimization process for the crack network growth dynamics model, which includes: after a set time interval, the latest image and point cloud data of the target jointed rock mass are re-acquired, and the latest freeze-thaw damage feature spectrum is obtained through the crack growth feature extraction process; from the latest freeze-thaw damage feature spectrum, the actual observed crack propagation path information is extracted; the actual observed crack propagation path information is matched and compared with the set of potential through paths generated by the crack network growth dynamics model; the difference between the actual observed length and the model predicted length, and the deviation angle between the actual observed direction and the model predicted direction of each matched path are calculated; generate a model growth probability rule correction parameter based on the difference and the average of the deviation angles; update iteratively the growth probability rule parameter in the fracture network growth kinetics model by using the model growth probability rule correction parameter, and obtain an optimized fracture network growth kinetics model, which is used for subsequent fracture path prediction processing.

Images