Intelligent analysis method and system for stability of tunnel surrounding rock in tunnel excavation

Abstract

The invention relates to the technical field of tunnel construction safety, in particular to an intelligent analysis method and system for tunnel surrounding rock stability in tunnel excavation. The method comprises the following steps: acquiring geological data of surrounding rock of a tunnel, and performing tunnel excavation simulation to obtain tunnel excavation data; joint fissure detection is conducted on the basis of the tunnel excavation data, and joint fissure data are obtained; performing stress disturbance analysis according to the tunnel excavation data to obtain stress disturbance data; identifying a tunnel instability area according to the stress disturbance data; surrounding rock convergence displacement detection is conducted on the tunnel instability area, and surrounding rock convergence displacement data are generated; and predicting the surrounding rock falling probability based on the surrounding rock convergence displacement data. According to the method, the stable state of the surrounding rock is dynamically and accurately judged on the basis of a tunnel construction safety technology, so that the accuracy and timeliness of tunnel surrounding rock instability risk identification and early warning are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of tunnel construction safety technology, and in particular to an intelligent analysis method and system for the stability of surrounding rock in tunnel excavation. Background Technology

[0002] Traditional tunnel excavation rock stability analysis relies heavily on manual experience and periodic on-site monitoring. Insufficient and continuous data collection leads to significant subjectivity and lag in rock stability assessment. The lack of refined identification and quantitative analysis of joints and fissures hinders accurate capture of rock fragmentation characteristics and structural surface activation, impacting the accuracy of instability risk assessment. Limited ability to analyze the dynamic response to stress disturbances and rock convergence displacement makes it difficult to achieve real-time, accurate location of instability zones and prediction of their evolution trends. Identification of potential unstable blocks and wedge-shaped blocks relies on simplified models, lacking multi-dimensional comprehensive judgment, resulting in unscientific instability risk level assessments. The absence of an intelligent integration platform leads to data silos, preventing efficient linkage between analysis modules and affecting the comprehensiveness and timeliness of rock stability analysis. Ultimately, this restricts the improvement of tunnel excavation safety management and risk control. Summary of the Invention

[0003] Therefore, it is necessary for the present invention to provide an intelligent analysis method and system for the stability of the surrounding rock of a tunnel during tunnel excavation, so as to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, an intelligent analysis method for the stability of surrounding rock in tunnel excavation includes the following steps: Step S1: Obtain geological data of the surrounding rock of the tunnel and perform tunnel excavation simulation to obtain tunnel excavation data; based on the tunnel excavation data, perform joint and fracture detection to obtain joint and fracture data; Step S2: Perform stress disturbance analysis based on tunnel excavation data to obtain stress disturbance data; identify unstable areas of the tunnel based on stress disturbance data; detect surrounding rock convergence displacement in unstable areas of the tunnel to generate surrounding rock convergence displacement data; predict the probability of surrounding rock detachment based on surrounding rock convergence displacement data. Step S3: Based on the surrounding rock spalling probability and joint and fracture data, perform structural surface activation detection to obtain structural surface activation data; based on the structural surface activation data, identify the structural surface shear slip trend to obtain structural surface shear slip data. Step S4: Identify potentially unstable blocks based on structural shear slip data to obtain potential unstable block data; perform wedge block detection based on potential unstable block data to obtain wedge block data; assess the instability risk level based on the wedge block data; perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

[0005] This invention utilizes tunnel excavation simulation combined with geological data for joint and fracture detection, enabling refined identification of surrounding rock structural characteristics and improving the quantitative analysis accuracy of joint and fracture distribution, orientation, and density, thus providing fundamental support for subsequent stability assessment. Secondly, through stress disturbance analysis and instability zone identification, the spatial changes in the stress state of the surrounding rock can be dynamically monitored, achieving early identification and regional location of instability risks and effectively improving early warning capabilities. Furthermore, surrounding rock convergence displacement detection and detachment probability prediction can reveal the deformation evolution trend of the surrounding rock, improving the accuracy of predicting deformation-induced instability. Structural surface activation detection and shear slip trend identification can comprehensively assess the response mechanism and potential slip risk of structural surfaces, enhancing the ability to identify instability causes. In the identification of potential unstable blocks and wedge-shaped block detection, by integrating slip data and spatial geometric features, the ability to identify and model complex block structures is improved, making the judgment of instability modes more scientific and reliable. Finally, by combining multi-source data to assess the instability risk level and conduct stability analysis, a closed-loop analysis of the entire process from data acquisition to risk assessment is achieved. This significantly overcomes the problems of strong subjectivity, poor timeliness, and simplified models in traditional methods, and improves the intelligence, systematization, and refinement of surrounding rock stability determination, providing strong support for risk prevention and safety decision-making in tunnel construction.

[0006] Preferably, step S1 specifically includes: Step S11: Obtain geological data of the surrounding rock of the tunnel and conduct tunnel excavation simulation to obtain tunnel excavation data; Step S12: Collect joint and fracture images based on tunnel excavation data; Step S13: Identify the layered joint structure based on the joint image; Step S14: Identify columnar joint fracture structures based on joint fracture images; Step S15: Integrate the layered joint fracture structure and the columnar joint fracture structure to obtain joint fracture data.

[0007] This invention, by introducing tunnel excavation simulation and joint and fracture image recognition technology, achieves high-precision automatic extraction and classification analysis of surrounding rock structural features, significantly improving the intelligence and scientific level of surrounding rock stability assessment. By acquiring and analyzing joint and fracture images, the subjective errors and delays caused by traditional reliance on manual on-site judgment can be eliminated, ensuring the continuity and comprehensiveness of joint and fracture data collection. At the image recognition level, a multi-feature recognition algorithm is used to extract layered and columnar joint and fracture structures separately, enabling accurate differentiation and quantification of various typical fracture types, providing a reliable structural foundation for subsequent slip trend analysis and block modeling. By integrating data from different fracture structure types and establishing a unified joint and fracture information model, not only is the data structuring and standardization improved, but the comprehensive perception of surrounding rock fracture characteristics and joint structure evolution trends is also effectively enhanced, thereby significantly improving the accuracy and reliability of structural surface activation prediction and instability risk assessment. Overall, this method enhances the ability to systematically identify and model the microscopic morphology and macroscopic distribution patterns of joints and fissures, promotes the transformation of tunnel surrounding rock stability analysis from experience-driven to data-driven, and from qualitative judgment to vectorized identification, lays a key data foundation for intelligent analysis of surrounding rock stability, and effectively solves the technical bottlenecks of low fissure identification accuracy, lagging data updates, and distorted risk assessment in traditional methods.

[0008] Preferably, step S13 specifically includes: Step S131: Perform image preprocessing based on the joint and fracture image to obtain the joint and fracture image to be processed; Step S132: Identify the fracture edges based on the joint fracture image to be processed, and obtain fracture edge data; Step S133: Perform linear detection based on the crack edge data to obtain linear crack data; Step S134: Calculate the fracture spacing based on the linear fracture data; Step S135: Calculate the parallelism of the fractures based on the linear fracture data; Step S136: Determine the layered joint structure based on the fracture spacing and the degree of fracture parallelism.

[0009] This invention significantly improves the accuracy and robustness of joint and fracture image analysis through image preprocessing and fracture edge recognition technology, ensuring a clearer, more continuous, and more interference-resistant data foundation for subsequent fracture structure identification. A linear detection mechanism is introduced based on edge recognition to effectively extract the geometric morphological features of joints and fractures, overcoming the subjective errors in fracture orientation identification inherent in traditional manual observation. Precise calculation of fracture spacing quantifies the degree of joint development, enabling a more detailed assessment of the integrity and fragmentation of the surrounding rock. Simultaneously, the calculation of fracture parallelism allows for the quantitative expression of the geometric patterns of fracture structure sequences, providing crucial information for the accurate identification of structural surface features. Identifying layered joint and fracture structures by integrating fracture spacing and parallelism not only achieves automatic determination of typical joint patterns but also enhances the analytical capabilities for the evolutionary characteristics of surrounding rock microstructures. Overall, this method overcomes the technical bottlenecks of traditional methods, such as coarse identification, unclear structural classification, and weak quantification capabilities, through structured processing and parameterized extraction of fracture images. It significantly improves the systematicness, accuracy, and practicality of joint identification, laying a solid structural data foundation for subsequent instability mechanism modeling and risk assessment, and helping to build an efficient, intelligent, and accurate surrounding rock stability analysis system.

[0010] Preferably, step S14 specifically includes: Step S141: Identify the vertical crack edge based on the joint and crack image to obtain the vertical crack contour data; Step S142: Identify polygonal cross sections based on vertical crack profile data to obtain polygonal crack profile data; Step S143: Calculate the profile closure degree based on the polygonal crack profile data; Step S144: Perform planar connectivity analysis based on the contour closure to obtain planar connectivity fracture data; Step S145: Calculate the spatial density based on the planar connected fracture data; Step S146: Determine and record the columnar joint fracture structure based on the spatial density.

[0011] This invention accurately captures the structural features of columnar joints in images through vertical fracture edge recognition and polygonal cross-section extraction, solving the problem of inaccurate understanding of two-dimensional projection of three-dimensional structures in traditional recognition processes. Polygonal fitting of fracture contours makes the fracture morphology structurally calculable, facilitating subsequent topological analysis and geometric feature quantification. The introduction of contour closure calculation effectively determines whether fractures form independent units, improving the accuracy of planar structure recognition. Planar connectivity analysis based on closure enables the identification of the continuity of columnar joint structures, thus supporting the three-dimensional reconstruction and propagation path determination of fracture connectivity networks. Calculating spatial density based on planar connected fracture data not only assesses the development intensity of columnar joints in space but also provides quantitative indicators for the mechanical strength and fracturing degree of surrounding rock. Finally, determining the columnar joint fracture structure based on spatial density overcomes the subjective limitations of traditional visual experience-based judgment, achieving automated, quantitative, and intelligent columnar structure recognition. Overall, this method establishes a high-precision, full-process analysis path for identifying columnar joint structures through step-by-step reasoning and parameter quantification of the fracture structure identification logic. This significantly enhances the scientific rigor and real-time performance of surrounding rock structure identification, providing solid data support and a structural foundation for subsequent structural surface activation analysis, instability risk identification, and intelligent decision-making.

