Rock slope support model construction scheme generation method and device, equipment and medium

Abstract

The invention relates to a rock slope support model construction scheme generation method and device, equipment and a medium. According to the method, a three-dimensional geological model fusing geological information and rock mass parameters is constructed, an initial damage field is obtained in combination with micro-seismic monitoring data and statistical learning inversion, and then a constitutive model capable of reflecting the damage and plastic coupling evolution law is established and calibrated; the model is used for dynamically predicting a spatio-temporal evolution path of a potential slip plane in the excavation process, the supporting opportunity and position are accurately judged based on the stress and damage state in the path, a spatio-temporal sequence scheme is generated, and finally optimal supporting parameters are solved through the multi-objective optimization model. And finally, a set of dynamic support construction scheme capable of actively controlling damage development and giving consideration to safety and economical efficiency is integrated and output, technical spanning from passive reinforcement to active intervention and from static design to dynamic optimization is achieved, and the accuracy and reliability of slope support are effectively improved.

Description

Technical Field

[0001] This invention belongs to the field of geotechnical engineering and digital twin technology, and in particular relates to a method, device, equipment and medium for generating construction schemes for rock slope support models. Background Technology

[0002] Rock slope engineering is a crucial component of infrastructure construction, including water conservancy, hydropower, and road transportation, and its stability directly impacts project safety and economic benefits. Traditional slope support design primarily relies on static stability analysis, employing limit equilibrium methods or strength reduction methods to calculate a single safety factor, neglecting the spatiotemporal evolution of damage during excavation. Practice shows that this static design approach often leads to two outcomes: firstly, premature or excessive support results in economic waste; secondly, delayed support fails to effectively suppress damage propagation, triggering slope instability.

[0003] With the development of monitoring technology and numerical simulation methods, the field of slope engineering has gradually recognized the importance of considering the damage evolution process. During slope deformation, the soil in the shear zone undergoes damage softening, manifested as a gradual decrease in shear strength with increasing deformation. The spatiotemporal distribution patterns of microfractures in rock mass during excavation, captured by microseismic monitoring technology, confirm that slope damage evolution exhibits significant nonlinear and time-dependent characteristics. This information provides a theoretical basis for establishing support design methods based on damage evolution laws.

[0004] However, existing technologies still face significant challenges: first, the quantitative characterization of the damage evolution process is insufficient, making it difficult to accurately predict the damage path; second, the quantitative relationship between damage state and support parameters is unclear, resulting in a lack of scientific basis for determining the timing and location of support. Therefore, an innovative method is urgently needed to predict the damage evolution path and optimize the spatiotemporal sequence of support. Summary of the Invention

[0005] Therefore, it is necessary to provide a method, device, equipment, and medium for generating construction schemes for rock slope support models to address the aforementioned technical problems.

[0006] Firstly, this application provides a method for generating a construction scheme for a rock slope support model, including:

[0007] S1. Construct a three-dimensional geological model of the slope based on geological survey data, and combine the Kriging space interpolation algorithm to allocate discrete rock mass mechanical parameters to each unit grid of the three-dimensional geological model of the slope to generate a rock mass parameter field.

[0008] S2. Collect rock mass microfracture data using a pre-deployed microseismic monitoring network. Based on the spatial distribution characteristics of the rock mass microfracture data and rock mass parameter field, invert the initial damage state of the rock mass using a statistical learning algorithm to generate the initial damage field.

[0009] S3. Based on the damage distribution characteristics of the initial damage field, establish an anisotropic damage constitutive model that considers the coupling between damage and plasticity; perform back-analysis calibration of the model parameters of the anisotropic damage constitutive model using experimental data to generate a damage evolution prediction model with calibrated parameters.

[0010] S4. Based on the damage evolution prediction model, the spatiotemporal evolution trajectory of the potential slip surface is predicted by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during the graded excavation process.

[0011] S5. Based on the stress state and damage degree in the spatiotemporal evolution trajectory, generate a spatiotemporal sequence scheme for support by determining the support timing and calculating the support location priority.

[0012] S6. With safety, economy and construction period as multiple objectives, and based on the support parameter constraints of the support spatiotemporal sequence scheme, establish a multi-objective optimization model; by solving the optimal support parameters of the multi-objective optimization model, a multi-objective optimized support scheme is formed.

[0013] S7 integrates spatiotemporal sequence schemes for support and multi-objective optimized support schemes, and outputs construction plans.

[0014] Secondly, this application also provides a device for generating a construction scheme for a rock slope support model, used to implement the method described in the first aspect, the device comprising:

[0015] The geological model construction and parameter allocation module is used to construct a three-dimensional geological model of the slope based on geological exploration data. It combines the Kriging space interpolation algorithm to allocate discrete rock mass mechanical parameters to each unit grid of the three-dimensional geological model of the slope, generating a rock mass parameter field.

[0016] The damage state dynamic inversion module is used to collect rock mass microfracture data using a pre-deployed microseismic monitoring network. Based on the spatial distribution characteristics of the rock mass microfracture data and rock mass parameter field, the initial damage state of the rock mass is inverted through statistical learning algorithms to generate the initial damage field.

[0017] The damage constitutive model calibration module is used to establish an anisotropic damage constitutive model that considers the coupling between damage and plasticity based on the damage distribution characteristics of the initial damage field; the model parameters of the anisotropic damage constitutive model are calibrated by back-analysis using experimental data, and a damage evolution prediction model with calibrated parameters is generated.

[0018] The slope stability prediction module is used to predict the spatiotemporal evolution trajectory of potential slip surfaces by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during the graded excavation process, based on the damage evolution prediction model.

[0019] The intelligent decision-making module for support timing is used to generate a spatiotemporal sequence scheme for support by determining the support timing and calculating the support location priority based on the stress state and damage degree in the spatiotemporal evolution trajectory.

[0020] The multi-objective optimization decision module is used to establish a multi-objective optimization model based on the support parameter constraints of the support spatiotemporal sequence scheme with safety, economy and construction period as multiple objectives; and to form a multi-objective optimized support scheme by solving the optimal support parameters of the multi-objective optimization model.

[0021] The intelligent integration module for construction schemes is used to integrate spatiotemporal sequence schemes for support and multi-objective optimized support schemes, and output construction schemes.

[0022] Thirdly, this application also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement a method for generating a construction scheme for a rock slope support model as described in the first aspect.

[0023] Fourthly, this application also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements a method for generating a construction scheme for a rock slope support model as described in the first aspect.

[0024] The aforementioned method, device, equipment, and medium for generating a construction scheme for rock slope support model, constructs a three-dimensional geological model by integrating geological survey data and spatial interpolation algorithms. It then inverts the initial damage field based on microseismic monitoring data and machine learning methods, establishes a constitutive model considering the plastic coupling effect of damage, and calibrates its parameters using experimental data. This forms a numerical model capable of accurately predicting the damage evolution path during excavation, dynamically identifying the development trajectory of potential slip surfaces, and formulating graded support strategies based on the spatiotemporal variation of stress state and damage degree. Finally, a multi-objective optimization algorithm balances economy and construction efficiency while ensuring safety, ultimately generating a forward-looking and adaptive dynamic support scheme. This represents a fundamental shift from traditional static experience-based design to dynamic and precise control based on damage evolution mechanisms. It significantly improves the scientific nature of support timing, the accuracy of support location, and the economic rationality of support parameters, effectively overcoming resource waste caused by premature support and control failure caused by late support. Therefore, it fundamentally improves the stability and controllability of slope engineering, construction safety, and the economic benefits throughout its entire life cycle. Attached Figure Description

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

[0026] Figure 1 A flowchart illustrating a method for generating a construction scheme for a rock slope support model provided by the present invention;

[0027] Figure 2 This is a schematic diagram illustrating the process of predicting the spatiotemporal evolution trajectory of a potential slip surface in one optional embodiment of the present invention.

[0028] Figure 3 This is a schematic diagram of a device for generating construction schemes for rock slope support models provided by the present invention. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0030] refer to Figure 1 The document presents a flowchart illustrating a method for generating a construction scheme for rock slope support model, as provided in this application. The method includes the following steps:

[0031] S1. Construct a three-dimensional geological model of the slope based on geological survey data, and combine the Kriging space interpolation algorithm to allocate discrete rock mass mechanical parameters to each unit grid of the three-dimensional geological model of the slope, thereby generating a rock mass parameter field.

[0032] Specifically, the first step is to collect geological survey data, covering the entire slope engineering area. The types of data collected include borehole data, geophysical data, and laboratory test data. Borehole data is used to record lithological stratification and rock quality-related indicators, geophysical data is used to obtain rock wave velocity and delineate lithological interfaces, and laboratory test data is used to test mechanical parameters of the rock mass, such as uniaxial compressive strength, elastic modulus, Poisson's ratio, cohesion, and internal friction angle.

[0033] When constructing the 3D geological model of the slope, professional geotechnical engineering numerical simulation software was selected. First, borehole data was imported to establish stratigraphic control points. Then, the Delaunay triangulation algorithm was used to generate stratigraphic interfaces, followed by unit mesh generation. The mesh size was determined according to the slope scale, and mesh refinement was performed at lithological stratification interfaces to accurately reflect the geological interface characteristics. The model boundary conditions were set to conform to the actual stress state of the project. The bottom boundary fixed the three-dimensional displacement, the two side boundaries restricted the horizontal displacement, and the top was set as a free surface. The initial geostress field was assigned values ​​according to a specific formula, which is: ,in The initial vertical ground stress, Let H be the unit weight of the rock mass and H be the burial depth of the calculation point; the initial horizontal in-situ stress is calculated as follows: Calculation, where denoted as the initial horizontal geostress, and K as the horizontal geostress coefficient, which is determined based on regional geological data.

[0034] The implementation of the Kriging space interpolation algorithm consists of four steps:

[0035] The first step is to preprocess the discrete rock mass mechanical parameters and use the Grubbs test to remove outliers to ensure the reliability of the parameters.

[0036] The second step is to select a variogram model. First, the variogram characteristics of the parameter space are determined by calculating the experimental variogram. The formula for calculating the experimental variogram is: ,in Here, N(h) represents the experimental variability function value, h is the distance between sample points, and N(h) is the number of sample pairs with a distance of h. For in position Rock mass mechanical parameter values ​​at the location, Let be the spatial coordinates of the i-th sample point. To and The spatial coordinates of sample points at a distance h are determined; subsequently, the goodness of fit of the spherical model, exponential model, and Gaussian model are compared, and the optimal model is selected based on the minimum root mean square error. The expression for the spherical model is: when... hour, ,when hour, In the formula The value of the variogram. C is the nugget value, reflecting measurement error and variations at minute scales; C is the partial sill value, reflecting variations at the regional scale. For variable range, it reflects the maximum distance related to the parameter space.

