Forest Aboveground Biomass Estimation Using Machine Learning Ensembles: Active Learning Strategies for Model Transfer and Field Sampling Reduction

Abstract

Biomass is a crucial indicator of the carbon sequestration capacity of a vegetation ecosystem. Its dynamic is of interest because it impacts on the carbon cycle, which plays an important role in the global climate and its changes. This work presents a novel technique, able to transfer a calibrated regression model between different areas by exploiting an active learning methodology and using Shannon’s entropy as a discriminator for sample selection. Model calibration is performed based on a reference area for which an extended ground truth is available and implemented via regression bootstrap. Then, re-calibration samples for model transfer are selected through active learning, allowing for choosing a limited number of points to be investigated for training data collection. Different sampling strategies and regression techniques have been tested to demonstrate that a significant reduction in the number of calibration samples does not affect the estimation performance. The proposed workflow has been tested on a dataset concerning Finnish forests. Experimental results show that the joint exploitation of regression ensembles and active learning dramatically reduces the amount of field sampling, providing aboveground biomass estimates comparable to those obtained using literature techniques, which need extended training sets to build reliable predictions.

1. Introduction

Forests cover about 40% of the Earth’s ice-free land surface, although, in the recent past, wide forested areas have been converted to agricultural or urban land uses. They constitute the vital core of the whole planet, because they are the sentinels of the carbon cycle and the custodians of its biodiversity [1]. The expected increase in the global population is a further threat, representing a potential trigger for increasing deforestation and overexploitation of resources. The challenge to be faced by both politicians and scientists is the definition of protocols and strategies allowing for the sustainable management of such critical areas [2]. In this context, it is important to support decision-makers with timely information about several biophysical parameters useful to implement adequate planning activities [3]. Among them, the aboveground biomass (AGB) is particularly interesting due to its impact on the carbon cycle and, consequently, on the global climate and its changes [4].

Traditionally, this parameter is estimated through field sampling [5] or the exploitation of remote sensing data [6]. The results of the first methodology are unfeasible for applications at forest scale. Conversely, remote sensing technologies are able to deliver synoptic pictures of large areas with high revisit time. Therefore, as widely demonstrated by the literature, they represent a unique data source to feed methodologies, allowing for AGB estimation at large scale [7,8,9].

To this end, several techniques have been proposed in the literature. Regression-based models are among the most exploited [10,11,12]. They have many advantages, like the possibility to combine heterogeneous data, the low computational complexity and the suitability to deal with sparse sampling, which is the typical case in forest studies [5]. Machine learning methodologies have been exploited to better deal with nonlinearities between measurements and remote sensing data. Previous studies demonstrated the effectiveness of solutions like random forests [13], gradient boost algorithms [14] and support vector machines [15]. However, some concerns have been raised against their parametrization, as tuning can significantly affect the performance of these models [16]. Moreover, most of the studies rely on sparse sampling over large areas for model calibration and validation, thus introducing issues related to generalization [17].

More recently, deep learning solutions are emerging to address AGB estimation problems, as proposed, as an example, in [18,19,20]. In this case, major constraints are related to the sparse sampling typically implemented in forestry applications, as this methodology is usually unsuitable for training neural networks [21]. In the past, SAR scattering models have demonstrated their reliability for this purpose [22,23,24]. However, these solutions require careful parametrization, including the specification of various stand properties to be retrieved through field campaigns [25].

Working at forest scale, the main issue related to the exploitation of algorithms requiring the acquisition of calibration points is probably generalization, as sampling is usually implemented in specific areas [17]. Therefore, a concept of interest is the introduction of methodologies allowing for the transfer of a model from one area to another, sharing similar characteristics [26]. This problem is particularly of interest, as an example, in precision agriculture, since the variations of the AGB are typically higher due to crop cycles and their sensitivity to environmental conditions. In this context, Reference [26] exploited active sampling to select significant samples to transfer a calibrated model for ryegrass AGB and nitrogen content estimation across different seasons. Wan et al. exploited the Kennard–Stone algorithm to select re-calibration samples acquired in different years, to extend the forecasting capabilities of a rice yield model [27]. Reference [28] proposed an active sampling framework to implement multi-year incremental monitoring of grassland biomass and nitrogen content.

