Global seasonal Sentinel-1 interferometric coherence and backscatter data set

Sentinel-1 data selection

The Copernicus Sentinel-1 mission was launched by the European Space Agency (ESA) in 2014 with the Sentinel-1A satellite, complemented with the second Sentinel-1B satellite in 2016. Each satellite has a 12-days repeat cycle. Continuity of the Sentinel-1 mission has been approved by ESA until 2030 and replacement satellites will be launched. The satellites operate in different acquisition modes over different parts of the globe. Land masses are covered primarily by the Interferometric Wide-Swath Mode (IW) with a 250 km swath width across-track. Single-look-complex (SLC) Level 1.1 data are required for interferometric processing. Along-track, Sentinel-1 data are sliced into consecutive frames (slices) of about 250 km length. Data are distributed via ESA’s Scientific Sentinel-1 Hub, which is mirrored at NASA’s Alaska Satellite Facility DAAC (ASF-DAAC). During production, Sentinel-1 SLC data were accessed on the ASF-DAAC data repository which resides on Amazon’s AWS S3 bucket in region us-west-2.

Sentinel-1 satellites cover various parts of Earth in ascending and descending flight direction in a total of 175 relative orbits. ESA’s flight plan has some areas covered every six days and in both flight directions, predominantly over Europe. For the production of this data set, Sentinel-1 SLC frames were selected from all available scenes acquired between December 1^st 2019 and November 30^th 2020. Over the one-year timeframe, a maximum of 30 to 31 acquisitions at 12-days repeat, and 60 to 61 acquisitions at 6-days repeat intervals can be expected. The following selection criteria were applied consecutively to achieve global coverage with uniform distribution of acquisitions across seasons (Fig. 1):

• Global descending data (Fig. 1a) were selected where the one-year stack size had at least 25 acquisitions.

• Spatial gaps were filled with ascending data (Fig. 1a) where the one-year stack size had at least 25 acquisitions.

• For spatial consistency, over conterminous North America north of Panama, preference was given to ascending data where both ascending and descending data existed with stack sizes over 25 acquisitions.

• For stack sizes less than 25 acquisitions, preference was given to the flight direction with the larger number of acquisitions.

• Remaining gaps were filled with data from the flight direction available.

Arctic and Antarctic regions are typically covered with polarization modes of horizontal transmit (HH single- or HH/HV dual-polarization). Figure 1b shows the global distribution of the processed data in horizontal and vertical polarization transmit modes, respectively. Table 1 summarizes the number of selected scenes in the two flight directions and various polarization modes. The total number of processed Sentinel-1 SLC frames came to ~205,000 scenes with a total raw input data volume of about 850 Terabytes. Figure 1c,d show the spatial distribution of the final scene selection with the number of 6- and 12-days repeat-pass image pairs. Consistent 6-days repeat coverage with about sixty image pairs from either ascending or descending orbits could be processed over Europe, the coastal areas of Greenland and Antarctica, and some smaller areas around the world; note that in some regions (e.g., India, interior Greenland, Northern Canada, Eastern China) 6-days repeat coverage was available in certain seasons only (Fig. 1c). A consistent coverage with 12-days repeat-pass imagery, instead, could be processed almost globally with the nominal maximum of about thirty repeat-pass pairs in areas where only one satellite, Sentinel-1A or Sentinel-1B, acquired data in all but few areas above 60° N in Canada, Greenland, or Russia (Fig. 1d). In some small areas in the Midwestern United States, the Khabarovsk region in Far-Eastern Russia, or in the Northern Sahara, neither Sentinel-1A nor Sentinel-1B acquire data in IW mode, leading to small gaps in the final data set.

Processing approach

The overall processing workflow was developed based on the interferometric processing software developed by GAMMA Remote Sensing and geared towards efficient processing in the Amazon Web Services (AWS) cloud environment utilizing Earth Big Data LLC’s cloud scaling solutions. The workflow is divided into three main blocks as illustrated in Fig. 2. Sentinel-1A and -1B acquire data along 175 relative orbits/orbital tracks. Blocks 1 and 2 were processed on a per relative orbit basis; block 3 was initiated after blocks 1 and 2 had been completed for all relative orbits.

Processing Block 1

For each SLC of a given relative orbit, processing block 1 entailed:

