The use of an unmanned aerial vehicle for fracture mapping within a marble quarry (Carrara, Italy): photogrammetry and discrete fracture network modelling

ABSTRACT This paper describes the use of a drone in collecting data for mapping discontinuities within a marble quarry. A topographic survey was carried out in order to guarantee high spatial accuracy in the exterior orientation of images. Photos were taken close to the slopes and at different angles, depending on the orientation of the quarry walls. This approach was used to overcome the problem of shadow areas and to obtain detailed information on any feature desired. Dense three-dimensional (3D) point clouds obtained through image processing were used to rebuild the quarry geometry. Discontinuities were then mapped deterministically in detail. Joint attitude interpretation was not always possible due to the regular shape of the cut walls; for every discontinuity set we therefore also mapped the uncertainty. This, together with additional fracture characteristics, was used to build 3D discrete fracture network models. Preliminary results reveal the advantage of modern photogrammetric systems in producing detailed orthophotos; the latter allow accurate mapping in areas difficult to access (one of the main limitations of traditional techniques). The results highlight the benefits of integrating photogrammetric data with those collected through classical methods: the resulting knowledge of the site is crucially important in instability analyses involving numerical modelling.


Introduction
The Carrara district in Italy is one of the most productive marble quarrying regions in Europe and represents an important economic resource. Many people work in the quarries, and safety in the workplace is an overriding priority. In this context, the use of new technologies such as unmanned aerial vehicle (UAV) systems to study and monitor mining and quarrying areas can be very important in risk assessment and management. Over the last decade, geomatic techniques such as terrestrial laser scanning (TLS) and digital terrestrial photogrammetry (DTP) have been used increasingly to characterize rock masses throughout the world (Abell an et al. 2006;Coggan et al. 2007;Ghirotti & Genevois 2007;Ferrero et al. 2009;Lato et al. 2009;Sturzenegger & Stead 2009a, 2009bSalvini et al. 2011;Jaboyedoff et al. 2012;Salvini et al. 2013;Francioni et al. 2015). One of the most common problems encountered in the use of these techniques, especially in the case of slopes with complex geometries, is the presence of occlusions. This study demonstrates how the use of UAV systems can overcome such limitations by acquiring data at different angles, depending on the orientation of walls. Photographs taken by UAV can be used to map all visible discontinuities irrespective of their position, attitude or elevation from the ground. Quarry walls in the Carrara marble district can be up to hundreds of metres high and are dominated by natural slopes with very complex morphologies. In these quarries, conventional structural and engineering geological survey data can be collected only at the foot of slopes or by climbing on high rock walls; these data often provide an incomplete understanding of the area.
The use of UAV can significantly improve the collection of geometric and structural data in such areas. The latest drone hardware and software developments allow their application in several fields of science, including the study of landslides and rockfall instability (Rau et al. 2011;Niethammer et al. 2012;Danzi et al. 2013;Salvini et al. 2014;Francioni et al. 2015;Giordan et al. 2015;Turner et al. 2015). Digital surface models (DSMs), orthophotos and three-dimensional (3D) models derived from processed UAV images allow accurate mapping of joint geometries and their degree of uncertainty. Such data can then be used both for the calculation of yield (in terms of commercial value) and for the deterministic analysis of stability and the study of stress distribution within slopes (in terms of safety). For these purposes, a case study was completed in the Lorano open pit (Apuan Alps, Italy). In particular, data from a UAV photogrammetric survey were used to map fractures in a buttress-shaped remnant of previous excavation activities. Data were integrated with information from traditional engineering geological surveys to derive important parameters such as fracture orientation, intensity, density, size and termination. This enabled the creation of a 3D discrete fracture network (DFN) model built on the basis of geomechanical data from traditional field surveys, the literature, laboratory tests, UAV image processing and geological photointerpretation. The DFN approach allows for an excellent use of collected fracture data and can be used to obtain a twodimensional (2D) or 3D synthetic rock mass (SRM) model. DFN is a stochastic representation of discontinuities derived from the analysis of fracture density and intensity parameters that are generally obtained by mapping line traces (discontinuities) on 2D or 3D rock exposures. The fracture pattern, which takes into account the natural variability of the rock mass, is therefore based on statistical inputs. The various stochastic models are expected to differ; consequently, several models are required to understand simulation outcomes (Weir & Fowler 2014).
The DFN approach is quite new in the field of engineering geology. Before 2010 there was almost no mention of DFNs in the proceedings of major conferences on rock mass stability analysis (Lorig 2014). However, in recent years the DFN approach has gained increased usage in numerical modelling for rock engineering purposes (Bakun-Mazor et al. 2009;Scholt es & Donz e 2012;Bahrani et al. 2014;Grenon et al. 2014;Lei et al. 2014;Borghi et al. 2015). This is because the method overcomes one of the principal limitations of traditional analyses: the unrealistic assumption of ubiquitous joints of infinite length (Weir & Fowler 2014). Ideally, each joint in the rock mass should be characterized; however, the limited visibility of discontinuities (outcrop observation) precludes such detailed study (Dershowitz & Einstein 1988), even using the most recent geomatic and geophysical techniques. DFN modelling is therefore the best available option. In this work, the DFN model was built to maximize the use of fracture data collected through exposure mapping and to build an SRM model. Since the accuracy of the DFN model depends directly on the quality and quantity of available data, high-definition images and 3D information derived from UAV data can make a strong contribution. In addition, both certain and uncertain joint attitudes were discerned in this work. This was done because on a flat surface like a quarry wall face it is difficult to correctly interpret the set to which a discontinuity belongs, an important point requiring further studies in the context of DFN model calibration and validation.

