Rock mass discontinuity automatic identification method and system based on density peak clustering

Abstract

This invention discloses an automatic identification method and system for rock mass discontinuities based on density peak clustering. The method includes: acquiring three-dimensional point cloud data of the rock mass; acquiring a normal vector field based on the three-dimensional point cloud data; automatically determining the number of clusters of dominant groups of discontinuities using a density peak clustering algorithm based on the normal vector field to obtain cluster centers; performing spatial-directional ensemble clustering on the three-dimensional point cloud data according to the cluster centers and the normal vector field to obtain an initial clustering result; optimizing the boundaries of the initial clustering result using a region growing algorithm to obtain an optimized clustering result; and performing multi-scale decomposition on the optimized clustering result to obtain the attitude parameters of the rock mass discontinuities.

Description

Technical Field

[0001] This invention belongs to the field of geotechnical engineering technology, and in particular relates to an automatic identification method and system for rock mass discontinuities based on density peak clustering. Background Technology

[0002] Discontinuities in rock mass (including joints, faults, bedding planes, and fractures) are fundamental factors controlling the mechanical, hydraulic, and deformation behavior of rock masses. These discontinuities form weak surfaces within the rock mass, dominating the rock mass response under load and controlling instability mechanisms, seepage channels, and stress distribution patterns. They are crucial for slope stability assessment, underground excavation design, and foundation analysis. The geometric properties of discontinuities (especially strike, dip, continuity, spacing, and roughness) directly determine the rock mass quality classification (RMR, Q-system) and control engineering response, making their accurate characterization fundamental to geotechnical engineering practice.

[0003] Traditional characterization of discontinuities relies on manual field measurements using geological compasses, survey lines, and window mapping techniques. Despite decades of procedural improvements, these methods suffer from five persistent limitations: (1) severely limited sampling coverage due to accessibility constraints; (2) safety concerns in hazardous terrain limiting measurement locations; (3) subjective measurement bias among surveyors, particularly when measuring irregular or weathered discontinuous surfaces; (4) considerable time commitment (2–4 hours per station); and (5) inability to capture spatial discontinuities beyond discrete measurement points. These limitations become particularly acute in complex geological environments requiring comprehensive three-dimensional characterization.

[0004] The advent of remote sensing technologies, particularly terrestrial laser scanning (TLS) and unmanned aerial vehicle (UAV) photogrammetry, has enabled the non-contact acquisition of dense 3D point clouds representing the entire geometry of a rock surface with sub-centimeter precision. This technological advancement theoretically overcomes the sampling and accessibility limitations of manual methods, capturing the complete spatial distribution of discontinuous surface networks. However, transforming these massive geometric datasets into meaningful geological interpretations presents significant computational and methodological challenges that current automated analysis methods have not yet fully addressed.

[0005] There are currently three fundamental gaps that hinder fully automated and geologically reliable characterization of discontinuities from point cloud data: First, computational and algorithmic efficiency remains an issue. Traditional clustering algorithms applied to high-dimensional, irregularly distributed point cloud datasets exhibit quadratic or higher computational complexity, making the analysis impractical for engineering workflows that require timely representation of discontinuous surfaces.

[0006] Second, the reliance on subjective parameters undermines reproducibility and objectivity. Most clustering-based methods require manually specifying the number of dominant groups on discontinuities through subjective interpretation of stereo projections, elbow coefficients, or profile coefficients. This requirement introduces analyst dependence and variability: analysts report variations of ±1-2 clusters in a significant proportion of complex geological environments, directly propagating uncertainty to subsequent stability analyses and rock mass ratings. A fundamental limitation is that traditional partitioning algorithms (K-means, fuzzy C-means) lack intrinsic validity measures for assessing cluster compactness at different scales, requiring external validation by geological experts, which contradicts the goal of automation.

[0007] Third, the methodological separation between directional analysis and spatial segmentation produces geologically inconsistent results. Existing methods overwhelmingly treat these as sequential, independent processes: first, clustering points based on normal vector similarity in directional space, then applying region growing in Euclidean space. This decoupling strategy fails to consider the physical reality that geological discontinuities exhibit both directional consistency and spatial continuity. Consequently, automated results often reveal boundary misclassifications, particularly at the intersections of multiple interacting discontinuities—precisely the most critical locations for wedge instability and block stability assessments. Furthermore, existing single-scale analyses either identify the dominant set of discontinuities (through clustering) or individual fracture surfaces (through segmentation), but rarely identify both simultaneously within an integrated framework. However, comprehensive geotechnical characterization requires multi-scale decomposition: the dominant set controls the instability pattern (wedge, planar, overturning), while fine-scale variations (spacing statistics, local directional discretization) determine block size and specific instability locations crucial for rock mass quality systems (RMR, Q-system, GSI).

[0008] Existing methodological approaches can be divided into four frameworks, each exhibiting different limitations relative to the gaps mentioned above: (1) Clustering-based methods use K-means, fuzzy C-Means or DBSCAN variants to partition the normal vectors in the three-dimensional space; (2) Geometric fitting methods use least squares, RANSAC or PCA-based plane fitting to identify planar surfaces; (3) The region growing method aggregates spatial connection points with similar normal vectors, maintaining spatial continuity but is highly sensitive to seed point selection and growth threshold; (4) Deep learning methods are promising for pattern recognition in noisy data, but require extensive labeled training datasets, have limited transferability between geological settings, and lack the interpretability necessary for engineering applications.

[0009] Therefore, there is an urgent need for an automatic identification method and system for rock mass discontinuities based on density peak clustering. This method can automatically determine the number of clusters, comprehensively consider spatial and directional attributes, and achieve multi-scale characterization. Summary of the Invention

[0010] To address the aforementioned technical problems, this invention proposes an automatic identification method and system for rock mass discontinuities based on density peak clustering. This method can automatically determine the number of clusters, comprehensively consider spatial and directional attributes, and achieve multi-scale characterization.

[0011] To achieve the above objectives, this invention provides an automatic identification method for rock mass discontinuities based on density peak clustering, comprising: Acquire 3D point cloud data of the rock mass; Based on the aforementioned 3D point cloud data, obtain the normal vector field; Based on the normal vector field, the number of clusters of discontinuous surfaces is automatically determined by the density peak clustering algorithm, and the cluster centers are obtained. Based on the cluster centers and the normal vector field, spatial-directional ensemble clustering is performed on the 3D point cloud data to obtain the initial clustering result; The initial clustering results are optimized by performing boundary optimization using a region growing algorithm to obtain optimized clustering results. The optimized clustering results are decomposed at multiple scales to obtain the attitude parameters of the rock mass discontinuities.

[0012] Optionally, obtaining the normal vector field based on the three-dimensional point cloud data includes: The three-dimensional point cloud data is preprocessed to obtain preprocessed point cloud data; Based on the preprocessed point cloud data, the normal vector field is obtained by calculating the normal vector through principal component analysis.

