Modelling Canopy Height of Forest-Savannah Mosaics in Togo Using ICESat-2 and GEDI Spaceborne LiDAR and Multisource Satellite Data

1. Introduction

The rapid rise in global greenhouse gas (GHG) emissions since the debut of the Industrial Revolution has led to considerable changes in the Earth’s climate, which has been the subject of much research over recent decades. The Intergovernmental Panel on Climate Change (IPCC) has determined that this increase in anthropogenic GHG emissions is the primary driver of climate change, which can push temperatures beyond the thermal tolerances of many species [1,2]. According to current trends, global warming that is linked to these gases will likely exceed 1.5 °C within several decades, despite aggressive emissions reduction strategies [3]. Removing CO2 from the atmosphere by favouring nature-based solutions (protected areas and forests) could contribute greatly to climate change mitigation [4,5,6].

Forests constitute one of the largest reservoirs of terrestrial carbon, thereby playing a vital role in offsetting the aforementioned climate changes and regulating global carbon balance. Annually, they contribute to about 50% of net terrestrial primary production, store about 45% of the planet's active carbon, and sequester around 33% of anthropogenic emissions [7,8,9]. Unfortunately, tropical forest ecosystems, which maintain the global ecological equilibrium, are constantly threatened by deforestation and degradation for economic purposes. Timber extraction and forest clearing for other land uses are among the largest sources of anthropogenic carbon emissions [10]. In this context, the implementation of the Paris Agreement on climate change and the 2030 Agenda for Sustainable Development (adopted in 2015) by UN member nations require quantitative studies that would provide essential data for monitoring forest dynamics [11,12,13].

Canopy height is one of several structural parameters that is used to monitor these forest dynamics. Its precise estimation is crucial for quantifying biophysical parameters, such as aboveground biomass, carbon storage, biodiversity, and many other parameters to which it is strongly linked [14,15,16]. Traditional methods of estimating forest height are based upon manual forest inventories, which are conducted in the field. While these inventories can provide accurate detailed information, they are very labour-intensive, time-consuming, and may cover only small spatio-temporal scales. Consequently, remote sensing has been employed in recent decades, in combination with field measurements, to estimate canopy height over large spatial extents [17,16]. The remote sensing observations are provided by various platforms. These include for instance multi-spectral optical data that are derived from Landsat [18,19], Sentinel 2 [20,21], or SPOT5 [22,23], as well as synthetic aperture radar (SAR) data from Sentinel 1 [24,25], TerraSAR-X [26,27], TanDEM-X [28,29] or ALOS PALSAR [30,31]. Regardless of whether the data are optical or obtained from radar, signal saturation (especially in dense forests) can constitute a substantial limitation [32,33,34,35]. The introduction of LiDAR (Light Detection And Ranging) has enabled notable advances in the estimation of canopy height, because of their capacity to detect vertical structure in the forest [16]. The most frequently used applications in forestry are based on telemetry from airborne laser scanning (ALS) [36,37] and terrestrial laser scanning (TLS) [38,39]. However, ALS and TLS exhibit spatio-temporal limitations, given that it is generally difficult to apply them over large areas and on a regular basis, due to their high costs of acquisition and to signal occultation, particularly in dense forests [40,41,42,43].

Airborne or satellite platforms have made it possible to extend canopy height estimation from local to global spatial scales [13]. The first platform, i.e., Ice, Cloud, and land Elevation Satellite (ICESat) carried the Geoscience Laser Altimeter System (GLAS). Between 2003 and 2009, this sensor made it possible to estimate the height of forests on a global scale on circular footprints with a diameter of about 60 m and a spacing of about 170 m along the transects [44,45,46,9]. The second ICESat-2 satellite was launched in September 2018 and carried the Advanced Topographic Laser Altimetry System (ATLAS), which uses photon-counting LiDAR technology. Between 88◦ S and 88◦ N, the laser produces three pairs of beams, thereby making it possible to obtain altimeter parameters on Earth’s surface for continuous monitoring of polar glaciers. In parallel with its main mission, this satellite also acquired terrestrial measurements of forest cover and vegetation. For the scientific community, this represents an important database for mapping plant biomass and for estimating carbon inventories at a global level [47,17]. The latest LiDAR instrument that was launched into space by NASA in December 2018 is the Global Ecosystem Dynamics Investigation (GEDI) system, which operates aboard the International Space Station (ISS) between 51.6◦ N and 51.6◦ S. GEDI is a multi-beam laser altimeter that measures parameters of vertical canopy structures at a very high sampling rate, thereby allowing forest height and wood volume estimation across different types of forest ecosystems, topography and latitudes [13]. Data that are acquired by GEDI have been increasingly used to estimate forest height and forest biomass [48,49,50]. These data consist of an impressive number of samples, which offer great potential for estimating canopy heights in complex savannah and forest mosaics.

ICESat-2 and GEDI LiDAR data are point clouds that permit the height of forest cover to be estimated, but only within ground acquisition footprints, rather than in a spatially continuous manner over large areas [51]. In contrast, optical or radar data provide continuous spatial coverage; but cannot, alone, allow the direct extraction of vertical profiles of the canopy. Therefore, the complementarity of different types of data can be exploited to map the height of the forest cover [52,53,54,55]. To accomplish this task, machine-learning models are being increasingly used to combine vertical LiDAR profiles with spectral or backscatter attributes [56,16]. For example, Li et al. [15] used variables derived from Sentinel 1&2 and Landsat-8 images over Northeast China to extrapolate ICESat-2 canopy height from the footprint-level to regional-level, using Deep Learning (DL) and Random Forest (RF) models, with correlations of (r =) 0.78 and 0.68, respectively. Zhu et al. [9] used stepwise regression and Random Forest (RF) approaches to estimate canopy heights in the United States. They obtained better results with RF using GEDI variables (R² = 0.93; RMSE = 2.99 m) compared to those of ICESat-2 (R² = 0.78; RMSE = 4.62 m). Sothe et al. [57] carried out continuous mapping of the forest cover height of Canada from the combination of GEDI and ICESat-2 data with PALSAR and Sentinel data. They found that both LiDAR products overestimated canopy height compared to ALS data, but GEDI outperformed ICESat-2, with an average difference of 0.9 m vs. 2.9 m and RMSE of 4.2 m vs. 5.2 m, respectively. To map China's forest canopy heights, Liu et al. [58] used neural network-guided interpolation to merge GEDI and ICESat-2 data. They then compared the height of the forest cover that was interpolated with the GEDI validation footprints (R² = 0.55; RMSE = 5.32 m), followed by drone-LiDAR validation data (R² = 0.58, RMSE = 4.93 m) and, finally, with field-collected data (R² = 0.60; RMSE = 4.88 m).