[0012] Preferably, step S2 specifically includes: Step S21: Perform surrounding rock evolution analysis based on tunnel excavation data to obtain surrounding rock evolution data; Step S22: Calculate the principal stress variation rate based on the surrounding rock evolution data; Step S23: Identify stress concentration areas in the surrounding rock based on the principal stress variation rate; Step S24: Calculate the stress disturbance intensity based on the stress concentration area of ​​the surrounding rock to obtain stress disturbance data; Step S25: Identify the unstable areas of the tunnel based on stress disturbance data; Step S26: Detect the surrounding rock convergence displacement in the unstable area of ​​the tunnel and generate surrounding rock convergence displacement data; Step S27: Predict the probability of rockfall based on the surrounding rock convergence displacement data.

[0013] This invention, by introducing surrounding rock evolution analysis, can dynamically track the deformation and mechanical response trends of the surrounding rock during excavation, providing fundamental evolutionary characteristics for subsequent stability analysis. The calculation of principal stress variation rate enhances the ability to identify areas sensitive to dynamic stress disturbances, offering stronger timeliness and early warning capabilities compared to traditional static stress field analysis. Identification of stress concentration areas allows for precise pinpointing of local high-risk areas, facilitating differentiated support and key protection designs. Further quantification of stress disturbance intensity based on stress concentration areas lays the foundation for quantitative analysis and dynamic monitoring of instability factors. The identification of tunnel instability areas overcomes the limitations of previous reliance on manual inspections and deformation threshold judgments, enabling automated calibration of high-risk areas. Combined with surrounding rock convergence displacement detection, a logical chain of excavation disturbance—deformation response—instability trend is constructed, enabling continuous acquisition and evaluation of displacement characteristics in deformation concentration areas. Finally, the probability of surrounding rock detachment is predicted based on surrounding rock convergence displacement data, enhancing the quantitative early warning capability for local surrounding rock instability risks and avoiding safety hazards caused by delayed deformation judgment. Overall, this method constructs a full-chain stability analysis path with excavation disturbance as the main line, evolution trend as the clue, and deformation response as the support. It overcomes the shortcomings of traditional methods in terms of timeliness, continuity, quantification, and intelligent linkage, and significantly improves the accuracy and effectiveness of identifying unstable areas of surrounding rock and providing risk warnings.

[0014] Preferably, step S25 specifically includes: Step S251: Calculate the peak shear stress based on the stress disturbance data; Step S252: Obtain the yield strength of the surrounding rock; identify potential unstable elements based on the yield strength and peak shear stress of the surrounding rock to obtain potential unstable element data; Step S253: Perform laser cross-section scanning based on the potential unstable element data to obtain laser cross-section data; Step S254: Calculate the deformation convergence based on the laser cross-section data; Step S255: Identify high-deformation section areas based on deformation convergence; perform rock mass microseismic monitoring based on the high-deformation section areas to obtain rock mass microseismic data; Step S256: Perform vibration frequency statistics based on the rock mass microseismic data to obtain high-frequency rock mass microseismic data; Step S257: Determine the tunnel instability zone based on high-frequency rock mass microseismic data.

[0015] This invention achieves highly sensitive identification of abnormal stress zones in surrounding rock through statistical analysis of peak shear stress, enhancing the ability to initially judge instability causes. Combining the yield strength of the surrounding rock with the identification of potential instability units enables precise targeting of high-risk areas based on mechanical critical states, significantly improving the physical and scientific rigor of the identification. The introduction of laser cross-section scanning technology acquires high-precision deformation information, providing data support for the continuous quantification of deformation trends. The calculation of deformation convergence reflects the aggregation characteristics of cross-sectional deformation, effectively identifying the weakest structural areas most affected by disturbances. Further identification of high-deformation cross-sectional areas enables spatial focusing of deformation, providing precise target areas for microseismic monitoring deployment and improving the efficiency of sensor resource allocation. Rock mass microseismic monitoring can reflect microscopic instability precursors such as joint slippage or structural loosening in real time. Statistical analysis of vibration frequencies further filters out high-frequency abnormal signals, effectively avoiding false alarms. Based on high-frequency rock mass microseismic data to determine instability areas, a multi-dimensional linkage identification logic of "mechanical criteria—deformation trend—dynamic disturbance—energy release" is constructed, realizing a fully coupled identification path from static discrimination to dynamic early warning. Overall, this method establishes an information channel between static indices, deformation information, and microseismic data, improving the accuracy, timeliness, and intelligence of tunnel instability zone identification, and providing efficient and reliable technical support for dynamic assessment of surrounding rock stability and risk prevention and control.

[0016] Preferably, step S26 specifically includes: Step S261: Perform crown region detection on the unstable area of ​​the tunnel to obtain crown region data; Step S262: Identify the initial profile of the cross-section based on the data of the arch area; Step S263: Continuously monitor the prism target based on the data of the dome area to obtain prism target monitoring data; Step S264: Calculate the radial displacement of the cross section based on the prism target monitoring data; calculate the convergence rate of the cross section based on the prism target monitoring data; perform crack propagation analysis based on the radial displacement and convergence rate of the cross section to obtain the crack propagation coordinates; Step S265: Map the fracture propagation coordinates to the initial profile of the cross section to determine the degree of surrounding rock convergence displacement and obtain surrounding rock convergence displacement data.

[0017] This invention achieves precise positioning of key structural components through automatic detection of the arch region within the tunnel instability area, enhancing the targeted nature of monitoring potential deformation and instability of the arch. Identifying the initial cross-sectional contour establishes an accurate benchmark for subsequent deformation analysis, facilitating precise measurement of deformation offset. Continuous monitoring with a prism target enables high-precision, continuous deformation data acquisition, significantly improving data timeliness and reliability. Calculating the radial displacement and convergence rate of the cross-section based on monitoring data forms a dynamic tracking capability for surrounding rock deformation trends, providing quantitative support for fracture propagation analysis. Extracting fracture propagation coordinates enables a refined characterization of the fracture evolution process, deepening the understanding of the surrounding rock failure mechanism. Mapping the propagation coordinates to the initial contour yields surrounding rock convergence displacement data, reflecting the spatial coupling relationship between actual deformation offset and fracture development, providing precise evidence for quantitative assessment of surrounding rock stability. Overall, this method achieves full-process linkage from structural identification, continuous monitoring, dynamic calculation to spatial mapping, enhancing the accuracy and real-time performance of surrounding rock deformation evolution identification, and significantly improving the scientific and intelligent level of tunnel instability early warning.

[0018] Preferably, step S3 specifically includes: Step S31: Calculate the joint density index based on the joint and fracture data; calculate the joint intersection angle based on the joint and fracture data; Step S32: Determine the highly intersecting cutting area based on the joint density index and joint intersection angle; Step S33: Map the probability of surrounding rock detachment to the highly interlaced cutting area to obtain the area with high probability of surrounding rock detachment; Step S34: Identify the characteristics of broken rock blocks based on areas with high probability of rock fragmentation, and obtain the characteristic data of broken rock blocks; Step S35: Determine the degree of structural surface activation based on the characteristic data of the fractured rock blocks, and obtain structural surface activation data; Step S36: Identify the shear slip trend of the structural surface based on the structural surface activation data to obtain the structural surface shear slip data.

[0019] This invention, through in-depth calculation of joint and fracture data, obtains the joint density index and intersection angle, which can quantitatively describe the distribution characteristics and structural complexity of joints, providing a precise basis for subsequent structural identification. Combining density and intersection angle to determine highly intersecting cutting areas helps to accurately identify potential structural weak zones. Mapping the probability of surrounding rock detachment to this area can effectively locate high-risk areas, achieving centralized early warning of unstable zones. Identifying the characteristics of broken rock blocks can further reveal the geometric morphology and mechanical weakening degree of local broken structures, providing a data foundation for structural activation degree analysis. Determining the activation degree of structural surfaces based on the characteristics of broken rock blocks can reveal the true response state and potential displacement trend of structural surfaces in the surrounding rock. Finally, by identifying the shear slip trend of structural surfaces through activation data, the path and direction of structural surface instability and slippage can be accurately predicted, thereby enhancing the understanding and control of the surrounding rock instability mechanism. The overall process realizes a full-chain intelligent analysis from structural complexity identification to risk mapping, and from structural activation judgment to shear slip trend identification, significantly improving the systematicness, scientificity, and forward-looking nature of tunnel surrounding rock stability analysis.

[0020] Preferably, step S4 specifically includes: Step S41: Calculate the gravity direction vector based on the structural shear slip data; calculate the slip direction vector based on the structural shear slip data; calculate the slip surface inclination angle based on the gravity direction vector and the slip direction vector; identify potential unstable blocks based on the slip surface inclination angle to obtain potential unstable block data; Step S42: Identify geometric closure based on potential unstable block data; determine wedge-shaped blocks based on geometric closure to obtain wedge-shaped block data; Step S43: Identify the degree of voiding of the structural face based on the wedge block data to obtain voiding data of the structural face; assess the instability risk level based on the voiding data of the structural face; Step S44: Perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

[0021] This invention, through the calculation of the vectors of slip direction and gravity direction, can accurately determine the dip angle of the slip surface, thereby identifying blocks with potential instability risks and achieving a quantitative assessment of the shear instability trend of the structural surface. Furthermore, by identifying wedge-shaped blocks through geometric closure judgment, it can pinpoint key blocks with poor structural integrity and prone to slippage based on geometric structural features, improving the accuracy of identifying unstable units in complex structures. Based on this, analyzing the degree of openness of the structural surface can reveal in depth the weakening of structural surface support and its actual impact on stability, helping to accurately locate areas with missing support or stress release in local structures. Subsequently, an instability risk level assessment can be conducted, systematically integrating multi-source information such as structural geometric features, mechanical state, and support conditions to form a comprehensive basis for judging the degree of risk. Finally, the surrounding rock stability data generated through stability analysis provides reliable quantitative support for tunnel excavation safety management and zonal reinforcement. The overall process realizes a closed-loop analysis chain from slip trend identification, block structure analysis, openness feature analysis to risk level assessment, comprehensively improving the scientific, systematic, and real-time nature of tunnel surrounding rock stability judgment.

