In this study, we used the DJI M300 RTK drone for aerial photography. The drone’s 3D flight path was designed to continuously adjust its altitude, maintaining a consistent distance from the ground to ensure a uniform ground sampling distance (GSD) in the captured images. The drone was equipped with high-precision cameras capable of multi-angle measurements, capturing images from four oblique angles and one vertical angle along each flight path. This configuration provided high-resolution texture data for both the top and side views of the target surface.

Control points were first established within the mining area to correct spatial positioning and attitude errors, improving the accuracy of the aerial data (Figure 2). A low-resolution flight was initially conducted to collect terrain data, which was used to create a preliminary model, providing basic information on coordinates and elevation changes. Based on this data, the drone’s flight paths were further refined, and slope-following flights were conducted. Due to the irregular development of rock masses, the camera’s distance to the surface was continuously adjusted to maintain a consistent relative position. The relative distance was calculated using Equation 1. Since the drone’s attitude could not be precisely determined, the efficiency and accuracy of image matching were reduced, requiring sufficient overlap between images. The forward and side overlap percentages were calculated using Equations 2, 3, respectively. H = G S D × f S p (1) P f = Q x Q y × 100 % (2) P l = L x L y × 100 % (3) where, H represents the distance from the camera’s perspective center to the surface of the object (in meters). The closer the distance, the higher the resolution. f is the camera focal length (in mm), and S p is the pixel size (in pixels). The resolution is also defined as the ground sampling distance (GSD). L x and L y represent the image width and height (in meters), while Q x and Q y are the overlap lengths in the horizontal and vertical directions.

Based on site conditions, the forward overlap was set to 85%, and the side overlap to 70%. The initial flight altitude was set at 300 m, while the slope-following flight altitude was set at 60 m, with an optimized area height of 5 m. In total, 100,000 images were captured, generating approximately 2 TB of data. The data was processed using a 10-node Smart3D software parallel computing setup, resulting in an ultra-high-resolution model with a precision of 10 mm, as shown in Figure 1.

The extraction of rock mass structural planes involves calculating the normal vectors of all point clouds and clustering these vectors to identify different groups of structural planes. However, a preliminary grouping of structural planes is insufficient for calculating parameters related to rock mass stability. To address this, this study utilizes the refined 3D point cloud data obtained in Section 2.1. The Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm is applied to further refine and classify the initially grouped structural planes.

First, using the oblique photogrammetry model, the point cloud orientation information of the rock mass is extracted. The K-means clustering algorithm is then employed to perform an initial grouping of structural plane orientations. Next, octree segmentation is used to determine the normal direction of the fitted plane within each cubic voxel. Finally, the Fuzzy C-Means (FCM) clustering algorithm is used to group point clouds with similar orientations, allowing for the extraction of different sets of structural planes (Eivazy et al., 2017). However, FCM is sensitive to the initial random selection of cluster centers, which can lead to different clustering results and cause the algorithm to converge to local minima instead of the global optimum. And FCM is highly sensitive to the choice of parameters, particularly the fuzziness coefficient, which requires experimental tuning to find the optimal setting, adding complexity to parameter selection. The results of the classification and extraction process are shown in Figure 3.

The refinement of structural plane classification relies on the preliminary classification and extraction of structural planes. After completing the initial classification, the DBSCAN algorithm (Schubert et al., 2017) is introduced to further refine the classification and extraction of structural planes. The DBSCAN algorithm is a density-based clustering method that defines clusters using a neighborhood radius (γ) and a minimum number of sample points (MinPts). Its key advantages include: no need to preset the number of clusters, the ability to detect clusters of arbitrary shapes, natural handling of noise points, and adaptability to uneven data density distributions. These features make DBSCAN particularly suitable for geological structural plane identification, as it effectively captures complex data distributions and excludes outliers.

The basic workflow of the algorithm is as follows:

Input the target point set, D = {P ₁,P ₂, … ,P _n};

