Next Article in Journal
Multi-Temporal and Spectral Analysis of High-Resolution Hyperspectral Airborne Imagery for Precision Agriculture: Assessment of Wheat Grain Yield and Grain Protein Content
Next Article in Special Issue
East Africa Rainfall Trends and Variability 1983–2015 Using Three Long-Term Satellite Products
Previous Article in Journal
Real-Time Precise Point Positioning Using Tomographic Wet Refractivity Fields
Previous Article in Special Issue
A Multisensor Approach to Global Retrievals of Land Surface Albedo
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Intercomparison and Validation of SAR-Based Ice Velocity Measurement Techniques within the Greenland Ice Sheet CCI Project

by John Peter Merryman Boncori 1,*, Morten Langer Andersen 2, Jørgen Dall 1, Anders Kusk 1, Martijn Kamstra 3, Signe Bech Andersen 2, Noa Bechor 4, Suzanne Bevan 5, Christian Bignami 6, Noel Gourmelen 7, Ian Joughin 8, Hyung-Sup Jung 9, Adrian Luckman 5, Jeremie Mouginot 10,11, Julia Neelmeijer 12, Eric Rignot 10, Kilian Scharrer 13, Thomas Nagler 13, Bernd Scheuchl 10 and Tazio Strozzi 14
1
DTU Space Institute, Technical University of Denmark (DTU), DK-2800 Kongens Lyngby, Denmark
2
Geological Survey of Denmark and Greenland (GEUS), Department of Glaciology and Climate, DK-1350 Copenhagen, Denmark
3
Science & Technology AS (S&T), N-0160 Oslo, Norway
4
Earth Resources Laboratory, Massachussets Institute of Technology (MIT), Boston, MA 02139, USA
5
Department of Geography, College of Science, Swansea University, Swansea SA28PP, UK
6
Istituto Nazionale di Geofisica e Vulcanologia (INGV), Centro Nazionale Terremoti, 00143 Rome, Italy
7
Institute of Geography, University of Edinburgh, Edinburgh EH89XP, UK
8
Applied Physics Laboratory, University of Washington, Seattle, WA 98105, USA
9
Department of Geoinformatics, University of Seoul, 02504 Seoul, Korea
10
Department of Earth System Science, University of California Irvine, Irvine, CA 92967, USA
11
CNRS, IRD, Grenoble INP, IGE, Université Grenoble Alpes, 38000 Grenoble, France
12
German Research Centre for Geosciences (GFZ Potsdam), 14473 Potsdam, Germany
13
Environmental Earth Observation IT GmbH (ENVEO), A-6020 Innsbruck, Austria
14
Gamma Remote Sensing AG, CH-3073 Gümligen, Switzerland
*
Author to whom correspondence should be addressed.
Remote Sens. 2018, 10(6), 929; https://doi.org/10.3390/rs10060929
Submission received: 20 April 2018 / Revised: 15 May 2018 / Accepted: 11 June 2018 / Published: 12 June 2018
(This article belongs to the Special Issue Remote Sensing of Essential Climate Variables and Their Applications)

Abstract

:
Ice velocity is one of the products associated with the Ice Sheets Essential Climate Variable. This paper describes the intercomparison and validation of ice-velocity measurements carried out by several international research groups within the European Space Agency Greenland Ice Sheet Climate Change Initiative project, based on space-borne Synthetic Aperture Radar (SAR) data. The goal of this activity was to survey the best SAR-based measurement and error characterization approaches currently in practice. To this end, four experiments were carried out, related to different processing techniques and scenarios, namely differential SAR interferometry, multi aperture SAR interferometry and offset-tracking of incoherent as well as of partially-coherent data. For each task, participants were provided with common datasets covering areas located on the Greenland ice-sheet margin and asked to provide mean velocity maps, quality characterization and a description of processing algorithms and parameters. The results were then intercompared and validated against GPS data, revealing in several cases significant differences in terms of coverage and accuracy. The algorithmic steps and parameters influencing the coverage, accuracy and spatial resolution of the measurements are discussed in detail for each technique, as well as the consistency between quality parameters and validation results. This allows several recommendations to be formulated, in particular concerning procedures which can reduce the impact of analyst decisions, and which are often found to be the cause of sub-optimal algorithm performance.

Graphical Abstract

1. Introduction

While the feasibility of measuring ice-motion from space-borne optical and Synthetic Aperture Radar (SAR) imagery has long been established ([1,2], and references therein), the routine generation of ice velocity maps of Greenland and Antarctica has become feasible only recently [3,4,5], due to several factors: the unprecedented availability of regularly acquired free and open access data, provided in particular by the European Space Agency (ESA) Sentinel-1a and 1b SARs and the U.S. Geological Survey (USGS)/National Aeronautics and Space Administration (NASA) Landsat-8 mission; the maturity of the data processing techniques; the increased performance and affordability of storage and computing resources; the establishment of long-term Earth Observations programs, such as the ESA Climate Change Initiative (CCI) and the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs).
The currently available ice-sheet velocity products consist of continental-scale seasonal mean-velocity maps and time-series of the horizontal ice-motion components for selected areas. Generation of the latter based on SAR data, which are the focus of this study, entails post-processing and mosaicking measurements of satellite line-of-sight (LoS) and flight-path (azimuth) velocity components, derived independently for a large number of radar tracks.
This paper describes an activity carried out within the ESA Greenland Ice Sheet CCI project, referred to as the “ice velocity round-robin algorithm intercomparison”, in which several international research groups were invited to take part in four experimental tasks, related to different measurement techniques/scenarios. Participants were provided with common SAR datasets, and asked to apply their processing algorithm of choice to derive ice-velocity measurements, which were then intercompared and validated against GPS data by the project consortium.
The measurement techniques considered in this study are reviewed in Section 2.1, whereas the SAR and the validation datasets are described in Section 2.2 and Section 2.3, respectively. In Section 3, the results provided by the participants, their inter-comparison and their validation are discussed for each experimental task. The impact of several algorithmic choices on the derived ice-velocity products is discussed in Section 4, and recommendations based on the best-performing approaches are provided. Conclusions are drawn in Section 5.

2. Methods

2.1. Sythetic Aperture Radar Ice-Motion Measurement Techniques

In this study, three core processing algorithms for the measurement of the SAR LoS and/or azimuth motion component are considered, namely Differential SAR Interferometry (DInSAR), offset-tracking, and Multi Aperture Interferometry (MAI).
Differential SAR Interferometry yields measurements of LoS displacement between two SAR acquisitions, exploiting the proportionality between LoS-distance (slant-range) differences and phase-delay differences of corresponding (co-registered) image pixels, provided the radar returns at the times of the two acquisitions are statistically similar (coherent). The contributions of topography and displacement can be separated either by using a single image pair and an external Digital Elevation Model (DEM), an approach known as DEM Elimination (DEME) [6], or by using two SAR image pairs, a method referred to in the following as double-difference (DD) [7]. The measurement accuracy provided by DInSAR is a fraction of the centimetric SAR wavelength, provided interferometric coherence between the two acquisitions is not lost, which can occur in polar environments due to weather-induced surface changes and/or high deformation gradients. The achievable spatial resolution of the measurements in each image dimension is in the order of several radar resolution cells, which in turn are between 2 m and 20 m in each image dimension for SAR acquisition modes relevant for ice-sheet studies. From the processing point of view, DInSAR requires image co-registration and 2D phase unwrapping, which are time-consuming and sometimes delicate steps for non-stationary and/or partially coherent scenes, particularly for imagery acquired in Terrain Observation by Progressive Scans (TOPS) acquisition mode ([8,9]).
Offset-tracking methods measure both the LoS and the azimuth components of the displacement between two SAR acquisitions by estimating the mis-registration of corresponding points on ground due to motion. This is achieved by tracking incoherent amplitude or intensity features (incoherent offset-tracking) and/or coherent SAR speckle (coherent offset-tracking) in several ways. Amplitude or intensity cross-correlation (ACC and ICC in the following) [10], as well as translation estimation in the Fourier domain via phase correlation [11], are able to exploit both the presence of features and of correlated speckle, typically yielding results on the featureless, but possibly coherent, ice-sheet interior as well as on the crevassed, but often incoherent, glacier tongues. Improvements in terms of accuracy and spatial resolution can in some cases be obtained using coherence maximization [12] or complex cross-correlation (CCC) [13,14], although these methods are only applicable in practice to areas/datasets in which a high level of phase coherence is retained. The achievable accuracy and spatial resolution of all offset-tracking techniques is typically an order of magnitude worse than DInSAR. From the processing point of view, although these methods are computationally intensive, they are also straightforward to parallelize and can be automated to a high degree, since image co-registration and most importantly phase unwrapping are not required.
Multi Aperture Interferometry [15] measures the azimuth displacement component between two SAR acquisitions, and can be proven to be an alternative formulation of the azimuth spectral-diversity technique originally proposed for image co-registration [16]. The same measurement principle, also referred to as split-bandwidth interferometry [14], can also be applied in the range dimension, where it was first proposed for absolute phase estimation [17]. To our best knowledge, however, the range variant of this method has never been applied to measure ice-motion. Split-bandwidth techniques exploit the fact that image mis-registrations in the azimuth (or LoS) dimension are necessarily coupled with phase variations, which are proportional to the registration offsets and to the central Doppler (or range) frequency of the processed data bandwidth. Such phase variations can thus be estimated by forming interferograms between upper and lower azimuth (or range) spectral sub-bands, from which displacement are then computed by differencing and scaling. From the performance point of view, these methods achieve a higher spatial resolution and accuracy compared to amplitude or intensity offset-tracking for a given level of coherence ([14,18]), and are more robust and computationally efficient than complex cross-correlation or coherence maximization. Although split-bandwidth techniques share the same phase-coherence requirements as DInSAR, and also require image co-registration, it is only for large displacements (comparable to the size of the SAR resolution cell) that phase unwrapping is required, so that this processing step can often be neglected.
Finally, for all of the above-mentioned techniques, it can be necessary to estimate and remove slowly-varying image-wide error trends due to satellite orbit and radar system uncertainties, and to the relative nature of the measurements in the DInSAR case [19,20]. This is typically done by estimating the parameters of suitable error models based on points of known position and velocity, referred to in the following as Ground Control Points (GCPs).

2.2. Experiment Description

The ice velocity round-robin algorithm inter-comparison experiment was organized into four different tasks, namely DInSAR (Task 1), MAI (Task 2), incoherent offset-tracking (Task 3), and partially-coherent offset-tracking (Task 4). Invitations to participate in one or more tasks were issued by the Greenland Ice Sheet CCI consortium to academic, research and business organizations, through announcements on the project website and on the CRYOLIST e-mailing list for snow and ice related research, as well as through personal communications. For each task participants were provided with the same input SAR and auxiliary datasets (Table 1 and Table 2) and were asked to apply their preferred implementation of one of the three processing algorithms described in Section 2.1. A set of deliverables was requested for each task (Table 2) and participants were requested to fill in a feedback form providing information on their processing algorithm and parameters.
The Differential SAR Interferometry dataset (Task 1) consisted of two European Remote Sensing (ERS)-Tandem pairs (I1 and I2 in Table 1) covering an area of about 200 km × 100 km (two standard frames) around the Storstrømmen glacier in northeastern Greenland (Figure 1b). Precise state vectors from Delft University of Technology and the 90 m posting Greenland Mapping Project (GIMP) DEM distributed by the Byrd Polar and Climate Research center [21] were also provided. Both SAR image pairs are highly coherent due to the small spatio-temporal baselines (Figure 2a,b), except for an unphysical band of lower coherence for the I2 pair, due to missing lines in the raw data. Participants were asked to generate a LoS velocity map using one of two applicable DInSAR techniques, namely DEME and DD, and to provide quality parameters if applicable as well as the position of any GCPs used in the processing.
The Multi Aperture Interferometry dataset (Task 2) consisted of the previously mentioned I2 pair in Table 1. Users were provided with the Level-0 (raw) data products and with the same ancillary data as for Task 1. Participants were asked to generate an azimuth velocity map and supply any applicable quality parameters and GCP positions.
For the offset-tracking experiments (Tasks 3 and 4) image pairs I3 and I4 in Table 1 were provided. The I3 dataset consisted of an Environmental Satellite Advanced Synthetic Aperture Radar (ENVISAT ASAR) C-band image pair spanning two standard frames (about 200 km × 100 km) around the Upernavik Isstrøm glaciers in northwestern Greenland (Figure 1c). The Synthetic Aperture Radar raw data were made available as well as Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) precise state vectors and the GIMP DEM. Due to the C-band radar frequency and the relatively long temporal baseline, the I3 pair has very low coherence values on the ice sheet (Figure 2c). The I4 dataset consisted of an Advanced Land Observing Satellite 1 Phased Array L-band Synthetic Aperture Radar (ALOS PALSAR) Fine Beam Single image pair, spanning 5 standard frames and covering an area of about 300 km × 70 km, featuring several large tidal outlet glaciers, such as Rink Isbræ, Store Gletscher and Sermeq Avannarleq (Figure 1d). The Synthetic Aperture Radar raw data and the GIMP DEM were provided to the participants. Due to the longer radar wavelength, the I4 pair shows several patches of relatively high coherence values, both on the bedrock outcrops and on slowly-moving areas on the ice sheet (Figure 2d). Coherence is, however, also lost for this pair closer to the ice margin, and in particular on the tongues of the above-mentioned glaciers. For both offset-tracking tasks, participants were asked to generate LoS and azimuth velocity maps using their offset-tracking method of choice, and to provide any applicable quality parameters and the position of any GCPs used during the processing.

