Intelligent Interpretation of the Geometric Properties of Rock Mass Discontinuities Based on an Unmanned Aerial Vehicle

The geometric properties of rock mass discontinuities are essential for the evaluation of the safety of rock masses. Numerous studies have recently been performed on the extraction of discontinuity information. However, most methods are characterized by poor data collection and processing efficiency. This paper presents a UAV-based methodology for the accurate and complete acquisition of rock surface data, as well as the automatic extraction of discontinuity information. Moreover, a program called Random Sample Consensus (RANSAC) Discontinuity Detection (RDD) is developed to extract discontinuity information based on the proposed method. The conclusions of this research are as follows. 1) RANSAC Discontinuity Detection (RDD) can identify the feature point set of discontinuities from a raw point cloud, and can calculate the discontinuity orientation. 2) The boundary of a discontinuity can be precisely depicted using the improved Graham scan algorithm. 3) The orientations of marked discontinuities extracted by RDD are compared with those extracted by the three-point method in CloudCompare. The differences in the orientations extracted by the two methods are found to be less than 3° for flat discontinuities and only about 4.87° for rough discontinuities, which are within a reasonable error range in practical engineering applications. Therefore, the feasibility of the proposed method is verified.


INTRODUCTION
Research indicates that discontinuities are an intrinsic characteristic of rock masses (Umili et al., 2013), and they have significant influences on rock mass deformation and stability (Kong et al., 2020). Therefore, the accurate and comprehensive extraction of rock mass discontinuity information is critical for the assessment of the safety of rock masses.
Traditional surveys are conducted via a window method or a line-scanning method (Gigli and Casagli, 2010;Kong et al., 2020), and require physical contact with the rock surface (Gigli and Casagli, 2010;Umili et al., 2013); however, this is time-consuming Kong et al., 2020) and subject to the expertise of the operator (Kong et al., 2020). With the advancement of measurement techniques, new non-contact surveying methods have been developed to acquire three-dimensional (3D) rock mass data, and include the total station method (Feng et al., 2001), close-range photogrammetry (De et al., 2012;Kaufmann, 2012;Francioni et al., 2019), and 3D laser scanning (Deliormanli et al., 2014;Monsalve et al., 2019;Wichmann et al., 2019;Jiang et al., 2020). These techniques have been rapidly utilized in slope monitoring (Kromer et al., 2019;Giacomini et al., 2020), rock mechanics and stability analysis (Firpo et al., 2011;Assali et al., 2014), and geomorphology (Brodie et al., 2015;Boothroyd et al., 2016), as well as the geological and geotechnical research fields (Giordan et al., 2018). While 3D laser scanning and close-range photogrammetry have made great development progress, these two survey means are characterized by the following disadvantages: 1) under unique and complex terrain conditions, it is difficult to find a suitable observation point at which to set up instruments; 2) scanning devices are expensive, and the prices of mainstream scanners on the market are more than 1 million yuan; 3) some rock mass data will inevitably be missing due to the scanning direction. Due to the conspicuous cost reduction of vehicles and sensors, as well as the recent progress in data processing software over the past decade (Manfreda et al., 2018), the application of unmanned aerial vehicles (UAVs) in the collection of rock mass information has been ensured. UAV photogrammetry has many benefits, such as light and flexible equipment, strong adaptability to various terrains, and wide coverage . In addition, a multi-rotor UAV can take images from different positions and in different directions, thereby avoiding possible shadows or vertical deviations in high and inaccessible rock surfaces (Salvini et al., 2020). Therefore, the use of a light and small UAV to acquire rock mass information is superior to 3D laser scanning and close-range photogrammetry in terms of the equipment cost, portability, efficiency, and integrity of data collection.
The UAV technique has been applied in slope monitoring (Rodriguez et al., 2020;, the failure mechanism analysis of landslides (Xu et al., 2017;, photogrammetric inspection , and topographic reconstruction (Agüera-Vega et al., 2018). However, it is rarely utilized in the extraction of rock mass discontinuity information. Yathunanthan et al. (2014) conducted the imaging analysis of a data set generated from UAV photography to map geological structures, after which the 3D feature coordinates corresponding to the pixel coordinates of two-dimensional (2D) feature points were calculated from the digital elevation model (DEM), and the best-fit plane coefficients were computed. Finally, the discontinuity orientation was extracted (Yathunanthan et al., 2014). However, in this method, the geological analysis of the images requires the acute intuition and deductive and inductive reasoning of interpreters. Jia et al. (2018) manually selected the exposed discontinuities from a 3D point cloud generated from UAV images, and then extracted discontinuity information using the least-squares plane-fitting algorithm (Jia et al., 2018). It is evident that the degree of automation of these methods remains to be increased.
To overcome the low efficiency of the existing methods for rock mass data acquisition and the low degree of automation in discontinuity identification, a UAV-based approach for the automatic identification of rock mass discontinuities is proposed in this paper. Efficient and comprehensive image acquisition can be realized by employing a multi-rotor UAV, and images collected by the UAV can be converted to a point cloud model of the rock mass via 3D model reconstruction. Moreover, rock mass discontinuities can be automatically recognized based on an improved random sample consensus (RANSAC) algorithm. The proposed method includes the following steps: 1) 3D model reconstruction; 2) Normal vector calculation; 3) Discontinuity extraction; 4) Boundary delineation; 5) Discontinuity orientation calculation.

