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

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.

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,4,5,6,7,8], and (2) smarter use of existing data through statistical or machine learning techniques that mutually enrich unrelated data sources (multi-source approaches [9,10,11,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,14,15,16,17,18,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,21,22,23,24,25,26,27,28,29,30,31,32,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,37,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,40,41,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,44,45,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,31,32,49,50,51,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.

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∗DBHi^b

(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. Range of the forest structure parameters for the 83 sample plots selected in our study.

	Minimum	Maximum	Mean
AGB (tons/ha)	1	136	70
BA (m2/ha)	0.8	44.1	21.2
Mean tree height (m)	2	24.7	14.4
Mean DBH (m)	0.03	0.57	0.22
Density (tree/ha)	87	2622	835
Age (year)	3	72	22.6

.

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. Coefficient of determination between the forest structure parameters (plot-level).

	r2
	AGB	BA	Height	DBH	Density	Age
AGB	1	-	-	-	-	-
BA	0.90	1	-	-	-	-
Mean tree height	0.71	0.46	1	-	-	-
Mean DBH	0.49	0.24	0.83	1	-	-
Density	0.19	0.04	0.55	0.64	1	-
Age	0.40	0.17	0.74	0.90	0.53	1

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² > 0.8) for diameters under 0.15 m. There is no AGB - DBH correlation for wider diameters (r² < 0.1). The DBH is correlated with tree height (Figure 2a) until a diameter of ~0.4 m (r² = 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² > 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.

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. List of the 36 selected features as input for the regression models.

Sensors	Feature Types		Features	Abbreviation
ALOS-PALSAR-2	L-band backscatter (L-band)		L-band HV	L-HV
			L-band HH	L-HH
			L-band HV/HH	L-Ratio
Sentinel-1	C-band backscatter (C-band) DES-33°IA orbit, annual mean values		C-band VH	C-VH
			C-band VV	C-VV
			C-band VH/VV	C-Ratio
	C-band textural indexes (C-TI) DES-33°IA orbit, annual mean values		VH-homogeneity	C-VH-hom
			VH-cluster shade	C-VH-clu
			VV-homogeneity	C-VV-hom
			VV-cluster shade	C-VV-clu
			VH-correlation	C-VH-cor
			VV-correlation	C-VV-cor
			VH-Haralick correlation	C-VH-Hcor
			VV-Haralick correlation	C-VV-Hcor
Sentinel-2	Spectral indexes (S2-SI)		NDVI July	S2-NDVI-sum
			NDVI December	S2-NDVI-win
			BI July	S2-BI-sum
			BI December	S2-BI-win
			NDWI July	S2-NDWI-sum
			NDWI December	S2-NDWI-win
	S2 textural indexes (S2-BI-TI)	Summer image	BI-homogeneity	S2-BI-sum-hom
			BI-cluster shade	S2-BI-sum-clu
			BI-correlation	S2-BI-sum-cor
			BI-Haralick correlation	S2-BI-sum-Hcor
		Winter image	BI-homogeneity	S2-BI-win-hom
			BI-cluster shade	S2-BI-win-clu
			BI-correlation	S2-BI-win-cor
			BI-Haralick correlation	S2-BI-win-Hcor
Spot-6	Spot-XS textural indexes (Spot-BI-TI)		BI-homogeneity	Spot-BI-hom
			BI-cluster shade	Spot-BI-clu
			BI-correlation	Spot-BI-cor
			BI-Haralick correlation	Spot-BI-Hcor
	Spot-PAN textural indexes (Spot-PAN-TI)		PAN-homogeneity	Spot-PAN-hom
			PAN-cluster shade	Spot-PAN-clu
			PAN-correlation	Spot-PAN-cor
			PAN-Haralick correlation	Spot-PAN-Hcor

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).

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.

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,64,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.

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.

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. 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.

	L	C	C-TI	S2-SI	S2-BI-TI	r2 and RMSE
AGB	L-HV (1) L-HH (5) L-Ratio (9)		C-VH-clu (8)	S2-NDVI-sum (2) S2-BI-sum (3) S2-NDWI-sum (4) S2-BI-win (7)	S2-BI-win-hom (6)	r2 = 0.76 19.5 tons/ha (28.0%)
BA	L-HV (1) L-HH (2)		C-VH-cor (5) C-VH-Hcor (6)	S2-BI-sum (3)	S2-BI-sum-hom (4) S2-BI-sum-Hcor (7)	r2 = 0.73 5.7 m2/ha (26.9%)
DBH	L-HV (2)	C-Ratio (6)	C-VV-cor (1) C-VH-clu (4) C-VH-cor (5)		S2-BI-sum-cor (3)	r2 = 0.88 0.04 m (19.8%)
Age	L-HV (2) L-HH (5)		C-VV-hom (1) C-VH-Hcor (3) C-VV-cor (4) C-VH-clu (8)	S2-BI-win (9) S2-BI-sum (10)	S2-BI-sum-hom (6) S2-BI-sum-Hcor (7) S2-BI-win-cor (11)	r2 = 0.93 3.95 years (17.4%)
Density	L-HV (2)	C-Ratio (5)	C-VV-Hcor (1) C-VH-hom (4) C-VH-Hcor (6)	S2-BI-sum (7)	S2-BI-sum-hom (3)	r2 = 0.86 204 trees/ha (24.4%)
Height	L-HV (4) L-HH (10)	C-VH (9)	C-VH-Hcor (1) C-VV-clu (8) C-VH-clu (11) C-VV-cor (12) C-VH-cor (13) C-VV-hom (14)	S2-NDWI-win (2)	S2-BI-sum-cor (3) S2-BI-sum-hom (5) S2-BI-win-hom (6) S2-BI-sum-Hcor (7)	r2 = 0.93 1.75 m (13.2%)

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.

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).

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.

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).

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.

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.

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. 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).

	AGB	BA	DBH	Age	Density	Height
Best feature type	L-band 32%	L-band 29%	C-TI 26%	C-TI 27%	C-band 37%	C-TI 18%
Best two types combination	L + S2-BI-TI	L + C	C-TI + S2-BI-TI	L + C-TI	L + S2-BI-TI	C-TI + S2-BI-TI
	28.0%	27.1%	21.2%	18.6%	30.7%	14.3%
Five types open access images	28.0%	26.9%	19.8%	17.4%	24.4%	13.2%
Seven types (Spot-6 included)	27.9%	27.0%	19.3%	18.6%	23.3%	13.7%

present the relative RMSE (%) obtained for the estimation of the different forest parameters according to the four types conditions.

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,41,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 forests 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.