3D GLACIER MAPPING BY MEANS OF SATELLITE STEREO IMAGES: THE BELVEDERE GLACIER CASE STUDY IN THE ITALIAN ALPS

The authors group is within the Glacier Lab of Politecnico di Torino (part of the CC-LAB, a laboratory for climate change monitoring), which is working on glacier monitoring since 2016, mainly exploiting Geomatics techniques to measure the extent and to model the surface of glaciers over the years. Measurement campaigns were carried out within the ASP (Alta Scuola Politecnica – Poliecnico di Torino e Milano) DREAM projects (Drone tEchnnology for wAter resources and hydrologic hazard Monitoring) The manuscript is focused on a specific case study related to the Belvedere glacier, a valley glacier located in northern Italy.In the framework of the Belvedere glacier monitoring, several Geomatics approaches have already been applied in the last four years by the cc-glacier-lab and DREAM Projects with the goal to monitor both the extent of the glacier and its surface. Such monitoring enables the multi-temporal comparison of the glacier digital surface model (DSM), highlighting areas of ice loss and gain. Considering the limitations of aerial surveys in high altitude environments, the authors started assessing the suitability of a satellite based approach, mainly focusing on positional accuracy assessment. The paper is focused on a monitoring based on a high resolution (0.5 m) satellite optical stereo pair. Several tests were carried out with the goal to test the 3D positional accuracies, assessing the impact of different configurations of Ground Control Point (GCP) in terms of numerosity and distribution and focusing on the DSM validation. The results demonstrated the fit-for-purpose of a satellite-based approach for glacier monitoring.


INTRODUCTION
The Department of Environment, Land and Infrastructure Engineering (DIATI) at Politecnico di Torino recently founded the CC-LAB, an integrated and multisite laboratory for climate change (CC) monitoring at national and international level. The CC_LAB includes an outfield laboratory to study technologies for carbon sequestration and greenhouse effect mitigation such as green walls and green roofs (cc-green-rooflab), a moving laboratory to control and monitor climate change effects (cc-moving-lab), a cave laboratory to study paleoclimate (cc-paleo-lab) and a laboratory on a glacier to study the evolution of glacier masses (cc-glacier-lab). From 2016 several surveys and activities were carried out within the "Dream" projects (DRone tEchnnology for wAter resources and hydrologic hazard Monitoring)) involving teachers and students from the ASP (Alta Scuola Politecnica -Politecnico di Torino e Milano). The authors group are within the cc-glacier-lab, which is working on glacier monitoring since 2016, mainly exploiting Geomatics techniques to measure the extent and to model the surface of glaciers over the years, during the same season. This activity is particularly important in the framework of climate change, since a tangible effect of global warming is glacier melting, retreating and shrinking: monitoring glacier extent and volumes over the time supports global warming related analyses. The manuscript is focused on a specific case study related to the Belvedere glacier, a valley glacier located in northern Italy (Piemonte Region). In the framework of the Belvedere glacier monitoring, several Geomatics approaches have already been applied in the last four years by the PoliTO cc-glacier-lab and DREAM ASP Projects including GNSS positioning, aerial (both * Corresponding author manned and unmanned) photogrammetry, LiDAR and Total Station surveying (Avanzi, 2018, Colombero, 2019. The main goal is to monitor both the extent of the glacier and its surface to enable the multi-temporal comparison of the glacier digital surface model (DSM), highlighting areas of ice loss and gain in different time periods. As far as 3D UAV mapping is concerned, the main lessons learnt is that the effort required to survey the glacier is really high, in terms of human resources, equipment, hardware and time. Generally, a team (4-5 skilled operators, including a certified UAV pilot) deployed in the field is needed for at least three days to complete the required tasks, i.e. positioning and measuring artificial markers by means of GNSS surveying, planning the flight plans, carrying out the flight missions, downloading and storing the acquired data (data processing in office hours is not considered, being common to all the mentioned approaches). Despite the quantitative results highlighted the possibility to achieve very high 3D positional accuracies (in the range of few cm, as confirmed by De Michele, 2016), only the lower part of the glacier could be safely mapped, regardless of the use of multi-rotor or fixed wing aerial platforms. This strong limitation is due to the approximate 2000 m elevation difference between the ablative tongue of the glacier and its highest areas. For this reason, in 2019 a manned aerial survey was carried out, exploiting a small aircraft that can complete the acquisition of a photogrammetric mission (equipped with a 150 Mpx medium format Phase One iXM-RS150F photogrammetric camera, 50 mm focal length and a dual frequency GNSS receiver) in about two hours, including landing and taking off from a small airport at about 100 km distance from the glacier. Also in this case, the highest part of the glacier could not be safely mapped due to technical limitations (maximum flight height is approximately 4000 m, also due to the payload weight) of the employed small aircraft. It was therefore agreed to exploit a satellite based approach to overcome the intrinsic limitations of mountain areas mapping by means of small manned and unmanned aerial surveys. To this purpose, a high resolution satellite optical stereo pair available in the commercial archives has been acquired.

Case study
The case study is related to the Belvedere glacier, a valley glacier located in north-western Italy (Piemonte Region), as shown in Figure 1. Belvedere Glacier constitutes the terminus of the glaciers descending the steep eastern slope of Monte Rosa (4633 m) in the Italian Alps. The glacier lies at the base of the east face of Monte Rosa reaching approximately 2,200 m a.m.s.l. at its highest point. Belvedere Glacier is a typical debris-covered glacier (Haeberli, 2002). The tongue is almost completely covered with debris. The glacier virtually lacks a true accumulation basin and is fed mostly by avalanches flowing down from the eastern face of Monte Rosa, the highest face in the Alps (Diolaiuti, 2003).

Reference data
The outputs from a previous (27/10/2019) aerial photogrammetric survey were selected as reference data, specifically an aerial orthoimage and a DSM characterized by a very high horizontal and vertical accuracy of about ± 0.02 m (GSD = 0.15 m). The products have been generated in Agisoft Metashape, processing 287 images (acquired with a PhaseOne sensor) and using 19 GCPs related to artificial marker measured on the ground with GNSS receivers). The aerial orthoimagery was used for the (automatic) identification of GCP and Check Point (CP) and the extraction of their 2D coordinates. The aerial DSM was used for the extraction of the height coordinate for GCP and CP as well as for the validation of the DSM generated by means of satellite data processing.

