Retrieval of digital elevation models from Sentinel-1 radar data – open applications, techniques, and limitations

2.1.1 Theoretic framework

The basic idea behind synthetic aperture radar interferometry (InSAR) is the combination of information provided by two different radar images, which have a known distance, the perpendicular baseline (Figure 1a). Based on the positions of both satellites and their perpendicular baseline, three-dimensional information of the Earth’s surface can be retrieved by the different path lengths of the radar echo signals from a surface as received by both antennas [41,42]. In the case of SRTM, these images were acquired at the same time (t ₁ = t ₂) by two antennas, which were separated by a mast of 60 m length; this is also referred to as bistatic or single-pass interferometry. As the scattering mechanisms of the Earth’s surface rarely do change within the time of data acquisition and the perpendicular baseline is constant, the retrieved signals are largely coherent (in-phase, Figure 1b), and the acquisition geometry can be precisely resolved for the topographic variations in the observed area [41]. The same bistatic principle is applied in the TanDEM-X mission, which consists of a constellation of two similar but independent satellites, which simultaneously acquire images within a helix-shaped orbit [32]. During its operation, their distance was changed in various phases to systematically analyze the impact of the perpendicular baseline on the retrieved information [43].

However, most radar satellites operate in a monostatic mode, which means that they do not have a tandem configuration, but utilize images acquired by the same antenna, but at two different times (t ₁ < t ₂). This approach is known as repeat-pass interferometry and is limited by the revisit rate of a satellite, which is defined at the interval at which the sensor recurs to the same position after orbiting the Earth [24]. For the Sentinel-1 satellite mission, this period is 12 days for each of the two satellites, but as they have the same configuration and a systematic offset, the revisit time can be reduced to 6 days when the images of S1A and S1B are combined. Repeat-pass interferometry enables the measuring of elevations based on the same principles as single-pass InSAR but is unavoidably susceptible to a number of error sources originating from potentially different conditions at the time of the acquisitions. These include receiver phase instabilities, incorrect estimates of the perpendicular baseline, atmospheric differences, as well as changes in altitude or land cover between the first and the second image [44]. While atmospheric and sensor-specific factors introduce patterns that can be easily misinterpreted as topographic variation [45], the problem of phase decorrelation leads to bad interferometric data quality in general, mostly over vegetated areas, which cannot be restored [24,46]. Short wavelengths are particularly prone to temporal decorrelation, such as those transmitted and received by the C-band of Sentinel-1 (5.405 GHz, corresponds to around 5.5 cm). This will be exemplified in the following sections.

2.1.2 Main processing steps

In the following, the derivation of DEMs from InSAR is briefly explained. However, as the well-established development of approaches for topographic mapping dates back into the 1980s, only the most important processing steps are explained, and the reader is advised to consult the given references for more detailed descriptions and mathematical elaborations. In turn, the presented methods and their parameters are calculated based on the data listed in Table 1 and illustrated by the maps shown in Figure 4.

In accord with the slant-looking geometry of every SAR measurement, signals transmitted by the antenna return to the sensor at different temporal delays, based on the distance R of the objects on the Earth’s surface to the satellite [44]. As this signal propagates in a sinusoidal nature (Figure 2), the measured delay τ is equivalent to a phase change φ between the transmitted and the received signal [47]. As indicated in equation 1, φ is proportional to the two-way travel distance 2R divided by the transmitted wavelength λ.

Figure 2

Amplitude and phase of a radar wave (adapted from Reference [197]).

Due to differences regarding the spatial resolution of Sentinel-1 SLC products (5 × 20 m) and the scale of the wavelength (5.55 cm), the phase change between neighboring pixels, representing only the last fraction of 2R, looks random at first (Figure 3, left). However, when cross-multiplying the phase of the first image with the complex conjugate of the second image in a pixel-by-pixel manner, an interferogram is retrieved, which consists of phase variations φ as a result of topography φ _topo, the Earth’s curvature φ _flat, surface deformation between both acquisitions φ _disp, atmospheric contributions φ _atm, and system-induced noise φ _noise (equation 2) [48] under the given acquisition geometry (Figure 3, right; Figure 4c) [25].

Figure 3

Combination of two complex phase products to an interferogram.

Figure 4

Steps of DEM generation with Sentinel-1: (a) Sentinel-2 image from 02.07.2019 (for visual reference), (b) Sentinel-1 image from 02.07.2019, (c) interferogram from 26.06.19 and 02.07.2019, (d) coherence image, (e) unwrapped interferogram, and (f) hillshade of the derived DEM. The extent is shown in Figure A1. This figure shows the modified Copernicus data.

A mathematical description of this process, including the subtraction of the flat-earth phase φ _flat component (interferogram flattening), is given in previous studies [22,44,47]. However, the formation of an interferogram relies on the accurate stacking of the first (reference) and the second (secondary) images, so that the phase difference is calculated for the correct pixel pairs. This can be achieved at sub-pixel accuracy by methods of coregistration, which are based on cross-correlation and spectral diversity approaches [49,50,51]. This grants that the information of both images refers to the same object at the ground and that no errors are introduced at an early stage.

As presented by Hanssen [24], the main reasons for noise induced in interferograms are as follows:

• Changes of scattering mechanisms between both images (temporal decorrelation).

• Changes in the incidence angle between both acquisitions (baseline decorrelation).

• Penetration of the radar wave into mediums (volume decorrelation).

• Thermal noise and other antenna-related sources of bias.

• Additional processing-induced errors, such as low coregistration quality.

To estimate the quality of the interferogram throughout the image, the local coherence γ is calculated between the reference and the secondary image. It is a cross-correlation measure based on a defined window size as demonstrated in equation 3 and serves as a spatial indicator of the reliability and robustness of the later products derived from the interferogram [42,52]. In this equation, the amplitude of each pixel of the interferogram |v| is proportional to the product of the amplitudes of the reference |u ₁| and the secondary image |u ₂|, and E is the expected value of the random variable x [47,53].

