Evaluating BFASTMonitor Algorithm in Monitoring Deforestation Dynamics in Coniferous and Deciduous Forests with LANDSAT Time Series: A Case Study on Marmara Region, Turkey

Abstract

Time series analysis combined with remote sensing data allows for the study of abrupt changes in the environment due to significant and severe disturbances such as deforestation, agricultural activities, fires, and urban expansion, as well as gradual changes such as climate variability and forest degradation in the ecosystem. The precision of any change detection analysis is highly dependent upon its ability to separate actual changes and fluctuations on a seasonal scale. One of the efficient methods in this context is using the Breaks for Additive Seasonal and Trend (BFAST) set of algorithms. This study aims to perform a comprehensive and comparative evaluation of different Vis’ performance in forest degradation with the Landsat 8 images and BFASTMonitor approach. Through evaluation, the study also considers the potential effects of different forest types and deforestation scales in the Marmara region of Turkey. For this purpose, the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Normalized Difference Moisture Index (NDMI), and Normalized Burn Ratio (NBR) vegetation indices (VI) were selected for a comparative evaluation. The overall accuracy of VIs in deciduous forests was around 85% for NDVI, NDMI, and NBR, and 78.80% for EVI, while in coniferous forests, the overall accuracy demonstrated higher values of about 88% for NDVI, NDMI, and EVI, and 87.28% for NBR. Consequently, water-sensitive VIs that utilize shortwave infrared bands proved to be slightly more sensitive in detecting forest disturbances while chlorophyll-sensitive VIs represented lower accuracy for both forest types. Overall, all VIs faced an underestimation error in deforested area detection that was observable through negative BIAS. The results illuminate that BFASTMonitor can be considered as a tool in monitoring forest environments due to its acceptable deforestation determination capability in deciduous and coniferous forests, with slightly higher performance for small-scale deforestation patterned regions.

1. Introduction

Global deforestation and biodiversity loss due to the significant industrial and economic development have become a major concern [1,2]. Furthermore, land cover changes caused by human and/or natural disasters result in significant costs of species extinction, which contributes to climate change and biodiversity loss [3,4]. Deforestation and forest degradation are the dominant results of climate change due to human activities and urban expansion. As forest biomass encompasses a myriad amount of carbon, deforestation and forest degradation contribute to considerable human-induced carbon emissions into the atmosphere, which is one of the reasons for climate change in this century [5].

One way to address such environmental issues is to perform change detection analyses with multi-temporal remote sensing data. Time series change detection analysis of satellite images from various platforms such as optical, synthetic aperture radar, and LIDAR sensors demonstrated promising results in detecting abrupt changes such as deforestation [6]. Image difference, tasseled cap transform, principal component analysis, and other commonly used methods for change detection analysis provide simple analysis structures and algorithms. However, they require images with ideal acquisition conditions such as having similar acquisition geometries, stable sensor characteristics, and minimum seasonal variations [7]. More importantly, the change detection accuracy is heavily dependent on the threshold definitions that describe the change [8]. These limitations of the traditional change detection methods push the scientific community to conduct research in time series analyses with dense data.

Recently, Breaks for Additive Seasonal and Trend (BFAST) with high-frequency data enabled the demonstration of changes along with detecting the time of the changes and more importantly, the magnitude of them [9]. These properties make BFAST an efficient method for determining the abrupt characteristics of changes such as deforestation and forest degradation [10]. The main challenge with time series analysis is differentiating the land-cover change from other phenological vegetation changes through change detection. At this point, BFAST decomposes the changes through trend, seasonal, and remainder components iteratively and tries to separate phenological changes from land-cover changes [11].

Verbesselt et al. demonstrated that BFAST indicates accurate results in significant gradual or abrupt changes over a long period and can eradicate noise or false positive breaks [12]. However, Schultz et al. indicated the enormous number of errors related to the BFAST algorithm, such as topography, atmosphere, edge effects, data availability, and variance [13]. These factors would affect commission errors; however, the availability of data during the monitoring period is significant, as the number of observations has a prominent effect on accuracy and omission errors. BFASTMonitor was developed as an alternative to the original BFAST, which does not try to detect seasonal changes and trend changes separately. It works in a near real-time strategy, in which the deviation of new observations from stable historic data can be identified as a change. This structure of the BFASTMonitor makes it more flexible to data gaps, thus providing an advantage over the initial BFAST structure [14].