Satellite data processing
A high resolution satellite optical stereo pair available in the commercial archives has been acquired. The in-track stereo images have been acquired by the Pleiades 1-A satellite on August, 20 th 2017 (GSD of about 0.72 m), embedding both panchromatic and multispectral (visible and near infrared) information.
The satellite data processing was carried out using the off-theshelves software package Geomatica Banff (for Education). A standard photogrammetric workflow (shown in Figure 3) was carried out, including GCP, CP and Tie Point (TP) collection, image orientation (with the calculation of the model consisting in a traditional block bundle adjustment), epipolar geometry reconstruction and generation of added value products such as orthoimagery and DSM. The Rational Function Model (RFM) was used for the model orientation, considering that RFM parameters are commonly available for most of satellite sensors and embedded in the image metadata as an additional file containing the rational polynomial coefficient estimated by the satellite data provider for that specific dataset (Tao, 2001, Boccardo, 2004. Different configurations of 3D Stereo GCPs (ranging from 0 to 15, where 0 GCP means that no adjustments are applied to the provided RFM model) were adopted, while a constant set of 17 CP was used to assess the 3D positional accuracy (details in section 3.1). The 2D position of GCPs was identified by means of an automatic procedure based on image autocorrelation (exploiting the aerial orthoimagery as reference data) and then manually checked and refined when needed (assuming a ± 0.5 pixel accuracy in point collimation). The elevation component was extracted from the aerial photogrammetric DSM. GCP were collected as Stereo GCP, i.e. the same point was identified in both images. The use of Stereo GCPs provides a stronger math model: they add redundancy and are weighted more heavily in the bundle adjustment (PCI Geomatics, 2019). The set of CPs was collected manually. As shown in Figure 3, the geometric distribution of GCPs and CPs was limited by two main factors: firstly the coverage of the aerial orthoimagery and DSM adopted as reference data (covering only the eastern part of the satellite footprints) and secondly the lack of time-invariant points over the glacier tongue and snow-covered areas.  Once the model is oriented, the stereo pair is reprojected into an epipolar pair in which the images have the same orientation, with corresponding features between the images laying on the same axis. Matching pixels are extracted throughout the employment of automatic image correlation algorithms, which calculate the 3D position of the points starting from the sensor geometry and the mathematical model (PCI Geomatics, 2019). It is thus possible to exploit the oriented model to interpolate a regular DSM starting from the denser 3D point cloud (Figure 4). The DSM can be eventually exploited as ancillary data source for the orthocorrection of one of the images of the oriented stereo pair ( Figure 5). All the outcomes have been referenced with respect to the ETRF2000 Datum (UTM projection, ellipsoidal height) using an UTM projection (Zone 32N). The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition)

