Automatic identification and characterization of discontinuities in rock masses from 3D point clouds

Journal Pre-proof Automatic identification and characterization of discontinuities in rock masses from 3D point clouds

Deheng Kong, Faquan Wu, Charalampos Saroglou PII:

S0013-7952(19)30584-8

DOI:

https://doi.org/10.1016/j.enggeo.2019.105442

Reference:

ENGEO 105442

To appear in:

Engineering Geology

Received date:

31 March 2019

Revised date:

21 November 2019

Accepted date:

29 November 2019

Please cite this article as: D. Kong, F. Wu and C. Saroglou, Automatic identification and characterization of discontinuities in rock masses from 3D point clouds, Engineering Geology (2019), https://doi.org/10.1016/j.enggeo.2019.105442

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier.

Journal Pre-proof

Automatic Identification and Characterization of Discontinuities in Rock Masses from 3D Point Clouds

Deheng Konga*, Faquan Wua, b, Charalampos Saroglouc

Department of Geotechnical Engineering, College of Civil Engineering, Tongji

oo

f

a

University, Shanghai 200092, China

Key Laboratory of Rock Mechanics and Geohazards of Zhejiang Province,

pr

b

Department of Geotechnics, School of Civil Engineering, National Technical

Pr

c

e-

Shaoxing University, Shaoxing 312000, China

rn

*Corresponding author.

al

University of Athens, Athens 15780, Greece

Jo u

E-mail addresses: [email protected] (D. Kong), [email protected] (F. Wu), [email protected] (C.Saroglou).

a b c

Deheng Kong developed the methodology, performed data analysis, and wrote the manuscript. Faquan Wu supervised the research and contributed to preparation of the manuscript. Charalampos Saroglou discussed the algorithms and results and critically reviewed the manuscript.

1

Journal Pre-proof

Highlights: 

An automatic method using several machine learning algorithms for identification and extraction of rock mass discontinuities is proposed.



Four discontinuity parameters, namely number of sets, orientation, spacing and trace length can be obtained. Discontinuity location, best fitting plane,

oo

This method is applied for two real cases and produces reliable and

rn

al

Pr

e-

pr

accuracy results.

Jo u



f

and 3D trace mapping can also be performed.

2

Journal Pre-proof

Abstract The routine application of remote surveying techniques which can quickly acquire 3D digital data with high resolution, in particular digital photogrammetry, light detection and ranging (LiDAR) and unmanned aerial vehicle (UAV) for rock mass characterization has rapidly grown over the past decade. In this paper, a new method for automatic identification and interpretation of rock mass and characterization of

f

of discontinuity sets

oo

discontinuities, clustering

discontinuity orientation, persistence and spacing using 3D point clouds, is

pr

presented. The proposed method is based on a four-stage procedure consisting

e-

of: (1) normal vector calculation using the iterative reweighted plane fitting (IRPF)

Pr

method, (2) discontinuity sets clustering by fast search and find of density peaks (CFSFDP) algorithm, and Fisher’s 𝐾 value iterative calculation to eliminate

al

noise points, (3) discontinuity segmentation using density-ratio based method,

rn

and discontinuity plane fitting using the random sample consensus (RANSAC)

Jo u

algorithm, (4) persistence and spacing calculation using the theory of analytic geometry. The method is applied to two case studies (i.e. rock slopes) and compared with the results from previous studies and from manual survey. It is concluded that the proposed method is reliable and yields a great accuracy for automatic identification of discontinuities in rock masses. Keywords: 3D point cloud, rock mass, automatic, discontinuity, LiDAR, UAV.

3

Journal Pre-proof

List of Symbols Point cloud data

𝑥, 𝑦, 𝑧

Three dimensional coordinates of a point

𝑁

Number of points

𝑝𝑖

The 𝑖th point at point cloud 𝑃

𝑁𝑘 (𝑝𝑖 )

The 𝑘 nearest neighbors of a point 𝑝𝑖

𝑃𝑙

The iterative reweighted fitting plane

𝑑

The distance from a point to the fitted plane

𝜎𝑑

Distance bandwidth

𝜎𝑟

Fitting residual bandwidth

𝜎𝑛

Normal difference bandwidth

𝐶𝑀

Positive semi-definite covariance matrix

𝑡

Number of iterations

𝒏

Normal vector

𝜃

Dip direction (°)

𝑙, 𝑚, 𝑛

oo

pr

e-

Pr

al

rn

Jo u

𝛿

f

𝑃

Dip angle (°)

Three components of a normal vector

𝑄

Angular converted value for dip direction calculation

𝑙𝑑

Local density of the pole data in CFSFDP algorithm

𝑑𝑚

The distance between pole data in CFSFDP algorithm

𝑚𝑑

The minimum distance between pole point and another higher

density pole point in CFSFDP algorithm 𝑑𝑐𝑓

Cutoff distance in CFSFDP algorithm

4

Journal Pre-proof

𝛼

The angular deviation from the mean orientation

𝐾

The Fisher Constant

𝑀

The total number of discontinuities for one set The mean orientation

𝑒𝑝𝑠

Radius of neighborhood

𝜂

Radius of the larger neighborhood

𝜏

Density ratio threshold

𝑚𝑖𝑛𝑠𝑖𝑧𝑒

Threshold fordiscarding small discontinuity clusters

𝑢𝑎,𝑢𝑏,𝑢𝑐

Three components of the unit normal vector to the fitted plane

𝐷𝐿

Discontinuity location

𝑠𝑧

The total number of one discontinuity cluster points

𝑇𝐿

Trace length

𝑒𝑑

The Euclidean distance of two points

𝐴, 𝐵, 𝐶, 𝐷

Normal unit vector (𝐴, 𝐵, 𝐶) and the constant parameter 𝐷 in

rn

al

Pr

e-

pr

oo

f

𝑟𝑀

the adjacent discontinuity plane equation Spacing value of two adjacent discontinuities

Jo u

𝑆𝑃

5

Journal Pre-proof

1. Introduction Discontinuities have an essential and pivotal role in engineering geology and rock mechanics, and have a significant influence on deformation, hydraulics and stability of rock masses. ISRM (1978) proposed a suggested method for the description of geometrical, mechanical and hydraulic features of discontinuities, consisting of the following parameters: orientation, spacing, persistence (or trace

oo

f

length), roughness, aperture, wall strength, filling, seepage, number of sets and block size. The accurate quantification of those discontinuity properties has

pr

always been considered as an initial key step for rock mechanics analyses. The

e-

traditional survey by using a geological compass and measuring tapes to apply

Pr

the scanline (one-dimensional) or window sampling method (bi-dimensional) is time-consuming and causes the user bias (Abellán et al., 2014; Gigli and

al

Casagli, 2011; Priest, 1993). The accuracy of statistical results by hand mapping

rn

relies on the selected rock mass area and the geological expertise of the user.

Jo u

Alternatively, remote sensing techniques such as digital photogrammetry, light detection and ranging (LiDAR, airborne and terrestrial lasering scanning) and unmanned aerial vehicle (UAV) with structure from motion (SfM) technique, which can acquire high-resolution 3D information (e.g. 3D point clouds), have been rapidly developed in geological and geotechnical research fields (Zekkos et al., 2018; Giordan et al., 2017), such as structural geology (Cawood et al., 2017; Telling et al., 2017), rock mechanics and stability analysis (Assali et al., 2014; Ferrero et al., 2009; Nguyen et al., 2011) as well as landslide monitoring (Lucieer et al., 2014). When compared with traditional survey, they have the following advantages (Buyer and Schubert, 2017; Haneberg, 2008):

6

Journal Pre-proof

1. High-precision, real-time, three-dimensional point cloud data acquisition; 2. Fast and efficient measuring feasibility for large-scale areas or time restricted areas; 3. Contactless and secure investigation for inaccessible and hazardous areas;

f

4. Reproducible and objective results acquisition;

oo

5. Variable degree of acquired data source from macro-scale (e.g. regional

pr

geological hazards survey) to laboratory-scale (e.g. joint roughness

These

techniques

are

now

e-

coefficient estimation using 3D laser device).

well-established

for

the

quantitative

Pr

characterization and mapping of rock mass discontinuities (Abellán et al., 2014;

al

Jaboyedoff et al., 2012). A literature review on the developments in recent years

rn

for the determination of discontinuity orientation, persistence and spacing, is presented in the next paragraphs. Additionally, several semi-automatic and

Jo u

automatic processes are briefly presented. 1.1. Review on discontinuity characterization using new techniques (a) Discontinuity orientation In order to define the geometry of the main discontinuities, sub-regions belonging to discontinuities in the point cloud are separately picked manually and the orientation of each discontinuity is calculated based o n normal vector of the fitting plane (Deliormanli et al., 2014; Sturzenegger and Stead, 2009). However, this method requires a tedious and time-consuming selection of sub-region points. Utilization of a series of triangulated irregular network (TIN)

7

Journal Pre-proof

surfaces generated by interpolated meshing has been developed by some researchers (Vöge et al., 2013; Mah et al., 2011; Gigli and Casagli, 2011). Unfortunately, TIN process for surface simplification will reduce the accuracy of normal vector calculation and sequentially affect orientation results. Other approaches have also been proposed, such as the region-growing approach (Wang et al., 2017), the volumetric pixels (voxels) method and the

oo

f

principal component analysis (Guo et al., 2017; Gomes et al., 2016; Riquelme et al., 2014). However, the basic computing cell (e.g. the dimension of voxels)

pr

which controls the calculation precision has a close correlation with the waviness

e-

and roughness of discontinuity surface (Chen et al., 2016). Hence the

Pr

implementation of new algorithms which accelerate the computation speed and maintain the precision is needed (Abellán et al., 2014).

al

(b) Discontinuity persistence and spacing

rn

In order to obtain discontinuity persistence, the discontinuity shape is

Jo u

assumed to be a circle and the diameter is equated with an “equivalent trace length” (Sturzenegger and Stead, 2009). Umili et al.(2013) detected the vertices curvature values of the surface break lines directly as traces. And Li et al. (2016) utilized the normal tensor voting theory to detect and extract traces. The circular windows method is adopted to statistically compute the mean trace length and estimate the trace intensity (Li et al., 2016; Sturzenegger et al., 2011). Nevertheless, those extraction methods may cause pointless traces existence and true discontinuity traces absence and segmentation.

8

Journal Pre-proof

There are three different types of discontinuity spacing, which can be summarized as (Priest, 1993): total spacing, set spacing, and normal set spacing. Discontinuity clusters and their spatial distributions are studied to calculate total spacing by adopting virtual scanline method (Laux and Henk, 2015; Slob, 2010). A reasonable assumption is made that each discontinuity cluster plane from the same set is parallel for normal set spacing calculation

oo

f

(Riquelme et al., 2015). 1.2. Semi-automatic and automatic processing

pr

Nowadays several semi-automatic or automatic programs have been

e-

developed to characterize rock mass structure. Such examples are ShapeMetrix

Pr