UNMANNED AERIAL VEHICLES MEASUREMENT SYSTEM AND WORKFLOW
The survey was conducted using a four-rotor DJI Phantom 3 Professional UAV (Figure 1), which is equipped with a 20-mm low-distortion wide-angle camera, a GPS/GLONASS dual-mode system, and an automatic return function. Moreover, the aircraft supports 4 K video capture at 30 frames per second, and the two photo formats of JPG and Raw. In addition, the aircraft has the ability to capture smooth and stable video pictures, and can actively record all the details of each flight, including the course, flight time, and other information. The specific parameters of the UAV are listed in Table 1.  A scene survey is the foremost step of, and lays a solid foundation for, data acquisition. After surveying the site conditions, the measurement range was determined, and the easily measured discontinuities were selected and marked. The DJI GO App, a mobile route-planning software, was used to plan the route. The flight path should exceed the survey range to shoot the whole area. The flying height was selected to be beyond the highest obstacle in the flying area to avoid a collision. It should be noted that it is necessary to balance the flying height with the flight time, as the battery consumption caused by an excessive flying height will exceed the battery capacity. Ultimately, the slope was photographed from left to right and from top to bottom. The flight speed was selected as the maximum value, and the camera angle was set to 45°. To maintain the integrity of the measured data, some overlaps were maintained between the images. An overlap of 80% and a side-lap of 50% were respectively used. After the field investigation, a 3D point cloud was generated from the photos taken by the UAV.

MATERIALS AND METHODS
The proposed method is divided into five steps, as illuminated in the flowchart in Figure 2.

Feature Point Extraction
Due to the strong distortion of photos taken by UAVs, it is difficult to effectively apply the traditional extraction method based on geometric and texture features to extract feature points. The scale-invariant feature transform (SIFT) algorithm is characterized by the three properties of scaling, rotation, and Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 711866 3 affine invariance, via which it is capable of resisting certain illumination changes and viewpoint transformation. Hence, the SIFT algorithm is adopted for feature point extraction in the proposed method. The main concept of the SIFT algorithm is that the scale-space representation of UAV aerial images is established, after which the extreme points of images are searched in the scale space and extracted as feature points.

Image Matching
Image matching is conducted to reconstruct 3D information from multiple 2D images. However, the process of image matching using only SIFT feature points is slow. The location data of GPS coordinates in the images collected by the UAV and the attitude angle data provided by the inertial measurement unit (IMU) can assist in the construction of the topological structure between images. Next, the nearest-neighbor method is utilized to find the corresponding relationships between the feature points of images and establish a set of matching feature points that meet the geometric constraints. A large number of coordinate points constitutes a 3D point cloud of the target object in space.

Structure From Motion
The image points in the photo are projected into spatial coordinates according to the principle of camera imaging. The error function is defined as the sum of squares of the reprojection errors. The objective function is defined as follows: where C p {C 1 , C 2 , C 3 , ...C n } are camera parameters, X {X 1 , X 2 , X 3 , ...X m } are the coordinates of space points, the variable v ij represents the visibility of space point X i in camera C i , n is the number of images, m is the number of feature points obtained by precise matching, and the function f (P(C i , X j ), q ij ) 2 represents the projection error of point X j in camera C i .
Finally, sparse beam adjustment is used for step-by-step iteration to minimize the reprojection error between the projected points and the points on the observed images, thereby calculating coordinates of the 3D point cloud in an optimal camera pose and camera scene.

Normal Vector Calculation
Considering that normal vector calculation is a necessary process for the extraction of discontinuities, the next step after obtaining the point cloud is to calculate the normal vector.
The normal vector calculation method consists of two key steps, namely 1) finding the k-nearest neighbors of each point P i and creating the point set Q i (Figure 3), and 2) conducting plane-fitting for each point set Q i and calculating the surface variation.