1.1. Related Works

BFAST has been widely used in the detection of biotic and abiotic environmental disturbances such as droughts, fires, and vegetation changes [15,16] in agriculture [17] and forested landscapes [18,19,20]; thus, it has been confirmed as a reliable tool [21]. Moreover, BFAST has demonstrated promising results in detecting changes over the different land-cover types, such as tropical dry forests [22,23,24], wetlands [25], wildlife nature reserves [26], cities [27], vegetation [12,22], agriculture [28], savannas [29], vegetation fire detection [30], and abandoned energy [31]. Several studies proved the capability of the BFAST algorithm and Landsat time series using a variety of spectral vegetation indices in change detection and deforestation monitoring [20,23,24,32,33]. Apart from the application of BFAST in forest change detection, Cai et al. indicated the capability of BFAST to measure the flooding of the Yangtze River [34]. Che et al. assessed the fluctuation of lake sizes in Tibet [35]. Potter used MODIS 250 m Normalize Difference Vegetation Index (NDVI) time series as an input for the BFAST method to evaluate vegetation changes specifically due to wildfires in Yellowstone National Park (USA) [36]. This study concluded that BFAST can indicate burned spots as well as represent the re-growing domains over the study area. Recently, Wu et al. implemented BFAST and Landsat NDVI time series in the West Dongting Lake region to investigate conversions of forests to other land-cover categories, from other land-cover categories to forests, and from forest to forest due to flooding or reforestation. They demonstrated that BFAST is capable of detecting multi-type forest changes with low data availability [7]. In a very current study, Gao et al. focused on the time series part of the BFAST package in detail over two different forest types with distinct seasonal characteristics and concluded that the historical period and amplitude are correlated with disturbance detection concerning tropical dry forests and temperate forests [23].

BFASTMonitor is mainly utilized in the field of deforestation and forest degradation, but there are several investigations that have applied this method in other applications, such as the detection of burned areas in a savanna landscape [37] and the estimation of forest canopy change [38]. Bueno et al. used BFASTMonitor in seven different vegetation indices, Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Soil Adjusted Vegetation Index (SAVI), Modified Soil Adjusted Vegetation Index (MSAVI), Normalized Burn Ratio (NBR), Normalized Burn Ration2 (NBR2), and Normalized Difference Moist Index (NDMI), and three distinct domains (Atlantic Forest, savanna, and semi-arid woodland) in Brazil to compare provided disturbance maps and analyze if there was any correlation among them. They concluded that having variation in input data leads to producing non-spatially correlated disturbance maps and representing site-specific sensitivity [32]. The aggregation of BFASTMonitor and Google Earth Engine was investigated by Hamunyela et al. [39]. They used a time series of Landsat NDMI dataset to evaluate their GEE BFASTMonitor implementation to detect forest disturbances in three distinct forest areas, including a humid tropical forest, dry tropical forest, and miombo woodland, and compared the results with the original BFASTMonitor with R programming language. They ascertained the high amounts of spatial and temporal agreements between the results stemmed from both the original BFASTMonitor and their GEE BFASTMonitor implementations for forest disturbance.

1.2. Paper Contributions

The studies provided in the related works section of this paper point out the effective use of the original BFAST and BFASTMonitor change detection methods, where BFAST is a good option in detecting seasonal changes but requires a stable and periodic image dataset, while BFASTMonitor is more flexible in terms of data gaps and noise but more efficient in permanent changes.

However, these studies mainly concentrated on a single VI-based comparison with other methods or used several VIs to limited geographic–climatologic extents. This study aims to perform a comprehensive and comparative evaluation of different Vis’ performances in forest degradation with the BFASTMonitor approach. Through evaluation, the study also considers the potential effects of different forest types and deforestation scales to fill this gap.