[0013] Optionally, the three-dimensional point cloud data is preprocessed to obtain preprocessed point cloud data, including: The three-dimensional point cloud data is downsampled using voxels to obtain the downsampled point cloud; Outlier points in the downsampled point cloud are identified and removed to obtain the preprocessed point cloud data.

[0014] Optionally, based on the preprocessed point cloud data, normal vectors are calculated using principal component analysis to obtain the normal vector field, which includes: An adaptive neighborhood is constructed for each point in the preprocessed point cloud data based on the KD tree structure. Principal component analysis is performed on the points within the adaptive neighborhood to obtain the covariance matrix; The covariance matrix is ​​decomposed into eigenvalues, and the eigenvectors corresponding to the smallest eigenvalues ​​are determined as the normal vectors of each point. The surface curvature of each point is calculated based on the magnitude of the minimum eigenvalue, and points with curvature higher than the threshold are marked as potential boundary points. The normal vector field is obtained based on the normal vector and the surface curvature.

[0015] Optionally, based on the normal vector field, the number of clusters of dominant groups of discontinuities is automatically determined by the density peak clustering algorithm, and the cluster centers are obtained, including: The normal vectors in the normal vector field are converted into spherical coordinates, and the similarity between the vectors is calculated using the dot product. The local density of each normal vector in the normal vector field is calculated based on the spherical coordinate system representation. Calculate the minimum distance from each normal vector to the nearest normal vector with higher local density; The γ value of each normal vector is determined based on the product of the local density and the minimum distance, where the γ value is used to determine whether the corresponding point is a cluster center; The γ values ​​are sorted from high to low, and the cluster centers are obtained based on the sorted γ values.

[0016] Optionally, based on the cluster centers and the normal vector field, spatial-directional ensemble clustering is performed on the 3D point cloud data to obtain initial clustering results, including: Calculate the weighted distance from each point to each cluster center, where the weighted distance is obtained by weighted summation of spatial distance and normal vector similarity terms; Each point is assigned to the nearest cluster center based on the weighted distance; Iteratively update the cluster centers until the change in the cluster centers is less than the convergence threshold, and obtain the initial clustering result.

[0017] Optionally, the initial clustering results can be optimized by using a region growing algorithm to obtain optimized clustering results, including: Based on the initial clustering results, high-confidence points are selected from each cluster as seed points; For each seed point's neighboring points, calculate the normal vector consistency weight and spatial distance weight between the neighboring points and the seed point; Based on the normal vector consistency weight and the spatial distance weight, the final clustering assignment of each neighboring point is determined through a weighted voting mechanism to obtain the optimized clustering result.

[0018] Optionally, the optimized clustering results can be decomposed at multiple scales to obtain the attitude parameters of the rock mass discontinuities, including: For each cluster in the optimized clustering results, a hierarchical density clustering algorithm is used to perform sub-clustering decomposition, identify fine-scale substructures in the main discontinuities, and obtain a multi-scale clustering structure. For each cluster and sub-cluster in the multi-scale clustering structure, a plane fitting is performed using the random sample consensus algorithm to obtain the fitted plane equation; Based on the normal vector of the fitted plane equation, the dip and dip angle parameters are calculated to obtain the attitude parameters of the rock mass discontinuity surface.

[0019] The automatic identification system for rock mass discontinuities based on density peak clustering includes: a data acquisition module, a normal vector estimation module, an adaptive clustering module, a region growth optimization module, and a multi-scale feature extraction and optimization module; The data acquisition module is used to acquire three-dimensional point cloud data of the rock mass; The normal vector estimation module is used to obtain the normal vector field based on the three-dimensional point cloud data; The adaptive clustering module is used to automatically determine the number of clusters of dominant groups of discontinuities based on the normal vector field and the density peak clustering algorithm, thereby obtaining the cluster centers. The region growing optimization module is used to perform spatial-directional ensemble clustering on the 3D point cloud data based on the cluster centers and the normal vector field to obtain an initial clustering result; and to perform boundary optimization on the initial clustering result using a region growing algorithm to obtain an optimized clustering result. The multi-scale feature extraction and optimization module is used to perform multi-scale decomposition on the optimized clustering results to obtain the attitude parameters of the rock mass discontinuity surface.

[0020] Compared with the prior art, the present invention has the following advantages and technical effects: 1) This invention eliminates the need for manually specifying the number of clusters through an automatic parameter determination method for Density Peak Clustering (DPCA), a common limitation of traditional clustering methods. By analyzing the density distribution of normal vectors in three-dimensional space, this method automatically identifies the number of dominant groups and their main directions on optimal discontinuities, enhancing objectivity and repeatability through density peak analysis. The decision map provides a quantitative assessment of clustering effectiveness, enabling consistent results across different datasets and analysts without prior knowledge of geological structures.

[0021] 2) This invention integrates the parallel normal vector criterion with spatial-direction weighted clustering, effectively solving the directional ambiguity problem that often leads to misclassification in traditional methods. The weighted distance metric jointly considers spatial approximation and normal vector similarity within a unified framework, enabling more accurate grouping of discontinuous surfaces with similar geometric properties even when normal vectors point in opposite directions. This innovation prevents the artificial segmentation observed in methods that rely solely on normal vector direction without considering geometric equivalence, resulting in improved boundary partitioning, particularly in complex transition regions where multiple discontinuous surfaces interact.

[0022] 3) This invention implements a region growing algorithm with a weighted voting mechanism. By considering the collective influence of neighborhood point sets rather than the binary inclusion criterion, it significantly improves the boundary partitioning accuracy between adjacent discontinuous surfaces. By simultaneously combining spatial distance and normal vector similarity, the weighted voting mechanism enhances the classification robustness of boundary regions, where orientation changes and measurement noise can introduce classification uncertainty.

[0023] 4) This invention applies HDBSCAN for hierarchical sub-clustering, which can identify fine-scale structural features within major discontinuity groups, overcoming a key limitation of single-scale methods. This multi-scale capability is particularly important for comprehensive rock mass characterization, where regional structural trends and local variations influence rock behavior and hydrogeological properties. Unlike traditional clustering methods that apply a uniform density threshold, HDBSCAN adaptively identifies clusters across multiple density scales, making it well-suited for detecting dominant and secondary structural features within structural groups.

[0024] 5) This invention has undergone comprehensive validation, demonstrating advantages in accuracy and efficiency. Regarding direction measurement accuracy, this method achieves results within 3° in dip direction and within 2° in tilt angle compared to established methods. Computational efficiency is improved by 72%, reducing the processing time for the benchmark dataset from 605 seconds to 170 seconds, and its overall complexity demonstrates scalability suitable for practical field applications. The method exhibits robustness across different rock types, data acquisition methods, and point cloud sizes (from 1.5 million to 6 million points), confirming its applicability to diverse geological environments and data characteristics.

[0025] In summary, this invention can automatically determine the optimal number of dominant groups for discontinuities, integrate spatial and directional attributes for classification, and achieve multi-scale characterization from dominant sets to fine-scale variations. This invention significantly improves the automation, accuracy, and computational efficiency of discontinuity identification, providing a powerful tool for slope stability assessment, tunnel design, and geotechnical risk assessment. Attached Figure Description