The coherence ranges between 0 and 1, where large values indicate a high correlation between the reference and the secondary image, which is the case over urban areas and bare soils, while low values occur in areas of vegetation and other rapidly changing surfaces (Figure 4d). It helps to decide whether an area produces coherent (in-phase, Figure 1b) scattering. Although there is no universal threshold above which the phase information is considered usable, values between 0.1 and 0.4 are often used as a pragmatic solution to mask out noisy parts from the interferogram and to avoid errors becoming introduced during later steps [54,55]. However, it is important for the subsequent unwrapping of the interferogram that the remaining parts are still connected.

Although isolated patches of high coherence can be technically unwrapped correctly, the result for the full scene fails if they are not sufficiently connected over the entire image [24]. This is especially valid for the derivation of digital surfaces, which form a continuum within the scene [56,57].

As shown in Figure 4c, each color cycle in the interferogram represents one phase difference between the reference and the secondary image. These patterns are also called fringes, where each fringe has a distance of 2π radians with respect to the employed wavelength of the system. With respect to equation 2 and the examples presented in Table 1 (A1), the amount of topographic height variations φ _topo represented by this 2π cycle is dependent on the wavelength λ (5.55 cm), the perpendicular baseline b _perp (−153 m), the range distance between the sensor and the Earth R (845 km), and the incidence angle θ (39°) as demonstrated in equation 4a [24,27]. If this equation is solved for the given example, one fringe cycle H _2π represents 147.7 m of topographic difference (equation 4b), which is also known as the height ambiguity [48]. For the mapping of topography, smaller height ambiguities allow a more precise description of terrain features.

To convert this cyclic measure (also called the principal or wrapped phase) into a continuous measure, a wide range of phase unwrapping approaches have been proposed [58,59], and these integrate the phase at a pixel-by-pixel level over the entire image along closed paths, considering that phase differences greater than the principal interval of π cannot occur [47,60]. The result is a full phase representation of topographic variations as shown in Figure 4e. However, as mentioned earlier, this process is prone to errors induced by complex topographies, phase noise in areas of low coherence, as well as by atmospheric artifacts [24,48,56]. Smaller amounts of noise can be tackled by specialized phase filters [61,62,63,64,65], but especially in interferograms with larger patches of low coherence, errors can be introduced in the form of phase jumps or entire trends superimposing the actual topographically induced phase variation [22]. Furthermore, based on the relative nature of InSAR measurements and the unavailability of a priori information about the terrain, it is worth mentioning that no ultimate solution to phase unwrapping exists, and it is always an approximation [44,60].

In a final step, the unwrapped phase is geocoded and translated into a geophysical unit, which represents the height of the calculated surface over the spherical Earth [24]. This can be done by integrating the height of either surveyed reference points or extracted from an external DEM [66], but across-track approaches also exist, which are independent of a reference [67]. Finally, the retrieved heights are referenced to an ellipsoid (e.g., WGS84), projected into a coordinate reference system based on either geographical or UTM coordinates, and resampled to a uniform pixel spacing [22,66,68]. A potential result is shown as a hillshade representation in Figure 4f.

2.1.3 Impact of temporal baseline

As described in the previous sections, one major reason for phase decorrelation is the change of scattering mechanisms between the acquisition of the reference and the secondary image. Currently, Sentinel-1 is a mission employing two satellites, which allow the generation of repeat-pass interferograms with a temporal baseline of 6 days under best conditions. However, as land surfaces can rapidly change, this period is not ideal compared to the bistatic missions of SRTM and TanDEM-X. Given that the wavelength of Sentinel-1 is comparably short, vegetation areas are especially likely to lose their coherent scattering properties within a couple of seconds [54,69]. Therefore, it is impossible to retrieve reliable phase information on the height of forest stands or agricultural fields based on a single Sentinel-1 image pair [70]. These shortcomings were addressed in previous times when the temporal baseline between acquisitions of ERS-1 and ERS-2 was systematically reduced to 1 day [71] or images of ERS and ENVISAT of the same day were combined for DEM generation [72,73,74]. Yet, the limitations of volume decorrelation and increased sensitivity to tropospheric disturbance and water vapor persist for C-band systems [75,76].

To illustrate the impact of the temporal baseline on the quality of interferometric analyses, five interferograms were computed for one reference image (July 2, 2019) and secondary images acquired with a separation of 6 days (A1), 18 days (A2), 24 days (A3), 54 days (A4), and 96 days (A5) (Table 1). For reasons of consistency, all image pairs have a similar perpendicular baseline (discussed in Section 2.1.4).

As shown in Figure 5, a loss of coherence is observed with growing temporal baseline. On average, an extension of the temporal baseline from 6 to 18 days leads to a decrease of coherence of −19.2%. To put this observation into an environmental perspective, coherence was analyzed with respect to the main land cover classes in the area as retrieved from the Corine Land Cover 2018 dataset [77]. The largest relative decrease is observed for the nonforest vegetation areas (−30.6% from 6 to 18 days and −56.1 from 6 to 96 days), while urban areas are comparably robust at a high level (above 0.5). In turn, agricultural areas and water bodies already have low coherence values for very short temporal baselines and show a smaller relative decrease of only −20% over the entire period. The figure underlines that InSAR quality is highly sensitive to temporal decorrelation over natural surfaces, especially over areas of vegetation and water areas. While this relationship can be effectively exploited for multitemporal land cover mapping [78,79], the capabilities of Sentinel-1 for topographic mapping are severely limited in densely vegetated areas.

Figure 5

Loss of coherence of selected land cover classes for the investigated temporal baselines (shown Table 1).

To visually underline the findings of Figure 5, the coherence products of different baselines, as well as the interferograms resulting from these combinations, are compared in Figure 6. It shows the city of Erzincan, Turkey, which is surrounded by agricultural and natural plains and framed by the Esence Mountains in the north and the Mercan Mountains in the south. Height differences between the city and the ridges in the selected extent are up to 1,000 m. The figure illustrates that the fringes are well aligned with the topographic conditions of the mountains and that the interferogram forms a consistent pattern in the plains. However, with increasing temporal baseline, these patterns suffer from decorrelation, caused by phenological changes due to the crop growth in the agricultural areas and seasonal changes over natural vegetation. This decorrelation strongly affects the interferogram quality, indicated by a higher amount of phase noise and less pronounced fringes. Most importantly, the different parts of the interferogram are no longer connected, which will lead to severe errors in the subsequent unwrapping, which relies on the steady integration of the continuous phase signal [60]. The results of the 6-day baseline image pair are presented in Section 2.1.4.

