Mapping Tree Height in Burkina Faso Parklands with TanDEM-X

: Mapping of tree height is of great importance for management, planning, and research related to agroforestry parklands in Africa. In this paper, we investigate the potential of spotlight-mode data from the interferometric synthetic aperture radar (InSAR) satellite system TanDEM-X (TDM) for mapping of tree height in Sapon é , Burkina Faso, a test site characterised by a low average canopy cover (~15%) and a mean tree height of 9.0 m. Seven TDM acquisitions from January–April 2018 are used jointly to create high-resolution (~3 m) maps of interferometric phase height and mean canopy elevation, the latter derived using a new, model-based processing approach compensating for some effects of the side-looking geometry of SAR. Compared with phase height, mean canopy elevation provides a more accurate representation of tree height variations, a better tree positioning accuracy, and better tree height estimation performance when assessed using 915 trees inventoried in situ and representing 15 different species/genera. We observe and discuss two bias effects, and we use empirical models to compensate for these effects. The best-performing model using only TDM data provides tree height estimates with a standard error (SE) of 2.8 m (31% of the average height) and a correlation coefﬁcient of 75%. The estimation performance is further improved when TDM height data are combined with in situ measurements; this is a promising result in view of future synergies with other remote sensing techniques or ground measurement-supported monitoring of well-known trees.


Introduction
Parklands are an important agroforestry system widespread across the Sudano-Sahelian zone of West Africa [1,2]. In these landscapes, agriculture and livestock production systems are integrated under a sparse cover of scattered trees. Parklands have been shaped monitoring of tree height, which is a proxy of many other tree properties. In this paper, we first introduce a novel, model-based processing approach that mitigates range offset and ground scattering effects, which are especially significant when studying individual trees with high-resolution TDM data. Using seven spotlight-mode TDM acquisitions with azimuth resolution of 1.1 m, we show that the new processing approach provides a mean canopy elevation estimate that has better potential for tree height mapping than the commonly used phase height, both in terms of top-of-canopy tree height estimation and tree positioning. We also compare the results for 915 trees from 15 different species/genera and we observe and discuss two bias effects in the data: one caused by vegetation bias in the DTM, and one caused by crown shape variability across different species. Finally, we investigate the potential improvement in tree height estimation accuracy when using TDM data with empirical models, both alone and in combination with in situ data of diameter at breast height, crown diameter, and species.
This paper begins with a description of the InSAR measurement method (Section 2) and the available experimental data (Section 3). Thereafter, the results are presented (Section 4) and discussed (Section 5). Finally, conclusions are given in Section 6.

Method
InSAR is an active, microwave remote sensing tool capable of high-resolution elevation measurements, independent of clouds and solar illumination. InSAR systems use the phase difference between two complex-valued images acquired from slightly different positions in space to measure the elevation of objects above a reference surface [22]. TanDEM-X (TDM) consists of two X-band (centre frequency: 9.65 GHz, wavelength: 3.1 cm) SAR satellites flying in a close tandem formation [35]. Because of the short wavelength and the small distance between the satellites (typically below 1000 m), the TDM DEM over vegetated areas often represents canopy height variations. However, there are several challenges associated with the estimation of vegetation height from a TDM DEM.
Firstly, information about ground topography is needed for estimation of vegetation height from an InSAR DEM [36], ideally in the form of a digital terrain model (DTM). Since most global, wall-to-wall DEMs have been acquired with either high-frequency InSAR or optical data [37], they contain vegetation bias and are not useful as ground reference in densely vegetated areas. Airborne lidar scanning (ALS) can be used in combination with InSAR data to provide accurate forest height measurements [26,28,29,32,38], but ALS data are costly to acquire and are not available for most parts of the world. Spaceborne lidar sensors, like the Global Ecosystem Dynamics Investigation (GEDI) mission from NASA [39], are a viable option for global mapping of topography, but the spatial resolution of the acquired DTMs is still relatively low (1 km in the case of GEDI).
Secondly, the elevation measured by an InSAR system is influenced by all objects within a resolution cell, which is slanted due to the side-looking geometry of SAR. Consequently, the measured elevation depends on InSAR system configuration and geometry, as well as the distribution of targets within the resolution cell. In sparsely forested areas, ground and vegetation targets can be observed within the same resolution cell. Additionally, geometric distortions prevail due to the projection of a three-dimensional scene onto a two-dimensional image plane [40].
Thirdly, the difference between the top-of-canopy tree height and the canopy elevation perceived by the radar depends on canopy density, shape, structure, phenology, and moisture, as well as radar polarisation and incidence angle. The relationship between elevation measured with TDM and top-of-canopy height varies in time and with acquisitions. Figure 1 shows the geometry of an InSAR measurement of a single tree with top-ofcanopy tree height h top , average tree crown diameter D avg , and stem diameter at breast height (1.3 m) d bh . Two InSAR quantities used in this paper are also indicated: h pha is the phase height, i.e., the difference between the DEM and the DTM, while h cnp is the mean canopy elevation above the DTM. The slanted geometry of the InSAR measurement is also Remote Sens. 2021, 13, 2747 4 of 25 shown; as a result, ground and canopy objects within the same resolution cell are located at a ground range offset ∆r gr . moisture, as well as radar polarisation and incidence angle. The relationship between elevation measured with TDM and top-of-canopy height varies in time and with acquisitions. Figure 1 shows the geometry of an InSAR measurement of a single tree with top-ofcanopy tree height ℎ , average tree crown diameter , and stem diameter at breast height (1.3 m) . Two InSAR quantities used in this paper are also indicated: ℎ is the phase height, i.e., the difference between the DEM and the DTM, while ℎ is the mean canopy elevation above the DTM. The slanted geometry of the InSAR measurement is also shown; as a result, ground and canopy objects within the same resolution cell are located at a ground range offset Δ . Figure 1. Illustration of some quantities used in this paper. The top-of-canopy tree height ℎ is a purely geometrical quantity. Two InSAR quantities related to tree height are also shown: (1) Phase height (ℎ ) is the estimated difference between a digital elevation model (DEM) and a digital terrain model (DTM). (2) Mean canopy elevation (ℎ ) is a model-based estimate of canopy height above the DTM. The ground range offset Δ is caused by geometrical distortion due to the slanted measurement geometry of SAR.
is the average diameter of the tree crown and is the stem diameter at breast height (1.3 m). Table 1 contains a summary of different metrics (mostly height-related) used throughout this paper, including some that will be introduced later.
In the next section, we define ℎ in terms of the InSAR data. Thereafter, we define the model used to compute ℎ from ℎ and coherence data, as well as the geometric correction Δ . Finally, we discuss the relationship between ℎ , ℎ , and ℎ . Figure 1. Illustration of some quantities used in this paper. The top-of-canopy tree height h top is a purely geometrical quantity. Two InSAR quantities related to tree height are also shown: (1) Phase height (h pha ) is the estimated difference between a digital elevation model (DEM) and a digital terrain model (DTM). (2) Mean canopy elevation (h cnp ) is a model-based estimate of canopy height above the DTM. The ground range offset ∆r gr is caused by geometrical distortion due to the slanted measurement geometry of SAR. D avg is the average diameter of the tree crown and d bh is the stem diameter at breast height (1.3 m). Table 1 contains a summary of different metrics (mostly height-related) used throughout this paper, including some that will be introduced later. Table 1. Summary of some metrics used throughout this paper.

Metric Explanation
In In the next section, we define h pha in terms of the InSAR data. Thereafter, we define the model used to compute h cnp from h pha and coherence data, as well as the geometric correction ∆r gr . Finally, we discuss the relationship between h cnp , h pha , and h top .

Phase Height
The main interferometric quantity is the complex correlation coefficient ( γ), defined as [22]: where s 1 and s 2 are two single-look complex (SLC) images acquired by an InSAR system; γ = | γ| is the coherence, i.e., the magnitude of the complex correlation coefficient; φ is the interferometric phase; φ 0 is a modelled interferometric phase for a reference surface; and E(·) is the expectation value operator. Coherence is a value between 0 and 1 representing the degree of similarity between s 1 and s 2 , while interferometric phase is the phase difference between s 1 and s 2 , corrected for the variability induced by the reference surface. In practical application, the expectation value in (1) is replaced by spatial averaging using a sliding window. The coherence and phase estimated using this approach are affected by wellknown errors [22]. In the case of coherence, coherence overestimation occurs for low numbers of samples and/or low coherence values. In the case of interferometric phase and disregarding 2π phase ambiguities, the estimator is unbiased, but the estimated phase has a zero-mean Gaussian error. Assume a well-designed InSAR system, an acquisition geometry (in terms of incidence angle and distance between the two satellites) providing good sensitivity to typical vegetation heights, negligible temporal change between the two acquisitions, and adequate signal processing including common-band and wavenumber shift filtering. Furthermore, assume that φ 0 in (1) represents the topographic phase modelled from a DTM. Under these assumptions, phase height can be estimated from (1) using: where unw(·) represents unwrapping, i.e., estimation of phase and removal of 2π ambiguities; and κ is the height-to-phase conversion factor, determining the sensitivity of the InSAR system to height variations. κ depends on system properties such as wavelength, separation between the satellites, and incidence angle [41]. κ is related to the height-ofambiguity (HOA), i.e., the height offset corresponding to a 2π shift of the interferometric phase, according to: Note that in (2), unwrapping is done after topographic phase (φ 0 ) removal in (1). It is also possible to estimate h pha by first creating a DEM and then subtracting a DTM. However, because the DEM is expected to have larger height variations than h pha , that approach is more susceptible to phase unwrapping errors.

Mean Canopy Elevation
Volume decorrelation is a loss of coherence caused by the distribution of targets in the direction perpendicular to the image plane. Several models for volume decorrelation in vegetation have been used with TDM data in the past [26,28,29,33,42]. In this paper, we will use the two-level model (TLM) [29,43], where the scattering is assumed to originate from only two levels. It directly separates ground and canopy contributions in InSAR data and does not require multi-polarised data or allometric equations. In the TLM, volume decorrelation is modelled using two discrete levels, ground and vegetation [44]: where ζ is the vegetation scattering fraction and h TLM is the elevation of the vegetation level above the ground level. The vegetation scattering fraction describes the distribution of scattering between the two levels, and it depends on the canopy cover η and ground-tovegetation backscatter ratio ρ as [44]: In this paper, we use the multiplicative coherence model from [22,45] to estimate volume decorrelation using: where γ is the complex correlation coefficient from (1) and γ 0 is an estimate of signal-tonoise decorrelation, i.e., the loss of similarity due to different noise representations in the two images [45]. Note that (4) neglects all decorrelation effects other than volume and signal-to-noise decorrelation. By fitting model (3) to the measured and calibrated data from (4), estimates of ζ and h TLM are obtained.
The h TLM estimated above contains the ground range offset ∆r gr described at the beginning of Section 2 and shown in Figure 1. If E and N are the respective easting and northing coordinates (in metres) of the measured γ and h pha data, then the value of h TLM at these coordinates is assumed to represent the mean canopy elevation (h cnp ) at a position shifted in ground range away from the radar by ∆r gr . For a right-looking system (like TanDEM-X) and disregarding topographic undulations, h cnp can be estimated from h TLM (in metres) through the following interpolation: where θ is the incidence angle and α is the flight heading angle (see Figure 2). The 1/ tan θ factor maps h TLM to the ground range offset ∆r gr , while factors cos α and − sin α project ∆r gr in the east and north directions, respectively.

Estimation of Tree Height
Both ℎ and ℎ provide biased estimates of ℎ due to the penetration of Xband radar signals into the canopy and residual effects of ground topography and inaccuracies in the DTM. Note that ℎ includes the ground scattering contribution, which ℎ aims to mitigate.
In this subsection, we use models to address two potential bias sources: crown shape variations and vegetation bias from the DTM. In Section 3.4 we introduce some empirical models used to compensate for these effects. Other effects, including varying canopy density, moisture, and phenology, are addressed further in Section 5.

Effect of Crown Shape
For the high-resolution InSAR data studied in this paper, tree crowns often occupy

Estimation of Tree Height
Both h cnp and h pha provide biased estimates of h top due to the penetration of X-band radar signals into the canopy and residual effects of ground topography and inaccuracies in the DTM. Note that h pha includes the ground scattering contribution, which h cnp aims to mitigate.
In this subsection, we use models to address two potential bias sources: crown shape variations and vegetation bias from the DTM. In Section 3.4 we introduce some empirical models used to compensate for these effects. Other effects, including varying canopy density, moisture, and phenology, are addressed further in Section 5.

Effect of Crown Shape
For the high-resolution InSAR data studied in this paper, tree crowns often occupy several grid cells. In this paper, we approximate a tree height measurement by taking the topmost pixel within the tree crown extent of the images of h pha and h cnp . The corresponding tree height proxies are referred to as h pha and h cnp , respectively, with the bar symbol indicating the selected topmost pixel.
Trees with the same top-of-canopy height h top , but with differently shaped crowns may have different h cnp values due to different distribution of canopy objects within the topmost pixel (and, indirectly also different h pha values). This is illustrated with a geometric model in Figure 3a. The model assumes ellipsoidal tree crowns with the horizontal semiaxes both equal to D avg 2 and the vertical semiaxis equal to h top 2 . Furthermore, we assume no penetration into the tree crown and that the resolution cell giving h cnp is centred on the tree trunk and has a width δ DEM < D avg . This crown shape bias can then be modelled as the maximal height difference within the resolution cell: with the negative sign indicating underestimation of tree height.

Effect of Vegetation Bias from the DTM
Vegetation bias is the residual effect of vegetation on a digital terrain model (DTM), typically resulting in an overestimation of topographic height in vegetated areas, see Figure 4a. This, in turn, manifests itself as underestimation of tree height if a biased DTM is used as reference. In this work, we use a low-resolution InSAR DEM as a DTM. Overall, this approach provides good results in sparsely forested areas (like the parklands studied here). However, vegetation bias from the DTM may still be noticeable around tall trees with wide canopies.
Vegetation bias from the DTM can be modelled using the TLM described in Section 2.2. We here assume that the average canopy cover within a DTM resolution cell is and the height of all trees is ℎ . Additionally, we assume that the backscattering coefficients for ground and vegetation are equal ( = 1) and that there is no penetration into tree canopies. Under these assumptions, the effect of vegetation bias from the DTM is modelled by the TLM expression (3) to: where the minus sign reflects the fact that a positive vegetation bias in the DTM causes a negative bias in the ℎ (and ℎ ). The resulting bias effect is shown in Figure 4b, for different forest height and canopy  Figure 3b shows expression (6) plotted against D avg /h top , which is an indicator of canopy shape, with lower values indicating narrow canopies and higher values indicating wide canopies. Three values for h top are used and δ DEM = 3 m is assumed. The effect of crown shape bias is largest for shorter trees and/or trees with narrow canopies (low D avg /h top ). For taller trees or trees with wider canopies, the bias is typically less than −1.0 m.

Effect of Vegetation Bias from the DTM
Vegetation bias is the residual effect of vegetation on a digital terrain model (DTM), typically resulting in an overestimation of topographic height in vegetated areas, see Figure 4a. This, in turn, manifests itself as underestimation of tree height if a biased DTM is used as reference. In this work, we use a low-resolution InSAR DEM as a DTM. Overall, this approach provides good results in sparsely forested areas (like the parklands studied here). However, vegetation bias from the DTM may still be noticeable around tall trees with wide canopies.

Test Site
Saponé (12.08° N, 1.57° W) is a rural commune located about 30 km south from Ouagadougou, the capital of Burkina Faso ( Figure 5). The test site is a 10 km × 10 km area dominated by parklands, but also featuring other land cover classes, such as patches of woodlands, small-scale plantations (Mangifera indica L., Eucalyptus camaldulensis Dehnh., and Tectona grandis L.f.), and riparian formations. The terrain is relatively flat, with altitudes varying between 293 and 363 m above sea level [46]. The climate is semi-arid, where the mean annual rainfall is around 800 mm with high inter-annual and inter-seasonal variability. The rainy season generally extends between May and October, with July and August producing the largest proportion of rainfall. The dry season features only sporadic and limited rainfall, and typically starts around November and lasts until May or June.
The tree cover in the parklands is actively shaped by the farmers practicing naturally assisted regeneration, where favoured tree species are protected when preparing the fields before sowing [5]. Crown pruning is also practiced, to revitalise fruit production, limit shading of crops, and provide fodder [47]. Mean tree canopy cover within the test site is around 15 percent. The tree cover is dominated by traditional agroforestry species, including the native V. paradoxa, P. biglobosa, and Lannea microcarpa Engl. and K. Krause, as well as M. indica [48], native to the Indian subcontinent but often cultivated in the area. Although these species are generally considered deciduous, they are seldom leafless due to a progressive leaf replacement throughout the year [49]. Notable exceptions include L. microcarpa, which loses all its leaves early in the dry season, and F. albida, which has reverse phenology and is foliated only during the dry season. Vegetation bias from the DTM can be modelled using the TLM described in Section 2.2. We here assume that the average canopy cover within a DTM resolution cell is η and the height of all trees is h TLM . Additionally, we assume that the backscattering coefficients for ground and vegetation are equal (ρ = 1) and that there is no penetration into tree canopies. Under these assumptions, the effect of vegetation bias from the DTM is modelled by the TLM expression (3) to: where the minus sign reflects the fact that a positive vegetation bias in the DTM causes a negative bias in the h cnp (and h pha ). The resulting bias effect is shown in Figure 4b, for different forest height and canopy cover values, and assuming κ = 0.063 (HOA ≈ 100 m). For the studied forest area, the average tree height is about 9 m and the average canopy cover at 1 ha scale is typically 15%, so the modelled effect of vegetation bias from the DTM is −1.3 m. For some 1 ha areas centred around tall trees with wide canopies, the average height at 1 ha-scale can reach 16 m and canopy cover can reach 25%; in that case, effect of vegetation bias in the DTM can cause a bias about −3.7 m. Furthermore, if penetration into vegetation canopies and crown shapes are also considered, the vegetation bias will be less significant. It is thus expected that only for some of the largest trees, the effect of vegetation bias will be noticeable.

Test Site
Saponé (12.08 • N, 1.57 • W) is a rural commune located about 30 km south from Ouagadougou, the capital of Burkina Faso ( Figure 5). The test site is a 10 km × 10 km area dominated by parklands, but also featuring other land cover classes, such as patches of woodlands, small-scale plantations (Mangifera indica L., Eucalyptus camaldulensis Dehnh., and Tectona grandis L.f.), and riparian formations. The terrain is relatively flat, with altitudes varying between 293 and 363 m above sea level [46]. The climate is semi-arid, where the mean annual rainfall is around 800 mm with high inter-annual and inter-seasonal variability. The rainy season generally extends between May and October, with July and August producing the largest proportion of rainfall. The dry season features only sporadic and limited rainfall, and typically starts around November and lasts until May or June.

Reference Data
Two orthorectified, high-resolution optical satellite images over Saponé were used in this study: one image acquired in December 2012 with the WorldView-2 satellite and one image acquired in October 2017 with a Pléiades satellite. Both images had a ground sampling distance of about 50 cm.
For initial geocoding of the TDM spotlight-mode data, we used the freely available, global 90 m TDM DEM [25]. We did not use the available higher-resolution TDM DEMs because we wanted to reduce vegetation bias, and the resolution of 90 m was found sufficient for initial geocoding in this relatively flat test site.
The in situ tree inventory data included three datasets collected in 2012, 2017 and 2018 (Table 1). For each dataset, trees with stem diameter at breast height ( , measured at 1.3 m) ≥ 5 cm were georeferenced using a handheld Garmin Oregon 550 GPS and tree species were recorded. Top-of-canopy tree height (ℎ ) was measured using a Haglöf electronic clinometer from a distance between 10-20 m, depending on the line of sight. The average crown diameter ( ) was determined by averaging two perpendicular crown diameter measurements. Positional uncertainties related to GPS accuracy were accounted for by manually matching easily identifiable trees with the high-resolution optical imagery using information on crown diameter, height, and species.
The in situ dataset from 2012 was collected between October and December within 76 plots (50 m × 50 m in dimensions), randomly distributed throughout the 100 km 2 test site and equally divided between three canopy cover classes, derived using the WorldView-2 image [15,20,48]. This resulted in a total of 1125 measured trees.
The in situ dataset from 2017 was collected in October using a sampling approach where plots were randomly placed in active parkland fields identified in the Pléiades image acquired two weeks before the inventory. This dataset consisted of 637 trees distributed over the entire 100 km 2 test site.
The in situ dataset from 2018 was acquired in June. Three large plots (about 500 m× 100 m in dimensions) were laid out in the central part of the test site, and a total of 321 trees were measured within these three plots. The three plots were chosen to cover areas with well-separated trees.
The three datasets were subsequently filtered to only include trees covered by all TDM acquisitions (see Section 3.3 and Figure 5). Trees that had undergone visible change The tree cover in the parklands is actively shaped by the farmers practicing naturally assisted regeneration, where favoured tree species are protected when preparing the fields before sowing [5]. Crown pruning is also practiced, to revitalise fruit production, limit shading of crops, and provide fodder [47]. Mean tree canopy cover within the test site is around 15 percent. The tree cover is dominated by traditional agroforestry species, including the native V. paradoxa, P. biglobosa, and Lannea microcarpa Engl. and K. Krause, as well as M. indica [48], native to the Indian subcontinent but often cultivated in the area. Although these species are generally considered deciduous, they are seldom leafless due to a progressive leaf replacement throughout the year [49]. Notable exceptions include L. microcarpa, which loses all its leaves early in the dry season, and F. albida, which has reverse phenology and is foliated only during the dry season.

Reference Data
Two orthorectified, high-resolution optical satellite images over Saponé were used in this study: one image acquired in December 2012 with the WorldView-2 satellite and one image acquired in October 2017 with a Pléiades satellite. Both images had a ground sampling distance of about 50 cm.
For initial geocoding of the TDM spotlight-mode data, we used the freely available, global 90 m TDM DEM [25]. We did not use the available higher-resolution TDM DEMs because we wanted to reduce vegetation bias, and the resolution of 90 m was found sufficient for initial geocoding in this relatively flat test site.
The in situ tree inventory data included three datasets collected in 2012, 2017 and 2018 (Table 1). For each dataset, trees with stem diameter at breast height (d bh , measured at 1.3 m) ≥ 5 cm were georeferenced using a handheld Garmin Oregon 550 GPS and tree species were recorded. Top-of-canopy tree height (h top ) was measured using a Haglöf electronic clinometer from a distance between 10-20 m, depending on the line of sight. The average crown diameter (D avg ) was determined by averaging two perpendicular crown diameter measurements. Positional uncertainties related to GPS accuracy were accounted for by manually matching easily identifiable trees with the high-resolution optical imagery using information on crown diameter, height, and species.
The in situ dataset from 2012 was collected between October and December within 76 plots (50 m × 50 m in dimensions), randomly distributed throughout the 100 km 2 test site and equally divided between three canopy cover classes, derived using the WorldView-2 image [15,20,48]. This resulted in a total of 1125 measured trees.
The in situ dataset from 2017 was collected in October using a sampling approach where plots were randomly placed in active parkland fields identified in the Pléiades image acquired two weeks before the inventory. This dataset consisted of 637 trees distributed over the entire 100 km 2 test site.
The in situ dataset from 2018 was acquired in June. Three large plots (about 500 m× 100 m in dimensions) were laid out in the central part of the test site, and a total of 321 trees were measured within these three plots. The three plots were chosen to cover areas with well-separated trees.
The three datasets were subsequently filtered to only include trees covered by all TDM acquisitions (see Section 3.3 and Figure 5). Trees that had undergone visible change (e.g., removal or significant pruning) between the December 2012 WorldView-2 image and the October 2017 Pléiades image were also excluded. The remaining trees were then sorted by species (or genus, if species could not be determined) and species/genera with less than five trees were excluded from the dataset. This resulted in a total of 915 trees left for this analysis, representing 15 different species/genera, see Table 2. Table 2. Summary of the in situ data from Saponé, Burkina Faso, used in this study. The second column from the left contains the number of trees fulfilling the conditions described in Section 3.2, as well as the total number of trees sampled in field.

Dates
Used Note that although the TDM data were acquired in early 2018, the temporal offsets between all three in situ datasets and the TDM data were ignored in this study. In particular, this concerns the dataset from 2012, which was included because it contained the largest number of in situ measured trees, allowing a more reliable statistical analysis and a more extensive study of the effect of tree species; also, it featured trees with the lowest average tree heights and crown diameters, thus improving the sampled interval of tree heights and crown diameters, see Table 2. Additionally, the growth rate for most of the tree species is low in the relevant conditions and our empirical investigations showed no noticeable effect of temporal difference between the three datasets.

TanDEM-X Data Processing
This study used a total of seven HH-polarised TDM spotlight-mode acquisitions made between January and April 2018. Three acquisitions were made from the descending orbit (flight heading around 191 • relative to true north), at an incidence angle of 32 degrees (at scene centre), with HOA values between 39 and 55 metres, and with approximate ground range and azimuth resolutions at scene centre of 2.1 and 1.1 metres, respectively. The remaining four acquisitions were made from the ascending orbit (flight heading about 349 • ), at an incidence angle of about 25 degrees, with HOAs between 49 and 87 metres, and with ground range and azimuth resolutions of 2.8 and 1.1 metres, respectively. Table 3 contains a summary of the acquisition parameters, together with air temperatures recorded at Ouagadougou Airport at the time of acquisition [50]. All acquisitions were made during the dry season and no precipitation was recorded at Ouagadougou Airport during the 72 h prior to any of the seven acquisitions [50]. Figure 5 shows the approximate outlines of the ascending and descending acquisitions. The joint area covered by all seven acquisitions was 4200 hectares. Table 3. Summary of the SAR acquisitions over Saponé, Burkina Faso used in this study. "No" refers to the relative orbit number, "Dir" refers to orbit direction ("dsc" is descending and "asc" is ascending), α and θ are the flight heading and incidence angles, respectively (see Figure 2), "Pol" refers to polarisation, "Res" refers to resolution ("grg" is ground range and "az" is azimuth), and "Temp" refers to temperature.

Date
Time ( All acquisitions from the same orbit were co-registered to one single master image using the GAMMA software package [51]. The remaining processing was then conducted using Python scripts based on [52]. Fine-resolution interferograms (images of γ) were formed using a sliding 1 × 2 window (range × azimuth), giving an approximate azimuth resolution of 2.2 m and ground range resolutions of 2.1 m for the descending data and 2.8 m for the ascending data. For coherence, a 5 × 5 window was used to reduce bias in coherence estimation [22].
Next, orbit state vectors and the free global 90 m TDM DEM [25] were used to create simulated DEM phase images, which were then subtracted from the fine-resolution interferograms, providing flattened interferograms. These were subsequently filtered using a sliding 45 × 45 averaging window, producing flattened interferograms with an approximate resolution of 100 m. A digital terrain model (DTM) was created from these coarse interferograms by unwrapping and scaling to height, then averaging across all seven acquisitions, and finally adding back to the 90 m TDM DEM. Phase height (h pha ) images were created by first subtracting the coarse-resolution flattened interferograms from the fine-resolution flattened interferograms, and then unwrapping and scaling to height.
The seven phase height images obtained were then geocoded and averaged into two images: h asc pha for the ascending orbit and h dsc pha for the descending orbit. Then, spatial cross-correlation was used to match h asc pha and h dsc pha to each other, and the final phase height image h pha was created by taking the maximum phase height value from h asc pha and h dsc pha , individually for each pixel. Coherence was calibrated using (4), where signal-to-noise decorrelation γ 0 was estimated using the SNR decorrelation model from [35,45] and the noise model provided with the TDM data [53].
TLM fitting was carried out using the principles of active surface modelling, by minimising the following cost function based on [54,55] with respect to h TLM and ζ: where w 1 and w 2 are inversion parameters; h E and h N are the first order spatial derivatives of h TLM in the east and north directions, respectively; h EE and h NN are the respective second order spatial derivatives; γ i is the volume decorrelation for acquisition i estimated using (4); κ i is the corresponding height-to-phase scaling factor; |·| is the magnitude operator; and N is the total number of acquisitions in each geometry (three for descending, four for ascending). A multi-temporal TLM inversion approach was used because it is less susceptible to phase unwrapping errors when used with data with different HOA values [44]. The minimisation of (8) was carried out using a gradient descent algorithm, with the smoothing parameters w 1 and w 2 determined empirically to provide smoothed height estimates and sharp canopy edges when evaluated against optical imagery. This procedure yielded two height images, h asc TLM and h dsc TLM , one for each flight heading. Each pixel was subsequently interpolated using (5), separately for h asc TLM and h dsc TLM . In case multiple grid cells were shifted to the same position, the maximal height value was selected. This procedure resulted in two images of mean canopy elevation: h asc cnp and h dsc cnp , one for each orbit direction, which were subsequently merged into the final mean canopy elevation image (h cnp ), by minimizing the following cost function based on [54,55] with respect to h cnp : where w 1 and w 2 are inversion parameters; h E and h N are the first order spatial derivatives of h cnp , in the east and north directions, respectively; and h EE and h NN are the respective second order spatial derivatives. The minimization was carried out in the same way as for (8).

Estimation of Tree Height from Phase Height and Mean Canopy Elevation
Tree-level estimates of h pha and h cnp were extracted for all reference trees described in Section 3.2 and Table 2 as the maximal pixel values found within the visible parts of the canopies: where T is the set of grid cells located within the extents of the tree crown. A simple geometric model was used to determine T: all trees were assumed to have circular crowns with crown diameter D avg and, in case of overlapping tree crowns, the taller trees were assumed to obscure the shorter trees within the area of crown overlap. Tree position was estimated from the spatial position of the maximal pixel value within the tree crown.
Due to the complex way in which h pha and h cnp depend on tree properties (see Section 2.3), we used empirical models to estimate h top from these quantities. All models were of the following form: where 1 and 2 are zero-mean Gaussian errors, f is one of the two data transform functions: and h is one of the two calibrated TDM-based proxies of tree height: (i) Calibrated phase height: (ii) Calibrated mean canopy elevation: where c 0 and c 1 are calibration constants estimated from the data. We investigated all possible variants of (10) with between one and seven parameters (i.e., up to six of the seven parameters p i set to zero), excluding the trivial case of a constant model ( f h top = p 0 + 2 ). For each of the models, model parameters were estimated twice: once jointly for all tree species/genera, and once individually for each species/genus. This resulted in a total of 760 models. Depending on the type of data used, these 760 models consisted of 16 models using TDM data only, 248 with in situ data only, and 496 models that combined in situ and TDM data.
Calibration constants c 0 and c 1 were estimated by fitting a linear model of h pha or h cnp to in situ measured data for h top by means of orthogonal distance regression (ODR), which is a regression technique assuming that some of the independent variables are also affected by a measurement error [56]. For models that included f (h ) (i.e., with non-zero p 1 ), we used ODR to estimate the parameters by fitting to the in situ measured data and assuming that f (h ) is the only independent variable with an associated uncertainty. Meanwhile, for models with p 1 = 0, ordinary least squares (OLS) regression was used. In the case when model parameters were fitted individually for each species, the number of trees within each species/genus was required to be at least ten times higher than the number of parameters to be estimated. For that reason, the number of used trees and species/genera was reduced than for the species-specific models.
The estimated model parameters, the inverse transforms of f , and (10) were used to predict an initial top-of-canopy heightĥ top . This estimate was subsequently corrected for bias using:ĥ This correction aimed to remove both the logarithmic bias occurring when using f (x) = ln(x) and the bias occurring when zero-intercept (p 0 = 0) models were used.

Results
This section consists of five parts. First, we study samples of the produced maps of h pha and h cnp , and assess their potential for mapping canopy height variations (Section 4.1). Then, in Section 4.2, we assess tree positioning accuracy. Subsequently, we study the correlation between TDM proxies of top-of-canopy height (h pha and h cnp ) and the in situ measured h top ; first, we do this for all 915 trees from 15 species/genera (Section 4.3), and then individually for each species/genus (Section 4.4). Finally, we use empirical models to evaluate the potential of TDM data for tree height estimation, with and without supporting in situ measurements (Section 4.5). Figure 6 shows seven different images for the same 20-hectare area in Saponé. The area features a typical parkland environment, with trees of various heights and crown diameters scattered among agricultural fields (Figure 6a). Backscatter intensity images for the descending and ascending orbit directions illustrate effects of the side-looking geometry of SAR (Figure 6b,c). Shadows occur in areas obscured by the trees, while the parts of tree crowns facing the radar are brighter than the rest of the image, indicating enhanced vegetation scattering. These effects are especially prominent for the row of large trees visible in the central part of the studied area (area A in Figure 6). Varying degrees of distortion and shadowing between the ascending and descending data is due to the different incidence angles for the two orbit directions (32 degrees for descending, 25 degrees for ascending). Additionally, the oblique side-looking geometry causes a tree height-dependent ground range offset. This becomes evident when images of h pha from descending and ascending directions are compared (Figure 6d,e): for large, tall trees, the difference between ascending and descending is larger than for smaller trees. Therefore, the combination of ascending and descending phase height images using a constant (space-invariant) offset produces an image that is blurred for some of the trees. This effect is clearly visible for areas B and C, where some trees are well focussed while others are blurred in Figure 6f, indicating different range offsets. The model-based approach proposed in Section 2.2 provides an image of h cnp , which is better focussed across different tree heights (Figure 6g). Again, this becomes obvious for areas B and C, where all trees are better focussed than in Figure 6f. Moreover, for the image of h cnp in Figure 6g, the outlines of the tree canopies are better matched with the reference satellite photo in Figure 6a, compared with the image shown in Figure 6f.  Table 4 evaluates tree positioning accuracy when using maps of ℎ and ℎ . The offset between tree positions from TDM and in situ data is quantified in terms of its mean  Table 4 evaluates tree positioning accuracy when using maps of h pha and h cnp . The offset between tree positions from TDM and in situ data is quantified in terms of its mean value and standard deviation in the east and north directions. Three tree height groups are used: short (up to 8 m), medium (8-16 m), and tall (above 16 m). For these three groups, the respective tree counts are: 418, 436, and 61. The results are provided separately for ascending and descending data, and for the final, combined images. Table 4. Tree positioning performance for phase height and mean canopy elevation images, individually for ascending and descending data and for the final, combined images. The results are provided separately for three tree height groups: short (below 8 m), medium (8-16 m), and tall (above 16 m), as well as for all trees. The biases (mean offsets) between measured and reference tree positions in the east and north directions are denoted with µ E and µ N , respectively, while the corresponding standard deviations are σ E and σ N . The highlighted values are discussed in Section 4.2. The tree height-dependent range offset primarily manifests itself as a bias (mean offset) in the eastern direction observed in phase height images. For tall trees, the treetops are shifted on average 6.0 m in the eastern direction for the ascending orbit and −3.9 m for the descending orbit. The bias is lower for shorter trees because the co-registration routine uses spatial cross-correlation, which matches the two images using the more abundant short trees. By combining the ascending and descending data, the average bias is decreased, but the uncertainty in location is large, with the standard deviation for the east direction being 7.3 m. The corresponding values for the mean canopy elevation images are better: the tall trees are shifted in the east direction by about −0.6 m for the ascending orbit and 2.8 m for the descending orbit, and the standard deviation for height positioning in the east direction is 3.5 m. In the north direction, the bias and standard error are typically much lower and similar for both phase height and mean canopy elevation. Figure 7 shows scatterplots comparing h pha and h cnp to h top for all 915 trees from 15 species/genera used in this study. Bias is indicated with a red line obtained by fitting a linear function to the data using orthogonal distance regression. The corresponding bias equation is also given, together with the coefficient of determination (R 2 ) for the regression. The respective estimation statistics are given in Table 5. ℎ ⋆ , while = −1.88 was used with (12) to obtain the calibrated mean canopy height estimate ℎ ⋆ .  Table 5. and ℎ (with a Pearson correlation coefficient, , between 66% and 75%) and the standard error (SE) is 20-35% of the average ℎ . For the four less abundant species (Azadirachta indica A. Juss., Entada africana Guill. and Perr., Tectona grandis, and Eucalyptus camaldulensis), the correlation between ℎ and ℎ is low (4-39%), and the relative SE values are higher (32-49%). Note that these groups contain some of the smallest trees in this study: Entada africana and Tectona grandis have the lowest average and ℎ , while Eucalyptus camaldulensis and A. indica have some of the lowest average . Finally, each of the remaining seven least abundant species/genera (Ficus spp., Sclerocarya birrea (A. Rich.) Hochst., Bombax costatum Pellegr. and Vuillet., Acacia spp., Diospyros mespiliformis Hochst. ex A. DC., Terminalia laxiflora Engl. and Diels, and F. albida) shows good correlation between ℎ and ℎ (65-96%), and a relative SE of 16-41%.  Table 5. Table 5. Tree height estimation statistics for h pha and h cnp . " Figure" refers to the figure with the corresponding scatterplot, N is the number of available tree measurements, and S is the number of species/genera represented in the available tree measurements. Bias and standard error metrics are given in three tree height categories: short trees (below 8 m), medium trees (8-16 m), and tall trees (above 16 m), as well as for all trees.   The results indicate that h pha is a biased estimator of tree height, and the bias varies with tree height. On average, for shorter trees, h pha underestimates the tree height with over 2 m, and the underestimation increases to about 5 m for taller trees, as indicated by the red line in Figure 7a. For h cnp , the bias is on average lower, and the underestimation is consistently below 2 m, as indicated by the red line in Figure 7b. Note that for some trees, the estimated phase height is negative.

Effect of Species/Genera on Mean Canopy Elevation
The intercept values shown in Figure 7 were subsequently used for calibration of phase height and mean canopy elevation, needed for the empirical models studied in Section 4.5: c 0 = −2.26 was used with (11) to obtain the calibrated phase height estimate h pha , while c 1 = −1.88 was used with (12) to obtain the calibrated mean canopy height estimate h cnp .

/ℎ
indicates a short and wide crown. Marker sizes are proportional to average ℎ and error bars are also shown, indicating the 25th and 75th percentiles of bias and /ℎ .
There is a correlation between bias and /ℎ . A larger underestimation of ℎ observed for species/genera with low /ℎ , i.e., with narrow crowns. One notable exception from the general trend is P. biglobosa, for which a negative bias is measured, compared with the positive bias that would be expected from the overall trend. For the four most abundant species/genera (V. paradoxa, Lannea spp., M. indica, and P. biglobosa), there is a clear correlation between h cnp and h top (with a Pearson correlation coefficient, r P , between 66% and 75%) and the standard error (SE) is 20-35% of the average h top . For the four less abundant species (Azadirachta indica A. Juss., Entada africana Guill. and Perr., Tectona grandis, and Eucalyptus camaldulensis), the correlation between h cnp and h top is low (4-39%), and the relative SE values are higher (32-49%). Note that these groups contain some of the smallest trees in this study: Entada africana and Tectona grandis have the lowest average D avg and h top , while Eucalyptus camaldulensis and A. indica have some of the lowest average D avg . Finally, each of the remaining seven least abundant species/genera (Ficus spp., Sclerocarya birrea (A. Rich.) Hochst., Bombax costatum Pellegr. and Vuillet., Acacia spp., Diospyros mespiliformis Hochst. ex A. DC., Terminalia laxiflora Engl. and Diels, and The dependence of bias on h top and D avg is investigated further in Figure 9, where the bias shown in Figure 8 is plotted against the D avg /h top ratio, which is an indicator of crown shape. A low D avg /h top indicates a tall and narrow tree crown, while a high D avg /h top indicates a short and wide crown. Marker sizes are proportional to average h top and error bars are also shown, indicating the 25th and 75th percentiles of bias and D avg /h top .  Figure 8. Dependence of the estimated mean canopy elevation (ℎ ) on the in situ-measured tree height (ℎ ) and species/genera for the 915 trees shown in Figure 6. ⟨ ⟩ and ⟨ℎ ⟩ are the average crown diameter and tree height for each species while is the number of trees. Marker sizes are proportional to .

Tree Height Estimation with Empirical Models
The potential of TDM height proxies and in situ data for tree height estimation was also investigated using empirical models. Although we investigated a total of 760 models, only a small subset of the best-performing models is presented here, while the reader is referred to Supplementary Materials for a comprehensive compilation of the results for all tested models. This evaluation serves three purposes: (i) models using only TDM data assess the potential of two different techniques for TDM-based tree height measurement, (ii) models using in situ data only provide performance metrics that can be used as benchmark when assessing the remote sensing-based methods, and (iii) models combining TDM proxies with in situ data provide an interesting operational alternative for long-term monitoring of trees with some, easy-to-measure quantities accessed from the ground and height monitored from satellite. Moreover, these models also assess the potential synergies of TDM-based height estimation with future techniques capable of providing reliable estimates of, e.g., crown diameter and species.
Overall, it was found that ℎ ⋆ provided better tree height estimation results than ℎ ⋆ , and logarithmic models typically performed better than linear. Of the two fieldmeasured metrics and , the latter provided better estimation results, both when used alone and in combination with the TDM-based estimates. Species-specific models provided significantly better results than species-independent models. Figure 10 shows scatterplots for six selected models, while Table 6 contains mathematical expressions and performance metrics for the respective models. Each of the six models was the best-performing in terms of standard error (SE, i.e., the standard deviation of residuals) out of all tested models with that input data. In case multiple models gave the same SE (to the first decimal), the model with the fewest parameters was selected. There is a correlation between bias and D avg /h top . A larger underestimation of h top observed for species/genera with low D avg /h top , i.e., with narrow crowns. One notable exception from the general trend is P. biglobosa, for which a negative bias is measured, compared with the positive bias that would be expected from the overall trend.

Tree Height Estimation with Empirical Models
The potential of TDM height proxies and in situ data for tree height estimation was also investigated using empirical models. Although we investigated a total of 760 models, only a small subset of the best-performing models is presented here, while the reader is referred to Supplementary Materials for a comprehensive compilation of the results for all tested models. This evaluation serves three purposes: (i) models using only TDM data assess the potential of two different techniques for TDM-based tree height measurement, (ii) models using in situ data only provide performance metrics that can be used as benchmark when assessing the remote sensing-based methods, and (iii) models combining TDM proxies with in situ data provide an interesting operational alternative for long-term monitoring of trees with some, easy-to-measure quantities accessed from the ground and height monitored from satellite. Moreover, these models also assess the potential synergies of TDM-based height estimation with future techniques capable of providing reliable estimates of, e.g., crown diameter and species.
Overall, it was found that h cnp provided better tree height estimation results than h pha , and logarithmic models typically performed better than linear. Of the two field-measured metrics D avg and d bh , the latter provided better estimation results, both when used alone and in combination with the TDM-based estimates. Species-specific models provided significantly better results than species-independent models. Figure 10 shows scatterplots for six selected models, while Table 6 contains mathematical expressions and performance metrics for the respective models. Each of the six models was the best-performing in terms of standard error (SE, i.e., the standard deviation of residuals) out of all tested models with that input data. In case multiple models gave the same SE (to the first decimal), the model with the fewest parameters was selected.
Using TDM data alone in an empirical model results in a correlation of 75%, an SE of 2.8 m and the error is similar for all three height groups (Figure 10a). For comparison, in situ-based measurements of average crown diameter (D avg ) and diameter at breast height (d bh ) give somewhat better overall results, although the performance varies more strongly with height: D avg seems to perform better for tall trees, while d bh is a better estimate of height for shorter trees (Figure 10b,c). Combination of D avg and d bh gives a correlation of 82% and an SE of 2.3 m (Figure 10d). Using species-specific models with TDM data clearly helps to mitigate some species-specific bias effects, improving the correlation to 79% and the SE to 2.6 m (Figure 10e). Finally, for comparison, the best results obtained using in situ data only with species-specific models is a correlation of 87% and an SE of 2.0 m (Figure 10f). The cases where TDM data are combined with in situ-measurements are also of interest. The results provided in Supplementary Materials indicate that h cnp together with d bh in the best-performing empirical model gives an SE of 2.3 m and a correlation of 82% if species-independent models are used, and 2.1 m and 85%, respectively, for species-specific models. Using TDM data alone in an empirical model results in a correlation of 75%, an SE of 2.8 m and the error is similar for all three height groups (Figure 10a). For comparison, in situ-based measurements of average crown diameter ( ) and diameter at breast height ( ) give somewhat better overall results, although the performance varies more strongly with height: seems to perform better for tall trees, while is a better estimate of height for shorter trees (Figure 10b,c). Combination of and gives a correlation of 82% and an SE of 2.3 m (Figure 10d). Using species-specific models with TDM data clearly helps to mitigate some species-specific bias effects, improving the correlation to 79% and the SE to 2.6 m (Figure 10e). Finally, for comparison, the best results obtained using in situ data only with species-specific models is a correlation of 87% and an SE of 2.0 m (Figure 10f). The cases where TDM data are combined with in situ-measurements are also of interest. The results provided in Supplementary Materials indicate that ℎ together with in the best-performing empirical model gives an SE of 2.3 m and a correlation of 82% if species-independent models are used, and 2.1 m and 85%, respectively, for species-specific models.  Table 6. For (a-d), the models were fitted to data from all 915 trees, disregarding the species/genus of the trees. For (e,f), the models were fitted individually for each species/genus with at least ten times more trees than model parameters. In each panel, the first scatterplot from the left shows the estimated tree height (ℎ ) on the -axis against the reference tree height (ℎ ) on the -axis, while the second scatterplot shows the obtained tree height residual (ℎ − ℎ ) against average canopy diameter ( ).

Tree Height Estimation Performance
This paper assessed the potential of spotlight-mode, interferometric TanDEM-X (TDM) data for mapping of tree height in the parklands of Burkina Faso. Two approaches were compared: one using phase height (ℎ ), i.e., the elevation of an InSAR digital elevation model (DEM) above a digital terrain model (DTM); and another using mean canopy elevation (ℎ ) derived using a novel, model-based processing approach correcting for the side-looking geometry of SAR. The latter, more complex processing approach provided a better geometric representation of canopy height variations, better tree positioning accuracy, and better tree height estimation performance, with a standard error (SE) of  Table 6. For (a-d), the models were fitted to data from all 915 trees, disregarding the species/genus of the trees. For (e,f), the models were fitted individually for each species/genus with at least ten times more trees than model parameters. In each panel, the first scatterplot from the left shows the estimated tree height (ĥ top ) on the y -axis against the reference tree height (h top ) on the x -axis, while the second scatterplot shows the obtained tree height residual (h top −ĥ top ) against average canopy diameter (D avg ). Table 6. Tree height estimation statistics for selected empirical models of TDM data and in situ-measurements. " Figure" refers to the figure with the corresponding scatterplot, P is the total number of estimated parameters, N is the number of available tree measurements, and S is the number of species/genera represented in the available tree measurements. For the species-specific models, the ratio between the number of trees within each species/genus and the number of model parameters to-be-estimated was required to be at least 10, thus reducing the number of included species and the total number of trees. Bias and standard error metrics are given in three tree height categories: short trees (below 8 m), medium trees (8-16 m), and tall trees (above 16 m), as well as for all trees. This table is an excerpt from the full results that can be found in Supplementary Materials.

Tree Height Estimation Performance
This paper assessed the potential of spotlight-mode, interferometric TanDEM-X (TDM) data for mapping of tree height in the parklands of Burkina Faso. Two approaches were compared: one using phase height (h pha ), i.e., the elevation of an InSAR digital elevation model (DEM) above a digital terrain model (DTM); and another using mean canopy elevation (h cnp ) derived using a novel, model-based processing approach correcting for the side-looking geometry of SAR. The latter, more complex processing approach provided a better geometric representation of canopy height variations, better tree positioning accuracy, and better tree height estimation performance, with a standard error (SE) of 2.8 m (31% of the average tree height of 9.0 m) and a small overall bias for most trees.
To the authors' current knowledge, no studies so far have evaluated satellite-based measurements of individual tree height in parkland areas, while only a few studies have evaluated satellite-based estimation of individual tree height or average tree height within small plots or sparsely forested areas. [31] used TDM data to estimate tree height in northwestern Canada. A mean absolute error (MAE) of 0.72 m was reported for 4185 trees with an average reference height from ALS data of 2.47 m. In our study, the best model using only TDM data and all 915 trees from 15 species/genera would give a MAE of 2.3 m for an average h top of 9.0 m. This would translate to a relative MAE of 25%, while for [31], the corresponding value would be 29%. [57] measured individual tree height in lichen woodlands in Canadian subarctic using WorldView-3 stereo-photogrammetry. A root-mean-square error (RMSE) of 1.27 m was reported for 96 trees with heights in the interval 2-12 m. [58] used full-waveform data from the ICESat GLAS spaceborne laser scanning system to estimate average tree height within 23 1.5-hectare footprints in a savanna landscape in Kruger National Park, South Africa. The best models used in the study provided an RMSE of 2.42 m for an average tree height range of 5-22 m. [59] used repeat-pass InSAR data from Sentinel-1 and RADARSAT-2 to estimate average tree height for 11 0.1-hectare plots with average canopy height of 12.7 m. The best obtained MAE and RMSE were 1.30 m and 1.34 m, respectively.
In this study, we showed that improvement of height estimation could be obtained by combining TDM measurements with selected in situ data. In particular, the use of d bh proved advantageous from the point of view of tree height estimation: while using d bh alone provided estimates with SE of 2.5 m, combination of d bh and h cnp further improved the estimation performance to an SE of 2.3 m. This observation has a practical implication: the measurement of d bh is easy to conduct with simple tools (e.g., through circumference measurement with a measuring tape) and is not affected by canopy pruning or season, so it is expected to be more stable over time. Meanwhile, D avg is more difficult to measure and is affected by pruning, moisture, and phenology. However, while d bh is difficult to measure with remote sensing methods, high-resolution optical satellite images can be used to estimate D avg [15]. For that reason, models combining TDM-based tree height metrics with D avg are also of interest for future applications.
Furthermore, if species information is also available, then species-specific models can be derived, giving an improvement of TDM-based estimation performance and an SE of 2.6 m. Although tree species determination from satellite data is a notoriously difficult task, especially in areas where species diversity is high and the geographical extent is large [60], species do not change over time which is useful for monitoring of existing trees or plantations. The development in spatiotemporal and spectral resolution of recent satellite systems and improvements in image classification methods may pave the way for accurate tree species mapping in the near future [61].
The results obtained in this study and in [31] show that TDM has good potential for mapping and monitoring of height for individual trees, in particular in remote and/or frequently cloud-covered areas, where other measurement methods are ineffective. In this study, to get a sufficiently large dataset of tree height measurements, we used in situ data acquired up to 6 years prior to the TDM measurements. Due to the lack of suitable information on growth and pruning activities, we neglected temporal changes occurring between the in situ and TDM measurements. The unaccounted temporal changes have certainly hampered the observed estimation performance.
The analysis revealed that tree height estimation performance varies across species/genera ( Figure 8). The observed underestimation was largest for tree species/genera with tall and narrow crowns, while most species with wide crown showed less bias (Figure 9). A notable exception was P. biglobosa, which showed a relatively large underestimation of tree height, despite being the tallest tree species in this comparison, with the largest tree crowns. Two potential explanations for these effects are (1) crown shape bias, caused by varying distribution of canopy objects within the topmost pixel and most prominent for trees with narrow canopies, and (2) vegetation bias from the DTM, most prominent for tall trees with wide canopies. In this paper, the observed systematic errors could be reduced with empirical models (Figure 10), but better understanding of the systematic errors is key for future large-scale use of the methods presented in this paper.
In this study, we did not observe any clear dependence of tree height estimation bias on canopy density, moisture, and phenology. However, this is most likely due to the limited temporal extent of the TDM data and the lack of reliable, quantitative information about tree canopies. Structural and moisture properties of the canopy are expected to have a significant effect on radar penetration, but dedicated follow up studies are needed before that impact can be measured.

Implementation Aspects and InSAR Data Considerations
In this study, phase height was estimated from TDM spotlight data using a lowresolution DTM derived from the same data, using a large averaging window. This approach was selected to reduce vegetation bias in the reference height model; it provided meaningful results thanks to the low canopy cover (~15%) of parklands and relatively flat topography. However, some vegetation bias could still be observed, especially for areas with tall trees with wide canopies (see Section 5.1), and the effect of topographic undulations was not studied at all. Future work should focus on improving the DTM estimation methodology used in this paper and/or synergy with topographic data provided by the current GEDI and future BIOMASS missions [39,62].
This study assessed the potential of TDM for tree height estimation, and some operational aspects were not addressed. This includes the delineation of trees in TDM data, which in this paper it was done using in situ measured position and crown diameter. Using remotely sensed estimates of tree position and crown diameter, e.g., from high resolution optical satellite data rather than in situ data is expected to generate additional uncertainties in the estimation, but it is outside the scope of this study. Furthermore, this paper disregarded geolocation and co-registration inaccuracies, shadowing of entire trees (e.g., small tree located underneath a larger tree), errors in in situ measurements, and numerical errors introduced during InSAR processing and modelling.
The proposed model-based approach to InSAR processing compensates partly for the side-looking geometry of SAR and provides an improvement in both tree positioning and tree height estimation performance, as compared with only using phase height. However, it requires complex processing and substantial InSAR data: multi-temporal, spotlight-mode acquisitions were used to ascertain high resolution and stable TLM inversion in sparsely forested areas, while the combination of ascending and descending data allowed height estimation in areas shadowed from one of the directions.
Depending on the application, phase height may be a sufficient proxy for tree height, in particular if adequate training data are available and if most trees are of similar shape, size, and structure, so that a constant ground range offset correction may be applied. However, phase height is affected by ground scattering to a larger degree than mean canopy height, which can introduce bias effects related to both ground properties and canopy cover. These, in turn, can lead to unexpected results in the data, like the negative phase height values observed in very sparsely areas in boreal forests [29]. The negative phase height values observed in Figure 7b are likely due to the combination of height calibration uncertainties and the aforementioned ground scattering effects.
Future work should address tree height estimation with the stripmap-mode TDM data used to create the global DEM [25]. These data are more abundant, and they provide a substantial advantage in terms of spatial coverage (typically 30 km × 50 km, as opposed to 10 km × 5 km for the spotlight-mode data used in this paper), albeit at the cost of azimuth resolution (typically 3.3 m, compared with 1.1 m for the spotlight-mode data).
Note that in this study, the aspect of polarimetry was ignored because only HHpolarised data were available. However, polarimetric data may provide significant additional information in parklands: the vertically oriented trunks are expected to be more exposed in these sparsely forested and relatively dry areas, thus potentially causing more polarimetric diversity at X-band than in more densely forested areas. This prospect should be addressed in follow-up studies.

Conclusions
This study investigated the potential of using single-polarised TanDEM-X spotlightmode data with 1.1 m azimuth resolution for mapping tree height and position of trees of multiple species in a sparse canopy cover (~15%) parkland environment in Burkina Faso. A new, model-based InSAR processing approach was developed for this study, providing a high spatial resolution (~3 m) mean canopy elevation map for a 4200 hectare test area. These data will be used in an ongoing study focussing on the effect of trees on agriculture in parkland environments.
Tree height was estimated with a standard error of 2.8 m (or 31% of the average tree height of 9.0 m), when evaluated against in situ data from 915 trees which were from 15 species/genera and using an empirical model. Systematic variations were studied across species/genera and explained with two effects: differences in crown shape and vegetation bias from the DTM. Further improvement in tree height estimation was obtained by combining TanDEM-X data with in situ data on crown diameter, diameter at breast height, and/or species, which significantly reduced the observed biases and promised synergies with other sensors (e.g., optical data with higher spatial resolution in the future and augmented bands) and/or long-term monitoring of well-known test sites with selected in situ measurements.
To fully explore the potential of the existing TanDEM-X data, future work should focus on adapting the methodology to polarimetric and/or stripmap-mode datasets. Time series spanning different seasons should provide important insights into the effect of phenology and moisture changes in trees. Future work should also consider synergies with other remote sensing data for tree crown delineation and more accurate removal of topographic effects. Funding: This work was funded by the Swedish National Space Agency project: "An integrated approach to explore the unknown role of trees in dryland crop production" (dnr: 112/16).