Validation of fracture data recognition in rock masses by automated plane detection in 3D point clouds

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

Contents lists available at ScienceDirect

International Journal of Rock Mechanics and Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms

Validation of fracture data recognition in rock masses by automated plane detection in 3D point clouds

T

⁎

T. Drewsa, , G. Miernika, K. Andersb, B. Höfleb, J. Profeb,c, A. Emmerichd, T. Bechstädta,d a

Institute of Earth Sciences, Heidelberg University, Im Neuenheimer Feld 234, Heidelberg, Germany Institute of Geography, 3DGeo Research Group, Heidelberg University, Im Neuenheimer Feld 368, Heidelberg, Germany c Department of Geography, Justus-Liebig-University, Senckenbergstr. 1, Gießen, Germany d GeoResources Steinbeis Transfer Center, Im Neuenheimer Feld 234, Heidelberg, Germany b

A R T I C LE I N FO

A B S T R A C T

Keywords: LiDAR 3D point cloud Geological outcrop acquisition Plane detection Fracture intensity Spherical-orientation statistics

This paper presents (1) an automated method to extract planes and their spatial orientation directly from 3D point clouds, followed by (2) extensive validation tests accompanied by thorough statistical analysis, and (3) a fracture intensity calculation on automatically segmented planes. For the plane extraction, a region growing segmentation algorithm controlled by several input parameters is applied to a point cloud of a granite outcrop. Within its complex surface shape, more than 1000 compass measurements were conducted for validation. In addition, digitally handpicked planes in the software Virtual Reality Geological Studio (VRGS) were used for single plane comparison. In a second test site, we performed fracture intensity calculation in Petrel based on results of the segmentation algorithm on mechanical layers of a clastic sedimentary succession. The comparison of automated segmentation results and compass measurements of three different plane sets shows a deviation of 0.70–2.00°, while the mean single plane divergence amounts to 4.97°. Hence, this study presents a fast, precise, and highly adaptable automated plane detection method, which is reproducible, transparent, objective, and provides increased accuracy in outcrops with rough and complex surfaces. Moreover, output formats of spatial orientation and location of planes are designed for simple handling in other workflows and software.

1. Introduction Various mining and civil engineering projects require extensive investigations of rock masses. To ensure safety and long-term reliability it is essential to characterize their mechanical behavior. While traditional data acquisition requires direct access to the rock face, terrestrial laser scanning (TLS) for generating digital outcrop models (DOMs) overcomes this necessity to some extent and has evolved into a standard method for describing and quantifying reservoir potential and structural characteristics of geological strata.1–4 To characterize the structural setting of rock masses the fracture pattern is a crucial property and is usually generated by statistical models such as discrete fracture networks (DFNs). The precision of these models is increased the more fracture data is available. Therefore, it is worthwhile to extract the inherent fracture data from TLS point clouds for such outcrop analog studies. Moreover, spatial information of other geological surfaces such as fault or bedding planes are of equal importance for rock mass characterization and can be extracted from 3D point cloud data as well. At many study sites, manual geological compass measurements and scan line surveys are not feasible or are restricted to small areas because ⁎

Corresponding author. E-mail address: [email protected] (T. Drews).

https://doi.org/10.1016/j.ijrmms.2018.06.023 Received 23 December 2016; Received in revised form 16 May 2018; Accepted 26 June 2018 1365-1609/ © 2018 Elsevier Ltd. All rights reserved.

of inaccessibility or rockfall hazard. In addition, manual measurements are time-consuming and error-prone, especially in complex settings. Plane analyses of 3D TLS data provide a tool to overcome all these problems. Another vital aspect is that each single orientation measurement obtained from TLS data can be easily associated with precise spatial information (x-, y-, z-coordinates), which is of great importance for straightforward import into modeling software and further data processing. Wide accessibility of close-range laser scanning techniques and its increased utilization has recently led to many studies dealing with discontinuity analysis from TLS-derived data.5–18 Despite the growing amount of publications, the authors of this study have been the first that focused explicitly on (semi-)automated discontinuity recognition through a 3D region growing segmentation algorithm based on a robust plane fit for normal estimation.19–21 Much too often a lack of extensive validation tests (“ground truth”) with conventional measurement techniques and other software can be observed in existing literature. Furthermore, little value has been attributed to passing on the results for subsequent DFN development in commercial software packages. Published methods regarding estimations of discontinuity spacing often lack the necessary precision for

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

point clouds has been widely accepted.15 Jaboyedoff et al.32 proposed the application of the principal component analysis (PCA) for the calculation of normal vector orientation of every point and its coplanar neighbors. This concept was adopted by several authors for plane identification in point clouds.15,23,25,26 Three recently published studies5,6,15 utilized normal vector calculation as a basis for subsequent discontinuity set recognition, performed directly on 3D point clouds. For this purpose Assali et al.5 applied a spherical k-mean algorithm on the normal vectors associated with each point to classify subsets of discontinuities. Riquelme et al.15 made use of a coplanarity test to remove anomalous points and applied a Kernel Density Estimation (KDE) to identify discontinuity sets. Most recently, Wang et al.6 published a region growing approach for an automated plane extraction, that is similar to the algorithm used in the present study. On the other hand, García-Sellés et al.8 performed a planar regression for normal vector estimation and subsequently applied a region growing segmentation algorithm, but their approach is not (semi-) automated. Other studies by5,27,33 recognized discontinuity spacing in point clouds. These studies apply a virtual scanline perpendicular to identified discontinuities belonging to one set. Gigli and Casagli9 applied a similar approach, but used an elongated cylinder instead of a scanline to cut the discontinuities.

significant DFN models and do not include all relevant discontinuities. Hence, these methods may not be sufficient to reflect the true fracture intensity. The originality of our study lies within several aspects. (1) We present the newest version of a constantly improved algorithm for automated detection of 3D planes, which has multiple advantages. The algorithm needs only a minimum of data preprocessing, if any. It works efficiently and fast without necessity of high processing power; hence, it is applicable on field laptops for in situ analysis. The algorithm is highly adaptable to specific applications or different outcrop conditions. In addition, the extensive output data is optimized for straightforward import into analytical and modeling software for further use. Finally, it delivers very precise and accurate orientation data in a short time. (2) Furthermore, we provide extensive validation tests for statistical significance by comparing our results with two alternative methods, which leads to new insights about ground truth fit and orientation data quality. (3) In addition, a workflow is proposed for precise fracture intensity (or spacing) calculations, which is based on orientation data extracted by the presented algorithm, applicable for fracture network modeling. Moreover, pitfalls of intensity sampling are highlighted. First, the automated plane detection algorithm is described in detail including data preprocessing steps. The algorithm performs a 3D region growing segmentation of point cloud data controlled by several parameters. To test the automated plane detection, an outcrop with intricately developed planes was chosen. More than 1000 geological compass measurements and orientation data, handpicked digitally with the software VRGS, were used to validate the results. By comparing the manual and digital methods, we reveal the statistical significance of plane set mean orientation calculations. Finally, we present a case study for digital fracture intensity estimation using the commercial modeling software Petrel.22

3. Study sites 3.1. A: Natural outcrop Wilkenfels This natural outcrop is a rock face in the nature sanctuary Russenstein in Heidelberg (located at UTM zone 32 U 479951E and 5474025N), which consists of biotite granite of Variscan age and is located close to the Upper Rhine Graben main boundary fault. It is intensely fractured with several visible minor faults orientated parallel to the main boundary fault. The most prominent surfaces are parallel quartz mineralized slickensides intersecting the whole outcrop with an average spacing of 1.2 m, partially giving the outcrop a layered appearance. The overall fracture pattern is very complex with several different fracture sets and sometimes irregularly orientated cracks. Blocks or surfaces often show spherical weathering. For this study, the westernmost part of the outcrop was investigated (Fig. 1). It is about 15 m long, 14 m wide and 8 m high. This site has several specific characteristics that are important for proper ground truth validation tests: (1) it provides full accessibility for compass measurements; (2) the TLS unit can be positioned at different levels around the target to minimize occlusion effects; (3) there are planes developed in almost all possible directions providing an optimal dataset to reveal systematic errors; (4) fractures are not concealed by quarry walls or other secondary surfaces in this natural outcrop; and (5) complex fracturing, rounded or curved planes and rough surfaces provide a challenging test for the segmentation algorithm. However, there