The aforementioned examples show real potential for using GEDI and ICESat-2 data alone or in combination with other spatial data. However, they also raise several questions, which depend upon the ecosystems that are being considered. Of particular interest is the following: Can the use of GEDI or ICESat-2 data (alone or in combination) with multisource optical or radar satellite observations make it possible to estimate satisfactorily the canopy heights of complex mosaics of forests and savannahs in a tropical environment? Our study attempts to answer the question through analyses of the forest-savannah mosaics of the Sudano-Guinean zone of West Africa, particularly in Togo, where research of this type is almost non-existent. Our main objective is to develop models for estimating height of the canopy in forest-savannah mosaics using a combination of space LiDAR, optical data and radar data. The specific objectives that we pursued are: 1) to analyze the performance of covariates and ICESat-2 and GEDI data in predicting canopy height in these forest ecosystems; 2) to develop canopy height prediction models that are adapted to forest-savannah mosaics; and 3) to establish continuous mapping of canopy height in these forest types from these discontinuous satellite LiDAR data. To accomplish these goals, optical and radar co-variables, such as spectral reflectances, vegetation indices, texture and backscatter variables, are derived from the spatially continuous satellite data, which were then integrated with those derived from ICESat-2 and GEDI, using machine-learning models.

2. Materials and Methods

2.1. Study Aera

Our study was conducted in Ecological Zone 4, southwest Togo. The Togolese Republic is a coastal nation in West Africa that is bordered on the north by Burkina Faso, by the Atlantic Ocean (Gulf of Guinea) to the south, by Benin to the east, and by Ghana to the west. The country is subject to a tropical Sudano-Guinean climate, with rainfall varying across four seasons from 1000 to 1600 mm/year in the southern regions. Average temperature is 27 °C [59]. Between 1939 and 1957, the classification of 14.2% of the country’s land area into protected areas (classified forests, national parks and reserves) served to preserve its forest cover. Today, many of these areas have been encroached by human populations seeking arable land and wood for energy. Vegetation types are composed of Sudano-Guinean forests, which are located in mountainous areas of the country, gallery forests along the main waterways, dry forests or dense tree savannahs in the arid northern half, and tree savannahs in the south and center. Five ecological zones have been designated, spanning the country.

Located in the southern Togo Mountains, Ecological Zone 4 (6397 km2), is one of such subdivision that is characterized by the landscape variability of its ecosystems [60]. Always known as the most heavily forested of the country’s ecological zones, Zone 4 is dominated by interspersed semi-deciduous forest and mosaics of Guinean savannah, with the latter having been degraded in recent years by the combined effects of slash-and-burn agriculture, vegetation fire and logging [61,62]. Over the past three decades, the zone has lost more than 27% of its forest cover, yet it remains the most forested and least degraded of the country's five ecological zones [63]. The study area is presented in Figure 1.

2.2. Methodology

To develop canopy height prediction models, multisource variables were extracted from the collected data, and made usable during models’ development process. Modelling involves the use of relative canopy heights and other variables that have been extracted from satellite LiDAR data (GEDI and ICESat-2), together with those that have been extracted from spatially continuous data. After the preprocessing stage, prediction variables were extracted from radar (Sentinel 1), optical (Sentinel 2) and topographical (SRTM) data. The former consist of native bands, vegetation indices, texture measurements and topographical variables. The methodological flowchart of this study using these data is illustrated in Figure 2.

2.3. Data Acquisition

The data that were collected and organized prior to use in this research come from seven sources (see Figure 2), separated into two main categories. The first category concerns remote sensing data, notably optical, radar, topographical and satellite LiDAR sources; the second category consists of dendrometric field data. The general parameters of these data sources data are summarized in Table 1, while the variables that are extracted and their extraction methods are described in Section 2.4.

Spatially discontinuous data from the Global Ecosystem Dynamics Investigation (GEDI) and Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) satellites have been downloaded from NASA's Land Processes Distributed Active Archive Center (LPDAAC) website (https://lpdaac.usgs.gov/). It should be noted that GEDI data are products of the International Space Station (ISS). The granules are GEDI data reduced size from one full ISS orbit to four segments per orbit [66]. The structures of the ground footprints of these satellite LiDAR data are illustrated in Figure 3a and Figure 3b, for ICESAt-2 and GEDI, respectively. Figure A1 (Appendix A) shows elevations of the ground surface (a) and canopy top (b) extracted from photon returns in ATL08 data acquired by ICESat-2 in the southwest of the city of Badou (7°35′8″N, 0°36′33″E). Continuous Sentinel 1, Sentinel 2 and SRTM data were collected from archives of the Google Earth Engine platform.

Regarding field data, a field campaign that was conducted between October 2020 and February 2021 made it possible to collect dendrometric field parameters within the rectangular footprints (17 m by 100 m) of the ICESat-2 data. Within these footprints, the total height was measured on the footprints using a Suunto clinometer and the diameter at breast height (DBH, 1.3 m) was measured using a direct diameter measurement tape, for all trees with a diameter greater than or equal to 10 cm. To supplement these data, we added measurements taken from the second National Forest Inventory (IFN2). An existing land use map from 2020 was also used in this study as auxiliary data to identify the most forested areas.

2.4. Variable Extraction

From Sentinel 1, Sentinel 2 and SRTM data, variables were extracted or calculated to serve as independent variables when predicting canopy height. The latter are made up of native bands, vegetation indices, texture measurements and topographical variables. They included 29 variables for radar data, 28 for optical data, and 3 for topographical data. We used JavaScript code when extracting these continuous variables. Table 2 contains the summary list of all variables that were resampled at 30 m resolution, particularly those variables which were not already measured at 30 m, so that they could be used with the other covariates.

Formulas used for calculating some variables, together with their respective references, are provided in Table A1 (Appendix B). From the GEDI L2A granules of the GEDI data, we extracted location parameters (latitude, longitude) and relative height values as percentiles (RH50, RH55, RH60, RH65, RH70, RH80, RH85, RH90, RH95, RH98 and RH100), which are frequently used for canopy height prediction. Other variables that are considered useful for modelling forest structure, such as beam type (coverage, full power), data quality indicator and sensitivity of the waveform to penetrate vegetation, were extracted. For the ICESat-2 data, we downloaded ATL08 products from which we extracted the powerful beams, representing the best option for detecting ground and canopy photons according to Neuenschwander and Pitts [67]. We used Matlab code when extracting these ICESat-2 and GEDI variables. Relative heights were also calculated for each of the plots, given that the field measurements were taken from individual trees.

2.5. Preparation of Variables for Modelling

Since the variables extracted from the data could not be used directly to develop the canopy height models, they were prepared for that purpose. This involved validating the relative heights, filtering the data and calculating zonal statistics.

2.5.1. Validation of Satellite LiDAR Data

