Brief communication: Remotely piloted aircraft systems for rapid emergency response: road exposure to rockfall in Villanova di Accumoli (central Italy)

Abstract

The use of remotely piloted aircraft systems (RPASs) in geosciences is often aimed at the acquisition of an image sequence to produce digital models and orthophotographs of the topographic surface. The technology can be applied for rockfall hazard and risk assessment. To study rockfalls, an approach consists in the application of numerical models for the computation of rockfall trajectories. Data required for such simulations include digital terrain models, location of the instability source areas, and the mechanical properties of the terrain. In this article, we present an analysis of the earthquake-triggered rockfall that occurred along the SP18 in Villanova di Accumoli (Lazio, central Italy) during the seismic sequence that started on 24 August 2016. A survey with a multicopter was carried out to obtain a surface model of the terrain and identify and characterize the source areas and other instable blocks in areas not accessible in the field. The investigated area extends for 6500 m${}^{\mathrm{2}}$ and was covered by 161 photographs that were used to obtain an orthophoto with a ground resolution of 2.5 cm and a digital surface model with a ground resolution of 20 cm $\times$ 20 cm, which was processed and fused with GNSS real-time kinematic data. To obtain a map of potential rockfall trajectories, we run the numerical model STONE, using as origin of the boulders both source areas mapped in the field and pixels with a slope angle above a selected threshold. Results showed that only the part of the road SP18 already affected by the rockfall was exposed to further rockfall impacts. In particular, it was observed that 29.2 % (i.e. 12 123) of the 41 500 simulated trajectories may potentially reach or cross this tract of the road. Based on these data, limited protection measures were suggested. The combined use of RPAS data, fused with ground GPS points, an accurate geomorphological survey, and terrain static and dynamic parameters from the literature allows fast, low-cost, and replicable rockfall numerical modelling useful for emergency response and adoption of proper protection measures.

Introduction

Rockfall is a widespread natural hazard that poses continuous risk to the population in mountain areas worldwide . Due to the prevalent type of movement (free falling, bouncing, rolling, sliding), rockfalls are among the fastest and deadliest landslide type . Rockfalls can be triggered by earthquakes, intense rainfall, frost weathering, wind, root growth, and even traffic . The elements most exposed to rockfall hazard are transport corridors , which often cross hazardous areas. In particular, rockfalls cause relevant damage to structures and infrastructures along secondary and minor transport networks, where adequate protection measures are not economically sustainable , and pose a severe risk to people.

During triggering events such as seismic sequences, areas exposed to rockfall hazard can be hit in several places , causing much damage and multiple interruptions of the infrastructure, particularly roads and rails. In such events, it is important to characterize the most critical situations along roads and rails to guarantee, in short time spans, the safe use of the transportation infrastructure. Operating workflows to better characterize rockfall source areas and trajectories could lead to identifying areas where rockfalls can interact with roads, allowing for tailored protection measures, and developing more rapid infrastructure protection activities.

