A high and steep slope dangerous rock mass instability movement risk evaluation method and application thereof

Abstract

The application relates to the technical field of geotechnical engineering and geological disaster prevention, in particular to a high and steep slope dangerous rock mass instability movement risk evaluation method and application thereof, which comprises the following steps: constructing a refined digital terrain model by using an unmanned aerial vehicle (UAV) oblique photography and cloth simulation filtering algorithm; performing probabilistic movement simulation based on a physical engine, and extracting characteristic parameters such as impact kinetic energy and bounce height; introducing a Wasserstein distance as a target function, and combining a grid self-adaptive direct search algorithm to iteratively optimize physical parameters; establishing a multi-field coupling correction model containing precipitation, underground water and seismic load, and dynamically correcting the movement process; and finally, constructing a harmfulness grade matrix through cross coupling of danger classification and harm object grade. Through quantitative simulation and multi-factor coupling, the application provides a quantitative index for risk evaluation under the condition of no historical rockfall data, and provides physical parameter support for engineering protection design.

Description

Technical Field

[0001] This invention belongs to the fields of geotechnical engineering, geological disaster prevention and control and numerical simulation technology. Specifically, it relates to a method for risk assessment of unstable rock mass movement on steep slopes and its application. Background Technology

[0002] Steep slope rock mass hazards pose a significant challenge to geotechnical engineering and geological disaster prevention. Their instability and collapse are highly sudden and destructive, severely threatening the safety of hydraulic structures, resettlement areas, transportation facilities, and construction camps. To effectively address this risk, existing research typically combines rock mass stability assessment with risk evaluation, quantifying instability hazard by constructing risk assessment models and indicator systems. These methods can establish assessment criteria and determine risk levels based on the main factors inducing collapse. However, in practical applications, these methods primarily rely on static ratings using statistical indicator systems and cannot provide detailed simulations of the dynamic movement process after rock mass instability.

[0003] Chinese patents CN119251714A and CN107194049A establish a risk assessment index system for rockfall based on the main factors inducing slope collapse and establish a classification standard for rockfall risk assessment, determining the risk assessment level. However, the above methods are mainly based on the risk rating of the statistical index system and lack detailed dynamic simulation of the instability process. In terms of rockfall movement simulation, more and more studies emphasize the necessity of fine reconstruction of three-dimensional terrain for simulating the movement trajectory of rockfall, and carry out dynamic simulation and risk analysis on this basis. For example, patents CN114964241A and CN119991994A first obtain slope point cloud data, reconstruct a real-world three-dimensional model of the slope, extract slope elevation points, plot the slope elevation points, establish a three-dimensional slope model and divide it into grids, and finally build a model in PFC for calculation. However, the above methods still rely on simplified rockfall particles.

[0004] These simulation methods still largely rely on simplified rockfall particle models and primarily focus on model inversion and parameter correction for rockfall events that have already occurred. Their core research objectives are often limited to improving the accuracy of trajectory prediction, failing to effectively address the risk prediction of potential unstable rock masses under conditions where instability has not yet occurred, and also lacking a risk assessment index system based on the dynamic characteristics of the movement process. Summary of the Invention

[0005] Existing technologies primarily focus on parameter calibration or predictive capability assessment of rockfall movement models, emphasizing inversion correction of models using past rockfall events. Their research objective is to improve the accuracy of trajectory prediction, but they do not address the risk prediction of potential unstable rock masses under conditions where instability has not yet occurred, nor do they construct a risk assessment index system based on the movement process. This invention, starting from engineering risk decision-making, constructs the spatial propagation characteristics of unstable rock mass movement through probabilistic simulation under conditions without historical rockfall events. Furthermore, it couples these characteristics with the threatened object to achieve quantitative evaluation of the risk of unstable rock mass movement and provide guidance for engineering protection.

[0006] To address the shortcomings of existing technologies in assessing the risk of unstable rock masses, such as relying heavily on static statistical index systems, lacking dynamic simulation of instability processes, excessively relying on historical rockfall event inversion for parameter acquisition, and failing to fully consider the coupling effects of multiple environmental factors, this invention provides a method for assessing the risk of unstable rock mass movement on steep slopes.

[0007] The technical solution adopted in this invention is as follows: This method achieves a quantitative assessment of the instability risk of potential unstable rock masses on steep slopes by constructing a refined digital terrain model, performing probabilistic motion simulation based on a physics engine, optimizing physical parameters using the Wasserstein distance and grid adaptive direct search algorithm, introducing a multi-field coupled dynamic correction mechanism, and constructing a hazard level matrix. Specific steps include: S1: Data Acquisition and Refined Construction of the Digital Terrain Model. Multi-angle imagery of the study area was acquired using UAV oblique photogrammetry, and an original digital terrain model was generated through 3D reconstruction technology. To eliminate the influence of surface vegetation and buildings on the true elevation of the terrain, a cloth simulation filtering algorithm was used to correct the original point cloud data. Specifically, the cloth simulation filtering algorithm is implemented as follows: the original point cloud is projected onto a horizontal plane, and a simulated cloth grid consisting of mass points and springs is constructed above the point cloud; the mass points fall under gravity until they contact ground points in the point cloud; the deformation of the cloth grid is controlled by constructing an objective function containing smoothness and friction terms, where the smoothness term constrains the elevation difference between adjacent mass points, and the friction term simulates the interaction between the cloth and the ground surface. By minimizing the objective function, the point cloud is automatically classified into ground points and non-ground points. After removing non-ground points, a high-precision digital terrain model with a resolution better than 2m×2m is output as the basic grid unit for subsequent simulation calculations.

