Large-Scale Forest Height Mapping by Combining TanDEM-X and GEDI Data

The present study addresses the development, implementation, and validation of a forest height mapping scheme based on the combination of TanDEM-X interferometric coherence and GEDI waveform measurements. The very general case where only a single polarisation TanDEM-X interferogram, a set of spatially discrete GEDI waveform measurements, and no DTM are available is assumed. The use of GEDI waveforms to invert the TanDEM-X interferometric measurements is described together with a set of performance criteria implemented to ensure a certain performance quality. The emphasis is set on developing a methodology able to invert forest height at large scales. Combining 595 TanDEM-X scenes and about 15 million GEDI waveforms, a spatially continuous 25-m resolution forest height map covering the whole of Tasmania Island is achieved. The derived forest height map is validated against an airborne lidar-derived canopy height map available across the whole island.


Large-Scale Forest Height Mapping by Combining
TanDEM-X and GEDI Data

I. INTRODUCTION
A CCURATE forest height measurements at subhectare scales are critical for characterizing the successional state of forests and/or their disturbance regime. At the same time, forest height can be related to forest biomass through allometric models, so used to initialize (or constrain) model estimates of above-ground biomass [1], [2], [3], [4]. Despite their importance for forest inventory and modeling, forest height measurements on the ground remain difficult, so the generation of accurate high spatial resolution forest height maps over large areas remains a remote sensing challenge.
The introduction of polarimetric synthetic aperture radar (SAR) interferometry (Pol-InSAR) at the end of the 90s was a decisive step toward measuring forest height accurately at large scales [5], [6], [7], [8], [9], [10]. Relying on the inherent sensitivity of the interferometric coherence to the vertical structure of volume scatterers, Pol-InSAR techniques have been established alongside lidar measurement techniques for accurate estimation of forest height in the context of airborne and spaceborne applications [11], [12], [13], [14]. From the proposed Pol-InSAR estimation algorithms, model-based ones have in general proved more robust and with a better performance. Such estimation algorithms model the vertical reflectivity in the forest by means of a two-layer model (accounting for vegetation and a ground scattering contribution) [5], [15], and forest height is then obtained by inverting the established model using interferometric measurements at different polarizations or vertical wavenumbers (e.g., spatial baselines) [10], [16]. Launched in 2010, TerraSAR-X add-on for Digital Elevation Measurement (TanDEM-X) introduced a new era in spaceborne radar remote sensing allowing single-pass interferometric measurements from space in a bistatic configuration [17], [18]. The sensitivity of the interferometric TanDEM-X coherence to the vertical forest structure and especially to forest height initiated a large number of studies on forest height estimation from TanDEM-X SAR interferometry (InSAR) data across all possible forest types and conditions [15], [19], [20], [21], [22], [23], [24]. However, since the conventional TanDEM-X global digital elevation model (DEM) observation mode is limited to a single polarization interferometric acquisition, only a highly simplified implementation of a single-layer model is possible, which often requires additional external information, such as an external digital terrain model (DTM), for its inversion [16], [25]. The achieved performance is, in general, remarkably good as long as the forest conditions allow sufficient penetration to ensure the "visibility" of the whole (vertical) forest extent, and the forest (and terrain) heterogeneity can be matched by the simplified single-layer parameterization of the vertical reflectivity. While the limited penetration is physical and must be accepted as long as the underlying topography is unknown, the limited ability of the inversion model to adapt to the local forest and terrain conditions can be addressed if a larger observation space is available [16].
The wide availability of airborne lidar data triggered several attempts to use lidar measurements to compensate for the underdetermination of the forest height inversion problem when addressed in terms of single-polarimetric TanDEM-X acquisitions. Besides using the lidar-derived DTM to directly This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ enable the forest height inversion [22], [23], [25], [26], [27], lidar data have also been used to constrain individual model parameters reducing, in this way, the dimensionality of the forest height inversion problem [28], [29], [30], [31], [32], [33]. More recently, the direct use of lidar waveforms to define the full X-band vertical reflectivity profile has been proposed [34]. This allows, depending on the spatial density of the available waveform measurements, adaptation to local forest and terrain conditions and improved performance, especially in spatially heterogeneous forest conditions.
The launch in late 2018 and subsequent operation of NASA's Global Ecosystem Dynamics Investigation (GEDI) Mission, has been a critical development for large-area forest height estimation [35], [36]. GEDI provides dense and well-distributed lidar waveform measurements across Earth's tropical and temperate forests and defines an ideal framework to explore the synergetic use of waveform lidar and interferometric X-band SAR measurements at local to global scales. Indeed, the synergies between TanDEM-X and waveform lidar data [34] have been confirmed by several studies using the TanDEM-X/GEDI framework [30], [31], [37], [38].
This study investigates the potential of this synergy for largescale forest height mapping for the case where only a single polarisation TanDEM-X interferogram, a set of spatially discrete GEDI waveform measurements, and no DTM are available. The emphasis is set on developing a methodology able to invert forest height at large scales. Accordingly, decisions and tradeoffs are taken in favor of optimizing the data handling and processing flow and the overall inversion performance accounting for minimizing the effect of systematic errors. This is done accepting that the performance achieved may be locally inferior compared to more sophisticated, but also more complex, inversion approaches. Section II describes the selected forest site for the experiments, i.e., Tasmania in Australia, the experimental data (TanDEM-X and GEDI), and the reference height measurements available. Section III addresses the estimation of forest height based on the combination of TanDEM-X and GEDI data. In Section IV, the forest height estimation performance is addressed. A performance model is introduced and used to define the optimum height inversion range. Section V describes the actual processing and main results over the test site. Finally, in Section VI, the key findings are summarized and the conclusions are drawn.

A. Tasmania
The island State of Tasmania, located in the south of the Australian mainland, has an area of 68401 km 2 with an extent of about 350 km in longitude and 330 km in latitude. The island is predominantly forested by temperate rainforests and has mountainous as well as flat regions. The forested area is very dynamic with fire events, logging activities, and regrowth, which makes the test site attractive for evaluating different forest conditions.

B. Lidar Canopy Height Model
From 2008 to 2019 the Government of Tasmania performed several airborne laser scanning campaigns to map a large part of the island. One of the products of these campaigns was a Canopy Height Model (CHM) of 2-m spatial resolution [39]. The CHM data, mosaicked together to a single raster dataset are shown in Fig. 1(a) together with a map that indicates the acquisition year for each area in Fig. 1(b). Although the CHM model results from data acquired over a period of about 10 years, the availability of such a large and continuous forest height dataset allows the evaluation of the proposed methodology not only at a local but also at larger (country-wide) scales.
For the validation of the TanDEM-X heights associated with the upper canopy height, the H100 height [16] (corresponding to the mean height of the 100 tallest trees within a hectare or equivalent to the mean of the 6 tallest trees within a 25 m × 25 m window) was derived from the CHM as the mean of the tallest 3 CHM values (accounting for their 2-m spatial resolution) within a 25 m × 25 m window. The obtained H100 values correspond well to GEDI's RH98 heights, as shown in Fig. 2.

C. Forest Non-Forest Map
To avoid inconsistent heights over areas where the (inversion) model is not applicable, the TanDEM-X-derived forest/nonforest (FNF) map was used to mask water and urban areas [40]. For all remaining areas, forested and non-forested height estimates are computed and accounted for in the validation.

D. TanDEM-X Dataset
In this study, 595 conventional Coregistered Single look Slant range Complex (CoSSC) data products acquired and processed in an HH polarized strip map mode during the global DEM [18] and change DEM [41]. TanDEM-X mission phases between 2011 and 2019 have been used. The dataset includes ascending and descending acquisitions with a typical resolution of about 3 m in range and azimuth and with vertical wavenumbers ranging from 0.05 to 0.15 rad/m (with a median of about 0.1 rad/m). Their coverage map is shown in Fig. 1(c).

E. GEDI Dataset
The GEDI instrument consists of three lasers (operating with a wavelength of 1064 nm) that are split into 4 beams and dithered along-track to produce a total of eight parallel acquisition tracks, spaced approximately 600 m apart on the Earth's surface in the cross-track direction. Along each acquisition transect, waveforms with approximately 25-m footprints are spaced every 60 m along the track. About 15 million GEDI footprints acquired in the 18 first months of the mission between April 2018 and October 2019 are used in this study. Their spatial distribution is shown in Fig. 1(d), covering an area of less than 1% of the whole island.
From the available GEDI footprints, only those with a signalto-noise ratio sufficient to penetrate through the forest with up to 95% canopy cover have been used [42]. The geolocation accuracy of the footprint data has been improved to ∼10 m at 1-sigma in Version 2 [43]. The GEDI waveforms and the Level 2A RH98 values [36] have been used to initialize the height estimation from the TanDEM-X coherences.

III. FOREST HEIGHT INVERSION
This section introduces the forest height inversion methodology based on the combination of TanDEM-X and GEDI data.

A. TanDEM-X Coherence Estimation
The (complex) interferometric coherenceγ Obs is obtained from the interferometric image pair s 1 and s 2 as where <·> denotes the expectation value. In the TanDEM-X global DEM and change DEM mission phases, s 1 and s 2 are acquired in a bistatic interferometric strip map mode where one of the two satellites transmit (in H polarization) and both satellites receive the scattered signal quasi simultaneously (in H polarisation). The absence of temporal decorrelationγ Obs comprises two main decorrelation contributions The first term, i.e.,γ Sys , includes the decorrelation effects induced by the nonideal SAR system. The most prominent system decorrelation contribution is the (additive) noise decorrelation γ SNR . Modeling the received signal to be composed by the scattering amplitude a and the noise amplitude n, i.e., s = a + n, γ SNR can be written as [16] where SNR = A/N is the signal-to-noise ratio, with P = A + N the received power, A = |a| 2 the scattered power, and N = |n| 2 the noise power.
A second contribution to the system decorrelationγ Sys is γ Quan , which is caused by the lossy raw data compression process [44], [45]. In the case of TanDEM-X, the received backscattered signal is first digitized by an 8-b analog-to-digital converter and then further compressed by a block adaptive quantizer with a compression rate of 8:3 (or 8:4). This is associated with a coherence loss of about 3.5% (or 1%) [45], [46]. The second term, i.e.,γ Scat , reflects the phase stability of the scatterer under the different incidence angles induced by the interferometric baseline. After range and azimuth spectral filtering [5]γ Scat reduces to the volume decorrelation contributionγ Vol [5], [6]γ where F(z) is the vertical reflectivity function (also referred to as the vertical reflectivity profile) expressing the vertical distribution of scatterers seen by the interferometer. The vertical direction is expressed by z and z 0 indicates the position of the underlying ground. Accordingly, F(z) depends on the frequency and polarisation of the interferometer as well as on the interferometric acquisition geometry. The upper bound of F(z) is given by z 0 + h V which in the case of a forest scatterer corresponds to the (top) forest height. The lower bound of F(z) is given by z 0 . For the bistatic TanDEM-X acquisition, the vertical (interferometric) wavenumber κ z is defined as [47] where θ is the local incidence angle given as the difference between the (nominal) incidence angle θ 0 and the terrain slope angle in the range direction α, λ the wavelength, and Δθ the change of the incidence angle induced by the spatial baseline (i.e., the spatial separation of the two satellites). The vertical wavenumber is often expressed by the so-called height of ambiguity HOA = 2π/κ z , i.e., the height corresponding to an interferometric phase of 2π. The local incidence angle is estimated using the terrain slope angle in the range direction α derived from the TanDEM-X DEM. Before using the interferometric coherences for inversion, γ SNR and γ Quan need to be compensated. The compensation of γ SNR is performed using the noise equivalent sigma zero (NESZ) patterns along the range provided for each of the two SLC images (i.e., TSX and TDX) in the CoSSC data product [16]. Accordingly, the SNR for each of the two SLCs is calculated from the NESZ pattern, and the individual backscattering coefficient σ 0 (6) and the SNR-induced decorrelation is then obtained as [16] γ SNR = 1 Using the γ SNR values obtained from (7) and assuming the γ Quan accounts constantly 3.5% of the decorrelation (γ Quan = 0.965), the volume decorrelation contributionγ Vol is obtained asγ After calibration, the interferometric coherences are ready to be used for inversion. Any volume decorrelation values greater than 1 (due to the inherent standard deviation of interferometric coherence) are set to 1.

B. Mean Vertical Reflectivity Profile
Following the approach proposed in [34], the GEDI waveforms P i (z) are used to approximate the vertical reflectivity profile F(z) in (4). Because GEDI's waveform measurements are spatially discrete a mean reflectivity profile, a representative for a larger area, has to be used. The generation of such a mean reflectivity profile from the available GEDI waveforms follows the approach proposed in [34]. First, the available waveforms P i (z) normalized to unit height and resampled to a common number of height samples are used as columns of the so-called profile matrix [P]. Accordingly, the number of rows of [P] is given by the sampling in the height dimension and the number of columns is the number of available GEDI waveforms n.
From the profile matrix [P], a covariance matrix [R] is formed (where [·] T indicates the transpose operator) and diagonalized where where N < M represents the number of eigenvectors used to compose F m (z). In the following, N = 1 is used so that F m (z) = P 1 (z). Considering only the dominant profile component and omitting higher-order eigenvectors correspond to an attempt to represent all forest states within the scene by a common lowfrequency structural component. This can be justified by the fact that higher-order structural components play often a secondary role as the structural dependency on the volume coherence is weaker compared to its dependency on height. Finally, the use of F m (z) in (4) allows the estimation of forest height h V for eachγ Vol (κ z ) sample. However, there are several points that have to be considered. The mean reflectivity profile is estimated for each TanDEM-X scene individually using all available GEDI waveforms within the scene. Fig. 3(a) shows an example of a mean reflectivity profile F m (z) representative for a whole TanDEM-X scene obtained by combining the more than 5000 GEDI waveforms located within the scene. Due to the height normalization and the subsequent resampling of the GEDI waveforms when forming profile matrix [P], the mean reflectivity profile becomes a long "tail," i.e., an asymptotic behavior toward higher heights, especially when the scene is dominated by short(er) forest stands. This, if not accounted for, biases the obtained forest height estimates later. To avoid this, the mean reflectivity profile is limited by an intensity threshold: the profile is cut at the height where the (profile) intensity decreases below a given (arbitrary) threshold (3 dB in the following) with respect to the highest located reflectivity maximum. After this, the profile is re-normalized between 0 and 1 (i.e., from 0 to 1) to obtain the final F m (z) [see Fig. 3 Furthermore, the use of a single F m (z) for all (forest) samples within the scene implies an intrinsic error depending on how well F m (z) represents the underlying reflectivity profile F(z) and its spatial variability within the scene. One way to address this problem is to reduce the area in which F m (z) is estimated. In this sense, one can segment a TanDEM-X scene in multiple areas estimating for each a mean reflectivity profile. However, there is a tradeoff between the number of GEDI waveforms available and the size of the area to be represented by the mean reflectivity profile: with decreasing area size, the number of GEDI waveforms within the area decreases making the estimation of the mean profile less robust. On the other hand, the fact that the mean reflectivity profile is estimated and used on a TanDEM-X scene basis can lead to inverting the same sample in different scenes with different mean reflectivity profiles. This can be the case for partially overlapping TanDEM-X scenes.

C. Modification of Profiles With Height
The high attenuation at lidar and X-band radar frequencies implies for both GEDI and TanDEM-X measurements a high sensitivity to the "visible" geometric architecture of the forest volume. This, together with the comparable spatial resolutions of the two measurements, makes the GEDI waveforms P(z) and the X-band vertical reflectivity profiles F(z) for many forest types similar. It is this similarity that justifies the use of P(z) as an approximation for the vertical reflection profile F(z). However, it is also clear that GEDI waveforms and the X-band vertical reflectivity profiles are not equal. The differences are induced primarily by the different acquisition geometries in which the two configurations operate.
The nadir-looking GEDI geometry causes a stronger ground contribution-especially in open canopy forest conditionsthan in the TanDEM-X reflectivity, where the side-looking geometry results in a larger path through the canopy and thus a stronger attenuation of the ground contribution. This leads to an overrepresentation of the ground contribution in the mean reflectivity profile that can bias the obtained height estimates.
The overrepresentation of the ground contribution can be clearly seen in Fig. 3(b) where the mean reflectivity profile has a much stronger ground than canopy contribution. While in open canopy forest conditions, the ground contribution is overrepresented, in denser and/or closed canopy forest conditions the ground contribution in the GEDI waveforms is weaker than the one expected at X-band (vertical) reflectivity. In this case, the ground contribution in the mean reflectivity profile is underrepresented. This makes the direct use of the mean profile F m (z) as derived from the GEDI waveforms suboptimum. Equally suboptimum, in terms of performance, is also the use of the same profile for the whole height range, i.e., the assumption of the same distribution of scatterers for all forest heights (from 5 to 70 m).
To account for both effects, the overrepresentation of the ground contribution and the variation of the distribution of scatterers with forest height, an empirical height-dependent attenuation correction on F m (z) is performed where ε(z) is a height-dependent attenuation factor The constant attenuation factor ε 0 expresses the difference in the effective attenuations of the side-looking X-band and the nadirlooking lidar propagation and is universally fixed at 0.1 dB/m for all scenes. This is, of course, a strong claim that cannot be confirmed by any theoretical or experimental means. However, the main function here is to attenuate the overrepresented ground scattering contribution in the mean profile F m (z) derived from the GEDI waveforms, especially for smaller stands, rather than to provide an accurate model. In this sense, any similarly chosen attenuation value would have a very similar effect.
One should note here that the difference between the two profiles, the TanDEM-X and the GEDI one, is more a question of forest density than of forest height. However, the absence of any large-scale structure or density information makes it impossible to account for this. On the other hand, the significant errors introduced by the overrepresented ground contribution make its generic attenuation [as performed by (12) and (13)] in terms of overall performance more effective even if for several cases this general correction leads to worse estimates.
Once the attenuation factor ε 0 and the reference height h Ref are defined, (12) provides for each mean reflectivity profile F m (z) a set of profiles F mh (z), one for each height. Fig. 3(c) shows an example of a set of 70 profiles each associated with a forest height from 1 m up to 70 m.
The use of F mh (z) in (4) allows now the estimation of forest height h V for eachγ Vol (κ z ) sample. The inversion can be implemented by means of a simple look-up table (LUT) that maps the absolute value of the volume decorrelation to a forest height accounting for the actual (local) vertical wavenumber κ z . Fig. 4(a) shows an example of such a volume decorrelation-forest height LUT for different vertical wavenumbers κ z .

IV. HEIGHT INVERSION PERFORMANCE ANALYSIS
To evaluate the quality of the obtained forest height estimates on a sample basis, a set of performance criteria are introduced that allow to mask suboptimum height estimates and compensate for systematic (height) offsets.

A. Lower Coherence Level
Lower coherence values are associated with a higher (amplitude and phase) variance that propagates into the height estimates. To avoid such estimates, a coherence threshold has been introduced to mask out the heights obtained from samples with a lower (absolute) coherence. A threshold value of 0.3 has proven to be appropriate and is used throughout the following.

B. Vertical Wavenumber Performance
Knowledge of the reflectivity profile allows a simplified forest height inversion performance evaluation accounting for two key performance effects: 1) inversion biases induced by a residual nonvolumetric decorrelation contribution affecting primarily the lower height ranges; and 2) the reduced sensitivity (or even saturation) to higher heights induced by too large or too small vertical wavenumbers. For a given forest height h V and reflectivity profile F mh (z), the expected interferometric volume coherenceγ Vol (κ z ) can be estimated using (4). The effect of residual (nonvolumetric) decorrelation contributions can be addressed by applying a residual decorrelation contribution γ Res to the volume decorrelatioñ γ Vol (κ z ) and by inverting the distorted coherence γ Res ·γ Vol (κ z ) (instead ofγ Vol (κ z )). The height estimates obtained for the whole height range are then compared to the original ones to estimate the height-dependent bias induced by γ Res . This estimation bias can be visualized by plotting the estimated heights against the reference heights as shown in Fig. 5  The obtained estimation bias can be used to define the optimum forest height performance range by means of an upper and lower limit for each individual height estimate using the reflectivity profiles F mh (z) and the associated vertical wavenumber κ z . The upper and lower limits of the optimum forest height performance range are defined by means of a performance threshold that can be the same or different for the upper and lower case. For example, in Fig. 5, a 20% error threshold for the lower heights is indicated by the vertical dark green line, and a 10% error threshold for the higher heights is indicated by the vertical light green line. According to the plots for a vertical wavenumber κ z ≈ 0.1 rad/m (middle), this leads to a lower height limit of about 12 m and an upper height limit of about 58m.
The reduced sensitivity at higher heights induced by too large vertical wavenumbers is accounted for by means of the height derivative |γ Vol (κ z , h V , F mh (z))| / h V for a given vertical wavenumber κ z . The height derivatives for the profile assumed in Fig. 3(c) are plotted in Fig. 4(b) clearly demonstrating the loss of sensitivity with increasing height and/or vertical wavenumber. The associated loss of performance can be accounted for by introducing a threshold on the derivative. A sensitive (yet conservative) way to do this is to accept height estimates only up to the height that corresponds to the minimum of the |γ Vol (κ z , h V , F mh (z))| / h V derivative. This height limit is indicated by the light blue vertical line in the performance plots of Fig. 5(b). For a vertical wavenumber of κ z ≈ 0.1 rad/m, this leads to an upper height limit of about 40 m. This is the height above that the inherent sensitivity of the observation configuration (expressed by the vertical wavenumber) degrades. The estimates mentioned above are rejected.
For defining the upper height limit of the validity range, the most restrictive height limit of the two performance criteria is selected. Of course, there is always a tradeoff between applying less restrictive thresholds to allow more valid pixels at the price of potentially higher height errors. Accordingly, for this case, only height estimates between 15 m and 40 m are accepted.

C. Global Bias Correction
The proposed performance model accounts for biases induced by residual nonvolumetric decorrelation contributions and/or suboptimal vertical wavenumbers assuming the validity of the vertical reflectivity profile(s) F mh (z). The bias (and the inaccuracy) induced by the profile mismatch remain unaccounted for. These can be significant where the assumed vertical reflectivity profile F mh (z) is very different from the (actual) underlying reflectivity.
To compensate for the height bias arising from a systematic profile mismatch across the whole dataset the GEDI RH98 heights are employed. Fig. 6(a) shows on the top the interferometric volume coherence |γ Vol (κ z )| plotted versus the product of the estimated forest height h V with the local vertical wavenumber κ z , i.e., h V κ z at the locations of the GEDI footprints. In Fig. 6(b), the same plot is shown, but this time the interferometric volume coherence |γ Vol (κ z )| is plotted versus the product of the GEDI RH98 height h RH98 with the local vertical wavenumber κ z , i.e., h RH98 κ z . The difference in the two plots [see Fig. 6(a) and (b)] reflects a systematic bias -independent of the local vertical wavenumber-primarily due to the profile mismatch and residual nonvolumetric decorrelation contributions. This bias can be compensated by means of ordinary least square bisector fitting where a 0 and a 1 are the fitting coefficients. Fig. 6(c) shows the interferometric volume coherence |γ Vol (κ z )| plotted versus the product of the corrected [by means of (13)] forest height estimate h V after with the local vertical wavenumber κ z , i.e., h V κ z .

V. DATA PROCESSING AND RESULTS
The proposed processing flow and their main outputs, as described in Sections III and IV, are shown in Fig. 7. For each of the 595 TanDEM-X scenes (CoSSC data products) first the volume decorrelation contributionγ Vol and the terraincorrected vertical wavenumber κ z are estimated, as described in Section III-A. For each scene, the available GEDI waveforms are used to calculate the mean profile following the procedure outlined in Section III-B. The attenuation correction, discussed in Section III-C, is applied to the mean profile to derive a set of reflectivity profiles that are used to estimate forest height from each volume decorrelation sample. For each obtained height estimate, the estimation bias is calculated using the set of profiles and the terrain-corrected vertical wavenumber as described in Section IV-B. A 20% height bias threshold (for the lower heights) and (upper) height threshold associated to the minimum of the |γ Vol (κ z , h V , F mh (z))| / h V derivative is used to derive the validity map to mask out suboptimum performance samples. Similarly, samples with (absolute) coherence lower than 0.3 are discarded (see Section IV-A). In total, about 2% of forest areas are excluded. After processing and inverting each of the TanDEM-X scenes, the forest height estimates along with an estimate of the height bias and the validity masks for the entire island are available and are georeferenced. The next step is to mosaic the individual scene-level maps into 1 degree × 1 degree (latitude, longitude) nonoverlapping tiles. A single location (within a tile) may be covered by multiple TanDEM-X scenes and it may be associated to several forest height estimates. In such a case, and to obtain a single height estimate, either the available estimates can be combined (for example on the basis of their individual estimation biases or by forming their median value) or a single one is selected to fulfill certain vertical wavenumber, acquisition timing or performance selection criteria. Each of these approaches has its own pros and cons and can lead under certain circumstances to outliers. In the following the approach of using the "mean κ z " (select the closest κ z to mean κ z value of all candidates) height value of the valid height estimates at each location is followed. The next step is the global bias correction, which is performed once as described (13) in Section IV. Finally, water and settlement bodies are masked by using the forest / nonforest map. Fig. 8(a) shows the final forest height map at 25-m resolution for the complete island of Tasmania. The heights are validated against the reference CHM. Fig. 8(b) shows the difference between the obtained forest heights and the reference CHM across the whole island. The majority of the samples have an error of ± 5 m with a mean value of -0.08 m. While over large parts of the island, the height difference is very homogenous, locally spatial patterns of over-and underestimation are clearly visible. The reason for these patterns is manifold. Many of these areas can be attributed to the extensive wildfires occurring in   the years between the lidar and the TanDEM-X acquisition, while the underestimation on the western part may be due to the dense forest conditions there. The overestimation of the eastern part is partially due to the forest change and partially due to the mismatch of the derived reflectivity profiles and the real underlying reflectivity. Finally, Fig. 8(c) shows the validation of the estimated forest heights against the CHM accounting for only the valid height estimates with an RMSE of 7.3 m with a Pearson coefficient of 0.66.

VI. PERFORMANCE DISCUSSION AND CONCLUSION
The proposed inversion scheme overcomes the limitations in large-scale forest height estimation caused by the limited dimensionality of TanDEM-X observations with single polarisation by using the waveform and height measurements from GEDI.
Besides the compensation of all nonvolumetric decorrelation contributions [by means of (8)] and the terrain correction of the vertical wavenumber [by means of (5)], the forest height inversion performance critically depends on the ability of the derived set of reflectivity profiles to match the real underlying reflectivity. Even if the effect of the reflectivity profile on the interferometric coherence is, when compared to the effect of the height itself, only secondary, it remains significant enough to have a decisive impact on the achieved performance. A mismatch between the assumed and the underlying reflectivity introduces a (positive or negative) bias on the estimated forest height. In contrast, a well-matching reflectivity profile is able to correctly interpret (4) for a wide range of vertical wavenumbers and provide the "same" height estimates for different vertical wavenumbers-of course, under the performance constraints imposed by the actual vertical wavenumber (as discussed in Section IV). In this sense, to test how well the derived reflectivity profiles match the real underlying reflectivity, forest height estimation performance and its consistency for different vertical wavenumbers are compared. Fig. 9 provides a representative overview of such a performance comparison on the basis of a single TanDEM-X scene. In Fig. 9(a), the reference CHM is shown for an area of 50 × 50 km located in the northern part of the island, whereas in Fig. 9(b), the GEDI RH 98 heights are shown. The forest height map (at 25-m resolution) obtained from a TanDEM-X scene acquired on Jun 19th, 2011 with a (mean) vertical wavenumber of κ z = 0.125 (HoA = 50 m) along a descending orbit is in Fig. 9(c). The validation plot against the CHM is shown in Fig. 10(a) with a performance described by an RMSE of 6.6 m and a Pearson coefficient of 0.63. A second forest height map from a TanDEM-X scene was acquired also along a descending orbit on December 23rd, 2011 with an approximately 40% smaller (mean) vertical wavenumber of κ z = 0.0757 (HoA = 83 m) is shown in Fig. 9(d). The corresponding validation plot against the CHM is shown in Fig. 10(b) with a very similar performance described by an RMSE of 7.2 m and a Pearson coefficient of 0.65. The correlation plot of the two derived height maps is shown in Fig. 11(a) confirming their relatively good agreement.
However, even if the derived set of reflectivity profiles is up to a degree able to describe the mean underlying reflectivity, the single set of profiles used for the entire extent of a TanDEM-X scene is not able to capture the spatial variability of the underlying reflectivity within the scene. For (spatially) uncorrelated mismatch between the assumed and the underlying reflectivity, the effect of the local height biases will appear across the whole scene in form of increased height variance. However, when the mismatch between the assumed and the underlying reflectivity becomes spatially correlated, for example, when the underlying reflectivity changes locally because of the (local) incidence, the height estimates will be (locally) biased. This can be the case in the presence of (positive or negative) terrain slopes or when scenes acquired along ascending and descending orbits are combined. To obtain an impression of the magnitude of this effect forest height obtained from TanDEM-X scenes acquired along ascending and descending orbits are compared.
The forest height map obtained from a TanDEM-X scene acquired on February 5th, 2019 with a (mean) vertical wavenumber of κ z = 0.133 (HoA = 47 m) along an ascending orbit is shown in Fig. 9(e). The validation plot against the CHM is shown in Fig. 10(c) indicating an RMSE of 7.8 m and a correlation 0.59. Finally, in Fig. 9(f), the forest height map obtained from a second TanDEM-X scene acquired along an ascending orbit on January 3rd, 2014 with a (mean) vertical wavenumber of κ z = 0.0938 (HoA = 67 m) that overlaps significantly with the scene shown in Fig. 9(d) (acquired about two years earlier along a descending orbit with a very similar vertical mean wavenumber). The validation plot against the CHM is shown in Fig. 10(d) indicating a similar performance with an RMSE of 7.4 m and a correlation of 0.65. While the agreement between the height estimates from the ascending orbits is similar to the correlation between the height estimates from the descending orbits, the agreement between the ascending and descending orbits is (as expected) somehow lower with an RMSE of about 5.9 and 7.3 m and correlation values of 0.82 and 0.61 [see Fig. 11(b) and (c)]. It is important to remark here that the lidar-derived CHM is by itself associated with a certain uncertainty as well as those different imaging geometries may cause an additional error when comparing the CHM with the radar-derived height maps at sample level.
Finally, the performance on terrain slopes is shown in Fig. 12. In Fig. 12(a), the TanDEM-X 90 m DEM of the site is shown, whereas in Fig. 12(b) and (c), the height validity maps for the two acquisitions on December 23rd, 2011 (descending orbit) and January 3rd, 2014 (ascending orbit). 20% and 13% of the forested areas respectively are masked out because of performance criteria, respectively. For both datasets, most of the positive slope areas are excluded. Most of the height estimates on slopes are on negative slopes allowing us to obtain valid estimates on 98% of the samples when combining both orbits. The common valid, gentle sloped areas, are characterized by a stronger height deviation due to the mismatch of the assumed and the underlying reflectivity profile.
Of course, the proposed combination of GEDI and TanDEM-X data cannot avoid significant absolute height errors, arising from fundamental limitations, such as the limited penetration at X-band. In this case, the interferometric volume coherence no longer represents the full vertical extent of the forest, and, regardless of the assumed vertical reflectivity profile, leads to underestimated forest heights. Accordingly, underestimation of high and/or dense stands because of insufficient X-band penetration cannot be avoided and, more importantly, it is very difficult if not impossible to be detected without additional information (for example, the availability of a DTM). Nevertheless, the achieved overall performance demonstrates the potential of combining lidar and interferometric SAR measurements for large-scale forest structure mapping.