Investigating the Impact of Digital Elevation Models on Sentinel-1 Backscatter and Coherence Observations

Abstract

Spaceborne remote sensing can track ecosystems changes thanks to continuous and systematic coverage at short revisit intervals. Active remote sensing from synthetic aperture radar (SAR) sensors allows day and night imaging as they are not affected by cloud cover and solar illumination and can capture unique information about its targets. However, SAR observations are affected by the coupled effect of viewing geometry and terrain topography. The study aims to assess the impact of global digital elevation models (DEMs) on the normalization of Sentinel-1 backscattered intensity and interferometric coherence. For each DEM, we analyzed the difference between orbit tracks, the difference with results obtained with a high-resolution local DEM, and the impact on land cover classification. Tests were carried out at two sites located in mountainous regions in Romania and Spain using the SRTM (Shuttle Radar Topography Mission, 30 m), AW3D (ALOS (Advanced Land Observation Satellite) World 3D, 30 m), TanDEM-X (12.5, 30, 90 m), and Spain national ALS (aerial laser scanning) based DEM (5 m resolution). The TanDEM-X DEM was the global DEM most suitable for topographic normalization, since it provided the smallest differences between orbital tracks, up to 3.5 dB smaller than with other DEMs for peak landform, and 1.4–1.9 dB for pit and valley landforms.

1. Introduction

Synthetic aperture radar (SAR) is an active imaging system with several advantages over optic sensors, such as Landsat OLI (Operational Land Imager) or Sentinel-2 MSI (Multi-Spectral Imager). SARs are independent of solar illumination and use wavelengths that can penetrate cloud cover and have unique interactions with ground targets. Furthermore, the capability to transmit and receive signals enables the use of both phase and polarization information, to monitor, among others, landslides, avalanches, snowmelt, and forests [1]. Terrain orientation affects the intensity of the backscattered signal based on lambert cosine law. The signal is further affected by the terrain scattering area. The SAR technique uses the return time to convert a table of recorded echoes into an image (focusing). These times are shortened in sensor-facing steep slopes, causing echoes to overlap and slopes to appear shortened in the focused image (a pixel represents more area). Both effects can be compensated for by using the acquisition geometry parameters and a digital elevation model (DEM). We refer to this compensation as topographic normalization. The backscattered intensity is normalized by accounting for the scattering area and the local incidence angle [2,3]. The coherence normalization is based on removing the topographic phase component from the interferogram before estimation [4]. Normalization results are heavily dependent on the DEM characteristics and quality [2,3] as they can be generated using different data sources, including remote sensing (optic, SAR, airborne laser scanning—ALS), and processing techniques (ALS point cloud, stereography, interferometry, and radargrammetry)

ALS uses the delay between emission and reception of light pulses to determine the 3D position of objects. When the pulses return, the instrument registers an intensity profile, which can be completely (full-waveform) or partially recorded (discrete return, i.e., the position of the leading edge before the peak) [5]. Airborne discrete return systems have become the source of national elevation datasets for many countries [6,7]. Spaceborne full-waveform Light Detection and Ranging (space LiDAR, such as the ICESat (Ice, Cloud, and land Elevation Satellite)) data have also been employed in the context of topographic mapping as a primary source for calibration or validation of global elevation datasets derived from other sensors [8,9]. However, global topographic mapping from space based on LiDAR is difficult, due to the sparse sampling, and the sensitivity to cloud cover.

Stereoscopic techniques are based on differences in the line of sight to objects (parallax) for common points (tie-points) in an overlapping set of images. Results are dependent on tie-point quantity, image contrast, noise, and features (such as shadows and homogeneous surfaces), which may pose problems [10,11,12]. These techniques have been applied over aerial [6,11] and satellite optical images [9,13], as well as SAR imagery (i.e., radargrammetry) [10,14]. Stereoscopic processing of optical images was used to generate the ASTER global DEM (ASTER GDEM) and the ALOS (Advanced Land Observation Satellite) World 3D Digital Surface Model (ALOS AW3D DSM). The ASTER GDEM was created using imagery from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) onboard the Terra satellite. ASTER stereo pairs were formed from two near-infrared images (nadir, backward) with 15m resolution. The AW3D DEM was based on data from the Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) onboard the ALOS. PRISM stereo acquisitions were formed with three panchromatic images (forward, nadir, backward) with 2.5 m resolution. In both cases, the most challenging task was masking clouds, snow, ice, or water on every acquisition [9] to avoid the introduction of outliers. Afterwards, images from each individual acquisition were matched, and elevation was calculated. All height estimates from individual acquisitions were stacked to ensure continuity and reduce noise [9,13]. The AW3D was corrected for biases using preexisting data, such as ICESat shots and the preexisting Shuttle Radar Topography Mission (SRTM) DEM (see next paragraphs) [9].

Even though the processing chains for ASTER GDEM and AW3D were similar, their accuracies are different. Studies comparing global DEMs based on ground control points (GCP) report that the ASTER GDEM has larger uncertainties and is affected by striping, hummock-like artifacts, and outliers [12,15,16,17]. These artifacts may stem from the tie point generation, the choice of band (NIR), its relatively low spatial resolution [16,17], or unremoved cloud patterns [15]. AW3D performed better, although hillslope and step-like artifacts (scene mismatch) have been found [12,16,17,18,19].

SAR interferometry uses two co-registered SAR images acquired from close orbits. The interferogram (i.e., the phase difference between the two SAR images) relates to the 3D position of each target on the ground. Thereof, an interferogram reproduces the topographic information, which appears in the form of fringes as phase is measured between 0 and 2π. To obtain absolute phase values, from which elevation can be estimated, the interferogram is unwrapped. Unwrapping may be hindered in areas of steep topography or areas affected by the lack of coherence between images in consequence of changes between acquisitions (wind-induced motion, precipitation, etc.) [1,14]. SAR interferometry was used to generate two global DEMs, the SRTM DEM, and the TanDEM-X DEM.

The SRTM acquired data over 80% of the Earth’s land surface (60° N–56° S) on an 11-day orbital flight in February 2000. SRTM operated two antennas physically separated in space by 60 m at C-band, as well as X-band. At C-band, a gap-free coverage was obtained with single-pass interferograms. The interferometric height was reconstructed (unwrapping) and re-gridded into map coordinates with variable-resolution smoothing. Data takes were combined using coincident tie points [20]. The main artifacts of the SRTM dataset were related to striping from uncompensated movements of the mast, voids in correspondence of steep slopes or for low coherence areas, or coarser than nominal spatial detail from the re-gridding step [17,20,21,22]. The global TanDEM-X DEM was generated using SAR data acquired during 2010–2015 by the TanDEM-X and TerraSAR-X satellites flying in formation. Individual scenes were focused, multi-looked (sample averaging) to 10–12 m pixel spacing and unwrapped. The DEM was generated in an iterative process with the first global coverage using data acquisition parameters (baseline) adequate for moderate terrain. The second global coverage was shifted half a swath, and its unwrapping was aided by the data from the first coverage. Over some areas, further coverages were acquired from a different viewing geometry to avoid errors caused by topographic distortions or volumetric scattering (e.g., forest, desert) [23,24]. Areas with height ambiguities were infilled using radargrammetric processing of the scenes. TanDEM-X DEM data are distributed with 12.5 m (0.4 arcsec, original), 30 m (1 arcsec), and 90 m resolution (3 arcsec). The latter was generated by the unweighted average of the overlapping 12.5 m pixels [17].

All DEMs are affected by contributions from elements covering the terrain, such as cities or vegetation, thus reporting elevations higher than those recorded for the ground control points [15,18]. For this reason, elevation refers to surface elevation rather than terrain elevation. For the specific case of vegetation, the main reason for this is the different penetration of each wavelength. On the one hand, the nanometric-scale wavelengths employed for generating photogrammetric DEMs have limited penetration, and thus the tie points and the generated surface tend to reflect canopy surface height. On the other hand, the centimetric wavelengths employed by SAR sensors are able to penetrate further, albeit the scattering center height depends on the frequency employed and the vegetation structure [20,25,26,27]. Furthermore, the quality of InSAR-based DEM also depends on the spatial and temporal baselines, or unaccounted sensor movements.

Despite the rather large range of DEMs available, their effect on terrain normalization is poorly understood. Hoekman [28] suggested that the SRTM 90 m spacing may not be adequate for the radiometric normalization of the SAR backscatter in complex terrains. Recently Truckenbrodt [29] compared and tested the SRTM (30 and 90 m), AW3D (30 m), and TanDEM-X (90 m) DEMs in the context of radiometric terrain normalization of Sentinel-1 data. They analyzed the deviation of each DEM from the pixel-wise median of all DEMs and performed a regression analysis between terrain flattened γ⁰ and the local incidence. The deviation analysis showed that the SRTM DEMs have the smallest difference from the median values, but the 30 m SRTM version had high deviation artifacts at one test site. The same errors were found for SRTM and AW3D, as the latter has been infilled with data from the former due persistent cloud cover. Both AW3D and TanDEM-X DEM contained outliers or noise over water areas. The 90 m TanDEM-X DEM was found to contain several large artifacts in mountainous areas. Regression analysis showed that all the analyzed DEMs largely removed the terrain influence, with complete removal (i.e., slope of 0) being observed in some experiments using higher resolution DEMs (SRTM 1 arcsec, AW3D) [29].

The objective of this study was to investigate the impact of global DEMs on the normalization of SAR backscatter and coherence observations by Sentinel-1 at two sites characterized by complex topography. To build on previous literature [29], we analyzed the performance of 12.5 (resampled to 20 m, see Section 3), 30, and 90 m pixel Tandem-X DEMs, along with AW3D (30 m), SRTM (30 m). A very high-resolution (5 m) ALS DEM was used to benchmark results. The impact of terrain normalization was assessed by investigating the inter-orbit variability of the observations by land cover and landforms. Then, we evaluated a land cover classification scheme based on the observations normalized for topography.

2. Study Area and Satellite Data

The study area consisted of an N-S transect over the Romanian Carpathians (11,700 km²) and the National Park of Sierra Nevada in Spain (2360 km²) (Figure 1). We could not add additional sites, due to TanDEM-X scientific proposal area limitations [30]. The sites were selected to account for the different vegetation types and structures encountered in the temperate and Mediterranean climates. Due to the more humid climatic conditions, the vegetation in the Carpathians is characterized by denser, taller, and more diverse forest types (broadleaf, needleleaf, mixed) when compared to the sparser and shorter forests dominated by pine species encountered in the Sierra Nevada.

Histograms of natural land covers for each site have been plotted to show their frequencies relative to the slopes they occupy (Figure 2). Needleleaf and mixed forests occupy steep slopes, whereas grasslands and broadleaf forests occupy moderate slopes at the Romanian site. At the Spanish site, bare and needleleaf forests occupy moderate slopes, although for the former, a significant fraction of the pixels occupy steep or very steep slopes.

For each site, we assembled the SRTM 1-arcsecond, i.e., 30 m, DEM [31] from the United States Geological Service (USGS) Earth Explorer [32], the AW3D 30 m DEM from the Japan Aerospace Exploration Agency (JAXA) Earth Observation Research Center (EORC) [33] and the TanDEM-X DEM (©DLR 2019) with a pixel spacing of 12.5 m (original resolution), 30 m and 90 m (resampled) from the German Aerospace Agency [34]. In addition, for the Spanish site, we used an ALS-based DEM to benchmark the results obtained from three global DEMs. The ALS DEM was available through the Spanish national plan of orthophotography (PNOA) from the National Center of Geographic Information of Spain (CNIG) [35,36]. The ALS DEM was created from ALS point clouds with a density of 0.5 returns/m². The ALS scan over our study area was performed in 2014 with the LEICA ALS60 sensor. The points were translated from ellipsoidal to ortho-metric heights, assigned color (RGB and near-infrared) from PNOA orthophotographs, and classified automatically using TerraScan [37,38]. Classification eliminates returns considered noise and filters point to avoid oversampling due to flight strip overlap. Ground points classification is based on slope, rugosity, and return count. Vegetation and Buildings are classified based on height, separating both based on NDVI values. The DEM is generated by calculating the mean value of all ground returns within a 5m pixel [39]. The reported accuracies of the DEM products are presented in

Table 1. Reported accuracies of the DEM used in this study.

DEM Product	Pixel Spacing	Accuracy		Relative Vertical Accuracy	Coverage	Reference
		Horizontal	Vertical
SRTM DEM	~30 m	≤12.6 m	≤9 m	≤9.8 m	Nearly global (60° N–56° S)	[40]
AW3D DSM	~30 m	-	<7 m	>3 m (slope ≤ 20%) >5m (slope > 20%)	Global	[41]
TanDEM-X DEM	~12.5 m	<10 m	<10 m	2 m (slope ≤ 20%) 4 m (slope > 20%)	Global	[23]
PNOA LiDAR DEM	~5 m	≤0.5m	≤0.5 m	-	Spain	[42]

.

The SAR dataset consisted of a time series of Sentinel-1 dual-polarized (VV and VH) images acquired in the Interferometric Wide Swath (IWS) mode. The images were obtained in Single Look Complex format (SLC) with a pixel spacing of 14.1 m in azimuth and 2.3 m in range. All SAR images were resampled to a pixel size matching the different DEMs used (see Section 3.2). For the Romanian site, 21 images acquired between 30 December 2016 and 6 February 2017 from three relative orbits (7, 29, 131) were used. For the Spanish site, 18 images acquired between 21 August 2018 and 8 October 2018 from relative orbits 1 and 81 were used.

3. Methods

The following analyses and processes where carried away using the software GDAL/OGR [43], GAMMA software [44], GRASS (Geographic Resources Analysis Support System) [45], Python [46], Rasterio [47], Pandas [48], Geopandas [49], Numpy [50], Scipy [51], and Matplotlib [52].

3.1. DEM Assembly

The global DEMs were provided in equiangular geographic coordinates. The SRTM and the AW3D DEMs height reference had to be shifted from geoidal to ellipsoidal heights without resampling. The TanDEM-X (TDX) DEMs were provided as height above the ellipsoid. The Tandem-X DEM at 30 m (TDX30) was used as provided. The Tandem-X 12.5 m DEM was resampled (bilinear interpolation) to 20 m pixel spacing (TDX20) to reduce pixel size difference with respect to the multi-looked Sentinel-1 image. The 90 m resolution Tandem-X DEM (TDX90) was resampled to 30 m (bilinear interpolation). The ALS DEM, originally projected to ETRS89 UTM zone 30N coordinate system, was translated to a height above the ellipsoid and resampled to a 20 m pixel size (bilinear interpolation). As the ALS DEM has not been re-projected, all products geocoded with it (i.e., geocoded Sentinel-1 backscatter) share the same projection.

3.2. SAR Data Preparation

For each SAR image, the SLC sub-swathes were mosaicked, and the resulting image was multi-looked by a factor of 7 in range and 2 in azimuth. The objective was to reduce noise and obtain the SAR backscattered intensity at a pixel spacing close to the target 20 m used for the analysis. For a given orbit, the first acquired image was used as a master. All remaining SAR images from the same orbit were co-registered to the master image using an iterative process based on intensity matching and spectral diversity aided by each DEM [53]. For each orbit, the master image was used to generate a lookup table (LUT) relating map and range doppler coordinates. The LUT was used to orthorectify the master and the co-registered images (interferograms and SAR backscatter) from the same orbit.

Interferograms were generated for each consecutive image pair (a-b, b-c, c-d, etc.), and the DEM-estimated topographic phase was subtracted from each. The interferometric coherence was estimated in a two-step adaptive approach [54,55]. The first estimate of coherence was obtained with a 3-by-3 window. To reduce the estimation bias due to the small window size [56], the coherence was then recomputed using a window size inversely proportional to the initial estimate of the coherence. As a trade-off between preserving spatial resolution and reducing the bias, the largest window size was set to 9-by-9 pixels. In addition, when the estimation window included scatterers with a coherence level different than the coherence of the target in the center of the window, the estimator masked out such features to preserve the true coherence of the latter target.

The backscatter coefficient was calibrated to terrain flattened ${\gamma }^0$, considering the scattering area on the ellipsoid and on DEM surfaces ($A_{flat}$ and $A_{slope}$) [3] and the incidence angle on the ellipsoid and on DEM surfaces (${\theta }_{ref}$ and ${\theta }_{loc}$), as reported in equation 1 [57]. The parameter $n$ can be employed to account for volume effects [57]. The parameter was set to 1, an adequate value for most land cover types, as dealing with volumetric effects was not the objective of this study. The backscatter intensity and coherence images were orthorectified using the LUT and an inverse distance resampling. LUT coordinates located more than two pixels (range) apart to its counterpart on the SAR image were masked as no data. To reduce speckle, the multi-temporal backscatter images were averaged, by polarization, in time. Seven images were averaged for each orbit for the Romanian site, while nine images were averaged for each orbit for the Spanish site. Similarly, the six and eight coherence images were averaged in time for each orbit and polarization for the Romanian and the Spanish site, respectively.

$${\gamma }^0={\sigma }^0\frac{A_{flat}}{A_{slope}}{\left(\frac{\cos {\theta }_{ref}}{\cos {\theta }_{loc}}\right)}^n$$

(1)

3.3. Auxiliary Datasets

In support of the analysis, a land cover dataset was created for each site, based on the agreement between the ESA CCI land cover map (2015) [58], the DLR’s global urban footprint (GUF) (2016) [59,60,61], the ALOS PALSAR forest map (ALOS FNF) (2017) [62], and either Corine land cover map 2012 (CLC) [63] (Romanian Site) or Spanish information system on soil occupation (SIOSE, 2014) [64] (Spanish site, more detailed and recent) (

Table 2. Composition of the analyzed land covers based on preexisting datasets. When a higher level of Corine Land Cover or SIOSE has been employed (CLC (Corine land cover map) Lvl.1), the rest have been filled as “x”. CCI LC forest types are further disaggregated by the fractional cover and were therefore aggregated. GUF, global urban footprint; ALOS FNF, Advanced Land Observation Satellite forest map; SIOSE, Spanish information system on soil occupation; CODIIGE, Board of directors of the geographic information infrastructure of Spain.

	CLC 2012 [63]	SIOSE 2014 [64] (CODIIGE)	CCI LC 2015 [66]	GUF 2016 [59,60,61]	ALOS FNF 2017 [62]
Urban	1xx: Artificial surfaces	1xx: Artificial surfaces	190: Urban areas	Urban	-
Crops	21x: Arable land	210: Crops (herbaceous)	10, 20: Cropland, 11: Herbaceous cover	Other	Other
Pasture	23x: Pastures	320: Pastures	11: Herbaceous 130: Grassland	Other	Other
Grassland	321: Grassland	320: Pastures	11: Herbaceous 130: Grassland	Other	Other
Bare	33x: Open spaces	354: Bare	200: Bare areas	Other	Other
Broadleaf forest	311: Forest (broadleaf)	311: Forest (broadleaf)	50–62: Tree cover, broadleaved	Other	Forest
Needleleaf forest	312: Forest (needleleaf)	312: Forest (needleleaf)	70–82: Tree cover, needle leaved	Other	Forest
Mixed forest	313: Forest (mixed)	313: Forest (mixed)	90: Tree cover, mixed leaf type	Other	Forest
Water	5xx: Water bodies	5xx: Water bodies	210: Water	Other	Other

). ESA CCI land cover maps are generated at 300 m resolution using optical imagery time series from AVHRR, MERIS, SPOT-VGT, and POBA-V imagery, as well as GlobCover unsupervised classification chain and machine learning [58]. GUF is generated from TanDEM-X imagery at 12 m resolution based on amplitude and texture [59,60,61]. ALOS forest map is generated based on local thresholding of annual composites of ALOS PALSAR 1/2 amplitude images (25 m pixel side) [62]. Both CLC and SIOSE are based on photointerpretation of satellite and aerial imagery, with minimum polygon surface of 25 Ha [63] and 1 Ha [65], respectively.

The ALOS FNF disagreed with the remaining data sets as part of the cities were classified as forest, and part of the water was classified as “other” (not water, nor forest). For cities, no ALOS FNF condition was applied, whereas the rest of the non-forest classes on other datasets were considered compatible with non-forest classes from ALOS FNF (non-forest, water). The polygons with the agreement were dissolved by the land cover to eliminate internal borders, and a negative buffer of 40 m was applied to avoid edge effects.

The analyses were undertaken in this study are also related to landforms (Figure 3), i.e., features of the terrain surface with a distinct and identifiable shape [67]. Landforms were labeled using the GRASS GIS add-on “r.geomorphon” [68], with a search window of 25 pixels and a “flatness” threshold of 5 degrees applied to the highest spatial resolution DEMs available for each site, i.e., the TanDEM-X DEM at 30 m for the Romanian site and the PNOA DEM aggregated to 30 m for the Spanish site.

3.4. Inter-Orbital Data Analysis

The topographic normalization (radiometric terrain normalization, topographic phase removal) and distortion masking (e.g., foreshortening, layover, shadows) of each SAR-derived variable (backscatter coefficients and coherence) were assessed using the inter-orbit range (IOR). IOR was calculated pixelwise for each SAR variable by subtracting the maximum and the minimum values available from all orbits. For example, for a pixel with data available from orbits a, b, and c, the IOR would be $max\left(a, b, c\right) – min\left(a, b, c\right)$. The inter-orbit range was plotted by land cover class (boxplots). The analysis was repeated by landforms for the needleleaf forest, the only common forest type between both sites, and classes appearing near the mountain tops (grassland, bare). Because the only difference in image processing is the DEM employed, low IOR values reflect improved topographic effects removal.

In addition, the scattering area estimates, backscatter coefficient, and the interferometric coherence obtained using the ALS DEM (the most detailed) were employed as a reference to assess the performance of the global DEMs at the Spanish site. Said assessment was based on the root mean square deviation (RMSD, Equation (2)), relative RMSD (RMSDrel, Equation (3)), mean absolute deviation (MAD, Equation (4)), and Offset (Equation (5)) for the products obtained from the global DEMs.

$$RMSD=\sqrt{\frac{1}{P}{{\displaystyle \sum }}_{p=1}^P{\left(v_p-r_p\right)}^2}$$

(2)

$$RMSDrel=\frac{RMSD}{\overline{r}}$$

(3)

$$MAD=\frac{1}{P}{{\displaystyle \sum }}_{p=1}^P\left|v_p-r_p\right|$$

(4)

$$Offset=\overline{v}-\overline{r}$$

(5)

where $P$ is the total number of pixels, $p$ is a specific pixel and $v_n$ and $r_n$ are the variable values (i.e., backscatter) obtained for said pixel using a global DEM ($v$) and the reference DEM ($r$), whereas $\overline{v}$ and $\overline{r}$ are the mean value for said variables.

3.5. Land Cover Classification

A Linear support vector machine (LinearSVM) classifier was selected for its robustness and short execution times. We employed the Scikit-Learn implementation [69] with default options with a regularization parameter of 1, primal problem optimization, 0.001 tolerance for stopping criteria, and 10,000 iterations maximum. The classifier was trained per orbit/DEM pair using 96,000 samples, using VV- and VH-polarized backscatter and co-pol coherence as features.

For each orbit and land cover class, 70% of the valid sample (foreshortened and shadowed pixels were masked during SAR processing) was used to calculate the median and the median absolute deviation (MAD) of each SAR variable. Median and MAD values were then employed to calculate the z-score for each predictor (VV- and VH-polarized backscatter and coherence) by land cover class. Only samples with an absolute z-score below three were retained. Depending on the land cover class, the number of pixels retained varied from several millions (forest and low vegetation) down to tens of thousands (urban and water). For each class, 12,000 pixels (the number of pixels available for the less extended class, i.e., water) were randomly selected and used for training (n). As low vegetation and forest classes were further split into three sub-classes each (crops, pastures, grasslands, and broadleaf, needleleaf, and mixed forests), the total numbers of training samples selected were thrice as much (3n) as for urban and water land cover classes.

The validation sample, formed by the remaining pixels (30%) of each class, was used to compute the confusion matrix and associated error metrics (i.e., user and producer accuracies, Cohen’s Kappa, as described below). Error metrics for valleys were calculated after resampling (nearest neighbor) the landform layer to match the spatial resolution of the DEMs.

The confusion matrix C (Equation (6)) represents the occurrences of the predicted (rows) against the actual land cover class (columns) ($r·r$ dimensions, where $r$ is the number of classes). Diagonal cells ($c_{ii}$) count pixels with the same class in the classification and the reference dataset (True Positive, TP). Cells over the diagonal count pixels of class $i$ that have received other class (False Negative, FN), whereas cells under the diagonal count pixels that have been classified as $i$, when they have other class in the reference dataset (False Positive, FP). Accuracies for a specific class $i$ are the count of correctly classified pixels for the class ($c_{ii}$) divided by the number of pixels classified as $i$ (count across columns, $c_{i+}$), in the case of user accuracy ($U{A}_i$, also called ‘precision’ in machine learning literature, Formula (7)), or by the number actual $i$ pixels (count across columns, $c_{+i}$), in the case of producer accuracy ($P{A}_i$, also called ‘recall’ in machine learning literature, Formula (8)).

$$C=\left[\begin{array}{ccc}{c}_{1,1}& c_{1,j}& c_{1,r}\\ c_{i,1}& c_{i,j}& c_{i,r}\\ c_{r,1}& c_{r,j}& c_{r,r}\end{array}\right]$$

(6)

$$U{A}_i=precisio{n}_i={\sum }_{i=1}^r\frac{c_{ii}}{c_{i+}} =\frac{TP}{TP+FP}$$

(7)

$$P{A}_i=recal{l}_i={\sum }_{i=1}^r\frac{c_{ii}}{c_{+i}}= \frac{TP}{TP+FN}$$

(8)

Cohen’s Kappa [70] is a measure of agreement between the predicted cover and the one appearing on the reference dataset.

$$\widehat{K}=\frac{N{\sum }_{i=1}^r{c}_{ii}-{\sum }_{i=1}^r\left(c_{i+}·c_{+i}\right)}{N^2-{\sum }_{i=1}^r\left(c_{i+}·c_{+i}\right)}$$

(9)

where $r$ is the number of rows, $c_{ii}$ is the number of pixels where there is agreement between the classification and the reference dataset (cells on the diagonal, with row $i$ and column $i$), $c_{i+}$ and $c_{+i}$ are the totals for row $i$ (count of pixels classified as $i$) and column $i$ (count of reference pixels with class $i$), and finally, $N$ is the total number of observations [71]. A Kappa value of 1 represents the complete agreement between both, 0 represents a classifier performance similar to random guessing, and values under 0 indicate results worse than random guessing.

4. Results

In Section 4.1 and Section 4.2, we show results of the IOR analysis by land cover and land form, respectively, Section 4.3 contains the ALS reference-based analysis, and Section 4.4 describes the results of the classification comparison. Through Section 4.1, Section 4.2 and Section 4.3, coherence showed very small differences based on the DEM employed, under 0.01 for the IOR based analyses, and under 0.02 for ALS reference-based analysis. For this reason, these tables have been omitted, as they carried little to no information.

4.1. Inter-Orbital Range by Land Cover

The inter-orbit range (IOR,

Table 3. Backscatter Inter-orbit range (IOR) by polarization and land cover class at each study site (GL, grassland; BLF, NLF, and MLF, are broadleaf, needleleaf, and mixed forest). Cell color shows the gradient between the lowest (green) and the highest value (yellow). “M.D.” column represents the maximum difference between global DEMs for each specific land cover.

	Romania							Spain
	AW	SR	TDX	TDX	TDX	M.		AW	SR	TDX	TDX	TDX	ALS	M.
	3D	TM	20	30	90	D.		3D	TM	20	30	90		D.
Urban	5.4	5.4	5.3	5.4	6.9	1.6
Crops	2.3	2.3	2.2	2.2	2.3	0.1		1.3	1.4	1.3	1.3	1.5	1.2	0.2
GL	3.4	3.4	2.9	3.1	3.5	0.6
Bare								3.1	3.5	2.9	3.1	3.9	2.7	1.0
BLF	2.0	2.0	1.7	1.8	2.0	0.3
NLF	2.5	2.5	1.8	1.9	2.5	0.7		1.2	1.5	1.1	1.2	1.5	1.0	0.4
MLF	2.4	2.4	1.9	2.0	2.4	0.5

) was analyzed as an indicator of the residual terrain influence on the normalized SAR metrics (higher IOR, higher influence). The mean IOR for urban cover varied very little between DEMs, except for TDX90. Crops and broadleaf forests presented little difference (up to 0.3 dB) depending on the DEM employed, but increased for mixed forests (0.5 dB, Carpathians), and classes appearing near mountain peaks, such as grassland (0.6 dB) and bare soil (1 dB). Needleleaf forests had similarly high differences at the Romanian site (0.7 dB), whereas they were lower at the Spanish site (0.4 dB).

TDX20/30 was the global DEM with the lowest IOR, whereas the highest IOR values were observed for the TDX90 DEM, followed by SRTM DEM. Results with AW3D over natural covers (bare, grasslands, forest) varied depending on the site: At the Romanian site, its IOR was close to SRTM or TDX90 values; whereas at the Spanish site, it was closer to the values obtained with TDX20/30 DEMs. ALS-based PNOA DEM had the lowest IOR, with a large improvement over results with SRTM or TDX90 (up to 1.2 dB for bare areas, up to 0.5 for needleleaf forest), and a slight improvement over TDX20 results (0.2 dB for bare areas, and 0.1 dB for needleleaf forest). IOR values for PNOA and TDX20/30 DEMs showed less spread (Figure 4), concentrating around lower values, whereas the spread was larger for SRTM and TDX90. The main differences between sites were the lower spread of IOR at the Spanish site, and the behavior of AW3D IOR, which was closer to the TDX20/30 values at the Spanish site, but closer to SRTM/TDX90 values at the Romanian site.

4.2. Inter-Orbital Ranges by Landform

Analyzing IOR values by land cover may have dampened potential differences as all pixels, regardless of landform, were averaged by class. A more detailed analysis was conducted by disaggregating mean IOR values by landforms for the grassland, bare land, and needleleaf forests (

Table 4. IOR values disaggregated by landform for classes on mountain tops (grasslands, Romania; bare, Spain). Cell color shows the gradient between the lowest (green) and the highest value (yellow) for each row, which represents the IOR value for a specific landform when a certain was DEM employed. “M.D.” column represents the maximum difference between global DEMs for each specific landform.

	Romania (Grasslands)							Spain (Bare)
	AW	SR	TDX	TDX	TDX	M.		AW	SR	TDX	TDX	TDX	ALS	M.
	3D	TM	20	30	90	D.		3D	TM	20	30	90		D.
peak	3.1	3.0	2.5	2.6	3.5	1.0		2.9	3.3	2.6	2.7	6.1	2.2	3.5
ridge	3.0	2.9	2.5	2.7	3.0	0.5		2.7	3.1	2.5	2.7	3.6	2.3	1.1
spur	3.2	3.2	2.7	2.9	3.1	0.5		3.1	3.4	2.9	3.0	3.9	2.7	1.0
slope	3.6	3.6	3.0	3.3	3.6	0.6		3.3	3.7	3.1	3.3	3.8	2.9	0.7
hollow	4.2	4.3	3.6	3.9	4.6	1.0		3.2	3.8	3.1	3.2	3.9	2.8	0.8
valley	4.5	4.7	4.1	4.4	5.5	1.4		3.0	3.4	2.9	3.1	4.7	2.6	1.8

). These classes were selected as they often occur on steep slopes (Figure 4). Only landforms with more than 1000 pixels have been analyzed to limit spurious results, due to small sample size.

Both grasslands and bare areas presented large differences between DEMs. In the case of grasslands, the landforms valley, hollow, and peak showed the largest differences between DEMs (1.0–1.4 dB). In the case of bare areas, peak and valley showed the largest differences (3.5 dB and 1.8 dB, respectively), followed by ridge and spur (1.1 and 1.0 dB). For needleleaf forest, the largest differences between DEMs were observed for concave landforms (

Table 5. IOR values disaggregated by landform for needleleaf forests. Cell color shows the gradient between the lowest (green) and the highest value (yellow) for each row, which represents the IOR value for a specific landform when a certain was DEM employed. “M.D.” column represents the maximum difference between global DEMs for each specific landform.

	Romania							Spain
	AW	SR	TDX	TDX	TDX	M.		AW	SR	TDX	TDX	TDX	ALS	M.
	3D	TM	20	30	90	D.		3D	TM	20	30	90		D.
peak	2.3	2.2	2.0	2.1	2.3	0.3
ridge	2.2	2.2	1.9	2.0	2.2	0.4		1.2	1.3	1.1	1.1	1.3	1.0	0.2
spur	2.2	2.2	1.7	1.8	2.2	0.5		1.2	1.3	1.1	1.1	1.4	1.0	0.3
slope	2.3	2.3	1.7	1.8	2.4	0.7		1.2	1.5	1.1	1.2	1.5	1.0	0.4
hollow	2.8	2.8	1.9	2.0	2.9	0.9		1.4	1.9	1.3	1.4	2.1	1.1	0.8
valley	3.7	3.5	2.4	2.5	3.8	1.4		1.6	2.1	1.5	1.6	2.3	1.2	0.8
pit	4.5	3.9	2.6	2.7	4.0	1.9

): Hollow (0.8–0.9 dB at either site), valleys (up to 0.8 dB at the Spanish site, and up to 1.4 dB at the Romanian site) and pits (up to 1.9 dB at the Romanian site).

Using the TDX20/30 generally yielded the lowest IOR among global DEMs, whereas using TDX90 yielded the highest. The absolute minimum was always observed when using the ALS DEM (0.1–0.4 dB lower when compared to TDX20/30). IOR values obtained with the ALS DEM for bare areas were up to 3.9 dB lower when compared to the TDX90 DEM and up to 1.1 dB lower when compared with SRTM. For needleleaf forest, IOR values for valleys using ALS DEM were up to 1.1 dB lower than those obtained with TDX90.

4.3. Differences with and ALS-Derived DEM

The PNOA ALS-derived DEM provided the lowest IOR in all previous analyses at the Spanish site and was used as a reference for quantitative analysis of the global DEMs (

Table 6. Quality assessment for needleleaf forests using PNOA as a reference. Cell color shows the gradient between the lowest (green) and the highest value (yellow).

		O001						O081
	Statistic	AW	SR	TDX	TDX	TDX		AW	SR	TDX	TDX	TDX
		3D	TM	20	30	90		3D	TM	20	30	90
Sc. area (m2)	Abs. RMSD	97	188	102	109	213		115	214	113	123	245
	Rel. RMSD	11%	21%	11%	12%	24%		12%	22%	11%	12%	25%
	MAD	38	81	32	36	87		47	96	38	43	104
	Offset	18	−4	−11	14	–21		25	0.4	−12	20	−18
VV (dB)	Abs. RMSD	0.54	0.96	0.63	0.54	2.89		0.56	0.96	0.64	0.55	3.35
	Rel. RMSD	5%	10%	6%	5%	29%		6%	10%	7%	6%	34%
	MAD	0.38	0.64	0.41	0.36	0.74		0.38	0.63	0.41	0.35	0.75
	Offset	−0.19	−0.27	−0.27	−0.17	−0.35		−0.24	−0.30	−0.35	−0.23	−0.41

). For the scattering area, the highest deviation (RMSD) was observed for the TDX90, followed by the SRTM DEMs. The relative RMSD obtained with these DEMs was at least 10% higher when compared to the remaining DEMs (TDX20, TDX30, and AW3D). The lowest RMSD were observed for AW3D (orbit 1) and TDX20 (orbit81), followed by TDX30. The MAD for the scattering area was higher when using the SRTM or the TDX90 DEMs when compared to AW3D, TDX30, and TDX20. For both orbits, the SRTM-derived scattering area was the least biased when compared to the ALS DEM (−4 m² for orbit 1 and 0.4 m² for orbit 81), followed by TDX20 (−11 m² for orbit 1 and −12 m² for orbit 81).

For both orbits (1 and 81), the use of AW3D, TDX30, and TDX20 DEMs resulted in backscatter coefficient (VV) values closest to those obtained using the ALS DEM, with relative RMSD under 8%. For both orbits, the smallest offset with respect to the ALS DEM was obtained using TDX30, followed by AW3D. In all cases, backscattering coefficient was underestimated when compared with ALS DEM results.

4.4. Land Cover Classification

The DEM used for radiometric terrain normalization and topographic phase removal showed little effect on the overall quality of the classification, regardless of the Sentinel-1 relative orbit (

Table A1. Confusion matrices by DEM and orbit (Reference > columns; Classified > Rows) at overall level. (Ur, urban; LV, low vegetation; Fo, forest; Wa, water. UA, user accuracy, PA, producer accuracy).

		O007						O029						O131
		Ur	LV	Fo	Wa	UA		Ur	LV	Fo	Wa	UA		Ur	LV	Fo	Wa	UA
AW3D	Ur	19147	6602	7310	43	57.84		18698	6465	19317	24	42.01		18749	5044	12423	44	51.71
	LV	2206	1548500	19344	271	98.61		2063	1544334	20493	321	98.54		2366	1545097	24852	358	98.25
	Fo	250	10075	811265	380	98.70		759	13244	787900	395	98.21		357	9861	739556	389	98.59
	Wa	0	506	506	5054	83.32		0	611	355	5007	83.83		0	688	468	4956	81.09
	PA	88.63	98.90	96.76	87.93			86.89	98.70	95.15	87.12			87.32	99.00	95.14	86.24
SRTM	Ur	18868	6245	7204	40	58.31		18439	6055	19377	24	42.01		18725	5048	12952	48	50.92
	LV	2278	1540604	19090	262	98.62		2108	1537069	20272	318	98.54		2283	1538477	23950	340	98.30
	Fo	265	9900	806632	379	98.71		774	12931	782981	386	98.23		285	9564	735257	371	98.63
	Wa	0	546	549	4944	81.87		0	566	283	4895	85.22		0	618	358	4861	83.28
	PA	88.12	98.93	96.78	87.89			86.48	98.74	95.15	87.05			87.94	99.02	95.18	86.49
TDX20	Ur	29990	9446	11850	73	58.39		28857	8597	30849	36	42.23		29327	7184	18985	73	52.78
	LV	3168	2411865	22525	363	98.93		3440	2410056	31421	475	98.55		3426	2407209	39924	475	98.21
	Fo	364	16915	1270072	619	98.61		1088	17561	1188573	552	98.41		439	12460	1060469	540	98.75
	Wa	1	1112	419	7772	83.53		0	960	251	7535	86.15		0	1063	311	7085	83.76
	PA	89.46	98.87	97.33	88.05			86.44	98.89	95.00	87.64			88.36	99.15	94.71	86.69
TDX30	Ur	19156	6430	8530	42	56.08		18677	6280	21653	24	40.05		18793	4847	13386	46	50.69
	LV	2034	1540695	16591	239	98.79		2010	1538410	20705	301	98.53		2205	1539752	27089	303	98.11
	Fo	221	9505	808033	390	98.76		642	11234	773058	350	98.44		276	8255	712706	356	98.77
	Wa	0	665	321	4954	83.40		0	718	220	4867	83.84		0	680	171	4678	84.61
	PA	89.47	98.93	96.95	88.07			87.57	98.83	94.78	87.82			88.34	99.11	94.60	86.90
TDX90	Ur	19027	6599	7449	46	57.45		18439	5990	18774	26	42.65		18727	5161	13647	53	49.82
	LV	2179	1542343	18571	245	98.66		2138	1538972	19112	309	98.62		2294	1540926	23537	344	98.33
	Fo	240	9780	808189	396	98.73		794	13421	789226	382	98.18		309	9608	743312	370	98.63
	Wa	9	666	921	4937	75.57		0	567	343	4905	84.35		0	641	357	4854	82.95
	PA	88.68	98.91	96.77	87.78			86.28	98.72	95.38	87.25			87.80	99.01	95.19	86.35

). The Cohen’s Kappa was between 0.94 and 0.96. Analyzing the confusion matrices and the associated error metrics (i.e., user accuracy, UE, producer accuracy, PE) showed that classes with a larger spatial extent (low vegetation or forest) had very high accuracies (>95%). However, most pixels misclassified as low vegetation (>85%) were, in fact, forest pixels according to the reference data. Water had a reasonable accuracy (>80%). Depending on the DEM, the “source class” of pixels misclassified water pixels varied. When using TDX20 for normalization, 21–27% of misclassified water pixels were recorded as forest pixels on the validation dataset, 20–33% with TDX30, 33–50% with SRTM or AW3D, and 36–58% with TDX90. The urban class had relatively small omission errors (10–15%), but the commission error was high (around 40% for orbit 7, around 60% for orbit 29, and around 50% for obit 31), mostly due to the misclassification of forests (52–78% of the pixels misclassified as urban were forest pixels in the validation dataset).

Classification results were also analyzed by landform (

Table A2. Confusion matrices by DEM and orbit (Reference > columns; Classified > Rows) for valley. (Ur, urban; LV, low vegetation; Fo, forest; Wa, water. UA, user accuracy, PA, producer accuracy).

		O007						O029						O131
		Ur	LV	Fo	Wa	UA		Ur	LV	Fo	Wa	UA		Ur	LV	Fo	Wa	UA
AW3D	Ur	304	188	827	2	23.01		305	225	2478	1	10.14		298	175	1292	2	16.86
	LV	25	5447	4117	26	56.65		20	5232	4292	30	54.65		27	5214	3984	33	56.32
	Fo	2	537	76000	19	99.27		1	697	73596	18	99.04		3	606	70952	20	99.12
	Wa	0	32	176	514	71.19		0	27	81	512	82.58		0	27	54	506	86.20
	PA	91.84	87.80	93.69	91.62			93.56	84.65	91.48	91.27			90.85	86.58	93.01	90.20
SRTM	Ur	319	239	816	1	23.20		318	193	1974	2	12.79		318	177	949	1	22.01
	LV	29	5195	3821	18	57.32		30	5137	4540	28	52.77		29	5102	4404	28	53.35
	Fo	3	598	72313	17	99.15		0	689	70070	14	99.01		1	583	68204	15	99.13
	Wa	0	29	128	429	73.21		0	27	84	421	79.14		0	29	58	420	82.84
	PA	90.88	85.71	93.82	92.26			91.38	84.97	91.39	90.54			91.38	86.61	92.65	90.52
TDX20	Ur	495	410	1795	1	18.33		492	371	4221	2	9.67		485	277	1923	1	18.06
	LV	35	8150	3874	29	67.42		35	8080	5511	47	59.09		33	7983	5496	44	58.89
	Fo	3	994	113029	34	99.10		1	1022	104714	21	99.01		4	844	95033	26	99.09
	Wa	0	38	108	724	83.22		0	31	36	673	90.95		0	27	48	596	88.82
	PA	92.87	84.97	95.14	91.88			93.18	85.02	91.47	90.58			92.91	87.43	92.72	89.36
TDX30	Ur	329	274	1267	1	17.58		323	258	2910	1	9.25		324	205	1377	1	16.99
	LV	19	5160	2846	15	64.18		25	5091	3604	24	58.22		19	5068	3732	23	57.32
	Fo	3	594	72894	21	99.16		0	651	68921	12	99.05		1	553	64973	18	99.13
	Wa	0	33	71	428	80.45		0	29	51	411	83.71		0	29	33	378	85.91
	PA	93.73	85.13	94.57	92.04			92.82	84.44	91.30	91.74			94.19	86.56	92.67	90.00
TDX90	Ur	305	258	1046	0	18.96		304	243	2341	0	10.53		309	187	1256	0	17.64
	LV	24	5299	3304	21	61.27		22	5183	3809	39	57.25		21	5216	3669	39	58.31
	Fo	3	606	70640	19	99.12		2	725	68702	14	98.93		2	610	67360	12	99.08
	Wa	0	26	259	451	61.28		0	28	116	437	75.22		0	25	106	439	77.02
	PA	91.87	85.62	93.88	91.85			92.68	83.88	91.64	89.18			93.07	86.39	93.05	89.59

), and in particular, for valleys. Valleys were selected as they showed the largest differences between DEMs in previous tests, and a reasonable number of samples were available (10 times more when compared to the number of samples available for pit landform). When only valley pixels were considered, Cohen Kappa was low (0.57–0.70). Using the TDX20 DEM allows for a marginal increase of the Cohen Kappa values over the value obtained using the rest of the DEMs (0.05 for orbits 7, 0.03 for orbit 131, and 0.02 for orbit 29). For valleys, user accuracy for low vegetation class dropped. In this context, the use of the TDX20 DEM reduced commission errors up to 11% for low vegetation and up to 21.9% for water when compared to the remaining DEMs (Figure 5). Using the TDX20 DEM also reduced the number of forest pixels misclassified as water, representing a smaller percentage of the pixels misclassified as such (17–27% less, depending on the orbit). Commission errors for urban cover increased for valley landform, especially when terrain normalization is performed with any of the TDX DEMs.

The backscatter coefficient (VV) for forests located on valleys was examined to better understand the results obtained with the TDX DEMs (Figure 6). The boxplot showed that products normalized using lower resolution DEMs (AW3D, SRTM, TDX90) had an increased frequency of low values on the valley when compared to TDX20/30.

5. Discussion

The influence of the DEM employed for terrain normalization of backscatter and coherence data variability was analyzed in three ways: (a) Comparing several orbital tracks (inter-orbit range, IOR); (b) using the results obtained with an ALS-derived DEM as a reference; and (c) assessing land cover classification results after a specific DEM is employed for normalization. Coherence varied very little with the DEM employed, whereas the effect was larger on the backscatter coefficient.

Terrain normalization was better served by high-resolution DEMs (i.e., TDX20, ALS DEM, AW3D at Spanish site), in agreement with prior research [17,29] (SRTM-1arcsec and AW3D outperformed SRTM-3arcsec and TanDEM-X 90m). Higher resolution DEMs have reduced vertical uncertainties (under 5 m over sloping terrain for TDX12.5 and AW3D) when compared to the SRTM DEM (Table 1), which may have contributed to reducing IOR values. Some DEMs (TDX20, TDX90, ALS) needed resampling prior, which may have impacted their performance. TDX20 and ALS DEMs, were down-sampled, which may have reduced the DEM detail with an associated increase of IOR values. However, the IOR values obtained after resampling was still smaller than those observed for the lower resolution DEMs, underlining the importance of the vertical uncertainty of the original DEM. TDX90 was resampled to a finer resolution. However, as resampling is a destructive operation, the only expected impact was the smooth interpolation of the original data to a denser grid, which does not provide additional information over the original DEM. The advantage provided by high-resolution DEMs was dependent on the specific land cover, and the landform it occupies. For instance, IOR for urban and crops showed little difference, as they occupy near-flat areas. Land cover classes occupying steeper slopes (i.e., forests, grasslands) received the largest benefits of using a more detailed DEM (minimum IOR). At the Romanian site, broadleaf forests showed smaller differences than mixed and needleleaf forest, as the latter grew on steeper slopes. Variability for needleleaf forests was smaller at the Spanish site, as it occupied milder slopes (Figure 3). Results were disaggregated by landform, as the “typical” slope of each land cover might obscure landform related effects. Peak, hollow, valley, and pit landforms showed the largest differences between the analyzed DEMs. Such differences can be attributed to the sensitivity of SAR to remote sensing artifacts, such as shadowing and foreshortening appearing with increasing slope. This affects DEM accuracy in sloped terrain, as shown in Table 1, and propagates into any analysis based on pixel neighborhood, such as slope, orientation [12], or terrain normalization. On these landforms, TDX20/30 clearly outperformed the rest of the DEMs, as it provided an improved characterization of smaller terrain forms supporting results reported by Grohmann et al. [17].

The use of an ALS-derived DEM resulted in the smallest IOR, pointing it as a suitable candidate to benchmark global DEMs at the Spanish site. Taking the ALS DEM as a reference, the lowest deviation and bias of the SAR metrics were observed for the AW3D and TDX20/30 DEMs. The results obtained by AW3D were explained by the combined effect of the resolution employed for its generation (5 m) and the low cover and height of Mediterranean forests. These factors may have eased the detection of vegetation-free pixels, “pushing” the reported data nearer to the true terrain surface once resampled to 30 m, as described by References [17,18]. TDX20/30 provided similar results thanks to the high spatial resolution of the X-band sensor employed for its creation. In addition, shorter wavelengths (X- as opposed to C-band) can capture finer spatial details [1,14,72]. In some cases, better results were observed for the TDX30 DEM when compared to the TDX20 DEM, pointing to a trade-off between the spatial detail and the effect of the DEM noise, with some improvement for slightly coarser resolutions (30 m instead of 20), but a high dispersion when resolution becomes too coarse (TDX 90). TDX90 and SRTM DEMs showed similar dispersion when compared to the ALS PNOA DEM reference values, possibly due to the variable-resolution smoothing [20] employed during the reprojection of SRTM to map coordinates, which may have decreased its detail, as shown in References [17,21,22].

DEM performance also varied across sites, with a larger inter-track variability being observed for needleleaf forests at the Romanian site. Such differences were explained by the steeper slopes this land cover class occupied in the Carpathians, as well as its characteristics (height and structure), which may have complicated DEM generation, due to volume decorrelation, as reported by References [20,24]: The Carpathians are covered by dense temperate forests, whereas sierra Nevada is populated by Mediterranean forests with lower tree height and canopy density. Among all DEMs, AW3D IOR had a distinct behavior. While at the Romanian site, the AW3D results were close to those observed to the SRTM DEM, at the Spanish, the results were closer to those observed when using the TDX20/30 DEMs. The quality mask layer showed that, at the Romanian site, the AW3D DEM has a large strip where missing data (due persistent cloud cover) have been infilled with SRTM data, a problem mentioned by Truckenbrodt [29]. For this reason, it is not possible to draw conclusions on the influence of canopy characteristics on the performance of AW3D.

The characteristics of each DEM propagated into the land cover classification results. Overall classification accuracy was similar regardless of the DEM, with reasonable accuracies for all classes except urban. The low accuracy of urban surfaces was attributed to (i) the prevalence of steep slopes, which difficult terrain normalization (not fully accounted scattering area) [29] and may introduce DEM artifacts (mischaracterizing terrain surface) [12,14,72], and (ii) the prevalence of discontinuous urban fabric at the Romanian site, which may cause confusion with ornamental and fruit tree cover present around residential areas, and genuine forested lands.

When analyzing classification results for specific landforms (i.e., valleys), the AW3D, SRTM and TDX90 showed larger commission errors (CE), due to the misclassification of the forest as low vegetation or water. Such errors can be explained by the mischaracterization of the thin crevices of the drainage network on the DEMs [12]. This is propagated to the lookup table and impacts both, the distortion masking process and the terrain normalization. Distortion masking is affected because the distance between the pixels in the range is altered. Therefore, such pixels are not marked as distorted or shadowed. Terrain normalization is affected as it uses the LUT, orientation, and slope layers to estimate the scattering area. In turn, the scattering area is overestimated, under-compensating the radiometric effects, keeping a pseudo-shadow with lowered values (as opposed to true shadow, caused by occlusion, which cannot be compensated). Therefore, forests are misclassified due to the apparently low backscatter.

Using the TDX20 reduced commission error for low vegetation (11%) and water (12–22%) on valleys, as well as the percentage of forest pixels among all pixels misclassified as water (decrease of 17–27%). This was explained by reduced pseudo-shadows when using the TDX20/30 DEMs. However, the improvement came at the cost of a 3–6% increase of commission error for urban cover located on valleys, which was caused by a slight increase in backscatter values over forests (overcompensation), and the large size difference between the validation sample for urban and forest (<10³ vs. 7.5·10⁴ at valleys), as misclassification of a small subset of the latter would be much larger when compared with the sample size for the former. Even with this trade-off, classification results using TDX20 were marginally better (0.01–0.05 higher Cohen’s Kappa). Furthermore, increased CE for urban cover on valleys did not affect the overall CE for the urban class, which was reduced by 2–3% when using the TDX20 DEM.

6. Conclusions

SAR observations are heavily affected by sensor-terrain geometry, which can be corrected using a DEM. Choosing a DEM for SAR data terrain normalization is not a trivial choice, as it affects backscattering coefficient variability, and mapping products generated downstream. High-resolution TanDEM-X DEM (20 or 30 m resolution) was the global DEM providing the largest reduction of terrain induced variability, followed by AW3D in sparse vegetation areas. Natural land covers (i.e., forest, bare areas, grasslands) occupying steeper slopes and complex landforms (i.e., peaks, pits, valleys) received the largest benefits. These benefits were felt on classification, where more forest pixels were classified correctly due to a better compensation of low values (valley pseudo-shadow). An ALS-based DEM was able to provide slightly better results (i.e., marginally reduced IOR) when compared to AW3D and TDX20/30 DEMs. However, AW3D and TDX20/30 DEMs seem suitable candidates to replace ALS-based local DEMs. However, AW3D should be checked for data infilling from older datasets (i.e., SRTM) as over such areas, its performance may be degraded.

This study showed the effect of several global DEMs on terrain normalization, highlighting their advantages and shortcomings when normalizing Sentinel-1 imagery. Further research should expand this analysis by including the recent NASADEM dataset (from re-processed SRTM) [73], using a reference ALS-based DEM for temperate forests, and studying the DEM-dependent normalization effects on SAR imagery acquired at different wavelengths. Finally, the effect of terrain normalization could be tested on downstream quantitative products, such as biomass estimates or canopy cover.