[0008] S2: Three-dimensional motion simulation and physical parameter setting of unstable rock masses. In a three-dimensional physical simulation environment, based on the distribution characteristics of unstable rock masses determined by geological surveys, a geometric model of the unstable rock mass is constructed, and its physical properties are set. These physical properties include rock mass mass, friction coefficient, normal restitution coefficient, tangential restitution coefficient, and energy dissipation coefficient. Furthermore, the initial motion conditions of the unstable rock mass are set, including three-dimensional spatial coordinates, geometric shape, initial linear velocity, and initial angular velocity. The geometric shape adopts a random polyhedron model or an equivalent sphere model to simulate the motion characteristics under different collapse and disintegration states. Based on this, probabilistic motion simulation is performed. For each typical unstable rock mass, at least 1000 collapse simulation calculations are conducted, and the spatial distribution data of motion trajectory, bounce height, motion velocity, and impact kinetic energy are recorded and extracted in real time.

[0009] S3: Parameter Optimization Based on Wasserstein Distance and Grid Adaptive Direct Search Algorithm. For potential unstable rock masses lacking historical rockfall data, an automatic parameter calibration technique is employed to reduce the bias of human-based parameter assignment. Specifically, the Wasserstein distance is introduced as the objective function to quantitatively measure the difference between the simulated probability distribution of rockfall stopping points and the preset actual observation distribution or the potential deposition distribution determined by geological experts. The Wasserstein distance is used to assess simulation error by calculating the minimum transportation cost required to transform one distribution into another. Further, an iterative parameter optimization algorithm is used, which establishes a dynamically adjusted grid in the parameter space. Through alternating execution of the search phase and the extreme point phase, it automatically adjusts the friction coefficient, normal restitution coefficient, and tangential restitution coefficient without requiring gradient information from the objective function. When the Wasserstein distance reaches a preset convergence threshold, the optimal physical parameter set is output, ensuring a high degree of consistency between the simulation results and actual physical characteristics.

[0010] S4: Dynamic Introduction of Multi-Field Coupling Factors. To reflect the influence of complex environments on the movement process of unstable rock masses, this scheme establishes a multi-field coupling correction model for precipitation, groundwater, and seismic loads. The precipitation coupling mechanism introduces a precipitation sensitivity coefficient to correct the friction coefficient of the contact surface in real time; the correction formula reflects a nonlinear decrease in the friction coefficient with increasing precipitation. The groundwater coupling mechanism establishes a functional relationship between groundwater level change and the normal restitution coefficient to simulate the influence of pore water pressure on the energy dissipation of rock mass collisions. The seismic load coupling mechanism applies a time-varying external acceleration field to the physical model to calculate the additional inertial force of the unstable rock mass under seismic forces, thereby correcting the instability criteria and acceleration of the unstable rock mass.

[0011] S5: Quantitative Assessment of the Risk of Rock Mass Instability and Movement. Based on the dynamic movement characteristic parameters obtained from the above simulations, and combined with the stability state of the rock mass itself and the importance of the threatened objects, a full-chain risk assessment system is constructed. First, based on the stability calculation results of the rock mass (e.g., unstable, poorly stable, basically stable) and its scale level (e.g., extra-large, large, medium, small), the hazard index is divided into three levels: large, medium, and small. Second, based on the level and importance of the threatened hydraulic structures, transportation facilities, or densely populated areas, they are classified into three levels of hazard objects: Level I to Level III. Finally, by constructing a hazard level matrix, the hazard level and the hazard object level are cross-coupled, and the risk is ultimately divided into four levels: Level I (extremely high risk), Level II (high risk), Level III (medium risk), and Level IV (low risk). Based on this, a comprehensive hazard zoning map is generated, providing quantitative parameter support for the deployment of engineering protection measures.

[0012] Furthermore, in step S1, to improve the accuracy of terrain reconstruction, the trajectory planning for UAV oblique photography adopts a cross-overlap mode, ensuring that the forward overlap rate is not less than 80% and the lateral overlap rate is not less than 70%. During point cloud processing, the displacement step size of the cloth simulation filtering algorithm is set to a value that matches the grid resolution. Through multiple iterations until the displacement increment of the mass point is less than the preset tolerance, the generated digital terrain model can accurately restore micro-topographic features such as gullies, steep slopes, and slope turning points, providing accurate terrain constraints for the turning simulation of rockfall trajectories.

[0013] Furthermore, in step S2, the collision mechanics model of the unstable rock mass adopts a nonlinear damped spring model. At the instant the rock hits the slope, the normal velocity component and tangential velocity component after the collision are calculated using the normal restitution coefficient and the tangential restitution coefficient, respectively. The physics engine incorporates the influence of air resistance on long-distance motion in real time during the calculation process and considers the rotational kinetic energy conversion caused by the irregular shape of the rock mass, thereby obtaining a more physically realistic impact kinetic energy distribution curve through simulation.