2.3. Validation Strategy

The SAR-based measurements provided by the participants were validated against in situ GPS measurements available from different sources. GPS coordinates and Cartesian velocity components are provided in Table S1, whereas their locations are shown in Figure 1 and Figure 2. Neither the GPS measurements nor their locations were made available to the participants beforehand.
For Tasks 1 and 2, velocities were derived from the positions of 10 stakes drilled into the ice surface, which were surveyed at least twice between 1992 and 1995 [22]. Since the last observations were in late June 1995, and thus between 6 and 8 months before the acquisition of the I1 and I2 SAR datasets (Table 1), the measurements were linearly extrapolated to a reference date (1 February 1996), to account for the post-surge, long-term adjustment of the ice dynamics, as discussed in [22].
For Task 3, three measurement sources were used, namely stations deployed during the Ice2sea project between August 2009 and July 2011 (UPE2, UPE3, UPE4) [23], Automatic Weather Stations deployed in August 2009 within the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) project (UPE_U, UPE_L) [24], and a GPS network deployed during a field campaign in summer 2010 by DTU (UPF1-3, UPM1-3).
For Task 4, GPS stations AVA1, RNK2, STO4 were used, which were deployed respectively on Sermeq Avannarleq, Rink Isbræ and Store Gletscher within the Ice2sea project between August 2009 and August 2010 [23].
To compare the GPS measurements, recorded in a local ENU Cartesian coordinate system (ve, vn, vu), with the LoS and azimuth velocity components provided by the SAR-based methods, the GPS velocities were projected onto the LoS and azimuth directions according to the following transformation:
( v r v a ) = ( cos φ cos θ sin φ cos θ sin θ sin φ cos φ 0 ) ( v e v n v u )
in which vr represents the SAR LoS velocity component (positive for motion towards the radar), va the SAR azimuth velocity component (positive for motion in the satellite flight path direction), φ the angle between the East direction and the ground projection of the SAR LoS vector, and θ the elevation angle from the horizontal plane to the SAR LoS vector. The values of these angles at scene center are listed in Table 1.

3. Results

In the following subsections, the results received for each task are described and intercompared. Contributions are anonymized and participants are referred to through group numbers, which are randomly assigned for each experiment.

3.1. Differential Synthetic Aperture Radar Interferometry (Task 1)

3.1.1. Measurements

For this task, results were provided by 6 groups, affiliated (in alphabetical order) to the Geological Survey of Denmark and Greenland (GEUS), the Italian National Institute for Geophysics and Volcanology (INGV) and the Technical University of Denmark (DTU). Some institutions provided more than one result, obtained with different processing algorithms. Image focusing was carried out by all groups with the software of [25], whereas the interferometric processing chains, which were used were based on different combinations of the software modules developed by [25] and on DTU’s IPP processor [26]. The LoS mean-velocity maps generated by the participants are shown in Figure 3.
Group 1 applied a DEME approach to the I1 dataset (Table 1), whereas groups 2, 3, and 5 applied the same approach to I2. Groups 4 and 6 processed both I1 and I2 using a DD algorithm. The main processing parameters are detailed in Table 3.
Co-registration was carried out by all groups according to a 2-step procedure, initiated by geometric calculations exploiting the available precise state-vectors and subsequently refined by fitting a mis-registration model to the results of an intensity image patch cross-correlation carried out on a coarse grid. Algorithms differed in the number of cross-correlation windows, their size, and in the order of the mis-registration polynomial model (Table 3).
Interferometric phase and coherence were estimated with similar methods, the main differences lying in the averaging factors (number of looks) and in the use of adaptive phase filtering [27], which was carried out by groups 1, 2, 3, and 6.
Concerning phase unwrapping group 1 used a branch-cut algorithm [28], whereas all other groups used different implementations of the Minimum Cost Flow (MCF) algorithm [29]. Weights for the latter were based on coherence for groups 2, 3 and 6 and on both coherence and intensity edge-detection for groups 4 and 5. All groups used a coherence threshold to generate a phase unwrapping mask, with thresholds selected either manually (group 2) or as a function of the number of interferogram looks (all other groups).
Absolute phase estimation and orbital refinement were carried out by all groups using GCPs manually selected on bedrock outcrops. The number and spatial distribution of the GCPs, however, varied significantly (Figure 4). In particular, group 2 used a single GCP to estimate only a constant phase calibration parameter, whereas all other groups estimated four parameters describing a linear baseline error model [19].
Quality parameters were provided by all groups (Figure 4). For groups 1 and 3 to 6, they consisted of maps of predicted LoS error standard deviation, estimated according to the framework of [30], whereas group 2 provided a coherence map.

3.1.2. Intercomparison

The received results were resampled to a common grid with a posting of 3 arcsec in latitude and 12 arcsec in longitude, corresponding to about 90 m × 90 m on ground. For comparison purposes, the differences of each LoS velocity map with respect to that provided by group 6 were computed (Figure 5). The measurements provided by group 1 (Figure 3a) show a more limited coverage compared to other groups, as well as significantly higher velocity magnitudes on ice and on bedrock outcrops, whereas the results provided by the remaining groups show local differences, in particular in the bottom-left and bottom-right corners of the image, as well as slightly varying image-wide trends (Figure 5b–e).
The differences observed for group 1 are most likely due to the more conservative branch-cut phase unwrapping algorithm, concerning coverage (i.e., the percentage of valid measurements), and to the selection of several GCPs in non-stationary areas (Figure 4a) as well as in areas affected by phase unwrapping errors, which cause a spatially varying bias in the measurements.
Concerning the results obtained by groups 2 to 6, some discontinuities are seen in the difference maps, e.g., in the lower-right of Figure 5b,d,e, which correspond to regions almost entirely delimited by areas of low-coherence in Figure 2a,b. These regions are prone to phase unwrapping errors, particularly if network approaches such as MCF are used [31], which could explain the observed differences. Results generated with a DEME approach applied to data pair I2 (Table 1), namely by groups 2, 3 and 5 (Figure 5b–d), show some common difference patterns with respect to the DD results based on the I1 and I2 pairs, obtained by groups 4 and 6. This could be due to the height estimation procedures, based respectively on an external DEM in the DEME case and on the joint inversion of height and displacement from two SAR data pairs in the DD case, and/or to the different sensitivity to GCP height errors, which affect only the DEME results [32].
Concerning the quality measurements, those provided by groups 3 to 6 are similar, except for the fact that the error predictions associated with the DD methods (Figure 4d,f) are slightly lower than those associated with the DEME ones (Figure 4c,e), due to the insensitivity to GCP height errors. Interestingly, the error predictions for group 1 are an order of magnitude higher than other groups. This is probably due to the fact that the phase unwrapping error modelling of Mohr et al. ([30]) identifies several GCPs to be in areas which are likely to be affected by single cycle phase unwrapping errors.

3.1.3. Validation

Measurement accuracy was assessed on exposed bedrock, identified using the ice-cover mask of Citterio et al. ([33]), and on ice-covered surfaces by comparison with the GPS positions of 10 stakes drilled into the ice surface [22]. The LoS-projected GPS velocities were up to 48 m/y in magnitude. Summary statistics are reported in Table 4, whereas the differences for each stake are detailed in (Table S2).
Table 4 confirms that the results of group 1 are affected by a significant image-wide error trend, whereas those of all other groups are similar and show a level of agreement with the GPS measurements, which is typically of 1 m/y or better in terms of median and Median Absolute Deviation (MAD). Interestingly, the group which achieves the lowest deviation with respect to the GPS-surveyed stakes, namely group 2, shows a slightly higher variability on bedrock areas compared to groups 3 to 6. This is probably due to the use of a single GCP for orbital refinement (Figure 4b), which leads to larger errors further away from this point, where some of the bedrock areas are located.
Finally, it is interesting to analyze to what extent the quality parameters provided by the participants predict the validation results. For groups 3 to 6, indeed, most observed errors fall within two predicted error standard deviations (Figure 4 and Table S2), whereas for group 1, although the predicted error standard deviations are higher by an order of magnitude, they still consistently underestimate the observed errors. This is due to the fact that the error prediction algorithm used by group 1 can account for errors in GCP velocity only if an appropriate uncertainty value is specified by the analyst. Furthermore, only single-cycle phase unwrapping errors are assumed to occur [30]. Concerning group 2, which provided a coherence map as a quality measurement, the values of this parameter are uniformly high and thus insensitive to the observed errors (Table S2).

3.2. Multi Aperture Interferometry (Task 2)

3.2.1. Measurements

For this task, the I2 dataset in Table 1 was processed by 4 groups, applying different processing chains listed in alphabetical order: DTU [34], Massachussets Institute of Technology [15], University of Edinburgh [18], and the University of Seoul [35]. Image focusing was carried out with the software of Werner et al. ([25]) and that of Rosen et al. ([36]). The resulting azimuth mean-velocity maps are shown in Figure 6.
Multi-aperture processing was carried out by all groups by creating two azimuth sub-apertures for each acquisition, with frequency separations between 0.4 and 0.5 of the SAR Pulse Repetition Frequency (Table 5). Different approaches were followed to generate the sub-apertures. Group 1 focused each image by processing the full azimuth bandwidth, and subsequently carried out spectral whitening combined with azimuth common-band filtering and range spectral-shift filtering. Co-registration was performed according to the two-step procedure described in Section 3.1.1. Finally, the two azimuth sub-apertures were formed by azimuth band-pass filtering. In contrast, groups 2 to 4 carried out the focusing step twice for each raw dataset, selecting different azimuth center frequencies and processing bandwidths to jointly carry out azimuth common-band filtering and sub-aperture formation. The corresponding sub-aperture images for each acquisition were then co-registered following a similar two-step procedure as for group 1.
All groups performed low-pass filtering and decimation (multi-looking), followed by a local mean filtering in the case of groups 2 and 3 and by a Goldstein adaptive filter [27] in the case of groups 1 and 4. Group 4 carried out multi-looking and adaptive filtering on each sub-band interferogram, prior to MAI (difference) interferogram formation, whereas group 1 applied both filtering steps to the final MAI interferogram alone. An additional processing step was carried out by group 4, to reduce the co-registration errors induced by ionospheric propagation also known as “azimuth streaks” [37]. The algorithm which was applied is similar to the one proposed in Raucoules et al. ([38]) for offset-tracking and consists in masking out areas with significant motion based on a user-defined threshold, rotating the image by a user-defined angle so that most ionospheric streaks are parallel to the slant-range direction (under the assumption that most streaks are almost straight lines), and finally isolating the streaks from residual motion and other error sources with a suitably designed spatial filter.
Phase unwrapping was carried out only by group 4, using an MCF approach [29]. MAI phase was then converted to azimuth displacement after performing an orbital refinement in the case of group 4, and applying different scaling factors, based either on the nominal antenna length (group 1) or on the sub-band frequency separation and the beam footprint (ground) velocity (groups 1,2,3).
Quality parameters were provided by all groups except group 1 and consisted of predicted error standard deviations (Figure 7), based on local moving-window estimates.

3.2.2. Intercomparison

