Flor Alvarez-Taboada

Miro Govedarica

David Morin
email: david.morin@univ-tlse3.fr

Milena Planells
email: milena.planells@cnes.fr

Nicolas Baghdadi
email: nicolas.baghdadi@teledetection.fr

Alexandre Bouvet
email: bouveta@cesbio.cnes.fr

Ibrahim Fayad
email: ibrahim.fayad@inrae.fr

Thuy Le Toan
email: thuy.letoan@cesbio.cnes.fr

Stéphane Mermoz
email: mermoz@globeo.net

Ludovic Villard
email: villardl@cesbio.cnes.fr

Improving Heterogeneous Forest Height Maps by Integrating GEDI-Based Forest Height Information in a Multi-Sensor Mapping Process

Keywords: forest inventory, Sentinel, ALOS PALSAR, GEDI, canopy height, support vector, temperate forest

Forests are one of the key elements in ecological transition policies in Europe. Sustainable forest management is needed in order to optimise wood harvesting, while preserving carbon storage, biodiversity and other ecological functions. Forest managers and public bodies need improved and cost-effective forest monitoring tools. Research studies have been carried out to assess the use of optical and radar images for producing forest height or biomass maps. The main limitations are the quantity, quality and representativeness of the reference data for model training. The Global Ecosystem Dynamics Investigation (GEDI) mission (full waveform LiDAR on board the International Space Station) has provided an unprecedented number of forest canopy height samples from 2019. These samples could be used to improve reference datasets. This paper aims to present and validate a method for estimating forest dominant height from open access optical and radar satellite images (Sentinel-1, Sentinel-2 and ALOS-2 PALSAR-2), and then to assess the use of GEDI samples to replace field height measurements in model calibration. Our approach combines satellite image features and dominant height measurements, or GEDI metrics, in a Support Vector Machine regression algorithm, with a feature selection process. The method is tested on mixed uneven-aged broadleaved and coniferous forests in France. Using dominant height measurements for model training, the crossvalidation shows 7.3 to 11.6% relative Root Mean Square Error (RMSE) depending on the forest class. When using GEDI height metrics instead of field measurements for model training, errors increase to 12.8-16.7% relative RMSE. This level of error remains satisfactory; the use of GEDI could allow the production of dominant height maps on large areas with better sample representativeness. Future work will focus on confirming these results on new study sites, improving the filtering and processing of GEDI data, and producing height maps at regional or national scale. The resulting maps will help forest managers and public bodies to optimise forest resource inventories, as well as allow scientists to integrate these cartographic data into climate models.

Introduction

Forests cover one third of Earth's land surface and play an essential role in stabilising climate, regulating the water cycle, and providing habitat to a multitude of life forms [START_REF] Fao; Unep | The State of the World's Forests 2020[END_REF][START_REF]Global Forest Resources Assessment 2020: Main Report[END_REF]. Besides that, forests are of great economic benefit by providing wood and timber, used as fuel or as raw material for numerous products [START_REF] Fao; Unep | The State of the World's Forests 2020[END_REF][START_REF]Global Forest Products Facts and Figures[END_REF]. However, tropical forests are threatened by deforestation [START_REF] Fao; Unep | The State of the World's Forests 2020[END_REF][START_REF]Global Forest Resources Assessment 2020: Main Report[END_REF], and the effects of climate change are more and more important for Remote Sens. 2022, 14, 2079 2 of 25 every region (tropical, boreal and temperate) with an increase in forest disturbances (fires, diseases, storms) [START_REF] Fao; Unep | The State of the World's Forests 2020[END_REF][START_REF] Dale | Climate Change and Forest Disturbances: Climate Change Can Affect Forests by Altering the Frequency, Intensity, Duration, and Timing of Fire, Drought, Introduced Species, Insect and Pathogen Outbreaks, Hurricanes, Windstorms, Ice Storms, or Landslides[END_REF] and long-term changes in the geographical distribution of forest species [START_REF] Fao; Unep | The State of the World's Forests 2020[END_REF][START_REF] Jump | Running to Stand Still: Adaptation and the Response of Plants to Rapid Climate Change[END_REF][START_REF] Yu | A Model-Based Method to Evaluate the Ability of Nature Reserves to Protect Endangered Tree Species in the Context of Climate Change[END_REF].

In France, as in western Europe, forests are one of the key elements in ecological transition policies [START_REF]European Forest Ecosystems: State and Trends[END_REF][START_REF] Pan | A Large and Persistent Carbon Sink in the World's Forests[END_REF][START_REF]National Low Carbon Strategy; French Ministry for the Ecological and Solidary Transition[END_REF]. Sustainable forest management is needed in order to optimise wood harvesting, while preserving carbon storage, biodiversity and other ecological functions of the forests [START_REF] Fao; Unep | The State of the World's Forests 2020[END_REF][START_REF]European Forest Ecosystems: State and Trends[END_REF][START_REF]National Low Carbon Strategy; French Ministry for the Ecological and Solidary Transition[END_REF]. Management techniques include timber extraction, planting and replanting of different species, and preventing and reducing the impacts of fires, diseases and storms.

For sustainable forest management and climate change mitigation measures, quantitative and dynamic information on 3D forest structure and biomass are needed by forest managers (local scale) and public bodies (local and regional scale). Furthermore, scientists need the information to understand the impacts of climate on forest health and growth in order to guide public policies for adaptation and mitigation measures.

Most countries in temperate regions conduct national forest inventories (NFI) on a five-year time scale, based on forest measurement plots sampled throughout the territory. The inventory database includes the following measurements (which are referred to in this paper as forest parameters): diameter at breast height (DBH), number of trees per hectare (tree number density or tree density), basal area (BA) and tree height. Known species-specific allometry equations allow the estimation of timber volume (i.e., wood resources) or above-ground biomass (AGB) (i.e., carbon storage) from these measurements. The inventory data are used to derive statistics of forest resources over large administrative regions. However, the forest inventory is not designed to produce geospatial distribution or high-resolution maps that are required by forest managers to carry out precise estimates at a local scale. For that purpose, forest managers usually conduct forest stands measurements themselves, that in most cases lack standardisation (due to estimation methods, spatial sampling and surveyor bias), require labor, and are less regularly updated than the NFIs performed by national bodies.

To produce forest parameter maps and to overcome the limitations of survey approaches, research studies have been carried out to assess the use of remote sensing data. The use of satellite images has been evaluated as a way to build operational and costeffective methods for mapping forest parameters over large areas. The focus is put on the existing open access images that are more and more diversified (radar, optical, high spatial and temporal resolution) and on the use of machine learning approaches [START_REF] Chen | Optimal Combination of Predictors and Algorithms for Forest Above-Ground Biomass Mapping from Sentinel and SRTM Data[END_REF][START_REF] Haywood | The Potential of Sentinel Satellites for Large Area Aboveground Forest Biomass Mapping[END_REF][START_REF] Laurin | Above-Ground Biomass Prediction by Sentinel-1 Multitemporal Data in Central Italy with Integration of ALOS2 and Sentinel-2 Data[END_REF][START_REF] Liu | Estimation of the Forest Stand Mean Height and Aboveground Biomass in Northeast China Using SAR Sentinel-1B, Multispectral Sentinel-2A, and DEM Imagery[END_REF][START_REF] Vafaei | Improving Accuracy Estimation of Forest Aboveground Biomass Based on Incorporation of ALOS-2 PALSAR-2 and Sentinel-2A Imagery and Machine Learning: A Case Study of the Hyrcanian Forest Area (Iran)[END_REF]. Several satellite programs currently provide global and open access high resolution images. In the optical domain, Sentinel-2 and Landsat-8 images are acquired, respectively, every 5 and 16 days, with a spatial resolution of 10-20 m and 30 m for visible, near infrared and short-wave infrared bands (vis, nir, swir). Optical reflectance can be related to the green leaf area of the canopy and its horizontal distribution [START_REF] Fassnacht | Estimating the Leaf Area Index of North Central Wisconsin Forests Using the Landsat Thematic Mapper[END_REF][START_REF] Melnikova | Estimation of Leaf Area Index in a Mountain Forest of Central Japan with a 30-m Spatial Resolution Based on Landsat Operational Land Imager Imagery: An Application of a Simple Model for Seasonal Monitoring[END_REF][START_REF] Turner | Relationships between Leaf Area Index and Landsat TM Spectral Vegetation Indices across Three Temperate Zone Sites[END_REF]. Studies have highlighted relationships between reflectance and forest parameters [START_REF] Castillo-Santiago | Estimation of Tropical Forest Structure from SPOT5 Satellite Images[END_REF][START_REF] Freitas | Relationships between Forest Structure and Vegetation Indices in Atlantic Rainforest[END_REF][START_REF] Woodcock | Mapping Forest Vegetation Using Landsat TM Imagery and a Canopy Reflectance Model[END_REF], although the signal quickly saturates [START_REF] Dong | Remote Sensing Estimates of Boreal and Temperate Forest Woody Biomass: Carbon Pools, Sources, and Sinks[END_REF][START_REF] Dube | Investigating the Robustness of the New Landsat-8 Operational Land Imager Derived Texture Metrics in Estimating Plantation Forest Aboveground Biomass in Resource Constrained Areas[END_REF][START_REF] Mutanga | Narrow Band Vegetation Indices Overcome the Saturation Problem in Biomass Estimation[END_REF]. Moreover, spatial texture indicators from optical sensors have been used with successful results [START_REF] Castillo-Santiago | Estimation of Tropical Forest Structure from SPOT5 Satellite Images[END_REF][START_REF] Dube | Investigating the Robustness of the New Landsat-8 Operational Land Imager Derived Texture Metrics in Estimating Plantation Forest Aboveground Biomass in Resource Constrained Areas[END_REF][START_REF] Proisy | Predicting and Mapping Mangrove Biomass from Canopy Grain Analysis Using Fourier-Based Textural Ordination of IKONOS Images[END_REF][START_REF] Sarker | Improved Forest Biomass Estimates Using ALOS AVNIR-2 Texture Indices[END_REF][START_REF] Wolter | Estimation of Forest Structural Parameters Using 5 and 10 Meter SPOT-5 Satellite Data[END_REF], and a few studies have shown that the temporal dimension can also improve forest parameter estimations [START_REF] Zhu | Improving Forest Aboveground Biomass Estimation Using Seasonal Landsat NDVI Time-Series[END_REF]. In the radar domain, Sentinel-1 carries a synthetic aperture radar (SAR), that provides data in C-band (wavelength = 5.6 cm) every 6 days at ~20 m spatial resolution, and annual global mosaics of L-band ALOS PALSAR (wavelength = 22.9 cm) imagery at ~25 m spatial resolution are also available. The SAR signal penetrates into the forest canopies, with deeper penetration at longer wavelengths. The backscatter energy is a function of the geometric and dielectric properties of the volume of scatterers in a forest, which interacts with the incoming waves. Indirectly, SAR backscatter could be related to the water content of the scattering volume, which in turn is related to biomass, as highlighted in large scale studies [START_REF] Austin | Estimating Forest Biomass Using Satellite Radar: An Exploratory Study in a Temperate Australian Eucalyptus Forest[END_REF][START_REF] Baghdadi | Evaluation of ALOS/PALSAR L-Band Data for the Estimation of Eucalyptus Plantations Aboveground Biomass in Brazil[END_REF][START_REF] Dobson | Estimation of Forest Biophysical Characteristics in Northem Michigan with SIR-C/X-SAR[END_REF][START_REF] Le Toan | Relating Forest Biomass to SAR Data[END_REF][START_REF] Mermoz | Biomass Assessment in the Cameroon Savanna Using ALOS PALSAR Data[END_REF][START_REF] Santoro | Reviewing ALOS PALSAR Backscatter Observations for Stem Volume Retrieval in Swedish Forest[END_REF]. However, for application at the local scale, relationships between SAR backscatter and forest biomass could be disturbed by other types of parameters (forest structure, soil and vegetation moisture, understory rebound, etc.), which need to be accounted for. A growing number of studies have shown that the combination of optical data, SAR data and textural data improves forest parameter estimations [START_REF] Chen | Optimal Combination of Predictors and Algorithms for Forest Above-Ground Biomass Mapping from Sentinel and SRTM Data[END_REF][START_REF] Vafaei | Improving Accuracy Estimation of Forest Aboveground Biomass Based on Incorporation of ALOS-2 PALSAR-2 and Sentinel-2A Imagery and Machine Learning: A Case Study of the Hyrcanian Forest Area (Iran)[END_REF][START_REF] Attarchi | Improving the Estimation of Above Ground Biomass Using Dual Polarimetric PALSAR and ETM+ Data in the Hyrcanian Mountain Forest (Iran)[END_REF][START_REF] Chen | Estimation of Forest Above-Ground Biomass by Geographically Weighted Regression and Machine Learning with Sentinel Imagery[END_REF][START_REF] Gao | Timber Production Assessment of a Plantation Forest: An Integrated Framework with Field-Based Inventory, Multi-Source Remote Sensing Data and Forest Management History[END_REF][START_REF] Morin | Estimation and Mapping of Forest Structure Parameters from Open Access Satellite Images: Development of a Generic Method with a Study Case on Coniferous Plantation[END_REF][START_REF] Morin | Estimation of Forest Parameters Combining Multisensor High Resolution Remote Sensing Data[END_REF]. Due to the diversity of remote sensing data and forests, most approaches use machine learning methods. In order to map large areas of forests with a diversity of tree species and management practices, these methods require statistical sampling to have a training dataset representative of forests [START_REF] Chave | Ground Data Are Essential for Biomass Remote Sensing Missions[END_REF][START_REF] Smith | Forest Inventory and Analysis: A National Inventory and Monitoring Program[END_REF][START_REF] Goedickemeier | Sampling for vegetation survey: Some properties of a gis-based stratification compared to other statistical sampling methods[END_REF]. However, field campaigns are very expensive and cannot be multiplied to represent forests with diverse and heterogeneous conditions, as found in Europe. For example, in France there are 11 large ecological regions (called GRECO) that represent the different topological and climate conditions and imply different silvicultural practices and forest growth. Because of the diversity of forest management practices (monospecific, mixed species with different compositions and age classes, different planting, thinning practices, etc.), studies often lack a sufficient volume and representativeness of data. Results rely on the forest region and forest type where training samples are measured. The main limitations of these methods are therefore the quality, quantity and representativeness of the reference field data.

Another approach to complement or replace ground survey is the use of light detection and ranging (LiDAR) methods. Airborne LiDAR (ALS) data provide highly accurate canopy height estimates, in addition to providing information on tree shape and understory. ALS data are the most widely used data operationally to assist in forest inventories [START_REF] Jochem | Estimation of Aboveground Biomass in Alpine Forests: A Semi-Empirical Approach Considering Canopy Transparency Derived from Airborne LiDAR Data[END_REF][START_REF] Yu | Comparison of Laser and Stereo Optical, SAR and InSAR Point Clouds from Air-and Space-Borne Sources in the Retrieval of Forest Inventory Attributes[END_REF]. Several studies used LiDAR data to estimate AGB, forest height and other forest parameters [START_REF] Jochem | Estimation of Aboveground Biomass in Alpine Forests: A Semi-Empirical Approach Considering Canopy Transparency Derived from Airborne LiDAR Data[END_REF][START_REF] Gobakken | Estimating Biomass in Hedmark County, Norway Using National Forest Inventory Field Plots and Airborne Laser Scanning[END_REF][START_REF] Oliveira | Estimating Forest Biomass and Identifying Low-Intensity Logging Areas Using Airborne Scanning Lidar in Antimary State Forest, Acre State, Western Brazilian Amazon[END_REF][START_REF] Bouvier | Generalizing Predictive Models of Forest Inventory Attributes Using an Area-Based Approach with Airborne LiDAR Data[END_REF][START_REF] Holopainen | Forest inventory attribute estimation using airborne laser scanning, aerial stereo imagery, radargrammetry and interferometry-finnish experiences of the 3d techniques[END_REF][START_REF] Vastaranta | TerraSAR-X Stereo Radargrammetry and Airborne Scanning LiDAR Height Metrics in Imputation of Forest Aboveground Biomass and Stem Volume[END_REF][START_REF] Vastaranta | Airborne Laser Scanning and Digital Stereo Imagery Measures of Forest Structure: Comparative Results and Implications to Forest Mapping and Inventory Update[END_REF], and showed improved results compared to photogrammetry, radargrammetry, and satellite optical and SAR images. Although ALS data are very accurate and useful to complement or replace ground survey for forest inventories, these data remain expensive. ALS campaigns are therefore often either spatially or temporally (repetitiveness) limited. Spaceborne LiDAR has also been used to improve the estimation of forest height. The Geoscience Laser Altimeter System (GLAS, ICESat satellite mission) has been used to improve the estimation of forest height [START_REF] Baghdadi | Testing Different Methods of Forest Height and Aboveground Biomass Estimations From ICESat/GLAS Data in Eucalyptus Plantations in Brazil[END_REF][START_REF] El Hajj | Mapping of Aboveground Biomass in Gabon[END_REF][START_REF] Fayad | Regional Scale Rain-Forest Height Mapping Using Regression-Kriging of Spaceborne and Airborne LiDAR Data: Application on French Guiana[END_REF][START_REF] Fayad | Canopy Height Estimation in French Guiana with LiDAR ICESat/GLAS Data Using Principal Component Analysis and Random Forest Regressions[END_REF][START_REF] Lefsky | Estimates of Forest Canopy Height and Aboveground Biomass Using ICESat[END_REF][START_REF] Pourrahmati | Capability of GLAS/ICESat Data to Estimate Forest Canopy Height and Volume in Mountainous Forests of Iran[END_REF][START_REF] Sun | Forest Vertical Structure from GLAS: An Evaluation Using LVIS and SRTM Data[END_REF][START_REF] Yu | Estimates of Forest Structure Parameters from GLAS Data and Multi-Angle Imaging Spectrometer Data[END_REF]. GLAS has provided transects of footprints between 2003 and 2009. The footprints have a nearly circular shape of ~60 m diameter and a footprint spacing of ~170 m along the track. More recently, the Global Ecosystem Dynamics Investigation mission (GEDI, full waveform LiDAR on board the ISS) has provided an unprecedented number of forest canopy height samples from 2019 and its mission is scheduled until beginning of 2023. The size of GEDI footprints (~25 m diameter) is smaller than ICESat, providing improved measurements over forest areas with relief higher than 10%, and is closer to the spatial resolution of satellite images from the Sentinel and ALOS PALSAR programs. GEDI forest height metrics have already been studied to estimate forest parameters [START_REF] Adam | Accuracy Assessment of GEDI Terrain Elevation and Canopy Height Estimates in European Temperate Forests: Influence of Environmental and Acquisition Parameters[END_REF][START_REF] Dorado-Roda | Assessing the Accuracy of GEDI Data for Canopy Height and Aboveground Biomass Estimates in Mediterranean Forests[END_REF][START_REF] Fayad | A CNN-Based Approach for the Estimation of Canopy Heights and Wood Volume from GEDI Waveforms[END_REF][START_REF] Fayad | Assessment of GEDI's LiDAR Data for the Estimation of Canopy Heights and Wood Volume of Eucalyptus Plantations in Brazil[END_REF][START_REF] Verhelst | Improving Forest Baseline Maps in Tropical Wetlands Using GEDI-Based Forest Height Information and Sentinel-1[END_REF].

The downside of the GLAS and GEDI is that the data are footprints (15-60 m radius) spaced along transects; thus, they are not spatially continuous (i.e., not a map with data on every forest pixel). On the other hand, the large amount of LiDAR footprints could be used as forest height reference data. Studies have shown that the combination of satellite images with spaceborne LiDAR data could be used to improve the production of forest height or AGB maps at national or continental scale. GLAS data have been used in combination with optical images such as MODIS [START_REF] Chi | National Forest Aboveground Biomass Mapping from ICESat/GLAS Data and MODIS Imagery in China[END_REF][START_REF] Nelson | Estimating Siberian Timber Volume Using MODIS and ICESat/GLAS[END_REF] or Landsat [START_REF] Chi | Estimation of Forest Aboveground Biomass in Changbai Mountain Region Using ICESat/GLAS and Landsat/TM Data[END_REF][START_REF] Duncanson | Integration of GLAS and Landsat TM Data for Aboveground Biomass Estimation[END_REF][START_REF] Hansen | Mapping Tree Height Distributions in Sub-Saharan Africa Using Landsat 7 and 8 Data[END_REF][START_REF] Liu | Comparison and Evaluation of Three Methods for Estimating Forest above Ground Biomass Using TM and GLAS Data[END_REF][START_REF] Wang | Estimation of Forest Canopy Height and Aboveground Biomass from Spaceborne LiDAR and Landsat Imageries in Maryland[END_REF][START_REF] Zhang | Estimation of Forest Aboveground Biomass in California Using Canopy Height and Leaf Area Index Estimated from Satellite Data[END_REF]. The newly available GEDI data have already been used to map forest structure [START_REF] Fayad | Estimating Canopy Height and Wood Volume of Eucalyptus Plantations in Brazil Using GEDI LiDAR Data[END_REF][START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF][START_REF] Chen | Improved Estimation of Forest Stand Volume by the Integration of GEDI LiDAR Data and Multi-Sensor Imagery in the Changbai Mountains Mixed Forests Ecoregion (CMMFE), Northeast China[END_REF]. A global forest canopy height map at 30 m spatial resolution was produced [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF] and is freely available [74] (hereafter called GLAD-2019). GEDI data offer an opportunity to improve representativeness in methods using multi-sensor satellite images to estimate and map forest height.

