Contributors Lin Gilbert (University of Leeds), Sebastian B. Simonsen (Technical University of Denmark), Jan Wuite (ENVEO IT GmbH)
Issued by: University of Leeds / Lin Gilbert
Issued Date: 27/04/2021
Ref: C3S_312b_Lot4.D1.IS.2_v3.0_IV_Algorithm_Theoretical_Basis_Document_i1.0
Official reference number service contract: 2018/C3S_312b_Lot4_EODC/SC2
History of modifications
List of datasets covered by this document
Related documents
Acronyms
Scope of the document
This document is the Algorithm Theoretical Basis Document for Ice Velocity (IV) as part of the Copernicus Ice Sheets and Ice Shelves service. It describes the satellite-mounted instruments, auxiliary datasets, auxiliary models and basic algorithms used to create the data products.
Executive summary
The service addresses three essential climate variables (ECVs) by providing four separate products.
- Ice velocity is given for Greenland in product D3.IS.4
- Gravimetric mass balance is given for Greenland and Antarctica in product D3.IS.5
- Surface elevation change is given for
- Antarctica in product D3.IS.6.1
- Greenland in product D3.IS.6.2
We describe here the inputs, processing stages and outputs for Ice Velocity.
1. Greenland Ice Sheet Velocity – D3.IS.4
1.1. Instruments
The primary satellite source for generation of the ice velocity (IV) fields is Copernicus Sentinel-1A and Sentinel-1B. The Sentinel-1 (S1) satellites carry a C-band synthetic aperture radar (SAR) instrument providing high resolution SAR images in different acquisition modes. The Interferometric Wide-Swath mode (IW) is the main acquisition mode over land areas, including Greenland. IW allows combining a large swath width (250 km) with a moderate geometric resolution (5 m x 20 m).
Sentinel-1A provides continuous coverage for the entire Greenland Ice Sheet margin since October 2014 with a 12-day repeat pass period. Sentinel-1B was launched in April 2016 and reduced the repeat pass period to 6 days. Besides the continuous coverage of the margin, the ice sheet is covered completely during dedicated annual (winter) campaigns, with 4 to 6 repeat observations per track. In 2019 additional tracks were added for continuous acquisition covering most of the interior of the ice sheet in ascending and descending direction, providing new opportunities for interferometric applications. The S1 acquisition plan for Greenland makes the constellation the primary source for year-round monitoring of IV. The main instrument parameters and IW mode are described in Table 1 and Table 2 .
Table 1: Sentinel-1 sensor specifications
Property | Value |
Sensor | C-SAR (C-band Synthetic Aperture Radar) |
Emitted frequency (GHz) | 5.4 |
Pulse repetition frequency (Hz) | 1 000 - 3 000 |
Pulse duration (microseconds) | 5-100 |
Bandwidth (MHz) | 0-100 |
Antenna length (m) | 12.3 |
Antenna beamwidth (degrees) | 0.23 |
Slant range resolution (m) | 5 |
Azimuth resolution (m) | 20 |
Table 2: Characteristics of Sentinel-1 IW mode
Parameter | Value |
Swath width (ground range) | 250 km |
Nom. L1b product length | 50 km |
Full performance incidence angle range | 20°- 46° |
Data access incidence angle range | 30° - 42° |
Azimuth resolution | 20 m |
Ground range resolution | 5 m |
Polarizations | HH+HV, VV+VH, VV, HH |
1.2. Input and auxiliary data
1.2.1. SAR data
Primary input data for IV retrieval are repeat pass Sentinel-1A and Sentinel-1B Level-1 Single Look Complex (SLC) images acquired in the Interferometric Wide (IW) swath mode. These products consist of focused SAR data geo-referenced using orbit and altitude data from the satellite and provided in zero-Doppler slant-range geometry. IW implements the Terrain Observation by Progressive Scans (TOPS) mode, consisting of 3 or 5 sub-swaths, each consisting of a series of bursts. Each burst has been processed as a separate SLC image. ESA provides Sentinel-1 data free of charge through the Copernicus Open Access Hub (SciHub) and various mirror sites.
1.2.2. Auxiliary data
Auxiliary data, used in the processing chain for IV production, include precise orbit files, a digital elevation model and ocean/land masks. Details are provided in the following:
1.2.2.1. Orbits
Orbit files (Precise Orbit Ephemerides - POE) are required for IV processing and enable accurate co-registration. Precise orbit files are available as separate products and are provided by the Precise Orbit Determination (POD) service for SENTINEL-1 approximately 3 weeks after acquisition. POE files cover approximately 28 hours and contain orbit state vectors at 10-second intervals.
1.2.2.2. DEM
A contemporary high-resolution digital elevation model (DEM) is required in the processing chain for various purposes, including co-registration, calculation of the vertical component of velocity (vz) and final geocoding of the product. The DEM used for the annual Greenland Ice Sheet velocity map is compiled from TanDEM-X 90m DEM tiles provided by DLR (Rizzoli et al., 2017). This is a product variant of the original 12m (0.4 arcsec) DEM product with reduced pixel spacing (Wessel et al., 2018). The extent of the DEM is equal to the IV product. Further information can be found in the TanDEM-X Ground Segment DEM Products Specification Document (Wessel, 2018).
1.2.2.3. Land/Ocean Mask
This is an ESRI-shapefile of the ocean/land boundary, based on recent optical imagery and used for masking the retrievals in the ocean (usually coming from tracked sea ice). The shapefiles are internal datasets which are updated as required.
1.3. Algorithms
1.3.1. Coherent and Incoherent Offset Tracking using Sentinel-1
Offset tracking (OT) is the principal method for generating velocity fields on glacier surfaces from Sentinel-1 operationally. OT requires less operator interaction than interferometric techniques (such as InSAR), i.e. can be more easily and reliably automated. Another advantage of the OT technique over InSAR is the provision of two components of the velocity vector (along track and line-of-sight, LOS) from data of a single track. Also, the Greenland Ice Sheet is subject to highly variable meteorological conditions, resulting in rapid temporal decorrelation of the radar signal phase impairing the comprehensive application of DInSAR. OT is less sensitive or insensitive to temporal decorrelation of the radar phase signal, albeit at a lower accuracy. In areas with distinct and stable surface features, as common on mountain glaciers and outlet glaciers of ice sheets, coherence of the repeat-pass SAR data is not required because the cross-correlation of image-templates is based on the amplitude signal. In the level interior sections of ice sheets, without distinct amplitude features, temporal stability of the target phase (coherence) is required for speckle tracking. Loss of stability of the target phase (temporal decorrelation) is less critical than for the InSAR method because a continuous path is not required for deriving velocities of an isolated patch in an image, as is the case with InSAR.
The ice velocity (IV) algorithm is developed as part of the ESA CCI and Austrian Space Applications (ASAP) Programmes. The IV system's core module performs coherent and incoherent offset tracking (OT), utilizing long stripes of Sentinel-1 data acquired in interferometric wide swath (IW) mode. OT refers to several related methods that include speckle tracking, coherence tracking and amplitude tracking or feature tracking. Feature tracking uses cross-correlation of image patches to find the displacement of surface features such as crevasses or rifts and edges, that move with the same speed as the ice and are identifiable on two co-registered amplitude images, and subsequently ice flow velocity. In coherence tracking, the offset which maximises the interferometric coherence within a certain window size is determined and used to derive the ice velocity.
For Sentinel-1 Terrain Observation by Progressive Scans (TOPS) mode acquisitions, OT algorithms need to support burst handling and stitching. Also, the azimuth-varying phase introduced by the TOPS beam steering must be accounted for, both for complex offset tracking and for incoherent offset-tracking, where a factor of two oversampling of the SLCs is required prior to detection.
OT accuracies depend on the level of coherence and on the correlation window size (Bamler and Eineder, 2005) as well as on the type of features that are being tracked. Ionospheric scintillations can have a large impact on offset-tracking by causing azimuth shifts that are introduced by the fluctuating electron density along the sensor path. Misregistration is not a problem, but the large azimuth shifts are interpreted as ice motion and observed offsets can exceed the azimuth pixel size. DEM errors have less impact on the algorithm performance compared to DInSAR, since there is no phase unwrapping. However, since offset-tracking can be applied with longer temporal baselines, and even though larger spatial baselines can be used for offset-tracking, in the typical scenario the sensitivity to DEM errors is much smaller.
For further details the reader is referred to the Algorithm Theoretical Basis Documents of the Greenland and Antarctic Ice Sheet CCI projects (Khvorostovsky et al. 2018; Nagler et al. 2018) and Nagler et al. (2015).
1.3.2. Generation of Regional Ice Velocity Maps from Multiple Tracks using OT
The generation of a regional ice velocity map requires the combination of offset tracking results, from several tracks, each providing slant range and azimuth displacements at local incidence and heading angle. Figure 1 illustrates the high-level processing line for IV production. The system includes 3 main modules, which are detailed in the System Quality Assurance Document, but briefly explained here.
|
Figure 1: High-level flow chart of the IV processing system. Green – input data, Blue – processing modules, Red - product and intermediate products, Yellow – product data base.
Within the IV module SAR data and orbit data are imported into the system and velocity maps are generated for pairs of repeat pass data of the same track. The long tracks are geocoded into the common map projection of the output grid, providing observations of slant range displacement and azimuth for each track separately. The local incidence and heading angle are calculated using the annotated image and orbit parameters and a DEM as input.
The Merge module combines all IV products from all tracks and image pairs over a specified time span (i.e. 1 year), applying a weighted average approach. The combination of displacement observations of multiple tracks applies a weighted least squares model to fit the multiple observed displacements according to:
where x is the horizontal velocity vector (in local Easting/Northing coordinates), y is a vector with the observed velocities (in slant range/azimuth), A is the matrix projecting the horizontal velocity to slant range/azimuth geometry, and 
is a noise vector. The matrices are defined as:
where f() is the master, t() is the template image at the pixel position (u,v). The peak in the correlation matrix corresponds to the best matching shift of the slave window with respect to the master window.
The weighted least squares estimate of the horizontal velocity vector is then:
An advantage of the approach above is that if crossing tracks (ascending/descending) are used, the displacement measurement in the direction with higher resolution dominates in the velocity retrieval. On the other hand, it allows to reduce e.g. ionospheric azimuthal artefacts by adjusting the weights for azimuth shifts. In the case of only one track (single range/azimuth offset-tracking pair), the proposed method provides an exact inversion, using the assumption of surface-parallel flow by incorporating the local height gradient derived from the DEM. The standard deviation error estimate is propagated in this case.
Lastly the Validation module facilitates quality assessment of the IV products, by automating various standard validation tests, including internal consistency checks and intercomparisons with independent data sets (e.g. in-situ GPS, or other published velocity maps). The output is statistical information on the intercomparisons compiled in a quality assessment report.
1.3.3. Filtering and interpolation
For various reasons the OT algorithm sometimes fails to find matching features or erroneous matches are returned leading to gaps and/or errors in the velocity fields. This can for example be caused by a lack of traceable surface features or low coherence due to melt. For smoothing and removing outliers a 3x3 simple median filter is applied, after which a 9x9 first order plane fit filter is applied to fill in small data gaps. Further filtering/gap filling is left to the user if required.
1.3.4. Uncertainty characterisation
For incoherent offset-tracking (the general case), the uncertainty is characterised for each pixel by calculating a local offset-map standard deviation in a 5x5 neighbourhood. A plane fit to the offset map in the 5x5 neighbourhood is subtracted prior to calculating the standard deviation, so that an actual velocity gradient is not interpreted as a noise signal. The standard deviation estimate is corrected for any averaging carried out, as well as correlation between neighbouring samples (i.e. if the radar data are oversampled). The Greenland Ice Sheet velocity map is accompanied by its associated error standard deviation. The latter is also a map, in the same geometry as the associated measurement, providing a measure of uncertainty on a per-pixel basis.
1.4. Output data
The ice velocity product covers the Greenland Ice Sheet and is provided as a NetCDF file with 7 separate layers/maps representing different variables: the 3 velocity components: vx, vy, vz as well as vv (the magnitude of the horizontal components), and gridded maps showing the valid pixel count and the uncertainty/standard deviation (easting & northing) (Table 3). The ice velocity map is annually averaged and provided at 250m grid spacing in North Polar Stereographic projection (EPSG: 3413). A new NetCDF file of annually averaged IV is produced each year. The horizontal velocity is provided in true meters per day, towards easting (vx) and northing (vy) direction of the grid, and the vertical displacement (vz), is derived from a digital elevation model. More details are given in Section 1.3 of the C3S Product User Guide for Ice Velocity [D2].
Table 3: IV main data variables.
Variable name | Variable description |
land_ice_surface_easting_velocity | Ice velocity East component [m/day] |
land_ice_surface_northing_velocity | Ice velocity North component [m/day] |
land_ice_surface_vertical_velocity | Ice velocity Vertical component [m/day] |
land_ice_surface_velocity_magnitude | Ice velocity magnitude [m/day] |
land_ice_surface_measurement_count | Valid pixel count [#] |
land_ice_surface_easting_stddev | Standard deviation easting [m/day] |
land_ice_surface_northing_stddev | Standard deviation northing [m/day] |
References
Bamler, R., M. Eineder (2005): Accuracy of differential shift estimation by correlation and split-bandwidth interferometry for wideband and delta-k SAR systems, IEEE Geosc. Rem Sens. Lett. 2, 151-155.
Khvorostovsky, K. et al., (2018): Algorithm Theoretical Baseline Document (ATBD) for the Greenland Ice Sheet CCI Project of ESA's Climate Change Initiative, version 3.2
Nagler, T, Rott, H., Hetzenecker, M., Wuite, J. and Potin, P., (2015): The Sentinel-1 Mission: New Opportunities for Ice Sheet Observations. Remote Sens., 7, 9371-9389.
Nagler, T. et al., (2018): Algorithm Theoretical Baseline Document (ATBD) for the Antarctic Ice Sheet CCI Project of ESA's Climate Change Initiative, version 3.0. Available from http://www.esa-icesheets-antarctica-cci.org/
Rizzoli, P., Martone, M., Gonzalez, C., Wecklich, C., Borla Tridon, D., Bräutigam, B., Bachmann, M., Schulze, D., Fritz, T., Huber, M., Wessel, B., Krieger, G., Zink, M., and Moreira, A. (2017): Generation and performance assessment of the global TanDEM-X digital elevation model. ISPRS Journal of Photogrammetry and Remote Sensing, Vol 132, pp. 119-139.
Wessel, B. (2018): TanDEM-X Ground Segment – DEM Products Specification Document", EOC, DLR, Oberpfaffenhofen, Germany, Public Document TD-GS-PS-0021, Issue 3.2. Available: https://geoservice.dlr.de/web/dataguide/tdm90/pdfs/TD-GS-PS-0021_DEM-Product-Specification.pdf
Wessel, B., Huber, M., Wohlfart, C., Marschalk, U., Kosmann, D., & Roth, A. (2018): Accuracy assessment of the global TanDEM-X Digital Elevation Model with GPS data. ISPRS Journal of Photogrammetry and Remote Sensing, 139, 171–182. doi:10.1016/j.isprsjprs.2018.02.01