The received results were resampled to a common grid with a posting of 3 arcsec in latitude and 12 arcsec in longitude, corresponding to about 90 m × 90 m on ground, and the difference of each azimuth velocity map with respect to that provided by group 4 was computed (Figure 8). The difference maps are dominated by ionospheric streaks, since group 4 is the only group, which implemented a dedicated filtering step for their estimation and removal (Table 5). This is apparent also from Figure 6d, in which the results of group 4 show significantly lesser streak artefacts.
Other differences are related to slowly-varying trends in both the range and azimuth dimensions and to the levels of spatially uncorrelated noise. Relative biases among the measurements can be due to different co-registration error models and parameter estimations, and to the orbital refinement phase calibration step in the case of group 4. Uncorrelated noise levels are instead influenced by the respective filtering approaches and parameters. In particular, groups 1 and 4 applied a higher degree of spatial filtering compared to groups 2 and 3. Concerning groups 1 and 4, although the same filtering algorithms were applied with similar parameters (Table 5), the results of group 4 contains a higher level of short-wavelength noise (Figure S1d). This could be due to the application of the phase filtering steps on the sub-band interferograms in the case of group 4, and on the differential MAI interferogram in the case of group 1, as discussed further in Section 4.2.
Finally, the results of group 1 contain an error due to the handling of missing-lines by the focusing software (red line crossing GPS stake 6 in Figure 8a), which is also visible in the coherence map generated from the same focused image pair (Figure 2b).
Concerning the error predictions in Figure 7, these are significantly lower for groups 2 and 3, compared to group 4, even though the local variations would actually seem lower for the latter (Figure S1). We speculate that this could be due to the use of smaller estimation windows for the computation of the local standard deviations on which the error estimates are based, although information on the window sizes was not provided by all participants and thus this hypothesis cannot be confirmed.

3.2.3. Validation

Measurement accuracy was assessed on exposed bedrock, identified using the ice-cover mask of Citterio et al. ([33]), and on ice-covered surfaces, by comparison with the same GPS data used for the validation of Task 1 [22]. The azimuth-projected GPS velocities were within 112 m/y in magnitude (Table S3). Summary statistics are reported in Table 6, whereas the differences for each validation location are detailed in Table S3.
Table 6 shows that group 4 achieved the highest level of agreement with the GPS measurements, as well as the lowest velocity variations on bedrock. Based on Figure 8, the main reason for this is the reduction of ionospheric streaks errors, which is also achieved to a lesser extent by the spatial filtering approaches applied by the remaining groups.
Concerning the quality parameters, the validation results of group 4 were essentially predicted by the error estimates (Figure 7), with all but one observed error lying between two predicted error standard deviations (Table S3). Groups 2 and 3 instead systematically underestimated the observed errors, since these are mainly due to ionospheric streaks, which are not accounted for by locally-estimated standard deviations. No quality parameter was provided by group 1.

3.3. Incoherent Offset-Tracking (Task 3)

3.3.1. Measurements

For this task the ENVISAT ASAR I3 dataset in Table 1 was processed by 6 groups, applying different processing chains, listed in alphabetical order: DTU [39], Environmental Earth Observation Information Technology (ENVEO, [3]), GAMMA-RS [25], German Research Centre for Geosciences (GFZ- Potsdam, [40]), GEUS [25], Swansea University [41]. The raw SAR data were focused with the software of [25,36,41]. LoS and azimuth mean-velocity maps were provided by all participants and are shown in Figure 9 and Figure 10 respectively.
The core processing algorithm was ICC for all groups, although several differences can be pointed out, as summarized in Table 7. Prior to offset-tracking, three groups carried out image co-registration, following the two-step approach described in Section 3.1.1, whereas the remaining groups skipped this step. Group 2 also carried out multi-looking and high-pass spatial filtering of the intensity images to enhance features. Concerning offset computation, groups 1, 4, and 6 chose the same oversampling factor, cross-correlation window size and a similar spacing between adjacent windows, leading to comparable spatial resolutions of the measurements. Group 2 chose window sizes corresponding to a larger ground coverage, whereas groups 3 and 5 applied larger oversampling factors prior to the cross-correlations, and chose smaller window sizes respectively in both image dimensions (group 3) and in the range dimension only (group 5).
Regarding outlier culling, all groups applied thresholds to one or more quality parameters, namely cross-correlation Signal-to-Noise-Ratio (SNR), as defined in [41], Normalized Cross-Correlation (NCC) coefficient, and maximum expected velocity magnitude. Group 2 additionally carried out a culling based on the local variability of flow magnitude and direction. The culled offsets were smoothed (filtered) by groups 1 and 6, which applied a moving average window, and by group 3, which used anisotropic diffusion filtering [42].
Estimation of image-wide error trends (calibration) was carried out by groups 1 and 6 using GCPs selected on bedrock, respectively by visual inspection (group 1) and based on the ice mask provided in Citterio et al. ([33]) (group 6). Group 4 implicitly performed a calibration step by fitting the coefficients of a 2D polynomial to high-SNR range and azimuth offsets computed on a coarse grid, and then co-registering the SAR images according to this model prior to offset-tracking. The positions of the GCPs used by each group are shown in Figure 11.
Quality parameters were not provided by all groups and were quite heterogeneous. They consisted of per-pixel LoS and azimuth error standard deviation estimates for groups 1 and 6, NCC coefficient amplitude for groups 3 and 5, and cross-correlation SNR [41] for group 4 (Figure 11). The predicted error standard deviations provided by groups 1 and 6 were generated by applying the framework of [30] and modelling measurement errors through the local variance, scaled to account for the partial overlap of adjacent cross-correlation windows [13].

3.3.2. Intercomparison

The coverage of the measurements is the first apparent difference between the results. Due to the low interferometric coherence (Figure 2c), measurements were expected to be feasible only on feature-rich areas, namely on bedrock outcrops and close to the ice-sheet margin. Indeed group 5 provided results also towards the interior of the ice-sheet, although these are clearly noisy (Figure 9e and Figure 10e) and correspond to very low NCC values (Figure 11e). Groups 2 and 3 (Figure 9b,c and Figure 10b,c) obtained almost complementary coverages, confined to glacier tongues in the case of the former and to rock and slow-moving ice-sheet areas for the latter. For group 2, this could be due to tight settings for the additional culling step based on local flow direction variability, which can be expected to be more sensitive to noise in areas where the velocity magnitudes are low. In contrast, the limited amount of good matches found by group 3 on glacier tongues seems due to the higher NCC threshold chosen for outlier culling, combined with the fact that the NCC values are reduced by the small cross-correlation window sizes (Figure 11c). Finally, groups 1, 6 and 4 all provided results on rocky surfaces and close to the ice-sheet margin as expected, although with a decreasing measurement coverage. Based on the insets in Figure 9 and Figure 10, group 1 appears to have used a less conservative hole-filling approach compared to both groups 4 and 6, although details of these procedures are not available to confirm this. Another lacking information is the culling procedure implemented by group 4, which is based on SNR thresholding and on a second unclearly described criterion, which is likely however to be the cause of the reduced measurement coverage compared to groups 1 and 6, since the latter used more conservative SNR thresholds.
To compare the location and spatial-correlation of the differences between groups, the received results were resampled to a common grid with a posting of 9 arcsec in latitude and 27 arcsec in longitude, corresponding to about 250 m × 250 m on ground. In Figure 12 and Figure 13 the differences between the LoS and azimuth velocities compared to group 6 are shown for the area in the insets in Figure 9 and Figure 10.
Differences are consistently higher for the azimuth measurements, probably due to the fact that for this dataset the highest velocity gradients, and hence challenges for the measurement algorithms, occur in this direction. Several spatially correlated differences are found in shear zones and more in general in areas adjacent to data voids. The fact that this holds also for participants, such as groups 1 and 6, whose measurements have the same spatial resolution and level of filtering, and also similar outlier culling thresholds (Table 7), suggests that this is mostly due to different hole-filling procedures. Spatially uncorrelated differences, such as those of group 5 compared to group 6, are instead more likely to be due to lower outlier culling thresholds and a lesser degree of spatial filtering.
Concerning the quality parameters (Figure 11), the predicted error standard deviations provided by groups 1 and 6 are similar, and mostly reflect the local spatial variability of the velocity measurements, which are higher for the azimuth component compared to the LoS one. The parameters provided by groups 3, 4, and 5 characterize the quality of the cross-correlation peaks, which are influenced by the strength (signal-to-clutter-ratio) of the intensity features and/or by the interferometric coherence.

3.3.3. Validation

Measurement accuracy was assessed on exposed bedrock, identified using the ice- mask of [33], and on ice-covered surfaces, by comparison with the GPS measurements. Summary statistics are reported in Table 8, whereas the differences for each GPS station are detailed in Table S4.
For all groups, the velocity errors with respect to the GPS measurements are larger in the azimuth component compared to the LoS one, except for group 5, whose statistics are significantly affected by outliers. In contrast, the LoS and azimuth variations on bedrock are comparable and, except for group 2, they are an order of magnitude less than the errors observed at the GPS locations. Thus, the worse agreement of the SAR azimuth measurements cannot be due to ionospheric propagation errors, i.e., the azimuth streaks discussed in Section 3.2, which would equally affect measurements on bedrock and at the GPS locations. A possible explanation could be that the limited (kilometric) spatial resolution of the SAR measurements (Table 7) has a greater impact on the accuracy of the azimuth component, since the velocity gradients are higher in this direction compared to the LoS.
Concerning the differences discussed in Section 3.2.1, a velocity profile is shown in Figure 14, which crosses the location of 4 GPS stations along the dashed line in the insets of Figure 9 and Figure 10. The negative differences between group 1 and most other groups in this area appear due to errors in the measurements of group 1, caused by excessive hole-filling as previously discussed. The higher levels of spatially uncorrelated noise of groups 4 and 5 compared to other groups are also apparent in Figure 14, as well as the unexpectedly high azimuth variability of group 2, which was also observed on bedrock as mentioned above, and the reasons of which are unclear.
Concerning the agreement between validation results and quality parameters, for groups 1 and 6, almost all observed errors fall within two predicted error standard deviations (Table S4). This is consistent with the previously mentioned hypothesis that the main error source for this dataset is the spatial variability of the underlying velocity field, particularly in the azimuth direction. The error predictions do in fact account for spatial variability, although just for the one between neighboring measurements, rather than for that within the coverage of the cross-correlation windows. For the groups which provided maps of NCC coefficient amplitude or cross-correlation SNR as quality parameters, there is not a clear relation between these and the observed errors (Table S4).

3.4. Partially-Coherent Offset-Tracking (Task 4)

3.4.1. Measurements

For this task, the ALOS-1 PALSAR I4 dataset in Table 1 was processed by 8 groups, applying different processing chains listed in alphabetical order: DTU ([25,39]), ENVEO [3], GAMMA-RS [25], GEUS [25], Swansea University [41], University of California, Irvine [43], and University of Washington [13]. The raw SAR data were focused with the software of Werner et al. ([25]) and with that of Rosen et al. ([36]).
The LoS and azimuth mean-velocity maps provided by the participants are shown in Figure 15 and Figure 16 respectively. Since group 7 provided the horizontal velocity components along the Cartesian axes of a National Snow and Ice Data Center (NSIDC) Polar Stereographic North projection (EPSG 3413), the LoS and azimuth components were derived for this group by using the surface parallel flow assumption and the GIMP DEM [21]. Based on a comparison with optical imagery, significant relative geocoding discrepancies were observed among groups, ranging between +/−300 m on ground, so that constant 2D shifts were estimated and applied to inter-compare and validate the results.
The participants applied several different offset-tracking techniques, of which the main processing steps and parameters are summarized in Table 9. Groups 1, 2, 4, 5, 6, and 8 all used ICC as a core calculation method, whereas Group 3 relied on ACC and Group 7 applied a combination of CCC and ACC.
Quality parameters, shown in Figure 17, consisted of per-pixel LoS and azimuth error standard deviation estimates for groups 1, 6 and 8, the amplitude of the NCC coefficient for group 5, and the cross-correlation SNR [41] for group 4. Groups 2 and 3 did not provide quality estimates, whereas LoS and azimuth quality parameters could not be derived for group 7, based on the information which was made available. The predicted error standard deviations provided by groups 1, 6 and 8 were generated by applying the framework of [30], adapted for offset-tracking as described in Section 3.2.1.

3.4.2. Intercomparison