RESULTS
Several model orientation tests were carried out with the goal to test the 3D positional accuracies, assessing the impact of different configurations of GCP in terms of numerosity and spatial distribution with respect to a constant set of 17 CPs. Additionally, the satellite-based DSM was validated with respect to the reference aerial DSM.
where n = n. of CP m = 2 for 2D RMSE, 3, for 3D RMSE The results are summarized in the Table 2 as well as in box plot graphics ( Figure 6 and Figure 7), where the limits of the boxes represent the first and third quartiles, the median is represented by the red line splitting the box and the whiskers display the minimum and the maximum value. Furthermore, some check points were excluded since considered outliers (red crosses).  The plot is integrated with the average RMSE (μ) and the interval [μ-σ; μ+σ], which provide a graphical indication of the accuracy and precision of the model. As shown in Figure 6, the main improvement in 2D RMSE can be noticed passing from 0 to 1 GCP, while a further increase of the number of GCP does not impact neither on the precision nor on the accuracy. The same behaviour can be also noticed in Figure 7, showing the DZ versus the number of GCP. Despite one single GCP could theoretically be the best choice in terms of effort, it is strongly recommended to use at least 3 GCPs for the model orientation step, to have redundancy for controllability and reliability and avoid the risk of gross errors when using just 1 GCP. In order to have a more realist measure of the planimetric accuracy of the final orthoimagery product, the 17 CPs were manually identified on the orthoimages generated with both the aerial DSM and the satellite DSM. A total of 16 satellite orthoimages (one for each GCP configuration shown in Table 2 using both satellite images) and 8 DSMs (one for each GCP configuration) have been generated and evaluated. The results are shown in Figure 8 and Figure 9. The results confirm that the orthoimage planimetric accuracy can be improved and stabilized with just 1 single GCP, while the impact of the input DSM is mainly on the precision (smaller standard deviations when using the more accurate aerial DEM) rather than on the accuracy. The results confirm the possibility to generate accurate orthoimagery using DSM derived by satellite stereopairs as elevation source.