Figure 6

Coherence (top) and interferograms (bottom) for selected temporal baselines (shown in Table 1). For reasons of visualization, the interferograms are combined with a hillshade representation retrieved from the SRTM data. The extent is shown in Figure A1. This figure shows the modified Copernicus data.

2.1.4 Impact of perpendicular baseline

From the perspective of the acquisition geometry, two points are crucial for the reliable retrieval of topographic information. First, the position of the satellites during the time of image acquisition must be known precisely. For Sentinel-1, this is granted by the provision of precise (or restituted) orbit information by the Copernicus Precise Orbit Determination (POD) Service [80]. These have an accuracy of 10 cm and are provided within a couple of days after image acquisition for non-time critical products [81]. Their high accuracy was also reported for combinations of S1A and S1B to facilitate the use of 6-day temporal baseline pairs [82]. Second, the perpendicular baseline (Figure 1a) of the repeat-pass acquisition should be between 150 and 400 m to allow a precise description of topographic variations by the observed fringes [59]. However, as reported in Section 1, the orbital tube of Sentinel-1 was designed for differential interferometry, which requires small perpendicular baselines [83]. For the majority of image pairs within intervals of 6 or 12 days, perpendicular baselines are below 150 m and even smaller than 25 m in most of all cases [84]. A detailed evaluation of perpendicular baseline lengths with respect to the different orbits of S1A and S1B is given in the Sentinel-1 performance report of 2019 [85]. Referring to equation 4, perpendicular baselines between 25 and 50 m generate fringes with height ambiguities of 500 and 1,000 m of altitude. Interferograms of such narrow perpendicular baselines are therefore less sensitive to the topographical contribution of the phase φ _topo and more suitable for the derivation of surface displacements between two passes [40]. The term of the critical baseline B _c was introduced by Zebker and Villasenor [46] as the upper limit of the perpendicular baseline above which the phase becomes completely decorrelated. It is dependent on the wavelength λ, the range distance R, the incidence angle θ, the inclination of the slope α, and the range resolution R _y and computed as shown in equation 5 (valid for flat terrains):

For Sentinel-1 (incidence angles between 29° and 46°), the critical baseline roughly ranges between approximately 5 and 10 km. Accordingly, it is not a limiting factor for topographic mapping under the acquisition geometries defined for S1A and S1B. This is also supported by the performance report, which denotes the maximum loss of coherence due to the length of the perpendicular baseline below 5% [85]. However, perpendicular baselines below 50 m are a problem for the retrieval of DEMs from Sentinel-1. It should be mentioned in this context that other C-band constellations employed significantly larger perpendicular baselines: ERS-1/ERS-2: 75–400 m [71,86,87,88,89]; ERS/ENVISAT: 1,650–2,050 m [72,90,91,92,93]; and Radarsat: 100–1,400 m [94,95], and also that the TerraSAR-X/TanDEM-X constellation allowed perpendicular baselines of up to 4 km [32,96].

The role of the perpendicular baseline for Sentinel-1 InSAR applications is illustrated in Figure 7. As also indicated by equation 4, the height difference is represented by one phase cycle increases with smaller perpendicular baselines. Accordingly, short baseline image pairs are less sensitive to elevation differences. The consequences on the quality of the derived DEM are displayed in the row at the bottom of Figure 7, which clearly indicates the loss of elevation information for smaller perpendicular baselines and the increasing amount of phase noise. This is also indicated by the drastically increasing RMSEs of the derived elevations, which were calculated in comparison with the available SRTM heights.

Figure 7

Interferograms (top) and hillshades of the derived DEMs (bottom) for selected perpendicular baselines (shown in Table 1). For reasons of visualization, the interferograms are combined with a hillshade representation retrieved from SRTM data. The extent is shown in Figure A1. This figure shows the modified Copernicus data.

2.1.5 Impact of other factors

Besides the temporal and perpendicular baselines, which are largely a matter of the selection of suitable images, there are geometric and topographic factors at the landscape scale, which also have an obvious impact on the quality of InSAR measurements. To demonstrate these effects, two profile lines were drawn across the study area as shown in Figure A1 (Supplementary data), one in azimuth direction across the valley (profile A, 40 km) and one in the range direction across the slopes of the Esence mountains in the north (profile B, 30 km). Elevations of the S1 DEM were extracted and analyzed for their differences to the reference SRTM heights. As shown by cross profile A (Figure 8, top left), there is a systematic overestimation of the overall elevation in the more level terrain (blue), while there is an underestimation of the heights of the mountain ridges (red). This is a typical sign for unwrapping errors accumulating across the image, slopes are especially underestimated then due to noise [97] and foreshortening [41]. A similar effect was reported by Singh et al. who created a DEM from ERS tandem pairs (temporal baseline of one day) for parts of the Swiss Alps and also observed that the underestimation of slopes led to a “smoothened” topography in comparison with a higher resolution reference DEM [98].

Figure 8

Impact of selected parameters on the height error (compared with SRTM). Top left: SRTM heights (red) vs S1 DEM heights (blue) along profile A (shown in Figure A1); top right to bottom left: height error in relation to coherence, local incidence angle, and slope for the pixels along profile B (shown in Figure A1).

Another reason can be orbital ramps, which are introduced by inaccurate coregistration and insufficient quality of orbit information [23,99]. One way to deal with these kinds of systematic errors is to model the overall spatial phase trend in the data by linear regression or polynomial approximation and remove it from the data before [86] or after unwrapping [43]. The impact of orbital errors on Sentinel-1 interferometry is comparably small because of the provision of precise orbit information [81] and high coregistration accuracy achieved through enhanced spectral diversity (ESD) coregistration, which corrects for the azimuth-dependent Doppler variation within the bursts [100]. Even if Sentinel-1 data can be processed with very high accuracies, already small azimuth offsets of a few centimeters can produce a phase ramp in the azimuth direction [101,102]. Together with the aforementioned slope underestimation, this explains most of the differences between SRTM and S1 elevations as observed along with profile A (Figure 8, top left).