[0014] Furthermore, in step S3, the grid size of the mesh adaptive direct search algorithm dynamically decreases during the iteration process. During the search phase, the algorithm randomly samples within the neighborhood of the current grid point. If a parameter combination with a smaller Wasserstein distance is found, the center point is updated and the grid is expanded; if no better solution is found, the algorithm enters the extreme point phase, performing a local search on a finer grid. Through this mechanism combining global and local searches, optimal approximation of physical parameters is achieved in a high-dimensional parameter space.

[0015] Furthermore, in step S4, the rainfall-based friction coefficient correction model considers changes in the saturation of the soil and rock mass. Specifically, differentiated rainfall correction weights are set for slopes with different geological compositions; for example, for slopes with a high clay content, the rainfall sensitivity coefficient is increased to reflect the slope lubrication effect caused by rainfall. For seismic loads, by inputting seismic wave sequences with specific frequencies and amplitudes, the dynamic tensile and shear stresses of the unstable rock mass during vibration are calculated, achieving a dynamic simulation of the entire earthquake-induced collapse process.

[0016] Furthermore, in step S5, the quantitative risk assessment system not only considers the probability of instability but also extracts the maximum bounce height, maximum impact kinetic energy, and maximum velocity from the simulated trajectory as key evaluation factors. Specifically, the maximum bounce height guides the design of the interception net's height, the maximum impact kinetic energy determines the structural strength requirements of the retaining wall, and the velocity distribution is used to define safe evacuation distances for personnel. By comparing these physical quantities with the disaster resistance level of the threatened object, the probability of disaster damage is calculated, thus achieving the transformation from qualitative description to quantitative assessment.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. Significantly improved the accuracy of terrain reconstruction and simulation. By employing a cloth-based simulation filtering algorithm to correct point cloud data and combining it with high-precision 2m×2m grid cells, interference from vegetation and other clutter on the terrain was eliminated, achieving accurate restoration of the micro-topographic features of steep slopes. This ensures that the simulated trajectory, turning characteristics, and stopping positions of unstable rock masses closely match the actual terrain constraints, providing a reliable geometric basis for subsequent risk assessment.

[0018] 2. This invention addresses the technical bottleneck of parameter determination for potential unstable rock masses when historical collapse events are lacking by introducing Wasserstein distance as the objective function for measuring error and combining it with a grid-adaptive direct search algorithm for automatic optimization. This method finds the optimal parameter combination through mathematical optimization, eliminating the subjectivity of human experience and thus enhancing the reliability of simulation results at the physical and mechanical levels.

[0019] 3. Enhanced adaptability of risk assessment to complex environments. By establishing a multi-field coupled dynamic correction mechanism for precipitation, groundwater, and seismic loads, this invention can simulate the instability and movement characteristics of unstable rock masses under different extreme weather and geological conditions. This dynamic correction mechanism makes the risk assessment results no longer a single static value, but a dynamic spectrum that changes with environmental factors, greatly improving the scientific rigor and relevance of early warning and forecasting.

[0020] 4. The risk assessment system constructed in this invention quantifies the entire process from rock mass instability and motion propagation to the response of the threatened object. By outputting specific engineering design parameters such as bounce height and impact kinetic energy, and deeply coupling them with the importance level of the threatened object, the generated risk zoning map can be directly used to guide the site selection and parameter design of protective projects such as interception nets and rockfall barriers, and has extremely strong engineering practical value.

[0021] 5. Improved efficiency in disaster prevention and control in complex geological environments. This invention, through automated data processing and probabilistic simulation algorithms, can quickly identify and classify the risks of unstable rock masses on large-scale, steep slopes. Compared to traditional on-site surveys and experience-based assessments, this method significantly reduces manpower and time costs while maintaining accuracy. Attached Figure Description

[0022] Figure 1 This is a schematic diagram of the UAV image acquisition steps proposed in this invention; Figure 2 This is a three-dimensional geological model of the study area proposed in this invention. Figure 3 This is a general distribution map of dangerous rocks in the study area according to an embodiment of the present invention; Figure 4 This invention is a simulation and assignment model for unstable rock mass and falling rock movement. Figure 5 This is a diagram showing the results of the three-dimensional numerical simulation proposed in this invention; Figure 6 This is the energy profile of a typical unstable rock mass collapse trajectory proposed in this invention; Figure 7 This is a schematic diagram of the energy change of a typical unstable rock mass collapse trajectory proposed in this invention; Figure 8 This is a typical velocity profile of a collapsed rock mass, as proposed in this invention. Figure 9 This is a schematic diagram of the velocity change of a typical unstable rock mass collapse trajectory proposed in this invention; Figure 10 This is a typical profile of the bounce height of a collapsed rock mass, as proposed in this invention. Figure 11 This is a schematic diagram illustrating the change in bounce height of a typical unstable rock mass collapse trajectory, as proposed in this invention. Figure 12 This is the comprehensive zoning map of the hazardous rock mass hazards in the study area proposed in this invention. Detailed Implementation

