Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Integration of a Mobile Laser Scanning System with a Forest Harvester for Accurate Localization and Tree Stem Measurements
Previous Article in Journal
Evaluation of Urban Microscopic Nighttime Light Environment Based on the Coupling Observation of Remote Sensing and UAV Observation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Characterization of Complex Rock Mass Discontinuities from LiDAR Point Clouds

School of Geography and Information Engineering, China University of Geosciences, Wuhan 430074, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2024, 16(17), 3291; https://doi.org/10.3390/rs16173291
Submission received: 2 August 2024 / Revised: 1 September 2024 / Accepted: 3 September 2024 / Published: 4 September 2024

Abstract

:
The distribution and development of rock mass discontinuities in 3D space control the deformation and failure characteristics of the rock mass, which in turn affect the strength, permeability, and stability of rock masses. Therefore, it is essential to accurately and efficiently characterize these discontinuities. Light Detection and Ranging (LiDAR) now allows for fast and precise 3D data collection, which supports the creation of new methods for characterizing rock mass discontinuities. However, uneven density distribution and local surface undulations can limit the accuracy of discontinuity characterization. To address this, we propose a method for characterizing complex rock mass discontinuities based on laser point cloud data. This method is capable of processing datasets with varying densities and can reduce over-segmentation in non-planar areas. The suggested approach involves a five-stage process that includes: (1) adaptive resampling of point cloud data based on density comparison; (2) normal vector calculation using Principal Component Analysis (PCA); (3) identifying non-planar areas using a watershed-like algorithm, and determine the main discontinuity sets using Multi-threshold Mean Shift (MTMS); (4) identify single discontinuity clusters using Density-Based Spatial Clustering of Applications with Noise (DBSCAN); (5) fitting discontinuity planes with Random Sample Consensus (RANSAC) and determining discontinuity orientations using analytic geometry. This method was applied to three rock slope datasets and compared with previous research results and manual measurement results. The results indicate that this method can effectively reduce over-segmentation and the characterization results have high accuracy.

1. Introduction

Rock slopes are geological entities with unique structural characteristics, displaying significant discontinuities, heterogeneity, and anisotropy. Due to complex tectonic movements at different historical periods, surface weathering, and stress release processes, rock masses often develop various discontinuities of different scales and properties, such as fractures, joints, bedding planes, and faults. The distribution and evolution of discontinuities in 3D space govern the deformation and failure characteristics of rock masses. It is one of the important factors that affect the strength, permeability, and stability of rock masses [1,2]. The characterization of rock mass discontinuities involves identifying the location, size, and some typical parameters of the discontinuities. Several geometric parameters can quantitatively describe rock discontinuities, such as orientation, roughness, spacing, persistence, and more [3]. Among these, the orientation (dip angle and dip direction) is particularly critical for analyzing the stability of jointed rock masses [4]. For instance, orientation is a critical input parameter in block theory [5], rock mass classification [6], and the use of the limited equilibrium method [7]. Hence, it is extremely important to accurately map these discontinuities and determine the orientation of rock discontinuities in rock engineering [1,8].
Traditionally, geologists use tools such as geological compasses and tape measures to manually measure the discontinuity plane of rock masses using methods like scanline surveys or window sampling to collect their geometric information [9,10]. Although this method is straightforward and user-friendly, it has several significant drawbacks: (1) The measurement process relies entirely on manual execution, which results in low efficiency [11]; (2) The working conditions are easily constrained by terrain and environmental factors [12]; (3) In complex terrains, especially in fragmented rock areas and steep slopes, surveyors may be unable to access specific areas for measurement [13]; (4) Investigators need to have relevant professional backgrounds, and the back of these makes it easy for human errors to occur [14,15]. These limitations constrain the ability to conduct an in-depth analysis of the structural characteristics of rock slopes.
In modern geotechnical engineering and geological exploration, Light Detection and Ranging (LiDAR) and digital photogrammetry are widely used remote sensing techniques for acquiring high-quality 3D point cloud data, capable of capturing high-resolution 3D information [16,17]. The origins of 3D laser scanning technology can be traced back to the mid-20th century. This technology uses 3D laser scanners to rapidly scan and measure target objects, directly capturing the distance, direction, and intensity values of the reflected laser from the object. Due to its ability to collect real-time, high-precision 3D point clouds of rock surfaces without physical contact, this technology eliminates the subjective factors associated with manual contact measurements. Therefore, compared to traditional manual measurement methods, 3D laser scanning offers a more objective and accurate means of data acquisition. These techniques have seen rapid advancement in geological and geotechnical research fields [18]. Rock mass discontinuities can be categorized into planar and linear types. Linear discontinuities primarily manifest as fractures, with their spatial distribution occurring within the rock mass. On the other hand, planar discontinuities are exposed fault surfaces that appear on the rock mass surface after sliding has occurred. This study focuses on the characterization of planar discontinuities. In recent years, many researchers have utilized 3D point cloud data to analyze discontinuities, including discontinuity recognition [19,20,21,22], slope stability analysis [23,24], discontinuity and rapid characterization [25,26]. The widespread use of remote sensing technology in the analysis of rock discontinuities has demonstrated its feasibility and effectiveness.
Currently, the main algorithms used in this field fall into three categories: model fitting, region-growing, and unsupervised clustering. The Random Sample Consensus (RANSAC) algorithm and Randomized Hough Transform (RHT) are two typical model-fitting methods that can detect planes in 3D point clouds. Ge et al. [27] employ the region-growing method to determine the geometric characteristics, including orientation, spacing, and roughness, of rock discontinuities. They perform growth through the differences in normal vectors, thereby obtaining sets of discontinuities with different orientations. The performance of this method is highly dependent on the selection of initial seed points. Unsupervised clustering classifies point clouds based on the similarity of their underlying structure, geometric features, or spatial distribution. Unsupervised clustering methods typically surpass the model-fitting approach in extracting discontinuities from point cloud data, especially when dealing with complex geometries and a wide range of shapes and sizes [22]. K-means [28,29] and its variants [30,31] are common distance-based clustering methods. These methods directly assign groups to each point, which can lead to inaccurate division of boundary points. Slob et al. [32] proposed the automatic identification of discontinuities using fuzzy k-means clustering. This method introduced fuzzy theory and addressed the boundary grouping problem more effectively by calculating membership degrees. However, a significant drawback of such algorithms is the need to preset the number of clusters, making the clustering results highly influenced by subjective factors. To address this issue, scholars have proposed several algorithms that can automatically determine the number of clusters, including the surface classification algorithm [33], improved K-means [34], firefly algorithm [35], unsupervised K-means [36], and optimized K-means [37]. However, partition-based clustering methods, like region-growing techniques, are affected by the choice of initial cluster centers, which affects the clustering results. Density-based clustering algorithms can avoid this issue. The Density-Based Spatial Clustering of Applications with Noise (DBSCAN) determines cluster boundaries based on the density around data points and is commonly used to further extract single discontinuity clusters from main discontinuity sets [38]. Mean Shift (MS) is a window-based method that enables samples to converge at density peaks. Compared to DBSCAN, this algorithm handles high-dimensional feature data more effectively and requires only one parameter, the bandwidth. It is robust against outliers and is frequently used for identifying main discontinuity sets.
Daghigh et al. [39] conducted a qualitative evaluation of RHT, region-growing, and RANSAC on a rock slope point cloud dataset. The results indicate that these methods exhibit poor performance when segmenting the discontinuity plane, often detecting spurious planes and over-segmenting. Due to the curvature and roughness of rock mass surfaces, point clouds often exhibit undulations [34]. Some discontinuities exhibit significant local undulations or even fragmentation, which can lead to over-segmentation. Therefore, recognizing discontinuities is a complex process that typically requires the use of multiple algorithms in combination. Leng et al. [40] combined Hough-Transform and Region-Growing (HT-RG) methods to segment the rock mass planes. This method can detect discontinuities of different sizes and has good segmentation accuracy. However, this approach may miss detecting small and low-density plane surfaces and is prone to over-segmentation. Hu et al. [41] combined RANSAC with Voxel-Based Region-Growing (RAN-RG) for plane extraction from a point cloud. This method effectively improved the issue of detecting spurious planes, but still experienced over-segmentation in areas with significant local surface undulations. Yu et al. [42] introduced the Supervoxel concept to improve RAN-RG (Supervoxel-RG), which to some extent improved over-segmentation. However, it sometimes results in holes within the segmented plane. To address the issue of over-segmentation, Daghigh et al. [22] refined the surface variation of the rock point cloud by modifying the search neighborhood size and performing downsampling. However, some discontinuities with undulations or small offsets still remain. Chen et al. [26] presented a novel approach that considers points with high curvature values as noise and proposed a modified MS algorithm that incorporates five denoising operations. While this method can improve over-segmentation to some extent, it inevitably overlooks some small, highly undulating discontinuities.
In summary, existing methods for characterizing rock discontinuities are not effective in improving the issue of over-segmentation, which impacts the accuracy of the results. To tackle this problem, we propose a method to improve over-segmentation and extract discontinuities along with their orientations from 3D point cloud data of rock surfaces.

2. Methods