But there are also effects in the range direction, which are discussed briefly in the following. It was already denoted in the study by Madsen et al. that height errors are strongly correlated with terrain modulation, especially with steep slopes [23]. Accordingly, the steepness of slopes is not only potentially underestimated over the entire image as mentioned earlier but also can lead to local height errors in extreme terrain. This effect is visible in Figure 8 (top right), where height errors are compared for all pixels along profile B in an aggregated form. The box plot shows that height errors, as the difference between the SRTM, are comparably small over slopes between 0 and 40°, but drastically increase for slopes above this value. The plot shows that underestimations of up to 150 m occur at all pixels along the profile, but with significantly higher interquartile ranges for slopes above 40°. Accordingly, in the investigated case study, the overall elevation estimates are good, but predominantly inaccurate for local extremes. This relationship between slope gradient and height errors was systematically analyzed and proven by Zhang et al. who compared a DEM derived from ENVISAT interferometry with the near-globally available SRTM [103].

To put these findings into another perspective, the second plot (Figure 8, bottom left) shows the height error in relation to the local incidence angle. This variable is different from the overall slope because it takes into consideration both the looking geometry of the sensor and the alignment of the landforms with respect to the flight direction of the satellite. As a consequence of the side-looking acquisition principle of spaceborne SAR sensors, geometric distortions occur in the image, often summarized as foreshortening, layover, and shadow [104]. Reasonably, these not only affect the retrieval of backscatter intensity but also result in errors in interferometric approaches, such as presented by Eineder and Holzner who demonstrate systematic effects for InSAR DEMs in areas of radar shadow [105]. In addition, Singh et al. and Breidenbach et al. report larger error rates for InSAR height estimates in areas facing away from the sensor [98,106]. This effect is also shown in the box plot, which confirms larger relative height errors for pixels with incidence angles below 30° and above 75°. The “optimal incidence angle for SAR images over alpine terrain” [105] of 45° shows the lowest error rates, furthermore underlining the local sensitivity of InSAR DEMs toward terrain aspects with respect to the satellite’s position. While the shadow areas are simply voids in the retrieved radar image, they cause phase jumps during the unwrapping, which, when untreated, cause false elevation estimates [22]. Different approaches are presented to overcome these voids in the phase image by finding optimum unwrapping paths [60], for example, by using coherence as a measure to connect the surrounding unwrapping paths [107], by restoring the phase in shadow areas using smooth phase gradients [108], or by the combination of interferograms of multiple looking directions and aspects [109,110].

Finally, sufficiently high coherence is required for accurate height estimates, which can strongly vary over the entire image. For instance, decorrelation can occur through variations in land cover and scattering mechanisms as discussed in Section 2.1.3, but also results from the aforementioned shadow or layover effects [97]. An extensive analysis of terrain effects on DEMs retrieved from ERS data is provided by Zebker et al. [111]. Unsurprisingly, InSAR height errors primarily occur in pixels with coherence smaller than 0.2 and continually decrease with higher coherence in our study as well (Figure 8, bottom right).

Another important source of error are atmospheric effects. Although the retrieved backscatter intensity of radar sensors is significantly less affected than the optical data, heterogenous refractivity distributions in the atmosphere can cause an inconsistent signal delay within an image [24]. This was already investigated by Goldstein in 1995 who identified patterns in an interferogram, which were not correlated with topographic features and could therefore be attributed to water vapor and turbulence in the troposphere [75]. However, their impact becomes visibly more obvious in DInSAR approaches where the topographic phase φ _topo is already removed [112]. Especially time-series approaches based on many interferograms allow to systematically screen for atmospheric effects that can then be subtracted from the analysis [47]. As removing the topographic phase is not possible within the scope of DEM generation, its contribution has to be assumed negligible [27] or targeted by more advanced methods, such as weather models, data fusion, or multi-baseline approaches [57,113,114,115].

Luckily, the unknown delay of tropospheric effects is independent of the wavelength, but inversely proportional to the length of the perpendicular baseline, so that it is minimized for configurations above 100 m as suggested in Section 2.1.4 [116]. However, the radar signal also undergoes a phase advance in the ionospheric layer of the atmosphere, leading to additional distortion, especially in the mid-latitudes [117]. Its effect on C-band Sentinel-1 data is lower than for long-wavelength systems, such as ALOS-2 (0.034 vs 0.13 ionospheric induced phase [118]), but especially for DInSAR approaches, it should not be neglected, where Gomba et al. propose a split spectrum–based approach for Sentinel-1 TOPS data [119].

3.2.1 Aims and data sources

As some studies were duplicate listings by both WoS and Google Scholar, the total number of papers relevant for this review was 21. A summary of all studies is given in Table 4. Most of these aim at the estimation of local or regional topography based on a pair of Sentinel-1 products, similar to the examples shown in Section 2. They are distributed over different continents, altitudes, and ecosystems: Northwestern Turkey [140,141], the Kavir Desert [142], and the cities of Tehran [143,144] and Jam [144] in Iran, the Cameron Highlands in Malaysia [145], the island of Mykonos, Greece [146], the Tonk district in India [147], and the municipality of Semarang, Indonesia [148]. To increase the quality of the elevation estimates, seven to eight images were used in studies on Ksiromero, Greece [149], Tsukuba, Japan [150], and the island of Ischia, Italy [151]. A study on the Greek islands of Mykonos and Lesbos takes a special role because 21 and 25 image products were combined to generate DEMs [152].

Three studies explicitly target the retrieval of surface heights of natural and man-made structures, which can be considered a special challenge because they are not part of the observed topography but rather have their own limited extent within the setting. These studies use different quantities of input images and deal with the estimation of canopy heights of the Tundi forest reserve, India [153] (n = 2), the surfaces of two glaciers (Videma and Frías) in Argentina [154] (n = 8 per glacier), and also building heights in the city of Larisa, Greece [155] (n = 6).

The third group of studies uses images of multiple dates to assess elevation changes related to geomorphologic processes. While this is usually done with differential interferometry (exclusion criterion for this review), three studies derived DSMs for pre and post landslide situations to calculate volume differences for the Hítardalur Valley, Iceland [156], and the landslides of Moira [157] and Panagopoula [158] in Greece. A fourth study investigated the movement of dunes around the city of Ardestan, Iran [159].

The study by Yang et al. [160] has a special role because of its experimental character. It does not aim at the construction of a DSM alone, but investigates if and how radar images of different pixel spacings can be combined to achieve a DSM at high spatial resolution. It was tested in the Kavir Desert in Iran based on o a 6 day image pair in a stripmap mode.