[0023] Example 1: This invention provides a method for assessing the risk of unstable rock mass movement on steep slopes. Through multi-source data fusion, refined terrain reconstruction, adaptive parameter optimization, and multi-field coupled numerical simulation, it achieves a quantitative risk assessment of potential collapse hazards on steep slopes. (See attached reference.) Figure 1 To be continued Figure 12 The specific implementation process of this invention is as follows.

[0024] Step 1: Data Acquisition and Digital Terrain Model (DTM) Construction: 1.1 High-resolution imagery data of the study area was acquired using UAVs, and 3D modeling was performed using Smart3D software to generate a high-precision digital terrain model (DTM). Each 2m×2m grid cell served as the basic unit for simulation calculations, ensuring high accuracy of the simulation results.

[0025] 1.2 The DTM is corrected based on actual terrain features (such as vegetation, soil layers, etc.). The CSF algorithm is used to process the point cloud and simulate the flow of cloth on the surface to distinguish between ground and vegetation points. The specific process of the DTM algorithm is as follows: (1) The input to the CSF algorithm is a set of coordinates of a three-dimensional point cloud. The coordinates of each point represent its position in space, typically as follows: ,in It is the index of the point. (2) For each point in the point cloud data First, the normal vector of the point is calculated through smoothing. Common smoothing methods include local neighborhood fitting or least squares: ; in, It is a point The set of neighborhood points, It represents the number of neighboring points. The normal vector reflects the local plane tilt of the point cloud data.

[0026] (3) In simulating the fabric, each point is imagined as a mass point on the fabric, and the surface of the fabric is formed by these mass points connected by springs. These springs simulate the elastic tension and friction of the fabric. The spring constant between each point (representing the tightness of the fabric) is a key parameter. According to the principles of physics, the fabric surface will deform under the action of external forces, and the friction of the fabric surface will vary according to the normal force at each point. The goal of the fabric simulation model is to minimize the bending and stretching of the fabric surface so that the fabric can conform to the ground.

[0027] (4) Based on the deformation of the fabric simulation model, the CSF algorithm classifies the point cloud in the following ways: Ground points: A point is classified as a ground point if its normal vector is parallel to the fabric surface and its location is close to the fabric surface.

[0028] Non-ground points: If the normal vector of a point has a large angle with the cloth surface, or if the point is far away from the cloth surface, then the point is classified as a non-ground point, usually vegetation or other obstacles.

[0029] (5) During fabric simulation, the CSF algorithm attempts to optimize the fabric shape by minimizing the objective function. The objective function typically consists of two parts: ① Smoothness term (height change of points within the smooth area): here, It is a point height, It is a set of ground points. It is a point The neighborhood set.

[0030] ② Friction term (simulating the friction of fabric): here, Point and The friction between them and These are their normal vectors.

[0031] (6) The CSF algorithm classifies the point cloud into ground points and non-ground points, and outputs a corrected DTM. Ground points will be used for subsequent high-precision 3D modeling and simulation calculations.

[0032] Specifically, the implementation process began with the use of drones to acquire high-resolution images of the study area. The drones were equipped with a five-lens oblique photography camera with a resolution of 20 megapixels and an integrated RTK positioning module, achieving a horizontal positioning accuracy of 1 cm ± 1 ppm and a vertical accuracy of 1.5 cm ± 1 ppm. During the aerial survey preparation phase, flight paths were planned based on the steep terrain features of the study area, ensuring a forward overlap of at least 80% and a lateral overlap of at least 70% to obtain complete three-dimensional surface features. The raw image data acquired from the aerial survey was imported into Smart3D software for aerial triangulation calculations, establishing a three-dimensional reality model with geographic coordinate reference, and outputting it in OSGB format. To obtain a digital terrain model that accurately reflects the actual surface undulations, vegetation removal and terrain correction were performed on the original point cloud.

[0033] This embodiment employs a cloth simulation filtering algorithm for ground point classification. The algorithm works by projecting the original point cloud onto a horizontal plane and constructing a virtual cloth mesh above the point cloud, consisting of mass points connected by springs. Each mass point falls under simulated gravity until it contacts a ground point in the point cloud. During the cloth's descent, the cloth shape is optimized by minimizing an objective function, which includes smoothness and friction terms.

[0034] The smoothness term constrains the height variation of adjacent mass points, ensuring terrain continuity. The friction term simulates the interaction force between the fabric and the ground by calculating the angle between the mass point's normal vector and the fabric surface. Through iterative calculation, when the fabric shape stabilizes and the mass point displacement increment is less than a preset tolerance, points in contact with the fabric are classified as ground points, and the remaining points are classified as non-ground points (such as vegetation or buildings). The final output is a high-precision digital terrain model with a resolution of 2m × 2m, providing fine geometric constraints for subsequent instability motion simulation, as shown in the attached figure. Figure 2 -Appendix Figure 4 As shown, the model accurately reproduces the micro-landforms such as gullies and steep slopes on the surface of high and steep slopes.