3D (3GSM GmbH), Sirovision (Datamine and CSIRO), Split-FX (Slob, 2010), Coltop-3D (Jaboyedoff et al., 2009), Plane Detect (Vöge et al., 2013),

al

Discontinuity Set Extractor (DSE, Riquelme et al., 2014) and RockScan (Ferrero

rn

et al., 2009). Some software mentioned earlier such as Coltop-3D, have the

Jo u

ability to conduct regional discontinuity orientation trend analysis, but have the limitation to delineate discrete fracture networks (Lato and Vöge, 2012). In the present paper, a

new method for automatic analysis and

characterization of a rock mass based on 3D point clouds is proposed. Several discontinuity parameters, i.e. 1) discontinuity orientation, 2) number of sets, 3) normal set spacing and trace length, and 4) the exact location of every discontinuity at the outcrop can be automatically obtained. 3D mapping of discontinuities and discrete fracture network (DFN) models can also be generated for subsequent discrete fracture analysis.

9

Journal Pre-proof

This study is structured as follows: an introduction of remote sensing techniques and their applications is presented in Chapter 1; the fundamental methodology is presented in Chapter 2; Chapter 3 presents two case studies where the proposed method is applied; application characteristics of the method are discussed in Chapter 4; and conclusions are summarized in Chapter 5.

f

2. Methodology

oo

In this study, several machine learning algorithms are developed to calculate

pr

automatically the discontinuity network of the rock face using 3D point clouds. Machine learning, which is seen as a subset of artificial intelligence, is the

e-

scientific study of algorithms and statistical models for data analysis. The

Pr

computer systems can learn from data, build analytical models, identify patterns and make decisions with minimal manual intervention.

rn

al

The proposed methodology is developed in four main steps (Fig. 1): Step 1: Calculation of normal vectors in order to acquire the major directions

Jo u

from the outcrop point cloud applying the adaptive algorithm (Wang et al., 2013). Step 2: Determination of principal discontinuity sets and mean set orientations. Statistical analysis of the poles (i.e. normal vectors of the discontinuities points) in stereographic projection is performed. Step 3: Cluster segmentation analysis. Every single discontinuity is isolated from its set family and the best fitting plane is calculated. Step 4: Calculation of geometric discontinuity parameters. The location, orientation, trace length of every discontinuity, and normal set spacing are acquired. The 3D map of discontinuities and DFN model can be established.

10

Journal Pre-proof

A typical section of a rock slope (Fig. 2a) as well as its point cloud data (Fig. 2c) obtained by terrestrial laser scanning are selected as the experimental dataset. In addition, three discontinuities which belong to different sets are labeled (Fig. 2b and d) for methodology explanation in this chapter. 2.1. Normal vector calculation There are several methods to estimate the normal vectors from a point cloud,

oo

f

some of which are the following:

pr

• based on local surface fitting (Hoppe et al., 1992),

e-

• based on Delaunay/Voronoi triangulation (Amenta et al., 1999),

Pr

• based on robust statistics (Fleishman et al., 2005). Those conventional methods commonly work accurately when the potential

al

surface is smooth, but show poor and puzzling results significantly near noise

rn

and sharp features (Huang et al., 2013). In order to define sharp features of a

Jo u

point cloud, an assumption is normally made that the underlying surface of a point cloud is piecewise smooth (Fleishman et al., 2005). A sharp feature can be regarded as a corner, an edge or an intersection where different surfaces abut (Wang et al., 2013; Lai et al., 2007). Generally, the rock outcrop surfaces are rough and curved, as well as they contain a great quantity of large-scale undulations and small-scale roughness (Li et al., 2016). Thus, the rock outcrop point cloud dataset is a complex model with high precision which typically contains numerous sharp features. In the present paper, the iterative reweighted plane fitting (IRPF) method proposed by Wang et al. (2013) based on initial local surface fitting studies is applied. IRPF method

11

Journal Pre-proof

provides a reliable and robust estimation in the case of large irregular and un-oriented point clouds containing noise, intersections, outliers and other sharp features. The point cloud 𝑃 (coordinates: 𝑥, 𝑦 and 𝑧, and number of points: 𝑁 ), taken

as

input in

the

proposed

method,

is

a

set

of 3D

points

𝑃 = {𝑝1 , 𝑝2 , … , 𝑝𝑖 , … , 𝑝𝑁 , }. 𝑁𝑘 (𝑝𝑖 ) is used to denote 𝑘 nearest neighbors of a

oo

f

point 𝑝𝑖, which collects 𝑘 neighbors nearest to 𝑝𝑖 regardless of the distance between the point and its neighbors.

pr

The iterative reweighted fitting plane (𝑃𝑙) can be formulated as:

e-

𝑃𝑙 𝑡 (𝑑 𝑡 , 𝒏𝑡 ) = argmin ∑𝑘𝑗=1 𝑟𝑗𝑡 𝑤𝑑 (𝑝𝑗 )𝑤𝑟 (𝑟𝑗𝑡 −1 )𝑤𝑛 (𝒏𝑗𝑡 −1 ) ∥ 𝒏𝑡 ∥= 1

(1)

Pr

where point 𝑝𝑗 belongs to the neighbor 𝑁𝑘 (𝑝𝑖 ), and 𝑑 is the distance from

al

𝑝𝑖 to the fitted plane 𝑃𝑙. 𝒏 is the fitted plane normal, which also represents the

Jo u

𝑡h iteration.

rn

normal of point 𝑝𝑖. 𝑟𝑗𝑡 = 𝑑 𝑡 + (𝑝𝑗 − 𝑝𝑖 )T 𝒏𝑡 is the point 𝑝𝑗 fitting residual at the

𝑤𝑑 (𝑝𝑗 ) = exp(−∥ 𝑝𝑗 − 𝑝𝑖 ∥/𝜎𝑑2 ), 𝑤𝑟 (𝑟𝑗 ) = exp(−(𝑟𝑗 /𝜎𝑟 )2 ) and 𝑤𝑛 (𝒏𝑗 ) = exp(−∥ 𝒏𝑗 − 𝒏𝑖 ∥/𝜎𝑛2 ) respectively represent the Gaussian weight function related to the distance from point 𝑝𝑖 to 𝑝𝑗 , the fitting residual of point 𝑝𝑗 , and the normal difference between 𝒏𝑗 and 𝒏𝑖 . 𝜎𝑑 , 𝜎𝑟 and 𝜎𝑛 are respectively the distance, fitting residual and normal difference bandwidths, which govern the relative contribution of each neighbor. Using the algorithms proposed by Wang et al. (2013), those bandwidth parameters are iteratively computed and determined automatically and adaptively.

12

Journal Pre-proof

The IRPF method can be interpreted as iteratively reweighted plane fitting at point 𝑝𝑖 from its neighbor incorporating distance, residual and normal difference weights. When those three weights are set at every iteration step, Eq. (1) can convert to a positive semi-definite covariance matrix 𝐶𝑀 for eigenvalue decomposition applying the Lagrange multiplier:

(2)

∑ 𝑤𝑑 (𝑝𝑗 )𝑤𝑟 (𝑟𝑗 )𝑤𝑛(𝒏𝑗 )(𝑝𝑗 −𝑝𝑖 ) ∑ 𝑤𝑑 (𝑝𝑗 )𝑤𝑟(𝑟𝑗 )𝑤𝑛 (𝒏𝑗 )

(3)

pr

𝑢=

oo

f

𝑝1 − 𝑝𝑖 − 𝑢 T 𝑝1 − 𝑝𝑖 − 𝑢 𝑝 − 𝑝𝑖 − 𝑢 𝑝 − 𝑝𝑖 − 𝑢 𝐶𝑀 = [ 2 ] [ 2 ] ⋯ ⋯ 𝑝 𝑘 − 𝑝𝑖 − 𝑢 𝑝𝑘 − 𝑝 𝑖 − 𝑢

e-

At the current iteration, the eigenvector corresponding to the smallest eigen

Pr

value of matrix 𝐶𝑀 is the normal vector of point 𝑝𝑖. As shown in Fig. 3b, the iteration terminates when the fitted plane (shown as

al

bold red line) is stabilized. the competent and adaptive neighbor of 𝑝𝑖 (in the

rn

shadowed ellipse) effectively discards the points belonging to the other surface

Jo u

(shown as blue points) due to their insignificant normal difference weights , as well as discards the spatial outlier points (shown as green points) due to their insignificant fitting residual weights during the iteration of Eq. (1). Finally, the neighbor of 𝑝𝑖 only includes the points belonging to the same surface (shown as red points) and replaces the original given neighbor (shown as the black dotted circle) for calculation. The normal vector calculation of point 𝑝𝑖 (the red arrow) shows a more reasonable and accurate outcome than the initial method. Eventually, the normal vectors are assigned on both sides of the edge correctly (Fig. 3a). As the rock outcrop surfaces are commonly rough and curved, IRPF method which applies automatic evaluation of the local adaptive parameters

13

Journal Pre-proof

without any manual intervention, is an optimal choice for precise and reliable normal vector calculation. Subsequently, the associated orientation (dip direction 𝜃, dip angle 𝛿) can be found according to the following equations (Priest, 1993): 𝑚

𝜃 = tan−1 ( 𝑙 ) + 𝑄

(4)

𝑛

𝛿 = tan−1 (√𝑙2 +𝑚2 )

oo

f

(5)

where 𝑙, 𝑚 and 𝑛 are the three components of a normal vector. 𝑄 is an

pr

angular converted value for dip direction calculation and is determined as follows:

Pr

and 𝑚 < 0 or if 𝑙 < 0, and 𝑚 ≥ 0.

e-

𝑄 = 0° , if 𝑙 ≥ 0, and 𝑚 ≥ 0; Q = 360° , if 𝑙 ≥ 0, and 𝑚 < 0; Q = 180° , if 𝑙 < 0,

In conclusion, normal vector calculation and assignment of points to planes

al

are obtained. The next step is the determination of the mean discontinuity set

rn

orientations and removal of non-discontinuity points.

Jo u

2.2. Determination of the main discontinuity sets The subsequent methodology is based on the similarity of the orientations of discontinuities at the same set. If the normal vector 𝒏𝑖 of point 𝑝𝑖 is parallel or approximately parallel to the normal vector 𝒏𝑗 of point 𝑝𝑗 , point 𝑝𝑖 and 𝑝𝑗 are expected to belong to the same discontinuity set. The statistical analysis for the determination of the main discontinuity sets is performed by the pole density contour map in the stereographic projection. Multiple data plotted in an equal area stereonet can be contoured with respect to density, which is useful when evaluating concentrations of poles for the main discontinuity sets.

14

Journal Pre-proof

2.2.1. Identification of the mean orientation of a discontinuity set In literature, the clustering algorithms for the identification of the main discontinuity sets, include: (a) the kernel density estimation (Riquelme et al., 2014), (b) the k-means algorithm (Chen et al., 2016) and (c) the fuzzy c-means algorithm (Guo et al., 2017). In the present methodology, clustering by fast search and find of density peaks (CFSFDP) algorithm is applied (Rodriguez and

oo

f

Laio, 2014). This algorithm is based on the local density and controlled distance of each discontinuity through measuring the similarity of all discontinuity data.

