# Assessment of GEDI's LiDAR Data for the Estimation of Canopy Heights and Wood Volume of Eucalyptus Plantations in Brazil

Ibrahim Fayad, Nicolas Baghdadi, Clayton Alcarde Alvares, Jose Luiz Stape, Jean-Stéphane Bailly, Henrique Ferraco Scolforo, Mehrez Zribi, Guerric Le Maire

<figure><img src="image_1.png" /></figure>

## To cite this version:

Ibrahim Fayad, Nicolas Baghdadi, Clayton Alcarde Alvares, Jose Luiz Stape, Jean-Stéphane Bailly, et al.. Assessment of GEDI's LiDAR Data for the Estimation of Canopy Heights and Wood Volume of Eucalyptus Plantations in Brazil. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14, pp.7095-7110. 10.1109/JSTARS.2021.3092836. hal-03318147

HAL Id: hal-03318147

https://hal.inrae.fr/hal-03318147

Submitted on 16 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assessment of GEDI's LiDAR Data for the Estimation of Canopy Heights and Wood Volume of Eucalyptus Plantations in Brazil

Ibrahim Fayad ID, Nicolas N. Baghdadi, Clayton Alcarde Alvares ID, Jose Luiz Stape, Jean Stéphane Bailly ID, Henrique Ferraço Scolforo, Mehrez Zribi, and Guerric Le Maire ID

**Abstract**—Over the past two decades spaceborne LiDAR systems have gained momentum in the remote sensing community with their ability to accurately estimate canopy heights and aboveground biomass. This article aims at using the most recent global ecosystem dynamics investigation (GEDI) LiDAR system data to estimate the stand-scale dominant heights ($H_{dom}$), and stand volume (V) of Eucalyptus plantations in Brazil. These plantations provide a valuable case study due to the homogenous canopy cover and the availability of precise field measurements. Several linear and non-linear regression models were used for the estimation of $H_{dom}$ and V based on several GEDI metrics. $H_{dom}$ and V estimation results showed that over low-slopped terrain the most accurate estimates of $H_{dom}$ and V were obtained using the stepwise regression, with an root-mean-square error (RMSE) of 1.33 m ($R^2$ of 0.93) and 24.39 $m^3 \cdot ha^{-1}$ ($R^2$ of 0.90) respectively. The principal metric explaining more than 87% and 84% of the variability ($R^2$) of $H_{dom}$ and V was the metric representing the height above the ground at which 90% of the waveform energy occurs. Testing the postprocessed GEDI metric values issued from six available different processing algorithms showed that the accuracy on $H_{dom}$ and V estimates is algorithm dependent, with a 16% observed increase in RMSE on both variables using algorithm a5 vs. a1. Finally, the choice to select the ground return from the last detected mode or the stronger of the last two modes could also affect the $H_{dom}$ estimation accuracy with 12 cm RMSE decrease using the latter.

Manuscript received November 13, 2020; revised January 11, 2021 and May 5, 2021; accepted June 15, 2021. Date of publication June 28, 2021; date of current version July 26, 2021. This work was supported in part by the French Space Study Center (CNES, TOSCA 2020 project), and in part by the National Research Institute for Agriculture, Food, and the Environment (INRAE). (Corresponding author: Ibrahim Fayad.)

Ibrahim Fayad and Nicolas N. Baghdadi are with the French National Research Institute for Agriculture, Food and the Environment (INRAE), CIRAD, CNRS, TETIS, AgroParisTech, Université de Montpellier, 34093 Montpellier, France (e-mail: ibrahim.fayad@inrae.fr; nicolas.baghdadi@inrae.fr).

Clayton Alcarde Alvares is with the UNESP, Faculdade de Ciências Agronômicas Botucatu 18610-034, Brazil, and also with the Suzano SA, Limeira 13465-970, Brazil (e-mail: calcarde@suzano.com.br).

Jose Luiz Stape is with the UNESP, Faculdade de Ciências Agronômicas, Botucatu 18610-034, Brazil (e-mail: jlstape@gmail.com).

Jean Stéphane Bailly is with the INRAE, IRD, Institut Agro, LISAH, Université de Montpellier, 34060 Montpellier, France, and also with the AgroParisTech, 75005 Paris, France (e-mail: bailly@agroparistech.fr).

Henrique Ferraço Scolforo is with the Suzano SA, Limeira 13465-970, Brazil (e-mail: hscolforo@suzano.com.br).

Mehrez Zribi is with the Center for the Study of the Biosphere from Space (CNRS/UPS/IRD/CNES/INRAE), 31401 Toulouse, France (e-mail: mehrez.zribi@ird.fr).

Guerric Le Maire is with the CIRAD, UMR Eco&Sols, 34398 Montpellier, France, and also with the Eco&Sols, CIRAD, INRA, IRD, Montpellier SupAgro, Université de Montpellier, 34060 Montpellier, France (e-mail: guerric.le_maire@cirad.fr).

Digital Object Identifier 10.1109/JSTARS.2021.3092836

**Index Terms**—Brazil, dominant heights, eucalyptus, global ecosystem dynamics investigation (GEDI), LiDAR, wood volume.

## I. INTRODUCTION

IN THE last couple of decades, global concerns on the increased atmospheric concentration of greenhouse gases, such as CO$_2$ has risen the interest in quantifying the state and change of forest resources due to the key role of forests in the global carbon cycle [1], [2]. Forests sequester a large quantity of carbon in their woody biomass where they store around 70% to 90% of the global terrestrial biomass ranging from 385x10$^9$ to 650x10$^9$ Mg [3]. Hence, the accurate estimation of forest biomass is needed to better determine its precise role in the global carbon cycle [4], [5]. Forest plantations represent a small fraction (6.9%) of the total forested land ([6]) but are becoming increasingly important around the world, economically, socially and environmentally ([7], [8]).

The primary source of above ground biomass (AGB) estimation in tropical forests at large scales came in the last years from observations and measurements from different satellite remote sensing platforms. Methods based on remotely sensed data are less accurate than field measurements, however, their major advantages are their global and frequent coverage and the low or free acquisition costs for the end user. Currently optical, radar, and LiDAR are the three main sources of remotely sensed data used in AGB estimation techniques. Nonetheless, current data sources are either limited to low AGB levels (<150 Mg/ha) (sensor saturation at certain biomass levels with radar and optical data) or have a limited spatial coverage (e.g., airborne LiDAR data). LiDAR systems either airborne or spaceborne have the capability to capture the horizontal and vertical structure of vegetation comprehensively [9], and can thus estimate biomass with better precision in comparison to the techniques using radar or optical data [10], [11]. To date, there have been only three satellite LiDAR missions. The first mission was the Ice, cloud, and land elevation satellite (ICESat-1) which carried the geoscience laser altimeter system (GLAS) from 2003 until 2009 [12]. Although GLAS's ~60 m diameter footprint was larger than the ideal resolution for forest observations [13], its capability to estimate forest parameters (e.g., canopy heights and biomass) has been exploited in numerous studies during its operational and post-operational periods [5], [14]–[20].
ICESat-1 was followed in 2018 by ICESat-2 that carried the advanced topographic laser altimeter system (ATLAS) with a goal to measure ice-sheet topography, cloud and atmospheric properties and global vegetation. However, the wavelength of the equipped laser (532 nm) has a spectral region of high radiation absorption by the vegetation. This results in a low number of reflected photons measured by ATLAS over vegetation [21], and limits its ability to estimate forest canopy heights [21].

The most recent spaceborne LiDAR system is GEDI on board the ISS, which was launched in December 2018 with on-orbit checkout in April 2019. GEDI's mission objective is to provide information about canopy structure, biomass and topography, and is estimated to acquire 10 billion cloud free shots in its two years mission [22]. GEDI measures vertical structures similar to ICESat-1 (i.e., waveforms). However, given GEDI's higher sampling rate (242 versus 40 Hz for ICESat-1), and the much smaller footprint size (~25 versus ~60 m for ICESat-1), GEDI provides a highly improved coverage and waveform precision.

GEDI's ability to estimate forest height and wood volume on different types of forest ecosystems, topography and latitudes is of paramount importance. GEDI datasets are organized in different levels of products, from raw acquisition data to more elaborated data obtained by performing signal analysis and metrics extraction from the waveforms. This results in a large number of metrics for each acquired footprint, from which, many different models could be used to retrieve canopy heights and wood volume. While direct metrics could be used as good proxies, it is however acknowledged that combining different metrics yields higher accuracies. For instance, such algorithms make use of linear or non-linear regression models applied on sets of metrics extracted from GEDI waveforms, and eventually combined with digital elevation models (DEMs). The full waveform LiDAR data can potentially give access to more information on canopy structure than the basic “top” and “bottom” return signals, being itself potentially informative for canopy height and volume prediction. Therefore, it is critical to explore which metrics, or combination of metrics, and with which type of models (e.g., linear versus nonlinear) provide the best forest parameter estimates. It is also important to evaluate the effect of the uncertainty of the metrics estimation themselves, which results from differences in preprocessing algorithms, as well as other acquisition characteristics that may influence the final models, such as beam acquisition angles.

Precise evaluation of forest height and volume is not an easy task. One of the main issues is that uncertainties in field measurements can propagate through the models and create larger uncertainties in the estimates [23]. For example, Saarela et al. [24] and Holm et al. [25] found that not accounting for errors in field measurements could underestimate the uncertainty in final satellite-based AGB maps by a factor of three or more. Feldpausch et al. [26] and Kearsley et al. [27] found that uncertainties in tree height measurements led to increased bias in the biomass and carbon stock estimates. Other obstacles include: the influence of tree growth during the timespan between the field measurements and satellite acquisitions which cannot be neglected [28]; the comprehensive model validation limited by the sparsity of *in situ* data [29]; and the method used to measure tree heights [30].

Consequently, in such evaluation and comparison of processing algorithms and models for forest height and volume estimation, the field dataset plays a critical role. In this article, we want to focus on the uncertainty coming from the GEDI metrics and models, minimizing the influence of the uncertainty on *in situ* measurements. To reach this objective, we analyzed a large dataset of forest plantations in Brazil, which has many advantages to serve as a test case: large number of sites, different climate and topographical environments, numerous and frequent measurements, precise measurements of tree heights, good allometric relationships for wood volume and homogeneous canopies, etc. (see description in Section II-C). And even though not representative of all forests, the results obtained on forest plantations can give a notion of the *reachable* precision on height and wood volume estimation using GEDI data on a more structurally simple forest than natural forests, while also removing part of the errors due to *in situ* measurements.

The main objectives of this article are therefore summarized in the following questions.

1. What are the more important GEDI metrics linked to canopy height and volume?

2. Are linear and nonlinear models using subsets of metrics more efficient in predicting height and volume?

3. What is the importance of the different pre-processing algorithms on the final uncertainties?

4. Is there an influence of other acquisition characteristics, such as viewing angle on the estimated forest characteristics?

5. Are other stand information, such as data from DEM or age of the stand relevant for the estimation of height and volume on forest plantations?

The manuscript presents first the GEDI dataset, followed by the processing of GEDI data and the main metrics that will be used for the estimation of canopy heights and wood volume. Next, a description of the used methods for the estimation of the forest characteristics is presented in Section II. Finally, the results, discussions, and main conclusions are presented in Sections IV, V, and VI, respectively

## II. STUDY SITE AND DATASETS

### A. Study Area

The study area is located in four regions in Brazil, (Bahia & Espírito Santo, Mato Grosso do Sul, São Paulo, and Maranhão) across a large latitudinal gradient (see Fig. 1) and covering different climate and soil types. The studied plantations are managed in order to produce high yield pulpwood growing at short rotations. Clonal seedlings of mainly *E. grandis* (W. Hill) and *E. urophylla* (S.T. Blake) and different types of hybrids are planted in rows at a density of 1000–1667 trees/ha, rationally fertilized with nitrogen, phosphorus, and potassium and micronutrients to alleviate any nutritional limitations. Harvest occurs every six to seven years, and very little tree mortality (under 7% from original plantation) is noticed. The annual productivity of the plantations was on average 40 m³/ha/year, with 80% of the stands being between 30–50 m³/ha/year and some stands could reach values as high as 60 m³/ha/year. At harvest time, the stand volume is between 180 and 300 m³/ha, with a dominant height
<figure><img src="image_3.png" /><figcaption>Fig. 1. (a) Location of the four study sites. (b) Example of GEDI tracks over some stands. (c) Eucalyptus stand during harvest (approx. 30 m high) illustrating the clearly separated crown and trunk strata.</figcaption></figure>

of 20 to 35 m range (for 80% of the stands). These plantations were managed locally by stand units, generally around 50 ha, where the same management is applied: planting, harvesting, weed control, genetic material, soil preparation and fertilization. There are generally sparse understory and herbaceous strata in these plantations, as result of chemical weeding the first year, the closing of the canopy, and the high competitive strength of Eucalyptus. Tree height is very homogeneous within a stand, with 95% of the trees having heights at +/- 10.5% around the average tree height in plot inventories. The plantations exhibit a simple structure, with a tree crown strata of 3 to 10 m in width above a “trunk strata” with few Eucalyptus leaves and few understories [see Fig. 1(c)]. The “soil strata” is mainly constituted of litter accumulation of branches and leaves, with some patches of herbaceous species.

## B. GEDI Data

1) Processing of GEDI Waveforms: GEDI uses three on-board lasers that produce eight parallel tracks of observations. GEDI lasers illuminate a surface or footprint on the ground with a 25 m diameter, at a frequency of 242 Hz, over which 3D structures are measured. The footprints are separated by ~60 m (center to center) along the track, and the tracks are separated by ~600 m across. Moreover, GEDI has the ability to rotate up to six degrees, allowing the lasers to be pointed as much as 40 km on either side of the ISS’s ground track [22]. GEDI measures vertical structures using a 1064-nm laser pulse, and

the echoed waveforms are digitized to a maximum of 1246 bins with a vertical resolution of 1 ns (15 cm), corresponding to a maximum of 186.9 m of height ranges, with a vertical accuracy over relatively flat, non-vegetated surfaces of ~3 cm [31].

As described in the algorithm theoretical basis document (ATBD) [32], [33], the received waveforms are first smoothed to reduce the noise in the signal, and thus permitting the determination of the useful part of the waveform within the corresponding footprint. Waveform smoothing is performed by means of a Gaussian filter with various widths. As mentioned in the ATBD, currently a width of 6.5 ns was used for the Gaussian filter (Smooth width). After smoothing, two locations in the waveform denoted as *search start* and *search end* are determined [see Fig. 2(a)]. *search start* and *search end* are, respectively, the first and last positions in the signal where the signal intensity is above the following threshold:

$$ threshold = mean + \sigma \cdot v \quad (1) $$

where “mean” is the mean noise level, “σ” is the standard deviation of noise of the smoothed waveform, and “v” is a constant currently set at 4. After determining the locations of *search start* and *search end*, the region between them, denoted as the waveform extent, is extended by a predetermined number of sample bins, currently set to 100 bins at both sides. Within the waveform extent, the highest (*toploc*) and lowest (*botloc*) detectable returns are determined [see Fig. 2(a)]. The metrics *toploc* and *botloc* respectively represent the highest and lowest locations within the waveform extent where two adjacent
<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) Example of an acquired GEDI waveform (Rw) over a Eucalyptus stand ($H_{dom}$ = 25.9m; V = 230.7 m³.ha⁻¹), its smoothing (Sw) and corresponding waveform metrics. (b) Cumulative energy of the waveform (CE) between *botloc* and *toploc* and the corresponding relative heights ($RH_n$) at different percentages "n" for the same waveform. One (1) ns corresponds to 15 cm sampling distance in the waveform. The waveform amplitudes are counts from the analog to digital converter on the instrument.</figcaption></figure>

intensities are above a threshold. The threshold equation used to determine *toploc* and *botloc* is the same as (1), with “v” an integer fixed at 2, 3, 4, or 6 (depending on the processing configuration). In the ATBD, the value of “v” used to determine *toploc* is named “Front_threshold” and “back_threshold” for *botloc*. Waveform metric values are extracted using thresholds on Smoothwidth_zcross, front_threshold, and back_threshold. Currently, there are six configurations (henceforth referred to as algorithms) of different thresholds on these variables, which are used to determine waveform metrics with high precision for a variety of acquisition scenarios (see Table I). Finally, the location of the distinctive peaks or modes in the waveform, such as the ground peak, or top of canopy peaks is determined using a second Gaussian filtering of the waveform section between *toploc* and *botloc*, and then finding all the zero crossings of the first derivative of the filtered waveform [see Fig. 2(a)]. The width of the second Gaussian filter (“Smoothwidth_zcross”) is fixed to either 3.5 or 6.5 ns (based on the algorithm used).

Finally, the position of the ground return within the waveform is determined using the position of the last detected mode. The six different algorithms, noted a1 to a6 correspond to different values of the above-mentioned parameters (see Table I) and lead to different estimates of the waveform metrics, and could in turn lead to six different canopy height estimates. Over forest stands, the recorded waveforms are multimodal in shape, with each mode representing a reflection from a distinct surface height. Fig. 2(a) shows a typical waveform over a Eucalyptus forest stand on relatively flat terrain. Over flat terrain, the first Gaussian corresponds to a reflection from the top of the canopy while the last Gaussian mostly refers to the lowest point in the footprint, i.e., the ground surface.

GEDI data used in this article have been already processed and published by the land processes distributed active archive center (LP DAAC). Currently, three products (L1B, L2A, and L2B) are available for download. The L1B data product [32] contains detailed information about the transmitted and received
<table><caption>TABLE I
DIFFERENT THRESHOLDS USED IN EACH OF THE SIX ALGORITHMS FOR THE ANALYSIS OF THE RECEIVED WAVEFORMS</caption><thead><tr><th>Algorithm</th><th>Smooth width</th><th>Smoothwidth</th><th>zcross</th><th>Front threshold</th><th>Back threshold</th></tr></thead><tbody><tr><td>a1</td><td>6.5</td><td>6.5</td><td></td><td>3</td><td>6</td></tr><tr><td>a2</td><td>6.5</td><td>3.5</td><td></td><td>3</td><td>3</td></tr><tr><td>a3</td><td>6.5</td><td>3.5</td><td></td><td>3</td><td>6</td></tr><tr><td>a4</td><td>6.5</td><td>6.5</td><td></td><td>6</td><td>6</td></tr><tr><td>a5</td><td>6.5</td><td>3.5</td><td></td><td>3</td><td>2</td></tr><tr><td>a6</td><td>6.5</td><td>3.5</td><td></td><td>3</td><td>4</td></tr></tbody></table>

waveforms, the location and elevation of each waveform footprint and other ancillary information, such as mean and standard deviation of the noise and acquisition time. The L2A product [33] contains data of elevation and height metrics of the vertical structures within the waveform. These height metrics are issued from the processing of the received waveforms from the L1B product. Finally, the L2B data product [34] provides footprint-level vegetation metrics, such as canopy cover, vertical profile metrics, plant area index and foliage height diversity.

In this article, the received waveforms, their geolocation (longitude, and latitude), as well as their acquisition times were extracted from the L1B data product. In the L2A data product, the derived metrics are also grouped by algorithm. Therefore, for each beam, the metrics derived from each of the six algorithms, as well as the parameters used for each algorithm are available. Therefore, we extracted from L2A for each beam and for each of the six algorithms (a1 through a6), the following variables.

1) The position within the waveform of *toploc* and *botloc*.

2) The amplitude of the smoothed waveforms lowest detected mode (zcrosssub).

3) The quality flag for each waveform (quality_flag).

4) The number of detected modes (num_detectedmodes).

5) The position and amplitude of each detected mode.

6) The Relative height metrics at 10% intervals from *botloc* (0%) to *toploc* (100%) (RHn, 10% ≤ n ≤ 100%, step 10%). RHn represents the height between *botloc* and the location at n% of cumulative energy [see Fig. 2(b)]. No metrics were extracted from the L2B product, as they were not relevant to this article.

2) *Calculation of Relevant GEDI and Terrain variables:* Several linear and non-linear regression models will be tested in order to estimate the stand dominant height $H_{dom}$ (m) and stand merchantable wood volume V (m³·ha⁻¹) from GEDI data. The models were tested with a *priori* variables that were extracted from GEDI waveforms. These variables represent canopy features, such as canopy top, canopy trunks, ground, or a mix of these elements. In addition to the available GEDI waveform metrics described in the previous section, several additional metrics were also extracted. The first is the waveform extent ($W_{ext}$) which is the height difference between *botloc* and *toploc*. Next, to remove the effects of canopy height variability and terrain slope, two indices relying on waveform structure were determined. The leading edge extent (Lead$_{ext}$) is, as defined by Hilbert and Schmullius [35], the difference between the position

<figure><img src="image_5.png" /><figcaption>Fig. 3. Allometric relation between *in situ* V and $H_{dom}$.</figcaption></figure>

of the first mode (*Vloc*) and *toploc*, while the trailing edge extent (Trail$_{ext}$) is the difference between *botloc* and the ground return (*Gloc*) [35]. Two methods will be used to determine the ground position: the position of the mode selected as the lowest non-noise mode from the L2A data product, and the position of the highest mode between the last two detected modes. The viewing angle of GEDI, which represents the angle between the looking direction of the instrument and nadir at acquisition time, and for each shot, has also been calculated using the geolocation of the GEDI instrument available from the L1B data product. The viewing angle has been demonstrated in Urban *et al.* [36] to increase elevation errors for ICESat-1 GLAS when the viewing angle deviates from nadir due to precision attitude determination.

Finally, as the wood volume (V) increases with canopy height in a nonlinear shape (see Fig. 3), we calculated RH$_n$ for several power values ($RH_n^p$, $1 < p \le 3$, step 0.2).

All the used variables for the estimation of the stand dominant height $H_{dom}$ (m) and stand merchantable wood volume V (m³·ha⁻¹) are given in Table II.

The values of the extracted waveform metrics vary by the algorithm used for the processing of the waveforms. Therefore, in this article, all GEDI metrics were determined for the six available algorithms. The variability of the metric values based on the processing algorithm is given in Table III.

3) *Filtering of GEDI Waveforms:* Not all GEDI acquisitions are viable, as atmospheric conditions (e.g., clouds) can affect them. Therefore, a waveform was not investigated further if it met any of the following criteria.

1) Waveforms with reported elevations that are significantly higher or lower than the corresponding elevations from the SRTM DEM [16]. In essence, we removed all waveforms where the absolute difference is higher than 100 m.
<table><caption>TABLE II
LIST OF ALL THE VARIABLES CALCULATED FROM GEDI WAVEFORMS</caption><thead><tr><th>Variable</th><th>Unit</th><th>Abbreviation</th></tr></thead><tbody><tr><td>Search Start</td><td>bin</td><td>-</td></tr><tr><td>Search End</td><td>bin</td><td>-</td></tr><tr><td>Vegetation mode location</td><td>bin</td><td>Vloc</td></tr><tr><td>Ground mode location</td><td>bin</td><td>Gloc</td></tr><tr><td>Sginal start (Top of canopy)</td><td>bin</td><td>toploc</td></tr><tr><td>Signal end</td><td>bin</td><td>botloc</td></tr><tr><td>Waveform extent</td><td>m</td><td>W<sub>ext</sub></td></tr><tr><td>Trailing edge extent [44]</td><td>m</td><td>Trail<sub>ext</sub></td></tr><tr><td>Leading edge extent [44]</td><td>m</td><td>Lead<sub>ext</sub></td></tr><tr><td>Relative canopy height (botloc (0%) to toploc (100%)) at n% of cumulative energy</td><td>m</td><td>(RH<sub>n</sub>, 10% ≤ n ≤ 100%, step 10%)</td></tr><tr><td>3x3 pixels terrain index</td><td>m</td><td>T<sub>i</sub></td></tr><tr><td>3x3 pixels terrain roughness</td><td>m</td><td>ROUG</td></tr><tr><td>Terrain slope</td><td>%</td><td>S</td></tr><tr><td>Viewing angle</td><td>°</td><td>VA</td></tr></tbody></table>

Variables to be used as predictor variables in the canopy height and wood volume estimation models are highlighted in gray.

<table><caption>TABLE III
MEAN AND STANDARD DEVIATION OF SOME GEDI METRIC VALUES FROM EACH OF THE SIX PROCESSING ALGORITHMS USING ALL GEDI SHOTS OVER THE 566 SELECTED EUCALYPTUS PLANTATIONS</caption><thead><tr><th rowspan="3">Variable</th><th rowspan="3">Unit</th><th colspan="12">Algorithm</th></tr><tr><th colspan="2">a1</th><th colspan="2">a2</th><th colspan="2">a3</th><th colspan="2">a4</th><th colspan="2">a5</th><th colspan="2">a6</th></tr><tr><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th></tr></thead><tbody><tr><td>W<sub>ext</sub></td><td>(m)</td><td>25.00</td><td>5.50</td><td>28.74</td><td>6.66</td><td>25.00</td><td>5.50</td><td>23.91</td><td>6.01</td><td>32.99</td><td>8.70</td><td>26.78</td><td>5.91</td></tr><tr><td>Lead<sub>ext</sub></td><td>(m)</td><td>4.89</td><td>2.23</td><td>4.34</td><td>1.08</td><td>4.61</td><td>2.05</td><td>3.80</td><td>1.01</td><td>4.34</td><td>1.08</td><td>4.38</td><td>1.21</td></tr><tr><td>Trail<sub>ext</sub></td><td>(m)</td><td>4.24</td><td>1.60</td><td>6.40</td><td>2.84</td><td>3.86</td><td>1.65</td><td>4.24</td><td>1.60</td><td>1.93</td><td>1.89</td><td>3.45</td><td>2.15</td></tr><tr><td>RH<sub>100</sub></td><td>(m)</td><td>20.73</td><td>5.32</td><td>22.31</td><td>5.97</td><td>21.11</td><td>5.40</td><td>19.64</td><td>5.72</td><td>31.02</td><td>9.43</td><td>23.30</td><td>6.42</td></tr></tbody></table>

<table><caption>TABLE IV
DISTRIBUTION OF GEDI SHOTS ACROSS THE FOUR STUDY REGIONS</caption><thead><tr><th>Site</th><th>GEDI count</th><th>in situ H<sub>dom</sub><br>(m)</th><th>in situ V<br>(m<sup>3</sup>.ha<sup>-1</sup>)</th></tr></thead><tbody><tr><td>BA&ES</td><td>1220</td><td>25.7</td><td>207.9</td></tr><tr><td>MA</td><td>808</td><td>27.2</td><td>208.0</td></tr><tr><td>MS</td><td>2586</td><td>29.4</td><td>300.0</td></tr><tr><td>SP</td><td>1068</td><td>32.4</td><td>380.6</td></tr></tbody></table>

in situ H<sub>dom</sub> and in situ V represent the 95th percentile of in situ values for each site.

2) Waveforms with a difference between waveform extent (W<sub>ext</sub>) and (Gloc-Vloc) higher than 400 bins (corresponding to 60 m)

A total of 6166 footprints were acquired over our reference stands between April 2019 and September 2019, with the majority of these footprints (92.15%) providing exploitable waveforms. Table IV gives the distribution of GEDI shots across the four regions.

GEDI data accessible through NASA's LP DAAC contain a quality flag (quality_flag) for each acquired waveform. A waveform with a quality flag set to “1” indicates that the waveform

meets certain criteria based on energy, sensitivity, amplitude, and real-time surface tracking quality, and thus can be processed further. However, in this article waveforms with either value of the quality_flag were analyzed.

### C. Inventory Measurements

A total of 566 Eucalyptus stands were selected, corresponding to stands where GEDI footprints acquired between April 20, 2019 and September 4, 2019 were totally included. An additional 50 m internal buffer strip from the stand borders was used to account for any footprint geolocation errors and to avoid footprints that match the boundary between the stand of interest and the surrounding medium. These 566 Eucalyptus stands were also selected because they had field inventories performed by the company close to GEDI's acquisition date (time difference fewer than two months). Field inventories are performed on several permanent inventory plots within each stand. These inventory plots are systematically distributed throughout the stand with a density of one plot per 10 ha (i.e., a 20 ha stand will have two inventory plots while an 80 ha plot will contain eight inventory plots). These permanent inventory plots had each an area of approximately 400 m² including 30 to 100 trees (average of 58
<figure><img src="image_6.png" /><figcaption>Fig. 4. Relationship between stand age and (a) $H_{dom}$ (m) and (b) Volume (m³.ha⁻¹).</figcaption></figure>

trees) in function of the inventory plot size and planting density. During a field inventory, the diameter at breast height (DBH, 1.3 m above the ground) of each tree in the inventory plot, the height of a central subsample of 10 trees and the height of the 10% largest trees in terms of DBH (dominant trees) were measured. The mean height of the 10% of the largest trees defined the dominant height of the plot ($H_{dom}$), while the mean height of all trees in the plot (measured + estimated) defined the average height of the plot ($H_{mean}$). $H_{dom}$, basal area and age on the inventory date were then used in local volume equations to estimate the plot total and merchantable volume (merchantable volume is the tree volume up to the diameter outside bark of 6 cm). Stem biomass was then estimated from the stem volume using age-dependent estimates of wood biomass density.

As the dates of the inventory measurements were different from GEDI acquisition dates, only data with a difference fewer than two months in the date between GEDI acquisitions and inventory were used. In fact, on these fast growing plantations, a two-month difference could result in an up to 50 cm growth in $H_{dom}$ [see Fig. 4(a)] and 10 m³.ha⁻¹ in V [see Fig. 4(b)]. However, this reasonable compromise allows keeping a large number of stands including a large variability of age and growing conditions. Fig. 5 shows the distribution of field measured $H_{dom}$ and volume.

## D. Digital Elevation Model

The DEM with a spatial resolution of 30 m, derived from the Shuttle Radar Topography Mission (SRTM), was used in this article. Three variables were derived from the DEM: slope

<figure><img src="image_7.png" /><figcaption>Fig. 5. Distribution of measurements of (a) dominant canopy heights and (b) wood volume from field inventories of the 566 Eucalyptus stands. The Y-axis represents the percentage of samples within each height range (a) and wood volume (b).</figcaption></figure>

(S), Terrain Index ($T_I$), and surface Roughness (Roug). The $T_I$ map was obtained by calculating the difference between the highest and lowest altitude in a 3 × 3 pixel-moving window. The surface roughness map was obtained by computing the standard deviation of the elevation in a 3 × 3 pixel-moving window.

## III. METHODOLOGY

### A. Forest Height Estimation

The simplest method to estimate $H_{dom}$ from a GEDI waveform over forest stands with a gently sloping terrain is the height difference between signal start (toploc) and the ground position (Gloc) [37]

$$H_{dom} = \text{toploc} - \text{Gloc.} \quad (2)$$

For the previous ICESat-1 GLAS waveforms, the ground return was assumed to be the stronger of the last two detected modes [15]. Therefore, in this article, the stronger between the last two detected modes, as well as the position of the mode in the field “selected_mode (SM)” from the L2A data product will be considered as the ground return separately.

Estimating canopy heights using (2) has several caveats. For example, over sloping terrain, the ground peak becomes wider, and the returns from ground and vegetation can be mixed in the case of large footprints, making the identification of the ground peak return difficult and the estimation of forest height inaccurate ([5], [37]). To remove or minimize the terrain slope effect on the waveforms, as well as the vegetation variability,
statistical approaches have been developed and used in several studies to predict canopy heights from GLAS data (e.g., [5], [15], [35], [38], [39]). These approaches proposed regression models based either on only waveform metrics or on both waveform metrics and terrain information derived from DEMs.

The first statistical model was developed by Lefsky et al. [5] to estimate the maximum canopy height ($H_{dom}$) from GLAS waveforms

$$H_{dom} = aW_{ext} - bT_{I}. \quad (3)$$

The coefficients a and b are fitted using least squares regression ($H_{dom}$ given by inventory measurements, $W_{ext}$ is derived from the GEDI waveform, and $T_{I}$ is calculated from the SRTM DEM, see Section II-D). For our dataset, $T_{I}$ values calculated from the SRTM DEM ranged from 1 to 46 m. The incorporation by Payn et al. [6] of the waveform leading edge extent in (4) showed a slight improvement on canopy height estimation

$$H_{dom} = aW_{ext} - bT_{I} + cLead_{ext}. \quad (4)$$

Over sloping terrain, Lefsky et al. [38] observed that the waveform extent is insufficient for estimating canopy heights. Hence, a new model based on the waveform extent, leading edge extent, and trailing edge extent was proposed. However, Pang et al. [39] observed inaccurate estimates of canopy heights with the improved model by Lefsky et al. [38], especially for small waveform extents, and thus proposed a simpler model to estimate canopy heights using the following equation:

$$H_{dom} = aW_{ext} - \{b(Lead_{ext} + Trail_{ext})\}^{c}. \quad (5)$$

The nonlinear model by Pang et al. [39] was further simplified by Chen [15]

$$H_{dom} = aW_{ext} - b(Lead_{ext} + Trail_{ext}) \quad (6)$$

Baghdadi et al. [16] tested additional models for the estimation of canopy heights using ICESat-1 GLAS waveforms, of which, two will be tested in this article. The first model uses the $Trail_{ext}$ and $T_{I}$

$$H_{dom} = aW_{ext} - bT_{I} + cTrail_{ext}. \quad (7)$$

The second model uses exclusively GEDI metrics

$$H_{dom} = aW_{ext} - bLead_{ext} - cTrail_{ext} + d. \quad (8)$$

In addition to the previously described models, a stepwise multilinear regression aimed at estimating canopy dominant heights (SRH) and automatically choosing a set of predictive explanatory variables among all possible variables presented in Table II. The choice of adding or removing a variable from the SRH model is based on the increase or decrease of the mean squared error (MSE).

We also estimated canopy dominant heights through nonlinear nonparametric regressions by means of a random forest regressor (RFH). Random forests are an ensemble of machine learning algorithms used for classification or regressing by fitting a number of decision trees on various sub-samples of the dataset, and use averaging to improve the predictive accuracy and control overfitting [40]. Compared to linear models, RF is advantageous for being able to model also nonlinear relationships (threshold

effect) between the explanatory variables. For this article, the number of trees in the RF were set to 100 trees (higher tree count slightly increased model accuracy), with a tree depth equal to the square root of the number of available factors.

Finally, since random forests are nonlinear and nonparametric, we only used the original relative heights without modification (i.e., $RH_{n}^{1}$, $10\% \le n \le 100\%$, step 10%.)

## B. Wood Volume Estimation

The estimation of aboveground biomass has been proven to be successful using ICESat-1 GLAS waveforms as demonstrated by several studies ([15]–[17]). In this article, four models were tested to estimate wood volume from GEDI waveforms based on Hdom estimates. The first model was adapted from Lefsky et al. [5] for the estimation of wood volume (instead of AGB in its original formulation), using the squared dominant canopy heights ($H_{dom}$)

$$V = a + bH_{dom}^{2}. \quad (9)$$

The second tested model was adapted from Saatchi et al. [41], and uses a power law relationship between the volume and Lorey's height

$$V = aH_{L}^{b} \quad (10)$$

where $H_{L}$ is Lorey's height which weighs the contribution of trees (all trees >10 cm in diameter) to the stand height by their basal area. In this article, the relationship defined in (10) was used by replacing Lorey's height with the dominant height as both height values were similar ($H_{L}$ was lower than $H_{dom}$ by a maximum of 0.9 m at the end of the rotation of the Eucalyptus plantation) [16]. For both models (9 and 10), the coefficients a and b were first fitted using in situ measurements of dominant height and wood volume (see Fig. 6), and then, the calibrated equations were used to estimate wood volume using the dominant height predicted from GEDI footprints (best model from Section III-A).

Similarly to Section III-A, a stepwise linear regression model (SRV) and a random forest regressor (RFV) were used to estimate the wood volume.

## C. Model Assessment

To assess how the tested models generalize to an independent data set, a five-fold cross validation was used. Large k-fold values mean less bias towards overestimating the true expected error (as training folds will be closer to the total dataset). Moreover, since there are several GEDI footprints inside each stand, and the stands are very homogeneous, the five-fold splitting was also done along the stands in order to reduce fitting bias. In essence, GEDI footprints inside the same stand were used for either training or validation. Finally, models' performance were assessed using the coefficient of determination ($R^2$), the bias (measured—estimated), the root-mean-square error (RMSE), the root mean squared percentage error (RMSPE), and the Akaike information criterion (AIC). $R^2$, RMSE, and RMSPE
<figure><img src="image_8.png" /><figcaption>Fig. 6. Comparison of measured vs. estimated $H_{dom}$ from the models presented in Section III-A using GEDI metrics extracted with algorithm a1 (see Table I). RMSE is expressed in meters (m).</figcaption></figure>

are defined as follows:

$$R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \quad (11)$$

$$\text{RMSE} = \sqrt{\frac{1}{n} \cdot \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \quad (12)$$

$$\text{RMSPE} = 100 \cdot \sqrt{\frac{1}{n} \cdot \sum_{i=1}^{n} \left( \frac{y_i - \hat{y}_i}{y_i} \right)^2} \quad (13)$$

where $y_i$ is the observed value, $\hat{y}_i$ the estimated value, $\bar{y}$ is the mean of all the observed values, and $n$ is the sample size.

The AIC proposed by Akaike [42] is a measurement of the relative goodness of fit of a statistical model to the true values. By calculating AIC values for each model, the most performant model based on the lowest AIC values can be identified.

## IV. RESULTS

### A. Canopy Height Estimation

We start our model performance analysis using GEDI metrics extracted from algorithm a1 (see Table I), and the ground location as determined from the SM field from the L2A dataset

(last detected nonnoise mode). The estimation of the canopy dominant heights ($H_{dom}$) using the linear regression models [(3) through 8] with five-fold cross validation shows an accuracy (RMSE) between 1.70 and 2.31 m with a coefficient of determination ($R^2$) between 0.80 and 0.89 (see Fig. 6). Moreover, the contribution of the trailing edge extent appeared to be higher than that of the leading edge extent [see (7) versus (4), Table V]. However, the best model between (3) through (8) was (8) (RMSE = 1.70 m and $R^2$ = 0.89) which uses both Leading and Trailing edge extents, with an independent coefficient fitted for each variable. The introduction of terrain information in the linear regression models did not show any significant improvements on the accuracy of the estimations.

The stepwise linear regression model (see Fig. 6, SRH) showed slightly better accuracy for the estimation of canopy heights (RMSE = 1.44 m, $R^2$ = 0.93) in comparison to Eq.8. However, unlike Eq.8 which relied on $W_{ext}$, $Lead_{ext}$, and $Trail_{ext}$, the most contributing variables for the estimation of the canopy heights using the SRH model were $RH_{90}$, followed by $RH_{10}$, $RH_{80}$, and $RH_{100}$. Meanwhile, the other metrics (e.g., $Lead_{ext}$, $Trail_{ext}$, TI, etc.,) were not necessary.

Furthermore the estimation of canopy heights using only $RH_{90}$ (by linear fitting) showed an RMSE of 1.63 m with an $R^2$ of 0.90, and this accuracy could be improved to an RMSE of 1.5 m ($R^2$ of 0.91) by only adding $RH_{10}^{1.8}$. The estimation of
<table><caption>TABLE V
MODELS' PERFORMANCE AND THE FITTED LINEAR EQUATIONS FOR ESTIMATING EUCALYPTUS STAND DOMINANT HEIGHTS</caption><thead><tr><th>Model</th><th>RMSE (m)</th><th>RMSPE (%)</th><th>Bias (cm)</th><th>R<sup>2</sup></th><th>AIC</th></tr></thead><tbody><tr><td>EQ. 3: H<sub>dom</sub> = 0.83W<sub>ext</sub> - 0.02T<sub>t</sub></td><td>2.30</td><td>12.4</td><td>1.7</td><td>0.80</td><td>9217.1</td></tr><tr><td>EQ. 4: H<sub>dom</sub> = 0.82W<sub>ext</sub> + 0.03T<sub>t</sub> - 0.03Lead<sub>ext</sub></td><td>2.31</td><td>12.5</td><td>1.9</td><td>0.80</td><td>9254.6</td></tr><tr><td>EQ. 5: H<sub>dom</sub> = 0.91W<sub>ext</sub> - {0.26(Lead<sub>ext</sub> + Trail<sub>ext</sub>)}<sup>0.90</sup></td><td>2.20</td><td>11.6</td><td>10.1</td><td>0.81</td><td>8710.3</td></tr><tr><td>EQ. 6: H<sub>dom</sub> = 0.90W<sub>ext</sub> - 0.21(Lead<sub>ext</sub> + Trail<sub>ext</sub>)</td><td>2.20</td><td>11.6</td><td>7.3</td><td>0.82</td><td>8683.9</td></tr><tr><td>EQ. 7: H<sub>dom</sub> = 0.97W<sub>ext</sub> + 0.041T<sub>t</sub> + 0.90Trail<sub>ext</sub></td><td>1.75</td><td>9.5</td><td>8.0</td><td>0.88</td><td>6174.2</td></tr><tr><td>EQ. 8: H<sub>dom</sub> = 0.91W<sub>ext</sub> + 0.06Lead<sub>ext</sub> - 0.98Trail<sub>ext</sub> + 1.83</td><td>1.70</td><td>9.2</td><td>1.7</td><td>0.89</td><td>5858.0</td></tr><tr><td>SRH: H<sub>dom</sub> = 0.64RH<sub>90</sub><sup>1.2</sup> - 0.18RH<sub>10</sub><sup>1.8</sup> - 0.75RH<sub>80</sub> + 0.9RH<sub>100</sub> + 6.42</td><td>1.44</td><td>7.1</td><td>0.4</td><td>0.92</td><td>4245.1</td></tr><tr><td>RFH</td><td>1.55</td><td>7.8</td><td>3.1</td><td>0.91</td><td>4887.6</td></tr></tbody></table>

The variables are described in Section II-B-2, with the models described in Section III-A

<figure><img src="image_9.png" /><figcaption>Fig. 7. Classification of the variable importance by decreasing order of importance in the RFH model for stand dominant height estimation. The importance is measured via the average percentage increase of MSE (%IncMSE) over 50 repetitions. The red bars indicates the standard deviation of %IncMSE.</figcaption></figure>

<table><caption>TABLE VI
ACCURACY (RMSE IN M) OF THE MODELS PRESENTED IN SECTION III-A FOR THE ESTIMATION OF H<sub>dom</sub> USING GEDI METRIC VALUES EXTRACTED USING THE SIX DIFFERENT ALGORITHMS (A1 THROUGH A6)</caption><thead><tr><th rowspan="3">Model</th><th colspan="12">Algorithm</th></tr><tr><th colspan="2">a1</th><th colspan="2">a2</th><th colspan="2">a3</th><th colspan="2">a4</th><th colspan="2">a5</th><th colspan="2">a6</th></tr><tr><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>Eq.3</td><td>2.30</td><td>0.80</td><td>3.82</td><td>0.44</td><td>2.53</td><td>0.75</td><td>3.35</td><td>0.57</td><td>4.96</td><td>0.06</td><td>3.09</td><td>0.64</td></tr><tr><td>Eq.4</td><td>2.31</td><td>0.80</td><td>3.77</td><td>0.46</td><td>2.54</td><td>0.75</td><td>3.32</td><td>0.58</td><td>5.01</td><td>0.04</td><td>3.00</td><td>0.66</td></tr><tr><td>Eq.5</td><td>2.20</td><td>0.81</td><td>3.32</td><td>0.58</td><td>2.47</td><td>0.77</td><td>3.11</td><td>0.63</td><td>5.06</td><td>0.02</td><td>3.07</td><td>0.64</td></tr><tr><td>Eq.6</td><td>2.20</td><td>0.82</td><td>3.52</td><td>0.52</td><td>2.46</td><td>0.77</td><td>3.24</td><td>0.60</td><td>4.69</td><td>0.16</td><td>3.08</td><td>0.64</td></tr><tr><td>Eq.7</td><td>1.75</td><td>0.88</td><td>3.50</td><td>0.53</td><td>2.40</td><td>0.78</td><td>3.13</td><td>0.63</td><td>4.49</td><td>0.23</td><td>2.98</td><td>0.66</td></tr><tr><td>Eq.8</td><td>1.70</td><td>0.89</td><td>3.08</td><td>0.64</td><td>2.50</td><td>0.76</td><td>2.67</td><td>0.73</td><td>4.06</td><td>0.37</td><td>2.79</td><td>0.70</td></tr><tr><td>SRH</td><td>1.44</td><td>0.92</td><td>1.46</td><td>0.92</td><td>1.48</td><td>0.92</td><td>1.60</td><td>0.90</td><td>1.67</td><td>0.88</td><td>1.49</td><td>0.92</td></tr><tr><td>RFH</td><td>1.55</td><td>0.91</td><td>1.58</td><td>0.90</td><td>1.72</td><td>0.89</td><td>1.75</td><td>0.88</td><td>1.80</td><td>0.88</td><td>1.61</td><td>0.90</td></tr></tbody></table>

H<sub>dom</sub> using the random forest regressor (RFH, Fig. 6) with the GEDI metrics in Table II (p in RH<sub>n</sub><sup>p</sup> was set to 1 for RFH) as the dependent variables showed an accuracy on the canopy height estimates similar to that of the SRH model.

The variable importance test of the metrics (see Fig. 7) used in RFH showed that the most contributing factors for the estimation of GEDI canopy heights is a combination of RH<sub>90</sub>, RH<sub>100</sub>, and to lesser extent RH<sub>80</sub>. These results show that in a low relief area the use of other metrics in addition to the RH<sub>90</sub> only slightly improved the precision of the estimation of canopy heights.

The estimation of H<sub>dom</sub> using the models described previously with GEDI metrics extracted from the five remaining algorithms (a2 to a6, see Table I) has been also tested. The results presented in Table VI show that for the linear regression models [(3) through (8)], canopy height estimation was worst with the metrics from algorithms a2 through a6 in comparison to the metrics from algorithm a1 with an RMSE on the canopy height estimates ranging from 2.40 m (R<sup>2</sup> of 0.78, a3) to 5.06 m (R<sup>2</sup> of 0.02, a5). This is to be expected given the low terrain relief in our study area (mean slope of 4.7 ± 3%).
<figure><img src="image_10.png" /><figcaption>Fig. 8. Comparison between Measured $H_{dom}$ and estimated $H_{dom}$ using only $W_{ext}$ values ($H_{dom} = \alpha \cdot W_{ext} + \beta$) from the six algorithms (a1 through a6). RMSE is expressed in meters (m).</figcaption></figure>

The low accuracy obtained with algorithm a5 is due to the low thresholds used for the front and back thresholds (3.σ and 2.σ, Table I), which result in larger waveform extents. This is evident when trying to estimate canopy heights based solely on the waveform extent ($H_{dom/insitu} = \alpha \cdot W_{ext} + \beta$), with the results in Fig. 8 showing that the metrics extracted using algorithm a5, especially the waveform extent ($W_{ext}$), were the least correlated to $H_{dom}$, with an RMSE of 4.38 m ($R^2$ of 0.26).

In contrast to the linear regression models, canopy height estimation using SRH or RFH with metrics from algorithms a2, a3, a5, and a6 showed accuracies similar to those obtained with algorithm 1 (see Table VI). In contrast, algorithm a5 was slightly less accurate with an RMSE of respectively 1.6 m ($R^2$ of 0.90) and 1.80 m ($R^2$ of 0.88) when using SRH and RFH.

Finally, the effects of the method to select the ground return has been studied. The results presented thus far have been based on detecting the ground mode from the SM provided in the L2A data product. SM detects the ground return as being the lowest, nonnoisy mode, which usually refers to the last detected mode. Previous studies that used GLAS waveforms suggested that the mode with the higher amplitude between the last two modes (HL2M) is a better indicator of the ground return [15], [16]. In this article, $H_{dom}$ estimation was also tested using the same models described in Section III-A, with the metrics calculated relatively to the mode with the higher amplitude between the last two modes, as the ground return. The results in Table VII show that the models relying mostly on the relative canopy heights ($RH_n$), such as SRH, or the trailing edge extent fitted separately [(7) and (8)] had an increase in accuracy on the canopy height estimates between 12 and 20 cm.

## B. Wood Volume Estimation

Four models were used to estimate the stand volume V. Two power functions as presented in Section III-B [9] and (10)], a stepwise multilinear regressing model (SRV) and a random forest based model (RFV). The coefficients of the models as

<table><caption>TABLE VII
DIFFERENCE IN ACCURACY (RMSE IN M) ON $H_{dom}$ BASED ON THE CHOICE OF THE SELECTED GROUND MODE FOR THE DIFFERENT MODELS DESCRIBED IN SECTION III-A, AND METRICS EXTRACTED USING ALGORITHM A1</caption><thead><tr><th rowspan="3">Model</th><th colspan="4">Ground Mode Selection</th></tr><tr><th colspan="2">SM</th><th colspan="2">HL2M</th></tr><tr><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>Eq.3</td><td>2.30</td><td>0.80</td><td>2.30</td><td>0.80</td></tr><tr><td>Eq.4</td><td>2.31</td><td>0.80</td><td>2.31</td><td>0.80</td></tr><tr><td>Eq.5</td><td>2.20</td><td>0.81</td><td>2.15</td><td>0.82</td></tr><tr><td>Eq.6</td><td>2.20</td><td>0.82</td><td>2.14</td><td>0.82</td></tr><tr><td>Eq.7</td><td>1.75</td><td>0.88</td><td><strong>1.57</strong></td><td><strong>0.91</strong></td></tr><tr><td>Eq.8</td><td>1.70</td><td>0.89</td><td><strong>1.50</strong></td><td><strong>0.91</strong></td></tr><tr><td>SRH</td><td>1.44</td><td>0.92</td><td><strong>1.33</strong></td><td><strong>0.93</strong></td></tr><tr><td>RFH</td><td>1.55</td><td>0.91</td><td>1.50</td><td>0.91</td></tr></tbody></table>

SM = Ground mode from SM (last detected mode) provided in the L2A data product. HL2M = ground mode corresponding to the higher amplitude between the last two modes.

presented in (9) and (10) were fitted using *in situ* V as the estimated variable and *in situ* $H_{dom}$ as the predictor on the 566 studied stands. Next, to estimate V, the fitted models in (9) and (10) used the estimated $H_{dom}$ values from SRH with GEDI metrics extracted using algorithm a1, while SRV and RFV used the GEDI waveform metrics from Table II extracted with algorithm a1 ($p$ in $RH_n^p$ was set to 1 for RFV). The estimation results of stand volume V (see Fig. 9, Table VIII) show that the four tested models produced similar accuracies on the estimation of V, with an RMSE between 24.39 and 27.45 $m^3 \cdot ha^{-1}$ and coefficient of determination ($R^2$) between 0.87 and 0.90. Moreover, the results in Fig. 9 also show that the estimations of V were close to the 1:1 for all values of V between 0 and 250 $m^3 \cdot ha^{-1}$, while they under estimated the volume for V values higher than 250 $m^3 \cdot ha^{-1}$. For
<table><caption>TABLE VIII
MODELS' PERFORMANCE AND THE FITTED LINEAR EQUATIONS FOR ESTIMATING EUCALYPTUS STAND WOOD VOLUME (V)</caption><thead><tr><th>Model</th><th>RMSE (m<sup>3</sup>.ha<sup>-1</sup>)</th><th>RMSPE (%)</th><th>Bias (m<sup>3</sup>.ha<sup>-1</sup>)</th><th>R<sup>2</sup></th><th>AIC</th></tr></thead><tbody><tr><td>Eq. 9: V = -12.44 + 0.34H<sub>dom</sub><sup>2</sup></td><td>28.19</td><td>21.8</td><td>0.64</td><td>0.87</td><td>36844.2</td></tr><tr><td>Eq. 10: V = 0.17H<sub>dom</sub><sup>2.2</sup></td><td>28.04</td><td>21.9</td><td>0.96</td><td>0.87</td><td>36786.8</td></tr><tr><td>SRV: V = 0.3RH<sub>90</sub><sup>2.2</sup> - 0.34RH<sub>10</sub><sup>2.8</sup> + 0.02RH<sub>30</sub><sup>2.2</sup> - 0.15RH<sub>60</sub><sup>2.2</sup> + 0.03RH<sub>20</sub><sup>2.6</sup></td><td>24.39</td><td>20.4</td><td>-0.03</td><td>0.90</td><td>35467.8</td></tr><tr><td>RFV</td><td>27.48</td><td>21.3</td><td>0.65</td><td>0.88</td><td>36789.0</td></tr></tbody></table>

The variables are described in section II.B.2, with the models described in section III.B

<table><caption>TABLE IX
ACCURACY (RMSE IN M<sup>3</sup>.HA<sup>-1</sup>) OF THE MODELS PRESENTED IN SECTION III-B FOR THE ESTIMATION OF V USING GEDI
METRIC VALUES EXTRACTED USING THE SIX DIFFERENT ALGORITHMS (A1 THROUGH A6)</caption><thead><tr><th rowspan="3">Model</th><th colspan="12">Algorithm</th></tr><tr><th colspan="2">a1</th><th colspan="2">a2</th><th colspan="2">a3</th><th colspan="2">a4</th><th colspan="2">a5</th><th colspan="2">a6</th></tr><tr><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th><th>RMSE</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>Eq. 9</td><td>28.19</td><td>0.87</td><td>28.72</td><td>0.87</td><td>28.59</td><td>0.87</td><td>29.79</td><td>0.86</td><td>29.83</td><td>0.86</td><td>28.12</td><td>0.87</td></tr><tr><td>Eq.10</td><td>28.04</td><td>0.87</td><td>28.53</td><td>0.87</td><td>28.42</td><td>0.87</td><td>29.69</td><td>0.86</td><td>29.82</td><td>0.86</td><td>27.92</td><td>0.88</td></tr><tr><td>SRV</td><td>24.39</td><td>0.90</td><td>27.25</td><td>0.88</td><td>26.32</td><td>0.89</td><td>27.66</td><td>0.88</td><td>28.96</td><td>0.87</td><td>26.85</td><td>0.88</td></tr><tr><td>RFV</td><td>27.48</td><td>0.88</td><td>29.01</td><td>0.87</td><td>31.59</td><td>0.84</td><td>30.87</td><td>0.85</td><td>33.36</td><td>0.82</td><td>29.65</td><td>0.86</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Comparison of measured vs. estimated wood volume from the models presented in section III.B using GEDI metrics extracted with algorithm a1. RMSE is expressed in m<sup>3</sup>.ha<sup>-1</sup>.</figcaption></figure>

the four models, the relative RMSE increased from ~18% for V less than 250 m<sup>3</sup>.ha<sup>-1</sup> to ~40% for V higher than 250 m<sup>3</sup>.ha<sup>-1</sup>. The bias (mean difference of *in situ* V and estimated V) for V higher than 250 m<sup>3</sup>.ha<sup>-1</sup> was also more apparent, and decreased from 2.2 m<sup>3</sup>.ha<sup>-1</sup> (average bias from all models) for V less than 250 m<sup>3</sup>.ha<sup>-1</sup> to 26.5 m<sup>3</sup>.ha<sup>-1</sup> for V higher than 250 m<sup>3</sup>.ha<sup>-1</sup>.

The variable importance test of the GEDI metrics (see Fig. 10) showed that the three most contributing factors on the estimation of V using the random forest regressor (RFV) were the same as those for the estimation of canopy heights, with the highest contributor being RH<sub>90</sub>, followed by RH<sub>80</sub> and RH<sub>100</sub>.

The estimation of V was also tested using GEDI metric values extracted using the remaining five algorithms (a2 through a6). The results presented in Table IX show that the estimates of V were mostly similar with (9) and (10) across all algorithms. On the other hand, SRV and RFV show less accurate estimates of V using GEDI metrics from algorithms a2 through a6 compared to a1.

Finally, the choice of the method for detecting the ground return was also studied. The results, unlike those obtained when estimating H<sub>dom</sub>, did not show any significant improvements when considering the last detected mode, or the stronger between the last two modes.

## V. DISCUSSION

The different tested models in this article showed that GEDI waveform metrics could be used to obtain good accuracies of canopy heights and wood volumes, with a RMSPE of 7.1% on canopy height estimation and 20.4% on wood volume estimation. Moreover, GEDI waveforms appear to be of high quality given the very little variability on the estimation of H<sub>dom</sub> and V from the individual footprints within a given stand. Indeed, the accuracy (RMSE) on the estimation of H<sub>dom</sub> using the mean estimates from SRH of the individual footprints was 1.33 m (R<sup>2</sup> of 0.93) versus 1.32 m (R<sup>2</sup> of 0.93) when averaging H<sub>dom</sub> estimates from GEDI for each stand. Similarly, the accuracy on the estimates of V using the mean estimates from SRV was 24.39 m<sup>3</sup>.ha<sup>-1</sup> (R<sup>2</sup> of 0.90) and 23.93 m<sup>3</sup>.ha<sup>-1</sup> (R<sup>2</sup> of 0.91) for the average of V from GEDI over each stand.

The most important GEDI variable for the estimation of H<sub>dom</sub> and V is RH<sub>90</sub>, which explained respectively more than 87% and 84% of the variability of H<sub>dom</sub> and V. Some of the remaining variability are explained by different GEDI metrics based on the
<figure><img src="image_12.png" /><figcaption>Fig. 10. Classification of the variable importance by decreasing order of importance in the RFV model for stand wood volume estimation. The importance is measured via the average percentage increase of RMSE (%IncMSE) over 50 repetitions. The red bars indicate the standard deviation of %IncMSE.</figcaption></figure>

model used. In the case of the stepwise linear regression models, RH₁₀, RH₈₀, and RH₁₀₀ are used for the estimation of Hₐ₉₆ and RH₁₀, RH₃₀, RH₈₀, and RH₂₀ for the estimation of V.

The success of the models in estimating H<sub>dom</sub> and V with good accuracy from GEDI metrics relies heavily on the accuracy of extracting such metrics from the raw waveform data, and on the precision of the field measurements. GEDI datasets provide metrics issued from six algorithms, and the same metric can differ in value from one algorithm to another. This was evident from our results where the accuracy (RMSE) in estimating both H<sub>dom</sub> and V was slightly less accurate with GEDI metrics generated from algorithms a2 through a6 in comparison to a1. Therefore, the choice of algorithm of which the metrics are calculated from, is important. Nonetheless, GEDI L2A datasets already provide a field for each footprint called “selected_algorithm” that recognizes the algorithm selected as identifying the lowest mode (last detected mode) with less noise. This field could be used as an indicator of the algorithm that provides the most accurate metrics. Indeed, for our dataset, the best accuracy on both H<sub>dom</sub> and V was observed with metrics extracted using algorithm a1, which was also the algorithm suggested by the “selected_algorithm” for more than 99% of the studied footprints.

Another variable that also affects the accuracy of the extracted metrics, such as RH<sub>n</sub>, Trail<sub>ext</sub>, or Lead<sub>ext</sub> is the location of the ground return. In this article, the ground return was identified either by the SM field in the L2a dataset, or by identifying the mode with the higher amplitude between the last two detected modes. Our results, which were obtained over study sites with mostly homogenous canopy cover and flat terrain, indicated that the second method improved the estimation accuracy (RMSE) on H<sub>dom</sub> by up to 20 cm. On the other hand, for the estimation of V, the choice of ground return did not have any effect on the accuracy. However, for highly dense vegetated areas as in the case of tropical forests, where the ground return is not easily identifiable, choosing the strongest mode between the last two as the ground return should provide better accuracies on H<sub>dom</sub> and V [15], [16].

Some of the uncertainties on the estimation of H<sub>dom</sub> and V can be attributed to some biophysical properties of the canopy that could not be quantified using GEDI data alone, or due to instrumental factors. To understand the effect of these biophysical properties, additional parameters could therefore be required.

<figure><img src="image_13.png" /><figcaption>Fig. 11. The effects of adding stand age (in situ information) as a predictor variable on the estimation of (a) stand dominant heights H<sub>dom</sub> (SRHa, m), and (b) stand volume V (SRVa, m<sup>3</sup>.ha<sup>-1</sup>) in a stepwise regression model.</figcaption></figure>

Moreover, since GEDI data were acquired up to two months before or after *in situ* measurements, stand growth that happened between this time gap added to the uncertainty of our estimations (up to 50 cm difference). Canopy age plays an important role in understanding the variability of both H<sub>dom</sub> and V. Fig. 4(a) and (b) show that both H<sub>dom</sub> and V are a log function of stand age. In fact, adding the log transform of stand age to the stepwise linear regression models enhances the accuracy (RMSE) on H<sub>dom</sub> by ~13% (RMSE of 1.33 m without age and 1.18 m with age) and the volume estimation by ~7% (24.39 versus 22.77 m<sup>3</sup>.ha<sup>-1</sup> with age) (see Fig. 11). Moreover, the addition of the age of canopies to the SRH and SRV models also helps reduce the difference between estimates and measurements of some outliers points (see Fig. 11). The interest of using the stand age has also been demonstrated by Le Maire *et al.* [43] in their study over Eucalyptus plantations with MODIS optical data. Other properties of the canopy, such as crown dimension, leaf area index, distribution of leaf angles, tree gaps due to mortality, could also add variability to GEDI waveforms. Therefore, by adequately filtering the stand, one could obtain more accurate estimates of H<sub>dom</sub> and V. On the opposite, the results presented in this article are specific to Eucalyptus plantations, since most sources of canopy variations are included.

The instrumental factors affecting the estimation accuracy include the viewing angle of GEDI at acquisition time. In fact, GEDI acquires data along eight beams with varying viewing angles (VA), which could affect the estimation accuracy. In this article, the analysis of estimation accuracy on both H<sub>dom</sub> and V according to the acquiring beam (see Table X) shows that
<table><caption>TABLE X
COMPARISON OF THE ESTIMATION ACCURACY ON $H_{dom}$ AND V BASED ON THE ACQUIRING BEAM</caption><thead><tr><th rowspan="2">Beam</th><th colspan="2">$H_{dom}$ (m)</th><th colspan="2">V (m³·ha⁻¹)</th></tr><tr><th>RMSE</th><th>R²</th><th>RMSE</th><th>R²</th></tr></thead><tbody><tr><td>1</td><td>1.34</td><td>0.92</td><td>22.92</td><td>0.91</td></tr><tr><td>2</td><td>1.41</td><td>0.93</td><td>25.99</td><td>0.91</td></tr><tr><td>3</td><td>1.22</td><td>0.95</td><td>21.92</td><td>0.91</td></tr><tr><td>4</td><td>1.28</td><td>0.94</td><td>22.08</td><td>0.90</td></tr><tr><td>5</td><td>1.68</td><td>0.90</td><td>30.65</td><td>0.86</td></tr><tr><td>6</td><td>1.39</td><td>0.93</td><td>22.73</td><td>0.91</td></tr><tr><td>7</td><td>1.42</td><td>0.91</td><td>25.30</td><td>0.91</td></tr><tr><td>8</td><td>1.25</td><td>0.93</td><td>22.50</td><td>0.91</td></tr></tbody></table>

Estimation results are produced using stepwise linear regression models with GEDI metric values extracted using algorithm a1.

the most accurate data were from beams 8, 3 and 4, while the least accurate $H_{dom}$ and V estimates were obtained using data from beam 5. Indeed, the difference between the most accurate estimates (beam 3) and the least accurate (beam 5) is 46 cm for $H_{dom}$ and of 9 m³·ha⁻¹ for V. A preliminary analysis of VA values shows that beam 5 acquisitions had on average higher VA values.

A part of the uncertainties on the estimation of $H_{dom}$ and V can be attributed to the heterogeneity of the Eucalyptus stands. In fact, the present study compared GEDI acquisitions with stand-scale averaged $H_{dom}$ and V obtained from 1 to up to 10 permanent inventory plots. These *in situ* measurements are therefore the result of only a few observations within the stand. The hypothesis that the stands are homogeneous enough for using stand-scale averages is sometimes challenged. Indeed, some stands have high variability between their inventory plots, which reaches more than 6 m in some cases. A more precise analysis of intra-stand variability of GEDI waveforms variations could help determine which Eucalyptus stand could be compared to stand-scale values of $H_{dom}$ and V.

Another part of the uncertainties stems from the generalized models used in this article to estimate $H_{dom}$ and V. Indeed, we used a single model across the entire dataset with disregard to the variability of growing conditions over each study site, leading to different canopy structures and stand-scale allometric relationships. A preliminary analysis of locally trained random forest models shows that for the estimation of $H_{dom}$, a locally trained model could slightly improve the estimation results (a maximum RMSE decrease of 16 cm was observed, Table S1). However, for the estimation of V, a locally trained model could reduce volume estimation errors by as much as 12.8 m³·ha⁻¹ (see Table S2). The difference in accuracies on the estimation of V are most probably due to the allometric relations between $H_{dom}$ and V, which could vary greatly between one region and another even for the same tree types, but further analysis is required to confirm these results.

Finally, after the filtering scheme applied to our dataset (see Section II-B-3), the remaining footprints with either value of

the quality_flag (either 0 or 1) showed the same accuracy on the estimation of both $H_{dom}$ and V for all the tested models.

## VI. CONCLUSION

In this article, we analyzed GEDI data in order to determine its accuracy in estimating stand-scale dominant heights ($H_{dom}$) and stand volume (V) of intensively managed Eucalyptus plantations in Brazil. $H_{dom}$ and V values have been estimated using the most accurate models used for estimating forest height and aboveground biomass from ICESat-1 waveforms. The GEDI waveform metrics used in the $H_{dom}$ and V estimation models were extracted using algorithms provided by the land processes distributed active archive center (LP DAAC), in addition to already established metrics for ICESat-1 waveforms. Overall, 5517 GEDI shots over 566 Eucalyptus stands were analyzed over our study area.

For our study site defined by flat and gently sloping terrains (average slope < 5°), six regression models, a stepwise linear regression model (SRH), and a random forest regressor (RFH) using GEDI waveform metrics were assessed on their accuracy to estimate $H_{dom}$. Results showed that the most accurate model was SRH with an RMSE on the $H_{dom}$ estimates of 1.33 m ($R^2$ of 0.93) using metrics extracted with the configuration of algorithm a1 and considering the higher mode between the last two modes as the ground return. For this model, the most relevant metrics for the estimation of $H_{dom}$ was the 90th percentile relative height (RH$_{90}$), followed by RH$_{10}$, RH$_{80}$, and RH$_{100}$.

Stand wood volume (V) was modeled following power law with the canopy height, a stepwise regression model (SRV), and a random forest regressor (RFV). The four tested models showed similar accuracies, with SRV being the most accurate one with an RMSE of 24.39 m³·ha⁻¹ (relative error ~20% of the wood volume average). Similar to SRH, the most relevant metrics for the estimation of V using SRV was RH$_{90}$, followed by RH$_{10}$, RH$_{30}$, RH$_{80}$, and RH$_{20}$.

The choice of the algorithm used to extract the waveform metrics affected sometimes the accuracy, as metrics extracted using algorithm a5 showed ~16% higher RMSE on the estimation of both $H_{dom}$ and V. Nonetheless, the field “selected_algorithm” from the L2A dataset provides a robust indicator of the algorithm that provides the most accurate metrics. In addition to the accuracy of the metrics on the accuracy of $H_{dom}$, selecting the higher of the last two modes as the ground return could potentially increase the accuracies on $H_{dom}$ from 12 to 20 cm.

Not all the variability on $H_{dom}$ and V could be explained by GEDI alone. For example, including the stand age in the stepwise regression models help decrease the RMSE on the estimation of $H_{dom}$ and V by, respectively, ~13% and ~7%. Nonetheless, and despite our efforts, underestimation of V using GEDI data was still observed for V greater than 250 m³·ha⁻¹.

Finally, given the high accuracy of GEDI data on the estimation of tree heights and volume, GEDI can provide an excellent source of information to calibrate and validate upcoming and future radar missions, such as the upcoming P-band BIOMASS mission. GEDI can also supplement radar data by means of data fusion models in order to obtain high resolution and very accurate wall-to-wall maps of forest properties.
## AUTHOR CONTRIBUTIONS

Ibrahim Fayad—Conceptualization, Methodology, software, validation, formal analysis, data curation, visualization, writing—original draft. Nicolas Baghdadi—conceptualization, methodology, validation, formal analysis, data curation, writing—original draft. Clayton Alcarde Alvares—Conceptualization, validation, writing—review and editing. Jose Luiz Stape—Conceptualization, validation, writing—review and editing. Jean Stéphane Bailly—Validation, writing—review and editing. Henrique Ferraço Scolforo—Conceptualization, validation, writing—review and editing. Mehrez Zribi—Validation, writing—review and editing. Gueric Le Maire—Conceptualization, validation, writing—review and editing.

## ACKNOWLEDGMENT

The authors would like to thank the GEDI team and the NASA LPDAAC (Land Processes Distributed Active Archive Center) for providing GEDI data. The authors acknowledge Suzano's researchers Italo Ramos Cegatta, Renan Tarenta Meirelles Brasil and Carla Foster Feria for their technical support and the CIRAD Suzano project. Suzano SA Company supported the forest-field data collection.

## REFERENCES

[1] M. Main-Knorn et al., "Monitoring coniferous forest biomass change using a Landsat trajectory-based approach," *Remote Sens. Environ.*, vol. 139, pp. 277-290, Dec. 2013.

[2] A. Peregon and Y. Yamagata, "The use of ALOS/PALSAR backscatter to estimate above-ground forest biomass: A case study in Western Siberia," *Remote Sens. Environ.*, vol. 137, pp. 139-146, Oct. 2013.

[3] R. Houghton, F. Hall, and S. J. Goetz, "Importance of biomass in the global carbon cycle," *J. Geophys. Res., Biogeosci.*, vol. 114, no. G2, 2009.

[4] T. E. Fatoyinbo and M. Simard, "Height and biomass of mangroves in Africa from ICESat/GLAS and SRTM," *Int. J. Remote Sens.*, vol. 34, no. 2, pp. 668-681, Jan. 2013.

[5] M. A. Lefsky et al., "Estimates of forest canopy height and aboveground biomass using ICESat: ICESAT Estimates of Canopy Height," *Geophys. Res. Lett.*, vol. 32, no. 22, Nov. 2005.

[6] T. Payn et al., "Changes in planted forests and future global implications," *Forest Ecol. Manage.*, vol. 352, pp. 57-67, 2015.

[7] P. Elias and D. Boucher, *Planting for the Future. How Demand for Wood Products Could Be Friendly to Tropical Forests*. Cambridge, MA, USA: Union Concerned Scientists, 2014.

[8] R. Pirard, L. Dal Secco, and R. Warman, "Do timber plantations contribute to forest conservation?," *Environ. Sci. Policy*, vol. 57, pp. 122-130, Mar. 2016.

[9] M. A. Lefsky, W. B. Cohen, G. G. Parker, and D. J. Harding, "Lidar remote sensing for ecosystem studies," *Bio. Sci.*, vol. 52, no. 1, pp. 19-22, 2002.

[10] R. Nelson, K. J. Ranson, G. Sun, D. S. Kimes, V. Kharuk, and P. Montesano, "Estimating Siberian timber volume using MODIS and ICESat/GLAS," *Remote Sens. Environ.*, vol. 113, no. 3, pp. 691-701, Mar. 2009.

[11] R. O. Dubayah et al., "Estimation of tropical forest height and biomass dynamics using Lidar remote sensing at La Selva, Costa Rica: Forest dynamics using LiDAR," *J. Geophys. Res.*, vol. 115, Jun. 2010.

[12] B. E. Schutz, H. J. Zwally, C. A. Shuman, D. Hancock, and J. P. DiMarzio, "Overview of the ICESat mission," *Geophys. Res. Lett.*, vol. 32, no. 21, 2005, Art. no. L21S01.

[13] Y. Pang, M. Lefsky, G. Sun, and J. Ranson, "Impact of footprint diameter and off-nadir pointing on the precision of canopy height estimates from spaceborne LiDAR," *Remote Sens. Environ.*, vol. 115, no. 11, pp. 2798-2809, Nov. 2011.

[14] I. Fayad et al., "Aboveground biomass mapping in French Guiana by combining remote sensing, forest inventories and environmental data," *Int. J. Appl. Earth Observ. Geoinf.*, vol. 52, pp. 502-514, Oct. 2016.

[15] Q. Chen, "Retrieving vegetation height of forests and woodlands over mountainous areas in the pacific coast region using satellite laser altimetry," *Remote Sens. Environ.*, vol. 114, no. 7, pp. 1610-1627, Jul. 2010.

[16] N. Baghdadi et al., "Testing different methods of forest height and aboveground biomass estimations from ICESat/GLAS data in eucalyptus plantations in Brazil," *IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens.*, vol. 7, no. 1, pp. 290-299, Jan. 2014, doi: 10.1109/JSTARS.2013.2261978.

[17] J. Boudreau, R. Nelson, H. Margolis, A. Beaudoin, L. Guindon, and D. Kimes, "Regional aboveground forest biomass using airborne and spaceborne LiDAR in Québec," *Remote Sens. Environ.*, vol. 112, no. 10, pp. 3876-3890, Oct. 2008.

[18] M. El Hajj, N. Baghdadi, N. Labrière, J.-S. Bailly, and L. Villard, "Mapping of aboveground biomass in Gabon," *Comptes Rendus Geosci.*, vol. 351, no. 4, pp. 321-331, Apr. 2019.

[19] M. R. Pourrahmati et al., "Capability of GLAS/ICESat data to estimate forest canopy height and volume in mountainous forests of Iran," *IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens.*, vol. 8, no. 11, pp. 5246-5261, Nov. 2015, doi: 10.1109/JSTARS.2015.2478478.

[20] M. R. Pourrahmati et al., "Mapping Lorey's height over Hyrcanian forests of Iran using synergy of ICESat/GLAS and optical images," *Eur. J. Remote Sens.*, vol. 51, no. 1, pp. 100-115, Jan. 2018.

[21] A. Neuenschwander and K. Pitts, "The ATL08 land and vegetation product for the ICESat-2 mission," *Remote Sens. Environ.*, vol. 221, pp. 247-259, Feb. 2019.

[22] R. Dubayah et al., "The global ecosystem dynamics investigation: High-resolution laser ranging of the earth's forests and topography," *Sci. Remote Sens.*, vol. 1, Jun. 2020, Art. no. 100002.

[23] A. Persson, J. Holmgren, and U. Soderman, "Detecting and measuring individual trees using an airborne laser scanner," *Photogramm. Eng. Remote Sens.*, vol. 68, no. 9, pp. 925-932, 2002.

[24] S. Saarela et al., "Hierarchical model-based inference for forest inventory utilizing three sources of information," *Ann. Forest Sci.*, vol. 73, no. 4, pp. 895-910, Dec. 2016.

[25] S. Holm, R. Nelson, and G. Ståhl, "Hybrid three-phase estimators for large-area forest inventory using ground plots, airborne LiDAR, and space LiDAR," *Remote Sens. Environ.*, vol. 197, pp. 85-97, Aug. 2017.

[26] T. R. Feldpausch et al., "Tree height integrated into pantropical forest biomass estimates," *Biogeosciences*, vol. 9, no. 8, pp. 3381-3403, Aug. 2012.

[27] E. Kearsley et al., "Conventional tree height-diameter relationships significantly overestimate aboveground carbon stocks in the central Congo Basin," *Nat. Commun.*, vol. 4, no. 1, Oct. 2013.

[28] Y. Su, Q. Ma, and Q. Guo, "Fine-resolution forest tree height estimation across the Sierra Nevada through the integration of spaceborne LiDAR, airborne LiDAR, and optical imagery," *Int. J. Digit. Earth*, vol. 10, no. 3, pp. 307-323, Mar. 2017.

[29] H. Tang et al., "Deriving and validating leaf area index (LAI) at multiple spatial scales through lidar remote sensing: A case study in Sierra National Forest, CA," *Remote Sens. Environ.*, vol. 143, pp. 131-141, Mar. 2014.

[30] Y. Wang et al., "Is field-measured tree height as reliable as believed – A comparison study of tree height estimates from field measurement, airborne laser scanning and terrestrial laser scanning in a Boreal forest," *ISPRS J. Photogramm. Remote Sens.*, vol. 147, pp. 132-145, Jan. 2019.

[31] R. Dubayah et al., "The global ecosystem dynamics investigation: High-resolution laser ranging of the earth's forests and topography," *Sci. Remote Sens.*, vol. 1, Jun. 2020, Art. no. 100002.

[32] S. L. R. Dubayah, "GEDI L1B geolocated waveform data global footprint level V001," NASA EOSDIS Land Processes DAAC, 2020. Accessed: Jul. 2021. [Online]. Avilable: https://doi.org/10.5067/GEDI012_B.001

[33] S. L. R. Dubayah, "GEDI L2A elevation and height metrics data global footprint level V001," NASA EOSDIS Land Processes DAAC, 2020. Accessed: Jul. 2021. [Online]. Avilable: https://doi.org/10.5067/GEDI/GEDI02_A.001

[34] S. L. R. Dubayah, "GEDI L2B canopy cover and vertical profile metrics data global footprint level V001," NASA EOSDIS Land Processes DAAC, 2020. Accessed: Jul. 2021. [Online]. Avilable: https://doi.org/10.5067/GEDI/GEDI02_B.001

[35] C. Hilbert and C. Schmullius, "Influence of surface topography on ICESat/GLAS forest height estimation and waveform shape," *Remote Sens.*, vol. 4, no. 8, pp. 2210-2235, Jul. 2012.

[36] T. J. Urban, B. E. Schutz, and A. L. Neuenschwander, "A survey of ICESat coastal altimetry applications: Continental coast, open ocean island, and inland river," *Terr. Atmospheric Ocean. Sci.*, vol. 19, pp. 1-19, 2008.
[37] D. J. Harding, “ICESat waveform measurements of within-footprint topographic relief and vegetation vertical structure,” *Geophys. Res. Lett.*, vol. 32, no. 21, 2005, Art. no. L21S10.

[38] M. A. Lefsky, M. Keller, Y. Pang, P. B. De Camargo, and M. O. Hunter, “Revised method for forest canopy height estimation from geoscience laser altimeter system waveforms,” *J. Appl. Remote Sens.*, vol. 1, no. 1, 2007, Art. no. 013537

[39] Y. Pang, M. Lefsky, H.-E. Andersen, M. E. Miller, and K. Sherrill, “Validation of the ICEsat vegetation product using crown-area-weighted mean height derived using crown delineation with discrete return lidar data,” *Can. J. Remote Sens.*, vol. 34, pp. S471–S484, 2008.

[40] L. Breiman, “Random forests,” *Mach. Learn.*, vol. 45, no. 1, pp. 5–32, 2001.

[41] S. S. Saatchi et al., “Benchmark map of forest carbon stocks in tropical regions across three continents,” *Proc. Nat. Acad. Sci.*, vol. 108, no. 24, pp. 9899–9904, Jun. 2011.

[42] H. Akaike, “Information theory and an extension of the maximum likelihood principle,” in *Selected Papers of Hirotugu Akaike*, E. Parzen, K. Tanabe, and G. Kitagawa, Eds. New York, NY, USA: Springer, 1998, pp. 199–213.

[43] G. le Maire et al., “MODIS NDVI time-series allow the monitoring of eucalyptus plantation biomass,” *Remote Sens. Environ.*, vol. 115, no. 10, pp. 2613–2625, Oct. 2011.

<figure><img src="image_14.png" /></figure>

**Ibrahim Fayad** received the Engineering degree in computer and telecommunications in 2011 and the Ph.D. degree in automatic and microelectronic systems both from the University of Montpellier, Montpellier, France in 2015.

He is currently a Research Engineer with the National Research Institute for Agriculture, Food and the Environment, Montpellier, France. His research interests include machine learning for the retrieval of environmental parameters using remote sensing data.

<figure><img src="image_15.png" /></figure>

**Nicolas N. Baghdadi** received the Ph.D. degree from the University of Toulon, Toulon, France in 1994.

From 1995 to 1997, he was a Postdoctoral Researcher with INRS Ete—Water Earth Environment Research Centre, Quebec University. From 1998 to 2008, he was with the French Geological Survey, Orléans, France. Since 2008, he has been a Senior Scientist with the National Research Institute for Agriculture, Food and the Environment, Montpellier, France. He is currently a Scientific Director with the THEIA Land Data Centre.

His research activities are in the areas of microwave remote sensing, image processing, and satellite and airborne remote sensing data analysis. His main field of interest is the analysis of remote sensing data (SAR, Lidar, optical) and the retrieval of environmental parameters (e.g., soil moisture content, surface roughness, biomass, etc.).

<figure><img src="image_16.png" /></figure>

**Clayton Alcarde Alvares** received the Ph.D. degree in forestry science from the University of São Paulo, São Paulo, Brazil, in 2012.

His research focused on mapping and edaphological modeling of productivity of Eucalyptus plantations in Brazil. Since 2011, he has been a Forest Scientist with Suzano SA, São Paulo, Brazil, where he has been developing in-depth analyzes in climatology, remote sensing, production environments zoning, and applied ecophysiology, and is strongly committed to delivering tools for operational uses in the short,

medium and long term. Since 2019, he has been a Professor and an Advisor of the Postgraduate Program in Forest Science, State University "Júlio Mesquita Filho"—UNESP, Faculty of Agronomic Sciences—Campus de Botucatu.

<figure><img src="image_17.png" /></figure>

**Jose Luiz Stape** received the Ph.D. degree in forest ecology from Colorado State University, Fort Collins, CO, USA, in 2002.

He is a Permanent Graduate Professor of Forest Ecophysiology with Sao Paulo State University (UNESP, Brazil), Sao Paulo, Brazil. He was with the University of Sao Paulo, and with the North Carolina State University and across many countries and companies, looking to improve silvicultural recommendations for the sustainability of forest plantations including: clonal deployment; site-preparation; nutrition and spacing. To better evaluate the factors limiting forest productivity and controlling C allocation, he coordinated the establishment, with other scientists, of four large Eucalyptus and Pine cooperative research programs in Brazil via IPEF (BEPP, Eucflux, TECHS and PPIB) and a research network at Suzano company (G2M2P2). Nowadays the use of remote sensing to improve monitoring, management and modeling planted forests has been a main focus of his research.

<figure><img src="image_18.png" /></figure>

**Jean Stéphane Bailly** received the Engineering degree in agronomy, the M.Sc. degree in biostatistics, and the Ph.D. degree in hydrology from the University of Montpellier, Montpellier, France, in 1990, 2003, and 2007.

He is currently a Senior Lecturer of physical geography and geostatistics in AgroParis-Tech, Paris, France. He is currently a Research Fellow with the UM LISAHLab, Montpellier, France. His research interests include spatial observations and modeling for hydrological issues.

<figure><img src="image_19.png" /></figure>

**Henrique Ferraço Scolforo** received the Ph.D. degree in forest biometrics from the North Carolina State University, Raleigh, NC, USA, in 2018.

His research focused on growth and yield modeling sensitive to climate and clonal variation applied to eucalypt stands in Brazil. Since 2018, he has been leading the biometrics, inventory, growth, and yield studies with Suzano SA, São Paulo, Brazil.

<figure><img src="image_20.png" /></figure>

**Mehrez Zribi** received the B.E. degree in signal processing from the Ecole Nationale Supérieure d'Ingénieurs en Constructions Aéronautiques, Toulouse, France, in 1995, and the Ph.D. degree in signal processing and remote sensing from the Université Paul Sabatier, Toulouse, France, in 1998.

He is currently a Research Director with the Centre National de Recherche Scientifique, Paris, France. In 1995, he was with the Centre d'Etude des Environnements Terrestre et Planétaires Laboratory/Institut Pierre Simon Laplace, Vélizy, France. Since October 2008, he has been with the Centre d'Etudes Spatiales de la Biosphère, Toulouse, France. He has authored or coauthored more than 140 articles in refereed journals. He is currently the Director of Centre d'Etudes Spatiales de la Biosphère, Toulouse, France. His research interests include microwave remote sensing applied to hydrology, microwave modeling for land surface parameters estimations and finally airborne microwave instrumentation.

<figure><img src="image_21.png" /></figure>

**Guerric le Maire** received the M.Sc. degree in agronomy, specialized in ecology, from the National Institute of Agronomy (INA P-G, Paris, France), Paris, France, in 2002 and the Ph.D. degree in plant ecology from Paris XI University, Orsay, France, in 2005.

During 2006–2007, he was a Postdoctoral Researcher with the Le Laboratoire des Sciences du Climat et de l’Environnement Laboratory (Saclay-France). He was a Researcher with the CIRAD, Montpellier, France, in 2008. His research interests include remote-sensing image processing/analysis and process-based ecophysiological forest models development applied to tropical forest plantations.