Defining the potential interaction of rockfalls with elements at risk, knowing their geographic location, is a challenging task. Empirical, process-based, and GIS-based software including STONE , HY-STONE , CONEFALL , and Rockyfor3D has been implemented and used to model the spatial pattern of rockfall trajectories. All the modelling approaches are based on at least a given distribution of rockfall source areas and a digital elevation model (DEM). Detection and mapping of unstable rock masses can be accomplished by field activities, interpretation of high- and very high-resolution images and DEMs, or using morphometric criteria . Depending on the scale of the modelling, elevation data can be coarse to very fine. For example, for regional-scale modelling, a ground sampling distance (GSD) of 25 m (Guzzetti, 2003) (http://damocles.irpi.cnr.it/docs/final_reports/Perugia-Detailed-Report-Feb2001-Mar2003.pdf, last access: 25 January 2019) or 10 m is adequate , whereas very high-resolution elevation data (GSD of 1 m or less) are suitable for slope-scale models.

High- and very high-resolution elevation data can be successfully acquired using cameras mounted on remotely piloted aircraft systems (RPASs). The acquired high(ultra)-resolution images with different orientations of the camera (nadir, oblique) can be used to obtain elevation data through the application of structure from motion (SfM) algorithms . SfM is a photogrammetric-range imaging technique capable of finding homologue points in a sequence of images acquired with a high overlap and transforming the homologue points in a georeferenced cloud of elevation points . The accuracy of the positioning of SfM products can be improved using ground control points (GCPs) whose position is measured using GNSS (Global Navigation Satellite System) . The main products of RPAS photogrammetry that can support the study of rockfalls are digital surface models (DSMs) and true orthophotographs .

High-resolution elevation data can also be collected by aerial and terrestrial laser scanners, which acquire a set of measures of distance between the laser scanner and the ground surface . A terrestrial laser scanner (TLS) is often adopted to study rocky slopes and other types of landslides . Use of a TLS often requires the integration of multiple surveys taken by different points of measure to ensure the coverage of complex morphologies. Arguably, the application of such a technique is limited by an overall benefit in terms of quality and time compared to other techniques. In particular, the local morphological setting could require multiple points of acquisition, which would imply a longer acquisition and processing time. An aerial laser scanner (ALS) has shown good performances in densely vegetated areas where photogrammetric techniques are more limited in obtaining ground elevation data. Conversely, ALS is not effective for the characterization of sub-vertical slopes, and rock cliffs, where the acquisition of oblique photo-sequences is the best solution for a more detailed representation of rockfall source areas . Use of RPASs for SfM is a solution that requires limited costs, short mission planning, and a brief setup time compared to ALS, particularly when the study area is small (up to a few square kilometres). Such advantages can be important, particularly where a multi-temporal approach has to be adopted to detect the evolution of the studied area or during emergencies, when support to decision makers for the identification and planning of first responses must be provided quickly . The capability of RPASs to provide very high-resolution georeferenced images in a short time can be effective to support fast emergency responses in areas affected by natural disasters . The major limitation of RPASs is the extension of the area that can be surveyed in a single flight .

According to , two approaches can be adopted for the use of RPASs depending on the slope steepness and the target resolution. For steep slopes, the use of multirotor RPASs, which allow us to acquire oblique photographic-sequences, is recommended. Conversely, on gentle slopes nadir acquisition carried out by fixed-wing RPASs can be considered the best solution to obtain a high-resolution coverage of large areas. In the study of areas affected by rockfall, a hybrid approach can be considered for the acquisition of a DSM that can be used for the definition of rock block trajectories and the identification of source areas. RPASs guarantee the possibility of also acquiring high-resolution DEMs in inaccessible areas . Several examples of the use of RPASs for the acquisition of oblique and nadiral image sequences for the identification of rock slope instabilities have been published in recent years.

Despite such an increasing attention to applications of RPAS images and photogrammetry to landslide risk assessment, only a limited number of articles focuses on applications during emergencies and on the assessment of the extent of damage caused by natural disasters . Furthermore, none of them, to the best of our knowledge, focus on testing a procedure that guarantees semi-quantitative information in a relatively short time to provide an evaluation of the residual rockfall risk during emergencies, when there are time and budget constraints.

In 2016, central Italy was affected by a very long and severe seismic sequence that began on 24 August with a $M_{\mathrm{w}}$ 6.0 earthquake, followed by a second main shock on 26 October ($M_{\mathrm{w}}$ 5.9), and third on 30 October ($M_{\mathrm{w}}$ 6.5). The seismic sequence, characterized by more than 50 000 aftershocks in 4 months, triggered numerous rockfalls that caused damage to structures and infrastructures. In this paper, we describe the use of RPASs to study an earthquake-triggered rockfall in the vicinity of the Villanova di Accumoli, Rieti, central Italy. The provincial road SP18 near Villanova di Accumoli was closed after a 1 m${}^{\mathrm{3}}$ boulder fell from a rock cliff and crossed the road (Fig. 1). During the seismic emergency, the Italian National Department of Civil Protection was supported by the Research Institute for Geo-Hydrological Protection of the Italian National Research Council (CNR-IRPI) in the definition of the conditions to safely reopen the road. To answer the request, a numerical model was applied to evaluate possible rockfall trajectories and to define locations where the road was exposed to further rockfall hazard.

Location of the study area. (a) Topographic map: Carta Tecnica Regionale of the Lazio region, at $\mathrm{1}:\mathrm{10}\mathrm{000}$ scale. Colour orthophoto map: Ortophoto 50 cm copyright 2012 AGEA – all rights reserved. Symbols named S1 and S2 indicate the source areas of the fallen boulders. The red star shows the location where a boulder hit and crossed the SP18 road. Red triangles indicate the boulders found at the farthest locations from the source areas along the talus slope during the field survey. (b) Map of surface deposits obtained through the visual interpretation of the available aerial photographs. (c) Three-dimensional perspective of the study area. Image base and elevation data were obtained ad hoc through a RPAS flight carried out on 10 October 2016. (d) S1 rockfall source area. Oblique photograph taken by the RPAS. (e) Detail of the orthophotograph obtained by the RPAS showing impact points of the boulder detached from the S2 source area. (f) Flight plan of the multicopter for the acquisition of aerial photographs.

[Figure omitted. See PDF]

Study area

The study area is located close to the village of Villanova di Accumoli, in the Accumoli municipality, central Italy (Fig. 1a). In the area, the main vulnerable element subject to rockfalls was a portion of the SP18 road connecting Villanova di Accumoli to the San Giovanni village (Fig. 1a).

The area consisted of a rocky hillslope from which rockfalls detached, reaching the SP18 road (Fig. 1a). Visual interpretation of stereoscopic black and white aerial photographs taken in 1954 at $\mathrm{1}:\mathrm{33}\mathrm{000}$ scale revealed that the cliff from where the rockfall originated is part of a partially dismantled scarp of a deep-seated, disrupted rockslide (Fig. 1b). Additional visual interpretation of stereoscopic colour aerial photographs taken in 1977 at $\mathrm{1}:\mathrm{13}\mathrm{000}$ scale revealed that the main landslide rocky escarpment is characterized by two areas that underwent progressive retreat, forming three cone-shaped, convex, sparsely vegetated talus deposits (Fig. 1b). The analysis of the two sets of aerial photographs revealed that in the central part of the escarpment the talus deposit was mined and partly removed between 1954 and 1977 (Fig. 1b). Outside of the talus deposits, the landslide scarp shows an average slope of around 40${}^{\circ }$, it is mainly bare and bedrock crops out. At the base of the slope, the top of the landslide deposit is covered partially by talus that supports a dense deciduous wood.

In the study area, a regional thrust brings the carbonatic units of the Umbria-Marche stratigraphic sequence over the Laga Formation . The boulders that hit the SP18, and the majority of the boulders found on the talus slope above the SP18 road, belong to the Maiolica Fm., a layered and highly fractured mudstone, upper Jurassic to lower Cretaceous in age.

In the following, we describe the two rockfall source areas, shown as “site 1” (S1) and “site 2” (S2) in Fig. 1a, c, recognized during the field survey.

Site 1 (S1 in Fig. 1c, d) is a 30 m high cliff in fractured limestone 110 m NE from the road. The talus downhill of the cliff has an average slope of about 31${}^{\circ }$. It was partly exploited in the 1950s, and in the quarry area, the average slope is 40${}^{\circ }$. Moving from the source area to the road, the slope shows a decreasing inclination. For a distance of 50 m uphill of the road, a system of scarps, counter scarps and rough and ruined embankments parallel to the SP18 was developed to protect the road during the quarry activities that removed a large part of the original talus. The quarry protection system is now completely wooded. Field observations revealed that the boulders detached from S1, moved primarily along ballistic trajectories (“ree fall”) as suggested by the impact points found on the ground and the scars left on the tree trunks. At impact points on bedrock and talus, the fallen blocks broke up in pieces, the larger of which 0.25–0.30 m${}^{\mathrm{3}}$ in volume, were found on the nearly flat area uphill of the road. Analysis of the scars left by the falling blocks on the trees revealed that the rockfall trajectories did not exceed an elevation above the local terrain of about 1.5 m. Based on the boulders found in the field, the area enclosing the trajectories of the rockfalls detached from S1 is estimated in about 4000 m${}^{\mathrm{2}}$. The pre-failure morphology of the S1 source area was estimated heuristically, by visual interpretation of the orthophotos, which then allowed to estimate the volume of the detached mass (about 7 m${}^{\mathrm{3}}$) from the DEM, where the height of the unstable wedge was measured. Based on the evidence of the fracturing of the remaining unstable rock mass (i.e. orthophotos and field observations), it is possible to hypothesize that the initial detached rock mass consisted of a few large boulders that then further fragmented. The rocky material detaching from S1 did not reach the SP18, but stopped on the talus deposit and on the quarry protection system. Visual analysis of the photographs taken by the RPAS (Fig. 1d) allowed identification of several unstable rock blocks present in the vicinity of S1. Such rock blocks were partially separated from the bedrock by open fractures, which could be further enlarged due to intense rainfall and by frost and thaw cycles, not infrequent in the area from late autumn to early spring. The overall volume of the unstable blocks was estimated in the range of 30–50 m${}^{\mathrm{3}}$ by combining a GIS measure of the area recognized on the orthophoto (10 m${}^{\mathrm{2}}$) and field observation together with 3-D measures obtained from the DEM, for which the height of the unstable rock blocks was estimated to range between 3 and 5 m.

Site 2 (S2 in Fig. 1c) is located along a slope covered by a talus deposit, with a mean slope angle of about 32${}^{\circ }$. From S2, a 1 m${}^{\mathrm{3}}$ boulder detached and reached the SP18 road (Fig. 1c, e). The analysis of the boulder impact points on the road (Fig. 1e) and its final position allowed us to reconstruct the boulder trajectory, from 70 to 90 m uphill of the road, where we were able to identify the most likely source area (Fig. 1a, c). Along the hillslope of S2, six boulders of similar size were found on the talus deposit. It is possible that the seismic shocks worsened the stability conditions of the boulders, increasing the rockfall hazard along the road.

Analysis of the distribution of the rockfall impact points and of the tree scars suggested that the height of the trajectory of the 1 m${}^{\mathrm{3}}$ boulder never exceeded 1 m above the ground. The boulder did not break along its path and stopped 15 m downslope of the SP18.

Methods

Elevation data acquisition

Numerical modelling of rockfalls requires the availability of an accurate digital model of the topographic surface. For this purpose, we performed a dedicated RPAS image sequence acquisition using a SenseFly^® Albris^® multicopter equipped with a 34 Mpixel RGB camera and an on-board GNSS system for the accurate geolocation of the acquired images. A total of 161 nadir photographs with a frontal overlap of 75 % and a side overlap of 60 % were taken at an altitude of 90 m above the ground (Fig. 1f). The altitude above the ground was set up using the dedicated mission planner using the SRTM DEM () for reference. Keeping constant the elevation of the flight with respect to the ground guaranteed a homogeneous GSD across the study area. The resulting photographs cover an area of about 65 000 m${}^{\mathrm{2}}$ with a GSD of about 1.5 cm $\times$ 1.5 cm. Using this ultra-resolution photographic sequence, we prepared an orthophotograph of the entire study area with a nominal ground resolution of about 2.5 cm $\times$ 2.5 cm and a DSM with a nominal ground resolution of about 20 cm $\times$ 20 cm.

Structure from motion applied to RPAS ultra-resolution images

We processed the collected images using Agisoft PhotoScan^® software (http://www.agisoft.com, last access: 25 January 2019) and produced the true orthophotograph and the DSM. A total of 10 GCPs were used to improve the accuracy of the geographic positioning of the DSM and the orthophotograph. GCPs consisted in 40 cm $\times$ 40 cm targets. The positioning accuracy of the resulting DSM and true orthophotograph was evaluated using 70 additional CPs. The geographical coordinates of both GCPs and CPs were obtained using a GNSS RTK VRS (Global Navigation Satellite System, real-time kinematic, virtual reference station) positioning technique. The estimated planimetric–altimetric accuracy of the absolute positioning of the GCPs and CPs was about 10 cm.

To produce a DEM, the fingerprint of vegetation and other elements that could hide the ground were filtered out. For this purpose, we first used a geometric and a radiometric filter available in Photoscan^®. Next, we performed a manual cleaning of the terrain elevation data. After the cleaning, Photoscan^® was used to interpolate the remaining ground points. In some sectors, the interpolation was based on a very limited number of data points due to the coverage of trees, which resulted in an over-smoothed DEM (Fig. 2a).

Shaded relief images of the elevation of the study area. (a) Image prepared using the DEM obtained with the photogrammetric procedure, after vegetation filtering. (b) Image prepared integrating the GNSS RTK data points (yellow squares in a) with the DEM. DEM shown in (b) was resampled at a resolution of 1 m $\times$ 1 m, for modelling purposes.

[Figure omitted. See PDF]

GNSS RTK VRS elevation data integration

In the areas where the elevation information was unavailable (due to the presence of trees), we obtained additional elevation information with a GNSS survey, applying an RTK VRS correction. For this purpose, we used a Leica^® Zeno 20 (https://leica-geosystems.com/, last access: 25 January 2019) device in conjunction with an Android smartphone with GPRS connection for the real-time correction of the positioning (RTK VRS). We acquired a total of 73 points (Fig. 2a), corresponding to an average density of one point every 100 m${}^{\mathrm{2}}$ (i.e. a GSD of approximately 10 m) with a nominal positional accuracy of less than 1 m.

To obtain a continuous DEM covering the entire study area, we merged the two elevation data sets obtained from the aerial and the ground-based surveys. For this purpose, a mask of the wooded area was drawn in the GIS, and a 5 m external buffer was applied to enclose a 5 m wide band containing elevation data from the DEM generated by SfM in the area. Elevation data contained in the buffer were then included in the interpolation of the 73 GNSS RTK VRS data points using the linear Delaunay interpolation method available in the v.surf.nnbathy GRASS GIS tool . Using the elevation data available in the buffer guaranteed that no step-like features occurred at the junction between the original DEM and the DEM obtained by the GNSS RTK VRS survey in the wooded area (Fig. 2b). The GSD of the merged DEM was set to 1 m, the highest resolution handled by the rockfall modelling software STONE .

We stress that the DEM accuracy is not uniform over the entire study area. The information is denser and more accurate in unwooded areas, where the elevation was obtained by the SfM, as opposed to areas where data were collected by the GNSS survey. Despite the information under the tree canopies being less dense than in the rest of the area, the GNSS points were acquired to best represent the ground morphology.

Numerical modelling of rockfalls

Numerical modelling of rockfall was performed using STONE, a software for the 3-D, kinematical simulation of rockfalls . STONE simulates the movement of a boulder along a slope within a three-dimensional approach, considering the moving boulder dimensionless and with all the mass concentrated in the centre of mass, and computing the trajectory at discrete time steps. Within each time step, the boulder can be in one of three “states”, i.e. free fall, rolling, or bouncing. The trajectory of a boulder is computed automatically from the DEM, and it depends on the starting point, the topography, and the coefficients used to simulate the loss of velocity at the impact points or during the rolling state. The coefficients can be obtained from existing thematic data or estimated analysing geological, geomorphological, and land cover maps (Table 1). For each simulation, and for each simulated trajectory, the program allows a random variation in the coefficients and in the initial direction of motion, resulting in an output with probabilistic content. The minimal input required to perform a simulation with STONE consists of four raster maps, containing (i) elevation information (i.e. the DEM), (ii) the location of the rockfall detachment areas (source grid cells), (iii) values of the parameters describing the tangential and the normal energy restitution at each impact point, and (iv) values of the parameter used to describe the dynamical friction coefficient where the rockfall is rolling .

Values of the dynamic rolling friction angle and of the normal and tangential energy restitution coefficients assigned to different terrain types.

Terrain type	Rolling	Normal	Tangential
	friction	restitution	restitution
Landslide scarp	0.30	65	80
Talus	0.65	35	55
Landslide deposit	0.40	45	55

The output of the software consists of three raster maps in which, for each grid cell, the following information is given: (i) the number of modelled rockfall trajectories, (ii) the maximum velocity of the simulated rockfalls, and (iii) the maximum height reached by the simulated trajectory. The three quantities can be used to evaluate rockfall hazard . In particular, we considered the count of the simulated rockfall trajectories as a proxy for the probability that a given grid cell is affected by a rockfall. The larger the number of trajectories, the larger the expected likelihood of rockfall occurrence in a given cell. The grids containing the maximum velocity and height of the possible rockfall were not used in this work. Analysis of the map of the trajectories allows us to distinguish segments of the road that are predicted to be safe from the ones with a non-negligible probability to be hit by a rockfall. Rockfall source areas were selected based on field observations and orthophoto interpretation. Additional source areas were singled out by selecting cells with a slope steeper than 60${}^{\circ }$ according to the measurement of the slope (derived from the DEM) of the potentially unstable rocky escarpments observed in the field.

Assigning the values of energy restitution and dynamic friction parameters required by the STONE model is a somewhat arbitrary, heuristic operation. In this work, parameters were assigned modifying the values selected in previous works according to the field observations.

We have prepared 100 different DEMs modifying the elevation values of the 73 GNSS RTK points by adding delta values to the original elevation data. Delta values were obtained randomly sampling from a Gaussian distribution that reproduces the error values declared by the instrumentation ($\mathit{\mu }=\mathrm{0}$, $\mathit{\sigma }=\mathrm{0.25}$). Each set of modified elevation points was interpolated following the approach described before. The set of 100 DEMs was used to evaluate the spatial distribution of rockfall trajectories considering the topographic uncertainty.

Results and discussion

From the 1938 cells identified as source areas (light blue polygons in Fig. 3a), STONE was run 100 times, one for each modified DEM (see Sect. 3.4), simulating a total of 193 800 rockfall trajectories per simulation (i.e. 100 launches per source cell). The raster map in Fig. 3a shows, for each cell, the mode of the number of trajectories obtained in the 100 simulations, whereas the plot in Fig. 3b shows the values obtained along the SP18 track. Figure 3c shows, for the pixel along the road, the number of simulations in which the trajectory count is greater than 0.

(a) The map shows the mode of count of trajectories of 100 STONE simulations (see text for explanation). Black circles indicate the location of large boulders observed during the field survey at their farthest distance from the source areas (light blue polygons). For each pixel along the road track, (b) shows the values of the map (a) and (c) shows the number of STONE simulations in which trajectory count is greater than 0. The star indicates the position of the impact point along the SP18; blue symbols show the end points of the road track.

[Figure omitted. See PDF]

Considering the map in Fig. 3a, the trajectories affect an area of 56 032 cells (m${}^{\mathrm{2}}$), with the count values ranging between 1 and 652. Furthermore, a total of 32 800 pixels assume a value greater than 1.

Inspection of Fig. 3b–c reveals that the tract of the road downhill of S2 is the most hazardous and that the location of the impact point matches the area where the STONE output assumes the highest values of both the number of trajectories (Fig. 3b) and the hit frequency (Fig. 3c). The portion of the SP18 closest to the source site S1 was not hit by the rockfalls during the seismic sequence. This is confirmed by the STONE simulations, which reveal a total of 957 pixels affected by possible trajectories with a mode value of 1 within the road or downhill of it. The 957 pixels represent 0.49 % of the total number of the simulated trajectories and 0.6 % of the trajectories simulated uphill from this tract of the road (152 300). Such figures can be considered negligible. The plot in Fig. 3c also confirms that in most of the simulations, the trajectories did not reach this tract of the SP18. Contrary, the portion of the SP18 downhill from S2, and which was hit by the rockfall during the seismic sequence, is correctly predicted by the model as subject to rockfall hazard. In detail, the model output identifies a total of 12 123 trajectories crossing the road, accounting for 6.2 % of the total number of simulated trajectories and 29.2 % of those simulated uphill from this portion of the road track. Furthermore, local values reaching 204 trajectories were identified downhill from S2.

Finally, at the pixel scale, we computed the coefficient of variation (CV $=\mathit{\sigma }/\mathit{\mu }$) of the number of trajectories computed in the 100 simulations. CV is a measure of the variability in a sample normalized by the average value of the distribution. The distribution of the CV values shows a standard deviation of 0.17 and a mean of 0.20, which indicates that the model can be considered stable despite the errors introduced by the GNSS RTK topographic measurements.

The workflow described in this paper for the assessment of the road exposure to rockfall hazard was carried out during a seismic emergency phase to support the Italian National Department of Civil Protection in the decision whether or not to safely reopen the SP18 road to traffic. The diffuse presence of trees in the study area would have required an ALS survey to collect ground elevation data under the tree canopy. However, since the seismic emergency framework imposed a strict time constraint, a lidar acquisition was ruled out since it would have required a long planning, preparation, and post-processing time compared to other techniques. Among the existing photogrammetric acquisition approaches, use of RPAS photogrammetry has a number of advantages in terms of feasibility, planning and management easiness, and good-quality data acquisition, which also make it one of the most widespread and used approaches in the private professional sector. From the methodological point of view, the application described in this study can be considered a procedure, applicable during emergencies, for the acquisition of all the required elements that can support the implementation of numerical models for rockfall hazard and risk assessment at a relatively low cost. For areas up to a few tens of square kilometres, the use of RPASs is the cheapest and fastest method for the acquisition of orthophotos and DSMs but, conversely , vegetation hampers the quality of the resulting DEMs, which is fundamental to obtaining reliable numerical rockfall simulations. To overcome this limitation, integration of aerial imagery with GNSS RTK points measured in the field is mandatory for accurate mapping applications, especially when the accuracy required as input by the numerical model is not very high. The availability of ultra-high-resolution images can be very useful for the delineation of rockfall source areas, one of the inputs required by numerical modelling.

The results of STONE are influenced by the values of the tangential and normal restitution and dynamical friction coefficient, which are assigned based on the characterization of the land cover, in most cases accomplished by geomorphological mapping (Fig. 1b). In this case, the detailed geomorphological mapping was carried out based on a rigid set of rules and criteria that make it repeatable and less subjective . To prove the influence of the detail of the geomorphological mapping on the results, we run STONE assigning the parameters by generalizing the initial geomorphological map, considering the entire scarp of the rockslide as a talus, instead of distinguishing the three smaller talus deposits in the same slope. Results showed that assigning the parameters based on a generalized version of the geomorphological map results in trajectories that do not reproduce the location of the boulders observed in the field. This is consistent with the relatively larger value of the friction coefficient ascribed to talus than landslide scarp (). Such evidence shows that the quality of the geomorphological mapping has a relevant effect on model results. The rockfall simulation prepared by STONE confirmed that the portion of the SP18 road affected by the rockfall is rather limited (Fig. 3c). During the emergency phase, this information was provided to the National Department of Civil Protection to support the decision about the road reopening conditions. It was advised that the road could be safely reopened to the traffic, provided that (i) protection measures were installed above the SP18 road where the model indicated the exposure of the road to rockfall trajectories and (ii) that the rockfall barriers installed should take into account the size of the boulders recognized in the field.

Modelling results show that outside the most hazardous part of the SP18 (Fig. 3), only a few locations are potentially affected by rockfalls. Here, the pixels that could potentially be reached by rockfalls along the road show count values of 1. It is worth noting that, over 100 trajectory simulations for each source pixel, a count value of 1 suggests a probability of occurrence that is equal to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$. It actually corresponds to probability values much smaller since most frequently a single pixel can be crossed by trajectories starting from different (even not so close) locations. In the case of the tract of the road threatened by the site S1, the 957 trajectories that could reach the road represented 0.6 % of the total number of simulated trajectories (152 300). This estimate was obtained by counting the number of pixels used as source area located uphill of the road in S1 and the total number of trajectories that reached the road in the mode map (Fig. 3a). Such a method represents a semi-quantitative estimate. Like other rockfall modelling software, STONE outputs a frequency map, whereas the evaluation of the probability of rockfall still represents an open problem due to the presence of the pixels where modelled trajectories overlap with others modelled from other source areas. Hence, a rigorous probabilistic estimation would require that each trajectory should be modelled keeping track of the source area. Possible future development of STONE or similar rockfall models should be able to define a probability value, which would allow better comparison of the rockfall hazard conditions in different areas.

STONE does not take into account the vegetation effect, although, in this case, the presence of a thick wooded area at the foot of the slope acts as a natural rockfall attenuator between the source areas and the road. As already pointed out by , despite the rockfall modelling being aimed at describing possible rockfall scenarios given the current state of the places, it appears poorly precautionary to consider the barrier effect of the trees in the final assessment. Future wildfires could leave the road unprotected, causing an increase in the overall rockfall hazard.

A non-negligible point of discussion consists in the caveats that should be taken into account when using model outputs that contain sources of uncertainty at some step of the procedure to support decision-making. In this case, for instance, a DSM at 20 cm GSD, which is inhomogeneous across the study area due to the presence of the tree canopy, was acquired. The entire DSM was resampled at 1 m GSD, which is the minimum working resolution of the model STONE and a reasonable resolution able to capture micro-morphological features and portray the general morphology reconstructed under the tree canopy by the GNSS RTK survey. Despite the average GNSS point density being one point per 100 m${}^{\mathrm{2}}$, the acquisition was planned to best represent the ground morphology, which is characterized by a series of two scarps and counter-escarpments probably built during the quarry exploitation. In our case, it was not possible to obtain a better acquisition because of the following practical limitations: (i) the GNSS signal under the tree canopy was unstable; hence almost every measure had to be repeated several times to obtain the minimum required accuracy of 1 m, and (ii) the GPRS signal was also poor, causing several interruptions in the RTK connection. For these reasons, when acquiring such data in the field, it is advised to collect as many data points as possible. Moreover, in places where the environmental conditions are unfavourable, acquisition should be performed to catch the most representative morphological features.

The overall cost of the materials (i.e. the RPAS, software licenses, GNSS receiver) used in our test case totalled about EUR 23 000, which is considerably less than lidar acquisition at a comparable resolution. The entire procedure was carried out in 4 working days (estimated to be about 15 working days for one person), including the initial survey, the RPAS and the GNSS acquisitions, the photogrammetric processing and the DSM filtering, and the integration of the RPAS DSM and the GNSS data points. The application of the presented procedure based on a RPAS acquisition has an overall cost that is lower compared to the same procedure involving an ALS survey.

Conclusions

The study described in the article was conducted during the seismic sequence that hit central Italy in the period between 24 August 2016 and January 2017 to support the National Department of Civil Protection in identifying the conditions to safely reopen the road exposed to rockfalls. The study consisted in the application of a rockfall numerical model (STONE) prepared using a digital terrain model generated by the integration of a DEM obtained by photogrammetry applied to remotely piloted aerial system (RPAS) acquisition and a GNSS RTK (Global Navigation Satellite System in real-time kinematic mode) survey in densely vegetated areas. The numerical model allowed us to perform a semi-quantitative evaluation of the residual rockfall risk posed to the road SP18. It was observed that the tract of the road that had been hit by a rockfall during the seismic emergency was predicted as unsafe by the model since 29.2 % of the total simulated trajectories uphill from the tract of SP18 reached the road. The remaining portion of the studied tract of the SP18 was reached by 0.6 % of the modelled trajectories uphill from the tract of SP18, and hence its exposure to rockfall was considered negligible. It is worth mentioning that results of this study were used to set up protection measures (i.e. elastic barriers) along the track of the road more exposed to rockfall impacts.

The described procedure is general and can be successfully applied during emergency phases and in mountainous regions like the study area. The complete photogrammetry-based procedure described in the paper was carried out in 4 working days by a multidisciplinary team, which correspond to a total of 15 working days for one person at an overall cost of EUR 23 000. Considering all the working phases, it is faster and cheaper than lidar surveys and has fewer logistic constraints. Conversely, it is limited by high vegetation cover. The GNSS RTK survey can partially solve this limitation, but it is important to know that the final DEM has a different quality under the tree canopies and on bare areas.

STONE fails to reproduce the position of the boulders actually fallen during the seismic sequence if parameters are assigned based on generalized geomorphological maps of the covers. Therefore, it is advised that the quality of the input geomorphological mapping be as high as possible and compatible with the resolution of the input DEM.

Data are available at 10.17605/OSF.IO/VKGW8 ().

MS contributed to the field survey and the geomorphological mapping, acquired the GNSS RTK data, analysed the results, and wrote the text. MA ran the model STONE and contributed to writing the paper. MB acquired and processed the RPAS data. MC contributed to the geomorphological mapping and reviewed the text. DG acquired and processed the RPAS data, analysed the results, and contributed to writing the text. FG contributed to the field survey and reviewed the text. IM acquired the GNSS RTK data, analysed the results, and contributed to writing the text. PR reviewed the text.

The authors declare that they have no conflict of interest.

In this work, use of copyright, brand, trade names, and logos is for descriptive and identification purposes only, and does not imply endorsement from the authors, or their institutions.

Acknowledgements

The work was partly funded by the Italian National Department of Civil Protection. Edited by: Norman Kerle Reviewed by: two anonymous referees