Monitoring Alpine glacier surface deformations with GB-SAR

ABSTRACT We present methods and results from interferometric data processing of a long-lasting survey campaign monitoring the Planpincieux glacier, located on the Italian side of the Mont Blanc, using a ground-based synthetic aperture radar (GB-SAR). Monitoring a European Alpine glacier during the winter, when the meteorological conditions are highly variable, presents some difficulties in radar data interpretation. The main issues to tackle in interferometric processing are unwrapping errors and high amplitude dispersion (DA), mainly due to the high velocity and dielectric heterogeneity of the backscattering surface. To improve the reliability of the results, a coherence-driven pixel-selection criterion for identifying the glacier area and a simple approach to reduce possible unwrapping errors in interferograms with low coherence are here proposed. The development of a new 2D polynomial regression model, as a function of elevation, for atmospheric phase screen (APS) estimation is also discussed. A comparison with the results obtained with a vision-based approach gave showed good agreement.


Introduction
In the past two decades, space-based synthetic aperture radar (SAR) systems have become valuable tools for monitoring glaciers when the satellite revisiting time, at best on the order of a few days, is adequate for the required temporal sampling (constrained by the process and contexts under investigation). Usually, when the survey needs shorter observation durations, as in case of short-live phenomena with significant variations on time scales as small as an hour, or in the presence of unfavourable geometries (e.g., steep or north-facing slopes), the installation of ground-based survey networks appears to unavoidable. The costs and difficulties related to the setup of onsite sensors suggest the use of terrestrial radar. The first applications on snow-covered and glacial surfaces (Luzi et al. 2007;Noferini et al. 2009;Martinez-Vazquez and Fortuny-Guasch. 2008) demonstrated the capabilities of terrestrial radar to monitor glacier and snow evolution. More recently, in Strozzi et al. (2012) and Riesen et al. (2011), the authors focused on monitoring a short-term surface ice motion linked to the fast drainage of an adjacent lake. In Allstadt et al. (2015) the results of a radar survey were used as validation of a glacier deformation model. Although this novel technique is able to provide only a partial estimate of the actual ice flow, displacement maps of glaciers of some square kilometres in size, from 3-4 km distances and entirely remotely sensed, can be obtained. Exhaustive reviews of this state of the art technique, the different instrumentation available, and the performances thereof, are available in Caduff et al. (2015) and Monserrat, Crosetto, and Luzi (2014). Especially in the case of alpine glaciers, which are particularly affected by cryosphere changes due to global warming (Deline et al. 2012), research and monitoring can largely benefit from the spatial and temporal sampling provided by terrestrial radar surveys.
In this paper, we discuss the results of a 40-day experimental campaign carried out with a commercial Ku-band terrestrial radar based on interferometric processing. During this campaign, we monitored the surface deformations of the Planpincieux glacier, which is found on the Italian side of the Mont Blanc massif.

Methods
The interferometric procedure used in this work to estimate the glacier surface deformation is based on the analysis of the differential phase between two consecutive groundbased SAR (GB-SAR) images. In the following, we briefly explain the main steps performed during data processing and present a novel approach to estimating the atmospheric phase screen (APS).
The ground-based interferometric procedure addresses some issues to obtain the differential phase due to the target displacement, φ displ . In particular, thermal and scattering noise, φ noise , phase wrapping, 2π multiple, and the phase change due to atmospheric effects, φ atm ; must be removed from the observed interferometric phase, Φ int (Monserrat, Crosetto, and Luzi 2014).
To estimate the contribution of these terms, some procedures are applied directly to the radar image. The first step of the processing chain consists of applying a 2D unwrapping algorithm (Costantini 1998) to all the images. The unwrapping is performed on pixels selected with a coherence-driven criterion (Berardino et al. 2002), with the condition γ>0:55; where γ is the temporal mean coherence. The high value of the coherence assures the statistical significance of the measured phase.
Second, the interferograms with high noise and significant unwrapping errors are excluded from the stack of the interferograms (Allstadt et al. 2015). In the case at hand, the huge number of radar images necessitates the definition of a criterion for identifying the interferograms to be rejected. The interferograms with mean spatial coherence, hγi<0:65 are excluded from the stack, and the related displacement is estimated with linear interpolation.
Third, the disturbance of the atmosphere caused by meteorological variables (i.e., humidity, temperature, and pressure, in order of relevance [Zebker, Rosen, and Hensley 1997]) is filtered from the interferometric data. The first study of the APS in GB-SAR applications (Luzi et al. 2004) found an empirical linear relationship between φ atm and range. Some studies applied polynomial (Noferini et al. 2005) or multilinear regressions (Iglesias et al. 2014;Caduff et al. 2014) to the phase observed in areas assumed to be stable, where Φ int ' φ atm , to estimate spatial maps of the APS.
Linear regression assumes homogeneous atmosphere conditions, but in our circumstance the line of sight (LOS) crosses areas with great elevation change and the assumption of a homogeneous atmosphere is no longer satisfied.
In the case at hand, the horizontal variability of the meteorological fields is considered homogeneous over distances of a few kilometres; hence, in the horizontal direction, the APS varies linearly with the range. In the vertical direction, the temperature lapse rate is linear, while humidity and pressure decrease exponentially with the height; nevertheless, it is possible to approximate with linear behaviour for limited differences in elevation (Iglesias et al. 2014). Therefore, it is necessary to introduce a term to describe the vertical atmospheric change.
According to what was just discussed, let us consider the general relationship between φ atm and range, r: Here, we described the term À r ð Þ, which represents the atmospheric disturbance as divided in a horizontally homogeneous term, hÀi ¼ const, and a term varying with the range, dÀ=dr. Considering that the atmospheric conditions are constant in the horizontal direction and vary linearly with elevation and writing dr ¼ dx þ dy þ dz, we obtain: (2) Therefore, substituting Equation (2) in Equation (1) we have: Where the scalar products z Á dx and z Á dy are null. Therefore, writing hÀir ¼ c 0 1 r; the φ atm can be estimated with the following regression of the interferometric phase over static areas (e.g., bedrock): Where a 0 ¼ c 0 is the residual offset of the integral and corresponds to the phase observed at the reference point, a 1 ¼ c 1 þ c 0 1 and a 2 ¼ c 2 . Equation (4) is similar to the APS models proposed by Noferini et al. (2005) and Iglesias et al. (2014), but here the quadratic term is a function of the elevation because horizontal and vertical atmospheric behaviours are treated independently.