[0037] The third step involves interpolating the mesh of each element in the 3D model, using the element center as the target point, setting a reasonable search radius to ensure sufficient coverage of sample points, and constructing the Kriging equation system. Where K is the spatial correlation matrix of the sample points, Let I be the Lagrange multiplier, and I be the identity matrix. This is a vector of weight coefficients. For a constant vector, after solving for the weight coefficients using the Cholesky decomposition method, according to... The rock mass mechanical parameter values ​​at the center of the calculation unit are as follows: For target point The interpolation parameter value at point n, where n is the number of sample points involved in the interpolation. Let be the weight coefficient of the i-th sample point.

[0038] After generating the rock mass parameter field in the fourth step, three-dimensional cloud maps of each parameter are output for verification. The verification includes the uniformity of parameter variation within the same lithological layer and the abrupt change amplitude of parameters between different lithological layers, ensuring that the parameter field can truly reflect the spatial distribution of rock mass mechanical properties.

[0039] S2. Collect rock mass microfracture data using a pre-deployed microseismic monitoring network. Based on the spatial distribution characteristics of the rock mass microfracture data and rock mass parameter field, invert the initial damage state of the rock mass using a statistical learning algorithm to generate the initial damage field.

[0040] Specifically, the number and distribution of sensors are determined based on the slope's planar dimensions. Three-component microseismic sensors are selected to capture microseismic signals from different directions. The sensors are buried inside the slope to avoid surface interference, and the spacing between adjacent sensors and the monitoring blind zone area are controlled within a reasonable range. The sensors are connected to the data acquisition instrument via shielded cables, and appropriate data storage intervals are set to continuously collect background microfracture data for a period of time prior to excavation.

[0041] Rock mass microfracture data processing includes three key steps: first, microseismic event localization, which involves establishing a set of localization equations based on the arrival time difference between P-waves and S-waves. The equations are expressed as follows: ,in Let be the wave propagation time when the i-th sensor receives the j-th microseismic event. , Let i be the planar coordinates of the i-th sensor. , Let J be the plane coordinates of the j-th microseismic event. The propagation speed of the P-wave. The propagation speed of the S-wave, First, to determine the initial time of the microseismic event, the least squares method is used to solve the system of equations to determine the spatial coordinates of the microseismic event. Second, the magnitude of the microseismic event is calculated using a simplified Richter scale method commonly used in engineering, with the following formula: ,in The Richter magnitude, representing a microseismic event, is an indicator used to quantify the energy of a microseismic event and has no unit. A represents the maximum amplitude of the microseismic signal recorded by the sensor, measured in millimeters. It is read through the waveform recording function of the microseismic data acquisition instrument, and is usually taken as the average amplitude of three consecutive valid waveforms to reduce errors. The distance between the microseismic event and the sensor, i.e., the epicentral distance, is expressed in kilometers and can be calculated using the coordinates of the previously deployed sensors and the location results of the microseismic event. C is a correction factor used to compensate for the attenuation effect of different rock media on wave propagation. For common sedimentary rock slopes, C is typically taken as 0.5~1.0, and for igneous rock slopes, it is taken as 1.0~1.5. Specific values ​​can be determined by referring to monitoring data from similar projects. Thirdly, the microseismic event fracture volume is estimated, calculated based on the correlation between microseismic energy and rock fracture volume, using the formula: Where V is the volume of rock mass fracture caused by the microseismic event, in cubic meters; E is the energy released by the microseismic event, in joules, which can be obtained from the amplitude A and epicentral distance monitored in the early stage. Derivation, specifically according to Calculation (this formula is used to convert magnitude into energy, which conforms to the conventional estimation logic of rock mass microseismic energy); k is the energy conversion coefficient, which reflects the proportion of microseismic energy converted into rock mass fracture energy, and is determined according to the integrity of the rock mass; The value represents the uniaxial compressive strength of the rock mass, measured in Pascals. It is determined through laboratory tests, and the value is consistent with the mechanical parameters in the previous rock mass parameter field to ensure the consistency of the calculation basis.

[0042] The initial damage field inversion uses the random forest algorithm, and the implementation steps are as follows:

[0043] The first step is to divide the analysis into units, with the unit size consistent with the unit size of the 3D geological model. Microfracture characteristic parameters and rock mass parameters are then statistically analyzed for each unit to form an input feature vector. The microfracture characteristic parameters include the number of microseismic events within the unit, the average magnitude, and the fracture volume density. The formula for calculating the fracture volume density is as follows: , The fracture volume density, This is the sum of the rupture volumes of all microseismic events within the element.

[0044] The second step is to obtain training samples. Representative specimens are selected from the borehole cores, and the initial damage value of the specimens is determined through uniaxial compression tests. The formula for calculating the initial damage value is: Where D is the initial damage value, The peak strength of the rock mass after damage. The peak strength of the intact rock mass is represented by the measured initial damage value as the output label.

[0045] The third step is to construct a random forest model, setting model parameters such as the number of decision trees and the maximum tree depth. Cross-validation is used to optimize these parameters. The model training objective is to minimize the mean squared error between the predicted damage value and the experimentally measured damage value. The formula for mean squared error is: ,in Let n be the mean squared error, and n be the sample size. Let i be the damage value predicted by the model for the i-th sample. The damage value of the i-th sample was measured in the experiment.

[0046] The fourth step involves substituting the input feature vectors of all analysis units into the trained model to invert and obtain the initial damage value for each unit. The initial damage value ranges from 0 to 1, where 0 represents intact rock mass and 1 represents completely destroyed rock mass. Based on this, an initial damage field is generated. After the initial damage field is generated, it is validated. The validation includes the spatial overlap between high-damage areas and densely micro-fractured areas, the correspondence between low-damage areas and areas without micro-fractures, and the matching between high-damage areas and low-strength rock mass areas, to ensure that the inversion results are reasonable and reliable.

[0047] S3. Based on the damage distribution characteristics of the initial damage field, establish an anisotropic damage constitutive model that considers the coupling between damage and plasticity; perform back-analysis calibration of the model parameters of the anisotropic damage constitutive model using experimental data to generate a damage evolution prediction model with calibrated parameters.

[0048] Specifically, the anisotropic damage constitutive model is established based on tensor theory to reflect the anisotropic characteristics of rock mass damage and the coupling effect between damage and plasticity. The damage state in the model is described using a second-order symmetric damage tensor, the expression of which is: ,in These are the components of the second-order damage tensor. This is the damage evolution coefficient, whose value varies with the equivalent plastic strain. For equivalent plastic strain, , All of these are components of the normal vector of the micro-fracture surface, which is determined statistically from microseismic monitoring data, taking the dominant fracture direction of the micro-fracture event.

[0049] Plastic deformation is described using the Drucker-Prager yield criterion, and the expression for the yield function is: Where f is the yield function value, k is the yield parameter. As the first invariant of stress, Let be the second invariant of the stress deviatoric tensor. The formula for calculating the first invariant of stress is: , , , These are the three principal stresses; the formula for calculating the second invariant of the stress deviator tensor is... To reflect the coupling between damage and plasticity, the yield parameter is correlated with the damage state, and the correlation formula is as follows: , ,in , Let D be the initial yield parameter (i.e., the yield parameter when the damage is 0), m be the damage variable, and m be the coupling coefficient. The damage variable is calculated from the trace of the damage tensor, using the following formula: , , , These are the components of the damage tensor along the three coordinate axes.

[0050] The elastic deformation characteristics of rock mass take into account the influence of damage, and the elastic matrix adopts a modified form, the modified formula is as follows: ,in These are the components of the corrected elasticity matrix. Let be the components of the initial elastic matrix (i.e., the elastic matrix when damage is 0), and n be the elastic damage coefficient. The expression for the initial elastic matrix is: Where E is the initial elastic modulus. Poisson's ratio, , , , , , All are Kronecker symbols. When the symbols have the same subscript, the value is 1; when the subscripts are different, the value is 0.

[0051] The calibration of model parameters combines indoor experiments and back analysis methods. The specific steps are as follows: First, prepare test samples. Select representative rock masses from the slope and process them into standard specimens. The number of specimens in each group must meet the experimental accuracy requirements. Conduct triaxial compression-unloading tests under different confining pressures, and collect stress-strain curves, volumetric strain curves, and acoustic emission data during the test. The acoustic emission data is used to simultaneously monitor the rock mass damage evolution process. Second, determine the value range of the parameters to be calibrated. These parameters include initial elastic modulus, Poisson's ratio, initial yield parameter, coupling coefficient, elastic damage coefficient, initial damage coefficient, and damage evolution index. The parameter range needs to be initially set based on the rock mass type and engineering experience. Third, use a particle swarm optimization algorithm for parameter back analysis, using the mean square error between the experimental curve and the model calculation curve as the objective function. The objective function formula is: ,in Let i be the test stress value at the i-th data point. The stress value is calculated for the model at the i-th data point, where n is the number of data points. Algorithm parameters such as the number of particles, number of iterations, inertia weight, and acceleration factor are set, and the parameter combination that minimizes the objective function is obtained through iterative calculation. The fourth step is to verify the parameters by selecting confining pressure conditions that were not included in the calibration to carry out verification tests, and calculating the relative error between the model prediction value and the experimental measurement value. If the relative error meets the requirements, the parameter calibration is completed, and a damage evolution prediction model that has been calibrated is generated accordingly.

[0052] S4. Based on the damage evolution prediction model, the spatiotemporal evolution trajectory of the potential slip surface is predicted by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during the graded excavation process.

[0053] Specifically, the staged excavation simulation is matched with the actual construction progress of the project. First, the number of excavation stages and the height of each stage are divided according to the total excavation height of the slope, and the time interval for each stage is set to ensure that the simulation process conforms to the actual construction rhythm. The three-dimensional geological model of the slope constructed in the early stage and the damage evolution prediction model after parameter calibration are imported into the numerical simulation software, and the calculation parameters are set, including the time step and convergence criteria. The time step needs to ensure the stability of the numerical calculation, and the convergence criterion uses the unbalanced force ratio as the control index.

[0054] The stress redistribution simulation is performed step-by-step according to the excavation stages. Before each excavation stage, an initial stress equilibrium calculation is performed on the model, and the calculation is iterated until the convergence criterion is met. The stress field after equilibrium is recorded, including the magnitude and direction of the three principal stresses of each element. During the excavation process, specific methods are used to delete elements in the excavation area or reduce the stiffness of the elements. After deletion or stiffness reduction, the stress calculation is re-performed to analyze the stress transfer law caused by excavation. The focus is on the stress release zone near the excavation slope surface and the stress concentration zone at the slope toe and inside the slope body. The stress difference and the change in the direction of the principal stress of each element after each excavation stage are recorded. The formula for calculating the stress difference is: , For stress difference, For the maximum principal stress, It is the minimum principal stress.

