UAV survey method to Monitoring and analysing geological hazards: The Case study of mud volcano of Villaggio Santa Barbara, Caltanissetta (Sicily)

All active geological processes determine effects on the soil due to different deformation processes: surface uplift and subsidence, shear lineaments with differential kinematics in relation to the source and the soils involved. Among all the active geological processes on Santa Barbara mud volcano (Caltanissetta town, Italy), represents a dangerous site because it caused, on 11 August 2008, a paroxysmal event, which determined severe damages to the infrastructures at around to 2 km the paroxysmal event. The remote sensing of surface deformation now represents a key tool for the evaluation and monitoring of the hazard. The use of unmanned aerial vehicles (UAVs) in contexts of natural danger presents three main steps for risk assessment and monitoring: pre-post event data acquisition, emergency support and monitoring. Here we present a methodology for monitoring deformation processes that may be precursors of paroxysmal events on the Santa Barbara mud volcano. Among the precursors, the lifting and development of structural features are the most important, with dimensions ranging from centimetre to decimetre. Therefore in relation to the magnitudes of the phenomena involved, the objective of this work is (going from the acquisition phase, to the SfM processing chain and the use of the M3C2-PM algorithm for the comparison between point clouds and uncertainty analysis with a statistical approach) the monitoring of deformation processes, with centimetre precision and a temporal frequency of 1 2 months, as precursor indications of hazard.

Applications include runoff laboratory trials (Morgan et al., 2017), applied geology (Niethammer et al., 2012;Russell, 2016;2 Saito et al., 2018), geomorphology (Bakker and Lane, 2016;Balaguer-Puig et al., 2017;Dietrich, 2015;Dietrich, 2016 a/b;Heindel et al., 2018;Javernick et al., 2014;Marteau et al., 2016;Mercer and Westbrook, 2016;Pearson et al., 2017;Prosdocimi et al., 2017;Seitz et al., 2018;Smith and Vericat, 2015;Snapir et al., 2014;Vinci et al., 2017), glaciology (Immerzeel et al., 2017;Piermattei et al., 2016), coastal morphology (Casella et al., 2016;Brunier et al., 2016;James and Robson, 2012), volcanology (Andaru and Rau, 2019; Bretar et al., 2013;Carr et al., 2018;De Beni et al. 2019;Favalli et al., 35 2018;Giordan et al., 2017Giordan et al., , 2018James and Robson, 2012;Müller et al. 2017;Witt et al. 2018 ) and geophysics (Amici et al., 2013a;Di Felice et al., 2018;Federico et al., 2019;Greco et al., 2016;Zahorec et al., 2018). SfM is commonly used in the cultural heritage field for the 3D reconstruction (Jalandoni et al., 2018;Sapirstein, 2016Sapirstein, , 2018Sapirstein and Murray, 2017). The monitoring of active geological processes is a preventive action in risk mitigation (Deng et al. 2019;Diefenbach et al., 2018, Rosa et al., 2018Stöcker et al., 2017;Turner et al., 2017b). Disasters occur when two factors -hazard and 40 vulnerability -coincide. The risk is proportional to the magnitude of the hazards and the vulnerability of the affected population. Among the deformation monitoring systems, the photogrammetry technique from Unmanned Aerial Vehicles (UAVs) is spreading thanks to the high efficiency in data acquisition, the low cost compared to traditional techniques and in the acquisition of high resolution images (Fonstad et al., 2013;Harwin and Lucieer, 2012;James and Robson, 2012;James et al., 2017/a/b;James et al., 2020;Javernick et al., 2013;Johnson et al., 2014;Westoby et al., 2012). This technique is 45 important to study the catastrophic natural events such as floods, earthquakes, landslides, subsidence, etc. Different acquisition methods and the ability to obtain high spatial and temporal resolutions (Boccardo et al. 2015) allow to acquire detailed information on the evolution of the landscape, therefore UAVs are an effective and complementary tool for field investigations. Furthermore, UAVs have other advantages including: (i) the ability to fly at low altitudes, (ii) the ability to reach remote locations, (iii) the ability to host multiple sensors (cameras, Lidar, thermal imaging cameras, navigation / 50 inertial sensors, etc.), (iv) the ability to capture images at different angles, and (v) the flexibility to carry out monitoring operations on a small, medium and large scale (Jordan et al. 2017). Ground Control Points (GCPs) are used to improve the accuracy of the resulting data. Therefore, recognisable points on the UAV imagery are measured with a high-precision surveying device to georeferenced the data. In this process a correct number of GCPs is required which lead to a greater accuracy of the resulting data (point clouds, 3D grid, orthomosaic or Digital Surface Model (DSM)). The precision of the 55 resulting data is also controlled by other variables, such as: the focal distance of the camera; flight path and flight altitude; the orientation of the camera; the picture quality; the processing chain and the category of UAV system (fixed or rotary wings).
In this paper, we present the results and analysis of the surface deformation monitoring of the mud volcano of Santa Barbara (Caltanissetta, central Sicily) (Fig.1). We have applied the statistical analysis of significant changes with 95% Level of 60 Detection (LoD 95%). In detail, we used precision maps and the M3C2-PM (Lague et al., 2013;James et al., 2017b) algorithm to determine the surface variations. The statistical analysis allows to verify i) the uncertainty between the different surveys, ii) the spatial variability of the accuracy in the surveys (James et al., 2017/b), iii) the quality of the georeferencing of the surveys based on the number of GCPs. The mud volcano of Santa Barbara is located within the Caltanissetta foredeep basin of the Apennines-Maghrebian collisional chain which developed from the Late Miocene to the Quaternary, along the border of the converging Eurasia-70 Nubia plate (Catalano et al., 2008;Dewey et al., 1989;Serpelloni et al., 2007). This structural domain is formed by a foreland fold and a thrust belt involving the deposition of clastic sediments which were gradually contracted and moved from the late Miocene to the Pleistocene (Monaco and Tortorici, 1996;Lickorish et al., 1999, and references therein).
According to Madonia et al., (2011) mud volcanoes are most of the time in stasis, but they are a preferential way for rising fluids rich in methane and sludge, therefore they can be considered a risk to urbanized areas or sites with an economy 75 dedicated to natural attractions.
On the mud volcanoes (e.g. Ayaz-Akhtarma and Khara Zira Island mud volcanoes in Indonesia), geomorphic/structural features have been observed from one year to shortly before the paroxysmal event (Antonelli et al., 2014;Madonia et al., 2011). Also the area of the Santa Barbara volcano was affected in 2008 by paroxysmal mud eruption which was preceded by deformation features (Fig.2). Moreover, the surface of the mud volcanic cone is incised by a drainage system ( Fig.1) 80 characterised by hydrographic basins with elongated dendritic geometry arranged to a centrifugal development from the areas of the summit craters towards the lower slopes of the volcano complex. The higher order of the thalwegs present deep recessed meanders (landscape rejuvenation process). This morphometric structure is typical of uplifting areas and therefore relative decrease of the base level. This suggests an inflection process of the volcano ground surface induced by increase of fluid pressure inside of shallower stagnation chamber which is located at a depth of about 30 m . The 85 stagnation chamber has a "sill-like" geometry, a radius of about 50 m and a thickness of about 30 m. This morphostructural configuration supported by geophysical data configures the active geological structure as a high potential geological hazard. On the surface, fractures and shear lineaments extend outside the areas of the volcano, (Madonia et al., 2011;Bonini et al. 2012;INGV, 2008a;Regione Siciliana, 2008) and they highlight the high stress and strain environment induced by the mud 95 volcano (Fig.2). Such structures have been detected in the past and we speculate that they are still active. This development has often been a precursor of paroxysmal events such as the 11th August 2008 event (INGV, 2008a;Regione Siciliana, 2008).