Site study and dataset
The Planpincieux glacier is located at approximately 45.85°N 6.97°E in the Aosta Valley Region (North-Western Italy). The glacier lies on the southern side of the Mont Blanc massif towards the Ferret Valley, and it is part of the composite Grandes Jorasses-Planpincieux glacier at an elevation between 2530-3700 m asl, covering approximately 1 km 2 (Figure 1(a)).
As described by (Giordan et al. 2016), the lower part of the Planpincieux glacier is intensely crevassed. The morphological analysis indicates that this part is separated into two different ice flows by a central ridge of bedrock (Figure 1(b)).
The Research Institute for Hydro-geological Protection of the National Council of Research of Italy (CNR IRPI) considers the Grandes Jorasses massif an open-air laboratory for the development of monitoring systems and techniques. Since 2013, CNR IRPI is monitoring the glacier with low-cost vision-based equipment, installed in the opposite side of the valley (for further details, see Giordan et al. 2016).
Within a common research activity, CNR IRPI, the University of Pavia and the Technological Centre of Telecommunications of Catalonia (CTTC), installed a GB-SAR in the Planpincieux hamlet in August 2015 (Figure 1(a)). The positioning of the apparatus aims at maximizing the parallel condition of the LOS to the estimated direction of the western flow. This choice affected the measurement of the eastern tongue that flows diagonally w.r.t. the LOS. The mean distance to the glacier is approximately 2700 m with differences in elevation of 1200-1400 m (Figure 2). Radar images are georeferenced on a 1 m-resolution digital surface model (DSM).
During this study, an Ibis-L IDS TM instrument is used with a revisiting time of 16 minutes; the total number of processed images is 3567. The survey lasted 40 days, from 4 th September to 14 th October 2015. Figure 3(a,b) represent the mean amplitude (MA) and the DA maps, respectively. The APS estimate is applied in the area within the frames.

Meteorological data
To discover the role of the APS, we analysed the meteorological dataset from the Ferrachet automatic weather station (AWS), which was provided by the Centro Funzionale della Valle d'Aosta. The AWS was placed at 2290 m asl on the same side of the valley as the Planpincieux glacier (Figure 1(a)). The following variables were collected: air temperature, air relative humidity (RH), cumulative rainfall in 30 minutes, and snow deposited on the ground. Figure 4 reports the acquired weather data.
The period of the survey coincided with the beginning of the Alpine cold season. The temperature was close to 0 C; and some snowfall events occurred. Cold temperatures and a high RH affected the radar signal with an increased APS. To take into account the effects of the altitude difference between the AWS location and the glacier, the air temperature of the glacier should be corrected according to the Standard Atmosphere lapse rate of À 0:007 C=m; therefore, the air temperature of the glacier was approximately 3°C lower, and RH and precipitation varied accordingly.