[0035] Step 2: In the simulation of unstable rock mass and rockfall movement on steep slopes, the three-dimensional rock mass movement simulation accurately predicts the instability process and movement trajectory of the rock mass under the action of external factors by constructing physical models and numerical methods.

[0036] 2.1 Set the physical parameters related to rock mass movement, including rock mass mass, friction coefficient, and coefficient of restitution. These parameters determine the behavior of the rock mass in the simulation. Among them, the friction coefficient... The coefficient of restitution is used to describe the frictional force between the rock mass and the slope surface, affecting the sliding and rolling behavior of the rock mass. and This describes the rebound behavior of a rock mass after colliding with a slope. Normal coefficient of restitution. and tangential restitution coefficient This determines the energy recovery efficiency during rock mass rebound.

[0037] Coefficient of Restoration Model: ; in, It is the velocity after the collision. It is the velocity before the collision. The coefficient of recovery.

[0038] Energy dissipation coefficient : Describes the energy loss of rock masses during movement, especially during collisions or friction.

[0039] 2.2 In RocPro3D, each simulated rock block needs initial conditions set, including position, velocity, mass, shape, etc. Based on the actual situation, initial parameters for multiple rock blocks are set to simulate their different motion trajectories.

[0040] The initial position of the rock mass is usually determined by field surveys or topographic data; it is set through field monitoring or by assuming an initial velocity. A common setting is based on historical data or experimental values ​​to set a standard initial velocity; and RocPro3D supports modeling rock masses of different shapes, including spheres, cuboids, or custom shapes. The choice of shape will affect the collision and bounce model.

[0041] Specifically, based on the constructed digital terrain model, the 3D physical simulation software RocPro3D was used to simulate the unstable movement of unstable rock masses. First, based on the results of the on-site geological survey, the distribution location, scale, and stability status of potential unstable rock masses were identified in the 3D model. A geometric model of the unstable rock mass was constructed in the simulation environment. This geometric model could be set to a sphere, cuboid, or a custom irregular polyhedron according to the actual shape of the collapsed blocks. For each unstable rock mass, its physical and mechanical parameters were set, including mass, friction coefficient, normal restitution coefficient, and tangential restitution coefficient. For example, for the dacite unstable rock mass in the study area, its normal restitution coefficient was set to 0.52, and its tangential restitution coefficient to 0.86; for the loose gravelly soil covering the slope, the normal restitution coefficient was set to 0.4, and its tangential restitution coefficient to 0.65. The restitution coefficient model determined the energy loss by calculating the ratio of velocities before and after the collision. Simultaneously, an energy dissipation coefficient was introduced to describe the nonlinear energy loss of the rock mass during complex movement. When setting initial conditions, the initial position and initial velocity of the unstable rock mass are given based on geomechanical analysis. The initial velocity is usually set to 0 m / s to simulate gravity-induced initial instability. Subsequently, probabilistic simulation calculations are performed, with no fewer than 1000 random path simulations for each unstable rock mass to cover uncertainties caused by minor terrain variations. During the simulation, the physics engine calculates in real time the sliding, rolling, and bouncing behavior of the unstable rock mass geometric model as it moves along the surface of a steep slope, as shown in the attached figure. Figure 4 As shown, a complete simulated motion trajectory is generated.

[0042] Step 3: To improve the accuracy and reliability of RocPro3D simulations, this approach optimizes key simulation parameters through algorithm adjustments, ensuring maximum matching with actual observation data. Key parameters include the normal recovery coefficient. Tangential coefficient of restitution and energy dissipation coefficient ; 3.1 Optimization Objective and Error Assessment: Using Wasserstein distance: Wasserstein distance is a way to measure the difference between two distributions. In simulations, it is used to compare the simulated distribution of rock block dwell points with the actual observed distribution of dwell points, calculating the minimum transport cost between the two distributions.

[0043] Wasserstein distance formula: in: This is the actual distribution of stop points in the data. It is the distribution of dwell points in the simulated data; and These are the actual and simulated locations of the rest stops, respectively. It is the transportation flow between simulated stop points and actual stop points. Is the actual distribution The Middle The location of each stop It is a simulated distribution The Middle The location of each stop.

[0044] The error between the simulated stop point and the actual stop point can be calculated using the Wasserstein distance, and this error can be used as the objective function for optimization.

[0045] 3.2 Minimize the Wasserstein distance as the objective function and iteratively update key parameters in the simulation (such as the normal restitution coefficient). Tangential coefficient of restitution and energy dissipation coefficient The specific steps are as follows: ① Select an initial parameter set This is obtained through experience or preliminary simulation data.

[0046] ② Use RocPro3D to simulate multiple rockfall paths. In each simulation, the initial position, velocity, and friction coefficient of the rocks will be different. The simulation results include multiple rockfall paths, distribution of stopping points, energy changes, etc.

[0047] ③ For each simulation, calculate the Wasserstein distance and evaluate the error between the simulated stop point distribution and the actual stop point distribution.

[0048] ④ Use the MADS algorithm to adjust the parameter set based on the Wasserstein distance obtained in each simulation. This minimizes the error in the simulation results.