[0055] Damage evolution simulation is conducted based on a damage evolution prediction model. The equivalent plastic strain of each element is calculated according to the stress state after each excavation stage. The formula for calculating the equivalent plastic strain is as follows: ,in The plastic strain tensor components are used; the equivalent plastic strain is substituted into the damage evolution formula to calculate the damage value of the element. The damage evolution formula is as follows: ,in The damage evolution index (determined through parameter calibration) is set; a damage threshold is set, and when the damage value of a unit reaches or exceeds the threshold, the unit is marked as a high-damage unit. The number and spatial distribution range of high-damage units after each level of excavation are counted, and the expansion pattern of high-damage units is analyzed.

[0056] The prediction of the spatiotemporal evolution trajectory of potential slip surfaces adopts a dual criterion of "damage connectivity + stress yielding," and the implementation steps are as follows: First, a connectivity analysis is performed on the high-damage units after each level of excavation to construct a unit adjacency matrix. Matrix analysis is used to determine whether high-damage units form continuous damage zones. The determination of continuous damage zones must meet the requirements of the number of units and the lateral penetration length. Second, a stress yielding determination is performed on the formed continuous damage zones (i.e., the initial slip surface shape). The yield function value of each unit in the initial slip surface shape is calculated, and the average yield function value of each unit is taken as the judgment index. The formula for calculating the average yield function value is as follows: ,in Let n be the average yield function, and n be the number of elements in the initial slip surface shape. Let be the yield function value of the i-th element. When the average yield function value satisfies the yield condition, the continuous damage zone is confirmed as a potential slip surface. The third step is to record the spatiotemporal parameters of the potential slip surface. The spatial parameters include the spatial coordinates of the slip surface's start and end points, the slip surface's dip angle, and its length. The formula for calculating the dip angle is... , The angle of inclination of the slip surface. The vertical height difference between the starting and ending points of the slip surface. The formula for calculating the length of the horizontal distance is: L is the length of the slip surface. The lateral distance between the start and end points of the slip surface; time parameters include the number of excavation stages and days corresponding to the formation of the slip surface; the fourth step uses the simplified Bishop method to calculate the safety factor of the potential slip surface, and the safety factor formula is: ,in For safety factor, c is the rock mass cohesion. The length of the slip surface of the segmented slip body. The weight of the segments of the sliding body. The angle of inclination of the slip surface. The friction angle within the rock mass is used. When the safety factor is lower than the allowable value specified in the code, the potential slip surface is marked as a high-risk slip surface, and its spatiotemporal evolution trajectory is closely tracked.

[0057] S5. Based on the stress state and damage degree in the spatiotemporal evolution trajectory, generate a spatiotemporal sequence scheme for support by determining the timing of support and calculating the priority of support location.

[0058] Specifically, the timing of support determination is based on the stress state and damage degree of the potential slip surface, establishing a three-level determination index system. The first level is the damage development rate index, calculated as the ratio of the average damage increment of the slip surface region between adjacent excavation stages to the excavation time interval. This ratio is the damage development rate, expressed by the formula: ,in The average damage increment, The first level is the excavation interval, which triggers an early warning when the damage development rate exceeds a set threshold; the second level is the stress yield degree index, which uses the ratio of the average yield function value of the slip surface to the yield critical value as the judgment basis. The yield critical value is 0. When the ratio reaches the set threshold, it indicates that the stress of the slip surface is close to the yield state; the third level is the safety factor index. When the safety factor of the potential slip surface drops to a specific range of the allowable value in the specification, support preparation needs to be initiated.

[0059] The weighted comprehensive judgment method is used to calculate the support timing judgment value. For example, the formula for calculating the judgment value can be:

[0060]

[0061] Where S is the value for determining the timing of support. As a warning threshold for the rate of damage development, This is the threshold for determining the degree of stress yielding. To standardize the permissible safety factor, The threshold for the safety factor is set; a critical value for the judgment value is set, and when the judgment value reaches or exceeds the critical value, it is determined to be the appropriate time for support, and the time limit for support construction is specified.

[0062] The priority calculation of support location adopts a multi-factor weighted method. First, the slope is divided into several support units, and the size of each unit is set. Three key factors influencing the priority of support location are determined: 1) the distance coefficient, i.e., the shortest distance from the center of the support unit to the potential slip surface; 2) the damage coefficient, i.e., the average damage value of the support unit; and 3) the stress concentration coefficient, i.e., the ratio of the maximum principal stress of the support unit to the compressive strength of the rock mass. Each factor is quantified: the quantification formula for the distance coefficient is... ,in This is the quantified value of the distance coefficient, where d is the shortest distance from the center of the support unit to the slip surface. The maximum burial depth of the slope; the quantitative formula for the damage coefficient is: ,in This is the quantified value of the damage coefficient. This represents the average damage value of the support unit. This is the critical damage value; the quantitative formula for the stress concentration factor is: ,in This is the quantified value of the stress concentration factor. ( For the maximum principal stress, (Rock mass compressive strength) This is the threshold for determining stress concentration.

[0063] The formula for calculating the priority of support location can be: P represents the priority of the support location; priority levels are divided according to the priority value, including emergency support level, routine support level and delayed support level, and different priority levels correspond to different support urgency.

[0064] The spatiotemporal sequence scheme for support integrates three-dimensional information: time, location, and type. The time dimension categorizes support into emergency, conventional, and delayed support based on the timing of support decisions, clearly defining the construction time windows for each type. The location dimension prioritizes support resources based on location priority, ensuring critical areas receive support first. Support type is determined by the stress state, damage level, and priority level of the support location: high-priority locations utilize a combination of anchor bolts and shotcrete, medium-priority locations use anchor bolts, and low-priority locations use shotcrete. The scheme clarifies the parameter relationships between different support types, such as the relationship between anchor bolt length and slip surface depth, the range of anchor bolt spacing, and the setting of shotcrete thickness. Simultaneously, the construction efficiency of each support procedure is defined to ensure that the support construction progress matches the excavation progress, avoiding delays.

[0065] S6. With safety, economy and construction period as multiple objectives, a multi-objective optimization model is established based on the support parameter constraints of the support spatiotemporal sequence scheme; by solving the optimal support parameters of the multi-objective optimization model, a multi-objective optimized support scheme is formed.

[0066] Specifically, regarding the safety objective, the goal is to maximize the minimum safety factor of the slope after support. The safety factor after support takes into account the contribution of the support structure, and the calculation formula is as follows:

[0067]

[0068] in, For the tensile strength of the anchor bolt, The angle between the anchor bolt and the sliding surface; the formula for calculating the tensile strength of the anchor bolt can be... Where d is the diameter of the anchor bolt. The anchor bolt yield strength is used as the constraint condition; the minimum safety factor after support must not be lower than the allowable value specified in the standard.

[0069] With an economic objective in mind, the goal is to minimize the total cost of the support structure. This total cost includes the cost of anchor bolts, shotcrete, and construction machinery and labor. The formula for the total cost is: ,in For anchor bolt costs, To reduce the cost of shotcrete, The costs include construction machinery and labor; the formula for calculating anchor bolt costs is:

[0070]

[0071] Where s is the anchor bolt spacing. Where L is the anchor bolt spacing, and L is the anchor bolt length. For the density of steel, This refers to the unit price of steel. The cost is for anchor bolt installation; the formula for calculating the cost of shotcrete is:

[0072]

[0073] Where h is the thickness of the shotcrete. For concrete density, This refers to the unit price of concrete. The cost of spraying construction; the formula for calculating the costs of construction machinery and labor is:

[0074]

[0075] in, For construction machinery rental fees, This refers to labor costs.

[0076] Regarding the project schedule objective, the goal is to minimize the total duration of the support structure project. The total duration is the maximum of the anchor bolt construction duration and the shotcrete construction duration, expressed by the formula: ,in For the anchor bolt construction period, The construction period for shotcrete is [specified]; the formula for calculating the construction period for anchor bolts is [formula]. The formula for calculating the construction period of shotcrete is: Set constraints for the project duration target, namely, the total project duration must not exceed the set upper limit.

[0077] The constraints on support parameters are set based on relevant standards and engineering practices, including:

[0078] 1) Anchor bolt parameter constraints: The length of the anchor bolt must be no less than a specific multiple of the slip surface depth and no more than the maximum allowable length to avoid construction difficulties; the anchor bolt spacing must be no less than the minimum allowable spacing to prevent rock splitting and no more than the maximum allowable spacing to prevent support failure; the anchor bolt diameter must be no less than the minimum allowable diameter to meet the tensile bearing capacity requirements and no more than the maximum allowable diameter to prevent increased drilling difficulty.

[0079] 2) Shotcrete parameter constraints: The thickness of shotcrete must be no less than the minimum allowable thickness to meet crack resistance requirements, and must not exceed the maximum allowable thickness to prevent detachment; the strength grade of shotcrete must be no less than the minimum strength grade requirement.

[0080] 3) Coupling constraints: When the anchor spacing exceeds a specific value, the thickness of the shotcrete must not be less than the corresponding minimum thickness; the total support period must not exceed the set upper limit of the period.

[0081] The NSGA-II algorithm is used to solve the multi-objective optimization model. The solution steps are as follows: First, the support parameters are encoded into real numbers to form chromosomes. Chromosomes include parameters such as anchor length, anchor spacing, anchor diameter, shotcrete thickness, and concrete strength grade. The population size is then set. Second, the population is initialized so that the parameter values ​​of individuals within the population are randomly distributed within the constraint range. Third, the fitness of each individual is calculated. Fitness is determined based on the objective function value and the satisfaction of the constraints. Individuals that satisfy the constraints proceed to the sorting stage, while individuals that do not satisfy the constraints have their fitness set to 0. Fourth, non-dominated sorting is performed, dividing the population into different non-dominated layers. The crowding distance between individuals within each non-dominated layer is calculated. The crowding distance measures the dispersion of individuals in the solution space. Fifth, selection, crossover, and mutation operations are used to generate... For the offspring population, the selection operation uses the roulette wheel method, prioritizing individuals with high non-dominated levels and large crowding distances. The crossover operation uses the SBX crossover method, setting the crossover probability and distribution index. The mutation operation uses the multinomial mutation method, setting the mutation probability and distribution index. In the sixth step, the parent and offspring populations are merged, and the non-dominated sorting and crowding distance calculations are performed again. A specific number of individuals are selected to form a new parent population. This iterative process is repeated until the convergence condition is met. The convergence condition is controlled by the change range of the optimal solution over several generations. In the seventh step, the ideal point method is used to select the optimal solution from the final set of non-dominated solutions. The standardized Euclidean distance from each non-dominated solution to the ideal point is calculated. The objective function values ​​of the ideal point are set as the theoretical optimal values. The solution with the smallest distance is selected as the optimal support parameter, thereby forming a multi-objective optimized support scheme.