3.2.2 Methods and outcomes

This section presents the contents of the studies (Table 4) in greater detail and reflect their approaches with respect to the principles and aspects of InSAR elevation extraction as outlined in Section 2.1. This section summarizes the current efforts and shortcomings in DEM generation from Sentinel-1 data to stimulate joint discussions on how to shape this field in the future. To grant an objective view on these studies, emphasis will be placed on the impact of different parameters and the reported form of validation of results as major indicators for transparent and good scientific practice.

Articles that were published by the same authors are evaluated together to embrace them as parts of a larger research aim, which was split into different aspects.

Against this background, Nikolakopoulos and Kyriou were among the first to approach the potential of Sentinel-1 data for the mapping of topographic heights of the Greek islands of Mykonos [146] and Lesbos [152] in 2015, shortly after the launch of Sentinel-1 in late 2014. In both studies, their approach was based on a series of interferograms, which were systematically divided into two groups of small (nine image pairs [b _perp: 175–350 m]) and large (21 image pairs [b _perp: 175–1,225 m]) perpendicular baselines to analyze the effect of this parameter on the DEM quality for Mykonos. For the island of Lesbos, results gained from ascending (nine image pairs) and descending (14 image pairs) orbits were compared to analyze the impact of look direction on the final result. For both islands, image pairs of perpendicular baselines larger than 175 m (up to 1,225 m) are reported, so that all of them fulfill the quality criteria suggested by the scientific literature as outlined in Section 2.1.4. Neither the interferograms falls below the minimum baseline of 150 m nor exceeds the critical baseline of around 5 km. The actual method of combining elevations derived from the different pairs into a single result is not described in the studies, but it is referred to in a previous article by Choussiafis et al. [161], which proposes a form of mosaicking where elevation values at single pixels are combined based on weights derived from a height accuracy map, which is produced along each DEM. This height accuracy map is computed based on an external reference DEM and allows the determination of the quality of a local height estimation. After selection of a threshold value, which was reported to lie around the vertical accuracy of the reference DEM, height pixels that exceed this threshold value are masked out and only pixels remaining from other interferograms are then used to compute a final height value weighted by coherence and baseline. The authors report that this method produces DEMs with significantly higher similarity to the reference DEM compared to the standard averaging methods. The authors’ studies evaluated in this article are based on this method and reference it in their methods section; yet, short discussions on the principles underlying the mosaicking process would increase transparency and help readers understand how the DEMs were generated in each respective study. After computation of DEMs of the two islands, the studies propose various forms of validation: (a) visual comparison with reference DEMs, (b) statistical comparison with reference DEMs, and (c) the evaluation of the derived heights against measured ground control points (GCPs). Splitting the accuracy assessment into several parts is advisable because each of them addresses different aspects: visual comparison helps to evaluate the overall data quality and search for artifacts (e.g., generated during the unwrapping process) and systematic errors, such as underlying ramps or general trends (e.g., underestimation of high values, like demonstrated in Figure 8). This was done by comparing color-coded elevation maps and cross-profiles from the generated DEMs of Mykonos and Lesbos with reference DEMs. The authors used data of different spatial resolutions obtained from SRTM (30 m), ALOS PRISM (2.5 m), and the Hellenic Cadastre authority (5 m), and also DEMs derived from contour lines topographic maps (7.5 m). This allowed obtaining a more objective and scale-independent view of the results without being limited to only one reference dataset, although they conclude that in most cases, reference data of the highest resolution and the lowest vertical accuracy should be preferred. In a second step, descriptive statistics on the elevation values of the produced and reference DEMs are compared. This does not give detailed insights on the local accuracy of the results, but does enable checking if the DEMs fall within the same value ranges as the reference data (minimum and maximum) and for an overall offset of the values (median). For example, the S1 DEMs produced for both Mykonos and Lesbos have higher maximum values (392 and 970 m) compared to their reference datasets (360 and 940 m), indicating potential overestimation of the highest peaks. For a more robust analysis of these indications, GCPs of certified elevation were used for the computation of the RMSE as a more specific accuracy measure, which compares the difference of the produced DEM elevations to that measured at the GCPs. As shown in Table 4, depending on the reference DEM, the RMSE lies between 5 and 22 m for Mykonos and between 10 and 85 for Lesbos, but it has to be noted that the overall number of GCPs is not reported in the articles, but is considerably lower for Mykonos (∼19) than for Lesbos (>100), based on a visual inspection of the presented GCP maps. As a second measure, the average difference (AD) was calculated for Mykonos between the S1 InSAR DEM as well as both the Hellenic Cadastre reference DEM (−3.2 m) and the available GCPs (−19.5 m). With respect to their study design, the studies found no difference in accuracy between the short and long perpendicular baseline groups of Sentinel-1 image pairs because both groups lie within the critical baseline of around 5 km. In addition, DEM mosaics computed from ascending and descending tracks led to similar accuracies, so that no benefit regarding the looking direction could be identified for this area as well.

In a later study, the authors applied the same methodology of mosaicking of Sentinel-1 derived elevations to characterize natural karst depressions in limestones of the Ionian geotectonic zone in Western Greece [149]. Besides the elevation model generated using eight Sentinel-1 products, openly available data of SRTM, ASTER GDEM, ALOS GDEM (all 30 m), and a DEM were generated from contour lines of a topographic map (20 m), as well an orthomosaic generated by photogrammetric processing of aerial images by the Hellenic Cadastre (5 m). Of all elevation products, the Sentinel-1 DEM showed the largest RMSE 38.9 m based on the elevation information of 38 GCPs collected in the field, which is high compared with the openly available DEMs (8–22 m) and the Hellenic Cadastre (2.4 m). In addition, the suitability of the Sentinel-1 DEM for the automated detection of karst features was determined with an identification rate of 69%, which lies above those of ALOS (66%) and the Greek Cadastre (22%), but below those of SRTM (73%) or the one derived from the topographic map (89%). The authors assume that the automated detection of karst features is sensitive to the scale of the input DEM, and therefore, products of higher resolution partly result in the lower quality. Yet, the DEM from Sentinel-1 data resulted in similar results as other freely available elevation products.