In order to use the relative heights extracted from ICESat-2 data as predictors of canopy height, we validated them with those that were collected from their ground footprints. To this end, they were compared with those calculated from field data. Given that the field data were collected on individual trees, we first aggregated them into data which could be easily compared to those extracted from ICESat-2, by calculating their statistics per plot. For each plot, the statistics calculated from the field data are: the minimum, the maximum, the first quartile, the mean, the median, and relative heights of 50%, 55%, 60%, 65%, 70%, 75 %, 80%, 85%, 90%, 95% and 98%.

The variable that was extracted from the ICESat-2 data and used for this comparison with field data is the relative height at 98%, designated h_canopy. Once these field variables had been calculated, the latter were matched with them by plot in an Excel table, which was used to establish a correlation matrix. The Pearson correlation values obtained between the field variables and relative height, h_canopy, allowed to use it as reference data for modelling canopy height. It should be noted that this validation using field data only concerns ICESat-2, given that we did not have field data from the GEDI footprints to conduct a similar analysis.

2.5.2. Data Filtering

During variable preparation, the data were filtered using various parameters internal to the dataset. During extraction of relative heights, other variables were extracted from satellite LiDAR and auxiliary data to use them as filtering parameters to clean the database used for modelling.

• For the ICESat-2 data, the number of photons classified as canopy in the segments (n_ca_photons) and the spacecraft orientation tracking parameter (sc_orient) were considered. During filtering, the n_ca_photons parameter was used to eliminate all footprints with a low number of canopy photons, while the sc_orient parameter was used to identify the strong beams from the weak beams: When sc_orient is equal to 0, the satellite is ascending and the strong beams are on the left on the ground tracks, while the weak beams are on the right, whereas the opposite occurs when it is equal to 1.
• For GEDI data, the parameters used are the identifiers for the type of beam which can be coverage or full power (Beam), data quality indicator (Quality_flag), sensitivity of the waveform to penetrate vegetation (Sensitivity) and firing time (delta_time). Filtering using these parameters made it possible to select footprints according to the nine configurations that are presented in Table 3.

• The geographical coordinates extracted from the ICESat-2 and GEDI data were used to display them on a GIS land cover map to identify the land cover classes into which each of their footprints fell. This action made it possible to select footprints that fell into the forest classes, and to remove from the database those that fell into other classes such as "crops and fallow land", "buildings and bare soil" and "grassy savannah".

The data filtering stage therefore made it possible to remove from the database any values that could have added noise to the modelling process or which reduced the quality or performance of the prediction models being developed. Filtering in relation to the location of LiDAR data footprints in forest areas was conducted in ArcGIS 10.8.1. Filtering of other parameters that were extracted from LiDAR data was conducted using a Python code developed for this purpose.

2.5.3. Calculation of Zonal Statistics

To properly integrate the data during modelling, zonal statistics were calculated on the continuous variables from the ICESat-2 and GEDI data footprints. For each variable extracted from the Sentinel 1, Sentinel 2 and SRTM data, the mean, median and standard deviation of the pixel values were calculated for each footprint of the satellite LiDAR data superimposed upon it. These three calculated statistics were evaluated at the start of modelling to select the one that provided the most accurate canopy height prediction results. Figure 4 shows an example of the GEDI and ICESat-2 footprints superposed over Sentinel 1 for the calculation of zonal statistics.

For each of these variables, the statistics were matched with relative heights and other variables extracted from the satellite LiDAR data to form a table for each variable. These data tables per variable were then grouped by data source (optical, radar, topographic) to provide the databases used in modelling.

2.6. Modelling

Once the variables had been extracted and prepared, we proceeded with modelling, which consisted of variable selections, predictive model development, and evaluation. Given that these continuous variables provide different information, we analyzed different scenarios for combining them, as presented in Table 4, to make the best possible choices.

To explore the influence of height classes on modelling, we divided the data into two groups of height classes. Because of their small number, the data with canopy heights ≤ 5 m (]2 - 5]) formed the first class in these two groups. Group 1 contained 7 classes increasing in 5 m increments from the upper bound of the first class (i.e., (]2 - 5]), (]5 - 10]), (]10 - 15]), etc.). Group 2 contained 10 height classes in which the step size or increment was set at 3 m relative to the upper bound of the first class (i.e. (]2 - 5]), (]5 - 8]), (]8 - 11]), etc.). Like the first class, canopy heights > 30 m (]30 – 50]) constituted a separate class in each of the two groups. When data were filtered (see Section 2.5.2), individuals with heights < 2 m or > 50 m were deleted from the database with reference to the field data.

The database created in the previous step is such that $D=\left\{\left(x_1,h_1\right),...,\left(x_N,h_N\right)\right\}$, where $x_i=\left(x_{i,1},...,x_{i,P}\right)$, $N$ is the set of sampled GEDI or ICESat-2 footprints, $P\in R$ represents the set of extracted attributes, and $h_N$ represents the canopy height value that was extracted from the $N$observations of GEDI or ICESat-2 footprints. With these input data ($D$), one part ($D_T$ = 80%) trained the models, and the other part ($D_t$ = 20%) tested these models to evaluate their performance.

2.6.1. Variable Selection

In this study, the importance of variables was determined using the RF variable importance approach. Decision trees (DTs) of the RF model and the nodes that link them provide information about the contributions of the input variables used in prediction, and it is this information that enables their importance to be assessed. When an upper node (parent node) is divided into two lower nodes (child nodes) by the variable i, the importance of the variable is defined as the difference between the root-mean-square error (RMSE) of the parent node and the sum of the RMSEs of the child nodes [68]. According to Hwang et al. [68], the importance of variable i for a given decision tree ${I\left(f_i\right)}_{DT}$ is defined by equation 1 below.

$${I\left(f_i\right)}_{DT}=\frac{{\sum }_{j:nodejsplitson{f}_i}{G}_j}{{\sum }_{a\in allnodes}{G}_a}$$

(1)

where $G_j$ is the information gain of node j branched by the variable i and $G_a$ of all nodes. The importance of each variable can then be normalized according to the following equation (2):

$${I\left(f_i\right)}_{norm}=\frac{I\left(f_i\right)}{{\sum }_{j\in allfeatures}I\left(f_j\right)}$$

(2)

Finally, the importance of each variable in the RF model is averaged over all decision trees in RF as in equation 3:

$${I\left(f_i\right)}_{RF}=\frac{{\sum }_{j\in alltrees}{I\left(f_j\right)}_{norm}}{N_T}$$

(3)

where ${I\left(f_i\right)}_{RF}$ is the importance of variable i in the RF model and $N_T$ is the total number of decision trees in RF. Ranking all the variables, therefore, allows us to retain only those with the highest values of ${I\left(f_i\right)}_{RF}$.

When testing the four modelling algorithms, we also used the SHAP (SHapley Additive exPlanations) method [69,70], which allowed us to select certain predictors with much greater effects on the performance of the models. This method is a game-theoretic approach that explains the predictions of machine-learning models and, thus, facilitates their interpretation [71,72].