[0049] ⑤ Multiple simulations were conducted. In each optimization process, multiple rockfall paths were simulated using RocPro3D, and the optimal simulation results were obtained through statistical analysis. Each simulation trajectory generated a set of stop points and energy data. Finally, the Wasserstein distance was used to evaluate the difference from the actual data, and the simulation parameters were further optimized.

[0050] Specifically, this invention dynamically adjusts key physical parameters through an optimization algorithm. The Wasserstein distance is introduced as the objective function to quantitatively evaluate the difference between the simulated rockfall deposition area distribution and the potential deposition distribution determined by actual observation or experts. The Wasserstein distance, by calculating the minimum transport cost between the two distributions, can more accurately reflect the morphological differences in spatial distribution than the traditional mean square error. During parameter optimization, the Grid Adaptive Direct Search (MADS) algorithm is used to perform parameter optimization. This algorithm establishes a dynamically adjusted grid within the parameter space (friction coefficient, coefficient of restitution, etc.) and seeks the global optimum through alternating search and extreme point phases. In the search phase, the algorithm randomly samples within the neighborhood of the current grid point. If a parameter combination with a smaller Wasserstein distance is found, the center point is updated; if no better solution is found, the algorithm enters the extreme point phase, reducing the grid size for a local fine-tuning search. Through multiple iterative calculations, the simulated deposition point distribution achieves the greatest possible match with the actual physical characteristics, thereby outputting a calibrated optimal parameter set, significantly reducing the subjective error of human-based parameter assignment.

[0051] Step 4: Implementation of Multi-Field Coupling Analysis in RocPro3D 4.1 The impact of simulated precipitation and groundwater flow The effect of precipitation on the friction coefficient: As precipitation increases, the friction coefficient of the rock surface may change. This can be addressed by adjusting the friction coefficient... This can be achieved. For example, as precipitation increases, the coefficient of friction may decrease, making sliding easier. This can be expressed by the following formula: ; in It is the initial coefficient of friction; It is the precipitation sensitivity coefficient, which indicates the degree to which precipitation affects the friction coefficient. It is usually determined based on experimental data. It refers to precipitation, which can be measured in millimeters per hour or millimeters per day, and 0 ≤ 1 - α·precipitation < 1.

[0052] The impact of groundwater level on rock mass stability: A rising groundwater level increases pore pressure in the rock mass, reducing its strength. This effect can be simulated using the following model. Specifically: ; in, It is the initial normal restitution coefficient. It is the influence coefficient of groundwater level on the normal recovery coefficient; It is the change in groundwater level.

[0053] 4.2 Introduction of Seismic Loads: Applying seismic loads Seismic loads are typically applied as external forces to the rock mass in a model. In RocPro3D, external loads can be applied using an acceleration field: ; in: It is in time The applied seismic load, It refers to the quality of the rock mass or rock block. It is in time The applied seismic acceleration.

[0054] Specifically, the stability and motion characteristics of unstable rock masses on steep slopes are significantly affected by environmental factors such as precipitation, groundwater, and earthquakes. This embodiment introduces a multi-field coupling correction mechanism into the physical model. For precipitation, a correction model is established to reflect the change in friction coefficient with precipitation. As precipitation increases, the initial friction coefficient is reduced using a precipitation sensitivity coefficient to simulate the lubrication effect of rainwater on the contact surface. The reduction formula ensures that the friction coefficient decreases non-linearly with increasing precipitation. For groundwater, a functional relationship is established between groundwater level change and normal restitution coefficient to simulate the change in collision energy recovery efficiency caused by increased pore water pressure. For seismic loads, an external acceleration field that varies with time is applied to the geometric model of the unstable rock mass to calculate the additional inertial force induced by the seismic force. The acceleration field is set according to the seismic intensity and dynamic response characteristics of the study area. By introducing these coupling factors, this method achieves dynamic prediction of the trajectory and energy characteristics of unstable rock masses under different working conditions, ensuring that the simulation results can realistically reflect the disaster evolution process under extreme environments.

[0055] Step 5: Substitute the parameters to perform calculations and dynamic feature extraction.

[0056] The optimal parameter set optimized in step S3 and the environmental correction factors determined in step S4 are substituted into the physics engine to perform the final numerical calculation. During the calculation, key physical quantities in the movement of the unstable rock mass, including movement velocity, impact kinetic energy, and bounce height, are monitored and extracted in real time. (See attached...) Figure 5 As shown, by performing profile analysis on the simulated motion trajectory, the curves of how each physical quantity changes with the distance of motion are obtained.

[0057] like Figures 6-11Calculations show that the bounce height and speed of falling rocks as they move downhill are controlled by both slope hardness and topographic gradient. In areas with exposed bedrock, the bounce height of falling rocks increases significantly; however, when entering gullies or areas covered with loose gravel and soil, the speed decreases rapidly due to increased frictional resistance. Furthermore, the volume of falling rocks is positively correlated with the length of their path and their impact kinetic energy; the larger the volume, the higher the destructive energy they carry and the stronger their potential destructive force on downstream threatened objects.