Nearest-Neighbor Search
A point cloud model commonly includes massive target points in a 3D region, and lacks topological information. Therefore, the principal problem of processing point cloud data is to establish a topological relationship among discrete points and realize the fast search of the nearest adjacent points.
A k-d tree (referred to as a k-dimensional tree) is a data structure that represents spatial partitions, and is mainly applied to search key data in multidimensional space (such as range searches and nearest-neighbor searches). In this study, the k-nearest neighbors are searched by the k-dimensional tree.

Surface Variation Calculation
The problem of determining the normal of a point on a surface is similar to the problem of estimating a section of a normal of a surface. Therefore, the problem can be transformed into one concerned with least-squares plane-fitting estimation. In this study, the surface normal is evaluated by analyzing the eigenvectors and eigenvalues (or principal component analysis, PCA) (Riquelme et al., 2014;Robson et al., 2016;Guo et al., 2017) of the covariance matrix created from the nearest points. The covariance matrix C corresponding to each point P i can be defined as where k is the number of the adjacent points of P i , P represents the 3D centroid of point sets, and λ j and V j → are the eigenvalue and eigenvector of the covariance matrix, respectively. The normal vector can be determined by the eigenvector that corresponds to its minimum eigenvalue. After all the normal vectors are obtained, the subsequent step is to extract rock mass discontinuities.

Rock Mass Discontinuity Extraction Based on the Improved Random Sample Consensus Algorithm
The RANSAC algorithm is an iterative computational algorithm, which determines the parameters of a predefined mathematical FIGURE 3 | Q i is the subset of P i , and α is the normal vector of P i . Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 711866 model by randomly selecting a subset sample, and then calculates the distance from all points to the model. A point is defined as an inlier if the distance between a point and the model is less than the threshold; otherwise, the point is regarded as an outlier. The number of inliers in each iteration is recorded, and the model with the largest number of inliers is considered as the optimal model (Uhercík et al., 2010). Figure 4 presents the simple application of the RANSAC algorithm to 2D data. Figure 4A displays a set of points that includes inliers and outliers, and in Figure 4B, the blue line is the generated mathematical model, and the inliers are indicated in red. For rock discontinuity extraction (the mathematical model is a plane), the RANSAC algorithm has two advantages: 1) it can be directly applied to raw point cloud data without triangulation gridding, and 2) it has strong robustness and can process more than 50% of the outliers. Based on these advantages, the RANSAC algorithm has been studied for the extraction of planes or discontinuities (Wang et al., 2019). However, due to the large number of points of a rock mass, most approaches are inefficient. Therefore, an improved RANSAC algorithm is proposed to greatly improve the accuracy and speed of the original algorithm.

Overview of the Improved Random Sample Consensus Algorithm
Given a point cloud set P {P 1 , ..., P N } and the normal vector N {n 1 , ...n N } of all points, the output result is a series of parameters ψ {ψ 1 , ...ψ N } of the plane model. In this paper, local sampling is proposed to acquire new candidate planes in each iteration. The RANSAC algorithm is then applied to determine the parameters of the plane model with the highest score (i.e., the largest number of inliers). All the candidate planes are placed in the set C, and a new evaluation method is used to calculate the score m of the best plane. Moreover, |m| is the number of points in a candidate plane, |c| is the number of candidate planes, and p(|m|, |c|) is the probability of ignoring the planes with a higher score. When p(|m|, |c|) is large enough, the extracted plane is the best, and the remaining points will be used for the subsequent iteration. When p(τ, |c|) is sufficiently large, the iteration process is terminated. Finally, τ (default) is the number of minimum points on a plane.

Probability Calculation
Consider a point cloud with N points and a plane with n points, and k is the number of points in the minimum point cloud set that determines a plane. Provided that any subset with k points will generate a plane model, then the probability of detecting the plane model in one iteration is as follows.
When s candidate planes are detected, the probability of detecting the plane Ψ is as follows. P(n, s) 1 − (1 − P(n)) s The threshold value p t is artificially set. The number of planes T that meet the requirement P(n, T) ≥ p t can then be obtained by solving s, as follows.
Because the value of P(n) is commonly small, its logarithm ln(P(n)) can be expanded by its Taylor series: ln(1 − P(n)) −P(n) + O(P(n) 2 ).
The substitution of ln(P(n)) into Eq. 5 yields the following.