[0022] Preferably, this specification also provides an intelligent analysis system for the stability of surrounding rock in tunnel excavation, used to execute the intelligent analysis method for the stability of surrounding rock in tunnel excavation as described above. The intelligent analysis system for the stability of surrounding rock in tunnel excavation includes: The joint and fracture detection module is used to acquire geological data of the surrounding rock of the tunnel and to simulate tunnel excavation to obtain tunnel excavation data; based on the tunnel excavation data, joint and fracture detection is performed to obtain joint and fracture data. The surrounding rock spalling probability prediction module is used to perform stress disturbance analysis based on tunnel excavation data to obtain stress disturbance data; identify unstable areas of the tunnel based on the stress disturbance data; detect the surrounding rock convergence displacement of the unstable areas of the tunnel to generate surrounding rock convergence displacement data; and predict the probability of surrounding rock spalling based on the surrounding rock convergence displacement data. The structural surface shear slip trend identification module is used to detect structural surface activation based on the surrounding rock detachment probability and joint and fracture data to obtain structural surface activation data; and to identify structural surface shear slip trend based on structural surface activation data to obtain structural surface shear slip data. The surrounding rock stability analysis module is used to identify potentially unstable blocks based on structural surface shear slip data, obtain potential unstable block data; detect wedge-shaped blocks based on the potential unstable block data, obtain wedge-shaped block data; assess the instability risk level based on the wedge-shaped block data; and perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

[0023] The present invention relates to an intelligent analysis system for the stability of surrounding rock in tunnel excavation. This system can implement any of the intelligent analysis methods for the stability of surrounding rock in tunnel excavation according to the present invention. It is used as a medium for combining the operation and signal transmission between various modules to complete the intelligent analysis method for the stability of surrounding rock in tunnel excavation. The internal modules of the system cooperate with each other to achieve dynamic and accurate determination of the stability state of the surrounding rock, thereby significantly improving the accuracy and timeliness of identification and early warning of tunnel surrounding rock instability risks. Attached Figure Description

[0024] Other features, objects, and advantages of the invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings: Figure 1 This is a schematic diagram of the steps of an intelligent analysis method for the stability of surrounding rock in tunnel excavation according to the present invention. Figure 2 This is a detailed flowchart of step S1 in the present invention; Figure 3 This is a detailed flowchart of step S13 in the present invention; The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

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