In addition, to the best of our knowledge, there is no specific study related to the performance investigation of BFASTMonitor-based analysis on forest degradation in Turkey, which exhibits a variable climatic characteristic that can be explained as a variant of the dry semi-continental Mediterranean structure.

The spotlights of this study are constructed according to three main questions:

• Does BFASTMonitor produce accurate results in the detection of disturbances over coniferous and deciduous forests in the Marmara region of Turkey?
• Does BFASTMonitor produce accurate results in both large- and small-scale deforestation?
• Among water-sensitive and chlorophyll-sensitive VIs, which ones selected in this study indicated highly accurate results?

2. Study Area and Data

2.1. Study Area

The forest types included in this study are Euxine–Colchic deciduous forests, which are categorized as temperate broadleaf and mixed forests ecoregion; and Anatolian coniferous and deciduous forests (Figure 1). Both regions are located within the Marmara region of Turkey. The topographic and climatologic characteristics of both regions are given in

Table 1. Geographical and climatic features of each study area.

Study Site	Coordinate (Lat/Long)	Climate
Site A (Selected part of Istanbul forest)	41°11′44.38″ N	Warm and hot weather during the summer, and winter is rainy and mild. Annual precipitation in the year 2020 was 820 mm and the average temperature was 14.95 °C.
	28°45′57.42″ E
Site B (Çanakkale province)	40°9′4.50″ N	Transition climate type, summer is hot and dry, winter is cold and rainy. July is the warmest month while January is the coldest month with long-term averages of 6.4 °C and 25 °C.
	27°8′58.88″ E

.

2.1.1. Euxine–Colchic Deciduous Forests

This ecoregion is located in the area south of the Black Sea, extending from the short region in southeastern Bulgaria to the northern parts of Turkey and to the east of Georgia. A specific part of this forest from the north of Turkey (site A) that includes large-scale human-induced deforestation due to the construction of the new airport in Istanbul was selected for this study [40]. According to the 2012 forest management plan, Istanbul’s total forest coverage is around 250,000 ha, with 84% of it being purely deciduous and the remaining 16% being coniferous and mixed types.

2.1.2. Anatolian Coniferous and Deciduous Forests

These forests are located in the southern part of the Marmara Sea region and the western parts of Turkey, which is surrounded by the Aegean Sea, the Marmara Sea, and Çanakkale Strait. The selected study site belongs to Çanakkale province (site B). This region suffers from small-scale disturbances in the forest due to mining activities [41]. The forests cover 54% of the whole province and mainly consist of Turkish fir (dominant), torch pine, red pine, and small forest patches of juniper and birch, especially in high altitudes. This region is selected as a representative of coniferous forests.

2.2. Landsat Data

In this study, Landsat 8 OLI, Collection 1 Level 2 satellite images were used as the input. The Level 2 processed Landsat 8 OLI image contains surface reflectance products that were geometrically and radiometrically corrected and quantized by a 10,000-scale factor. The satellite images from the Landsat 8 mission that was launched in February 2013 have a 30 m spatial resolution with a 16-day temporal resolution. In Landsat 8, the sensor was updated to Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS), which together contained eleven spectral bands and improved radiometric precision [42]. Within the context of this study, all available Landsat 8 (OLI) data from the two study areas (Path/Row, 180/031 and 181/032), that had cloud cover of less than 10% between the period of February 2013–December 2020 were obtained. Overall, 108 images that met the criteria were used for this study.

3. Methodology

The applied methodology begins with the pre-processing of the Landsat 8 OLI satellite images through cloud masking and VI calculation, then producing the forest masks for the study sites, applying the BFASTMonitor steps, and lastly, performing the accuracy assessment. A detailed flowchart of the applied methodology is presented in Figure 2.

3.1. Defining Forest and Deforestation Scale