To validate the reliability of the automatic structural plane recognition algorithm, point cloud data from the eastern slope of the Yanshan mine (validation area in Figure 1) was selected for structural plane identification. The area of the validation region is 142.3 m × 30.2 m. The identification results in Figure 5 depict the same set of structural planes, with different colors used for clear differentiation. The automatic structural plane recognition program identified a main set of 78 joints, with an average orientation of 245.4°∠45°. The average spacing and trace length were 0.78 m and 1.54 m, respectively. A comparative validation with previous literature (Deng et al., 2021) was performed, where data was manually extracted from UAV models, accurately representing the structural plane features. The comparison results, shown in Table 1, indicate that the method performs well in clustering dip direction, dip angle, and spacing. However, in some cases, a single structural plane was split into two clusters, shortening the trace length and increasing deviation from the original data. This issue will require further improvement in future work. Using the UAV oblique photogrammetry model and the automatic structural plane recognition system, a total of 2,895 joints were identified in this study, as shown in Figure 6.

Figure 7 shows the statistical distribution of the identified structural plane information. The dip angle and dip direction follow a normal distribution, while the trace length follows a log-normal distribution with a base of 10. The mean values of the dip angle, dip direction, and trace length are 44.5°, 249.9°, and 0.78 m, respectively. The coefficient of variation (CV), defined as the ratio of the standard deviation to the mean, was calculated to measure the dispersion of the data. The CVs for dip direction, dip angle, and trace length are 11.6%, 15.7%, and 131.4%, respectively.

Both dip direction and dip angle exhibit low coefficients of variation, each below 20%, indicating minimal fluctuation in these parameters. However, the CV for trace length is as high as 131%, significantly exceeding that of dip direction and dip angle. This indicates that the trace length distribution is highly uneven, with notable differences between long and short joints or the strong influence of a few long-trace joints on the overall statistics. This high variability likely reflects the strong influence of geological structures and stress fields on joint development, leading to the pronounced heterogeneity of trace length distribution within the rock mass.

Statistical analysis of the structural plane proportions in different areas of the open-pit mine (Figure 8) shows that the northern slope and northeastern slope account for 29.8% and 48.2% of the joints, respectively, while the eastern slope accounts for 20%. The eastern slope has two main sets of structural planes (266°∠44° and 293°∠44°), with a joint density of 0.83 joints/m³. The northern slope has two main sets of structural planes (206°∠45° and 237°∠42°), with a joint density of 1.71 joints/m³. The northeastern slope has two main sets of structural planes (240°∠44° and 266°∠42°), with a joint density of 2.68 joints/m³. Calculations of the joint connectivity rate in different areas reveal that the eastern slope has the lowest connectivity rate at 0.53, followed by the northern slope at 0.61, and the northeastern slope with the highest rate at 0.90. Based on the evaluation of rock mass fragmentation and connectivity rate, the northeastern slope and northern slope are identified as critical areas requiring focused attention.

To determine the slope failure patterns, the slope orientations and structural plane information were analyzed using the Dips software for kinematic analysis (Smith, 2015). Figure 9 illustrates the failure probabilities in different areas as the slope angle increases. As the slope angle increases from 35° to 49°, the failure probability consistently rises across all regions. The critical point of significant failure increase is at 39°, where the failure probability exceeds 20% across all regions. After 43°, the failure probability accelerates rapidly, peaking at 49° (with wedge failure nearing 70% and planar sliding reaching 55%). Overall, maintaining the slope angle below 43° is a key measure to reduce slope failure risks.

In terms of failure types, wedge failure consistently dominates, with probabilities higher than planar sliding at all slope angles. This gap widens further when the slope angle exceeds 43°. In high slope angle regions (47°–49°), wedge failure probabilities reach 65%–70%, while planar sliding probabilities are between 50% and 55%. This indicates that increased slope angles have a more significant impact on wedge failure. Regionally, the eastern slope exhibits the highest failure probability at all slope angles, followed by the northeastern slope, with the northern slope having the lowest probability. However, considering the total number of failure events, the northeastern slope has the highest number of events, followed by the northern slope, with the eastern slope having the least.

RocSlope uses a test matrix to determine the geometric mobility of blocks (Goodman and Shi, 1985). This is calculated based on the relationship between the normal direction of the joint planes and the block position, as shown in Equation 4. T i j k = sign n i ^ × n j ^ · n k ^ I a k (4) where, n i ^ is the upward unit normal of joint face i; I a k represents the value corresponding to joint i in the signed block code; For a joint k, its value in the signed block code is defined as Equation 5: I a k = + 1 , if the block lies above the plane defining joint face  k − 1 , if the block lies below the plane defining joint face  k ± 1 , if theblock lies between  2  planes parallel to joint face  k (5)