[0026] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings: Figure 1 This is a flowchart of the automatic identification method for rock mass discontinuities based on density peak clustering according to an embodiment of the present invention; Figure 2 This is a geometric representation of the relationship between the normal vector and the angle in an embodiment of the present invention; Figure 3 This is a schematic diagram of the DPCA decision diagram according to an embodiment of the present invention; Figure 4This is a parameter sensitivity analysis diagram of an embodiment of the present invention; Figure 5 This is a visualization of the main discontinuities automatically identified in an embodiment of the present invention (Case A); Figure 6 This is a visualization diagram of multi-scale substructure recognition in an embodiment of the present invention (Case A), wherein (a) is a top view of the multi-scale substructure, (b) is a front view of the multi-scale substructure, and (c) is a side view of the multi-scale substructure. Figure 7 This is a visualization of the main discontinuities automatically identified in an embodiment of the present invention (Case B); Figure 8 This is a visualization diagram of multi-scale substructure recognition in an embodiment of the present invention (Case B), wherein (a) is a top view of the multi-scale substructure, (b) is a front view of the multi-scale substructure, and (c) is a side view of the multi-scale substructure. Detailed Implementation

[0027] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0028] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0029] This embodiment proposes an automatic identification method for rock mass discontinuities based on density peak clustering, such as... Figure 1 As shown, the specific steps include: Based on the aforementioned 3D point cloud data, obtain the normal vector field; Based on the normal vector field, the number of clusters of the dominant group of discontinuities is automatically determined by the density peak clustering algorithm, and the cluster centers are obtained. Based on the cluster centers and the normal vector field, spatial-directional ensemble clustering is performed on the 3D point cloud data to obtain the initial clustering result; The initial clustering results are optimized by performing boundary optimization using a region growing algorithm to obtain optimized clustering results. The optimized clustering results are decomposed at multiple scales to obtain the attitude parameters of the rock mass discontinuities.

[0030] Specifically, S1, Data Acquisition and Preprocessing: Acquire three-dimensional point cloud data of the rock mass, and preprocess the raw point cloud data through voxel downsampling, statistical outlier removal and quality assessment to ensure geometric integrity while reducing computational burden; S2, Normal Vector Estimation: Based on the adaptive neighborhood selection of KD tree structure, robust normal vector calculation is performed through principal component analysis (PCA), and consistent direction propagation is performed through minimum spanning tree; S3. Density Peak Clustering: Automatically Determines the Number of Clusters by Analyzing the Normal Vector Density Distribution in Three-Dimensional Space and Automatically Identifying the Number of Dominant Groups on the Optimal Discontinuity Surface by Detecting Density Peaks, thereby Eliminating the Need for Subjective Parameter Specification. S4. Adaptive Spatial-Directional Integrated Clustering: Enhanced K-means clustering using weighted distance metric, jointly considering spatial coordinates and normal vector similarity, and then improving boundary partitioning accuracy through region growing optimization with weighted voting mechanism; S5. Multiscale analysis and parameter extraction: Fine-scale substructures within the main discontinuities are identified through hierarchical HDBSCAN decomposition, and then geological occurrence parameters (dip and dip angle) for engineering applications are extracted through RANSAC plane fitting and SVD optimization.

[0031] Furthermore, based on the three-dimensional point cloud data, obtaining the normal vector field includes: The three-dimensional point cloud data is preprocessed to obtain preprocessed point cloud data; Based on the preprocessed point cloud data, the normal vector field is obtained by calculating the normal vector through principal component analysis.

[0032] Furthermore, the three-dimensional point cloud data is preprocessed to obtain preprocessed point cloud data, including: The three-dimensional point cloud data is downsampled using voxels to obtain the downsampled point cloud; Outlier points in the downsampled point cloud are identified and removed to obtain the preprocessed point cloud data.

[0033] Specifically, in step S1, the preprocessing includes: S11. Voxel downsampling: The voxel size is set according to the point cloud density and target resolution. The original point cloud space is divided into voxel grids of equal size. The points in each non-empty voxel are averaged to obtain the downsampled point cloud. S12. Statistical outlier removal: Based on local neighborhood analysis, the average distance to the m nearest neighbors (m=30) of each point is calculated. Points exceeding the threshold μ+2σ (where μ is the average distance and σ is the standard deviation) are classified as outliers and removed. This threshold can adapt to the density changes in the dataset and effectively eliminate measurement noise while maintaining the surface characteristics of discontinuous surfaces.

[0034] Furthermore, based on the preprocessed point cloud data, normal vectors are calculated using principal component analysis to obtain the normal vector field, which includes: An adaptive neighborhood is constructed for each point in the preprocessed point cloud data based on the KD tree structure. Principal component analysis is performed on the points within the adaptive neighborhood to obtain the covariance matrix; The covariance matrix is ​​decomposed into eigenvalues, and the normal vector of each point is determined based on the eigenvector corresponding to the smallest eigenvalue, wherein the normal vector is obtained based on the direction of the smallest eigenvalue. The surface curvature of each point is calculated based on the magnitude of the minimum eigenvalue, and points with curvature higher than the threshold are marked as potential boundary points. The normal vector field is obtained based on the normal vector and the surface curvature.

[0035] Specifically, in step S2, the normal vector estimation includes: S21. Adaptive Neighborhood Construction: Use the KD-tree data structure to construct a local neighborhood with adaptive search criteria for each point p; S22. Calculate the covariance matrix C: ; in, p is the centroid of the neighborhood points. i Let k represent each neighboring point, k be the number of neighboring points, and T be the transpose of the matrix. S23. Eigenvalue decomposition: Perform eigenvalue decomposition on the covariance matrix to obtain three eigenvalues ​​(λ1≥λ2≥λ3) and corresponding eigenvectors (v1,v2,v3). The eigenvector v3 corresponding to the minimum eigenvalue λ3 represents the estimated normal vector because it is aligned with the direction of minimum variance of the local point distribution. S24. Normal vector direction consistency processing: The ambiguity of normal vector direction is solved by using a consistency direction procedure based on minimum spanning tree (MST), ensuring that the normal vector follows a consistent pattern throughout the point cloud; S25. Local Surface Curvature Calculation: Calculate the local surface curvature at each point using the eigenvalue ratio. Curvature = λ3 / (λ1+λ2+λ3); This index quantifies the degree of surface variation at each point. Higher values ​​indicate potential edges or corners, and points with high curvature values ​​(>0.1) are marked as potential boundary regions between discontinuous surfaces, which are given special consideration in subsequent processing stages.

[0036] Furthermore, based on the normal vector field, the number of clusters of dominant groups of discontinuities is automatically determined by the density peak clustering algorithm, resulting in cluster centers including: The normal vectors in the normal vector field are converted into spherical coordinates, and the similarity between the vectors is calculated using the dot product. The local density of each normal vector in the normal vector field is calculated based on the spherical coordinate system representation. Calculate the minimum distance from each normal vector to the nearest normal vector with higher local density; The γ value of each normal vector is determined based on the product of the local density and the minimum distance, where the γ value is used to determine whether the corresponding point is a cluster center; The γ values ​​are sorted from high to low, and the cluster centers are obtained based on the sorted γ values.

