Accounting for Deciduous Forest Structure and Viewing Geometry Effects Improves Sentinel-1 Time Series Image Consistency

Microwave scattering from forests generates pixel geolocation shifts in synthetic aperture radar (SAR) data that require an adequate representation within digital elevation models (DEMs) for preprocessing. We analyze the impact of DEM properties on the radiometry and geolocation of radiometric terrain corrected Copernicus Sentinel-1 imagery of forests to improve consistency in backscatter intensities for time series analyses. To account for the penetration depth of the C-band sensor, we approximate the structure of stands in a temperate deciduous forest using height percentiles from aerial laser scanning (ALS) point clouds in the Hainich National Park (HNP), Germany. Comparing the RTC results obtained using DEMs of SRTM, Copernicus, and ALS DEMs, the latter reduces topographically induced errors, resulting in visibly smaller effects from topography and spatially shifted information. Based on the P50 ALS vegetation elevation, the results show homogeneous intensities within the same orbit and reduce variance from 2.4 to 1.2 dB2 in the difference in mid-range data from ascending and descending azimuth directions. Over forest, we observe lower intensities on sensor-facing and increased intensities on away-facing slopes and correlations with the illuminated pixel area (IPA) and local incidence angle (LIA). We reduce this bias with linear regressions of intensity on IPA. ALS DEMs in RTC and the proposed regression correction increase the consistency of images across orbits, measured by the inter-orbit range (IOR), throughout the selected year at our study site. We suggest the proposed method applies to other areas, requiring further testing under different forest types and topography.


I. INTRODUCTION
R EMOTE sensing of forests relies on precise and recurrent measurements to link the observations by satellite-based Earth observation with the developments on the surface.The C-band synthetic aperture radar (SAR) of the Copernicus Sentinel-1 (S-1) satellites (S-1A, S-1B), which is able to penetrate clouds, as well as parts of the canopy, is widely used, for example, in mapping disturbances [1], forest types [2], response to drought [3], and forest biodiversity [4] among others.
Since 2015, S-1A, joined by S-1B from mid-2016 until the end of 2021, has been delivering global SAR coverage at high temporal and spatial resolution free of charge [5].Over Central Europe, up to four measurements are available within a six-day period due to the overlap of relative orbits.The shutdown of S-1B halved the temporal sampling rate, rendering consistent data from relative orbits important for denser time series.

A. SAR Preprocessing Challenges in Forested Areas
Especially under forest land cover, the path of the C-band SAR signal is influenced by a complex combination of structure and dielectric properties of the canopy [6], [7].With most of the signal reflected in the canopy, the scattering center's elevation ideally references the measured intensities the best [8], [9].Since S-1 is a monostatic sensor, the scattering center cannot be derived directly from S-1 data in densely forested areas due to temporal and volume decorrelation [10]; therefore, a precise external digital elevation model (DEM) is required for the transformation of S-1 data to map geometry for land cover analysis, for which previous research efforts in terrain correction of SAR [11], [12], [13], [14] provide the basis for the automatic geocoding of SAR imagery, further explained in [15] and [16].The available ground range detected (GRD) product with ∼20 × 22-m spatial resolution already accounts for the ellipsoid height in the geocoding process [17].Further SAR data processing to reduce distortion from the slanted data acquisition usually includes radiometric and geometric corrections [18].
It is proposed that vegetation height offsets of the used DEM used in preprocessing lead to an inaccurate geometric positioning and radiometric correction of signal intensities from forests.
1) Geometric correction references the signal by the time it travels between the sensor and the scattering center.Most remotely sensed DEMs account indirectly for the vegetation with their sensors' penetration depth which depends on the considered wavelength.During georeferencing, a mismatch This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. of the forest scattering center and an elevation model likely attributes the intensities to an incorrect surface area, toward the sensor if underestimated and vice versa.The information is shifted toward the sensor in the range direction from its actual geolocation.This shift, along with layover and foreshortening effects, is most prominent along the forest borders (Fig. 1).Opposing viewing geometries of ascending and descending orbits cause shifts in their respective range directions, doubling the offset that may be encountered in time series.With the out-of-the-box error of up to 4 m for S-1 interferometric wide-swath (IW) in the slant range due to atmospheric path delays [19], the analysis of pixel-level time series becomes challenging.Hence, we propose using more accurate information on the forest height in the form of DEMs that purposefully include the height structure of vegetation.Throughout this article we use the abbreviation DEM for terrain and surface models at several heights in the forest.
2) Radiometric corrections of the measured intensities, either angular-based [20] or area-based [18], rely on the angle or area of the DEM relative to the sensor's viewing angle and will likewise use an inaccurate surface angle or area in the calculation if the DEM underestimates the scattering center.In addition, volume scattering is prevalent in forested areas.Assuming a solid surface or its angle toward the sensor neglects the actual process and overcorrects the resulting backscatter intensity.

B. DEMs Used in SAR Preprocessing for Forest Studies
Many global elevation models generated to measure topographic height, such as SRTM-1 HGT (SRTM DEM) [21], are affected by the forest structure and height.Offsets between these elevation models and the more recent vegetation structure encountered by the S-1 SAR sensors affect data processing.However, a survey of recent studies using S-1 shows that SRTM-based corrections are commonplace for forest-aimed studies in on-site processing and data provided by Google Earth Engine (GEE) [1], [2], [22], [23], [24].
Few studies use more detailed national DEM, digital terrain models (DTMs), or digital surface models (DSMs), but they do not compare geocoding results from different DEMs over the same area [25], [26], [27].Like the SRTM DEM, DTMs and DSMs may under-or overestimate the canopy height perceived by the sensor, introducing deviations between the DEM and the acquired signal (Fig. 1).Recently, the use of top-canopy elevation models proved essential analyzing drought at pixel level [28].
Direct comparisons of DEM resolution and error have been performed with Seasat, detecting the largest displacement of measured intensities at topographic peaks [29].More recently, Truckenbrodt et al. [30] highlighted inconsistencies in RTC backscatter generated with DEMs of ASTER, SRTM, and TanDEM-X (TDX) missions with 30-and 90-m spatial resolutions processed with different software.Using the DEMs for radiometric terrain correction, they found the best compensation for the incidence angle when using the SRTM DEM.Borlaf-Mena et al. [31] compare the similarity of S-1 scenes from different viewing angles for different land cover using DEMs of SRTM, Advanced Land Observing Satellite (ALOS) World 3D, and TDX.They used the inter-orbit range (IOR, explained in Section III-B3), measuring the range of intensities from different relative orbits per pixel, using a DTM generated from aerial laser scanning (ALS) data as a reference and noticed a considerable dependence of backscatter intensity on the choice of the DEM.However, the ground elevation was used without focusing on vegetation, omitting overlying structures and their impact on backscatter.Also, the recently published Copernicus GLO-30 DEM (COP DEM) [32], generated based on TDX data [33], has not yet been used in comparative studies of S-1.

C. Radiometric Corrections in Forests
Radiometric corrections are performed to reduce topographically induced SAR artifacts and to harmonize data in various ways.However, the imagery of S-1 is not handled uniformly.For example, orbits are kept as separated data variables of the same area without correction for terrain effects [22].S-1 σ 0 intensities are normalized to a joint viewing angle [2], [34], which is used to combine different orbits of similar incidence angles [23], [35], or within images from one azimuth direction [36].Hoekman and Reiche [37] use different angular normalization functions depending on vegetation type.In GEE, angular corrections are also available [20], [38].The radiometric terrain correction (RTC) [18] has gained much interest in recent times for reducing topographic effects using the illuminated pixel area (IPA) to calculate γ 0 RTC , with more notable effect in mountainous areas [39].For example, the RTC renders intensity information more comparable to combining overlapping relative orbits as temporal composites [40].Following the improved correction with area-based methods [41] and the implementation within the planned S-1 analysis-ready data (ARD) processing chain of the normalized radar backscatter product [42], the RTC is preferred as the preprocessing method.

D. Research Approach
We propose a method of including the forest structure based on ALS-derived forest height percentiles to reduce displacement and radiometric artifacts in the S-1 time series for forest analyses.Starting from the observation that currently used global DEMs introduce systematic errors in preprocessing due to an inaccurate representation of the forest's scattering center, we hypothesize that backscatter differences between relative orbits depend largely on DEM impact within the correction.We compare the resulting intensity based on ALS-derived processing to datasets (DSs) obtained using DEMs from SRTM and Copernicus.First, we address DEM-related variations by analyzing S-1 data processed to γ 0 RTC while exchanging the DEM used for radiometric and geometric correction.After visual exploration, we show the relationships of backscatter intensities with properties of the viewing geometries from different DEMs.Our proposed regression normalization uses partial residuals from ordinary least squares (OLSs) to remove the bias explained by the viewing geometry impacting γ 0 RTC backscatter intensities to obtain normalized very dense γ 0 RTC,norm SAR time series suitable for seasonal analysis.The consistency of images from different relative orbits over the year is measured using the IOR on a monthly basis.With a better representation of the forest structure, we expect to be able to reduce the shift in the intensity information of the canopies and improve the overall radiometric detail in pixel-level dense time series analyses.

A. Study Area
To assess the impact of viewing geometry on backscatter intensity, we chose an area of interest (AOI) in the Hainich National Park (HNP), Germany (Fig. 2).The HNP comprises mostly unmanaged old-growth temperate forests of deciduous broadleaved beech (Fagus sylvatica) and ash (Fraxinus excelsior) trees [43].The 2 × 1-km AOI contains mainly forest and some grassland.Within the AOI, the mean forest height was 23.6 m (derived from ALS) in 2017, with a distinct new-growth forest parcel at the center, directly southwest of the northeastern meadow.The new-growth stand has a lower mean forest height of 19.2 m, visible in the COP DEM and the P50 DEM in Fig. 3.The area is chosen for its relatively flat ground topography to reduce the impact of topography on the SAR signal and focus on forest height differences instead.Since the most recent ALS measurements in the AOI were acquired in early 2017, which coincides with the first full year of S-1 A and B data, we limit our analysis to that year.The annual mean temperature and precipitation height are 9.5 • C and 774.2 mm, respectively, measured at the nearest DWD weather station "Eisenach," 10-km south of the AOI [44].

B. Elevation Data
As elevation data (Table I), we select SRTM DEM [45] as a reference for its widespread use in recent studies.In addition, we include the COP DEM [32] since its global coverage, open data availability, and more recent acquisition date make it a suitable alternative to the SRTM DEM.It is worth noting that the COP DEM is based on TDX data with 12-m pixel spacing, which has a smaller penetration depth and therefore results in higher elevation estimates in forested areas than SRTM's C-band.To account for the vegetation structure, we use the point clouds of open-access ALS surveys (GDI-Th, Freistaat Thueringen, TLVermGeo).The ALS surveys were carried out during the winter of 2017 in leaf-off conditions for deciduous trees and preclassified for general land cover types containing a vegetation class.The measured ALS points attributed to vegetation, corresponding to the trees' trunks and branches, are later used for the elevation percentile calculation.
The three DEMs differ strongly in their ability to represent the forest (Fig. 3).The SRTM DEM shows a blurry image, making it difficult to visually separate forest from grassland.The COP DEM has a much better-defined forest area, but the ALS-derived P50 DEM provides the most detail.The SRTM DEM generally estimates a lower forest elevation, falling between the DTM and P10 from the ALS point cloud.This places SRTM for most of the AOI below the canopy (Fig. 1), which the acquisition can explain during the leaf-off period with steep sensor incidence, a large temporal gap concerning possible forest change, and a longer wavelength of the sensor.The COP DEM generally lies between P30 and P70 of the ALS vegetation elevation.

C. SAR Data
Copernicus S-1 C-band dual-polarimetric IW SAR data of all the relative orbits over the HNP for 2017 is downloaded as GRD [17] from the Alaska Satellite Facility [46].Two ascending (44 and 117) and two descending (66 and 168) orbits are available, acquired at approximately 18:00 and 6:30 CET, respectively for the AOI, with an incidence angle difference of about 8 • between the acquisition directions (Table II), resulting in a combined mean temporal sample rate of 1.5 days.

D. Ancillary Data
Ancillary vector data from the digital base land cover map (DLM, derived from the German land cover catalog, ATKIS) are used.Specifically, nonforest areas are excluded using vegetation layer 2, containing deciduous, coniferous, and mixed forest types with a geometric precision of 15 m (GeoBasis-DE/BKG 2022) [47].According to DLM classification, the study area comprises only deciduous forests.

A. Data Preparation
The overall workflow of the study is given in Fig. 4. First, all the DEMs are resampled to a common 10 × 10-m grid.This spacing was chosen as a suitable denominator between the ∼20-m S-1 GRD product and the 30-m global elevation products, as well as being on line with our interest in analyzing highly resolved time series, where knowledge of information shift may be crucial [28].Then, the S-1 GRD scenes are geocoded to γ 0 RTC with the GAMMA software [48] using the different DEMs.Forest areas are then clipped.A regression of intensity on IPA is fit per scene, and its coefficient is applied to the intensity in the proposed regression normalization.IORs are finally calculated from monthly averages for both γ 0 RTC and γ 0 RTC,norm . 1) DEM Creation: To consider and analyze the impact of vegetation height in RTC and geocoding, we use various DEMs.The SRTM and COP DEMs are processed with the pyroSAR (v.0.10.1)function dem_autocreate in Python [49], which mosaics, corrects from geoid to ellipsoidal height, crops, and resamples the elevation data to a chosen area, resolution, and coordinate reference system.The elevation data are resampled bilinearly and projected to WGS84 UTM32N at a 10 × 10-m pixel spacing, also used for ALS data aggregation for consistency.
The elevation models derived from the preclassified ALS point cloud are created for ground elevation and vegetation elevation percentiles from 0 to 100 in steps of 10% using the point data abstraction library (PDAL) [50].The aggregation process includes limiting the point cloud to one flight path Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.with the lowest deviation from the nadir per cell of the aforementioned 10 × 10-m grid to reduce the impact from the ALS's sensor incidence angle of ±22 • and higher point densities from overlapping swaths.Then we calculate the elevation percentiles weighted by the intensity of the discrete returns per cell.With a small AOI selected for low topography covered in one ALS survey flight, the intensities were comparable, which may not hold for larger areas.The weighting favors more strongly reflecting elements in the forest structure, resulting in slightly increased elevations.To avoid data gaps in areas with low or no vegetation, grid cells with fewer than 50 points attributed to vegetation are filled with the mean ground elevation from the point cloud.Finally, the elevation data are corrected to the ellipsoidal height, and a 3 × 3 mean filter is applied, as the variability in the ALS DEMs leads to missing values in the RTC otherwise.Then, the ALS-derived elevations are converted to rasters matching the SRTM and COP DEMs.In the following, the ALS DEMs are referred to by their respective percentile, e.g., P50 DEM.
2) SAR Processing: The S-1 IW GRD SAR DSs are processed with each of the DEMs used for both the area-based radiometric and geometric corrections in the RTC in GAMMA (GAMMA RS, v. 20210701) [48] to γ 0 RTC with 10 × 10-m pixel spacing [18] using the Python module pyroSAR [49].The detailed workflow is provided in the pyroSAR documentation [51].The resulting intensity of VV and VH polarizations in decibels and the associated viewing geometry properties, local incidence angle (LIA), and IPA are converted into a map geometry matching the DEM grid.LIA describes the angle (in degree) between the pixel's surface normal and the sensor look direction in range.IPA additionally considers the angle in the azimuth direction and provides the illuminated area of the pixel in square meters.
In the following, the SAR DSs are named after the DEM used for preprocessing, i.e., SRTM DS, COP DS, Ground DS, and P0-P100 DS for preprocessing with the various ALS DEMs.
3) Filtering to Homogeneous Forest: Focusing on the forested areas, we mask out nonforest areas using the DLM land cover vectors.As the masked areas still include nonforest pixels at the border and small clearings within the forest, density-based spatial clustering of applications with noise (DBSCAN) [52] is used for further classification.This is performed in the 4-D feature space spanned by polarizations VV and VH, IPA, and LIA.Using this method with a temporal one-year mean, the largest resulting cluster provides the prevalent land cover pattern, which is, in our case, closed forest cover.Other smaller clusters and noise are excluded from the subsequent analysis, as they are assumed to belong to other land cover classes.
The filtering is performed per relative orbit to arrive at one common set of pixels for each DS, retaining only those pixels that belong to the largest cluster in all the orbits.The filter is again combined for comparisons between DSs, representing ∼29% of the AOI.
4) Proposed Regression Normalization: We observe linear correlations in both the viewing geometry variables LIA and IPA.IPA is more fitting as a predictor variable in linear regression models of γ 0 RTC as the area is a physically more meaningful variable, and better results are reported with area-based correction measures [18], [41].Then the intensity values are corrected by eliminating the effect of IPA using where s is the estimated slope coefficient for a given S-1 scene calculated by the regression between IPA and the backscatter.γ 0 RTC refers to the RTC backscatter intensity in either VV or VH polarization.The linear regressions were fit with an intercept term using the implementation in scikit-learn (1.0.2) [53].
In (1), we normalize the intensities to each DS's mean IPA using IPA − IPA for the adjustment.This allows us to keep the mean intensity response of the AOI and only eliminate the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
effect of deviating from that mean.Note that IPA is specific to the AOI within a scene.

B. Analysis
A seasonal average of the DSs from July to September is created for visual analysis (more in detail in Section III-B1).Then, the observed intensity bias and its reduction are given by comparing the seasonal average's γ 0 RTC and γ 0 RTC,norm relationships to LIA and IPA, as well as the aspect and slope of the respective DEM (more in detail in Section III-B2).We calculate the IOR each month across the four orbits for all the DSs.Finally, the lowest IOR per month is calculated by interchanging the DEM in preprocessing for each relative orbit (more in detail in Section III-B3).
1) Visual Assessment of DEM Sensitivity: We first analyze the spatial patterns of the sensitivity of the preprocessed intensity to DEM choice.We compare γ 0 RTC backscatter intensities in VH polarization, seasonally averaged from July to September (leaf-on) for the SRTM, COP, and P50 DS of orbit 117 (asc) and their difference with orbit 66 (desc).Orbits 117 and 66 were chosen because of their similar mid-range incidence angle for both the directions, enabling a comparison with minimized influence of the LIA.
2) Relationships of Intensity and Viewing Geometry Properties: Second, we investigate dependencies in the DS on the viewing geometry properties of the SAR data and the DEMs.γ 0 RTC and γ 0 RTC,norm are analyzed to determine whether the proposed regression normalization reduces the existing dependencies.Density scatterplots of the relationships between intensities and IPA, LIA, aspect (relative to north, 0 • -360 • clockwise), and slope are compared, focusing on the forest area.Aspect and slope are calculated from the DS's respective DEM using gdal.DEMProcessing with its default "Horn" algorithm.To summarize general trends, median intensities are calculated within bins along the x-axes.
3) Across-Orbit Consistency Analysis: Finally, we investigate which DEM improves the consistency of the SAR data across different orbits by calculating the monthly averaged IOR with, on average, five scenes per orbit [31].We examine the usually applied approach of using a single DEM to preprocess all the relative orbits.Hypothesizing that the scattering center depends on viewing geometry, acquisition time, and vegetation density, the best representation to reference S-1 intensity might change over the year.To investigate whether the consistency is further increased when using forest height from different DEMs over the year and for different orbits, we calculate the lowest IOR per month by ALS DS combinations.Therefore, we 1) calculate monthly backscatter intensity averages of each orbit and ALS DS; 2) iterate over all the combinations of monthly averages; and 3) calculate the average IOR over the AOI for each combination.
The IOR gives the intensity range per pixel by subtracting the minimum from the maximum across the considered orbits.It is formally defined as where γ (x, y, o) is the backscatter intensity at the location x, y in the image from orbits o.
A lower IOR is, here, associated with a better-fitting DEM, producing similar backscatter values regardless of the orbit.We calculated the IOR per pixel on monthly backscatter intensity averages for each DS with all the four available relative orbits as input.The monthly IORs are finally spatially averaged within the masked forest area.The number of orbits compared and the period for averaging the backscatter intensities impact the IOR and must be consistent for direct comparisons.With differing pixel spacing, the number of compared orbits, AOI size, and averaging period, our IORs are not directly comparable to the results presented by [31].

A. Visual Assessment of DEM Sensitivity
Fig. 5 shows clearly that seasonal mean intensity images were most consistent within one and between different mid-range orbits [117 (asc) and 66 (desc)] in VH polarization when preprocessed with the P50 DEM and least consistent when using the less detailed SRTM DEM.This tendency is visible both within forest stands and along their edges: Intensity differences between ascending and descending orbits are reduced in magnitude and affected area [Fig.5(b)].The bright sensor-facing forest edge, most visible in the SRTM DS, can be attributed to the smoothed height change in the SRTM DEM from meadow to forest.With a shallow slope at the forest edge in the DEM used for radiometric correction, the intensity is not corrected for layover and foreshortening effects and therefore shows greater intensity values.Similar results were obtained in VV polarization (not shown).
Within the forest area [Fig.5(a)], evident features are marked with white arrows.Central in the north of the AOI, marked by the upper arrow, a V-shaped feature is visible in the SRTM DS.The two diagonal low-and high-intensity ridges translate to areas of ≥5 dB difference between the ascending and descending azimuth directions.Directly below (south), marked by the second arrow, another ridge is visible in the SRTM and COP DSs at a new-to-old-growth forest boundary.In the COP DS, this appears as a broad blue and slim red ridge in the intensity difference image [Fig.5(b)].Both the features are visible as small mounds in the respective elevation data (Fig. 3) and are greatly reduced in the P50 DS in Fig. 5(a) and (b).Prominent in SRTM and COP DEM are small clearings in the right (east) of the AOI, which are much more present in the intensity differences in Fig. 5(b) than the P50 DS.The lowest row in Fig. 5 shows the distribution of intensity values and across-orbit differences in the above images.For one orbit, the minimum and maximum of SRTM (short dashed) and COP DSs (dashed) extend over the P50 DS (full line), which is, in turn, slightly more centered at −14 dB.The intensity differences in the ascending and descending orbits are more closely centered at 0 for P50 DS compared with SRTM and COP DSs as well.The variance in the difference image decreases, with 2.4, 2.0, and 1.2 dB 2 for SRTM, COP, and P50 DS, respectively, which confirms the visual impression that preprocessing with the P50 DEM improves the consistency between different orbits.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Relationships of Intensity and Viewing Geometry Properties
Before applying the proposed regression normalization, the average July-September VH-polarized backscatter intensities show strong linear trends with the properties describing the viewing geometry (Fig. 6).LIA shows a positive linear trend and a mean increase of 1 dB for the three depicted DSs, processed with the SRTM DEM, the COP DEM, and the ALS-derived P50 DEM.IPA shows a negative linear trend with a mean intensity decrease of 1 dB over the value range.The relationship for the aspect variables resembles a sine curve, peaking at 90 • and reaching its minimum at 250 • -270 • .With the sensor's viewing direction toward 75 • in relative orbit 117, we observe a higher backscatter intensity over slopes facing in the look direction than toward the sensor.The backscatter intensities do not show a consistent and substantial slope dependency.Note that aspect and slope angles refer to elevations within the vegetation structure, not generally to the ground level.The results obtained for the VV polarization and other relative orbits are similar and are therefore not shown.
After applying the proposed regression normalization, calculated per scene from the correlation of intensity with IPA, γ 0 RTC,norm backscatter intensities no longer display IPA-related trends (Fig. 7).The trend in LIA is also successfully removed in the P50 DS, but a weak positive linear trend is still present for LIA in the SRTM and COP DSs.The sine-shaped pattern concerning the aspect has also vanished, perhaps except for the SRTM DS.

C. Across-Orbit Consistency Analysis
The monthly IOR of the DSs shows that most ALS DSs have higher consistency during the summer period, further expanding over the year after the proposed regression normalization Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.(Fig. 8).Overall, P50 DS provides the best results before the proposed regression normalization, whereas Ground, COP, and SRTM DSs only surpass it in IOR for some winter months.The proposed regression normalization improves the IOR overall, most visibly in the winter months of P30-P90.Within all the DSs, a similar temporal pattern is visible: high consistency with three minima in January, May, and October, and maxima in February, July, and November with lower consistency between the orbits.
Comparing the IOR of P50 DS to the better IOR using data processed with varying ALS-derived DEMs shows that there is only marginal further improvement (Fig. 9).The lowest IOR values are often found with DEMs P20-P60 over the observation period.VH polarization shows higher consistency than VV polarization, except in February and July.

A. Impact of Forest Representation in Preprocessing of SAR Data
We show that the S-1 backscatter is impacted by the DEM choice, introducing visible features with differences of more than 5 dB [Fig.5(b)] between viewing geometries, and correlates with the viewing geometry by ± 1 dB (Fig. 6).Using more detailed and recent DEMs increases the consistency in SAR imagery from one orbit and reduces the differences between ascending and descending orbits.Lessened over-and underestimations of the forest height of the P50 instead of the SRTM DEM reduce the observed patterns [Fig.5(b)] which comply with [34] of different azimuth angles [34] and enhance them on a finer scale.We interpret the visually flatter intensities over the forest area of the P50 DS to a better result in the radiometric correction and geocoding of the SAR data through more fitting elevation estimates of the SARs' scattering center within the forest area, refining the area-based radiometric normalization [18], [41].
The patterns observed between the viewing geometry properties, the aspect of the DEM, and γ 0 RTC hint at an overcorrection of radiometry within vegetation volume.Frey et al. [41] show a negative trend between European Remote Sensing Satellite (ERS) intensity σ 0 and LIA, no trend of γ 0 with LIA, and Truckenbrodt et al. [30] show negative trends between S-1 γ 0 and LIA.The sine-shaped pattern over the aspect with a minimum for sensor-facing slopes within the continuous forest is contrary to the results of [20], [37].Typically, higher intensity values are observed at sensor-facing slopes [18], [41].Focusing on a small subset of a homogeneous forest area in the leaf-on period with flat ground elevation, the trend is inverse: larger intensities with increasing LIA and lower intensities with increasing IPA, and the intensities are generally lower Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.for sensor-facing slopes.The RTC process uses IPA without accounting for the land cover in its calculation.However, the change rate concerning IPA might differ in surface or volume scattering.The mismatch of these interactions might explain the inverted trends, with IPA in vegetation volume with less impact on the measurements than anticipated in the RTC.
The increased consistency between viewing geometries of ascending and descending orbits from ALS DEMs supports the hypothesis that deviations in intensity can be reduced if the scattering center is well-approximated.The best results in the IOR are provided using the P50 DEM (Fig. 8), which fits centrally within the canopy.With similar median height, the difference in across-orbit consistency between COP and P50 DS can be explained by the representation of the forest height variability.With a mean elevation difference of ∼9 m (Table I), elevation estimates from SRTM fall between DTM and P10 DEM and are for most of the AOI below the canopy and produce the highest IOR overall.Similarly, an increase in IOR can be seen in using the P90 DEM (Fig. 8), which deviates toward the upper boundary of the canopy.
The lowest IOR from different ALS DEMs per orbit (Fig. 9) shows that while the P50 DS performs almost identically, the choice in the DEM combinations empirically shows the hypothesized shift of the scattering center during the year.Vegetation density is increased during summer by the presence of leaves and diurnally changing vegetation water content, which increases signal attenuation.This is reflected in the DEM combinations, with an arc during the leaf-on period and the tendency for similar DEM choice for the azimuth directions and acquisition time of day.The descending orbits match better with higher ALS percentiles compared with the ascending orbits.This tendency could be verified against in situ information on water allocation in trees along the diurnal cycle.Although at high variance, the remaining months tend to lower percentiles, implying a good scattering center approximation in the lower canopy under leaf-off conditions.
Compared with the findings in [31], we do not see an improvement in IOR when using the DTM over SRTM in preprocessing for the deciduous forest.The temporal pattern for all the considered DEMs may limit the comparison of IOR to the same season.Lower IOR values achieved with the ALS percentiles and the proposed regression normalization show the improvement by our approach.

B. Proposed Regression Normalization for Forested Areas
The proposed regression normalization improves the IOR with a larger impact on the DSs preprocessed with ALS DEMs.The pattern of IPA and backscatter intensities (Fig. 6) indicate a radiometric overcorrection of the RTC instead of  the normalization over viewing geometries shown in other studies [20], [30], [37], [41].For robust measurements over large areas or time steps from various viewing geometries, we show that the removal of systematic variations in the backscatter (Fig. 7) increases across-orbit consistency (Fig. 8) [18].By considering the target land cover and normalizing intensities to the predominant IPA, we consider different scattering responses instead of normalizing a scene to a Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
predefined incidence angle.While the yearly temporal average generates the land cover filter per orbit here, coefficients are processed individually for each scene.The proposed regression normalization is applicable for single scenes to reduce DEM artifacts in vegetation and multiple scenes in time series.
With the preservation of the mean overall intensity in (1) by normalizing around the backscatter response of IPA, the impact of large area factors such as rainfall and forest characteristics remain in the data.This also includes offsets between sensors noted in [54], which can be tackled by coupling the mean backscatter of several scenes of the same relative orbit over time at the cost of temporal detail.Differences induced by vegetation water content due to different acquisition times of ascending and descending orbits have been shown to impact backscatter intensities systematically [55], [56], if present should also remain within the backscatter information.Effects of precipitation or temperature were not investigated but reduced by taking temporal averages.Finally, we cannot account for the path taken through the vegetation canopy by the C-band signal, which is different for each viewing angle and acquisition due to the randomness of volume scattering [7].As such, measurements will display some degree of difference for each relative orbit per pixel that geocoding cannot resolve.
Application on larger areas using ALS-percentile DEMs with the proposed regression normalization must be extended for other forest types.Similar to [37], the best fitting percentile is dependent on forest structure, and the proposed regression normalization needs to be adjusted for areas with different trends of intensity to IPA.An alternative to the ALS percentiles could be modeling the scattering volume with the available ALS point cloud, but it was outside this analysis's scope.Because of the observed linear trend of IPA to backscatter intensity in Fig. 6 within the Hainich AOI, a linear regression was deemed fit for the proposed regression normalization.Adjustments for larger areas and varying topography could be a local average instead of IPA in (1), masking of pixels displaying layover, most easily detected by unrealistically large IPA, or nonlinear regression.

VI. CONCLUSION AND OUTLOOK
In this work, we have done the below.1) Halve the across-orbit variance in S-1 γ 0 RTC intensities using DEMs derived from ALS point clouds instead of SRTM-1HGT in deciduous forest area.2) Show improved preprocessing results with Copernicus GLO-30 compared with SRTM-1HGT.3) Observe linear correlations between backscatter intensity and IPA, LIA, and a sine-shaped pattern with aspect with all the DEMs for preprocessing.4) Improve consistency among intensity images with different viewing geometries by the proposed regression normalization.We conclude the following.1) A viewing geometry normalization is beneficial for time series involving different orbits to reduce systematic radiometric overcorrection in areas with strong volume scattering.
2) The more detailed ALS DEM at the 50th elevation percentile improves spatial fit, reduces topography-induced artifacts, and improves internal radiometric consistency of forest within a single orbit.
3) The best across-orbit consistency is achieved using RTC with ALS DEMs and the proposed regression normalization.In forestry applications, the proposed approach enables better interpretation with reduced SAR artifacts, increased data consistency, and, consequently, a higher temporal sample rate in time series using combinations of relative orbits.This outcome is of critical importance due to the recent loss of Sentinel-1B.However, the proposed bias corrections should be extended and tested on a more extensive range of forest types with different volume scattering behaviors for a wider application.

Fig. 1 .
Fig. 1. (Left) DEM cross sections overlaid on the ALS point cloud.(Right) seasonal mean backscatter intensity from RTC with (red) SRTM and (green, blue) ALS P50 DEMs in a small portion of the AOI.Based on data GDI-Th, Freistaat Thueringen, TLVermGeo, NASA JPL, Copernicus.

Fig. 3 .
Fig. 3. Elevation models of the AOI.(Top to bottom) ALS DTM (mean ground-point elevation aggregated to 10-m resolution); SRTM DEM resampled to 10 m; COP DEM resampled to 10 m; and P50, the median elevation of ALS vegetation points aggregated to 10 m, with gaps from low or no vegetation cover filled with ground elevation.Forest structure is better represented by the COP and ALS DEMs than the SRTM DEM.Based on data from GDI-Th, Freistaat Thueringen, TLBG, NASA JPL, DLR e.V. 2010-2014 and Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by the European Union and ESA; all rights reserved.

Fig. 4 .
Fig. 4. Schematic workflow of the study.Parallelograms: data; rectangles: processing steps; pill shapes: analysis steps.Data are one year of S-1, covering the AOI.Except for the spatial filter, processes are performed per scene.Analyses are performed on seasonal (B1, 2) and monthly (B3) data aggregates.

Fig. 5 .
Fig. 5. (a) Sentinel-1 seasonal mean γ 0 RTC intensity image from orbit 117 (asc) in VH polarization, July-September.The red arrow shows the range direction; white arrows point to features discussed in Section IV-A.(b) Difference in Sentinel-1 seasonal mean γ 0 RTC intensities between orbit 117 (asc) and orbit 66 (desc) in VH polarization, July-September (positive: intensities of orbit 117 greater than intensities of orbit 66).The fourth row shows the empirical distribution of values in the images above.ASF DAAC 2023 contains modified Copernicus Sentinel data 2017, processed by ESA.

Fig. 6 .
Fig. 6.Relationships between VH γ 0 RTC orbit 117 July-September (leaf-on) mean intensity data and viewing geometry properties related to SAR processing (LIA, IPA) and DEM (aspect, slope) for SRTM, COP, and P50 DS.Markers (horizontal lines) show mean intensity values within intervals of the variable shown on the x-axis.Note that ranges of slope angles vary.

Fig. 7 .
Fig. 7. Relationships between VH γ 0 RTC,norm orbit 117 July-September (leaf-on) mean intensity data, corrected with the proposed regression normalization, and viewing geometry properties related to SAR processing (LIA, IPA) and DEM (aspect, slope) for SRTM, COP, and P50 DS.Markers (horizontal lines) show mean intensity values within intervals of the variable shown on the x-axis.Note that ranges of slope angles vary.

Fig. 8 .
Fig. 8. Monthly IORs for S-1 VH polarization of deciduous forest in the AOI.(a) γ 0 RTC and (b) γ 0 RTC,norm processed with the SRTM DEM, COP DEM, and the ALS-derived DEMs for different elevation percentiles and the ground level.

Fig. 9 .
Fig. 9. Lowest IOR for monthly VH and VV γ 0 RTC,norm intensity averages of deciduous forest in the AOI.Each month, the lowest IOR is determined across all the orbits processed using ALS DEMs.(a) IOR values of backscatter intensities in VH and VV polarizations to the monthly IOR of P50 DS.(b) and (c) Indicate which DEM achieves the best IOR in each polarization and month across all the obits.For example, in January the smallest IOR in VH polarization is achieved using P40 DS for orbits 44 and 168, P50 DS for orbit 117, and P60 DS for orbit 66.

TABLE I ELEVATION
DATA FOR THE AOI

TABLE II SENTINEL
-1 DATA FROM THE YEAR 2017 OVER THE AOI