In this section, we introduce a workflow that integrates unsupervised clustering with RANSAC model fitting, utilizing a watershed-like algorithm to effectively address the issue of over-segmentation in rock discontinuities. The methodology is structured into five main steps:
Step 1: Adaptive resampling of the point cloud data based on density ratio, transforming the dataset through adaptive scaling according to density ratio. This task is detailed in Section 2.2.
Step 2: Estimating normal vector using Principal Component Analysis (PCA) and orienting them in the positive direction. This task is outlined in Section 2.3.
Step 3: Determination of the main discontinuity sets, which includes: (1) Applying a watershed-like algorithm to identify the non-planar areas of the discontinuities. This part is described in Section 2.4.2. (2) Segmenting the dominant discontinuity sets using the Multi-threshold Mean Shift (MTMS) algorithm. This part is covered in Section 2.4.3.
Step 4: Recognizing individual discontinuity clusters within each set using the DBSCAN algorithm. This task is explained in Section 2.5.
Step 5: Fitting a plane to each cluster from Step 4 using the RANSAC algorithm and calculating the orientation of each discontinuity. This task is described in Section 2.6.
A summary of this procedure is illustrated in Figure 1.

2.1. Description of the Experimental Datasets

The experimental dataset used in this study originates from a quarzitic rock mass at a road cut slope in Ouray, CO, USA (Figure 2a). It was acquired from the Rock bench repository [43] using an Optech Ilris 3D Laser Scanner. It consists of 1,515,722 points with an average point spacing of approximately 2 cm. Figure 2b shows its corresponding trimmed point cloud, which contains 1,024,521 points. This dataset is a popular resource in the fields of geotechnical and geological engineering for estimating the orientations of structural surfaces. In previous studies, two sets of ground truth orientations have been documented. Numerous researchers have created both semi-automatic and fully automatic methodologies to determine the orientations of structural surfaces within this point cloud.

2.2. Density-Ratio-Based Resampling Method

Point cloud data obtained from laser scanning often suffer from issues like uneven data density distribution due to factors such as vegetation cover, occlusion, and noise. Density-based clustering algorithms, such as DBSCAN, are highly sensitive to variations in point cloud density and struggle to detect all clusters in datasets with significantly varying densities. Therefore, in our method, we draw inspiration from an image processing technique known as histogram equalization. This technique applies nonlinear stretching to an image, redistributing the grayscale values so that the grayscale levels within a certain range become approximately equal. The most common method for equalization is the Cumulative Distribution Function (CDF).
Based on this principle, we substitute the grayscale values with the local point density of the point cloud and use a density estimation algorithm as the density estimator to estimate the local density of each sampled point. We then resample the data according to the CDF estimated over a neighborhood ( η ). The effect of this process is to achieve a roughly uniform spatial distribution of the data [44]. This step enhances the accuracy of the DBSCAN algorithm in recognizing single discontinuity clusters, as described in Section 2.5. The key symbols and notations utilized in this study are provided in Table 1.
The specific process is as follows.
(1)
The original dataset is D . When it is divided into ψ equal intervals, the boundary positions are { x 1 , x 2 , , x ψ + 1 } ;
(2)
At each boundary position x i , calculate p d f η ^ ( x i ) = 1 2 n η j Ι ( x i x j η ) ;
(3)
Calculate the resampled dataset: x D , y S y = x i x p d f η ^ ( x i ) .
The theoretical basis for steps (2) and (3) is as follows.
Let the CDF at the sampling point x be c d f η ^ ( x ) = x p d f η ^ ( x ) d x . We define the cumulative distribution function. We resample the dataset D into S using the follow mapping function:
x D , y S   y = c d f η ^ ( x )
This is a probability integral transform, where x D , y S   p d f η ^ ( y ) = p d f η ^ ( c d f η ^ ( x ) ) = 1 / n , representing a uniform distribution function.
Specifically, let p d f ε ^ ( y ) be the average density of y = c d f η ^ ( x ) within a neighborhood of size ε < η , then its density in the complete space is expressed as:
p d f ε ^ ( y ) = 1 2 n ε j Ι ( y y j ε ) = 1 2 n ε j Ι ( c d f η ^ ( x ) c d f η ^ ( x j ) ε ) 1 2 n ε j Ι ( P X p d f η ^ ( X [ min ( x , x j ) , max ( x , x j ) ] ) ε )
When n = | D | is large and varies slowly around x , this expression can be approximated as:
p d f ε ^ ( y ) = 1 2 n ε j Ι ( p d f ε ^ ( x ) x x j ε ) 1 2 n ε j Ι ( x x j ε ) p d f η ^ ( x ) = p d f ε ^ ( x ) p d f η ^ ( x )
Therefore, p d f ε ^ ( y ) r p d f ^ ε , η ( x ) , i.e., when resampling with a smaller ε , the resampling result is scaled by different degrees according to the density ratio between the small neighborhood ε and the large neighborhood around the sampling point. In other words, when the sampling point x is in a low-density region, the points in that area are less thinned out, while in a high-density region, the points are thinned out more significantly. This process is summarized in Table 2.

2.3. Normal Vector Calculation

Accurate normal estimation in a point cloud is essential for effectively segmenting the desired plane features. Principal Component Analysis (PCA) is a statistical method that focuses on variation to uncover significant patterns in a dataset through an orthogonal transformation. PCA is widely used to identify local planar surfaces and linear points within a 3D point cloud. Using KD-tree nearest neighbor search, extract the neighborhood point sets ( D 1 , D 2 , D 3 D n ) for each sample point within the search radius ( R 1 ,   R 2 ,   R 3 R n ), respectively. Calculate the centroids of each neighborhood point set ( S 1 ,   S 2 ,   S 3 S n ). Construct the covariance matrices ( C 1 ,   C 2 ,   C 3 C n ) relative to the centroids, which reflect the correlation of each coordinate component in the point sets. The 3×3 covariance matrix for a sample point p i ( x , y , z ) can be written as:
C = 1 k i = 1 k ( p i p ¯ ) T ( p i p ¯ ) ; p ¯ = 1 k i = 1 k p i
Utilize Singular Value Decomposition (SVD) on the covariance matrix of each point. The covariance matrix C can be represented as:
C = [ e 1 e 2 e 3 ] [ λ 1 0 0 0 λ 2 0 0 0 λ 3 ] [ e 1 T e 2 T e 3 T ]
Among them, λ 1 , λ 2 , λ 3 are the eigenvalues of C , and e 1 , e 2 , e 3 are the corresponding eigenvectors. After normalizing the eigenvalues and sorting them by size, we have λ 1 λ 2 λ 3 > 0 , where λ 1 , λ 2 and λ 3 represent the fitting deviations in the directions of the corresponding eigenvectors of the local plane. The eigenvector corresponding to the smallest eigenvalue, λ 3 , indicates the direction of the least variation in the point set distribution. This direction can be considered the normal direction of the local surface.
Since PCA can only provide a direction without distinguishing between “inward” and “outward”, an additional post-processing step is required to correct any normal vectors that are flipped by 180°. This study uses a method based on the continuity of normal directions to adjust and ensure the consistency of the normal vector orientations.

2.4. Determination of the Main Discontinuity Sets

Identifying group discontinuities is fundamental for pinpointing individual discontinuities, which is essential for accurately identifying rock discontinuities. This phase focuses on categorizing each point into a discontinuity group based on distribution patterns.
MS is a feature space analysis method that estimates the probability density distribution of feature data through Kernel Density Estimation (KDE). It clusters data points by moving them in the direction of increasing probability density until they converge at the local density maximum. The method can process point clouds of any shape without needing a predefined number of clusters. It requires only one parameter, the bandwidth size. The bandwidth determines the step size for each adjustment of the probability density function. For MS, tuning the bandwidth is the most critical hyperparameter.

2.4.1. The Impact of Bandwidth on MS

The bandwidth, as the sole input parameter for MS, has a crucial impact on the segmentation effect. Figure 3a illustrates the kernel density of a one-dimensional dataset. The dashed lines indicate the kernel function for each sample, while the solid line depicts the probability density of the dataset, which is the aggregate of the sample kernel functions. It can be observed that the probability density of this dataset has only one local maximum, so the current dataset is segmented into one cluster using MS with a bandwidth of 6. If we set the bandwidth to approach zero infinitely, the sample kernel functions will be as shown in Figure 3b, approximating impulse functions, and the probability density will be as shown in Figure 3c. The dataset will have five local maxima and be segmented into five clusters. Therefore, the bandwidth determines the probability density of the dataset. If the bandwidth is too large, MS may overly smooth the data, failing to capture local details and features, resulting in under-segmentation. Conversely, if the bandwidth is too small, it will include fewer points, leading to the amplification of some local details and features, causing over-segmentation.
Figure 4 illustrates the effect of bandwidth on different types of discontinuities. Figure 4a shows a representative area of the experimental dataset with flat discontinuities. The spatial orientation of plane 1, plane 2, and plane 3 is illustrated in Figure 4b, where each plane belongs to a different set of discontinuities. By comparing segmentation effects with different bandwidths in this area, we find that the segmentation effect in Figure 4c with a bandwidth of 0.25 is ideal because it effectively segments the area into three discontinuities while preserving edge details. For the non-planar discontinuities targeted in this study (Figure 4d), Figure 4e simulates their planes in 3D space. It is clear that the non-planar discontinuities fit well with the plane model, suggesting that most points on the non-planar discontinuities should theoretically be classified into the same discontinuities set. Figure 4f shows the segmentation effect of using different bandwidths for MS in this area. When the bandwidth is 0.25, the points are over-segmented. To avoid over-segmentation, the bandwidth is increased, and the result obtained with a bandwidth of 0.35 is considered the most suitable representation because most points in this area are assigned to the same discontinuities set. However, for flat discontinuities, when the bandwidth is 0.35 (see Figure 4c), plane1 and plane3 are assigned to the same discontinuities set, resulting in under-segmentation.
Based on the comprehensive analysis above, it is evident that the optimal bandwidth for non-planar discontinuities and flat discontinuities differs. However, most existing MS methods use a global fixed bandwidth [22,26]. To ensure accuracy, this often leads to over-segmentation in areas with non-planar discontinuities. On the other hand, expanding the threshold range can reduce accuracy when segmenting flat discontinuities. This makes it difficult to balance both types of discontinuities within the same dataset, leading to over-segmentation or under-segmentation. Therefore, it is necessary to identify non-planar areas and set different bandwidths for points located in these areas and those in flat discontinuities. To tackle this issue, we propose a new method for detecting individual discontinuities, which primarily includes: (1) applying a watershed-like algorithm to identify non-planar areas of discontinuity. (2) using the Multi-threshold Mean Shift (MTMS) algorithm to segment the advantageous discontinuity sets. Detailed descriptions of these improvements are provided in the following subsections.

