Spatially explicit estimation of forest structural compartments has become a main target of studies and forest monitoring during the last decades (Lim
Landsat image time series with simultaneous pre-processing of spectral data that are freely available on the cloud-computing platform Google Earth Engine have created new capabilities for mapping of forests and their biophysical parameters and features (Cohen & Goward, 2004). New temporally continuous and spatially explicit approaches have significantly empowered a robust estimation of deadwood (Gonzalez
Herewith, precise estimation of forested areas and quantification of ecosystems’ biophysical values is linked to the various statistical classification and regression approaches. For the last decades, non-parametric machine learning methods have shown a strong capability to be applied for forest remote sensing data. A number of studies refer to successful implementation of numerous classifiers (e.g., decision trees (Ruefenacht
To date, ensemble machine learning methods are successfully used in remote-sensing applications (Belgiu & Dragut, 2016). Using a set of so-called weak learners (mainly decision trees), these models apply either a bootstrap or bagging process (well-known Random Forest, RF), or boosting fitting (e.g., gradient boosting machines, GBM). The main challenges those ensemble and above-mentioned methods face are the quality of reference and ancillary data, while tuning of modelling hyper-parameters (e.g., learning rate of number of iterations) was reported as an issue more dependent on dataset features (Maxwell
In the twenty-first century, Ukraine is facing both an increasing interest in a full carbon accounting of forest ecosystems, including estimation of deadwood stocks, and substantial challenges due to the absence of the main sources of reference data (e.g., sampling plots-based National Forest Inventory, LiDAR applications, Shvidenko
Here we aim to (
A study site in Ukrainian Polissya was chosen with the aim to fully capture regional geographical and forest composition features. The study was carried out for a site in the Chernihiv region (Northern Ukraine, Figure 1) with an area of
The local climate represents the continental climate subtype, with a mean annual temperature of 7.2 °C and an average total precipitation of 702 mm. Soils are mainly sod-podzolic, with alluvial soils in river banks. Considering the forest type classification implemented in Ukraine, Scots pine grows there in A2 and B2 types (i.e., poor and semi-poor level of soil fertility and fresh sites in terms of humidity), silver birch – B2 and C3, and black alder – C4 and C3 (soils with a rich nutrient content and high humidity level; for more detailed classification of forest types see Bilous
Forest inventory maps delineating forest polygons have been incorporated in the study as a source of forest attributes derived from the forest planning and management (FPAM) data (available for 2011). We considered utilizing FPAM data instead of statistical forest inventory datasets, since these are currently limited to several regions (Myroniuk, 2018) which do not include the study area. FPAM data for these stands contain information on dominant tree species, mean diameter at breast height, mean height, age, site index, relative stocking, forest type, and GSV estimated during the last forest accounting. In total, we used data from 2,125 primary units available from the FPAM dataset for the study area. These polygons were given as ESRI shapefiles.
Stocks of CWD biomass were estimated using allometric linear models developed for the Ukrainian Polissya region (Lakyda
Accuracy and reliability of produced maps is strictly dependent on the agreement between the used terrestrial-based and remote sensing data. Considering that FPAM data for the study area was available only for 2011, the same point of time was chosen for obtaining remote sensing data. We used the open United States Geological Service archive presented in the Google Earth Engine (GEE, Gorelick
The main workflow of this study is illustrated in Figure 2. Further, the methods of classification and data imputation are described as given in this flowchart. For the land cover classification task, the RF classifier was chosen as the most popular among the models used in the various scientific papers (e.g., Jhonnerie
We provided a two-step classification of Landsat mosaics. First, land cover classes were predicted using the RF machine learning algorithm to produce a forest cover mask. Then, the classification on dominant tree species across forest stands within the produced mask was performed. For both steps we used the same set of 6 predictors: Normalized Vegetation Difference Index (NDVI), Normalized Water Difference Index (NDWI), Tasselled Cap Transformation (TCT, namely brightness TCB, greenness TCG and wetness TCW), and elevation from the Digital Elevation Model (DEM) downloaded from the GEE for the study area. We chose NDVI and TCT as predictors because of their ability to capture GSV and biomass stock changes (Kauth & Thomas, 1976). In order to better classify water bodies, we utilized the NDWI predictor here. As an auxiliary non-spectral variable, elevation from the DEM is useful for distinguishing black alder from other broad-leaved species, since these stands are more likely to grow in lower sites regarding elevation above sea level (Lakyda
Multispectral data for land cover classification was extracted for a set of 1,010 sample points, randomly distributed across the study area. Sample plots were generated in the Quantum Geoinformatics System (QGIS 3.4) software, and land cover classes were then visually interpreted using freely available software from Google Inc. (Google Earth for visual interpretation). Each sampling point was visually inspected using high-resolution images available on Google Earth as a reference to assign one of eight thematic classes: bog, croplands (including gardens), forests (including the stands of linear spatial shape in the fields, across roads and rivers), grasslands (including unforested sites surrounded by forests, but without signs of clear-cuts being identified with the time series of high-resolution images available on Google Earth), settlements, shrublands, water bodies (rivers and lakes) and other (sands, rocks).
The validation dataset for the land cover classification was also compiled independently using the Collect Earth plugin and included 964 points, different from points used in the training dataset.
The dataset for tree species classification was extracted from the FPAM database for 742 points that fell within forested areas. All points were carefully interpreted on Google Earth to ensure they reliably represented characteristics of the inventory polygons which they were linked to. For the training model we used 557 points (75%) and the rest (185 points) were used for the validation task. We considered four dominating tree species in the study: PISY (acronym for Scots pine), BEPE (acronym for silver birch), ALGL (acronym for black alder), and POTR (acronym for common aspen). Other species were grouped in the thematic class OTHER.
RF is an ensemble method that decreases the variance of output, using a set of simple learners (here Classification and Regression Trees, CARTs) that preferably give outputs with low bias. The training dataset is split into small subsets used by each learner in the RF model in an independent way, then the results are averaged (Breiman, 2001) to get the best performance (a process called “bootstrap”, or “bagging”) and the lowest variance. We used the “randomForest” R-package (Liaw & Wiener, 2002) to classify land cover. Default parameters (number of decision trees
Performance of the RF models was assessed in terms of the following metrics: overall accuracy and the kappa (Congalton, 1991) index for all datasets, and user accuracy and producer accuracy for the class “forest”. Confusion matrixes were built in the “caret” R-package (Kuhn, 2008).
We used a training dataset that consisted of 439 sample stands with information on deadwood biomass stocks calculated for each forest stand using the abovementioned allometric models. Such a number of stands and points was chosen because of full spatial consistency between FPAM-based polygon map data and the forest mask obtained by the machine classifier for these stands. Imputation of continuous values of deadwood biomass stocks was performed using
Since the allometric models of deadwood biomass are consistent with the parameters of a growing forest stand, we directly model this compartment instead of linking to GSV values (Pflugmacher
For imputation purposes, the performance of the
GBM was implemented in the “gbm” (Greenwell
Hyperparameter grid for tuning GBM models in this study.
Parameter to tune | Hyperparameter values | |||||
---|---|---|---|---|---|---|
Shrinkage | 0.01 | 0.05 | 0.02 | 0.001 | ||
Interaction depth | 2 | 4 | 6 | 8 | 10 | |
Number of trees | 1000 | 5000 | 10000 | |||
Number of observations in the terminal nodes | 5 | 10 | 15 | |||
Bagging fraction | 0.5 | 0.75 | 1.0 |
The shrinkage parameter (Table 1) defines a learning rate applied to each CART in the expansion. Interaction depth shows the highest level of variable interactions allowed in each tree. The number of trees shows how many CARTs were utilized in the GBM model. Bagging fraction is a proportion of observations randomly chosen to propose the next CART in the expansion, i.e., making the given GBM a stochastic boosting model (if this fraction is less than 1.0).
We assessed the produced GBM models, considering their root mean square error (RMSE) produced by a random inner validation subset (25% of the whole training dataset, as the training fraction for all models was 0.75). Performance of the GBM and
Additionally, we compared model performance on a test landscape which was not covered by the forest stands from the training dataset. This landscape (1,020 ha) represents forest communities typical of the entire study area: Scots-pine dominated stands with admixtures of silver birch and rare deciduous forests. Mean and total
The RF model produced has shown quite good overall accuracy (81.9%) and the kappa (74.6%) index. Land cover types “water”, “forest” and “croplands” on the validation dataset were interpreted with a high level of user accuracy (93.1%, 92.2%, and 90.0%, respectively), while the largest misclassifications were met for the types “bog”, “other” and “shrublands” (Table 2).
Confusion matrix of RF land cover classification.
Predicted | Observed | |||||||
---|---|---|---|---|---|---|---|---|
Bog | Croplands | Forest | Grasslands | Others | Settlement | Shrubland | Water | |
Bog | 0 | 1 | 0 | 0 | 0 | 1 | 0 | |
Croplands | 0 | 11 | 10 | 19 | 12 | 7 | 0 | |
Forest | 6 | 5 | 6 | 0 | 1 | 11 | 0 | |
Grasslands | 1 | 22 | 2 | 0 | 1 | 19 | 1 | |
Others | 0 | 0 | 0 | 0 | 1 | 0 | 0 | |
Settlement | 0 | 4 | 2 | 0 | 2 | 1 | 0 | |
Shrubland | 1 | 6 | 6 | 14 | 0 | 0 | 1 | |
Water | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
The tree classification model has shown relatively good performance (overall accuracy is 78.9%, the kappa index is 56.4%), with user accuracy ranging from 90.0% for Scots pine to 60.0% for silver birch and 80.0% for black alder. The class “OTHER” and common aspen achieved worse accuracies (23% and 29%, respectively).
Elevation was the most important variable to predict forests of black alder (Figure 5), while TCG and TCW contributed the most to the prediction of Scots pine stands. All predictors rather weakly classified less abundant broadleaved tree species (silver birch, common aspen and other), which resulted in higher errors (Table 3). Figure 6 illustrates the outputs of both RF classification models.
Confusion matrix of RF tree species classification.
Predicted | Observed | ||||
---|---|---|---|---|---|
ALGL | BEPE | OTHER | PISY | POTR | |
ALGL | 1 | 5 | 4 | 3 | |
BEPE | 0 | 3 | 7 | 0 | |
OTHERS | 0 | 0 | 1 | 2 | |
PISY | 3 | 7 | 2 | 0 | |
POTR | 0 | 0 | 0 | 1 |
According to the map provided in Figure 6a, forest cover within the study area in 2011 was 40,367 ha, or 32.0% of the total area. This forested area was covered (Figure 6b) by stands of the following species (with 95% confidence interval, CI): 3320 ± 1328 ha of black alder forests (8.2 ± 3.3%), 3958 ± 1459 ha of silver birch forests (9.8 ± 3.6%), 1893 ± 1310 ha of common aspen forests (4.7 ± 3.2%), 27862 ± 1979 ha of Scots pine forests (69.0 ± 4.9%), and 3334 ± 1634 ha of stands where other species grew (8.2 ± 4.0%).
Tuning the GBM algorithm has shown that the lowest RMSE (2.09 t·ha−1, or 26% of the mean deadwood biomass stock) on the random 25% inner validation subset was achieved by a model with the following hyperparameters: number of trees = 1,000, shrinkage = 0.01, minimum number of observations in the terminal node = 15, bagging fraction = 0.5, interaction depth = 10. The other nine best GBM models produced quite similar errors (ranged from 2.087 to 2.093 t·ha−1). According to their hyper-parameters, the best learning rate for training the GBM in this case is quite high – 0.01. As well, all the best models used a bagging fraction value of 0.5 and a minimum number of observations in the terminal node of 15. Interaction depths there are mostly high (> 6), i.e., decision trees with a low number of “branches” did not produce the lowest error. In addition, models with various numbers of CARTs utilized were almost equal in terms of RMSE on the inner validation subset (25% of the total dataset).
When predicting deadwood biomass on the training dataset without an inner validation subset, there was a lower RMSE – 1.3 t·ha−1, or 16% of the mean deadwood biomass stock. According to the validation results, all 10
Importantly, low deadwood biomass values (2–3 t·ha−1) were predicted by the GBM model (Figure 7b) as two or three times higher. Therefore, the
Comparison of mean and total values of deadwood biomass stock within the study area produced by the best performing models.
Type of model | Tree species | ||||
---|---|---|---|---|---|
ALGL | BEPE | OTHER | PISY | POTR | |
Mean ± | |||||
7.0 ± 3.4 | 8.5 ± 1.9 | 8.1 ± 2.8 | 8.6 ± 1.9 | 8.1 ± 2.1 | |
GBM | 6.8 ± 2.3 | 8.0 ± 1.7 | 8.1 ± 1.7 | 8.4 ± 1.5 | 8.1 ± 2.1 |
Total, thousands t | |||||
39.0 | 34.8 | 19.9 | 237.6 | 4.9 | |
GBM | 37.9 | 32.8 | 19.9 | 232.0 | 4.9 |
The total deadwood biomass stock within the study area produced by the
Maps of predicted dead biomass stock for test landscapes (with 95% CI) are given in Figure 8, with a respective comparison in Table 5.
Comparison of mean and total values of deadwood biomass stock for the test landscape at polygon level.
Type of model | Mean ± CI, t·ha−1 | Total stock, t | Mean difference with reference per polygon, t·ha−1 |
---|---|---|---|
7.8 ± 0.3 | 8011 | – | |
9.0 ± 0.1 | 9107 | +4.5 | |
GBM | 8.7 ± 0.1 | 8996 | +4.1 |
While there were no likely large differences between model outputs for the test landscape (Figure 8), there was at least 13% total overestimation of deadwood biomass stock compared to reference FPAM data (Table 5). However, mean differences at polygon (stand) level could reach half of the mean dead biomass stock.
To date, Ukraine is facing increasing challenges due to the lack of precise, statistically-based NFI data (Lesiv
The spatial resolution of Landsat-5 TM imagery can still be considered as sufficient for the classification models while using the tools for validation purposes like Google Earth and Collect Earth. The cloud-masked mosaic, even with a limited period (one vegetation season in 2011), reliably produces a multispectral dataset without noises. The forest RF-based mask produced in this study may be compared with the forest cover map proposed by Global Forest Change (GFC, Hansen
Some challenges linked to local geographical and biophysical features were met in this study. Linear and relatively narrow forest stands are common in Ukraine (shelterbelts, across roads and rivers, etc.). Myroniuk (2018) reported that the spatial resolution of Landsat bands is too coarse for proper classification of such forests surrounded by fields (Figure 9b). However, utilization of the NDWI as a predictor can be sufficient to provide reliable distinction of forest patches in river banks from water bodies (Figure 9a) and thus create a more precise forest mask. This task can be solved by using a dense Landsat time series (Cohen & Goward, 2004).
Because of the sparse distribution of silver birch stands within the study area, tree species classification (Table 3) has shown strong biases in distinction of this species from Scots pine forests. Lakyda
Successful use of the
In this study we used linear allometric models to obtain data on deadwood biomass stocks. These models were developed based on the biophysical parameters of a growing stand and thus are strictly linked to natural mortality processes in different successional stages and reflect local site productivity. However, such models cannot predict deadwood stocks after abrupt events (natural disturbances at different scales) and different management considerations common in Ukraine (Bilous
Deadwood stocks are strictly related to various issues ranging from biodiversity conservation to carbon management, but hitherto are not reflected in operational forest management data. However, those can be considered in the forest inventory agenda of countries with developing or transition economy utilizing free available remote sensing data. To meet climate mitigation and ecosystem diversity maintenance targets, there is credible capability to apply a bundle of novel ensemble machine learning methods, e.g. Random Forest and Gradient Boosting Machines, with satellite multispectral data. These non-parametric approaches, together with the traditional
Assessment and mapping tasks of deadwood in Ukraine require urgent implementation of a national statistical inventory, preferably with utilization of diverse sources of data (like point clouds provided by LiDAR). Also, the proposed approach matches with the necessity to estimate dead biomass, while this ecosystem compartment is crucial to maintain forest resilience under the changing global climate.