Debris Flow and Rockslide Analysis with Advanced Photogrammetry Techniques Based on High-Resolution RPAS Data. Ponte Formazza Case Study (NW Alps)

: The use of a Remotely Piloted Aircraft System (RPAS) for the characterization and monitoring of landslides has been widely improved in the last decade. In particular, the use of this system is particularly effective for the study of areas prone to geohazards. Zones affected by landslides, such as rock slides and debris ﬂows, are often quite critical in terms of accessibility due to unstable blocs that can strongly limit the direct access to the studied area. In this paper, we present the case study of Ponte Formazza in NW Italian Alps. In June 2019, a massive and complex debris ﬂow re-mobilized about 300,000 m 3 of a rockslide deposit that occurred in 2009. In this particular environment, we tested traditional, direct and mixed photogrammetric approaches using various conﬁgurations of Ground Control Points (GCPs) of the photogrammetric block and by calculating the relative errors. The minimum conﬁguration of GCPs was established to reduce in situ measurements without de-grading the accuracy of the cartographic products. The images of three RPAS campaigns (2017, 2018 and 2019), processed with a Structure from Motion (SfM) technique, allowed us to obtain very high-resolution orthophoto and digital surface models (DSMs) before and after the 2019 event. A few GCPs, geolocated with a Global Navigation Satellite System (GNSS), improved the orthophoto and DSM quality (Root Mean Squared Error RMSE 5 cm) even in the areas far from the drone deployment. The availability of high-resolution models has been fundamental for the identiﬁcation of the volume changes. Furthermore, the 3D view supported and completed the geomorphological mapping of affected areas, particularly in the areas where the ﬁeld survey is dangerous. The use of ancillary meteorological data and Sentinel-2 satellite images allows for a better deﬁnition of the kinematics and the predisposal and triggering factors of the 2019 debris ﬂow.


Introduction
In the last few years, the use of unmanned aerial vehicles (UAVs) or Remotely Piloted Aircraft Systems (RPASs) in the study of natural hazards has significantly increased [1,2].
The improvement of autopilot or semi-autopilot systems, high-resolution digital cameras, and GNSS and inertial systems has allowed RPAS to increase the precision and accuracy of the data. Captured images are geocoded by various approaches requiring onboard and/or ground-based solutions. UAVs are usually equipped with positioning apparatuses. To improve the final positioning accuracy, ground control points (GCPs) are positioned, surveyed and used in the image post-processing. A possible technical solution adopted for improving the quality of the GNSS position is the use of Real-Time Kinematic (RTK) receivers [3,4] onboard UAVs. This approach, called direct photogrammetry, allows RPAS surveys without using traditional GCPs on the ground [5,6]. However, such a technique is a recent technology and needs to be tested and compared with traditional photogrammetry approaches.
The Structure from Motion SfM [7][8][9] is a codified process that can generate a 3-D point cloud from the RPAS geocoded images. Starting from a densified pointcloud, SfM software can also generate a Digital Elevation Model (DEM) and orthoimage. These map layers can be used to extract image features (e.g., colors, texture) or geomorphological features from DEM or 3-D image analysis (e.g., slope, aspect, hydrographic network).
Among the natural processes, the landslides are suitable for UAV analysis, and in the literature, it is possible to find several examples of this [10][11][12][13][14][15]. The multi-temporal analysis of UAV data allows us to detect and study the morphological changes that have occurred to landslides using images [16,17]. The DEM time series provides detailed multi-temporal sets of 3D surfaces useful for vertical displacements and volume change investigation [7,18,19].
The debris flows in mountain areas are among the main risks for human life and cause damage to infrastructure and buildings [20]. Like many other mountainous areas, the Alps are affected by debris flow, especially in the summer and autumn seasons [21][22][23].
Several studies observed an increase in the frequency and intensity of debris flow in the Alps that could be related to climate change that provides more extreme rainfall events [24][25][26][27]. Other events that could activate debris flows are the sudden outburst of glacial lakes [28] or the increase in source materials like rockfall accumulation [29] that creates unstable deposits that could be activated during extreme rainfall events.
An important element that should be considered in the study of debris flow is the sequence of events that generates the critical condition for the activation of the gravitational process. In particular, the identification and the characterization of the deposits that the debris flow has mobilized are fundamental elements that can be evaluated for a correct evaluation of the possible repetition of the process. The identification of the critical sequence of the critical events and the evaluation of the possibility that a new debris flow could occur are the basis for the improvement in early warning systems [30]. To reach this goal, very high-resolution orthophoto and digital surface model (DSM) obtained by RPAS could be a precious help [31].
In this work, we present the application of direct and mixed photogrammetry, SfM based on RPAS data to study the debris flow that hit the village of Ponte Formazza (NW Alps, Italy) on 10 June 2019. The availability of pre-an post-event images made with RPAS data allowed detailed mapping of debris flow. We tested and compared the accuracy and the volume estimation of different photogrammetry approaches: traditional, direct, and mixed. We applied the SfM technique to extract pre-and post-event DSMs and calculate the topography and volume variation. We also used low-resolution satellite data of Sentinel-2 to make a rapid mapping of the affected area and estimate the snow coverage. Meteorological and snow data allowed us to reconstruct the predisposing and the triggering factor of this debris flow. The potentialities and limitations of the RPAS applied to large a landslide study are discussed here.