The coverage of the measurements differs in slow-moving areas towards the interior of the ice-sheet, in fast-moving areas close to the glacier termini and in terms of the overall azimuth extent.
Towards the interior of the ice-sheet some results show significantly more data voids than others (e.g., Figure 15b,d,f). Several processing steps and parameter trade-offs can contribute to this. For instance, the conservative coverage of group 2 (Figure 15b and Figure 16b) could be influenced by the high-pass pre-filtering of the image intensities, since the ice-sheet interior is poor of intensity features, and/or by outlier culling based on the local variability of the flow direction, which was not applied by any other group. The limited coverage of group 4 compared to group 8, despite the choice of a lower cross-correlation SNR threshold (Table 9), is most likely due to the application of an additional culling procedure, which was however unclearly described. Interestingly, group 6, which used the same SNR threshold as group 8, also achieved a poorer coverage compared to the latter. This could be due to the coarser azimuth posting of the cross-correlations (150 m vs. 60 m), which reduces the chances of finding good cross-correlation peaks. The most comprehensive coverage was achieved by group 7 (Figure 15g and Figure 16g) and is most likely due to the combination of a relatively dense posting for the cross-correlation computations and to the use of variable cross-correlation window sizes (Table 9).
In fast moving areas the coverage differences are apparent from the insets in Figure 15 and Figure 16. Different strategies were applied by the participants to accommodate the full motion range. Group 2 used a large (2.5–3 km) search window size, which would be sufficient to track a 10 km/y motion. Groups 3, 7 and 8 used a priori information on the velocity magnitude to either increase the search window size (group 3), to selectively increase the cross-correlation window size and apply multi-looking (group 7) or to select the position of the image patches to cross-correlate (group 8).
Groups 1, 4 and 6, which did not implement any of such strategies, achieved a limited coverage close to the termini of the outlet glaciers (panels a, d and f in Figure 15 and Figure 16).
Finally, coverage differences can be noticed in the overall azimuth length of the ice velocity maps in Figure 15 and Figure 16. Besides group 4 (panel d), which unintentionally processed only part of the provided dataset, these differences are due to the processed azimuth bandwidth at the image borders, which is typically a user-configurable parameter during the image focusing step.
Concerning the properties of the velocity measurements, the LoS maps differ primarily due to image-wide trends in the east-west direction. To highlight such differences all results were resampled to a common grid with a posting of 4.5 arcsec in latitude and 13.5 arcsec in longitude, corresponding to about 125 m × 125 m on ground, and differenced with respect to the results of group 8. Figure 18 shows the results for an area of interest around Store Gletscher and Sermeq Avannarleq (STO4 and AVA1 GPS stations respectively). Smaller scale differences are also apparent, particularly in the Rink Isbræ area (Figure 19), and could be caused by geocoding differences, particularly in the north-south direction, which were not corrected by the rigid 2D shift applied in the pre-processing procedure.
More subtle differences in the LoS measurements can be inferred from Figure 20, in which the velocities derived by all groups are shown relative to a common reference point and are wrapped so that each color cycle represents a 200 m/y variation. For instance, the results from group 2 (Figure 20b) show some anomalous bands parallel to the azimuth direction, whereas the measurements of group 7 show a greater fringe detail (spatial resolution) particularly in the glacier shear zones.
The azimuth velocities also differ primarily due to East-West trends (Figure 21). Furthermore, the results of group 2 show a difference pattern with respect to other groups, which is correlated with the underlying azimuth velocity values (panel b in Figure 21 and Figure 22). The measurements of all groups are affected by ionospheric streaks [37], as seen in Figure 23, although for group 4 these artifacts are reduced due to a dedicated filtering step (Table 9). Finally, differences are observed concerning the spatial resolution of the measurements, which can be inferred from the clarity of the fringes on Store Gletscher and on Sermeq Avannarleq in Figure 23.
Concerning the quality parameters (Figure 17), the predicted error standard deviations provided by groups 1 and 6 are very similar, whereas those provided by group 8 are slightly lower, although the reason for this is not apparent based on the processing parameters. Since for this track the GCPs are located within 50 km of any measurement, the error predictions mainly reflect the local variability of the measurements. In contrast, the indices provided by groups 4, and 5 (Figure 17c,d), characterize the quality of the cross-correlation peaks, and are thus influenced by the signal-to-clutter ratio of amplitude/intensity features and by the level of phase coherence.

3.4.3. Validation

Measurement accuracy was assessed on exposed bedrock, identified using the ice-mask of [33], and on ice-covered surfaces, by comparison with GPS measurements. Summary statistics are reported in Table 10, whereas the differences for each GPS station are detailed in Table S5.
On bedrock outcrops, assumed to be stationary, median LoS velocities smaller than 1 m/y are measured by virtually all groups, whereas most groups show a negative bias in the azimuth component. Statistics on bedrock could not be computed reliably for group 5, since most of the variations found on bedrock fall within the quantization step of the measurements determined by the oversampling factor in Table 9, which corresponds to 3 m/y and 4.6 m/y respectively in the azimuth and LoS velocity components. For the remaining groups median absolute deviations lie between 1 and 4 m/y in LoS and 3 to 8 m/y in azimuth. A greater azimuth variability is expected due to the presence of azimuth streaks caused by ionospheric propagation [37]. Using the simplified approach of Raucoules et al. ([38]), the contribution of this error source is found to be approximately within +/− 50 m/y (Figure S2). For group 4 this error contribution is lower (Figure 23d), since an ionospheric filtering procedure was applied (Table 9).
On ice surfaces, in the area of the STO4 station, all groups underestimated the LoS-projected GPS measurement by an amount varying between 7.7 m/y (group 7) to 40 m/y (group 1) (Table S5). The profile in Figure 24c, taken from points β to β’ in Figure 21, shows that the differences between groups are spatially correlated and not confined only to the vicinity of the STO4 station. In particular, around the peak of the profile, at about 20 km from point β, groups 5, 7, and 8 are in good agreement with each other, whereas the results from other groups differ significantly. This is very likely due to a higher azimuth resolution of the measurements in the shear zone of the glacier, which can be inferred from Figure 23e,g,h. In turn this could be caused by the use of smaller azimuth cross-correlation window sizes (Table 9). We speculate that azimuth spatial resolution is also the main cause of the LoS differences seen at STO4. Concerning the azimuth measurements, ionospheric streaks explain most of the errors for groups 3, 5, 7, and 8 (Tables S5 and S6). For groups 1 and 6 an additional error contribution is provided by long-wavelength errors (Figure 21a,f), whereas group 2 (Figure 21b) shows spatially correlated differences with respect to all other groups, which seem correlated with the underlying azimuth velocities, and could thus depend on the use of a different scale-factor in the conversion from pixel offsets to velocity values.
In the area of the AVA1 station, concerning the LoS velocity component, group 7 shows quite a different profile compared to most other groups (Figure 24a) and also the highest disagreement with the GPS measurement (Table S5). Based on a visual inspection this is due to a correlation of group 7′s results with topography, most likely due to the pre-processing step used to convert Cartesian velocities to LoS and azimuth components. Groups 1, 4, and 6, show LoS gradients which differ from the other groups (Figure 18a,d,f) and seem incorrect, based on the profile through GPS station AVA1 in Figure 24a, in which point α’, located on a bedrock outcrop, has negative velocity values only for these groups. Such image-wide error trends are due to different strategies in the selection of the GCPs used for orbital refinement. Groups 1 and 6 manually selected points on bedrock outcrops, assuming these to have a zero motion, whereas group 4 selected points with a high cross-correlation SNR, although some of these do not fall in stationary areas (Figure 17c). Groups 3, 7 and 8 selected zero-motion GCPs automatically based on an ice-cover mask, whereas group 2 did not carry out orbital refinement. Concerning the azimuth measurements, ionospheric propagation errors (Table S6) account for the differences observed for group 5, and only partially for group 2 (Table S5), since the latter also shows an error pattern which is correlated with the velocity field (Figure 21b). Groups 1, 3, 4, 6, and to a lesser degree groups 7 and 8, underestimate the GPS azimuth velocity at AVA1, even more if the ionospheric errors are accounted for, as seen in Figure 24b and Table S5. As for the STO4 case, the smaller errors of groups 5, 7, and 8 could be due to the higher azimuth resolution of the measurements.
In the area of the RNK2 station, concerning the LoS velocity component, all groups except group 4 overestimate the GPS measurement (Table S5 and Figure 24c). The reason for this is not clear. A speculation is that this could be due to orbital refinement inaccuracies due to a less favourable distribution of rock surfaces and thus of GCPs, compared to the STO4 and AVA1 validation sites. From Figure 19, it is seen that the results of groups 1, 2, 7, and 8 differ only locally, whereas groups 3, 4 and to a lesser extent also 5 show differences which are spatially correlated with the underlying LoS velocity pattern (Figure 15). This could be due to relative geocoding shifts, which were not corrected by the constant shift model used in the pre-processing step. Groups 4, 5 and 6 also show differences with respect to image-wide East-West trends, due to the coregistration step and/or the GCP-based calibration of the offsets. Concerning the azimuth differences, group 2 (Figure 22b) shows a bias in the measurements towards lower values, which is also visible in the profile in Figure 24f, whereas group 4 shows spatially correlated differences with respect to other groups, which could be due to relative geocoding shifts and/or to the coregistration step.
Concerning the agreement between validation results and quality parameters, the observations are too few to draw reliable conclusions (Table S5). Two shortcomings of the approaches which were used can, however, be pointed out. Firstly, none of the quality parameters account for ionospheric errors, particularly those affecting the azimuth velocities, although this is the dominant error source on bedrock outcrops, where the GCPs used for calibration are located. Secondly, quality parameters, which reflect the strength of the cross-correlation peaks, e.g., cross-correlation SNR provided by group 4 (Figure 17), fail to characterize spatially-correlated errors, due for instance to the GCP-based calibration step (Figure 18d).

4. Discussion

The experiments discussed in Section 3 provide a comprehensive picture of the properties of the different SAR-based ice velocity measurements techniques. At the same time, it was unexpected for the results provided by the participants to differ often significantly in terms of accuracy and/or coverage. This is due to the variety of processing algorithms applied and to the influence of analyst choices on some critical processing steps.

4.1. Differential Synthetic Aperture Radar Interferometry

Concerning the DInSAR task (Section 3.1), most groups achieved LoS MADs better than 1 m/y (Table 4) and spatial resolutions of 45 m × 45 m on ground, using two- or four-pass approaches, applied to highly coherent ERS Tandem data. The most significant differences among groups consisted of image/wide trends due to the manual selection of GCPs for orbital refinement and absolute phase estimation, as well as to large-scale phase unwrapping errors (Figure 5).
Errors in the GCP selection process could be reduced by automated procedures, such as those implemented by some of the tasks 3 and 4 participants, in which an ice-cover mask was used to avoid points in non-stationary areas. Methods to ensure a homogeneous spatial distribution of the GCPs throughout the scene have also been proposed in literature [44]. Moreover, for several missions which became operational in the last decade, the orbit control and overall geolocation accuracy are significantly improved compared to the sensors considered in this study (e.g., [45,46,47]), which allows to simplify GCP-based calibration approaches by reducing the number of error-model parameters, and thus the estimation errors.
In contrast, no general automated method is available for detecting phase unwrapping errors in single interferograms, although synergies with range split-bandwidth or offset tracking techniques could be exploited for sensors with sufficiently high range resolutions [14]. Since phase unwrapping errors are due to local violations of the assumption that interferometric phase gradients are smaller than +/−π radians, their occurrence can be greatly reduced by considering only acquisitions with small temporal and spatial baselines, thereby limiting phase noise due to temporal decorrelation and sensitivity to topography, and by applying these methods towards the interior of the ice-sheet, to limit the phase contribution of velocity gradients.

4.2. Multi Aperture Interferometry