2. State of the art Several authors used the acquired TLS data for manual plane extraction.12–14,23,24 However, an automated identification of planes is desirable. Hence others introduced (semi-) automated methods on 3D point cloud derived 3D meshes10,18,23 and used a fuzzy k-mean clustering algorithm for plane set classification. However, due to their nature, 3D meshes are accompanied by smoothing and might therefore lead to an unwanted alteration of the original data and loss of detail. Because of the constraints mentioned above, the number of studies regarding the (semi-)automated identification of discontinuities in 3D point clouds has increased.5–7,9,11,15,17,23,25–27 Some authors proposed the use of least squares adjustment on a subset of points through planar regression28,29 or through a sampling cube in the point cloud,9,27 while others30 applied a procedure based on the Random Sample Consensus (RANSAC)31 algorithm by setting up an iterative process, in which the inliers of the best plane are removed from the data at each iteration until no more planes are found. The normal vector calculation on 3D

Fig. 1. (a) Upper section of the natural outcrop Wilkenfels. View in SW direction. Meditating geologist for scale. (b) Quarry Rockenau. View in W direction. 20

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

are some disadvantages: (1) the orientation of inconsistent or very rough surfaces is hard to determine precisely by compass measurements, and a time-consuming record of large statistically significant datasets is therefore necessary; (2) intense fracturing with many orientations leads to a wide divergence of the measurement results that impedes the definition of single plane sets.

plane detection were imported into Petrel. The format of a XYZ coordinate of the centroid for every detected plane was used to define fracture sets. In addition, the complete point cloud was imported containing the segment ID for every point defining its affiliation to a certain plane. These formats make it possible to define fracture sets, triangulate fracture planes and compute fracture set intensities.

3.2. B: Artificial quarry Rockenau

4.4. Manual orientation measurements

The still operating quarry (located at UTM zone 32 U 499833E and 5476264N) of the Lower Triassic (Bunter) consists primarily of sandstones with a lesser amount of silt and clay intercalations. The quarry walls predominantly consist of massive sandstone beds often showing cross bedding. The dimensions of the investigated area are: length ca. 80 m, width ca. 40 m and height ca. 38 m. Fractures frequently terminate at interbedded argillaceous and silty layers, representing incompetent beds and hence defining a distinct mechanical layering. Being a main geological control factor for fracture intensity, the mechanical layer thickness is of crucial importance for extrapolating fracture statistics.34–36 Therefore, the quarry serves as an ideal object for outcrop analog fracture spacing analysis. As an example, for the present study a sandstone layer between two siltstone beds was chosen whose outcropping fractures are almost all represented by fracture faces. For this reason, no additional examination of fracture traces is needed in the presented workflow.

We took 1196 individual geological compass measurements of 259 planes at the Wilkenfels outcrop for ground truth validation. The number of measurements per plane depends on its size, continuity and surface roughness. All measured planes were registered and their exact positions were marked in digital outcrop photos on a laptop during the field work. The measured dips and dip directions of each plane were converted into vectors to calculate the mean orientation, mean orientations resultant length, normalized resultant length, and apical angle of confidence cones for each plane. Subsequently, mean vectors were converted back into dips and dip directions. Errors resulting from bipolar (i.e. partly overtilted planes) or widely scattered (i.e. rough surfaces or curved planes) orientation sets were detected by a threshold for minimum resultant length of 0.95·n to 0.985·n, depending on size of n, where n is the number of orientations. For these data sets, the mean orientations were determined by the eigenvectors of orientation-distribution matrices from dips and dip directions.37–40 Steeply bipolar dipping planes were averaged by the orientation-distribution matrices of their normals to avoid errors of the calculated dip direction.

4. Methodology 4.1. Terrestrial laser scanning

4.5. Digitally handpicked orientation measurements Point cloud data is acquired by the TLS system Optech ILRIS-HD with a laser pulse repetition rate of 10 kHz and a measurement range from 3 m to 1800 m within a 40° x 40° field of view. For each point cloud, a corresponding high-resolution digital image was taken for subsequent photo-texturing. The study site “A: Wilkenfels” was covered by 15 point clouds from four different positions, with a mean scan distance between 3 and 25 m. The average point spacing was 12 mm at the mean distance of each scan. Full accessibility of the outcrop allowed for a good vertical and horizontal distribution of TLS positions in order to minimize occlusion effects. The test area “B: Rockenau” is built up of 5 point clouds taken from a single position with an average point spacing of 10 mm at distance of 15 m. Due to its advantageous morphology occlusion effects are negligible. Further processing of the datasets comprises point cloud co-registration (mean registration error of 0.4 mm), merging of point clouds and photo-texturing. To ensure best comparability of compass clinometer and (semi-)automated measurements, the merged point clouds were cropped according to the area covered by field work.

A number of 122 planes were digitally handpicked in the point cloud of the Wilkenfels outcrop for single plane comparison using the VRGS software published in.41 Depending on the respective size, varying amounts of points (a minimum of three) were picked on a section of the point cloud representing a plane. The software, subsequently, uses these points to compute a plane polygon, which in turn provides spatial orientation information for comparison with other results. 4.6. Statistical methods for comparing rotationally symmetric orientation distributions To test if two corresponding plane sets of compass survey and automated plane detection data have a common mean direction, a graphical method presented in37 is adopted. This method uses a quantilequantile uniform probability plot, where quantiles of the pooled samples are plotted against the uniform quantile. For this purpose the pooled weighted mean directions of the surface normals (θw , ϕw ) are calculated in polar coordinates (colatitude θ, longitude ϕ). In this case, where the spherical standard errors of the mean directions σi and the numbers of orientations ni for the ith compared plane set are not very different, weights based on sample size only, can be used: wi = ni/N, where N is the number of orientation of the pooled sample. Afterwards the deviations of all normals (θ*, ϕ*) from (θw , ϕw ) are obtained by multiplying the 3 × 3 rotation matrix:

4.2. Georeferencing of digital outcrop models The Wilkenfels outcrop data was georeferenced by determining the coordinates of fifteen reference points distributed over the study site via a total station (Leica TPS 700) to ensure high accuracy. In the Rockenau quarry the position of the TLS unit was determined by dGPS and used as reference. Its spatial orientation was assessed by clinometer measurements. In addition, two boards (50 × 50 cm) were integrated into the point clouds. The spatial orientation of those two boards was obtained by geological compass measurements for more precise alignment of the reference scan.

⎛ cos θ w cos ϕw cos θw sin ϕw − sin θw ⎞ ⎟ A (θw , ϕw ,0)=⎜ − sin ϕw cos ϕw 0 ⎜ ⎟ ⎝ sin θwcosϕw sin θw cos ϕw cos θw ⎠

(1)

with each normal vector. Under the hypothesis that all mean directions are the same, the N azimuths ϕ1*…ϕN* should be approximately uniformly distributed on (0, 2π). This means that the sample quantiles ϕn*/2 π plotted in ascending order against their uniform quantiles of the standardized uniform distribution on [0, 1] Ui = (i-0.5)/n for ith

4.3. Data processing The main point cloud processing task is the 3D plane detection, which is described in detail in Section 5. The results obtained from 21

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