Several studies defined a forest as an area with more than 30% of canopy cover and with a patch size greater than 0.5 ha [13]. However, there are no specific definitions of small- and large-scale deforestation. In Kalamandeen et al. study, the discussion was performed on the PRODES (Monitoramento do Desmatamento na Amazônia Legal por Satélite) program in which deforestation in the Brazilian Amazon was monitored since 1988 with Landsat images. The authors stated that the minimum mapping unit of 6.25 ha was considered for small-scale deforestation in PRODES data; thus, it does not present small-scale disturbances due to mining or agricultural practices if they are not accumulated or expanded through years. Their study also demonstrated that small forest-loss patches (smaller than 1 ha) increased by over 300% in French Guiana and Suriname in one year. They defined deforestation patches bigger than 50 ha as large-scale deforestation in their study [43]. In light of these size considerations and the potential of Landsat imagery, any deforestation mapped by the “BFASTMonitor” algorithm greater than or equal to the Landsat minimum mapping unit, which is 0.09 ha or 1 pixel, was considered for validation of the final estimation of forest cover changes in this study. The use of this approach ensures the assessment of even very small changes in the forest cover due to the expansion of agriculture and logging.

3.2. Spectral Vegetation Indices

This study evaluates the performances of two chlorophyll-sensitive VIs and two VIs belonging to water-sensitive parts of the electromagnetic spectrum (EMS) in deforestation detection. Among them, NDVI and EVI are related to healthy vegetation as they belong to the chlorophyll-absorption part of the electromagnetic spectrum, which is highly correlated with photosynthesis and plant transpiration [44]. While NDVI mostly responds to variation in the red band, EVI is responsive to NIR. Moreover, EVI tends to reduce the saturation issue in NDVI over dense canopy covers or a high amount of leaf area index [45]. NDMI is mostly used to define vegetation water content, water stress, and plant biomass changes, calculated by the application of NIR and SWIR. In this regard, it is sensitive to canopy cover and absorption by leaf moisture; therefore, it has the potential of detecting deforestation. Recently, NBR and its second version, NBR2, have been widely used in the Landsat time series due to their sensitivity to detect forest fires and burn severity. Moreover, it is proven that in the context of characterizing forest dynamics, these spectral indices have considerable abilities for various forest areas [46]. The formulas of the VIs are provided as follows:

NDVI = (NIR − Red)/(NIR + Red),

(1)

EVI = G × ((NIR − Red)/(NIR + C₁ × Red − C₂ × Blue + L)),

(2)

NDMI = (NIR − SWIR₁)/(NIR + SWIR₁),

(3)

NBR = (NIR − SWIR2)/(NIR + SWIR2),

(4)

where NIR—near-infrared (Band 5 in Landsat 8), Red—red band (band 4 in Landsat 8), Blue—blue band (band 2 in Landsat 8), SWIR1—short-wave infrared 1 band (band 6 in Landsat 8), SWIR2—short-wave infrared 2 band (band 7 in Landsat 8), G—gain factor for correction (2.5 for Landsat 8), C₁ and C₂—Coefficients of aerosol resistance term Blue (6 and 7.5 for Landsat 8, respectively).

Before VI calculation, the cloud masking process was performed with pixel QA band by using the “processLandsatBatch” function of R. Then, the above-mentioned VIs were calculated with the function “vi” that is inside the “processLandsatBatch” module of the BFAST Spatial algorithm.

3.3. Forest Mask

The primary goal of this research is to detect disturbances in forest cover. In this regard, the forest mask was created using the Landsat satellite image from the first year of the time series of this study, 2013, without the inclusion of any clouds and/or cloud shadows, as described by Schultz et al. [13]. To create the forest mask, initially, the Landsat Vegetation Continuous Field (VCF) product, which provides the regions with a canopy cover of more than 30% and a forest area of at least 0.5 ha, was used as a base map [47]. However, Landsat VCF products are available at 5-year intervals (e.g., 2005, 2010, 2015); thus, the VCF map of 2010 needed to be updated for the year 2013. Initially, according to the VCF map of 2010, the 2013-dated Landsat images were classified into the forest and non-forest classes with a supervised Random Forest (RF) classifier by using the threshold of “canopy cover greater than 30%”. Then, the forest-class patches that covered less than 0.5 ha area were considered as no class.

3.4. BFASTMonitor Implementation