Concerning the MAI task (Section 3.2), the best MAD accuracy achieved with respect to the GPS surveys was of 7 m/y (Table 6), and azimuth velocities maps with a spatial resolution of about 90 m × 90 m were generated. The algorithmic differences among groups, which most significantly affected performance, were the handling of ionospheric propagation errors [10] and the reduction of spatially uncorrelated noise due to loss of interferometric coherence.
The mitigation of ionospheric streak errors was successfully carried out by one of the participants, although the procedure which was used cannot be considered of general applicability, since moving areas were excluded from the estimation, an empirically assigned threshold on the MAI phase. For systems with a sufficiently high range bandwidth, and datasets/areas meeting the above-mentioned conditions to reduce the occurrence of phase unwrapping errors, general ionospheric-streak correction approaches could be investigated based on the split-spectrum method [48,49,50], in which motion and ionospheric phase contributions are separated exploiting the frequency dependence of the latter.
Regarding the filtering of spatially uncorrelated noise, a difference among participants concerned the application of spatial filtering methods to each sub-band interferogram as opposed to their application to the final differential MAI interferogram. The superiority of the former approach was established for the case in which multi-looking is the only filtering method applied, based on analyses carried out on real datasets [35], as well as on theoretical derivations supported by numerical simulations [51]. However, the intercomparison between groups 1 and 4 in Section 3.2.2 suggests that these conclusions should be re-examined for the common case in which adaptive filtering methods, such as [27], are used. Unlike multi-looking, these techniques can in fact be expected to perform better on the differential MAI interferogram, in which fringes common to both sub-band interferograms (e.g., due to LoS motion) have been cancelled out.
The Multi Aperture Interferometry quality parameters provided by the participants were based on the local variability of the measurements. While this can provide a reasonable characterization for spatially uncorrelated error sources, such as decorrelation, other spatially-correlated error sources such as the above mentioned ionospheric streak errors and image-wide co-registration errors cannot be accounted for based on a single acquisition pair.

4.3. Offset Tracking

The results of both offset-tracking experiments (Section 3.3 and Section 3.4) were the most heterogeneous. The highest levels of agreement with the GPS measurements yielded LoS and azimuth MADs of 20 m/y and 40 m/y respectively for the ENVISAT ASAR dataset (task 3), and of 10 m/y in LoS and 20 m/y in azimuth for the ALOS PALSAR dataset (task 4). The spatial resolutions can be considered at best to approach the cross-correlation window spacing and at worst twice the size of the correlation window size [52], which on average amounted to values between 200 m × 170 m and 2.7 km × 1.7 km (ground range × azimuth) for task 3, and between 160 m × 120 m and 640 m × 1100 m for task 4. For both processing tasks the limited spatial-resolution of the measurements was one of the major error sources for areas with high velocity gradients, in agreement with the SAR vs. GPS comparisons reported in [23]. Furthermore, sub-optimal outlier culling and hole-filling procedures were also relevant in particular for the void-rich task 3 dataset, whereas GCP calibration errors as well as ionospheric streaks were noticeable for task 4, due to the larger spatial extent of the measurements and the lower (L-band) radar frequency. Finally, both offset-tracking tasks intentionally provided a challenge regarding the measurement of very high velocities (i.e., >2 km/y) close to the glacier termini.
Concerning the choice of the cross-correlation window sizes, which typically determine also the spacing at which the cross-correlations are computed, several groups chose large values, assigned based on analyst experience, and constant for the processing of the whole dataset. Although this approach typically succeeds in providing good-quality cross-correlations in ice-sheet regions characterized by different levels of coherence and intensity features, the optimal window sizes depend on surface and radar properties, which can vary significantly. Within the validation area of the task 3 dataset (Section 3.2.3), the smallest window sizes which provide good quality measurements lie between 1.0 km and 1.5 km in both ground-range and azimuth, whereas for the task 4 dataset (Section 3.4.2), good cross-correlation peaks at the GPS locations are found already for window sizes between 300 m and 500 m in each dimension. Since these values cannot be known beforehand, a good compromise to preserve spatial resolution would seem to be the hierarchal approach of Joughin ([13]), in which cross-correlations are attempted for each measurement with increasing window sizes until quality criteria are satisfied. Concerning the latter, the use of an NCC threshold, as proposed by Joughin ([13]), could be complemented by a second threshold on cross-correlation SNR [41], to account also for off-peak cross-correlation values in the quality characterization.
Measurements in fast-moving areas were provided by groups which applied processing approaches capable of accommodating the full range of motion. This was done by adjusting either the search window size or its position, or both, possibly guided by a priori velocity information. Increasing the search window size has typically a significant computational cost and can lead to a larger percentage of false-detections due to the presence of competing cross-correlation peaks. In contrast the positioning of search windows based on a priori information, such as a seasonal velocity map, depends on the velocity differences compared to the latter, as well as on the pixel spacing and temporal baseline of the SAR acquisitions. This study does not provide clear evidence in favor of one approach or the other; however, the results of task 4 (insets in Figure 15 and Figure 16) in particular show that in general one of these approaches should be implemented to derive measurements close to the glacier termini. Furthermore, since interferometric coherence is lost in fast moving areas, processing steps which enhance features with respect to radar speckle, such as the low- or high-pass filtering implemented by some of the participants, can be expected to improve the quality of the cross-correlation peaks in these areas. However, clearer indications on this can be provided by dedicated studies, such as [41].
Regarding outlier culling, several groups obtained good validation results on static bedrock with approaches based on the thresholding of NCC magnitude and cross-correlation SNR, using the same thresholds for the datasets in tasks 3 and 4. Further culling steps based on the analysis of local velocity magnitude and/or flow-direction variations proved to be effective, but also very sensitive to analyst choices. Concerning the hole-filling procedures, the results of task 3 suggest that conservative approaches are preferable, to avoid unphysical interpolations in areas with strong velocity gradients such as shear margins.
Spatially-correlated measurement biases were observed, related to the selection of GCPs and to ionospheric propagation errors in the case of the azimuth measurements. Concerning the latter, the same considerations discussed above for MAI hold also for offset-tracking. However, a general ionospheric-streak correction procedure based on the split-spectrum principle, as discussed in Section 4.2 for MAI, would only be applicable to azimuth offset-tracking measurements in coherent areas, where interferometry and phase unwrapping can also be carried out. Concerning GCPs, the considerations discussed in Section 4.1 for DInSAR hold also for offset-tracking.
Finally, concerning the quality parameters associated with the offset-tracking measurements, error estimates were either provided in the form of intermediate quality parameters, such as phase coherence, NCC magnitude or cross-correlation SNR, or in the form of 1σ uncertainties, accounting for the local variability of the measurements and in some cases for the position of the GCPs used for image-wide velocity calibration. The latter provided in several cases the correct order of magnitude for the observed errors, although it can be considered only a partial error characterization, since uncertainties due to the limited spatial resolution in both LoS and azimuth, to image-wide co-registration errors and to ionospheric propagation remain unaccounted for. In contrast intermediate parameters, although useful for outlier culling, provided a very poor indication of end-to-end measurement errors, besides being of difficult interpretation for an end-user unfamiliar with the processing techniques.

5. Conclusions

In the current context in which SAR-based ice-velocity products are being made available by an increasing number of independent providers, it would be desirable from an end-user point of view, for these products to share similar state-of-the-art properties in terms of measurement spatial resolution, accuracy, coverage, and error characterization. However, the experiments described in this paper show that single-track LoS and azimuth ice-velocity measurements, which represent the building block of higher-level products distributed to the end-users, can have quite heterogeneous properties when generated by different processing chains and analysts.
The contribution of this study is two-fold. Firstly, a total of 24 SAR-based ice-velocity processing results are described, generated by 13 international research groups, and a total of 24 GPS measurements used for validation are provided. The results cover all SAR processing techniques currently used to carry out ice-motion measurements, namely DInSAR, MAI (azimuth split-bandwidth) and offset-tracking. Both the SAR-based results and the GPS data provide a unique benchmark to assess the performance of processing chains implementing these measurement techniques.
Secondly, the algorithmic steps and/or parameters responsible for performance variability are analyzed in detail for each technique. Interestingly, the causes for sub-optimal performance are most often related to the application of procedures, which are strongly influenced by analyst choices, e.g., concerning processing parameters. Useful indications to automate some key processing steps, while at the same time yielding improved performances, are obtained by analyzing the available processing results, as discussed in Section 4 and detailed in the inter-comparison and validation sections of each experiment (Section 3).
Finally, although this study was carried out on archive data from ESA and ESA Third Party Missions, the considerations presented in Section 4 concerning the processing algorithms are not sensor-specific and are thus equally applicable to current (e.g., Sentinel-1a,b) and forthcoming (e.g., NASA-ISRO) SAR missions. In particular, since the latter by design support the systematic application of interferometric techniques, this will allow the synergies between different measurement algorithms to be exploited to a much higher degree, thus contributing to define more robust and analyst-independent procedures for several critical processing steps, besides of course noticeably improving the spatial resolution and the accuracy of the final ice-velocity products.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-4292/10/6/929/s1, Figure S1: Task 2 wrapped azimuth velocity maps, Figure S2: Task 4 estimated azimuth velocity ionospheric propagation errors, Table S1: GPS data overview, Table S2: Task 1 validation results, Table S3: Task2 validation results, Table S4: Task 3 validation results, Table S5: Task 4 validation results.

Author Contributions

Conceptualization, J.P.M.B.; Formal analysis, J.P.M.B., M.L.A., A.K., S.B.A., N.B., S.B., C.B., N.G., I.J., H.-S.J., A.L., J.M., J.N., E.R., K.S., T.N., B.S. and T.S.; Funding acquisition, J.D.; Investigation, J.P.M.B., M.L.A.; Software, M.K.; Validation, M.L.A.; Writing, J.P.M.B., J.D. and A.K.

Funding

This research as well as the APC were funded by European Space Agency contract 4000104815/11/I-NB. The SAR data was obtained through CAT-1P proposal 11339.

Acknowledgments