[0058] Step Six: Risk Assessment of Rock Mass Instability and Movement The risk is assessed based on the simulation results, and the criteria are shown in the table below: Table 1 Hazard Classification of Rock Mass Table 2 Classification of Hazardous Rock Masses Table 3 Classification of Hazard Levels of Dangerous Rock Masses Specifically, based on the dynamic characteristic data obtained from the above simulation, and combined with the inherent properties of the unstable rock mass and its importance to the surrounding environment, a risk assessment is conducted. The risk assessment system consists of three parts: hazard classification, hazard object classification, and hazard level matrix.

[0059] First, based on the stability calculation results (e.g., unstable, basically stable) and size (e.g., extra-large, large, medium, small) of the unstable rock mass, the hazard level is classified into three levels: large, medium, and small. For example, an unstable rock mass with poor stability and an extra-large size is classified as high hazard.

[0060] Secondly, the threats are classified according to their importance. Level 1 threats include key parts of Level 1 permanent hydraulic structures, resettlement areas and camps; Level 2 threats include key parts of Level 2 and Level 3 hydraulic structures; and Level 3 threats include Level 4 and Level 5 hydraulic structures and permanent bridges in the work area.

[0061] Finally, the hazard level matrix was used to cross-couple the hazard classification with the hazard object classification, ultimately determining the risk into four levels: I, II, III, and IV. Level I represents extremely high risk, and Level IV represents low risk. Based on the evaluation results, a comprehensive hazard zoning map of the study area was generated, marking the boundaries of the hazard zones and labeling the hazard objects with different colors.

[0062] Example 2: Taking a steep slope as an example, the risk of unstable movement of the unstable rock mass is evaluated.

[0063] Implementation steps: 1. Data Acquisition: Extract high-precision terrain data based on basic geological data, UAV imagery, etc.

[0064] Aerial surveys of the study area were conducted using a DJI Mavic 3D drone equipped with a five-lens oblique photography camera, obtaining a complete 3D model, stereoscopic images, and a digital elevation model of the entire survey area. An engineering project was created in Smart3D software, importing the original images captured by the drone. The CGCS2000 coordinate system was set as the spatial reference for the output model. Finally, the OSGB format was selected for output, generating a three-dimensional geological model of the study area. Parameters are shown in the table below: Table 4 Technical parameters of UAVs Using modern surveying methods such as manual field investigations and drone aerial photography, the distribution range of unstable rock masses in the area is determined. Under the condition of potential rockfall events, probabilistic simulations are used to construct the spatial propagation characteristics of unstable rock mass movement, and further coupled with the threatened objects. For example... Figure 2 As shown.

[0065] 2. Simulation of Rock Mass Instability and Rockfall Movement: A three-dimensional model of a steep, unstable rock mass is established based on acquired terrain data. The CSF algorithm is used to optimize the surface model within the three-dimensional model, and physical parameters related to rock mass movement are set, including rock mass mass, friction coefficient, and coefficient of restitution. For example... Figure 3 As shown.

[0066] 3. Minimize the Wasserstein distance as the objective function and iteratively update key parameters in the simulation (such as the normal restitution coefficient). Tangential coefficient of restitution and energy dissipation coefficient The initial values ​​of the specific parameters are as follows: Table 5 Initial Parameter Values 4. Simulation results are as follows Figures 5-11As shown, in terms of motion patterns, after a rockfall occurs in the study area, the bounce height, speed, and impact kinetic energy of the rockfall moving down the slope are controlled by changes in slope surface hardness caused by whether the bedrock is exposed or whether it is covered by colluvial rock debris and loose gravel. The trajectory is controlled by micro-geomorphic conditions such as the slope gradient. For example, after a rockfall, unstable rocks near the source of a gully or on the lateral slope of the gully will enter the gully within a certain period of time, collide with the bottom of the gully, and then move down the gully. When a gully changes direction, some falling rocks may bounce and tumble over obstacles, but most will follow the gully's turn. The larger the volume of the falling rocks, the farther their path and the greater their kinetic energy, but their speed and bounce height are not closely related to their volume. The steeper the slope of the terrain where the falling rocks are located, the greater their bounce height, speed, and corresponding impact kinetic energy. If the gully contains colluvial rock debris and loose gravel after the falling rocks enter, their speed will decrease significantly, and the probability of them accumulating in the gully will increase considerably. 5. Risk Assessment of Rock Mass Instability and Movement: Based on the three-dimensional simulation analysis of movement characteristics, the hazards of unstable rocks, rock clusters, rock masses, isolated boulders, and isolated rock clusters in the study area can be classified into four levels: I, II, III, and IV. The summary statistical results of the hazards of various types of unstable rock masses in the study area are shown below. Figure 12 As shown in the figure, the hazardous rock masses in the hub area are mainly classified as Class II and III, followed by Class I, with some dangerous rocks and rock masses classified as Class IV.

[0067] In summary, this invention integrates advanced surveying and mapping technology, physical simulation technology, and mathematical optimization algorithms to construct a complete and scientific risk assessment process for steep slopes and unstable rock masses. This method not only solves the problems of insufficient simulation accuracy and difficulty in parameter acquisition in traditional assessments, but also improves the applicability of risk assessment under complex working conditions through multi-field coupling analysis. The generated risk zoning map and dynamic physical parameter curves provide intuitive and accurate decision-making basis for geological disaster prevention and control, significantly improving the disaster reduction and prevention level of major projects in mountainous areas. The operating principle and process demonstrated in this embodiment fully prove the feasibility and advancement of this method in actual complex geological environments, realizing scientific management of the entire life cycle of potential collapse hazards.