Verbesselt et al. developed BFAST, which stands for Breaks for Additive Seasonal Trends, as a statistical package in R programming language [9]. It is based on a harmonic analysis model for detection changes in time series; for each pixel in the Landsat scene, it fits the best seasonal regression model with a trend component. The main objective of BFAST is to determine the moment and location of changes during the time series. In this regard, BFAST decomposes time series into three main components: harmonic, trend, and remainder or error term. Verbesselt et al. revealed the new domain of BFAST, which is BFAST Monitor [12]. The novel method is more flexible and robust, especially when data are irregular and scarce [48]. BFASTMonitor performs based on historical data observations, in which it defines the model and predicts a monitoring period, and new observations are compared to this model output. Change magnitude values are defined when there is an extreme deviation between the observed and modeled values during the monitoring period. According to the previous studies, the first-order harmonic model is fitted to the Landsat pixels observations by Equation (5) [49]:

$${\mathrm{y}}_{\mathrm{t}}=\mathsf{\alpha }+\mathsf{\gamma }\sin \left(\frac{2\mathsf{\pi }\mathrm{t}}{\mathrm{f}}+\mathsf{\delta }\right)+{\mathsf{\epsilon }}_{\mathrm{t}}$$

(5)

where y_t is the dependent variable, t is the independent variable, f is temporal frequency, α is the intercept, γ and δ are the amplitude and phase of the components, respectively, and ε_t is the noise components (residual). For detecting significant deviations among pixels during the monitoring period from the historical period, Ordinary Least Square (OLS) moving sum of residuals (MOSUM) was used (Equation (6)). The fraction of observation, defined by user-defined lag value (h) during the historical period, is as follows:

$${\mathrm{MO}}_{\mathrm{t}}=\frac{1}{\widehat{\mathsf{\sigma }}\sqrt{\mathrm{n}}}{\displaystyle \sum }_{\mathrm{s}=\mathrm{t}-\mathrm{h}+1}^{\mathrm{t}} \left({\mathrm{y}}_{\mathrm{s}}-{\hat{\mathrm{y}}}_{\mathrm{s}}\right)$$

(6)

where respectively y and ŷ are real and estimated observations, n is the number of sample observations, h is the fraction of the number of observations known as the bandwidth of MOSUM during the history period (n) [12], and $\widehat{\mathsf{\sigma }}$ is the estimator of the variance [49,50]. The signal of breakpoint would be defined as deviations from zero to beyond the 95% significance boundary.

BFASTMonitor also provides the opportunity to learn the magnitude of change M during the monitoring period. This model assumes each new pixel in a time series as a historical period and a monitoring period, which is defined by t_i ε [t₁, t_N]. The monitoring period starting time at t_n, the historical period, would be defined as t₁ < t_i < t_n, and the monitoring period as the remainder of the time, which is t_n < t_i < t_N, by calculating the median of residuals:

$$\mathrm{M}=\mathrm{median} \left\{{\mathrm{y}}_{\mathrm{t}}-{\hat{\mathrm{y}}}_{\mathrm{t}}\right\}$$

(7)

where y_t and ŷ_t are real and estimated observations, respectively [49].

The vitally important parameters that can be modified by BFASTMonitor are listed below. Although there are other parameters, these are the most significant ones [49,50]:

• Formula—regression model formula (harmonic and/or trend component);
• Order—order of the harmonic term;
• Start—starting date of the monitoring period;
• History—specification of the stable history period;
• h—bandwidth relative to the sample size in the MOSUM monitoring process, the numeric between 0–1.

Previous research demonstrated that using a first-order harmonic model with an h value of 0.25 provided the most accurate result. The h = 0.25 means that a 4-year window of data is considered in the computing of “OLS-MOSUM” statistics. In other words, by this amount of the “h-value”, only one break could be detected every 4 years [12,49].

After finishing the process over all the pixels of the images, the output of BFASTMonitor is produced as a raster brick object with three different layers: named breakpoint, magnitude, and error, respectively. The breakpoint layer is essentially the time of breaks that are detected for each pixel, the magnitude defines the median of the residuals through the monitoring period, and the error layer provides the value of “1” for pixels when an error is detected or “NA” when the method triumphs [50]. To obtain better manipulation of output layers, there is a function named the “changeMonth” used for the separation of breakpoints by year and month; moreover, a map could be created to just represent the magnitude of breaking points.