Geological and Geomorphological Settings
The Ponte Formazza is a small village located in the upper Formazza Valley, Piemonte region (NW Italian Alps), here several hotels and touristic facilities were built on the east side of the Toce river.
The study area ( Figure 1) is inside Lepotine nappe of the Pennidic domain, and it is located at the tectonic boundaries between Teggiolo formation (made by calcschist, metapelite lithology) and the Antigorio formation (orthogneiss lithotype) [32].
The main tectonic lineaments have directions SW-NE, E-W and NW-SE, these lineaments drive at the local scale the fracture systems of the rock mass. The upper Formazza Valley shows a North-South direction, and its morphology is driven by the glacial erosion with a typical U-shape. Near the village of Ponte Formazza, on the left flank of the valley, there are the basins of some small creeks (Enni and Rich) that contributed to the debris flow. We can identify three main sectors along the slope: (i) the low sector with a gentle gradient made by an alluvial fan and talus and rockfall accumulation; (ii) the middle sector 1500-2000 m located on the tectonic contact between Teggiolo formation and Antigorio formation, here the creek erosion isolated a cliff of a calcschist where a system of fractures separated several blocks prone to fall; (iii) the upper part of the slope (2000-2800 m) shows a glacial landform with the talus and periglacial deposits that cover the bedrock. The area involved in the June 2019 debris flow comprises the rockslide deposit located in the lower part of the slope and the bottom of the valley.