neighbors of the seed points for the defined homogeneity criteria and starts to grow with the seed points having lowest curvature. These criteria comprise (1) the deviation of a point's normal vector to the adjusted plane's normal, (2) the distance of candidate points to the plane as well as to segmented points, and (3) the curvature and plane sigma values of candidate points. Both are assessed as percentiles and prohibit the exceedance of specified thresholds. Plane segments are completed once all neighboring points have been checked. Subsequently, the principal orientation of segments is determined via an adjusted plane through all points belonging to a segment and is then assigned to each point of the respective segment. Points are discarded if they have not been classified as part of a plane or if they are part of a plane that does not meet the requirement of minimum segment size. Parameters used for plane detection in this study are specified in Table 1.

Table 1 Specification of parameters for plane segmentation used in this study. Parameter description

Value

Maximum angle between the normal vectors [°] Maximum distance of point to the plane [m] Maximum distance of point to the segmented points [m] Maximum allowed curvature as percentile [%] Point selection threshold as percentile of roughness [%] Minimum segment size for output [no. of points]

Wilkenfels

Wilkenfels adjusted

Rockenau

20

5

35

0.07

0.05

0.08

0.05

0.05

0.5

90

75

85

90

75

90

800

600

1200

6. Results 6.1. Validation of orientation measurements and statistical significance

sample should result in approximately linear plots along x = y if the hypothesis is true. Furthermore, possible deviations in the overall precision were checked by comparing the concentration parameters κ of the Fisher distribution. The concentration parameter κ is a shape parameter of the Fisher distribution. It is a measure of the degree of preferred orientation within a rotational symmetric population of orientations. κ is zero for a uniform distribution and increases without limits if the clustering of orientations improves. In this study, two different approximations for κ are used. Given that R/n ≥ 0.95, κ is approximated by see42,43:

n−1 n−R

κ=

To compare plane orientations from automated plane detection with ground truth data and data generated with the VRGS software, respectively, two different procedures were implemented. First, the orientations of individual planes are compared. Afterwards, sets of congeneric planes are determined within the orientation data of the compass survey and automated plane detection. Identified corresponding plane sets are extracted and used for comparative statistical analyses. 6.1.1. Individual plane equations By marking the outlines of all measured planes in digital outcrop pictures on a laptop during the compass survey at outcrop Wilkenfels, these planes can be compared with individual digitally generated segments. The segmentation algorithm was deliberately carried out with broad parameter minima and maxima thresholds to capture all plane sets in one analysis (Fig. 2). Therefore, some unwanted planes as vegetation-covered sections were not completely excluded by the algorithm. Such planes mostly correspond to the approximate overall orientation of the slope, their appearance can be observed in Fig. 3 by a slightly overrepresentation of poles to planes in the fourth quadrant, compared to the compass measurements displayed in Fig. 4. A total of 122 corresponding planes were identified and their orientations compared. Results are listed in the Appendix A. In some cases, two or more planes or segments were merged to achieve a better match of outlines for its equivalent counterpart (hereafter referred to as congruency). Seven orientations of the automated plane detection and fifteen orientations of the compass survey were obtained by merging. The varying congruency of corresponding planes was classified into 5 ranks through visual comparison (Table 2). The location of planes picked in the point cloud with the VRGS method was defined with the help of its corresponding segments generated by the algorithm. However, the extents of planes were not adjusted exactly to the segmented counterpart but, instead, interpreted independently. Nonetheless, this suggests a higher match of orientations between these two methods. Yet, only a slight correlation between the ranks of congruency and the match of orientations between compass and automated plane segmentation measurements could be observed, and 80% of planes are among the first three ranks with similar deviation angles (Table 2). Therefore, this effect seems not to be substantial. The calculated mean deviation angle αAC between the compared orientation of the plane detection and the compass survey is 4.97°, angle αAV between plane detection and in VRGS picked planes is 1.79°, and αCV spanned by the means of compass survey planes and VRGS planes is 5.22° (see Appendix A). The relatively small angle αAV in comparison with the significant larger and approximately equal

(2)

and for sample sizes n < 16 by37,44:

(n−1)2 n (n − R)

κ=

(3)

where n is the number of orientations and R is the resultant length, or the magnitude of the mean vector. 5. Automated plane detection algorithm A segmentation algorithm for the application on geological outcrops has been developed for the detection of planes in the acquired point clouds.19,21 The algorithm segments the point cloud by performing 3D region growing regarding specified homogeneity criteria in the neighborhood of seed points. Based on the assumption that low curvature neighborhoods are most probable to make up part of a plane, seed points are sorted in order of increasing curvatures in respect of their 3D neighborhood. The growing of segments is limited mainly by a maximum angular deviation of surface normals. As input for the plane detection algorithm, the data is preliminarily processed by computing the local plane information, i.e. surface normals, for every point considering their full 3D neighborhood using the modular program system OPALS (Orientation and Processing of Airborne Laser Scanning data).45 In this process, local normals are estimated through a robust plane fit to the nearest neighbors of each point in a specified 3D radius. An output of the normals estimation further used in the plane segmentation is the standard deviation for the plane fit, which can be understood as the roughness of the local surface and is subsequently addressed as plane sigma. The sorting of seed points in the algorithm is performed by initially computing the curvature (Ce) for each point according to its respective eigenvalue (e1, e2, e3) ratio46:

Ce =

e1 , e1 + e2 + e3

e1 < e2 < e3

(4)

Subsequently, the algorithm checks a specified number of nearest 22

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

Fig. 2. (a) Initial point cloud of the natural outcrop Wilkenfels. (b) Segmentation of the Wilkenfels point cloud with “Wilkenfels” parameters (see Table 1). Individual segments are distinguished by random colors. Segments of same color are not related. View in N direction.

Fig. 4. Poles to planes in equal area stereographic projection of compass survey results at the outcrop Wilkenfels. Blue poles: Fracture set 1 with mean direction (black square) and confidence cone (black circle); Green poles: Fracture set 2; Purple poles: Set 3.

Fig. 3. Poles to planes in equal area stereographic projection of segmentation algorithm results at the outcrop Wilkenfels. Blue poles: Fracture set 1 with mean direction (black square) and confidence cone (black circle).

compass-dependent angles αAC and αCV indicates that the two digital methods have a lower random error. The fit of orientations, generated with the VRGS software, is also visually checked by fitting discs to the point cloud. In contrast, discs created with the divergent orientations of the compass survey often fail this visual test. This affirms the previous assumption that the digital methods have a lower random error. The requirement for this conclusion is that no significant spatial orientation error of the used point clouds exists. Since no systematic error was noticed (see Section 6.1.2), this requirement is considered as fulfilled. A few outliers of αAV indicate that the spatial orientation of some planes calculated by the algorithm may differ slightly from the true mean orientation. This is the case for curved planes with a variable point density, which leads to a concentration of points in one specific part of the plane. Through equal weighting of every point, the part with higher point density will dominate and therefore shifts the planes’ main orientation. Although the

Table 2 Congruency of algorithm segments compared to compass survey planes. Ranks of congruency: (1) very high besides contours; (2) high except marginal sections; (3) high at central sections; (4) moderate but still more than 50% of total planes area; (5) partial. Rank of congruency

1

2

3

4

5

Total

Number of planes Mean deviation angle [°]

28 3.82

42 4.74

27 4.89

18 6.83

7 6.43

122 4.97

algorithm generally works well with heterogeneous point densities, this rare case can be avoided by re-sampling the point cloud to harmonize its spatial point density, using such methods as an Octree filter for choosing a random point per voxel or the center of gravity for each voxel. 23

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