pr

The number of potential discontinuity sets and mean orientation values can be

e-

obtained by the decision graph approach (Rodriguez and Laio, 2014; Gao et al.

Pr

2019).

According to previous studies (Kulatilake et al., 1990; Priest, 1993), the

al

orientations belonging to one discontinuity set basically conform to a Fisher

rn

distribution. The density of the orientations near the mean orientation is very

Jo u

large, and decreases if the distance from the mean orientation increases. This feature coincides well with CFSFDP method and CFSFDP method has been verified more reliable and faster than traditional methods for clustering and identifying discontinuity sets by Gao et al. (2019). The local density 𝑙𝑑𝑖 of the pole data 𝑖 in a stereographic projection is calculated by the following equations: 𝑙𝑑𝑖 = ∑𝑗 𝜒(𝑑𝑚 𝑖𝑗 − 𝑑𝑐𝑓 )

(6)

1, 𝑑𝑚 < 0 𝜒(𝑑 ) = { 0, 𝑑𝑚 ≥ 0

(7)

15

Journal Pre-proof

where 𝑑𝑚 𝑖𝑗 is the distance between pole data points 𝑖 and 𝑗, and 𝑑𝑐𝑓 is a cutoff distance. The local density 𝑙𝑑𝑖 represents the number of discontinuities that are closer than 𝑑𝑐𝑓 to the discontinuity 𝑖. When a discontinuity pole is farther than 𝑑𝑐𝑓 away from all other poles, its local density is defined as zero. The determination of a measured distance (or similarity) between discontinuity orientations is a core issue for subdividing data into sets successfully

oo

f

(Jimenez-Rodriguez and Sitar, 2006). In the original CFSFDP algorithm, the distance measured between every data is set equal to the Euclidean distance.

sets

identification.

Based

e-

discontinuity

pr

However, this setting cannot meet exactly the clustering requirement for on

former

researches

Pr

(Jimenez-Rodriguez and Sitar, 2006; Jimenez, 2008; Xu et al., 2013), the sine or sine-squared of the acute angle between discontinuity normal vectors is used as

al

a measure of their similarity. In this paper, the sine-squared similarity measure is

rn

adopted.

Jo u

Parameter 𝑚𝑑𝑖 is computed as the minimum distance between data point 𝑖 and another higher density pole point: 𝑚𝑑𝑖 = 𝑚𝑖𝑛𝑗:𝜌𝑗> 𝜌𝑖 (𝑑𝑖𝑗 )

(8)

According to 𝑙𝑑 and 𝑚𝑑 values of every discontinuity, CFSFDP algorithm is able to automatically detect the number of clusters. The cluster centers which are measured as local density peaks represent the main orientations in the 3D point cloud. The determination of cutoff distance 𝑑𝑐𝑓 value is discussed in Chapter 4.3.

16

Journal Pre-proof

As shown in Fig. 4a and b, four discontinuity sets are obtained (𝑑𝑐𝑓 = 0.02) for the chosen rock slope area. Comparison of mean orientations with field survey is listed in Table 1. For dip direction values, the maximum difference is 5.5° and the mean difference is 3.6°, and for dip angle values, maximum difference is 3.2° and the mean difference is 2.0°. This result shows that this algorithm is reliable and of high accuracy.

oo

f

2.2.2. Elimination of noise points

Since non-discontinuity points exist in the point cloud, the proposed method

pr

detects if each point belongs to a given discontinuity set or not. If a point is not

e-

assigned to a main set, it is considered as a noise point and Fisher’s 𝐾 value

Pr

iterative calculation is used to eliminate it (Kulatilake, 1985; Priest, 1993; slob, 2010).

al

The variable 𝛼 that obeys the Fisher distribution has the following density: 𝐾 𝑠𝑖𝑛 𝛼𝑒𝐾 𝑐𝑜𝑠 𝛼 𝑒𝐾 −𝑒 −𝐾

(9)

Jo u

rn

𝑓( 𝛼 ) =

where 𝛼 (0 < 𝛼 <

𝜋 2

) is the angular deviation from the mean orientation, and

𝐾 is the Fisher Constant, which is a measure of dispersion of the cluster. When 𝐾 > 3, 𝐾 can be estimated by the following probability density distribution (Kulatilake, 1985; Priest, 1993; Baghbanan and Jing, 2007): 𝑀 −1

𝐾 = 𝑀 −| 𝑟

𝑀|

(10)

where 𝑀 is the total number of discontinuities for one set, and 𝑟𝑀 is the mean orientation. During Fisher’s 𝐾 value iterative calculation, boundary poles (low density points), which are far away from the mean orientation poles (high

17

Journal Pre-proof

density peaks), are removed iteratively until the boundary poles are cleared by plotting on the stereographic projection (Slob et al., 2007). It proved to be hard to define a fixed or standard number of iterations to control this removal process. The number of iterations of every set is varied which depends strongly on the number and the scatter of poles in each set (Slob, 2010). At this point, the corresponding points assigned to discontinuity sets from the

oo

f

3D point cloud are obtained. Four discontinuity sets represented by different colors of the selected rock slope area are presented in Fig.4c. The points that do

pr

not belong to discontinuity sets are discarded.

e-

2.3. Cluster segmentation analysis

Pr

2.3.1. Recognition of single discontinuity cluster

In order to automatically extract points belonging to a single discontinuity

al

plane from one discontinuity set, the algorithm called “Density-Based Spatial

rn

Clustering of Applications with Noise” (DBSCAN) was applied in previous

Jo u

studies (Riquelme et al., 2014; Tonini and Abellan, 2014). However, a disadvantage of DBSCAN algorithm is that it has some problems when the dataset density varies significantly (Ester et al., 1996; Kriegel et al., 2011; Riquelme et al., 2014). The point clouds acquired by LiDAR always have significantly varying densities after pre-processing (e.g. registration and stitching from multi-station scanning data, noise reduction and filtering). Hence in the present methodology, a density-ratio based method proposed by Zhu et al. (2016) is applied, which has the capacity to identify all the clusters from a point cloud dataset of largely varying densities.

18

Journal Pre-proof

The parameters required at this algorithm are: (a) 𝑒𝑝𝑠, which is the distance threshold to determine two points as neighbors (radius of neighborhood); (b) 𝜂, which is defined as radius of the larger neighborhood; (c) 𝜏, which represents the threshold on density ratio. In order to avoid identifying and outputting excessive small clusters, a threshold value 𝑚𝑖𝑛𝑠𝑖𝑧𝑒 is used for discarding small discontinuity clusters. The structure of this algorithm is shown in Appendix A.

oo

f

The main steps of the process and more details can be found in former studies (Ester et al., 1996; Zhu et al., 2016).

pr

Ester et al. (1996) and Riquelme et al. (2014) recommended that the 𝑒𝑝𝑠

e-

value is the 4th neighbor distances mean value plus two standard deviations. Zhu

Pr

et al. (2016) stated that 𝜂 only needs to be slightly larger than 𝑒𝑝𝑠 value (generally 𝜂 is set to 𝜂 = 1.1 × 𝑒𝑝𝑠), and the best value range of 𝜏 is usually

al

between 0.7 and 0.9 based on datasets from real cases. In this paper, the

rn

density ratio threshold 𝜏 is set to 0.7, and 𝜂 is set to 𝜂 = 1.1 × 𝑒𝑝𝑠.

Jo u

Three discontinuities (shown and labeled in Fig.2b), which are individually subordinate to set 1, set 3, and set 4 are selected to verify the quality and effectiveness of this algorithm. As shown in Fig. 4d, the points belonging to those three discontinuities are assigned successfully (threshold 𝑚𝑖𝑛𝑠𝑖𝑧𝑒 is set to 50). In conclusion, at this step we applied the density-ratio based method to extract every single discontinuity cluster from its related set. This is an important step to determine the geometrical parameters of each discontinuity, namely location, orientation, trace length and spacing, as presented later.

19

Journal Pre-proof

2.3.2. Discontinuity plane fitting Considering that a rock discontinuity surface is usually rough and undulating, the orientation of a discontinuity is defined as the orientation of its best fitting plane. Since the spatial points cluster of each discontinuity surface have been identified in the last step, the random sample consensus (RANSAC) algorithm (Fischler and Bolles, 1981) which is used for plane fitting, takes the volatility of

oo

f

the points into consideration and pro vides a more realistic estimation. RANSAC

(Ferrero et al., 2009; Chen et al., 2016).

pr

algorithm has been successfully implemented in previous studies for plane fitting

e-

After utilizing the RANSAC algorithm, the orientations of all discontinuity

Pr

planes can be determined according to Eqs. (11) and (12), where 𝑢𝑎, 𝑢𝑏, and 𝑢𝑐 are the three components of the unit normal vector to the fitted plane.

Jo u

rn

al

90° − 𝑡𝑎𝑛−1 (𝑢𝑏⁄𝑢𝑎) 0° 𝐷𝑖𝑝𝑑𝑖𝑟 𝜃 = 180° −1 𝑢𝑏 {270° − 𝑡𝑎𝑛 ( ⁄𝑢𝑎 ) 0° 𝐷𝑖𝑝 𝛿 = { 90° − 𝑡𝑎𝑛−1 (

|𝑢𝑐| √𝑢𝑎2 +𝑢𝑏 2

)

𝑢𝑎 > 0 𝑢𝑎 = 0＆𝑢𝑏 ≥ 0 𝑢𝑎 = 0＆𝑢𝑏 < 0 𝑢𝑎 < 0

(11)

𝑢𝑎2 + 𝑢𝑏2 = 0 𝑢𝑎2 + 𝑢𝑏2 ≠ 0

(12)

2.4. Geometric parameters calculation Additional geometric parameters, such as location, trace length and spacing are needed when developing a three-dimensional DFN model and performing numerical analysis. 2.4.1. Location

20

Journal Pre-proof

Discontinuity location (𝐷𝐿) is defined as the coordinates of the centroid of all spatial points in one discontinuity cluster, according to Eq. (13) (Sturzenegger and Stead, 2009): 𝐷𝐿 =

∑𝑠𝑧 𝑖=1 𝑝𝑖 𝑠𝑧

𝐷𝐿 ∈ ℝ3 , 𝑠𝑧 ∈ 𝑁

(13)

where 𝑠𝑧 is the total number of the discontinuity cluster points.

oo

f

2.4.2. Trace Length Persistence implies the 3D size of a discontinuity and can be crudely

pr

represented in 2D by the discontinuity trace lengths on a surface exposure. As

e-

shown in Fig. 5a, Laux and Henk (2015) and Riquelme et al. (2018) considered

Pr

the exposed discontinuity surface as a polygon. They calculated all distances between the vertices and assigned the maximum length as the representative

al

trace length of the discontinuity plane. In the present study, the trace length (𝑇𝐿)

rn

is calculated according to this method by Eq. (14): 𝑇𝐿 = 𝑚𝑎𝑥: 𝑒𝑑(𝑝𝑖 , 𝑝𝑗 ) = ||𝑝𝑖 − 𝑝𝑗 || ∀ 𝑖, 𝑗 ∈ {1, … , 𝑠𝑧} 𝑇𝐿, 𝑒𝑑 ∈ ℝ, 𝑠𝑧 ∈ 𝑁

(14)

Jo u

2

where 𝑒𝑑(𝑝𝑖 , 𝑝𝑗 ) is the Euclidean distance of two points 𝑝𝑖 and 𝑝𝑗 , and 𝑠𝑧 is the total number of the discontinuity cluster points. 2.4.3. Spacing 3D spatial discontinuity planes offer the possibility of an accurate and objective measurement of normal set spacing values within each set in three dimensions. In general, the orientation of each discontinuity from the same set is similar but not identical, so an assumption is made that every discontinuity orientation from the same set is assigned to the same mean orientation. As

21

Journal Pre-proof

illustrated in Fig. 5b, the normal set spacing value is defined as the normal distance between two adjacent discontinuity surfaces in one set (Wichmann et al., 2018; Riquelme et al., 2015). When the discontinuity planes are ideally parallel to each other, the spacing can be calculated by Eq.(15): 𝐷𝑖𝑠𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑡𝑦1 : 𝐴𝑥 + 𝐵𝑦 + 𝐶𝑧 + 𝐷1 = 0

f

𝐷𝑖𝑠𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑡𝑦2 : 𝐴𝑥 + 𝐵𝑦 + 𝐶𝑧 + 𝐷2 = 0

oo

|𝐷 −𝐷 |

1 2 𝑆𝑃12 = √𝐴2 +𝐵 2 +𝐶 2

(15)

pr

where 𝐷𝑖𝑠𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑡𝑦1 and 𝐷𝑖𝑠𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑡𝑦2 are two adjacent discontinuity plane

e-

polynomial equations defined by the normal unit vector ( 𝐴, 𝐵, 𝐶) and the

Pr

constant parameter 𝐷1 and 𝐷2 . 𝑆𝑃12 is the spacing value of these two discontinuities.

al

Several parameters (e.g. orientation, and trace length) and fitting plane

rn

equations of the three selected joints are presented in Fig.6. The comparison of

Jo u

orientation and trace length results obtained automatically with those from manual measurements is listed in Table 2. The error values of joint 1 and joint 2 are small, and the error values of joint 3 are slightly high.

3. Application of the proposed method In order to demonstrate the proposed methodology, two rock slope case studies were analyzed. The dataset of case study 1 is used to clarify the workflow, show the results of the proposed method and make a comparison with the methods from previous studies of Chen et al. (2016) and Riquelme et al.

22

Journal Pre-proof

(2014). Case study 2 is used to show the applicability of the proposed method when compared with manual measurements in the field. By using a common laptop (2.30_GHz(R) Intel Core i5-6300Q, 4 GB RAM memory), the total calculation time for case study 1 was 5.5 hours (the number of points is approximately 1,500,000) and for case study 2, 1.5 hours (the number of points is approximately 500,000), respectively. It was found that

oo

f

computational efficiency of this proposed method is acceptable for practical engineering applications.

pr

3.1. Case study 1

e-

The proposed method is applied to a real case study from the Rockbench

Pr

repository (Lato et al., 2013), and the dataset is from a quartzitic rock mass exposed at a road cut slope located in Ouray, Colorado, USA (Fig. 7a). The

al

database is based on a LiDAR point cloud that consists of 4 scan positions

rn

collected by an Optech Ilris 3D Laser Scanner in 2004. The number of data

Jo u

points is 1,515,722, and the point spacing is approximately 2 cm. Using the proposed method, the orientation, number of sets, trace length and spacing values of the identified discontinuities are determined. Orientation data obtained by this method is compared with the results of Chen et al. (2016) and Riquelme et al. (2014). Based on automatic analysis ( 𝑑𝑐𝑓 = 0.007 ), five discontinuity sets are identified (Fig. 8). A summary of the mean orientations, proportions, and number of discontinuities is presented in Table 3. Based on the automatic identification analysis, every single discontinuity surface is acquired. The main discontinuity

23

Journal Pre-proof

sets on the point cloud denoted with different colors are presented in Fig. 12a. Once the RANSAC algorithm is applied, the location, orie ntation, and trace length information of each discontinuity can be determined. Four calculation examples are shown in Fig. 9 and are summarized in Table 4. The 3D trace map of every set from the slope using different colors is presented in Fig. 10a. Additionally, the histograms of trace lengths and the fitted log -normal distribution

oo

f

for each set are presented in Fig. 10b. When calculating the spacing of discontinuities, the discontinuity persistence

pr

(full-persistent or non-persistent) should be considered. Existing methods

e-

presume either full-persistence (e.g., Slob, 2010) or non-persistence (e.g.,

Pr

Riquelme et al., 2015), which add an additional complexity to the discontinuity spacing calculation. In this case, the normal set spacing values of the sets are

al

calculated considering fully persistent discontinuities. For example, the

rn

schematic diagram of the spacing calculation of set 𝐽2 , and set 𝐽4 is shown in Fig. 11a. The histograms of spacing and fitted negative exponential distribution

Jo u

for each set are presented in Fig. 11b. Riquelme et al. (2014) and Chen et al. (2016) analyzed the discontinuity orientations within the same rock mass. The mean orientations of the five sets obtained by Riquelme et al. (2014) and by our proposed method are summarized in Table 3. It is evident that the results agree very well as the maximum deviation of dip direction is 2.6° and the maximum deviation of the dip angle is 4.5°. The orientations of several discontinuities were also calculated by Riquelme et al. (2014) and Chen et al. (2016). These discontinuities are labeled using

24

Journal Pre-proof

different colors in Fig. 12. A comparison of orientation measurements based on the classical approach (Riquelme et al., 2014), our proposed method, and the other two identification methods introduced by Riquelme et al. (2014) and Chen et al. (2016) is presented in Table 5. It is found that the results are similar, and the absolute values of deviation of dip direction and dip angle between the proposed method and the classical approach is approximately between 1° and

oo

f

5°, except for a few discontinuities, such as No. 13, No. 31 and No. 32. The larger deviation of these discontinuities from our results could be attributed to the

pr

surface roughness and waviness or the poor quality and small size of assigned

e-

points.

Pr

3.2. Case study 2

In this case study, a rock slope at Laoyingshan quarry in the City of Shaoxing,

al

China is analyzed (Fig. 7b). A representative and accessible rock mass outcrop

rn

of this slope is chosen for field investigation including both window sampling

Jo u

measurements and terrestrial LiDAR technique. The approximate size of the selected area is 12 m × 2.2 m. This database consists of 3 scan positions collected by a Topcon GLS-20003D Laser Scanner. A total of 504,831 points of the selected area were scanned from a distance of approximately 15 m. Using the proposed method, the orientation, sets and trace length are obtained and the results are compared with the manual measurements. During the automatic analysis ( 𝑑𝑐𝑓 = 0.01), four main sets are detected (Fig. 13a). A summary of the mean orientations, and number of discontinuities is presented in Table 6.

25

Journal Pre-proof

In order to verify the correctness and effectiveness of the proposed method, a careful manual investigation by compass and tape at the selected site is performed. Fig. 13b shows 72 poles obtained manually at the site, and the four main sets, which are identified, are presented in Table 6. The maximum deviation of dip direction is 5.8° and the maximum deviation of the dip angle is 4.7°. The average differences of all trace lengths acquired from the manual

oo

f

survey and the proposed method are presented in Table 6. The average differences of every set are in the range of 0.15 m to 0.23 m, which are

pr

acceptable.

e-

The precise fracture trace mapping by the proposed method and the manual

Pr

survey are shown in Fig.14. It can be seen that the results are very similar with some minor differences as follows:

al

(a) The total number of discontinuities by the proposed method is 80, and the

rn

number by the manual survey is 72. More small discontinuities are extracted by

survey.

Jo u

using the automatic method as these may often be ignored during manual

(b) As shown in the dotted area (1) of Fig.14, due to the presence of vegetation, there are no discontinuities identified by the proposed method. By the influence of the environment, the view of the laser scanner was blocked, and that caused the great difficulty in data collection at this area. (c) As shown in the dotted area (2) of Fig.14, due to the inaccessibility in this area, there are only a few discontinuities measured at this area by the manual survey.

26

Journal Pre-proof

4. Discussion 4.1. Point cloud generation and required quality for the analysis Terrestrial laser scanning and digital photogrammetry with SfM algorithm are two popular remote sensing technologies which are very suitable for close-range dense point cloud data acquisition. The quality of the 3D point clouds generated

f

from on-site surveys has a major influence for accuracy of the final results. The

oo

main conclusions are as follows:

pr

a) Reference system. Registration of the 3D point cloud model within a

e-

precise geodetic coordinate system is essential for the orientation calculation. It is preferred to use geodetic technologies, such as total

Pr

station or RTK GPS device as auxiliary methods, to establish non-collinear ground control points (GCPs) during field survey (Armesto et al., 2009).

al

Based on the known coordinates of GCPs, the coordinates of the 3D point

rn

clouds can be transferred from the local reference system to the global

Jo u

reference system.

b) Ambient occlusions. Some environmental factors (e.g. moving objects, vegetation and debris, obstructions, reflections from wet surfaces, humidity in the air) affect the quality of 3D point clouds during the data collection process (Olsen et al., 2015). For terrestrial LiDAR equipment, selecting a suitable position and performing a multi-site scan can be an appropriate way to avoid these negative effects in the point cloud. Filtering technologies are also recommended by different authors to delete extraneous points, remove vegetation, and reduce instrumental errors

27

Journal Pre-proof

(Brodu and Lague, 2012; Abellán et al., 2014). c) Resolution. One of the most important parameters that represent the quality of the 3D point cloud is the spatial resolution. For terrestrial laser scanning, point cloud resolution is determined by sampling interval variability, beam divergence, and angle of observation (Lichti and Jamtsho, 2006). For digital photogrammetry, the spatial density is chosen