3.5. Reference Data and Validation

The accuracy assessment of BFASTMonitor to detect deforestation over both study areas was performed according to Congalton and Olofsson et al. (1991), whose method is broadly used in the concept of map validation in scientific studies [51]. The error matrix was used to compare map values associated with the breakpoint magnitude from BFAST to demonstrate the deforestation with the ground-truth information provided from reference data. Several metrics of accuracy can be calculated, including producer’s accuracy (PA), user’s accuracy (UA), overall accuracy (OA), and bias from the error matrix [52]. For the accuracy assessment, stratified random points were produced according to the breakpoint magnitudes of vegetation indices for each study area. For this purpose, 500 points for study area A and 700 points for study area B were selected, which were defined according to the total forest area of each region as one point per square kilometer.

In the second stage, these point data were labeled according to whether they represented deforestation (D) or stable (S) land cover/land use [49,53,54,55]. This ground truth information was extracted from multi-temporal, very-high-spatial-resolution (VHR) imagery by comparing the image of 2013 and late 2020 through the Google Earth platform. This platform provides satellite images with less than 5 m spatial resolution for both test sites. During the monitoring period, there were some regrowth areas with sufficient biomass to be recognized from VHR images. These areas were picked up by BFASTMonitor as deforested areas (D). However, these areas are collected as stable (S) spots for validation purposes.

In general, negative magnitude values correspond to the disturbances. However, it is worth mentioning that not all nonnegative magnitudes indicate disturbances due to forest regrowth [49]. Thus, defining the thresholds accordingly to reduce false positives is essential. In this study, the threshold for change magnitude was selected as <−0.05, where the probability of true positive disturbances was almost 50% for moderate to negative magnitudes of BFASTMonitor for all VIs [11,49,54]. The accuracy assessment was performed for all VIs (NDVI, NDMI, NBR, and EVI) according to this threshold (magnitude of change <−0.05).

4. Results

4.1. Breakpoints and Magnitude

The distribution maps of breakpoints during 2015–2020 (monitoring time) for each site are given in Figure 3 and Figure 4. In these maps, the breakpoints are labeled according to their magnitude values using a red > yellow > green color gradient scheme. Red color corresponds to slight to extreme negative breakpoint magnitudes, while yellow and green correspond to slight to moderate positive breakpoint magnitudes (Figure 3 and Figure 4). It should be considered that in the magnitude layer, all pixels are assigned a value regardless of whether any breakpoints have been detected. The yellow to green magnitude values represented extreme positive values of breaks, which are considered to be the sudden increases in the values of the vegetative areas.

According to the magnitude maps of each study area, it is observed that BFASTMonitor can detect disturbances where there is a conversion from forest to another land cover. These conversations are related to the construction of the new airport in Istanbul in study area A and the clear-cutting of forests due to the mining activities in study area B.

The breakpoints detection summary for both regions is represented in

Table 2. Percentage of all detected breakpoints in deciduous and coniferous forests.

Vegetation Index		NDVI	NDMI	NBR	EVI
Total Breakpoints	Deciduous forest	124,822	141,229	173,126	199,275
Total Pixels		275,913	275,913	275,913	275,913
Percentage		45.24%	51.19%	62.74%	72.22%
Total Breakpoints	Coniferous forest	375,605	388,647	420,210	249,464
Total Pixels		471,281	471,281	471,281	471,281
Percentage		79.70%	82.47%	89.16%	52.93%

. Based on this summary, the deforestation/breakpoints in study area A are higher than in study area B. There are approximately 557 km² of breakpoints in study area A, which is considered to be a large-scale deforestation due to the construction of a new airport in that area of northern Istanbul. Meanwhile, in study area B, small-scale deforestation occurred due to the mining activities and expansion of agriculture. Figure 5 and Figure 6 represent the time series of forest disturbance detection for all vegetation indices in study areas A and B, respectively. The solid red line indicates the time that a break was detected by the BFASTMonitor algorithm, which was on the 137th day in 2015 in study area A and the 101st day in 2017 in site B.