6.1.2. Comparing plane sets statistics For all three measuring methods discussed above the dimension of segmentation, varying deviation tolerances, and ambiguous plane boundaries can cause different individual orientations, which does not necessarily imply that one of the results is less correct. This applies in particular to outcrops with complex surfaces where several uncertainties in plane interpretation are common. Therefore, besides comparing individual orientation data of compass and automated point cloud based measurements, it is essential to test the overall consistency of the two methods. For this purpose, we clustered the measurements into congeneric plane sets, which can be statistically compared. Additionally, in most practical cases, mean orientations of plane sets are more relevant than the orientation of individual surfaces, making these tests even more relevant. Overlapping orientations of different plane sets and arbitrary orientated weathering surfaces impede the identification of plane sets. Nevertheless, three usable plane sets were identified in both datasets for statistical analysis. Set 1 consists of the most prominent surfaces of the outcrop Wilkenfels, formed by consistent slickensides described in Section 3.1. Because of their wide spacing and high consistency only 11 surfaces and 12 segments, respectively, of these slickensides were recorded with each method (Figs. 3 and 4). Set 2 and set 3 are both steep dipping fracture sets. For enhanced capturing of these sets and to minimize interference with other planes, the automated plane detection was carried out with adjusted parameters (Table 1 and Fig. 5). Due to the bipolar distribution of both sets, the data was rotated for better visualization in the stereographic projections to allow correct vector analyses (Figs. 6 and 7). For a first comparison, the mean orientations were calculated and their semi-apical angle θα of the conical confidence interval with a confidence level of 95%, corresponding to a significance level of α = 0.05 was approximated with see37,43,47:

Fig. 6. Rotated poles to planes in equal area stereographic projection of two fracture sets from the segmentation algorithm at the outcrop Wilkenfels. Red: Set 2; blue: Set 3. Rotation axis 85°, magnitude of rotation 60°.

1

⎞ ⎛ n−R 1 ⎞ × ⎛ ⎛ ⎞ (n −1) −1⎞ θα = cos−1 ⎜1−⎛ ⎟⎟ ⎝ R ⎠ ⎜⎝ ⎝ α ⎠ ⎠⎠ ⎝

(5)

Lastly, the angles spanned by the mean orientations of corresponding sets were calculated. As can be seen in Table 3, all angles

Fig. 7. Rotated poles to planes in equal area stereographic projection of two fracture sets from the compass survey at the outcrop Wilkenfels. Red: Set 2; blue: Set 3. Rotation axis 85°, magnitude of rotation 60°.

spanned by corresponding means are very small (0.7° to max. 2°) and all associated θα values exceed them by far. Accordingly, confidence regions do not just overlap but also mutually enclose each other's mean directions (Figs. 7 and 8). Hence, it can be concluded that the mean directions are similar at a 95% confidence level. Even though our first test already delivered very high congruency between the mean orientations of compass- and LiDAR-based plane sets, a graphical test was additionally performed to check whether the

Fig. 5. Poles to planes in equal area stereographic projection of segmentation algorithm results at the outcrop Wilkenfels with adjusted parameters. Green poles: Fracture set 2; Purple poles: Set 3. 24

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

Table 3 Comparative chart of same plane sets mean orientation and their deviation, obtained by automated plane detection algorithm and compass measurements. Plane set number

1 2 3

compass survey

automatic analysis

angle between mean orientations [°]

mean orientation [°]

R/n

confidence cone [°]

mean orientation [°]

R/n

confidence cone [°]

237.16/63.27 346.8/1.0 136.8/4.1

0.9895 0.9841 0.9874

3.1 3.1 3.8

235.75/63.05 167.0/0.6 138.5/5.1

0.9962 0.9862 0.9854

2.8 2.7 4.1

0.7 1.6 2

by,37 p. 114). The graphical method also uses quantile-quantile uniform probability plots and the assumption of rotational symmetry is implausible if the resulting plot is not approximately linear along x = y (Fig. 8d). Nevertheless, Fig. 8c corroborates the statistical match of mean orientations with good approximation. All statistical tests carried out in order to compare the mean orientations of single clusters show no significant variances, although the precision or number of single measurements may vary due to errors in measurements or more likely differences in segmentation. This is also valid for clusters with small sample sizes as shown for set 1 with down to 11 orientations for one sample. Furthermore, the tests reveal high accuracy for those clusters; i.e. there were no significant systematic errors occurring from spatial calibration of the point cloud data. Possible deviations in the overall precision can also be checked by comparison of the concentration parameters κ. The summarized results in Table 4 show that associated κ values are of similar magnitude and thus the precision can be assumed to be similar for both methods.