[0037] Specifically, in step S3, the automatic determination of the number of clusters for density peak clustering includes: S31, convert the normal vector to spherical coordinates for analysis, and use the dot product to calculate the similarity between vectors: ; Where, θ ij Represents the normal vector n i and n j The angular separation between them, S(n) i , n j Let θ be a function of the similarity between two normal vectors, and let θ be the angle between the two normal vectors. S32, Calculate local density : ; Where, d ij =arccos(n i ·n j ) represents the angular distance between the normal vectors, d c The cutoff distance is determined as the average distance r (where r is 2% of the total number of points), χ(x) = 1 (if x < 0) and 0 (otherwise). This formula ensures that the cutoff distance is only valid at angular distance d. c Only points within the local area contribute to the local density metric; S33, calculate the minimum distance δ to any point with higher density. i : δ i =min(d ij ), j: > ; Where, δ i This represents the distance from point i to any point with a higher density (i.e., satisfying the condition that...). > The minimum distance between the points ( ), Let j be the local density. Let be the local density at point i.

[0038] For the point with the highest density, δ i Let max(d) ij To ensure proper boundary handling; S34, by calculating γ i = The values ​​of γ are then selected, and the points with the highest γ values ​​in the top 10% are chosen as cluster centers to identify potential cluster centers. γ i = ; Where, γ i A comprehensive decision value is calculated for each data point.

[0039] The threshold was optimized through systematic testing on multiple datasets to balance undersegmentation and oversegmentation of discontinuous surfaces.

[0040] Further, based on the cluster centers and the normal vector field, spatial-directional ensemble clustering is performed on the 3D point cloud data to obtain initial clustering results, including: Calculate the weighted distance from each point to each cluster center, where the weighted distance is obtained by weighted summation of spatial distance and normal vector similarity terms; Each point is assigned to the nearest cluster center based on the weighted distance; Iteratively update the cluster centers until the change in the cluster centers is less than the convergence threshold, and obtain the initial clustering result.

[0041] Specifically, in step S4, the adaptive spatial-orientation integrated clustering includes: S41, Weighted distance calculation: For each vector with normal vector n p The point is calculated to have a normal vector n. ci Cluster center c i The weighted distance D(p, c) i ): ; Where w s and w n These are the weighting factors for spatial distance and normal vector similarity, d max Normalization is provided. The optimal weight w was determined through system parameter optimization across multiple datasets. s = 0.3 and w n = 0.7, emphasizing directional consistency while maintaining spatial coherence; S42, Parallel Normal Vector Criterion Implementation: Implements the parallel normal vector criterion to handle discontinuous surfaces with similar directions but opposite normal vector directions. ; This criterion identifies the geometric equivalence of planes with opposite normal vectors, preventing the artificial segmentation in traditional clustering methods that treats relative normal vectors as different clusters. S43, Iterative Clustering Update: Each point is assigned to its nearest cluster center using a weighted distance metric, and then the cluster centers are updated based on the average coordinates and average normalized normal vector of the assigned points. This process continues until convergence (typically 15-25 iterations), defined as the maximum change in cluster centers being less than ε = 0.01.

[0042] Furthermore, the initial clustering results are optimized by performing boundary optimization using a region growing algorithm, resulting in optimized clustering results including: Based on the initial clustering results, high-confidence points are selected from each cluster as seed points; For each seed point's neighboring points, calculate the normal vector consistency weight and spatial distance weight between the neighboring points and the seed point; Based on the normal vector consistency weight and the spatial distance weight, the final clustering assignment of each neighboring point is determined through a weighted voting mechanism to obtain the optimized clustering result.

[0043] Specifically, in step S4, the region growth optimization includes: S44, selects high-confidence seed points from each cluster based on local density and distance to the cluster center; S45, using KD-tree search to identify radius r max The spatial domain of each seed point within; S46, Normal Vector Consistency Criterion: It must be proven that the normal vector is consistent with the seed point, determined by the angle threshold θ. max definition: ; Where θ max = π / 12 (15 degrees), providing an appropriate balance between inclusiveness and accuracy based on empirical testing; S47, Spatial Continuity Criterion: Spatial continuity must be maintained within the dynamic adjustment radius. r max = β×d avg ; Where d avg It is the average local point spacing, and β = 5 is a scaling factor determined experimentally to adapt to different point densities in the dataset; S48, Weighted Voting Mechanism Implementation: By considering the collective influence of neighborhood point sets instead of the binary inclusion criterion, a weighted voting mechanism is implemented to improve boundary accuracy. ; Among them, W(p, C) k Let p be the point and C be the cluster C. kThe association weights between points are N(p), which is the neighborhood of point p, and w(p, q) are the weights based on distance and normal vector similarity, where I(q∈C) is the weight of the association between points. k ) is when point q belongs to cluster C k An indicator function that equals 1 when the value is true and 0 otherwise. The weighting function incorporates spatial and directional factors: ; in, Let be the variance parameter of the spatial distance. The variance parameter of the normal vector direction gives higher influence to domains that are both spatially close and have similar normal vector directions, thus enhancing the robustness of boundary region assignment. S49, Final Cluster Assignment: The final cluster assignment is determined by maximizing weighted voting. ; Where C(p) is the cluster category corresponding to point p (i.e., the cluster to which point p belongs).

[0044] Furthermore, the optimized clustering results are decomposed at multiple scales to obtain the attitude parameters of the rock mass discontinuities, including: For each cluster in the optimized clustering results, a hierarchical density clustering algorithm is used to perform sub-clustering decomposition, identify fine-scale substructures in the main discontinuities, and obtain a multi-scale clustering structure. For each cluster and sub-cluster in the multi-scale clustering structure, a plane fitting is performed using the random sample consensus algorithm to obtain the fitted plane equation; Based on the normal vector of the fitted plane equation, the dip and dip angle parameters are calculated to obtain the attitude parameters of the rock mass discontinuity surface.

[0045] Specifically, in step S5, the multi-scale analysis and parameter extraction include: S51, Hierarchical Sub-clustering: The Hierarchical Density Clustering (HDBSCAN) algorithm is applied to each identified discontinuity to identify substructures within the main discontinuity group. Unlike traditional clustering methods, HDBSCAN constructs a distance-based minimum spanning tree and selects stable clusters across multiple density levels, making it well-suited for detecting fine-scale variations within a set of structures. S52, RANSAC Plane Fitting: For each identified cluster and sub-cluster, a robust plane fitting is performed using the Random Sample Consensus (RANSAC) algorithm. RANSAC iteratively samples three points to define a candidate plane and calculates the distance d from each point to the plane. i : ; Where, a is the coefficient of the plane equation (x-component of the plane normal vector); b is the coefficient of the plane equation (y-component of the plane normal vector); c is the coefficient of the plane equation (z-component of the plane normal vector); d is the constant term of the plane equation; x i The x-coordinate of a point in space; the y-coordinate i The z-coordinate is the y-coordinate of a point in space; i Let z be the z-coordinate of the spatial point, and select the plane that maximizes the interior point count. Set the distance threshold to 1.5 times the average point spacing, and the number of iterations is adaptively determined according to the cluster size (usually 1000-5000 iterations). S53, Calculation of geological attitude parameters: For each fitted plane a x +b y +c z +d=0, calculate geological occurrence parameters according to standard practice: Tendency = ; Inclination angle = ; Among them, u y For the unit normal vector y-component, u x For the x-component of the unit normal vector, u z For the z-component of the unit normal vector, (n x ,n y , n z ) represents the unit normal vector of the fitted plane.