Dealing with forests, the problem of the update of an existing model to make it work when applied to other areas, or to the same area in different times, has been scarcely addressed, as the dominant approach is the validation of new techniques fed by variously collected data [4,7,10,13,29]. The most relevant literature on the topic concerns the processing of airborne laser scanning (ALS) data. Reference [30] exploited repeated ALS acquisitions and ground measures to update models predicting aboveground carbon, basal area, stand density index and total stem volume within a study area in northern Idaho (USA). Tompalski et al. demonstrated the feasibility of transferring between different areas of forest attributes derived using ALS data [31]. Reference [32] proposed a methodology to generate new maps of forest structural attributes, integrating ALS samples within pre-fitted models, pointing out that the direct transfer of an available model elsewhere may involve extrapolations or its application to conditions not included in the training dataset.

In this regard, the purpose of this paper is to introduce a strategy to tackle the typical generalization issues found in the past AGB estimation literature. In particular, active learning principles [33] are exploited to build a novel workflow suitable for transferring a calibrated regression model to new areas, through the selection of relevant re-calibration samples. To this end, a novel technique for the selecting of the most informative samples, based on their entropy rank, is introduced. Different from commonly adopted solutions, in which training data are distributed over large areas and used to train a single model, the principal novelty here is the exploitation of several regression models, generated from a reference one, retrained with few samples collected on specific areas. The aim is to extend the forecasting capability of an existing training dataset through the addition of a limited number of samples, made available by new acquisitions, through the implementation of active learning [26,34]. Finally, particular attention is paid to validation. This is one of the weaknesses found in the literature, which tends to assess the performance of the estimation over very large areas using a few tens of plots. In this work, the obtained results will be evaluated against an extended ground truth, accounting for more than 850,000 validation points, derived from the LiDAR processing of data acquired over Finnish forests.

The work is organized as follows. Materials and methods are introduced in Section 2. Experimental results are presented in Section 3 and discussed in Section 4. Conclusions are drawn at the end of the work.

2. Materials and Methods

2.1. Data

In Finland, ALS data are systematically collected for digital terrain mapping and forest parameters estimation [35], including AGB. These data, opportunely integrated with satellite images and historical field sampling, have been exploited by the Finnish Forest Centre to produce the reference AGB data used in this study [36], which analyzes an area of about 85 km², partitioned in 13 clips of 6.55 km² extension each. The clips, whose global statistics are reported in

Table 1. AGB statistics for the area considered in this study.

Min AGB Clip (t/ha)	Max AGB Clip (t/ha)	Mean AGB (t/ha)	Std AGB (t/ha)	Total Area (km2)
22.1	81.3	69.5	43.1	85.1

, are quite dishomogeneous in terms of the average AGB amount and its distribution (see detailed data in Table S1, provided as supplementary material). The average AGB ranges between 37.1 t/ha and 81.3 t/ha. Its standard deviation ranges between 22.2 t/ha and 63.2 t/ha. The overall AGB mean and standard deviation are in the order of 69.5 t/ha and 43.1 t/ha. The area is characterized by almost flat topography, with elevation values ranging approximately 100–200 m above sea level, with a lack of significant urban areas and agricultural lands. It belongs predominantly to the middle boreal vegetation zone, with a dominant presence of coniferous [37].

Satellite data, used to predict the AGB, belong to a publicly available dataset (see the data availability statement at the end of the work) composed of several clips. It is presented as a benchmarking tool for AGB studies. Therefore, data are provided without any geographic and/or temporal attributes. For each clip, the dataset includes, beyond AGB data, one year of Sentinel-1 SAR and Sentinel-2 MS observations, made approximately on a monthly basis (see Figure 1). In particular, SAR data have been acquired regularly once per month, due to their insensitivity to atmospheric conditions. MS acquisitions are mostly concentrated during spring and summer in order to cope with massive cloud coverage typically registered in the other periods of the year.

2.2. Methodology