2.6.2. Development of Prediction Models

We evaluated four machine learning algorithms, namely Random Forest (RF), Support Vector Machine (SVM), Extreme Gradient Boosting (XGBoost), and Deep Neural Network (DNN). The goal was to assess their relevance in order to select the best model, which was optimized using an automated learning tool. The objective is to construct a model $M\left(\theta \right)$ according to the equation:

$$M\left(\theta \right):\widehat{h}=f_{\theta }\left(x_N\right)$$

(4)

with $\theta$ representing parameters and hyperparameters of the model to be learned and $\widehat{h}$ is predicted height. This model should be capable of predicting the height that minimizes the cost function:

$$argmin{\sum }_{x\in D_T}l\left(\widehat{h},h\right)$$

(5)

where $l$ represents the loss function. Equation (5) is a general form of the loss function, but here we have used root-mean-square error (RMSE) as the measure of loss. The four tested algorithms are briefly presented below.

• Random Forest

RF is a learning algorithm that constructs a series of decision trees generated by training samples taken at random with or without replacement. It uses the decision trees $t_j\left(P,{\Theta }_j\right)$ as the training base, where j is the number of base trees. Given the training database $D$ that is defined above for a particular realization ${\theta }_k$of ${\Theta }_j$ (with $k=\left\{1,...,j\right\}$), the trained decision tree is defined by $\widehat{t}\left(p,{\theta }_k,D\right)$. Although this formulation follows Breiman [73], the random realization ${\theta }_k$ is implicitly used to introduce randomness in two ways. First, "bagging" fits each tree to a randomly drawn sample of the original database. Second, when splitting a node, the best split is found on a randomly selected subset of p predictors rather than all predictors, independently at each node.

Decision trees are then constructed without being pruned; the resulting trees are combined as a weighted average, which is presented in equation 6.

$$f\left(x\right)=\frac{1}{J}{\sum }_{j=1}^{J}t\left(p\right)$$

(6)

Although the RF algorithm has demonstrated its prediction ability, three parameters can be adjusted to improve its performance, depending on the situations and applications. These are the number of predictive variables that are randomly selected at each node (p), the number of trees in the random forest (J) and the size of the tree [74].

• Support Vector Machine

The SVM algorithm is a supervised, non-parametric statistical learning technique, the initial purpose of which was to solve binary classification problems, which was later extended to regression problems [75]. The objective here is to model a function that would allow us to predict canopy height $\widehat{h}$ while maximizing the hyperplane, i.e., the margin between the predicted value ($\widehat{h}$) and the actual value (h), for all the training data. In the case of a non-linear support vector regression approach, this involves applying a transformation ($\varphi :{\mathbb{R}}^n\to {\mathbb{R}}^{\alpha }\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}x\to \varphi \left(x\right)$) from the space of the input data to a higher dimensional space in order to better predict the input data. Note that $n$ and $\alpha$ are respectively dimensions of the space of the input data and transformed with $\alpha >n$. A kernel function such as a radial basis function is generally used to solve the non-linearity problem [76,77]. From then on, linear regression in $\alpha$ dimensional space will be written as follows:

$$f\left(x\right)=<w,\varphi \left(x\right)>+b$$

(7)

where $w=\left(w_1,......,w_N\right)\in {\mathbb{R}}^n$ is the hyperplane coefficient vector; $b\in \mathbb{R}$ is a scalar denoting bias.

• Extreme Gradient Boosting

The XGBoost algorithm is a machine-learning ensemble model, which is an efficient and scalable implementation of the gradient boosting machine algorithm [78]. As a learning set approach XGBoost uses multiple decision trees to achieve optimal prediction performance. The output of the model predicted by this approach will have the same decision rules as a classical decision tree model. Let K be the number of trees used, and the output prediction result is the sum of all the scores predicted by K trees, as shown in equation (8):

$$\widehat{h}={\sum }_{k=1}^K{f}_k\left(x_m\right),f_kϵF$$

(8)

XGBoost adopts the same gradient boosting as the Gradient Boosting Machine (GBM) algorithm [79], but provides a small improvement to the objective function by regularizing it, as presented in equation (9):

$$\mathcal{L}\left(\theta \right)={\sum }_i\mathcal{l}\left({\widehat{h}}_i,h_i\right)+{\sum }_k\Omega \left(f_k\right)$$

(9)

Here, $ℒ$ is the total objective function, $\theta$ represents the hyperparameters of the model, ℓ is a differentiable convex loss function which measures the distance between the prediction ${\widehat{h}}_i$ and the true value $h_i$, the second term representing the regularization which reduces the variation in output of the new tree. Detailed information regarding the XGBoost model can be found in the literature [80,81].

In order for the model to work efficiently, XGBoost also stores data in in-memory units for parallel learning, thereby allowing it to handle larger datasets and to run much more rapidly [80].

• Deep Neural Network

DNN is a neural network that is organized into several hidden and densely connected layers, which characterize the input-output parameters of the network [82]. DNNs were selected because this architecture guarantees a high capacity for finding relationships between variables and for generating machine learning based on data representations [83]. The capacity that characterizes DNN is based upon the fact that each layer is continuously updated by repetitive learning, which is referred to as "backpropagation", to find the appropriate weights and biases [83]. Backpropagation is carried out until the difference between the predicted value and real value (the error) is optimal. In this perspective, the output (Oj) of the DNN layer j is defined according to equation 10, considering the input X (X = {x1, . . ., xN}) the activation function (σ), the weight matrix (w) and the bias vector (b)

$$O_j=\sigma \left({\sum }_{N}X.w_N+b_N\right)$$

(10)

This output is accurately predicted through careful adjustment of parameters, such as the activation function, the learning rate, number of neurons in each hidden layer, the number of hidden layers, batch size, and number of epochs, i.e., the number of forward-backward passes through the dataset [47], [84,85,86].

When using Machine Learning (ML) algorithms, it is sometimes difficult to manually identify the best hyperparameters for good model training or generalization [87,88]. To achieve this goal, the ML community has used approaches such as grid search and random search [89] or Bayesian optimization [90]. Yet, these ML models are configured by a set of hyperparameters with values that can substantially affect their performance, which means that we cannot know whether a given technique is truly better or simply better tuned [91,92]. Automated machine learning (AutoML) approaches automate the selection of algorithms and the implementation of other parallel operations to efficiently optimize hyperparameters, while taking into account the particularities of the input data. For this purpose, we employed the AutoML TPOT (Tree-based Pipeline Optimization Tool) and AutoGluon to compare their performance with that of the algorithm, which offered the most accurate results among the four being tested [93,94].

2.6.3. Performance Evaluation of the Developed Models

ICESat-2 canopy height was validated using Pearson's Correlation Coefficient (r) with reference to the measurements, before being used as reference data for canopy height prediction modelling. The prediction models and canopy height maps that were subsequently produced were evaluated, comparing one to another, and with existing models. The comparisons were made using traditional performance indicators, i.e., correlation r, RMSE (Root-Mean-Square Error) and MAE (Mean Absolute Error) to retain the models that stood out in terms of their performance. Equations 11 to 13 show the mathematical expressions of the three indicators:

$$r=\frac{{\sum }_{i=1}^n\left(x_i-\widehat{x}\right)\left(y_i-\widehat{y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(x_i-\widehat{x}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(y_i-\widehat{y}\right)}^2}}$$

(11)

$$RMSE=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(x_i-y_i\right)}^2}$$

(12)

$$MAE=\frac{1}{n}{\sum }_{i=1}^n\left|x_i-y_i\right|$$

(13)

where $x_i$ represents the ith observed value of the canopy height data extracted from ICESat-2 or GEDI, and $y_i$ is the ith predicted value; $\widehat{x}$ and $\widehat{y}$ are respectively the means of all $x_i$ and $y_i$, while $n$ represents the total number of canopy height samples from ICESat-2 or GEDI.

These parameters were used to validate models developed from ICESat-2 data on the basis of data collected in the field from their footprints. The same parameters were calculated to evaluate model performance by comparing predicted canopy heights on the plots of the second National Forest Inventory (NFI2) with those measured in the field during the inventory. They also enabled us to compare the cartographic products of this study with those of other authors.

2.7. Forest Height Mapping

Following evaluation of the models developed, those obtaining the best performance in predicting canopy height on the satellite LiDAR footprints were retained for cartographic inference. Indeed, the execution of these models yielded predicted canopy heights in only LiDAR footprints, leaving blank areas between them. To move from spatially discontinuous to continuous data, we first created a stacked multi-band image, using the variables that had contributed the most to prediction, according to variable importance defined in Section 2.6.1. Given that each band in this stacked image was considered a predictor, execution of the trained model uses the content of each pixel and their corresponding values in the underlying bands, to produce a new height pixel. Several iterations of the process on all pixels of the stacked image complete the image of the heights that is produced. Figure 5 illustrates the cartographic inference that is performed with the models, which were developed from GEDI or ICESat-2 data.

The resulting canopy height maps were formatted in ArcGIS. To analyze differences between types of LiDAR data that were used, some statistical metrics of the map resulting from cartographic inferences with the GEDI-based model were compared with those obtained with the ICESat-2-based model. These two maps were compared to similar data existing locally or globally to obtain an understanding of the particularities that are related to our study area. The performance evaluation parameters that are presented in Section 2.6.3 were also used during these comparisons.

3. Results

3.1. Validation of the Reference Data

The correlation matrix presented in Table 5 compared 98% relative heights extracted from ATL08 data of ICESat-2, with those calculated from field data. In Table 5, the values that were derived from the field measurements are shown, ranging from the minimum value to the 98% relative heights. They are compared with one another, and were compared respectively with h_canopy, which is RH98 height that was extracted from the ICESat-2 data.

The last row of the correlation matrix summarizes the correlations (r) between h_canopy and the field data (Table 5). These positive correlations range from weak to moderate associations, i.e., the last being 98% relative height (r = 0.53; RMSE = 4.85; MAE = 3.84). We used these ICESat-2 data as reference data and this relative height as a predictor of canopy height in the different modelling scenarios examined in this study.

3.2. Selection and Combination of Multisource Variables

The scenarios where different combinations of variables were applied in preliminary modelling allowed not only the importance of variables to be evaluated, but also allowed four different algorithms to be tested: RF, SVM, XGBoost, and DNN. The methodology and the different tested scenarios S1 to S7 are explained in Section 2.6 (see Table 4). Table 6 indicates, in the form of a heatmap, the performance metrics of the canopy height prediction models for the seven scenarios applied to the four algorithms. At this stage, only ICESat-2 data were used to detect the scenario and to assess the algorithm being used at the actual modelling stage.

The RF approach obtained the best performance metrics among all tested combinations with S7 scenario. If we rank each of the performance metrics for S7 across the four models, their agreement is weakly consistent (Kendall’s W = 0.48, p = 0.23; W = 1, where there is complete agreement among rankings). As the best model, RF ranked highest, while the poorest-performing model was DNN, with the lowest ranks. In comparing performance among the seven RF scenarios, the three metrics (R, RMSE, MAE) were strongly consistent (W = 0.865, c²_r = 15.57, df = 3, p = 0.0014), which is reflected in the heatmap. We ordered RF scenarios from worst to best, as: S3 (Topographic) ≤ S2 (Radar) ≤ S6 (Radar-Topographic) < S1 (Optical) = S4 (Optical-Radar) < S5 (Optical-Topographic) ≤ S7 (Optical-Radar-Topographic). We further analyzed the importance of variable contributions to the RF/S7 combination. These results are depicted in Figure 6. Optical variables apparently contribute much more to modelling than do those derived from other data sources.

Over the course of modelling, the SHAP method allowed us to select about 20 predictors that exerted the strongest effects on model performance. Indeed, the importance of these features in estimating canopy height was further evaluated. For example, Figure 7 depicts the results for RF and XGBoost algorithms.

Fourteen importance features were common to the two algorithms. Yet, in the absence of unique features, those features that were shared by RF and XGBoost were very consistent in their respective rankings (W = 0.953, c²_r = 24.77, df = 13, p = 0.025). Similar depictions of feature importance in the development of the remaining models (SVM, DNN) using SHAP are presented in Figure A2 (Appendix C).

3.3. Modelling Canopy Height Using ICESat-2 Data

Variables selected during the preliminary modelling stage and RF algorithm that resulted in remarkable performances allowed us to obtain the following results. It must be emphasized that only footprints containing more than 50 photons were considered to ensure that representative data were analyzed. Ultimately, 9781 ICESat-2 footprints were distributed in two groups of height classes that were defined in Section 2.6. Slightly higher metrics were obtained with Group 2 (r = 0.58; RMSE = 5.33; MAE = 3.92), compared to Group 1 (r = 0.51; RMSE = 5.46; MAE = 4.05).

Subsequent use of AutoML TPOT (see Section 2.6.2) made it possible to optimize the learning architecture by intelligently exploring thousands of pipelines (i.e., processing chains from preprocessing to modelling). This exploration made it possible to find the appropriate pipeline and automatically choose the appropriate hyperparameters for the model that best suited our data. The processing chain that led to these results is represented in Figure 8.

