An indirect approach to predict deadwood biomass in forests of Ukrainian Polissya using Landsat images and terrestrial data

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 ca 1,192 km², between N 50°54′1″, E30°44′60″ and N51°14′15″, E31°12′16″. The dominated land cover types (according to the local operational forest management data) are croplands, grasslands in Desna river bank and homogeneous coniferous replanted forests (with Scots pine (Pinus sylvestris L.) as the main tree species). Some forest stands are formed by silver birch (Betula pendula Roth.) and common aspen (Populus tremula L.); natural forests of black alder are also presented. Pine-dominated forests are mainly pure, with insignificant admixtures (i.e., the share of mean GSV is ca less than 10%) of silver birch. Forests within the study area are managed by two main authorities: State Forest Enterprise “Osterskiy” and Regional Landscape Park “Mezhyrichensky”.

Figure 1

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 A₂ and B₂ types (i.e., poor and semi-poor level of soil fertility and fresh sites in terms of humidity), silver birch – B₂ and C₃, and black alder – C₄ and C₃ (soils with a rich nutrient content and high humidity level; for more detailed classification of forest types see Bilous et al., 2019).

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 et al., 2019). According to the classification implemented in Ukraine (Bilous et al., 2019), CWD is divided into such well-common deadwood components as snags (i.e., dead standing trees), logs (i.e., downed stems and stumps) and litter of coarse (d > 1 cm) branches. In this study, only biomass stocks of snags and logs were estimated: according to Lakyda et al. (2019), the litter of coarse branches weakly correlates with spectral data and mainly depends on the soil and fine litter features. The sum of the studied stocks further is considered as deadwood biomass stock. The allometric model used is given as (1): (1)M=a1⋅Da2⋅Ha3⋅Sa4M = {a_1} \cdot {D^{{a_2}}} \cdot {H^{{a_3}}} \cdot {S^{{a_4}}} where M is the biomass stock of snags or logs (t·ha⁻¹), a₁, a₂, a₃ and a₄ are regression parameters, D is the diameter at breast height (cm), H is the mean height (m), S is the relative stocking (ranging from 0 to 1). According to the source of these models (Lakyda et al., 2019), there was a relatively small root mean square error (3% of the estimated mean). This error was ignored in the following modelling and uncertainties estimations. For more detailed information on the regression parameters and predictor sets across compartments and tree species, see the respective paper (Lakyda et al., 2019).

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 et al., 2017) to download Landsat-5 Thematic Mapper (TM) images (with 30 m spatial resolution) and create a cloud-free mosaic for the given site and time frame from 22 May 2011 to 22 September 2011 (according to the duration of the vegetation season in Ukraine). Landsat-5 TM imagery was already geometrically corrected. Radiometric correction to the surface reflectance and mosaicking with the aim of masking pixels with clouds were conducted using the GEE. The mosaics were then downloaded from the GEE and processed in the R environment (R Core Team, 2019).

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 et al., 2015; Thanh Noi & Kappas, 2018). Then, this type of a model was utilized for the tree species classification task (e.g., Lang et al., 2018). The relatively novel approach of gradient boosting (e.g., Sachdeva et al., 2018) was compared to the k-NN, as the last method is commonly used for regression and imputation purposes (e.g. Lang et al., 2014; Latifi et al., 2015; Fu et al., 2019).

Figure 2

Classification reference data

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 et al., 2019). Ancillary datasets are illustrated in Figure 3.

Figure 3

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.

Imputation of deadwood biomass

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 k-NN and GBM models based on the extracted mean spectral values of all pixels within the polygon of each stand from the sample dataset.

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 et al., 2012). This decision was made based on an apparently strong relationship between GSV and deadwood biomass values obtained from FMAP data for the study area (Figure 4).

Figure 4

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

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.

Figure 5

Figure 6

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⁻¹, 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⁻¹). 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⁻¹, or 16% of the mean deadwood biomass stock. According to the validation results, all 10 k-NN models have shown rather similar performance: RMSE ranges from 1.42 to 1.51·ha⁻¹. The lowest RMSE was achieved by the model with k = 1 (1.42 t·ha⁻¹, or 17 % of the mean deadwood biomass stock). The best GBM and k-NN models thus were used to impute on the spatial dataset (Figure 7).