The method of mosaicking several Sentinel-1 DEMs was also applied by this group of authors to an area in the northern Peloponnese, which is reported as being especially susceptible to landslides, and DSMs of different dates could be used for the estimation of surface changes [158]. This study does not analyze a special landslide event but can be seen as a preparatory work for the development of a monitoring strategy based on freely available InSAR data. Eleven Sentinel-1 products were used to generate a DEM as proposed by Choussiafis et al. [161], and outcomes were visually and statistically compared with a reference DEM of the Hellenic Cadastre (5 m), an unknown number of GCPs, and a DEM generated from a series of TerraSAR-X Spotlight products with a spatial resolution of 1 m. All three DEMs were statistically validated against the GCPs (but not against each other), reporting an AD of 9.1 m and an RMSE of 14.7 m for the Sentinel-1 DEM. The authors conclude that neither the Sentinel-1 nor TerraSAR-X DEMs provide the necessary accuracy to detect and measure small-scaled landslides in this area. This study was continued by Kyriou and Nikolakopoulos in 2018 [157] who undertook extensive surveys in the area with Global Navigation Satellite System (GNSS) measurements and flight campaigns of unmanned aerial vehicles (UAV) to measure topographic heights of the Moira landslide in January 2016 and had an extent of 300 m by 330 m. A total number of 30 interferograms were computed with perpendicular baselines between 97 and 172 m of which 15 describe the surface before the landslide event and 15 afterward. The study does not give information on how these 15 interferograms were merged to one final DSM. As the authors coincide with the previously presented studies, it is assumed that the same method of mosaicking was undertaken, but still, this is a drawback because it prevents reproducibility of the results. The elevation profiles of the pre-landslide and post-landslide DEMs were compared along three different transects, clearly showing a reduction of surface heights between 5 and 30 m, varying along the distance between the landslide crown at the top and the accumulation zone at the bottom. Unfortunately, thousands of GNSS measurements, which were collected along tracks in the study area after the landslide, were used only to estimate the volume of the displaced material by comparing them with the heights of the Hellenic Cadastre and not for the validation of the post-landslide DEM of Sentinel-1. Yet, the authors conclude from this assessment that the difference of surface heights retrieved by the two S1 DEMs are overestimated and do not provide a reliable estimation of the actual displacement.

The second group of articles by Ghannadi et al. who tested the interferometric capabilities of Sentinel-1 for different locations in Iran is presented later. Large parts of their study areas are characterized by sparse vegetation cover, so that the impact of temporal decorrelation is low compared to other studies. The first article targeted the derivation of surface topography around the city of Tehran based on a single image pair with a temporal baseline of 72 days [142,143]. No information is given on the perpendicular baseline, but the large distance of the fringes in the interferogram presented in this study indicates a baseline smaller than 100 m (comparable to that of the shortest perpendicular baseline shown in Figure 7). The results were validated against two sets of height reference points extracted from a DEM with a spatial resolution of 1 m (source not specified): the standard deviations (StDs) between both datasets were 1.26 m (flat areas, based on 138,761 reference points) and 10.32 m (mountainous areas, based on 78,196 reference points), with large proportions of errors smaller than 10 m, but outliers reaching up to 20 m.

To tackle these outliers more effectively, the findings of Ghannadi et al. were refined in another study in Iran [144]. It uses two image pairs of different regions around the cities of Tehran and Jam with temporal baselines of 180 and 24 days, respectively. It is not justified in the article why these image pairs were selected, yet perpendicular baselines are not presented. The approach of the study is of rather experimental nature and analyzed how DEM accuracies change after the application of a median filter and a two-dimensional Kalman filter [162]. They also used a reference DEM generated by photogrammetric processing of GeoEye-2 image products with a spatial resolution of 1 m and a vertical accuracy of 0.5 m. Presumably, the same reference DEM was used in the studies mentioned earlier as well although its accuracy is not evaluated at a greater detail. Findings of both regions show that the RMSE of 24.3 m (Tehran) and 20.7 m (Jam) can be reduced slightly by the median filters (23.1 m and 20.4 m, respectively) and even more so by the proposed 2D Kalman filter (20.1 m and 18.5 m, respectively). The same is reported for the mean absolute error (MAE), which also slightly decreases after postprocessing. However, the authors conclude that the effect of filtering on the DEM quality is mostly limited to outliers, and it can be assumed that image pairs with more suitable temporal baselines (6–12 days) and larger perpendicular baselines can provide DEMs with lower error rates in general.

A similar study design has been presented by Atalay and Sefercik who estimate topography based on a single pair of Sentinel-1 images with a temporal baseline of 12 days in Zonguldak in Turkey [140] and validate the results based on a reference DSM with a spatial resolution of 10 m provided by local authorities. They performed a co-registration of the elevation models to minimize systematic errors between both datasets to correct for horizontal shifts. This step resulted in horizontal shifts of 39.4 m and −23.8 m in x and y directions, which were corrected before the accuracy assessment. Afterward, they compared the StD and the normalized median absolute deviation (NMAD) between the reference DEM and the S1 DEM. The calculated StD was 4.2 m on the flat terrain and 5.2 m for the entire area. To put these values into perspective, StD was also assessed between the reference DEM and the AW3D30 model, attributing similar error rates (StDs of 4.0 and 5.5 m). The calculated NMAD was 3.8 m (flat) and 4.0 m (entire area) for Sentinel-1, and 3.2 m (flat) and 4.7 m (entire area) for AW3D30. Finally, the morphologic quality of the DSMs was assessed by visual comparison of color-coded contour levels. The authors conclude that Sentinel-1 data can deliver DEMs of similar quality compared to AW3D30, but the latter is more precise on flat terrain in the study area.

The comparison of DSMs from SRTM, AW3D30, and a Sentinel-1 image pair (b _temp = 12 days, b _perp = 86 m) was extended by Sefercik et al. for the metropolitan area of Istanbul [141]. A reference DSM with a spatial resolution of 5 m was used to calculate StD and NMAD of the height difference for different surface types. The study showed that StD of Sentinel-1 (5.6 m) is lower than that of SRTM (6.2 m) and AW3D30 (6.3 m) over the entire area (water bodies excluded), but especially AW3D30 performs better on flat terrain (4.9 m). With respect to different types of land cover, open areas were found to have the lowest NMAD (3.5 m), followed by forest (4.3 m) and built-up areas (3.8 m) for the Sentinel-1 DSM. Additional height error maps, contour lines, and cross sections were presented to illustrate the differences between the datasets. Especially the profile cross sections showed that some heights of building blocks were visible in the AW3D30 because it derived from a product with initially 5 m spatial resolution, while SRTM and Sentinel-1 were of comparably poor quality over the urban area.