Sampling Method
The complexity of an algorithm is closely related to the sampling method. The sampling method used in this study is detailed as follows. Shape is a local phenomenon, and the closer two points are, the more likely they belong to the same plane. The sampling efficiency can be greatly improved by utilizing this characteristic.
Research has revealed that it may be effective to increase the number of inliers within the model by utilizing the locality of the shape to the sample (Myatt et al., 2002). Generally, in random sampling, a circle with a given radius is used to randomly select sample points, but the radius needs to be determined in advance according to the density and distribution of the points. However, the density and distribution of outliers vary greatly for different models; even at different locations on the same model, the density of outliers can vary dramatically. Therefore, a method is presented in this paper to adapt to the density of outliers. The octree structure is an effective method by which to establish spatial proximity between sampling points. First, the point p 1 is chosen without restriction to create the candidate plane, and then a set C that includes p 1 is randomly selected from the structural layers of the octree structure. Finally, the remaining K − 1 sample points are selected from the set C. The probability of finding the plane ψ containing n points in this manner can be calculated as follows: where n/N is the first probability value, and the second probability value relies on the choice of the elements in the set C. The set C is considered to be an optimization if abundant points on the plane ψ are included in the set. Most points on a plane, excluding boundary points and edge points, have neighbors that belong to the plane. Generally speaking, although the adjacent points on a plane cannot be determined by the element set of the octree structure, the candidate plane model includes a large number of points to ensure that it is more representative of the actual data. Therefore, the number of these neighbors must be as large as the number of elements in the octree structure, excluding a few points. To facilitate the analysis, it is assumed that the set C embodies all points p i on the plane ψ (p i ⊂ ψ), the number of all points on the plane ψ is half of the number of points in the set C, and the other half of the points in the set C contains outliers or noise points. The probability of Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 711866 choosing a large set C is conservatively evaluated to be 1/d, where d refers to the depth of the octree structure. Therefore, the conditional probability of choosing points p 2 and p 3 (p 2 , p 3 ⊂ ψ) from the set C set can be calculated by Eq. 8. Then, by introducing Eq. 8 into Eq. 7 and Eq. 9 can be obtained.

Evaluation Method
The evaluation function σ p is utilized to evaluate the extracted candidate plane model. This evaluation function mainly includes the following three aspects. 1) After the selected plane is created, the points whose distance to the plane is less than the distance threshold value ε are regarded as the points on the candidate plane.
2) The points that meet the distance requirement will be further filtered. When the angles between the normal vectors of the detected points and the normal vectors of the candidate planes are less than the angle threshold value α, these points will be selected as the points in the plane. 3) A new threshold value β, which represents continuity, is added to the proposed method. For the points that have met the first two requirements, only those satisfying the continuity requirement can be selected as inliers on the plane.
In short, for a plane ψ, its evaluation function σ p can be expressed as follows.
For example, points P ψ on a plane model ψ can be defined in two steps: where d(ψ, p) is the Euclidean distance from point p to plane ψ, n(p) is the normal vector of point p, n(ψ, p) is the relationship between the normal vector of the plane model ψ and the projection of the normal vector of point p on the plane ψ, and maxcomponent(ψ, p ψ ) refers to the point set where the projected points on the plane ψ can form the largest connected part.

Modified Graham Scan Algorithm
Rock mass discontinuities are generally not standard, and even very irregular planes. Hence, a new approach is presented to accurately depict the boundary of the discontinuities. Feature points of discontinuities extracted by the method in this paper are normally distributed on two sides of the vertical direction of extracted discontinuities. Thus, this 3D issue can be transformed into a problem of searching the optimal contour of a 2D point cloud by projecting points onto the fitting plane. The Graham scan algorithm is a straightforward and efficient convex hull algorithm (Ferrada et al., 2020), the general concept of which is to remove points that are not part of the convex hull. Given a point cloud set S, the point p 0 is obtained from S with the minimum y-coordinate, and the points are then sorted counterclockwise. By scanning from p 0 , if p 0 , p 1 , and p 2 are on the convex surfaces, they must meet the following property: p 2 is on the left of the vector 〈p 1 , p 2 〉; otherwise, p 1 should be removed.
However, the disadvantage of the Graham scan algorithm is that the first boundary obtained by the algorithm is convex. Therefore, optimization should be conducted to tackle this issue. Given the boundary points M k and nonboundary points N j , and provided that there are the two new lines M 1 − N 1 and N 1 − M 2 , N 1 will be regarded as a new boundary point if the triangle formed by the two new lines and M 1 − M 2 does not include the new boundary point and the angle formed by M 1 − N 1 , and N 1 − M 2 conforms to the tolerance requirement of a concave angle. In addition, the two new sides can be used to calculate the new boundary points by utilizing a recursive algorithm. The results of the Graham scan algorithm and the modified Graham scan algorithm are respectively presented in Figures 5A,B.
After calculation, all the boundary points with concaveconvex features are acquired. Figure 6 displays the extraction results of a set of points; Figure 6A presents the result of the fitting plane, and Figure 6B depicts the boundary of the fitting plane detected by the improved Graham scan algorithm.