Results and discussion
As previously mentioned, the first step of the processing consisted of applying a 2D unwrapping algorithm (Costantini 1998).
In the next step, the interferograms with hγi<0:65were excluded from the stack and the displacement occurred during the rejected data was estimated with linear interpolation of the contiguous interferograms. The temporal distribution of the rejected data (i.e., hγi<0:65Þis well correlated with adverse meteorological conditions (Figure 4) that can provoke the change of the dielectric properties due to snowdrift, snow melting and metamorphism (Luzi et al. 2009). The discarded interferograms account for approximately 16% of the total. The comparison between the 'cleaned' and 'raw' results exhibited a difference of 2-5%; therefore, probable unwrapping errors due to noise effects are strongly reduced.
The APS is evaluated in the stable areas in the surroundings of the glacier. We considered static points those pixels withDA<0:35. Figure 5(a) represents the cumulative APS in radians in the considered area, computed with Equation (4). Figure 5(b) shows the cumulative surface deformation map of the glacier. Negative data indicate points approaching the radar along the LOS, while positive data are moving away.
The deformation results in Figure 5(b) agree with the morphological observation of the glacier, highlighting two main streams with different kinematics (Giordan et al. 2016). In fact, the western flow moves downstream almost parallel to the LOS, while the eastern part flows with a remarkable horizontal component towards east. Therefore, the LOS-aligned motion component of this region is moving away.
The eastern flow presents a very low amplitude. Moreover, the radar signal passes through an atmospheric layer with strong turbulence and heat fluxes that affect the In the upper box, RH (in blue) and temperature (in red) are displayed. In the middle box, the snow deposit at ground level is drawn in grey (mm), the rainfall over 30 minutes is drawn in black (mm/30ʹ). In the lower box, the time series of the mean spatial coherence is shown.
wave phase (Iannini and Monti Guarnieri 2011;Barucci et al. 2010). For these reasons, the reliability of the results of this area is low.
Bedrock areas display almost no movement, as expected, but in the bedrock ridge on the right side (downstream) of the glacier, a sharp difference is present between the upper and the lower part. This result is probably due to unwrapping errors between non-continuous areas.
In Figure 6, time series of cumulative displacement of a set of selected points are displayed. The point locations are shown in Figure 5(b). The picture shows the stability of the two fixed points (black and grey lines) as expected, apart from a low drift that had occurred since the end of September. Points belonging to the same glacier sectors have similar behaviours. Points 1 and 2 display probable errors in the final part of the period due to their location in the front sector of the glacier, which was particularly affected by ice collapses and loss of coherence. We note a slowing down of the deformation ratio of the points in the western flow some days after the two coldest periods on 23 rd September and 29 th September. The delay is probably due to the thermal inertia of the ice.
Finally, the accuracy of the results was estimated by analysing the distribution of the interferometric phase in stable areas (Tarchi et al. 2003;Luzi et al. 2006). The mean value of the distribution is 1°, which is very close to the theoretical value of 0°, and the standard deviation converted to metric displacement is approximately 0.75 mm.

Comparison with photographic data
In the following, we compare the main radar data outcomes with some results of a contemporary monitoring campaign based on a completely independent methodology, performed by vision-based sensors (Giordan et al. 2016) during the period from 4-27 September 2015. Optical photographs estimate the surface deformation through normalized cross-correlation of consecutive co-registered images. The two systems measure different components of movement: the GB-SAR detects the component parallel to the LOS, while the vision-based system detects the component orthogonal to the LOS. The actual displacement is estimated by correcting the measured deformations with the known slope angle (Caduff et al. 2014).
In Table 1, we report the results of the comparison between the two methods, considering the mean value of the pixels within three glacier sectors of the western flow and the related standard deviation as reported in Giordan et al. (2016). The data are adjusted with a local mean slope of 45°. The results are compatible within the error interval, showing a good reliability of the GB-SAR processed data.

Conclusions
We presented the results of a survey campaign aimed at estimating the surface deformation of an Alpine glacier using a GB-SAR. A new method capable of correcting atmospheric artefacts in this specific case has been described and applied to the acquired data. The measurement campaign, which lasted forty days, was carried out within a project monitoring the lower part of the Planpincieux glacier. Critical issues in Table 1. Comparison between two different systems for evaluating the glacier surface deformation (from Giordan et al. 2016).
In the presented method, we reduced possible unwrapping errors by discarding interferograms with low coherence, and the displacement of the rejected interferograms is estimated with linear interpolation. Moreover, we developed a new APS-filtering model based on a 2D polynomial regression explicitly taking elevation into account .
The results of the data processing were compared with the results of a recent study based on the cross-correlation of optical images (Giordan et al. 2016), obtaining good agreement.
We believe that our study might be useful for developing a better methodology for GB-SAR data processing for use in monitoring the fast surface deformation processes of glaciers.