Using Spatial Features to Reduce the Impact of Seasonality for Detecting Tropical Forest Changes from Landsat Time Series

In forested areas that experience strong seasonality and are undergoing rapid land cover conversion (e.g., Brazilian savannas), the accuracy of remote sensing change detection is affected by seasonal changes that are erroneously classified as having changed. To improve the quality and consistency of regionally important forest change maps, we aim to separate process related change (for example, spectral variability due to phenology) from changes related to deforestations or fires. Seasonal models are typically used to account for seasonality, but fitting a model is difficult when there are insufficient data points in the time series. In this research, we utilize remotely sensed data and related spectral trends and the spatial context at the object level to evaluate the performance of geostatistical features to reduce the impact of seasonality from the NDVI (Normalized Difference Vegetation Index) of Landsat time series. The study area is the São Romão municipality, totaling 2440 km2, and is part of the Brazilian savannas biome. We first create image objects via multiresolution segmentation, basing the objects on the characteristics found in the first image (2003) of the 13-year time series. We intersected the objects with the NDVI images in order to extract semivariogram indices, the RVF (Ratio Variance—First lag) and AFM (Area First lag—First Maximum), and spectral information (average and standard deviation of NDVI values) to generate the time series from these features and to derive Spatio-Temporal Metrics (change and trend) to train a Random Forest (RF) algorithm. The NDVI spatial variability, captured by the AFM semivariogram index time series produced the best result, reaching 96.53% of the overall accuracy (OA) to separate no-change from forest change, while the greatest inter-class confusion occurred using the average of the NDVI values time series (OA = 63.72%). The spatial context approach we presented is a novel approach for the detection of forest change events that are subject to seasonality (and possible miss-classification of change) and mitigating the effects of forest phenology without the need for specific de-seasoning models.


Introduction
Savannas are the dominant biome in South Hemisphere, covering approximately 45% of the area of South America [1]. In Brazil, this biome consists of a mosaic of land cover types, undergoing a strong seasonality in climate, accompanied with a widespread occurrence of fires that impose environmental detection results. We used an object-based approach associated with a set of Spatio-Temporal Metrics by computing geostatistical features, specifically the semivariogram indices developed by Reference [35]. We also compared the change detection results to those provided by traditional spectral features (average and standard deviation of NDVI values). Implementation of this approach enables us to enrich object-based remote sensing methods using the spatial information provided by semivariograms that capture the spatial variability of images which are directly related to land cover changes, improving the change detection classification performance in highly seasonal and ecologically important environments.

Study Area
The study area is the São Romão municipality in the state of Minas Gerais (MG), Brazil (Figure 1a). Totaling 2440 km 2 , approximately 84% of this area is occupied by land cover classes comprising Brazilian savannas [39] (Figure 1b). The Brazilian savanna's natural vegetation consists of mixtures of grasslands, shrublands, and woodlands which are found over a gradient of vegetation density, growth form, vertical structure, and biomass [40] (Figure 1c). Riparian forests are found along the rivers and streams of the study area. As vegetation progresses from grasslands to woodlands and forests, there is an increase in cover, height, density, and the basal area of the trees present [41]. Our analysis focused on the main vegetation types found in the study area (Table 1): open grassland, open grassland with sparse shrubs, mixed grassland, shrublands, trees up to seven meters, and riparian vegetation (palm swamps and semideciduous forests) [42,43]. The climate is typically tropical (Aw) with rains concentrated in the October-May period, which characterizes the high seasonality in the region [44].  [39]; (c) Savannas natural vegetation gradients (adapted from Reference [40]).  [39]; (c) Savannas natural vegetation gradients (adapted from Reference [40]).

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 22 Table 1. The examples of the savannas land cover classes found in the study area. Landsat TM images Red-NIR, Green-SWIR 1 and Blue-Red composite.

Land Cover Classes Description
Grasslands Shrublands Savanna Woodlands

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive

Shrublands
Open grassland with sparse shrubs

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents (Figure 3):


No-change: seasonal changes due to phenological effects (natural events)  Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive

Methodology
We introduce a method for object-based change detection using geostatistical features to derive Spatio-Temporal Metrics from the time series, eliminating the effect of forest phenology. The method is divided into three main steps ( Figure 2).
We chose as the change target the Brazilian savannas land-cover classes and the following change agents ( Figure 3):

•
No-change: seasonal changes due to phenological effects (natural events) • Forest change: land cover changes such as deforestation (human-induced) and fires (both natural events and human induces)

Satellite Data and Processing
We downloaded all available cloud-free Landsat TM (Thematic Mapper) and OLI (Operational Land Imager) images between 2003 to 2016 (total of 66 images) from the free and open access USGS Earth Resources Observation and Science (EROS) Center archive ( Table 2). The data as downloaded is already processed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Processing Remote Sens. 2018, 10, 808 5 of 21 System (LEDAPS) and L8SR algorithm, respectively. Clouds and cloud shadows have been masked with the CFmask function [45].
We used the NDVI vegetation index [46] for further analysis. The advantages of using indices rather than the original spectral band observations include minimizing the soil and other background effects, reducing data dimensionality, providing a degree of standardization for comparison, and enhancing the vegetation signal [47].
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 22 Processing System (LEDAPS) and L8SR algorithm, respectively. Clouds and cloud shadows have been masked with the CFmask function [45]. We used the NDVI vegetation index [46] for further analysis. The advantages of using indices rather than the original spectral band observations include minimizing the soil and other background effects, reducing data dimensionality, providing a degree of standardization for comparison, and enhancing the vegetation signal [47].   Processing System (LEDAPS) and L8SR algorithm, respectively. Clouds and cloud shadows have been masked with the CFmask function [45]. We used the NDVI vegetation index [46] for further analysis. The advantages of using indices rather than the original spectral band observations include minimizing the soil and other background effects, reducing data dimensionality, providing a degree of standardization for comparison, and enhancing the vegetation signal [47].

Object-Based Image Analysis
Image segmentation is the division of the satellite image into groups of pixels spatially continuous and spectrally homogeneous, minimizing the within-object variability compared to the between-object variability [48]. We used the multiresolution segmentation algorithm [49] implemented in the eCognition software. eCognition is a powerful development environment for object-based image analysis designed to produce homogeneous image objects based on spatial and spectral homogeneity criteria. The eCognition technology examines image pixels not in isolation, but in the spectral and spatial contexts. It builds up a given segmentation iteratively, recognizing groups of pixels as objects, using color, shape, texture, shape, and the size of objects, as well as their context and relationships [50].
Three different choices are typically considered when defining the analysis unit (image objects) for multitemporal object-based image analysis [51]: (i) image-object overlay: image-objects are generated by segmenting one of the images in the time series. A comparison against other images is then made by simple overlay; (ii) image-object comparison: image-objects are generated by segmenting each image Remote Sens. 2018, 10, 808 7 of 21 in the time series independently; and (iii) multitemporal image-object: image-objects are generated by segmenting the entire time series together. Our study focused on the image-object overlay approach, as expected to capture the intra-object heterogeneity well.
Three key segmentation parameters (shape, compactness, and scale) control the size, shape, and spectral variation of segmented image objects [52]. We set the shape and compactness parameter to 0.1 and 0.5, respectively. The scale parameter controls the size of the image objects and sets a homogeneity threshold that determines the number of neighboring pixels that can be merged together to form an image object [53]. The scale of segmentation directly influences the size of the objects connected to the semivariogram predefined criteria (lag distance) and the minimum number of pixels inside each object necessary to generate the semivariogram makes this segmentation parameter crucial at this point. Based upon starting values found in the literature [48], we then used a "trial-and-error" approach [54] to find the appropriate scale parameter [55] that was set to 150 [37]. The image segmentation outputs were assessed based on a visual assessment of segmentation suitability. The segmentation parameters were also defined to generate sufficiently small objects to not be spatially or spectrally confused, yet large enough to allow for calculation of the semivariogram and related geostatistical features.