Geographic and geological setting
The Lorano open pit is located in the Province of Massa and Carrara, north-western Tuscany (Italy). In 1997 a rockfall in the quarry interrupted excavation activities for a few weeks. Remediation works were subsequently carried out with the aim of ensuring safety, and a marble buttress accessible from three sides was emplaced (Figure 1). Due to the ongoing quarrying activity the buttress is now about 150 m high, 30 m wide and 40 m thick.
The quarry is located in a fold-and-thrust belt of the Northern Apennines, which formed through the Tertiary collision (66 Ma) between the Sardinia-Corsica block and the Adria plate (Boccaletti et al. 1971;Scandone 1979;Dercourt et al. 1986). The study area belongs to the Apuan Alps metamorphic complex (Figure 2), first described by Zaccagna (1932); deep levels of the belt are exposed in this, the largest tectonic window in the inner Northern Apennines (Elter 1975;Carmignani & Kligfield 1990;Molli 2008). The metamorphic complex consists of two main tectono-metamorphic units, the 'Apuan unit' and the overlying 'Massa unit.' According to the classical interpretation, the regional tectonic setting of the Apuan Alps is the result of two main tectono-metamorphic events, D1 and D2 (Carmignani & Kligfield 1990). D1, the ductile compressional event, was determined by the Tertiary continental collision between the Sardinia-Corsica block and the Adria plate; it was followed by the D2 extensional event that led to isostatic re-equilibration (Carmignani & Kligfield 1990). During the D1 event, tectonic units of the Tuscan and Ligurian domains were stacked, and deformation occurred progressively in two stages (Molli & Meccheri 2000). The greenschist foliation (Sp) in the axial plane of isoclinal micro-to kilometre-scale folds represents the main stage. During the D2 event, the previously formed structures were deformed and developed different generations of folds and local high-strain zones associated with exhumation and vertical shearing (Molli 2012). The late stages of D2 are characterized by the development of brittle structures (low-and high-angle faults and joint systems) associated with the final exhumation and uplift of the metamorphic units in the context of the late-to post-orogenic regional extension of the inner Northern Apennines (Molli et al. 2010).
According to Carmignani et al. (2002), the brittle stages of the D2 phase led to the development of the three main discontinuity systems that characterize the marble district. The first system (J1), shows an average anti-Apennine direction ranging from N 20 À30 E to N 80 À90 E and a general vertical or subvertical inclination (dipping up to 50 À60 both to the NW and the SE). The second system (J2) ranges in direction from N 120 E to N 150 E with a mediumÀhigh dip generally to the E and the NE. The third and last system (J3) is often pervasive and almost parallel to the foliation (Sp); it shows a direction similar to J2 and dips steeply to the SW.
The Lorano marble quarry is located in the normal limb of the 'Pianza anticline' which, together with the 'Vallini syncline,' represents a coupled antiformÀsynform with a core of Jurassic marbles and cherty meta-limestone. These are minor folds (hectometre-scale) between two D1 structures known as the 'Carrara syncline' and 'Vinca anticline' located to the NW and the SE, respectively (Molli & Meccheri 2012). Most of the quarried marble belongs to the white marble group, characterized by homogeneous marbles of mediumÀfine grain size (about 100À200 mm), and ordinary marble (Meccheri 1996), with a medium grain size (about 200 mm). Two subordinate types of marble in the quarry are veined grey marble and breached marble (Carmignani et al. 2007).