DSM Validation
The satellite DSM was therefore compared to the reference aerial DSM, with the main goal to assess the DSM elevation accuracy The evaluation of the vertical precision was carried out calculating the elevation differences between the satellite DSM (generated with model oriented with 3 GCPs, according to outcomes of the model accuracy assessment) and the aerial DSM.
Considering that the glacier surface is continuously changing over the time, the comparison needs to be limited to areas considered stable over the years, mainly areas around the glacier ablative tongue. For this reason also vegetated and permanent snow-covered areas were excluded (Aguilar, 2012); the latter represents an additional limitation in the comparison of the DSMs, allowing the assessment of DSM elevation accuracy to be carried out only in areas at lower elevation ( Figure 10). As for the planimetric accuracy assessment, the elevation DSM differences are shown as a box plot and are charted with respect to the number of GCPs (Poli, 2014). The differences (Figure 11) have been evaluated over more than 21 million points randomly sampled. In line with the outcomes of the planimetric accuracy, it can be noticed the main enhancement in term of accuracy and precision can be achieved with just 1 single GCP.  Table 2 To evaluate the possible relation between slope and vertical accuracy (as highlighted in previous studies as Dowman et al., 2012), the DSM accuracy versus the local slope of the terrain was mapped and analysed. In particular, 5 categories with different slope ranges were used, as shown in Table 3. Table 3. Slope ranges adopted for the assessment of the DSM accuracy vs slope The results were reported in a chart representing the empirical cumulative distribution function of |ΔZ|. For each |ΔZ|, the value of the function is calculated as: where n = n. of points where |ΔZ| has been evaluated The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 12 -Empirical cumulative distribution function of |DZ| for different categories of slope (m) calculated over more than 21M points As shown in Figure 12, the main decrease in DSM accuracy can be noticed for the two categories with slope values greater than 70%. More specifically, in areas with a limited slope (< 30%) 95% of |Z| are lower than 1.5 m, while this value is in the range 2-3 m for slopes > 70%.
In order to perform also a qualitative evaluation of the satellite DSM and to identify possible outliers, elevation profiles have been plotted and compared. Two examples are provided in Figure  13, with respect to the transects shown in Figure 10. It can be highlighted an overall overlapping of the aerial and satellite elevation profiles, except for the part related to the glacier ablative tongue where high height differences are present. Considering the 2 years' time gap among the aerial and satellite acquisitions, it was indeed expected a loss of glacier volume (Donizetti, 2017). Lastly, it is also possible to generate realistic 3D scenes for qualitative assessment (e.g. to spot errors in the DSM deriving by the interpolation of areas with a low density of 3D points or with outliers) and dissemination purposes, as in the example shown in Figure 14.

Multi-temporal analysis
As highlighted, the availability of a multi-temporal dataset is posing some limitations in the accuracy assessment stage. Nevertheless, multi-temporal data enable the possibility of 4D glacier monitoring. To this purpose pixel-by-pixel DSM differences in the time period 2017-2019 have been calculated, providing information for glacier evolution (including volumes) monitoring (and for a possible global warming monitoring). Figure 15 shows the elevation DSM differences limited to the area covered by the glacier ablative tongue. Additionally, the comparison of elevation profiles extracted from DSMs related to different years allows surface changes to be clearly detected and measured ( Figure 13).

CONCLUSIONS
The main focus of the paper was on the assessment of the planimetric and vertical accuracies of DSM and orthoimagery derived by satellite stereo pairs (for glacier monitoring purposes).
One of the main outcome is that satellite stereopair based analyses allows the limitations of aerial surveys in mountain regions to be overcome. The whole glacier can be mapped (including areas that can't be surveyed with UAV or with small aircraft with suitable accuracies and a very high level of detail (GSD = 0.5 m). Specifically, 2D positional accuracy lower than 1 m and vertical accuracy lower than 1.5 m (with standard deviations lower than 1 m) can be achieved with one single GCP when suitable reference data are available. Nevertheless, it is recommended to use 3 GCPs for a more robust solution, especially to detect potential gross errors. A qualitative analysis of elevation profiles highlights a very good match between reference and satellite DSM, clearly showing areas of glacier volume loss. Nevertheless the DSM should be properly inspected to detect potential outliers in areas characterized by poor image matching performances. The aforementioned outcomes confirms the possibility to successfully exploit satellite remote sensing for 4D glacier monitoring when proper reference data are available for the identification of GCPs (e.g. an aerial orthoimagery and related DSM). Indeed, a satellite-based monitoring is nowadays characterized by a high level of detail (maximum commercial GSD = 0.3 m as of May 2020), a high revisiting time (theoretically up to 1 image/day, mainly limited by cloud coverage conditions) and with 3D accuracies lower than 2 m even without any GCP. The accuracy assessment was mainly limited by the lack of reference data at higher altitudes: a field survey exploiting as base camp the hut "Capanna Regina Margherita" (the highest building in Europe at 4,554 m a.m.s.l) is already planned. Further developments will be focused on: i) validating the whole satellite dataset with additional reference dataset (e.g. regional LiDAR or DTM from National Mapping and Cadastre authorities) ii) using different satellite sensors to verify if the proposed methodology can be successfully generalized and iii) check the fit for purpose of a satellite-based direct georeferencing approach, especially in areas where no reference data are available or can be properly measured.