[0082] S7 integrates spatiotemporal sequence schemes for support and multi-objective optimized support schemes, and outputs construction plans.

[0083] Specifically, the first step is to match the spatiotemporal parameters. Using the "time-location" unit in the spatiotemporal sequence of the support scheme as an index, the optimal support parameters in the multi-objective optimized support scheme are embedded into the corresponding unit to form a "time-location-parameter" correspondence table. The parameter matching results are then verified, with a focus on verifying the compatibility of the support parameters with the geological conditions, damage state, and stress requirements of the support location. For example, whether the anchorage length of the anchor bolt meets the specifications, and whether the thickness of the shotcrete matches the degree of stress concentration, etc., to ensure the rationality of parameter application.

[0084] The second step involves collaborative process design. Based on a "time-location-parameter" correspondence table, a corresponding construction process flow is developed. This flow includes procedures such as surveying and setting out, drilling, hole cleaning, anchor installation, grouting, and shotcrete. The operational requirements for each procedure are clearly defined: surveying and setting out requires the use of professional surveying equipment to ensure hole accuracy; drilling requires the selection of appropriate drilling equipment to control the hole diameter and inclination angle; hole cleaning requires specific methods to remove impurities and ensure cleanliness; anchor installation requires mechanical pushing to control the exposed length of the anchors; grouting requires the selection of appropriate grouting materials and proportions, and control of grouting pressure and fullness; shotcrete requires specific processes to control the thickness and interval of each layer. Simultaneously, the time nodes for each procedure are clearly defined, and the construction time for each procedure is calculated to ensure the overall project duration meets requirements.

[0085] The third step is to verify the scheme. Numerical simulation is used to verify the safety of the support scheme. The integrated support parameters are substituted into the three-dimensional geological model of the slope to simulate the stress distribution and damage evolution process after support. The minimum safety factor and damage growth rate after support are calculated to verify whether the safety target is met. Typical areas are selected for on-site trial construction to test the performance indicators of the support structure, such as anchor pull-out force and shotcrete strength, to verify the actual effect of the support parameters and the feasibility of the construction technology.

[0086] The final construction plan document contains complete project information, including: ① Project overview, introducing the spatial coordinates, scale, and geological parameters of the slope; ② Overall support design, presenting a time-location-parameter correspondence table of the support spatiotemporal sequence and a table of support parameters; ③ Detailed construction rules for each area, clearly defining the location coordinates, support parameters, and corresponding construction process flowcharts for each area; ④ Quality control requirements, specifying the quality control indicators and testing methods for each process, such as the allowable range of hole position error, grouting fullness requirements, and pull-out force testing frequency; ⑤ Safety measures, including plans for high-altitude operation protection, machinery operation safety, and emergency response; ⑥ Schedule and cost table, detailing the schedule and cost composition of each process; ⑦ Appendices, including technical charts such as the 3D model of the slope, damage field cloud map, and optimized convergence curve. The document must comply with relevant national standards and specifications to ensure that construction personnel can implement the support construction according to the plan without any creative work.

[0087] The aforementioned method for generating a construction scheme for rock slope support model constructs a three-dimensional geological model by integrating geological survey data and spatial interpolation algorithms. It then inverts the initial damage field based on microseismic monitoring data and machine learning methods. Subsequently, a constitutive model considering the plastic coupling effect of damage is established, and its parameters are calibrated using experimental data. This forms a numerical model capable of accurately predicting the damage evolution path during excavation. Based on this, the development trajectory of potential slip surfaces is dynamically identified. A graded support strategy is formulated based on the spatiotemporal variation of stress state and damage degree. Finally, a multi-objective optimization algorithm balances economy and construction efficiency while ensuring safety, ultimately generating a forward-looking and adaptive dynamic support scheme. This represents a fundamental shift from traditional static experience-based design to dynamic and precise control based on damage evolution mechanisms. It significantly improves the scientific nature of support timing, the accuracy of support location, and the economic rationality of support parameters. It effectively overcomes the resource waste caused by premature support and the control failure caused by late support, thereby fundamentally improving the stability and controllability, construction safety, and life-cycle economic benefits of slope engineering.

[0088] In one optional embodiment, rock mass microfracture data are collected using a pre-deployed microseismic monitoring network. Based on the spatial distribution characteristics of the rock mass microfracture data and rock mass parameter field, the initial damage state of the rock mass is inverted using a statistical learning algorithm to generate an initial damage field, including the following steps:

[0089] S11. Through the microseismic monitoring network, monitor the spatiotemporal coordinates, energy release, and magnitude parameters of rock microfractures during the excavation process, and generate the original dataset of microseismic events.

[0090] Specifically, the microseismic monitoring network consists of several three-component microseismic sensors, a data acquisition instrument, and a data transmission module. The sensors are deployed in a preset array inside and on the surface of the slope rock mass to ensure full coverage monitoring of micro-fracture activity in the rock mass within the excavation influence area. During the monitoring process, the sensors capture elastic wave signals generated by micro-fractures in real time, convert the signals into electrical signals, and transmit them to the data acquisition instrument. The acquisition instrument performs analog-to-digital conversion and stores the signals according to a set sampling frequency.

[0091] The acquisition of spatiotemporal coordinates is based on the principle of elastic wave propagation time difference positioning: Multiple sensors record the arrival times of the P-wave and S-wave of the same micro-fracture event, establishing a correlation equation between the time difference and the spatial position of the sensors. Solving this equation yields the three-dimensional spatial coordinates (x, y, z) and the time t of the micro-fracture event, forming the event's spatiotemporal information. Energy release is calculated by analyzing the amplitude, frequency, and duration of the microseismic waveform. Based on the elastic wave energy propagation law, the waveform characteristic parameters are converted into the energy value released during the rock mass micro-fracture process. This energy value reflects the scale of the micro-fracture. Magnitude parameters are calculated using the standard magnitude calculation method based on waveform amplitude. By measuring the maximum amplitude of the waveform within a specific frequency range and correcting for the distance between the sensor and the micro-fracture event (epochal distance), a magnitude value characterizing the intensity of the micro-fracture is obtained.

[0092] The raw dataset of microseismic events contains complete information for each micro-fracture event, including: event number, occurrence time, three-dimensional spatial coordinates, energy release, magnitude, corresponding sensor array response waveform segment, and environmental parameters (such as temperature and humidity) during data acquisition. The dataset is stored in a standardized format to ensure the accessibility and consistency of parameters during subsequent data processing, providing fundamental data support for subsequent spatial clustering analysis and damage variable calculation.

[0093] S12. Based on the original dataset of microseismic events and the rock mass parameter field, rock mass type partitioning is performed. The spatial clustering characteristics of microseismic events are analyzed by geostatistical variation function to generate regional division schemes for different damage levels.

[0094] Specifically, the rock mass type is first divided into zones based on the rock mass parameter field. According to the spatial distribution differences of rock mass mechanical parameters (such as elastic modulus, cohesion, and internal friction angle), the three-dimensional geological model of the slope is divided into several rock mass type zones. The mechanical parameters of the same rock mass are similar, while the parameters of different rock masses are significantly different. The zoning results must be consistent with the actual geological stratification characteristics of the slope to ensure that subsequent analysis can reflect the differences in the response of different rock masses to micro-fracture activities.

[0095] Geostatistical variogram analysis is used to quantify the spatial distribution characteristics of microseismic events. The core function of the variogram is to describe the degree of variation in the distribution of microseismic events at different spatial distances. Its basic expression is: ,in Let N(h) be the variogram value at a distance of h, where h is the distance between two points in space, and N(h) is the number of microseismic event pairs at a distance equal to h. For position Microseismic event density at a location (number of microseismic events per unit volume). Let i be the spatial location of the i-th microseismic event. To and Spatial location of the microseismic event at a distance h.

[0096] During the calculation process, the variation function is calculated separately for each rock mass type partition: the microseismic event density data in each partition is gridded, an appropriate distance interval h is determined (set according to the partition size and event density), the number of event point pairs N(h) at each distance interval is counted, and the corresponding values ​​for each distance interval are calculated by substituting them into the formula. Values; by fitting the variogram curve, three key parameters are obtained: range, sill value, and nugget value. The range reflects the maximum distance of spatial correlation of microseismic events, the sill value reflects the total variability of microseismic event density, and the nugget value reflects measurement error and microscale variation.

[0097] The spatial clustering characteristics of microseismic events are determined based on the parameters of the variogram: when the range is small and the sill value is large, it indicates that the microseismic events are concentrated in a small area, forming clustered regions. Combining the microseismic event density, total energy release, and magnitude distribution of the clustered regions, regions with different damage levels are divided: regions with high microseismic event density, large total energy release, and a high proportion of high-magnitude events correspond to high damage levels; regions with low microseismic event density, small total energy release, and a high proportion of low-magnitude events correspond to low damage levels; and regions in between correspond to medium damage levels. The final region division scheme needs to clearly define the spatial extent, boundary coordinates, and corresponding microseismic activity characteristic indicators of each damage level region, providing a basis for region classification for subsequent calculation of grid unit damage variables.

[0098] S13. Map the regional division scheme onto the grid cells of the three-dimensional geological model of the slope. Calculate the initial damage variable for each grid cell based on the number and energy amplitude of microseismic events falling into each grid cell. The formula for calculating the initial damage variable is:

[0099]

[0100] in, is the initial damage variable of the grid cell, characterizing the initial damage degree of the rock mass within the grid cell; N is the number of microseismic events within the grid cell, obtained by statistically analyzing the microseismic events falling into the grid cell from the original microseismic event dataset; The number of failure threshold events for a grid cell is determined based on the rock mass type at the corresponding location of the grid cell in the rock mass parameter field. The energy amplitude of the microseismic event is obtained from the original microseismic event dataset; The reference energy amplitude is determined based on the statistical characteristics of the energy amplitudes of all microseismic events in the original microseismic event dataset.

[0101] Specifically, the different damage level area division schemes generated by S12 are first mapped onto the grid cells of the three-dimensional geological model of the slope according to the spatial coordinate correspondence, ensuring that each grid cell can clearly correspond to a certain damage level area. At the same time, the spatial location, rock mass type and corresponding damage level area category of each grid cell are recorded to establish the relationship between area division and grid cell.

[0102] For each grid cell, the number N of microseismic events falling within its spatial range is counted. The statistical process is based on matching the three-dimensional coordinates of the events in the original dataset of microseismic events with the spatial boundary coordinates of the grid cell. If the coordinates of a certain microseismic event are within the boundary range of a certain grid cell, then the event is determined to fall within that grid cell. The statistical results need to be verified twice to avoid duplicate counting or omission of events.