According to

Table 3. Percentage of the deforested area detected by each VI with respect to deciduous and coniferous forests.

Vegetation Index		NDVI	NDMI	NBR	EVI
Total Deforested Pixels	Deciduous forest	45,195	35,515	49,510	61,679
Percentage		16.38%	12.87%	17.94%	22.35%
Total Deforested Pixels	Coniferous forest	24,537	14,243	34,921	13,239
Percentage		5.20%	3.02%	7.40%	2.80%

, the detected forest disturbance percentage in the mixed deciduous forest (site A) varied between 13% and 23% for different VIs. For the coniferous forest (site B), the variation of the detection was between 3% and 7.5%. The deforested areas for both study sites are represented in Figure 7 and Figure 8, respectively.

4.2. Accuracy Assessment of Detected Changes

In this study, the OA, UA, and PA metrics, which were calculated with respect to the labeled reference points, were used to demonstrate the accuracy of the BFASTMonitor method in the detection of changes over deciduous and coniferous forest types (

Table 4. Confusion matrix for accuracy assessment. This matrix shows Vis’ accuracy for detecting deforestation by implementing the BFASTMonitor algorithm over deciduous and coniferous forests.

Deciduous Forest (A)					Coniferous Forest (B)
NDVI	D	S	Total	UA	NDVI	D	S	Total	UA
D	101	14	115	87.82%	D	126	10	136	92.64%
S	62	323	385	83.89%	S	71	493	564	87.41%
Total	163	337	500		Total	197	503	700
PA	61.96%	95.84%	OA=	84.80%	PA	63.95%	98.01%	OA=	88.42%
NDMI	D	S	Total	UA	NDMI	D	S	Total	UA
D	99	13	115	86.08%	D	119	1	120	99.16%
S	64	324	385	84.15%	S	78	502	580	86.55%
Total	163	337	500		Total	197	503	700
PA	60.73%	96.14%	OA=	84.60%	PA	60.40%	99.80%	OA=	88.71%
NBR	D	S	Total	UA	NBR	D	S	Total	UA
D	103	19	122	84.42%	D	131	23	154	85.06%
S	60	318	378	84.12%	S	66	480	546	87.91%
Total	163	337	500		Total	197	503	700
PA	63.19%	94.36%	OA=	84.20%	PA	66.49%	95.42%	OA=	87.28%
EVI	D	S	Total	UA	EVI	D	S	Total	UA
D	101	44	142	71.12%	D	118	4	122	96.72%
S	62	293	355	82.53%	S	79	499	578	86.33%
Total	163	337	500		Total	197	503	700
PA	61.96%	86.94%	OA=	78.80%	PA	59.89%	99.20%	OA=	88.14%

). For the mixed deciduous forest, apart from massive deforestation due to the construction of the new airport in Istanbul, the OA of VIs tolerated from 79% for EVI to almost 85% for NDVI. Moreover, coniferous forest regions with smaller areas of disturbances indicated an OA at about 89% for NDVI and around 88% for others. In the confusion matrix, it is indicated that the amount of PA for the deforestation detection class over the deciduous forest varied between VIs, which was around 62% for EVI and NDVI, while NBR showed 63% of PA. Meanwhile, the coniferous forest site’s percentage of PA related to the deforestation detection class was 64% for NDVI, 66% for NBR, and roughly 60% for EVI and NDMI.

In addition, the BIAS metric is calculated as the difference between the Commission Error (CE) and Omission Error (OE), where the CE is calculated as 100—UA and the OEs are calculated as 100—PA. A higher OE represents a failure to detect deforestation occurrences and is considered an underestimation of change, which is observable with negative BIAS. On the other hand, a higher CE represents false detection of deforestation and is considered an overestimation of change, which is observable with positive BIAS [56]. Calculated OE, CE, and BIAS values of VIs for each study area are represented in Figure 9. It is obvious from this figure that BIAS has raised with negative values, which means that errors produced by underestimation are higher than by overestimation of deforestation.