Local Network
In order to monitor active deformation in the mud volcano area, a local GNSS network was created according to the criteria According to these criteria two GNSS benchmarks were created: CTN0 and CTN1, located on the roof of a building on the northern sector of the studied area (Fig.3). The benchmarks were surveyed using double frequency (L1/L2) receivers (TOPCON Hiper V and HiPer SR) in static mode.
Once the stability of the benchmarks has been assessed, we surveyed CTN0 and CTN1 5 and 2 times respectively. 115 Post processing of GNSS data was carried out by AUSPOS Online service (Geoscience Australia, 2011;Jia, M et al., 2014). Finally, the optimal ITRF2014-UTM33N coordinates have been definitively assigned to CTN0 and CTN1. 125

Ground Control Points (GCPs)
Various acquisition methods of Ground Control Points (GCPs) were tested in order to define the most suitable one. The main function of GCPs is to geo-reference outcomes from SfM.

Initially (2016-2018) we used only GNSS receivers in different configurations: Real Time Kinematics (RTK) and Static 130
Ultra Rapid. The GPCs were made of 50 cm x 50 cm alveolar polypropylene square targets. Using these configurations (Tab. 1), errors ranging from centimetre to decimetre were recorded. In these early phases, errors were only computed by SfM software PhotoScan (v 1.4.5.7554).
PhotoScan provides different types of error estimation: XY error (m) -root mean square error for horizontal coordinates for a GCP location; Z error (m) -error for elevation coordinate for a GCP location; Error X, Y and Z (m) -root mean square 135 error for X, Y, Z coordinates for a GCP location; Error Img (pix) -root mean square error for X, Y coordinates on an image for a GCP location averaged over all the images; Total Error (m) -implies averaging over all the GCP locations.
Using Static Ultra Rapid mode, the total error was about 18 cm, whereas using the RTK configuration the total error were reduced up to about 4 cm. .
From 2018 the use of the TopCon DS-103 TST was introduced, obtaining lower values on the GCP errors than those 140 detected with the GNSS technique (Tab.1). With this measurement technique the total error has been reduced to about 3 cm.
On the first 3 campaigns, we used to georeferenced the cloud points 6 GCPs (Fig.4). In the following section, these first 3 campaigns will be not considered in the computation of the significant changes for their incomparable to the last campaigns.
Since 2019 only TST has been used for the GCP survey and according to Tahar et al. (2013) the number of GCPs has been increased (Tab.1). Considering these two improvements, the total error has been reduced to about 1.4 cm and in the last 145 campaign about 0.7 cm.

155
In the last two campaigns we used TST to obtain the coordinates of the GCPs. The TST was positioned on the CTN0 point of the local network, (Fig.3) which coincides with the roof of the nearby private houses in the northern part ( Fig.5 A). A classic celerimetric survey was carried out.
The measurements of the GCPs were carried out with a surveying ranging rod equipped with a reflecting prism (offset of -30 160 mm) assisted by a tripod with a bubble level (Fig.5 B), to ensure the verticality and stability of the measurement.

165
Based on the CTN0 point and directing on the CTN1 point, we assumed an instrumental error of 5 mm (named in PhotoScan "Marker Accuracy"). This value has been assumed due to the uncertainty of the CTN0 and CTN1 point coordinates, obtained through GNSS measurements, considering that the uncertainty derived from the TST is negligible.
To validate GCPs data we have done an analysis (compared with the SfM software data) using the python script 170 "Monte_Carlo_BA.py", with a statistical iterative approach (Monte Carlo approach) (James et al. 2017/a). For clarity, the following terminology will be used in the text: To be more precise values ranging from 10% up to 80% of GCPs have been set.

Photo Acquisition
We have performed five measurement campaigns in approximately one year. The same flight plan was used for all five of them. The AOI (Area Of Interest) was captured by a DJI Phantom 4 Standard, a quadcopter UAV, at a flight height of 33 m above the ground. The sensor size of the UAV's digital camera is 6.17 mm by 4.55 mm, capable of shooting images with a resolution of 12 MP (4000 × 3000 pixels) with a mechanical shutter. Each flight planning was carried out with the Pix4D 185 Mapper software, adopting a frontal and side overlapping of 80% and 70% respectively. The camera was set-up in a nadiral orientation to capture vertical imagery. The flight was carried out in a single grid (simple geometric flight patterns without intersections). An average of 280 images for each survey were acquired with about 1.1 cm Ground Sampling Distance (GSD).

Data Processing SfM 190
The photogrammetric processing was performed using the commercial software Agisoft PhotoScan. The photogrammetric processing is based on the workflow formulated by the USGS (2017)   The procedure above is conventional, but we modified same parameters during the cleaning procedure of sparse point cloud.
We adopted the optimal values of "Reconstruction Uncertain" and "Projection Accuracy" in the "Gradual Selection" option.
Finally, after that the model has been georeferenced, we adopted the optimal value of the "Reprojection Error" (USGS, 2017). 200 The first parameter, "Reconstruction Uncertain", improves the geometric reconstruction of the cloud point.
The second one, the "Projection Accuracy", improves pixel mismatches in images.
After these operations, the GCPs can be tied to the sparse point cloud.
After the georeferencing procedure, we have to perform the "Gradual Selection" option, adjusting the "Reprojection Error" After that, we applied the "precision_estimate.py" script (James, et al. 2017 / b). The script operates, by an interactive Monte Carlo approach, to estimate the accuracy of SfM surveys through photogrammetric and georeferencing parameters, which are then used to provide spatially variable confidence limits for the detection of surface variations. 210 The precision estimates are calculated through multiple "bundle adjustments" ("Optimize Camera Alignment" in PhotoScan) with different pseudo-random offsets (in this case 4000 pseudo-random offsets) (Fig.7)  The 3D topographic change is usually detected from sparse cloud points that have been cleaned to exclude the vegetation that interferes with the comparison techniques. Next step is to link those points to the precision estimates of the sparse cloud; 230 this has been done in CloudCompare v. 2.11. Through CloudCompare the sparse cloud is interpolated (with the relative precision values on the three-dimensional components) with the dense cloud point. In this phase we decide which interpolation technique is the more suitable, in this case the best is "nearest neighbors" that consider the three closest points to the sparse cloud, using the median value (better outliner mitigation) to assign the error values to the dense cloud points (Fig.9). This methodology has been chosen due to the 235 heterogeneous distribution of the points in the sparse cloud, avoiding point with null values. • significant change (Fig.10 A). 250 • M3C2 distance (Fig.10 B).
The first output (Fig.10 A) shows the changes which exceeds the uncertainty values in both point clouds. It represents a confidence interval constrained by values of Level of Detection 95% (LoD95%) which are spatially variable. This is applicable in any morphological setting, providing a reliable 3D analysis of topographical change. 255 The second output (Fig.10 B) shows the calculated distances between the two clouds.
The third output (Fig.10 C) shows the uncertainty values of the distances and their spatial variation between the two clouds.
Once the uncertainty value is defined, the changes are significant when they overtake the value of the uncertainty. In the 2020/06/15 survey, we used callipers to test and validate the method. Five numbered callipers, with steps of increasing 265 height of about 2 cm, were positioned on the mud dome. The heights are 2, 4, 6, 8 and 10 cm.
These were used to obtain an instrumental sensitivity of the measures (Fig.11). All callipers are detected, we show as example the smallest (Fig.11).