data pairs have a common mean direction. The method of this test is described in Section 4.6 and presented in.37 If the hypothesis is true that the compared plane sets have a common mean direction., the resulting quantile-quantile uniform probability plots should be approximately linear along x = y. At first, set 1 was tested. Despite the relatively small size of the pooled sample (n = 23), the uniform probability plot (Fig. 8a) illustrates a very clear result. All data points are scattered closely around x = y without any outliers and the regression line shows only a negligible deviation. Plotted data of set 2 (Fig. 8b) with a much larger pooled sample size (n = 72) also correlates well with the expected distribution. The probability plot of set 3 (Fig. 8c) with n = 38 appears again reasonably linear, but apparently most of the data points lie above the x = y line. This is probably due to the fact that this statistical test is designed for orientation data with rotationally symmetric (or Fisher) distributions. The orientation distribution of compass survey data of set 3 does not fit these criteria perfectly as a performed test for rotational symmetry about the mean directions indicated (test described

Fig. 8. (a), (b) & (c) Quantile-quantile uniform probability plots to test whether the same fracture sets obtained with compass measurements and the automated plane detection have common mean directions for set 1 (a), set 2 (b) and set 3 (c). (d) Quantile-quantile uniform probability plot to test compass survey data of set 3 for rotational symmetry about the mean direction. Dotted black lines; regression lines; red dashed lines: 45° lines through origin. 25

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

Table 4 Comparative chart of concentration parameters for plane sets obtained from compass measurements (κcompass) and automated plane detection (κauto).

κcompass κauto

Set 1

Set 2

Set 3

220 187

61 71

75 65

6.2. Layered fracture intensity estimation for DFN modeling In this study, the fracture intensity is defined by the number of fractures of one set encountered per unit length along a sample line (scanline) perpendicular to the mean orientation of this set. This crucial parameter for DFN modeling is also referred to as P10 value.48 In literature, the commonly reported average spacing is the inverse of the fracture intensity.36 For intensity estimation, series of strata representing mechanical layers were chosen as targets and extracted. The mechanical layer thickness is a main controlling factor for fracture intensity. Since these two parameters often show a proportional correlation,34,35,49 it is possible to extrapolate and compare results within the surveyed lithology. The examples presented below are excerpts of a regional study on the Triassic Bunter Formation, which demonstrate a technique for gathering accurate fracture intensities solely based on point cloud data and digital outcrop images. All presented intensity calculations arise from the same mechanical layer at outcrop B “Rockenau”. The first analyzed section is 16 m long and every fracture is developed as an outcropping plane (Fig. 9a). For this reason, no fracture traces had to be included for intensity estimations. Section 2 covers 21 m of the example layer and lies 14 m north of Section 1 (Figs. 9a and 9b). For one fracture set an additional interpretation of possible fracture traces was necessary. This was accomplished by means of high-resolution photos projected on meshed point clouds. Scanline orientations of the two presented sections form an angle of approximately 60°. Greater angles between scanlines enhance the chance to capture all present fracture sets and minimize bias. The intensity analysis is based on the TLS orientation data and carried out with the commercial modeling software Petrel. Hence, imported and generated data in Petrel can subsequently be used for straightforward DFN modeling. For example, fracture intensity logs can be created on the basis of scanlines with precise spatial information. Tests performed with automated approaches such as clustering algorithms for plane set identifications led to incorrect or overly imprecise results. This is shown in the plot of set classification (Fig. 10). Here it can be clearly seen that clustering algorithms based on fracture orientation only lead to erroneous results. Although filter criteria for false plane elimination were determined (described below), an expert's interpretation is still needed to guarantee accurate results. The filtering process of false planes begins with choosing ideal parameters for the segmentation algorithm. A step-by-step workflow

Fig. 10. Poles to planes in equal area stereographic projection of two fracture sets from the automated plane detection of Section 1 at outcrop Rockenau. Red and blue poles belong to two differed fracture sets and green poles to undefined or “false” planes.

was designed to tackle the false planes during the process. In the example presented a sufficiently high minimum segment size was of particular importance. Secondly, rough contours of clusters are interpreted manually in the stereographic projections. Those clusters still contain false planes that tamper fracture intensity calculations. These are primarily edges or marginal fractions and do not transect the whole layer, in contrast to relevant fractures. Therefore, the centroids of these planes are usually significantly above or below the fracture sets centroids and have extreme values regarding the z-coordinate. By means of a logical filter determining z-coordinate outliers, most of the remaining false planes can be eliminated. For both scanlines of Section 1, at this step, all false planes were filtered except for one. This remaining false plane results from a fracture divided into two segments through lateral displacement along a clay lens (on the far right of Fig. 9b). In this case, one of the two redundant planes had to be eliminated manually. Furthermore, care must previously have been taken to ensure not to completely eliminate such partially displaced fractures with the z-value filter or the minimum segment size parameter. (Fig. 11) In Section 2 the same two fracture sets are present as in Section 1,

Fig. 9. Steps for estimation of layered fracture intensities at the outcrop Rockenau. (a) Point cloud with photo projection of investigated horizons Section 1 (between red brackets). (b) Segmented point cloud with blue and red planes belonging to two different fracture sets and undefined green planes visualized in Petrel. The undefined planes represent edges or marginal fractions except the extreme right one, which is a part of a fracture divided through lateral displacement. (c) and (d) extracted fracture sets visualized as discs around plane centroids with white scanline. 26

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

Fig. 11. Example Section 2 of investigated horizon at outcrop Rockenau. (a) Point cloud with photo projection (Section 2 between red brackets). (b) Segmented point cloud with blue and red planes belonging to two different fracture sets and undefined or redundant green planes visualized in Petrel.

even though the mean orientations seem to be a little bit shifted (Fig. 12). One set is very well clustered without any false planes, which makes filtering unnecessary. However, in addition five fracture traces belonging to the same set were identified and picked in the DOM. Fig. 13 shows that the fracture traces lie within the range of plane orientations. The second set has a relatively low incident angle (about 19°) to the rock face. Therefore, several redundant planes are within the cluster of the set (Fig. 12) but with predominantly peripheral positions. As in Section 1 all false planes were filtered by means of filter determining z-coordinate outliers (6 planes) or a combination of z- and yposition (3 planes) except for one. The latter is a two times outcropping fracture and one of its segments had to be eliminated manually. The endpoints of each scanline are defined by the centroids of the two outermost fractures of each set (Figs. 9c and 9d). Since for intensity estimation the true scanline length L90 must be measured perpendicular to the mean orientation of the fracture sets,50,51 the calculated length is corrected by:

L90 = LS cos α

(6)

where α is the intersection angle between scanline and mean normal orientation of the fracture set and LS is the measured scanline length. For the sake of completeness, the following fracture intensities were calculated for the example layers’ two sections: 1.03 and 0.94 fracture per meter for set 1 and 1.25 and 1.23 for set 2. Fig. 12. Poles to planes in equal area stereographic projection of two fracture sets from the automated plane detection of Section 2 at outcrop Rockenau. Red and blue poles belong to two differed fracture sets and green poles to undefined or “false” planes.

6.3. Processing time for plane detection To run the segmentation algorithm no special processing power is required. The processing speed primarily depends on the size of the point cloud used and secondarily on the parametrization. The most time-consuming operation of the whole plane detection workflow is the normal calculation in OPALS. For the Wilkenfels point cloud with 2.03 million points the normal calculation took 34 min on an Intel Xeon X 5647, 24 GB RAM workstation. The following entire segmentation process took only 120–130 s (depending on chosen parameters) on the same device. 7. Conclusions The interest in extracting orientation information from TLS 3D point clouds has increased rapidly over the past few years. For our study multiple advantages can be pointed out. A novel solid 3D plane detection algorithm that is highly adjustable through several freely selectable input parameters is presented. This novel algorithm can be adapted to multiple conditions like weathering phenomena, varying curvature, roughness, dimensions of fractures, point densities and noise of point clouds. This is shown by the precise results of the presented segmentations applied to rock masses with complex fracture patterns. In addition, a first filtering of unwanted planes can be realized through parameter adaption. The segmentation algorithm works efficiently and fast without the necessity of high processing power; hence, it is also applicable on field laptops for in situ analyses. Absence of extensive preprocessing like meshing prevents distortion or loss of information. The initial point cloud information is kept during plane detection but to each point additional attributes are added (segment ID, orientation, roughness, curvature, and number of points used for normal

Fig. 13. Poles to planes (red dots), great circles of planes (red circular arcs) and fracture traces (blue triangles) of fracture set 1 of Section 2 at outcrop Rockenau in equal area stereographic projection.

27

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

Finally, although the applied ground truth tests delivered statistically clear and consistent results, additional tests are needed for confirmation. Unfortunately, statistically significant tests are very timeconsuming. Mapping precise compass measurements specifically is the hardest and most time-consuming task, since angles within complex planes can vary for several tens of degrees and therefore even numerous measurements may not ensure high precision of orientation data. This also points out that, besides the usually mentioned advantages of high sampling speed and no limitation in accessibility of TLS based plane analysis, increased accuracy is a strong argument for its application in outcrops with rough and complex surfaces. Expert knowledge is still needed for interpreting the origin of segments or plane sets and to ensure that false planes like weathering surfaces or grassy slopes are not included into fracture analysis. When it comes to DFN modeling, fieldwork is still needed in most cases, since some fracture attributes such as aperture cannot be extracted satisfactorily with digital methods. Therefore, the presented method and workflows should not be seen as standalone techniques but as a chance to significantly facilitate fracture sampling and increase data quality.

estimation). This process minimizes unwanted bias in the dataset. Another key finding of this study is the excellent ground truth of the presented algorithm validated by extensive statistical comparison with compass data, even though there are some outliers. Orientation data gathered by picking planes with the VRGS software showed an even better agreement to the algorithm-based orientations, indicating that compass measurements were less precise than the two digital methods. There are several reasons, such as roughness, curvature or accessibility, that explain why complex surfaces are sometimes difficult to measure precisely with the geologic compass, even though time-consuming multiple measurements are carried out. On the other hand, the algorithm sometimes captures unwanted sections, e.g. areas covered with vegetation. Inhomogeneous point density can affect curved plane orientations negatively in rare cases. For outcrops with considerably curved planes, a re-sampling of the point cloud is recommended to harmonize its spatial point density in order to avoid shifts of calculated orientations. The fracture set statistics performed reveals that the differences of single measurements are balanced for whole sets. All mean orientations of a certain set calculated from compass and automated plane detection data lie within each other's confidence limits. This result is of even greater importance for the proposed application of DFN modeling, where the mean orientations of fracture sets are of primary relevance. The comparison of set statistics also reveals high precision and accuracy of the TLS based orientation data. Significant systematic errors can therefore be excluded. In recent years, several automatic or semi-automatic approaches for intensity estimations have been published. Tested automated methods greatly risk incorporating false planes into the analysis and do not consider the mechanical layer thickness. The presented layer-wise analysis is carried out on a sedimentary sequence showing a defined mechanical layering. It combines manual fracture interpretation within point cloud data with digital plane filtering tools. Except for very few cases, in this study, false planes were able to be filtered by adjusting algorithm parameters; spatial outliers of fracture centroids were also able to be filtered out. Precision can be maximized and fast straightforward DFN modeling is possible with this method.

Acknowledgements We appreciate the support and partially funding of the German Federal Ministry for Economic Affairs and Energy (BMWi), previously of the German Federal Ministry for the Environment, Nature Conservation and Nuclear Safety (BMU) as part of the Research and Development Projetct AuGE (Outcrop Analog Studies in Geothermal Exploration) within the framework of the 5th Energy Research Program (grant number: 0325302). The authors thank Schlumberger for the support of academic Petrel licenses and Janine Lange for total station measurements. We also would like to thank the Regierungspräsidium Karlsruhe for the authorization to enter the nature sanctuary Russenstein in Heidelberg and Hermann Graser from Bamberger Natursteinwerk for the permission to work in the Rockenau quarry. Last, the authors would like to acknowledge the reviewers’ input and Annette Stemmerich for providing English language revision.

Appendix A See Appendix Table A1

Table A1 Comparative table of orientations of a single plane gathered with automated plane detection, compass measurements and handpicking tool in VRGS software. Deviation angles are angles spanned by αAC: automated plane detection and compass survey planes; αAV: automated plane detection and VRGS plane; αCV: compass survey and VRGS plane. Automated (LiDAR)

Compass survey

VRGS

deviation angles [°]

Dip dir. [°]

Dip [°]

n

Dip dir. [°]

Dip [°]

Con.

Dip dir. [°]

Dip [°]

αAC

αAV

αCV

273.14 81.91 257.33 321.84 54.97 76.42 324.43 51.22 274.2 108.47 138.46 111.02 302.78 268.42 137.82 348.38 94.76 300.54

75.3 85.17 78.62 73.05 82.24 84.33 75.77 25.91 84.34 32.18 81.45 69.24 85.29 83.36 38.87 62.98 87.66 68.44

8 4 3 8 4 2 6 6 12 3 5 3 5 4 5 3 1 1

267.02 80.67 249.41 321.79 57.08 79.56 324.5 48.82 272.85 104.7 140.26 104.28 292.94 263.9 129.63 335.67 275 300

77.08 86.264 76.72 75.55 81.65 85.5 81.06 27.23 77.65 35.03 80.64 70.55 83.41 82.52 33.05 55.68 84 64

1 1 1 1 2 1 3 1 2 2 3 2 2 2 2 4 1 2

272.77 80.98 255.71 323.03 56.46 74.22 323.07 53.82 276.47 109.82 140.1 110.89 303.45 269.5 139.28 349.23 93.1 299.02

76.42 84.17 77.03 72.8 83.39 83.12 74.54 24.37 84.9 34.03 79.25 68.95 85.14 85.4 34.4 63.19 87 68.67

6.20 1.65 7.97 2.50 2.17 3.34 5.29 1.70 6.82 3.53 1.95 6.46 9.97 4.56 7.54 13.13 8.34 4.47

1.18 1.36 2.24 1.16 1.87 2.50 1.80 1.89 2.33 1.99 2.73 0.31 0.68 2.31 4.55 0.79 1.78 1.43

5.64 2.12 6.14 3.00 1.85 5.82 6.67 3.59 8.08 3.07 1.40 6.40 10.60 6.27 5.52 13.86 9.20 4.76

(continued on next page) 28

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

Table A1 (continued) Automated (LiDAR)

Compass survey

VRGS

deviation angles [°]

Dip dir. [°]

Dip [°]

n

Dip dir. [°]

Dip [°]

Con.

Dip dir. [°]

Dip [°]

αAC

αAV

αCV

57.61 117.74 212.77 100.51 96.66 95.86 108.68 256.71 249.96 65.9 27.41 105 187.2 90.62 84.43 316.25 51.08 271.78 28.27 276.7 338.58 166.9 167.33 18.98 180.21 357.13 193.73 356.04 346.27 332.05 20.13 171.53 140.15 142 134.93 187.57 179.21 139.83 180.12 188.37 162.45 160.53 189.55 184.86 138.19 228.2 180 164.86 182.66 200.41 177.01 257.83 53.98 75.55 65.52 78.25 133.47 137.79 186.1 99.13 148.28 166.63 315.52 135.85 128.07 330.73 165.08 159.21 105.25 56.82 96.13 92.42 90.87

23.32 79.17 10.17 79.89 85.66 86.51 74.93 71.58 89.89 23.94 68.96 46.02 29.84 76.18 47.48 84.58 27.76 72.06 85.58 72.61 84.42 24.9 17.64 84.01 89.99 88.8 80.12 87.67 85.64 83.94 86.44 32.01 83.49 30.25 20.74 22.02 31.55 87.37 34.31 26.16 20.29 16.39 24.38 30.36 16.22 27.9 21.93 84.27 55.02 76.98 80.14 79.08 30.3 59.9 32.7 80.55 21.94 24.76 26.26 17.08 22.53 18.1 84.17 89.81 86.72 86.07 89.26 88.42 85.23 27.06 85.63 87.39 89.68

2 3 7 8 3 3 7 6 6 3 5 5 6 3 2 26 12 4 5 3 10 15 4 3 4 4 4 6 4 4 2 3 5 8 8 18 10 3 10 14 3 10 8 3 5 3 13 5 4 3 6 2 9 6 18 4 5 20 11 3 4 14 2 7 11 10 28 6 5 7 3 4 3

52.86 112.1 244.91 97.8 281 94.14 102.47 253.88 76.1 64.65 25.02 101.38 202.58 85.33 85.5 318.68 49.18 267.55 25.74 277.86 339.37 181 189.49 19.5 1.5 357.3 194.73 352.49 344.31 333.11 198.2 172.55 318.43 155.57 125.41 189.9 188.79 319.8 184.89 201.6 161.55 156.9 200.61 192.33 134.26 230.82 170.6 169.37 180.86 201.4 171.27 257 54.75 73.25 62.75 77.14 130.3 142.9 174.5 95.92 143.95 173.2 317.56 311.5 130 327.13 344.7 162.7 102.4 56.98 100 271.24 89

25.91 83.01 11.25 85.59 88.67 85.34 80.75 70.88 87.9 26.69 70.21 51.62 33.55 76.67 49.01 84.05 29.61 70.01 82.04 67.33 86.21 25.78 18.84 82 86.22 87.5 82.51 80.45 86.76 84.06 87.5 37.37 87.2 32.25 20.58 28.83 34.08 89.83 42.03 26.56 29.94 15.46 26.7 28.7 23.09 34.62 21.76 86.31 47.11 78.35 86.34 80 28.49 64.01 35.07 82.95 22.43 23.5 24.85 20.11 21.35 20.05 85.5 87.57 89.73 83.16 82.87 89.5 86.38 26.9 90 88.5 89.33

2 1 4 2 2 1 2 2 2 3 1 2 4 5 1 1 1 1 1 2 1 2 3 3 2 2 3 4 2 1 1 3 5 3 2 3 1 3 5 2 5 2 3 4 4 2 4 5 3 1 2 1 1 1 3 3 2 2 1 2 3 2 1 3 4 5 3 2 1 2 2 2 3

58.01 118.76 217.34 101.36 97.31 97.22 108.84 257.61 69.04 65.02 27.89 104.93 195.68 91.08 85.34 316.43 52.55 272.08 28.62 275.76 338.12 152.92 165.14 17.83 181.87 359.19 193.67 354.9 345.69 332.89 201.35 175.26 139.74 142.58 137.16 186.36 178.96 319.94 181.92 187.82 168.51 165.85 187.2 182.87 138.96 227 167.92 166.53 180.07 198.9 176.96 258.28 54.16 76.8 64.2 79.2 128.32 137.17 190.26 96.73 147.78 167.6 314.71 135.77 127.26 331.24 164.13 159.57 104.28 57.23 96.32 92.14 91.08

23.7 80.32 10.95 80.5 86.79 85.08 77.4 71.79 89.78 23.39 68.06 49.92 30.91 75.56 48.66 83.18 27.9 69.66 84.96 71.31 86.25 24.33 20.02 86.18 89.46 88.62 80.81 88.66 87.06 83.19 85.5 33.75 83.27 27.04 19.65 22.78 30.49 89.23 34.41 26.2 20.02 17.72 25.13 28.67 15.74 25.96 24.76 85.81 57.64 79.05 80.61 80.7 30.31 58.92 31.86 78.49 20.96 25.04 17.88 16.06 22.42 17.94 83.54 88.43 84.34 84.8 87.68 89.49 85.48 26.94 84.05 86.81 89.63

3.26 6.77 5.99 6.30 7.14 2.08 8.41 2.77 6.52 2.80 2.56 6.23 8.86 5.17 1.73 2.48 2.06 4.49 4.34 5.39 1.96 6.09 7.00 2.08 4.00 1.31 2.59 8.04 2.25 1.06 6.36 5.39 9.47 7.30 3.36 6.88 5.77 2.80 8.26 5.88 9.66 1.36 5.29 4.04 6.99 6.85 3.50 4.94 8.03 1.68 8.42 1.23 1.85 4.58 2.83 2.64 1.29 2.44 5.19 3.20 2.00 2.90 2.43 5.08 3.57 4.62 7.88 3.65 3.07 0.18 5.83 4.28 1.90

0.41 1.53 1.14 1.04 1.30 1.97 2.47 0.88 0.98 0.65 1.00 3.90 4.42 0.76 1.36 1.41 0.70 2.42 0.71 1.58 1.89 5.84 2.48 2.45 1.74 2.07 0.69 1.51 1.53 1.12 8.15 2.67 0.46 3.22 1.33 0.89 1.07 3.40 1.02 0.25 2.10 2.05 1.24 1.95 0.52 2.01 5.55 2.27 3.39 2.54 0.47 1.68 0.09 1.46 1.10 2.26 2.12 0.38 8.52 1.23 0.22 0.34 1.02 1.38 2.51 1.37 1.84 1.13 1.00 0.22 1.59 0.64 0.22

3.09 7.12 5.27 6.20 5.85 3.08 7.09 3.65 7.30 3.30 3.44 3.23 4.53 5.69 0.37 2.40 2.36 4.27 4.09 4.44 1.25 11.88 8.13 4.50 4.34 2.20 2.00 8.55 1.41 0.90 3.73 3.95 9.62 8.25 4.14 6.24 6.35 0.62 7.84 6.12 10.33 3.41 6.05 4.54 7.51 8.87 3.18 2.88 10.55 2.55 8.05 1.44 1.84 5.97 3.31 4.90 1.64 2.81 8.97 4.06 1.78 2.79 3.45 5.85 6.04 4.40 9.47 3.13 2.08 0.12 6.99 4.78 2.10

(continued on next page) 29

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

Table A1 (continued) Automated (LiDAR)

Compass survey

VRGS

deviation angles [°]

Dip dir. [°]

Dip [°]

n

Dip dir. [°]

Dip [°]

Con.

Dip dir. [°]

Dip [°]

αAC

164.86 289.58 289.14 153.74 110.21 116.07 113.95 179.4 130.35 154.29 317.6 178.17 78.25 129.1 154.66 152.02 196.85 156.41 300.61 294.52 114.67 161.72 101.93 166.14 82.61 93.71 84.45 117.95 140.2 137.74 112.08

84.27 87.8 86.13 83.82 84.64 86.88 74.71 80.22 83.87 85.75 81.86 24.84 80.55 35.23 14.17 77.07 54.43 26.12 88.7 88.57 68.64 79.18 78.62 82.06 72.03 46.14 47.88 52.11 46.25 31.68 32.9

7 3 6 9 7 14 6 3 2 7 7 14 4 10 7 6 4 7 10 15 10 2 4 4 2 5 6 4 3 10 14

344.7 107.7 117.8 160.9 287.2 113.2 108.9 175.8 129.61 165.7 315.09 180.3 77.14 134.93 151.97 159.77 198.16 162.31 297.4 289.8 108.33 166.11 99.03 168.07 79 92.27 80.88 122.2 132.94 132.5 103.8

89.85 89.67 89.49 89.33 86 89.96 76.57 83.33 78.51 89.88 81.9 27.18 82.95 42.67 17 82.07 63.24 30.35 81.95 78.67 80.36 80.5 81.26 86.26 75 47.03 50.69 52.79 46.38 32.54 33.4

4 2 3 4 3 4 3 2 2 4 4 3 3 4 2 1 4 3 3 2 4 1 2 3 2 4 2 4 3 5 2

166.53 289.52 290.99 154.84 110.23 115.27 113.68 177.72 127.28 152.82 316.77 176.65 79.2 131.32 152.45 150.85 195.04 158.8 300.35 294.62 114.73 157.06 101.29 166.11 81.71 92.16 83.06 121.89 141.23 136.67 112.36

85.81 87.38 87.88 84.92 85.4 88.5 74.09 78.3 83.82 83.15 82.32 24.22 78.49 28.63 14.95 76.91 52.36 26.01 89.18 88.81 68.38 76.95 79.82 82.91 73.71 45.88 48.44 51.9 48.58 31.63 33.11

5.88 3.15 9.70 9.02 9.83 4.21 5.23 4.73 5.41 12.12 2.49 2.52 2.64 8.29 2.92 9.11 8.88 5.06 7.47 10.95 13.21 4.52 3.89 4.62 4.56 1.37 3.90 3.44 5.25 2.91 4.55 Mean deviation 4.967

αAV 2.27 0.42 2.54 1.55 0.76 1.81 0.67 2.53 3.05 2.98 0.94 0.88 2.26 6.70 0.96 1.15 2.53 1.06 0.55 0.26 0.27 5.08 1.35 0.85 1.89 1.15 1.18 3.11 2.45 0.56 0.26 angles [°] 1.785

αCV 4.71 3.47 7.30 7.49 9.12 2.53 5.25 5.38 5.79 14.50 1.72 3.36 4.90 14.19 2.05 10.17 11.19 4.64 7.80 11.21 13.47 9.56 2.65 3.88 2.91 1.15 2.80 0.92 6.49 2.39 4.70 5.215

15. Riquelme AJ, Abelian A, Tomas R, Jaboyedoff M. A new approach for semi-automatic rock mass joints recognition from 3D point clouds. Comput Geosci. 2014;68:38–52. 16. Slob S, van Knapen B, Hack R, Turner K, Kemeny J. Method for automated discontinuity analysis of rock slopes with three-dimensional laser scanning. Transp Res Rec: J Transp Res Board. 2005;1913:187–194. 17. Roncella R, Forlani G. Extraction of planar patches from point clouds to retrieve dip and dip direction of rock discontinuities. Proceedings of Laser Scanning. Enschede, The Netherlands; 2005:162-167. 18. Vöge M, Lato MJ, Diederichs MS. Automated rockmass discontinuity mapping from 3-dimensional surface data. Eng Geol. 2013;164:155–162. 19. Anders K, Hämmerle M, Miernik G, et al. 3D geological outcrop characterization: automatic detection of 3d planes (azimuth and dip) using lidar point clouds. ISPRS Ann Photogramm, Remote Sens Spat Inform Sci. 2016;III-5:105–112. 20. Bechstädt T, Drews TM, Miernik G, Soyk D, Zühlke R, Stinnesbeck W. AuGE Aufschlussanalogstudien und ihre Anwendbarkeit in der geothermischen Exploration : Diagenese von Geothermie-Reservoiren im Oberrheingraben und deren Aufschlussanaloga & 3D Reservoirmodellierung mit terrestrischem Laser-Scanning. 2016; 2016:64. 21. Miernik G, Profe J, Höfle B. et al. Modelling fractured reservoirs from LiDAR derived digital outcrop models (DOMs). In: Proceedings of the 30th IAS Meeting of Sedimentology. Manchester, 2; 2013. 22. Schlumberger. Petrel Fundamentals Course. Houston, TX, USA; 2012. 23. Slob S, Hack H, Feng Q, Röshoff K, Turner A. Fracture mapping using 3D laser scanning techniques. Proceedings of the 11th Congress of the International Society for Rock Mechanics, Lisbon, Portugal, 1; 2007:299-302. 24. Pearce MA, Jones RR, Smith SA, McCaffrey KJ. Quantification of fold curvature and fracturing using terrestrial laser scanning. AAPG Bull. 2011;95(5):771–794. 25. Duan Y, Li X, Maerz N, Otoo J. Automatic 3D facet orientations estimation from LiDAR imaging. Proceedings Engineering Research and Innovation Conference, Atlanta, Georgia, USA, 2011:4–7. 26. Otoo J, Maerz N, Duan Y, Xiaoling L. LiDAR and optical imaging for 3-D fracture orientations. Proceedings Engineering Research and Innovation Conference, Atlanta, Georgia, USA, 8, 2011. 27. Laux D, Henk A. Terrestrial laser scanning and fracture network characterisation–perspectives for a (semi-) automatic analysis of point cloud data from outcrops. Z der Dtsch Ges für Geowiss. 2015;166(1):99–118. 28. Fernández O. Obtaining a best fitting plane through 3D georeferenced data. J Struct Geol. 2005;27(5):855–858. 29. Abellán A, Vilaplana J, Martínez J. Application of a long-range terrestrial laser scanner to a detailed rockfall study at Vall de Núria (Eastern Pyrenees, Spain). Eng Geol. 2006;88(3):136–148. 30. Roncella R, Forlani G, Remondino F. Photogrammetry for geological applications: automatic retrieval of discontinuity orientation in rock slopes. Proceedings of SPIE-IS

References 1. Barton N. Suggested methods for the quantitative description of discontinuities in rock masses. ISRM, Int J Rock Mech Min Sci Geomech Abstr. 1978;15(6). 2. Hodgetts D. Laser scanning and digital outcrop geology in the petroleum industry: a review. Mar Pet Geol. 2013;46:335–354. 3. Fabuel-Perez I, Hodgetts D, Redfern J. A new approach for outcrop characterization and geostatistical analysis of a low-sinuosity fluvial-dominated succession using digital outcrop models: upper Triassic Oukaimeden sandstone formation, central High Atlas, Morocco. AAPG Bull. 2009;93(6):795–827. 4. Rarity F, Van Lanen X, Hodgetts D, et al. LiDAR-Based Digital Outcrops for Sedimentological Analysis: Workflows and Techniques. 387. Geological Society of London Special Publications; 2014:153–183. 5. Assali P, Grussenmeyer P, Villemin T, Pollet N, Viguier F. Surveying and modeling of rock discontinuities by terrestrial laser scanning and photogrammetry: semi-automatic approaches for linear outcrop inspection. J Struct Geol. 2014;66:102–114. 6. Wang X, Zou L, Shen X, Ren Y, Qin Y. A region-growing approach for automatic outcrop fracture extraction from a three-dimensional point cloud. Comput Geosci. 2017;99:100–106. 7. Assali P, Grussenmeyer P, Villemin T, Pollet N, Viguier F. Solid images for geostructural mapping and key block modeling of rock discontinuities. Comput Geosci. 2016;89:21–31. 8. García-Sellés D, Falivene O, Arbués P, Gratacos O, Tavani S, Muñoz JA. Supervised identification and reconstruction of near-planar geological surfaces from terrestrial laser scanning. Comput Geosci. 2011;37(10):1584–1594. 9. Gigli G, Casagli N. Semi-automatic extraction of rock mass structural data from high resolution LIDAR point clouds. Int J Rock Mech Min Sci. 2011;48(2):187–198. 10. Lato MJ, Vöge M. Automated mapping of rock discontinuities in 3D lidar and photogrammetry models. Int J Rock Mech Min Sci. 2012;54:150–158. 11. Olariu MI, Ferguson JF, Aiken CL, Xu X. Outcrop fracture characterization using terrestrial laser scanners: deep-water Jackfork sandstone at Big Rock Quarry, Arkansas. Geosphere. 2008;4(1):247–259. 12. Sturzenegger M, Stead D. Close-range terrestrial digital photogrammetry and terrestrial laser scanning for discontinuity characterization on rock cuts. Eng Geol. 2009;106(3):163–182. 13. Sturzenegger M, Stead D. Quantifying discontinuity orientation and persistence on high mountain rock slopes and large landslides using terrestrial remote sensing techniques. Nat Hazards Earth Syst Sci. 2009;9(2):267–287. 14. Sturzenegger M, Stead D, Elmo D. Terrestrial remote sensing-based estimation of mean trace length, trace intensity and block size/shape. Eng Geol. 2011;119(3):96–111.

30

International Journal of Rock Mechanics and Mining Sciences 109 (2018) 19–31

T. Drews et al.

31.

32.

33. 34. 35. 36. 37. 38. 39. 40.

41. Hodgetts D, Gawthorpe R, Wilson P, Rarity F. Integrating digital and traditional field techniques using Virtual Reality Geological Studio (VRGS). Proceedings of the 69th EAGE Conference and Exhibition Incorporating SPE EUROPEC 2007; 2007. 42. Fisher R. Dispersion on a sphere. Proc R Soc Lond A. 1953;217(1130):295–305. 43. Mardia K. Probability and Mathematical Statistics: Statistics of Directional Data. London: Academy Press; 1972. 44. Best DJ, Fisher NI. The bias of the maximum likelihood estimators of the von misesfisher concentration parameters. Commun Stat-Simul Comput. 1981;10(5):493–502. 45. Pfeifer N, Mandlburger G, Otepka J, Karel W. OPALS–a framework for Airborne Laser Scanning data analysis. Comput, Environ Urban Syst. 2014;45:125–136. 46. Weinmann M, Jutzi B, Mallet C. Semantic 3D scene interpretation: a framework combining optimal neighborhood size selection with relevant features. ISPRS Ann Photogramm, Remote Sens Spat Inform Sci. 2014;2(3):181–188. 47. Watson G. Analysis of dispersion on a sphere. Geophys Suppl the Mon Not R Astron Soc. 1956;7(4):153–159. 48. Dershowitz WS, Herda HH. Interpretation of fracture spacing and intensity. Proceedings of the 33th US Symposium on Rock Mechanics (USRMS). American Rock Mechanics Association; 1992. 49. Narr W, Suppe J. Joint spacing in sedimentary rocks. J Struct Geol. 1991;13(9):1037–1048. 50. Terzaghi RD. Sources of error in joint surveys. Geotechnique. 1965;15(3):287–304. 51. Priest S. Discontinuity Analysis for Rock Engineering. New York: Chapman & Hall; 1993.

&T Electronic Imaging. 5665. San Jose, CA, USA: International Society for Optics and Photonics; 2005:17–27. Fischler MA, Bolles RC. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. Commun ACM. 1981;24(6):381–395. Jaboyedoff M, Metzger R, Oppikofer T, et al. New insight techniques to analyze rockslope relief using DEM and 3D-imaging cloud points: COLTOP-3D software. Proceedings of the Rock Mechanics: Meeting Society’s Challenges and Demands. 1; 2007:61–68. Riquelme AJ, Abelian A, Tomas R. Discontinuity spacing analysis in rock masses using 3D point clouds. Eng Geol. 2015;195:185–195. Bogdanov A. The intensity of cleavage as related to the thickness of beds. Sov Geol. 1947;16(000). Ladeira F, Price N. Relationship between fracture spacing and bed thickness. J Struct Geol. 1981;3(2):179–183. Ortega OJ, Marrett RA, Laubach SE. A scale-independent approach to fracture intensity and average spacing measurement. AAPG Bull. 2006;90(2):193–208. Fisher NI, Lewis T, Embleton BJ. Statistical Analysis of Spherical Data. Cambridge: Cambridge University Press; 1987. Hall GG. Matrices and Tensors. Oxford: Pergamon Press; 1963. Scheidegger A. On the statistics of the orientation of bedding planes, grain axes, and similar sedimentological data. US Geol Surv Prof Pap. 1965;525:164–167. Woodcock N. Specification of fabric shapes using an eigenvalue method. Geol Soc Am Bull. 1977;88(9):1231–1236.

31