GLACIER VOLUME CHANGE MONITORING FROM UAV OBSERVATIONS: ISSUES AND POTENTIALS OF STATE-OF-THE-ART TECHNIQUES

: Alpine glaciers play a key role in our society through the production of freshwater for domestic, industrial and agricultural use. As they are severely affected by climate change, it is of crucial importance to understand their behaviour and monitor their morphological evolution, with the primary aims to estimate ice volume and mass changes. However, the accurate retrieval of glacier morphology changes over time is not an easy task. In this context, the use of Unmanned Aerial Vehicles (UAVs) is of interest to the glaciological community because of their flexibility, fine spatial detail and ease of processing with state-of-the-art software packages, which makes them an ideal candidate to investigate glacier changes. The goal of this work is to assess the accuracy that can be obtained with UAVs observations when comparing volume changes computed from multi-temporal acquisitions on an Alpine glacier, on the basis of a photogrammetric pipeline implemented in Leica Infinity software. The study area is Forni Glacier in Raethian Alps, Italy. Two photogrammetric blocks were acquired in 2014 and 2016 using different UAVs: a fixed-wing drone in 2014 and an in-house multicopter in 2016. Ground Control Points (GCPs) were established only during the 2016 survey which was used to establish the reference datum. Different techniques to co-register the 2014 dataset to the 2016 dataset were applied and compared: 1) using points extracted from the 2016 Dense Point Cloud (DPC) as GCPs for the 2014 DPC generation; 2) shifting and rotating the raw 2014 DPC, using manually digitised common points from the 2014 and 2016 DPCs in Leica Infinity; 3) first manually shifting, then automatically roto-translating with the Iterative Closest Point (ICP) algorithm the raw 2014 DPC in CloudCompare. The investigation shows a good agreement of the three co-registration methods in terms of height and ice volume changes and the potential of UAV data processing with Leica Infinity for glacier monitoring even when the acquisition conditions are problematic (lack of ground control points, sub-optimal image quality).


