1. Introduction
Radiometrically terrain-flattened synthetic aperture radar (SAR) imagery is widely considered to be the first SAR-derived analysis ready dataset (ARD) family [
1] and is of considerable interest to the sustainability, ecosystem and agriculture science communities. The Sentinel-1 SAR constellation [
2], part of Europe’s Copernicus Earth Observation Programme, has significantly increased adoption of SAR data within geospatial data frameworks, including data cubes, that were originally developed for optical imagery. Terrain-flattened SAR ARD datasets from other SAR missions like ALOS-2 and NovaSAR are also under active development, see
https://ceos.org/ard. A large number of national agencies and commercial entities now offer Sentinel-1 backscatter imagery as a foundational service on which to build applications, e.g., [
3]. These ARD datasets are expected to support a wide number of applications including land cover classification, multi-temporal change detection etc.
Most terrain flattening processors implement the Gamma Flattening approach [
4] or related variants that have been optimized for performance, e.g., [
5]. The Gamma Flattening approach operates on Level-1 SAR imagery in slant-range or ground-range coordinate systems and generates geocoded, radiometrically terrain-flattened imagery. The rest of this manuscript assumes familiarity with the Gamma Flattening formulation [
4], as we rely heavily on associated terminology and algorithm description.
1.1. Terminology
Before proceeding, we describe some of the key terms used in this manuscript. We attempt to stay consistent with [
4] and with terminology used by the Sentinel-1 toolbox [
6] user community.
A geocoded product or
Geocoded Terrain Corrected (GTC) product, as they are referred to in SAR mission product user guides, is derived by precisely geolocating SAR imagery using a Digital Elevation Model (DEM) [
7]. Alternately, the GTC products themselves could be generated directly by focusing raw radar pulses onto a regular map grid [
8]. When starting from Level-1 products, SAR imagery is usually calibrated to
before it is interpolated on to a regular map grid. GTC products calibrated to
or
are also common. GTC products have been geolocated precisely [
9] but have not been corrected for terrain-related radiometric effects. While we provide mathematical expression for working with different calibration levels of GTC products, we specifically focus on
GTC products, which we process at Descartes Labs at a global scale [
10].
A terrain flattened or
Radiometrically Terrain Corrected (RTC) product [
4] or normalized radar backscatter (NRB) [
11] product, is a special type of GTC product where the imagery has been corrected for terrain-related radiometric effects. In the context of this manuscript, we always assume that an RTC product has been calibrated to
(Table I of [
4]). RTC products are widely considered to be the most ready-for-analysis product derived from SAR imagery and most similar to optical imagery for developing similar applications [
1,
11].
In general, a collection of GTC products generated on the same map grid is referred to as a
geocoded stack. In the context of this manuscript, we specifically refer to GTC products generated on a common grid from interferometrically compliant acquisitions as a geocoded stack, unless mentioned otherwise. Such products are usually labelled with a common Path-Frame identifier (ERS, ALOS etc) or unique burst identifiers (Sentinel-1) [
10,
12]. These identifiers represent unique imaging geometry configurations, i.e., all images in the collection share baselines of less than a few kilometers with respect to each other and are acquired at similar incidence angles.
1.2. Manuscript Organization
This manuscript is organized as follows:
Section 2 describes our approach to transforming GTC SAR imagery to RTC SAR imagery.
Section 3 describes our time-series simulation experiments using actual SAR metadata and provides justification for the use of static radiometric terrain flattening factors for transforming stacks of GTC products to RTC products. In
Section 4, we discuss layover in SAR imagery in the context of radiometric terrain corrections and its implications for the proposed approach. In
Section 5, we describe a simple method to compute global scale terrain flattening factors using the Sentinel-1 toolbox [
6]. In
Section 6, we discuss the implications of applying radiometric terrain corrections to geocoded SAR data instead of applying radiometric terrain corrections to SAR data in its original radar geometry.
2. Revisiting the Gamma Flattening Formulation
In this section, we reinterpret the Gamma Flattening formulation [
4] in the context of adopting the method to apply it directly to GTC products. We try to use the same terminology as presented in [
4] to allow readers to relate to our interpretation more easily. We ignore the issue of heteromorphism (layover) in this section, as we will address this in great detail in
Section 4.
2.1. Single DEM Facet
We will consider a single triangular facet of a DEM represented by points
,
and
, with their projections onto the reference ellipsoid represented by points
,
and
. Let
represent the local incidence angle of the facet at
, the centroid of
.
is a point on the reference ellipsoid that projects to the same point in the slant range geometry as
and
represents the nominal incidence angle at
.
,
and
represents the projection of the DEM facet on to the slant range plane and
represents the projection angle [
13], i.e, the angle between the normal to the DEM facet and the normal to the slant range plane.
Figure 1 depicts the geometry and is similar to a combination of Figures 2 and 5 of [
4], except that the points
,
and
in Figure 5 of [
4] represent regular grid points on a map for GTC products.
For the case where there exists a one-to-one mapping between pixels in the original radar image projection and a well known map projection, the conservation of energy argument, as presented in [
4], and the area relationships from
Figure 1, show that
We make the following observations about the formulation presented above:
Equation (
1) can be used to flatten GTC products corresponding to any of the standard calibration levels -
,
and
.
The formulation can be applied to GTC products in any well known map projection system [
14] as long as the actual area computations are performed in a geocentric cartesian projection system, e.g., EPSG:4978, to avoid projection system related distortions.
Since the transformation of GTC products according to Equation (
1) only involves computation of simple facet-by-facet area normalization factors (assuming no layover), we can significantly speed up processing and circumvent the use of large radar image index lookup tables.
If points
,
and
lie on the reference ellipsoid,
and
are complementary (see Equation (15) of [
5]) and Equation (
1) reduces to the classic Equations (2) and (3) of [
4], i.e.,
The Calculation of area normalization factors in Equation (
1) can be easily implemented within any open SAR processing software [
6,
15,
16] that includes functions for interpreting SAR imaging geometry, in combination with a map projection transformation library like PROJ [
14].
2.2. Extension to Rectangular Pixels
Extending Equation (
1) for use with a rectangular GTC pixel using two-facets [
4] or four-facets including the center [
5] is straightforward using area summation of facets visible to the radar (Equation 11 of [
4]).
where
represents all facets associated with the GTC pixel of interest and
represents a subset of
that is visible to the imaging SAR sensor. A facet is considered visible to the SAR sensor if it is not in radar shadow and
.
is usually set to
, consistent with the area ratio threshold in Section II-H of [
4] or more conservatively to
, e.g., in [
17]. Note that this threshold is to avoid amplification of the noise due to the
term in the denominator. We don’t need similar thresholds for the
term in the numerator but we do use it’s absolute value as projection angles can be greater than
when the facet itself is impacted by foreshortening. In general,
and the look vector to the center of the pixel could be reused with all the contributing facets to further reduce the number of computations involved. Ref. [
5,
18] described a similar formulation but used a constant
and
for the entire rectangular pixel in slant range coordinates. The formulation presented here operates on geocoded pixels and accounts for contribution from individual DEM facets.
3. Terrain Flattening of Geocoded Stacks
In
Section 2, we described an approach to terrain-flatten a single pixel of a GTC product. In this section, we analyze the sensitivity of the area flattening factor from Equation (
1) to the variation of imaging geometry within a stack for actual SAR acquisitions from sensors like Sentinel-1 and ALOS-1. In this section, we again consider a single triangular DEM facet and study the effect of variations in
,
and
. The backscatter term in Equation (
1) is obtained from the source GTC product that we want to transform and the area term
is independent of the imaging geometry. For all examples presented in this section, we consider a facet located at near-range of the imaged swath as the impact of change in imaging geometry decreases with slant range [
19]. We also specifically consider the transformation of
GTC products to
as this is most relevant for use with our global scale radar backscatter product [
10]. This analysis can be easily extended in the same framework to study transformation of
and
GTC products.
3.1. Sentinel-1
We consider a stack of 58 acquisitions corresponding to a single Sentinel-1 burst footprint with European Space Agency (ESA) identifier IW1-0151226 (Descartes Labs identifier [
10] 071-2637-IW1-VV-RD) over Big Bear, California in the USA and spanning the time period of January 1, 2020 to January 1, 2021. The perpendicular baseline variation of this stack is about 200 meters. This footprint is imaged from a right looking descending geometry and we set up a single DEM facet as follows:
where
and
represent the Easting and Northing coordinates of a nominal near-range pixel of the GTC product in a Universal Transverse Mercator (UTM) system. We vary
in the interval
meters and
in the interval
meters, to study the impact of imaging geometry in the relationship between
and
. We note that the altitude of point
is fixed to zero for the simulations presented here and moving the entire facet up or down by a constant within the limits of earth’s topography (-500 to 9000 meters) only modifies the mean incidence angles for space-borne missions, and does not affect the interpretation of our results.
Figure 2 shows that, for this stack, for most local incidence angles the variation is well below 0.01 dB and we only start approaching 0.02 dB variation for local incidence angles greater than 85 degrees. The bright strip in
Figure 2(Left) corresponds to
and is an numerical artifact from estimating a ratio of two small numbers. When interpreted in the context of Equation (
3), facets with
would contribute very little to the sum in the numerator. The perpendicular baseline spread in this stack is comparable to other Sentinel-1 stacks due to the mission’s narrow orbital tube [
2]. It is clear that the observed variation is much smaller than the typical radiometric requirement of 0.1 dB associated with change detection applications, even for fairly steep terrain, i.e over a wide range of
and
. Consequently, we can conclude that a constant pixel-by-pixel flattening factor (per imaging geometry) is sufficient to transform
GTC products to
for Sentinel-1 stacks. We extend this argument to suggest that constant pixel-by-pixel flattening factor (per imaging geometry) is sufficient to efficiently transform GTC products to RTC products from other SAR missions with narrow orbital tubes like ALOS-2 and NISAR as well.
Another interpretation of the result above is that the area projection of DEM facets do not vary significantly between passes for missons with narrow orbit tubes. But the fact that each Level 1 SAR product is distributed in its own projection system [
10] requires users to expend computational resources and build elaborate lookup tables, from these projection systems to well known map coordinate systems, to carefully account for area contribution to each SAR product pixel. Our proposed approach of starting from a GTC product eliminates the need for mapping area projections for each product. A static terrain flattening factor layer can be built using a reference orbit and this correction factor can be reused for all GTC SAR products acquired from a similar imaging geometry. Additionally this framework allows us to estimate error introduced by using static terrain flattening factors using actual orbit data.
Navacchi
et al. [
19] also observed that the variation in projected incidence angle is on order of 0.005 degrees (standard deviation) for Sentinel-1, which is consistent with our simulations, and proposed the use of static correction factors to transform
to projected local incidence angle (PLIA) normalized backscatter product. With our simulations, we are able to show that the same approach can be used for generating
RTC products on the fly as well.
3.2. ALOS-1
While we considered a narrow orbital tube mission in
Section 3.1, here we consider a stack of 34 ALOS-1 acquisitions corresponding to Path 216, Frame 740 over Long Valley, California in USA and spanning the time period of June 1, 2006 to March 1, 2011 and exhibiting a perpendicular baseline variation about 6500 meters. This frame is imaged from a right looking ascending geometry and we set up a single DEM facet as follows:
where
and
represent the Easting and Northing coordinates of a nominal near-range pixel of the GTC product in a Universal Transverse Mercator (UTM) system. Note the change in Easting of
to accommodate the change in the pass direction.
Figure 3 shows that for this stack, the variation in area normalization factor in Equation (
1) can be as high as 0.3 dB which is larger than the typical radiometric requirement of 0.1 dB associated with change detection applications. Even in this case, for most incidence angles the observed variation is below 0.15 dB but we see significantly higher variation for incidence angles greater than 80 degrees. The band corresponding to
is also clearly visible and is broader than the Sentinel-1 case due to larger variation in imaging geometry. In general, we can deduce that a constant pixel-by-pixel correction factor per imaging geometry is insufficient to transform
GTC products to
for ALOS-1 stacks or missions with wider orbital tubes in general. In
Section 3.3 we present a solution for these types of missions.
3.3. Generalized Formulation
In this section, we present a generalized formulation, inspired by SAR interferometry (InSAR), that allows us to exploit interferometric baseline information with static terrain flattening terms to efficiently flatten GTC products from InSAR-capable SAR missions. Since, this manuscript focuses on modern SAR missions with narrow orbit tubes like Sentinel-1, ALOS-2 and NISAR, we only present the motivation for the formulation without delving into greater detail and analysis.
While we showed in
Section 3.2 that a constant pixel-by-pixel flattening factor is insufficient for flattening stacks with large baseline variations, we also observe that the form of Equation (
1) suggests that there might a relationship between interferometric baselines and
, akin to the topography phase term used in differential InSAR analysis, e.g., [
20], due to the dependence on the incidence angle term
and the projection angle term
. Following a similar Taylor-series expansion as Equations (8)–(10) in [
20] we can generalize that for each GTC pixel
where
represents the ratio computed using a reference imaging geometry (could be a reference orbit or a reference scene in the same stack) in decibels,
and
represent the perpendicular and along track baselines of the GTC pixel in product of interest w.r.t the reference imaging geometry, and
and
represent a constant scaling factor associated with the GTC pixel. The parallel baseline (
) does not contribute to Equation (
4) as this term does not modify the line-of-sight vector and its impact on
is minimal. The scaling factors,
and
, for each GTC pixel depend on the slope of facets associated with the pixel and can be easily computed numerically, contemporaneously with the computation of
, by recomputing
and
after perturbing the estimated reference satellite location by unit perpendicular and along-track baseline vectors. For modern sensors with Doppler steering capability, e.g., Sentinel-1, ALOS-2, TERRASAR-X etc, the along track baseline is on the order of few meters and the associated term (
) can be ignored in Equation (
4).
We demonstrate the linear relationship between
and perpendicular baseline (
) in
Figure 4 for 3-different sets of combinations of
and
in the ALOS-1 simulation from
Section 3.2. The minor deviations from the linear trend that we observe are due to the orientation of the facets and the resulting sensitivity to the along-track baseline, which can be as high as 600 meters for ALOS-1. While we demonstrate this relationship with a single triangular DEM facet, the idea can be extended to area summation of facets in Equation (
3). In general, we can deduce that in addition to a constant flattening factor computed using a reference imaging geometry, similar to
Section 3.1, an additional pixel-by-pixel perpendicular baseline-related scale factor (
in Equation (
4)) should be sufficient to transform
GTC products to
for ALOS-1 stacks. Using the perpendicular baseline-related scale factor with Sentinel-1 (
Section 3.1) or ALOS-1 (
Section 3.2) stacks increases the precision of radiometric correction to well below 0.005 dB, even for incidence angles greater than 85 degrees in our simulations. Similar baseline terms can also be used to extend the approach proposed in [
19] for use with other SAR missions. The rest of the manuscript will focus on using static terrain flattening factors with the Sentinel-1 mission.
4. Impact of Layover
Radar shadow and layover are two inherent limitations of side-looking SAR imaging systems, e.g., [
21]. Regions impacted by radar shadow have to be masked out from SAR imagery as no energy is scattered back to the SAR sensor from these areas and the values for these pixels in GTC products are essentially noise. As far as we can tell, there is no disagreement amongst the user community regarding the masking approach to radar shadow.
The original Gamma Flattening approach [
4] proposes a solution, for regions impacted by layover or heteromorphism. We attempt to layout the inherent assumptions of this approach and describe the impact of correcting for heteromorphic effects on an image-by-image basis on a stack of geocoded imagery on multi-temporal change detection applications. Specifically, we will make the case for complete masking of pixels impacted by layover compared to using the computationally expensive approach of tracking all facets in each individual SAR image.
4.1. Single SAR Image
Layover is a bigger challenge with airborne SAR data rather than spaceborne SAR data, which is acquired with nominal incidence angles of 23 to 45 degrees (in most modern spaceborne SAR sensors [
22]), due to the comparatively low altitude of the airborne sensing platform. For Sentinel-1, the fraction of data affected by layover for a region like the United Kingdom is below one percent [
23] and could be as high as ten percent when specifically looking at targeted areas of interest (AOI) with steep topography [
24]. The fraction of pixels that are impacted by layover or shadow in all corresponding Sentinel-1 imaging geometries is very low, e.g., [
25,
26].
The Gamma Flattening approach [
4,
5] rigorously tracks the topological relationship between map coordinates and slant range radar geometry to account for the many-to-one mapping of energy in map coordinates to slant range geometry. The total area of all contributing facets are accumulated before imagery in
is normalized to
. However, the starting point of our proposed approach is a GTC product and not an image in slant range geometry implying that we suffer from the issue clearly indicated in Section II-D of [
4]. Pixels impacted by layover clearly stand out as too bright in GTC products as observed in number of previous works, e.g., [
4,
27].
The conventional terrain flattening approach corrects for the layover effect under a very specific assumption of a distributed scattering mechanism. It is assumed that the scattering mechanism in all pixels on the map contributing to the same slant range pixel is the same and hence their relative contributions are proportional to the projected area. This assumption is not necessarily true in most real life scenarios. For example, consider the textbook case of layover in SAR imagery where the top of a hill and its base are at same slant range from the imaging platform. The conventional terrain flattening approach redistributes energy under the assumption that the material and scattering mechanisms are the same at all the geographic locations mapping to the same pixel in slant range geometry. While this approach definitely reduces the number of bright outliers and improves the histogram, as shown in [
4], it may not be correct.
4.2. Stack of SAR Images
Stacks of SAR imagery are often used in the context of multi-temporal change detection, e.g., [
28]. In such scenarios, careful handling of regions affected by layover, which cannot be corrected, become important to avoid false detections and to improve robustness of any analysis building on SAR backscatter imagery [
3,
27]. In
Section 4.1, we highlighted that energy is redistributed proportional to the area contributions for layover regions in the Gamma Flattening approach. While this approach improves the histogram and visualization of a single image, it does not guarantee consistent redistribution of energy over time in a stack of RTC products. Our simulations in
Section 3 show that if the same facets contributed to a layover affected pixel in every image in the stack, their relative area contributions will be consistent. However, slight variations in satellite position between passes results in different facets contributing to layover affected regions at different time epochs, resulting in inconsistent redistribution of energy over time. In fact, SAR tomography techniques have the potential to coherently redistribute energy in layover regions by exploiting the variation in imaging geometry in a stack of SAR images but this topic is beyond the scope of this manuscript. The need for masking out layover regions has also been empirically observed by a number of research groups, across a spectrum of applications—e.g., [
29,
30].
For regions not impacted by layover, the assumption of distributed scattering mechanism could still be physically wrong, particularly in heterogeneous terrain like in urban areas. However, our simulations in
Section 3 show that terrain flattening can be interpreted as a pixel-by-pixel scaling within a stack and relative changes in time are preserved for such regions.
4.3. Shadow-Layover Mask
A number of shadow-layover mask generation methodologies have been developed, in both slant range geometry [
21] as well as in map coordinates directly [
27]. It is worth noting that the shadow-layover mask varies slightly between acquisitions in a Sentinel-1 stack due to variation in platform imaging positions. Having collectively analyzed a number of shadow-layover masks corresponding to the same imaging geometry for Sentinel-1, in our experience, it suffices to use the shadow-layover mask of a single reference acquisition and buffer it by 150–200m before applying it to the whole stack during change detection analysis. This buffering also accounts for imperfections in the source shadow-layover masks which are typically generated using a simple ray tracing method. Our approach to buffering the shadow-layover mask is similar to [
27], but is meant to compensate for variations of imaging geometry within the orbit tube in the stack of Sentinel-1 images and the original mask is derived in radar-coordinates using traditional approach. In this work, we do not delve deeper into the accuracy of different approaches of generating a shadow-layover mask but emphasize that we can significantly reduce amount of computation needed by exploiting the narrow orbital tube of Sentinel-1 and using a buffered shadow-layover mask from a reference acquisition or orbit for the whole stack. In the context of our proposed flattening approach in
Section 2 and Equation (
3), a GTC pixel should be considered as impacted by shadow or layover even if a single contributing facet is affected.
5. Experiments with Sentinel-1 Toolbox
In
Section 2 and
Section 3, we presented simulations using individual DEM facets and the results justify the use of static factors for terrain flattening Sentinel-1 imagery. In this section, we describe experiments we conducted with the Sentinel-1 toolbox [
6] to validate our simulations. We emphasize that we chose the Sentinel-1 toolbox because of its accessibility and the fact that it is widely used, allowing any interested reader to replicate these experiments. We believe that results from similar experiments with other open or commercial software will be of interest of the user community, especially if a large volume of RTC data is going to be generated in a systematic fashion with the software.
We used the following experimental setup with the Sentinel-1 toolbox (SNAP version 9.0.0). We staged Sentinel-1 Single Look Complex (SLC) products corresponding to the same imaging geometry (burst footprints) and replaced the imagery arrays in the TIFF files with a constant Digital Number (DN) value of 8000. We also staged the same global DEM that we use for processing our global backscatter and InSAR products [
10] and used it as an input to Sentinel-1 toolbox workflows. We used standard workflows to generate geocoded layers for calibrated
, terrain flattened
and shadow-layover masks for each of these products. The geocoded products were generated on a predefined 10 meter regular grid in a Universal Transverse Mercator (UTM) coordinate system. No thermal noise correction was applied as these experiments are designed to purely capture the impact of variation of imaging geometry within a stack. We then analyzed the statistics of the ratio of
to
using the geocoded products across the footprint and over time. We repeated the experiment with different values of the
oversamplingMultiple parameter ranging from one to four in the terrain flattening module to understand its impact on the processing errors, following the observations of [
5]. We expect to see variations in this ratio across pixels in line with our simulations in
Section 3.1, but we expect this ratio to be nearly constant in time.
As a constant DN was used instead of real SAR imagery, the resulting ratios should represent the consistency between SNAP and the method proposed here in the computation of projected facet areas and their summation. The other benefits of using constant DN imagery for this experiment include:
It eliminates the effect of differences introduced by InSAR-grade interpolators [
31] used in complex-value interpolation of SLC data and noisier bilinear or bicubic interpolators used with real valued intensity data in the workflow, and lets us focus on geometric inconsistencies.
A GTC
product derived from a constant DN image in slant range coordinates will also be constant valued image, thus allowing us to compare outputs with terrain flattened products generated from GTC products as described in
Section 2.
It eliminates the effects introduced by inconsistent spatial averaging due to the use of multilooking operator in slant range coordinates, as multilooked products of constant DN images are also constant valued. This effect is similar to phase closure artifacts observed in pair-by-pair InSAR analysis as described in [
10].
These different effects also contribute in their own way to the overall processing error budget but quantifying them is beyond the scope of this manuscript.
5.1. Open Ocean
The first example is a test over open ocean off the coast of California, USA - ESA identifier IW3-0151231 (Descartes Labs identifier 071-2655-IW3-VV-RD) with a stack of 38 acquisitions spanning the time period of Jan 1, 2020 to Jan 1, 2021. We picked an open ocean region as this would not be impacted by shadow-layover effects or the quality of DEM being used. We would expect our results to mimic our simulations in
Section 3.1. We observe a variation of ∼20 meters in the reference terrain height values in the individual burst metadata, which is negligible.
Figure 5 (Left) shows a clean peak around 0.04 dB for the histogram of pixel-by-pixel standard deviation of the ratio between
and
. This number is larger than the 0.01 dB (peak-to-peak) variation from the simulations in
Section 3.1. However, Sentinel-1 toolbox performance significantly improves when we increasing the DEM oversampling by a multiple of 2 in
Figure 5 (Right). The peak is observed to be around ∼0.005 dB and is well within acceptable limits for change detection applications. This experiment confirms the fact that using static flattening factors should be more than sufficient (< 0.1 dB) for regions not impacted by shadow-layover. This type of metric would be useful for the end user community to understand change detection sensitivity of operational processing workflows that will be deployed for generating various CARD4L normalized radar backscatter or ocean radar backscatter products [
11].
5.2. Rugged Terrain
We repeated the same experiment with our original burst footprint from
Section 3.1 over Big Bear, California. This region is characterized by steep terrain as well as some flat regions. The Sentinel-1 toolbox estimates only 1 percent of this footprint to be impacted by layover and that a negligible area is impacted by shadow. We masked out the estimated ratios with shadow-layover mask, using a buffer of 150m as described in
Section 4, before generating the multi-temporal statistics.
Based on our simulations in
Section 3.1, we expect the histogram from this region to look similar to that of the open ocean region above but with a slightly broader distribution as there are additional DEM interpolation operations involved which could introduce numerical noise. While we observe a nice sharp peak around 0.04 dB just like for open ocean, we also observe a long tail extending to about 0.3 dB in
Figure 6 (Left) with default settings. The standard deviation of the ratio for the majority of the image (83 percent) is below 0.1 dB. At an DEM oversampling multiple of 4,
Figure 6 (Right), we observe that the peak is centered around ∼0.025 dB and the standard deviation of 87 percent of the image is below 0.1 dB. However, the tail of the distribution still extends to 0.3 dB.
The pixels contributing to the long tail are correlated with rugged terrain. We include
Figure 7 over a particularly rugged region of the burst footprint to illustrate the correlation between observed standard deviation and the DEM, which indicates that this is likely a systematic processing artifact in Sentinel-1 toolbox. Another possible reason for the observed correlation is the underestimation of the areas impacted by shadow-layover. A number of additional software implementation factors can contribute to this observation including thresholds chosen for convergence of geometry computations, possible use of single precision representations for certain intermediate computations, handling of map projections etc.
It is clear from these experiments that oversampling of the DEM is critical for reliable terrain flattening, as observed by [
5]. The better terrain flattening results using higher DEM oversampling multiples are achieved at the expense of longer processing times and more compute resources. The area projection approach [
5] has been shown to have better properties than the bilinear weighting approach of [
4] in terms of compute resources and efficiency. In the future, we also plan to study the area projection implementation in ISCE3 [
16] in the presented framework. We also note that one of the advantages of our proposed approach in
Section 2 is that by using the DEM and GTC pixels on the same standard grid, it eliminates the need for additional DEM oversampling and produces consistent terrain flattening factors, as long as the GTC pixel is not impacted by layover.
5.3. Global Terrain Flattening Product
In [
10], we described our pipeline for generating near real-time SAR backscatter
GTC product from Sentinel-1 data. Following the approach described above and in
Section 3, one can also generate a global terrain flattening factor product at a granularity of a single Sentinel-1 burst with the following layers using Sentinel-1 toolbox:
Such a global terrain flattening factor product can be used to transform GTC products to RTC products on-the-fly using simple band math in standard geospatial data frameworks. Additional useful layers like local incidence angle, projection angle, line-of-sight angles and nominal incidence angles can also be included for use with other workflows on-the-fly. For consistency, the computation of these layers should be done with the same processing engine that was used to geocode or terrain correct the imagery. One of the big advantages of the proposed approach of using static flattening factors is that these factors can be easily re-estimated efficiently and cost-effectively at a future date with the same underlying SAR metadata and DEM if a better method were to be developed, without having to reprocess all the backscatter imagery data globally. Decoupling terrain flattening from terrain correction or geocoding allows us to support a wider-range of applications with SAR data more efficiently and cost-effectively. These global layers are significantly smaller in volume than the imagery archives and can easily be distributed openly for efficient use in cloud-enabled geospatial frameworks.
6. Discussion
6.1. Applicability of Terrain Flattening
While terrain flattening of SAR imagery is widely considered to be a prerequisite for a number of applications, there are a number of applications that don’t necessarily require it. We discuss some cases below:
Equations (
1) and (
3) clearly show that terrain flattening can be considered as a correction of a pixel-by-pixel bias term. Consequently, if analysis of individual SAR backscatter products can be reformulated as ratio of polarization channels—e.g, radar vegetation indices, the terrain flattening effects cancel out. Such analysis can be directly performed on GTC products.
Section 3.1 shows that the pixel-by-pixel bias is consistent for narrow orbital tube missions. Consequently, if multi-temporal backscatter analysis can be reformulated to work with relative changes w.r.t a reference epoch or a temporal average, terrain flattening effects cancel out. This is similar to using a reference epoch in InSAR time-series analysis.
Terrain flattened products from different imaging geometries—e.g., ascending vs descending passes, are not necessarily comparable over heterogeneous terrain like urban areas where the scattering mechanism is not necessarily distributed in nature. Comparing GTC products acquired from similar imaging geometries would allow for more sensitive change detections.
Multi-temporal, multi-modal change detection frameworks are increasingly becoming popular for wide area monitoring and change detection applications, e.g., [
28,
32]. These frameworks are designed to analyze time-series from multiple types of sensors and combine change detections. Sensitivity of change detection from SAR data can be improved by just considering different imaging geometries as a different sensor in such frameworks.
6.2. Efficient Processing
As of January 1, 2023 Sentinel-1 has imaged over 50 million bursts in Interferometric Wide Swath (IW) mode. However, the number of unique footprints that have been imaged is on the order of 340 thousand. We can reduce computation required for terrain flattening by a factor of ∼150, using static flattening factors. A new flattening factor for a burst footprint only needs to be computed when a previously unimaged burst is imaged by Sentinel-1. Use of static flattening factors allows terrain flattening to be implemented on-the-fly using simple band math within standard geospatial frameworks. Computational resources needed to generated a GTC product are also significantly less than those needed to generate a RTC product, thus reducing costs and resources needed to keep up with the live stream of Sentinel-1 data.
We also observe that using static flattening factors is the same as creating a stack of coregistered data with the Sentinel-1 toolbox before terrain flattening. This method uses the slant range geometry of the reference image of the stack to terrain flatten all the images after aligning them. The stack coregistration process is performed in slant range geometry in the traditional implementations and is computationally more expensive than the geocoding approach, which results in a coregistered stack in map coordinates as described in [
10]. Stacking is the approach recommended by Sentinel-1 toolbox developers to minimize geolocation inconsistencies and our proposed approach produces similar results while requiring lot fewer resources and computations.
6.3. Validation of Terrain Flattening Processors
Software implementation differences in terrain flattening processors, particularly over sloped terrain, have also been observed in previous comparative studies [
1,
33]. Similar comparative and validation studies are also needed for shadow-layover mask generation approaches, as masking is an important aspect for robust multi-temporal change detection applications with SAR data in general. The Gamma Flattening approach [
4] involves many more interpolation operations, both geometrically and over imagery, than our proposed approach (
Section 2) which only involves geometric interpolation. Identifying the exact source of possible discrepancies observed in
Section 5 will require more detailed comparative studies and development of more comprehensive synthetic tests and metrics. These topics are beyond the scope of this manuscript. We again emphasize that metrics similar to the ones presented in
Section 5 will be useful for the end user community to understand the limitations of the underlying terrain flattening implementation, open or commercial, that might be used to generate large datasets.
6.4. Analysis Ready Data Interoperability
Committee of Earth Observation Satellites (CEOS) is currently working on standardization of five different families of Analysis Ready Datasets (ARD) derived from SAR imagery (
https://ceos.org/ard/).
Normalized Radar Backscatter (NRB)
Interferometric Radar (InSAR)
Geocoded Single-Look Complex (GSLC)
Polarimetric Radar (POL)
Ocean Radar Backscatter (ORB)
These different families of products are often generated by independent processing chains and different Level-1 data as sources, e.g., NRB from Ground Range Detected (GRD) products and InSAR from SLC products. Lack of synchronization or cross-validation of these independent processing paths can lead to issues in interoperability of these datasets not limited to geolocation offsets. Our proposed approach of using static terrain flattening factors within the framework, implemented in [
10], ensures that all these families of products listed above are efficiently and consistently processed. In addition to these, our GTC
products [
10] can directly support sea ice, soil moisture and oceanography science users who do not typically work with terrain flattened SAR data. Our proposed approach of using static factors works equally well for transforming GTC products to
and projected local incidence angle (PLIA) normalized products, e.g., [
19,
34].
6.5. Common Framework with InSAR
ALOS-1 is considered an extreme case in terms of baseline variation for modern InSAR-capable SAR sensors. Since most modern day SAR sensors, e.g., TerraSAR-X, COSMO-SkyMed, ERS, EnviSAT etc, all exhibit much smaller baseline variation than ALOS-1, we can argue that the presented framework for flattening GTC products to RTC products using static flattening factors will work with all of these InSAR-capable sensors. We also observe that use of baseline information is a standard feature in InSAR time-series analysis, e.g., [
35], and the presented approach of transforming GTC products to RTC products can be easily incorporated into the same InSAR time-series analysis tools. We can further simplify the transformation process and reduce it to band-math operations in standard geospatial frameworks once baselines and incidence angles are represented as coarse three dimensional cubes [
36].
7. Conclusions
In this manuscript, we have described the design principles and implementation details of a method to transform GTC SAR products to RTC SAR products. We demonstrate that static flattening factors, one-per-imaging-geometry, is more than sufficient for efficiently flattening SAR imagery for terrain effects from missions characterized by a narrow orbit tube like Sentinel-1 on-the-fly. The presented approach is efficient, cost-effective, and highly scalable; and is suited for handling, in near-realtime, large volumes of SAR data that are expected to be acquired by missions such as Sentinel-1, NISAR, ALOS, and other commercial providers in the near future.
Author Contributions
Conceptualization, S.A.A., P.S.A., M.S.W. and M.T.C.; methodology, P.S.A. and M.S.W.; investigation, P.S.A.; data curation, P.S.A. and M.S.W.; software, P.S.A., M.S.W. and M.T.C.; validation, P.S.A., M.S.W. and S.A.A.; writing—original draft preparation, P.S.A., M.T.C. and M.S.W.; writing—review and editing, P.S.A., M.T.C., M.S.W. and S.S.A.; visualization, P.S.A. and M.T.C.; All authors have read and agreed to the published version of the manuscript.
Funding
This research received no external funding.
Data Availability Statement
The single-look complex (SLC) Sentinel-1 and ALOS-1 imagery and associated metadata are available at the Alaska Satellite Facility’s Vertex Portal here:
https://search.asf.alaska.edu/#/.
Acknowledgments
The authors would like to thank members of Descartes Labs Applied Science group for providing useful feedback during pipeline development. We thank John Truckenbrodt and Marcus Engdahl for providing useful feedback on our Sentinel-1 toolbox related experiments.
Conflicts of Interest
The authors assert that they have no conflicts of interest.
References
- Truckenbrodt, J.; Freemantle, T.; Williams, C.; Jones, T.; Small, D.; Dubois, C.; Thiel, C.; Rossi, C.; Syriou, A.; Giuliani, G. Towards Sentinel-1 SAR analysis-ready data: A best practices assessment on preparing backscatter data for the cube. Data 2019, 4, 93. [Google Scholar] [CrossRef]
- Torres, R.; Snoeij, P.; Geudtner, D.; Bibby, D.; Davidson, M.; Attema, E.; Potin, P.; Rommen, B.; Floury, N.; Brown, M.; Traver, I.N.; Deghaye, P.; Duesmann, B.; Rosich, B.; Miranda, N.; Bruno, C.; L’Abbate, M.; Croci, R.; Pietropaolo, A.; Huchler, M.; Rostan, F. GMES Sentinel-1 mission. Remote Sensing of Environment 2012, 120, 9–24. [Google Scholar] [CrossRef]
- Dyke, G.; Rosenqvist, A.; Killough, B.; Yuan, F. Intercomparison of Sentinel-1 Datasets from Google Earth Engine and the Sinergise Sentinel Hub Card4L Tool. 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS, 2021, pp. 1796–1799. [CrossRef]
- Small, D. Flattening Gamma: Radiometric Terrain Correction for SAR Imagery. IEEE Trans. Geosci. Remote. Sens. 2011, 49, 3081–3093. [Google Scholar] [CrossRef]
- Shiroma, G.H.; Lavalle, M.; Buckley, S.M. An area-based projection algorithm for SAR radiometric terrain correction and geocoding. IEEE Trans. Geosci. Remote. Sens. 2022, 60, 1–23. [Google Scholar] [CrossRef]
- Veci, L.; Lu, J.; Foumelis, M.; Engdahl, M. ESA’s Multi-mission Sentinel-1 Toolbox. EGU General Assembly Conference Abstracts, 2017, p. 19398.
- Meier, E.; Frei, U.; Nüesch, D. Precise terrain corrected geocoded images. SAR Geocoding: Data and Systems 1993, pp. 173–185.
- Zebker, H.A. User-Friendly InSAR Data Products: Fast and Simple Timeseries Processing. IEEE Geosci. Remote. Sens. Lett. 2017, 14, 2122–2126. [Google Scholar] [CrossRef]
- Small, D.; Schubert, A. Guide to ASAR geocoding. ESA-ESRIN Technical Note RSL-ASAR-GC-AD 2008, 1, 36. [Google Scholar]
- Agram, P.S.; Warren, M.S.; Calef, M.T.; Arko, S.A. An Efficient Global Scale Sentinel-1 Radar Backscatter and Interferometric Processing System. Remote. Sens. 2022, 14. [Google Scholar] [CrossRef]
- CEOS CARD4L Normalized Radar Backscatter Committee. Analysis Ready Data For Land: Normalized Radar Backscatter, Version 5.5. Technical Report CARD4L-PFS_NRB_v5.5, Committee on Earth Observation Satellites, 2021.
- Sentinel-1 Mission Performance Cluster. Sentinel-1 Burst Id map, version 20220530. Generated by Sentinel-1 SAR MPC.2022. Available online: https://sar-mpc.eu/test-data-sets/ (accessed on).
- Ulander, L. Radiometric slope correction of synthetic-aperture radar images. IEEE Trans. Geosci. Remote. Sens. 1996, 34, 1115–1122. [Google Scholar] [CrossRef]
- PROJ contributors. PROJ coordinate transformation software library. Open Source Geospatial Foundation, 2022. 2022. [CrossRef]
- Sandwell, D.; Mellors, R.; Tong, X.; Wei, M.; Wessel, P. Open radar interferometry
software for mapping surface Deformation. EOS Trans. Am. Geophys. Union 2011,
92, 234–234. Available online:
https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2011EO280002 (accessed on). [CrossRef]
- Rosen, P.A.; Gurrola, E.M.; Agram, P.; Cohen, J.; Lavalle, M.; Riel, B.V.; Fattahi, H.; Aivazis, M.A.; Simons, M.; Buckley, S.M. The InSAR Scientific Computing Environment 3.0: A Flexible Framework for NISAR Operational and User-Led Science Processing. IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium, 2018, pp. 4897–4900. [CrossRef]
- Samuele, D.P.; Filippo, S.; Orusa, T.; Enrico, B.M. Mapping SAR geometric distortions and their stability along time: A new tool in Google Earth Engine based on Sentinel-1 image time series. Int. J. Remote. Sens. 2021, 42, 9135–9154. [Google Scholar] [CrossRef]
- Lavalle, M.; Shiroma, G.H.; Agram, P.; Gurrola, E.; Sacco, G.F.; Rosen, P. Plant: Polarimetric-interferometric lab and analysis tools for ecosystem and land-cover science and applications. 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). IEEE, 2016, pp. 5354–5357. [CrossRef]
- Navacchi, C.; Cao, S.; Bauer-Marschallinger, B.; Snoeij, P.; Small, D.; Wagner, W. Utilising Sentinel-1’s orbital stability for efficient pre-processing of sigma nought backscatter. ISPRS J. Photogramm. Remote. Sens. 2022, 192, 130–141. [Google Scholar] [CrossRef]
- Rosen, P.; Hensley, S.; Joughin, I.; Li, F.; Madsen, S.; Rodriguez, E.; Goldstein, R. Synthetic aperture radar interferometry. Proc. IEEE 2000, 88, 333–382. [Google Scholar] [CrossRef]
- Kropatsch, W.G.; Strobl, D. The generation of SAR layover and shadow maps from digital elevation models. IEEE Trans. Geosci. Remote. Sens. 1990, 28, 98–107. [Google Scholar] [CrossRef]
- Simard, M.; Riel, B.V.; Denbina, M.; Hensley, S. Radiometric Correction of Airborne Radar Images Over Forested Terrain With Topography. IEEE Trans. Geosci. Remote. Sens. 2016, 54, 4488–4500. [Google Scholar] [CrossRef]
- Novellino, A.; Cigna, F.; Brahmi, M.; Sowter, A.; Bateson, L.; Marsh, S. Assessing the Feasibility of a National InSAR Ground Deformation Map of Great Britain with Sentinel-1. Geosciences 2017, 7. [Google Scholar] [CrossRef]
- Eckerstorfer, M.; Malnes, E.; Müller, K. A complete snow avalanche activity record from a Norwegian forecasting region using Sentinel-1 satellite-radar data. Cold Reg. Sci. Technol. 2017, 144, 39–51. [Google Scholar] [CrossRef]
- Karbou, F.; Veyssière, G.; Coleou, C.; Dufour, A.; Gouttevin, I.; Durand, P.; Gascoin, S.; Grizonnet, M. Monitoring Wet Snow Over an Alpine Region Using Sentinel-1 Observations. Remote. Sens. 2021, 13. [Google Scholar] [CrossRef]
- Burrows, K.; Marc, O.; Remy, D. Using Sentinel-1 radar amplitude time series to constrain the timings of individual landslides: a step towards understanding the controls on monsoon-triggered landsliding. Nat. Hazards Earth Syst. Sci. 2022, 22, 2637–2653. [Google Scholar] [CrossRef]
- Vollrath, A.; Mullissa, A.; Reiche, J. Angular-Based Radiometric Slope Correction for Sentinel-1 on Google Earth Engine. Remote. Sens. 2020, 12. [Google Scholar] [CrossRef]
- Durieux, A.M.; Rustowicz, R.; Sharma, N.; Schatz, J.; Calef, M.T.; Ren, C.X. Expanding SAR-based probabilistic deforestation detections using deep learning. Applications of Machine Learning 2021. International Society for Optics and Photonics, 2021; Volume 11843, p. 1184307. [CrossRef]
- Nagler, T.; Rott, H.; Ripper, E.; Bippus, G.; Hetzenecker, M. Advancements for Snowmelt Monitoring by Means of Sentinel-1 SAR. Remote. Sens. 2016, 8. [Google Scholar] [CrossRef]
- Zhao, J.; Pelich, R.; Hostache, R.; Matgen, P.; Cao, S.; Wagner, W.; Chini, M. Deriving exclusion maps from C-band SAR time-series in support of floodwater mapping. Remote. Sens. Environ. 2021, 265, 112668. [Google Scholar] [CrossRef]
- Hanssen, R.; Bamler, R. Evaluation of interpolation kernels for SAR interferometry. IEEE Trans. Geosci. Remote. Sens. 1999, 37, 318–321. [Google Scholar] [CrossRef]
- Durieux, A.M.; Ren, C.X.; Calef, M.T.; Chartrand, R.; Warren, M.S. Budd: Multi-modal bayesian updating deforestation detections. IGARSS 2020-2020 IEEE International Geoscience and Remote Sensing Symposium. IEEE, 2020, pp. 6638–6641. [CrossRef]
- Ticehurst, C.; Zhou, Z.S.; Lehmann, E.; Yuan, F.; Thankappan, M.; Rosenqvist, A.; Lewis, B.; Paget, M. Building a SAR-Enabled Data Cube Capability in Australia Using SAR Analysis Ready Data. Data 2019, 4. [Google Scholar] [CrossRef]
- Bauer-Marschallinger, B.; Cao, S.; Navacchi, C.; Freeman, V.; Reuß, F.; Geudtner, D.; Rommen, B.; Vega, F.C.; Snoeij, P.; Attema, E.; others. The normalised Sentinel-1 Global Backscatter Model, mapping Earth’s land surface with C-band microwaves. Sci. Data 2021, 8, 277. [CrossRef]
- Fattahi, H.; Amelung, F. DEM Error Correction in InSAR Time Series. IEEE Trans. Geosci. Remote. Sens. 2013, 51, 4249–4259. [Google Scholar] [CrossRef]
- Eineder, M. Efficient simulation of SAR interferograms of large areas and of rugged terrain. IEEE Trans. Geosci. Remote. Sens. 2003, 41, 1415–1427. [Google Scholar] [CrossRef]
|
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).