[0026] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0027] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0028] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides an intelligent analysis method for the stability of surrounding rock in tunnel excavation, the method comprising the following steps: Step S1: Obtain geological data of the surrounding rock of the tunnel and perform tunnel excavation simulation to obtain tunnel excavation data; based on the tunnel excavation data, perform joint and fracture detection to obtain joint and fracture data; In this embodiment, geological parameters of the tunnel surrounding rock are collected using on-site surveying methods such as drilling, seismic exploration, and ground-penetrating radar. These parameters include, but are not limited to, the bedding structure, spatial distribution of joints and fractures, rock type, and mechanical properties (such as elastic modulus, Poisson's ratio, and compressive strength). Based on the collected geological data, a three-dimensional digital geological model is constructed, and CAD software or specialized geological modeling software is used to digitally represent the tunnel structure and geological characteristics of the surrounding rock. Subsequently, based on numerical calculation methods, finite element analysis software (such as ABAQUS or ANSYS) is used to simulate tunnel excavation. The excavation cross-section dimensions, excavation sequence, and geostress field parameters are set. These parameters include the magnitude of the initial principal stress (e.g., the maximum principal stress is 20 MPa, and the minimum principal stress is 10 MPa) and its direction (angle relative to the surface direction), simulating the disturbance of the surrounding rock stress field caused by tunnel excavation. During the excavation simulation, the elements at the tunnel boundary are refined according to the mesh division to ensure that the spatial resolution of stress changes is within 0.1 meters. Based on the stress and deformation field data output by simulation, image processing algorithms and fault recognition technology are used to automatically detect joints and fractures in the surrounding rock. The detection methods include the Canny edge detection algorithm and morphological processing. Combined with neighborhood connectivity analysis, the spatial distribution characteristics of the fractures are determined, and the length, orientation and density of the joints and fractures are calculated. Finally, the joint and fracture data are output, including information such as joint location coordinates, length range (e.g., 0.5 meters to 3 meters), and cracking angle range (0°-180°).

[0029] Step S2: Perform stress disturbance analysis based on tunnel excavation data to obtain stress disturbance data; identify unstable areas of the tunnel based on stress disturbance data; detect surrounding rock convergence displacement in unstable areas of the tunnel to generate surrounding rock convergence displacement data; predict the probability of surrounding rock detachment based on surrounding rock convergence displacement data. In this embodiment, based on the tunnel excavation data obtained in step S1, a detailed analysis of stress disturbance within the excavation area is performed. Numerical simulation software is used to extract the principal stress and shear stress changes caused by excavation. When calculating stress disturbance data, the difference between the baseline stress field and the post-disturbance stress field is set. The finite difference method is used to calculate the stress change rate at each point, and a critical disturbance threshold of 0.2 MPa is set. Areas exceeding this value are marked as areas with significant stress disturbance. Based on this data, a spatial clustering algorithm (such as DBSCAN) is used to identify continuous areas with large stress disturbances, thus identifying potential tunnel instability areas. For the identified instability areas, a high-precision laser rangefinder and displacement sensor are installed to detect the surrounding rock convergence displacement. The sensor sensitivity reaches 0.01 mm, and the sampling frequency is set to once per hour to ensure continuous acquisition of dynamic displacement data. The collected displacement data is processed using a time-series analysis method to calculate the surrounding rock convergence rate. The rate threshold is set to 0.05 mm / hour for subsequent statistical modeling of the surrounding rock detachment probability. Probabilistic statistical models (such as Poisson process or Bayesian inference method) are used to analyze the surrounding rock convergence displacement data. Combined with convergence rate and historical rockfall event data, the probability value of surrounding rock falloff is output. The value range is limited to 0 to 1, and a value greater than 0.7 is considered high risk.

[0030] Step S3: Based on the surrounding rock spalling probability and joint and fracture data, perform structural surface activation detection to obtain structural surface activation data; based on the structural surface activation data, identify the structural surface shear slip trend to obtain structural surface shear slip data. In this embodiment, the activation status of the structural surface is detected by combining the surrounding rock detachment probability obtained in step S2 and the joint and fracture data obtained in step S1. The activation detection is based on structural mechanics analysis. The stress concentration factor on the structural surface is calculated, and mechanical parameters such as shear modulus (G) and normal stress (σn) are selected. The specific calculation formula is τ / (c+σntanφ), where τ is the shear stress, c is the cohesion (taken as 0.5MPa), and φ is the internal friction angle (taken as 35°). The activation degree of the structural surface is numerically assigned, ranging from 0 to 1. A value greater than 0.6 indicates that the structural surface is in an activated state. Based on the activation data, the finite difference method is used to identify the shear slip trend of the structural surface, analyzing the slip direction vector and velocity. Regions with a slip velocity exceeding 0.1 mm / day are determined to have a shear slip trend, and the slip direction is calculated based on the principal stress field direction. Finally, the structural surface activation data and shear slip data are output, including the slip velocity vector and activation degree score.

[0031] Step S4: Identify potentially unstable blocks based on structural shear slip data to obtain potential unstable block data; perform wedge block detection based on potential unstable block data to obtain wedge block data; assess the instability risk level based on the wedge block data; perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

[0032] In this embodiment, the structural shear slip data obtained in step S3 is used to identify potentially unstable blocks using a three-dimensional geometric analysis method. The block stability is calculated based on the slip surface inclination angle, and the safety factor of the unstable block is calculated using limit equilibrium theory. A threshold of 1.2 is set, and blocks with a safety factor lower than this value are marked as potentially unstable blocks. The block geometry is obtained using point cloud scanning technology, and geometric closure is assessed in conjunction with structural surface data. Blocks with a closure gap width of less than 5 mm are considered well-closed, while those with a gap width greater than 5 mm are considered wedge-shaped blocks. The degree of openness of the structural surfaces of wedge-shaped blocks is quantitatively assessed using ultrasonic testing and void scanning technology. Areas with a porosity greater than 10% are identified as high-risk areas for openness. Based on the block safety factor and degree of openness, a standard for classifying instability risk levels is established, divided into low, medium, and high levels, with critical values ​​of 1.5 and 1.0, respectively. Based on the risk level data, a comprehensive surrounding rock stability analysis is performed using statistical methods, outputting a final surrounding rock stability data report, including a block distribution map, a risk level heat map, and a stability score.

[0033] Preferably, step S1 specifically includes: Step S11: Obtain geological data of the surrounding rock of the tunnel and conduct tunnel excavation simulation to obtain tunnel excavation data; In this embodiment, geological parameters of the tunnel surrounding rock were systematically collected through various means, including borehole sampling, ground-penetrating radar scanning, and on-site observation. These parameters included rock type, mineral composition, stratum thickness, and the distribution characteristics of joints and fractures. A high-precision 3D laser scanner was deployed on-site to scan the tunnel excavation face and acquire high-resolution point cloud data. The point cloud spacing was controlled within 1 mm to ensure the capture of subtle geological features. Based on the on-site geological parameters and point cloud data, the tunnel excavation process was simulated using the finite element method. The excavation cross-section was set to be 10 meters wide and 8 meters high, with an excavation depth of 30 meters. The initial ground stress values ​​were based on on-site measured data, with a maximum principal stress of 25 MPa and a directional deviation angle of 30°, and a minimum principal stress of 12 MPa and a directional deviation angle of 120°. The static loading method was used during the excavation simulation, gradually removing the mesh elements of the tunnel excavation face with a time step controlled within 0.01 seconds to calculate the stress distribution and deformation response of the surrounding rock. All calculation results are output in spatial coordinates, with data accuracy controlled at the level of 0.001 MPa and 0.01 mm to ensure the accuracy of subsequent crack detection.

[0034] Step S12: Collect joint and fracture images based on tunnel excavation data; In this embodiment, based on the tunnel excavation simulation data obtained in step S11, a high-resolution industrial camera array is set up on the tunnel excavation face, containing at least six 12-megapixel cameras evenly distributed to cover the excavation face, with shooting angles ranging from 0° to 90°, and an image resolution of 4000×3000 pixels. Uniform lighting is ensured during shooting, and a ring-shaped LED fill light is used to ensure no shadows or reflections. After acquisition, the images are transmitted to an image processing workstation for preprocessing operations, including noise reduction, color correction, and distortion correction. Specifically, a median filter is used to remove salt-and-pepper noise, with the filter window set to 3×3 pixels. Color correction is performed based on a standard color chart for linear RGB channel correction. The preprocessed images serve as the basic data for joint and fracture identification, used for subsequent structural identification processing.

[0035] Step S13: Identify the layered joint structure based on the joint image; In this embodiment, image segmentation technology is applied to identify layered joint fracture structures based on preprocessed joint fracture images. First, the color image is converted to a grayscale image, and the Otsu thresholding method is used to automatically calculate the binarization threshold, typically between 120 and 150 grayscale values. After binarization, morphological closing operations are applied to fill small pores within the fractures, using a 3×3 square kernel as the structural element, with two iterations. Subsequently, a directional filter is used to extract linear features from the image, specifically a Gabor filter with a directional angle range of 0° to 180° and a step angle of 15°, to detect the orientation of the layered joint fractures. The extracted linear features are further enhanced using Hough transform to detect straight segments of the main layered joints, with a length threshold set to at least 50 pixels. Segments with a fracture spacing of less than 10 pixels are merged into continuous layered joints. Finally, detailed parameters such as the orientation, spacing, and distribution location of the layered joint fractures are output.

[0036] Step S14: Identify columnar joint fracture structures based on joint fracture images; In this embodiment, columnar joint fracture structures are identified in joint fracture images using texture analysis. Texture feature parameters, such as contrast, correlation, and energy, are obtained by calculating the gray-level co-occurrence matrix (GLCM) of the image. The texture window size is set to 15×15 pixels, with a sliding step of 5 pixels, covering the entire image area. Columnar joints often appear as mutually perpendicular or intersecting linear textures. Fourier transform based on texture directionality is used to analyze the main texture direction and identify typical intersection angles of columnar joints, with an angle threshold set between 70° and 110°. Combining texture features and geometric morphology, connected component analysis is used to segment the columnar joint regions, with a minimum region area threshold of 500 pixels to filter noise. The spatial distribution, orientation, and density data of columnar joint structures are systematically recorded, forming a detailed dataset of columnar joint fracture structures.

[0037] Step S15: Integrate the layered joint fracture structure and the columnar joint fracture structure to obtain joint fracture data.

[0038] In this embodiment, spatial data fusion is performed on the identified layered joint fracture structures and columnar joint fracture structures. A three-dimensional coordinate registration method is used to map the two-dimensional identification data of the two types of structures into a three-dimensional tunnel excavation face model. First, rigid body transformation is used to unify the coordinate systems of different structural datasets. The transformation parameters include translation vectors (Tx, Ty, Tz) and rotation matrices (rotation angles around the X, Y, and Z axes). The translation error is controlled within 1 mm, and the rotation error is controlled within 0.1 degrees. During the fusion process, the spatial proximity principle is used, setting the neighborhood radius to 0.5 meters, and merging the overlapping areas of layered and columnar joint data. A weighted average method is used to determine the attribute parameters of the overlapping parts. The final generated joint fracture data contains complete spatial distribution information, strike, dip angle, and fracture width (measurement accuracy 0.1 mm), and is labeled as layered and columnar joints by category. The data format conforms to the Geological Information System (GIS) standard, facilitating subsequent analysis module calls and processing.

[0039] Preferably, step S13 specifically includes: Step S131: Perform image preprocessing based on the joint and fracture image to obtain the joint and fracture image to be processed; In this embodiment, after importing the original joint and fracture image into the image processing system, image preprocessing is performed first. Preprocessing includes image denoising, enhancement, and geometric correction. Denoising uses a median filter with a 3×3 pixel filter window to effectively remove salt-and-pepper noise while preserving edge information. Image enhancement uses histogram equalization to adjust the image grayscale distribution, increasing contrast by at least 20% to facilitate clear extraction of fracture edges. Geometric correction uses camera intrinsic parameters and distortion coefficients for distortion correction. The distortion coefficients are calculated using standard images acquired by a calibration board; typical radial distortion parameters k1 are -0.2, k2 is 0.03, and tangential distortion parameters p1 and p2 are both less than 0.01. The corrected image resolution remains at 4000×3000 pixels. After the above preprocessing, the resulting joint and fracture image has low noise, clear edges, and no distortion, meeting the requirements for fracture edge recognition.

[0040] Step S132: Identify the fracture edges based on the joint fracture image to be processed, and obtain fracture edge data; In this embodiment, crack edge recognition is performed on the joint crack image to be processed. The Canny edge detection algorithm is used, with specific parameter settings of a low threshold of 50, a high threshold of 150, and an operator size of 3×3 pixels to ensure the accuracy and detail integrity of edge detection. Before edge detection, a Gaussian blur filter with a kernel size of 3×3 and a standard deviation of 1.4 is applied to the image to reduce the interference of high-frequency noise on edge recognition. The algorithm flow includes calculating the image gradient magnitude and direction, non-maximum suppression, double threshold detection, and edge connection. The edge output data is a binary image, where a pixel value of 1 represents a detected crack edge, and 0 represents the background. The edge data undergoes morphological thinning processing using a 3×3 structuring element and one iteration to ensure that the crack edge line width is uniformly 1 pixel.

[0041] Step S133: Perform linear detection based on the crack edge data to obtain linear crack data; In this embodiment, linear crack detection is performed based on crack edge data. A probabilistic Hough transform-based method is used to extract straight line features from the image. Specifically, the edge image is used as input, with a Hough spatial resolution of 1 pixel and an angular resolution of 1 degree. The detection threshold is set to 50 pixels to ensure that only line segments longer than 50 pixels are detected, excluding noisy short line segments. The detected lines are stored as parameters, including polar coordinate parameters ρ (distance from the origin, in pixels) and θ (angle, in degrees). To eliminate duplicate line segments, the polar coordinate parameters are clustered, with a clustering distance threshold of 10 pixels and an angular difference threshold of 5 degrees, merging line segments with similar angles and distances. Finally, linear crack data is obtained, clearly recording the start and end coordinates and direction angle of each crack.

[0042] Step S134: Calculate the fracture spacing based on the linear fracture data; In this embodiment, the fracture spacing is calculated based on linear fracture data. First, all linear fractures are classified according to their orientation angle, with a classification threshold of ±5 degrees considered as the same group. For each group of linear fractures, they are sorted based on the ρ value in their polar coordinate parameters, and the distance difference between adjacent fractures is calculated to obtain a spacing sequence. The spacing unit is converted to actual length according to the image resolution conversion ratio (e.g., 1 pixel corresponds to 0.5 mm). The mean and standard deviation of the spacing sequence are calculated, and outliers (exceeding ±2 times the standard deviation of the mean) are removed. Finally, the average spacing of each group of fractures is output, with the spacing accuracy controlled within 0.1 mm.

[0043] Step S135: Calculate the parallelism of the fractures based on the linear fracture data; In this embodiment, the parallelism of fractures is calculated based on linear fracture data. Using the orientation angle of linear fractures within the same group as a basis, the deviation angle between all linear fractures in the group and the average orientation of that group is calculated. The standard deviation is calculated as a quantitative indicator of parallelism, and parallelism is defined as 1 minus the standard deviation and the maximum permissible deviation angle (set to 15 degrees). During the calculation, the deviation angle is in degrees, and the parallelism result is limited to between 0 and 1. Parallelism ≥ 0.85 is considered highly parallel, and parallelism < 0.5 is considered low parallelism. The parallelism of each group is used as an important parameter for subsequent judgment of layered joint fracture structure.

[0044] Step S136: Determine the layered joint structure based on the fracture spacing and the degree of fracture parallelism.

[0045] In this embodiment, layered joint fracture structures are determined based on fracture spacing and fracture parallelism. Combining fracture spacing data and parallelism indicators, a judgment threshold is set: linear fracture groups with an average fracture spacing of less than 500 mm and a parallelism greater than 0.85 are judged as layered joint fracture structures. Combined with spatial distribution information from images, the layered joints are confirmed to have obvious layered arrangement characteristics. For fracture groups meeting the judgment criteria, the strike, dip angle (calculated by image-to-spatial coordinate transformation, dip angle range 0°~90°, error not exceeding 2°), and spacing distribution characteristics of the layered joints are recorded to construct a layered joint fracture structure dataset. The data format conforms to the Geological Information System (GIS) standard for easy subsequent analysis and retrieval.

[0046] Preferably, step S14 specifically includes: Step S141: Identify the vertical crack edge based on the joint and crack image to obtain the vertical crack contour data; In this embodiment, an input joint fracture image is used to first identify the vertical fracture edges. A gradient-based edge detection algorithm, specifically the Sobel operator, is employed. Gradient values ​​are calculated in both the horizontal and vertical directions, with the operator size being 3×3 pixels. After calculating the gradient magnitude and direction for each pixel, edge points consistent with the vertical direction (i.e., close to 90° or 270°, with a tolerance of ±10°) are extracted. To improve detection accuracy, a threshold of 50 (grayscale level) is set for the gradient magnitude; edge points below this value are discarded to prevent weak edges from affecting identification. Connectivity analysis is performed on the extracted vertical edge points, using 8-neighborhood connectivity to filter out isolated point groups smaller than 100 pixels. The final vertical fracture edge data is saved as a binary image with an edge line width of 1 pixel for subsequent contour analysis.

[0047] Step S142: Identify polygonal cross sections based on vertical crack profile data to obtain polygonal crack profile data; In this embodiment, polygonal cross-section recognition is performed based on the vertical crack edge data obtained in step S141. A contour extraction algorithm (such as the `findContours` function in OpenCV) is used to detect closed edge contours. The contour point set is approximated as a polygon using the Douglas-Peucker algorithm, with a tolerance parameter set to 5 pixels to reduce the number of points while preserving contour shape features. The identified polygonal contours must have at least 4 vertices and side lengths greater than 20 pixels, excluding contours that are too small or too fragmented. The polygonal crack contour data is stored as a vertex coordinate sequence, with the coordinate system consistent with the input image, facilitating spatial relationship calculations.

[0048] Step S143: Calculate the profile closure degree based on the polygonal crack profile data; In this embodiment, the closure degree of the polygonal crack profile is calculated based on the profile data. The closure degree is defined as the ratio of the actual perimeter of the profile to the perimeter of the smallest bounding rectangle of the polygon. Specifically, the perimeter is first calculated by summing the lengths of all sides of the polygon, and then the perimeter of the rectangle is obtained by summing the lengths of the four sides of the smallest bounding rectangle. The closure degree value ranges from 0 to 1; the closer the value is to 1, the closer the profile is to a regular closed structure. For each polygonal profile, its closure degree is calculated, and a closure degree threshold is set to 0.8. Profiles below this threshold are excluded and not included in subsequent planar connectivity analysis.

[0049] Step S144: Perform planar connectivity analysis based on the contour closure to obtain planar connectivity fracture data; In this embodiment, planar connectivity analysis is performed on polygonal fracture contour data with acceptable contour closure. First, a spatial adjacency matrix of the polygonal contours is constructed. A boundary intersection determination method is used: if two polygonal contours have intersecting boundaries or touching vertices, the corresponding element in the adjacency matrix is ​​set to 1; otherwise, it is 0. Based on the adjacency matrix, a connected subgraph detection algorithm in graph theory is executed to extract all planar connected fracture groups. Each connected group is a set of spatially connected polygonal fracture contours. This data is stored in the form of connected group numbers and their contained contour vertex data for easy subsequent spatial statistical analysis.

[0050] Step S145: Calculate the spatial density based on the planar connected fracture data; In this embodiment, spatial density is calculated for surface-level connected fracture data. The tunnel surrounding rock area is divided into a 3D grid using a voxel mesh method, with each grid cell having a side length of 1 meter. The number of fracture contours contained within each grid cell is counted to obtain the spatial density distribution. Spatial density is defined as the ratio of the total fracture length to the volume per unit volume, expressed in meters per cubic meter. When calculating the spatial density of each connected fracture group, the perimeters of all contours in the group are summed and divided by its envelope volume (the convex hull algorithm is used to calculate the 3D envelope volume, with an accuracy error of less than 5%). Connected fracture groups with a spatial density higher than 1.5 meters per cubic meter are recorded as high-density regions.

[0051] Step S146: Determine and record the columnar joint fracture structure based on the spatial density.

[0052] In this embodiment, columnar joint fracture structures are determined and recorded based on the spatial density calculated in step S145. A columnar joint fracture structure is defined as a region containing planar connected fractures with a spatial density greater than 1.5 m / m³. Connected fracture groups meeting the criteria are marked as columnar joint fracture structures, and the spatial extent of the structure (envelope volume coordinate range), the number of polygonal outlines within the fracture group, and the average closure degree are recorded. The structural information is stored in a database, with fields including structure ID, spatial boundary box, density value, and fracture characteristic statistics. This data is used for subsequent stability analysis and risk assessment.

[0053] Preferably, step S2 specifically includes: Step S21: Perform surrounding rock evolution analysis based on tunnel excavation data to obtain surrounding rock evolution data; In this embodiment, tunnel excavation data is used to analyze the evolution of the surrounding rock. First, based on the three-dimensional surrounding rock structure data and physical and mechanical parameters collected during tunnel excavation, a time-series simulation of the deformation and failure process of the surrounding rock is conducted using numerical calculation methods. Input parameters include the initial stress state of the surrounding rock, an elastic modulus typically ranging from 30 to 50 GPa, a Poisson's ratio ranging from 0.25 to 0.30, the specific geometric dimensions of the excavation, and the excavation advance rate. By analyzing the stress, strain, and fracture conditions of the surrounding rock at different time points, the deformation process of the surrounding rock is meticulously recorded. A three-dimensional unit with a side length of 0.5 meters is used for mesh generation, and the time step is set to 0.1 days to ensure accurate capture of the spatiotemporal changes in the surrounding rock evolution. All simulation results are saved in a three-dimensional data format containing displacement and stress changes for further processing.

[0054] Step S22: Calculate the principal stress variation rate based on the surrounding rock evolution data; In this embodiment, based on the surrounding rock evolution data obtained in step S21, the principal stress values ​​of each three-dimensional grid cell are extracted, including the maximum, intermediate, and minimum principal stresses. The principal stress values ​​at adjacent time points are differentially calculated to determine the change in principal stress per unit time, with a time interval of 0.1 days. The rate of change of each principal stress component is calculated separately, and the maximum rate of change of principal stress is used as the representative value for that grid cell. The unit of change rate is megapascals per day. All calculation results are summarized into a three-dimensional distribution map of the principal stress change rate to accurately locate the area of ​​fastest stress change in the surrounding rock.

[0055] Step S23: Identify stress concentration areas in the surrounding rock based on the principal stress variation rate; In this embodiment, the principal stress change rate data calculated in step S22 is used as the basis. A threshold is set to filter out grid cells with a change rate greater than 0.5 MPa per day, and these grid points are considered potential stress concentration areas. A spatial clustering algorithm is used to aggregate adjacent grid points with high change rates. The neighborhood search radius is set to 1 meter, and the minimum number of points in the cluster is 10. Spatial clusters smaller than 5 cubic meters are filtered out to avoid noise interference. The finally determined stress concentration areas include their spatial extent, volume, and average principal stress change rate, providing a location basis for subsequent analysis.

[0056] Step S24: Calculate the stress disturbance intensity based on the stress concentration area of ​​the surrounding rock to obtain stress disturbance data; In this embodiment, for the identified stress concentration areas, a weighted average is performed based on the principal stress change rate of each grid cell and the grid volume to calculate the stress disturbance intensity of that area. Specifically, the principal stress change rate of each grid cell is multiplied by its corresponding volume, the products of all grid cells are summed, and then divided by the total volume of the entire area. The result is taken as the stress disturbance intensity of that area, expressed in megapascals per day. This value is used to quantitatively describe the stress disturbance level of the area.

[0057] Step S25: Identify the unstable areas of the tunnel based on stress disturbance data; In this embodiment, the stress disturbance intensity calculated in step S24 is used to determine the areas surrounding the tunnel, and areas with a disturbance intensity greater than 0.6 MPa per day are defined as unstable areas. Through three-dimensional connectivity analysis, the spatial continuity of these areas is confirmed, excluding isolated volumes smaller than 2 cubic meters to ensure the spatial integrity of the unstable areas. The spatial boundaries of the unstable areas are recorded in the form of maximum and minimum coordinates in a three-dimensional coordinate system for easy location and monitoring.

[0058] Step S26: Detect the surrounding rock convergence displacement in the unstable area of ​​the tunnel and generate surrounding rock convergence displacement data; In this embodiment, for the unstable area identified in step S25, convergence displacement detection of the surrounding rock is carried out. A combination of ground-based laser scanning technology and multi-temporal oblique photogrammetry is used to periodically collect surface point cloud data of the unstable area, with the point spacing controlled within 5 mm, maintaining measurement accuracy at the millimeter level. The multi-temporal point cloud data is registered to achieve high-precision displacement calculation, obtaining the three-dimensional displacement vector and time series of each monitoring point. A convergence displacement threshold of 3 mm is set; displacements below this threshold are considered measurement noise. The detection results are saved as the displacement magnitude and direction of the three-dimensional spatial points for subsequent dynamic monitoring and analysis.

[0059] Step S27: Predict the probability of rockfall based on the surrounding rock convergence displacement data.

[0060] In this embodiment, the surrounding rock convergence displacement data obtained in step S26 is used to predict the probability of surrounding rock detachment. The maximum convergence displacement value, convergence velocity per unit time, and spatial gradient of the convergence displacement are selected as input parameters. The maximum convergence displacement value is directly extracted from the maximum displacement of all monitoring points within the instability area. The convergence velocity is calculated using the displacement difference between adjacent time points, in millimeters per hour. The spatial gradient is obtained as the ratio of the displacement difference between neighboring monitoring points to their distance, in millimeters per meter. Statistical regression analysis is performed on these parameters using historical data to determine the weight coefficients of each parameter, establishing a multiple regression equation. The predicted result is a probability value between 0 and 1, with values ​​exceeding 1 being assigned as 1 and values ​​below 0 as 0. The prediction results are summarized and recorded according to spatial regions for risk assessment purposes.

[0061] Preferably, step S25 specifically includes: Step S251: Calculate the peak shear stress based on the stress disturbance data; In this embodiment, shear stress in each grid cell within the tunnel surrounding rock area is extracted and statistically analyzed based on stress disturbance data. The shear stress value is calculated from the three-dimensional stress tensor data collected in real-time by the surrounding rock mechanical monitoring system; specifically, the principal stress components are transformed into shear stress components using stress tensor transformation. During the statistical process, the shear stress values ​​are discretized according to a time series, and the maximum shear stress value within a certain time window is selected as the peak shear stress. The time window is set to 24 hours, and the statistical accuracy is 0.1 MPa. The peak shear stress is expressed in MPa, and the statistical results are stored in a spatial coordinate system to ensure that each grid cell corresponds to a unique peak shear stress.

[0062] Step S252: Obtain the yield strength of the surrounding rock; identify potential unstable elements based on the yield strength and peak shear stress of the surrounding rock to obtain potential unstable element data; In this embodiment, the yield strength data of the surrounding rock is obtained based on field rock test data. Yield strength is typically determined through indoor triaxial compression tests of rock. Test parameters include the peak compressive strength and friction coefficient of the rock. Typical peak compressive strength ranges from 50 to 150 MPa, and the friction coefficient ranges from 0.6 to 0.85. Using the Mohr-Coulomb yield criterion combined with peak shear stress, each grid cell is identified as a potential instability cell. Specifically, when the peak shear stress exceeds the critical value of the corresponding yield strength of the surrounding rock, the grid cell is marked as a potential instability cell. The data for potential instability cells is stored in three-dimensional coordinates and instability levels (divided into three levels according to the degree of shear stress exceeding the yield strength, corresponding to shear stress exceeding the yield strength by 10%, 20%, and 30%, respectively).

[0063] Step S253: Perform laser cross-section scanning based on the potential unstable element data to obtain laser cross-section data; In this embodiment, a 3D laser cross-section scanner is used to scan the area where potential unstable elements are located. The scanner uses a laser wavelength of 1550 nanometers, with a measurement accuracy of ±0.5 millimeters. The scanning range covers the surrounding rock cross-section where the potential unstable elements are located, and the cross-section size is not less than 5 meters × 5 meters. During the scanning process, the spacing between scanning points is controlled at 1 millimeter to ensure complete acquisition of cross-sectional details. The point cloud data obtained from the scan is registered and filtered to remove environmental noise and data redundancy, generating high-density laser cross-section data. The data format adopts a standard XYZ point cloud file for easy subsequent calculations.

[0064] Step S254: Calculate the deformation convergence based on the laser cross-section data; In this embodiment, laser cross-sectional data is used to calculate the deformation convergence of the cross-section. Specifically, based on point cloud data from multiple time points before and after the cross-section, rigid body registration technology is employed to align the cross-section positions, with the registration error controlled within 0.2 mm. Convergence is calculated by measuring the displacement vectors of adjacent points on the cross-section to obtain the average displacement value of the overall convergence direction of the cross-section. Deformation convergence is defined as the number of millimeters the cross-section converges per unit time, with the time interval fixed at 12 hours. By statistically analyzing the average and maximum convergence values ​​at each point on the cross-section, a detailed spatial distribution map of the cross-sectional deformation convergence is obtained.

[0065] Step S255: Identify high-deformation section areas based on deformation convergence; perform rock mass microseismic monitoring based on the high-deformation section areas to obtain rock mass microseismic data; In this embodiment, based on the deformation convergence data obtained in step S254, a deformation convergence threshold of 2 mm every 12 hours is set, and cross-sectional areas exceeding this threshold are defined as high-deformation cross-sectional areas. Microseismic monitoring points are located within the spatial range of these high-deformation cross-sectional areas, and a multi-channel seismic sensor array is deployed. The sensor sampling frequency is set to 1000 Hz, and the sensitivity is -110 dB. Through continuous real-time monitoring, rock mass microseismic data, including source location, magnitude, and time, is collected. The microseismic data uses an event-triggered storage method, with a trigger threshold of microseismic magnitude ≥ -3, and the data format is the standard SEED format.

[0066] Step S256: Perform vibration frequency statistics based on the rock mass microseismic data to obtain high-frequency rock mass microseismic data; In this embodiment, vibration frequency statistical analysis is performed on the rock mass microseismic data. First, time-frequency analysis is conducted on the acquired microseismic signals, using short-time Fourier transform to decompose the signal frequency components, with a frequency resolution set to 1 Hz. Statistical analysis is performed in one-hour time windows, counting the number of microseismic events in each frequency band. High-frequency rock mass microseismic events are defined as events with frequencies higher than 500 Hz. When the number of high-frequency events exceeds 10 within the statistical window, the microseismic data for that time period is marked as high-frequency rock mass microseismic data. The statistical results are saved as a frequency-time two-dimensional distribution matrix for easy dynamic monitoring.

[0067] Step S257: Determine the tunnel instability zone based on high-frequency rock mass microseismic data.

[0068] In this embodiment, based on the high-frequency rock mass microseismic data obtained in step S256 and combined with sensor spatial layout information, a three-dimensional positioning algorithm is used to determine the spatial distribution of high-frequency microseismic events. Spatial density clustering is used, with a neighborhood radius of 1.5 meters, to identify spatial clusters of highly concentrated microseismic events. The tunnel instability area is determined to be a region where the spatial density of microseismic events exceeds 50 events per cubic meter. The spatial boundary of the instability area is described using a three-dimensional polyhedral model, and the coordinate range and volume information are recorded and stored in a GIS-compatible format for easy visualization analysis and risk management.

[0069] Preferably, step S26 specifically includes: Step S261: Perform crown region detection on the unstable area of ​​the tunnel to obtain crown region data; In this embodiment, the location of the crown region is determined using a combined laser scanning and ground-penetrating radar (GPR) detection technique, based on the three-dimensional spatial coordinate data of the tunnel instability area. The laser scanning equipment uses LiDAR, with a scanning accuracy controlled within ±1 mm, and the scanning range covers the top of the tunnel instability area and an adjacent 5-meter radius. The GPR uses a 100MHz antenna with a penetration depth of 5 meters to detect changes in the internal structure of the surrounding rock. The laser scanning point cloud and radar reflection data are fused, and the crown region is extracted based on terrain features and rock mass geometry. Specifically, the extraction method involves calculating the geometric curvature of the crown region, setting a curvature threshold of 0.05 / mm, and defining the set of points with curvature greater than this threshold as the crown region. Finally, point cloud data of the crown region is generated in XYZ coordinate set format for easy subsequent processing.

[0070] Step S262: Identify the initial profile of the cross-section based on the data of the arch area; In this embodiment, the initial profile of the cross-section is identified based on the point cloud data of the arch area obtained in step S261. A cross-section cutting technique is used to cut the arch point cloud data into multiple cross-sections at 1-cm intervals along the longitudinal direction of the tunnel. Within each cross-section, edge detection algorithms are used to identify contour edge points. The edge detection employs the gradient-based Canny operator, with thresholds set to a low threshold of 30 and a high threshold of 90 to ensure complete edge point capture. The identified edge point set is fitted with polygons, and outliers are removed using the RANSAC algorithm to obtain a closed polygon of the initial profile of the cross-section. The contour polygon data includes the three-dimensional coordinates of each vertex, with the vertex spacing controlled within 5 mm to ensure the accuracy of the cross-section profile.

[0071] Step S263: Continuously monitor the prism target based on the data of the dome area to obtain prism target monitoring data; In this embodiment, optical prism targets, each 50 mm × 50 mm in size, are installed in the arch area as deformation monitoring targets, with an installation density of 4 targets per square meter. A total station is used to continuously and automatically measure the prism targets at a frequency of once per hour, with a measurement accuracy of ±0.1 mm. The total station utilizes laser ranging and angle measurement technology, combined with trigonometric leveling, to determine the spatial coordinates of each prism target. The monitoring data includes a time series of the target point's three-dimensional coordinates, stored in a CSV file format, with timestamps accurate to the second. During monitoring, the data is calibrated in real time to remove outliers, ensuring the reliability and consistency of the monitoring data.

[0072] Step S264: Calculate the radial displacement of the cross section based on the prism target monitoring data; calculate the convergence rate of the cross section based on the prism target monitoring data; perform crack propagation analysis based on the radial displacement and convergence rate of the cross section to obtain the crack propagation coordinates; In this embodiment, based on the prism target monitoring data obtained in step S263, the radial displacement and convergence rate of the cross-section are calculated. Radial displacement is obtained by comparing the changes in the radial coordinates of the corresponding prism targets at continuous measurement time points; the radial direction is defined as the perpendicular direction from the tunnel centerline to the target point. The convergence rate is the change in radial displacement per unit time, with a fixed time interval of 1 hour. The radial displacement and convergence rate of all target points on the cross-section are statistically analyzed, and a weighted average method is used to calculate the overall cross-section convergence rate, with the weight being the reciprocal of the distance from each target point to the tunnel center. Based on the radial displacement and convergence rate data, the fracture propagation criterion in fracture mechanics theory is used to determine the fracture propagation coordinates. The specific calculation steps are as follows: when the prism target displacement in a certain area is greater than 2 mm and the convergence rate exceeds 0.5 mm / hour, the location is determined to be a fracture propagation zone. The coordinates of the target points that meet the conditions are summarized to form a set of fracture propagation coordinates.

[0073] Step S265: Map the fracture propagation coordinates to the initial profile of the cross section to determine the degree of surrounding rock convergence displacement and obtain surrounding rock convergence displacement data.

[0074] In this embodiment, the fracture propagation coordinates obtained in step S264 are mapped to the initial cross-sectional contour determined in step S262. The mapping process first performs a two-dimensional planar projection on the cross-sectional contour polygon to establish a local coordinate system, and then transforms the three-dimensional fracture propagation points into this plane. Based on the distance between the spatial point and the edge of the cross-sectional contour, it is determined whether the fracture propagation point falls within the internal region of the cross-section; internal points are then included in the surrounding rock convergence displacement range. Subsequently, the degree of surrounding rock convergence displacement is calculated, using the average and maximum radial displacements of the fracture propagation points as convergence displacement indicators, with units in millimeters. The convergence displacement data structure includes the cross-sectional number, the number of fracture propagation points, the average convergence displacement, and the maximum convergence displacement, and the storage format adopts a database table structure for easy subsequent querying and analysis.

[0075] Preferably, step S3 specifically includes: Step S31: Calculate the joint density index based on the joint and fracture data; calculate the joint intersection angle based on the joint and fracture data; In this embodiment, based on the collected images and 3D scan data of joints and fissures in the tunnel surrounding rock, the spatial distribution information of the joints and fissures is first extracted. The joint density index is calculated using the total length of joint surfaces per unit volume. Specifically, within the sampling volume of the tunnel surrounding rock (standard volume set at 1 cubic meter), all joint and fissure surfaces are identified and their lengths are calculated using laser scanning point cloud data combined with a fault plane detection algorithm. The sum of all joint lengths is divided by the sampling volume to obtain the joint density index, expressed in meters per cubic meter. The calculation of joint intersection angles is based on the angle between the normal vectors of the joint surfaces. The normal vectors of each joint surface are first extracted, and the angle between any two joint surfaces is calculated using the vector dot product formula. After statistically analyzing all angle data, an angle interval is divided using a frequency distribution histogram, ranging from 0° to 90°, into nine 10° intervals. By calculating the joint intersection frequency in each interval, the distribution characteristics of the joint intersection angles are obtained. In the above calculation process, data preprocessing includes removing joint surfaces with a length of less than 0.1 meters to eliminate interference from minor fractures.

[0076] Step S32: Determine the highly intersecting cutting area based on the joint density index and joint intersection angle; In this embodiment, the highly intersecting cutting regions in the tunnel surrounding rock are determined using the joint density index and joint intersection angle distribution calculated in step S31. The specific criteria are: a joint density index exceeding 15 m / m³ and a joint intersection angle with a frequency exceeding 40% in the 15° to 45° range. First, the surrounding rock space is divided into 1-meter square grid cells. For each cell, the joint density index and intersection angle frequency are calculated, and grid cells meeting the above numerical conditions are selected. This process is implemented using spatial data processing software, with data in a three-dimensional grid format. Grid numbers, coordinates, and corresponding values ​​are all recorded. All grid cells meeting the conditions constitute a spatial set of highly intersecting cutting regions. For regions with ambiguous boundaries, a dilation operator is used to expand the boundary by 3 cm to ensure regional continuity.

[0077] Step S33: Map the probability of surrounding rock detachment to the highly interlaced cutting area to obtain the area with high probability of surrounding rock detachment; In this embodiment, the rockfall probability data obtained in previous steps is mapped to the highly staggered cutting regions identified in step S32. The rockfall probability data comes from the rock stability analysis database, and the data format is a three-dimensional coordinate point cloud with a point spacing of approximately 5 cm. The probability value is labeled at each point. The mapping process first spatially overlays the three-dimensional grid cells of the highly staggered cutting regions with the rockfall probability point cloud, and uses spatial indexing technology to quickly match grid cells with probability points. For all probability points contained in a grid, the average probability is calculated as the rockfall probability of that grid cell. A rockfall probability threshold of 0.7 is set, and grid cells exceeding this threshold are classified as high rockfall probability regions. Finally, spatial distribution data of three-dimensional high rockfall probability regions is formed and stored as a structured data table, containing grid coordinates, grid number, and rockfall probability.

[0078] Step S34: Identify the characteristics of broken rock blocks based on areas with high probability of rock fragmentation, and obtain the characteristic data of broken rock blocks; In this embodiment, the characteristics of broken rock blocks within the high-probability rockfall area determined in step S33 are further identified. The broken rock block feature data is acquired based on point cloud data collected by a high-precision 3D laser scanner on-site, with scanning accuracy controlled within ±2 mm. A point cloud segmentation algorithm is used to segment the point cloud within the high-probability rockfall area according to the rock block boundaries. Boundary identification is based on point cloud density and surface curvature characteristics, with a density threshold of no less than 5000 points per cubic meter and a curvature threshold of 0.03 / mm. The volume, maximum size, and surface irregularity of each segmented rock block are measured. Rock blocks with a volume less than 0.1 cubic meters or a maximum size less than 0.2 meters are discarded, and the remaining data constitute the broken rock block feature dataset. Each rock block feature includes the rock block number, spatial coordinate range, volume, size, and surface roughness.

[0079] Step S35: Determine the degree of structural surface activation based on the characteristic data of the fractured rock blocks, and obtain structural surface activation data; In this embodiment, the activation degree of the structural surface is calculated based on the characteristic data of the broken rock blocks obtained in step S34. The activation degree of the structural surface is defined as a function of the number of broken rock blocks and their total volume per unit volume. The surrounding rock space is divided into 1 cubic meter grid cells, and the number and total volume of broken rock blocks in each cell are counted to calculate the activation index. The activation index is calculated as follows: Activation index = Total volume of broken rock blocks (cubic meters) × Number of broken rock blocks, and the activation index threshold is set to 5. Cells with an activation index exceeding this value are marked as highly activated areas of the structural surface. All grid cell activation index data are stored in a structured database, including fields such as grid number, activation index, number of rock blocks, and volume.

[0080] Step S36: Identify the shear slip trend of the structural surface based on the structural surface activation data to obtain the structural surface shear slip data.

[0081] In this embodiment, based on the structural surface activation data obtained in step S35, shear slip trend identification is performed on the surrounding rock structural surfaces. Geomechanical analysis methods are employed, combined with on-site stress measurement data and activation index data. Regions with a structural surface activation index greater than 5 are selected, and the corresponding structural surface normal direction and shear force direction data are extracted. The normal direction is determined by fitting a plane to the structural surface point cloud extracted by laser scanning, and the shear force direction is calculated from stress measurement point data. The angle between the structural surface normal and the direction of the maximum principal stress is calculated; structural surfaces with an angle less than 30° are considered shear surfaces. Combining displacement monitoring data on the shear surfaces, shear slip velocity is calculated through time series analysis, with a slip velocity threshold set to 0.1 mm / hour. Finally, structural surfaces that meet the shear surface condition and whose slip velocity exceeds the threshold are marked as having a shear slip trend. Shear slip data is stored with structural surface number, spatial location, slip velocity, and direction to support subsequent analysis and early warning.

[0082] Preferably, step S4 specifically includes: Step S41: Calculate the gravity direction vector based on the structural shear slip data; calculate the slip direction vector based on the structural shear slip data; calculate the slip surface inclination angle based on the gravity direction vector and the slip direction vector; identify potential unstable blocks based on the slip surface inclination angle to obtain potential unstable block data; In this embodiment, for the shear slip data of the structural surfaces, the spatial location information and slip vector data of each unstable structural surface are first extracted. The gravity direction vector is fixed as a vertically downward unit vector, with its direction perpendicular to the geographic coordinate system of the tunnel excavation area. The slip direction vector is calculated from the slip direction and magnitude data collected by the three-dimensional displacement sensors installed on site. The direction angle (azimuth angle) and inclination angle of the structural surface slip vector are calculated and converted into a three-dimensional spatial vector representation. The inclination angle of the slip surface is calculated by the angle between the structural surface normal and the gravity direction vector, using the vector dot product formula to calculate the angle, with the inclination angle threshold set within the range of 20° to 70°. Structural surfaces with inclination angles exceeding 20° but less than 70° are identified as slip surfaces. All slip surfaces are screened to identify the spatial range of potential unstable blocks. Potential unstable blocks are formed by the intersection of multiple adjacent slip surfaces. The spatial range is generated using three-dimensional geometric construction technology and stored as a three-dimensional polyhedral data structure, containing the normals of each surface, boundary points, and vertex coordinate information.

[0083] Step S42: Identify geometric closure based on potential unstable block data; determine wedge-shaped blocks based on geometric closure to obtain wedge-shaped block data; In this embodiment, the geometric closure of the potentially unstable block data obtained in step S41 is calculated. Geometric closure is determined based on the spatial overlap and boundary conformity between all structural surfaces of the potentially unstable block. Specifically, the ratio of the length of the intersection line of each structural surface to the theoretical intersection line length is calculated. Boundary conformity is defined as a distance error between boundary points not exceeding 5 cm. By calculating the spatial error of the boundary points of each surface of the block, blocks with an error less than 5 cm and no obvious cracks or faults on the structural surfaces are considered to have good geometric closure. The portion of the block with a geometric closure greater than 90% is identified as a wedge-shaped block. The wedge-shaped block data includes the block number, the number of structural surfaces, the block volume (cubic meters), the orientation of each structural surface, and the length of the intersection line. The data is stored in a structured table format for easy subsequent querying and analysis.

[0084] Step S43: Identify the degree of voiding of the structural face based on the wedge block data to obtain voiding data of the structural face; assess the instability risk level based on the voiding data of the structural face; In this embodiment, the degree of free space on the structural face is assessed based on the wedge-shaped block data identified in step S42. The degree of free space on the structural face is defined as the ratio of the actual gap size between the block and the surrounding rock to the block's surface area. A 3D laser scanner is used to acquire boundary point cloud data between the surrounding rock and the wedge-shaped block, with an accuracy controlled within ±2 mm. Point cloud registration technology is used to align the block's point cloud with the surrounding rock's point cloud, calculate the distance between the block's surface and the surrounding rock's surface, and statistically analyze the average gap width, setting a gap threshold of 5 mm. When the distance between the block and the surrounding rock exceeds this threshold, a free space is considered to exist. The volume of the free space and its proportion in the block's surface area are recorded as structural face free space data. Structural face free space data is stored in the form of block number, free space volume (cubic meters), and free space ratio (percentage). Based on the free space ratio classification, blocks with a free space ratio exceeding 10% are rated as high free space level, 5%-10% as medium free space level, and below 5% as low free space level. Free space level data is stored in the structural database to support subsequent analysis.

[0085] Step S44: Perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

[0086] In this embodiment, the stability analysis of the tunnel surrounding rock is performed in conjunction with the instability risk level determined in step S43. The instability risk level is based on comprehensive indicators such as the openness level of the structural face, the geometric closure of the block, the dip angle of the slip surface, and the volume of the block. A weighted scoring method is used to calculate the total risk score, with the weights set as follows: openness of the structural face 40%, geometric closure 30%, dip angle of the slip surface 20%, and volume of the block 10%. Each indicator is mapped to a score from 0 to 1 according to a preset standardized interval. A risk score exceeding 0.7 is considered high risk, 0.4-0.7 is medium risk, and below 0.4 is low risk. The surrounding rock stability analysis process is based on on-site sensor data and a three-dimensional geological model. The mechanical state of the surrounding rock is calculated using numerical simulation and spatial data analysis techniques, and the output includes the stability level, the location of potential unstable blocks, and the risk level. The final generated surrounding rock stability data includes spatial coordinates, risk level, and detailed values ​​of risk factors, and is stored in a structured database to support tunnel excavation safety management and early warning applications.

[0087] Preferably, this specification also provides an intelligent analysis system for the stability of surrounding rock in tunnel excavation, used to execute the intelligent analysis method for the stability of surrounding rock in tunnel excavation as described above. The intelligent analysis system for the stability of surrounding rock in tunnel excavation includes: The joint and fracture detection module is used to acquire geological data of the surrounding rock of the tunnel and to simulate tunnel excavation to obtain tunnel excavation data; based on the tunnel excavation data, joint and fracture detection is performed to obtain joint and fracture data. The surrounding rock spalling probability prediction module is used to perform stress disturbance analysis based on tunnel excavation data to obtain stress disturbance data; identify unstable areas of the tunnel based on the stress disturbance data; detect the surrounding rock convergence displacement of the unstable areas of the tunnel to generate surrounding rock convergence displacement data; and predict the probability of surrounding rock spalling based on the surrounding rock convergence displacement data. The structural surface shear slip trend identification module is used to detect structural surface activation based on the surrounding rock detachment probability and joint and fracture data to obtain structural surface activation data; and to identify structural surface shear slip trend based on structural surface activation data to obtain structural surface shear slip data. The surrounding rock stability analysis module is used to identify potentially unstable blocks based on structural surface shear slip data, obtain potential unstable block data; detect wedge-shaped blocks based on the potential unstable block data, obtain wedge-shaped block data; assess the instability risk level based on the wedge-shaped block data; and perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

[0088] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0089] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. An intelligent analysis method for the stability of surrounding rock in tunnel excavation, characterized in that, Includes the following steps: Step S1: Obtain geological data of the surrounding rock of the tunnel and perform tunnel excavation simulation to obtain tunnel excavation data; perform joint and fracture detection based on the tunnel excavation data to obtain joint and fracture data; Step S2: Perform stress disturbance analysis based on tunnel excavation data to obtain stress disturbance data; The unstable areas of the tunnel are identified based on stress disturbance data; the surrounding rock convergence displacement is detected in the unstable areas of the tunnel to generate surrounding rock convergence displacement data; and the probability of surrounding rock detachment is predicted based on the surrounding rock convergence displacement data. Step S3: Based on the surrounding rock spalling probability and joint and fracture data, perform structural surface activation detection to obtain structural surface activation data; based on the structural surface activation data, identify the structural surface shear slip trend to obtain structural surface shear slip data. Step S4: Identify potentially unstable blocks based on structural shear slip data to obtain data on potentially unstable blocks; Wedge block detection is performed based on potential unstable block data to obtain wedge block data; the instability risk level is assessed based on the wedge block data; Based on the instability risk level, the stability of the surrounding rock is analyzed, and the stability data of the surrounding rock is generated.

2. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 1, characterized in that, Step S1 is as follows: Step S11: Obtain geological data of the surrounding rock of the tunnel and conduct tunnel excavation simulation to obtain tunnel excavation data; Step S12: Collect joint and fracture images based on tunnel excavation data; Step S13: Identify the layered joint structure based on the joint image; Step S14: Identify columnar joint fracture structures based on joint fracture images; Step S15: Integrate the layered joint fracture structure and the columnar joint fracture structure to obtain joint fracture data.

3. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 2, characterized in that, Step S13 is as follows: Step S131: Perform image preprocessing based on the joint and fracture image to obtain the joint and fracture image to be processed; Step S132: Identify the fracture edges based on the joint fracture image to be processed, and obtain fracture edge data; Step S133: Perform linear detection based on the crack edge data to obtain linear crack data; Step S134: Calculate the fracture spacing based on the linear fracture data; Step S135: Calculate the parallelism of the fractures based on the linear fracture data; Step S136: Determine the layered joint structure based on the fracture spacing and the degree of fracture parallelism.

4. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 2, characterized in that, Step S14 is as follows: Step S141: Identify the vertical crack edge based on the joint and crack image to obtain the vertical crack contour data; Step S142: Identify polygonal cross sections based on vertical crack profile data to obtain polygonal crack profile data; Step S143: Calculate the profile closure degree based on the polygonal crack profile data; Step S144: Perform planar connectivity analysis based on the contour closure to obtain planar connectivity fracture data; Step S145: Calculate the spatial density based on the planar connected fracture data; Step S146: Determine and record the columnar joint fracture structure based on the spatial density.

5. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 1, characterized in that, Step S2 is as follows: Step S21: Perform surrounding rock evolution analysis based on tunnel excavation data to obtain surrounding rock evolution data; Step S22: Calculate the principal stress variation rate based on the surrounding rock evolution data; Step S23: Identify stress concentration areas in the surrounding rock based on the principal stress variation rate; Step S24: Calculate the stress disturbance intensity based on the stress concentration area of ​​the surrounding rock to obtain stress disturbance data; Step S25: Identify the unstable areas of the tunnel based on stress disturbance data; Step S26: Detect the surrounding rock convergence displacement in the unstable area of ​​the tunnel and generate surrounding rock convergence displacement data; Step S27: Predict the probability of rockfall based on the surrounding rock convergence displacement data.

6. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 5, characterized in that, Step S25 is as follows: Step S251: Calculate the peak shear stress based on the stress disturbance data; Step S252: Obtain the yield strength of the surrounding rock; identify potential unstable elements based on the yield strength and peak shear stress of the surrounding rock to obtain potential unstable element data; Step S253: Perform laser cross-section scanning based on the potential unstable element data to obtain laser cross-section data; Step S254: Calculate the deformation convergence based on the laser cross-section data; Step S255: Identify high-deformation section regions based on deformation convergence; Rock mass microseismic monitoring was conducted in the high-deformation section area to obtain rock mass microseismic data; Step S256: Perform vibration frequency statistics based on the rock mass microseismic data to obtain high-frequency rock mass microseismic data; Step S257: Determine the tunnel instability zone based on high-frequency rock mass microseismic data.

7. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 5, characterized in that, Step S26 is as follows: Step S261: Perform crown region detection on the unstable area of ​​the tunnel to obtain crown region data; Step S262: Identify the initial profile of the cross-section based on the data of the arch area; Step S263: Continuously monitor the prism target based on the data of the dome area to obtain prism target monitoring data; Step S264: Calculate the radial displacement of the cross section based on the prism target monitoring data; calculate the convergence rate of the cross section based on the prism target monitoring data; perform crack propagation analysis based on the radial displacement and convergence rate of the cross section to obtain the crack propagation coordinates; Step S265: Map the fracture propagation coordinates to the initial profile of the cross section to determine the degree of surrounding rock convergence displacement and obtain surrounding rock convergence displacement data.

8. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 1, characterized in that, Step S3 is as follows: Step S31: Calculate the joint density index based on the joint and fracture data; calculate the joint intersection angle based on the joint and fracture data; Step S32: Determine the highly intersecting cutting area based on the joint density index and joint intersection angle; Step S33: Map the probability of surrounding rock detachment to the highly interlaced cutting area to obtain the area with high probability of surrounding rock detachment; Step S34: Identify the characteristics of broken rock blocks based on areas with high probability of rock fragmentation, and obtain the characteristic data of broken rock blocks; Step S35: Determine the degree of structural surface activation based on the characteristic data of the fractured rock blocks, and obtain structural surface activation data; Step S36: Identify the shear slip trend of the structural surface based on the structural surface activation data to obtain the structural surface shear slip data.