Geostatistical and Spectral Features
We used two semivariogram indices [35] as geostatistical features (Table 3). Inside each object, we extracted the NDVI values to generate the experimental semivariogram and calculate the semivariogram indices. The experimental semivariogram is defined as half of the average squared difference between values separated by a given lag, where this lag is a vector in both, distance, and directions [56]. It was estimated using Equation (1): where γ (h) is the estimator of the semivariance for each distance h, N(h) is the number of pairs of points separated by the distance h, Z(x) is the value of the regionalized variable in the point x and Z(x + h) is the value of the point (x + h). The graphic for spatial variance versus distance h represents the semivariogram, which allows obtaining the estimate of the variance value for different combinations of pairs of points. The semivariance functions are characterized by three parameters: sill (σ 2 ), range (ϕ), and nugget effect (τ 2 ). The sill parameter is the plateau reached by semivariance values and shows the quantity of variation explained by the spatial structure of the data. The range parameter is the distance where the semivariogram reaches the sill, showing the distance until where the data are correlated. The nugget effect is the combination of sampling errors and variations that happen in scales smaller than the distance between the sampled points [57].
When the semivariogram is calculated at the object level, an important factor for consideration is the lag distance. It should not be larger than the spatial extent of the object; on the other hand, an exceedingly small distance fails to provide a complete description of textural features. We attempted to find an optimal lag distance to ensure that sill values would provide a concise description of data variability. We fixed the number of lags as 20 pixels and the lag size equivalent to image spatial resolution (30 m), resulting in a lag distance of 600 m [37,58].
Using FETEX 2.0 [59], we calculated the semivariogram indices that describe the shape of the experimental semivariograms and, therefore, the properties that characterize the spatial patterns of each image object. In addition, fitting a model is unnecessary, which reduces both processing time and possible errors from choosing a sub-optimal model. FETEX 2.0 is a software application enabling automatic descriptive feature extraction from image-objects. The input data include a multispectral digital image or single band, such as NDVI, and a vector file in shapefile format containing the boundaries of the objects. The output file is a table containing a vector of features for every object processed. This table of numeric values describing the objects from different points of view can be then used as the input data to other classification software.
RVF is the ratio between the values of the total variance γ max_1 and the semivariance at first lag ( γ1). Since the semivariogram tends to reach the sill near the total variance, this parameter is an indicator of the relationship between the spatial correlation at long and short distances. Its value increases when the high variability at long distances and low variability at short distances occurs. This feature also presents high values for images with large primitives or periodic patterns [35]. AFM is the area between the semivariogram value in the first lag and the semivariogram. This parameter provides information about the semivariogram curvature and is also related to the variability of the data (Figure 4). processed. This table of numeric values describing the objects from different points of view can be then used as the input data to other classification software.
RVF is the ratio between the values of the total variance γ max_1 and the semivariance at first lag (γ1). Since the semivariogram tends to reach the sill near the total variance, this parameter is an indicator of the relationship between the spatial correlation at long and short distances. Its value increases when the high variability at long distances and low variability at short distances occurs. This feature also presents high values for images with large primitives or periodic patterns [35]. AFM is the area between the semivariogram value in the first lag and the semivariogram. This parameter provides information about the semivariogram curvature and is also related to the variability of the data ( Figure 4). To verify the potential of the semivariogram indices, we also used spectral features for comparative purpose. We extracted the average and standard deviation (STDV) of the NDVI values inside the objects.