Engineering geological survey
An engineering geological survey was carried out to characterize the geomechanical properties of discontinuities. Data were collected in accessible areas using traditional scan-line mapping techniques whereby all the discontinuities that intersect a defined scan line are recorded. About 100 discontinuities longer than 10 m were identified in seven scan lines. Collected data were compared with those of Profeti and Cella (2010) for the same area and integrated with DTP and TLS data from inaccessible zones, as described in Francioni et al. (2015). The attitude of 301 joints was calculated from TLS and an additional 236 discontinuities were manually measured from the photointerpretation of UAV imagery . Data processing revealed the following four sets of discontinuities (named K1, K2, K3 and S1 in accordance with Profeti & Cella 2010) describing the current state of the buttress ( Figure 3): S1 (foliation), SW dipping with average dip of about 50 , corresponding to Sp of Carmignani et al. (2002); K1, SE dipping, subvertical, corresponding to J1 of Carmignani et al. (2002); K2, NE dipping with average dip of about 50 , corresponding to J2 of Carmignani et al. (2002); and K3, SW dipping, subvertical, corresponding to J3 of Carmignani et al. (2002).
The K1, K2 and K3 systems are characterized by metre spacing, millimetre to centimetre apertures, moderate roughness and no infill. According to the rock mass rating (RMR, Bieniawski 1989) the rock mass is of good quality (basic RMR or RMRb D 76). Geological strength index (GSI, Hoek & Brown 1997) values within engineered areas (i.e. open pit) are variable and differ with respect to those of the natural rock outcrops (i.e. the slope above). Considering the characteristics of the marble, the GSI was estimated to vary between 50 and 60 on the natural slope and between 70 and 80 on the face of the open pit. This variation is probably due to the degree of weathering of the natural slope surface and to the fact that the joint spacing usually increases with depth .

