1. Introduction
Forest and its ecosystems play a vital role in the overall health and sustainability of the planet earth, providing indispensable services such as carbon sequestration, biodiversity conservation, and water regulation [
1,
2]. These ecosystems support a vast array of life forming processes for both plant and animal species, contributing to ecological functionality, stability, and resilience [
3,
4]. In addition, they offer economic benefits through the provision of resources such as timber and non-timber forest products (NTFPs), highlighting their multifaceted importance [
5]. Forest growing stock volume (GSV) is a critical metric in forestry resources assessments (FRA), representing the total volume of all living trees within a given area [
6,
7], serving as an indicator for measuring the forest management capability, forest quality, and carbon sequestration potential [
8]. The accurate estimation of this essential variable is pivotal for sustainable forest management, conservation practices, quantification of ecosystem services, and aid in providing detailed information to support decision making mechanism [
9,
10]. Advances in remote sensing data capturing technology and machine learning (ML) algorithms have significantly enhanced the precision and efficiency of GSV assessments, casing a fundamental change in estimation thereby providing valuable insights for resource management [
11,
12].
Traditional and Remote Sensing (RS) based estimation methods in assessing GSV offer distinct advantages and present challenges, shaping their complementary roles in forest inventory and management. Traditional methods, such as field sampling and harvesting data, provide accurate measurements at localized scales by directly observing tree dimensions and timber volume [
6]. However, the fieldwork for acquiring data for forest growing stock volume at a larger scale is challenging [
13] and holistically, these approaches are labour-intensive, costly, and may lack continuous coverage over large, forested area [
14]. Conversely, the RS method offers significant advancements by providing spatially explicit and continuous assessments of forest related variables. With a long history of remote sensing data availability, different data sources have been utilized to estimate forest variables out of which optical satellite imagery such as Landsat has been widely used [
15]. The synoptic coverage, moderate ground cover resolution, high repetitive coverage, and free access policy makes some optical RS the most widely used [
15,
16,
17] meanwhile, While RS methods enhance efficiency and coverage, it is expedient to note that they cannot emphatically replace field sampling methods. They require careful calibration and validation with ground-truth data from traditional methods to ensure accuracy and precision [
18]. Integration of both approaches through hybrid methodologies facilitates comprehensive GSV assessments [
19], leveraging the strengths of RS for broad-scale monitoring and the precision of traditional methods for validation and localized accuracy.
Optical sensors usually provide detailed information on the horizontal structure of forest vegetation based on the configuration, meanwhile, information about the vertical structure is not presented. On the contrary, LiDAR remote sensing can provide both horizontal and vertical information, depending on the configuration of the LiDAR system used and therefore, it is well suited for describing forest structure [
20,
21]. The use of optical imagery with high resolution for estimating GSV has proven to be effective for example, Mura et al. [
22] explored the capabilities of the sentinel-2 for predicting growing stock volume in Italy. Zhou & Feng [
8] compared the precision of sentinel-2 and Landsat 8 OLI in estimating forest stock volume and concluded that sentinel-2 based model achieved higher accuracy. Hence the choice of integrating LiDAR metric and Sentinel-2 data in our current study. Research efforts geared towards the usage of Airborne laser scanning (ALS) in the estimation of GSV has proven its resourcefulness as a valuable data source [
13,
23,
24,
25,
26]. Lindgren et al [
27] emphasized that in addition to ALS, there are many other sources of remote sensing data that could be used for automated prediction of forest variables like GSV However, lack of good prediction accuracy has reduced their use in operational forest management.
In providing more accurate, non-destructive frequent updates, remotely sensed data can improve the quality of forest inventory database, thus, the synergy between traditional and RS methods in GSV estimation not only enhances forest management practices but also supports sustainable resource utilization and conservation efforts on a global scale thereby improving the resources management activities [
19,
28,
29,
30,
31]. Integrating LiDAR and optical remote sensing data for estimating forest biophysical parameters presents some associated challenges necessitating careful attention, robust methodologies, and many a time a hybrid method that leverages the strengths of both LiDAR and optical sensors while mitigating their limitations. Researchers have addressed these challenges through various innovative approaches for example, Hudak et al. [
32] proved that the integration of LiDAR and optical remote sensing data significantly improves forest inventory accuracy, most especially in predicting forest structure attributes when LiDAR is combined with multispectral imagery. Asner et al. [
33] stressed the importance of ground-based calibration and validation to achieve accurate carbon stock estimates with airborne LiDAR. Fernández-Landa et al. [
34] analyzed the effectiveness of combining national wide LiDAR information and optical imagery to estimate forest variables and concluded that this integration may lead to an increase in estimation accuracy most especially when the species are homogenized.
Several modelling techniques either as parametric (e.g., multiple linear regression, logistic regression, geographically weighted regression) and non-parametric (e.g., k-nearest neighbors (k-NN), random forests, artificial neural networks, gradient nearest neighbor (GNN), support vector machines, and decision trees) have been used to establish the relationship between the response variables (forest variables) and predictor variables (remote sensing derived variables) [
8,
23,
35,
36,
37,
38,
39]. Machine learning algorithms such as Random Forests, Support Vector Machines, Extreme Gradient Boosting, and Artificial Neural Networks (ANN) are some of the commonest ML algorithms used for forest variables estimation and they have demonstrated superiority and proven efficient in estimating GSV, outperforming traditional methods, highlighting their robustness across diverse datasets and varying forest ecological conditions [
8,
11,
36,
37,
40,
41].
The aim of our research was to assess the predictability of forest growing stock volume across various Cartesian planes in Estonia, corresponding to the flight schedules of ALS based on forestry purpose, by utilizing three ML algorithms to analyse NFI plot data alongside multimodal remote sensing data. We integrated ALS height percentile with multispectral data, to create a comprehensive dataset for the analysis and incorporates ground-truth data from NFI plots to enhance the accuracy of remote sensing data interpretations and model training. We applied and compared the performance of three distinct machine learning algorithms such as RF, SVR, and XGBoost to predict forest growing stock volume, focusing on their accuracy and reliability at the plot level to ensure localized accuracy in forest growing stock volume. By achieving these objectives, the research aimed to advance the understanding of forest growing stock volume estimation using multimodal remote sensing data and machine learning techniques, ultimately contributing to more effective and sustainable forest management in Estonia.
3. Results
A total of 14880 sample plots were selected from the NFI data following the removal of outliers through PauTa criteria. The share of each quadrant is represented in
Table 1 where the largest share (3852) is on the NE zone followed by the NW with 3712 sample plots while the SE and SW has a share of 3665 and 3651 respectively. Four (4) variables combinations as discussed in section 2 were used in training the three predictive models i.e., RF, SVR and XGBoost. A maximum of a 10-fold cross validation was considered sufficient enough to optimize the model according to minimum RMSE and a considerable coefficient of determination. The sample plots in each share were used to generate a corresponding model and each model’s performance metrics is explained as follows.
3.1. Performance Assessment Based on Random Forest Predictive Model
The performance metrics for the RF-based model for each combination of variables at the quadrant level are illustrated in
Table 3. In the NW zone, RMSE values ranged from 125.39 m³/plot to 141.84 m³/plot for CO3 and CO1 categories respectively, with CO4 and CO3 showing relatively close performance R² values of 0.64 (CO4) and 0.63 (CO3). The CO4 seems to have a relatively higher RMSE (126.74 m³/plot) and a higher RMSE% (13.09) as compared to the CO3 (
Table 3) with this, preference was given to the combination with the lowest RMSE and RMSE%. The highest R² (0.78) also corresponds to a larger RMSE (141.84) as with the CO1 model, this necessitates the choice of CO3 variable combination for the best performance at this zone. In the SW zone, RMSE varied from 126.85 m³/plot to 139.21 m³/plot for CO3 and CO1 respectively. In this zone, the CO4 variable combination exhibited a good fit with an RMSE of 129.22 m³/plot, R² of 0.74 and RMSE% of 14.04. Although, the CO1 and CO2 models shows a higher R² of 0.79 and 0.77 yet, their associated RMSE is larger as compare to the CO4 combination (
Table 3).
For the NE zone, CO4 yielded the best results with an R² of 0.64 and an RMSE of 133.77 m³/plot, followed by CO3 model with an R² of 0.61 and an RMSE of 132.17 m³/plot. The highest R² and RMSE under this zone also correspond to the CO2 and CO1 which is accompanied by a higher RMSE. Similarly, in the SE zone, CO4 variables performed the best with an R² of 0.70 and an RMSE of 120.56 m³/plot, followed by the CO3 variables with an R² 0.68 and an RMSE of 119.35 m³/plot. Although, CO1 had the highest R² (0.81) and RMSE (135.49), the emphasis on minimizing the RMSE favored the CO4 variables as the variability explain under this combination is on an average good.
3.2. Performance Assessment Based on Support Vector Regression
Support vector regression model carried out at quadrant level using each variable combinations is as presented in
Table 3. In the NW zone, the best performance in terms of the lowest RMSE and RMSE% is the CO3-generated model having an R² of 0.52, RMSE of 123.41 m³/plot and an RMSE% of 12.74%. However, this model captured less variability compared to its RF counterpart (
Table 3). A notable trend of a consistent highest R² value across each quadrant is exhibited by CO1 model and also, it is characterized by the highest RMSE. In the NW zone, for instance, the CO1 variable achieves an R² of 0.81 and an RMSE of 150.01 m³/plot, indicating significant error when applied to the 3712 plots within the zone.
In the SW zone, CO1 variable produced an R² value of 0.79 and an RMSE of 141.35 m³/plot, while CO4 and CO3 both achieved an R² value of 0.64. However, the CO3 model has the lowest RMSE (124.53 m³/plot) compared to CO4 (129.57 m³/plot). Despite CO2's relatively average performance with an R² of 0.74 and an RMSE of 136.17 m³/plot, the CO3 model is deemed more suitable with an RMSE of 124.53 m³/plot. The range of coefficient of determination in the NE zone is between 0.43 – 0.65 where the CO1 and CO2 explain the same variability with an R² of 0.65. However, the RMSE in CO1 is 4.37 m³/plot higher than CO2. The CO3 model exhibits the lowest RMSE (129.83 m³/plot) and RMSE% (12.14), although its R² (0.43) is the lowest among all variable combinations. On average, CO4 performs notably better with an R² of 0.60 and an RMSE of 138.98 m³/plot. In the SE zone, CO1 demonstrates the highest R² (0.86) and RMSE (144.46 m³/plot), while CO3 exhibits the lowest R² (0.63) and corresponding lowest RMSE (116.43 m³/plot). Although the CO3 and CO4 models shows a relatively good fit, selecting CO4 (R² = 0.68, RMSE = 122.34 m³/plot) is justifiable. In comparison with the RF model, a better precision is obtained with the RF using CO4 model.
Table 4.
Support Vector Regression performance metric.
Table 4.
Support Vector Regression performance metric.
Variables |
|
R² |
RMSE |
RMSE% |
CO1 |
NW |
0.81 |
150.01 |
15.49 |
CO2 |
|
0.53 |
130.53 |
13.48 |
CO3 |
|
0.52 |
123.41 |
12.74 |
CO4 |
|
0.51 |
125.05 |
12.91 |
CO1 |
SW |
0.79 |
141.35 |
15.36 |
CO2 |
|
0.74 |
136.17 |
14.79 |
CO3 |
|
0.64 |
124.53 |
13.53 |
CO4 |
|
0.64 |
129.57 |
14.08 |
CO1 |
NE |
0.65 |
153.53 |
14.06 |
CO2 |
|
0.65 |
149.16 |
13.95 |
CO3 |
|
0.43 |
129.83 |
12.14 |
CO4 |
|
0.6 |
138.98 |
13 |
CO1 |
SE |
0.86 |
144.46 |
14.76 |
CO2 |
|
0.79 |
132.45 |
13.52 |
CO3 |
|
0.63 |
116.43 |
11.89 |
CO4 |
|
0.68 |
122.34 |
12.5 |
3.3. Performance assessment based on Extreme gradient boosting
The XGBoost outperform the SVR in the NW zone meanwhile, a similar performance was achieved within the variable combination of RF and XGBoost models. Notably, for the CO4 variable, both models (XGBoost and RF) had the same R² value of 0.64 and closely correlated RMSE values, differing by only 0.59 m³/plot. For the CO1 model, both XGBoost and RF also achieved an R² of 0.78, though XGBoost had a higher RMSE of 1.8 m³/plot compared to the RF model. An R² of 0.61 for XGBoost was found in the CO3 model, with an RMSE of 124.63 m³/plot, highlighting the robustness of CO3 variables for growing stock volume estimation.
In the SW zone, we compared the unit differences in performance metrics and found that the CO4 model-based algorithm met our goals best, with an R² of 0.69 and the lowest RMSE of 128.29 m³/plot. This was followed by the CO3 model with an R² of 0.68 and an RMSE of 129.32 m³/plot. As with the other models, the CO1 variable yielded the highest R² of 0.77, with an RMSE of 139.39 m³/plot. At the NE zone, the RMSE are in the range of 132.58 m3/plot to 148.16 m3/plot which is similar to the range obtained with RF. The CO3 produced a more accurate estimate considering its R² and RMSE of 0.60 and 132.58 m3/plot respectively. In the NE zone, the RMSE values range from 132.63 m³/plot to 148.16 m³/plot, similar to the range obtained with RF. The CO3 model provided a more accurate estimate, with an R² of 0.60 and an RMSE of 132.63 m³/plot.
Comparing the results obtained in the SE zone based on the four variable combinations, we found the CO3 model to be the most precise, with an R² of 0.65 and an RMSE of 119.48 m³/plot. The CO4 model also achieved an R² of 0.65, but its associated RMSE is 123.36 m³/plot, which is 4.18 m³/plot higher than CO3, making CO3 the preferred model. The CO2 model shows an increased R² of 0.75 but a decreased accuracy, with an RMSE of 131.40 m³/plot. Although the CO1 model has the highest R² of 0.80, its associated RMSE of 138.01 m³/plot indicates a notable decrease in predictive power.
Table 5.
XGBoost performance metric.
Table 5.
XGBoost performance metric.
Variables |
|
R² |
RMSE |
RMSE% |
CO1 |
NW |
0.78 |
143.67 |
14.84 |
CO2 |
|
0.68 |
133.77 |
13.81 |
CO3 |
|
0.61 |
124.46 |
12.85 |
CO4 |
|
0.64 |
126.15 |
13.03 |
CO1 |
SW |
0.77 |
139.37 |
15.19 |
CO2 |
|
0.74 |
136.17 |
14.79 |
CO3 |
|
0.68 |
129.32 |
14.05 |
CO4 |
|
0.69 |
128.29 |
13.94 |
CO1 |
NE |
0.69 |
139.63 |
13.09 |
CO2 |
|
0.72 |
148.16 |
13.86 |
CO3 |
|
0.6 |
132.63 |
12.4 |
CO4 |
|
0.59 |
134.25 |
12.56 |
CO1 |
SE |
0.8 |
138.01 |
14.1 |
CO2 |
|
0.75 |
131.4 |
13.4 |
CO3 |
|
0.65 |
119.48 |
12.17 |
CO4 |
|
0.65 |
123.36 |
12.6 |
3.4. Model Predictive Power at Each Quadrant/ Comparative Analysis of Predictive Model
In all quadrants, the RF algorithms proves to be the most effective. For instance, RF algorithm with the CO3 variable combination proves to be the best predictive model in the NW zone, with the SVR-based model achieving the lowest RMSE of 123.41 m³/plot, followed by the CO3-XGBoost model with an RMSE of 124.63 m³/plot; however, considering the differences in R² values and associated RMSEs, the RF-based model with an R² value of 0.63 and an RMSE value of 125.39 m³/plot is more preferable (
Figure 3 NW_RF), and the CO3 variable combination demonstrates robust performance, consistently yielding the lowest RMSE values across all three predictive model. In the SW zone, the RF algorithm again outperforms other models in estimating GSV based on the CO4 variable combination (
Figure 3 SW_RF), although the SVR model achieved the lowest RMSE, it also showed more scattering with an R² of 0.6; shifting focus to the NE zone, the RF algorithms using the CO4 variable combination yield more preferable results (Figure NE_FR), followed by the XGBoost model using the CO3 variable combination, with the least effective results found in the SVR model (
Table 6); turning to the SE zone, the use of RF-based algorithms with the CO4 variable combination also shows appreciable performance (Figure 4 SE_RF), outperforming the SVR and XGBoost models (
Table 6), yet as observed in other zones, the lowest RMSE does not necessarily correspond to the RF-based model; nonetheless, the minimal differences in performance metrics among these predictive models favour the RF model, and generally, the SVR yields the lowest RMSE across all quadrants and predictive models, but its CO1 and CO2-based models are characterized by relatively high RMSE values.
4. Discussion
In this study, we demonstrated the prediction of forest tree growing stock volume using three predictive algorithms such as RF, SVM and XGBoost. In other to enrich our prediction, we made use of Estonian National Forest Inventory data, Sentinel 2, and LiDAR height percentile data as the use of multi-source data has the potential of improving the prediction of growing stock volume [
37,
40,
59]. Our study presented a pioneering model estimation of plot GSV using a complete cycle NFI over Estonia forest and at such provide an architectural framework for comparative study in the country. We used four variable combinations derived from the remotely sensed data as explained in section 2 of this article. These includes; LiDAR and Vegetation Indices (CO1), Lidar and band reflectance (CO2), Vegetation Indices and band reflectance value (CO3), band reflectance value, vegetation indices and LiDAR height percentile (CO4). In our characterization, we examined which of the machine learning algorithms performs best and also delved into variable importance and their associated sources.
Notably, in all the quadrant, the RF based algorithm outperformed other predictive algorithms (
Table 6). Although the lowest RMSE does not necessarily correspond to the RF model, nonetheless, the associated unit differences between the performance metric favors the RF model. For example, in the NW, the RF achieved an R² of 0.63 and an RMSE of 125.39 m³/plot, followed by the XGBoost with an R² of 0.64 and an RMSE of 126.15 m³/plot, the lowest RMSE correspond to the SVR algorithm with an R² value of 0.52 and an RMSE value of 123.41 m³/plot. The lowest RMSE of 124.53 m³/plot in the SW also corresponds to the SVR whereas, the best predictive model was found in the RF models with an R² 0.74 and an RMSE of 129.22 m³/plot. This further establish the strength of RF algorithms in estimating forest variables.
The use of RF to estimate GSV has been demonstrated overtime, in some cases, it is used exclusively [
14,
37,
60] and in some, it was used comparatively with other ML algorithms [
36,
37,
40]. For example, Chirici et al. [
40] used two non-parametric and two parametric prediction methods (RF, K-NN, MLR and GWR) to produce a wall – wall spatial prediction of growing stock volume using Italian NFI plot and remotely sensed data and concluded that RF forest algorithm produced the best accuracy. Although the variable used in their prediction is not the same as ours. In our study, we do not consider the environmental factors, nor distinguish the forest type and species at each quadrant level this of course we believe would have help to further homogenized our findings. For example, Estonia is characterized by both coniferous and deciduous tree species. Lang et al. [
23] in their research to facilitate the construction of maps for forest height, standing wood volume and tree species composition in Estonia suggested that for a more reliable prediction, empirical data can be divided into subset that account for variability before fitting the model.
The use of optical remote sensing data for the prediction of forest variable is usually accompanied with issue relating to saturation and this is believed to have effect on the accuracy of the predicted variables. Areas within the forest with dense canopy closure are highly susceptible to saturation problem as the optical sensor only capture horizontal information about the forest. One of the ways to possibly reduce this impact is the inclusion of vegetation indices derived from red edge band [
61,
62] and also through the inclusion of satellite data that captures the vertical information about the forest structure. Hence, we incorporated the height percentile data from LiDAR metric and vegetation indices derived from red edge bands (table 2). As seen from the results of our findings, combination that incorporated vegetation indices, band reflectance values and height percentile data consistently yielded the best results in each quadrant except in the NW quadrant where combination of bands reflectance number and height percentile is considered to be more accurate (table 3). Although the unit difference between the performance metrics of the CO3 and CO4 variable combinations at this quadrant is not well pronounced having an RMSE of 125.39 m³/plot and an R² value of 0.63 in the CO3 and an RMSE 126.7 m³/plot and an R² value of 0.64 in the CO3 (table 3).
When band reflectance and height percentiles were used, this approach yields the lowest RMSE across all quadrants and the three predictive models (
Table 3,
Table 4 and
Table 5). Additionally, incorporating vegetation indices further improves the coefficient of determination (R²). For example, in the SW zone, adding vegetation indices improves the R² by 4.11% and decreases the RMSE by 1.56%. Similarly, in the NE zone, there is a 4.7% improvement in R² and a 1.19% decrease in RMSE. In the SE zone, the R² increases by 2.86% and the RMSE decreases by 1.13%. The CO1 and CO2 combination consistently yields significantly higher R² values and an erroneous RMSEs (
Table 3,
Table 4 and
Table 5).
Lahssini et al. [
59] was of the opinion that the availability of multimodal remote sensing data alone does not necessarily guarantee better performance; instead, how these data are integrated is crucial.
Appendix A shows how the various variable were integrated and ranked in order of their importance. Based on the prediction in the NW quadrant, LiDAR height percentile and band reflectance number shows a good fit in GSV prediction with height percentile at p95, p90 and p50 the most important height percentile variable this had been established by other studies [
63,
64,
65] and the red edge band 7 and 5 are the most important band reflectance number this correspond to other research where the red edge bands had been established to outperformed other bands in forest variable prediction [
8,
66,
67].
In the SW, NE and SE quadrat where the combination of height percentile, individual band reflectance value and vegetation indices proved to be the most efficient combination in GSV prediction, the height percentile at p95, p90, p80 and p50 are crucial among the height percentile variables, followed by individual band reflectance values. This correspond to [
63,
64,
68] where these height percentiles were also reported to contribute effective to forest variables prediction. Among the individual band reflectance number, the red edge and SWIR are the most commonly reported bands that correlate with GSV. Our findings further establish this where the SWIR1, SWIR2, red edge bands were part of the most important variable that contribute significantly to the GSV prediction except in the NE quadrant where red edge band 6 was removed.
In our research, we only examine the effect of the inclusion of various band reflective value, vegetation indices and height percentile data on the overall performance and reliability of GSV prediction meanwhile the initial performance of ML algorithm can be improved with hybrid model [
37,
69]. For instance, in Hu et al. [
37], accuracy of three machine learning including RF, SVR and an artificial neural network (ANN) were improved by building a hybrid model as an extension of these ML these later yielded an improvement by reducing the initial RMSE of RF by 11.45 %, and an increase in the associated R² by 2.17% for the SVR, the RMSE was reduced by 20.37% and the R² increased by 7% these are called ordinary kriging model hybrid.
5. Conclusions
Estimating and quantifying forest related variables such as aboveground biomass and growing stock volume are becoming increasingly important and there already exist considerable research finding on this, nonetheless, there exist a dearth of comprehensive research effort in Estonia. The primary aim of this study was to enhance the estimation of growing stock volume by integrating data from National Forest Inventory, Sentinel 2, and LiDAR, employing three predictive algorithms such as RF, SVR, and XGBoost. Our findings reveal that the combined use of height percentile, band reflectance value, and vegetation indices significantly improves the accuracy of GSV predictions. The combination of band reflectance value and height percentile also exhibited a promising trend. Among the predictive algorithms tested, RF demonstrated the best performance, yielding the most considerable accuracy at each quadrant. The most occurring height percentile variable are P95, p90 and p50 while the occurrence of the red edge and SWIR bands are significant notable. Despite the promising outcomes of our study, certain research gap exists that need addressing. For example, the variability of forest types and conditions across different geographical regions suggests that further testing and refinement of these methods are necessary. Efforts should attempt to homogenized the tree species in each quadrant into individual strata before modelling and additionally, integrating other data sources, such ecological, climatological, and experimenting with emerging machine learning techniques could further enhance prediction accuracy. In conclusion, the enriched estimation of growing stock volume using NFI, Sentinel 2, and LiDAR data, combined with advanced predictive algorithms, offers a robust tool for forest management. This approach promises significant benefits for sustainable forestry practices and ecosystem monitoring, paving the way for more informed and effective climate smart decision-making in forest conservation and management.