[0046] The following examples will illustrate this embodiment in detail: like Figure 1 As shown, the embodiments of the present invention demonstrate a complete processing flow including five stages: data acquisition and preprocessing, normal vector estimation, density peak clustering to automatically determine the number of clusters, adaptive spatial-directional ensemble clustering, and multi-scale analysis and parameter extraction.

[0047] Phase 1 - Data Acquisition and Preprocessing: Raw point cloud data is preprocessed through voxel downsampling, statistical outlier removal, and quality assessment to ensure geometric integrity while reducing computational burden.

[0048] Phase 2 - Normal Vector Estimation: Robust normal vector computation is achieved through adaptive neighborhood selection based on KD-tree structure (via PCA), and consistent directional propagation is performed through minimum spanning tree.

[0049] Phase 3 - DPCA clustering automatically determines the number of clusters: It analyzes the normal vector density distribution in three-dimensional space and automatically identifies the number of dominant groups on the best discontinuity surface by detecting density peaks, thereby eliminating the need for subjective parameter specification.

[0050] Phase 4 - Adaptive Clustering and Spatial-Oriented Integration: Enhanced K-means clustering uses a weighted distance metric to jointly consider the similarity of spatial coordinates and normal vectors, and then improves the boundary partitioning accuracy through region growing optimization with a weighted voting mechanism.

[0051] Phase 5 - Multiscale Analysis and Parameter Extraction: Hierarchical HDBSCAN decomposition identifies fine-scale substructures within the main discontinuities, followed by RANSAC plane fitting and SVD optimization to extract geological attitude parameters (dip and dip angle) for engineering applications.

[0052] This embodiment includes the following steps: S1, Data Acquisition and Preprocessing: Two complementary datasets were used to validate the framework under different geological environments and data acquisition methods. Dataset A (Case A) includes Rockbench benchmark data: 1,515,722 points acquired by an Optech LiDAR scanner with a resolution of less than 2 cm. Dataset B (Case B) represents a more challenging scenario: 5,953,005 points of sedimentary rock slopes captured by UAV photogrammetry.

[0053] Preprocessing included voxel downsampling and outlier removal. Statistical outlier removal based on local neighborhood analysis was implemented to improve data quality while preserving geometric integrity. For each point, the average distance to the m nearest neighbors (m=30) was calculated, and points exceeding a threshold μ+2σ (where μ is the average distance and σ is the standard deviation) were classified as outliers and removed. This threshold adapts to density variations in the dataset, effectively eliminating measurement noise while preserving the surface characteristics of discontinuous surfaces.

[0054] S2, Normal Vector Estimation: Accurate normal vector estimation provides the foundation for subsequent clustering and segmentation operations. Unlike fixed-radius methods, which may exhibit inconsistent performance under different point densities, an adaptive KD-tree method is implemented to adapt to local density variations.

[0055] like Figure 2 The diagram shown is a geometric representation of the relationship between the normal vector and the angle. For each point p, a local neighborhood is constructed using a KD-tree data structure. The covariance matrix C is calculated as follows: ; in p is the centroid of the neighborhood points. i Let represent each neighborhood point, and k be the number of neighborhood points. Eigenvalue decomposition produces three eigenvalues ​​(λ1 ≥ λ2 ≥ λ3) and corresponding eigenvectors (v1, v2, v3). The eigenvector v3 corresponding to the smallest eigenvalue λ3 represents the estimated normal vector because it is aligned with the direction of minimum variance of the local point distribution.

[0056] To address the ambiguity of normal vector direction (the vector may point inside or outside the surface), a consistent direction procedure based on minimum spanning tree (MST) propagation is implemented. This method ensures that normal vectors follow a consistent pattern throughout the point cloud, facilitating more accurate clustering in subsequent stages.

[0057] Quantize the local surface curvature at each point using the eigenvalue ratio: Curvature = λ3 / (λ1+λ2+λ3); This metric quantifies the degree of surface variation at each point, with higher values ​​indicating potential edges or corners. Points with high curvature values ​​(>0.1) are marked as potential boundary regions between discontinuous surfaces and receive special consideration in subsequent processing stages.

[0058] By conducting system parameter sensitivity tests on three datasets with different point densities, the optimal neighborhood size k=60 was determined, providing the best balance between noise reduction and feature preservation. Comparison of normal vector estimation results using different neighborhood parameters (k=20, k=60, k=100) shows that k=60 offers superior performance in terms of accuracy and computational efficiency, with a smoother normal vector field and better preservation of discontinuous surface boundaries.

[0059] S3, Density Peak Clustering automatically determines the number of clusters: The fundamental innovation of this invention is the use of the Density Peak Clustering (DPCA) algorithm to automatically determine the optimal number of dominant groups on discontinuities. Unlike traditional methods that require manually specifying the number of clusters, DPCA identifies cluster centers based on density distribution patterns in the normal vector space, eliminating subjective bias in parameter selection.

[0060] First, the normal vectors are converted to spherical coordinates for analysis, and their dot product is used to calculate the similarity between the vectors: ; Where θ ij Represents the normal vector n i and n j The angle between them.

[0061] DPCA identifies cluster centers based on two key metrics: local density. and the minimum distance δ to higher density points i The local density is calculated as follows: ; Where d ij = arccos(n i ·n j ) represents the angular distance between the normal vectors, d cIt is the cutoff distance determined as the average r-distance (r is 2% of the total number of points), χ(x) = 1 (if x < 0) and 0 (otherwise). This formula ensures that only the angular distance d is considered as the cutoff distance. c Only points within the local area contribute to the local density metric.

[0062] Minimum distance δ to any point with higher density i The calculation is as follows: δ i = min(d ij ), j: > ; For the point with the highest density, δ i Let max(d) ij To ensure proper boundary handling.

[0063] like Figure 3 The diagram shown is a schematic of the DPCA decision diagram. Potential cluster centers are determined by calculating γ. i = The values ​​of γ are then used to select the top 10% of cluster centers for identification. γ i = ; The threshold was optimized through systematic testing on multiple datasets to balance undersegmentation and oversegmentation. The method effectively identifies clusters across different scales, capturing dominant and secondary structural features without requiring prior knowledge.