Ahmadabadi et al. used a Sentinel-1 image pair separated by 12 days and a perpendicular baseline of 77 m to create a DEM of dunes near the city of Ardestan in the Iranian desert. It was validated against 110 GCPs with height information and achieved an overall agreement of R ² = 0.73, an RMSE of 2.42 m, and an MAE of 2.1 m. They used it to estimate the direction and the horizontal movement of sand dunes based on a reference DEM and additional Sentinel-2 data, which were in line with the wind directions of a nearby meteorologic station. The detected horizontal displacements of the dunes showed an RMSE of 4.7 m and were found to be more frequent with increasing elevation, although this could also be attributed to a computational bias. The authors conclude that the looking direction of the satellites should be taken into consideration for future studies to increase accuracy.

In contrast to the sparsely vegetated areas, which were presented in the previous studies, Mohammadi et al. present a study based on an image pair of the Cameroon Highlands in Malaysia. which is characterized by dense forests, farmlands. and orchards and altitude differences of more than 1,000 m [145]. These challenging conditions, combined with a very large temporal baseline of 602 days, led to an entirely decorrelated interferogram without any fringe patterns as shown in the figures of the article. The unwrapped phase is therefore characterized by a strong ramp in the azimuth direction superimposing any topographically induced phase variations. As indicated by the graphical workflow, the topographic phase removal was applied to the data, and hence, the topographic phase estimation based on the SRTM data was then wrongfully interpreted as the actual output DEM. Accordingly, the accuracy assessment led to a coefficient of determination R ² of 0.996 and a standard error of 1.9 m, which was basically just the SRTM data validated against itself at a resampled spatial resolution of 10 m. The authors conclude that the method was highly effective, but based on the presented data sources and the outlined workflow, the results clearly underlie false conclusions and were based on technically incorrect processing.

A short study was presented by Soni et al. who used a 12-day image pair to derive elevations in the district of Tonk, India [147]. No information on the perpendicular baseline is given, and the quality of the resulting DEM was only visually validated and discussed based on the coherence layer, which was reported low for croplands and water bodies. The authors suggest addressing this issue by combining images of S1A and S1B to reduce the temporal baseline to 6 days.

While the previously evaluated studies did not present or discuss the role of perpendicular baseline, all the following studies placed more emphasis on the selection of suitable image pairs with favorable configurations. Sunu et al. found an image pair separated by 18 days but with a perpendicular baseline of 146 m, which comes close to the ideal value range suggested in the previous studies mentioned in Section 2.1.4. Their study area is the region around Semarang, the capital of Central Java province in Indonesia. Altitudes roughly range between 0 and 500 m, and large parts of the study area are covered by swampy forests and cultivations. Therefore, the fringes in the interferogram are visible, but are interrupted by low coherence areas at several points. Their phase information could not be unwrapped correctly leading to artifacts in the final DEM. Validated against 137 GCPs with elevation information, the AD was 38.7 m and reached values up to 206 m. Also, the RMSE of 52.3 m is comparably high, and according to the authors, the largest errors are concentrated in swamp areas and dense vegetation, while the AD was lowest over urban areas.

Similar challenges were reported by Ullo et al. who generated a DEM of the island of Ischia, Italy, in the context of an earthquake event, which was also assessed by differential interferometry (DInSAR). Hence, they used four different image pairs of varying look directions (ascending and descending), temporal baselines (6–12 days), and perpendicular baselines (between 7 and 108 m). They carefully discuss the potential impact of different error sources (atmospheric noise, temporal decorrelation, look angle, and surface orientation) and also underlining that pairs with the largest perpendicular baseline and highest coherence should bring the most promising results for topographic mapping. The resulting DEM contained smaller errors in the areas of dense forests, but according to the authors, the overall elevation patterns and value ranges were in accordance with those official topographic maps. No statistical evaluation based on independent reference data has been carried out for the generated DEM. Similar to the study by Mohammadi et al. [145], the information content of the topographic phase, which was estimated for the estimation of displacement of the study, was partly misinterpreted within the scope of DEM generation.

To study the impacts of the perpendicular baseline on DEM generation, Nonaka et al. analyzed six interferograms of Sentinel-1 image pairs of the city of Tsukuba, Japan [150]. Perpendicular baselines ranged between 46 and 137 m with temporal baselines between 6 and 24 days. They derived the height ambiguities for these pairs according to equation 4a and confirmed that larger height ambiguities caused by shorter perpendicular baselines lead to lower DEM accuracies. These were assessed by using a quantity of 61 height reference points and were reported the lowest (RMSE of 6.3 to 7.7 m) for the pairs with perpendicular baselines larger than 100 m and reached 31 m for the pair with 72 m and 24 days. The study furthermore showed that errors were larger for reference points collected over parking lots and roads than those from other ground surfaces. As the temporal baselines were not constant through the image pairs, the interpretability of the results is limited.

While the previous study intentionally limited the GCPs to ground surfaces, Letsios et al. tested the capabilities of Sentinel-1 interferometry for building height estimations in Larisa, Greece [155]. Six image pairs with temporal baselines between 6 and 28 days and perpendicular baselines between 114–171 m were used to generate individual DSMs, which were stacked and averaged in a second step to form a final urban DSM. Validation was based on 35 elevation measures collected at rooftops to obtain the topographic height including the surface of the building objects. Absolute errors ranged between 0.01 and 4.2 m and resulted in a mean average error (MAE) of 1.69 m. The authors conclude that the results can be useful at the block level, but the spatial resolution of Sentinel-1 is not sufficient to reliably predict single building heights. Also, the difference between the ground height and the object height is not assessed in this study.