1. 1.

   Conversion of SLC image files to a GAMMA specific format. Each Sentinel-1 SLC, covering an area of ~250 × 250 km, consists of six SLC image files (one SLC image file for each of the three sub-swaths in co- (VV or HH) and cross-polarizations (VH or HV).

2. 2.

   Compensation of the SLC amplitudes for the noise equivalent sigma zero (NESZ).

3. 3.

   The orbit state vectors provided with the original Sentinel-1 SLCs were updated with the precision state vectors (AUX_POEORB) distributed by the Sentinel-1 payload data ground segment 20 days after data take with a precision (3σ) generally of the order of 1 cm (target requirement < 5 cm).

4. 4.

   Each Sentinel-1 sub-swath SLC typically comprises nine to ten (sometimes more) individual bursts separated by a few no-data lines, a format that we will from now on refer to simply as “burst SLC”. For each Sentinel-1 burst a unique burst ID number was defined that permits identifying the corresponding bursts in all repeat acquisitions. We then extracted individual bursts in each sub-swath and wrote each individual burst image and associated parameter files to a temporary AWS S3 bucket.

Processing Block 2

Once all SLCs of a given relative orbit have been imported and individual bursts have been stored to AWS S3, the pre-processing block for generating terrain-corrected and geocoded coherence and backscatter images was initiated. The processing block started with grouping all bursts in a selected relative orbit into 1 × 3 burst groups/segments (3 sub-swaths comprising 1 burst each) and the creation of consistent stacks of burst SLCs for each 1 × 3 burst segment across all multi-temporal Sentinel-1 observations of the same relative orbit per segment. All processing steps described below were then performed on a segment-by-segment basis. The main processing parameters have been summarized in Table 2.

Having created consistent stacks of burst SLCs for all repeated observations from the same relative orbit, we defined for each burst SLC segment a reference (primary) acquisition to which all other (secondary) burst SLCs of the same segment were co-registered. In the next step, transformation functions (so called geocoding look-up tables) for resampling between slant-range and map geometry were calculated based on the available orbit information and the 3 arcsecond Copernicus DEM³⁴ for the burst SLC segment selected as reference. Alongside the geocoding look-up table, maps of the local incidence angle and layover/shadow masks were produced. An intensity cross-correlation method was subsequently applied to check the geometric accuracy of the transformation functions³⁵. This accuracy check uses a simulated SAR backscatter image, calculated by means of the DEM and an assumed topography dependence of the backscatter³⁶, as reference geometry against which the matching of the Sentinel-1 intensity images in small image chips was performed. The standard deviation of offsets in range and azimuth with respect to a polynomial regression fit to the offset estimates served as quality control. Throughout the entire processing chain, the geocoding look-up table served two purposes:

1. 1.

   Resampling of the DEM in map projection (geographic, EPSG: 4326) to the slant-range geometry of the reference SLC image to be able to i) consider topography related offsets between reference and secondary SLCs in the co-registration step, and ii) simulate the topographic phase when computing the differential interferograms,

2. 2.

   terrain-corrected geocoding of all imagery that matches the geometry of the selected reference acquisition/segment (i.e., all coherence and backscatter intensity images derived from SLCs that were co-registered to the selected reference).

The co-registration of the burst SLCs was done based on transformation functions (referred to as co-registration look-up tables) derived from the orbit information and the DEM resampled to the slant-range geometry of the selected reference burst SLC. The accuracy of the orbit information available for Sentinel-1 suggests that these transformation functions allow for resampling the secondary burst SLCs to the geometry of the reference with a co-registration accuracy of few hundredth of a pixel, i.e., sufficient to preserve coherence⁵. However, to ensure accurate and reliable co-registration, we also applied a single matching refinement (intensity cross-correlation) that served, if necessary, as correction of the co-registration look-up tables and to document the quality of the co-registration.

For simulating the baseline-dependent interferometric phase introduced by topography, the DEM in slant-range geometry and an orbit-based phase model was used. In the calculation of the differential interferometric phase (in SAR geometry), the topographic phase was subtracted and common band filtering in range was applied. Differential interferograms were calculated for all 6-, 12-, 18-, 24-, 36-, and 48-days repeat-pass image pairs in the co-registered stack of images. This was done in a multi-look burst geometry, meaning that bursts were kept separate so that no data from different bursts were combined in the estimation. The differential interferograms were calculated using 12 range and 3 azimuth looks resulting in interferograms with a ground-range resolution of about 40 m. The number of interferograms varied dependent on whether 6- or 12-days repeat-pass imagery was available. In case of a complete 6-days coverage throughout the observation period (Fig. 1c), a total of the order of 340 interferograms could be computed. In case of a consistent 12-days coverage (Fig. 1d), the number of interferograms was reduced to about 110.