As stated in the Introduction, the purpose of the paper is to introduce a methodology allowing for the transfer of a calibrated partial least squares regression (PLSR) [38] model to different areas sharing similar characteristics. In other words, it is assumed that a model, calibrated using an extended ground truth, is available. The question to be addressed is about how to select relevant samples for fitting this model with new scenes.

To this end, the workflow depicted in Figure 2 is adopted. The inputs are a calibrated regression model and two time series of data, one MS and one SAR, which are exploited to calculate a set of predictors, a total of 159, to be exploited for regression. MS predictors are mostly constituted by vegetation indices. SAR ones are constituted by both intensity-based and texture variables derived from the Haralick co-occurrence matrix [39].

The predictors are calculated for each available acquisition and then averaged on a yearly basis. Following this operation, the processing forks in two branches. The first one is dedicated to the retrieval of the model to be transferred, which is obtained via bootstrap, i.e., a regression calibrated with a random sample selection. Bootstrapping is implemented as follows. It is assumed that the AGB is known in a specific area (Clip 0 in Table S1). Then, among 1000 PLSR experiments calibrated with random sample selection, the one providing the lowest RMSE is selected as reference for transferring.

Model transferring requires that the same predictors are used for each regression. However, the literature highlighted that the large amount of data usually ingested into PLSR can contain irrelevant information that can cause the reduction in the model performance [40]. Therefore, a set of potentially informative predictors is selected through the calculation of the variable importance in prediction (VIP) score [41].

As the VIP score is an output of PLSR, a tentative regression is implemented with this purpose. In other words, this will not be used for AGB estimation. In this phase, the number of latent variables, i.e., the set of orthogonal factors having the best predictive power, is set to 10.

The VIP score is defined, for each variable j, as the sum, over latent variables f, of its PLS-weight value $w_j$ weighted by the percentage of explained variance of the specific latent variable $SS{Y}_f$. For the j-th variable, it holds [41]:

$$VI{P}_j=\sqrt{\frac{J\cdot {\sum }_{f=1}^F{w}_{jf}^2\cdot SS{Y}_f}{F\cdot {\sum }_{f=1 }^{F}SS{Y}_f}},$$

(1)

where F is the total number of latent variables and J the total number of variables. The relation for the calculation of $SS{Y}_f$ is given by [41]:

$$SS{Y}_f={\mathit{b}}_f^2{\mathit{t}}_f^T{\mathit{t}}_f,$$

(2)

where ${\mathit{b}}_f$ and ${\mathit{t}}_f$ are the PLS inner relation vector of coefficients and the score vector relevant with the f-th latent variable, respectively. It is common practice in the literature to assume that informative variables are associated with a VIP score larger than one. This means that those variables have an above-average influence on the building of the model explaining the observations [42].

According to Figure 2, following the definition of the model to be transferred, active sampling is implemented to select re-calibration data. To this end, the stack composed of all the predictors for the new area is processed with a principal component analysis (PCA) for dimensionality reduction. Then, an ISODATA clustering is implemented to obtain labelled regions. In this case, the maximum number of allowed classes is set to 10.

The classified image provided by ISODATA clustering is treated with connected component labelling [43] to retrieve elementary regions, i.e., the image segments. They are used to calculate the spatial average of the regression variables, which constitute the input for active learning.

As described in [34], active learning techniques can be categorized based on uncertainty [44] or diversity criteria [45]. In the first case, samples are ranked according to their uncertainty. The higher the uncertainty, the better their rank. Among these techniques, those based on variance-based pool of regressors (PAL) are probably the most interesting [46]. Active learning techniques based on diversity criteria, instead, select samples based on the dissimilarities they introduce in the training datasets [34]. Different metrics can be used to assess such dissimilarity like the Euclidean distance [47], the cosine angle distance [48] or the entropy [49], which is the one used in this work according to the definition provided by Shannon [50]:

$$H=-{\displaystyle \sum }_{i=1}^N{P}_n{\log }_2{P}_n,$$

(3)

where $P_n$ is the normalized probability of the n-th histogram quantization level and N is the total number of bins. The entropy gives information about the shape of the histogram and measures the quantity of information carried by a signal. The higher the entropy, i.e., the flatter the histogram, the higher its informative content.

According to Reference [28], the histogram is calculated on the vector constituted by the region-wise average of the predictors. In other words, given the N predictors for the generic orange region provided by the ISODATA clustering depicted in Figure 3, the reference vector for histogram calculation is composed of the N average values of the predictors. The entropy $H_1$ is calculated according to Equation (3) and stored as metrics for active sampling selection.

The second step of the active sampling process, the vector representative of the green region of Figure 3, is appended to the one previously considered. The entropy H₂ of this new vector, constituted by 2N elements is calculated. If the condition $H_2>H_1$ is verified, it is retained. Otherwise, the vector corresponding to the green region is discarded. In any case, the entropy value used for the acceptance test is the last one computed. The process continues until all the regions have been tested. Those determining an increase in the entropy against the last available one are appended to the selection and the corresponding regions marked as informative.

In summary, the adopted active sampling procedure can be described as follows:

a.

Consider each retrieved cluster and make the spatial average for each spectral band;

b.

Starting from a randomly selected cluster, calculate its histogram considering all the elements of the average spectral response and its associated entropy $H_1$ according to Equation (3);

c.

Add another cluster to the dataset by appending its (average) spectral response to the vector constituted by the one previously considered. Calculate the new histogram and its entropy $H_2$;

d.

If $H_2>H_1$ mark the cluster as informative and continue the process by adding new clusters to the dataset. Do not delete those marked as not informative. Clusters marked as informative are those to be sampled to retrieve model calibration data.

The average number of samples, i.e., regions, selected by the adopted active sampling methodology is around 301 within areas of about 6.5 km². In an operational environment, this number can be considered high. Therefore, a simple strategy for reducing the number of samples to be considered for model transfer has been implemented.

The output of the above-discussed procedure is the list of the regions marked as informative. Each of them is associated with an entropy value, calculated by adding the vector representative of the i-th region to the one collecting those already accepted. Moreover, a class ID coming from the ISODATA clustering is available.

Informative regions are then separated, based on the ISODATA class ID, and sorted in descending entropy order. The first k samples of the chart are retained for model calibration. In such a way, the maximum number of calibration samples is fixed to $n_{max}=kt$, where t is the maximum number of allowed classes in ISODATA clustering. In the following experiments, t has been set to 10. Two different settings for k have been tested, i.e., k = 10, yielding $n_{max}=100$, and k = 5, yielding $n_{max}=50$. In the following, these samples will be referred to as very informative samples (VIS). The adoption of this strategy allowed for the reduction in the areas to be sampled to 73 per clip, in the case of k = 10, and to 36 per clip, in the case of k = 5.

Regressions for test areas are implemented using different strategies. The first one is a PLSR. It is a statistical method, aiming at finding a linear regression model between the observations and the independent variables, by projecting the predicted variables and the observed ones to a new space [38]. This technique is particularly suited when the matrix of predictors has more variables than observations and there is multicollinearity among predictors, as in the case study at hand.

The second regression technique used for AGB predictions Is gradient boosting (GB). It consists of a finite set of weak learners (typically decision trees) that minimize an arbitrary differentiable loss function. A meta-learner assigns weights to each learner and combines their predictive results through voting methods to provide a better predictive performance for regression problems. In particular, the least-squares boosting algorithm has been used [51]. This means that the loss function to be minimized is the mean squared error of the learners. The ensemble of GB and PLSR represents the third regression strategy. In this case, the final result is given by the average of the two predictions.

3. Results

The aggregated results of the proposed AGB estimation methodology for the 12 considered clips are reported in

Table 2. Aggregated AGB estimation results for the different tested sampling settings and benchmarking techniques. The column “Inc” refers to the incremental sampling setting.

Proposed							Benchmark
	PLSR		GB		Ensemble		Bootstrap		NN
Samples	RMSE		RMSE		RMSE		RMSE*	$\overline{RMSE}$
	Area	Inc	Area	Inc	Area	Inc
All	28.8	31.3	33.7	36.3	28.8	30.6	26.2	46.8	30.4
k = 10	30.7	32.1	34.2	36.7	30.0	31.5	27.7	47.7
k = 5	31.6	32.2	34.0	37.9	30.7	32.3	28.5	49.7

RMSE*: minimum bootstrap RMSE; $\overline{\mathrm{RMSE}}$: average bootstrap RMSE.

, together with benchmarking results. The breakdown for each single clip is reported in Table S2 as supplementary material.

Using all the samples selected by the adopted active sampling technique (see Table S2), the obtained RMSE, using PLSR, ranges between 20.8 t/ha and 47.2 t/ha with an average of 28.8 t/ha. GB provided a minimum RMSE of 24.6 t/ha and a maximum one of 48.5 t/ha with an average of 33.7 t/ha. The ensemble resulted in RMSE values ranging between 21.4 t/ha and 44.7 t/ha with an average of 28.8 t/ha.

Reducing the samples to VIS does not significantly affect the estimate (see also Table S2). In the case of k = 10, the RMSE ranges between 21.5 t/ha and 51.6 t/ha with an average of 30.7 t/ha, using PLSR. GB estimates provided RMSE values ranging between 25.1 t/ha and 52.1 t/ha, with an average of 34.2 t/ha. The estimates obtained using the ensemble range between 21.4 t/ha and 46.2 t/ha, with an average of 30.0 t/ha.

When the setting is k = 5 (see Table S2), the obtained RMSE ranges between 21.7 t/ha and 50.3 t/ha, with an average of 31.6 t/ha, using PLSR. GB exploitation provided RMSE values ranging between 24.8 t/ha and 51.4 t/ha with an average of 34.0 t/ha. Using the ensemble, the RMSE ranges between 21.4 t/ha and 49.3 t/ha with an average 30.7 t/ha.

The column named “Inc” reports the results obtained using an incremental model. In this configuration, all the points used to calibrate clips previously analyzed are included within the calibration set for the current regression. In other words, as an example, the calibration set for the third regression includes the data used to train the models used for the first and the second regression. Indeed, the exploitation of more calibration points, referring to different areas, seems to not be beneficial for the estimate, as the average RMSE is higher than the one retrieved using the standard setting in most of the experiments (see Table S2).

Some scatterplots relevant to selected experiments (see Table S2 for clip indexing) are depicted in Figure 4. In particular, the first row of the graphics reports the results of AGB estimation for Clip 0 (Figure 4a) and Clip 4, relative to ensemble regression with the setting k = 5 (Figure 4b) and the neural network output (Figure 4c). The regressions shown in the second row of the graphics concern Clip 7. The one reported in Figure 4d has been obtained using ensemble regression calibrated with all the samples selected by active learning. The plot shown in Figure 4e concerns ensemble regression with k = 5. All the other plots are omitted for brevity, but aggregated results in the plane (R², RMSE) are shown in Figure 4f, relative to ensemble regression.

The variables used for regression, determined through VIP score thresholding, are all derived from the MS dataset. In particular, the 17 indices reported in

Table 3. VIP variables extracted from the regression of the reference area.

Name	Formula	Ref
Band 2	${\mathrm{R}}_{492}$
Band 5	${\mathrm{R}}_{704}$
Band 8	${\mathrm{R}}_{832}$
Aerosol free vegetation index	${\mathrm{AFRI}}_{1600}={\mathrm{R}}_{830}-0.66\frac{{\mathrm{R}}_{1613}}{{\mathrm{R}}_{832}+0.66{\mathrm{R}}_{1613}}$	[52]
Aerosol free vegetation index	${\mathrm{AFRI}}_{2100}={\mathrm{R}}_{832}-0.5\frac{{\mathrm{R}}_{2202}}{{\mathrm{R}}_{832}+0.5{\mathrm{R}}_{2202}}$	[52]
Ashburn vegetation index	$\mathrm{AVI}=2{\mathrm{R}}_{832}-{\mathrm{R}}_{560}$	[53]
Chlorophyll absorption ratio index	$\mathrm{CARI}=\frac{{\mathrm{R}}_{704}}{{\mathrm{R}}_{664}}\sqrt{\frac{{\left(670\mathrm{a}+{\mathrm{R}}_{664}+\mathrm{b}\right)}^2}{{\left({\mathrm{a}}^2+1\right)}^{0.5}}}$	[54]
Difference 800/550	${\mathrm{D}}_{800/550 }={\mathrm{R}}_{832}-{\mathrm{R}}_{560}$	[55]
Green difference vegetation index	$\mathrm{GDVI}={\mathrm{R}}_{864}-{\mathrm{R}}_{560}$	[56]
Differenced vegetation index MSS	$\mathrm{DVIMSS}=2.4{\mathrm{R}}_{832}-{\mathrm{R}}_{664}$	[53]
Global environment monitoring index	$\mathrm{GEMI}=\mathrm{n}\left(1-0.25\mathrm{n}\right)-\frac{{\mathrm{R}}_{664}-0.125}{1-{\mathrm{R}}_{664}}$	[53]
Misra soil brightness index	$\mathrm{MSBI}=0.406{\mathrm{R}}_{492}+0.600{\mathrm{R}}_{560}+0.645{\mathrm{R}}_{704}+0.243{\mathrm{R}}_{832}$	[53]
Misra yellow vegetation index	$\mathrm{MYVI}=0.723{\mathrm{R}}_{492}+0.597{\mathrm{R}}_{560}+0.206{\mathrm{R}}_{704}+0.278{\mathrm{R}}_{832}$	[53]
Modified chlorophyll absorption in reflectance index	$\mathrm{MCARI}=\left[\left({\mathrm{R}}_{704}-{\mathrm{R}}_{832}\right)-0.2\left({\mathrm{R}}_{704}-{\mathrm{R}}_{560}\right)\right]\frac{{\mathrm{R}}_{704}}{{\mathrm{R}}_{664}}$	[57]
Nonlinear vegetation index	$\mathrm{NLI}={\mathrm{R}}_{782}^2-{\mathrm{R}}_{664}$	[58]
Reflectance at the inflexion point	$\mathrm{RRE}=\frac{{\mathrm{R}}_{664}+{\mathrm{R}}_{782}}{2}$	[59]
Simple ratio 833/1649	${\mathrm{SR}}_{800/1649}=\frac{{\mathrm{R}}_{832}}{{\mathrm{R}}_{1613}}$	[60]

were informative for model training. Remarkably, none of them are derived from the SAR dataset.

The adopted benchmarking techniques are randomly calibrated PLSR and NN-based estimation. In the first case, a number of samples equal to the one determined via active sampling has been adopted for a fair comparison. Bootstrapping has been implemented to retrieve the minimum RMSE and the average one, for a total of 1000 experiments.

Bootstrapping returned a minimum RMSE (indicated as RMSE* in Table 2) ranging between 21.0 t/ha and 39.4 t/ha (see Table S2). The average of the minimum RMSE is 26.2 t/ha. The column reporting the mean RMSE values ($\overline{\mathrm{RMSE}}$ in Table 2), obtained via bootstrap, is useful to aid in pointing out that this method can fail, as not all the experiments can be considered successful whatever the combination of samples is used for model calibration (see Table S2).

NN-based estimation has been implemented using ReUse [17]. This is a UNet network trained to perform a pixel-wise regression task, mapping Sentinel-1 and Sentinel-2 images, concatenated along the spectral axis into a single tensor, into AGB raster. The architecture consists of two pathways. The first is a contracting pathway (the encoder) that captures the image context. The second pathway is a symmetric expanding pathway (the decoder) designed to produce pixel-wise predictions using transposed convolutions. Several skip connections connect the two pathways.

It is worthwhile to remark that the deep learning benchmark was carried out without the usage of training samples belonging to the clips under investigation for AGB estimation. The training phase has been implemented using the same number of clips exploited for validation. As for the bootstrap benchmark, which is typically based on sparse sampling, the training on points belonging to the clips under testing allows the model to learn the characteristics of the specific clip effectively. This approach is unsuitable for deep learning, which is typically trained using patches, here retrieved from clips not considered for AGB estimates.

The network generates a pixel-wise regression map and can extract spatial and spectral features from satellite images using an end-to-end paradigm. In particular, a patch-wise approach is adopted [57]. Each input image and ground truth AGB map is divided into nonoverlapping 32 × 32 pixel patches. The training has been implemented using a ground truth with comparable dimensions of the area to be estimated. It stops when the monitored validation loss stops to improve after 35 epochs. The maximum number of epochs is set to 200. The optimizer used is Adam [58] with default parameters. The learning rate is reduced by a factor of 0.2 if no improvement in validation loss is observed for 25 epochs. The adopted loss function is the mean absolute error. The RMSE values obtained using this technique range between 20.6 t/ha and 43.9 t/ha (see Table S2), with an average of 30.4 t/ha.

4. Discussion

Biomass is a crucial indicator of the carbon sequestration capacity of a vegetation ecosystem [61]. Its dynamic is of interest as it impacts the carbon cycle, which plays an important role in the global climate and its changes. This work introduces a novel framework for AGB estimation based on the transferring of a calibrated regression model between different scenes through active selection re-calibration samples. In this regard, a simple, fully unsupervised technique, exploiting Shannon’s entropy as a discriminator, has been adopted, to select a number of calibration samples compatible with operational environments.

The obtained results testified the reliability of the workflow, which restituted an average RMSE of 28.8 t/ha, in case all the samples selected by the adopted active sampling methodology, on average 46/km², are exploited. The adoption of the ISODATA-based selection allowed for the reduction in the samples to be collected to about 10/km² (when the setting was k = 10) and 6/km² (when the setting was k = 5) with negligible impact on the AGB estimation, whose RMSE increased to 30.0 t/ha when the setting was k = 10 and to 30.7 t/ha when the setting was k = 5. The registered values of the coefficient of determination, R², as expected, follow the trend of the RMSE, i.e., they tend to increase as the RMSE decreases, as shown in Figure 4f. As a general comment, it is remarkable that the incremental framework, in which all the samples collected for previous estimates are involved in the current regression, did not add value to the estimates.

The aforementioned RMSE values refer to the most performing regression configuration which, interestingly, is provided by the ensemble of PLSR and GB. This configuration outperforms the single regression technique in 18 experiments out of 32. This trend increases as the number of calibration points is reduced. The analysis of the errors reveals that they are strongly correlated with the AGB standard deviation of the clip. The higher the standard deviation, the higher the RMSE. This result is in line with the literature, as reported, as an example, in [62].

As shown in the previous section, the variables selected for regression through VIP scoring belong to the MS dataset. Despite the well-known relation between SAR backscattering and the AGB [8], spectral information in the near-infrared and red edge frequencies, which appear in almost all the indices listed in Table 3, resulted dominant. This is due to the scarce sensitivity of C-band frequencies to AGB within the resolution cell [63]. However, the contribution of SAR features is expected to be more significant using L-band data, as this wavelength typically shows higher correlation with the AGB [8]. This result can be found, as an example, in [62], where the authors highlighted how AGB estimates can benefit from the presence of L-band SAR features among regression variables.

The RMSE values obtained adopting the proposed methodology are comparable with those declared in the literature for the specific application. As an example, Reference [7] claimed an RMSE of 18.9 t/ha on an area having an average AGB of 49.4 t/ha, with a standard deviation of 33.8 t/ha using extreme gradient boost algorithm. Reference [14] reported an RMSE of about 30.0 t/ha in the study of an area with an average AGB of 55.5 t/ha and a standard deviation of 18.2 t/ha, using extreme gradient boosting as well. Vafaei et al. [62] reported an RMSE of 38.7 t/ha, using support vector regression in an area with an average AGB of 206.5 t/ha and a standard deviation of 78.0 t/ha.

All these studies were characterized by the same approach, consisting of sparse sampling of a few plots (149 in Reference [62] and 367 in Reference [7]) distributed over large areas. These data are used in part to train the model and in part to validate it. In other words, the authors used a few samples to calibrate a single model to fit a large area. The validation was performed on a small subset of the ground truth, as the training phase tends to absorb the majority of the collected samples.

The approach proposed here is different. First, the AGB estimation problem is tackled via active learning, which allows for the selection of highly informative calibration samples, thus introducing a further guide for field operations which can be accordingly significantly reduced. It requires extended reference data about a small area, which can be retrieved with a single LiDAR flight, useful to calibrate a machine learning model via bootstrapping, which, as demonstrated in the previous section, is the most effective solution in the presence of extended ground data. This model can be effectively transferred elsewhere through the collection of a few new calibration samples, according to active learning principles. This means that the objective was to demonstrate that a reliable selection of re-calibration samples feeding simple regression models can provide results comparable with methodologies exploiting more sophisticated solutions tuned on the specific case study and/or requiring extended ground data.

In this regard, the proposed active learning technique dramatically reduces the necessity of field sampling, which is the principal limitation in large-scale AGB mapping [5], up to about 6 samples/km². The concept of active sampling for model transfer also opens the possibility of the exploitation of historical field data in combination with newly acquired ones with the purpose of updating and/or expanding past AGB maps [37]. Differently from the literature, using, as aforementioned, a few plots to retrieve estimation quality metrics, the validation has been performed on an extended ground truth for a total of more than 850,000 validation points.

The reliability of the proposed methodology has been tested through comparison against some of the most popular techniques for AGB estimation, like bootstrap and NN-based estimation. The benchmark revealed that random sampling is not a viable option if an extended ground truth is not available. Using this methodology, in most cases, the average RMSEs are much higher than the typical ones indicated in the literature. On the other hand, in the presence of extended ground data, it can be considered as an option, as the provided minimum RMSEs constitute, in five cases out of seven, the best overall estimation performance (see Table S2).

NNs are popular tools for AGB estimation allowing for reliable estimates, provided that an extended ground truth is available for model training [18,64]. Dealing with AGB estimates, this represents a significant limitation, since reference data are usually retrieved via sparse plot sampling. The obtained RMSE values, using ReUse, testify that NNs are suitable for model transfer, i.e., it is possible to train the net using reference data acquired on a test area and then use the weights for making inference elsewhere. However, the amount of reference data should be comparable with the extension of the area to be predicted.

The performed experiments revealed that the average RMSE values obtained using ReUse are comparable with those provided by the proposed framework and by the best bootstrap performance. This means that the choice of the best solution for the specific AGB estimation problem to be solved can be a function of the available data and computational power. In other words, the usage of machine learning ensembles can be a viable solution to cope with the sparsity of reference data. On the other hand, when extended ground truths can be collected, NNs and bootstrapping, or even a combination of the two, could represent an option, as they can slightly reduce the estimation error.

5. Conclusions

Forests cover about 40% of the Earth’s ice-free land surface and constitute the vital core of the whole planet since they represent the sentinels of the carbon cycle. Monitoring their biomass is key, as this parameter is the principal indicator of the carbon sequestration capacity of a vegetation ecosystem.

Remote sensing literature about forest aboveground biomass estimation is principally oriented towards machine learning, due to its ability to deal with sparse samples acquired during field campaigns. This also allows researchers to cope with large-scale problems, as samples are tendentially distributed all over the area to be investigated.

In this work, the problem of aboveground estimation is tackled using active learning principles, allowing for the transfer of a calibrated regression model through a data-driven selection of the most informative re-calibration samples. In this regard, a technique exploiting Shannon’s entropy as a diversity criterion has been discussed, with a particular focus on how it allows for the dramatic reduction in the number of samples to be collected, to a few samples per square kilometer, which is a quantity suitable for operational environments. The proposed framework requires one reference area, for which an extended ground truth should be available, to calibrate the model to be transferred to the rest of the forest, assuming that its structure does not exhibit significant variations. Reference AGB data can be obtained, as an example, with a LiDAR acquisition.

The estimates provided by the proposed methodology have been compared against techniques set out in the literature that exploit massive reference data for model training, including bootstrapping and neural networks. This investigation revealed that the exploitation of machine learning ensembles, in combination with active learning, provides results, measured in terms of the root mean square errors, fully comparable with the literature. This suggests that the proposed methodology is a viable solution to cope with the sparsity of reference data. On the other hand, in case an extended ground truth is available, both neural networks and bootstrapping represent feasible solutions as they allow for a slight reduction in the estimation error.