GPS campaign measurements over Upernavik Isstrøm were provided by Abbas Khan, DTU. All figures were generated with the Generic Mapping Tools software Version 5.2.1 [53].

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Scambos, T.A.; Dutkiewicz, M.J.; Wilson, J.C.; Bindschadler, R.A. Application of image cross-correlation to the measurement of glacier velocity using satellite image data. Remote Sens. Environ. 1992, 42, 177–186. [Google Scholar] [CrossRef]
  2. Goldstein, R.M.; Engelhardt, H.; Kamb, B.; Frolich, R.M. Satellite radar interferometry for monitoring ice sheet motion: Application to an Antarctic ice stream. Science 1993, 262, 1525–1530. [Google Scholar] [CrossRef] [PubMed]
  3. Nagler, T.; Rott, H.; Hetzenecker, M.; Wuite, J.; Potin, P. The Sentinel-1 mission: New opportunities for ice sheet observations. Remote Sens. 2015, 7, 9271–9389. [Google Scholar] [CrossRef]
  4. Mouginot, J.; Rignot, E.; Scheuchl, B.; Millan, R. Comprehensive annual ice sheet velocity mapping using Landsat-8, Sentinel-1, and RADARSAT-2 data. Remote Sens. 2017, 9, 364. [Google Scholar] [CrossRef]
  5. Joughin, I.; Smith, B.E.; Howat, I. A complete map of Greenland ice velocity from satellite data collected over 20 years. J. Glaciol. 2017, 64, 1–11. [Google Scholar] [CrossRef]
  6. Massonnet, D.; Feigl, K. Discrimination of geophysical phenomena in satellite radar interferograms. Geophys. Res. Lett. 1995, 22, 1537–1540. [Google Scholar] [CrossRef]
  7. Kwok, R.; Fahnestock, M.A. Ice sheet motion and topography from radar interferometry. IEEE Trans. Geosci. Remote Sens. 1996, 34, 189–200. [Google Scholar] [CrossRef]
  8. De Zan, F.; Monti Guarnieri, A. TOPSAR: Terrain observation by progressive scans. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2352–2360. [Google Scholar] [CrossRef]
  9. Scheiber, R.; Jäger, M.; Prats-Iraola, P.; De Zan, F.; Geudtner, D. Speckle tracking and interferometric processing of TerraSAR-X TOPS data for mapping nonstationary scenarios. IEEE JSTARS 2015, 8, 1709–1720. [Google Scholar] [CrossRef]
  10. Gray, A.L.; Mattar, K.E.; Vachon, P.W. InSAR results from the RADARSAT Antarctic mapping mission data: Estimation of glacier motion using a simple registration procedure. In Proceedings of the IGARSS, Seattle, WA, USA, 6–10 July 1998. [Google Scholar]
  11. Michel, R.; Rignot, E. Flow of glaciar Moreno, Argentina, from repeat-pass shuttle imaging radar images: Comparison of the phase correlation method with radar interferometry. J. Glaciol. 1999, 45, 93–100. [Google Scholar] [CrossRef]
  12. Derauw, D. DInSAR and coherence tracking applied to glaciology: The example of Shirase glacier. In Proceedings of the FRINGE’99, Liège, Belgium, 10–12 November 1999. [Google Scholar]
  13. Joughin, I. Ice-sheet velocity mapping: A combined interferometric and speckle-tracking approach. Ann. Glaciol. 2002, 34, 195–201. [Google Scholar] [CrossRef]
  14. Bamler, R.; Eineder, M. Accuracy of differential shift estimation by correlation and split-bandwidth interferometry for wideband and delta-k systems. IEEE Geosci. Remote Sens. Lett. 2005, 2, 151–155. [Google Scholar] [CrossRef]
  15. Bechor, N.B.D.; Zebker, H.A. Measuring two-dimensional movements using a single InSAR pair. Geopyhs. Res. Lett. 2006, 33, L16311. [Google Scholar] [CrossRef]
  16. Scheiber, R.; Moreira, A. Coregistration of interferometric SAR images using spectral diversity. IEEE Trans. Geosci. Remote Sens. 2000, 38, 2179–2191. [Google Scholar] [CrossRef]
  17. Madsen, S.N.; Zebker, H.A. Automated absolute phase retrieval in across-track interferometry. In Proceedings of the IGARSS, Houston, TX, USA, 26–29 May 1992. [Google Scholar]
  18. Gourmelen, N.; Kim, S.W.; Shepherd, A.; Park, J.W.; Sundal, A.V.; Bjornsson, H.; Palsson, F. Ice velocity determined using conventional and multiple-aperture InSAR. Earth Planet. Sci. Lett. 2011, 307, 156–160. [Google Scholar] [CrossRef] [Green Version]
  19. Joughin, I.; Kwok, R.; Fahnestock, M. Estimation of ice-sheet motion using satellite radar interferometry: Method and error analysis with application to Humboldt glacier, Greenland. J. Glaciol. 1996, 42, 564–575. [Google Scholar] [CrossRef]
  20. Gray, A.L.; Short, N.; Mattar, K.E.; Jezek, K.C. Velocities and flux of the Filchner ice shelf and its tributaries determined from speckle tracking interferometry. Can. J. Remote Sens. 2001, 27, 193–206. [Google Scholar] [CrossRef]
  21. Howat, I.M.; Negrete, A.; Smith, B.E. The Greenland Ice Mapping Project (GIMP) land classification and surface elevation datasets. Cryosphere 2014, 8, 1509–1518. [Google Scholar] [CrossRef]
  22. Reeh, N.; Mohr, J.J.; Madsen, S.N.; Oerter, H.; Gundestrup, N.S. Three-dimensional surface velocities of Storstrømmen glacier, Greenland, derived from radar interferometry and ice-sounding radar measurements. J. Glaciol. 2003, 49, 201–209. [Google Scholar] [CrossRef]
  23. Ahlstrøm, A.P.; Andersen, S.B.; Andersen, M.L.; Machguth, H.; Nick, F.M.; Joughin, I.; Reijmer, C.H.; van de Wal, R.S.W.; Boncori, J.P.M. Seasonal velocities of eight major marine-terminating outlet glaciers of the Greenland ice sheet from continuous in situ GPS instruments. Earth Syst. Sci. Data 2013, 5, 277–287. [Google Scholar] [CrossRef] [Green Version]
  24. Citterio, M.; van As, D.; Ahlstrøm, A.P.; Veicherts, M. Automatic weather stations for basic and applied glaciological research. Geol. Surv. Den. Greenl. Bull. 2015, 33, 69–72. [Google Scholar]
  25. Werner, C.; Wegmuller, U.; Strozzi, T.; Wiesmann, A. Gamma SAR and interferometric processing software. In Proceedings of the ERS-ENVISAT Symposium, Gotheburg, Sweden, 15–20 October 2000. [Google Scholar]
  26. Mohr, J.J.; Merryman Boncori, J.P.; Dall, J. Requirements and Interface Specification (RIS) for Preparation of Interferometric SAR Processor (COISP); Ver. 7.0.1., R 734; Technical University of Denmark: Kongens Lyngby, Denmark, 2008. [Google Scholar]
  27. Goldstein, R.M.; Werner, C.L. Radar interferogram filtering for geophysical applications. Geophys. Res. Lett. 1998, 25, 4035–4038. [Google Scholar] [CrossRef] [Green Version]
  28. Goldstein, R.M.; Zebker, H.A.; Werner, C.L. Satellite radar interferometry: Two-dimensional phase unwrapping. Radio Sci. 1988, 23, 713–720. [Google Scholar] [CrossRef]
  29. Costantini, M. A novel phase unwrapping method based on network programming. IEEE Trans. Geosci. Remote Sens. 1998, 36, 813–821. [Google Scholar] [CrossRef]
  30. Mohr, J.J.; Merryman Boncori, J.P. An error prediction framework for interferometric SAR data. IEEE Trans. Geosci. Remote Sens. 2008, 46, 1600–1613. [Google Scholar] [CrossRef] [Green Version]
  31. Chen, C.W.; Zebker, H.A. Network approaches to two-dimensional phase unwrapping: Intractability and two new algorithms. J. Opt. Soc. Am. 2000, 17, 401–414. [Google Scholar] [CrossRef]
  32. Mohr, J.J.; Reeh, N.; Madsen, S.N. Accuracy of three-dimensional glacier surface velocities derived from radar interferometry and ice-sounding radar measurements. J. Glaciol. 2003, 49, 210–222. [Google Scholar] [CrossRef]
  33. Citterio, M.; Ahlstrøm, A. Brief communication “The aerophotogrammetric map of Greenland ice masses”. Cryosphere 2013, 7, 445–449. [Google Scholar] [CrossRef]
  34. Merryman Boncori, J.P.; Pezzo, G. Measuring the north-south coseismic displacement component with high-resolution multi-aperture InSAR. Terra Nova 2015, 27, 28–35. [Google Scholar] [CrossRef]
  35. Jung, H.S.; Won, J.S.; Kim, S.W. An improvement of the performance of multi-aperture SAR interferometry (MAI). IEEE Trans. Geosci. Remote Sens. 2009, 47, 2859–2869. [Google Scholar] [CrossRef]
  36. Rosen, P.A.; Hensley, S.; Peltzer, G.; Simons, M. Updated repeat orbit interferometry package released. EOS Trans. AGU 2004, 85, 47. [Google Scholar] [CrossRef]
  37. Gray, A.L.; Mattar, K.E. Influence of ionospheric electron density fluctuations on satellite radar interferometry. Geophys. Res. Lett. 2000, 27, 1451–1454. [Google Scholar] [CrossRef] [Green Version]
  38. Raucoules, D.; De Michele, M. Assessing ionospheric influence on L-band SAR data: Implications on coseismic displacement measurements of the 2008 Sichuan earthquake. IEEE Geosci. Remote Sens. Lett. 2010, 7, 286–290. [Google Scholar] [CrossRef]
  39. Dall, J.; Kusk, A.; Nielsen, U.; Merryman Boncori, J.P. Ice velocity mapping using TOPS SAR data and offset tracking. In Proceedings of the Fringe 2015 Workshop, Frascati, Italy, 23–27 March 2015. [Google Scholar]
  40. Sarmap. Available online: http://www.sarmap.ch (accessed on 30 March 2018).
  41. De Lange, R.; Luckman, A.; Murray, T. Improvement of satellite radar feature tracking for ice velocity derivation by spatial frequency filtering. IEEE Trans. Geosci. Remote Sens. 2007, 45, 2309–2317. [Google Scholar] [CrossRef]
  42. Weickert, J. Anisotropic Diffusion in Image Processing, 1st ed.; Teubner: Stuttgart, Germany, 1998; ISBN 3-519-02606-6. [Google Scholar]
  43. Mouginot, J.; Scheuchl, B.; Rignot, E. Mapping of ice motion in Antarctica using synthetic-aperture radar data. Remote Sens. 2012, 4, 2753–2767. [Google Scholar] [CrossRef]
  44. Bähr, H. Orbital Effects in Spaceborne Synthetic Aperture Radar Interfermoetry. Ph.D. Thesis, Karlsruher Institut für Technologie, Karlsruhe, Germany, 2013; p. 66. [Google Scholar]
  45. Eineder, M.; Minet, C.; Steigenberger, P.; Cong, X.; Fritz, T. Imaging geodesy—Toward centimeter-level ranging accuracy with TerraSAR-X. IEEE Trans. Geosci. Remote Sens. 2011, 49, 661–671. [Google Scholar] [CrossRef]
  46. Peter, H.; Jäggi, A.; Fernández, J.; Escobar, D.; Ayuga, F.; Arnold, D.; Wermuth, M.; Hackel, S.; Otten, M.; Simons, W.; et al. Sentinel-1A—First precise orbit determination results. Adv. Space Res. 2017, 60, 879–892. [Google Scholar] [CrossRef]
  47. Schubert, A.; Miranda, N.; Geudtner, D.; Small, D. Sentinel-1a/b combined product geolocation accuracy. Remote Sens. 2017, 9, 607. [Google Scholar] [CrossRef]
  48. Rosen, P.; Hensley, S.; Chen, C. Measurement and mitigation of the ionosphere in L-band interferometric SAR data. In Proceedings of the IEEE Radar Conference, Arlington, TX, USA, 10–14 May 2010. [Google Scholar]
  49. Brcic, R.; Parizzi, A.; Eineder, M. Estimation and compensation of ionospheric delay for SAR interferometry. In Proceedings of the IGARSS, Honolulu, HI, USA, 25–30 July 2010. [Google Scholar]
  50. Gomba, G.; Parizzi, A.; De Zan, F.; Eineder, M.; Bamler, R. Toward operational compensation of ionospheric effects in SAR interferograms: The split-spectrum method. IEEE Trans. Geosci. Remote Sens. 2016, 54, 1446–1461. [Google Scholar] [CrossRef]
  51. De Zan, F.; Prats-Iraola, P.; Rodriguez-Cassola, M. On the dependence of Delta-k efficiency on multilooking. IEEE Geosci. Remote Sens. Lett. 2015, 8, 1745–1749. [Google Scholar] [CrossRef]
  52. Pritchard, H.; Murray, T.; Luckman, A.; Strozzi, T.; Barr, S. Glacier surge dynamics of Sortebræ, east Greenland, from synthetic aperture radar feature tracking. J. Geophys. Res. 2005, 110, F03005. [Google Scholar] [CrossRef]
  53. Wessel, P.; Smith, W.H.F. New, improved version of the generic mapping tools released. EOS Trans. AGU 1998, 79, 579. [Google Scholar] [CrossRef]