oo

f

depending on the pixel size when doing the SfM calculations (Assali et al., 2014). Abellán et al. (2011) selected the optimal point spacing that ranges

pr

from 4.5 up to 5.5 cm in the study area for rockfall monitoring by terrestrial

e-

laser scanning. Olariu et al. (2008) and Assali et al. (2014) discussed the

Pr

effect of point cloud spatial density on the extraction and interpretation of discontinuity data. Through our attempts of the case studies, the

al

orientation results that plot in stereographic projection for resolutions of 1

rn

cm, 2 cm, 3 cm, and 4 cm were quite similar. The resolution change will cause centimeter-level errors for trace length and spacing measurements

Jo u

which can be acceptable for practical engineering application. As most mainstream instrumentations can reach this resolution under normal use, the authors recommend the point cloud spatial resolution can be in the range from 1 cm to 4 cm. 4.2. Advantages and limitations of the proposed method Based on the analyses and comparisons of the case studies, the proposed method showed a good accuracy for automatic identification of discontinuities in rock mass. In relation to other current discontinuity extraction methods which are based on 2.5D interpolated approach (Slob and Hack, 2004) or TIN to simplify

28

Journal Pre-proof

the surfaces (Gigli and Casagli, 2011), a main advantage of the proposed method herein is that the real point cloud information attributed to each point can be used directly. This approach improves greatly the accuracy of the normal vector calculation. For identification and clustering of rock discontinuity sets, some general algorithms, e.g. k-means algorithm, have the disadvantage that an initial guess

oo

f

of the number of clustering groups a nd clustering centers must be made, which introduces human bias. Based on the CFSFDP algorithm, the number of clusters

pr

and clustering centers doesn’t need to be specified in advance.

e-

It should be noted that this method aims to determine the exposed plane-type

Pr

discontinuity at outcrop and is not capable of determining the line-type discontinuity (i.e. traces). Additionally, other aspects such as scale effect should

al

be considered when using these techniques for rock mass characterization.

rn

4.3. Relevant parameters needed for the analyses

Jo u

For normal vector calculation (at Step 1), IRPF method is able to automatically determine the parameter values without manual data input. For determination of mean orientations and main discontinuity sets based on the CFSFDP algorithm (at Step 2), the cutoff distance 𝑑𝑐𝑓 needs to be defined by the user. The cutoff distance 𝑑𝑐𝑓 influences the local density of the poles. If 𝑑𝑐𝑓 is larger than an appropriate value, it will decrease the differences of local densities and the cluster centers (local density peaks) may not be identified effectively. If 𝑑𝑐𝑓 is smaller than an appropriate value, it will result in a very large number of cluster centers. A detailed sensitivity analysis of 𝑑𝑐𝑓 has been

29

Journal Pre-proof

carried out by Gao et al. (2019), and he proposes a lower and upper limit for 𝑑𝑐𝑓 of 0.005 and 0.12, respectively. Normally, a value of 𝑑𝑐𝑓 larger than 0.12 in the process of clustering for identifying discontinuity sets is inappropriate, as 𝑑𝑐𝑓 of 0.12 leads to an acute angle between two discontinuity normal vectors of approximately 20°. Consequently, a 𝑑𝑐𝑓 of 0.12 is considered as the upper limit of the appropriate range of 𝑑𝑐𝑓 . A 𝑑𝑐𝑓 of 0.005 is selected as the lower limit, as

oo

f

the acute angle between two discontinuity normal vectors is approximately 4°. It should be noted that the determination of a proper 𝑑𝑐𝑓 value depends on the

pr

dispersion of the data.

e-

In this study, the Fisher’s 𝐾 iterative calculation method was employed for the

Pr

elimination of non-discontinuity points. In principle, this method yields good results, when the variation of data points is small. An alternative way to filter out

al

non-discontinuity points is by defining an angle threshold as the maximum

rn

allowed angle difference between the normal vector of one point and the normal

Jo u

vector of the associated mean orientation (Riquelme et al., 2014). During the discontinuity extraction process (at Step 3), the parameters required for the density-ratio based method have been discussed at chapter 2.3.1. The threshold 𝑚𝑖𝑛𝑠𝑖𝑧𝑒 which determines the discontinuity number of each set has an influence on cluster segmentation and scale parameters calculation. When the parameter 𝑚𝑖𝑛𝑠𝑖𝑧𝑒 is small, many minor clusters will be generated, and when the 𝑚𝑖𝑛𝑠𝑖𝑧𝑒 is large, several real clusters will be discarded. For this reason, the threshold 𝑚𝑖𝑛𝑠𝑖𝑧𝑒 for each set should be

30

Journal Pre-proof

selected carefully and it is usually based o n the resolution and size of every set (Table 3 and 6). For geometric parameters calculation (at Step 4), there are no additional parameters needed which have to be defined manually. Futhermore, due to the intersection angle between the discontinuity surface and the outcrop face, acquision and anlysis of discontinuity parameters using 3D

oo

f

point clouds can be affected by sample bias similar to traditional methods, e.g. scanline method (Lato et al., 2010). The Terzaghi bias correction (Terzaghi,

pr

1965) is adopted to reduce sample bias in our research.

e-

4.4. Comparison with manual measurements

Pr

Comparing with manual measurements, the proposed method integrated with close-range remote sensing technologies has a great advantage for field data

al

collection and analysis. The manual field survey is often slow, and dangerous.

rn

For high steep slope or open pit mine, it is always difficult or impossible to

Jo u

access to rock faces directly. Additionally, manual measurements are restricted to the slope foot area, as the height which can be reached without rock climbing equipment is about 2 m. Moreover, the properties of discontinuity usually vary a lot throughout the rock mass, even within a small volume. A large amount of field data is required to accurately describe 3D geometry of a fractured rock mass. Therefore, these remote sensing technologies overcome the problems with the manual measurements (Deliormanli et al., 2014). Considering that an exposed discontinuity surface is usually rough and undulating, through the case studies above, the differences of orientations

31

Journal Pre-proof

calculated by automatic method and manual measurements occur mainly due to the roughness and waviness of discontinuity surfaces. For rough exposed surfaces, manual measurements by compass just make use of an accessible part of the joint surface without considering the whole surface variation. In our proposed method, RANSAC fitting plane algorithm which takes this roughness and waviness into consideration can produce more objective estimation. The

oo

f

discontinuity orientation is related on the best fitting plane. In this respect, we think that the proposed method provides more accurate and reliable orientation

pr

results than manual measurements.

e-

4.5. Application and purpose of the method for practical questions

Pr

Probability statistical distributions for various discontinuity parameters can be derived based on the large amount of field data collected by the proposed

al

method. The comprehensive database including orientation, trace length and

rn

spacing distributions can be applied for establishing a DFN model. DFN model is

Jo u

very useful when analyzing the stability of rock mass and when performing numerical modeling using the discrete element method. For example, a DFN model of the outcrop rock mass in case study 1, is generated using thin circular shape (Fig.15).

In addition, the proposed method can be well applied to perform rock mass classification, fluid flow simulation, or solve other practical engineering problems.

5. Conclusions In the present paper, a new approach for automatic identification and extraction of rock mass discontinuities, clustering of discontinuity sets and

32

Journal Pre-proof

characterization of discontinuity orientation, trace length and spacing using 3D point cloud data is presented. The proposed method employs several machine learning algorithms to perform automatic calculation. The main innovative contributions of the proposed work are: (1) precise normal vector calculation using adaptive IRPF algorithm without simplification or TIN surface; (2) automatic determination of main discontinuity sets and mean

oo

f

orientations; (3) cluster segmentation analysis using density-ratio based method which has the capacity to process a dataset of largely varying densities; (4)

pr

plane fitting using RANSAC algorithm for every single discontinuity; (5)

e-

calculation of geometric parameters and automatic generation of trace maps.

Pr

Using this method, four discontinuity parameters, i.e. number of sets, orientation, spacing and trace length (as suggested by ISRM, 1978) can be

al

obtained. Additionally, every discontinuity location, best fitting plane, trace

rn

mapping and DFN model can also be performed. The method is applied to two

Jo u

real cases and the results are compared with previous studies and manual measurements. Based on the comparisons, it is evident that the proposed method yields accurate and reliable results and can satisfy the practical needs of engineering,

especially

in

a

complex

and

inaccessible

environment.

Furthermore, the proposed method can achieve automatic calculation and reduce manual intervention, as only few parameters have to be defined manually. Future research work will be focused on the extraction of the other discontinuity parameters from 3D point clouds, such as roughness, block size, and improvement of the computational efficiency.

33

Journal Pre-proof

Acknowledgements This research was supported by the National Natural Science Foundation of China (Grant No. 41831290) and the Zhejiang Provincial Natural Science Foundation of China (Grant No. LGF18D020002). The authors gratefully acknowledge Chen et al. (2016) and Riquelme et al. (2014) for their previous studies. The authors also thank Rockbench repository for providing data of the

oo

f

study area. And the authors appreciate the editor and reviewers whose detailed comments and suggestions helped improve this manuscript significantly.

e-

pr

References

1. Abellán, A., Oppikofer, T., Jaboyedoff, M., Rosser, N.J., Lim, M., Lato, M.J.,

Pr

2014. Terrestrial laser scanning of rock slope instabilities. Earth Surf. Process. Landforms 39, 80–97. doi:10.1002/esp.3493

al

2. Abellán, A., Vilaplana, J.M., Calvet, J., García -Sellés, D., Asensio, E., 2011.

rn

Rockfall monitoring by Terrestrial Laser Scanning - Case study of the basaltic

Jo u

rock face at Castellfollit de la Roca (Catalonia, Spain). Nat. Hazards Earth Syst. Sci. 11, 829–841. doi:10.5194/nhess-11-829-2011 3. Amenta, N., Bern, M., Kamvysselis, M., 1999. Surface reconstruction by Voronoi Filtering. Discret. Comput. Geom. 22, 481–504. 4. Armesto, J., Ordóñez, C., Alejano, L., Arias, P., 2009. Terrestrial laser scanning used to determine the geometry of a granite boulder for stability analysis

purposes.

Geomorphology

106,

271–277.

doi:10.1016/j.geomorph.2008.11.005 5. Assali, P., Grussenmeyer, P., Villemin, T., Pollet, N., Viguier, F., 2014. Surveying and modeling of rock discontinuities by terrestrial laser scanning

34

Journal Pre-proof

and photogrammetry: Semi-automatic approaches

for linear outcrop

inspection. J. Struct. Geol. 66, 102–114. doi:10.1016/j.jsg.2014.05.014 6. Baghbanan, A., Jing, L., 2007. Hydraulic properties of fractured rock masses with correlated fracture length and aperture. Int. J. Rock Mech. Min. Sci. doi:10.1016/j.ijrmms.2006.11.001 7. Brodu, N., Lague, D., 2012. 3D terrestrial LiDAR data classification of

oo

f

complex natural scenes using a multi-scale dimensionality criterion: Applications in geomorphology. ISPRS J. Photogramm. Remote Sens.

pr

doi:10.1016/j.isprsjprs.2012.01.006

from

3D

Point

Clouds.

Procedia

Eng.

191,

270–278.

Pr

doi:10.1016/j.proeng.2017.05.181

e-

8. Buyer, A., Schubert, W., 2017. Calculation the Spacing of Discontinuities

al

9. Cawood, A.J., Bond, C.E., Howell, J.A., Butler, R.W.H., Totake, Y., 2017.

rn

LiDAR, UAV or compass-clinometer? Accuracy, coverage and the effects on structural models. J. Struct. Geol. 98. doi:10.1016/j.jsg.2017.04.004

Jo u

10. Chen, J., Zhu, H., Li, X., 2016. Automatic extraction of discontinuity orientation from rock mass surface 3D point cloud. Comput. Geosci. 95, 18–31. doi:10.1016/j.cageo.2016.06.015 11. Deliormanli, A.H., Maerz, N.H., Otoo, J., 2014. Using terrestrial 3D laser scanning and optical methods to determine orientations of discontinuities at a granite

quarry.

Int.

J.

Rock

Mech.

Min.

Sci.

66,

41–48.

doi:10.1016/j.ijrmms.2013.12.007 12. Ester, M., Kriegel, H.P., Sander, J., Xu, X., 1996. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. Proc. 2nd Int.

35

Journal Pre-proof

Conf. Knowl. Discov. Data Min. 226–231. 13. Ferrero, A.M., Forlani, G., Roncella, R., Voyat, H.I., 2009. Advanced geostructural survey methods applied to rock mass characterization. Rock Mech. Rock Eng. 42, 631–665. doi:10.1007/s00603-008-0010-4 14. Fischler, M.A., Bolles, R.C., 1981. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated

oo

f

cartography. Commun. ACM 24, 381–395. doi:10.1145/358669.358692 15. Fleishman, S., Cohen-Or, D., Silva, C.T., 2005. Robust moving least-squares with

sharp

features.

ACM

Trans.

Graph.

24,

544.

pr

fitting

e-

doi:10.1145/1073204.1073227

16. Gao, F., Chen, D., Zhou, K., Niu, W., & Liu, H., 2019. A fast clustering method

Pr

for identifying rock discontinuity sets. KSCE Journal of Civil Engineering,

al

23(2), 556-566. doi:10.1007/s12205-018-1244-7

rn

17. Gigli, G., Casagli, N., 2011. Semi-automatic extraction of rock mass structural data from high resolution LIDAR point clouds. Int. J. Rock Mech.

Jo u

Min. Sci. 48, 187–198. doi:10.1016/j.ijrmms.2010.11.009 18. Giordan, D., Hayakawa, Y., Nex, F., Remondino, F., Tarolli, P., 2017. Review article: The use of remotely piloted aircraft systems (RPAS) for natural hazards monitoring and management. Nat. Hazards Earth Syst. Sci. Discuss. doi:10.5194/nhess-2017-339 19. Gomes, R.K., De Oliveira, L.P.L., Gonzaga, L., Tognoli, F.M.W., Veronez, M.R., De Souza, M.K., 2016. An algorithm for automatic detection and orientation estimation of planar structures in LiDAR-scanned outcrops. Comput. Geosci. 90, 170–178. doi:10.1016/j.cageo.2016.02.011

36

Journal Pre-proof

20. Guo, J., Liu, S., Zhang, P., Wu, L., Zhou, W., Yu, Y., 2017. Towards semi-automatic rock mass discontinuity orientation and set analysis from 3D point

clouds.

Comput.

Geosci.

103,

164–172.

doi:10.1016/j.cageo.2017.03.017 21. Haneberg, W.C., 2008. Using close range terrestrial digital photogrammetry for 3-D rock slope modeling and discontinuity mapping in the United States.

oo

f

Bull. Eng. Geol. Environ. 67, 457–469. doi:10.1007/s10064-008-0157-y 22. Hoppe, H., DeRose, T., Duchamp, T., McDonald, J., Stuetzle, W., 1992.

pr

Surface reconstruction from unorganized points. ACM SIGGRAPH Comput.

e-

Graph. 26, 71–78. doi:10.1145/142920.134011

23. Huang, H., Wu, S., Gong, M., Cohen-Or, D., Ascher, U., Zhang, H. R., 2013. point

set

resampling.

Pr

Edge-aware

ACM

Trans.

Graph.

al

doi:10.1145/2421636.2421645

rn

24. ISRM, 1978. Suggested methods for the quantitative description of discontinuities in rock masses. Int. J. Rock Mech. Min. Sci. 15, 319–368.

Jo u

25. Jaboyedoff, M., Couture, R., Locat, P., 2009. Structural analysis of Turtle Mountain (Alberta) using digital elevation model: Toward a progressive failure. Geomorphology. doi:10.1016/j.geomorph.2008.04.012 26. Jaboyedoff, M., Oppikofer, T., Abellán, A., Derron, M.H., Loye, A., Metzger, R., Pedrazzini, A., 2012. Use of LIDAR in landslide investigations: A review. Nat. Hazards 61, 5–28. doi:10.1007/s11069-010-9634-2 27. Jimenez, R., 2008. Fuzzy spectral clustering for identification of rock discontinuity sets. Rock Mech. Rock Eng. doi:10.1007/s00603-007-0155-6 28. Jimenez-Rodriguez, R., Sitar, N., 2006. A spectral method for clustering of

37

Journal Pre-proof

rock

discontinuity

sets.

Int.

J.

Rock

Mech.

Min.

Sci.

doi:10.1016/j.ijrmms.2006.02.003 29. Kriegel, H., Kröger, P., Sander, J., Zimek, A., 2011. Density‐based clustering. WIREs Data Mining & Know. Discovery, 1: 231-240. doi:10.1002/widm.30 30. Kulatilake, P.H.S.W., 1985. Fitting Fisher Distributions to Discontinuity Orientation Data, Journal of Geological Education, 33:5, 266-269, DOI:

oo

f

10.5408/0022-1368-33.5.266

31. Kulatilake, P.H.S.W., Wu, T.H., Wathugala, D.N., 1990. Probabilistic

pr

modelling of joint orientation. Int. J. Numer. Anal. Methods Geomech. 14,

e-

325–350. doi:10.1002/nag.1610140503

Pr

32. Lai YK, Zhou QY, Hu SM, Wallner J, Pottmann H. Robust feature classification and editing. IEEE Transactions on Visualization and Computer

al

Graphics 2007; 13(1):34–45

rn

33. Lato, M., Kemeny, J., Harrap, R.M., Bevan, G., 2013. Rock bench:

Jo u

Establishing a common repository and standards for assessing rockmass characteristics using LiDAR and photogrammetry. Comput. Geosci. 50. doi:10.1016/j.cageo.2012.06.014 34. Lato, M.J., Vöge, M., 2012. Automated mapping of rock discontinuities in 3D lidar and photogrammetry models. Int. J. Rock Mech. Min. Sci. 54, 150–158. doi:10.1016/j.ijrmms.2012.06.003 35. Laux, D., Henk, A., 2015. Terrestrial laser scanning and fracture network characterisation – perspectives for a (semi-) automatic analysis of point cloud data from outcrops. [Terrestrisches Laserscanning zur Charakterisierung von Kluftnetz- werken – Perspektivenfüreine (halb-) automatische Analyse von

38

Journal Pre-proof

Punktwolken aus Oberflächen aufschlüssen.] Z. Dt. Ges. Geowiss. 166: 99–118. doi:10.1127/1860-1804/2015/0089 36. Li, X., Chen, J., Zhu, H., 2016. A new method for automated discontinuity trace mapping on rock mass 3D surface model. Comput. Geosci. 89, 118–131. doi:10.1016/j.cageo.2015.12.010 37. Lichti, D.D., Jamtsho, S., 2006. Angular resolution of terrestrial laser

oo

f

scanners. Photogramm. Rec. doi:10.1111/j.1477-9730.2006.00367.x 38. Lucieer, A., Jong, S.M. d., Turner, D., 2014. Mapping landslide displacements

pr

using Structure from Motion (SfM) and image correlation of multi-temporal

e-

UAV photography. Prog. Phys. Geogr. 38. doi:10.1177/0309133313515293 39. Mah, J., Samson, C., McKinnon, S.D., 2011. 3D laser imaging for joint Int. J.

Rock

Pr

orientation analysis.

Mech. Min. Sci. 48, 932–941.

al

doi:10.1016/j.ijrmms.2011.04.010

rn

40. Nguyen, H.T., Fernandez-Steeger, T.M., Wiatr, T., Rodrigues, D., Azzam, R., 2011. Use of terrestrial laser scanning for engineering geological applications

Jo u

on volcanic rock slopes - An example from Madeira Island (Portugal). Nat. Hazards Earth Syst. Sci. 11, 807–817. doi:10.5194/nhess-11-807-2011 41. Olariu, M.I., Ferguson, J.F., Aiken, C.L.V., Xu, X., 2008. Outcrop fracture characterization using terrestrial laser scanners: Deep-water Jackfork sandstone

at big rock quarry, Arkansas. Geosphere 4, 247–259.

doi:10.1130/GES00139.1 42. Olsen, M.J., Wartman, J., McAlister, M., Mahmoudabadi, H., O’Banion, M.S., Dunham, L., Cunningham, K., 2015. To fill or not to fill: Sensitivity analysis of the influence of resolution and hole filling on point cloud surface modeling

39

Journal Pre-proof

and individual rockfall event detection. Remote Sens. 7, 12103–12134. doi:10.3390/rs70912103 43. Priest, S.D., 1993. Discontinuity analysis for rock engineering, Chapman & Hall. doi:10.1016/0926-9851(93)90044-Y 44. Riquelme, A.J., Abellán, A., Tomás, R., 2015. Discontinuity spacing analysis in rock masses using 3D point clouds. Eng. Geol. 195, 185–195. doi:

oo

f

10.1016/j.enggeo.2015.06.009 45. Riquelme, A.J., Abellán, A., Tomás, R., Jaboyedoff, M., 2014. A new

pr

approach for semi-automatic rock mass joints recognition from 3D point

e-

clouds. Comput. Geosci. 68, 38–52. doi:10.1016/j.cageo.2014.03.014 46. Riquelme, A., Tomás, R., Cano, M., Pastor, J.L., Abellán, A., 2018. Automatic

Rock

Mech.

Rock

Eng.

51,

3005–3028.

al

Clouds.

Pr

Mapping of Discontinuity Persistence on Rock Masses Using 3D Point

rn

doi:10.1007/s00603-018-1519-9

47. Rodriguez, A., o, A., 2014. Clustering by fast search and find of density

Jo u

peaks. Science (80-.). doi:10.1126/science.1242072 48. Slob, S., 2010. Automated rock mass characterisation using 3 -D terrestrial laser scanning. (PhD thesis) Delft University of Technology. 49. Slob, S., Hack, R., 2004. 3D terrestrial laser scanning as a new field measurement and monitoring technique. Lect. Notes Earth Sci. 104, 179–189. 50. Slob, S., van Knapen, B., Hack, R., Turner, K., Kemeny, J., 2007. Method for Automated Discontinuity Analysis of Rock Slopes with Three-Dimensional Laser Scanning. Transp. Res. Rec. J. Transp. Res. Board 1913, 187–194.

40

Journal Pre-proof

doi:10.3141/1913-18 51. Sturzenegger, M., Stead, D., 2009. Quantifying discontinuity orientation and persistence on high mountain rock slopes and large landslides using terrestrial remote sensing techniques. Nat. Hazards Earth Syst. Sci. 9, 267–287. doi:10.5194/nhess-9-267-2009 52. Sturzenegger,

M., Stead,

D.,

Elmo, D.,

2011.

Terrestrial

remote

oo

f

sensing-based estimation of mean trace length, trace intensity and block size/shape. Eng. Geol. 119, 96–111. doi:10.1016/j.enggeo.2011.02.005

using

terrestrial

laser

scanning.

Earth-Science

Rev.

e-

research

pr

53. Telling, J., Lyda, A., Hartzell, P., Glennie, C., 2017. Review of Earth science

doi:10.1016/j.earscirev.2017.04.007

Pr

54. Terzaghi, R.D., 1965. Sources of Error in Joint Surveys. Géotechnique, 15,

al

287-304.

clouds:

A

rn

55. Tonini, M., Abellan, A., 2014. Rockfall detection from terrestrial LiDAR point clustering

approach

using

R.

J.

Spat.

Inf.

Sci.

Jo u

doi:10.5311/JOSIS.2014.8.123 56. Umili, G., Ferrero, A., Einstein, H.H., 2013. A new method for automatic discontinuity traces sampling on rock mass 3D model. Comput. Geosci. 51, 182–192. doi:10.1016/j.cageo.2012.07.026 57. Vöge, M., Lato, M.J., Diederichs, M.S., 2013. Automated rockmass discontinuity mapping from 3-dimensional surface data. Eng. Geol. 164, 155–162. doi:10.1016/j.enggeo.2013.07.008 58. Wang, X., Zou, L., Shen, X., Ren, Y., Qin, Y., 2017. A region-growing approach for automatic outcrop fracture extraction from a three-dimensional

41

Journal Pre-proof

point cloud. Comput. Geosci. 99. doi:10.1016/j.cageo.2016.11.002 59. Wang, Y., Feng, H.Y., Delorme, F.É., Engin, S., 2013. An adaptive normal estimation method for scanned point clouds with sharp features. CAD Comput. Aided Des. 45, 1333–1348. doi:10.1016/j.cad.2013.06.003 60. Wichmann, V., Strauhal, T., Fey, C., Perzlmaier, S., 2018. Derivation of space-resolved normal joint spacing and in situ block size distribution data

oo

f

from terrestrial LIDAR point clouds in a rugged Alpine relief (Kühtai, Austria). Bull. Eng. Geol. Environ. doi:10.1007/s10064-018-1374-7

pr

61. Xu, L.M., Chen, J.P., Wang, Q., Zhou, F.J., 2013. Fuzzy C-Means Cluster

Grouping

of

Discontinuity

D.,

Greenwood,

W.,

Rock

Lynch,

Mech.

J.,

Rock

Manousakis,

Eng.

J.,

al

62. Zekkos,

Sets.

Pr

doi:10.1007/s00603-012-0244-z

e-

Analysis Based on Mutative Scale Chaos Optimization Algorithm for the

rn

Athanasopoulos-Zekkos, A., Clark, M., Cook, K.L., Saroglou, C., 2018. Lessons learned from the application of UAV-enabled Structure-from-Motion

Jo u

photogrammetry in geotechnical engineering. International Journal of Geoengineering Case Histories, Vol.4, Issue 4, pp. 254-274. doi: 10.4417/IJGCH-04-04-03 63. Zhu, Y., Ming, K., Carman, M.J., 2016. Density-ratio based clustering for discovering clusters with varying densities. Pattern Recognit. 60, 983–997. doi:10.1016/j.patcog.2016.07.007

42

Journal Pre-proof

Figure Captions Fig. 1. Work flow chart of the proposed methodology. Fig. 2. A rock slope as well as its point cloud data obtained by terrestrial laser scanning are selected as experimental dataset for methodology explanation. (a) A photo of the rock slope; (b) Three discontinuities which belong to different sets; (c) The point cloud of the selected slope part; (d) Three labeled discontinuities in point cloud.

oo

f

Fig. 3. (a) Normal calculation using the IRPF method at a sharp position from a point cloud (the green dots represent the points and the short blue lines represent normal vectors); (b) Detailed schematic diagram of normal calculation (the red arrow) of one point 𝑝𝑖 combining three weights.

e-

pr

Fig. 4. Some results are obtained using the dataset of the selected slope area (Fig. 2). The number of data points is 88,075. The pole plots (a) and the contour plots (b) in the equal area stereonet (lower hemisphere); (c) Four sets are identified, and every point is assigned to the corresponding set represented by different colors ( 𝐽1 − blue, 𝐽2 − magenta, 𝐽3 − green, 𝐽4 − red ); (d) Automatic identification of three discontinuities (labeled in Fig. 2b).

Pr

Fig. 5. (a) Schematic view of the 3D point cloud for trace length calculation; (b) Schematic diagram of normal set spacing calculation.

al

Fig. 6. Calculation of the parameters: orientation, trace length and fitting plane equation of the three selected joint planes.

Jo u

rn

Fig. 7. (a) View of the road cut slope in Ouray, Colorado, USA (from Rockbench repository); (b) A rock slope of the Laoyingshan quarry in the City of Shaoxing, China. The area selected in the red rectangle is used for analysis. Fig. 8. Clustering of poles using stereographic projection (equal area net, lower hemisphere) of case study 1. Fig. 9. Calculation of location, dip direction, dip angle, trace length and fitting plane equation of each single discontinuity. Four examples of case study 1 are shown: (a) one discontinuity from set 𝐽1 , (b) one discontinuity from set 𝐽2 , (c) one discontinuity from set 𝐽3 , and (d) one discontinuity from set 𝐽5 . Fig. 10. (a) Trace mapping of each set using different colors (𝐽1 − red, 𝐽2 − orange, 𝐽3 − blue, 𝐽4 − yellow, and 𝐽5 − green ); (b) Trace length histogram with fitted log-normal distribution (probability density function, PDF). Fig. 11. Discontinuity normal set spacing calculation considering full persistence. (a) Normal set spacing results of set 𝐽2 and set 𝐽4 ; (b) Spacing histogram with fitted negative exponential distribution.

43

Journal Pre-proof

Fig. 12. (a) Discontinuity sets represented by different colors (𝐽1 − red, 𝐽2 − orange, 𝐽3 − blue, 𝐽4 − yellow, and 𝐽5 − green). (b)- (d) Several labeled discontinuities which are used for calculation. Fig. 13. Stereographic projection (equal area net, lower hemisphere) of case study 2 using automatic method (a) and manual measurements (b). Fig. 14. Trace mapping by the proposed method (a) and by manual investigation (b).

Jo u

rn

al

Pr

e-

pr

oo

f

Fig. 15. DFN model of the slope of case study 1 (𝐽1 − red, 𝐽2 − orange, 𝐽3 − blue, 𝐽4 − yellow, and 𝐽5 − green).

44

Jo u

rn

al

Pr

e-

pr

oo

f

Journal Pre-proof

Fig. 1. Work flow chart of the proposed methodology.

45

Pr

e-

pr

oo

f

Journal Pre-proof

Jo u

rn

al

Fig. 2. A rock slope as well as its point cloud data obtained by terrestrial laser scanning are selected as experimental dataset for methodology explanation. (a) A photo of the rock slope; (b) Three discontinuities which belong to different sets; (c) The point cloud of the selected slope part; (d) Three labeled discontinuities in point cloud.

46

f

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

oo

Fig. 3. (a) Normal calculation using the IRPF method at a sharp position from a point cloud (the green dots represent the points and the short blue lines represent normal vectors); (b) Detailed schematic diagram of normal calculation (the red arrow) of one point 𝑝𝑖 combining three weights.

47

Pr

e-

pr

oo

f

Journal Pre-proof

Jo u

rn

al

Fig. 4. Some results are obtained using the dataset of the selected slope area (Fig. 2). The number of data points is 88,075. The pole plots (a) and the contour plots (b) in the equal area stereonet (lower hemisphere); (c) Four sets are identified, and every point is assigned to the corresponding set represented by different colors ( 𝐽1 − blue, 𝐽2 − magenta, 𝐽3 − green, 𝐽4 − red ); (d) Automatic identification of three discontinuities (labeled in Fig. 2b).

48

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

oo

f

Fig. 5. (a) Schematic view of the 3D point cloud for trace length calculation; (b) Schematic diagram of normal set spacing calculation.

49

Jo u

rn

al

Pr

e-

pr

oo

f

Journal Pre-proof

Fig. 6. Calculation of the parameters: orientation, trace length and fitting plane equation of the three selected joint planes.

50

oo

f

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

Fig. 7. (a) View of the road cut slope in Ouray, Colorado, USA (from Rockbench repository); (b) A rock slope of the Laoyingshan quarry in the City of Shaoxing, China. The area selected in the red rectangle is used for analysis.

51

oo

f

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

Fig. 8. Clustering of poles using stereographic projection (equal area net, lower hemisphere) of case study 1.

52

Pr

e-

pr

oo

f

Journal Pre-proof

Jo u

rn

al

Fig. 9. Calculation of location, dip direction, dip angle, trace length and fitting plane equation of each single discontinuity. Four examples of case study 1 are shown: (a) one discontinuity from set 𝐽1 , (b) one discontinuity from set 𝐽2 , (c) one discontinuity from set 𝐽3 , and (d) one discontinuity from set 𝐽5 .

53

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

oo

f

Fig. 10. (a) Trace mapping of each set using different colors (𝐽1 − red, 𝐽2 − orange, 𝐽3 − blue, 𝐽4 − yellow, and 𝐽5 − green ); (b) Trace length histogram with fitted log-normal distribution (probability density function, PDF).

54

oo

f

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

Fig. 11. Discontinuity normal set spacing calculation considering full persistence. (a) Normal set spacing results for set 𝐽2 and set 𝐽4 ; (b) Spacing histogram with fitted negative exponential distribution for each set.

55

e-

pr

oo

f

Journal Pre-proof

Jo u

rn

al

Pr

Fig. 12. (a) Discontinuity sets represented by different colors (𝐽1 − red, 𝐽2 − orange, 𝐽3 − blue, 𝐽4 − yellow, and 𝐽5 − green). (b)- (d) Several labeled discontinuities which are used for calculation.

56

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

oo

f

Fig. 13. Stereographic projection (equal area net, lower hemisphere) of case study 2 using automatic method (a) and manual measurements (b).

57

oo

f

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

Fig. 14. Trace mapping by the proposed method (a) and by manual investigation (b).

58

f

Journal Pre-proof

Jo u

rn

al

Pr

e-

pr

oo

Fig. 15. DFN model of the slope of case study 1 (𝐽1 − red, 𝐽2 − orange, 𝐽3 − blue, 𝐽4 − yellow, and 𝐽5 − green).

59

Journal Pre-proof

Table 1 Comparison of discontinuity orientations using the proposed method with manual field survey. Proposed method

Field survey

Differences *

Dip dir /Dip (°)

Dip dir /Dip (°)

∆|𝐃𝐃|

∆|𝐃𝐀|

𝑱𝟏

181.4/68.4

186/70

4.6

1.6

𝑱𝟐

316.5/82.2

311/79

5.5

3.2

𝑱𝟑

353.1/55.0

355/53

1.9

2.0

𝑱𝟒

33.6/68.2

36/67

2.4

1.2

oo

f

Set

Jo u

rn

al

Pr

e-

pr

*∆|DD| and ∆|DA| respectively represent the dip direction and dip angle differences.

60

Journal Pre-proof

Table 2 Comparison of orientation and trace length results with manual measurements. Automatic results Joint

Set

Manual survey

Differences *

Dip dir /Dip (°)

𝑻𝑳 (m)

Dip dir /Dip (°)

𝑻𝑳 (m)

∆|𝐃𝐃|

∆|𝐃𝐀|

∆|𝐃𝐋|

Joint 1

No.4

216.5/66.4

2.31

214/68

2.20

2.5

1.6

0.11

Joint 2

No.3

169.2/54.4

1.69

166/50

1.60

3.2

4.4

0.09

Joint 3

No.1

183.6/66.7

3.79

190/61

3.63

6.4

5.7

0.16

oo

f

*∆|DL| (m) represents the trace length difference between automatic recognized

Jo u

rn

al

Pr

e-

pr

results and manual survey.

61

Journal Pre-proof

Table 3 Results of grouping discontinuity sets of case study 1. Se

Number of

Dip dir /Dip (°) Proportion *

t

Dip dir /Dip (°)

𝒎𝒊𝒏𝒔𝒊𝒛𝒆 value

Number of points

discontinuities

by proposed method

𝑱𝟏

246.4/34.7

46.62%

538,720

300

47

𝑱𝟐

172.5/81.3

2.88%

33,246

100

14

𝑱𝟑

135.5/81.9

19.19%

221,722

300

30

𝑱𝟒

90.8/51.1

8.96%

103,560

200

𝑱𝟓

290.2/72.7

13.47%

155,668

l a n 200

2.6

2.0

172.3/83.2

0.2

1.9

137.3/77.9

1.8

4.0

32

93.0/48.7

2.2

2.4

47

288.5/68.2

1.7

4.5

r P

r u o

J

62

∆|𝐃𝐀|

249.0/36.7

o r p

e

*Proportion represents the percentage of set points of all used data points.

f o

∆|𝐃𝐃|

by Riquelme (2014)

Journal Pre-proof Table 4 List of geometric parameters of four example discontinuities (shown in Fig. 9). No.

S et

Location*

Orientation(°)

x

y

z

Fitting plane equation (𝒁 = 𝑨𝑿 + 𝑩𝒀 + 𝑪)

𝑻𝑳 (m)

𝑨

𝑩

𝐂

a

𝑱𝟏

27.6720

6.2147

4.0846

256.3/35.9

10.53

0.7032

0.1709

-16.4366

b

𝑱𝟐

26.0937

-6.7581

6.0475

176.9/74.5

11.06

-0.1959

3.5965

8.3818

c

𝑱𝟑

27.0294

8.0901

-0.1330

308.9/54.2

8.20

1.0810

0.8716

-27.4773

d

𝑱𝟓

22.6070

1.3540

-2.3046

286.5/72.5

4.35

3.0354

-08967

-69.4553

Jo u

rn

al

Pr

e-

pr

oo

f

*The location here means the discontinuity position in the computing coordinate system, and can transfer to the geodetic coordinate system.

63

Journal Pre-proof

Table 5 Comparison of orientation results using different methods. No.

Classical approach

Riquelme (2014)

Chen (2016)

New proposed method

Riquelme (2014) ∆|𝐃𝐃|

11

249.2/40.2

246.2/39.0

244.6/38.4

247.5/39.0

3.0

12

264.2/57.0

256.9/52.3

256.2/52.2

258.4/55.0

7.3

13

264.0/41.9

70.3/35.8

251.0/36.2

74.8/38.5

14

252.6/36.5

252.7/35.5

251.4/33.9

256.3/35.9

15

248.7/37.0

249.7/35.9

250.8/36.8

249.6/36.4

16

254.8/29.8

70.5/35.9

250.5/35.9

251.8/33.3

17

249.9/35.9

255.1/32.7

253.2/33.5

21

338.7/82.4

339.5/83.3

157.6/83.8

22

347.5/79.0

166.3/76.6

23

341.0/89.5

160.2/89.9

24

353.5/76.4

173.6/76.9

31

314.1/77.2

32 33

Chen (2016)

f o

∆|𝐃𝐀|

o r p

∆|𝐃𝐀|

∆|𝐃𝐃|

∆|𝐃𝐀|

4.6

1.8

1.7

1.2

4.7

8.0

4.8

5.8

2.0

6.1

13.0

5.7

9.2

3.4

0.1

1.0

1.2

2.6

3.7

0.6

1.0

1.1

2.1

0.2

0.9

0.6

4.3

6.1

4.3

6.1

3.0

3.5

250.6/33.6

5.2

3.2

3.3

2.4

0.7

2.3

338.0/82.5

0.8

0.9

1.1

1.4

0.7

0.1

166.3/78.7

345.2/77.2

1.2

2.4

1.2

0.3

2.3

1.8

157.5/86.9

342.4/88.4

0.8

0.4

3.5

2.6

1.4

1.1

353.1/77.8

176.9/74.5

0.1

0.5

0.4

1.4

3.4

1.9

136.6/82.6

314.7/80.0

320.5/81.3

2.5

5.4

0.6

2.8

6.4

4.1

302.4/75.9

131.2/82.7

136.5/89.9

307.8/80.2

8.8

6.8

14.1

14.0

5.4

4.3

330.2/83.0

143.9/89.7

145.6/89.9

331.9/79.6

6.3

6.7

4.6

6.9

1.7

3.4

al

e 13.7

r P

n r u

o J

64

1.2

∆|𝐃𝐃|

New proposed method

Journal Pre-proof

41

286.1/58.9

97.6/63.2

286.0/59.8

283.7/60.2

8.5

4.3

0.1

0.9

2.4

1.3

42

274.2/51.1

91.1/50.2

272.6/47.6

268.3/50.7

3.1

0.9

1.6

3.5

5.9

0.4

43

277.2/46.4

96.6/48.0

277.3/49.3

276.0/44.3

0.6

1.6

0.1

2.9

1.2

2.1

51

305.0/77.6

123.4/76.2

305.0/77.6*

305.6/79.1

1.6

1.4

16.3

4.4

0.6

1.5

52

290.2/67.0

105.8/69.9

109.3/76.6

286.5/70.5

4.4

2.9

0.9

9.6

3.7

3.5

16.3

14.0

9.2

4.3

4.3

3.9

3.2

2.1

Maximum deviation

13.7

Average deviation

3.9

* It maybe a data printing error in the original text (Chen et al., 2016).

l a n

e

r P

r u o

J

65

f o

o r p 6.8 3.0

Journal Pre-proof

Table 6 Results of grouping and extraction discontinuity sets from case study 2. Automatic method S et

Dip dir /Dip (°)

𝒎𝒊𝒏𝒔𝒊𝒛𝒆

Number of

Dip dir /Dip (°)

discontinuities

Number of

Differences ∆|𝐃𝐃| ∆|𝐃𝐀| ∆|𝐃𝐋|

discontinuities

𝑱𝟏

28.9/78.8

100

15

30/83

16

1.1

4.2

0.21

𝑱𝟐

61.5/37.7

100

25

58/39

21

3.5

1.3

0.15

𝑱𝟑

270.1/68.3

150

21

269/73

19

1.1

4.7

0.18

𝑱𝟒

341.2/67.4

100

19

347/63

f

value

Manual measurements

5.8

4.4

0.23

Jo u

rn

al

Pr

e-

pr

oo

16

66

Journal Pre-proof

Appendix A. The structure of the density-ratio based method Input: 𝐷𝑠 : 3D points dataset of one discontinuity set; 𝑒𝑝𝑠 : radius of neighborhood; 𝜂: radius of the larger neighborhood; 𝜏: threshold on density ratio. Output: 𝐶𝑙 – Class labels of every points in dataset. Tag all points as uncalled status.

2.

While exist uncalled points in 𝐷𝑠 do

3.

Choose one uncalled points 𝑝 and tag 𝑝 as called status;

4.

𝑁𝐵𝑝𝑠 (𝑝) ← { 𝑞 ∈ 𝐷𝑠 | ∥ 𝑝 − 𝑞 ∥ ≤ 𝑒𝑝𝑠} ;

5.

𝑁𝐵𝑡𝑎 (𝑝) ← { 𝑞 ∈ 𝐷𝑠 | ∥ 𝑝 − 𝑞 ∥ ≤ 𝜂} ;

6.

If

7.

Create a new cluster 𝐶𝑆, and add 𝑝 into 𝐶𝑆;

8.

𝑁𝐵 ← 𝑁𝐵𝑝𝑠 (𝑝);

9.

For each point 𝑝 ′ in 𝑁𝐵 do If 𝑝 ′ is uncalled then

11.

Tag 𝑝 ′ as called status;

12.

𝑁𝐵𝑝𝑠 (𝑝 ′ ) ← { 𝑞 ∈ 𝐷𝑠 | ∥ 𝑝 ′ − 𝑞 ∥≤ 𝑒𝑝𝑠};

13.

𝑁𝐵𝑡𝑎 (𝑝 ′ ) ← { 𝑞 ∈ 𝐷𝑠 | ∥ 𝑝 ′ − 𝑞 ∥≤ 𝜂};

14.

If

15.

𝑁𝐵 ← 𝑁𝐵 ∪ 𝑁𝐵𝑝𝑠 (𝑝 ′ );

16.

End If

17.

End If

18.

If 𝑝 ′ is not yet a member of any cluster then

19.

Add 𝑝 ′ into 𝐶𝑆;

20.

End If

21.

End for

22.

Output 𝐶𝑆;

23.

Else

24.

Tag 𝑝 as NOISE;

25.

End If

26.

End While

al

≥ 𝜏 then

Jo u

rn

|𝑁𝐵𝑝𝑠 (𝑝′ ) |

Pr

10.

|𝑁𝐵𝑝𝑠 (𝑝′ ) |

oo

≥ 𝜏 then

pr

|𝑁𝐵𝑝𝑠 (𝑝) |

e-

|𝑁𝐵𝑝𝑠 (𝑝) |

f

1.

67

Journal Pre-proof

Highlights: 

An automatic method using several machine learning algorithms for identification and extraction of rock mass discontinuities is proposed.



Four discontinuity parameters, namely number of sets, orientation, spacing and trace length can be obtained. Discontinuity location, best fitting plane,

This method is applied for two real cases and produces reliable and

oo



f

and 3D trace mapping can also be performed.

pr

accuracy results.

Pr

e-

Author Statement: Deheng Kong developed the methodology, performed data analysis, and wrote the manuscript. Faquan Wu supervised the research and contributed to preparation of the manuscript. Charalampos Saroglou discussed the algorithms and results and critically reviewed the manuscript.

Jo u

rn

al

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

68