Time Series
To generate the time series, we included objects intersecting savanna land cover classes in the Landsat TM image acquired in July 2003 (first image of time series) and excluded classes as water, pasture, bare soil, crops, and urban areas. We used a land cover map [39] and subsequently, we performed a post-classification editing of the imagery carried out by a skilled human interpreter using the false color composite of the corresponding image (see Figure 1b). We overlapped the savanna objects with the remaining NDVI images in order to extract the geostatistical (RVF and AFM indices) and spectral features (average and STDV) and then generated the time series ( Figure 5). To verify the potential of the semivariogram indices, we also used spectral features for comparative purpose. We extracted the average and standard deviation (STDV) of the NDVI values inside the objects.

Time Series
To generate the time series, we included objects intersecting savanna land cover classes in the Landsat TM image acquired in July 2003 (first image of time series) and excluded classes as water, pasture, bare soil, crops, and urban areas. We used a land cover map [39] and subsequently, we performed a post-classification editing of the imagery carried out by a skilled human interpreter using the false color composite of the corresponding image (see Figure 1b). We overlapped the savanna objects with the remaining NDVI images in order to extract the geostatistical (RVF and AFM indices) and spectral features (average and STDV) and then generated the time series ( Figure 5).

Spatio-Temporal Metrics (STM)
We transformed the geostatistical and spectral features time series into two summary variables: Change and Trend, adapted from Reference [16]. As described by Reference [51], change is the maximum interannual absolute difference in the time series data within the period 2003-2016 (Equation (2)). The trend is the slope of linear regression applied to the full time series (Equation (3)).
The Change variable contains information on the greatest positive or negative change between consecutive data. A high value indicates that a high magnitude change occurred or there is an influence of seasonal variation of the feature response. The Trend variable indicates the overall (increasing, decreasing, or none) trend in the amount of green vegetation across the entire period [16].

Spatio-Temporal Metrics (STM)
We transformed the geostatistical and spectral features time series into two summary variables: Change and Trend, adapted from Reference [16]. As described by Reference [51], change is the maximum interannual absolute difference in the time series data within the period 2003-2016 (Equation (2)). The trend is the slope of linear regression applied to the full time series (Equation (3)).
The Change variable contains information on the greatest positive or negative change between consecutive data. A high value indicates that a high magnitude change occurred or there is an influence of seasonal variation of the feature response. The Trend variable indicates the overall (increasing, decreasing, or none) trend in the amount of green vegetation across the entire period [16].

Random Forest Algorithm
We used the Spatio-Temporal Metrics (STM) derived from geostatistical and spectral features to train the Random Forest (RF) algorithm [60] (Figure 6). These metrics were obtained from the time series of four features (RVF, AFM, Average, and STDV) extracted from NDVI values inside objects, forming four groups of input files as described in Table 4.

Random Forest Algorithm
We used the Spatio-Temporal Metrics (STM) derived from geostatistical and spectral features to train the Random Forest (RF) algorithm [60] (Figure 6). These metrics were obtained from the time series of four features (RVF, AFM, Average, and STDV) extracted from NDVI values inside objects, forming four groups of input files as described in Table 4.   RF is an ensemble method which generates a set of individually trained decision trees and combines their results. RF is a robust non-parametric classifier and can accommodate many predictor variables [19]. The advantages of RF include effectiveness in prediction, efficient implementation on large datasets [60], and insensitivity to noise in the training data [61]. We used the Random Forest package [60] available as an R package [62].
Two parameters need to be set in order to produce the random forest trees: the number of decision trees to be generated (Ntree) and the number of variables to be selected and tested for the best split when growing the trees (Mtry) [63]. The Gini criterion is used to select the split with the lowest impurity (highest homogeneity) at each node. For each tree, the predicted class for each observation is obtained and the class with the maximum number of votes among the trees is the predicted class of an observation [64].
The Ntree required to maintain a given level of accuracy has been assessed by several authors, and the minimum number of trees for optimal classification appears to be between 100 [65] and 300 [66], with the majority of studies setting the Ntree value to 500 because the errors stabilize before this number of classification trees is achieved [65]. Five hundred trees were grown for each classification and the number of features (Mtry = 1) was left at its default values.