Then, according to the Rule for Testing Finiteness in the referenced text, if every row of the testing matrix includes both positive and negative terms, the block is not removable; otherwise, the block is geometrically removable.

Once it is determined that the block is capable of movement, all individual force vectors acting on the sliding block are identified to calculate the driving and resisting forces. Generally, the driving force represents the motive force in the safety factor calculation, while the resisting force represents resistance, as defined by Equations 6, 7. A = W + C + X + U + E + B a + W s (6) where, A is the resultant active force vector; W is the block weight vector; C is the load force vector; X is the active pressure force vector; U is the water force vector; E is the seismic force vector; B a is the active bolt force vector; W s is the shotcrete weight vector. P = H + Y + B p (7) where, P is the resultant passive force vector; H is the shotcrete shear resistance force vector; Y is the passive pressure force vector; B p is the resultant passive bolt force vector.

The sliding direction of the block is determined by the driving force and is not influenced by the resisting force. When considering sliding along multiple joints, the direction must satisfy the following inequality as defined by Equations 8, 9. s i j ^ · n k i ^ = 0 (8) A · n k i ^ = 0 (9) where, s i j ^ is the sliding direction on joints i and j (along the line of intersection).

Once the sliding direction is determined, the normal force and shear strength on each joint plane can be calculated based on the sliding direction to evaluate the resistance of the joint planes (Mauldon and Ureta, 1996). The normal force is decomposed along the sliding direction, while the shear strength is determined by the joint plane’s frictional force and cohesion. Finally, the safety factor (FoS) is calculated using the limit equilibrium principle to assess the block’s stability.

The high-precision terrain model obtained from UAV oblique photogrammetry was imported into the RocSlope software. Through boundary reconstruction, a basic geometric model for block kinematic analysis was generated. The model dimensions are 1650 m × 1450 m × 400 m, as shown in Figure 10. To construct the slope geometric block failure model, the previously identified structural plane information was incorporated into the model, generating structural plane stereonets with orientation and trace length data. The calculations primarily focused on the safety factors under three conditions: natural conditions, water-saturated conditions, and water-saturated plus blasting-induced vibration.

During block failure, the resisting force of the sliding block is primarily provided by the shear strength of the joint planes. Among them, four failures occurred under natural conditions, while six failures developed after rainfall. For the natural-condition cases, the geometric characteristics of each failed block were obtained through detailed field investigation. We then performed sensitivity analyses in RocSlope by assigning different combinations of cohesion and friction angle to the joint planes. The parameters that yielded a factor of safety equal to 1 were taken as the calibrated values for each case. This process resulted in several feasible friction–cohesion pairs, from which the representative values of 25.1° (friction angle) and 10 kPa (cohesion) were determined.

Under water-saturated conditions, since the structural planes primarily exist in the surface layer of the slope, water pressure contributes minimally to block failure. However, water influences failure in two key ways: increasing the unit weight of the rock mass and reducing the shear strength of the structural planes. To account for the effects of water saturation, we referred to laboratory tests reported in Zhang Z. et al. (2025), which show that long-term immersion causes a 0.1% increase in rock unit weight, a 1%–5% reduction in friction angle, and a 15%–25% reduction in cohesion. These degradation trends were combined with back-analysis of five rainfall-induced failures. Based on these rainfall-related cases, the representative mechanical parameters under saturated conditions were determined to be 24.6° for friction angle and 8 kPa for cohesion. Under blasting-induced vibration conditions, the impact of the blasting influence factor on failure results was analyzed. Based on on-site blasting vibration tests, the failure influence factor was calculated to be 0.015. The mechanical parameters used in this calculation are summarized in Table 2.

Geometric mobility calculations identified 317 potentially movable rock blocks, most of which are lower than the height of a single bench (30 m). Among these, 61.1% are located in the northeastern slope, 38.6% in the northern slope, and only 0.3% in the eastern slope. The areas with a higher concentration of movable blocks closely align with the zones of highly fractured rock observed in field investigations, as shown in Figure 11a. Based on the exposed area and depth of the blocks (Figures 11b,c), most blocks have an area of 1–20 m², accounting for 84.2%, and a thickness of 0–3 m, representing 84.9%. These findings suggest that wedge-type failures in the Yanshan open-pit mine are likely to occur in blocks with heights below the bench height, exposed areas smaller than 20 m², and depths under 3 m. Such failure modes are concentrated in zones with fractured rock and warrant particular attention.