[0064] S4, Adaptive Spatial-Oriented Integrated Clustering: Following cluster center identification, a unique enhanced K-means algorithm combining spatial approximation and normal vector similarity is implemented to optimize discontinuous surface grouping. Traditional methods treat directional clustering and spatial continuity as independent processes, which can lead to spatially dispersed results or geometrically inconsistent grouping.

[0065] For each having a normal vector n p For point p, calculate the normal vector n. ci Cluster center c i The weighted distance is calculated as follows: ; Where w s and w n These are the weighting factors for spatial distance and normal vector similarity, d max Normalization is provided. The optimal weight w was determined through system parameter optimization across multiple datasets. s = 0.3 and w n = 0.7, emphasizing directional consistency while maintaining spatial coherence.

[0066] A key innovation of this method is the implementation of a parallel normal vector criterion to handle discontinuities with similar orientations but opposite normal vector directions: ; This criterion identifies the geometric equivalence of planes with opposite normal vectors, preventing the artificial segmentation of relative normal vectors as different clusters in traditional clustering methods.

[0067] The algorithm iteratively assigns each point to the nearest cluster center using a weighted distance metric, and then updates the cluster centers based on the average coordinates and average normalized normal vector of the assigned points. This process continues until convergence (typically 15-25 iterations), defined as the maximum change in cluster centers being less than ε = 0.01°.

[0068] The region growing algorithm optimizes clustering results by ensuring consistency of normal vectors and spatial continuity on discontinuous face sets, thus addressing boundary uncertainties. The process begins with high-confidence seed points selected for each cluster, based on local density and distance metrics to the cluster centers.

[0069] For each seed point, use a KD-tree search to identify the radius r. max The spatial domain within. First, it must be proven that the normal vector is consistent with the seed point, by the angle threshold θ. max definition: ; Where θ max =π / 12 (15 degrees), providing an appropriate balance between inclusiveness and accuracy based on empirical testing.

[0070] Secondly, spatial continuity must be maintained within the dynamically adjusted radius: r max =β×d avg ; Where d avg It is the average local point spacing, and β=5 is a scaling factor determined experimentally to adapt to different point densities in the dataset.

[0071] A significant enhancement of this method is the implementation of a weighted voting mechanism, which improves boundary accuracy by considering the collective influence of the neighborhood point set rather than the binary inclusion criterion. ; Where N(p) is the neighborhood of point p, w(p, q) is the weight based on distance and normal vector similarity, and I(q∈C) k ) is when point q belongs to cluster C k An indicator function that equals 1 when the value is true and 0 otherwise. The weighting function incorporates spatial and directional factors: ; This formula gives greater influence to domains that are both spatially close and have similar normal vector directions, enhancing the robustness of boundary region assignment.

[0072] The final cluster assignment is determined by the following: ; S5, Multiscale Analysis and Parameter Extraction: To identify substructures within a group of major discontinuities, a hierarchical density clustering algorithm (HDBSCAN) is applied to each identified discontinuity. Unlike traditional clustering methods that apply a uniform density threshold, HDBSCAN constructs a distance-based minimum spanning tree and selects stable clusters across multiple density levels, making it well-suited for detecting fine-scale variations within groups of structures.

[0073] Key parameters were optimized through sensitivity analysis: minimum cluster size = 100 points (approximately 0.5% of the typical discontinuity size), minimum number of samples = 10 (determined through sensitivity analysis to balance noise tolerance and cluster detection), and cluster selection epsilon = 0.5 (optimized through testing on multiple datasets). These values ​​ensure the detection of meaningful substructures while avoiding over-fragmentation that would otherwise complicate interpretation.

[0074] The final stage uses Random Sample Consensus (RANSAC) to perform robust plane fitting on each identified cluster and sub-cluster. RANSAC iteratively samples three points to define candidate planes and calculates the distance from each point to the plane: ; The plane that maximizes the interior point count is selected. The distance threshold is set to 1.5 times the average point spacing, and the number of iterations is adaptively determined based on the cluster size (usually 1000-5000 iterations).

[0075] For each fitted plane a x +b y +c z +d=0, calculate geological occurrence parameters according to standard practice: Tendency = arctan(n) y / n x )mod360°; Inclination angle = arccos(n) z ); Where (n x , n y , n z ) represents the unit normal vector of the fitted plane. Example 1

[0076] This embodiment provides a validation example of analyzing dataset A. For example... Figure 5 As shown, the method of this invention was applied to analyze dataset A, automatically identifying four main discontinuities (colored as follows: J1-red, J2-blue, J3-green, J4-yellow). The density peak clustering algorithm successfully determined the optimal number of clusters by analyzing the distribution of normal vectors in three-dimensional space, eliminating the inherent subjective interpretation in traditional methods.

[0077] Figure 3 For the decision graph of the density peak clustering algorithm, the relationship between the local density (ρ) and minimum distance (δ) of all normal vector data points is plotted. The point size is proportional to the γ value (γ = ρ × δ), showing the characteristic distribution pattern. Most data points are concentrated in the lower left region, with both low local density and small minimum distance, representing points inside the cluster. In contrast, four distinct cluster centers (C1-C4) appear as outliers with high γ values, demonstrating the characteristic distribution: - C1(ρ=16, δ=0.206): Medium density and relatively high minimum distance, indicating a well-separated cluster that represents the main structural features.

[0078] - C2(ρ=6, δ=0.347): Located in the region of highest minimum distance, representing isolated clusters with moderate internal density.

[0079] - C3 (ρ=20, δ=0.059): The region with the highest density, moderate separation, and characterized by the most prominent discontinuities.

[0080] - C4(ρ=8, δ=0.323): High minimum distance and medium density, indicating another well-separated cluster.

[0081] The four identified discontinuities exhibit distinct directions with clear geological significance: - J1 (dip / dipping angle: 246.12° / 35.01°): Represents a moderately dipping bedding plane consistent with the regional structural trend.

[0082] - J2 (172.42° / 81.78°) and J3 (136.46° / 76.04°): represent steeply dipping joints formed under different stress states, which may reflect tectonic compression from multiple deformation events.

[0083] - J4 (93.06° / 47.79°): Represents a moderately inclined joint that reflects local stress variations.

[0084] These results are consistent with the known geological features of Case A, where orthogonal joint systems developed during regional uplift and erosion are frequently observed.

[0085] Hierarchical HDBSCAN sub-clustering identified 13 distinct substructures, providing insights into local structural variations in impact hydrogeological behavior and hydrogeological characteristics. Figure 6 As shown in (a)-(c), the identified substructures exhibit slight but systematic variations in orientation (standard deviation: dip 2.78°, dip angle 1.92°), which may reflect local stress variations or sequential formation processes in the geological history of the rock mass. These fine-scale variations are crucial for comprehensive rock mass characterization because they affect the distribution of discontinuity spacing, block size estimation, and local stability conditions. Example 2