Orientation Calculation
The discontinuity of a plane can be represented as follows.
Ax + By + Cz + D 0 The normal vector of the plane is as follows.
The orientation of the discontinuity can be calculated by the following equations.When N z > 0, If N x ≥ 0, When N z < 0, CASE STUDY

Data Description
The study site was a rock mass slope located in Guishan Park in Wuhan City, China. The first task was to investigate the field   Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 711866 7 conditions to determine the survey range and mark the easily measured discontinuities. The course was then planned, and the slope was photographed. Ultimately, twenty-six consecutive images were collected. Each image taken by the UAV was attached with geotags to provide a geographic reference for the point cloud, thereby producing highly detailed landform information (Rodriguez et al., 2020).
Subsequently, images from the survey were used to generate a point cloud. The 3D model ( Figure 7) and point cloud data of the slope were automatically generated via 3D reconstruction with Smart3D software. There were 52,639,008 points in the point cloud data, which rendered it difficult to extract discontinuities in the later step. The point cloud model of the slope was acquired (Figure 8) by setting the point spacing of 16 mm in CloudCompare to down-sample the numerous point clouds.

Discontinuity Information Results
A new extraction procedure called RANSAC Discontinuity Detection (RDD) was developed to automatically extract rock mass discontinuity information. The dialog box of RDD is presented in Figure 9, where ε is the distance from the point to the discontinuity, α is the angle between the normal vector of a point and the normal vector of the plane, and β is the distance between the points that make up a continuous plane.
The slope point cloud presented in Data description Section was employed to analyze the extraction effect of RDD. Figure 9 exhibits the parameter settings when RDD is used to extract discontinuities information of rock mass from the slope point cloud. As shown in Figure 10, most relatively smooth rock mass discontinuities were entirely detected by RDD, as were some small and fragmented discontinuities, even when the surface was rough or there were too many crushed pieces. It should be noted that the colors of these discontinuities were randomly assigned by the plug-in.
The results of the three-point method were compared with those of the proposed RDD to quantitatively analyze the effect of RDD in identifying discontinuities. Thirteen planes from the slope model were marked (Figure 11). Similarly, these 13 planes   were marked in the calculation results of RDD ( Figure 12). Table 2 reports the comparative results of the orientations extracted by the three-point method and RDD. As shown in the table, when the discontinuities were flat enough, the difference between the two methods was less than 3°. Moreover, even for the rough discontinuities (such as plane 12), the difference in orientation was only 4.87°, which is within the reasonable error range in practical engineering applications.

CONCLUSION
In this paper, a UAV-based approach was proposed for the data acquisition of rock masses and the automatic extraction of discontinuity information, which increases the automation level of discontinuity extraction and overcomes the disadvantages of close-range photogrammetry and 3D laser scanning, such as insufficient data.
Via the proposed method, discontinuity information can be extracted from the raw point cloud data, and the discontinuity boundary can be described with high precision.
A new procedure, RDD, was also developed to realize the automatic extraction of rock mass discontinuities. The orientations of the marked discontinuities calculated by the three-point method and the proposed RDD were compared, and the experimental results demonstrate the following: 1) RDD can completely detect most relatively smooth discontinuities, and therefore exhibits a good discontinuity detection effect; 2) the difference in the values of orientations calculated by the RDD and the manual three-point method was less than 3°when the discontinuity was smooth enough. Even for relatively rough discontinuities, the error was within an acceptable range for practical engineering applications. Therefore, the practicability of the proposed method was proven.
Nonetheless, this research was characterized by some disadvantages. For example, the extraction results of rock mass discontinuity information do not include the discontinuity spacing, roughness, degree of weathering, etc. Therefore, future research will be conducted to systematize the information extraction of rock mass discontinuities.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
NC proposed the algorithm and developed the RDD procedure. CD was responsible for data collection and curation, as well as wrote the manuscript and edited the charts. XD verified the experimental design, supervised and guided the process of data acquisition, as well as reviewed and revised the manuscript.