UAV survey
Due to the considerable height of the marble buttress, a photogrammetric survey was carried out using UAV technology to obtain high-resolution digital products (DSMs and orthophotos). Given the morphology of the area, a multi-rotor Falcon 8 UAV vehicle (Asctec TM Gmbh) equipped with a Sony TM NEX-5N digital camera was used (Table 1). Falcon 8 consists of eight electric rotors, a remote control system (mobile ground station) and software for flight plan management. The UAV is equipped with a GPS and an inertial navigation system that can record 3D spatial coordinates and the orientation of the camera at every shot (description of the on-board sensors is reported in Table 2).  Five vertical flights with 80% overlap and 40% sidelap ensured coverage of the entire area; the photos were taken from an average distance of 30 m from the slope surface (Figure 4), yielding an estimated Ground Sample Distance of about 1 cm.
With the aim of improving the external orientation and level of accuracy of the entire photogrammetric model, a topographic survey was carried out using a Leica TM TCRP 1203 C R1000 reflectorless total station (TS) and two geodetic Leica TM 1200 GPS receivers. Due to the complex morphology of the buttress and the large area, a large number of ground control points (GCPs) were used to orient the photogrammetric model. The TS was used to acquire 54 evenly spaced GCPs within the marble buttress. Two geodetic devices operating in static mode were used to obtain the geographic coordinates of two points: the origin of the survey and its zero-Azimuth direction. The GPS data acquired were corrected using contemporary data recorded by two permanent GPS stations located in Lucca and Borgo a Mozzano .
Specific post-processing was used to determine the absolute coordinates of all GCPs. The 391 photos acquired through the UAV survey were processed using the Agisoft TM PhotoScan Professional software, version 1.2.3 (Agisoft 2016). This software, based on specific feature extraction algorithms, was used to solve the camera's interior and external orientation parameters and extract georeferenced products such as 3D point clouds, DSMs, orthophotos and textured models. The first task carried out in PhotoScan was alignment processing. This was used to solve the interior and relative orientation parameters and to extract millions of tie points. Although the software allows for a tie point extraction limit, this parameter was set to 0 (no limit) in order to improve alignment and  obtain a lower reprojection error. At the end of the process PhotoScan was able to align all the images and extract more than 2 million tie points with an average reprojection error of about 0.2 pixels. Twenty-nine GCPs measured with the TS during the topographic survey were used as control points for georeferencing the model and solving the exterior orientation parameters. The remaining 25 GCPs were used as check points to assess the accuracy of the photogrammetric model. This operation was carried out by identifying the features measured through the topographic survey directly on the photos. Each feature was also associated with its relative 3D geographic coordinates.
PhotoScan's optimize tool (Agisoft 2016) was used to adjust the estimated camera positions and to remove possible non-linear deformations, minimizing the errors due to reprojection and misalignment of the photos. At the end of the process the estimated root mean square error (RMSE) of GCPs was 0.029 m, and that of the check points was 0.033 m.
In order to assess the reliability of the obtained RMSEs, the selection and number of control and check points on the buttress were varied. Table 3 shows the results of three tests completed by varying the number of control and check points. Considering 54 GCPs, when the number of control and check points is varied from 29 to 16 and from 25 to 38, respectively, the resulting RMSE does not change (only 1 mm for check points). In order to validate the RMSE achieved with 29 control points, the accuracy of 10 check points was determined: 0.025 m (Table 3). From this it follows that the availability of a large number of spatially well-distributed GCPs allows accuracy assessment thanks to the numerous check points and the correct orientation of images. Figure 5 shows the spatial distribution and error of GCPs in the final configuration (29 control points and 25 check points).  Once the resulting RMSEs were checked and accepted, a dense cloud consisting of more than 20 million points was created (using the settings Medium quality and Aggressive depth filtering). This dense cloud was used to create a mesh of about 4 million facets. Lastly, the 3D model was used to create the orthophotos of the buttress by removing image distortion due to camera characteristics (i.e. lens distortions), camera tilt and topographic relief displacement. Unlike an uncorrected aerial photograph with a perspective projection, an orthophoto is geometrically corrected ('orthorectified') and can be used to measure true distances since it is 'scale-corrected.' Four different planes of reference, one for each side of the buttress (eastern, southern and western) and an overhead view, were selected for measurements. The resulting orthophotos with a spatial resolution of 1 cm/pixel were exported in the local coordinate system; only the top view was projected into the Italian National GaussÀBoaga system (Figure 6).