Figure 1. (a) SAR data frames (rectangles) and GPS stations/stakes used for validation (triangles) plotted on a horizontal velocity magnitude mosaic generated from Sentinel-1a/b SAR data acquired between October 2015 and October 2016 (©ENVEO). Insets represent the Storstrømmen (b); Upernavik (c) and central west ice margin (d) test sites referred to in Table 1.
Figure 1. (a) SAR data frames (rectangles) and GPS stations/stakes used for validation (triangles) plotted on a horizontal velocity magnitude mosaic generated from Sentinel-1a/b SAR data acquired between October 2015 and October 2016 (©ENVEO). Insets represent the Storstrømmen (b); Upernavik (c) and central west ice margin (d) test sites referred to in Table 1.
Remotesensing 10 00929 g001
Figure 2. Interferometric coherence of the I1 (a); I2 (b); I3 (c) and I4 (d) SAR datasets (Table 1).
Figure 2. Interferometric coherence of the I1 (a); I2 (b); I3 (c) and I4 (d) SAR datasets (Table 1).
Remotesensing 10 00929 g002
Figure 3. Task 1 LoS velocity maps provided by groups 1 to 6 (af). Positive values correspond to motion towards the radar along the LoS vector (from ground to satellite), the ground projection of which at scene center is shown by the arrow in (a). Triangles represent the GPS stations used for validation.
Figure 3. Task 1 LoS velocity maps provided by groups 1 to 6 (af). Positive values correspond to motion towards the radar along the LoS vector (from ground to satellite), the ground projection of which at scene center is shown by the arrow in (a). Triangles represent the GPS stations used for validation.
Remotesensing 10 00929 g003
Figure 4. Task 1 quality parameters provided by groups 1 to 6 (af). For group 2 (b) the quality parameter is interferometric coherence, whereas for all other groups it is the predicted LoS velocity error standard deviation in m/y. Triangles represent the GPS stations used for validation. Crosses represent GCPs used for orbital refinement and absolute phase estimation.
Figure 4. Task 1 quality parameters provided by groups 1 to 6 (af). For group 2 (b) the quality parameter is interferometric coherence, whereas for all other groups it is the predicted LoS velocity error standard deviation in m/y. Triangles represent the GPS stations used for validation. Crosses represent GCPs used for orbital refinement and absolute phase estimation.
Remotesensing 10 00929 g004
Figure 5. LoS velocity differences between groups 1 to 5 (ae) with respect to group 6.
Figure 5. LoS velocity differences between groups 1 to 5 (ae) with respect to group 6.
Remotesensing 10 00929 g005
Figure 6. Task 2 azimuth velocity maps provided by groups 1 to 4 (ad). Positive values correspond to motion along the satellite flight path, the ground projection of which is shown with an arrow in (a). The triangles represent the GPS stations used for validation.
Figure 6. Task 2 azimuth velocity maps provided by groups 1 to 4 (ad). Positive values correspond to motion along the satellite flight path, the ground projection of which is shown with an arrow in (a). The triangles represent the GPS stations used for validation.
Remotesensing 10 00929 g006
Figure 7. Task 2 azimuth quality parameters provided by groups 2 (a); 3 (b); and 4 (c). For all groups the quality parameter is the predicted azimuth velocity error standard deviation in m/y. Triangles represent the GPS stations used for validation.
Figure 7. Task 2 azimuth quality parameters provided by groups 2 (a); 3 (b); and 4 (c). For all groups the quality parameter is the predicted azimuth velocity error standard deviation in m/y. Triangles represent the GPS stations used for validation.
Remotesensing 10 00929 g007
Figure 8. Azimuth velocity differences between groups 1 to 3 with respect to group 4.
Figure 8. Azimuth velocity differences between groups 1 to 3 with respect to group 4.
Remotesensing 10 00929 g008
Figure 9. Task 3 LoS velocity maps provided by groups 1 to 6 (af). Positive values correspond to motion towards the radar along the LoS vector (from ground to satellite), the ground projection of which at scene center is shown by the arrow in (a). The insets correspond to the area in the dashed rectangle and show the GPS stations used for validation.
Figure 9. Task 3 LoS velocity maps provided by groups 1 to 6 (af). Positive values correspond to motion towards the radar along the LoS vector (from ground to satellite), the ground projection of which at scene center is shown by the arrow in (a). The insets correspond to the area in the dashed rectangle and show the GPS stations used for validation.
Remotesensing 10 00929 g009
Figure 10. Task 3 azimuth velocity maps provided by groups 1 to 6 (af). Positive values correspond to motion along the satellite flight path, the ground projection of which is shown with an arrow in (a).
Figure 10. Task 3 azimuth velocity maps provided by groups 1 to 6 (af). Positive values correspond to motion along the satellite flight path, the ground projection of which is shown with an arrow in (a).
Remotesensing 10 00929 g010
Figure 11. Task 3 quality parameters. Groups 1 slant-range (a) and azimuth error standard deviation (b); group 3 NCC amplitude (c); group 4 cross-correlation SNR (d); group 5 NCC amplitude (e); and group 6 slant-range (f) and azimuth error standard deviation (g). Triangles represent the GPS stations used for validation. Crosses represent GCPs used for calibration.
Figure 11. Task 3 quality parameters. Groups 1 slant-range (a) and azimuth error standard deviation (b); group 3 NCC amplitude (c); group 4 cross-correlation SNR (d); group 5 NCC amplitude (e); and group 6 slant-range (f) and azimuth error standard deviation (g). Triangles represent the GPS stations used for validation. Crosses represent GCPs used for calibration.
Remotesensing 10 00929 g011
Figure 12. Task 3 LoS velocity differences of groups 1 to 5 (ae) with respect to group 6. The dashed line in the inset represent the trace of the profiles shown in Figure 14.
Figure 12. Task 3 LoS velocity differences of groups 1 to 5 (ae) with respect to group 6. The dashed line in the inset represent the trace of the profiles shown in Figure 14.
Remotesensing 10 00929 g012
Figure 13. Task 3 azimuth velocity differences of groups 1 to 5 (ae) with respect to group 6.
Figure 13. Task 3 azimuth velocity differences of groups 1 to 5 (ae) with respect to group 6.
Remotesensing 10 00929 g013
Figure 14. Task 3 velocity profiles along the dashed line in Figure 12 and Figure 13. Crosses represent the projected GPS LoS and azimuth velocities for the stations intersected by the profile. Colors represent the results of groups 1 to 6.
Figure 14. Task 3 velocity profiles along the dashed line in Figure 12 and Figure 13. Crosses represent the projected GPS LoS and azimuth velocities for the stations intersected by the profile. Colors represent the results of groups 1 to 6.
Remotesensing 10 00929 g014
Figure 15. Task 4 LoS velocity maps provided by groups 1 to 8 (ah). Positive values correspond to motion towards the radar along the LoS vector (from ground to satellite), the ground projection of which at scene center is shown by the arrow in (a). The insets show the GPS stations used for validation.
Figure 15. Task 4 LoS velocity maps provided by groups 1 to 8 (ah). Positive values correspond to motion towards the radar along the LoS vector (from ground to satellite), the ground projection of which at scene center is shown by the arrow in (a). The insets show the GPS stations used for validation.
Remotesensing 10 00929 g015
Figure 16. Task 4 azimuth velocity maps provided by groups 1 to 8 (ah). Positive values correspond to motion along the satellite flight path, the ground projection of which is shown by the arrow in (a).
Figure 16. Task 4 azimuth velocity maps provided by groups 1 to 8 (ah). Positive values correspond to motion along the satellite flight path, the ground projection of which is shown by the arrow in (a).
Remotesensing 10 00929 g016
Figure 17. Task 4 quality parameters. Groups 1 slant-range (a) and azimuth error standard deviation (b); group 4 cross-correlation signal to noise ratio (c); group 5 normalized cross-correlation amplitude (d); group 6 slant-range (e) and azimuth error standard deviation (f); group 8 slant-range (g) and azimuth error standard deviation (h). Triangles represent the GPS stations used for validation. Crosses represent GCPs used for velocity calibration.
Figure 17. Task 4 quality parameters. Groups 1 slant-range (a) and azimuth error standard deviation (b); group 4 cross-correlation signal to noise ratio (c); group 5 normalized cross-correlation amplitude (d); group 6 slant-range (e) and azimuth error standard deviation (f); group 8 slant-range (g) and azimuth error standard deviation (h). Triangles represent the GPS stations used for validation. Crosses represent GCPs used for velocity calibration.
Remotesensing 10 00929 g017
Figure 18. Task 4 LoS velocity differences between groups 1 to 7 (ag) with respect to group 8 in the area of GPS stations STO4 and AVA1. The dashed lines represent the traces of the profiles in Figure 24a,c.
Figure 18. Task 4 LoS velocity differences between groups 1 to 7 (ag) with respect to group 8 in the area of GPS stations STO4 and AVA1. The dashed lines represent the traces of the profiles in Figure 24a,c.
Remotesensing 10 00929 g018
Figure 19. Task 4 LoS velocity differences between groups 1 to 7 (ag) with respect to group 8 in the area of RNK2 GPS station. The dashed lines represent the trace of the profiles in Figure 24e.
Figure 19. Task 4 LoS velocity differences between groups 1 to 7 (ag) with respect to group 8 in the area of RNK2 GPS station. The dashed lines represent the trace of the profiles in Figure 24e.
Remotesensing 10 00929 g019
Figure 20. Task 4 wrapped LoS velocity maps provided by groups 1 to 8 (ah) in the area of GPS stations STO4 and AVA1. Each color cycle represents a 200 m/y variation. Velocities are referred to the star (zero velocity).
Figure 20. Task 4 wrapped LoS velocity maps provided by groups 1 to 8 (ah) in the area of GPS stations STO4 and AVA1. Each color cycle represents a 200 m/y variation. Velocities are referred to the star (zero velocity).
Remotesensing 10 00929 g020
Figure 21. Task 4 azimuth velocity differences between groups 1 to 7 (ag) and group 8 in the area of GPS stations STO4 and AVA1. The dashed lines represent the traces of the profiles in Figure 24b,d.
Figure 21. Task 4 azimuth velocity differences between groups 1 to 7 (ag) and group 8 in the area of GPS stations STO4 and AVA1. The dashed lines represent the traces of the profiles in Figure 24b,d.
Remotesensing 10 00929 g021
Figure 22. Task 4 azimuth velocity differences between groups 1 to 7 (ag) and group 8 in the area of RNK2 GPS station. The dashed line represents the trace of the profiles in Figure 24f.
Figure 22. Task 4 azimuth velocity differences between groups 1 to 7 (ag) and group 8 in the area of RNK2 GPS station. The dashed line represents the trace of the profiles in Figure 24f.
Remotesensing 10 00929 g022
Figure 23. Task 4 wrapped azimuth velocity maps provided by groups 1 to 8 (ah) in the area of GPS stations STO4 and AVA1. Each color cycle represents a 200 m/y variation. Velocities are referred to the star (zero velocity).
Figure 23. Task 4 wrapped azimuth velocity maps provided by groups 1 to 8 (ah) in the area of GPS stations STO4 and AVA1. Each color cycle represents a 200 m/y variation. Velocities are referred to the star (zero velocity).
Remotesensing 10 00929 g023
Figure 24. Task 4 LoS (a,c,e) and azimuth (b,d,f) velocity profiles for groups 1–8 along the dashed lines shown in Figure 18, Figure 19, Figure 21 and Figure 22. Crosses represent the projected GPS velocity components.
Figure 24. Task 4 LoS (a,c,e) and azimuth (b,d,f) velocity profiles for groups 1–8 along the dashed lines shown in Figure 18, Figure 19, Figure 21 and Figure 22. Crosses represent the projected GPS velocity components.
Remotesensing 10 00929 g024
Table 1. SAR data overview. Locations are shown in Figure 1.
Table 1. SAR data overview. Locations are shown in Figure 1.
IdLocationSensorAcquisitionBt aBp bDf cΘ dΦ e∆r f∆a g
Date(Days)(m)(Hz)(°)(°)(m)(m)
I1StorstrømmenERS-11996-01-311140067−1517.903.94
ERS-21996-02-01
I2StorstrømmenERS-11996-04-101−2024567−1517.903.94
ERS-21996-04-11
I3UpernavikASAR2010-07-1135306−4667−1587.804.01
(ENVISAT)2010-08-19
I4Central westPALSAR2009-11-20465836152−1664.683.14
ice margin(ALOS-1)2010-01-05
a Temporal baseline; b Perpendicular baseline at scene center; c Doppler centroid difference at scene center; d,e Elevation and orientation angles at scene center in Equation (1); f,g Slant-range and azimuth pixel spacing. Acronyms: European Remote Sensing satellites 1 and 2 (ERS-1,2), Advanced Synthetic Aperture Radar (ASAR), Environmental Satellite (ENVISAT), Phased Array L-band Synthetic Aperture Radar (PALSAR), Advanced Land Observing Satellite 1 (ALOS-1).
Table 2. Participant tasks and deliverables.
Table 2. Participant tasks and deliverables.
TaskSAR DatasetsProcessing LevelAuxiliary DataProcessing AlgorithmDeliverables
1I1 and I2SLCDEM a, Precise SVs bDInSARLoS velocity map, quality parameters, GCPs
2I2RAWDEM a, Precise SVs bMAIAzimuth velocity map, quality parameters, GCPs
3I3RAWDEM a, Precise SVs cOffset TrackingLoS and azimuth velocity maps, quality parameters, GCPs
4I4RAWDEM aOffset TrackingLoS and azimuth velocity maps, quality parameters, GCPs
a 90 m posting DEM of Greenland (Byrd Polar Research Center, Ohio State University); b Precise State Vectors from Delft University of Technology; c DORIS Precise State Vectors. Acronyms: Single Look Complex (SLC), Differential Synthetic Aperture Radar Interferometry (DInSAR), Multi Aperture Interferometry (MAI), Ground Control Points (GCPs).
Table 3. Task 1 main processing parameters.
Table 3. Task 1 main processing parameters.
GroupG1G2G3G4G5G6
MethodDEMEDEMEDEMEDDDEMEDD
Quality parameter aσγσσσσ
Coregistration bNw: 20 × 20,
Ws: 64 × 256,
Npoly: 2
Nw: 64 × 128,
Ws: 64 × 64,
Npoly: 2.
Nw: 20 × 20,
Ws: 64 × 256,
Npoly: 2
Nw: 8 × 16,
Ws: 128 × 128,
Npoly: 0
Nw: 8 × 16,
Ws: 128 × 128,
Npoly: 0
Nw: 20 × 20,
Ws: 64 × 256,
Npoly: 2
Interferogram filtering cNL: 2 × 10, Goldstein
(Nfft: 32, α: 0.5)
NL: 2 × 10, Goldstein
(Nfft: 32, α: 0.8)
NL: 2 × 10, Goldstein
(Nfft: 32, α: 0.5)
NL: 3 × 15NL: 3 × 15NL: 2 × 10, Goldstein
(Nfft: 32, α: 0.5)
Phase unwrapping dBranch cut
thr = 0.41)
MCF
thr = 0.3)
MCF
thr = 0.41)
MCF
thr = 0.3)
MCF
thr = 0.3)
MCF
thr = 0.41)
Geophysical inversion49 GCPs1 GCP20 GCPs20 GCPs 20 GCPs20 GCPs
Geocoding eEQA,
280 m × 250 m
PS,
45 m × 45 m
EQA,
45 m × 45 m
EQA,
45 m × 45 m
EQA,
45 m × 45 m
EQA,
45 m × 45 m
a σ: per-pixel error standard deviation (m/y), γ: interferometric phase coherence; b Nw: number of cross-correlation windows (range × azimuth), Ws: cross-correlation window size (range × azimuth), Npoly: polynomial degree for co-registration refinement; c NL: number of interferogram looks (range × azimuth), Nfft: Goldstein filter 2D FFT size and α parameter [27]; d γthr: coherence threshold; e EQA: lat/lon grid, PS: NSIDC Polar Stereographic North (EPSG 3413). Acronyms: Digital Elevation Model Elimination (DEME), Double Difference (DD). Minimum Cost Flow (MCF). Ground Control Points (GCPs), Equal Area (EQA), Polar Stereographic (PS).
Table 4. Task 1 validation statistics.
Table 4. Task 1 validation statistics.
ParameterComponentGroup
G1G2G3G4G5G6
SAR–GPS
N aLoS101010101010
Median bLoS40.48−0.711.130.551.070.93
MAD cLoS13.700.550.580.590.580.81
Bedrock
MedianLoS3.1−1.410.36−0.210.240.30
MADLoS29.70.990.820.740.890.71
a Number of co-located SAR and GPS measurements; b Median of SAR differences with respect to GPS in m/y; c Median Absolute Deviation (MAD) of SAR differences with respect to GPS in m/y.
Table 5. Task 2 main processing parameters.
Table 5. Task 2 main processing parameters.
GroupG1G2G3G4
Quality parameter a-σσσ
Sub-aperture formation bn: 0.4n: 0.4n: 0.4n: 0.5
Interferogram filtering cNL: 4 × 20,
Goldstein (Nfft: 64, α: 0.8)
NL: 4 × 20,
Mean filter (12 pixel radius)
NL: 4 × 20,
Mean filter (18 pixel radius)
NL: 4 × 20,
Goldstein (Nfft: 64, α: 0.5),
Ionospheric filter
Phase unwrapping d---MCF
Geophysical inversionΔf, Vg scalingΔf, Vg scalingΔf, Vg scalingOrbital refinement,
La scaling
Geocoding ePS,
90 m × 90 m
EQA,
90 m × 90 m
EQA,
90 m × 90 m
PS,
90 m × 90 m
a σ: per-pixel error standard deviation (m/y); b n: normalized Doppler frequency separation between sub-apertures (as a fraction of the 1680 Hz PRF); c NL: number of interferogram looks (range × azimuth), Nfft: Goldstein filter 2D FFT size, and α parameter [27], r: mean filter radius in pixels; d La: nominal antenna length, Δf: sub-aperture frequency separation, Vg: ground (beam footprint) velocity; e EQA: lat/lon grid, PS: NSIDC Polar Stereographic North (EPSG 3413).
Table 6. Task 2 validation statistics.
Table 6. Task 2 validation statistics.
ParameterComponentGroup
G1G2G3G4
SAR–GPS
N aAzimuth9101010
Median bAzimuth4.9711.64.783.58
MAD cAzimuth10.3517.9913.596.80
Bedrock
MedianAzimuth6.83−3.82−3.240.12
MADAzimuth12.8919.4617.083.39
a Number of co-located SAR and GPS measurements; b Median of SAR differences with respect to GPS in m/y; c Median Absolute Deviation (MAD) of SAR differences with respect to GPS in m/y.
Table 7. Task 3 main processing parameters.
Table 7. Task 3 main processing parameters.
GroupG1G2G3G4G5G6
Method aICCICCICCICCICCICC
Quality parameter bσ-NCCSNRNCCσ
Pre-processingCoregistration1 × 5 multi-look; High-passCoregistrationCoregistration--
Offset computation cOsf: 2,
Ws: 1300 × 1000,
Wp: 200 × 200
Osf: 2,
Ws: 1600 × 1600,
Wp: 100 × 100,
Ss: 2300 × 2300
Osf: 16,
Ws: 650 × 256,
Wp: 100 × 100
Osf: 2,
Ws: 1300 × 1000,
Wp: 240 × 200
Osf: 8
Ws: 2000 × 400,
Wp: 400 × 200,
Ss: 650 × 130
Osf: 2,
Ws: 1300 × 1000,
Wp: 200 × 200
Outlier ulling dSNRthr = 7.0,
Max. velocity
SNRthr = 2.5,
Max. velocity,
Local variability (magnitude and direction)
NCCthr = 0.15,
Local variance
SNRthr = 4.0,
Unclear
NCCthr = 0.05SNRthr = 7.0,
NCCthr = 0.05
FilteringMoving average (5 × 5 window)-Anisotropic diffusion --Moving average (5 × 5 window)
Geophysical inversionStationary GCPs--High-SNR GCPs-Stationary GCPs
Geocoding eEQA,
280 m × 250 m
PS,
100 m × 100 m
PS,
100 m × 100 m
PS,
90 m × 90 m
PS,
90 m × 90 m
PS,
250 m × 250 m
a ICC: Intensity cross-correlation; b σ: per-pixel error standard deviation (m/y), NCC: normalized cross-correlation, SNR: cross-correlation peak Signal to Noise Ratio; c Osf: oversampling factor, Ws, Wp: cross-correlation window size and spacing in m (ground range × azimuth), Ss: search window size in m (ground range × azimuth); d SNRthr: SNR threshold, NCCthr: NCC threshold; e EQA: lat/lon grid, PS: NSIDC Polar Stereographic North (EPSG 3413).
Table 8. Task 3 validation statistics.
Table 8. Task 3 validation statistics.
ParameterComponentGroup
G1G2G3G4G5G6
SAR–GPS
N aLoS743486
Azimuth7434106
Median bLoS−12.3−16−5.45−15.8−15.61.01
Azimuth−34.9−11.4−67.04.7212.3−14.1
MAD cLoS27.623.521.150.5145021.2
Azimuth52.356.170.090.977242.6
Bedrock
Median bLoS−1.100.40−2.201.1010.20.60
Azimuth0.20−39.3−8.10−1.505.200.0
MAD cLoS2.436.745.134.923.262.16
Azimuth1.1529.955.679.2410.51.82
a Number of co-located SAR and GPS measurements; b Median of SAR-GPS differences in m/y; c Median absolute deviation of SAR differences with respect to GPS in m/y.
Table 9. Task 4 main processing parameters.
Table 9. Task 4 main processing parameters.
GroupG1G2G3G4G5G6G7G8
Method aICCICCACCICCICCICCCCC and ACCICC
Quality parameter bσ--SNRNCCσσσ
Pre-processingCoregistration1 × 2 multi-look, High-pass filter-Coregistration-CoregistrationCoarse ACC
Interf. formation
-
Offset computation cOsf: 2,
Ws: 960 × 800,
Wp: 150 × 150
Osf: 2,
Ws: 1120 × 940,
Wp: 150 × 125,
Ss: 3000 × 2500
Osf: none,
Ws: 1000 × 1000,
Wp: 240 × 200,
Flow-dependent search window size
Osf: 2,
Ws: 480 × 600,
Wp: 90 × 112
Osf: 8,
Ws: 720 × 300,
Wp: 150 × 60,
Ss: 240 × 100
Osf: 2,
Ws: 960 × 800,
Wp: 150 × 150
Osf: 2,
Ws: 180 × 80,
480 × 200,
1440 × 600,
3 × 3 multi-look Wp: 180 × 80
Osf: 2,
Ws: 960 × 400,
Wp: 150 × 60
Flow dependent window position
Outlier culling dSNRthr = 7.0,
Max. velocity
SNRthr = 1.8,
Max. velocity,
Local variability (mag. and dir.)
Median filter,
Coastline mask
SNRthr = 4.0,
Unclear
NCCthr = 0.04SNRthr = 7.0,
Max. velocity
NCCthr = 0.07, 0.18,
Median filter,
Coastline mask
SNRthr = 7.0,
NCCthr = 0.05,
FilteringMoving average --Ionospheric --Moving averageMoving average
Geophysical inversionStationary GCPs-Stationary GCPsHigh-SNR GCPs-Stationary GCPsStationary GCPsStationary GCPs
Geocoding eEQA,
280 m × 250 m
PS,
100 m × 100 m
PS,
150 m × 150 m
PS,
90 m × 90 m
PS,
90 m × 90 m
PS,
250 m × 250 m
PS,
200 m × 200 m
PS,
250 m × 250 m
a ACC: Amplitude cross-correlation, CCC: Complex cross-correlation, ICC: Intensity cross-correlation; b σ: per-pixel error standard deviation (m/y), NCC: normalized cross-correlation, SNR: cross-correlation peak Signal to Noise Ratio; c Osf: oversampling factor, Ws, Wp: cross-correlation window size and spacing respectively in m (ground range × azimuth), Ss: search window size in m (ground range × azimuth); d SNRthr: SNR threshold, NCCthr: NCC threshold. e EQA: lat/lon grid, PS: NSIDC Polar Stereographic North (EPSG 3413).
Table 10. Task 4 validation statistics.
Table 10. Task 4 validation statistics.
StationComponentGroup
G1G2G3G4G5G6G7G8
SAR–GPS
N aLoS23313132
Azimuth23313132
Median bLoS−3.22−1.90−16.1−9.84−14.127.834.42.4
Azimuth−21.2−5.4−23.7−10.613.620.12.63−8.39
MAD cLoS36.317.011.6-23.4-22.022.2
Azimuth43.923.319.1-22.7-17.641.2
Bedrock
Median bLoS−0.71−0.28−0.26−0.91-−1.050.6−0.22
Azimuth−4.86−8.13−1.0−0.63-−2.460.56−1.13
MAD cLoS1.833.571.214.01-2.361.570.92
Azimuth8.437.386.553.13-6.665.596.43
a Number of co-located SAR and GPS measurements; b Median of SAR-GPS differences in m/y; c Median absolute deviation of SAR-GPS differences in m/y.