[0103] Number of events that violate the threshold The determination of the number of microseismic events required for the grid cells to be damaged depends on the rock mass type. Different rock mass types have different damage resistance capabilities, and the number of microseismic events required for significant damage varies. Damage evolution curves of different rock mass types under microseismic action are obtained through indoor rock mass damage tests. The number of microseismic events corresponding to when the rock mass exhibits significant plastic deformation or a strength reduction of a preset percentage (e.g., a 10% strength reduction) is then used as the failure threshold event number for that rock mass type. If a grid cell spans multiple rock mass types, a weighted calculation is performed based on the volume percentage of each rock mass type within the grid cell to obtain the value of that grid cell. value.

[0104] Microseismic event energy amplitude Obtained from the original microseismic event dataset, each microseismic event falling into a grid cell corresponds to an energy amplitude. This value characterizes the amount of energy released by a single microseismic event. Its acquisition requires energy calculation and calibration of the waveform signal to ensure accuracy for different events. They are comparable. The sum of the energy amplitudes of all N microseismic events falling into the grid cell is calculated for each... Validity checks are performed to remove abnormal energy amplitude data caused by signal interference, ensuring the accuracy of the total calculation.

[0105] Reference energy amplitude The reference energy amplitude is determined based on the statistical characteristics of the energy amplitudes of all microseismic events in the original microseismic event dataset. First, a statistical analysis is performed on the energy amplitudes of the entire dataset, calculating the mean, median, mode, and standard deviation. Considering the influence of abnormally high or low energy events on the statistical results, the mean after removing extreme values ​​(e.g., removing the upper and lower 5% extreme values) is used as the reference energy amplitude. ,make sure It can reflect the overall level of energy amplitude in the dataset, avoiding the influence of individual extreme events. It deviates from the actual average level.

[0106] Initial damage variables In the calculation formula, The term reflects the proportion of the number of microseismic events within a grid cell relative to the failure threshold. The larger this ratio, the closer the number of events is to the critical value required for rock mass failure, and the greater their contribution to the damage. The term reflects the ratio of the average energy amplitude of microseismic events within a grid cell to the reference energy amplitude. A larger ratio indicates a higher average energy of the event and more severe damage to the rock mass; the product of the two terms yields... Taking into account the synergistic effect of the number of events and the energy of events on damage, its value ranges from 0 to 1. The closer the value is to 1, the more severe the initial damage to the rock mass within the grid cell. The closer it is to 0, the closer the rock mass is to a complete state.

[0107] S14. Based on the initial damage variables of all grid cells, a spatially continuous damage distribution field is generated using the kernel density estimation algorithm to obtain the initial damage field.

[0108] Specifically, the initial damage variables of all mesh elements are first organized. The corresponding grid cell center coordinates form a discrete set of damage data points, each data point containing three-dimensional coordinates. (Grid cell center coordinates) and corresponding The data point set must ensure that it covers all grid cells of the 3D geological model of the slope, with no missing or duplicate data points.

[0109] Kernel density estimation algorithms are used to convert discrete... The value is transformed into a spatially continuous damage distribution field. Its core principle is to use a kernel function to transform each discrete data point... The values ​​are spatially smoothed to ensure a continuous transition of damage values ​​between adjacent grid cells, avoiding abrupt changes in damage values ​​caused by mesh generation. The basic expression for kernel density estimation is:

[0110]

[0111] Where f(x) is the continuous damage value at spatial location x, n is the number of discrete data points (i.e., the total number of grid cells), h is the bandwidth parameter (controlling the smoothness), d is the spatial dimension (here d=3, corresponding to three-dimensional space), and K(·) is the kernel function. Let i be the center coordinates of the i-th discrete data point. Let be the initial damage variable for the i-th data point.

[0112] The kernel function K(·) must satisfy the requirements of nonnegativity, normalization, and symmetry. Here, the Gaussian kernel function is used, and its expression is:

[0113]

[0114] Where u is the standardized spatial distance vector. Let u be the Euclidean norm; the Gaussian kernel function can achieve a smooth transition of damage values ​​and has a stable calculation process, making it suitable for generating continuous fields in three-dimensional space.

[0115] The bandwidth parameter h needs to be determined while balancing smoothing and detail preservation. If h is too small, the generated damage distribution field will have obvious grid boundary marks and poor continuity; if h is too large, it will mask the detailed features of local high-damage or low-damage areas. The optimal h value is determined using cross-validation: the discrete data point set is randomly divided into a training set and a validation set, and kernel density estimation is performed for different h values. The actual value in the validation set is then calculated. The mean square error between the measured value and the estimated value is used to select the h value with the smallest mean square error as the optimal bandwidth parameter.

[0116] During the calculation process, the three-dimensional geological model of the slope is divided into calculation grids according to the set fineness (the calculation grid can be consistent with the original model grid or finer). For the center position x of each calculation grid, the corresponding continuous damage value f(x) is calculated by substituting it into the kernel density estimation expression. After the calculation is completed, the f(x) values ​​of all calculation grids are spatially interpolated to ensure the continuity and smoothness of the damage values ​​in three-dimensional space.

[0117] Finally, based on the distribution characteristics of continuous damage values ​​f(x), the damage distribution field is visualized to generate a damage distribution cloud map, clarifying the spatial distribution range, shape, and interrelationships of high damage regions (f(x) close to 1), medium damage regions (f(x) between 0.3 and 0.7), and low damage regions (f(x) close to 0). Simultaneously, the rationality of the damage distribution field is verified, checking the consistency between continuous and discrete damage values. The degree of agreement between the values ​​and the spatial correspondence between high-damage areas and microseismic event clustering areas ensure that the generated damage distribution field can truly reflect the spatial distribution characteristics of the initial damage of the rock mass, and finally obtain the initial damage field.

[0118] In one optional embodiment, based on the damage distribution characteristics of the initial damage field, an anisotropic damage constitutive model considering the coupling between damage and plasticity is established; the model parameters of the anisotropic damage constitutive model are calibrated through back-analysis using experimental data to generate a parameter-calibrated damage evolution prediction model, including the following steps:

[0119] S21. Based on the anisotropic damage tensor theory and damage distribution characteristics, a damage evolution equation incorporating the parameters of the material to be calibrated is established by defining the functional relationship between damage energy release rate and plastic strain rate; the damage evolution equation is:

[0120]

[0121] Where D is the damage variable, representing the current degree of damage to the rock mass; Y is the damage energy release rate, calculated based on the current stress state and damage state. and For the parameters of the material to be calibrated; t represents the plastic strain rate, calculated using the plastic flow rule; t represents time.

[0122] Specifically, the anisotropic damage tensor theory describes the direction dependence of rock mass damage using a second-order symmetric tensor. The components of the damage tensor are directly related to the spatial distribution characteristics of microfractures within the rock mass. By analyzing the damage distribution characteristics of the initial damage field, the principal direction and magnitude of the damage tensor at different spatial locations can be determined, providing a basis for the initial damage state in subsequent damage evolution analysis.

[0123] The damage energy release rate Y is a key mechanical quantity characterizing the development of rock mass damage. Its calculation combines the current stress state and the damage state, specifically derived through the coupling relationship between the stress tensor and the damage tensor. For anisotropic damaged rock masses, the stress tensor needs to be corrected for damage; the expression for the corrected effective stress tensor is as follows: ,in For the nominal stress tensor (a second-order tensor, components) In this context, i and j represent the direction of stress and the normal to the surface, respectively. The values ​​1, 2, and 3 correspond to the three-dimensional directions x, y, and z. When i = j, it is normal stress; when i ≠ j, it is shear stress. The effective stress tensor (considering the actual stress borne by the rock mass after damage, also a second-order tensor, with component meanings similar to...) (Consistent), D is the damage variable (scalar, value 0-1, 0 represents intact rock mass, 1 represents complete rock mass failure). Based on the effective stress tensor, the formula for calculating the damage energy release rate Y can be derived through the energy principle. ,in The components of the strain tensor (second-order tensor, and stress tensor) Correspondingly, the values ​​of i and j, 1, 2, and 3, represent the three-dimensional spatial directions of x, y, and z. When i = j, it is normal strain, reflecting the elongation or shortening deformation of the rock mass in that direction; when i ≠ j, it is shear strain, reflecting the shear deformation of the rock mass in the plane formed by the two directions. Represents the components of the strain tensor of the damage variable D. The partial derivatives reflect the changes in damage variables caused by changes in unit strain components. The formula as a whole embodies the energy released by unit damage changes and is directly related to the energy driving relationship between stress state and damage evolution.

[0124] Plastic strain rate The plastic flow law describes the relationship between the increment of plastic strain and the yield function. For the Drucker-Plag yield criterion commonly used in rock masses, the yield function is... ,in The first invariant of stress (scalar, calculated by the formula is) (reflecting the volume change effect of the stress field). The second invariant of the stress deviatoric tensor (scalar, calculated using the formula is) (reflects the shape change effect of the stress field). k is the yield parameter (both are scalars). The stress threshold for plastic yielding of the rock mass is determined by the friction angle within the rock mass and the cohesion of the rock mass. According to the plastic flow law, the plastic strain rate... It is proportional to the partial derivative of the yield function, and the expression is: ,in The plasticity multiplier (a scalar, non-negative, reflecting the degree of development of plastic deformation) The larger the value, the more significant the plastic deformation per unit time. The components of the plastic strain rate tensor (second-order tensor, i, j have the same meaning as...) Consistent, reflecting the change in plastic strain of the rock mass per unit time, i=j is the plastic normal strain rate, and i≠j is the plastic shear strain rate. for The equivalent plastic strain rate (scalar, through) The calculation shows that by converting the three-dimensional plastic strain rate tensor into a single scalar (which facilitates the quantification of the overall development rate of plastic deformation), this parameter reflects the development rate of plastic deformation in the rock mass and is an important driving factor for damage evolution.