9. The intelligent analysis method for the stability of surrounding rock in tunnel excavation according to claim 1, characterized in that, Step S4 is as follows: Step S41: Calculate the gravity direction vector based on the structural surface shear slip data; calculate the slip direction vector based on the structural surface shear slip data; calculate the slip surface inclination angle based on the gravity direction vector and the slip direction vector; Potentially unstable blocks are identified based on the dip angle of the slip surface, and data on potentially unstable blocks are obtained. Step S42: Identify geometric closure based on potential unstable block data; determine wedge-shaped blocks based on geometric closure to obtain wedge-shaped block data; Step S43: Identify the degree of void in the structural face based on the wedge block data to obtain the void data of the structural face; Assess the level of instability risk based on the structural face-space data; Step S44: Perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

10. An intelligent analysis system for the stability of surrounding rock in tunnel excavation, characterized in that, For executing the intelligent analysis method for tunnel surrounding rock stability during tunnel excavation as described in claim 1, the intelligent analysis system for tunnel surrounding rock stability during tunnel excavation includes: The joint and fracture detection module is used to acquire geological data of the surrounding rock of the tunnel and to simulate tunnel excavation to obtain tunnel excavation data; based on the tunnel excavation data, joint and fracture detection is performed to obtain joint and fracture data. The surrounding rock spalling probability prediction module is used to perform stress disturbance analysis based on tunnel excavation data to obtain stress disturbance data; identify unstable areas of the tunnel based on the stress disturbance data; detect the surrounding rock convergence displacement of the unstable areas of the tunnel to generate surrounding rock convergence displacement data; and predict the probability of surrounding rock spalling based on the surrounding rock convergence displacement data. The structural surface shear slip trend identification module is used to detect structural surface activation based on the surrounding rock detachment probability and joint and fracture data to obtain structural surface activation data; and to identify structural surface shear slip trend based on structural surface activation data to obtain structural surface shear slip data. The surrounding rock stability analysis module is used to identify potentially unstable blocks based on structural surface shear slip data, obtain potential unstable block data; detect wedge-shaped blocks based on the potential unstable block data, obtain wedge-shaped block data; assess the instability risk level based on the wedge-shaped block data; and perform surrounding rock stability analysis based on the instability risk level to generate surrounding rock stability data.

Images