[0086] This embodiment provides a verification example for the analysis of dataset B. For example... Figure 7 and Figure 8 As shown in (a)-(c), Case B presents a more challenging scenario with a significantly larger dataset (5.95 million points), representing a sedimentary rock slope with complex structural geometry. To manage computational constraints while maintaining analytical accuracy, a localized sampling method was implemented, including spatial segmentation into 1m... 3 Voxels, density sampling based on geometric importance (retaining 10% of points, prioritizing high curvature or boundary points), and normal vector consistency verification to ensure preservation of key structural features.

[0087] Applying the framework of this embodiment, three dominant groups of major discontinuities were automatically identified, each group exhibiting different spatial patterns and geometric features reflecting the potential geological structure of the sedimentary body. DPCA clusters were automatically determined through low-density (ρ) and minimum distance (δ) analysis, and cluster centers with high gamma values ​​were extracted using optimized thresholds (ρ>12, δ>0.1). This automatic determination represents a significant improvement over traditional methods that require subjective interpretation through stereo projection.

[0088] The three identified dominant groups exhibit distinct directional patterns, which have clear geological significance: - Dominant group J1 (76.14° / 61.85°): Shows a highly concentrated distribution in the normal vector density map, indicating strong directional consistency and surface integrity. This dominant group may correspond to oblique shear joints formed under regional tectonic stress, and the moderate dip angle is conducive to the formation of potential groundwater flow paths along discontinuous surfaces.

[0089] - Dominant group J2 (261.85° / 18.47°): Characterized by extensive lateral continuity in platy structure, typical of sedimentary bedding or stratified joints. Low dip angle (<20°) indicates stratification control in the rock mass structure, potentially representing the main weak surface of the planar sliding instability mechanism. Southwest dip suggests favorable slope stability in eastward excavation, but potential instability in westward slopes.

[0090] - Dominant group J3 (182.66° / 80.51°): Displays a spatial pattern of linear alignment and marginal clustering, indicating tectonic compression characteristics. The near-vertical orientation (dip > 80°) and southward dip suggest formation as tensional fractures or as conjugate shear fractures during regional extension. This group creates potential boundaries for overturning or wedge-shaped instability mechanisms when intersecting with other discontinuities, and is particularly critical for steep slope configurations.

[0091] Figure 4 This is a parameter sensitivity analysis diagram for the present invention. The impact of various key parameters on the results was evaluated through systematic parameter sensitivity testing. Figure 4 The results show that the number of clusters identified remains stable at the correct value under different angle threshold conditions (6°-24°), indicating the robustness of the determined cluster number. When the threshold is in the range of 10°-20°, the directional consistency remains at a high level (>0.75), indicating that this is the optimal parameter range.

[0092] For the region growing algorithm, the effect of segmentation quality on the angle threshold (θ) was also analyzed. max The sensitivity to the neighborhood scaling factor (β) and the angle threshold exhibits a clear optimal value of approximately θ. max = 15° (π / 12), with relatively stable performance in the range of 10°–20°. Values ​​below 10° lead to overly strict normal vector matching, causing fragmentation of continuous planes; values ​​above 20° lead to excessively inclusive grouping, merging different adjacent discontinuities. The selected value of 15° provides a suitable balance, consistent with the typical orientation measurement uncertainties in geological fieldwork.

[0093] The neighborhood scaling factor β shows less sensitivity, with values ​​between 3 and 7 producing similar results (F1 score variation < 0.87). The chosen value β = 5 ensures sufficient neighborhood coverage for reliable voting while avoiding excessive computational burden from very large neighborhoods. This stability indicates that the weighted voting mechanism is robust to moderate variations in the spatial search radius.

[0094] The computational performance analysis of the method in this invention shows that, for different point cloud sizes, the processing time increases with the increase of point cloud size according to the complexity of O(nlogn), which is consistent with the theoretical complexity analysis. Normal vector estimation and DPCA clustering are the most computationally intensive stages, accounting for approximately 80% of the total processing time. Region growing optimization contributes 17.4%, while adaptive clustering and structure analysis show relatively low computational costs (total <5%).

[0095] Benchmark tests compared to existing methods show significant performance improvements. For the Case A dataset (881,552 points), our method completed the analysis in 169.52 seconds, a reduction of 60% to 72% compared to existing methods. This improvement is primarily due to: an optimized KD-tree implementation for spatial operations with adaptive neighborhood size adjustment; an efficient DPCA algorithm for automatic cluster determination without iterative optimization; and a memory-efficient data structure that minimizes redundant computation.

[0096] These performance improvements make the framework particularly well-suited for handling large-scale point cloud datasets common in modern remote sensing applications. Computational efficiency enables practical field applications required for engineering decision-making, supporting time-sensitive scenarios such as emergency slope stability assessments after earthquakes or rapid characterization for excavation planning.

[0097] This embodiment also provides an automatic identification system for rock mass discontinuities based on density peak clustering, including: a data acquisition module, a normal vector estimation module, an adaptive clustering module, a region growth optimization module, and a multi-scale feature extraction and optimization module; The data acquisition module is used to acquire three-dimensional point cloud data of the rock mass; The normal vector estimation module is used to obtain the normal vector field based on the three-dimensional point cloud data; The adaptive clustering module is used to automatically determine the number of clusters for discontinuous surfaces based on the normal vector field using the density peak clustering algorithm, and obtain the cluster centers. The region growing optimization module is used to perform spatial-directional ensemble clustering on the 3D point cloud data based on the cluster centers and the normal vector field to obtain an initial clustering result; and to perform boundary optimization on the initial clustering result using a region growing algorithm to obtain an optimized clustering result. The multi-scale feature extraction and optimization module is used to perform multi-scale decomposition on the optimized clustering results to obtain the attitude parameters of the rock mass discontinuity surface.

[0098] Specifically, the system in this embodiment includes: Data Acquisition and Preprocessing Module: Acquires 3D point cloud data of the rock mass and performs voxel downsampling and outlier removal on the raw data. Voxel downsampling sets the voxel size based on the point cloud density and target resolution, dividing the original point cloud space into a uniformly sized voxel grid. Points in each non-empty voxel are averaged. Outlier removal is based on local neighborhood analysis, calculating the average distance from each point to its m nearest neighbors (m=30). Points exceeding the threshold μ+2σ are classified as outliers and removed.

[0099] Normal vector estimation module: Uses KD tree data structure to build local neighborhood with adaptive search criteria for each point, determines normal vector by calculating covariance matrix and its eigenvalue decomposition, handles consistency of normal vector direction by minimum spanning tree (MST), and calculates surface curvature to mark potential boundary regions.

[0100] Density Peak Clustering Module: Analyzes the normal vector density distribution in 3D space, calculates the local density (ρ) of each point and the minimum distance (δ) to higher density points, and automatically determines the optimal number of discontinuous surfaces by calculating the γ value (γ = ρ × δ) and selecting points with high γ values ​​as cluster centers.

[0101] Adaptive clustering module: Employs enhanced K-means clustering with weighted distance metric, jointly considering spatial coordinates and normal vector similarity, and simultaneously implementing the parallel normal vector criterion to handle discontinuous surfaces with similar directions but opposite normal vector directions, iteratively updating cluster centers until convergence.