Share and Cite

MDPI and ACS Style

Merryman Boncori, J.P.; Langer Andersen, M.; Dall, J.; Kusk, A.; Kamstra, M.; Bech Andersen, S.; Bechor, N.; Bevan, S.; Bignami, C.; Gourmelen, N.; et al. Intercomparison and Validation of SAR-Based Ice Velocity Measurement Techniques within the Greenland Ice Sheet CCI Project. Remote Sens. 2018, 10, 929. https://doi.org/10.3390/rs10060929

AMA Style

Merryman Boncori JP, Langer Andersen M, Dall J, Kusk A, Kamstra M, Bech Andersen S, Bechor N, Bevan S, Bignami C, Gourmelen N, et al. Intercomparison and Validation of SAR-Based Ice Velocity Measurement Techniques within the Greenland Ice Sheet CCI Project. Remote Sensing. 2018; 10(6):929. https://doi.org/10.3390/rs10060929

Chicago/Turabian Style

Merryman Boncori, John Peter, Morten Langer Andersen, Jørgen Dall, Anders Kusk, Martijn Kamstra, Signe Bech Andersen, Noa Bechor, Suzanne Bevan, Christian Bignami, Noel Gourmelen, and et al. 2018. "Intercomparison and Validation of SAR-Based Ice Velocity Measurement Techniques within the Greenland Ice Sheet CCI Project" Remote Sensing 10, no. 6: 929. https://doi.org/10.3390/rs10060929

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop