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

: 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 a ﬀ ected by cloud cover and solar illumination and can capture unique information about its targets. However, SAR observations are a ﬀ ected by the coupled e ﬀ ect 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 di ﬀ erence between orbit tracks, the di ﬀ erence with results obtained with a high-resolution local DEM, and the impact on land cover classiﬁcation. 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 di ﬀ erences 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.


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 Figure 1. Extent of the study areas (a, Romania, Ro; c, Spain, Sp) and digital elevation models (DEMs) used for synthetic aperture radar (SAR) data processing. The yellow box indicates the location of the area covered by the Sentinel-1 dataset used for the analysis. The red box indicates the extent of the subset shown in the right-hand side panels (b, Romania; d, Spain).
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 Figure 1. Extent of the study areas (a, Romania, Ro; c, Spain, Sp) and digital elevation models (DEMs) used for synthetic aperture radar (SAR) data processing. The yellow box indicates the location of the area covered by the Sentinel-1 dataset used for the analysis. The red box indicates the extent of the subset shown in the right-hand side panels (b, Romania; d, Spain).
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.
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. Figure 2. Distribution of the slopes for grassland, bare soils, and forest land. The latter has been separated according to leaf type. 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 2016/12/30 and 2017/02/06 from three relative orbits (7,29,131) were used. For the Spanish site, 18 images acquired between 2018/08/21 and 2018/10/08 from relative orbits 1 and 81 were used.

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 multilooked Sentinel-1 image. The 90 m resolution Tandem-X DEM (TDX90) was resampled to 30 m Figure 2. Distribution of the slopes for grassland, bare soils, and forest land. The latter has been separated according to leaf type.
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 2 . 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. 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.

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.

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 γ 0 , considering the scattering area on the ellipsoid and on DEM surfaces (A f lat and A slope ) [3] and the incidence angle on the ellipsoid and on DEM surfaces (θ re f and θ 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.

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 Remote Sens. 2020, 12, 3016 7 of 23 (CLC) [63] (Romanian Site) or Spanish information system on soil occupation (SIOSE, 2014) [64] (Spanish site, more detailed and recent) ( Table 2). 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. 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. 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.

CLC
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.
(a) (b) Figure 3. Landform classification. The general shape of the landforms (a), modified from GRASS documentation, based on [68]) and landforms over the shaded relief for a subset of the Spanish study area (b).
where is the total number of pixels, is a specific pixel and and are the variable values (i.e., backscatter) obtained for said pixel using a global DEM ( ) and the reference DEM ( ), whereas ̅ and ̅ are the mean value for said variables.

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

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(a, b, c) -min(a, b, c). 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.
Remote Sens. 2020, 12, 3016 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 v and r are the mean value for said variables.

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 (UA 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 (PA i , also called 'recall' in machine learning literature, Formula (8)).
Cohen's Kappa [70] is a measure of agreement between the predicted cover and the one appearing on the reference dataset.K 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.

Results
In Sections 4.1 and 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 Sections 4.1-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.

Inter-Orbital Range by Land Cover
The inter-orbit range (IOR, Table 3) 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). 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.  (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.

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

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). 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. 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 ( 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): 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). 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. 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.

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). 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 2 for orbit 1 and 0.4 m 2 for orbit 81), followed by TDX20 (−11 m 2 for orbit 1 and −12 m 2 for orbit 81). 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). Statistic  AW  SR  TDX  TDX  TDX  AW  SR  TDX  TDX  TDX   3D  TM  20  30  90  3D  TM  20  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.

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). 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), 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. 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. Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 24 Figure 6. Boxplot representing the VV backscatter coefficient for forests located on valleys by Sentinel-1 relative orbit: Mean value (triangle) median value (orange) and inter quantile-ranges (whiskers) for 5-95%.

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 highresolution 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].

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 (Xas 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 3 vs. 7.5·10 4 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.

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. Funding: This research was funded by the Romanian National Authority for Scientific Research and Innovation and the European Regional Development Fund under the project "Prototyping an Earth-Observation based monitoring and forecasting system for the Romanian forests" (EO-ROFORMON, project ID P_37_651/105058).
Acknowledgments: TanDEM-X DEM data were provided by the German Aerospace Center (DLR), under the proposal DEM_FOREST2614. We also would like to thank the publishers of the open datasets employed in this study: DLR for making available the TanDEM-X Global Urban Footprint dataset, European Environment Agency for providing the Corine Land Cover Dataset, European Space Agency for providing Sentinel-1 data and CCI land cover map, National Geographic Institute of Spain for providing PNOA LiDAR DEM and SIOSE land cover dataset, JAXA, for providing both AW3D DEM and ALOS forest/non-forest map, and USGS, for providing SRTM. Finally, we would like to thank INCDS Marin Dracea and University of Alcalá for hosting the research activities.

Conflicts of Interest:
The authors declare no conflict of interest.