Change Detection Evaluation
Evaluation of change detection outcomes is communicated using a confusion matrix [67] and its accuracy measures: overall accuracy (OA), producer accuracy (PA), user accuracy (UA), and the kappa coefficient (k). A dataset of 1000 objects well-distributed over the study area, with 500 samples per class (Table 5), was assigned based on visual interpretation (Figure 7a). The objects of no-change were those where no changes were found during the entire analyzed period 2003-2016 ( Figure 7b). Otherwise, regardless of the date when changes have occurred, the objects undergoing a deforestation or fires were selected and identified as forest change (Figure 7c). Seventy percent of the samples were randomly chosen as training samples, while the rest were used as validation samples. We used the class-stratified random sampling design that provides the option to increase the sample size in classes that occupy a small proportion of area and is one of the easier designs to implement [68]. The evaluated performances were compared among the four input features (average, standard deviation, RVF, and AFM) through the confusion matrix measures. We also generated a validation map to identify the savannas land cover changes (deforestation and fires) by visual interpretation to spatialize the areas classified as correct and incorrect. We used as reference data, a land cover map [39] which described qualitative and quantitative information of savannas vegetation type in the state of Minas Gerais, Brazil. The evaluated performances were compared among the four input features (average, standard deviation, RVF, and AFM) through the confusion matrix measures. We also generated a validation map to identify the savannas land cover changes (deforestation and fires) by visual interpretation to spatialize the areas classified as correct and incorrect. We used as reference data, a land cover map [39] which described qualitative and quantitative information of savannas vegetation type in the state of Minas Gerais, Brazil.

Semivariogram Analysis
We analyzed the temporal-spatial variability of NDVI data comparing the semivariogram shape and indices. The semivariograms reached the sill (σ 2 ) within the calculated distance (600 m), indicating that their spatial extents were sufficiently large to encompass the entire spatial variability. From no-change transitions (landscape vegetation cover that had the same land cover classes in both times, encompassing seasonal variations due to vegetation phenology), NDVI semivariograms had similar behaviors, presenting sill values around 0.0010 (Figure 8a,b). In the presence of changes (that is, fires), NDVI semivariograms considerably increased, reaching sill values around 0.0060 ( Figure  8c). Note how the shape and sill parameter of no-change areas remained constant, but changed markedly in areas with human activities (Figure 8d).

Semivariogram Analysis
We analyzed the temporal-spatial variability of NDVI data comparing the semivariogram shape and indices. The semivariograms reached the sill (σ 2 ) within the calculated distance (600 m), indicating that their spatial extents were sufficiently large to encompass the entire spatial variability. From no-change transitions (landscape vegetation cover that had the same land cover classes in both times, encompassing seasonal variations due to vegetation phenology), NDVI semivariograms had similar behaviors, presenting sill values around 0.0010 (Figure 8a,b). In the presence of changes (that is, fires), NDVI semivariograms considerably increased, reaching sill values around 0.0060 (Figure 8c). Note how the shape and sill parameter of no-change areas remained constant, but changed markedly in areas with human activities (Figure 8d).
These results indicate a very clear trend whereby the overall variability and shape, the semivariogram increased when land cover changes occurred, and were similar when the area had not undergone any land cover change. These results demonstrate that the semivariograms derived from NDVI images, is not affected by phenological changes (remain constant), while also demonstrating the potential of detecting land cover changes (notably increase) due to the mosaic of types of landscape with low NDVI and high NDVI, increases the NDVI spatial variability.
indicating that their spatial extents were sufficiently large to encompass the entire spatial variability. From no-change transitions (landscape vegetation cover that had the same land cover classes in both times, encompassing seasonal variations due to vegetation phenology), NDVI semivariograms had similar behaviors, presenting sill values around 0.0010 (Figure 8a,b). In the presence of changes (that is, fires), NDVI semivariograms considerably increased, reaching sill values around 0.0060 ( Figure  8c). Note how the shape and sill parameter of no-change areas remained constant, but changed markedly in areas with human activities (Figure 8d).