We built a processing chain [START_REF] Morin | Estimation and Mapping of Forest Structure Parameters from Open Access Satellite Images: Development of a Generic Method with a Study Case on Coniferous Plantation[END_REF][START_REF] Morin | Estimation et Suivi de la Ressource en Bois en France Métropolitaine par Valorisation des Séries Multi-Temporelles à Haute Résolution Spatiale D'images Optiques[END_REF] for the estimation and mapping of forest structure parameters using field measurements and freely available spaceborne images (Sentinel-1, Sentinel-2 and ALOS-2 PALSAR-2 annual mosaics) in supervised learning techniques. The chain was first tested on a case study with single species and even-aged maritime pine (Pinus pinaster) plantation forests, for the estimation of AGB, BA, DBH, mean height (Hmean), dominant height (Hdom) and density. In this study, we assess the method on mixed broadleaved and coniferous forests in France, and we test the possibility of integrating GEDI height metrics as reference data for forest dominant height estimation (i.e., replacing field measurements by GEDI height metrics in model learning). In the pine forest plantations, Hdom is roughly similar to Hmean since all the trees are of the same age and species. However, in uneven-aged forests and with several species, the height of trees is often more difficult to estimate in the field (lack of a direct view of the top of the canopy due to more vertical complexity), but also the dominant or average height no longer mean the same things for forest management. The dominant height has more significance in these uneven-aged forests since it takes into account the largest trees, which are those which have the most value for logging. From the perspective of satellite data, it is also harder to build models for estimating the forest parameters of uneven-aged and mixed forests. Satellite remote sensing of forests is partly based on allometric relationships, but these relationships change from one species to another and from one forest management activity to another. In addition, in pixel-based mapping procedures, the allometric equations need to account for the mixture of species and age classes in a pixel. Working with machine learning, the most direct way to improve Hdom estimates on these highly diverse forests is to increase the number and representativeness of reference samples. Since field measurements are time-consuming and not always representative, the GEDI data could constitute a large reference dataset that can help to overcome the lack of sufficient reference data to cover the European temperate forest diversity.

In this study, we assess the use of GEDI height metrics as training data for the estimation of forest dominant height (Hdom) with multi-sensor satellite imagery. We first conduct the estimation of forest parameters (BA, DBH, Density, Hdom) using field measurements for training and validation (Section 3.2, Annex A). Secondly, we conduct the estimation of Hdom using GEDI RH95% for training and validation (Section 3.3.2). Finally, we conduct the estimation of Hdom using GEDI RH95% for training and field measurement for validation (Section 3.3.2). Results are compared to the existing canopy height map, GLAD-2019 (Section 3.1).

Materials and Methods

Study Sites and In Situ Data

The study was conducted over two test sites in northern France (Figure 1). Both sites contain public forests managed by the Office National des Forêts (ONF). The Orléans forest contains coniferous high stands and broadleaved uneven-aged high stands. Both can contain understory and coppice, but coniferous stands are generally more regular. Predominant species are oaks (Quercus robur and Quercus petraea) and Scots pines (Pinus sylvestris), with about ten other secondary species. The climate is temperate with ~640 mm annual rain and a mean annual temperature of 11.2 • C. The soil is poor and acid, and the topography is flat. The Saint-Gobain Forest (hereafter called StGobain) contains mainly broadleaved uneven-aged high stands. Predominant species are oaks (Quercus robur and Quercus petraea), beech (Fagus sylvatica) and hornbeam (Carpinus betulus), with about 15 other secondary broadleaved species. The climate is also temperate with ~700 mm annual rain and a mean annual temperature of 10.4 • C. The soil is more nutrient-rich than in Orléans, so the forests are denser and more productive. The topography is generally flat, but some slopes can reach 20 • locally.

Data collection was conducted by the ONF in early 2015 in StGobain and during the 2016-2017 winter in the Orléans forests. Tree count and diameter at breast height (DBH) measurements were conducted on circular plots from 15 to 17 m radius; only trees with a DBH larger than 17.5 cm were considered. On part of the plots, the height of the largest trees was measured to calculate the dominant height (height of the 100 largest trees per hectare). Additionally, airborne LiDAR (ALS) data were acquired on each site. LiDAR parameters for Orléans and StGobain are, respectively, 37 and 18 points/m 2 scanning density; 6 and 15 cm vertical accuracy; and 20 and 30 cm absolute positional accuracy. ALS data enabled a canopy height model (CHM) to be obtained at 1 m spatial resolution. The dominant height of all the samples was estimated by the ONF at plot level with Computree software (v5.0) [START_REF] Krebs | [END_REF][START_REF] Piboule | Computree: A Collaborative Platform for Use of Terrestrial Lidar in Dendrometry[END_REF] using a crown segmentation approach on the CHM, and then, the heights of the 6 highest trees in the plot were extracted to obtain the dominant height. Results have been validated and corrected by the ONF with the samples where in situ height measurements were made. In the end, at plot level, the values of the basal area (BA), mean diameter at breast height (DBH), tree density (Density) and dominant height (Hdom) are available.

density; 6 and 15 cm vertical accuracy; and 20 and 30 cm absolute positional accuracy. ALS data enabled a canopy height model (CHM) to be obtained at 1 m spatial resolution. The dominant height of all the samples was estimated by the ONF at plot level with Computree software (v5.0) [START_REF] Krebs | [END_REF][START_REF] Piboule | Computree: A Collaborative Platform for Use of Terrestrial Lidar in Dendrometry[END_REF] using a crown segmentation approach on the CHM, and then, the heights of the 6 highest trees in the plot were extracted to obtain the dominant height. Results have been validated and corrected by the ONF with the samples where in situ height measurements were made. In the end, at plot level, the values of the basal area (BA), mean diameter at breast height (DBH), tree density (Density) and dominant height (Hdom) are available. We selected two classes for each site in order to work on consistent sets of samples (i.e., similar allometric relationships within each class). The classes are based on the dominant species that can be found in the forest cover cartographic database. At the Orléans test site, the "oak" class contains 75 plots with at least 75% of the basal area being oaks; secondary species are mainly hornbeams and pines. The "Scots pine" class contains 28 plots with at least 75% of the basal area being Scots pines; secondary species are mainly oaks and other pines. At the StGobain test site, the "oak" class contains 21 plots (still 75% basal area threshold); the secondary species is manly hornbeam. The "beech" class contains 25 plots (50% basal area threshold); secondary species are mainly oaks, hornbeams and ash (Fraxinus excelsior). Field inventory values are summarised in Table 1. The forest structure statistics for each class are broadly similar except for the dominant height which is higher at the StGobain site. We selected two classes for each site in order to work on consistent sets of samples (i.e., similar allometric relationships within each class). The classes are based on the dominant species that can be found in the forest cover cartographic database. At the Orléans test site, the "oak" class contains 75 plots with at least 75% of the basal area being oaks; secondary species are mainly hornbeams and pines. The "Scots pine" class contains 28 plots with at least 75% of the basal area being Scots pines; secondary species are mainly oaks and other pines. At the StGobain test site, the "oak" class contains 21 plots (still 75% basal area threshold); the secondary species is manly hornbeam. The "beech" class contains 25 plots (50% basal area threshold); secondary species are mainly oaks, hornbeams and ash (Fraxinus excelsior). Field inventory values are summarised in Table 1. The forest structure statistics for each class are broadly similar except for the dominant height which is higher at the StGobain site. Sentinel-2 images: Sentinel-2 data consist of optical time series acquired in the visible, near infrared (nir) and short-wave infrared (swir), with a spatial resolution of 10 to 20 m. The raw images are provided by the European Space Agency (ESA). The Theia Land Data Center [78] produces the level 2A images (top of canopy reflectance: TOC) with orthorectification, atmospheric corrections and cloud mask, using MAJA processing [START_REF] Hagolle | A Multi-Temporal Method for Cloud Detection, Applied to FORMOSAT-2, VENµS, LANDSAT and SENTINEL-2 Images[END_REF][START_REF] Hagolle | A Multi-Temporal and Multi-Spectral Method to Estimate Aerosol Optical Thickness over Land, for the Atmospheric Correction of FormoSat-2, LandSat, VENµS and Sentinel-2 Images[END_REF]. For both sites, we selected one clear acquisition in winter (March 2016) and one clear acquisition in summer (early September 2016). Spectral indices were extracted from Sentinel-2 bands (Equations ( 1)-(3): Normalized Difference Vegetation Index (NDVI), Brightness index (BI), and Normalised Difference Water Index (NDWI).

NDVI = B8 -B4 B8 + B4 (1)
BI = B2 2 + B3 2 + B4 2 (2) NDWI = B8A -B11 B8A + B11 (3)
where B2, B3, B4 and B8 are, respectively, blue, green, red and nir spectral bands at 10 m spatial resolution; B8A and B11 are, respectively, nir and swir at 20 m spatial resolution. NDWI is then resampled to 10 m spatial resolution. These indices from winter and summer acquisitions were used in the following processing.

Sentinel-1: C-band SAR time series resampled at 10 m spatial resolution are provided by the ESA and downloaded from the Sentinel Product Exploitation Platform [START_REF] Peps | Plateforme d'exploitation des produits Sentinel[END_REF] (PEPS). We used the Level-1 Ground Range Detected (GRD) products on the descending orbit (DES, morning acquisition) with ~33 • incidence angle (IA) on the test sites. The 2016 Sentinel-1 time series were calibrated, orthorectified and corrected from slope effect, and then, we applied multi-image filtering [START_REF] Bruniquel | Multi-Variate Optimal Speckle Reduction in SAR Imagery[END_REF][START_REF] Quegan | Filtering of Multichannel SAR Images[END_REF] on the 2016 time series, leading to an equivalent number of looks (ENL) of ~92 for each acquisition. Then, average images for VH and VV polarisations were calculated using all the dates included in an interval of 2 months, centred on the Sentinel-2 winter and summer dates. VH and VV winter and summer averages were used in the following processing.

ALOS-2 PALSAR-2: L-band SAR annual mosaics at 25 m spatial resolution are produced by the Japan Aerospace Exploration Agency (JAXA) through large scale mosaicking [START_REF] Shimada | Generating Large-Scale High-Quality SAR Mosaic Datasets: Application to PALSAR Data for Global Monitoring[END_REF] that includes orthorectification, slope correction and radiometric calibration between neighbouring strips. Five annual mosaics were sourced from JAXA for the years 2015 through 2019. We converted digital numbers (DN) to gamma naught values. Then, we applied multi-image filtering [START_REF] Bruniquel | Multi-Variate Optimal Speckle Reduction in SAR Imagery[END_REF][START_REF] Quegan | Filtering of Multichannel SAR Images[END_REF] using the 2 polarisations and the 5 years, leading to an ENL of 135 for each image. In the following, we used resulting HV and HH images from 2016.

Textural indices: spatial texture metrics have been generated from both Sentinel-2 10 m spatial resolution indices (NDVI, BI) and Sentinel-1 polarisations (VH, VV), using a Grey Level Co-occurrence Matrix (GLCM) with the Orfeo Toolbox library (OTB, [START_REF] Grizonnet | Orfeo ToolBox: Open Source Processing of Remote Sensing Images[END_REF]). Three second order Haralick metrics were calculated on GLCM: homogeneity (also called Inverse Difference Moment), contrast (also called Inertia), and Haralick correlation. Based on a previous sensitivity analysis [START_REF] Morin | Estimation et Suivi de la Ressource en Bois en France Métropolitaine par Valorisation des Séries Multi-Temporelles à Haute Résolution Spatiale D'images Optiques[END_REF], the following parameters were used to obtain the textural metrics for further processing: offset set to 1 pixel; 3 pixels radius; 100 grey levels (number of bins); minimum and maximum calculated separately between Sentinel-2 indices and dates, and between Sentinel-1 polarisations; and orientations 0 • , 45 • , 90 • and 135 • were computed then averaged.

Final set of features: after the processing steps explained above, we obtained a total of 36 features: 16 Sentinel-1 features (VH and VV for winter and summer dates and derived textural metrics), 18 Sentinel-2 features (BI, NDVI and NDWI for winter and summer dates, BI and NDVI textural metrics), and 2 ALOS-2 PALSAR-2 metrics (HV and HH polarisations).

GEDI Data

The Global Ecosystem Dynamics Investigation (GEDI) LiDAR is a laser instrument on board the International Space Station. It has three lasers; two of the lasers are full power and the other one is split into two beams (called coverage beams), producing a total of four beams and eight ground tracks (four power and four cover tracks). Footprints average 25 m in diameter and they are separated by 60 m along track and 600 across track.

The GEDI L2A data from May 2019 to August 2020 were downloaded [START_REF] Dubayah | The Global Ecosystem Dynamics Investigation: High-Resolution Laser Ranging of the Earth's Forests and Topography[END_REF] and preprocessed as presented in [START_REF] Fayad | Estimating Canopy Height and Wood Volume of Eucalyptus Plantations in Brazil Using GEDI LiDAR Data[END_REF]. On this data, the first settings group is available to derive relative height metrics (RH): Smoothing Width (Noise) = 6.5σ; Smoothing Width (Signal) = 6.5σ; Waveform Signal Start Threshold = 3σ; Waveform Signal End Threshold = 6σ. σ represents the standard deviation of the background noise level. Relative height metrics (RH n ) were available with n from 0 to 100% (with a step of 1%). In order to select the best RH n , we compared different GEDI RH percentiles to the airborne LiDAR (ALS) data available on the two test sites. Based on the results (see Section 3.2) and previous studies, [START_REF] Dorado-Roda | Assessing the Accuracy of GEDI Data for Canopy Height and Aboveground Biomass Estimates in Mediterranean Forests[END_REF][START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF][START_REF] Khati | The Role of Time-Series L-Band SAR and GEDI in Mapping Sub-Tropical Above-Ground Biomass[END_REF] we chose to use RH95 as the GEDI dominant height estimates. We also tested the relevance of using all the beams or only the full power beams.

After filtering the data as performed in [START_REF] Fayad | Estimating Canopy Height and Wood Volume of Eucalyptus Plantations in Brazil Using GEDI LiDAR Data[END_REF] (heights not consistent with tree height), we considered different options to reduce GEDI height metric errors. We kept only waveforms where at least 2 modes were detected to ensure a top of canopy and soil return signal. Finally, we tested to remove the noisiest waveforms, keeping only the 25, 50 or 75% waveforms with the best signal-noise ratio.

Methodology 2.3.1. Processing of Ancillary Products

The forest species on the GEDI samples were unknown; therefore, we used a combination of two datasets for tree species classification: the OSO land cover map (website: https://www.theia-land.fr/en/product/land-cover-map/, accessed on 16 March 2022) and BDForêtV2 (website: https://geoservices.ign.fr, accessed on 16 March 2022).

OSO: annual land cover maps are available for France on the Theia platform and produced with IOTA2 software (website: https://framagit.org/iota2-project/iota2, accessed on 16 March 2022); the main stages are described in [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. The nomenclature includes two classes of forests: broadleaved and coniferous. This nomenclature is not sufficient for our study, but this map offers a high spatial accuracy (10 × 10 m pixels) and up-to-date land cover information (annual product). The 2018 OSO map is used.

BDForetV2: the French national institute of geographic and forest information (IGN) has produced forest species maps using airborne optical images from 2004 to 2016, depend-ing on localisation. Over the two sites, Orléans and StGobain, forest maps were produced, respectively, with 2010 and 2016 images. The nomenclature provides knowledge of the location of the pure species (>75% of the surface is occupied by a given species), which is enough to assign tree species classes to GEDI footprints. However, the spatial accuracy, the updating and the minimum mapping unit are somewhat limited for our use.

Forest maps combination: for each GEDI footprint, we assigned a species (or group of species) class when the OSO and BDForetV2 agreed (e.g., OSO indicates broadleaved and BDForetV2 indicates oak). We kept only footprints with pure species.

In the end, we had 826 GEDI oak dominant footprints and 1701 Scots pine dominant footprints in Orléans, and 1195 oak dominant footprints and 1090 beech dominant footprints in StGobain.

Estimation of Forest Parameters

Our processing chain is based on a machine learning approach; the full methodology is described in detail in [START_REF] Morin | Estimation and Mapping of Forest Structure Parameters from Open Access Satellite Images: Development of a Generic Method with a Study Case on Coniferous Plantation[END_REF][START_REF] Morin | Estimation et Suivi de la Ressource en Bois en France Métropolitaine par Valorisation des Séries Multi-Temporelles à Haute Résolution Spatiale D'images Optiques[END_REF]. An example with our number of features and a forward approach is schematised in Figure 2. The key points are:

•

ONF Hdom measurements or GEDI RH95% are used as reference data, and the remote sensing features from Section 2.3 are used as predictive variables.

•

Reference data and remote sensing features are fed to a Support Vector Machine (SVR) algorithm for regression; the SVR has 3 key parameters that need to be optimised: cost parameter, gamma and epsilon.

•

Validation is made with a leave-one-out (LOO) cross-validation method when using the ONF measurement data because the number of samples is very low for these data (21 to 75 samples); when the reference data is GEDI RH95, we use a stratified K-fold cross-validation method.

•

We used a stepwise approach for feature selection. It is an iterative process that starts with 1 feature or N features and iteratively adds or removes features, respectively, for forward and backward selection. Within the 1 to N feature combinations, the one with the lowest RMSE is considered the best predictive set we can obtain on the dataset.

BDForetV2: the French national institute of geographic and forest information (IGN) has produced forest species maps using airborne optical images from 2004 to 2016, depending on localisation. Over the two sites, Orléans and StGobain, forest maps were produced, respectively, with 2010 and 2016 images. The nomenclature provides knowledge of the location of the pure species (>75% of the surface is occupied by a given species), which is enough to assign tree species classes to GEDI footprints. However, the spatial accuracy, the updating and the minimum mapping unit are somewhat limited for our use.

Forest maps combination: for each GEDI footprint, we assigned a species (or group of species) class when the OSO and BDForetV2 agreed (e.g., OSO indicates broadleaved and BDForetV2 indicates oak). We kept only footprints with pure species.

In the end, we had 826 GEDI oak dominant footprints and 1701 Scots pine dominant footprints in Orléans, and 1195 oak dominant footprints and 1090 beech dominant footprints in StGobain.

Estimation of Forest Parameters

Our processing chain is based on a machine learning approach; the full methodology is described in detail in [START_REF] Morin | Estimation and Mapping of Forest Structure Parameters from Open Access Satellite Images: Development of a Generic Method with a Study Case on Coniferous Plantation[END_REF][START_REF] Morin | Estimation et Suivi de la Ressource en Bois en France Métropolitaine par Valorisation des Séries Multi-Temporelles à Haute Résolution Spatiale D'images Optiques[END_REF]. An example with our number of features and a forward approach is schematised in Figure 2. The key points are:

•
ONF Hdom measurements or GEDI RH95% are used as reference data, and the remote sensing features from Section 2.3 are used as predictive variables.

•

Reference data and remote sensing features are fed to a Support Vector Machine (SVR) algorithm for regression; the SVR has 3 key parameters that need to be optimised: cost parameter, gamma and epsilon.

•

Validation is made with a leave-one-out (LOO) cross-validation method when using the ONF measurement data because the number of samples is very low for these data (21 to 75 samples); when the reference data is GEDI RH95, we use a stratified K-fold cross-validation method.

•

We used a stepwise approach for feature selection. It is an iterative process that starts with 1 feature or N features and iteratively adds or removes features, respectively, for forward and backward selection. Within the 1 to N feature combinations, the one with the lowest RMSE is considered the best predictive set we can obtain on the dataset. In this study, three ways to estimate forest parameters from satellite imagery are tested. First, is the estimation of BA, DBH, Density and Hdom using ONF field In this study, three ways to estimate forest parameters from satellite imagery are tested. First, is the estimation of BA, DBH, Density and Hdom using ONF field measurement for learning and validation (LOO cross-validation). Second, is the estimation of Hdom using GEDI RH95 for learning and validation (K-fold cross-validation). Finally, the estimation of Hdom using GEDI RH95 for learning and ONF field measurements for validation is conducted.

The statistics used to evaluate performances are Pearson's correlation (r, r 2 and p0), RMSE, relative RMSE, Mean Absolute Error (MAE) and relative MAE.

RMSE = 1:N ∑ (reference -prediction) 2 N (4) MAE = 1:N ∑ |reference -prediction| N ( 5
)
where N is the number of samples. Relative scores (in percentages) are these statistics divided by the mean of reference data values. RMSE gives more weight to outliers compared to MAE.

Results

GLAD-2019 Canopy Height Map Evaluation

A global map of forest canopy height at 30 m spatial resolution was produced by [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF], based on GEDI RH95% and Landsat-8 2019 times series. We evaluated the accuracy of this map over our test sites in temperate forests (France).

Figure 3 shows the scatterplot between the ONF field measurements (x-axis) based on dominant trees and the GLAD-2019 forest canopy height map (y-axis) over the four forest classes. The scatterplots show poor agreement between the two. The dynamic range of the GLAD-2019 values are restricted compared with the field data in the StGobain forest, and multiple outliers are found in the Orléans test site. Several possible causes can explain the poor agreement for these mixed forests (to silvicultural interventions at Orléans or to atypical forest stands, etc.). More importantly, the GLAD-2019 height estimates were not able to retrieve the forest height dynamic range of the reference data, as indicated by the low correlation coefficients. For the Orléans oak class, reference values range from 14 to 30 m, while the GLAD-2019 height values range from ~19 to ~25 m. The dynamics of the GLAD-2019 values are even more restricted for other classes compared to reference value dynamics. This is due both to the type of forests that we observed which are mixed (different ages, diameters and species within the plots) and spatially complex, and to the satellite image data used for height map production: Landsat-8 time series are optical data; several studies showed that they have a significant but limited relationship with forest parameters compared to radar and textural metrics [START_REF] Freitas | Relationships between Forest Structure and Vegetation Indices in Atlantic Rainforest[END_REF][START_REF] Dube | Investigating the Robustness of the New Landsat-8 Operational Land Imager Derived Texture Metrics in Estimating Plantation Forest Aboveground Biomass in Resource Constrained Areas[END_REF][START_REF] Zhu | Improving Forest Aboveground Biomass Estimation Using Seasonal Landsat NDVI Time-Series[END_REF].

measurement for learning and validation (LOO cross-validation). Second, is the estimation of Hdom using GEDI RH95 for learning and validation (K-fold cross-validation). Finally, the estimation of Hdom using GEDI RH95 for learning and ONF field measurements for validation is conducted.

The statistics used to evaluate performances are Pearson's correlation (r, r 2 and p0), RMSE, relative RMSE, Mean Absolute Error (MAE) and relative MAE.

1:N 2 (reference prediction) RMSE N - =  (4) 1:N | reference prediction | MAE N - =  (5)
where N is the number of samples. Relative scores (in percentages) are these statistics divided by the mean of reference data values. RMSE gives more weight to outliers compared to MAE.

Results

GLAD-2019 Canopy Height Map Evaluation

A global map of forest canopy height at 30 m spatial resolution was produced by [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF], based on GEDI RH95% and Landsat-8 2019 times series. We evaluated the accuracy of this map over our test sites in temperate forests (France).

Figure 3 shows the scatterplot between the ONF field measurements (x-axis) based on dominant trees and the GLAD-2019 forest canopy height map (y-axis) over the four forest classes. The scatterplots show poor agreement between the two. The dynamic range of the GLAD-2019 values are restricted compared with the field data in the StGobain forest, and multiple outliers are found in the Orléans test site. Several possible causes can explain the poor agreement for these mixed forests (to silvicultural interventions at Orléans or to atypical forest stands, etc.). More importantly, the GLAD-2019 height estimates were not able to retrieve the forest height dynamic range of the reference data, as indicated by the low correlation coefficients. For the Orléans oak class, reference values range from 14 to 30 m, while the GLAD-2019 height values range from ~19 to ~25 m. The dynamics of the GLAD-2019 values are even more restricted for other classes compared to reference value dynamics. This is due both to the type of forests that we observed which are mixed (different ages, diameters and species within the plots) and spatially complex, and to the satellite image data used for height map production: Landsat-8 time series are optical data; several studies showed that they have a significant but limited relationship with forest parameters compared to radar and textural metrics [START_REF] Freitas | Relationships between Forest Structure and Vegetation Indices in Atlantic Rainforest[END_REF][START_REF] Dube | Investigating the Robustness of the New Landsat-8 Operational Land Imager Derived Texture Metrics in Estimating Plantation Forest Aboveground Biomass in Resource Constrained Areas[END_REF][START_REF] Zhu | Improving Forest Aboveground Biomass Estimation Using Seasonal Landsat NDVI Time-Series[END_REF].

Dominant Height Estimation Using Field Campaign as Reference Data

The processing chain described in Section 2.3.2 and [START_REF] Morin | Estimation and Mapping of Forest Structure Parameters from Open Access Satellite Images: Development of a Generic Method with a Study Case on Coniferous Plantation[END_REF] was used to estimate the dominant height (Hdom), basal area (BA), mean diameter at breast height (DBH) and tree density (Density) on the ONF dataset for the forest classes in the Orléans and StGobain test sites. In this section, we focus on the Hdom results. BA, DBH and Density estimations are presented in Appendix A (Figures A1 andA2, Table A1) for informational purposes: the overall results on the broadleaved and Scots pines forests are good, with RMSE from 11% for dominant height to 30% for tree Density. The results by class show slightly less accurate predictions on the Orleans oaks for all the forest parameter estimates; the other classes have better and fairly similar results.

Using only in situ data (ONF field measurements for model training and validation with a LOO cross-validation), a forward stepwise feature selection process was applied, since the number of reference samples was low compared to the number of input features in the learning model. Figure 4 presents the evolution of the dominant height validation statistics (r 2 ) during the forward feature selection process. The r 2 first increases for all the four forest classes until it reaches a maximum (i.e., best combination) of between 5 and 15 features. The features selected are different for each class. At least two textural indices from Sentinel-1 and two others from Sentinel-2 are selected in the optimal combination, even if they differ between classes (sensor indices, textural indices or acquisition date). Spectral indices and L-or C-band polarisations are used even more variably. After reaching the maximum, precision decreases and shows poor estimation quality with all features in the model. This may be due to the low number of samples compared to the maximum number of features; there is over-fitting with 36 features and therefore the cross-validation results show large errors. This result shows the importance of feature selection when using small numbers of samples; forward stepwise selection allows very good results to be obtained where the total feature set causes over-fitting. However, the selection process can also cause excessive method specialisation and therefore difficulties in generalising the results.

Figure 5 presents the validation scatterplots of the four classes for dominant height estimation with the ONF samples as references for training and testing using a LOO crossvalidation. Good results are obtained for StGobain and for Scots pine in Orléans with 0.74-0.78 r 2 and low errors (1.9-3.1 m RMSE, 7.3-11.4% relative RMSE). For Orléans oak forest, the relative RMSE is still satisfactory (2.6 RMSE, 11.6% relative RMSE), but the r 2 (0.48) is not as good as for the other classes. There might be a saturation effect; dominant height beyond 25 m is underestimated.

Dominant Height Estimation Using Field Campaign as Reference Data

The processing chain described in Section 2.3.2 and [START_REF] Morin | Estimation and Mapping of Forest Structure Parameters from Open Access Satellite Images: Development of a Generic Method with a Study Case on Coniferous Plantation[END_REF] was used to estimate the dominant height (Hdom), basal area (BA), mean diameter at breast height (DBH) and tree density (Density) on the ONF dataset for the forest classes in the Orléans and StGobain test sites. In this section, we focus on the Hdom results. BA, DBH and Density estimations are presented in Appendix A (Figures A1 andA2, Table A1) for informational purposes: the overall results on the broadleaved and Scots pines forests are good, with RMSE from 11% for dominant height to 30% for tree Density. The results by class show slightly less accurate predictions on the Orleans oaks for all the forest parameter estimates; the other classes have better and fairly similar results.

Using only in situ data (ONF field measurements for model training and validation with a LOO cross-validation), a forward stepwise feature selection process was applied, since the number of reference samples was low compared to the number of input features in the learning model. Figure 4 presents the evolution of the dominant height validation statistics (r 2 ) during the forward feature selection process. The r 2 first increases for all the four forest classes until it reaches a maximum (i.e., best combination) of between 5 and 15 features. The features selected are different for each class. At least two textural indices from Sentinel-1 and two others from Sentinel-2 are selected in the optimal combination, even if they differ between classes (sensor indices, textural indices or acquisition date). Spectral indices and L-or C-band polarisations are used even more variably. After reaching the maximum, precision decreases and shows poor estimation quality with all features in the model. This may be due to the low number of samples compared to the maximum number of features; there is over-fitting with 36 features and therefore the cross-validation results show large errors. This result shows the importance of feature selection when using small numbers of samples; forward stepwise selection allows very good results to be obtained where the total feature set causes over-fitting. However, the selection process can also cause excessive method specialisation and therefore difficulties in generalising the results.

Figure 5 presents the validation scatterplots of the four classes for dominant height estimation with the ONF samples as references for training and testing using a LOO crossvalidation. Good results are obtained for StGobain and for Scots pine in Orléans with 0.74-0.78 r 2 and low errors (1.9-3.1 m RMSE, 7.3-11.4% relative RMSE). For Orléans oak forest, the relative RMSE is still satisfactory (2.6 RMSE, 11.6% relative RMSE), but the r 2 (0.48) is not as good as for the other classes. There might be a saturation effect; dominant height beyond 25 m is underestimated.

Dominant Height Estimation Using GEDI Height Metrics as Reference Data

As it has been discussed, the processing chain is able to estimate forest dominant height from a set of satellite features. But the main limitation in the performance results is the small number of samples for model learning. We are thus interested in the use of GEDI data to increase the number of reference samples for the height predictor. In this section, we first present the results of our analysis of the GEDI height metrics compared to ALS data. Secondly, we present the application of our processing chain to estimate dominant height using GEDI as reference data and the satellite image features.

GEDI RH95% Comparison with ALS Data

GEDI RH percentiles were compared to tree height metrics extracted from the ALS canopy height model (CHM) at 1 m spatial resolution. CHM metrics were calculated on the pixels inside the footprints of the GEDI sample (circles 25 m in diameter). Among different CHM metrics, we selected the 95% percentile (H95) as a proxy for dominant height [START_REF] López | Validation de L'estimation de la Hauteur des Jeunes Peuplements à Partir de Données LiDAR Aériennes[END_REF][START_REF] Lucie | Comparaison de Modèles Numériques de Surface photogrammétriques de différentes résolutions en forêt mixte. estimation d'une variable dendrométrique simple: La hauteur dominante[END_REF], and compared it to GEDI RH metrics. The analysis of GEDI RHn with n from 50 to 100% (steps (50,75,85,90,95,99,100)) revealed that RH95% was the best metric overall (according to r 2 and RMSEs of the four classes studied).

Figure 6 shows the comparison between CHM-H95 (x-axis) and GEDI RH95% (yaxis) for the four classes. The best result was obtained on the only coniferous class (Scots pine), with r 2 equal to 0.65 and 3.7 m RMSE. Broadleaved classes have lower coefficients of determination of around 0.48. In France, coniferous forests are generally managed as uniform high stands and thus easier to characterise. We can observe outliers overestimated in the low range and underestimated in the high range. This can be caused partly by growth or silvicultural interventions in the period between the ALS and GEDI acquisitions, or the acquisition date (GEDI winter dates can underestimate height due to leaf loss), or by not enough filtering of GEDI footprints. For the purposes of the study and to use the GEDI samples in the best conditions, outliers were removed before the application of our processing chain for forest parameter estimation with satellite image features. The threshold was set at 5 m away from the linear regression line between CHM-H95 and GEDI RH95% for each class. In the end, we had 479 GEDI oak dominant footprints and 958 Scots pine dominant footprints in Orléans, and 808 oak dominant footprints and 723 beech dominant footprints in StGobain. Apart from outliers, there is generally a bias or saturation for high values, more visible in StGobain because the heights are higher.

Dominant Height Estimation Using GEDI Height Metrics as Reference Data

As it has been discussed, the processing chain is able to estimate forest dominant height from a set of satellite features. But the main limitation in the performance results is the small number of samples for model learning. We are thus interested in the use of GEDI data to increase the number of reference samples for the height predictor. In this section, we first present the results of our analysis of the GEDI height metrics compared to ALS data. Secondly, we present the application of our processing chain to estimate dominant height using GEDI as reference data and the satellite image features.

GEDI RH95% Comparison with ALS Data

GEDI RH percentiles were compared to tree height metrics extracted from the ALS canopy height model (CHM) at 1 m spatial resolution. CHM metrics were calculated on the pixels inside the footprints of the GEDI sample (circles 25 m in diameter). Among different CHM metrics, we selected the 95% percentile (H95) as a proxy for dominant height [START_REF] López | Validation de L'estimation de la Hauteur des Jeunes Peuplements à Partir de Données LiDAR Aériennes[END_REF][START_REF] Lucie | Comparaison de Modèles Numériques de Surface photogrammétriques de différentes résolutions en forêt mixte. estimation d'une variable dendrométrique simple: La hauteur dominante[END_REF], and compared it to GEDI RH metrics. The analysis of GEDI RHn with n from 50 to 100% (steps (50,75,85,90,95,99,100)) revealed that RH95% was the best metric overall (according to r 2 and RMSEs of the four classes studied).

Figure 6 shows the comparison between CHM-H95 (x-axis) and GEDI RH95% (yaxis) for the four classes. The best result was obtained on the only coniferous class (Scots pine), with r 2 equal to 0.65 and 3.7 m RMSE. Broadleaved classes have lower coefficients of determination of around 0.48. In France, coniferous forests are generally managed as uniform high stands and thus easier to characterise. We can observe outliers overestimated in the low range and underestimated in the high range. This can be caused partly by growth or silvicultural interventions in the period between the ALS and GEDI acquisitions, or the acquisition date (GEDI winter dates can underestimate height due to leaf loss), or by not enough filtering of GEDI footprints. For the purposes of the study and to use the GEDI samples in the best conditions, outliers were removed before the application of our processing chain for forest parameter estimation with satellite image features. The threshold was set at 5 m away from the linear regression line between CHM-H95 and GEDI RH95% for each class. In the end, we had 479 GEDI oak dominant footprints and 958 Scots pine dominant footprints in Orléans, and 808 oak dominant footprints and 723 beech dominant footprints in StGobain. Apart from outliers, there is generally a bias or saturation for high values, more visible in StGobain because the heights are higher.

Using GEDI RH95% Instead of Field Data

In this section, we present the results of the processing chain presented in Section 2.3.2 to estimate the dominant height using GEDI RH95% as reference samples in Orléans and StGobain test sites (GEDI RH95% for model training and validation using a K-fold cross-validation). Figure 7 shows the evolution of r 2 on the four forest classes with the feature selection process. The r 2 is stable from 36 features (maximum) to ~15 features. Having a large number of samples makes it possible to have "unnecessary" features in the model without affecting the results. Among these ~15 features, the variability between classes is lower than when working with the small ONF datasets (all feature types are always present in the 15 features). However, differences remain when looking at precisely which features are selected for each class. When this threshold of 15 features is exceeded, the precision starts to fall faster and faster when there are a very few numbers of features in the model. The models no longer have enough information to allow a good estimation of the dominant height. The feature selection step might not be necessary when working with GEDI samples. Moreover, directly using all the features without specific selection for each class is better in order to be more easily generalisable. Indeed, working by class implies having cartographic information on these classes to assign them to the GEDI footprints, then to produce maps with the prediction models. This cartographic information is not always accurate, it is therefore preferable to have identical feature sets between classes to avoid having artefacts when there are species errors.

Figures 891011present the scatterplots of the estimation of Hdom using one model per forest class (algorithm parametrisation and model learning are made separately for each forest class). The graphs on the left side show GEDI RH95% cross-validation results, meaning the ability of satellite images to retrieve GEDI height metrics. The graphs on the right side show the application of GEDI RH95% based models on ONF measurement plots, meaning an independent validation of the dominant height predictions using GEDI data and satellite image features. GEDI RH95% K-fold cross-validation results are roughly similar for the four classes, just slightly worse on Orléans oak. Statistics are good, showing r 2 from 0.54 to 0.65, and RMSEs from 2.7 to 3.4 m (~12.9 to 16.7%). Despite the results being scattered, there is no significant bias and the satellite image features are able to retrieve the range of GEDI RH95% value dynamics. Learning on GEDI data and testing on the ONF measurement plots also provide good overall results (right graph on Figures 891011). Validation statistics are however more variable, partly because of the lower numbers of samples and, less importantly, the larger variability in the ONF value dynamics between

Using GEDI RH95% Instead of Field Data

In this section, we present the results of the processing chain presented in Section 2.3.2 to estimate the dominant height using GEDI RH95% as reference samples in Orléans and StGobain test sites (GEDI RH95% for model training and validation using a K-fold crossvalidation). Figure 7 shows the evolution of r 2 on the four forest classes with the feature selection process. The r 2 is stable from 36 features (maximum) to ~15 features. Having a large number of samples makes it possible to have "unnecessary" features in the model without affecting the results. Among these ~15 features, the variability between classes is lower than when working with the small ONF datasets (all feature types are always present in the 15 features). However, differences remain when looking at precisely which features are selected for each class. When this threshold of 15 features is exceeded, the precision starts to fall faster and faster when there are a very few numbers of features in the model. The models no longer have enough information to allow a good estimation of the dominant height. The feature selection step might not be necessary when working with GEDI samples. Moreover, directly using all the features without specific selection for each class is better in order to be more easily generalisable. Indeed, working by class implies having cartographic information on these classes to assign them to the GEDI footprints, then to produce maps with the prediction models. This cartographic information is not always accurate, it is therefore preferable to have identical feature sets between classes to avoid having artefacts when there are species errors.

Figures 891011present the scatterplots of the estimation of Hdom using one model per forest class (algorithm parametrisation and model learning are made separately for each forest class). The graphs on the left side show GEDI RH95% cross-validation results, meaning the ability of satellite images to retrieve GEDI height metrics. The graphs on the right side show the application of GEDI RH95% based models on ONF measurement plots, meaning an independent validation of the dominant height predictions using GEDI data and satellite image features. GEDI RH95% K-fold cross-validation results are roughly similar for the four classes, just slightly worse on Orléans oak. Statistics are good, showing r 2 from 0.54 to 0.65, and RMSEs from 2.7 to 3.4 m (~12.9 to 16.7%). Despite the results being scattered, there is no significant bias and the satellite image features are able to retrieve the range of GEDI RH95% value dynamics. Learning on GEDI data and testing on the ONF measurement plots also provide good overall results (right graph on Figures 891011). Validation statistics are however more variable, partly because of the lower numbers of samples and, less importantly, the larger variability in the ONF value dynamics between classes. The Orléans oak class has the worst r 2 (0.23), but the error and bias are still satisfactory. The best result is on the StGobain oak class, with 0.64 r 2 and 12.8% relative RMSE.

classes. The Orléans oak class has the worst r 2 (0.23), but the error and bias are still satisfactory. The best result is on the StGobain oak class, with 0.64 r 2 and 12.8% relative RMSE. classes. The Orléans oak class has the worst r 2 (0.23), but the error and bias are still satisfactory. The best result is on the StGobain oak class, with 0.64 r 2 and 12.8% relative RMSE.

Discussion

In a previous study, we built and validated a processing chain on forest coniferous plantations. In the study presented in this paper, our processing chain was able to accurately estimate forest parameters on more complex temperate semi-natural forests (uneven ages, mixed species, spatial variability). All the results obtained with the ONF data, GEDI data and the comparison with GLAD-2019 height map are summarised in Table 2.

Discussion

In a previous study, we built and validated a processing chain on forest coniferous plantations. In the study presented in this paper, our processing chain was able to accurately estimate forest parameters on more complex temperate semi-natural forests (uneven ages, mixed species, spatial variability). All the results obtained with the ONF data, GEDI data and the comparison with GLAD-2019 height map are summarised in Table 2.

The best results we obtained are relative RMSEs from 7.3% to 11.6% for dominant height estimation using the ONF measurements as a reference dataset in a LOO crossvalidation. The error rates obtained in this study with satellite imagery are similar to the results seen in the literature with photogrammetry and radargrammetry (6.6% to 14% relative RMSE depending on test sites), and slightly higher than the results seen with ALS data (4.6% to 7.8% relative RMSE) [START_REF] Yu | Comparison of Laser and Stereo Optical, SAR and InSAR Point Clouds from Air-and Space-Borne Sources in the Retrieval of Forest Inventory Attributes[END_REF][START_REF] Holopainen | Forest inventory attribute estimation using airborne laser scanning, aerial stereo imagery, radargrammetry and interferometry-finnish experiences of the 3d techniques[END_REF][START_REF] Vastaranta | Airborne Laser Scanning and Digital Stereo Imagery Measures of Forest Structure: Comparative Results and Implications to Forest Mapping and Inventory Update[END_REF][START_REF] Pearse | Comparison of High-Density LiDAR and Satellite Photogrammetry for Forest Inventory[END_REF]. Thereby, the combination of optical and SAR sensors with textural indices at 10 m spatial resolution, in a forward stepwise feature selection process, allowed accurate models for different coniferous and broadleaved forests to be built. Satellite images used in this study are freely available within long-term programs, acquired globally and regularly at high spatial resolution (10-25 m). It is very promising to be able to use these data to estimate forest parameters that can be used for forest management or carbon modelling. However, there are limitations to generalising the method. Field data that meet the needs of mapping by satellite imagery (precise geolocation, the distribution of values, the size of the sampled surfaces, etc.) are often of limited number and concentrated in small areas. The accuracy and the ability to generalise the results are therefore limited by field measurement availability. In this study, with small measurement datasets, we could produce accurate estimations by using an automated feature selection process (forward stepwise selection) and showed that this selection was important to improve the results. However, a selection of features, whether automatic or manual (operator analysis), also increases the risk of being unable to obtain high accuracy on other types of forest or when dealing with large areas (regional or country scale), with forests not represented in the reference data. GEDI data availability therefore presents a huge opportunity for methods based on machine learning. The data acquisition will be carried out only until 2023, but the coverage and the large number of samples have already the potential to be integrated into learning approaches for the mapping of local forests worldwide. We have shown that the heights provided by GEDI are in good agreement with ALS heights. When GEDI RH95% is compared with ALS-CHM 95% (Figure 6, Section 3.3.1), r 2 ~0.48 and RMSE ranges from 5.1 to 5.6 m for the broadleaved classes, and r 2 is 0.65 and RMSE = 3.7 m for the coniferous class (Orléans Scots pines). These statistics are improved when removing outliers that may be due to time shift, silvicultural interventions and acquisition conditions (r 2 = 0.78 and RMSE = 4.5 m for the Orléans oak class, r 2 = 0.84-0.86 and RMSE = 2.8-2.9 m for the other classes). The lower accuracy on the Orléans oak class can be explained by a larger presence of coppice forests and open forests in this class, while other classes are more dense forests.

A global map of forest canopy height has been produced by [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF] using GEDI RH95% and Landsat-8 time series. Yet, the validation on our test sites showed very poor results, making it necessary to produce new estimates of forest height that are accurate enough at local scale to allow forest management and monitoring of these temperate complex forests. We have shown that our processing chain using more diversified satellite images and textural indices could produce good estimates of the dominant height with GEDI data as reference samples for model learning. Moreover, the quantity of GEDI samples allows the use of all input features together without a feature selection process, which will facilitate the generalisation of the method. The GEDI RH95% cross-validation showed no saturation; r 2 from 0.54 to 0.65 and RMSE from 2.7 m to 3.4 m (12.9 to 16.7% relative RMSE). The validation of the GEDI's predictions on the measurement plots showed roughly equivalent errors. This is a very good result given that GEDI RH95% comparison to ALS data had RMSE from 3.7 to 5.6 m (Figure 6). Furthermore, the comparison of the map produced by our method and the CHM-H95 at 10 m spatial resolution (Table 2, bottom column) shows that the prediction values are similar or even better than the GEDI RH95% values. Bias and saturation trends observed between GEDI RH95% and CHM-H95 on footprints (Figure 6) are still observed between the prediction map and CHM-H95 at the pixel level. The validation results obtained in this study for height estimation combining GEDI and remote sensing images (r 2 from 0.54 to 0.65) are slightly better than the results seen in the literature combining GLAS and Landsat (r 2 ~0.46) [START_REF] Hansen | Mapping Tree Height Distributions in Sub-Saharan Africa Using Landsat 7 and 8 Data[END_REF][START_REF] Wang | Estimation of Forest Canopy Height and Aboveground Biomass from Spaceborne LiDAR and Landsat Imageries in Maryland[END_REF], and similar to [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF] which uses GEDI and Landsat-8 (r 2 ~0.61), although the GLAD-2019 map shows poor results on our validation data (r 2 from 0.01 to 0.21). More investigation on GEDI height metrics and their relationships with reference data (measurements or ALS data) could make it possible to apply a correction of the GEDI data on these temperate forests and to use more accurate values to train the models. Alongside this, the filtering of GEDI samples could be improved (outliers, noise, acquisition dates). Thereby, GEDI data could replace or complete field measurements in machine learning methods. The accuracy of the predictions slightly decreases using GEDI as training data instead of field measurements (relative RMSE 7.3-11.4% with plot measurements against 12.9-16.7% with GEDI on the same validation plots). Nonetheless, the quantity and spatial distribution of GEDI data can compensate for the lower accuracy than that of field measurements in large areas where several forest types are present.

Once the Hdom predictor model is built, the next step is to apply it to the whole satellite image to generate forest height maps. Figure 12 shows Hdom maps over Orléans (top) and StGobain (bottom). The graphs on the left side are extracted from ALS data, the graphs in the centre are produced with GEDI and satellite imagery (methodology presented in this paper), and the graphs on the right side are the GLAD-2019 maps [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF]. The forest mask presented in Section 2.3.1 has been used to select only pixels from forest classes. We see that the forest height map in StGobain shows a more mature and dense forest, with higher trees and not many gaps, whereas the Orleans forest height map shows a more open forest, with higher diversity and height levels. Predictions from models trained on GEDI and satellite images allow the retrieval of general patterns of dominant height in forest stands, although the dynamic range and spatial precision are not as good as those of the map from the CHM. Looking at the GLAD-2019 maps, we see, as in Figure 3 (Section 3.1), that the dynamic range is very low compared to the ALS data and even compared to the maps produced in this paper. This shows that it is possible to obtain better results from GEDI and satellite images when adapting the method to these managed temperate forests and with a larger diversity of satellite sensors. Beyond that, several issues must be taken into account to produce accurate maps over larger territories: standardise the satellite image data (temporal shifts between orbits, different cloud frequencies among regions, etc.) and split the territory between eco-climatic regions, improve the quality of land cover maps (forest localisation and species), and design a nomenclature that allows the best precision of the models by class.

Remote Sens. 2022, 14, x FOR PEER REVIEW 18 of 26 remote sensing images (r 2 from 0.54 to 0.65) are slightly better than the results seen in the literature combining GLAS and Landsat (r 2 ~0.46) [START_REF] Hansen | Mapping Tree Height Distributions in Sub-Saharan Africa Using Landsat 7 and 8 Data[END_REF][START_REF] Wang | Estimation of Forest Canopy Height and Aboveground Biomass from Spaceborne LiDAR and Landsat Imageries in Maryland[END_REF], and similar to [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF] which uses GEDI and Landsat-8 (r 2 ~0.61), although the GLAD-2019 map shows poor results on our validation data (r 2 from 0.01 to 0.21). More investigation on GEDI height metrics and their relationships with reference data (measurements or ALS data) could make it possible to apply a correction of the GEDI data on these temperate forests and to use more accurate values to train the models. Alongside this, the filtering of GEDI samples could be improved (outliers, noise, acquisition dates). Thereby, GEDI data could replace or complete field measurements in machine learning methods. The accuracy of the predictions slightly decreases using GEDI as training data instead of field measurements (relative RMSE 7.3-11.4% with plot measurements against 12.9-16.7% with GEDI on the same validation plots). Nonetheless, the quantity and spatial distribution of GEDI data can compensate for the lower accuracy than that of field measurements in large areas where several forest types are present.

Once the Hdom predictor model is built, the next step is to apply it to the whole satellite image to generate forest height maps. Figure 12 shows Hdom maps over Orléans (top) and StGobain (bottom). The graphs on the left side are extracted from ALS data, the graphs in the centre are produced with GEDI and satellite imagery (methodology presented in this paper), and the graphs on the right side are the GLAD-2019 maps [START_REF] Potapov | Mapping Global Forest Canopy Height through Integration of GEDI and Landsat Data[END_REF]. The forest mask presented in Section 2.3.1 has been used to select only pixels from forest classes. We see that the forest height map in StGobain shows a more mature and dense forest, with higher trees and not many gaps, whereas the Orleans forest height map shows a more open forest, with higher diversity and height levels. Predictions from models trained on GEDI and satellite images allow the retrieval of general patterns of dominant height in forest stands, although the dynamic range and spatial precision are not as good as those of the map from the CHM. Looking at the GLAD-2019 maps, we see, as in Figure 3 (Section 3.1), that the dynamic range is very low compared to the ALS data and even compared to the maps produced in this paper. This shows that it is possible to obtain better results from GEDI and satellite images when adapting the method to these managed temperate forests and with a larger diversity of satellite sensors. Beyond that, several issues must be taken into account to produce accurate maps over larger territories: standardise the satellite image data (temporal shifts between orbits, different cloud frequencies among regions, etc.) and split the territory between eco-climatic regions, improve the quality of land cover maps (forest localisation and species), and design a nomenclature that allows the best precision of the models by class.

(a) The fact that GEDI data are a one-time operation (acquired only over a period of 3 years) is a problem for continuing medium-or long-term forest monitoring. However, the Multi-footprint Observation LiDAR and Imager system (MOLI) is planned by JAXA to provide continuity in the 2022-2027 period.

The amount of data and the large spatial coverage will help to better understand forest parameters estimation using satellite images. Moreover, models learned from a large representativeness of forests should be able to be applied over the following years without new learning processes and thus allow annual map production with our processing chain. Future missions providing global scale open access images, such as NISAR (L-band time series at 10 m spatial resolution) or Sentinel-HR (optical images at very high spatial resolution), will help to provide accurate forest information maps and make sustainable forest monitoring possible with spaceborne images.

Conclusions

This study confirms the interest of machine learning techniques for forest parameters estimation in European temperate forests. These forests are highly variable (heterogeneous stands, uneven age, mixed species, management, environment), making it difficult to predict forest parameters on a large scale by simple laws relating remote sensing observables to a specific forest parameter or physical based models. There is a need for more complex systems such as machine learning methods. Our processing chain was validated on three semi-natural managed broadleaved forests and one pine forest in temperate areas. We demonstrated that the combination of Sentinel-1 backscatters, Sentinel-2 indices and textural metrics at 10 m spatial resolution from both sensors, together with ALOS-2 PALSAR-2 backscatter, was able to provide accurate estimations of forest dominant height with RMSE from 1.9 to 3.1 m (7.3 to 11.6% relative RMSE). Replacing field data with GEDI RH95% data for model learning provided satisfactory results, with RMSE from 3.0 to 4.6 m (12.8 to 16.7% relative RMSE). The study also shows the valuable contribution of GEDI data by using the GEDI RH95% as a reference dataset in the machine learning approach. This contribution is twofold: on one hand, the large amount of reference data using GEDI data allows the stepwise selection of input features in the processing chain to be bypassed, reducing the computational time and complexity. On the other hand, the use of GEDI data as reference data allows more representativeness and forest height mapping to be provided in areas where little or no in situ data are available. The fact that GEDI data are a one-time operation (acquired only over a period of 3 years) is a problem for continuing medium-or long-term forest monitoring. However, the Multi-footprint Observation LiDAR and Imager system (MOLI) is planned by JAXA to provide continuity in the 2022-2027 period.

The amount of data and the large spatial coverage will help to better understand forest parameters estimation using satellite images. Moreover, models learned from a large representativeness of forests should be able to be applied over the following years without new learning processes and thus allow annual map production with our processing chain. Future missions providing global scale open access images, such as NISAR (L-band time series at 10 m spatial resolution) or Sentinel-HR (optical images at very high spatial resolution), will help to provide accurate forest information maps and make sustainable forest monitoring possible with spaceborne images.

Conclusions

This study confirms the interest of machine learning techniques for forest parameters estimation in European temperate forests. These forests are highly variable (heterogeneous stands, uneven age, mixed species, management, environment), making it difficult to predict forest parameters on a large scale by simple laws relating remote sensing observables to a specific forest parameter or physical based models. There is a need for more complex systems such as machine learning methods. Our processing chain was validated on three semi-natural managed broadleaved forests and one pine forest in temperate areas. We demonstrated that the combination of Sentinel-1 backscatters, Sentinel-2 indices and textural metrics at 10 m spatial resolution from both sensors, together with ALOS-2 PALSAR-2 backscatter, was able to provide accurate estimations of forest dominant height with RMSE from 1.9 to 3.1 m (7.3 to 11.6% relative RMSE). Replacing field data with GEDI RH95% data for model learning provided satisfactory results, with RMSE from 3.0 to 4.6 m (12.8 to 16.7% relative RMSE). The study also shows the valuable contribution of GEDI data by using the GEDI RH95% as a reference dataset in the machine learning approach. This contribution is twofold: on one hand, the large amount of reference data using GEDI data allows the stepwise selection of input features in the processing chain to be bypassed, reducing the computational time and complexity. On the other hand, the use of GEDI data as reference data allows more representativeness and forest height mapping to be provided in areas where little or no in situ data are available.

Building on the results shown in this study, the next step is to apply our processing chain to larger areas such as regions or at country scale and to produce new and accurate

Figure 1 .

 1 Figure 1. Study sites in the north of France. Measurement plots are shown on both sites, background images are Sentinel-2 2016 summer dates in false colours (active vegetation appears in red).

Figure 1 .

 1 Figure 1. Study sites in the north of France. Measurement plots are shown on both sites, background images are Sentinel-2 2016 summer dates in false colours (active vegetation appears in red).

Figure 2 .

 2 Figure 2. Overall processing chain with the example of using a forward stepwise selection on 36 features.

Figure 2 .

 2 Figure 2. Overall processing chain with the example of using a forward stepwise selection on 36 features.

Figure 3 .Figure 3 .

 33 Figure 3. Cont.

Figure 3 .

 3 Figure 3. GLAD-2019 height map validation using plot measurement. (a) Orléans, oak; (b) Orléans, Scots pine; (c) StGobain, oak; (d) StGobain, beech. Statistical signicicance: * p ≤ 0.05.

Figure 4 .

 4 Figure 4. Coefficient of determination (r 2 ) of dominant height estimation (ONF field measurements for model training and validation with a LOO cross-validation), evolution with the number of features selected with a forward stepwise selection.

Figure 4 .

 4 Figure 4. Coefficient of determination (r 2 ) of dominant height estimation (ONF field measurements for model training and validation with a LOO cross-validation), evolution with the number of features selected with a forward stepwise selection.

Figure 4 .Figure 5 .

 45 Figure 4. Coefficient of determination (r 2 ) of dominant height estimation (ONF field measurements for model training and validation with a LOO cross-validation), evolution with the number of features selected with a forward stepwise selection.

Figure 5 .

 5 Figure 5. Validation scatterplot of dominant height (Hdom) estimations with the processing chain (forward selection), using ONF field measurements and satellite images data for model training and a LOO cross-validation. (a) Orléans, oak; (b) Orléans, Scots pine; (c) StGobain, oak; (d) StGobain, beech. Statistical signicicance: *** p ≤ 0.001.

Figure 6 .Figure 6 .

 66 Figure 6. Cont.

Figure 6 .

 6 Figure 6. Comparison of the height 95th percentile (H95, proxy for dominant height) from ALS data and GEDI relative height metric at 95% (RH95%). (a) Orléans, oak; (b) Orléans, Scots pine; (c) StGobain, oak; (d) StGobain, beech. Statistical signicicance: *** p ≤ 0.001.

Figure 7 .

 7 Figure 7. Coefficient of determination (r 2 ) of dominant height estimation (using GEDI RH95% for model training and validation using a K-fold cross-validation), evolution with the number of features selected with a backward stepwise selection (starts with all 36 features and ends with 1 remaining feature).

Figure 8 .

 8 Figure 8. Oak samples (Orléans). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (479 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (75 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 7 .

 7 Figure 7. Coefficient of determination (r 2 ) of dominant height estimation (using GEDI RH95% for model training and validation using a K-fold cross-validation), evolution with the number of features selected with a backward stepwise selection (starts with all 36 features and ends with 1 remaining feature).

Figure 7 .

 7 Figure 7. Coefficient of determination (r 2 of dominant height estimation (using GEDI RH95% for model training and validation using a K-fold cross-validation), evolution with the number of features selected with a backward stepwise selection (starts with all 36 features and ends with 1 remaining feature).

Figure 8 .

 8 Figure 8. Oak samples (Orléans). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (479 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (75 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 8 .

 8 Figure 8. Oak samples (Orléans). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (479 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (75 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 9 .

 9 Figure 9. Scots pine samples (Orléans). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (958 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (28 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 10 .

 10 Figure 10. Oak samples (StGobain). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (808 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (21 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 9 . 26 Figure 9 .

 9269 Figure 9. Scots pine samples (Orléans). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (958 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (28 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 10 .

 10 Figure 10. Oak samples (StGobain). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (808 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (21 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 10 .

 10 Figure 10. Oak samples (StGobain). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (808 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (21 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 11 .

 11 Figure 11. Beech samples (StGobain). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (723 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (25 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 11 .

 11 Figure 11. Beech samples (StGobain). Validation scatterplots of dominant height estimation using GEDI RH95% metric for the training model (723 samples). (a) GEDI's K-fold cross-validation; (b) application of GEDI model on ONF measurement plots (25 samples). Statistical signicicance: *** p ≤ 0.001.

Figure 12 .Figure 12 .

 1212 Figure 12. Cont.

Figure 12 .

 12 Figure 12. Local examples of Hdom maps of Orleans ((a), top) and StGobain ((b), bottom). Left maps: CHM-H95 computed at 10 m resolution; centre: prediction maps using the processing chain and the GEDI RH95% as reference for model learning; right: GLAD-2019 maps.

Figure A1 .Figure A2 .

 A1A2 Figure A1. Validation scatterplot forest parameters estimation with ONF field measurements and satellite remote sensing data. (a) dominant height (Hdom); (b) mean diameter at breast height (DBH). Statistical signicicance: *** p ≤ 0.001.

Figure A1 .Figure A1 .Figure A2 .

 A1A1A2 Figure A1. Validation scatterplot forest parameters estimation with ONF field measurements and satellite remote sensing data. (a) basal area (BA); (b) tree density (Density). Statistical signicicance: *** p ≤ 0.001.

Figure A2 .

 A2 Figure A2. Validation scatterplot forest parameters estimation with ONF field measurements and satellite remote sensing data. (a) dominant height (Hdom); (b) mean diameter at breast height (DBH). Statistical signicicance: *** p ≤ 0.001.

Table 1 .

 1 Mean and standard deviation of field inventory for basal area (BA), mean DBH (DBH), dominant height (Hdom) and tree density (Density).

		Orléans		StGobain
		Oak	Scots Pine	Oak	Beech
	Plot number	75	28	21	25
	Hdom (m

) 22.1 ± 3.6 20.8 ± 4.0 26.2 ± 3.7 27.3 ± 6.2 BA (m 2 /ha)

		18.0 ± 5.8	20.8 ± 8.3	21.0 ± 8.8	19.1 ± 9.2
	DBH (cm)	38.8 ± 12.0	35.8 ± 8.8	37.7 ± 7.2	36.8 ± 8.8
	Density (tree/ha)	164 ± 85	237 ± 133	155 ± 68	153 ± 66
	2.2. Remote Sensing Data
	2.2.1. Satellite Imagery
	This study uses three different satellite sensors: Sentinel-1, Sentinel-2 and ALOS-2
	PALSAR-2 (all 3 open access).

Table 2 .

 2 Validation statistics of the estimation of dominant height (Hdom) of the four classes. Training and test/validation data are indicated in the left column, RSdata: satellite image features from Sentinel-1, Sentinel-2 and ALOS-2 PALSAR-2.

		Site		Orléans	StGobain
		Class (Plots)	Oak (75)	Scots Pine (28)	Oak (21)	Beech (25)
	GLAD-2019 Test ONF	r 2 RMSE MAE	0.01 6.9 (31.0%) 4.9 (22.2%)	0.21 8.6 (41.6%) 7.1 (34.2%)	0.21 3.9 (14.8%) 3.3 (12.8%)	0.1 6.6 (24.3%) 5.9 (21.5%)
	Train ONF/RSdata Test ONF	r 2 RMSE MAE	0.48 2.6 (11.6%) 2.0 (9.1%)	0.78 1.9 (9.2%) 1.3 (6.2%)	0.74 1.9 (7.3%) 1.4 (5.4%)	0.78 3.1 (11.4%) 2.2 (7.9%)
	Train GEDI/RSdata Test ONF	r 2 RMSE MAE	0.23 3.3 (14.7%) 2.7 (12.4%)	0.55 3.0 (14.3%) 2.4 (11.3%)	0.64 3.3 (12.8%) 2.9 (11.1%)	0.54 4.6 (16.7%) 3.6 (13.3%)
		Class (footprints)	Oak (479)	Scots pine (958)	Oak (808)	Beech (723)
	Train GEDI/RSdata	r 2	0.54	0.64	0.65	0.59
	Test GEDI	RMSE	3.4 (16.7%)	2.7 (15.5%)	2.9 (12.9%)	3.4 (13.9%)
		MAE	2.6 (12.9%)	2.2 (12.2%)	2.4 (10.3%)	2.7 (11.1%)
	Maps comparison:	Class (pixels)	Oak (149,165) Scots pine (168,031) Oak (113,916)	Beech (107,213)
	Train GEDI/RSdata	r 2	0.49	0.57	0.56	0.47
	Test CHM-H95	RMSE	4.5 (22.9%)	3.5 (19.2%)	4.0 (16.9%)	4.8 (19.0%)

Table 2 .

 2 Validation statistics of the estimation of dominant height (Hdom) of the four classes. Training and test/validation data are indicated in the left column, RSdata: satellite image features from Sentinel-1, Sentinel-2 and ALOS-2 PALSAR-2.

		Site		Orléans	StGobain
		Class (Plots)	Oak (75)	Scots Pine (28)	Oak (21)	Beech (25)
		r 2	0.01	0.21	0.21	0.1
	GLAD-2019 Test ONF	RMSE	6.9 (31.0%)	8.6 (41.6%)	3.9 (14.8%)	6.6 (24.3%)
		MAE	4.9 (22.2%)	7.1 (34.2%)	3.3 (12.8%)	5.9 (21.5%)
		r 2	0.48	0.78	0.74	0.78
	Train ONF/RSdata Test ONF	RMSE	2.6 (11.6%)	1.9 (9.2%)	1.9 (7.3%)	3.1 (11.4%)
		MAE	2.0 (9.1%)	1.3 (6.2%)	1.4 (5.4%)	2.2 (7.9%)
		r 2	0.23	0.55	0.64	0.54
	Train GEDI/RSdata Test ONF	RMSE	3.3 (14.7%)	3.0 (14.3%)	3.3 (12.8%)	4.6 (16.7%)
		MAE	2.7 (12.4%)	2.4 (11.3%)	2.9 (11.1%)	3.6 (13.3%)
		Class (footprints)	Oak (479)	Scots pine (958)	Oak (808)	Beech (723)
	Train GEDI/RSdata	r 2	0.54	0.64	0.65	0.59
	Test GEDI	RMSE	3.4 (16.7%)	2.7 (15.5%)	2.9 (12.9%)	3.4 (13.9%)
		MAE	2.6 (12.9%)	2.2 (12.2%)	2.4 (10.3%)	2.7 (11.1%)
		Class (pixels)	Oak (149,165)	Scots pine (168,031)	Oak (113,916)	Beech (107,213)
	Maps comparison: Train GEDI/RSdata	r 2	0.49	0.57	0.56	0.47
	Test CHM-H95	RMSE	4.5 (22.9%)	3.5 (19.2%)	4.0 (16.9%)	4.8 (19.0%)
		MAE	3.6 (18.3%)	2.8 (15.1%)	3.2 (13.4%)	3.8 (14.9%)

Table A1 .

 A1 Cross-validation statistics of the estimation of basal area (BA), tree density (Density), dominant height (Hdom) and mean DBH (DBH).

			Orléans	StGobain
			Oak	Scots Pine	Oak	Beech
		r 2	0.37	0.62	0.64	0.77
	BA (m 2 /ha)	RMSE	4.62 (25.7%)	5.18 (24.9%)	5.38 (25.7%)	4.57 (23.9%)
		MAE	3.94 (22.0%)	4.17 (20.1%)	4.29 (20.4%)	3.14 (16.5%)
		r 2	0.61	0.78	0.64	0.56
	Density (tree/ha)	RMSE	53.4 (32.6%)	66.2 (28.0%)	42.3 (27.3%)	44.3 (28.9%)
		MAE	36.5 (22.3%)	47.3 (20.0%)	35.5 (22.9%)	34.2 (22.3%)
		r 2	0.48	0.78	0.74	0.78
	Hdom (m)	RMSE	2.56 (11.6%)	1.92 (9.2%)	1.90 (7.3%)	3.10 (11.4%)
		MAE	2.01 (9.1%)	1.29 (6.2%)	1.40 (5.4%)	2.16 (7.9%)
		r 2	0.56	0.88	0.89	0.78
	DBH (cm)	RMSE	7.98 (20.5%)	3.34 (9.4%)	2.43 (6.4%)	4.17 (11.3%)
		MAE	6.56 (16.9%)	2.56 (7.1%)	2.06 (5.5%)	3.41 (9.3%)

Remote Sens. 2022, 14, x FOR PEER REVIEW

forest dominant height maps for foresters, public bodies and climate modellers working on areas with mixed heterogeneous forests.

Acknowledgments:

The authors would like to thank the Office Nationale des Forêts (ONF), Jérôme Bock and Anne Jolly for providing the field measurement data and the ALS data used in this study.

Funding: This research was funded by the Centre National d'Etudes Spaciales (CNES, TOSCA 2021 project) and the Agence de l'Environnement et de la Maitrise de l'Energie (ADEME, GRAINE 2019 project).

Data Availability Statement: GEDI data are available from the NASA/USGS Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/, accessed on 16 March 2022). Sentinel-1 data are available from the PEPS platform (https://peps.cnes.fr, accessed on 16 March 2022). Sentinel-2 data are available from the Theia platform (https://theia.cnes.fr, accessed on 16 March 2022), and ALOS PALSAR mosaics are available from (https://www.eorc.jaxa.jp/ALOS/en/dataset/fnf_e.htm, accessed on 16 March 2022). Field measurements and ALS data are the property of the Office Nationale des Forêts (ONF).

 A1.

The overall results are good. Merged statistics show RMSE of 4.8 m 2 /ha (25%) for basal area, 53 trees/ha (30%) for tree density, 2.5 m (11%) for dominant height and 6.2 cm (16%) for mean DBH. Orléans oak parameters are slightly harder to estimate. The canopy is more open in these oaks stands; moreover, it is partly coppice forest and the threshold of 17.5 cm DBH for measurement probably introduces a bias. BA and Density results are globally more scattered, which is because these parameters are computed on surfaces and thus are more impacted by the variability of tree positioning, while Hdom and DBH are averaged variables and therefore more stable.