2.4.2. Watershed-like Algorithm for Identifying Non-Planar Areas

The Watershed algorithm is an image segmentation technique that separates and labels different regions within an image. It is based on the grayscale values or gradient information of the image, treating the depth image as a topographic map [45]. The algorithm assigns a gradient value to each pixel by calculating gradient operators, and simulating peaks and valleys (basins). As illustrated in Figure 5, assuming there are small holes on the surface of the minimal regions of the terrain, water seeps in through these holes and gradually floods the surrounding areas at a uniform rate, forming corresponding catchment basins. The boundary line of each catchment basin is known as the watershed line [46,47,48]. Figure 6 demonstrates the Watershed algorithm’s segmentation process in an image [49]. Figure 6a shows a limestone-blasted grayscale image, treated as a gradient image, where the grayscale values at corresponding coordinates represent elevation. Figure 6b depicts the 3D simulation of the water catchment basin in this image, assuming a gray threshold of 50 as the watershed line. Figure 6c shows the 3D simulation post-flooding, and Figure 6d presents the segmented 2D cross-section.
Current researchers are exploring the application of watershed-like concepts or modified algorithms to point cloud data, as the density or other attributes of point clouds may, to some extent, resemble “height” or “intensity”. By simulating this analogy, “terrain” in 3D space can be defined. Based on this idea, this paper proposes a watershed-like algorithm that uses non-planar indicators as gradient operators to restrict the flooding process of the watershed algorithm, thereby better identifying non-planar regions of discontinuities.

Calculate the Non-Planar Indicator

In Section 2.3, we used PCA to compute the local covariance matrix eigenvalues λ 1 , λ 2 , λ 3 as shown in Equation (2). According to the theory of Demantké et al., as shown in Figure 7, the eigenvalues can represent the dimensional properties of the local neighborhood within the point cloud [50]. When λ 1 λ 2 λ 3 , the local neighborhood of the point cloud tends to be rod-shaped; when λ 1 λ 2 λ 3 , the local neighborhood of the point cloud tends to be planar; when λ 1 λ 2 λ 3 , the local neighborhood of the point cloud tends to be a three-dimensional surface.
Therefore, we construct a local neighborhood dimension feature model based on these dimensional characteristics [51]:
L λ = λ 1 λ 2 λ 1 ,   P λ = λ 2 λ 3 λ 1 ,   S λ = λ 3 λ 1
In the study by Demantké et al. [50], the dimensional features L λ , P λ and S λ can, respectively, represent the probability that the local neighborhood of the point is a one-dimensional line, a two-dimensional plane, or a three-dimensional surface [50]. The planar dimension feature P λ indicates how well the local neighborhood of the point cloud aligns with a plane. The closer the point is to a planar structure, the larger its corresponding dimension feature P λ . Therefore, this paper defines an indicator to describe the curvature of the local neighborhood:
ω = ( 1 P λ ) 2
From Equation (7), it can be observed that the smaller the value of ω , the closer the local neighborhood shape is to a plane. When the local neighborhood is a curved surface, the value corresponding to ω is a fixed value related to the curvature of the surface.
To distinguish between non-planar and flat areas, as shown in Figure 8, surfaces with different degrees of curvature have varying normal direction deviations within the same neighborhood range. The more curved the surface, the greater the deviation in the normal direction. Therefore, we use the angle between the normal vectors of the sampling point and its neighboring points as a reference factor.
To simplify the complex curvature calculation process, we establish a local neighborhood V r , centered on the sampling point. The size of V r is customized according to the scale of the rock mass. Given that the number of neighborhood points within V r is k , the curvature at the sampling point is defined as follows:
C = 1 k i = 1 k λ i
where λ i represents the deviation between the sampling point and the normal vectors of the neighboring points within the local neighborhood. It can be observed that as λ i increases, indicating a greater degree of bending in the local structure, the curvature value correspondingly increases.
Combining the two-dimensional feature ω , this paper defines a non-planar indicator:
K = C ω
where C is an indicator that measures the degree of bending within the local neighborhood, and ω describes the curvature of the local neighborhood. The shape characteristics of the point cloud’s local neighborhood are determined based on the non-planar indicator K . When the local neighborhood is close to being planar, the non-planar indicator approaches 0. As K increases, it indicates a greater degree of bending in the local neighborhood of the point cloud.

Fitting Spatial Gradient Fields