Fracture mapping and DFN modelling
Different software allows DFN model simulations by defining the statistical distribution of joint length and intensity. In this case study the FracSim3D software (Xu & Dowd 2010) was used to create the 3D DFN model. This software can simulate 2D and 3D stochastic rock fracture networks using marked point processes. The geometric characteristics supported by this software are fracture location, size (diameter) and orientation. These properties can be defined using different processes and distribution functions (Table 4).
In this case the homogeneous model was adopted for fracture locations, whereas exponential and Fisher distributions were used for the size and orientation of the discontinuities, respectively. The first step is to simulate the fracture location by generating a number of points. For any joint set, the number of points falling inside a certain area (A), N(A), follows a Poisson distribution with mean g Ã v(A), where v(A) is the volumetric measure of A and g is the intensity of the distribution. The second step is to simulate the geometry of each individual fracture. A stochastic network of discontinuities can be generated in this way.
The 'P' system, introduced by Dershowitz and Herda (1992), was adopted in this case study to define discontinuity density (P 10 and P 20 ) and intensity (P 21 ). It is a common method for defining fracture density/intensity parameters starting from line traces mapping on 2D survey windows, and one that provides a convenient framework for moving between differing scales and dimensions. The mentioned parameters can be calculated as follows: The method was chosen for its simplicity and for the availability of high spatial resolution orthophotos and 3D information on the buttress. Three high-resolution orthophotos (eastern, southern and western sides of the buttress) were used to derive the 'P' system parameters. In particular, a detailed cartographic approach within a geographic information system (GIS) was adopted. Six different surface mapping windows, two for each side of the marble buttress, were drawn on the orthophotos using the ESRI TM ArcMap software. The survey windows, each measuring 10 m £ 15 m, were spaced a few metres apart in the central parts of the buttress. All discontinuities observed within the windows were then measured in order to calculate the P 10 , P 20 and P 30 values. In this case the flatness of the quarry wall face represented in a 2D orthophoto complicates things because it makes it difficult to assign a discontinuity to the correct set identified through the engineering geological and geomatic surveys (only line traces are visible). The apparent inclination of every discontinuity system on each side of the buttress was therefore calculated. Because the buttress faces are irregular due to quarrying, the difficulties relating to fracture characterization partially persisted. With the aim of facilitating this operation, the top view orthophoto was used to observe the quarry benches from above (so as to 'see' a second orthogonal surface). The 3D model formed by the mesh built in PhotoScan from the dense cloud and the orthophotos of the buttress were used within the ESRI TM ArcGIS Pro software to have an additional control of joint attitude in 3D visualization. No additional fractures were mapped in the 3D view since the rock fracture network was simulated using a stochastic, not a deterministic approach.
By combining 2D and 3D views with data collected during in situ surveys, discontinuities on all three sides of the marble buttress were identified with a high degree of certainty. Despite efforts to minimize uncertainty, discontinuities with uncertain attitude were also mapped because some highly dipping sets (e.g. K1 and K3) produce similar line traces on the buttress surface. The deterministic maps of discontinuities (with length greater than 2 m) are shown in Figures 7 and 8.
After identifying discontinuities and using the P system, it was possible to measure and collect all data required to produce several 3D DFN models using the FracSim3D code. In this context, P 10 , P 20 and P 21 values were calculated using the six windows located on all three sides of the buttress. In particular, P values were calculated for both the certain (P certain , including only joints assigned with certainty to a specific joint set; see Figures 7 and 8) and uncertain attitudes (P certainCuncertain , including also those joints of uncertain attitude; see Figures 7 and 8) of joint systems. Tables 5À7 show the results in terms of the calculated P values.
In addition, information about the mean (M) and standard deviation (SD) values of P as computed for every discontinuity system in all the mapped windows are presented in Table 8. The values calculated on different survey windows show good correlation and, as mentioned earlier, can be used to generate a number of DFN models by varying the distribution functions that characterize the discontinuity sets. Figure 9 shows an example of a DFN model generated using FracSim3D. In this case certain P 10 values obtained on windows 3 and 4 were used to define joint Figure 8. Survey windows 3À6 used to map discontinuity sets K1, K2 and K3, as well as foliation S1, on the southern and eastern sides of the buttress. density, a Fisher distribution was used for the orientation of discontinuities and an exponential distribution was used to assess discontinuity lengths. The software can also be used to associate additional geomechanical properties to the discontinuities using a uniform, normal, exponential, lognormal or non-parametric distribution.  Figure 9. Generation of a DFN model for the marble buttress using FracSim3D. K1, K2 and K3 joint sets included.
The S1 foliation was not included in the shown model because it does not represent a real plane of discontinuity in the studied quarry. Nevertheless, it is a plane of weakness that could affect the stability of the buttress and will be the subject of future research.
Several models can be created to assess the influence of uncertainty on the SRM model and improve the reliability of results. This is particularly important because several trial runs can be completed and analyzed in order to obtain reliable results for a stochastic representation of a natural discontinuity network. For example, Figures 10 and 11 show 2D representations (including stereographic projection and normal probability plots for fracture lengths) of two DFN models generated using respectively certain and certainCuncertain P 10 values (as average among all the six mapped windows).