TPOT contributed to an improved performance of the canopy height prediction model that was developed. We refer to this model as rf_icesat-2_rh98, a designation that refers to the RF algorithm, ICESat-2 data and RH98 from which it was developed. The metrics obtained with this model are: r = 0.65; RMSE = 5.10; and MAE = 3.80. The regression curve, indicating the dispersion of heights extracted from ICESat-2 compared to predicted values, is shown in Figure 9. An initial observation of the distribution of predictions on this graph compared with the 1:1 line shows that the model tends to overestimate small canopies (< 10 m), to correctly estimate medium canopy heights ([10 m -20 m]) and to underestimate large canopies (> 20 m).

AutoGluon was applied to improve modelling performance. Given that it is a very complex ML system, AutoGluon is computationally very intensive, resource-intensive, difficult to debug and may make inappropriate assumptions regarding both parameters and data types [95,96,58]. The results that we subsequently obtained (r = 0.64; RMSE = 5.12; MAE = 3.83) were slightly higher than those of a simple RF (see Table 5), yet they remain very close to those obtained with TPOT.

3.4. Modelling Canopy Height from GEDI Data

The same tools and methods that were used with the ICESat-2 data to inform the choice of hyper-parameters and selection of variables with strong contributions to modelling were applied to predictions from the GEDI data. Table 7 summarizes performance metrics for nine configurations of canopy prediction models defined during data filtering (Section 2.5.2) and developed from GEDI.