INTRODUCTION
Alpine glaciers are one of the key indicators of climate change and play an important role in our society by providing freshwater, which can be used for domestic, industrial and agricultural applications. In addition, Alpine glaciers are important accumulators of water, which marginally contribute to limit sea-level rise, while they guarantee the survival of specific mountain ecosystems. As they are severely affected by climate change, it is of crucial importance to understand their behaviour and monitor their morphological evolution, with the primary aims to estimate ice volume and mass changes. Observations show that during the 21 st century the rates of glacier retreat and mass loss have been accelerating globally (Zemp et al., 2015), and continuous worldwide monitoring of glaciers is therefore necessary to better understand their evolution and project future freshwater availability. In the Alps, retreat rates exceed 1.2% per year (Paul et al., 2020) and most glaciers are smaller than 1 km 2 (Smiraglia et al., 2015), which puts them at greater risk of disappearance. Overall, Alpine glaciers are projected to lose between approximately 35% and 90% of their area and volume by the end of the century, depending on different levels of warming in climate scenarios (Zekollari et al., 2019), with serious impacts on mountain ecosystems and the tourism industry.
Among the most important parameters to investigate the glacier state and climatic response are the area and mass balance, which relates to the mass of ice that is gained or lost during a hydrological year. The direct approach to observing mass balance through stakes drilled into the ice requires access to the glacier and is very time consuming, thus only 37 glaciers worldwide currently have a long record (> 30 years) of such direct observations (Zemp et al., 2015). The remote sensing, or "geodetic" approach, is instead based on the analysis of elevation changes of the surface, and also allows a better understanding of the glacier morphological changes. Even with remote sensing, however, the accurate retrieval of the glacier surface and its changes over time is not an easy task. With a GNSS-only survey it is complex to cover wide areas and the spatial resolution of the acquired data, being limited to the path that can be walked by a surveyor, does not provide a global insight of the studied glacial morphology. Terrestrial Laser Scanning (TLS) is affected by occlusions and may be logistically complex to operate in glacial and periglacial areas . Optical remote sensing based on satellite imagery is independent from logistic constraints, presenting however various issues when it comes to finding cloudless images due to quick changing weather conditions, as well as to the not easily tuneable temporal resolution. Synthetic Aperture Radar (SAR) images are independent from cloud coverage but often costly or unsuitable to study glacier volume variations in case of large temporal baselines. Compared to these techniques, UAV (Unmanned Aerial Vehicles) photogrammetry (Granshaw, 2018) has already demonstrated several advantages including flexibility, fine spatial detail and ease of processing with state-of-the-art software packages, which makes them an ideal candidate to investigate glacier changes. The coupling of fixed-wing or rotary-wing UAV platforms with Structure-from-Motion (SfM) (Westoby et al., 2012) and Dense Matching photogrammetric techniques probably represent the most efficient way to deal with this task today. The necessity of collecting data at multiple epochs to estimate topographic changes however introduces some challenges that make the photogrammetric approach the most suitable technique for this task. First of all, the ground control points (GCPs) may not be stable due to the difficulty of establishing permanent targets in the glacier areas, or to the limited GSM coverage in order to exploit real-time GNSS positioning services. In addition, the quality of the individual photogrammetric blocks may be uneven, due to the different drone adopted, the meteorological conditions and the flight plan (O'Connor et al., 2017;Pepe et al., 2018). The goal of this work is to assess the accuracy of Dense Point Clouds (DPCs) and Digital Surface Models (DSMs) obtained from photogrammetric blocks captured by UAVs with the aim of computing volume changes from multi-temporal acquisitions. The study area and focus of this research is Forni Glacier in Stelvio National Park (Raethian Alps, Italy). The image processing module of Leica Infinity is used to compute the image orientation, to derive DPCs, DSMs, and to compute volumes by differencing digital models.

DATA COLLECTION AND PROCESSING
In this research, we focused on Forni Glacier in Stelvio National Park (Raethian Alps, Italy). Two separate photogrammetric blocks were acquired in 2014 and 2016 using different UAVs: a fixed-wing drone in 2014 and an in-house multicopter in 2016 (see Fugazza et al., 2018 for further details). Both drones were equipped with compact digital cameras integrated with navigation-grade GNSS sensors. No GCPs were deployed during the first survey, while 12 markers were positioned and measured with GNSS in 2016, although 4 were later discarded owing to sub-optimal visibility in the images. The 2016 mission was used to establish the reference frame to develop all the analysis for volume comparison. Image blocks are processed with Leica Infinity, a geospatial surveying software focused on workflows to easily process, combine and integrate data collected from different sensors such as GNSS receivers, total stations, laser scanners and UAVs. Among other things, it offers some specific features dedicated to achieving accurate and reliable results with integrated digital data standards. Thanks to this option, for example, it is possible to assign a class of theoretical accuracy to the output DPC, which provides an easy way to analyse data. The image processing module of Leica Infinity is used to compute the image orientation (including camera self-calibration) on the basis of SfM integrated with bundle block adjustment, to derive DPCs, DSMs, and to compute volume variation by differencing DSMs. Concerning the overall photogrammetric processing, the 2016 dataset was processed using the acquired GCPs, to establish an accurate and reliable reference frame for photogrammetric block image processing. In fact, this reference frame might not be guaranteed even using surveying-grade GNSS sensors on board UAV, due to possible outages both for RTK positioning and signal acquisition. The 2014 dataset was processed using only the GNSS position of camera centres, since no GCPs had been collected due to logistic difficulties. With respect to the DPC generation, for both sets of data an upper threshold to the theoretical precision of points equal to 1 meter has been set. In fact, in Leica Infinity the quality of the generated DPC is controlled by the so-called "Filtering Threshold", which is the upper limit for the quality of the points to be generated to create the DPC; in this case 1 meter is a good trade-off between completeness and noise in the reconstruction. Three different co-registering approaches were investigated to compare the finally generated DPCs, to assess the potential of different methods which can be used when a reliable reference frame common to all the surveying epochs is not available . From each approach, a DPC, a DSM and an orthophoto were generated for the 2014 survey.
The first method (called CP2016 hereafter) consisted in digitising in Leica Infinity well visible features on the 2016 DPC to create 3D points, and marking them on the 2014 images to use as GCPs in the photogrammetric block orientation of the 2014 dataset. Five points (termed as DPC-GCPs hereafter) were manually selected outside the glacier (Figure 1).  The second method (called CONFBLOCKS hereafter) consisted in applying a conformal transformation to the DPC generated from the orientation of the 2014 dataset based only on the images positions. Leica Infinity has a function to shift, rotate and scale objects, with the option to apply a 5 or a 7 parameters transformation: for this case the 7 parameters transformation (or 3D Helmert transformation) was applied, to estimate the complete rigid-body transformation with scale. To apply this transformation, a set of common points is needed, so the same 5 feature points digitised on the 2016 DPC were identified and digitised on the 2014 DPC generated without GCPs. A 3D shift and a 3D rotation were estimated based on the common points and a roto-translation was then applied. For the specific case the scale has been kept fixed to 1: in fact, having the images a position and being the orientation computed, the expectation is that the scale is estimated when images are oriented with SfM. The 3D shift applied is the following: with standard deviations (SDs): The components of the rotation are the following: with standard deviation (SDs): Overall, only shifts in Northing and Height (for geoidal undulation) are significant with respect to their standard deviations.
The third method (called ICPDPC hereafter) is simply based on the Iterative Closest Point (ICP -see Pomerleau et al. 2013 for a review) algorithm implemented in CloudCompare (www.cloudcompare.org). The used DPCs are the same as in the previous step (DPC from 2016 generated with the acquired GCPs and DPC from 2014 generated without GCPs). As a first step, the 2014 DPC was manually shifted to roughly align to the 2016 DPC, in order to grossly remove the vertical shift of about 50 meters estimated between the two DPCs with the second method, since with the ICP algorithm only small transformations can be accurately estimated. Then the glacier was masked in both DPCs to exclude its area while performing ICP. The computed 3D shift is the following (5): with 3D standard deviations (SDs): The components of the rotation are the following: Standard deviations for rotations are not available.
It is interesting to note that, considering the respective standard deviations, the shifts estimated in CONFBLOCKS and ICPDPC approaches are equivalent.

Operational Burden of Three Co-registration Methods
One of the most important tasks in the process aimed at comparing DPCs obtained from photogrammetric datasets captured in 2014 and 2016 was related to the co-registration phase. While the UAV mission accomplished in 2016 included the deployment and measurement of some GCPs, these were not available during the 2014 campaign. The need to co-register both DPCs allowed us to check the differences in terms of final accuracy of the height comparison, as described in the previous subsection. Here we would like to analyse the adopted coregistration methods under the operational point-of-view. For the CP2016 method, it was necessary to identify and manually digitise on the 2016 DPC some well visible points in the stable regions outside the glacier. Not many well visible features are present on the glacier, so this step was the most time consuming and prone to errors. Once the DPC-GCPs were created, marking them on the 2014 dataset was easy and fast. In total the full registration process with the CP2016 method took 30 minutes. As regards the CONFBLOCKS, to estimate the transformation between the two DPCs in Infinity a set of common points was needed, so the same points digitised for the previous method were also identified on the 2014 DPC generated without using GCPs. These two sets of points were selected and used to estimate the transformation needed to co-register the DPCs. The full registration process took 40 minutes. In the case of ICPDPC method, it was necessary to carry out a selection of some stable regions outside the glacier body in order to apply the ICP registration by considering only those portions. Two regions at both rocky flanks of the glacier were selected for this purpose. The successive application of the ICP algorithm ran automatically. In the end, the transformation computed on stable areas was applied to the whole point clouds. In total, the full registration process by the ICPDPC method took 30 minutes.

RESULTS AND DISCUSSION
The 2014 DPC, originally geolocated with the GNSS image positions only, was registered to the 2016 DPC and georeferenced into its reference frame, following the three different approaches described in Section 2. This led to the generation of three different DPCs obtained from the 2014 flight (CP2016, CONFBLOCKS, ICPDPC). All these DPCs have been compared to the 2016 DPC to evaluate the glacier morphological changes.

Height Comparison Maps
The height differences between the 2016 DPC and the three registered versions of the 2014 DPCs were computed with the "Comparison Map" tool available in Leica Infinity. In all the comparisons, the DPC from 2016 was subtracted from each DPC generated from the 2014 photogrammetric dataset, which was used as reference to evaluate the evolution and changes of the glacier body.    As a general comment, the three comparison maps show a good level of matching, and broadly agree in the general pattern of topographic change on the glacier surface, although the first approach (CP2016) shows slightly larger differences (-15 m) on the southern and south-eastern parts of the glacier tongue. All methods capture well the areas with the largest changes, which exceed -20 m. These are located at the glacier terminus, southwestern and central part of the glacier tongue and correspond to the collapse of ice cavities caused by subglacial water, as well as the collapse of a portion of the glacier medial moraine owing to the melt-out of the ice core.
For a more detailed comparison between the three methods used for registration, a focus on smaller areas was carried out. In particular, the area outside the glacier (Figure 1) and Area 3 inside the glacier ( Figure 5) were analysed. In Figures 6, 7, 8, the height differences for the three DPCs for the area outside the glacier are reported. The comparison outside the glacier reveals some differences between the three methods: in general, areas outside the glacier do not show large differences but are mostly close to zero, with the exception of a small area to the hydrographic left of the glacier, subject to sediment reworking by a proglacial stream (red color in Figure 6, Figure 7 and Figure 8). The first approach applied for co-registration (CP2016 method) performed slightly worse than the others, showing evidence of positive and negative biases in the DPC comparison outside the glacier: positive biases occur in the western side of the glacier, while negative biases are located in the eastern side ( Figure 6). The other two approaches are in closer agreement (Figure 7 and Figure 8), with values generally very close to 0.   The good agreement seen on the glacier surface is confirmed by analysis of focus areas (Figure 9, Figure 10 and Figure 11); here, the CONFBLOCKS approach appeared to detect slightly smaller differences compared to the other two methods, but the general pattern is well detected by all three techniques, particularly in the areas with the largest differences (red color in Figure 9, Figure 10 and Figure 11), corresponding to the collapse of the terminal portion of the glacier medial moraine.

Volume Comparisons
Volume was computed on some focus areas on the glacier for the 3 DPCs generated from the 2014 dataset and for the DPC generated from the 2016 dataset. The areas of interest are reported in Figure 5; Area 1 and 4 represent collapsing glacier cavities, Area 2 a North-South transect along the glacier eastern medial moraine while Area 3 is the terminal part of the medial moraine also subject to collapse. Volume differences between three meshes created from the 2014 DPC, used as reference, and the mesh from the 2016 DPC were computed. Meshes were created by interpolating the DPCs. 23,05*10 3 120,15*10 3 4,540*10 3 23,74*10 3 Table 3. Size in m 2 of the four focus areas.

Cross-Sections
For a local comparison between the three 2014 DPCs and the 2016 DPC, four sections were defined inside the area of interest using CloudCompare (Figure 12). The sections were chosen to provide both longitudinal and transversal transects which cross some of the main areas of interest on the glacier, including the medial moraine and collapsing areas at the terminus. The comparisons along the longitudinal section and the three cross-sections are shown in Figure 13, Figure 14, Figure 15 and Figure 16. Furthermore, to better visualize the changes between 2014 and 2016, a focus on Section 2 was performed. All sections show good agreement between the three 2014 DPCs; differences of 15-20 m between the 2014 DPCs and the 2016 DPC were reported on the glacier surface while outside the glacier the differences are close to 0. For Section 1, differences are larger at the glacier terminus and collapsing parts of the medial moraine, gradually decreasing up valley ( Figure 13). Section 2 shows larger differences in the middle of the glacier, caused by the collapse of the medial moraine ( Figure 14 and Figure 17); Section 3 is rather homogeneous with respect to topographic changes ( Figure 15) while Section 4 reveals the largest differences on the western side of the glacier, caused by the collapse of an ice cavity ( Figure 16).

CONCLUSIONS
In this study, the DPCs obtained from images acquired by UAVs over two years were compared to estimate height and ice volume changes. The two acquisitions were captured in 2014 and 2016 using different drones and different techniques: in fact, ground control points were established only during the 2016 survey, which was used to establish the reference datum. The study area is Forni Glacier in Stelvio National Park (Raethian Alps, Italy). The photogrammetric blocks were processed with Leica Infinity. To assess the potential of different co-registering methods which can be used when a reliable reference frame common to all the surveying epochs is not available, three different coregistering approaches were investigated. The approaches analysed (CP2016, CONFBLOCKS, ICPDPC) show a good agreement in terms of morphological changes, both when height differences and volumes are compared. This study shows the potential of the photogrammetric pipeline and the analysing tools implemented in Leica Infinity software for glacier monitoring and confirms the potential of UAV data processing for glacier morphology changes assessment even when the acquisition conditions are problematic (lack of ground control points, sub-optimal image quality).