Result
As illustrate on the section 2.2, the results show that using a percentage of GCPs (considered as control points) between 40 and 60, the RMSE value has a slight minimal variation. Thus, the optimal minimum number of GPCs is between 12 and 18.
Exceeded the threshold of 60% there are no significant improvements (Fig.12). This result has been confirmed in all campaigns. 280 The results obtained from the photogrammetric comparisons supported by geodetic topographic survey, have an average uncertainty of about 6.4 cm, with a minimum of about 2 cm and a maximum of about 12 cm, relative to an area of 42700 m² 290 ( Fig.13 A). The uncertainty of the central area has an average value of about 3.9 cm, with a minimum of about 2 cm and maximum of 10 cm, on an extension of 360 m² (Fig.13 B).

cm. (B) the Gaussian distribution of the distance uncertainty of the central emission area is shown, the mean is about 3.9 cm and the standard deviation is about 1.4 cm. The reference system is WGS 84 / UTM zone 33N.
During the monitoring period, in addition to natural changes, we recorded anthropogenic action, due to the relocation of objects (garbage) in the area after unauthorized access. These changes have a decimetric order of magnitude and are easily 300 detected by the technique used (Fig.14). Two types of analysis were carried out: semi-quantitative and quantitative. The surveys 2019/07/29 and 2020/06/15 (time interval of about 1 year) were chosen to perform the semi-quantitative analysis. This analysis was carried out on the whole of 315 the mud volcano area, the objective is to detect important deformations (in the order of decimetres). To visualized important deformation (Fig.15), the values of M3C2 distance ranged between -2 and 2 cm were excluded, according to the minimum value of distance uncertainty. This range was verified instrumentally with the use of callipers. Considering the 2019/07/29 and 2020/06/15 surveys (Fig. 15) we observe that the total surface of the volcanic cone can be 325 considered not affected by deformations. Only morphological changes in the volcanic structures can be underlined, such as small eruptive cones, griphons, sauces and mud pools.
The quantitative analysis was performed on small central portions of the mud volcano. The aim of the quantitative analysis is to estimate the trend and evolution of the deformation. To assess the deformation and the local morpho-structural evolution, two temporal series have been developed in two areas: Z1 (collapse zone) and Z2 (uplifting zone) ( Fig.16 and Fig.17). The 330 zones were chosen related to the significant change. On the selected zone, the average distances between points with significant change were computed by M3C2-PM. Z1 has an extension of 1 m² (Fig.16) and Z2 has an extension of 0.84 m² (Fig.17).
As master was chosen the campaign of 19/07/29 (T0) (   In figure 16 the subsidence (about 10 cm / 60 days) of gryphon is represented by the southwest migration of the mud pool 345 edge. In figure 17, the uplift (about 4 cm / 60 days) of a blind gryphon is represented by the development of radial fractures on the cone surface. Close to the growing gryphon, another indication of deformation is the deviation of a mudslide flowing in the N-S direction towards the southern sector of the analysed area. Finally, the sudden appearance of the gryphon is well highlighted in the ortophoto of the last survey (Fig. 17).

355
The reference system is WGS 84 / UTM zone 33N.

Discussion
UAV technology combined with SfM is a valuable tool in geological risk assessment and monitoring, however, some questions must be considered. 360 The first one is the quantity and quality of the acquired GCPs. According to the results, the minimum optimal number for geo-referencing is 12 GCPs. A slight improvement is observable up to 18 GCPs, while no appraisable improvement is detected for a higher number of GCPs (Fig.12). In addition, the method of acquiring GCPs reduces their error by using the Total Station Theodolite (Tab.1). The combined use of high-precision topographic instruments with an optimal number of GCPs improves the reliability of the datasets. 365 The second question is the evaluation of distance uncertainty when two surveys are compared. The distance uncertainty between two data sets (surveys) can be considered as an estimate of the sensitivity of the methods to detect measurable topographic changes. The results of the M3C2-PM show that the average distance uncertainty between the first and the last surveys is about 6.5 cm over the whole area at the 95% confidence level (Fig.13 A). Furthermore, considering a smaller portion of the area, the uncertainty decreases to about 3.9 cm at the 95% confidence level (Fig.13 B). This allows us to 370 analyse certain morphological changes and anthropogenic activity on the volcano's shell (Fig.14). The anthropogenic activity determined a height decrease of about 16 cm where a wood platform was previously located. As well as, a height increase of about 14 cm is recorded at the place where the wood platform has been relocated. The values recorded are consistent with the real thickness of the object detected, differing by 2 cm (Fig.14). The callipers show that it is possible to measure changes at least 2 cm (Fig.11). 375 After that the uncertainty and sensitivity of the surveys were computed, the time series were made. Considering the significant changes and their errors, we computed the trend of deformation of two areas in order to reconstruct the evolution of the phenomena which generated these changes ( Fig.16 and Fig.17).
Considering the results obtained by SfM, we propose a monitoring system of Santa Barbara mud volcano based on a semiquantitative approach. We defined three states by setting a space-time range: 380 • Normal state if the deformation value does not exceed 5 cm in the range between 1 to 2 months.
• Pre-Alert state if the deformation value is between 5 cm and 10 cm during a period of between 1 to 2 months.
• Alert state if the deformation value exceeds 10 cm in the range between 1 to 2 months.
Obviously, these values must be considered in relation to the surface involved by deformation. To define the state of the activity of the mud volcano, the area affected by deformation must be a significant area of the total surface of the volcano. 385 The definition of the activity is an important aspect of the monitoring of hazard, because there are many natural phenomena that can produce deformation in a little area (like a gryphon) or in a big area (like the sedimentation occurred at the border of the mud volcano), so the significant area affected by deformation must be correctly defined by an operator (public agency, university, institute of research or Dipartimento della Protezione Civile) which really knows the natural phenomena that can occur on the Santa Barbara mud volcano. 390 The methodological approach is valuable and efficient from the point of view of quantity and quality of data collected in relation to the work and time spent. This monitoring system of the deformations is a useful tool to detect the early unrest phase of the mud volcano usually induced by changes in pressure and volume of fluid raising from the stagnation chamber.

Conclusion
In the management of the hazard, the SfM technique, (Gomez et al., 2016, Kaab, 2000, Fugazza et al., 2018Giordan et al., 2017Giordan et al., , 2018 starts to be largely used by the scientific community. In the monitoring of potentially dangerous active sites, the UAVs are very advantageous because they are not used only as a support in the post-disaster events (Rokhmana and Andaru, 400 2017;Hisbaron et al., 2018), but also for the pre-event monitoring.
According to Kopf (2002), Antonelli et al., (2014), Madonia et al., (2011), INGV (2008a and Regione Siciliana (2008), the deformations of the surface shell of mud volcanoes can occur to a year before the paroxysmal event with doming and developing of structural lineaments with order of magnitude from centimetres to decimetres.
The results allow us to define the criteria for monitoring and analysing the study area. For the mud volcano of Santa Barbara, 405 the monitoring criteria are: • Monitoring interval between 1 and 2 months.
• The optimal number of GCPs is between 12 and 18.
• The processing chain of the sparse point cloud according to workflow of USGS (2017) was enhanced by the correct 410 value of "Tie Point Accuracy" and "Marker Accuracy" as suggested by .
• An assessment of the state of activity of the mud volcano based on a semi-quantitative approach.
The frequencies of the campaigns depend on the status of activity, while the other criteria depend on the object / structure of the monitoring.
These criteria allow us to detect events with deformation of at least 2 cm. In case of anomalous values were detected the 415 monitoring campaigns must be improve. This involves extensive monitoring, such as: i) developing time series localized in key areas, ii) combining different methodologies, e.g. micro-seismicity monitoring and three-dimensional geophysical prospecting (Imposa et al., 2016) to improve the monitoring system of the active geological process.
at Sinabung Volcano, Sumatra, Indonesia, documented by structure-from-motion photogrammetry, Journal of Volcanology