Due to the similarity in curvature among closely spaced points, we sample the point cloud data to reduce computational load and improve efficiency. A voxelized grid method is used to create a voxel grid, which consists of 3D boxes in space, to encompass the input point cloud. The voxel size is determined by the user based on specific requirements. To sample feature points, select point p c e n t r o i d  located at the centroid of each voxel, and calculate the non-planar indicator K for each feature point in the experimental dataset. Figure 9 is an HSV simulation diagram of the experimental dataset colored according to the non-planar indicator K . It can be observed that the values of K in this dataset primarily range 0 0.2 , The smaller the K value, the flatter the region where the point is located; the larger the K value, the more curved the point’s location is. We use the ascending K values to simulate the elevation from valleys (basins) to peaks in the high-dimensional features, with the peaks representing the non-planar areas we aim to identify.
The watershed algorithm typically processes gradient data [52]. In this study, the Sobel operator is used to calculate the gradient data of the voxels V . Assuming V L T , V R T , V L U , V R U , V F T , V B T , V F U , V B U , V T , V U , V L , V R , V F , V B represent the points in the 14-neighborhood of voxel V in the directions: top-left, top-right, bottom-left, bottom-right, front-top, back-top, front-bottom, back-bottom, top, bottom, left, right, front, and back, respectively. The gradient values of voxel V in the x, y and z directions, denoted as d x ( V ) , d y ( V ) and d z ( V ) , are given by:
{ d x ( V ) = [ K ( p R T ) + 2 K ( p R ) + K ( p R U ) ] [ K ( p L T ) + 2 K ( p L ) + K ( p L U ) ] d y ( V ) = [ K ( p L T ) + 2 K ( p T ) + K ( p R T ) ] [ K ( p L U ) + 2 K ( p U ) + K ( p R U ) ] d z ( V ) = [ K ( p F T ) + 2 K ( p F ) + K ( p F U ) ] [ K ( p B T ) + 2 K ( p B ) + K ( p B U ) ]
The gradient magnitude d ( V ) is:
d ( V ) = [ d x ( V ) ] 2 + [ d y ( V ) ] 2 + [ d z ( V ) ] 2
For each point p i , in the point cloud data, assuming the index of its voxel is V p i , we approximate the gradient value d ( p i ) of the point p i as being equal to the gradient value of the voxel it resides in, that is
d ( p i ) = d ( V p )
Since point cloud data exist in 3D space, it is not feasible to simulate ridges and valleys on the same horizontal plane as one would with images. To facilitate the simulation of the watershed inundation process, we use the normal direction of each point as the vertical direction of the peaks. Through this approach, we define the “terrain” in the multidimensional space, allowing for more convenient simulation and analysis.
To determine the frequency distribution of each gradient level, we need to utilize a sorting algorithm to arrange the experimental dataset by gradient magnitude. Each point is then assigned to an ordered queue corresponding to its gradient level. Assume the number of gradient levels is N :
N = INT ( d max d min d g a p )
where d max is the maximum gradient value, d min is the minimum gradient value, d g a p is the interval of the gradient levels, and INT denotes the rounding operation.
The gradient level K p at point p i is shown as
K p = INT [ d ( p ) d min d g a p ] + 1
From Equation (13), it can be seen that the gradient level interval d g a p determines the number of gradient levels N, which in turn determines the number of points submerged during each step of the flooding process. According to Equation (14), the gradient level for each point can be calculated. The flooding process submerges points sequentially according to their gradient levels.

Flooding

First, let us clarify some basic concepts in this method:
Minimum Region: The minimum region M is an independent connected area where the gradient values within the region are lower than the gradient values of the surrounding neighboring points.
Water Catchment Basin: The water catchment basin B ( M ) is an area composed of a collection of points within the experimental dataset that form the minimum region M . Within the water catchment basin, water flow from any location will slide down along a certain path to reach the minimum region M .
Geodesic Distance: Assuming S is a connected region, the geodesic distance D S ( j , k ) is along the geodesic path from j to k within S .
Geodesic Influence Region: Let set S contain set W , where W consists of multiple connected regions W 1 , W 2 , , W k . Within set S , if there exists a point set Q whose geodesic distance to connected region W i is shorter than its distance to any other connected region within set W , then Q is considered the geodesic influence region Z S ( W i ) of connected region W i .
The flooding process is essentially a recursive process. Suppose the experimental dataset is E , and T f ( p ) is the threshold set for point p at level f ,
T f ( p ) = { p E , K p f }
The flooding process starts from the threshold set T 1 ( p ) , which is the set of points with the smallest gradients. Let the collection of first-level basins be:
U 1 = T 1 ( p )
Let n be the gradient level, and the recursive process of the f -th level basin U f is
f [ 1 , n 1 ] , U f + 1 Z T f + 1 ( p ) ( U f )
The set of catchment basins with w as the watershed line is
U w = 1 w 1 { U i + 1 Z T i + 1 ( p ) ( U i ) }
For the experimental dataset E , the set of non-planar areas H w (slope peaks) is
H w = E U w
According to Equation (18), the watershed line w determines the range of slope peaks and valleys. For the experimental dataset, we set the number of gradient levels to N = 15 , Figure 10 shows a comparison of the flooding results for different watershed lines, where the red area represents the non-planar area (slope peaks). It can be observed in Figure 10a–f that as the value of w decreases, the range of the identified non-planar area gradually increases. However, in Figure 10g,h, even though w continues to decrease, the range of the identified non-planar area remains almost unchanged. The study found that there exists a critical value for the watershed line, which distinguishes between planar areas and non-planar areas—essentially, the cliff edge of the peak. Due to the significant gradient difference between planar areas and non-planar areas, the identification results remain nearly constant when w continues to decrease. In the experimental dataset, w = 8 represents the critical value of the watershed line, where the division at the peak of the slope accurately reflects the location of non-planar area. Therefore, we use w = 8 as the watershed line for the experimental dataset.

2.4.3. The Multi-Threshold Mean Shift (MTMS) Algorithm

Currently, most discontinuity recognition methods based on MS employ a fixed bandwidth. However, using a single global bandwidth for rock slopes with numerous non-planar areas is not optimal, as a single global value may result in over-segmentation or under-segmentation. Section 2.4.1 provides a detailed explanation of the impact of bandwidth on segmentation performance. MTMS is an improvement over the traditional MS method. This approach assigns different bandwidths to points in planar and non-planar areas during the shifting process. By doing so, it controls the step size of the probability density function adjustment for points in different regions. This method achieves a better balance between local and global features. The specific process is shown in Figure 11.
Map the estimated normal vectors from Section 2.3 to the feature space, extending from a 2D Gaussian surface to a 3D Gaussian sphere (Figure 11a), as the input for MTMS. Define a flag value to mark whether a point is in a non-planar area. Search for the points located in the non-planar area defined in Section 2.4.2 in spherical coordinates, setting their flag to 1, while setting the flag to 0 for the remaining points. The resulting binary map is shown in Figure 11b.
Based on the flag, construct a trigger-based bandwidth allocation mechanism as detailed in Figure 11c. Randomly select a point from all unprocessed points as the center point, check its flag label, and if flag = 1, set its bandwidth to h f l a g = 1 ; if flag = 0, set its bandwidth to h f l a g = 0 . The kernel function of MS assigns certain weights to points around each point. In the feature space, it weights the points within a specified bandwidth region around the current center point to estimate the local density. The density calculation formula is
f ( x ) = 1 n h d i = 1 n K ( x x i h ) ,
where K ( x ) is the kernel function weight formula, n is the number of features, h is the bandwidth of the current search point, d is a normalization constant that ensures K ( x ) = 1.
Calculate the gradient direction of the dataset density, which is the direction of the fastest increase in local density, to obtain the drift vector. The gradient of the density function is:
f ( x ) = 2 C k , d n h d + 2 i = 1 n ( x i x ) [ i = 1 n x i g ( x x i h 2 ) i = 1 n g ( x x i h 2 ) x ] ,
where 2 C k , d n h d + 2 i = 1 n ( x i x ) is a real number, so the drift vector is
m h ( x ) = i = 1 n x i g ( x x i h 2 ) i = 1 n g ( x x i h 2 ) x
Translate the sample point x iteratively according to the drift vector. The center point for the next iteration is
M ( x ) = x + m h ( x i )
Repeat the iteration several times until m h ( x i ) = 0 , so that the samples converge to the local density peak, resulting in multiple sets of dominant points with similar principal directions.

2.5. Recognition of Single Discontinuity Cluster

To distinguish each individual planar surface among many planes with similar orientations in the primary cluster, secondary clustering is required. Within each discontinuity set, each single discontinuity creates a dense region in space, resulting in a density difference between the discontinuities. This characteristic makes them ideal for separation using density-based methods [26]. DBSCAN is a typical density-based clustering algorithm. It has performed well in previous studies for identifying individual discontinuities [21,22,26,53]. In this study, the DBSCAN algorithm is utilized, with parameters configured for the neighborhood search radius ( e p s ) and the minimum number of points needed to form a cluster ( M i n P t s ). The algorithm starts with a randomly chosen starting point p and searches within a radius of e p s . If there are at least M i n P t s points within this radius, the point ( p ) is classified as a core point. If the number of points is fewer than M i n P t s , the point is considered a boundary point. Core points, along with all points within the radius e p s , are grouped into a cluster. This process is repeated for each unprocessed point, ultimately yielding the second-stage clusters that belong to each structural plane. the results of identifying the first set of discontinuities are illustrated in Figure 12.

2.6. Discontinuity Plane Fitting

Given that rock discontinuity surfaces are typically rough and uneven, the orientation of a discontinuity is determined by the orientation of its best-fitting plane. RANSAC is an iterative method designed to estimate surface model parameters from a dataset. It identifies the parameters of a geometric shape by utilizing the smallest number of points needed to define that shape within a point cloud [54,55,56].
It accounts for the variability of the points, offering a more accurate estimation. The RANSAC algorithm is widely applied in shape detection [57,58] and has been effectively used in previous research for plane fitting [34,59].
The RANSAC algorithm repeatedly applies a plane equation to fit the discontinuity, aiming to find the best parameters. The equation for the discontinuity plane is expressed as:
A x + B y + C z + D = 0
In this equation, A, B, C, and D are the coefficients, with A, B, and C representing the three components of the plane’s normal vector. D denotes the perpendicular distance from the origin to the plane. The dip angle and dip direction of the planes are determined using Equations (25) and (26).
Dip   direction = { 90 t a n 1 ( B A ) A > 0 0 A = 0 & B 0 180 A = 0 & B < 0 270 t a n 1 ( B A ) A < 0
Dip   angle = { 0 A 2 + B 2 = 0 90 t a n 1 ( | C | A 2 + B 2 ) A 2 + B 2 0

2.7. Qualitative Segmentation Results

Daghigh et al. [39] reviewed the experimental findings of reference [60] and reference [25], highlighting various examples of deficiencies associated with rectangles and identifying the type of deficiency (OS: over-segmentation; US: under-segmentation) (Figure 13a,b). To optimize segmentation results, Daghigh et al. [22] modified the KNN value and utilized VG downsampling to reduce the surface variation in the point cloud. Nevertheless, certain discontinuities with undulations or minor offsets remain (as shown in Figure 13c, highlighted by the circled and rectangular areas). The visual segmentation results from these studies, along with the current study, are displayed in Figure 13. It can be seen that the results of the present study (Figure 13d) show an ideal segmentation effect in the areas marked with deficiencies, demonstrating that our method for determining the main discontinuity sets can effectively avoid over-segmentation.
For the extraction results of discontinuity orientations, the experimental dataset used in this study has several sets of true orientation values that have been published by other researchers. Riquelme et al. [53] and Chen et al. [60] manually measured the orientations of several typical discontinuities in this dataset. Daghigh et al. [39] conducted a critical evaluation of these measurements and identified significant errors in the data obtained by Riquelme et al. [53] and Chen et al. [60]. As a result, Daghigh et al. [39] re-evaluated the discontinuity orientations by employing the Segment tool in CloudCompare. This method entails drawing a closed polygon around the points related to a particular discontinuity plane and then using the best-fit plane tool to ascertain the plane’s orientation based on the chosen points. This paper uses these results as the true values for discontinuity orientations. Figure 14 shows the numbering of the discontinuities used for comparison, and Table 3 compares the measurement results of our method with those of reference [22] and the true values. It can be observed that our method has smaller errors compared to reference [22], validating that our method provides satisfactory characterization results.

3. Case Studies and Discussion

3.1. Case Dataset

To illustrate the proposed methodology, two case studies of rock slopes were analyzed. In Case A, data were gathered in Kingston, Ontario, Canada, using a Leica HDS600 scanner (Figure 15a). Figure 15b presents the corresponding trimmed point cloud (C1), which comprises 387,610 points with an average spacing of about 4 cm. Case B: A rock slope in a quarry in Wuhan, China (Figure 15c). Figure 15d shows the corresponding trimmed point cloud (C2), consisting of 266,294 points with an average point spacing of approximately 3 cm.

3.2. Case A

For C1, to effectively evaluate the segmentation results, Hu et al. [41] generated a colored ground truth point cloud. This ground truth was later revised by Daghigh et al. [22] and utilized as a benchmark for evaluating the segmentation outcomes (Figure 16). Therefore, we use this colored ground truth point cloud as the standard for assessing the segmentation results obtained by our method.
Figure 17 illustrates the results of employing the watershed-like algorithm to identify non-planar areas in C1. It is evident that the non-planar regions detected by the proposed method encompass those identified in Figure 17, while also revealing additional details. This indicates that the method is effective in accurately identifying non-planar areas.
Figure 18 is the comparison of discontinuous surface segmentation results. Figure 18a shows the color-coded ground truth point cloud. Figure 18b–g illustrate the segmentation results for several typical methods from Liu et al. [61], with red polygons used to emphasize the shortcomings in each outcome. The ellipse marks areas of over-segmentation, the rectangle indicates discontinuities, the parallelogram highlights incorrect results, the rhombus shows poor boundary definitions, and the trapezoid identifies missed detections [61]. Figure 18h,i is the segmentation results of MOE [61] and Supervoxel-RG [42] on C1, as summarized in reference [39]. It can be observed that MOE shows gaps in non-planar areas (marked with black rectangles), while Supervoxel-RG exhibits over-segmentation in non-planar areas (indicated by red ellipses). Reference [22] proposed an improved method and applied it to C1. As depicted in Figure 18j, this approach somewhat mitigates over-segmentation, proving most effective for large-area planes (such as blue and some red planes). For non-planar areas with high curvature, this method can still result in incorrect segmentation (marked with white rectangles). Figure 18k shows the segmentation results of this method on C1. It can be observed that the results align well with the color ground truth point cloud shown in Figure 18a. Additionally, in areas marked by red ellipses, where over-segmentation was prone to occur, the segmentation performance is improved, demonstrating the method’s high segmentation accuracy and its effectiveness in mitigating over-segmentation.

3.3. Case B

Figure 19 shows the recognition results of the watershed-like algorithm. We computed segmentation results for different bandwidths (Figure 20). Analysis reveals that for flat discontinuities (the blue areas in the figure), the best segmentation result is achieved when bandwidth = 0.24. For non-planar areas (the red areas in the figure), a bandwidth of 0.3 results in clearly defined and complete boundaries of the identified point sets with almost no over-segmentation. Therefore, we use 0.24 and 0.3 as the bandwidths for flat discontinuities and non-planar discontinuities, respectively.
Figure 21a shows the final discontinuity surface recognition results. By analyzing the 2D profiles of two typical non-planar discontinuities (as shown in Figure 21c), it can be observed that their spatial distribution approaches a plane, and the segmentation results meet expectations. Figure 21b provides a spherical projection of the characterization results, which includes four discontinuity groups and a total of 51 discontinuities.
To validate the accuracy of the results, we measured the directions using CloudCompare. Figure 22 shows the eight sample planes used for verification (two from each discontinuity concentration), and Table 4 compares our method with CloudCompare measurements. The average errors in dip direction and dip angle for the eight sample planes are 1.2 and 1.1, respectively, indicating that our method has good applicability.

3.4. Computation Time Comparison

Numerous scholars have conducted experiments on the dataset used in this study and recorded the runtime. Table 5 compares the computation time of the proposed method with the runtime of methods proposed by other researchers. The results indicate that the proposed method has a faster runtime, which meets the needs of practical engineering applications. It is important to note that the runtime of a method can be influenced by various factors, such as computer hardware, the programming tools employed, and the complexity of the method. Therefore, the runtimes summarized in Table 5 should not be considered the sole criterion for evaluating the efficiency of a method.
As observed in Table 5, our method takes two minutes longer than the method proposed by Daghigh et al. [3]. Upon analysis, it was found that the most time-consuming part of our method is the process of identifying the main discontinuity sets. The identification of non-planar areas has somewhat reduced efficiency; however, considering its contribution to optimizing the segmentation results, the additional time consumed by this step is acceptable.

4. Conclusions

This paper primarily investigates the problem of discontinuity characterization in unstructured point clouds of rock masses. To address the issue of over-segmentation that often arises when characterizing complex rock mass discontinuous surfaces, a new method based on laser point cloud data is proposed. The method is divided into five main stages: (1) adaptive resampling; (2) normal vector calculation; (3) determination of the main discontinuity sets; (4) recognition of single discontinuity cluster; and (5) discontinuity plane fitting.
The primary novel contributions of the proposed work are: (1) adaptive resampling of point cloud data based on density ratios, which has the capacity to process datasets with varying densities; (2) the watershed-like algorithm can identify the non-planar areas on the rock mass surface; (3) proposed MTMS, which includes a trigger-based bandwidth allocation mechanism, effectively avoiding over-segmentation in non-planar areas (4) accurate calculation of the orientation of discontinuities; (5) fewer critical hyperparameters require adjustment compared to other methods.
To verify the reliability of the results, we applied this method to three real rock mass datasets and compared it with previous research results and manual measurements. For the experimental dataset, it can be seen that the results of the present study achieve satisfactory segmentation in the set of marked deficiency areas, demonstrating that our method for determining the main discontinuity sets effectively avoids over-segmentation. When comparing the extracted discontinuity orientations from this dataset with the ground truth published by other researchers, our method shows smaller errors. For Case A, the identification results from our method align closely with the color-coded ground truth point cloud. Moreover, in areas previously recognized by scholars as prone to over-segmentation, our method provides better segmentation results, proving its high accuracy and effectiveness in mitigating over-segmentation. Furthermore, we validated our method on Case B by comparing the computed discontinuity orientations with measurements from CloudCompare. The average errors for the dip direction and dip angle across eight sample surfaces were all within 2°, indicating good generalizability of our method. We also analyzed the computational efficiency of our method by comparing its computation time with several representative methods on the same dataset. The results demonstrate that our method has a satisfactory computational efficiency, meeting the practical needs of engineering applications, especially for complex rock masses with many non-planar areas.

Author Contributions

Conceptualization, Y.L.; visualization, W.H.; validation, Y.L.; writing—original draft preparation, Y.L.; writing—review and editing, W.H., Q.C. and X.L.; supervision, X.L.; project administration, X.L.; funding acquisition, X.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Key R&D Program of Hubei (Grant No. 2022BCA080).

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

The authors would like to thank the Rock bench repository for providing the 3D point cloud data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hudson, J.; Harrison, J.; Popescu, M. Engineering Rock Mechanics: An Introduction to the Principles. Appl. Mech. Rev. 2002, 55, B30. [Google Scholar] [CrossRef]
  2. Zhang, F.; Damjanac, B.; Maxwell, S. Investigating Hydraulic Fracturing Complexity in Naturally Fractured Rock Masses Using Fully Coupled Multiscale Numerical Modeling. Rock. Mech. Rock. Eng. 2019, 52, 5137–5160. [Google Scholar] [CrossRef]
  3. Barton, N.R. Suggested methods for the quantitative description of discontinuities in rock masses. ISRM Congr. 1978, 15, 319–368. [Google Scholar]
  4. Sturzenegger, M.; Stead, D. Close-range terrestrial digital photogrammetry and terrestrial laser scanning for discontinuity characterization on rock cuts. Eng. Geol. 2009, 106, 163–182. [Google Scholar] [CrossRef]
  5. Goodman, R.E.; Shi, G.-H. Block Theory and Its Application to Rock Engineering; Springer: New York, NY, USA, 1985. [Google Scholar]
  6. Pantelidis, L. Rock slope stability assessment through rock mass classification systems. Int. J. Rock. Mech. Min. Sci. 2009, 46, 315–325. [Google Scholar] [CrossRef]
  7. Alejano, L.R.; Ferrero, A.M.; Ramírez-Oyanguren, P.; Álvarez Fernández, M.I. Comparison of limit-equilibrium, numerical and physical models of wall slope stability. Int. J. Rock. Mech. Min. Sci. 2011, 48, 16–26. [Google Scholar] [CrossRef]
  8. Mah, J.; Samson, C.; McKinnon, S.D. 3D laser imaging for joint orientation analysis. Int. J. Rock. Mech. Min. Sci. 2011, 48, 932–941. [Google Scholar] [CrossRef]
  9. Abellán, A.; Oppikofer, T.; Jaboyedoff, M.; Rosser, N.J.; Lim, M.; Lato, M.J. Terrestrial laser scanning of rock slope instabilities. Earth Surf. Process. Landf. 2014, 39, 80–97. [Google Scholar] [CrossRef]
  10. Andrews, B.J.; Roberts, J.J.; Shipton, Z.K.; Bigi, S.; Tartarello, M.C.; Johnson, G. How do we see fractures? Quantifying subjective bias in fracture data collection. Solid. Earth 2019, 10, 487–516. [Google Scholar] [CrossRef]
  11. Pagano, M.; Palma, B.; Ruocco, A.; Parise, M. Discontinuity Characterization of Rock Masses through Terrestrial Laser Scanner and Unmanned Aerial Vehicle Techniques Aimed at Slope Stability Assessment. Appl. Sci. 2020, 10, 2960. [Google Scholar] [CrossRef]
  12. Riquelme, A.; Tomás, R.; Cano, M.; Pastor, J.L.; Abellán, A. Automatic Mapping of Discontinuity Persistence on Rock Masses Using 3D Point Clouds. Rock Mech. Rock. Eng. 2018, 51, 3005–3028. [Google Scholar] [CrossRef]
  13. Chen, J.; Huang, H.; Zhou, M.; Chaiyasarn, K. Towards semi-automatic discontinuity characterization in rock tunnel faces using 3D point clouds. Eng. Geol. 2021, 291, 106232. [Google Scholar] [CrossRef]
  14. Zhang, K.; Wu, W.; Zhu, H.; Zhang, L.; Li, X.; Zhang, H. A modified method of discontinuity trace mapping using three-dimensional point clouds of rock mass surfaces. J. Rock. Mech. Geotech. Eng. 2020, 12, 571–586. [Google Scholar] [CrossRef]
  15. Cai, X.h.; Lü, Q.; Zheng, J.; Liao, K.W.; Liu, J. An efficient adaptive approach to automatically identify rock discontinuity parameters using 3D point cloud model from outcrops. Geol. J. 2023, 58, 2195–2210. [Google Scholar] [CrossRef]
  16. Janowski, L.; Wroblewski, R.; Rucinska, M.; Kubowicz-Grajewska, A.; Tysiac, P. Automatic classification and mapping of the seabed using airborne LiDAR bathymetry. Eng. Geol. 2022, 301, 106615. [Google Scholar] [CrossRef]
  17. Liu, Y.; Chen, J.; Tan, C.; Zhan, J.; Song, S.; Xu, W.; Yan, J.; Zhang, Y.; Zhao, M.; Wang, Q. Intelligent scanning for optimal rock discontinuity sets considering multiple parameters based on manifold learning combined with UAV photogrammetry. Eng. Geol. 2022, 309, 106851. [Google Scholar] [CrossRef]
  18. Zekkos, D.; Greenwood, W.; Lynch, J.; Manousakis, J.; Athanasopoulos-Zekkos, A.; Clark, M.; Saroglou, C. Lessons learned from the application of UAV-enabled structure-from-motion photogrammetry in geotechnical engineering. Int. J. Geoengin. Case Hist. 2018, 4, 254–274. [Google Scholar]
  19. Rodriguez, J.; Macciotta, R.; Hendry, M.; Roustaei, M.; Gräpel, C.; Skirrow, R. UAVs for monitoring, investigation, and mitigation design of a rock slope with multiple failure mechanisms—A case study. Landslides 2020, 17, 2027–2040. [Google Scholar] [CrossRef]
  20. Xu, Z.-h.; Guo, G.; Sun, Q.-c.; Wang, Q.; Zhang, G.-d.; Ye, R.-q. Structural plane recognition from three-dimensional laser scanning points using an improved region-growing algorithm based on the robust randomized Hough transform. J. Mt. Sci. 2023, 20, 3376–3391. [Google Scholar] [CrossRef]
  21. Ge, Y.; Cao, B.; Tang, H. Rock Discontinuities Identification from 3D Point Clouds Using Artificial Neural Network. Rock Mech. Rock. Eng. 2022, 55, 1705–1720. [Google Scholar] [CrossRef]
  22. Daghigh, H.; Tannant, D.D.; Jaberipour, M. A computationally efficient approach to automatically extract rock mass discontinuities from 3D point cloud data. Int. J. Rock. Mech. Min. Sci. 2023, 172, 105603. [Google Scholar] [CrossRef]
  23. Liu, H.; Jisen, S.; Tovele, G.S.V.; Tao, C.; Shuzhao, C.; Mawugnon, B.K.; Peng, L. Application of photogrammetry and in-situ test technology in the stability evaluation of gangue dump slope. Bull. Eng. Geol. Environ. 2022, 82, 2. [Google Scholar] [CrossRef]
  24. Cirillo, D.; Zappa, M.; Tangari, A.C.; Brozzetti, F.; Ietto, F. Rockfall Analysis from UAV-Based Photogrammetry and 3D Models of a Cliff Area. Drones 2024, 8, 31. [Google Scholar] [CrossRef]
  25. Kong, D.; Wu, F.; Saroglou, C. Automatic identification and characterization of discontinuities in rock masses from 3D point clouds. Eng. Geol. 2020, 265, 105442. [Google Scholar] [CrossRef]
  26. Chen, Q.; Ge, Y.; Tang, H. An unsupervised method for rock discontinuities rapid characterization from 3D point clouds under noise. Gondwana Res. 2024, 132, 287–308. [Google Scholar] [CrossRef]
  27. Ge, Y.; Tang, H.; Xia, D.; Wang, L.; Zhao, B.; Teaway, J.W.; Chen, H.; Zhou, T. Automated measurements of discontinuity geometric properties from a 3D-point cloud based on a modified region growing algorithm. Eng. Geol. 2018, 242, 44–54. [Google Scholar] [CrossRef]
  28. Jimenez, R. Fuzzy spectral clustering for identification of rock discontinuity sets. Rock. Mech. Rock. Eng. 2008, 41, 929–939. [Google Scholar] [CrossRef]
  29. Assali, P.; Grussenmeyer, P.; Villemin, T.; Pollet, N.; Viguier, F. Surveying and modeling of rock discontinuities by terrestrial laser scanning and photogrammetry: Semi-automatic approaches for linear outcrop inspection. J. Struct. Geol. 2014, 66, 102–114. [Google Scholar] [CrossRef]
  30. Lato, M.J.; Vöge, M. Automated mapping of rock discontinuities in 3D lidar and photogrammetry models. Int. J. Rock. Mech. Min. Sci. 2012, 54, 150–158. [Google Scholar] [CrossRef]
  31. Vöge, M.; Lato, M.J.; Diederichs, M.S. Automated rockmass discontinuity mapping from 3-dimensional surface data. Eng. Geol. 2013, 164, 155–162. [Google Scholar] [CrossRef]
  32. Slob, S.; van Knapen, B.; Hack, R.; Turner, K.; Kemeny, J. Method for automated discontinuity analysis of rock slopes with three-dimensional laser scanning. Transp. Res. Rec. 2005, 1913, 187–194. [Google Scholar] [CrossRef]
  33. Olariu, M.I.; Ferguson, J.F.; Aiken, C.L.; Xu, X. Outcrop fracture characterization using terrestrial laser scanners: Deep-water Jackfork sandstone at Big Rock Quarry, Arkansas. Geosphere 2008, 4, 247–259. [Google Scholar] [CrossRef]
  34. Chen, J.; Zhu, H.; Li, X. Automatic extraction of discontinuity orientation from rock mass surface 3D point cloud. Comput. Geosci. 2016, 95, 18–31. [Google Scholar] [CrossRef]
  35. Guo, J.; Liu, S.; Zhang, P.; Wu, L.; Zhou, W.; Yu, Y. Towards semi-automatic rock mass discontinuity orientation and set analysis from 3D point clouds. Comput. Geosci. 2017, 103, 164–172. [Google Scholar] [CrossRef]
  36. Sinaga, K.P.; Yang, M.S. Unsupervised K-Means Clustering Algorithm. IEEE Access 2020, 8, 80716–80727. [Google Scholar] [CrossRef]
  37. Yi, X.; Feng, W.; Wang, D.; Yang, R.; Hu, Y.; Zhou, Y. An efficient method for extracting and clustering rock mass discontinuities from 3D point clouds. Acta Geotech. 2023, 18, 3485–3503. [Google Scholar] [CrossRef]
  38. Wu, W.; Zhang, K.; Zhu, H. A fast automatic extraction method for rock mass discontinuity orientation using fast k-means++ and fast silhouette based on 3D point cloud. In IOP Conference Series: Earth and Environmental Science; Springer: New York, NY, USA, 2020; p. 570. [Google Scholar] [CrossRef]
  39. Daghigh, H.; Tannant, D.D.; Daghigh, V.; Lichti, D.D.; Lindenbergh, R. A critical review of discontinuity plane extraction from 3D point cloud data of rock mass surfaces. Comput. Geosci. 2022, 169, 105241. [Google Scholar] [CrossRef]
  40. Leng, X.; Xiao, J.; Wang, Y.J.T.P.R. A multi-scale plane-detection method based on the Hough transform and region growing. Photogramm. Rec. 2016, 31, 166–192. [Google Scholar] [CrossRef]
  41. Hu, L.; Xiao, J.; Wang, Y. Efficient and automatic plane detection approach for 3-D rock mass point clouds. Multimed. Tools Appl. 2020, 79, 839–864. [Google Scholar] [CrossRef]
  42. Yu, D.; Xiao, J.; Wang, Y. High-Precision Plane Detection Method for Rock-Mass Point Clouds Based on Supervoxel. Sensors 2020, 20, 4209. [Google Scholar] [CrossRef]
  43. Lato, M.; Kemeny, J.; Harrap, R.M.; Bevan, G. Rock bench: Establishing a common repository and standards for assessing rockmass characteristics using LiDAR and photogrammetry. Comput. Geosci. 2013, 50, 106–114. [Google Scholar] [CrossRef]
  44. Zhu, Y.; Ting, K.M.; Carman, M.J. Density-ratio based clustering for discovering clusters with varying densities. Pattern Recognit. 2016, 60, 983–997. [Google Scholar] [CrossRef]
  45. Haris, K.; Efstratiadis, S.N.; Maglaveras, N.; Katsaggelos, A.K. Hybrid image segmentation using watersheds and fast region merging. IEEE Trans. Image Process. 1998, 7, 1684–1699. [Google Scholar] [CrossRef] [PubMed]
  46. Beucher, S. Use of watersheds in contour detection. In Proceedings of the Proceedings International Workshop on Image Processing, Rennes, France, 17–21 September 1979; pp. 17–21. [Google Scholar]
  47. Beucher, S. Watersheds of functions and picture segmentation. In Proceedings of the ICASSP’82. IEEE International Conference on Acoustics, Speech, and Signal Processing, Paris, France, 3–5 May 1982; pp. 1928–1931. [Google Scholar]
  48. Chien, S.-Y.; Huang, Y.-W.; Chen, L.-G. Predictive watershed: A fast watershed algorithm for video segmentation. IEEE Trans. Circuits Syst. Video Technol. 2003, 13, 453–461. [Google Scholar] [CrossRef]
  49. Guo, Q.; Wang, Y.; Yang, S.; Xiang, Z. A method of blasted rock image segmentation based on improved watershed algorithm. Sci. Rep. 2022, 12, 7143. [Google Scholar] [CrossRef]
  50. Demantké, J.; Mallet, C.; David, N.; Vallet, B. Dimensionality based scale selection in 3D lidar point clouds. Remote Sens. Spat. Inf. Sci. 2012, 38, 97–102. [Google Scholar] [CrossRef]
  51. Weinmann, M.; Jutzi, B.; Hinz, S.; Mallet, C. Semantic point cloud interpretation based on optimal neighborhoods, relevant features and efficient classifiers. ISPRS J. Photogramm. Remote Sens. 2015, 105, 286–304. [Google Scholar] [CrossRef]
  52. Vincent, L.; Soille, P. Watersheds in digital spaces: An efficient algorithm based on immersion simulations. IEEE Trans. Pattern Anal. Mach. Intell. 1991, 13, 583–598. [Google Scholar] [CrossRef]
  53. Riquelme, A.J.; Abellán, A.; Tomás, R.; Jaboyedoff, M. A new approach for semi-automatic rock mass joints recognition from 3D point clouds. Comput. Geosci. 2014, 68, 38–52. [Google Scholar] [CrossRef]
  54. Fischer, M. Randam sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. Commun. ACM 1981, 24, 381–385. [Google Scholar] [CrossRef]
  55. Raguram, R.; Chum, O.; Pollefeys, M.; Matas, J.; Frahm, J.-M. USAC: A universal framework for random sample consensus. IEEE Trans. Pattern Anal. Mach. Intell. 2012, 35, 2022–2038. [Google Scholar] [CrossRef] [PubMed]
  56. Raguram, R.; Frahm, J.-M.; Pollefeys, M. A comparative analysis of RANSAC techniques leading to adaptive real-time random sample consensus. In Proceedings of the Computer Vision–ECCV 2008: 10th European Conference on Computer Vision, Marseille, France, 12–18 October 2008; pp. 500–513. [Google Scholar]
  57. Nguyen, A.; Le, B. 3D point cloud segmentation: A survey. In Proceedings of the 2013 6th IEEE Conference on Robotics, Automation and Mechatronics (RAM), Manila, Philippines, 12–15 November 2013; pp. 225–230. [Google Scholar]
  58. Xu, B.; Jiang, W.; Shan, J.; Zhang, J.; Li, L. Investigation on the Weighted RANSAC Approaches for Building Roof Plane Segmentation from LiDAR Point Clouds. Remote Sens. 2016, 8, 5. [Google Scholar] [CrossRef]
  59. Ferrero, A.M.; Forlani, G.; Roncella, R.; Voyat, H.I. Advanced Geostructural Survey Methods Applied to Rock Mass Characterization. Rock Mech. Rock. Eng. 2009, 42, 631–665. [Google Scholar] [CrossRef]
  60. Chen, N.; Cai, X.; Li, S.; Zhang, X.; Jiang, Q. Automatic extraction of rock mass discontinuity based on 3D laser scanning. Q. J. Eng. Geol. Hydrogeol. 2020, 54, qjegh2020–qjegh2054. [Google Scholar] [CrossRef]
  61. Liu, L.; Xiao, J.; Wang, Y. Major Orientation Estimation-Based Rock Surface Extraction for 3D Rock-Mass Point Clouds. Remote Sens. 2019, 11, 635. [Google Scholar] [CrossRef]
Figure 1. Flowchart of the presented method.
Figure 1. Flowchart of the presented method.
Remotesensing 16 03291 g001
Figure 2. Overview of rock mass photos and the corresponding point clouds: (a) photo of rock surface in Ouray, CO, USA (from Rock bench repository); (b) corresponding point clouds of the area of study.
Figure 2. Overview of rock mass photos and the corresponding point clouds: (a) photo of rock surface in Ouray, CO, USA (from Rock bench repository); (b) corresponding point clouds of the area of study.
Remotesensing 16 03291 g002
Figure 3. The impact of different bandwidths on MS: (a) bandwidth = 6; (b,c) bandwidth ≈ 0.
Figure 3. The impact of different bandwidths on MS: (a) bandwidth = 6; (b,c) bandwidth ≈ 0.
Remotesensing 16 03291 g003
Figure 4. The effect of bandwidth on different types of discontinuities. (a) schematic representation illustrating flat discontinuities; (b) model of fitted planes for flat discontinuities; (c) The segmentation effect of flat discontinuities by MS with different bandwidths. (d) schematic representation illustrating non-planar discontinuities; (e) model of fitted planes for non-planar discontinuities (Side view); (f) The segmentation effect of non-planar discontinuities by MS with different bandwidths. ①, ②, and ③ represent three discontinuous surfaces in the real scenario, with different colors representing the positions of the identified discontinuous surfaces.
Figure 4. The effect of bandwidth on different types of discontinuities. (a) schematic representation illustrating flat discontinuities; (b) model of fitted planes for flat discontinuities; (c) The segmentation effect of flat discontinuities by MS with different bandwidths. (d) schematic representation illustrating non-planar discontinuities; (e) model of fitted planes for non-planar discontinuities (Side view); (f) The segmentation effect of non-planar discontinuities by MS with different bandwidths. ①, ②, and ③ represent three discontinuous surfaces in the real scenario, with different colors representing the positions of the identified discontinuous surfaces.
Remotesensing 16 03291 g004
Figure 5. The principle of the Watershed algorithm.
Figure 5. The principle of the Watershed algorithm.
Remotesensing 16 03291 g005
Figure 6. The segmentation process of the Watershed algorithm [49]: (a) a limestone blasted grayscale image; (b) 3D simulation of water catchment basin; (c) 3D simulation after flooding; (d) 2D cross-section after segmentation.
Figure 6. The segmentation process of the Watershed algorithm [49]: (a) a limestone blasted grayscale image; (b) 3D simulation of water catchment basin; (c) 3D simulation after flooding; (d) 2D cross-section after segmentation.
Remotesensing 16 03291 g006
Figure 7. Local features of the point cloud: (a) rod-shaped; (b) planar; (c) three-dimensional surface.
Figure 7. Local features of the point cloud: (a) rod-shaped; (b) planar; (c) three-dimensional surface.
Remotesensing 16 03291 g007
Figure 8. The normal vectors of surfaces with different degrees of curvature: (a) smooth area; (b) sharp area. The dashed circles represent the local neighborhood of the sampling points, and the arrows represent the normal directions.
Figure 8. The normal vectors of surfaces with different degrees of curvature: (a) smooth area; (b) sharp area. The dashed circles represent the local neighborhood of the sampling points, and the arrows represent the normal directions.
Remotesensing 16 03291 g008
Figure 9. An HSV simulation diagram of K .
Figure 9. An HSV simulation diagram of K .
Remotesensing 16 03291 g009
Figure 10. Comparison of flooding results for different watershed lines: (a) w = 13; (b) w = 12; (c) w = 11; (d) w = 10; (e) w = 9; (f) w = 8; (g) w = 7; (h) w = 6. Red color represents the slope peaks, and blue color represents the slope valleys.
Figure 10. Comparison of flooding results for different watershed lines: (a) w = 13; (b) w = 12; (c) w = 11; (d) w = 10; (e) w = 9; (f) w = 8; (g) w = 7; (h) w = 6. Red color represents the slope peaks, and blue color represents the slope valleys.
Remotesensing 16 03291 g010
Figure 11. Schematic overview of MTMS: (a) spatial distribution of normal vectors; (b) mark the non-planar area; (c) drift process.
Figure 11. Schematic overview of MTMS: (a) spatial distribution of normal vectors; (b) mark the non-planar area; (c) drift process.
Remotesensing 16 03291 g011
Figure 12. Results of recognizing a single discontinuity cluster: (a) first set of discontinuities; (b) grouping data points that form an individual plane. Different colors represent the identified discontinuous surfaces.
Figure 12. Results of recognizing a single discontinuity cluster: (a) first set of discontinuities; (b) grouping data points that form an individual plane. Different colors represent the identified discontinuous surfaces.
Remotesensing 16 03291 g012
Figure 13. Segmentation comparison with examples of over-segmentation for experimental dataset (the raw point cloud was provided by Lato et al. [43]): (a) Kong et al. [25]; (b) Chen et al. [60]; (c) Daghigh et al. [22] (the circled and rectangular areas represent undulations or minor offsets that remain); (d) The results of the present study for determining the main discontinuity sets.
Figure 13. Segmentation comparison with examples of over-segmentation for experimental dataset (the raw point cloud was provided by Lato et al. [43]): (a) Kong et al. [25]; (b) Chen et al. [60]; (c) Daghigh et al. [22] (the circled and rectangular areas represent undulations or minor offsets that remain); (d) The results of the present study for determining the main discontinuity sets.
Remotesensing 16 03291 g013
Figure 14. Discontinuity plane numbers used for comparison and validation (different colors represent different groups).
Figure 14. Discontinuity plane numbers used for comparison and validation (different colors represent different groups).
Remotesensing 16 03291 g014
Figure 15. Overview of rock mass photos and the corresponding point clouds: (a) photo of the study area for the first case; (b) point cloud data C1 for the first case; (c) photo of the study area for the second case; (d) point cloud data C2 for the second case.
Figure 15. Overview of rock mass photos and the corresponding point clouds: (a) photo of the study area for the first case; (b) point cloud data C1 for the first case; (c) photo of the study area for the second case; (d) point cloud data C2 for the second case.
Remotesensing 16 03291 g015
Figure 16. Ground truth of planar surfaces for C1 (each number represents a discontinuous surface, and non-planar areas are shown in black).
Figure 16. Ground truth of planar surfaces for C1 (each number represents a discontinuous surface, and non-planar areas are shown in black).
Remotesensing 16 03291 g016
Figure 17. Results of identifying non-planar areas in C1 (red regions represent non-planar areas).
Figure 17. Results of identifying non-planar areas in C1 (red regions represent non-planar areas).
Remotesensing 16 03291 g017
Figure 18. Comparison of segmentation results for C1 using different methods: (a) Ground truth of planar surfaces for C1; (b) RHT; (c) Region-growing: (d) RANSAC; (e) Optimized RANSAC; (f) HT-RG; (g) DSE; (h) MOE; (i) Supervoxel-RG; (j) reference [22]; (k) our method. The red circles and red polygons represent the shortcomings in each outcome, the black rectangles represent gaps in non-planar areas, and the white rectangles represent incorrect segmentation results.
Figure 18. Comparison of segmentation results for C1 using different methods: (a) Ground truth of planar surfaces for C1; (b) RHT; (c) Region-growing: (d) RANSAC; (e) Optimized RANSAC; (f) HT-RG; (g) DSE; (h) MOE; (i) Supervoxel-RG; (j) reference [22]; (k) our method. The red circles and red polygons represent the shortcomings in each outcome, the black rectangles represent gaps in non-planar areas, and the white rectangles represent incorrect segmentation results.
Remotesensing 16 03291 g018
Figure 19. The recognition results for non-planar areas.
Figure 19. The recognition results for non-planar areas.
Remotesensing 16 03291 g019
Figure 20. Segmentation results for C2 with different bandwidths (different colors represent the distinct discontinuous surfaces identified).
Figure 20. Segmentation results for C2 with different bandwidths (different colors represent the distinct discontinuous surfaces identified).
Remotesensing 16 03291 g020
Figure 21. The characterization results of C2: (a) discontinuity surface recognition results; (b) spherical projection of the characterization results; (c) the 2D profiles of the two non-planar discontinuities.
Figure 21. The characterization results of C2: (a) discontinuity surface recognition results; (b) spherical projection of the characterization results; (c) the 2D profiles of the two non-planar discontinuities.
Remotesensing 16 03291 g021
Figure 22. The sample plane marked with different colors for verification and their corresponding numbers (the first digit of the serial number represents the associated discontinuous group).
Figure 22. The sample plane marked with different colors for verification and their corresponding numbers (the first digit of the serial number represents the associated discontinuous group).
Remotesensing 16 03291 g022
Table 1. Symbols and notations.
Table 1. Symbols and notations.
x The current sampling point
D , S a d-dimensional dataset with n instances
ε , η radius of neighborhood
ψ number of intervals for histogram
p d f ( · ) probability density function
p d f η ^ ( · ) the average density of the sampling point within the neighborhood η
p d f ε ^ ( · ) the average density of the sampling point within the neighborhood ε
c d f η ^ ( · ) The CDF of the sampling point within the neighborhood η
r p d f ^ ε , η ( · ) density-ratio estimation
Table 2. Resample the dataset D into S using the cumulative distribution function CDF.
Table 2. Resample the dataset D into S using the cumulative distribution function CDF.
Dataset D (Original Data) S (Resampled Data)
Instances x y = c d f η ^ ( x )
Density estimator p d f ε ^ ( x ) p d f ε ^ ( y ) = r p d f ^ ε , η ( x )
Table 3. Comparison of measurement results.
Table 3. Comparison of measurement results.
Plane No.Dip Direction/Dip Angle (°)
Truth Value (Daghigh et al. (2022)) [39]Daghigh et al. (2023) [22]Our MethodError Degrees
TM1M2T vs. M1T vs.M2
1246.1/39.1244.93/39.39244.07/38.541.17/0.290.86/0.85
2248.1/49.5252.77/50.32243.77/50.454.67/0.824.33/0.95
3249.8/35.7248.95/35.58249.90/36.060.85/0.120.10/0.36
4251.8/34.7251.85/34.58251.03/34.070.05/0.120.77/0.63
5249.3/35.4252.89/36.99248.77/35.233.59/1.590.53/0.17
6250.1/35.7250.45/35.85250.81/34.950.35/0.150.71/0.75
7253.4/33.3251.22/35.30253.47/33.364.18/2.000.07/0.06
8339.5/83.1339.67/83.27338.2/83.610.17/0.171.30/0.50
9345.7/73.2345.19/76.07163.66/76.120.51/2.872.04/2.92
10338.2/88.1341.02/87.59338.05/86.832.82/0.510.15/1.27
11173.7/76.7176.47/75.15174.18/76.232.77/1.550.48/0.47
12136.7/77.0134.40/79.82315.12/78.542.30/2.821.65/1.54
13147.2/89.1143.66/89.51147.19/89.213.54/0.410.01/0.11
1497.2/62.299.55/61.63101.01/64.732.35/0.573.81/2.53
1592.5/48.692.28/48.9993.13/47.520.22/0.390.63/1.08
1697.7/48.897.40/48.1797.57/47.820.30/0.630.13/0.98
17304.4/77.9307.35/73.42304.33/79.632.95/4.481.22/1.73
18286.7/70.7286.72/69.86286.41/69.300.02/0.840.29/1.40
Average 1.82/1.691.06/1.01
Table 4. Comparison of the measurement results between our method and CloudCompare.
Table 4. Comparison of the measurement results between our method and CloudCompare.
Plane No.Dip Direction/Dip Angle (°)
Our MethodCloudCompareError Degrees
1114.9/37.515.3/37.10.4/0.4
1212.2/40.112.1/40.50.1/0.4
21331.1/61.7329.2/59.51.9/2.2
22318.6/66.3320.1/67.91.5/1.6
31178.5/77.8176.5/76.32.0/1.5
32187.1/78.1185.9/78.91.2/0.8
41358.0/85.5358.3/84.10.3/1.4
42318.4/80.6316.6/79.91.8/0.7
Average 1.2/1.1
Table 5. Run-time comparison.
Table 5. Run-time comparison.
ReferencesInputProgrammingComputer SystemTime (min)
Riquelme et al. (2014) [54]3D pointsMATLABIntel Core i3-350 M and 8 GB DDR3 RAM88
Chen et al. (2016) [36]2.5D pointsNot reportedIntel Core i7-2600 CPU and 16 GB RAM150
Kong et al. (2020) [1]3D pointsNot reportedIntel Core i5 6300Q and 4 GB RAM330
Daghigh et al. (2023) [3]3D pointsPythonIntel Core i7-8565U CPU and 16 GB RAM4
Our method3D pointsPythonIntel Core i7-1195G7 CPU and 16 GB RAM6
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Liu, Y.; Hua, W.; Chen, Q.; Liu, X. Characterization of Complex Rock Mass Discontinuities from LiDAR Point Clouds. Remote Sens. 2024, 16, 3291. https://doi.org/10.3390/rs16173291

AMA Style

Liu Y, Hua W, Chen Q, Liu X. Characterization of Complex Rock Mass Discontinuities from LiDAR Point Clouds. Remote Sensing. 2024; 16(17):3291. https://doi.org/10.3390/rs16173291

Chicago/Turabian Style

Liu, Yanan, Weihua Hua, Qihao Chen, and Xiuguo Liu. 2024. "Characterization of Complex Rock Mass Discontinuities from LiDAR Point Clouds" Remote Sensing 16, no. 17: 3291. https://doi.org/10.3390/rs16173291

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop