# Estimation and mapping of forest structure parameters from open access satellite images: development of a generic method with a study case on coniferous plantation

David Morin, Milena Planells, Dominique Guyon, Ludovic Villard, Stéphane Mermoz, Alexandre Bouvet, Hervé Thevenon, Jean-François Dejoux, Thuy Le Toan, Gérard Dedieu

► To cite this version:

David Morin, Milena Planells, Dominique Guyon, Ludovic Villard, Stéphane Mermoz, et al. Estimation and mapping of forest structure parameters from open access satellite images: development of a generic method with a study case on coniferous plantation. Remote Sensing, 2019, 11 (11), pp.1-25. 10.3390/rs1111275. hal-02629195

HAL Id: hal-02629195

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

Submitted on 27 May 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Estimation and Mapping of Forest Structure Parameters from Open Access Satellite Images: Development of a Generic Method with a Study Case on Coniferous Plantation

David Morin<sup>1,*</sup>, Milena Planells<sup>1</sup>, Dominique Guyon<sup>2</sup>, Ludovic Villard<sup>1</sup>,
Stéphane Mermoz<sup>1,3 ID</sup>, Alexandre Bouvet<sup>1 ID</sup>, Hervé Thevenon<sup>1 ID</sup>, Jean-François Dejoux<sup>1</sup>,
Thuy Le Toan<sup>1</sup> and Gérard Dedieu<sup>1 ID</sup>

Received: 24 April 2019; Accepted: 25 May 2019; Published: 29 May 2019

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

check for
updates

**Abstract:** Temperate forests are under climatic and economic pressures. Public bodies, NGOs and the wood industry are looking for accurate, current and affordable data driven solutions to intensify wood production while maintaining or improving long term sustainability of the production, biodiversity, and carbon sequestration. Free tools and open access data have already been exploited to produce accurate quantitative forest parameters maps suitable for policy and operational purposes. These efforts have relied on different data sources, tools, and methods that are tailored for specific forest types and climatic conditions. We hypothesized we could build on these efforts in order to produce a generic method suitable to perform as well or better in a larger range of settings. In this study we focus on building a generic approach to create forest parameters maps and confirm its performance on a test site: a maritime pine (*Pinus pinaster*) forest located in south west of France. We investigated and assessed options related with the integration of multiple data sources (SAR L- and C-band, optical indexes and spatial texture indexes from Sentinel-1, Sentinel-2 and ALOS-PALSAR-2), feature extraction, feature selection and machine learning techniques. On our test case, we found that the combination of multiple open access data sources has synergistic benefits on the forest parameters estimates. The sensibility analysis shows that all the data participate to the improvements, that reach up to 13.7% when compared to single source estimates. Accuracy of the estimates is as follows: aboveground biomass (AGB) 28% relative RMSE, basal area (BA) 27%, diameter at breast height (DBH) 20%, age 17%, tree density 24%, and height 13%. Forward feature selection and SVR provided the best estimates. Future work will focus on validating this generic approach in different settings. It may prove beneficial to package the method, the tools, and the integration of open access data in order to make spatially accurate and regularly updated forest structure parameters maps effortlessly available to national bodies and forest organizations.

**Keywords:** forest structure parameters; aboveground biomass; map; forest plantation; Sentinel; ALOS-PALSAR; support vector regression; data combinations; management tools

<footnote n="1">¹ Centre Détudes Spatiales De La Biosphère (CESBIO), Université de Toulouse, CNES/CNRS/IRD/UT3, UMR 5126, 31401 Toulouse, France; milena.planells@cnes.fr (M.P.); villardl@cesbio.cnes.fr (L.V.); mermoz@globeo.net (S.M.); bouveta@cesbio.cnes.fr (A.B.); thevenonh@cesbio.cnes.fr (H.T.); dejouxjf@cesbio.cnes.fr (J.-F.D.); letoant@cesbio.cnes.fr (T.L.T.); dedieug@cesbio.cnes.fr (G.D.)</footnote>

<footnote n="2">² Intéraction Sol Plante Atmosphère (ISPA), INRA, UMR 1391, 33882 Villenave d'Ornon, France; dominique.guyon@inra.fr</footnote>

<footnote n="3">³ Global Earth Observation (GlobEO), 31400 Toulouse, France</footnote>

<footnote>* Correspondence: david.morin@cesbio.cnes.fr</footnote>
# 1. Introduction

Forests provide several ecosystem services such as carbon storage, climate regulation, biodiversity and wood production. Over the past 30 years, the extent of the world's forests has declined as human population continues to grow and its demand for food and land increases [1]. Forest changes and trends are not the same worldwide. Losses of forest area are mostly located in tropical zones, while the area of temperate forests is stable or sometimes even increasing in private and unmanaged forests [1]. The primary regions featuring temperate forests are Europe, North America, northeast Asia, Patagonia and New Zealand. With 25% of the total forest's areas, temperate forests are a net carbon sink that hold 16% of the global plant biomass [2].

While threats to tropical and boreal forests are widely acknowledged, temperate forests are also under climatic and economical pressures. Although temperate forests do not suffer from large-scale deforestation, they experience hotter droughts with higher frequency, longer heat waves, severe wildfires and pathogen outbreaks. The frequency of these events may affect their resilience in the long term. Beside these climatic pressures, European temperate and boreal forests have become a target for renewable energy (2020–2030 European energy and environmental policies) and the primary resource for new commercial endeavors. As a consequence, local, national and European bodies are looking forward to intensify already managed forest areas and to bring unmanaged forest areas under intensive management. The attention paid to sustainable forest management has never been higher for maintaining biodiversity, increasing carbon uptake and forest resilience.

Temperate forests are managed both globally and locally. Private forest owners are required to manage their forests according to legislation focused on sustainability and risks prevention. Policy makers and foresters have a growing need for regular estimation and mapping of forest resources indicators with high spatial resolution. Tree height, diameter at breast height (DBH) and basal area (BA) are essential parameters for foresters to assess the current or prospective commercial value of forest stands and to plan thinning and cutting. Tree density and age complement these forest structure parameters and allow to build allometric relationships in order to estimate aboveground biomass (AGB). AGB is indeed one of the essential climate variables (ECVs) identified by the Global Climate Observing System (GCOS). At a national scale, most countries in temperate region conduct forest inventories that produce normalized and actualized data on a yearly basis. The resulting statistics are accurate over large administrative areas and provide useful information such as changes in the national forest estate and the main tree species distributions. However, this approach is neither designed for mapping spatially continuous large area, nor to carry out precise estimates at a local scale that could be of interest to foresters. In order to satisfy their operational needs, foresters perform forest stand measurements that usually lack normalization (due to estimation methods, spatial sampling and surveyor bias) and are less often updated than the measurements performed by national bodies.

To overcome the limitations and related uncertainties of their surveys, national bodies have researched two main areas of improvement: (1) additional airborne/spaceborne remote sensing instruments (e.g., lidar, tomography or photogrammetry) with the view to generate better data [3–8], and (2) smarter use of existing data through statistical or machine learning techniques that mutually enrich unrelated data sources (multi-source approaches [9–12]). Our study fits into this second field, considering the use of freely available satellite images as a promising way to build a cost-effective method for mapping forest structure parameters [13–19].

The current availability of open access synthetic aperture radar (SAR) and optical data, as well as open source image processing software, makes the use of spaceborne remote sensing data more easily accessible and suitable for forest mapping and monitoring over large areas with high spatial resolution. Effectively this state of affairs provides an opportunity to develop better products for both communities. SAR sensors [20–33] can operate day and night and penetrate clouds. The transmitted microwaves penetrate into the forest canopies in more or less degree depending on the wavelength. The reflected energy is proportional to the volume of scatterers and the geometric and electromagnetic properties of the landscape. Shorter wavelengths X- and C-band are sensitive to smaller canopy elements (leaves
and small branches), while longer wavelengths L- and P-band are sensitive to large branches and trunks and have higher sensitivity with forest biomass [34]. C-band and L-band images are currently available, free and worldwide thanks to the Sentinel-1 and ALOS-PALSAR-2 satellites. The launch of the P-band radar of BIOMASS satellite in 2022 will further improve forests biomass estimations, in particular for dense tropical forests [35].

In the optical domain, there are direct relationships between the reflectance (visible and near to medium infrared) and the green leaf area of the whole canopy (tree and understory vegetation) and its horizontal and vertical distribution. However, the signal often saturates at high leaf area values. Optical data is also used to monitor the phenology and stress conditions of vegetation. Although optical data are not directly linked to forest parameters such as height or DBH, some studies show possible relationships between reflectance data and forest structure parameters [36–38]. In addition to the intensity of the reflectance signal, it is possible to use the spatial indicators and the temporal information of satellite image time series. Spatial indicators derived from very high spatial resolution (VHRS, <1 m) optical images, using Fourier-based or Grey Level Spatial co-occurrence Matrix (GLSM) textural metrics, have been used successfully to estimate forest structure parameters and AGB [39–42]. This approach allows to retrieve vertical and horizontal structure information of forests based on the distribution of the reflectance along adjacent pixels. VHRS allows to observe more precise and more specific patterns depending on tree crown sizes and tree spacing. Currently, Landsat and Sentinel-2 time-series are freely available worldwide with 10–30 m spatial resolution. A few studies successfully used images at this spatial resolution for retrieving forest structure parameters with textural indexes [43–46]. Temporal dimension is useful for tree species classification [47], to observe forest phenology and to monitor forest health status. Temporal dimension does also improve forest parameters estimations [38,44,48].

The value of remote sensing data to estimate forest parameters such as DBH, AGB, tree height and density has been demonstrated. Recent studies suggest that different combinations of satellite data, i.e., multi-frequency SAR data, optical and textural data improve forest parameters estimation [30–32,49–52]. Methodology and data vary greatly across these studies. Parametric semi-empirical and physically-based models provide accurate estimations and are easy to interpret. Non-parametric approaches are more flexible in combining different types of data and applying the method to other types of forests, making them potentially more suitable for building a generic method. In order to be useful and relevant, forest parameters maps should have the following characteristics: high spatial resolution (~10 to 50 m), coverage of large territories, ideally worldwide so it can be produced on any territory, regular updates (1~2 years maximum delay), and be cost effective.

In this context, our study aims at building an automatic and generic method using freely available SAR and optical data at high resolution, with open source software, to produce quantitative maps of forest structure parameters and AGB with regular updates. We use Sentinel-1, Sentinel-2 and ALOS-PALSAR mosaics acquired worldwide with high spatial resolution (10–25 m). We extract commonly accepted features for forest parameters estimation, then we analyze dimensionality reduction methods and algorithm performance for the estimation of forest parameters. The objective is to develop a flexible method applicable to a wide range of forests and regions, that can integrate several types of remote sensing features. This paper presents the method we developed and its validation on maritime pine plantations within the largest European plantation forest known as the *Landes of Gascogne*. In Section 2, we introduce the study site, the satellite sensors and image processing steps, and the selection of a set of remote sensing features used in our study. Section 3 describes the estimation methods and the machine learning algorithms we assessed. Section 4 presents the results obtained using all feature types (SAR L- and C-band, optical and texture metrics from both SAR and optical images) and the contribution of each one of them for the estimation of AGB, BA, DBH, age, density and height.
## 2. Materials

### 2.1. Study Site

The site is part of the largest European plantation forest (nearly 1 million hectares of maritime pine, *Pinus pinaster*) known as *Landes de Gascogne* forest, here referred as Landes forest (Figure 1). The site's climate is oceanic, the land is rather flat, the soil consists mainly of sandy podzols. The stands are even-aged and intensively managed. Thinning generally occurs every 5 years when the plantation is 10 to 15 years old. Stands are clear cut when they are approximately 50 years old, and replanted within 2 to 3 years. These management practices give homogeneity in the within-stand structure.

<figure><img src="image_4.png" /><figcaption>Figure 1. Study site location in Southwest France. The measured stands are in red.</figcaption></figure>

### 2.2. In Situ Data

We use the forest database produced in the framework of the GLORI project [53,54], from an in situ measurement campaign that was carried out in 2016 before the beginning of the tree growth. 83 stands of maritime pine were selected in an effort to cover all the range of forest structure parameters specific to maritime pine in the Landes forest. As in [41], one sample plot is selected in a rather homogeneous area within each stand. The plot size varies from 0.05 to 0.2 hectares. The center of each sample plot is geolocated using global positioning system (GPS). All trees are counted and the tree density is calculated. DBH and height are measured for the trees closest to the center (10 and 5 trees respectively). BA is estimated from the individual trunk cross-section derived from DBH and the tree density. AGB is estimated from the individual biomass calculated with a tree level allometric equation using DBH as input parameter [55] and the tree density:

$$AGBi = a \cdot DBHi^b \qquad (1)$$

where a = 5.013 tons/m and b = 2.48. AGBi is the individual aboveground biomass in ton, and DBHi is the individual DBH in meter. Stand age is provided by forest managers or estimated by counting the apparent number of annual growth cycles in height, by locating branch whorls.

The minimum, maximum, and mean values of the measured forest structure parameters in these 83 plots are presented in Table 1.
<table><caption>**Table 1.** Range of the forest structure parameters for the 83 sample plots selected in our study.</caption><thead><tr><th></th><th>Minimum</th><th>Maximum</th><th>Mean</th></tr></thead><tbody><tr><td>AGB (tons/ha)</td><td>1</td><td>136</td><td>70</td></tr><tr><td>BA (m²/ha)</td><td>0.8</td><td>44.1</td><td>21.2</td></tr><tr><td>Mean tree height (m)</td><td>2</td><td>24.7</td><td>14.4</td></tr><tr><td>Mean DBH (m)</td><td>0.03</td><td>0.57</td><td>0.22</td></tr><tr><td>Density (tree/ha)</td><td>87</td><td>2622</td><td>835</td></tr><tr><td>Age (year)</td><td>3</td><td>72</td><td>22.6</td></tr></tbody></table>

The 83 stands have an average size of 20.9 ha. Visual analysis of aerial photographs revealed local heterogeneities within the stands. We digitalized smaller sub-stands centered on the measured samples in order to have homogeneous sub-stands. The new 83 sub-stands digitalized have an average size of 2.7 ha, the minimum is 1 ha. They are used to extract values from remote sensing images. Figure 2 and Table 2 illustrate the plot-level relationships of forest parameters. They depend on allometry properties specific to maritime pine, soil fertility and history of forest management specific to each stand. The correlation between AGB and mean DBH is strong ($r^2$ > 0.8) for diameters under 0.15 m. There is no AGB - DBH correlation for wider diameters ($r^2$ < 0.1). The DBH is correlated with tree height (Figure 2a) until a diameter of ~0.4 m ($r^2$ = 0.94), above 0.4m DBH the height is stable but the DBH continues to grow. Finally, the DBH and density (Figure 2b) show a linear relationship for diameter in the range 0.2 to 0.4 m ($r^2$ > 0.7). As the DBH falls below 0.2 m the density can vary from 750 to 1750 trees/ha. Thus, the variety of silvicultural practices (such as density of plantation, age at first thinning, intensity and frequency of thinning) affects the allometry relationships between different forest parameters and must be taken into account to interpret the results.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Plot-level relationships of forest structure parameters: (a) mean DBH vs. mean height, (b) mean DBH vs. Density.</figcaption></figure>

<table><caption>**Table 2.** Coefficient of determination between the forest structure parameters (plot-level).</caption><thead><tr><th></th><th colspan="6">r<sup>2</sup></th></tr><tr><th></th><th>AGB</th><th>BA</th><th>Height</th><th>DBH</th><th>Density</th><th>Age</th></tr></thead><tbody><tr><td>AGB</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>BA</td><td>0.90</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Mean tree height</td><td>0.71</td><td>0.46</td><td>1</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Mean DBH</td><td>0.49</td><td>0.24</td><td>0.83</td><td>1</td><td>-</td><td>-</td></tr><tr><td>Density</td><td>0.19</td><td>0.04</td><td>0.55</td><td>0.64</td><td>1</td><td>-</td></tr><tr><td>Age</td><td>0.40</td><td>0.17</td><td>0.74</td><td>0.90</td><td>0.53</td><td>1</td></tr></tbody></table>
## 2.3. Remote Sensing Data

This study uses four different satellite sensors: Sentinel-1, Sentinel-2, ALOS-PALSAR-2 mosaics (all three open access), and SPOT (freely available over France for French public bodies).

### 2.3.1. Optical Data

Sentinel-2 images: Sentinel-2 data consist in optical time series acquired in the visible (blue, green, red), near infrared (nir) and short waves infrared (swir) with a spatial resolution from 10 to 20 m. The raw images are provided by the European Space Agency (ESA). The Theia Land Data Center (website: theia.cnes.fr) produces the level 2A images (top of canopy reflectance: TOC) with ortho-rectification, atmospheric corrections and cloud mask (MAJA processing, [56,57]). Eighteen dates are available in 2016, and seven dates have all the reference plots cloud free in May, August, September, November and December.

Spot-6 images: Spot-6 data consists of annual mosaics with a spatial resolution from 1.5 m (panchromatic image) to 6 m (multi-spectral image in the visible and near infrared). Annual cloud free mosaics are provided by Airbus and available from the Geosud platform (website: ids.equipex-geosud.fr) with ortho-rectification and top of atmosphere (TOA) reflectance. The acquisition date for our test site is 2016-03-22. Panchromatic (PAN) and multi-spectral bands are available.

### 2.3.2. SAR Data

ALOS-PALSAR: L-band SAR annual mosaics at 25 m spatial resolution are produced by the Japan Aerospace Exploration Agency (JAXA) through large scale mosaicking [58] that includes ortho-rectification, slope correction and radiometric calibration between neighboring strips. Seven annual mosaics are sourced from JAXA for years 2007 through 2010, 2015 and 2016. We converted digital numbers (DN) to gamma naught values (calibration following [58]). Then we applied a multi-image filtering [59,60] using 12 channels (six dates at HV and HH polarizations) with a 7 × 7 window size, leading to an equivalent number of looks (ENL) of 156. In the following, we use processed HV and HH channels from 2016.

Sentinel-1: C-band SAR time series at 10 m spatial resolution are provided by ESA and downloaded from the Sentinel Product Exploitation Platform (CNES, website: peps.cnes.fr). We used the Level-1 Ground Range Detected (GRD) products. Three orbits cover our test site, one ascendant (ASC) with ~37° incidence angle (IA) at the center of the test site, and two descendant (DES) with ~33° and ~43° IA. As acquisition occurs every 12 days, there are approximately 30 dates available in 2016 for each orbit. Sentinel-1 time series are calibrated, ortho-rectified and filtered using the same multi-image filter as for ALOS-PALSAR data. Each orbit is filtered separately using 60 channels (30 dates in 2016 at VH and VV polarization) and a 5 × 5 window size, leading to an ENL of 87.

### 2.3.3. Textural Metrics

Textural metrics have been generated from both optical and radar using a Grey Level Spatial co-occurence Matrix (GLSM) [61] with the Orfeo Toolbox library (OTB, open source [62]). Eight features were extracted: energy, entropy, correlation, homogeneity, inertia, cluster shade, cluster prominence, and Haralick correlation. We performed a sensitivity analysis of Haralick texture extraction. The offset, orientation, window size, and grey levels parameters were tested. With offset set to 1 pixel, textural indexes values were not significantly different for orientations 0°, 45°, 90° and 135° even for Spot's 1.5 m spatial resolution. Consequently we set the offset to 1 pixel orientation to 45°. We tested various window sizes and number of grey levels. We found that too small values of window size and number of grey levels produce binary outputs while a large window size is detrimental to the accuracy of dendrometry and land cover outputs. The best compromise was found with a window size of 7 pixels and 30 grey levels for the 10 m resolution images (Sentinel-1 and Sentinel-2 indexes),
9 pixels and 40 grey levels for the 6 m resolution images (Spot-6 multi-spectral indexes), 35 pixels and 100 grey levels for the 1.5 m resolution image (Spot-6 PAN indexes).

### 2.3.4. Selection of Remote Sensing Features

Almost 300 remote sensing features were extracted from the four satellite sensors (not including textural index parameters). When using statistical methods, the number of features must be less than the number of samples. In this case study, we use 83 samples, therefore it is necessary to reduce the number of features. This process is called feature selection. In addition to this requirement we wish features to be both generic and adaptable to different types of forest. Table 3 presents the 36 selected features based on regression tests and correlation analysis. We categorize them into seven feature types (Features types column) according to sensors and the type of processing (raw data, spectral or textural indexes).

<table><caption>Table 3. List of the 36 selected features as input for the regression models.</caption><thead><tr><th>Sensors</th><th>Feature Types</th><th>Features</th><th>Abbreviation</th></tr></thead><tbody><tr><td rowspan="3">ALOS-PALSAR-2</td><td rowspan="3">L-band backscatter (L-band)</td><td>L-band HV</td><td>L-HV</td></tr><tr><td>L-band HH</td><td>L-HH</td></tr><tr><td>L-band HV/HH</td><td>L-Ratio</td></tr><tr><td rowspan="12">Sentinel-1</td><td rowspan="4">C-band backscatter (C-band)<br>DES-33°IA orbit, annual mean values</td><td>C-band VH</td><td>C-VH</td></tr><tr><td>C-band VV</td><td>C-VV</td></tr><tr><td>C-band VH/VV</td><td>C-Ratio</td></tr><tr><td>VH-homogeneity</td><td>C-VH-hom</td></tr><tr><td rowspan="6">C-band textural indexes (C-TI)<br>DES-33°IA orbit, annual mean values</td><td>VH-cluster shade</td><td>C-VH-clu</td></tr><tr><td>VV-homogeneity</td><td>C-VV-hom</td></tr><tr><td>VV-cluster shade</td><td>C-VV-clu</td></tr><tr><td>VH-correlation</td><td>C-VH-cor</td></tr><tr><td>VV-correlation</td><td>C-VV-cor</td></tr><tr><td>VH-Haralick correlation</td><td>C-VH-Hcor</td></tr><tr><td>VV-Haralick correlation</td><td>C-VV-Hcor</td></tr><tr><td rowspan="12">Sentinel-2</td><td rowspan="6">Spectral indexes (S2-SI)</td><td>NDVI July</td><td>S2-NDVI-sum</td></tr><tr><td>NDVI December</td><td>S2-NDVI-win</td></tr><tr><td>BI July</td><td>S2-BI-sum</td></tr><tr><td>BI December</td><td>S2-BI-win</td></tr><tr><td>NDWI July</td><td>S2-NDWI-sum</td></tr><tr><td>NDWI December</td><td>S2-NDWI-win</td></tr><tr><td rowspan="6">S2 textural indexes (S2-BI-TI)</td><td>BI-homogeneity</td><td>S2-BI-sum-hom</td></tr><tr><td>BI-cluster shade</td><td>S2-BI-sum-clu</td></tr><tr><td>BI-correlation</td><td>S2-BI-sum-cor</td></tr><tr><td>BI-Haralick correlation</td><td>S2-BI-sum-Hcor</td></tr><tr><td>BI-homogeneity</td><td>S2-BI-win-hom</td></tr><tr><td>BI-cluster shade</td><td>S2-BI-win-clu</td></tr><tr><td rowspan="6">Winter image</td><td>BI-correlation</td><td>S2-BI-win-cor</td></tr><tr><td>BI-Haralick correlation</td><td>S2-BI-win-Hcor</td></tr><tr><td rowspan="6">Spot-XS textural indexes (Spot-BI-TI)</td><td>BI-homogeneity</td><td>Spot-BI-hom</td></tr><tr><td>BI-cluster shade</td><td>Spot-BI-clu</td></tr><tr><td>BI-correlation</td><td>Spot-BI-cor</td></tr><tr><td>BI-Haralick correlation</td><td>Spot-BI-Hcor</td></tr><tr><td rowspan="4">Spot-PAN textural indexes (Spot-PAN-TI)</td><td>PAN-homogeneity</td><td>Spot-PAN-hom</td></tr><tr><td>PAN-cluster shade</td><td>Spot-PAN-clu</td></tr><tr><td>PAN-correlation</td><td>Spot-PAN-cor</td></tr><tr><td>PAN-Haralick correlation</td><td>Spot-PAN-Hcor</td></tr></tbody></table>
ALOS-PALSAR-2: L-band backscatter data (~23 cm wavelength) are considered to be more sensitive to the woody components of the canopy, and thus to forest height and AGB. We use HV, HH and the HV/HH ratio from 2016. No texture metrics are extracted from ALOS-PALSAR because the spatial resolution is too coarse (~25 m): the window size would be too large compared to the size of some stands in these forests. Sentinel-1 and 2 provide a resolution more suitable for texture extraction.

Sentinel-1: although more limited than L-band, the penetration of C-band microwaves can be significant for temperate forests. C-band backscatter is sensitive to thorns, leaves and small branches. It provides information about vegetation structure and can be related to forest volume. SAR backscatter is sensitive to humidity, rain events can cause occasional increases in the signal regardless of vegetation evolution. For the purpose of ensuring that the results are reproducible, the selected dates meet two criteria: (1) one date every two months (six dates/year) and (2) no precipitation within the three previous days. This reduce the uncertainties caused by the rain, and increase the independence of the dates with climatic conditions and yearly variations. We assessed the suitability of Sentinel-1 temporal information for the estimation of forest structure parameters through correlation and regression analysis. The annual mean of the six dates provides results equivalent to the six dates together. As a consequence we keep the annual mean of VH, VV and VH/VV ratio. Details can be found in Appendix B. The three orbits (33°IA-DES, 43°IA-DES, 37°IA-ASC) are equally pertinent for the estimation of forest parameters, we retain the smallest incidence angle (33°). Consequently, textural indexes are calculated from the 33° IA orbit for the six dates, then averaged into annual values.

Sentinel-2: spectral reflectance is sensitive to green leaf area, vegetation cover fraction (trees and understory) and soil properties. Sentinel-2's observations provide information about vegetation development and photosynthetic activity. With view to build a generic method suitable for other sites and years (different species and climates), it is better to use winter and summer dates because the phenology of the vegetation is more stable during these periods (Appendix C). As a results we use 2016-08-12 for the summer and 2016-12-10 for winter period. We computed three spectral indexes from summer and winter dates: Normalized Difference Vegetation Index (NDVI), an indicator of vegetation activity; Normalized Difference Water Index (NDWI, using nir and swir bands), an indicator of leaves water content; and Brightness Index (BI, using visible bands), an indicator of global reflectance. Textural indexes were extracted from the 10 m resolution indexes (NDVI and BI) of summer and winter dates. Our tests show the textural indexes from BI are more efficient than the textural indexes from NDVI for the estimation of forest parameters (–10 to –20% relative RMSE depending on forest parameter studied). As a result, we keep summer and winter BI textural indexes for the following study.

Spot-6: Spot-6 mosaic images are used to test if a better spatial resolution for optical textural indexes can improve the results. As for Sentinel-2 images, we use the BI index to extract textural indexes on the multi-spectral bands at 6 m resolution. Additionally, the panchromatic band (sensible to the visible spectrum) is used to extract textural indexes at 1.5 m resolution.

Textural indexes in general: we computed eight features: energy, entropy, correlation, homogeneity, inertia, cluster shade, cluster prominence, and Haralick correlation. Energy and entropy indexes measure the level of "order" in the image (repetitive patterns). Homogeneity and inertia indexes give an indication of the local changes in intensity. Cluster shade and cluster prominence indexes indicate the tendency of clustering of the pixels in the window. Correlation indexes measure the dependency of grey levels on those of neighboring pixels (it means there is a predictable relationship between two neighboring pixels within the window). Energy, entropy, homogeneity and inertia are often strongly correlated. We chose to keep only the homogeneity index. Correlation and Haralick correlation indexes have no collinearity and shows significant correlations with forest parameters. The final textural indexes we consider for the following study are therefore homogeneity, correlation, cluster shade and Haralick correlation.
# 3. Methods and Analysis

## 3.1. Relationships Between Remote Sensing Features and Forest Structure Parameters

The linear correlations between the 36 remote sensing features (Table 3) and the six forest structure parameters can be divided in two groups of similar patterns: (1) AGB and BA, and (2) height, tree density, DBH and age. Figure 3 shows AGB and height absolute correlations with remote sensing features.

L-band wavelength penetrates tree crowns and is sensitive to trunk and large branches volume, and allows to retrieve AGB. L-HV polarization is more sensitive to wood volume, and L-HH to humidity. The results showed in Figure 3 confirm that our maritime pine plantations case study verifies these principles: $r > 0.7$ for both (L-HV, AGB) and (HV/HH-ratio, AGB) correlations. C-band wavelength is more sensible to crown elements (thorns and small branches). These elements are related to other forest structure parameters because allometric relationships are strong in these forests; this explains the high correlation of C-band VH polarization and textural indexes (-cor, -hom and -Hcor) with AGB. Moreover, C-band VH, VV and extracted textural indexes (VH and VV -cor, -hom and -Hcor) are highly correlated with height ($0.75 < r < 0.85$). Optical indexes from Sentinel-2 are poorly correlated with AGB and height. In contrast, some textural indexes extracted from BI show significant correlations with height ($r > 0.6$ for summer-hom, winter-hom, and winter-cor) and AGB ($r > 5$ for winter-hom). Springtime optical texture indexes from Spot-6 at 1.5 m spatial resolution show significant correlation with height ($r > 0.6$ for -cor and -hom) whereas the correlation is poor for both AGB and height at 6m spatial resolution.

This conclusion based on the correlations between remote sensing features and forest structure parameters highlights the capability of C-band and various textural indexes to describe linearly forest parameters on these plantations. Several studies showed high correlations between L-band SAR and AGB on different forest [25,26,34]; these correlations apply to our study site. The comparison between Sentinel-2 and Spot-6 shows that spring time 6 m spatial resolution image is less adequate than winter and summer 10 m spatial resolution images, whereas Spot's spring time 1.5 m spatial resolution image adds useful information.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** AGB and height Pearson's correlations with the remote sensing features presented in Table 3.</figcaption></figure>

## 3.2. Choice of Regression Algorithms and Parametrization

Statistical regressions have been used to establish the connections between the remote sensing features and AGB, BA, DBH, age, density and height variables. The choice of the regression algorithm is very important for the quality of the results. We investigated non parametric and multiple regression methods that could perform with several input predictors and a small number of samples. Deep learning approaches like neural network for regression are not relevant for a small dataset. Multi-linear regressions (MLRs) are very simple and frequently used in papers related to forest variables
estimation using remote sensing. MLRs assume that predictive variables have linear relationships with predicted variables, although it has been shown it is not always an acceptable approximation when applied to remote sensing features and forest structure parameters. Machine learning approaches like support vector machines and random forest algorithms are increasingly and successfully used in remote sensing domain [63–65]. In this study we tested MLR, random forest regression (RF) and support vector regression (SVR) with a Gaussian kernel (radial basis function, RBF), as implemented in the Python scikit-learn library [66].

SVR requires three key parameters: (1) the cost parameter (C) trades off misclassification of training samples against simplicity of the decision surface. Tested values are {0.1, 1, 5, 10, 50, 500, 1000, 1500}. (2) The gamma coefficient (G) defines how much influence a single training sample has. Tested values are {0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 20} and (3) epsilon (E) specifies the epsilon-tube within which no penalty is associated in the training loss function. Tested values are {0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 20}.

RF regression parameters are (1) number of trees {50, 100, 250}, (2) depth of trees {10, 25, None}, and (3) minimum impurity decrease {0, 0.1, 0.2}. The number of features considered at each split was set to 30%.

MLR has no parameter to be set.

### 3.3. Feature Selection Process (Dimensionality Reduction)

Although we selected 28 out of 215 features, there are still features inter-correlated and thus not necessarily useful. Moreover, too many features could lead to a degradation of results or over-fitting. To avoid that, we tested three approaches:

**Principal Component Analysis (PCA):** the first solution is to reduce the dimension of the data using PCA method. We test the PCA on all the 28 features together or on the seven feature types separately. We keep the number of principal components (PC) which correspond to 95% of the explained variance.

**Forward selection:** we start with the best efficient feature using the regression algorithm, and then for each iteration we choose which feature could be added providing the lower RMSE, until all features are selected. Finally, within these feature combinations, we select the one that provides the best result.

**Backward selection:** it is the reverse process, we start with all the features as inputs for the regression algorithm, and then we choose which feature could be removed while keeping the RMSE minimal, until only one feature remains. Again, within these feature combinations, we select the one that provides the best result.

### 3.4. Validation

The statistical significance of the tests stems from a large enough validation subset while the pertinence of the solutions is built on a learning subset that is large enough with regard to the features. In our case study, the validation procedure was sensitive to the small number of samples. In order to maximize the number of learning samples and to minimize statistical variance, we used the leave-one-out cross validation method (LOO) where each sample is estimated using all the other samples. We assessed the quality of the results through the following statistics: (a) the coefficient of determination (r²) reflects the variance explained by the regression. (b) The root mean squared error (RMSE) provides an indicator of the estimation errors, giving significant weight to large errors. (c) The relative RMSE is the division by the mean of reference samples. (d) The mean absolute error (MAE) is also an indicator of estimation errors but gives the same weight to all errors, the relative MAE is the division by the mean of reference samples.
# 4. Results

## 4.1. Method Selection

### 4.1.1. Machine Learning Algorithm

MLR assumes that remote sensing features and forest variables have linear relationships, it is very fast and simple to use. Random forest (RF) regression is a non-linear and non-parametric algorithm, easy to interpret. However, the small number of learning samples is challenging for the construction of independent trees in the random forest. SVR (RBF kernel) is also a non-linear and non-parametric algorithm. This algorithm can deal with small numbers of samples. Figure 4a shows the relative LOO RMSE (smaller values are better) obtained for MLR, RF and SVR using as input the 28 features from open access images with a forward feature selection. These three algorithms produce equivalent results for AGB and BA estimations. Differences are more pronounced for the estimation of DBH, age, density and height. The best estimates are obtained using SVR, RF and MLR in this order. The results confirm the benefits of using non-linear and non-parametric algorithms, and SVR is better at handling small datasets.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparison of algorithms and dimensionality reduction for the estimation of AGB, DBH, Density and Height. (a) Comparison between MLR, RF and SVR using the forward feature selection method of the 28 features from open access images. (b) Comparison between PCA, backward selection, forward selection or no selection using SVR on the 28 features from open access images.</figcaption></figure>

### 4.1.2. Feature Selection Approaches

We study the effectiveness of dimensionality reduction for improving the estimation of forest parameters from the perspective of four approaches: PCA, forward selection, backward selection, and no selection. The evaluation is performed with SVR on 28 features from Sentinel-1, Sentinel-2 and ALOS-PALSAR. Results are presented in Figure 4b. Backward selection was expected to perform best for it maintains groups of features and draws on their synergies. Expectations toward forward selection were lower since the approach adds one feature at a time and is known to dismiss prematurely features that have a beneficial contribution. Yet forward selection performed better and is our recommendation when dealing with small reference dataset.

Figure 5 shows how the relative LOO RMSE evolves given the number of features selected by a forward algorithm in order to estimate forest parameters.

Figure 5 shows that seven features are the “sweet spot” to minimize the RMSE. Additional features marginally reduce the RMSE or decrease it. We note that both estimates of AGB and BA are less dependent on the number of features than the other variables.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Evolution of relative RMSE with the forward feature selection method for AGB, BA, DBH, age, Density and Height. Vertical lines correspond to the lowest RMSE for each forest parameters.</figcaption></figure>

## 4.2. Forest Parameters Estimation

This section presents the estimates of AGB, BA, DBH, age, tree density and height computed using 28 features belonging to the five feature types (L-band, C-band, C-TI, S2-SI, S2-BI-TI) and both SVR regression and forward selection methods. Table 4 indicates which features are chosen by the algorithm to obtain the results for the six forest parameters. Validation scatterplots are presented in the following subsections.

<table><caption>**Table 4.** Summary of features selected for the estimation of forest parameters. The processing chain uses 5 feature types as input (L-band, C-band, C-TI, S2-SI, S2-BI-TI) and forward selection method with SVR regression. Numbers in parenthesis indicate the position order in the forward selection.</caption><thead><tr><th></th><th>L</th><th>C</th><th>C-TI</th><th>S2-SI</th><th>S2-BI-TI</th><th>r<sup>2</sup> and RMSE</th></tr></thead><tbody><tr><td rowspan="3">AGB</td><td>L-HV (1)</td><td></td><td rowspan="3">C-VH-clu (8)</td><td>S2-NDVI-sum (2)</td><td rowspan="3">S2-BI-win-hom (6)</td><td>r<sup>2</sup> = 0.76</td></tr><tr><td>L-HH (5)</td><td></td><td>S2-BI-sum (3)</td><td>19.5 tons/ha</td></tr><tr><td>L-Ratio (9)</td><td></td><td>S2-NDWI-sum (4)</td><td>(28.0%)</td></tr><tr><td rowspan="2">BA</td><td>L-HV (1)</td><td></td><td rowspan="2">C-VH-cor (5)</td><td rowspan="2">S2-BI-sum (3)</td><td>S2-BI-sum-hom (4)</td><td>r<sup>2</sup> = 0.73</td></tr><tr><td>L-HH (2)</td><td></td><td>S2-BI-sum-Hcor (7)</td><td>5.7 m<sup>2</sup>/ha</td></tr><tr><td rowspan="3">DBH</td><td rowspan="3">L-HV (2)</td><td rowspan="3">C-Ratio (6)</td><td>C-VV-cor (1)</td><td></td><td rowspan="3">S2-BI-sum-cor (3)</td><td>r<sup>2</sup> = 0.88</td></tr><tr><td>C-VH-clu (4)</td><td></td><td>0.04 m</td></tr><tr><td>C-VH-cor (5)</td><td></td><td>(19.8%)</td></tr><tr><td rowspan="4">Age</td><td rowspan="4">L-HV (2)</td><td rowspan="4">L-HH (5)</td><td>C-VV-hom (1)</td><td></td><td rowspan="4">S2-BI-sum-hom (6)</td><td>r<sup>2</sup> = 0.93</td></tr><tr><td>C-VH-Hcor (3)</td><td>S2-BI-win (9)</td><td>3.95 years</td></tr><tr><td>C-VV-cor (4)</td><td>S2-BI-sum (10)</td><td>(17.4%)</td></tr><tr><td>C-VH-clu (8)</td><td></td><td></td></tr><tr><td rowspan="3">Density</td><td rowspan="3">L-HV (2)</td><td rowspan="3">C-Ratio (5)</td><td>C-VV-Hcor (1)</td><td></td><td rowspan="3">S2-BI-sum-hom (3)</td><td>r<sup>2</sup> = 0.86</td></tr><tr><td>C-VH-hom (4)</td><td></td><td>204 trees/ha</td></tr><tr><td>C-VH-Hcor (6)</td><td></td><td>(24.4%)</td></tr><tr><td rowspan="6">Height</td><td rowspan="6">L-HV (4)</td><td rowspan="6">L-HH (10)</td><td>C-VH-Hcor (1)</td><td></td><td rowspan="6">S2-NDWI-win (2)</td><td>r<sup>2</sup> = 0.93</td></tr><tr><td>C-VV-clu (8)</td><td></td><td>S2-BI-sum-cor (3)</td><td>1.75 m</td></tr><tr><td>C-VH-clu (11)</td><td></td><td>S2-BI-sum-hom (5)</td><td>(13.2%)</td></tr><tr><td>C-VV-cor (12)</td><td></td><td>S2-BI-win-hom (6)</td><td></td></tr><tr><td>C-VH-cor (13)</td><td></td><td>S2-BI-sum-Hcor (7)</td><td></td></tr><tr><td>C-VV-hom (14)</td><td></td><td></td><td></td></tr></tbody></table>

The features presented in Table 4 are generated by the SVR algorithm with the features obtained by applying the forward method. The forward selection singles six to 14 features out of 28 (five types) for the estimation of AGB, BA, DBH, age, density and height. The evolution of relative RMSE with forward selection presented in Figure 5 reveals that results marginally improve for any feature added to a core of seven, and eventually degrade the performance. At least four of the five feature types are represented for each forest parameter estimate. Comparison with Table 4 and Figure 3 shows that the features presenting higher linear correlation with a given forest parameter are not necessarily leveraged by the algorithm. For AGB estimates, the algorithm exploits four spectral indexes (S2-SI, low correlations with AGB) while only L-HV and L-Ratio within the group of features whose r > 0.5. The same can be
observed among the features selected for height estimation. We see on second and third position in the selection order in Table 4 two features weakly correlated with height (S2-NDWI-win and S2-BI-sum-cor, respectively 0.1 and 0.3).

As it is well known, multiple independent data sources relevant for the problem at hand are a key factor to obtain a good model; this is what we observe here. In general, for the whole forest variables, L-band and textural metrics from C-band (C-TI) are very important, often present in the first selected features. C-band backscatter is less useful, maybe because VH and VV backscatters are highly correlated with their textural indexes (-Hcor, r > 0.95). Sentinel-2 spectral indexes and textural indexes are also often selected by the forward method although they have fewer correlations with forest parameters (particularly S2-SI). This further illustrates the importance of taking into account the synergy between features and not just individual performance to build the best models for forest parameters estimation.

### 4.2.1. Estimation of Aboveground Biomass (AGB)

Figure 6a shows the AGB predictions versus AGB references. The RMSE is 19.5 tons/ha (28%) and the absolute error is 22.7%. Figure 6c shows that the variance of the errors is approximately constant with the increase of the AGB, but the positive or negative orientation of these errors is not random. Indeed, Figure 6c shows the errors are negative from 100 tons/ha. This is due to signal saturation for AGB estimation. The analysis of the SVR model and the high Cost parameter (C = 500) selected suggest that the model gives high weights to each sample, and has trouble generalizing. This could be improved either by a larger reference dataset or in the future by adding new features that have stronger relationships with the AGB (P-band radar for example).

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Estimation of AGB using SVR (C = 500, G = 0.001, E = 0.001) and forward selection method from L-band, C-band, C-TI, S2-SI, S2-BI-TI (selected features presented in Table 4). (a) Validation scatterplot; (b) AGB reference samples distribution; (c) Estimation errors depending AGB on reference values.</figcaption></figure>

### 4.2.2. Estimation of Basal Area (BA)

The estimation of basal area shows very similar results to AGB estimation. The RMSE is 5.7 m²/ha (26.9%) and the absolute error is 20.3%. Figure A1 in Appendix A shows that medium values (15–20 m²/ha) are overestimated and capped under 30–35 m²/ha.

### 4.2.3. Estimation of Diameter at Breast Height (DBH)

Figure 7 shows very good results for DBH predictions with a RMSE of 4 cm (19.8%) and a r² of 0.88. Most of the reference samples are below 40 cm DBH (Figure 7b). The errors are well distributed along the DBH values, there is no saturation. Forest volume information (L-band features and C-Ratio), combined to forest structure information (S1 and S2 textural metrics) are able to provide very good mean DBH prediction at high spatial resolution in the coniferous plantations of our test study.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Estimation of mean DBH using SVR (C = 10, G = 0.01, E = 0.1) and forward selection method from L-band, C-band, C-TI, S2-SI, S2-BI-TI (selected features presented in Table 4). (a) Validation scatterplot; (b) DBH reference samples distribution; (c) Estimation errors depending on DBH reference values.</figcaption></figure>

### 4.2.4. Stand Age Estimation

Age estimation results (Figure A2 in Appendix A) are very similar to the DBH estimation results. This was expected since both variables are usually highly correlated in plantation forests. We obtained a RMSE of ~4 years (MAE = 2.7 years). There is no saturation. The oldest plantations (three samples > 65 years old) are underestimated, but most of the maritime pine plantation are usually clear cut when they get 45–55 years old. The samples above this age could have different management practices, and the age estimation method in reference data (counting the apparent number of annual growth cycles with branch whorls) could introduce some approximation for very high ages. Table 4 and Figure 5 show more features are needed to obtain the best predictions compared to DBH.

### 4.2.5. Tree Density Estimation

Results related to tree density estimates are presented in Figure 8. The RMSE is 204 trees/ha (24.4%), the relative MAE is 17.1%. The errors distribution shows no saturation. Results are therefore very good while this is an example of forest parameter for which there is no obvious correlation with other forest parameters (r² from 0.04 to 0.64).

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Estimation of tree density using SVR (C = 10, G = 0.1, E = 0.01) and forward selection method from L-band, C-band, C-TI, S2-SI, S2-BI-TI (selected features presented in Table 4). (a) Validation scatterplot; (b) Density reference samples distribution; (c) Estimation errors depending on density reference values.</figcaption></figure>
### 4.2.6. Dominant Height Estimation

Results related to forest height estimates are presented in Figure 9. The RMSE is 1.75 m (13.2%), and the MAE is 1.4 m. Figure 9c shows that the errors are well distributed and there is no outliers.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Estimation of mean tree height using SVR (C = 10, G = 0.01, E = 0.1) and forward selection method from L-band, C-band, C-TI, S2-SI, S2-BI-TI (selected features presented in Table 5). (a) Validation scatterplot; (b) Height reference samples distribution; (c) Estimation errors depending on height reference values.</figcaption></figure>

### 4.3. Mapping

Remote sensing images allow to spatialize forest parameters estimates. Figure 10 shows the AGB map of the studied area. We stack the features used for AGB prediction (Table 4) and we make a pixel-level application of the model learned (validation scatterplot presented in Figure 6). In order to have proper predictions, we mask non coniferous forested areas with the land cover map produced in 2016 with Sentinel-2 time series [67]. Within the forest classes, the nomenclature differentiates coniferous and broadleaved forests. With the help of ancillary land cover data, our method is able to produce high spatial resolution maps of forest parameters. We note that textural indexes that are sensitive to boundary conditions create unwanted effects. For example, AGB should decrease where there are pathways between the stands, but does rise instead.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Aboveground biomass (AGB) map covering the study site area. Pixel level application of the model using 5 feature types (L-band, C-band, C-TI, S2-SI, S2-BI-TI). Selected features are presented in Table 4 and validation scatterplot in Figure 6. Non coniferous pixels are masked using the 2016 land cover map [67].</figcaption></figure>
## 4.4. Other Feature Types and Combination Analysis

The goal of this section is to assess the robustness of our method with different sets of feature types. We test four types conditions: (1) single feature type, (2) combination of 2 feature types, (3) combination of open access data types only (whose results are presented Section 4.2), and (4) the addition of Spot-6 multi-spectral and panchromatic textural indexes (respectively 6m and 1.5m spatial resolution) to the open access feature types. Figure 11 and Table 5 present the relative RMSE (%) obtained for the estimation of the different forest parameters according to the four types conditions.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Relative RMSE (%) of the estimations of the six forest parameters using feature types individually (the seven first columns) and combined: two feature types (best efficient pairs, details in Table 4), the five feature types from open access data (L-band, C-band, C-TI, S2-SI, S2-TI), and 6 feature types (all feature types except S2-TI).</figcaption></figure>

<table><caption>**Table 5.** Best Relative RMSE (%) obtained for the estimations of the six forest parameters. Feature types individually (first line) and combined: 2 feature types (second line, best pair indicated), the 5 feature types from open access data (third line), and addition of Spot-6 textural indexes (last line).</caption><thead><tr><th></th><th>AGB</th><th>BA</th><th>DBH</th><th>Age</th><th>Density</th><th>Height</th></tr></thead><tbody><tr><th>Best feature type</th><td>L-band 32%</td><td>L-band 29%</td><td>C-TI 26%</td><td>C-TI 27%</td><td>C-band 37%</td><td>C-TI 18%</td></tr><tr><th>Best two types combination</th><td>L + S2-BI-TI<br>28.0%</td><td>L + C<br>27.1%</td><td>C-TI + S2-BI-TI<br>21.2%</td><td>L + C-TI<br>18.6%</td><td>L + S2-BI-TI<br>30.7%</td><td>C-TI + S2-BI-TI<br>14.3%</td></tr><tr><th>Five types open access images</th><td>28.0%</td><td>26.9%</td><td>19.8%</td><td>17.4%</td><td>24.4%</td><td>13.2%</td></tr><tr><th>Seven types (Spot-6 included)</th><td>27.9%</td><td>27.0%</td><td>19.3%</td><td>18.6%</td><td>23.3%</td><td>13.7%</td></tr></tbody></table>

### 4.4.1. Best Single Feature Type

For AGB and BA, L-band backscatter is the best feature to use, but only slightly more efficient than the Sentinel-1 C-band backscatter or textural indexes. Indeed, L-band wavelength (~27 cm) is able to penetrate tree crowns so it is sensible to trunks and large branches, allowing to estimate forest basal area and biomass. The performances of the L-band decrease strongly on DBH, age, density and height. For these forest parameters, Sentinel-1 C-band textural indexes are the best indicators, followed by C-band backscatter, Sentinel-2 and Spot panchromatic textural indexes. Sentinel-1 SAR C-band data is therefore the best pick when using only one data type. C-band (~5 cm wavelength) penetrates the canopy and is sensible to the volume of crown and small branches. As a consequence, it is generally well studied to estimate biomass on young forests below 50 tons/ha, but the signal quickly saturates in older forests. On these coniferous plantations of our study site however, allometric relationships between crown biomass and total tree biomass are strong and Sentinel-1 C-band data is pertinent to estimate forest biomass through the crown structure. We note that spatial information provided by textural indexes is a good proxy for the estimation of DBH, height, and for AGB as well. We also note that Sentinel-2 spectral indexes and Spot-6 multi-spectral textural indexes show poorer performances
if used alone. According to the results of these feature types, Sentinel-1 VH and VV textural indexes (10 m) are the most efficient data for the description of forest spatial structure.

### 4.4.2. Best Combination of Two Feature Types (Results of the Best Pairs)

The combination of two feature types improves clearly the quality of the estimations for all forest parameters over a single type scenario. The results are close to the best obtained when all the feature types are combined, but for tree density estimates. The best feature type (first line of Table 5) is generally present in the best pair selected. The second type is not necessarily within the best individual feature types; there are examples where the worst efficient type (the one at the end of the ranking list) is used: when used alone, L-band backscatter gives the worst results for age, DBH, height and density estimations, however the SVR algorithm combines it with the C-TI features and yields noticeably improved estimates. L-band and C-band are complementary when we are concerned with building biomass estimates: L-band wavelength saturates at high biomass values (~100–150 tons/ha) and C-band wavelength saturates from 50 tons/ha, but their combined use improves the quality of the estimates for the whole range 50-150 tons/ha.

### 4.4.3. Addition of Spot-6 Textural Metrics

The use of optical textural metrics with a better spatial resolution is not relevant in the context of our study site. Either this better resolution does not bring additional information on these forests, or the date of the Spot-6 image (early spring) is not suited. However, further studies on other forest types would be needed to conclude on these intermediate spatial resolutions between HR (~10–20 m) and VHR (~0.5–1 m).

## 5. Discussion

Previous studies demonstrated the usefulness of open access spaceborne data such as Sentinel-1, Sentinel-2 and ALOS-PALSAR mosaics when applied to different forest types and data combinations. In this paper, we build on these previous efforts in order to harmonize the methods and build a generic approach to produce forest parameters maps. We gather all open access data shown to be useful across different forest types and parameter estimations (SAR L- and C-band, optical indexes, and textural indexes) and develop an automated processing chain that allows to retrieve and spatialize forest key parameters accurately. The method is modular and can be easily tailored to different contexts. The processing chain has been developed with Python (especially scikit-learn library [66]) and Orfeo Toolbox (OTB) [62] open source software for image processing and machine learning algorithms. Performance has been evaluated on maritime pine forests within the Landes de Gascogne in Southwest France.

We estimate six forest parameters: AGB (28% relative RMSE), BA (27%), DBH (20%), age (17%), density (24%) and height (13%). Height and DBH and age are key variables for forest monitoring, highly related with thinning and harvesting operations. They are predicted with the best accuracy. BA and AGB are important for resources and carbon budget assessment. Their accuracy is lower. One reason can be that the remote sensing data we use are known to saturate for high biomass. Moreover, considering the high spatial resolution of the images, the biomass is more variable than DBH or height from one pixel to another because of tree cover gaps, especially for sowing afforestation. Image segmentation may assist with dealing with such heterogeneous areas.

As mentioned in the literature specialized remote sensing of forest parameters, the best accurate estimations are provided by airborne laser scanning (ALS) and SAR radargrammetry [5,68,69] with about 15–20% RMSE for AGB estimation, and 5–10% RMSE for height estimation on tropical, boreal and temperate forests. Estimations based on photogrammetry and textural metrics from very high resolution optical images (0.5–1 m) also provide very accurate results [40–42,68] with 20–25% and 10–15% relative RMSE for AGB and height. The estimates we produce on the basis of global coverage and freely available images, with open source software, are close to the best results published so far.
Out of the feature selection and machine learning algorithms we tested, SVR algorithm and forward feature selection are recommended when reference samples are scarce (<100). RF algorithm can be a good option for larger datasets. On these maritime pine forests, the analysis of feature types contributions reveal that the Sentinel-1 textural metrics are the most useful feature type when used alone for the estimation of the ensemble of forest variables. However, L-band is still the better to estimate specifically AGB and BA. With respect to the synergy, the best feature types combination is L-band SAR backscatter with textural indexes from Sentinel-1 or Sentinel-2. Nevertheless, other features can be useful for generalization and it is better to keep a wide range of feature types in order to apply the method efficiently to other forests. In addition, other predictive features can be used when available in the study area. We tested Spot-6 annual mosaics that are freely available on France territory. It would be possible to feed the processing chain with more precise data (airborne, drones) or information on soils and slopes, according to what is available and relevant on the targeted forests.

The method is designed to be generic, taking advantage of a large variety of remote sensing data that have been used in several studies on different types of forests. Although this first validation was made on a single type of forest, the method can be applied to other forest types and perform as well as the studies that have used the same data on other forests. We are currently testing the application on more complex forests (uneven-aged, mixed broadleaved and coniferous). Preliminary results confirm the approach is robust. Nevertheless, difficulties could appear in dense forests (due to the signal saturation of current satellite sensors) and in mountainous areas with slopes above 20° (due to inherent radar geometry distortions). Furthermore, the relationships between remote sensing data and forest parameters are expected to change from one forest to another according to the management, species, environmental conditions or even the fertility of the soil. Temperate forests can be very different from one stand to another in a neighboring area. Information on the forest types can therefore be important for learning and applying models. In order to produce accurate maps of forest parameters over large areas in temperate regions, the differentiation between models should be related to current work on land cover and forest species classification. The versatility of the processing chain is ground to continue validation effort on other forests, and to analyze its strengths and limitations.

## 6. Conclusions

The results of this study confirm that open access spaceborne data, such as Sentinel-1, Sentinel-2 and ALOS-PALSAR mosaics (all three acquired worldwide at high spatial resolution), are able to provide accurate estimations of forest structure parameters and AGB. Based on open source software, we built a processing chain that uses SAR L- and C-band, optical spectral indexes, and textural indexes, all together, in a synergistic way to produce quantitative maps of forest parameters. The method is intended to be generic. A quality assessment was done on maritime pine stands within the largest European plantation forest. We are currently applying and validating the method on other temperate forests (uneven-aged, broadleaved and mixed forests). Once we have thoroughly analyzed the strengths and limitations of this method on different forest types, the key tools of the processing chain could be packaged and made available for forest organizations, national bodies or other researchers to test the method anywhere. In recent decades many forest organizations (e.g., in France the Office Nationale des Forêts, ONF) have established networks of permanent plots in order to improve forest knowledge and evaluate silvicultural practices. In addition, regular field measurements are performed to assess the financial value of forest stands. These measurements could provide a solid learning database for the method presented in this paper. In return, the possibility of mapping forest structure parameters such as DBH, height and BA with remote sensing data would be useful to obtain quantitative information on stands without measurements, and to improve the distribution and representativeness of permanent plots. The methodology we presented could also be used directly on the geolocated measurements performed by the national forest inventories in order to improve the sampling design and the spatial resolution of their maps. Quantitative mapping could then be used at both the stand level and over large areas to monitor and better assess forest management strategies. In addition, other information
could be extracted from these maps such as ecological indicators on the landscapes, a better estimate
of the consequences of fires, storms and other sanitary disasters, and a finer spatialization of carbon
budget assessment.

**Author Contributions:** D.M. wrote the manuscript, prepared and accomplished the study. M.P. supervised the study and contributed to data processing, results analysis and review. D.G. allowed the acquisition of the reference dataset and contributed to data analysis and review. L.V., S.M. and A.B. contributed to data processing, results analysis and review. H.T. contributed to conceptualization and review. J.-F.D. helped for funding acquisition and contributed to conceptualization. T.L.T. contributed to results analysis, conceptualization and review. G.D. supervised the research activity and contributed to conceptualization and review.

**Funding:** This research is part of a PhD funded by The French spatial agency (CNES) and the Agence de l'Environnement et de la Maîtrise de l'Énergie (ADEME).

**Acknowledgments:** Thanks to INRA for making available the *in situ* forest database. We are grateful for data quality to Sylvia Dayau (INRA, UMR ISPA) and Bernard Issenhut (INRA, Unité expérimentale Forêt Pierroton, 69 route d'Arcachon, 33612 CESTAS, France).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

**Appendix A Supplementary Scatterplot Results**

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

**Figure A1.** Estimation of BA using SVR (C = 500, G = 0.001, E = 0.001) and forward selection method from L-band, C-band, C-TI, S2-SI, S2-BI-TI (selected features presented in Table 5). (a) Validation scatterplot; (b) BA reference samples distribution; (c) Estimation errors depending on reference values.

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

**Figure A2.** Estimation of stand age using SVR (C = 100, G = 0.01, E = 0.001) and forward selection method from L-band, C-band, C-TI, S2-SI, S2-BI-TI (selected features presented in Table 5). (a) Validation scatterplot; (b) Age reference samples distribution; (c) Estimation errors depending on age reference values.
## Appendix B Analysis on the Sentinel-1 Temporal Information

SAR backscatter is affected by the variations of environmental effects such as soil and vegetation moisture. We have analyzed the impact of rainfalls (using Global Satellite Mapping of Precipitation, GSMaP, provided by JAXA) and wetness (using Soil Water Index, SWI, provided by Copernicus Global Land Service) on the backscatter temporal evolution. We have found that the signal variations are mostly concomitant with rain events observed on GSMaP. Furthermore, significant positive correlations are found between Sentinel-1 time series and the SWI on forest stands with different conditions of tree density and biomass range. As we want to study the evolution of the signal with regard to vegetation, the impact of rain and soil moisture might be a nuisance. This is the rationale for being selective about the dates used in this context.

Thirty Sentinel-1 dates are available in 2016 for each orbit; we choose one date every 2 months, and with no precipitation within the three previous days (as observed on GSMaP for each descendant orbit). We study whereas the temporal information of the time series is useful in the context of coniferous study site. Figure A3 shows the correlations of AGB and DBH with Sentinel-1 SAR backscatter (VH, VV polarizations and VH/VV ratio; orbit with 33° of the mean incidence angle) for the six dates and the annual mean. VH backscatter sensitivity to AGB and particularly to DBH is rather stable in time. VV backscatter is more sensitive to soil conditions and more correlated to AGB and DBH in the summer (July–August), in sharp contrast with the VH/VV ratio. We average the 6 dates in an annual mean (showed on the right of Figure A3a,b). When averaging, VV correlation with AGB (Figure A3a) is decreased by 0.1, and the correlation of the annual mean of VH/VV with DBH (Figure A3b) is also decreased by 0.1. Overall the annual mean shows quite similar correlations to the dates within the year with the best scores. Moreover, we test the estimation of forest variables with Sentinel-1 dates or annual mean using SVR and forward selection method. The use of six dates (18 features for one orbit, instead of three with the annual mean) improves the results by ~1% for the AGB and density, but not for the other variables. Given these very slight improvements, we choose to keep only the annual mean values of the S1 radar backscatter for this study.

<figure><img src="image_17.png" /><figcaption>**Figure A3.** Correlation coefficients between Sentinel-1 backscatters (VH, VV, VH/VV; orbit 33°IA) and forest parameters for the 6 dates in 2016 and the annual mean. (a) correlations with AGB; (b) correlations with DBH.</figcaption></figure>
## Appendix C Analysis on the Sentinel-2 Temporal Information

Optical images are affected by cloud cover; the signal does not pass through clouds. Sentinel-2A acquires one image every 10 days, but when there is more than 90% of cloud cover, image processing can't be applied and the image is not delivered (ortho-rectification for level 1A cannot be applied). Images may still show clouds that cover part or all of the study area. For the Landes forest test site, 18 images are available in 2016 and only seven are cloud free on the reference samples. The dates and the number of usable images may change between years (see Figure A4) and between regions. We want to select dates that are representative of different vegetation stages, while being the most stable (based on the phenology) between years, regions or forest species. Figure A4 shows the temporal evolution of NDVI, BI and NDWI. Broadleaved and maritime pine forests are selected in order to compare different phenology of the vegetation. The analysis of time series reveals that vegetation cycles are repeated over the years. However, non-cloudy image dates may be very different between years. Rapid changes occur in the spring and autumn, while spectral indexes are more stable in the summer and winter. The period when coniferous and broadleaved forests differ the most is in winter. Therefore we choose to keep two cloud free images, one close to January (winter season) and another close to July 2016 (summer season).

<figure><img src="image_18.png" /><figcaption>**Figure A4.** Sentinel-2 spectral indexes times series on test site broadleaved (green lines) and maritime pine (black lines) closed forests. (a) NDVI time series; (b) BI time series; (c) NDWI time series. Vertical lines represent non cloudy acquisition dates (on our reference samples).</figcaption></figure>

## References

1. Food and Agriculture Organization of the United Nations. *Global Forest Resources Assessment*; FAO: Rome, Italy, 2015; ISBN 978-92-5-108826-5.

2. Dixon, R.K.; Brown, S.; Houghton, R.A.; Solomon, A.M.; Trexler, M.C.; Wisniewski, J. Carbon Pools and Flux of Global Forest Ecosystems. *Sci. New Ser.* **1994**, *263*, 185–190. [CrossRef] [PubMed]

3. Jochem, A.; Hollaus, M.; Rutzinger, M.; Höfle, B. Estimation of Aboveground Biomass in Alpine Forests: A Semi-Empirical Approach Considering Canopy Transparency Derived from Airborne LiDAR Data. *Sensors* **2010**, *11*, 278–295. [CrossRef] [PubMed]
4. Gobakken, T.; Næsset, E.; Nelson, R.; Bollandsås, O.M.; Gregoire, T.G.; Ståhl, G.; Holm, S.; Ørka, H.O.; Astrup, R. Estimating biomass in Hedmark County, Norway using national forest inventory field plots and airborne laser scanning. *Remote Sens. Environ.* **2012**, *123*, 443–456. [CrossRef]

5. D'Oliveira, M.V.N.; Reutebuch, S.E.; McGaughey, R.J.; Andersen, H.-E. Estimating forest biomass and identifying low-intensity logging areas using airborne scanning lidar in Antimary State Forest, Acre State, Western Brazilian Amazon. *Remote Sens. Environ.* **2012**, *124*, 479–491. [CrossRef]

6. Bouvier, M.; Durrieu, S.; Fournier, R.A.; Renaud, J.-P. Generalizing predictive models of forest inventory attributes using an area-based approach with airborne LiDAR data. *Remote Sens. Environ.* **2015**, *156*, 322–334. [CrossRef]

7. Ho Tong Minh, D.; Le Toan, T.; Rocca, F.; Tewaldini, S.; Villard, L.; Réjou-Méchain, M.; Phillips, O.L.; Feldpausch, T.R.; Dubois-Fernandez, P.; Scipal, K.; et al. SAR tomography for the retrieval of forest biomass and height: Cross-validation at two tropical forest sites in French Guiana. *Remote Sens. Environ.* **2016**, *175*, 138–147. [CrossRef]

8. Wulder, M.A.; White, J.C.; Nelson, R.F.; Næsset, E.; Ørka, H.O.; Coops, N.C.; Hilker, T.; Bater, C.W.; Gobakken, T. Lidar sampling for large-area forest characterization: A review. *Remote Sens. Environ.* **2012**, *121*, 196–209. [CrossRef]

9. Tomppo, E. Using coarse scale forest variables as ancillary information and weighting of variables in k-NN estimation: A genetic algorithm approach. *Remote Sens. Environ.* **2004**, *92*, 1–20. [CrossRef]

10. Mcroberts, R.; Tomppo, E. Remote sensing support for national forest inventories. *Remote Sens. Environ.* **2007**, *110*, 412–419. [CrossRef]

11. Koch, B. *Remote Sensing Supporting National Forest Inventories*; NFA: Chicago, IL, USA, 2013; p. 15.

12. Fernandez-Ordonez, Y.; Soria-Ruiz, J.; Leblon, B. Forest Inventory using Optical and Radar Remote Sensing. *Adv. Geosci. Remote Sens.* **2009**. [CrossRef]

13. Goetz, S.J.; Baccini, A.; Laporte, N.T.; Johns, T.; Walker, W.; Kellndorfer, J.; Houghton, R.A.; Sun, M. Mapping and monitoring carbon stocks with satellite observations: A comparison of methods. *Carbon Balance Manag.* **2009**, *4*. [CrossRef] [PubMed]

14. Dittmann, S.; Thiessen, E.; Hartung, E. Applicability of different non-invasive methods for tree mass estimation: A review. *For. Ecol. Manag.* **2017**, *398*, 208–215. [CrossRef]

15. Hansen, M.C.; Potapov, P.V.; Goetz, S.J.; Turubanova, S.; Tyukavina, A.; Krylov, A.; Kommareddy, A.; Egorov, A. Mapping tree height distributions in Sub-Saharan Africa using Landsat 7 and 8 data. *Remote Sens. Environ.* **2016**, *185*, 221–232. [CrossRef]

16. Chi, H.; Sun, G.; Huang, J.; Li, R.; Ren, X.; Ni, W.; Fu, A. Estimation of Forest Aboveground Biomass in Changbai Mountain Region Using ICESat/GLAS and Landsat/TM Data. *Remote Sens.* **2017**, *9*, 707. [CrossRef]

17. Deo, R.; Russell, M.; Domke, G.; Andersen, H.-E.; Cohen, W.; Woodall, C. Evaluating Site-Specific and Generic Spatial Models of Aboveground Forest Biomass Based on Landsat Time-Series and LiDAR Strip Samples in the Eastern USA. *Remote Sens.* **2017**, *9*, 598. [CrossRef]

18. Liu, K.; Wang, J.; Zeng, W.; Song, J. Comparison and Evaluation of Three Methods for Estimating Forest above Ground Biomass Using TM and GLAS Data. *Remote Sens.* **2017**, *9*, 341. [CrossRef]

19. Bolton, D.K.; White, J.C.; Wulder, M.A.; Coops, N.C.; Hermosilla, T.; Yuan, X. Updating stand-level forest inventories using airborne laser scanning and Landsat time series data. *Int. J. Appl. Earth Obs. Geoinf.* **2018**, *66*, 174–183. [CrossRef]

20. Luckman, A.; Baker, J.; Honzák, M.; Lucas, R. Tropical Forest Biomass Density Estimation Using JERS-1 SAR: Seasonal Variation, Confidence Limits, and Application to Image Mosaics. *Remote Sens. Environ.* **1998**, *63*, 126–139. [CrossRef]

21. Castel, T. Retrieval biomass of a large Venezuelan pine plantation using JERS-1 SAR data. Analysis of forest structure impact on radar signature. *Remote Sens. Environ.* **2002**, *79*, 30–41. [CrossRef]

22. Mitchard, E.T.A.; Saatchi, S.S.; Lewis, S.L.; Feldpausch, T.R.; Woodhouse, I.H.; Sonké, B.; Rowland, C.; Meir, P. Measuring biomass changes due to woody encroachment and deforestation/degradation in a forest-savanna boundary region of central Africa using multi-temporal L-band radar backscatter. *Remote Sens. Environ.* **2011**, *115*, 2861–2873. [CrossRef]

23. Carreiras, J.M.B.; Vasconcelos, M.J.; Lucas, R.M. Understanding the relationship between aboveground biomass and ALOS PALSAR data in the forests of Guinea-Bissau (West Africa). *Remote Sens. Environ.* **2012**, *121*, 426–442. [CrossRef]
24. Austin, J.M.; Mackey, B.G.; Van Niel, K.P. Estimating forest biomass using satellite radar: An exploratory study in a temperate Australian Eucalyptus forest. *For. Ecol. Manag.* **2003**, *176*, 575-583. [CrossRef]

25. Mermoz, S.; Le Toan, T.; Villard, L.; Réjou-Méchain, M.; Seifert-Granzin, J. Biomass assessment in the Cameroon savanna using ALOS PALSAR data. *Remote Sens. Environ.* **2014**, *155*, 109-119. [CrossRef]

26. Joshi, N.; Mitchard, E.; Schumacher, J.; Johannsen, V.; Saatchi, S.; Fensholt, R. L-Band SAR Backscatter Related to Forest Cover, Height and Aboveground Biomass at Multiple Spatial Scales across Denmark. *Remote Sens.* **2015**, *7*, 4442-4472. [CrossRef]

27. Wijaya, A.; Liesenberg, V.; Susanti, A.; Karyanto, O.; Verchot, L.V. Estimation of Biomass Carbon Stocks over Peat Swamp Forests using Multi-Temporal and Multi-Polratizations SAR Data. *Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.* **2015**, *XL-7/W3*, 551-556. [CrossRef]

28. Urbazaev, M.; Cremer, F.; Migliavacca, M.; Reichstein, M.; Schmullius, C.; Thiel, C. Potential of Multi-Temporal ALOS-2 PALSAR-2 ScanSAR Data for Vegetation Height Estimation in Tropical Forests of Mexico. *Remote Sens.* **2018**, *10*, 1277. [CrossRef]

29. Bouvet, A.; Mermoz, S.; Le Toan, T.; Villard, L.; Mathieu, R.; Naidoo, L.; Asner, G.P. An above-ground biomass map of African savannahs and woodlands at 25 m resolution derived from ALOS PALSAR. *Remote Sens. Environ.* **2018**, *206*, 156-173. [CrossRef]

30. Dobson, M.C.; Ulaby, F.T.; Pierce, L.E.; Sharik, T.L.; Lin, C.; Sarabandi, K. Estimation of Forest Biophysical Characteristics in Northern Michigan with SIR-C/X-SAR. *IEEE Trans. Geosci. Remote Sens.* **1995**, *33*, 19. [CrossRef]

31. Harrell, P.A.; Kasischke, E.S.; Bourgeau-Chavez, L.L.; Haney, E.M.; Christensen, N.L. Evaluation of approaches to estimating aboveground biomass in Southern pine forests using SIR-C data. *Remote Sens. Environ.* **1997**, *59*, 223-233. [CrossRef]

32. Berninger, A.; Lohberger, S.; Stängel, M.; Siegert, F. SAR-Based Estimation of Above-Ground Biomass and Its Changes in Tropical Forests of Kalimantan Using L- and C-Band. *Remote Sens.* **2018**, *10*, 831. [CrossRef]

33. Baghdadi, N.; Maire, G.L.; Bailly, J.; Osé, K.; Nouvellon, Y.; Zribi, M.; Lemos, C.; Hakamada, R. Evaluation of ALOS/PALSAR L-Band Data for the Estimation of Eucalyptus Plantations Aboveground Biomass in Brazil. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2015**, *8*, 3802-3811. [CrossRef]

34. Le Toan, T.; Beaudoin, A.; Riom, J.; Guyon, D. Relating forest biomass to SAR data. *IEEE Trans. Geosci. Remote Sens.* **1992**, *30*, 403-411. [CrossRef]

35. Le Toan, T.; Quegan, S.; Davidson, M.W.J.; Balzter, H.; Paillou, P.; Papathanassiou, K.; Plummer, S.; Rocca, F.; Saatchi, S.; Shugart, H.; et al. The BIOMASS mission: Mapping global forest biomass to better understand the terrestrial carbon cycle. *Remote Sens. Environ.* **2011**, *115*, 2850-2860. [CrossRef]

36. Freitas, S.R.; Mello, M.C.S.; Cruz, C.B.M. Relationships between forest structure and vegetation indices in Atlantic Rainforest. *For. Ecol. Manag.* **2005**, *218*, 353-362. [CrossRef]

37. Dube, T.; Mutanga, O. Evaluating the utility of the medium-spatial resolution Landsat 8 multispectral sensor in quantifying aboveground biomass in uMgeni catchment, South Africa. *ISPRS J. Photogramm. Remote Sens.* **2015**, *101*, 36-46. [CrossRef]

38. Zhu, X.; Liu, D. Improving forest aboveground biomass estimation using seasonal Landsat NDVI time-series. *ISPRS J. Photogramm. Remote Sens.* **2015**, *102*, 222-231. [CrossRef]

39. Proisy, C.; Couteron, P.; Fromard, F. Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images. *Remote Sens. Environ.* **2007**, *109*, 379-392. [CrossRef]

40. Singh, M.; Evans, D.; Friess, D.; Tan, B.; Nin, C. Mapping Above-Ground Biomass in a Tropical Forest in Cambodia Using Canopy Textures Derived from Google Earth. *Remote Sens.* **2015**, *7*, 5057-5076. [CrossRef]

41. Beguet, B.; Guyon, D.; Boukir, S.; Chehata, N. Automated retrieval of forest structure variables based on multi-scale texture analysis of VHR satellite imagery. *ISPRS J. Photogramm. Remote Sens.* **2014**, *96*, 164-178. [CrossRef]

42. Meng, S.; Pang, Y.; Zhang, Z.; Jia, W.; Li, Z. Mapping Aboveground Biomass using Texture Indices from Aerial Photos in a Temperate Forest of Northeastern China. *Remote Sens.* **2016**, *8*, 230. [CrossRef]

43. Sarker, L.R.; Nichol, J.E. Improved forest biomass estimates using ALOS AVNIR-2 texture indices. *Remote Sens. Environ.* **2011**, *115*, 968-977. [CrossRef]
44. Sarker, M.L.R.; Nichol, J.; Ahmad, B.; Busu, I.; Rahman, A.A. Potential of texture measurements of two-date dual polarization PALSAR data for the improvement of forest biomass estimation. *ISPRS J. Photogramm. Remote Sens.* **2012**, *69*, 146–166. [CrossRef]

45. Dube, T.; Mutanga, O. Investigating the robustness of the new Landsat-8 Operational Land Imager derived texture metrics in estimating plantation forest aboveground biomass in resource constrained areas. *ISPRS J. Photogramm. Remote Sens.* **2015**, *108*, 12–32. [CrossRef]

46. Attarchi, S.; Gloaguen, R. Improving the Estimation of Above Ground Biomass Using Dual Polarimetric PALSAR and ETM+ Data in the Hyrcanian Mountain Forest (Iran). *Remote Sens.* **2014**, *6*, 3693–3715. [CrossRef]

47. Sheeren, D.; Fauvel, M.; Josipović, V.; Lopes, M.; Planque, C.; Willm, J.; Dejoux, J.-F. Tree Species Classification in Temperate Forests Using Formosat-2 Satellite Image Time Series. *Remote Sens.* **2016**, *8*, 734. [CrossRef]

48. Laurin, G.V.; Balling, J.; Corona, P.; Mattioli, W.; Papale, D.; Puletti, N.; Rizzo, M.; Truckenbrodt, J.; Urban, M. Above-ground biomass prediction by Sentinel-1 multitemporal data in central Italy with integration of ALOS2 and Sentinel-2 data. *J. Appl. Remote Sens.* **2018**, *12*, 1. [CrossRef]

49. Santoro, M.; Cartus, O. Research Pathways of Forest Above-Ground Biomass Estimation Based on SAR Backscatter and Interferometric SAR Observations. *Remote Sens.* **2018**, *10*, 608. [CrossRef]

50. Castillo, J.A.A.; Apan, A.A.; Maraseni, T.N.; Salmo, S.G. Estimation and mapping of above-ground biomass of mangrove forests and their replacement land uses in the Philippines using Sentinel imagery. *ISPRS J. Photogramm. Remote Sens.* **2017**, *134*, 70–85. [CrossRef]

51. Gao, T.; Zhu, J.; Deng, S.; Zheng, X.; Zhang, J.; Shang, G.; Huang, L. Timber production assessment of a plantation forest: An integrated framework with field-based inventory, multi-source remote sensing data and forest management history. *Int. J. Appl. Earth Obs. Geoinf.* **2016**, *52*, 155–165. [CrossRef]

52. Vafaei, S.; Soosani, J.; Adeli, K.; Fadaei, H.; Naghavi, H.; Pham, T.; Tien Bui, D. 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). *Remote Sens.* **2018**, *10*, 172. [CrossRef]

53. Motte, E.; Zribi, M.; Fanise, P.; Egido, A.; Darrozes, J.; Al-Yaari, A.; Baghdadi, N.; Baup, F.; Dayau, S.; Fieuzal, R.; et al. GLORI: A GNSS-R Dual Polarization Airborne Instrument for Land Surface Monitoring. *Sensors* **2016**, *16*, 732. [CrossRef]

54. Zribi, M.; Guyon, D.; Motte, E.; Dayau, S.; Wigneron, J.P.; Baghdadi, N.; Pierdicca, N. Performance of GNSS-R GLORI data for biomass estimation over the Landes forest. *Int. J. Appl. Earth Obs. Geoinformation* **2019**, *74*, 150–158. [CrossRef]

55. Shaiek, O.; Loustau, D.; Trichet, P.; Meredieu, C.; Bachtobji, B.; Garchi, S.; EL Aouni, M.H. Generalized biomass equations for the main aboveground biomass components of maritime pine across contrasting environments. *Ann. For. Sci.* **2011**, *68*, 443–452. [CrossRef]

56. Hagolle, O.; Huc, M.; Pascual, D.V.; Dedieu, G. A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENUS, LANDSAT and SENTINEL-2 images. *Remote Sens. Environ.* **2010**, *114*, 1747–1755. [CrossRef]

57. Hagolle, O.; Huc, M.; Villa Pascual, D.; Dedieu, G. A Multi-Temporal and Multi-Spectral Method to Estimate Aerosol Optical Thickness over Land, for the Atmospheric Correction of FormoSat-2, LandSat, VENUS and Sentinel-2 Images. *Remote Sens.* **2015**, *7*, 2668–2691. [CrossRef]

58. Shimada, M.; Ohtaki, T. Generating Large-Scale High-Quality SAR Mosaic Datasets: Application to PALSAR Data for Global Monitoring. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2010**, *3*, 637–656. [CrossRef]

59. Bruniquel, J.; Lopes, A. Multi-variate optimal speckle reduction in SAR imagery. *Int. J. Remote Sens.* **1997**, *18*, 603–627. [CrossRef]

60. Quegan, S.; Yu, J.J. Filtering of multichannel SAR images. *IEEE Trans. Geosci. Remote Sens.* **2001**, *39*, 2373–2379. [CrossRef]

61. Haralick, R.M.; Shanmugam, K.; Dinstein, I. Textural Features for Image Classification. *IEEE Trans. Syst. Man Cybern.* **1973**, *3*, 610–621. [CrossRef]

62. Grizonnet, M.; Michel, J.; Poughon, V.; Inglada, J.; Savinaud, M.; Cresson, R. Orfeo Toolbox: Open source processing of remote sensing images. *Open Geospatial Data Softw. Stand.* **2017**, *2*. [CrossRef]

63. Mountrakis, G.; Im, J.; Ogole, C. Support vector machines in remote sensing: A review. *ISPRS J. Photogramm. Remote Sens.* **2011**, *66*, 247–259. [CrossRef]
64. Pal, M.; Mather, P.M. Support vector machines for classification in remote sensing. *Int. J. Remote Sens.* **2005**, *26*, 1007–1011. [CrossRef]

65. Belgiu, M.; Drăguț, L. Random forest in remote sensing: A review of applications and future directions. *ISPRS J. Photogramm. Remote Sens.* **2016**, *114*, 24–31. [CrossRef]

66. Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; et al. Scikit-learn: Machine Learning in Python. *J. Mach. Learn. Res.* **2011**, *12*, 2825–2830.

67. Inglada, J.; Vincent, A.; Arias, M.; Tardy, B.; Morin, D.; Rodes, I. Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series. *Remote Sens.* **2017**, *9*, 95. [CrossRef]

68. Vastaranta, M.; Wulder, M.A.; White, J.C.; Pekkarinen, A.; Tuominen, S.; Ginzler, C.; Kankare, V.; Holopainen, M.; Hyyppä, J.; Hyyppä, H. Airborne laser scanning and digital stereo imagery measures of forest structure: Comparative results and implications to forest mapping and inventory update. *Can. J. Remote Sens.* **2013**, *39*, 382–395. [CrossRef]

69. Vastaranta, M.; Holopainen, M.; Karjalainen, M.; Kankare, V.; Hyyppä, J.; Kaasalainen, S. TerraSAR-X Stereo Radargrammetry and Airborne Scanning LiDAR Height Metrics in Imputation of Forest Aboveground Biomass and Stem Volume. *IEEE Trans. Geosci. Remote Sens.* **2014**, *52*, 1197–1204. [CrossRef]

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

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).