[0102] The region growth optimization module selects high-confidence seed points based on local density and distance to the cluster center. It constructs a spatial neighborhood for each seed point and applies the normal vector consistency and spatial continuity criteria. It considers the collective influence of the neighborhood point set through a weighted voting mechanism and finally determines the optimal clustering affiliation for each point.

[0103] Multi-scale analysis module: The hierarchical density clustering algorithm (HDBSCAN) is used to perform sub-clustering on each major discontinuity, identify fine-scale structural changes, construct distance-based minimum spanning trees, and select stable clusters across multiple density levels.

[0104] Parameter extraction module: Apply RANSAC plane fitting to each cluster and sub-cluster, iteratively sample to define candidate planes and calculate the distance from points to the planes, select the plane that maximizes the interior points, and then calculate geological attitude parameters (dip and dip angle) to provide standardized structural data for engineering applications.

[0105] The modules are tightly integrated through data flow, forming a complete automated processing flow. The system interface design allows users to input point cloud data in various formats, including standard formats such as .las, .ply, and .pcd, and provides result visualization and export functions, supporting seamless integration with CAD and GIS software.

[0106] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. An automatic identification method for rock mass discontinuities based on density peak clustering, characterized in that, include: Acquire 3D point cloud data of the rock mass; Based on the aforementioned 3D point cloud data, obtain the normal vector field; Based on the normal vector field, the number of clusters of the dominant group of discontinuities is automatically determined by the density peak clustering algorithm, and the cluster centers are obtained. Based on the cluster centers and the normal vector field, spatial-directional ensemble clustering is performed on the 3D point cloud data to obtain the initial clustering result; The initial clustering results are optimized by performing boundary optimization using a region growing algorithm to obtain optimized clustering results. The optimized clustering results are decomposed at multiple scales to obtain the attitude parameters of the rock mass discontinuities.

2. The method for automatic identification of rock mass discontinuities based on density peak clustering according to claim 1, characterized in that, Based on the aforementioned 3D point cloud data, obtaining the normal vector field includes: The three-dimensional point cloud data is preprocessed to obtain preprocessed point cloud data; Based on the preprocessed point cloud data, the normal vector field is obtained by calculating the normal vector through principal component analysis.

3. The method for automatic identification of rock mass discontinuities based on density peak clustering according to claim 2, characterized in that, The three-dimensional point cloud data is preprocessed to obtain preprocessed point cloud data, including: The three-dimensional point cloud data is downsampled using voxels to obtain the downsampled point cloud; Outlier points in the downsampled point cloud are identified and removed to obtain the preprocessed point cloud data.

4. The method for automatic identification of rock mass discontinuities based on density peak clustering according to claim 3, characterized in that, Based on the preprocessed point cloud data, normal vectors are calculated using principal component analysis to obtain the normal vector field, which includes: An adaptive neighborhood is constructed for each point in the preprocessed point cloud data based on the KD tree structure. Principal component analysis is performed on the points within the adaptive neighborhood to obtain the covariance matrix; The covariance matrix is ​​decomposed into eigenvalues, and the normal vector of each point is determined based on the eigenvector corresponding to the smallest eigenvalue, wherein the normal vector is obtained based on the direction of the smallest eigenvalue. The surface curvature of each point is calculated based on the magnitude of the minimum eigenvalue, and points with curvature higher than the threshold are marked as potential boundary points. The normal vector field is obtained based on the normal vector and the surface curvature.

5. The method for automatic identification of rock mass discontinuities based on density peak clustering according to claim 1, characterized in that, Based on the normal vector field, the number of clusters in the dominant group of discontinuities is automatically determined using the density peak clustering algorithm, resulting in cluster centers including: The normal vectors in the normal vector field are converted into spherical coordinates, and the similarity between the vectors is calculated using the dot product. The local density of each normal vector in the normal vector field is calculated based on the spherical coordinate system representation. Calculate the minimum distance from each normal vector to the nearest normal vector with higher local density; The γ value of each normal vector is determined based on the product of the local density and the minimum distance, where the γ value is used to determine whether the corresponding point is a cluster center; The γ values ​​are sorted from high to low, and the cluster centers are obtained based on the sorted γ values.

6. The method for automatic identification of rock mass discontinuities based on density peak clustering according to claim 1, characterized in that, Based on the cluster centers and the normal vector field, spatial-directional ensemble clustering is performed on the 3D point cloud data to obtain the initial clustering results, including: Calculate the weighted distance from each point to each cluster center, where the weighted distance is obtained by weighted summation of spatial distance and normal vector similarity terms; Each point is assigned to the nearest cluster center based on the weighted distance; Iteratively update the cluster centers until the change in the cluster centers is less than the convergence threshold, and obtain the initial clustering result.

7. The method for automatic identification of rock mass discontinuities based on density peak clustering according to claim 6, characterized in that, The initial clustering results are optimized by performing boundary optimization using a region growing algorithm, resulting in the following optimized clustering results: Based on the initial clustering results, high-confidence points are selected from each cluster as seed points; For each seed point's neighboring points, calculate the normal vector consistency weight and spatial distance weight between the neighboring points and the seed point; Based on the normal vector consistency weight and the spatial distance weight, the final clustering assignment of each neighboring point is determined through a weighted voting mechanism to obtain the optimized clustering result.

8. The method for automatic identification of rock mass discontinuities based on density peak clustering according to claim 7, characterized in that, The optimized clustering results are decomposed at multiple scales to obtain the attitude parameters of the rock mass discontinuities, including: For each cluster in the optimized clustering results, a hierarchical density clustering algorithm is used to perform sub-clustering decomposition, identify fine-scale substructures in the main discontinuities, and obtain a multi-scale clustering structure. For each cluster and sub-cluster in the multi-scale clustering structure, a plane fitting is performed using the random sample consensus algorithm to obtain the fitted plane equation; Based on the normal vector of the fitted plane equation, the dip and dip angle parameters are calculated to obtain the attitude parameters of the rock mass discontinuity surface.

9. An automatic identification system for rock mass discontinuities based on density peak clustering, used to implement the method as described in any one of claims 1-8, characterized in that, include: The module includes a data acquisition module, a normal vector estimation module, an adaptive clustering module, a region growing optimization module, and a multi-scale feature extraction and optimization module. The data acquisition module is used to acquire three-dimensional point cloud data of the rock mass; The normal vector estimation module is used to obtain the normal vector field based on the three-dimensional point cloud data; The adaptive clustering module is used to automatically determine the number of clusters of dominant groups of discontinuities based on the normal vector field and the density peak clustering algorithm, thereby obtaining the cluster centers. The region growing optimization module is used to perform spatial-directional ensemble clustering on the 3D point cloud data based on the cluster centers and the normal vector field to obtain an initial clustering result; and to perform boundary optimization on the initial clustering result using a region growing algorithm to obtain an optimized clustering result. The multi-scale feature extraction and optimization module is used to perform multi-scale decomposition on the optimized clustering results to obtain the attitude parameters of the rock mass discontinuity surface.

Images