5. Discussions

To summarize, in this BFASTMonitor-based forest disturbance analysis, all VIs could detect the small- and large-scale deforested areas for both forest types with overall accuracies around 80% to 90%. However, there was an underestimation error associated with both deciduous and coniferous forests that was observable with a decreased PA and negative BIAS in all study site–VI combinations. The negative BIAS values could be related to the regrowth of vegetation or the succession level in the forest. It is also worth mentioning that the underestimation error of deforestation with large patches such as in site A is lower than the underestimation error of deforestation with small patches such as in site B, which can be explained by the permanent deforestation characteristics in site A due to the airport construction. In addition, the reason behind the false positive break detections could be explained by the increasing annual precipitation in the study areas, which can have an explicit impact in slight increments on VIs such as NBR and NDMI. Therefore, defining an efficient value for the threshold is essential and may be assisted with meteorological data.

In the study of Schultz et al. [46], NDMI proved to be improving the PA and balancing the errors of underestimation and overestimation in degraded areas; however, this study provided lower PA values of NDMI that are similar to other VIs [46]. The same study also suggested that the fusion of VIs can address the problem of imbalanced overestimation and underestimation errors, which is an open research area, and further studies are planned for the continuation of this study.

According to the comparison of different VIs, it has been found that the vegetation indices related to water absorption portions of the electromagnetic spectrum, such as NBR, are slightly more sensitive to the changes in forest canopy compared to vegetation indices associated with chlorophyll absorption, such as EVI. Previous research reported a lower performance of NDVI when compared to NDMI [57,58,59]. However, NDVI demonstrated reliable and similar results in both forest types in this study. A possible reason for this difference may be related to the climatologic differences between the study regions. Most of the previous studies focused on tropical and savanna ecosystems, where the amount of and changes in vegetation moisture content may be a dominant factor throughout the season, while the experimental sites of this study show significantly different climatic characteristics. It is worth mentioning that BFASTMonitor illuminates the capability for not only detection of small-scale deforestation, but also the large-scale example where several square kilometers of vegetation were cleared due to the construction of the airport. At this point, it is important to state that this study used mapping accuracy metrics; thus, results are dependent on the number of reference points, their distribution, and their accuracy, which does not directly yield to statistical significance.

There are several sources of errors related to the BFAST-based algorithms, which were investigated by Schultz et al. [13]. In this study, the primary source of error was the lack of data due to the cloud shadow and cloud cover. Cloud cover is known as a significant source of error in the investigation of time series analysis with BFAST-based algorithms. This effect was more significant in study area A; the total number of pixels for the 48 images from 2013–2020 after forest masking was initially almost 18,195,552 and 13,243,824 of them remained after cloud masking, which means that almost 27% of the input pixels were flagged as NA (Not Applicable). Nevertheless, BFASTMonitor seems to handle the irregular structure of the time series due to a lack of data, as previously stated by Awty-Carroll et al. [14].

6. Conclusions

This study evaluated the efficiency and accuracy of BFASTMonitor in detecting forest disturbances over coniferous and deciduous forest types in two study sites that represent different scales (small-scale and large-scale) of deforestation. In this study, we assessed the accuracy of VIs for coniferous and deciduous forests, which indicated that among VIs that belong to the chlorophyll-sensitive part of the electromagnetic spectrum (EMS), NDVI is more accurate than EVI in the deciduous forest region. For the coniferous forest region, overall accuracies were similar, but again, NDVI provided higher producer accuracy than the EVI. Regarding the VIs in the water absorption portion of EMS, NBR provided higher producer accuracy than NDMI for both regions while similar overall accuracies were achieved for both indices. Overall, BFASTMonitor demonstrated efficient results in characterizing deforestation and forest disturbances in both areas with different small- and large-scale deforestation, with slightly higher overall accuracy in the small-scale deforestation region. However, all VIs faced an underestimation problem in detecting deforested areas that were observed with lower producers’ accuracy and negative BIAS. Further studies are planned for implementing the fuse of several VIs to the BFASTMonitor approach with the intention to balance the overestimation and underestimation errors.