According to the regulations of the Ministry of Housing and Urban-Rural Development of the People’s Republic of China, slopes higher than 200 m must have a safety factor greater than 1.2. Under natural conditions, 17 blocks (5.4% of the total) had safety factors below 1.2, failing to meet the design requirements. All these blocks were located in the northeastern and northern slopes. Under water-saturated conditions, rainfall infiltration increases the saturation degree along persistent joint planes, and even a slight rise in pore water pressure reduces the effective normal stress acting on these planes, thereby decreasing their shear strength according to the Mohr–Coulomb criterion. This hydro-mechanical effect accounts for the experimentally observed reductions in friction angle (−2%) and cohesion (−20%). In addition, long-term water immersion slightly increases the unit weight of the fractured rock mass and accelerates weathering and micro-cracking along discontinuities, further weakening joint shear resistance and promoting block mobility. As a combined consequence of these rainfall-induced hydro-mechanical processes, the average safety factor of the blocks decreased by 13.6% compared with natural conditions, increasing the number of blocks that failed to meet the design standard to 23 (7.3% of the total). Under combined “water-saturation + blasting” conditions, the safety factor decreased further by 15.5% compared to natural conditions. The number of blocks failing to meet design requirements remained at 23, accounting for 7.3% of the total.

Based on the comprehensive analysis, the northeastern and northern slopes are high-risk areas for slope failure, with their stability significantly reduced under water-saturated and blasting conditions. In future slope design and stability monitoring efforts, these areas should be closely monitored, and slope designs optimized with appropriate support measures to enhance overall slope stability. Additionally, for potential wedge failures, precise treatment of fractured rock zones should be implemented to effectively reduce the risk of landslides.

In August 2023, following several rainfall events, a small landslide occurred in the reinforced test area shown in Figure 11a after normal blasting operations. The landslide, measuring 30 m in length, 15 m in height, and 3 m in thickness, did not cause any casualties or equipment damage, as shown in Figure 12. This area is cut by multiple joint sets, which caused the rock to experience sequential slab peeling. Additionally, this landslide confirmed the reliability of the intelligent hazard zone identification method proposed in this study.

Figure 13 shows the simulated sequence of rock slippage for this landslide. The entire model is 100 m long, 20 m wide, and 30 m high. During the entire failure process, region 1 experienced the initial slip. The slip in region 1 reduced the constraint on the right side and bottom of region 2, leading to a slip in region 2. The failure of regions 1 and 2 similarly reduced the constraint on the right side of region 3, creating a free surface on the right, which caused region 3 to slip as well. Subsequently, regions 4 and 5 also experienced sliding. A safety of factor analysis under various operating conditions showed that the average factor of safety under normal conditions is 1.15. Under the influence of water, the factor of safety decreased to 0.98, and under both water and blasting vibrations, it further dropped to 0.95. This indicates that the landslide was significantly influenced by blasting and groundwater.

Previous studies have shown that the thickness of landslides controlled by structural surfaces is generally small, and the area of a single landslide is also small. Therefore, high-ductility concrete was used for reinforcement. The principle of shotcrete is to use high-pressure equipment to spray concrete repair materials at high speed onto the slope surface and interior. This increases the cohesion and internal friction angle of the slope structure to some extent, enhancing its resistance to sliding. It also forms a thick concrete layer on the slope surface, making the blocks cohesive as a whole, improving the overall shear resistance and increasing the connectivity between individual block. Moreover, shotcrete effectively prevents water penetration, weathering, and freeze-thaw damage. The concrete parameters are shown in Table 3.

The simulation results are presented in Figure 14. In the simulation, after reinforcement, the average factor of safety under normal conditions increased from 1.05 to 1.35. Under the influence of water, the average factor of safety rose from 0.98 to 1.30, and when both blasting vibrations and water effects were considered, the factor of safety increased from 0.95 to 1.28. This demonstrates that shotcrete reinforcement has a significant impact on the stability of this type of slope. Following this, the management team began large-scale shotcreting on the slope, as shown in Figure 15. As of now, the treated slope has not experienced any further landslides.