Claims

1. A method for assessing the risk of instability and movement of unstable rock masses on steep slopes, characterized in that, Includes the following steps: S1: High-resolution image data of the study area is acquired by drones, and a real-scene 3D model is generated using modeling software. Point cloud data is then extracted. The cloth simulation filtering (CSF) algorithm is used to classify ground points and vegetation points in the point cloud, remove non-ground obstacles, and construct a digital terrain model (DTM) with a grid accuracy of 2m×2m. S2: Based on the DTM, establish a three-dimensional unstable rock mass motion simulation model in a three-dimensional physical simulation environment; set the physical properties and initial motion conditions of the three-dimensional unstable rock mass motion simulation model; simulate the motion trajectory of the rock mass after instability through a collision rebound model and record the spatial distribution data of the motion trajectory, bounce height, motion speed and impact kinetic energy. S3: Introduce Wasserstein distance as the objective function to measure the difference between the simulated rockfall accumulation area distribution and the preset distribution. Use a grid adaptive direct search algorithm to iteratively optimize the parameters. When the Wasserstein distance reaches the preset convergence threshold, output the optimal set of physical parameters. S4: Establish a multi-field coupled correction model that includes precipitation, groundwater and seismic load, and dynamically correct the physical parameters and instability criteria by using the precipitation sensitivity coefficient, groundwater level function relationship and external acceleration field; S5: Substitute the optimized parameters and multi-field coupling conditions to perform risk calculation. Based on the simulated rockfall trajectory, energy distribution and dwell range, combined with the importance level of the threatened object, output the risk assessment result of unstable rock mass movement.

2. The method according to claim 1, characterized in that, The CSF algorithm is implemented as follows: the original point cloud is projected onto a horizontal plane, and a simulated cloth mesh consisting of mass points and springs is constructed above the point cloud. The mass points fall to contact the ground points in the point cloud under the action of gravity. The deformation of the cloth mesh is controlled by minimizing the objective function containing smoothness and friction terms. The displacement step size of the CSF algorithm is set to a value that matches the mesh resolution. Through multiple iterations, the displacement increment of the mass points is less than the preset tolerance, and the point cloud is classified into ground points and non-ground points.

3. The method according to claim 2, characterized in that, The smoothness term: ； In the formula, It is a point height, It is a set of ground points. It is a point The neighborhood set; It is a point a certain neighborhood point Height; Friction term: ； In the formula, Point and The friction between them and It is the normal vector between points.

4. The method according to claim 1, characterized in that, In step S2, the physical properties include rock mass and friction coefficient. Normal recovery coefficient Tangential recovery coefficient and energy dissipation coefficient The initial motion conditions include three-dimensional spatial coordinates, geometry, initial linear velocity, and initial angular velocity; the geometry adopts a random polyhedron model.

5. The method according to claim 1, characterized in that, The Wasserstein distance formula is as follows: ； in: This is the actual distribution of stop points in the data. It is the distribution of dwell points in the simulated data; and These are the actual and simulated locations of the rest stops, respectively. It is the transportation flow between simulated stop points and actual stop points. It is the actual distribution The Middle The location of each stop It is a simulated distribution The Middle The location of each stop.

6. The method according to claim 1, characterized in that, The specific steps of step S3 are as follows: S31. Select the initial parameter set ; S32. Use three-dimensional physical simulation to simulate multiple sets of random rockfall paths and generate statistical samples of stopping points; S33. Calculate the Wasserstein distance for each group of samples, which serves as the search guide for the grid adaptive direct search algorithm, and perform grid partitioning and polling search in the parameter space. S34. Continuously adjust the values ​​of the normal restitution coefficient, tangential restitution coefficient, and energy dissipation coefficient until the global optimal solution that minimizes the Wasserstein distance is found.

7. The method according to claim 1, characterized in that, In step S4, the precipitation coupling mechanism corrects the friction coefficient of the contact surface through the precipitation sensitivity coefficient, and the friction coefficient decreases nonlinearly with the increase of precipitation; the groundwater coupling mechanism simulates the influence of pore water pressure on collision energy dissipation through the functional relationship between groundwater level change and normal recovery coefficient.

8. The method according to claim 1, characterized in that, The seismic load coupling mechanism calculates the additional inertial force of the unstable rock mass under seismic force by applying a time-varying external acceleration field to the physical model, and calculates the dynamic tensile stress and dynamic shear stress of the unstable rock mass during the vibration process based on the input seismic wave sequence.

9. The method according to claim 1, characterized in that, In step S5, the risk assessment system extracts the maximum bounce height, maximum impact kinetic energy, and maximum speed from the simulated trajectory as evaluation factors, compares the physical quantities with the disaster resistance level of the threatened object to calculate the probability of damage, and determines the interception and protection location, high-risk area, and medium- and low-risk area according to the risk level.

Images