Discussion
The adopted approach, based on the combined use of high-resolution UAV images and engineering geological data on a marble buttress, allowed the construction of a reliable 3D DFN rock mass model. Thanks to the acquisition of images from different angles it was possible to overcome data acquisition difficulties for high steep quarry walls and to create high-resolution orthophotos from different perspectives. GCPs measured using a TS and GPS receivers allowed for a high level of accuracy in the external orientation of photographs, which is particularly important for joint measurements. The methodology was used to obtain accurate deterministic measurements of discontinuity intensity and density for every joint set in the marble buttress. This was particularly important because the density and intensity of discontinuities can vary at different heights of the rock mass due to stress relaxation induced by excavation activity. An in-depth examination of the rock mass revealed differences in the discontinuities on different zones of the buttress. For example, the bottom of the buttress shows fewer discontinuities, probably because it has been excavated more recently and brittle fracturing associated with stress release has not yet propagated through the intact rock. The upper portions of the slope, instead, are more fractured. This difference is particularly evident in the S1 system (foliation): S1 discontinuities are numerous in the upper parts of the slope, whereas they are hard to find at the bottom of the buttress. In the lower buttress S1 is present only as a plane of weakness that may later evolve into a real discontinuity. This difference is also highlighted by P 10 , P 20 and P 21 results. Considering the buttress front, which is perpendicular to the supposed direction of relaxation of the structure, the density/intensity values of window 3 (upper window, 0.467) are substantially higher than those of window 4 (lower window, 0.133). Such information is very useful in creating reliable DFN models to be calibrated for stability analysis purposes. Therefore, the possibility to obtain information at different heights can allow the generation of realistic DFN models achieved by creating sub-domains where different P values can be used. Areas with high concentrations of S1 discontinuities are currently those most likely to experience problems due to stress relaxation and possibly failure. In future numerical modelling stability analyses it will be important to compare the obtained stress values with results from in situ stress measurements (already planned in collaboration with the Massa and Carrara USL (local health authority), Mining Engineering Operative Unit À Department of Prevention Hygiene Safety Workplace). This will allow better assessment of the possible effects of density/intensity values on DFN modelling and the stability of the buttress.
In addition, both certain and uncertain joint attitudes were discerned in this work. This was done because on a flat surface like a quarry wall face it is difficult to correctly interpret the set to which a discontinuity belongs, an important point requiring further studies in the context of DFN model calibration and validation. Although the geometric parameters required to build a DFN (fracture location, size, shape, orientation, density, flow property and number of sets) can all be derived deterministically through photointerpretation, additional statistical input is important for generating the discontinuity sets. The presence of joints with an uncertain attitude results in a model with a higher number of discontinuities showed, as an example, by comparing Figures 10 and 11. In fact, the DFN model obtained using P certainCuncertain values, resulted in a 32% increase of the number of Figure 11. Generation of a DFN model for the marble buttress using FracSim3D with P 10 certainCuncertain values. K1, K2 and K3 joint sets included. The location of the defined plane is the same as that of Figure 9. discontinuities, from 306 to 404 joints. On the other hand, uncertain discontinuities may not be representative of specific sets but of random fractures. The key aspect in this case is that the amount of data collected thanks to the adopted UAV-based approach can possibly lead to the development of a large number of different DFN models to be used for stability analysis purposes, which are more likely to satisfy the statistical criteria.
In this case study the particular shape of the buttress was an advantage: as the buttress was accessible from three sides, it was possible to create a complete 3D model of the rock mass. By measuring the joints from different perspectives, orientation bias was minimized. Moreover, it was possible to check whether the 3D DFN model obtained from the observation of 2D rock exposures was representative of the actual buttress condition. This is well illustrated by observing the obtained DFN stereographic projections (Figures 9(b), 10(a) and 11(a)) and the 2D sections of Figures 9(d), 10(c) and 11(c) which were created from the 3D model by a defined plane always in the same position. Results are respectively very similar to the stereographic projection from the engineering geological survey ( Figure 3) and to the fracture traces observable on the deterministically mapped windows (Figures 7  and 8). In this regard, useful information can also be obtained by analyzing the fracture trace lengths. Francioni et al. (2015) report that fracture trace lengths on the buttress range from 10 to 30 m for K1, K2 and K3. This work well correlates with these findings, as illustrated in the two normal probability plots (Figures 10(b) and 11(b)), that show how the majority of the fracture traces have lengths ranging between 10 and 35 m.
This validation of the model confirmed the quality of results and the great potential of UAV systems in detailed surveys, even in inaccessible areas. Conversely, possible limitations in the use of UAV can be linked to the need of a UAV pilot license and user experience on topographic survey. In fact, the accuracy of the output can be affected by the quality of data collected (photos and GCPs), hardware and software capability and, last but not least, user expertise. When performing manual photointerpretation, the precision and accuracy of measurements can vary, depending on the skill of the interpreter. Although several studies have demonstrated the use and reliability of automatic and semi-automatic processing of imagery and 3D point clouds for fracture mapping (Mah et al. 2011;V€ oge et al. 2013;Assali et al. 2014;Vasuki et al. 2014), a fully manual approach was here adopted to guarantee consistency and quality in the interpretation of discontinuities. Note that in most cases the flat and regular morphology of quarry walls only allowed photointerpretation of discontinuity traces. Final visual inspection of outputs is always required, even when using codes for the semi-automatic extraction of joints.