Kumar and Krishna also investigated the interferometric potential for the retrieval of canopy heights of the Tundi forest reserve (India) comparing image pairs from Radarsat-1 and Sentinel-1. Their comparability was limited because of the different polarization (HH and VV), direction (ascending and descending), temporal baseline (24 and 12 days), and perpendicular baseline (491 and 90 m). After interferometric retrieval of surface heights, including the forest canopy, ground heights have been assigned to each forest area based on the nearest available vegetation gap to derive the local height of each forest stand. Different accuracy measures were calculated based on 11 height measurements collected in the field, each consisting of an average of 25 trees: R ² between measured and predicted tree height was found to be 0.89 for both sensors, and RMSE (1.3 m) and MAE (1.34) were slightly lower for the Sentinel-1 canopy model, despite visually reported noise within dense forest areas caused by phase decorrelation. These low measures must be interpreted with care because of the low number of reference points (n = 11) and the long time between the Radarsat-1 DSM from 2004 and the Sentinel-1 DSM from 2017 which includes an expected height increase of around 4 m according to the authors. This is underlined by a comparably low R ² between both DSMs of only 0.48. Still, the use of forest gaps to determine the difference between ground and surface heights was proven to be effective in this study and expected to be more efficient when more image pairs were combined.

A study on glacier heights in contrast to terrain heights was proposed by Ortone Lois for two glaciers in the Los Glaciares National Park in Argentina, with altitude differences of over 2,000 m [154]. Four Sentinel-1 image pairs of different years were used for each of them, with temporal baselines between 6 and 24 days and perpendicular baselines between 24 and 180 m. No statistical accuracy assessment was conducted, but the impact of rainfall events on the phase information as well as the loss of coherence with the increasing temporal baseline was presented in the study. The quality of results was furthermore assessed by visual inspection of hillshade representations of Sentinel-1 DSMs compared to the SRTM data. The author concludes that the impact of vegetation, steep slopes, and moving ice masses was severely degrading the phase quality, and pairs with perpendicular baselines below 70 m were not usable in this case.

An approach comparing DSMs of three different acquisition periods was pursued by Dabiri et al. who used 6 day image pairs at perpendicular baselines between 110 and 159 m to assess the volume changes of a landslide in the Hítardalur Valley, Iceland [156]. As this event exceeded the maximum detectable displacement (half the size of the wavelength) and included non-coherent surface changes, a DInSAR approach was no longer applicable [24]. After masking the landslide area, volumes of 7, 12, and 109 million m³ were estimated by the three different post-event DSMs. To put these large variations into perspective, the ArcticDEM (2 m) [163] was used to validate the DEM areas outside the landslide event. Besides a statistical comparison, the vertical accuracy of the three post-landslide DSMs was assessed by the RMSE, while their horizontal accuracy was tested for autocorrelation by Moran’s I index. It was applied to the difference image between the S1 DSM and the ArcticDEM and tests for local homogeneity within a kernel of eight neighboring pixels [164]. Results show that the DEM calculated from the ascending data with a perpendicular baseline of 134 m has the lowest RMSE (30 m compared to 37 and 41 m), but the DEM with a perpendicular baseline of 142 m produced the best volume estimation (7 million cubic meters, also confirmed by independent literature [165]), and has the highest Moran's I (0.97). The authors conclude that besides the factor of the perpendicular baseline, which has been carefully selected in the study, temporal decorrelation and weather conditions have a significant impact on the quality of the results, and therefore, the suitability of Sentinel-1 for landslide volume estimation has to be further investigated.

The last study presented in this work by Yang et al. is of an experimental nature and tested the hypothesis regarding if Sentinel-1 products could be combined with images of higher resolution to generate interferograms with higher details for DEM generation [160]. Their test was based on an image pair acquired in a stripmap mode (5 m spatial resolution) with a temporal baseline of 6 days acquired over the Kavir Desert in Iran, of which one was downsampled to 10 and 25 m for testing reasons. In a second step, an interferogram was computed using a compressive sensing approach developed by the author based on the product of high resolution (5 m) and each of the lower resolution images (10 and 25 m). The filtered and unwrapped phases were then compared to the high-resolution outcome by the relative RMSE, resulting in a relative phase difference of −28.7 (10 m) and −25.1 (25 m). No conversion to metric elevations was carried out, but the authors conclude that this method is suitable to increase the spatial resolution and degree of detail of interferograms. Yet, it remains unclear why the relative RMSE used for the quality assessment is negative and given in the decibel unit. Accordingly, generating interferograms from both stripmap and TOPS products should be further investigated.

Figure 11 presents a short summary of the literature review: around one half of the 21 analyzed studies dealt with the analysis of a single image pair, while the others used more advanced approaches. Most of the studies underlined the importance of short temporal baselines, but only half of them explicitly worked with pairs of suitable perpendicular baseline. It is unclear if the others simply did not mention this parameter or if they were not aware of its impact. Sixty-eight percent of all studies used SNAP for the data processing, which is not surprising, because it is the main software provided by ESA when making use of their satellite constellations. Besides this, ERDAS and Matlab were used as commercial alternatives. Half of all studies used field measurements of heights to validate their DSMs, which is considered the most precise validation, but collecting such data is time consuming, and many studies only used less than 50 GCPs for the calculation of specific error measures. Accordingly, these points were not necessarily representative of the entire study area. In cases where the field data were not available, validation has been conducted based on a reference DEM of higher resolution. This approach allows a considerably larger number of values for statistical comparison, but the accuracy metrics can be biased by the quality of the reference DEM. An advantage of this approach is that it allows the identification of the spatial distribution of errors, for instance, in relation to different types of land cover or with respect to topographic characteristics (e.g., larger errors at slopes). Based on either point-wise measurements or a reference DEM, most studies were able to conduct a statistical accuracy assessment. RMSE and MAE were the most popular measures because they are resolved in the same unit as the investigated DEM (m), but also the coefficient of determination R ² was frequently used to quantify the overall agreement. Only few studies distinctively separated horizontal and vertical accuracy or checked for spatial autocorrelation with Moran’s I. Finally, the main reported reason for low quality was temporal decorrelation, causing noise in interferograms because of low coherence. But it is obvious that some studies were not fully exploiting the potential range of the imaging geometry by selecting pairs with higher perpendicular baselines.

Figure 11

Overview of the study area. Extents and profiles used for demonstration in Section 2.1. This figure shows the modified Copernicus data.