Change Detection
We used the Change and Trend Spatio-Temporal Metrics derived from geostatistical and spectral features time series as input data to train the random forest (RF) algorithm. The classification results ( Figure 9) achieved overall accuracies from 63.7% to 96.5%, with user's accuracies ranging from 61.5% to 98.2% and producer's accuracies from 58.9% to 97.7% ( Table 6). The AFM index provided the best results, reaching 96.5% of the overall accuracy and a 0.9 kappa index. The standard deviation also produced good results, with a kappa index reaching 0.8. The greatest inter-class confusion occurred using the average of the NDVI values and RVF index, revealing low values of overall, producer and user accuracies and kappa index. As a result, the feature average and the semivariogram index RVF were not able to accurately capture the differences between forest change and seasonal changes. These features are sensitive to phenological variations, resulting in the false detection of LULC changes, overestimating the areas undergoing fire or deforestation (Figure 10a,c). The standard deviation spectral feature presented good accuracy measures (overall accuracy reaching 88.3%), however, the map generated was not more accurate (Figure 10b). The AFM semivariogram index provided the best result. This index is directly influenced by the sill (σ 2 ) semivariogram parameter (overall variability), which is sensitive to land cover changes and is not affected by changes caused by vegetation seasonality, providing good classification results (Figure 10d).  As a result, the feature average and the semivariogram index RVF were not able to accurately capture the differences between forest change and seasonal changes. These features are sensitive to phenological variations, resulting in the false detection of LULC changes, overestimating the areas undergoing fire or deforestation (Figure 10a,c). The standard deviation spectral feature presented good accuracy measures (overall accuracy reaching 88.3%), however, the map generated was not more accurate (Figure 10b). The AFM semivariogram index provided the best result. This index is directly influenced by the sill (σ 2 ) semivariogram parameter (overall variability), which is sensitive to land cover changes and is not affected by changes caused by vegetation seasonality, providing good classification results (Figure 10d).

NDVI Spatial Variability
To capture the NDVI variability inside the objects we used the semivariogram, that is a graphical representation of the spatial variability in a given dataset [69], provided by the sill (σ 2 ), range (φ), and

NDVI Spatial Variability
To capture the NDVI variability inside the objects we used the semivariogram, that is a graphical representation of the spatial variability in a given dataset [69], provided by the sill (σ 2 ), range (ϕ), and nugget effect (τ 2 ) parameters. The semivariogram presents the ability to depict landscape spatial heterogeneity [70]. The differences of NDVI spatial variability between landscapes are mainly explained by the type of landscape, where the image spatial variability (sill-σ 2 ) increases considerably from homogeneous to heterogeneous areas.
The semivariogram shape and metrics generated from NDVI values, derived from Landsat TM images are useful in detecting deforestation in areas covered by savanna vegetation [62]. In deforested areas, the landscape change caused spatial variations that were quantified by the semivariogram metrics of range, sill, and shape. The range and sill were the two most important and complementary metrics. Both metrics increased their values after deforestation and remain similar if the land cover had not been changed.
Here, instead of using the semivariogram parameters to quantify the NDVI spatial variability we used the RVF and AFM semivariogram indices as they synthesize the most relevant information about the shape of the semivariogram. In addition, fitting a model is unnecessary, which shortens the processing time and reduces the errors generated by choosing a sub-optimal model [35]. From the wet to dry season, the NDVI spatial variability is almost constant and in the presence of land cover changes, such as deforestation and fires, the NDVI variability highly increased ( Figure 11). The preserved savannas land cover classes (both in wet and dry seasons) have low NDVI variability compared to disturbed savannas land cover classes. These high values are explained by the presence of high and low NDVI values in the same object, as land cover changes exceed seasonal variability. The total variance of NDVI values increases considerably from homogeneous to heterogeneous landscapes. In homogeneous objects as preserved savannas, the NDVI semivariance among the pixels is low, thus, the RVF index is also low because the ymax (~0.0013) is almost similar with the semivariance in the first lag (~0.0005). The AFM is also small (0.011) (Figure 12a). In heterogeneous objects, such as savannas undergoing a deforestation or fires, the RVF and AFM indices considerably increase, ranging from 3.33 to 8.38 and 0.011 to 0.13, respectively (Figure 12b). The total variance of NDVI values increases considerably from homogeneous to heterogeneous landscapes. In homogeneous objects as preserved savannas, the NDVI semivariance among the pixels is low, thus, the RVF index is also low because the ymax (~0.0013) is almost similar with the semivariance in the first lag (~0.0005). The AFM is also small (0.011) (Figure 12a). In heterogeneous objects, such as savannas undergoing a deforestation or fires, the RVF and AFM indices considerably increase, ranging from 3.33 to 8.38 and 0.011 to 0.13, respectively (Figure 12b). Figure 11. The example of NDVI spatial variability of objects under a seasonal change and deforestation land cover change.
The total variance of NDVI values increases considerably from homogeneous to heterogeneous landscapes. In homogeneous objects as preserved savannas, the NDVI semivariance among the pixels is low, thus, the RVF index is also low because the ymax (~0.0013) is almost similar with the semivariance in the first lag (~0.0005). The AFM is also small (0.011) (Figure 12a). In heterogeneous objects, such as savannas undergoing a deforestation or fires, the RVF and AFM indices considerably increase, ranging from 3.33 to 8.38 and 0.011 to 0.13, respectively (Figure 12b). Similar patterns were found by Reference [37] exploring and evaluating the performance of semivariogram indices derived from NDVI images in an object-based approach to detect land cover changes. The image's spatial variability changed considerably from native vegetation (pre-disaster image) to flooded areas (post-disaster image). The flooded areas had a low overall variability compared to native vegetation, making the detection of change using the semivariogram possible.

Seasonal Effects on Change Detection
We hypothesized that the semivariogram indices (RVF and AFM) are not affected by vegetation seasonality due to the NDVI variability remaining constant in the presence of seasonal changes and substantially increases in the presence of land cover changes, such as deforestation and fires (both natural events and human-induced).
Remote sensing of vegetation has an implied link to vegetation phenology. Phenology is an environmental process that operates independently of land cover or land use change. The ability to use remote sensing images to detect changes in seasonal areas relies on the capacity to remove the Similar patterns were found by Reference [37] exploring and evaluating the performance of semivariogram indices derived from NDVI images in an object-based approach to detect land cover changes. The image's spatial variability changed considerably from native vegetation (pre-disaster image) to flooded areas (post-disaster image). The flooded areas had a low overall variability compared to native vegetation, making the detection of change using the semivariogram possible.