An adaptive coherence estimator was used to produce coherence maps from all differential interferograms in multi-look burst geometry. The adaptive moving window estimator derives an initial estimate using a fixed 3 × 3 pixels estimation window size. Based on the initial estimate, a spatially adaptive estimator was applied so that the estimation included 5 × 5 or 7 × 7 pixels in multi-look geometry for low coherence levels. In the final coherence estimation, inverts of the deviations of the 3 × 3 pixels coherence estimates located within the larger estimation window served as weights to preserve spatial detail in the final estimate. Given the multi-looking and the size of the adaptive estimation window, coherence was estimated with a nominal number of looks of at least 324 and up to 1764. The large number of looks ensured a minimization of the bias inherent to coherence estimation². In areas for which decorrelation is expected to be complete (e.g., water surfaces, layover), the average coherence was of the order of 0.03. Finally, to produce seamless coherence images in slant-range geometry without no-data lines between subsequent bursts, coherence estimates for the three bursts per segment were mosaicked (referred to as de-bursting in ESA terminology). Only one of the available estimates in regions of overlapping bursts/sub-swaths was used. Averaging these would have resulted in different statistics between areas with or without burst overlap. All coherence images were eventually geocoded to the required map geometry with a 3 arcsecond sampling using the geocoding look-up table calculated for the reference acquisition.

Coherences were processed at co-polarization only (VV or HH dependent on data availability). The processing of coherence for cross-polarization data was not in the scope of the project. In the case of backscatter intensity, both polarizations acquired by Sentinel-1 in IW mode (VV/VH or HH/HV) were instead considered. Alongside the coherence estimation, we radiometrically terrain corrected (RTC) backscatter images acquired at co- and cross- polarizations. The pre-processing of backscatter intensity images deviated from the processing of coherence following the co-registration. The co- and cross-polarization burst SLCs were detected and multi-looked with multi-looking factors of 36 in range and 9 in azimuth. The multi-look backscatter images in burst geometry (calibrated to σ⁰) were subsequently de-bursted to create seamless images in slant-range geometry combining bursts from all three sub-swaths. In order to calibrate the backscatter images over sloped terrain, the varying extent of the ground area covered by image pixels dependent on topography was calculated as with Frey et al.³⁶. For each image, a normalization factor for the topographic correction of backscatter was calculated with the ratio of the pixel area referring to flat terrain (ellipsoid referenced σ⁰ areas) and the “true” pixel area considering topography (γ⁰ normalization area to produce “terrain-flattened” γ⁰ backscatter images³⁷. All backscatter images in co- and cross-polarizations were geocoded using the same geocoding look-up table used for geocoding the coherence imagery.

All coherence, backscatter, and ancillary imagery (local incidence angle maps, layover/shadow maps) were fully pre-processed for all 1 × 3 burst segments in a relative orbit, all resulting images were tiled to a fixed global 1° × 1° grid and stored as GeoTIFF files in a temporary AWS S3 storage bucket together with quality reports generated during processing. The quality/ancillary reports associated with each image comprised information on i) perpendicular baselines, ii) quality indicators for the co-registration, and iii) quality indicators for the geocoding accuracy.

Processing Block 3

Having pre-processed all coherence and backscatter imagery and created a complete global tile data base, the final processing block comprised i) seasonal compositing of coherence and intensity imagery, ii) fitting of a coherence decay model to seasonal stacks of coherence at pixel level, and iii) creation of global mosaics.

For each tile, seasonal composites of the coherence at different repeat intervals and backscatter imagery were calculated. We calculated the median coherence based on all coherence estimates per tile of a given repeat interval (6, 12, 18, 24, 36, and 48) per three-month period: 1) December, January, February 2) March, April, May 3) June, July, August, and 4) September, October, November. We chose the median operation to account for outliers. In the case of the backscatter intensity products, we calculated per three-month period the average backscatter intensity in VV and VH, or HH and HV, polarization.

The decay of coherence with increasing repeat interval was modelled for each season at pixel-level with the exponential model³⁸

where ρ_∞ and τ denote the long-term coherence and rate of coherence decay with increasing repeat interval, respectively. The estimation of the coherence model parameters ρ_∞ and τ was done based on the median-aggregated seasonal coherence estimates. The pixel-level estimates for ρ_∞ and τ, together with a quality indicator of the model fit, i.e., root mean square difference between modelled coherence and individual non-aggregated coherence observations, were stored in form of images in the same geometry as the tiled Sentinel-1 datasets. Curve fitting was achieved with the standard Levenberg-Marquardt least-squares regression unless the initial model fit resulted in non-physical solutions for the model parameters, i.e., negative estimates for ρ_∞ (cf., Technical Validation). In such cases, we adopted a different regression technique with the Trust Region Reflective regression³⁹, which is significantly slower (and thus not applied globally by default) but allows for constraining the potential range of parameter estimates, i.e., 0 ≤ρ_∞≤1.

Finally, we created global mosaics of each layer, i.e., seasonal 6-, 12-, 18-, 24-, 36-, and 48-days coherence at VV- and HH polarizations, seasonal backscatter in co- (VV and HH) and cross-polarizations (VH and HV), as well as ρ_∞, τ, and the corresponding model RMSE of the model fit (see Data Records).

Source: Ecology - nature.com