Previous Events: April 2009 Rockslide
The area of the 2019 debris flow was previously affected by several instability events; debris flows have already hit this area, as in August 1987 and in October 2000 [33], as reported in the Arpa Piemonte database, even if these were much less massive compared to the 2009 event. The geomorphological evidence from pre-2009 aerial photos suggests that the area was also periodically affected by large rockfalls.
The most massive event in recent time occurred on 9 April 2009, when a massive rockslide ( Figure 2) created a large deposit that can be considered the main predisposing factor for the event that occurred in June 2019.
In 2009, the rock mass release area (R.S.A. in Figure 2) was located on a ridge at 1950 m elevation and oriented parallel to NE-SW striking steeply dipping faults. The initial failure likely occurred along such a fault, then the detached mass fell in a steep gorge on the north side of the cliff [34][35][36].
The rockslide mass moved quickly along the slope, which has an inclination of about 40 • . The direction of the detached mass was about 210 • with a drop height of 530 m over a horizontal distance of 600 m. The estimated total volume was 700,000 m 3 . On 1 May 2009, a few days after the rockslide, a small shallow landslide affected the colluvial deposit where vegetation was eradicated (E.T. in Figure 2) [34]. The unstable rockslide accumulation (U.R.A. in Figure 2A), formed by some boulders of up to tens of metres in size, was partially blocked in a steep channel at the cliff base. This meta-stable deposit was massively moved by debris flow in June 2019. The field survey of 2017 also allowed us to detect an unstable block (a' in Figure 2B) that collapsed during the 2019 event.
From 2009 to 2019, the studied area was characterized only by limited rockfall detached from the unstable cliff or some movement of blocks of the rockfall deposit.
In November 2017 and May 2018, within a project funded by the Agency for the Informative System of Piemonte Region (CSI is the Italian acronym), we mapped the rockslide accumulation and part of the gorge close to the cliff using RPAS. These two surveys allowed us to acquire high-resolution DSMs and orthophotos to compare the pre-event morphology with the post-debris flow topography acquired in the 2019 survey.

Materials and Methods
The correct description of the 2019 debris flow requires the identification of the triggering and predisposing factors and the characterization of the effects of the slope instability. For this reason, we collected rainfall and snow depth data to assess the triggering factors. We used satellite data to map the limit of the occurred debris flow and the known past events and the snow coverage. RPAS data of pre-event (2017 and 2018) and post-event (2019) were used to create high-resolution orthophoto and DSM of the studied area. We use these data to detect the geomorphological feature of debris flow and quantify the changes.

Meteorological Data
We collected meteorological data (rainfall and snow depth) from the Agency for Environmental Protection of Piemonte Region (ARPA is the Italian acronym) meteorological database [37], to understand the triggering factors of debris flow. The data come from the closest stations to Formazza ( Figure 1A), located at a different altitude (Table 1). We first calculated the snowpack thickness in the contributing basin on the basis of: (i) the snow gauge data; (ii) Sentinel-2 images; (iii) the elevation from the Digital Terrain Model (DTM) of ARPA. Snow melting in early summer is an important parameter that has to be carefully considered because it could contribute to soil saturation. We also considered the nearest station data (Formazza Bruggi) to check the hourly peak rainfall that triggered the debris flow. Table 1. The meteorological stations of ARPA Piemonte used for this work were also located in Figure 1.

Low-resolution Mapping and Snow Coverage Estimation
Before planning a detailed RPAS survey, we made a rapid mapping of the affected area using the Sentinel-2 images. In particular, we made a difference in the pre-and post-NDVI spectral index (NDVI var ) ( Table 2). This approach can identify areas where vegetation was eradicated or covered by debris; using an approach similar to the one adopted to identify flooded areas [38,39], we identified the area covered by the debris flow deposit. The data of Sentinel-2 were also used to map the snow coverage before the debris flow event. We mapped the area affected by the rockslide in 2009 using the middle-high resolution orthophoto provided by Regione Piemonte. The data allowed mapping the detachment sector and part of the unstable rockslide accumulation (C.B. and U.R.A in Figure 2) that was not covered by RPAS flight ( Table 2).

High-Resolution Mapping Strategy
The area affected by a rockslide in 2009 was mapped at very high resolution in November 2017, July 2018, and July 2019 (Table 3) with a Phantom 4 Pro UAV. The flights have been planned to acquire a detailed DSM of the accumulation sector and part of the gorge occupied by debris; due to its vertical position and the distance from the take-off point, it was impossible to reach and map the cliff source of the collapsed block. Considering these logistic limits, the area affected by 2019 debris flow was divided into two main parts (S1 and S2, in Figure 3) that partly overlap (S1+S2, in Figure 3).
The upper part of the slope (S1 and S1+S2 in Figure 3) was mapped by traditional RPAS and processed with SfM (Table 3). This sector was too dangerous for installing new GCPs, and we used only the two survived and five of the new GCPs located at the bottom of the surveyed area, in safe sectors. From processing, we obtained a DSM for most of the debris flow area used for the comparison with the pre-event dataset of 2017 and 2018.
The survey made in 2019 in the lower part of the slope (sectors S2 and S1+S2, in Figure 3) was done using a DJI PH4 UAV equipped with double-frequency GNSS receivers.
It was possible to use all the new 13 GCPs made of 50×50 cm 2 plastic square sharply visible in the RPAS photo ( Figure 3C), which allowed us to test three different photogrammetry approaches (direct, traditional, and mixed). In direct photogrammetry ( Figure 3A"), the coordinates of the Projection Center (PCPs) are used, these are obtained from data of onboard GNSS and the camera parameters ( Table 4). The traditional photogrammetry uses the GCPs measured on the field ( Figure 3A'). The mixed approach uses both types of coordinates. We obtained the DSM for the lower sector of the slope (S1 and S+2 in Figure 3) from these processings, which only partly overlap with the pre-event DSMs.

Mapping with RPAS: Processing Strategies in Emergency Cases
The dangerous and inaccessible conditions and the urgency to survey areas affected by landslide and debris flow phenomena led to using RPAS photogrammetry mapping. In these emergency cases, the most time-consuming and dangerous activity in the field is the positioning and measure of GCPs.
In 2017, we installed GCPs measured with GNSS on the area up to the cliff base for photogrammetric orientation of images.
Unfortunately, most of the GCPs installed for 2017-2018 flights were removed or covered by the debris flow, and only two markers installed on the stable cliff on the flank were still visible in July 2019.
Due to the greater extent of the debris flow in 2019 compared to 2018, the higher danger and the work in situ to make the safety of homes and the road, it was impossible to create and measure GCPs inside the landslide.  Even if they are quick and precisely re-measurable with GNSS, these points can be located only at the landslide area edges and not where the photogrammetry technique requires it. On the other hand, with the so-called "Direct Photogrammetry" theoretically, the ground-based measurements are unnecessary because the GNSS receiver installed onboard the aircraft provides the Projection Centers Point (PCP) positions and the inertial measurement unit (IMU) provides the camera's angular assets. In both cases, however, these measurements require complex post-processing.
Unfortunately, the inertial instruments of adequate precision have weights and costs incompatible with small RPAS. For this reason, usually the information derived from the GNSS position of the PCPs is the only one used for the post-processing procedure. Adopting a high percentage of photographic coverage, the set of frames (called "block") can be oriented in the cartographic reference system even without the precise but expensive inertial measurements.
The coordinates of the PCPs enter into a software called "aerial triangulation" (AT), which allows us to determine the angular asset of all the frames and provides statistical parameters to evaluate the precision of the procedure. Thus, it is possible to improve this accuracy by adding GCP measurements like in classical photogrammetry. The approach used a 2019 survey to combine the two techniques ("mixed approach") to improve the accuracy and the reliability of the results.
The research activity aimed to find what could be a good compromise in these emergency cases. This compromise between precision, safety and speed consists in the use of all PCPs and a limited number of GCPs. For this purpose, the accuracies of the aerial triangulation calculations in some typical scenarios were evaluated and compared.
To achieve the survey to be sufficiently precise, i.e., a few cm, it is necessary that the coordinates of the PCPs also have similar accuracy and a mass-market GNSS receiver is insufficient for this purpose. It is mandatory to use a receiver that records the code and phase data to be processed in post-processing, such as the uBlox-F9P receiver. After the flight, the data must be treated with differential processing, i.e., also using the data of a second receiver close to the survey area. The post-processing software must also achieve the fixing of phase ambiguities to integer numbers.
The photogrammetric block was compensated with a reference both with GCPs and with PCs, analyzing different configurations.
The RPAS used, the DJI PH4 model, is equipped with a dual-frequency receiver to achieve this positioning accuracy. The GNSS antenna's position was obtained with a onesecond rate of acquisition with post-processed kinematic (PPK) methods [40], starting from a base station near the rockslide. Fixing of the integer phase ambiguity has been achieved with a formal precision always better than 1 cm in the three coordinates. In turn, the base station position coordinates have been obtained with a static relative positioning in the ETRF2000 system, assured by the correct use of GNSS permanent stations network "SPIN." (https://www.spingnss.it/spiderweb/frmIndex.aspx accessed on 2 May 2021).
The positions of the PCPs at the instant of the shooting were obtained with two subsequent calculations. The phase center of the GNSS antenna does not coincide with the center of the camera, and the opening instant does not coincide with the exact second position of the antenna already computed First, a temporal interpolation of the phase center's position was done to bring it backward or forward to the instant of shooting. The camera center position was then obtained from the phase center by correcting the known spatial eccentricity between the two centers. The eccentricity vector, called lever arm, was then rotated in the ground reference system, knowing the aircraft attitude's angular values. For this purpose, the low precision of the drone's onboard IMU is sufficient due to the distances of a few centimeters. Even an angular error of 3 • makes it possible to carry out this operation with centimeter accuracy. Therefore, the camera centers' position and accuracy do not depend solely on the accuracy of the GNSS antenna phase centers, and coordinates and accuracy enter together in the aerial triangulation procedure.

Traditional, Direct and Mixed Photogrammetry Techniques Using an SfM Software
The adopted aerial triangulation software used is part of the Agisoft Metashape™ suite that also performs the dense images matching. Agisoft Metashape™ is probably one of the most famous Structure for Motion software (SfM). Usually, the products of this software are the digital surface model (DSM) and orthophoto. These results are obtained with a so-called dense matching, i.e., the output of correlation, pixel by pixel of digital images.
All the SfM software are working with a series of routines that, in order, perform these tasks:

•
Finding homologous points, named key points, in all photograms.

•
Performing the relative orientation of the images (also named camera frames orientation).

•
Performing the absolute orientation of the block of photograms with some measures (GCPs and/or PCs). This part is classically named "Aerial Triangulation" (AT).

•
Building the dense cloud.

•
Building of mesh and texture. • Building of DSM, DTM, and Orthophoto.
In this paper, we present several tests conducted in the AT phase of the adopted SfM software. Tests aim to identify the most accurate results that we could obtain in the AT phase limiting the number of GCPs. We defined different scenarios, and we evaluated the accuracy obtained in AT that determine the accuracy of the final products (orthophotos and DSM). The AT procedure also provides the estimate of the metric characteristics of the camera (internal orientation of camera and distortions of the optic). This estimation is necessary, in particular with cameras like the one placed on the used RPAS that cannot be considered a professional camera, and therefore, it has non-negligible optical distortions. If it is not possible to evaluate accurately the mathematical parameters that model these distortions, the frames would be unusable for metric purposes. The estimation occurs together with the camera orientation parameters during the AT. The final result is usually conditioned by the presence of a high number of constraints (known coordinate points) in the photogrammetric block. The known coordinates points are made of GCPs and the PCPs of the images sequence.
To evaluate the amount of degradation of the accuracy due to the scarcity of measurement points on the ground of the AT process, we considered three different scenarios:

2.
Direct photogrammetry: the use of Projection Centre Points (PCPs).

3.
Mixed approach: the use of both PCPs and a variable number of GCPs (from one up to 13).
To compare the results and evaluate the volumetric and precision differences, we considered the third scenario (with all PCPs and GCPs) as the more precise reference. In intermediate scenarios, the differences between the unused points' coordinates and those obtained from AT have been evaluated (Check Points-CKP).
Both GCPs and PCPs enter in AT with their weight, i.e., their accuracy. Planimetric accuracy of 1 cm and an altimetric accuracy of 3 cm were achieved and established for the GCPs. For PC, a planimetric accuracy of 3 cm and an altimetric accuracy of 5 cm were used.
The configurations are summarised in Table 5. In Case 1, the AT uses only the GCPs. In Case 2, only the PCPs are used. The other cases are mixed solutions with all the PCPs and with a variable number of GCPs.

Volume Estimation and Geomorphological Mapping of the Debris Flow
The volume changes for the sectors S1 and S1+S2 ( Figure 3) were obtained by a DSM of difference computed in QGIS software. Specifically, we first divided the debris flow into different sectors avoiding the areas covered by vegetation. In a second step, we multiplied the area of these sectors for the average variation of elevation (DEM of difference) to obtain volume change. We calculated the volume variation using the 2019 and 2017 DSM, which shows a more overlapping area than the 2018 DSM.
For the lower part of debris flow (sector S2 in Figure 3) covered only by the 2019 RPAS RtK survey (A), we evaluated the volume with a different methodology that is made of three steps, made with Agisoft Metashape™ software: • we first cut out from DSM the area covered by debris ( Figure 4B); • we proceeded to interpolate the clipped part to reconstruct the surface before the event ( Figure 4C); • we calculated the difference between the original raw DSM and the reconstructed DSM pre-event. This allowed us to estimate the volume also of this sector (see Figure  8A in the results). A more precise volume estimation was also made by removing the vegetation on both pre-and post-event DSM, using a manual clipping. Then, using interpolation, the data gaps were filled.
Using the results obtained from SfM and volume change estimation, we used QGIS to map debris flow boundary and geomorphological sectors manually. The high-resolution DSM and orthophoto were also used to define the shape of the 2009 rockslide better and estimate the size of the most massive boulders.
The pre-and post-event DSM difference was used to compute the volume changes in different sectors and measure the size of some large boulders mobilized by debris flow. Using Qgis2threejs Plugin [41] for QGIS, we create a 3-D view that helps the detailed analysis of the deposits that we could not reach with field surveys. The use of 3D also helped reconstruct the event dynamic and trajectories of some most massive boulders.

Results
The main results of the 2019 debris flow analysis are shown in the following paragraphs.

Meteorological Analysis of the 2019 Event
The debris flow was triggered by intense rainfall on 10-11 June 2019, with accumulated rainfall of 100 mm in 24 h and a peak of 25 mm/h ( Figure 5B). These values are not extreme for the area; however, they occurred in the period of intense snow-melting that probably saturated the soil ( Figure 5A). From the data of snow depth of the nearby stations and the Sentinel-2 images of 01 June 2019 ( Figure 6B), it is possible to estimate snow depth and coverage at the moment of triggering rainstorm ( Figure 6A). Half of the basin that contributed to debris flow was covered by snowpack that in the highest sector (>2500 m) reached about 200 cm of thickness. From the end of May 2019, the snowpack starts a rapid melting at 8 cm/days. All these factors contributed to soil saturation at the moment of the storm.  Data from ARPA Piemonte [42] show that a peak of rainfall of about 15 mm in 10 min observed from RADAR data could be the final trigger of the debris flow. The rainfall that triggered the debris flow is coherent with several intensity/duration thresholds in the literature [43][44][45].
The main difference with previous debris flow events in this area, as the events of 1987 and 2000, is the massive amount of new debris released by the 2009 rockslide that became part of the debris flow mass. In particular, the rockslide accumulation that holds the canyon of the Enni stream represented an important source of material, and it could have temporally dammed the creek. The creation of a small temporary lake and its sudden collapse could amplify the magnitude of the debris flow and the volume of deposits involved in the event [42]. Besides, the vibration and the movement caused by debris flow contributed to the detachment of the large unstable block (a' in Figure 2) that rolled and bounced until a distance of 30 m to the base of the wall.

Low-Resolution Mapping
Low and medium spatial resolution images characterize the rockslide that occurred in 2009 and the 2019 debris flow. The combined study of these two events is fundamental as the rockslide created favorable conditions for the subsequent debris flow activation.
We first mapped the rockslide (Figure 2) using aerial images of ARPA Piemonte (Table 2) and on the basis of field surveys made by [34,35].
The processing of Sentinel-2 images ( Figure 7A,B) allowed us to map the area affected by debris flow in 2019 quickly. Thanks to the first post-event images available a few days after the event, the affected area map helped plan RPAS flights made on 2 July 2019. Moreover, the upper part of the debris flow was mapped only with Sentinel-2 and ground-based photos because this sector is too far from the RPAS take-off point.  Figure 7C shows the NDVI var perfectly matching the boundary of the debris flow area furtherly mapped with RPAS, especially in the bottom deposits sector. It is also possible to see that the decrease of NDVI is higher in the lower sector of the slope where vegetation was erased. The areas already occupied by debris and boulder of the 2009 event show a negative NDVI var because the debris flow buried pioneer vegetation of the 2009 rockslide deposits.
As reported in the previous paragraph, the Sentinel-2 images also allowed mapping the snow coverage area.

High-Resolution Mapping Results
In this section is presented the results of the high-resolution mapping of debris flow based on RPAS data.

SfM for DSM and Orthophoto Production
For the upper part of the slope (S1 and S1+S2 in Figure 3), we processed images of 2017, 2018, and 2019 using the SfM technique (using Agisoft Metashape™) to obtain point clouds DSMs, and orthophotos. The software also gives a report about the precision and accuracy of elaboration, time, and processing parameters. We tested several processes in order to obtain the most accurate results. Table 6 resumes the main characteristic of the RPAS acquisitions. It is possible to see that RMSE is within 2 cm on the plane and within 5 cm, including the vertical component. These values are accurate enough to map the topography variations (in the order of several meters) after the 2019 debris flow. Table 6. SfM processing of RPAS data: resume of the orthophotos and DSM spatial resolution and accuracy. The obtained results from the flight of 2019 show a higher RMSE related to the GCP distribution. As mentioned before, for the 2017 and 2018 flights, the GCPs used were almost the same and equally distributed in the affected areas, while for the 2019 flight, only two were surveyed, and the other new seven are located on the lower part.

Accuracy Results in SfM Software for the Aerial Triangulation Cases
In the S1+S2 sector, we tested different photogrammetric approaches using the 2019 survey. To verify the accuracy of the AT, even with a very low number of ground measurements, some elaborations have been carried out corresponding to three cases. These cases are indicated in Table 5. The results obtained are reported in Table 7. Recall that the accuracy of the digital terrain and volume calculations largely depend on these preliminary operations. The obtained results show that the direct photogrammetry solution (Case 2), if not supported by precise CPs measurements, is less accurate than the traditional (Case 1) and mixed approaches. In particular, the mixed one requires only 3 GCPs (Case 3b) for obtaining a good accuracy (planar residuals about 2 cm and vertical 5 cm). The best result is obtained using the whole 13 GCPs (Case 3c).
We propose a detailed description of each case and the obtained residuals with graphical representation in the Supplementary Material (Figures S1 and S2).

DSM and Volume Accuracy in the Lower Sector of Landslides
As shown in the previous section, case 3c is the solution we used for the official DSM creation and volume calculation. However, we also used solution 3b to evaluate how the higher residual effectively affected the volume estimation.
We first computed the volume difference between the 2019 DSM created with cases 3c and 3b for the whole areas covered by the RPAS survey. The results show a difference of 3112 m 3 over an area of 114,226 m 2 , corresponding to an average elevation difference of 2.7 cm. This difference is less than the expected bias, usually 2-3 times the ground sampling distance, GSD.
After we computed the effective volume changes caused by the landslides using the DSM of 2018 for comparison, we removed vegetation and buildings from both DSMs to have an effective digital terrain model (DTM) of difference and a more realistic volume change.
The results Figure 8A show that the volume of accumulated material is respectively 30,909 m 3 (Case 3b) and 31,064 m 3 (Case 3c). The difference between direct photogrammetry with only 3 GCPs and full solution (3c) is about 720 m 3 (2.4%); this is excellent accuracy results for geomorphological analysis. It is possible to appreciate that the excavators were also removed from the volume computation. Using the methodology described in paragraph 3.4, we calculated the DTM of change and the volume for the low-energy sector of the debris flow ( Figure 8B) where there is no available 2018 data (12s in Figure 13). In this case, the estimated cumulated volume is about 155 m 3 , a modest amount compared with the rest of the debris flow volume (described in Section 4.3.4).

Rockslide and Debris Flow Volume Estimation and Mapping Based on SfM
The High-resolution RPAS surveys carried out in 2017 and 2018 allowed to obtain an extremely detailed orthophoto ( Figure 9A,B) and DSM of the rockslide accumulation ( Figure 9C). The 3-D views allowed us to measure and estimated the volumes of large boulders as on the field. It is interesting to note how most of the GCPs installed in 2017 are still well visible in the 2018 orthophoto (e.g., Figure 9A',B').  Table 8). This is in agreement with the fact that no significant event occurred during the period. The main noise is related to vegetation or shadow areas created by some blocks. The low 2018-2017 DSM difference a confirmation of the excellent quality of geocoding obtained with SfM. The orthophoto and DSM obtained from the SfM allowed us to compare the change that occurred after the June 2019 debris flow.
The orthophoto allowed a precise mapping of most debris flows ( Figure 10A), especially in the lower part, and at the downslope of the wall. Here the relatively low-energy flow created fine-sized debris deposits (from the decimetric block to the sand matrix). Some UAV footage available on the Local Team web video platform [46] shows the debris flow deposit just a few hours after the event. The videos helped to understand some depositional patterns of the restored area at the moment of our RPAS flight 15 days later. The difference between DSM 2019 and 2017 allowed us to understand the complexity of the processes during the 2019 event ( Figure 10B). The main events and their geomorphological evidence can be resumed as follows: (i) The canyon of Enni Stream was utterly washed out of the 30 m thickness of debris and boulders, and now the bedrock is not covered by rockslide deposits; (ii) the giant block a' that was positioned at the top of the fan fall down during the 2019 event following a linear trajectory independent from the rest of debris flow; (iii) most of the debris mass was deposited in the middle sector of the slope, here it is possible to see how the Enni stream migrated from the central part of the fan to the left sector joining the Rich stream; (iv) a sector characterized by the presence of a heterogeneous deposit composed by large blocks in a coarse-grained matrix due to the debris flow transport effect and the rolling component of the movement of unstable rock blocks; (v) the contribution of Rich stream energy and water to the lower sector of the debris flow that accumulated 10 m of materials behind the wall built after the 2009 rockslide; (vi) part of debris flow overpassed the wall, and the rest get around the wall and reached the main Toce river, with two separate branches. These flows hit several buildings, roads, and some sports facilities. The relatively low-energy debris flow avoided more severe damage to buildings, and the few residents were evacuated some minutes before, avoiding life risk.
In the general 3D view (Figure 11), it is possible to detect the area of erosion in the upper part of the channel, the detachment on the large boulders (a'), and its trajectory until its final position (a"). The main accumulation (m.a.) is located at an altitude of 1500 m. In the lower part of the slope, the debris flow shows the fine-size of the blocks (i.a.), and the restoration work moved part of this debris (l.a.e) to make a new section of the wall (a.a.). In Figure 12, a 3D view of the areas with the largest boulders is represented with more detail. It is interesting to note how the DSM difference allowed us to distinguish the boulder "c" that is part of the 2009 rockslide, from the block a" derived from the new detachment and the rest of the boulders like "b" that is part of the debris flow.

Geomorphological Interpretation
We finally produced a geomorphological map of the 2019 event ( Figure 13) based on the orthophoto, 3D view, and 2019-2017 DSM differences (Table 9). We mapped four main geomorphological groups of sectors:

•
Erosion sectors (red/violet colors in Figure 13). Here it is possible to find the eroded channel (1) (source of most of the material mobilized; the area of the detachment of the block a') (2) and the erosion caused by its rolling and bounce (4 a/b); a small area of erosion (3) on the left side of the rockslide deposit related to debris flow action increased by the Bich stream. • A small transit sector (5) where erosion and deposits are balanced (orange color in Figure 13). Another transit sector is the lateral debris flow (12) that increased the energy for the middle section of the slope. • Deposition sectors (green colors in Figure 13). Most of the debris was deposits in the middle sector of the slope just after the channel (6) at 1400 m a.s.l. Just below this deposit, we can find a dispersed accumulation where most of the boulders stopped (about 211,000 m 3 ) (7); the lower sector of the slope (under 1350 m) is characterized by middle energy deposits blocked by the wall (8) and partly removed by remedial work in July 2019 (8b). Part of the removed material was used to build a new section of the wall (10). On the right side of the accumulation, we can also find the collapsed boulder a" (2000 m 3 ) that stopped about 30 m from the wall (9). • Low-energy sectors (cyan colors in Figure 13). A small portion of debris flow overpassed the wall for these deposits (11). We calculated the volume for the 's' portion (155 m 3 ) with the method described in Section 3.5 and showed in Figure 8B.

Discussion
In this paper, we analyzed the evolution of the slope located between the Enni and Bich streams in the Ponte Val Formazza municipality. This area has been affected by different gravitational processes in the last decade. In particular, the rockslide that occurred in 2009 created a large chaotic deposit mobilized during the intense rainfall event in 2019 that caused a large debris flow. Different remote sensing solutions (satellite and RPAS datasets) supported a detailed analysis of the evolution that characterized the slope evolution. The particular condition of the studied area after the event occurred in 2019 makes the presented case study a good example of the strategy that could be adopted for the acquisition of a high-resolution dataset in a dangerous environment. After the event, the direct accessibility to a large part of the slope was limited by unstable blocks, and the RPAS survey has been planned considering these limitations. The availability of an RPAS equipped with double-frequency GNSS receivers allowed us to test and compare different post-processing techniques.
In particular, we found that the direct photogrammetry solution, if not supported by precise PCPs measurements, is less accurate than the mixed one, but the accuracies are acceptable with minimum support of only 3 GCPs. The internal calibration parameters of the camera are: (i) the focal distance f (more correctly called "main distance"), (ii) the position of the main point within the frame (cx, cy in Table 10), and (iii) a variable number of parameters of a polynomial indicated by Ki, through which the radial distortion is modeled. Other distortion parameters can also be calculated, which are not essential for the explanation here. The result of the calculation of these parameters is visible in Table 10 for cases 3-c and 3-a.
From the comparison of the two sets of parameters, we can draw some considerations. The precision with which they are estimated is very similar, except for the principal distance, which is more affected by the different configuration of the constraint points. In the case of complete support, it goes from 0.31 µm to 0.11µm, and the two values of the principal distance are significantly different. The differences in the position of the principal point are within 0.26 µm (about 1/9 of a pixel). Even calculating the radial distortion with the Ki coefficients in the two cases shows that the difference between the configuration with full and minimum support of GCP does not exceed 0.27 µm (about 1/9 of a pixel).
Instead, we can observe a significant variation ∆f = 15 µm of the principal distance. It involves a change in altitude in altitude ∆H = n∆f = 17 cm, where n indicates the inverse of the frame's average scale. On the ground level (Table 10), a variation ∆f implies a planimetric variation of 2.6 cm, comparable to the GSD value.
Therefore, the number and configuration of GCPs are helpful to improve the estimation of the reference system scale, which is the main effect due to the variation of the principal distance. Thus, a single GCP is not sufficient, but at least three are required for a sufficiently correct estimate of the focal length.
The calculation scenarios described were also used to obtain various DSMs, which are the most important product obtained from dense photogrammetry. Therefore, the real comparison concerns the differences found in the landslide volumes obtained in the various calculation cases already discussed. The most accurate reference used is the digital surface model obtained from Case 3-c. The precision reached by Case 3-c could be used to detect the slight ground variation and estimate the volume and debris removal cost.
For the geomorphological study of such large events, the accuracy and precision reached by the photogrammetry techniques and SfM are very good both for DSM creation and volume estimation.
The use of RPAS data associated with the SfM technique has shown many advantages and some limitations for the study of this complex and large debris flow.
The availability of pre-and post-event RPAS data allowed reconstructing the elevation and volume variations for most of the affected area.
The very high resolution obtained from orthophoto and DSM allowed the study of the debris flow deposits like in a field survey, without safety risk (area with unstable blocks and affected by rockfall). Moreover, the 3D view and DSM change were useful to understand the dynamics of some blocks that created a path with bounce and rolling.
On the other hand, the main limitation of this RPAS survey is the limited area that could be covered (mainly for legal limitation to their distance and flight height). These regulations limited the survey of the upper part of the slope (>500 m of vertical distance) and vertical cliff.
The gap in acquisition could be covered by a flight made from the upper part for the slope, but it would have required a long time hiking to reach the point of RPAS deployment. A terrestrial laser scanner (TLS) or helicopter LIDAR to supply this lack of data needs a budget far higher than RPAS surveys. In our case study, we used the satellite data of Sentinel-2 and the comparison of ground photos taken from the bottom slope to fill the gap of RPAS flight.
Focusing on the event of 10-11 June 2019, it is not ordinary debris flow for its dynamics but a more complicated process that involved several mechanisms and kinematic not so commonly found in the literature. The rockslide deposits in the channel were triggered by water flow, and then it started to move under gravitational kinematics with a roll and bounce movement. These also triggered a rockfall of a block (a' > a") with its independent trajectory. Finally, the lower part of the mass moved to show a more flow-based mechanism also implemented by a liquid contribution of Bich and other small creeks. The wall projected to protect the village of Ponte Formazza from another rockslide failed to contain this type of landslide. The debris flow deposit reached the same height of the wall (about 10 m) and, in some points, overpassed it. The older structures built to protect from debris flow (e.g., after the 1987 event) are several orders of volume below the capacity to retain the volume (about 50,000 m 3 ) of material that reached the wall in June 2019.

Conclusions
We tested the potentiality of advanced photogrammetry techniques and structure from motion (SfM) based on RPAS data in this work. We applied these techniques to study the debris flow that affected Ponte Formazza village (NW Alps) in June 2019. The collected data supported a better definition of the predisposal, preparatory, and trigger causes of the occurred event.
The meteorological and Sentinel-2 data analysis shows that the debris flow was triggered by an intense rainfall event (the registered peak of precipitation was 25 mm/h) that occurred during the most intense snow melting period of the year. The massive amount of debris and boulders inside a steep gorge is derived from a rockslide deposit. The rockslide occurred in 2009 and its deposit can be considered the most crucial predisposal factor of this event.
After the debris flow, we made a RPAS survey and we used the acquired material for a detailed mapping of the complex debris flow and rockfall that occurred in June 2019 in upper Formazza Valley (NW Alps). The availability of two surveys done in 2017 and 2018 allowed a detailed comparison of pre-and post-event topography.
The event affected an area of 0.4 km 2 , with an elevation range of 600 m representing a technical limit for many RPAS flights. Unfortunately, these limitations did not allow us to map the highest altitude sector study areas. The low-resolution mapping made with Sentinel-2 data and ground photos allowed us to supply the area of debris flow (7%) not reached by RPASS.
The post-event photogrammetric RPAS survey allowed us to make a precise and rapid mapping of hazard areas and test different photogrammetric approaches. We found that a mixed approach, based on the direct photogrammetry (PCPs) and at least 3 GCPs, allowed us to reach a centimetric accuracy almost the same as traditional photogrammetry that uses a dense network of GCPs. This approach also allowed us to obtain precise orthophotos and DSM volume variations (2.4% of difference compared to traditional photogrammetry). RPAS data processed with SfM allowed us to obtain very-high resolution orthophoto and DSMs for most of the debris flow, with a satisfactory accuracy (5 cm RMSE). This result is significant for the use of RPAS in emergency conditions when the use of GCPs could often represent a critical element of the survey.
The obtained data supported the precise mapping of 2019 debris flow deposits and helped us to understand the dynamics, volume and trajectory followed by massive boulders (up to 2000 m 3 ). The good precision of pre-and post-DSM allowed us to precisely estimate the volume changes (about 300,000 m 3 ), including the anthropic restoration that occurred between the event and the flight (July 2019). Using the very high-resolution data obtained with mixed photogrammetry methods, we also calculated the volume of the low part of the debris flow note covered by pre-event DSM by applying a clipping and interpolation.
We used the 3D views created with GIS to make a virtual field trip on the deposits in a potentially dangerous area and measure the most massive boulders included some blocks of the 2009 rockslide.
These results confirm that the RPAS surveys could be an accurate, safe, and affordable solution to map and quantify changes in small areas and limited relief energy affected by geological processes. The areas with a high elevation range, however, still need a different mapping methodology.
Simultaneously, the integration of RPAS data with low-resolution images and ancillary data is still fundamental to understand the kinematic of the processes in a multidisciplinary approach.