[0125] Damage evolution equation middle, The differential increment of the damage variable (a scalar, reflected in a small time interval) The change in the internal injury variable D, This indicates that the damage has worsened. (Indicates damage stability), characterizing the amount of damage change within a small time interval; and The material parameters to be calibrated (all are positive scalars) The initial threshold reflecting the resistance of rock mass to damage development. The larger the value, the more damage energy the rock mass needs to accumulate before significant damage occurs. It reflects the sensitivity of damage evolution to the damage energy release rate. The larger the value, the more likely a small change in Y will cause [a certain effect]. Significant fluctuations); t is a time variable (a scalar that reflects the time process of construction or deformation development). For a small time interval (a scalar, the derivative of the time variable t, used to describe the instantaneous changes in damage evolution), the equation is derived from the damage energy release rate Y and the plastic strain rate. The coupling established a quantitative relationship between damage and time development, while also incorporating the inherent properties of the material ( , As a parameter to be calibrated, it provides a mathematical basis for subsequent parameter calibration based on experimental data.

[0126] S22. Obtain stress-strain curve data from triaxial compression and direct shear tests of the slope rock mass, and generate the test dataset required for model parameter calibration.

[0127] Specifically, the triaxial compression test was conducted according to the rock mechanics test specifications. During the test, the confining pressure was kept constant, and axial loads were applied in stages, with axial stress, axial strain, and lateral strain data recorded simultaneously until the rock mass failed. The test was conducted under different confining pressure conditions to obtain the rock mass mechanical response under different stress states. At least three parallel tests were performed under each confining pressure condition to ensure data reliability. The axial stress-axial strain curves obtained from the test contain information about the entire process of the rock mass from elastic deformation and plastic yielding to damage and failure, and are the core data source for calibrating the parameters of the damage evolution model.

[0128] Direct shear tests were conducted using a rock mass structural plane direct shear apparatus. Rock mass specimens containing natural structural planes were selected. By applying normal stress and shear load, the relationship curves between shear stress and shear displacement were recorded, and the volume change during the shearing process was monitored. Direct shear tests were carried out under different normal stresses to obtain the shear strength parameters and shear deformation characteristics of the structural planes. The test data can supplement the damage evolution information of the rock mass under shear load, complementing the triaxial compression test data and ensuring that the test dataset can cover the mechanical behavior of the rock mass under different stress modes.

[0129] The raw data obtained from the experiment were preprocessed. First, outlier data points were removed, and experimental noise was eliminated through data smoothing. Then, the stress-strain data was converted into a format that the model could use, including stress tensor components, strain tensor components, plastic strain components, and corresponding time series data. The experimental data under different confining pressures and normal stresses were integrated to form an experimental dataset containing multiple sets of stress-strain curves. Each set of data needed to clearly define the corresponding initial stress state, experimental loading path, and final failure state, providing comprehensive sample information for the subsequent training of the backpropagation neural network.

[0130] S23. Employing a backpropagation neural network algorithm, using stress-strain data from the experimental dataset as training samples and the damage evolution equation as physical constraints, the parameters of the material to be calibrated are iteratively optimized by minimizing the error between the model predictions and experimental measurements. and The value of is used to obtain the calibrated damage evolution prediction model.

[0131] Specifically, the backpropagation neural network algorithm consists of an input layer, hidden layers, and an output layer. The number of neurons in the input layer is determined based on the feature dimensions of the training samples, and strain components and time variables from the experimental data are selected as input features. Two to three hidden layers are set, and the number of neurons in each layer is determined through trial and error to ensure that the network has sufficient fitting ability and avoids overfitting. The neurons in the output layer output the stress components predicted by the model, which are compared with the stress components measured in the experiment.

[0132] The stress-strain data in the experimental dataset were divided into training and validation sets in a 7:3 ratio. The training set was used for iterative updates of the network weights and biases, while the validation set was used to monitor the network's generalization ability and prevent overfitting. During network training, the damage evolution equation was used as a physical constraint, and the parameters in the damage evolution equation were... , As an optimizable parameter of the network, it is embedded into the network's computation process to ensure that the network prediction process meets the physical laws of rock mass damage evolution, rather than simply mathematical fitting.

[0133] The error between model predictions and experimental measurements is evaluated using the mean squared error (MSE). The formula for calculating the MSE is: Where n is the number of data points, The stress value predicted by the model for the i-th data point. This represents the stress value at the i-th data point measured in the experiment. The error is calculated using the backpropagation algorithm based on the network parameters (weights, biases) and the parameters of the material to be calibrated. , The gradient of the validation set is used to iteratively update the parameter values ​​using the gradient descent method. After each iteration, the mean square error of the validation set is calculated. The iteration stops when the mean square error of the validation set remains stable for multiple consecutive iterations and is less than a set threshold.

[0134] After iterative optimization, extract the final result. , The values ​​are selected and substituted into the damage evolution equation. Combined with other related equations of the anisotropic damage constitutive model (such as the effective stress calculation formula and the plastic strain rate calculation formula), a complete damage evolution prediction model is formed. The model is then validated by inputting the stress-strain data of the validation set into the model and calculating the relative error between the model's predicted values ​​and the measured values. If the relative error meets the engineering accuracy requirements, it indicates that the model parameter calibration is effective, and the obtained damage evolution prediction model can be used for damage evolution simulation in the subsequent slope excavation process.

[0135] refer to Figure 2 In one optional embodiment, based on a damage evolution prediction model, the spatiotemporal evolution trajectory of the potential slip surface is predicted by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during graded excavation, including the following steps:

[0136] S31. Based on the three-dimensional geological model of the slope, the excavation and unloading process is simulated through the element birth and death technique to generate a graded excavation numerical model.

[0137] Specifically, the core of the element birth and death technique is to correct the simulated mechanical response after the rock mass is removed by using the stiffness matrix, as shown in the formula: .in The stiffness matrix of the element after "killing" it represents the stiffness state after the element's mechanical contribution is reduced. This is the initial stiffness matrix of the element, reflecting the original stiffness characteristics of the element before excavation. This is the stiffness reduction factor, used to control the degree of stiffness reduction in order to match the rock mass removal effect.

[0138] The time interval between excavation steps is matched with the construction rhythm, and the calculation formula is as follows: .in The time interval between adjacent excavation steps; h is the single-stage excavation height; This represents the average daily excavation speed.

[0139] The initial stress field equilibrium must satisfy the convergence criterion. .in These are the unbalanced forces in the model, reflecting the degree of imbalance in the model's mechanical state. The total external forces on the model include gravity, ground stress, etc. The convergence threshold is used to determine whether the stress field has reached equilibrium. The stress is updated after each excavation stage, using the following formula: .in This refers to the stress that is renewed after excavation; The initial stress before excavation; This refers to the stress increment caused by excavation and unloading.

[0140] S32. In the process of solving the numerical model of graded excavation, the changes in stress field and damage field at each excavation time step are calculated by the damage evolution prediction model, and a spatiotemporal evolution dataset containing stress variables, strain variables and damage variables is generated.

[0141] Specifically, the core of the calculation is the coupled calculation of stress and damage, and the coupling formula is as follows: .in [D] represents the stress after damage correction; [D] is the elastic matrix, reflecting the elastic deformation characteristics of the rock mass. The total strain includes both elastic and plastic strain. D represents plastic strain, reflecting the degree of plastic deformation of the rock mass; D is the damage variable, characterizing the degree of damage to the rock mass.

[0142] The damage evolution process is described by the damage evolution equation, the formula of which is: .in is the increment of the damage variable; Y is the damage energy release rate, reflecting the driving effect of stress state on damage development; and These are the parameters of the material to be calibrated, which are related to the damage characteristics of the rock mass itself. t represents the plastic strain rate, reflecting the rate of change of plastic strain; t represents time.

[0143] The formula for calculating the damage energy release rate Y is as follows: .in This is the transpose matrix of the total strain; This is the derivative of the effect of damage on mechanical properties.

[0144] Plastic strain rate Calculated using the plastic flow rule, the formula is as follows: .in is the plasticity multiplier, which controls the development range of plastic strain; f is the yield function, used to determine whether the rock mass has entered the plastic state; The partial derivative of the yield function with respect to pressure reflects the direction of plastic strain development.

[0145] S33. Based on the strain field data in the spatiotemporal evolution dataset, the maximum shear strain rate criterion is used to identify the set of units where the difference between the shear strain rate and the background value is greater than the first threshold, and a spatial distribution map of the damage localization region is generated.

[0146] Specifically, calculating the maximum shear strain rate is the core of identifying damage localization, and the formula is as follows: .in The maximum shear strain rate reflects the severity of the element's shear deformation. The maximum principal strain rate characterizes the rate of strain change of the element in the most dominant deformation direction; The minimum principal strain rate characterizes the strain rate of an element in the secondary deformation direction.

[0147] The background shear strain rate was obtained through statistical calculation, and the formula is as follows: .in The background shear strain rate reflects the overall normal deformation level of the slope; N is the total number of elements after excluding abnormal elements. is the maximum shear strain rate of the i-th element.

[0148] The formula for setting the first threshold is: .in is the first threshold used to filter elements with abnormal shear strain rate; k is a coefficient related to rock mass stability requirements; The standard deviation of the maximum shear strain rate is calculated using the following formula: This reflects the degree of dispersion of the maximum shear strain rate.

[0149] S34. Perform cluster analysis on the spatial distribution map, and aggregate the damage localization units that are spatially continuous and have strain characteristic differences less than the second threshold into the same potential slip surface, generating a set of spatial locations of potential slip surfaces and their development trends at each excavation time step.

[0150] Specifically, spatial continuity is determined based on the Euclidean distance formula. .in The spatial distance between unit i and unit j; , , The center coordinates of element i; , , Let be the center coordinates of element j.

[0151] The formula for calculating the difference in strain characteristics is as follows: .in The difference in strain characteristics between the k-th unit within the cluster and the cluster mean; The maximum shear strain rate of the k-th element within the cluster; The mean of the maximum shear strain rate within the cluster is calculated using the following formula: M is the number of units within the cluster.

[0152] The formula for setting the second threshold is as follows: .in The second threshold is used to determine the consistency of strain characteristics within a cluster; These are the discrete coefficients; The standard deviation of the maximum shear strain rate within the cluster reflects the degree of dispersion of the strain characteristics within the cluster.

[0153] The formula for calculating the center coordinates of the potential slip surface is as follows: .in , , The center coordinates of the potential slip surface; , , The coordinates are the center coordinates of the k-th unit within the cluster.

[0154] S35. Based on the spatial location set of potential slip surfaces and their development trend, calculate the penetration rate of each potential slip surface at different time steps using the vector integration method, and calculate the expansion rate of each potential slip surface by the change in the geometric shape of the slip surface between time steps.

[0155] Specifically, the breakthrough rate is calculated using the vector integral method, with the core formula being: Where P is the penetration rate of the potential slip surface, characterizing the integrity of the slip surface; The total length of the through path is calculated using the following formula: , , , The coordinates of the endpoints of the free surface. , , The coordinates of the endpoints of the internally stable rock mass; The contribution coefficient for the breakthrough is correlated with the damage variable; For a small vector segment of the path, the expression is: ; The magnitude of a small vector segment is calculated using the following formula: .

[0156] The formula for calculating the scalar expansion rate is as follows: .in The velocity of the slip surface along its length; Let t be the length of the sliding surface at time step t; The length of the sliding surface at time step t+1; This represents the time interval between adjacent time steps.

[0157] The formula for calculating the vector spread velocity is: .in This represents the migration velocity vector of the feature points on the slip surface. , , Here are the coordinates of the feature point at time step t; , , Here are the coordinates of the feature point at time step t+1. The formula for calculating the magnitude of the vector velocity is: , , , These are the components of the velocity vector in the x, y, and z directions, respectively.

[0158] S36. Based on the penetration rate and expansion speed of all potential slip surfaces in all excavation time steps, construct the spatiotemporal evolution trajectory of potential slip surfaces from initiation to penetration through the entire process using a spatiotemporal interpolation algorithm.

[0159] Specifically, spatial interpolation uses the Kriging interpolation method, with the core formula being: .in Points to be interpolated The damage or strain value; n is the number of known data points; The weight coefficient for the i-th known data point; This represents the damage or strain value at the i-th known data point. The weighting coefficients are obtained by solving the Kriging equations. Confirmed, among which For known data points and The variability function between them For known data points Interpolation point The variability function between them It is a Lagrange multiplier.

[0160] Time interpolation employs cubic spline interpolation to construct a continuous function for the throughput. ( Where P(t) is a continuous function of the penetration rate over time; t is time. , , , These are the parameters of the interpolation function; This refers to the moment at the i-th time step; This represents the time step at time i+1. The parameters are obtained by solving a system of equations. Confirmed, among which for The throughput rate of the time step; for The throughput rate of the time step.

[0161] The integrated formula for the spatiotemporal evolution trajectory is: .in The spatiotemporal evolution trajectory of the potential slip surface; is the damage or strain value after spatiotemporal interpolation; P(t) is the penetration rate after time interpolation; v(t) is the function of the expansion rate as a function of time.

[0162] In one optional embodiment, based on the stress state and damage degree in the spatiotemporal evolution trajectory, a spatiotemporal support sequence scheme is generated by determining the support timing and calculating the support location priority, including the following steps:

[0163] S41. Based on the damage variables and plastic zone penetration rate of the potential slip surface in its spatiotemporal evolution trajectory, the timing of support is determined using a dual-threshold judgment rule; wherein, the criterion for determining the timing of support is:

[0164] When damage variables or plastic zone penetration rate Immediate support should be provided when necessary. and Phased support, when and Temporarily no support;

[0165] in, The first critical damage threshold, The second critical damage threshold is determined based on the rock mass type at the corresponding location in the rock mass parameter field. The first breakthrough rate threshold, This is the second breakthrough rate threshold, and .

[0166] Specifically, the dual-threshold judgment method achieves precise determination of support timing through the coordinated judgment of damage variables and the penetration rate of the plastic zone. The damage variable D characterizes the degree of damage to the rock mass in the potential slip surface area, and its value ranges from 0 to 1. The larger the value, the more severe the rock mass damage. The first critical damage threshold, The second critical damage threshold is determined based on the rock mass type at the location of the potential slip surface. Different rock mass types have different critical damage tolerance values ​​due to differences in mechanical properties.

[0167] The plastic zone penetration rate β characterizes the degree of penetration of the plastic deformation region in the potential slip surface area. Its value ranges from 0 to 1. The larger the value, the closer the plastic zone is to complete penetration. The first breakthrough rate threshold, It is the second breakthrough rate threshold, and satisfies The relationship ensures that the two thresholds form a reasonable gradient division within the range of the pass rate, avoiding overlap or gaps in the judgment intervals.

[0168] Immediate support is suitable for scenarios where rock mass damage or plastic zone penetration is close to a critical state, requiring rapid intervention to prevent further development of the slip surface; phased support is suitable for scenarios where the rock mass is moderately damaged or the plastic zone is partially penetrated, and support can be implemented in stages according to the construction progress; temporary non-support is suitable for scenarios where the rock mass damage is minor and the plastic zone has not developed significantly, where monitoring can continue and support resources can be temporarily withheld, achieving a precise match between the timing of support and the state of the rock mass.

[0169] S42. Based on the difference between the current stress and the initial stress of the potential slip surface in its spatiotemporal evolution trajectory, and combining damage variables and displacement increments, the support location priority is calculated by weighted summation; the formula for calculating the support location priority is:

[0170]

[0171] in, σ represents the priority of support locations, used to quantify the urgency of support required at each location in the 3D geological model of the slope; σ is the current stress. The initial stress is extracted from the rock mass parameter field; ΔL is the displacement increment, obtained by comparing the current displacement of the potential slip surface with the initial displacement. The reference displacement is determined based on the geometric dimensions of the three-dimensional geological model of the slope; , , This is a weighting coefficient, determined based on the project's safety level.

[0172] Specifically, the priority of support location is the core indicator for quantifying the urgency of support at various locations on a slope; the higher the value, the more urgently support is needed at that location. The formula includes three key calculation terms, reflecting support requirements from three dimensions: stress change, damage state, and displacement development. A comprehensive assessment of multiple factors is achieved through weighted summation.

[0173] The first calculation term is the stress difference coefficient, which is composed of the current stress σ and the initial stress. and weighting coefficients Composition. The current stress σ is the real-time stress value of the corresponding position of the potential slip surface in the spatiotemporal evolution trajectory, and the initial stress. Extracted from the rock mass parameter field, it represents the original stress value at that location when it was not affected by excavation. This represents the absolute difference between the current stress and the initial stress, divided by the initial stress. Then the relative degree of stress difference is obtained, avoiding judgment bias caused by the magnitude of the absolute value of stress; This is the weighting coefficient for this calculation item, used to adjust the importance of stress variation factors in priority assessment.

[0174] The second calculation term is the damage contribution term, which consists of the damage variable D and the weighting coefficients. Composition. The definition of damage variable D is consistent with that in S41, directly reflecting the degree of damage to the rock mass; This is the weighting coefficient for this calculation item, used to adjust the importance of damage status factors in priority assessment.

[0175] The third calculation item is the displacement increment coefficient, which is composed of the displacement increment ΔL and the reference displacement. and weighting coefficients Composition. The displacement increment ΔL is obtained by comparing the current displacement with the initial displacement at the corresponding location of the potential slip surface, reflecting the displacement change at that location during the excavation process; the reference displacement... The geometric dimensions of the slope are determined based on the three-dimensional geological model, ensuring the relativization of displacement increments. This is the weighting coefficient for this calculation item, used to adjust the importance of displacement development factors in priority assessment.

[0176] , , The sum of the three is 1. The specific value is determined according to the safety level of the project. The higher the safety level, the greater the weight coefficient of key factors such as stress change and damage state, so as to ensure that higher risk locations are given priority in high safety requirements.

[0177] S43. Based on the timing and location priority of support, a topological sorting algorithm is used to generate a spatiotemporal sequence of support operations, resulting in a support spatiotemporal sequence scheme.

[0178] Specifically, the core of the topological sorting algorithm is to determine a reasonable execution order of support operations based on the sequential constraints between them, while also considering the time requirements of support timing to form a three-dimensional spatiotemporal sequence of "time-location-operation order". First, a directed acyclic graph of support operations is constructed. The nodes in the graph represent specific support locations (divided by priority), and the directed edges represent the sequential constraints between operations. Nodes with higher priority support locations point to nodes with lower priority, indicating that support for higher priority locations must be completed first, followed by support for lower priority locations.

[0179] Based on the directed acyclic graph (DAG), a time window is assigned to each node according to the time requirements of support. For nodes requiring immediate support, the time window is set to a short-term interval starting from the current moment; for nodes requiring phased support, the time window is divided into multiple consecutive intervals according to the phased requirements; for nodes not requiring support immediately, the time window is set to a long-term interval. The topological sorting algorithm traverses the DAG and extracts nodes sequentially according to the principle of "earlier time window, higher priority," forming the execution order of support operations.

[0180] The generated spatiotemporal sequence plan for support must specify the specific operation time (time window), operation sequence (position in the overall sequence), and corresponding support type recommendations (based on priority and rock mass condition matching) for each support location. The plan must ensure that support operations at high-priority locations are executed first within their corresponding time windows to avoid reduced support effectiveness due to improper sequence or time delays. At the same time, it should achieve reasonable allocation of support resources and improve the efficiency and safety of support projects.

[0181] The aforementioned method for generating a construction scheme for rock slope support model constructs a three-dimensional geological model by integrating geological survey data and spatial interpolation algorithms. It then inverts the initial damage field based on microseismic monitoring data and machine learning methods. Subsequently, a constitutive model considering the plastic coupling effect of damage is established, and its parameters are calibrated using experimental data. This forms a numerical model capable of accurately predicting the damage evolution path during excavation. Based on this, the development trajectory of potential slip surfaces is dynamically identified. A graded support strategy is formulated based on the spatiotemporal variation of stress state and damage degree. Finally, a multi-objective optimization algorithm balances economy and construction efficiency while ensuring safety, ultimately generating a forward-looking and adaptive dynamic support scheme. This represents a fundamental shift from traditional static experience-based design to dynamic and precise control based on damage evolution mechanisms. It significantly improves the scientific nature of support timing, the accuracy of support location, and the economic rationality of support parameters. It effectively overcomes the resource waste caused by premature support and the control failure caused by late support, thereby fundamentally improving the stability and controllability, construction safety, and life-cycle economic benefits of slope engineering.

[0182] It should be understood that although the steps in the flowcharts of the embodiments described above are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages of other steps.

[0183] Based on the same inventive concept, this application also provides an apparatus for implementing the above-mentioned method for generating construction schemes for rock slope support models. The solution provided by this apparatus is similar to the solution described in the above method. Therefore, the specific limitations of one or more embodiments of the rock slope support model construction scheme generation apparatus provided below can be found in the limitations of the rock slope support model construction scheme generation method described above, and will not be repeated here.

[0184] In one exemplary embodiment, such as Figure 3 As shown, a rock slope support model construction scheme generation device 30 is provided to implement the methods in the above-described method embodiments. The device includes:

[0185] The geological model construction and parameter allocation module 31 is used to construct a three-dimensional geological model of the slope based on geological exploration data. It combines the Kriging space interpolation algorithm to allocate discrete rock mass mechanical parameters to each unit grid of the three-dimensional geological model of the slope, thereby generating a rock mass parameter field.

[0186] The damage state dynamic inversion module 32 is used to collect rock mass microfracture data using a pre-deployed microseismic monitoring network, and based on the spatial distribution characteristics of the rock mass microfracture data and rock mass parameter field, to invert the initial damage state of the rock mass through a statistical learning algorithm and generate an initial damage field.

[0187] The damage constitutive model calibration module 33 is used to establish an anisotropic damage constitutive model that considers the coupling between damage and plasticity based on the damage distribution characteristics of the initial damage field; and to perform back-analysis calibration of the model parameters of the anisotropic damage constitutive model using experimental data to generate a damage evolution prediction model with calibrated parameters.

[0188] The slope stability prediction module 34 is used to predict the spatiotemporal evolution trajectory of potential slip surfaces by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during the graded excavation process, based on the damage evolution prediction model.

[0189] The intelligent decision-making module 35 for support timing is used to generate a spatiotemporal sequence scheme for support by determining the support timing and calculating the support location priority based on the stress state and damage degree in the spatiotemporal evolution trajectory.

[0190] The multi-objective optimization decision module 36 is used to establish a multi-objective optimization model based on the support parameter constraints of the support spatiotemporal sequence scheme with safety, economy and construction period as multiple objectives; and to form a multi-objective optimized support scheme by solving the optimal support parameters of the multi-objective optimization model.

[0191] The intelligent integration module 37 for construction schemes is used to integrate spatiotemporal sequence schemes for support and multi-objective optimized support schemes, and output construction schemes.

[0192] Embodiments of this application also provide a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the aforementioned method embodiments.

[0193] Embodiments of this application also provide a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps in the above-described method embodiments.

[0194] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to in the description of the method embodiments. The device embodiments described above are merely illustrative. The components described as separate parts may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this disclosure according to actual needs. Those skilled in the art can understand and implement this without creative effort.

[0195] The above-described embodiments are merely illustrative of several implementation methods of the embodiments of this application, and their descriptions are relatively specific and detailed. However, they should not be construed as limiting the scope of the patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the embodiments of this application, and these modifications and improvements all fall within the protection scope of the embodiments of this application.

Claims

1. A method for generating a construction scheme for a rock slope support model, characterized in that, The method includes: S1. Construct a three-dimensional geological model of the slope based on geological survey data, and combine the Kriging space interpolation algorithm to allocate discrete rock mass mechanical parameters to each unit grid of the three-dimensional geological model of the slope to generate a rock mass parameter field. S2. Collect rock mass microfracture data using a pre-deployed microseismic monitoring network. Based on the spatial distribution characteristics of the rock mass microfracture data and the rock mass parameter field, invert the initial damage state of the rock mass using a statistical learning algorithm to generate an initial damage field. S3. Based on the damage distribution characteristics of the initial damage field, establish an anisotropic damage constitutive model that considers the coupling between damage and plasticity; perform back-analysis calibration on the model parameters of the anisotropic damage constitutive model using experimental data to generate a damage evolution prediction model with calibrated parameters. S4. Based on the damage evolution prediction model, the spatiotemporal evolution trajectory of the potential slip surface is predicted by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during the graded excavation process. S5. Based on the stress state and damage degree in the spatiotemporal evolution trajectory, generate a spatiotemporal support sequence scheme by determining the support timing and calculating the support location priority. S6. Taking safety, economy, and construction period as multiple objectives, and based on the support parameter constraints of the support spatiotemporal sequence scheme, establish a multi-objective optimization model; by solving the optimal support parameters of the multi-objective optimization model, a multi-objective optimized support scheme is formed. S7. Integrate the spatiotemporal sequence scheme of the support and the multi-objective optimized support scheme, and output the construction scheme.

2. The method according to claim 1, characterized in that, The process involves collecting rock mass microfracture data using a pre-deployed microseismic monitoring network, and based on the spatial distribution characteristics of the rock mass microfracture data and the rock mass parameter field, inverting the initial damage state of the rock mass using a statistical learning algorithm to generate an initial damage field, including: S11. Through the microseismic monitoring network, monitor the spatiotemporal coordinates, energy release, and magnitude parameters of rock microfractures during the excavation process, and generate a raw dataset of microseismic events. S12. Based on the original dataset of the microseismic events and the rock mass type partitioning of the rock mass parameter field, the spatial clustering characteristics of the microseismic events are analyzed by geostatistical variability function to generate regional division schemes for different damage levels. S13. Map the region division scheme onto the grid cells of the three-dimensional geological model of the slope, and calculate the initial damage variable for each grid cell based on the number and energy amplitude of microseismic events falling into each grid cell; the formula for calculating the initial damage variable is: in, The initial damage variable of the grid cell characterizes the initial damage degree of the rock mass within the grid cell; N is the number of microseismic events within the grid cell, obtained by statistically analyzing the microseismic events falling into the grid cell in the original dataset of microseismic events. The number of failure threshold events for the grid cell is determined based on the rock mass type at the corresponding location of the grid cell in the rock mass parameter field; The energy amplitude of the microseismic event is obtained from the original dataset of the microseismic event. The reference energy amplitude is determined based on the statistical characteristics of the energy amplitudes of all microseismic events in the original microseismic event dataset. S14. Based on the initial damage variables of all the grid cells, a spatially continuous damage distribution field is generated by a kernel density estimation algorithm to obtain the initial damage field.

3. The method according to claim 1, characterized in that, Based on the damage distribution characteristics of the initial damage field, an anisotropic damage constitutive model considering damage and plastic coupling is established; the model parameters of the anisotropic damage constitutive model are calibrated through back-analysis using experimental data to generate a parameter-calibrated damage evolution prediction model, including: S21. Based on the anisotropic damage tensor theory and the aforementioned damage distribution characteristics, a damage evolution equation incorporating the parameters of the material to be calibrated is established by defining a functional relationship between the damage energy release rate and the plastic strain rate; the damage evolution equation is as follows: Where D is the damage variable, representing the current degree of damage to the rock mass; Y is the damage energy release rate, calculated based on the current stress state and damage state. and The parameters of the material to be calibrated; t represents the plastic strain rate, calculated using the plastic flow law; t represents time. S22. Obtain stress-strain curve data from triaxial compression and direct shear tests of the slope rock mass, and generate the test dataset required for model parameter calibration; S23. Using a backpropagation neural network algorithm, with the stress-strain data in the experimental dataset as training samples and the damage evolution equation as physical constraints, the parameters of the material to be calibrated are iteratively optimized by minimizing the error between the model predictions and the experimental measurements. and The value of is used to obtain the calibrated damage evolution prediction model.

4. The method according to claim 1, characterized in that, The method based on the damage evolution prediction model, by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during graded excavation, predicts the spatiotemporal evolution trajectory of the potential slip surface, including: S31. Based on the three-dimensional geological model of the slope, the excavation and unloading process is simulated by the element birth and death technique to generate a graded excavation numerical model. S32. In the process of solving the numerical model of graded excavation, the stress field and damage field changes at each excavation time step are calculated by the damage evolution prediction model to generate a spatiotemporal evolution dataset containing stress variables, strain variables and damage variables. S33. Based on the strain field data in the spatiotemporal evolution dataset, the maximum shear strain rate criterion is used to identify the set of units where the difference between the shear strain rate and the background value is greater than a first threshold, and a spatial distribution map of the damage localization region is generated. S34. Perform cluster analysis on the spatial distribution map, and aggregate the damage localization units that are spatially continuous and have strain characteristic differences less than the second threshold into the same potential slip surface, generating a set of spatial locations of potential slip surfaces and their development trends at each excavation time step. S35. Based on the set of spatial locations of the potential slip surfaces and their development trends, calculate the penetration rate of each potential slip surface at different time steps using the vector integration method, and calculate the expansion rate of each potential slip surface by the change in the geometric shape of the slip surface between consecutive time steps. S36. Based on the penetration rate and expansion speed of all potential slip surfaces at all excavation time steps, construct the spatiotemporal evolution trajectory of the potential slip surfaces from initiation to penetration through a spatiotemporal interpolation algorithm.

5. The method according to any one of claims 1 to 4, characterized in that, The step of generating a spatiotemporal support sequence scheme based on the stress state and damage degree in the spatiotemporal evolution trajectory by determining the support timing and calculating the support location priority includes: S41. Based on the damage variables and plastic zone penetration rate of the potential slip surface in the spatiotemporal evolution trajectory, the support timing is determined using a dual-threshold judgment rule; wherein, the criterion for determining the support timing is: When the damage variable Or the penetration rate of the plastic zone Immediate support should be provided when necessary. and Phased support, when and Temporarily no support; in, The first critical damage threshold, The second damage critical threshold is determined based on the rock mass type at the corresponding location in the rock mass parameter field; The first breakthrough rate threshold, This is the second breakthrough rate threshold, and ; S42. Based on the degree of difference between the current stress and the initial stress of the potential slip surface in the spatiotemporal evolution trajectory, and in conjunction with the damage variable and displacement increment, the support location priority is calculated by weighted summation; the formula for calculating the support location priority is: in, The support location priority is used to quantify the urgency of support required at each location in the three-dimensional geological model of the slope; σ is the current stress. The initial stress is extracted from the rock mass parameter field; ΔL is the displacement increment, obtained by comparing the current displacement of the potential slip surface with the initial displacement. The reference displacement is determined based on the geometric dimensions of the three-dimensional geological model of the slope. , , The weighting coefficient is determined based on the project's safety level. S43. Based on the support timing and the support location priority, a topological sorting algorithm is used to generate a spatiotemporal sequence of support operations, thereby obtaining the support spatiotemporal sequence scheme.

6. A device for generating construction schemes for rock slope support models, used to implement the method described in any one of claims 1 to 5, characterized in that, The device includes: The geological model construction and parameter allocation module is used to construct a three-dimensional geological model of the slope based on geological exploration data. It combines the Kriging space interpolation algorithm to allocate discrete rock mass mechanical parameters to each cell grid of the three-dimensional geological model of the slope, thereby generating a rock mass parameter field. The damage state dynamic inversion module is used to collect rock mass microfracture data using a pre-deployed microseismic monitoring network, and based on the spatial distribution characteristics of the rock mass microfracture data and the rock mass parameter field, invert the initial damage state of the rock mass through a statistical learning algorithm to generate an initial damage field. The damage constitutive model calibration module is used to establish an anisotropic damage constitutive model that considers the coupling between damage and plasticity based on the damage distribution characteristics of the initial damage field; and to perform back-analysis calibration on the model parameters of the anisotropic damage constitutive model using experimental data to generate a damage evolution prediction model with calibrated parameters. The slope stability prediction module is used to predict the spatiotemporal evolution trajectory of potential slip surfaces by simulating the stress redistribution and damage evolution process of the three-dimensional geological model of the slope during the graded excavation process, based on the damage evolution prediction model. The intelligent decision-making module for support timing is used to generate a spatiotemporal sequence scheme for support by determining the support timing and calculating the support position priority based on the stress state and damage degree in the spatiotemporal evolution trajectory. The multi-objective optimization decision module is used to establish a multi-objective optimization model based on the support parameter constraints of the spatiotemporal sequence of the support scheme, with safety, economy and construction period as multiple objectives; and to form a multi-objective optimized support scheme by solving the optimal support parameters of the multi-objective optimization model. The intelligent integration module for construction schemes is used to integrate the spatiotemporal sequence scheme of the support and the multi-objective optimized support scheme, and output the construction scheme.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 5.

Images