Conclusions
One of the main advantages of using a UAV for rock slope stability studies is the ability to acquire data in inaccessible areas where foliation and fractures may vary their attitude and other fundamental geometric characteristics. Flight plans can be customized to reach any height above the ground so that observations are optimized and occlusions are avoided. UAV data collection is non-invasive, safe and inexpensive compared to TLS and helicopter surveys (no crew required). UAV surveys can be completed quickly. The produced DSM from 3D point clouds and orthophotos are easily interpreted and can be managed in a GIS environment. A large number of features can be accurately mapped both in 2D and 3D, with great flexibility in data editing.
This study shows how data from UAV image photointerpretation can be successfully used to study the geological structure of a rock mass and build a 3D DFN model of its fracture system. The importance of mapping uncertain joints has been highlighted, as well as the benefits of accurately measuring fracture density/intensity in different portions of the rock mass. Nevertheless, model calibration and validation using either in situ stress measurement or displacement data are required to obtain reliable results.
Other advantages of using UAV systems in the study of rock slope stability are linked to the possible acquisition of multi-temporal data, and the ability to update data and customize the nominal scale of photo shoots. The study of rock mass fracturing can greatly benefit from this type of approach in determining both yield (in terms of commercial value) and stability (in terms of safety). Within this context, GIS environments allow for the simple, comprehensive management, processing and analysis of georeferenced data.