Table 7 presents a heatmap of accuracy metrics for 63 prediction models that had been established, with seven relative heights under the nine configurations (see Section 2.5.2). For reference purposes, these models are designated as rf_gedi_configx_rhy, i.e., the model that was established with the Random Forest algorithm, based on GEDI data in configuration x with relative height y, (where x is the configuration number ranging from 1 to 9, and y is the percentage of relative height increasing from 75 to 100% in 5% increments. Through use of the simple RF algorithm and considering only Pearson coefficients, the rf_gedi_config9_rh98 model attained relatively high performance (r = 0.80), compared to the other models.

The heatmaps that were created for Pearson’s r, RMSE and MAE metrics exhibit similar overall trends, i.e., Configuration 9 emerges as the best performer, while Configurations 6, 2 and 1 are the worst. The remaining configurations are intermediate between the two extremes. In ranking correlations across configurations (columns) for each relative height category (rows) and subsequently determining their consistency, we found strong concordance among relative heights for Pearson’s correlations (Table 7). Visual assessments of r colour-ratings were consistent with rank ordering of the RF configurations. Strong ordering from worst to best performance was likewise noted for RMSE (ranked from highest to lowest error). MAE ranks (highest to lowest error) were also very consistent. Despite high levels of concordance, rank ordering was not identical among the performance metrics, given that error estimates also appeared to increase with increasing relative height percentiles, particularly for the worst performance scenarios (Table 7).

It should be noted that results of modelling with the RF algorithm from RH98 in group 1 (r = 0.73; RMSE = 4.95; MAE = 3.65) and in group 2 (r = 0.74; RMSE = 4.93; MAE = 3.66) also remain lower than those obtained for the GEDI data without the height classes. With the use of AutoML TPOT, the optimization of hyperparameters allowed us to produce a model with much higher performance (r = 0.84; RMSE = 4.15; MAE = 2.36) that was based on RF. The processing chain that yielded better results is represented in Figure 10.

In total, 28478 GEDI footprints that met the filtering criteria of Configuration 9 produced the more efficient model rf_gedi_config9_rh98. The regression curve in Figure 11 indicates strong dispersion of GEDI-based relative heights compared to the predicted values. Here again, according to the distribution of predictions in relation to the 1:1 line, it can be seen that the model slightly overestimates small canopies (< 10 m), makes a relatively better estimate of medium canopy heights ([10m -25 m]) and underestimates larger canopies (> 25 m). But this distribution is much more aligned or oriented more closely with the 1:1 line than does the scatterplot depicted in Figure 9.

Application of AutoML AutoGluon to the GEDI data resulted in a model with good performance (r = 0.83; RMSE = 4.16; MAE = 2.65) compared to using the RF algorithm alone. The results are comparable to those obtained with Auto-ML TPOT. Table 8 shows the effects of using Auto-ML TPOT and AutoGluon in improving the performance of models developed with both types of satellite LiDAR data compared to models developed simply with the RF algorithm.

Both sources of data showed consistent improvement in performance metrics with the application of AutoGluon to the RF model, which was then followed by an improvement with TPOT, with progressively increasing r and progressively decreasing error values with each improvement (W = 0.975, c²_r = 14.62, df = 5, p = 0.012). Mean (± SD) Pearson coefficient for GEDI was 0.823 (± 0.021), while that of ICESat-2 was 0.633 (± 0.021). RMSE and MAE estimates were consistently lower for GEDI compared to ICESat-2. The three performance metrics obviously differed between the two datasets, significantly so (1-df directed contrast: Z = 11.57). The application of AutoGluon produced metrics that were slightly lower than those obtained with TPOT, which exhibited the best performance. We then continued analyses of the results of this study with the models developed using the TPOT method.

3.5. Forest Canopy Height Mapping from Developed Models

3.5.1. Forest Canopy Height Map Created from the ICESat-2 Based Model

The rf_icesat-2_rh98 is the model that was selected from the analysis of the ICESat-2 data. It was used to produce the continuous canopy height map of the study area at a spatial resolution of 30 m. The map encompassing Ecological zone 4 is presented in Figure 12.

Regarding the performance of the model used to produce this map, the predicted minimum and maximum canopy heights are respectively 4.20 m and 38.75 m, while predicted mean height (± SD) is 14.26 m (± 4.24 m). Figure 12 contains two zoom rectangles at different scales (Fig. 12b and Fig 12c), which present patterns that were very similar to those that are observed in high-resolution Google Earth images corresponding to the extents of these zooms (Figs. 12d and 12e, respectively).

3.5.2. Forest Canopy Height Map from GEDI-Based Model

After all analyses were completed and improvements were made to the developed models (see Section 3.4), the rf_gedi_config9_rh98 model was retained for the GEDI data. It was used to produce a continuous canopy height map of the study area at a spatial resolution of 30 m. This map is presented in Figure 13 and follows the same format as Figure 12.

This continuous map is likewise the result of predicting canopy height using the model that learned best from GEDI (spatially discontinuous) data, relying upon a combination of continuous multisource variables, which provided the most accurate results when developing the model. Given the performance of the model used to produce this map, predicted minimum and maximum canopy heights are respectively 2.56 m and 44.20 m, while mean (± SD) height is 11.23 m (± 5.17 m). Figure 13 contains two zooms at different scales (Figs. 13b and 13c), which also show patterns almost similar to their corresponding high-resolution Google Earth images for the same footprints (Figs. 13d and 13e).

3.6. Comparative Analysis of Developed Models with Existing Products

No similar products exist for comparable studies that have been conducted locally in our study area to which our results could be compared. Regression was conducted on the field data collected from the ICESat-2 footprints, and those values predicted from these same plots made it possible to validate the model that was developed with these data (r = 0.54; RMSE = 3.11; MAE = 2.54). Given that we did not have field data for footprints of the GEDI products to which we could compare predictions of the optimal model established from these data, we used relative heights derived from data collected in the second National Forest Inventory (IFN2) to perform a similar regression. The only maps of canopy height available for the area are global maps that had been created by Lang et al. [97] and Potapov et al. [10]. From these maps, canopy heights were extracted from footprints of the ICESat-2, GEDI and IFN2 plots to compare them with their predicted values. Various linear regressions conducted with these data made it possible to estimate correlations between the extracted or existing data and their predicted values (Table 9). For the sake of simplicity, we refer to Lang et al. [97] and Potapov et al. [10] as “Lang" and "Potapov," respectively, in the following tables and figures, and in their canopy height mapping products.

This table shows that the predictions of canopy height made with the model based on data extracted from ICESat-2 are more closely correlated with those of the Lang map (r = 0.71, RMSE = 3.38, MAE = 2.55) than they are with those of Potapov (r = 0.62, RMSE= 3.80, MAE = 2.93). Similarly, predictions of canopy height using the model based on data extracted from GEDI are more consistent with Lang's map (r = 0.65, RMSE = 5.50, MAE = 4.17) than with Potapov's (r = 0.55, RMSE= 6.04, MAE = 4.64). On the other hand, the heights measured during NFI2 are closer to the canopy height estimates of the model based on data extracted from GEDI (r = 0.63, RMSE= 3.40, MAE = 2.65), than they are to those predicted with the model based on ICESat-2 (r = 0.55, RMSE = 3.65, MAE = 2.98).

In order to better compare the canopy height maps produced during this study, both with each other and with those of other authors, most notably Lang and Potapov, we performed image subtractions and then analyzed the results. Histograms of the maps resulting from the models developed by this study and those of their differences are presented in this section. Those maps resulting from the subtraction between the GEDI-based map and those of Lang and Potapov are presented and analyzed in the Discussion (see Section 4.3). Figure 14 presents the histograms of the GEDI-based map, the ICESat-2-based map and the map resulting from differences between these two maps.

The analysis of the histogram of the map that is based on GEDI (Figure 14a) reveals that more than 300 000 pixels have heights less than 5 m, while ICESat-2 (Figure 14b) shows less than 50 000 pixels in that range. This implies that the prediction model based on GEDI is much more sensitive to shorter canopies compared to that of ICESat-2. This response is also observed in the histogram of the resulting difference map (Figure 14c), which indicates most pixels are negative. The analysis of this last histogram further reveals that the model developed from ICESat-2 data overestimates heights most of the time, compared to that developed from GEDI data. Nevertheless, the average of the deviations is relatively small, and it should be noted in this histogram (Figure 14c) that the curve of the normal distribution is completely flattened in the tails, tending towards 0. This means that pixels with large positive or negative deviations are very few and that the source maps from which these difference maps are derived are relatively close to one another. The maps relating to these different image subtractions can be consulted in Figure A3, Figure A4, and Figure A5 in Appendix D, while the related discussions are presented in session 4.3.

4. Discussions

4.1. Performance of Multisource Satellite Variables in Estimating Forest Height

In this study, variables from optical, radar and topographic data were integrated with ICESat-2 and GEDI satellite LiDAR data. These combinations were used to develop models for estimating canopy height of forest-savanna mosaics in the Sudano-Guinean zone of West Africa. The importance of these variables in contributing to model development depends upon the method that is used for the assessment, and the algorithm that is used for this purpose. For example, Figure 7 reveals that in predicting responses with the RF model (Figure 7a), the ten most important covariates (in order) are swir2, swir1, slope, elevation, ndbi, ndii, vari, s1vv, mndwi and blue bands. In contrast, predictions made with the XGBoost model (Figure 7b) consisted (in order) of 10 covariates: swir1, slope, ndbi, swir2, elevation, vari, nirnarrow, rededge2, blue and arvi.

From literature reports, it is generally accepted that radar data would exhibit increased sensitivity to the vertical structure of the forest and the variables derived from them would allow for perfect estimation of the parameters of forest vertical structure, including canopy height [97,98,99] . Contrary to this assertion, which is widely accepted by several other studies, our results illustrated in Figure 6 indicate a strong contribution of variables from optical and topographic data sources to canopy height estimation, while those derived from radar data contribute very little. Our finding is consistent with Xi et al. [17], who also noted that vegetation indices and topographic information from Sentinel-2 and SRTM optical data respectively contributed much more effectively to the establishment of canopy height prediction models compared to texture measurements and backscatter variables from Sentinel-1 radar data. This result of our research is supported by the study by Luo et al. [100], who concluded that Sentinel-2-derived variables significantly contributed to the canopy height estimation model, unlike backscatter coefficients and textural parameters derived from Sentinel-1. In future studies, further investigations of contributions of variables by source would be relevant, together with the effects of their combinations in estimating forest structural parameters, particularly in complex tropical ecosystems of forest-savannah mosaics.

4.2. Comparative Analysis of ICESat-2 and GEDI Data Performance

In combination with multi-spectral, radar and topographic data, ICESat-2 and GEDI data have made it possible to develop models for predicting canopy height. The metrics of the best model that was established from ICESat-2 data are: r = 0.65; RMSE = 5.10 m; and MAE = 3.80 m. In contrast, metrics of the model established with GEDI data appear to be much improved (r = 0.84; RMSE = 4.15 m; MAE = 2.36 m). Therefore, GEDI emerges as providing better estimates of heights compared to ICESat-2. Moreover, these different results had been obtained by considering relative height RH98. Furthermore, the in situ heights measured during NFI2 national inventory are more consistent with the canopy heights estimated by the GEDI model (r = 0.63; RMSE = 3.40; MAE = 2.65) than with those estimated by the ICESat-2 model (r = 0.55; RMSE = 3.65; MAE = 2.98). The superiority of GEDI over ICESat-2 has been reported in previous work in other regions. This is the case, for example, of Zhu et al. [9], whose models with relative height RH98 obtained RMSE error parameters ranging from 3.61 m to 4.23 m with GEDI data, and from 4.76 m to 10.23 m with ICESat-2. The performance superiority of GEDI versus ICESat-2 is consistent with several other literature reports, including those of Liu et al. [13], Liu et al. [58] and Zhu et al. [101], who likewise had reported better estimates with GEDI.

One possible reason for this response is that relative heights contained in the ICESat-2 data, each representing an average value over 100 m, consist of several laser pulses that are spaced 0.7 m apart (Figure 3a) and, thus, contain a certain level of imprecision. What results in a low level of correlation of tree heights with ground surface data collected in their footprints. The other reason is that the GEDI data, when acquired, are more densely sampled than the ICESat-2 data and have individual and direct ground footprints of their emitted laser beams (Figure 3b). In this study, the ground data were collected in the footprints of the ICESat-2 data only. This does not permit direct validation of the estimates that are made using GEDI data. Collection of ground data in GEDI footprints, therefore, would be of great importance for future work in our study area, or more generally, in forest-savannah mosaics.

4.3. Comparison of Map Products with Similar Recent Work

To properly compare the cartographic products of this study with existing ones, we established histograms from the maps resulting from the differences between the GEDI-based map and that of Lang on the one hand, and the GEDI-based map and that of Potapov on the other hand. Given that these two researchers did not use the ICESat-2 data, this comparison could only be done with our GEDI-based model. Figure 15 below presents the distribution of these differences in canopy heights between these products.

It is evident from Figure 15a that the mean of the deviations is negative, which could mean that the Lang model overestimates some canopy heights compared to our GEDI-based model. This is opposite to the phenomenon in Figure 15b, in which the histogram exhibits a positive mean of the deviations. The Potapov model, therefore, underestimates some canopy heights compared to our model. In this situation, where the models of the two researchers appear to overestimate or underestimate canopy height compared to our model, could be related to differences in algorithms that are used by each of these studies regarding the development of these models. Indeed, while we used the RF algorithm, Lang et al. [97] considered a probabilistic deep-learning approach based upon an ensemble of deep convolutional neural networks (CNNs), while Potapov et al. [10] employed a machine-learning algorithm based on an ensemble of bagging regression trees. In this case, it would be interesting for future studies to further investigate effects of the choice of algorithms on the quality of the developed models. .

The standard deviation of the difference map between the GEDI model versus Lang appears to be lower than that depicted in the difference map between the GEDI model versus Potapov. This implies that the distribution of these deviations is much more tightly grouped about their mean in the first case, while in the second case, these deviations are more dispersed and relatively far from their mean. This observation indicates that the GEDI-based model of this study would be closer to that of Lang than to Potapov. Nevertheless, the mean of the deviations is relatively low, and we also note that in these two histograms, the curve of the normal distribution is completely flattened in the tails, tending towards 0. Furthermore, the number of pixels with large positive or negative deviations is quite limited and would be more likely linked to an edge-effect of the images. As in the case of Figure 14c, the source maps from which each of these difference maps were derived are therefore relatively close to each other.

4.4. Important Factors and Limitations in Estimating Canopy Height

Histogram analysis had demonstrated that the more sensitive the model is to shorter vegetation, the better it would be able to correctly estimate the canopy heights of our study area. Our study area consists of forest-savannah mosaics, which are ecosystems with mostly small and often scattered trees. In these ecosystems, the GEDI-based model therefore performed better than the ICESat-2-based model, given its ability to better estimate the heights of both tall and short canopies. This performance would also be linked to the high density of GEDI footprints compared to ICESat-2. Similar canopy height estimation studies by Liu et al. [58] in China, Zhu et al. [9] in the United States, and Sothe et al. [57] in Canada have shown remarkable performance, partly due to the continuity and homogeneity of the forest landscape in which their models were applied. Therefore, further investigation is warranted regarding the effect of different ecosystems and forest types on performance levels of canopy height estimation models developed from GEDI and ICESat-2 data, especially in forest-savannah mosaics of the Sudano-Guinean zone.

Furthermore, the choice of the dependent variable in such predictions greatly influences the accuracies of the models that are being developed. In this study, we evaluated the relative heights that were extracted from satellite LiDAR data before selecting RH98 as the optimal dependent variable in model development. Our results from regressions used to assess correlations between extracted or existing data and predicted values (Table 9) also indicate that the choice of a prediction parameter is very critical to the modelling process. Indeed, Lang et al. [97] likewise used RH98 as a dependent variable, while Potapov et al. [10] used RH95. The NFI2 data were slightly further away from the canopy heights extracted from the Potopov map (r = 0.46; RMSE = 4.21; MAE = 3.28), but closer to the canopy heights extracted from the Lang map (r = 0.64; RMSE = 3.96; MAE = 3.09), which are most similar to our own (r = 0.63; RMSE = 3.40; MAE = 2.65). Therefore, our study joins those of Lang et al. ([97]) and Ngo et al. [102], as well as many others, in confirming that RH98 would be the most appropriate RH to be used in best estimates of canopy height based upon GEDI data.

Ultimately, the development of canopy height prediction models from satellite LiDAR data is strongly dependent upon several parameters. In order to optimize the parameters that would make these models more accurate and geographically adapted, the spatio-temporal coherence of the independent covariates that are used must be taken into account, together with the characteristics of the ecosystems to which these models are being applied, and the algorithms or learning models that are used. To the best of our knowledge, such studies involving the combination of optical, radar, topographic and satellite LiDAR data to develop models for estimating canopy heights in the forest-savanna mosaics of the Sudano-Guinean zone did not exist prior to ours. Taking into account the elements of discussions and perspectives of this study would help improve the quality of canopy height or forest biomass models in future studies.

5. Conclusions

In our search for greater accuracy in estimating the canopy height of forest-savannah mosaics in the Sudano-Guinean zone, the research relied upon a combination of multisource and multisensor data. Four machine-learning algorithms (RF, SVM, XGBoost, DNN) were evaluated before selecting Random Forest, which was demonstrably more efficient in predicting canopy height across the study area than the remaining models. Our results indicate that the final model developed from GEDI data is more efficient than that derived from ICESat-2. Investigations carried out during this research also reveal that the prediction models based on grouping data by height classes did not provide any improvement compared to those where the height classes were not defined, using both GEDI and ICESat-2 data.

Estimation of canopy height was best achieved using combinations of several types of remote sensing data rather than using each one in isolation. Yet, we found that variables derived from optical and topographic data contributed much more to the development of better performing models that did those derived from radar data, which showed very little sensitivity. Furthermore, future studies should be able to assess the effects of different variables extracted from the radar data for the estimation of forest and vegetation structural parameters. The return signal that is received by satellite LiDAR sensors depends upon the characteristics of the land cover (height, density, canopy closure). Future studies, therefore, could also analyze the effects of ecosystem characteristics on the quality of GEDI-based models, especially in the ecosystems within our study area, which are characterized by sparse and relatively small vegetation types or patches.

In the short-term, our results of canopy height modelling could be used by local decision-makers for forest management in the study area by favouring the use of GEDI data in estimating canopy height. In the long-term, the estimation of field dendrometric parameters within GEDI data footprints is more necessary than ever in order to better validate the models that have been developed. Validation of these GEDI data and models in this particular eco-climatic zone would also provide forest managers with an appropriate tool for estimating aboveground biomass to better understand forest dynamics and carbon fluxes and, thus, adapt their practices and management methods to the requirements of REDD+.