Figure 7

Importantly, low deadwood biomass values (2–3 t·ha⁻¹) were predicted by the GBM model (Figure 7b) as two or three times higher. Therefore, the k-NN model (Figure 7a) visually has produced more “outliers” throughout the graph. A comparison of mean (with standard deviations, SD) and total deadwood biomass outputs for the entire study area is given in Table 4.

The total deadwood biomass stock within the study area produced by the k-NN model is 0.34 Mt, while the best GBM resulted in similar stock (0.33 Mt, Table 4). The highest differences in the mean stock between model outputs were found for silver birch, and, with less extent, for Scots pine and black alder.

Maps of predicted dead biomass stock for test landscapes (with 95% CI) are given in Figure 8, with a respective comparison in Table 5.

Figure 8

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 et al., 2018). A pilot project run in two regions (Sumy and Ivano-Frankivsk) a decade ago showed both new evidence (forest biometric parameters of stands located outside state-owned forests, deadwood stocks) on the state of local Ukrainian forests and a strong demand for further expanding such inventory (Myroniuk, 2018). However, Ukraine does not have a full national forest inventory due to a number of reasons (Bilous et al., 2019). Since the quality of ancillary data strictly defines the performance of any classification of remote sensing data, uncertainties and mistakes always occur while working with imprecise data. In this study, FPAM data also had a limited temporal frame (2010–2011), so only Landsat-5 TM imagery could be used for land cover classification and further processing, while the Landsat-8 Operational Land Imager archive has been available for years. FPAM data in Ukraine is also commonly recognized as biased and limited in precision, strictly focusing on management purposes and neglecting the ecological functions of forests (Lesiv et al., 2018). One possible suggestion is to carry out independent research of forest inventory (Bilous et al., 2017). However, the use of such an approach for larger areas as proposed in this study requires substantial time and funds.

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 et al., 2013). The GFC mask was utilized with a 40% canopy closure threshold, which is recognized as a suitable value for distinguishing forests in Ukraine (Myroniuk, 2018). The forested area according to this data was 31.2%, thus closely matching the value produced by our RF model (32.0%).

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

Figure 9

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 et al. (2019) indicated that there was uncertainty in distinguishing between black alder and silver birch, since these species had similar spectral reflectance. This issue was solved by adding soil maps as an ancillary dataset, while in this study the non-spectral predictor DEM was important for delineating stands of black alder.

Successful use of the k-NN method for forest inventory purposes throughout the world (McRoberts, 2012) has proven its capability to provide robust data imputation. However, this approach is sensitive to the quality of field-sampled data (Tomppo et al., 2017). Gradient boosting machines are commonly used for classification purposes and often compared to another ensemble method – RF, so there is some lack of information for tests of the imputation performance (Dube et al., 2014). In this study we have observed that GBM can slightly outperform the k-NN model in terms of accuracy (Figure 7). On the other hand, k-NN does not require control of as many modelling parameters, so it is more user-friendly. The GBM tuning algorithm needs several hours to run depending on the number of hyperparameter combinations and computation power. However, it allows inner validation using a training fraction parameter, k-fold cross-validation and stochastic nature through the bagging fraction hyperparameter. Because of these advantages, XGBoost (Chen & Guestrin, 2016), the next extension of the AdaBoost (and thus GBM) framework, has shown the best accurate results in various data science competitions working with structured datasets. Relatively low attention to GBM capabilities in forest imputation tasks requires further utilization of this method to obtain more robust woody stock and biomass predictions. Possible limitations to such use, compared to the simple and classic k-NN method, could be associated with local spectral features, the size and quality of the training dataset, nature of the satellite spatial and temporal resolution.

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 et al., 2019). A network of sample plots with a reasonable allocation scheme can be harnessed for this issue but requires data from statistical national inventory (currently available only for two regions in Ukraine), global datasets, such as Schepashenko et al. (2019), or substantial costs to carry out the LiDAR measurements (Gonzalez et al., 2018). While we estimated deadwood indirectly through satellite images, more precise data could be directly obtained with the use of laser scanning (Zald et al., 2016).