Seasonal Effects on Change Detection
We hypothesized that the semivariogram indices (RVF and AFM) are not affected by vegetation seasonality due to the NDVI variability remaining constant in the presence of seasonal changes and substantially increases in the presence of land cover changes, such as deforestation and fires (both natural events and human-induced).
Remote sensing of vegetation has an implied link to vegetation phenology. Phenology is an environmental process that operates independently of land cover or land use change. The ability to use remote sensing images to detect changes in seasonal areas relies on the capacity to remove the effects of seasonal variations, while still identifying changes caused by deforestation, urbanization, floods, and fires. To avoid the classification of seasonal effects as land cover change, time series have been used to derive the seasonal patterns when mapping deforestation. For instance, Verbesselt et.al developed a powerful change detection approach for time series, using 16-day MODIS NDVI images, involving the detection and characterization of Breaks For Additive Seasonal and Trend (BFAST) that integrates the iterative decomposition of time series into trend, seasonal, and noise components with methods for detecting changes, without the need to select a reference period, set a threshold, or define a change trajectory [15].
Rather than using the BFAST seasonal model to account for seasonality in the image time series, Hamunyela et.al researches used the spatial context in a pixel-based approach using Landsat data assessing how the size of spatial window influences the deforestation detection and related computational time [12]. They demonstrated that, in dry tropical forests, deforestation events are detected much earlier when the spatial context approach is used than when a seasonal model is used. Smaller spatial windows achieved higher overall spatial accuracy than larger windows and increasing the size of the spatial window resulted in a linear increase in the computational time per pixel.
The window size is an important aspect of the spatial context method since it influences the resultant spatial information due to the amount of variance included [71,72]. The spatial moving window needs to be sufficiently larger than the size of the deforestation event to avoid smoothing out the disturbance impact from the time series, and also for it to capture local variations [12].
Here, we used the spatial context (semivariogram indices) in an object-based approach, without the need of use of the spatial window, basing the objects on the characteristics found in the first image of the time series. The time series of the AFM semivariogram index provided the best separability (global accuracy = 96.53%) between objects with seasonal changes and objects undergoing a forest change. Similar results were found by Reference [73], assessing the potential of individual geostatistical features derived from bitemporal NDVI images to accurately detect land cover changes. They demonstrated that the accuracies of the sill parameter and AFM semivariogram indices were higher than those of the other indices, which were not affected by vegetation seasonality because they represent the information provided by the entire semivariogram rather than specific parts of it.

Method Limitations
Despite the good performance of the proposed method, in this section, two limitations to the approach used in this study are elaborated upon, including the quantification of areas changed and the timing of forest changes: (1) Forest change quantification: The method can identify the objects undergoing land cover changes without being impacted by seasonal variations; however, it is not possible to accurately calculate the changed area. Due to use of an object-based approach, some changes do not cover the entire object because the forest change class comprised both small changes reaching up to 20% of the objects and large changes ranging from 50% to 100% of the object area. This limitation is related to the initial image segmentation performed on the first image of the time series (image object overlay) rather than segmenting the entire time series (multitemporal image object). We used this technique because it provides a more meaningful framework for spatial measures [51], capable of capturing the intra-object heterogeneity [74]. Further research into how the proportion of pixels changed inside the objects (real area of forest change) affect the resultant accuracies remain outstanding. (2) Timing of changes: The method allows for the discrimination between forest change and no-change classes, however, it is not possible to indicate when exactly the forest changes have occurred. Although the timing of the change is an important feature of forest monitoring, we did not attempt to model the time of the change, as our main focus is not hinged on the capacity to date forest changes, rather we focus on the advantages of using variability measures as input data to train the RF algorithm.

Conclusions
We have developed a method to mitigate the presence of phenological effects from time series using NDVI trends to detect changes over a savannah forest ecosystem. The approach combining spatial and spectral information and allowed for the detection of changes without being impacted by seasonal variations. We demonstrated that NDVI spatial variability, captured by AFM semivariogram index, is not affected by vegetation seasonality and, therefore, can produce time series that accurately differentiate forest changes from seasonal changes, resulting in fewer classification errors. Our proposed method is simple and accurate, efficiently detecting changes while eliminating the effects of phenology.
The spatial context approach we presented in this paper is a novel and useful object-based remote sensing method for the detection of forest change events using Landsat NDVI data in areas where forests exhibit strong seasonality (that is, in Brazilian seasonal biomes). Our results demonstrate that the semivariograms derived from NDVI images are not affected by phenological changes (remain largely constant). The spectral values may differ over time within a given year, but absent actual change (for example, fire or deforestation), the amount and nature of the variability present at the object level remains similar. This method addresses the challenge of accurately detecting deforestation and fires without erroneously labeling seasonal phenological differences as changes without the need to use de-seasoning models.