Mapping Burned Areas with Multitemporal–Multispectral Data and Probabilistic Unsupervised Learning

Abstract

The occurrence of forest fires has increased significantly in recent years across the planet. Events of this nature have resulted in the leveraging of new automated methodologies to identify and map burned areas. In this paper, we introduce a unified data-driven framework capable of mapping areas damaged by fire by integrating time series of remotely sensed multispectral images, statistical modeling, and unsupervised classification. We collect and analyze multiple remote-sensing images acquired by the Landsat-8, Sentinel-2, and Terra satellites between August–October 2020, validating our proposal with three case studies in Brazil and Bolivia whose affected regions have suffered from recurrent forest fires. Besides providing less noisy mappings, our methodology outperforms other evaluated methods in terms of average scores of 90%, 0.71, and 0.65 for overall accuracy, F1-score, and kappa coefficient, respectively. The proposed method provides spatial-adherence mappings of the burned areas whose segments match the estimates reported by the MODIS Burn Area product.

1. Introduction

Climate change is one of the most significant environmental challenges facing humanity. Occasionally, these changes in climate increase the probability of extreme events that alter rainfall and drought patterns [1,2] in several regions of the planet, besides triggering the process of natural disasters arising mainly from floods, landslides, droughts, and forest fires [3].

Extreme droughts are among the worst climatic disasters since they produce large economic and social losses [4,5]. Moreover, these phenomena also cause forest fires [6], an environmental issue that has been observed over the years, leading to severe ecological damages [7]. Frequent fires are hazardous because they affect forest regeneration, making fire a disruptive agent of this environment’s dynamics [8]. Consequently, fire mapping and its monitoring are essential to detect, mitigate, and even prevent harmful impacts [9,10].

Aiming at addressing this issue, one can apply remote sensing (RS) technologies, which gather a comprehensive collection of multitemporal data and effective imaging products. In fact, new RS-driven computational apparatuses have been proposed in the last years, enabling the extraction of high-quality data over vast monitored regions with high temporal frequency and reasonable operational costs [11]. Moreover, the development of new and more accurate data-driven tools for tracking, monitoring, and predicting fire events is leveraged in both the amount and quality of RS data, as continuously provided by a variety of sensors [12,13].

Considering the particular application of fire detection using RS data, several techniques have been created to obtain data from the Earth’s surface while mapping burned areas [14,15]. According to Birch et al. [16], the use of time series of images obtained by remote sensing is essential for studies about monitoring, management, and recovery of burned areas. In more computational terms, spectral indices as well as supervised machine-learning methods can be properly combined to produce representative detection maps for the fire. For instance, Gibson et al. [17] proposed a methodology that unified spectral indices taken from Sentinel-2, and used the random forest algorithm to map the burning severity in fire-affected areas. Similarly, Mohajane et al. [18] employed a historic forest fire inventory, composed of both local and environmental variables, to formulate distinct frequency-ratio machine learning-based methods for generating maps of fire risk areas. The creation of fire susceptibility maps was also the goal of Tien Bui et al. [19], whose authors introduced a kernel logistic regression model by taking in situ data. An extensive study by Collins et al. [20] discussed the training data requirements for adequately modeling the random forest method applied to distinguish burned areas. Likewise, Pacheco et al. [21] compared the K-nearest neighbor and random forest methods concerning the resulting burning maps using as input data images recorded by the sensor Operational Land Imager (OLI), Multispectral Instrument (MSI), and Moderate-Resolution Imaging Spectroradiometer (MODIS) on-board the satellites Landsat-8, Sentinel-2, and Terra/Aqua, respectively.

Concerning proposals that have applied spectral indices for fire mapping, recently Lasaponara et al. [10] proposed a self-organizing map-driven approach that relies on popular spectral indices, which include Normalized Difference Vegetation Index (NDVI), Normalized Burn Ratio (NBR), and Burned Area Index for Sentinel (BAIS), to enhance and identify burned areas. Mallinis et al. [22] analyzed the correlation between field survey measures and spectral indices computed from images acquired by the MSI and OLI sensors. By taking the difference Normalized Burn Ratio ($\mathsf{\Delta }$NBR) [23] and its relative $\mathsf{\Delta }$NBR (R$\mathsf{\Delta }$NBR) [24] for fire mapping applications, Cai and Wang [25] performed a statistical analysis to demonstrate that the $\mathsf{\Delta }$NBR index is more accurate and reliable than the R$\mathsf{\Delta }$NBR for burn severity classification. Gholinejad and Khesali [26] employed the NDVI and Normalized Difference Water Index (NDWI) to first identify regions subjected to burning events and then run a change-point procedure based on $\mathsf{\Delta }$NBR values to find a convenient threshold to distinguish burned areas. Mapping burned areas has not only been restricted to orbital sensors. For instance, McKenna et al. [27] adopted images recorded by unmanned aerial vehicle and the band indices excess green index, excess green index ratio, and modified excess green index to train a decision model for identifying distinct burn severity classes.

Fire detection has also been performed from the deep-learning perspective. More recently, Hu et al. [28] took deep-learning strategies for the identification of regions impacted by fire by handling uni-temporal OLI and MSI images. By employing a deep learning-based formulation, Ba et al. [29] took an empirical spectral index based-rule to automatically select samples and model a back-propagation neural network to distinguish burned areas from native vegetation on MODIS images. Similarly, Ban et al. [30] took neural networks as a primary tool to perform the mapping of wildfires. In a simple but efficient way, Matricardi et al. [31] proposed the use of probabilistic models to estimate canopy cover and fire potential in forest areas.

Despite the development of new data-driven algorithms for the detection and quantification of wildfires, the literature still lacks methods designed to classify wildfires automatically while ensuring high accuracy in the punctual mapping [30]. Indeed, for most supervised learning methods, obtaining a representative set of ground-truth data samples or a sizable number of annotated images for training are difficult tasks to accomplish in practice [32,33,34]. Another critical issue frequently found in most fire detection algorithms is that they are prone to failure when coping concurrently with temporal data, lighting conditions, and spectral bands’ complexity [35], a condition also present and addressed in other related remote sensing applications [36]. Finally, as highlighted by Khosravi et al. [37] and Barmpoutis et al. [32], there is no definitive data-driven approach that is optimal for hazard mappings, as the dependence on the study area and the availability of data can make most fire detection methods inoperative for hazard assessment.

Aiming at addressing most of the issues raised above, this paper introduces a novel data-driven multitemporal methodology for mapping burned areas by combining remote-sensing data and probabilistic models, identifying and classifying targets according to their changes measured in terms of the NBR index. These changes are turned into a mapping of burned areas by a fully unsupervised approach that applies logistic regression and spatial smoothing with Markov random fields modeling.

In summary, the introduced method comprises:

• A fully unsupervised methodology that unifies spectral indices, logistic regression, and Markov fandom fields into a spatial–temporal representation of burned areas;
• The proposed approach is capable of dealing with data acquired by sensors with distinct spatial radiometric resolutions;
• The proposed methodology retains some desirable features such as classification accuracy for images composed of different spatial resolutions, capability to handle severe imaging conditions such as cloud/shadow occurrences, and stability w.r.t. the systematic fire mapping task.
• Our computational approach is modular and, thus, flexible enough to be integrated with other learning models.

In our assessments, three case studies were carried out in regions located in Brazil and Bolivia using OLI, MSI, and MODIS images. Our choices were motivated by the fact that these study regions have suffered frequent fire events in recent years.

This article is organized as follows. Section 2 provides a brief review of fundamental concepts and techniques required for the proposed method, which is formalized in Section 3. Section 4 presents the application of the proposed method to map burned areas in three distinct areas using different sensors. Finally, Section 5 concludes this article.

2. Background

2.1. Definitions and Notation

Let ${\mathcal{I}}^{\left(t\right)}$ be an image defined on the support $\mathcal{S}\subset {\mathbb{N}}^2$ acquired at an instant t, $\mathcal{X}$ be the set of possible values at each position, and ${\mathit{x}}_i^{\left(t\right)}={\mathcal{I}}^{\left(t\right)}\left(s_i\right)$ the observation ${\mathit{x}}_i^{\left(t\right)}\in \mathcal{X}$ at position $s_i\in \mathcal{S}$ of ${\mathcal{I}}^{\left(t\right)}$. In the remote-sensing context, the components of ${\mathit{x}}_i^{\left(t\right)}$ are the values measured by the sensor, or derived features, over a specific Earth surface position.

2.2. Spectral Indexes for Burnt Area Detection

Spectral indices allow the extraction and analysis of remotely sensed data. In the face of a feature of interest, a spectral index can assist in its identification. This approach is also essential due to the impossibility of modifying orbiting imaging sensors and the difficulty of obtaining field data with the same spatial and temporal resolution [38]. Generally speaking, spectral indices are derived from algebraic operations on the attributes of $\mathcal{X}$ that characterize the behavior of ${\mathit{x}}_i^{\left(t\right)}$ at every pixel $s_i$ of ${\mathcal{I}}^{\left(t\right)}$.

Among dozens of spectral indices in the literature [39], typical examples are the NDVI [40], the NDWI [41], and the Soil Adjusted Vegetation Index (SAVI) [42].

Usually, such indices combine the information from visible (VIS), near-infrared (NIR), and shortwave-infrared (SWIR) spectral bands, being sensitive to variations in color, composition, soil moisture, and vegetation chlorophyll. The NBR [23] is a convenient index to identify burned areas.

Consider an image ${\mathcal{I}}^{\left(t\right)}$ with ${\mathrm{NIR}}_i^{\left(t\right)}$ and ${\mathrm{SWIR}}_i^{\left(t\right)}$ as components of the attribute vector ${\mathit{x}}_i^{\left(t\right)}$, which represents the target behavior at position $s_i$ with respect to the NIR and SWIR wavelengths; the image of NBR values is defined as:

$${\mathcal{I}}_{\mathrm{NBR}}^{\left(t\right)}\left(s_i\right)=\frac{{\mathrm{NIR}}_i^{\left(t\right)}-{\mathrm{SWIR}}_i^{\left(t\right)}}{{\mathrm{NIR}}_i^{\left(t\right)}+{\mathrm{SWIR}}_i^{\left(t\right)}},\mathrm{for}\mathrm{every}{s}_i\in \mathcal{S}$$

(1)

Spectral indices may be used at a single time t or as a time series. In the latter case, the use of differences between instants may help with identifying and mapping specific events. When the interest is mapping burned areas through remote-sensing images, the approaches found in the literature use only post-fire images or combine pre- and post-fire data [43]. In particular, Key and Benson [23] propose the “Delta-NBR” feature ($\mathsf{\Delta }\mathrm{NBR}$), aimed at measuring the severity of fire events on the vegetation. Such methodology comprises a characterization of fire severity through the following difference:

$$\mathsf{\Delta }\mathrm{NBR}={\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{pre}\right)}-{\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{post}\right)}$$

(2)

where ${\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{pre}\right)}$ and ${\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{post}\right)}$ are the NBR index computed at two distinct instants, before and after the fire event. Burning signs are usually characterized by values above $0.1$.

Table 1. Burning severity categorization.

Severity	$\mathsf{\Delta }$ NBR ($\times {10}^3$)
Enhanced Regrowth, High	−500 to −251
Enhanced Regrowth, Low	−250 to −101
Unburned	−100 to 99
Low Severity	100 to 269
Moderate–Low Severity	270 to 439
Moderate–High Severity	440 to 659
High Severity	660 to 1300

displays a fire severity categorization proposed by Key and Benson [23].

Although the $\mathsf{\Delta }\mathrm{NBR}$ offers good perspectives regarding the fire severity effects on the landscapes, it is important to highlight that atmospheric interference, or even variations in the limits of NIR and SWIR spectral bands, may require adjustments in the burn severity categories (Table 1). Another limitation stems from the intrinsic features of each ecosystem, which may impose significant mismatches between the severity classes and those intervals [44].

We avoid the aforementioned limitations by building a binary map that incorporates both radiometric-derived and contextual information. The former enters through a logistic regression (Section 2.3); the latter uses the Iterated Conditional Modes (ICM) algorithm (Section 2.4).

2.3. Logistic Model

Logistic regression has been widely used to estimate the probability of membership in “element-classes” in cases of binary association [45]. It acts as a class indicator provided a set of observations $\mathcal{D}=\left\{({\mathit{x}}_i,y_i)\in \mathcal{X}\times \mathcal{Y}:i=1,\dots ,m\right\}$, where $\mathcal{Y}=\left\{0,1\right\}$.

This regression belongs to the class of generalized linear models [46], in which the response $y_i$ follows a Bernoulli law with probability of success

$$f(\mathit{x};\mathbf{\theta })=\frac{1}{1+exp\{-[\mathbf{1},\mathit{x}]·{\mathbf{\theta }}^T\}}$$

(3)

where $\mathbf{\theta }=\left[{\theta }_0,{\theta }_1,\dots ,{\theta }_n\right]$ is the parameter vector with one dimension more than the input vector $\mathit{x}$ and $\mathbf{1}$ is a column of ones. The behavior of $f(\mathit{x};\mathbf{\theta })$ due to its sigmoid structure has an image defined in the range $[0,1]$.

Beyond regression purposes, the model $f(\mathit{x};\mathbf{\theta })$ may also be applied for binary classification tasks. Equation (4) represents the classification result $\mathcal{M}$ resulting from applying f on each position $s_i$ of $\mathcal{I}$:

$$\mathcal{M}\left(s_i\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}f({\mathit{x}}_i;\mathbf{\theta })<0.5\hfill \\ 1\hfill & \mathrm{if}f({\mathit{x}}_i;\mathbf{\theta })\ge 0.5\hfill \end{array}\right.$$

(4)

where the values 0 and 1 are class indicators.

The parameter $\mathbf{\theta }$ that characterizes f can be estimated by minimizing the following objective function [47]:

$$J\left(\mathbf{\theta }\right)=-{\frac{1m}{\sum }}_{i=1}^m\left[y_{i}logf({\mathit{x}}_i;\mathbf{\theta })+(1-y_i)log\left(1-f({\mathit{x}}_i;\mathbf{\theta })\right)\right]$$

(5)

Although the derivative of (5) with respect to $\mathbf{\theta }$ does not have an explicit form, its convexity grants that optimization algorithms such as the descending gradient converge to a global minimum [48].

2.4. Iterated Conditional Modes

The ICM algorithm is a classification technique that incorporates both the site evidence and the local information [49]. It is based on modeling the former with a suitable marginal distribution and the latter with a Markov random field for the underlying (unobserved) classes. It consists of solving the following maximization problem [50]:

$${\mathit{x}}_i\mathrm{is}\mathrm{assigned}\mathrm{to}{\omega }_j\iff j=\underset{k=1,\dots ,c}{\arg \max }\left\{P({\mathit{x}}_i\mid {\omega }_k)+\beta \sum _{\eta \in W}\sum _{t\in \mathcal{N}\left(i\right)}\left(1-\delta (y_i,y_t)\right)\right\}$$

(6)

where ${\omega }_k$ denotes $y_i=k$, $\beta \ge 0$ weights the influence of neighboring pixels (the context), $\eta$ is a clique (a clique of an undirected graph $\mathcal{G}=(\mathcal{A},\mathcal{B})$ is a set of vertices $\overline{\mathcal{A}}\subseteq \mathcal{A}$ such for all $i,j\in \overline{\mathcal{A}}$, there is an edge between i and j) that contains i, W is the set of cliques, $\mathcal{N}\left(i\right)$ is the coordinate set of neighboring pixels of $s_i$, and

$${\displaystyle \delta (y_a,y_b)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{y}_a=y_b\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.}$$

is the Kronecker delta function. Equation (6) is solved iteratively until a convergence criterion is met. The parameter $\beta$ is estimated, at each iteration, by maximum pseudolikelihood using the previous classification as input; see details in Frery et al. [50]. In the limit $\beta =0$, the solution of (6) is the maximum likelihood classification map, while when $\beta \to \infty$, the observations are discarded in favor of the most frequent class in the neighborhood.

3. Method Proposal

3.1. Conceptual Formalization

Figure 1 depicts a high-level overview of the proposed method for burned area detection.

3.1.1. Multitemporal NBR Computing

The red block outlines how the NBR is computed. Initially, assume there is a reference image ${\mathcal{I}}^{\left(\mathrm{Ref}\right)}$ expressing the usual behavior of targets regarding a spatial domain $\mathcal{S}$. Furthermore, let ${\mathcal{I}}^{\left(t\right)}$, for $t=1,2,\dots ,k$, be a multitemporal image series also regarding $\mathcal{S}$. Given $s_i\in \mathcal{S}$, denote ${\mathit{x}}_i^{\left(t\right)}=[x_{i1}^{\left(t\right)},\dots ,x_{in}^{\left(t\right)}]\in \mathcal{X}$, for $t\in \left\{\mathrm{Ref},1,2,\dots ,k\right\}$ as an attribute vector with measures registered at n spectral wavelengths. Using the appropriate components of ${\mathit{x}}_i^{\left(t\right)}$, it is possible to compute the NBR index at each position of $\mathcal{S}$ and then determine ${\mathcal{I}}_{\mathrm{NBR}}^{\left(t\right)}$.

3.1.2. Spatio-Temporal Representation for $\mathsf{\Delta }\mathrm{NBR}$

The $\mathsf{\Delta }\mathrm{NBR}$ comprises a derived spectral index obtained by the difference of NBR at “pre” and “post” instants, in order to distinguish burn-affected or recovered burnt vegetation from other targets. In this sense, fixing ${\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{Ref}\right)}$ as a “pre-fire event” image, we define ${\mathsf{\Delta }}_{\mathrm{NBR}}^{\left(t\right)}={\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{Ref}\right)}-{\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{t}\right)}$ for $t=1,2,\dots ,k$. We denote ${\mathsf{\Delta }}_{\mathrm{NBR}}^{\left(t\right)}\left(s_i\right)={\delta }_i^{\left(t\right)}\in \mathbb{R}$. Following the temporal order $t=1,2,\dots ,k$, the computed $\mathsf{\Delta }{\mathrm{NBR}}^{\left(t\right)}$ are combined to determine ${\mathcal{I}}_{\mathsf{\Delta }}$. In this representation, for a given $s_i$, we have ${\mathcal{I}}_{\mathsf{\Delta }}\left(s_i\right)={\mathbf{\delta }}_i=[{\delta }_i^{\left(1\right)},\dots ,{\delta }_i^{\left(k\right)}]$. These operations are represented in the “Spatio-temporal representation for $\mathsf{\Delta }\mathrm{NBR}$” block (Figure 1, green box).

3.1.3. Deviation Analysis

It is worth observing that ${\mathcal{I}}_{\mathsf{\Delta }}$ is a spatio-temporal representation for the $\mathsf{\Delta }\mathrm{NBR}$ values. Locations $s_i$, without changes due to burns or recovering, tends to show similar ${\delta }_i^{\left(t\right)}$ values independently of t. Inversely, the absence of a significant correlation between two or more components of ${\mathbf{\delta }}_i$ is evidence of a temporal change.

The “point-to-line” distance is a way of quantifying the correlation between the components of ${\mathbf{\delta }}_i$, where ${\mathbf{\delta }}_i$ is the “point” and the “line” is the geometric place where ${\delta }_i^{\left(t\right)}={\delta }_i^{(\ell )}$ with $t,\ell =1,\dots ,k$. This process is summarized in the “Deviation analysis” block (Figure 1, blue box).

Using linear algebra concepts, for a generic vector ${\mathbf{\delta }}_i$, the mentioned distance, herein denoted as $d_i$, is computed by:

$$d_i=∥{\mathbf{\delta }}_i-\frac{〈{\mathbf{\delta }}_i,\mathbf{w}〉\mathbf{w}}{{\parallel \mathbf{w}\parallel }^2}∥$$

(7)

where $\mathbf{w}=[1,1,\dots ,1]\in {\mathbb{R}}^k$ is the vector director for the line used as reference (i.e., where ${\delta }_i^{\left(t\right)}={\delta }_i^{(\ell )}$); $\parallel ·\parallel$ and $〈·,·〉$ are the norm and the inner product in an Euclidean space. The submission of ${\mathbf{\delta }}_i$ assigned to each $s_i\in \mathcal{S}$ into (7) gives place to a “deviation image” ${\mathcal{I}}_D$.

In the context of this proposal, assuming that small distances for $d_i$ represent no or irrelevant changes at $s_i$ over time, such positions should represent non-burned or non-recently-burned areas. However, as consequence of the Hughes’ phenomenon [51], the dimension of ${\mathbf{\delta }}_i\in {\mathbb{R}}^k$ plays a direct influence on the values of $d_i$, since it expresses a distance measure.

Based on this premise, an initial upper-bound value $\mathsf{\Xi }$ is defined according to $d_i$ values related to ${\mathbf{\delta }}_i$ vectors whose components are limited by a tolerance $\xi$. Such a margin acts as a parameter in the proposed method (and has a physical meaning according to the NBR spectral index, as discussed in Key and Benson [23]):

$$\mathsf{\Xi }=max\left\{d_i:\underset{j=1,\dots ,k}{max}\left\{{\delta }_i^{\left(j\right)}\right\}\le \xi \right\}$$

(8)

For convenience, let us define ${\displaystyle D\left(\xi \right)=\left\{d_i:\underset{j=1,\dots ,k}{max}\left\{{\delta }_i^{\left(j\right)}\right\}\le \xi \right\}}$.

Similarly, it is possible to define an initial lower-bound $\mathsf{\Lambda }$ for the $d_i$ values related to positions that may represent a relevant change as consequence of a burn event. A convenient definition for $\mathsf{\Lambda }$ is:

$$\mathsf{\Lambda }=\mathsf{\Xi }+\alpha gm{a}_{\mathsf{\Xi }}$$

(9)

where $gm{a}_{\mathsf{\Xi }}$ is the observed standard deviation of $D\left(\xi \right)$ and $\alpha \in {\mathbb{R}}_{+}^{*}$ is a parameter that controls the separation/transition from irrelevant-to-relevant changes.

According to the above-presented discussions, the derived information from the original multitemporal image series allows discriminating between affected (burned) and unaffected locations over the analyzed period. The upper ($\mathsf{\Xi }$)- and lower ($\mathsf{\Lambda }$)-bound values are critical elements for building a model able to express the probability of temporal changes as a consequence of a burning event.

$\mathcal{D}=\left\{(d_i,y_i)\in \mathbb{R}\times \left\{0,1\right\}:i=1,\dots ,m\right\}$ is the reference set, where $y_i=0$ or 1 when $d_i\ge \mathsf{\Lambda }$ or $d_i\le \mathsf{\Xi }$, respectively. It is worth observing that, while $y_i=1$ represents a non-burning event at $s_i$ in the analyzed period (i.e., irrelevant changes), $y_i=0$ stands for relevant changes due to burning. Consequently, we may observe that $\mathsf{\Xi }$ and $\mathsf{\Lambda }$ define a margin between irrelevant and relevant changes. The intuitive assignment “change-is-burn” is a consequence of the initial stages where the NBR and $\mathsf{\Delta }\mathrm{NBR}$ values give place to ${\mathcal{I}}_{\mathsf{\Delta }}$ and then $\mathcal{D}$.

3.1.4. Logistic Regression and MRF Modeling

With its basis in $\mathcal{D}$, a logistic regression model (Section 2.3) $f(d;\mathbf{\theta })$ is built. In the light of (4), this model classifies each position $s_i\in \mathcal{S}$ and returns a map $\mathcal{M}$ of either irrelevant or relevant changes. In the sequence, we incorporate the spatial context by using the map $\mathcal{M}$ and the probabilities interpreted from $f(d;\mathbf{\theta })$ as input for the ICM algorithm (Section 2.4). With this, we generate the regularized map $\tilde{\mathcal{M}}$. All the presented discussions, from the definition of the $\mathcal{D}$ dataset to the logistic regression modeling and ICM application, are encompassed by the “Logistic and Markov Random Field Modeling” block (Figure 1, yellow box).

In the following experiments and discussions, the proposed method is referenced as MRF+UFD, since it combines Markov Random Field (MRF) concepts into an unsupervised fire detection (UFD) approach.

3.2. Implementation Details

The previous formalization is implemented in two stages, namely: (i) data acquisition and pre-processing; and (ii) model application and processing. The code of the proposed framework is freely available at https://github.com/rogerionegri/ufd (accessed on 20 October 2022).

In the first stage, we used Python 3.8 [52] and the Application Programming Interface for Google Earth Engine (API-GEE) [53] to access and obtain the multitemporal image series. Additionally, we relied on functions provided in the Pandas [54] and GDAL [55] libraries for data manipulation. More specifically, the remote-sensing images (OLI, MSI and MODIS data; see Section 4.2) were obtained after defining a given period and region of interest. Only images with less than 50% of cloud/shadow cover were considered. Furthermore, inspired by Gómez-Chova et al. [56] and Basso et al. [57], the cloud/shadow-affected regions were removed by taking a simple masking operation, i.e., by replacing the missing regions with the median information computed from the immediately past period. Particularly, regions affected by clouds and shadows were replaced with the local-spectral trends as selected from the immediate previous period of the time series. In more mathematical terms, let ${\mathcal{C}}^{\left(t\right)}\subset \mathcal{S}$ be a boolean segmentation mask w.r.t. image ${\mathcal{I}}^{\left(t\right)}$, created to distinguish between the occurrence of cloud/shadow (i.e., True) and non-occurrence (i.e., False). From the immediate previous period, which comprises the image time-series ${\mathcal{I}}^{(t-1)},{\mathcal{I}}^{(t-2)},\dots ,{\mathcal{I}}^{(t-p)}$, the median image $\tilde{\mathcal{I}}$ is computed so as to embed the local-spectral trend into each position $s\in \mathcal{S}$. The boolean image ${\mathcal{C}}^{\left(t\right)}$ is defined by using the band/pixel quality control flags (e.g., qa_pixel; qa60; and qa for OLI, MSI and MODIS sensors, respectively). The algebraic operation $({\mathcal{I}}^{\left(t\right)}\ominus {\mathcal{C}}^{\left(t\right)})\oplus (\tilde{\mathcal{I}}\ominus \neg {\mathcal{C}}^{\left(t\right)})$ provides the overlapped segments of cloud/shadow-affected regions w.r.t. the spectral-temporal trends, expressed by $\tilde{\mathcal{I}}$. The symbols ⊖ and ⊕ act as element-wise subtraction and addition, and the “not operator’’ ¬ is applied to invert the masking structure given by ${\mathcal{C}}^{\left(t\right)}$. As a result, the cloud/shadow gaps are then filled by a regular local historic tendency without impairing the ${\mathcal{I}}_D$ representation. Preliminary tests revealed that taking p as a six-month period comprises a large enough data history to obtain a representative image collection ${\mathcal{I}}^{(t-1)},{\mathcal{I}}^{(t-2)},\dots ,{\mathcal{I}}^{(t-p)}$. From the collected image time series, we computed the median profile (i.e., $\tilde{\mathcal{I}}$), which expresses the target’s radiometric response related to the target area. Inversely, in situations of dense and frequent cloud cover, values of p assigned to short periods tend to be insufficient for extracting the proper target tendencies, and thus may impair the proposed interpolation scheme.

The reference image ${\mathcal{I}}^{\left(\mathrm{Ref}\right)}$ is the median image computed in a second past period regarding the same region of interest. We then computed the NBR spectral index for each image (i.e., ${\mathcal{I}}_{\mathrm{NBR}}^{\left(t\right)}$, $t\in \{\mathrm{Ref},1,\dots ,k\}$), followed by the calculation of the $\mathsf{\Delta }$NBR index at instant t as ${\mathsf{\Delta }}_{\mathrm{NBR}}^{\left(t\right)}={\mathcal{I}}_{\mathrm{NBR}}^{\left(\mathrm{Ref}\right)}-{\mathcal{I}}_{\mathrm{NBR}}^{\left(t\right)}$, for $t=1,\dots ,k$.

Regarding the second stage, we implemented the processes comprised of Equations (7) to (9), including the Logit model and ICM algorithm, in Interactive Data Language (IDL) v. 8.8 [58]. We verified the overall performance with values of $\alpha$ (Equation (9)) in the grid $\alpha \in \{0.01,0.05,0.1,0.25,0.5,1.0,1.5,2.0\}$ for different datasets. Although the results did not change dramatically, $\alpha =0.5$ produced the best outcomes.

4. Experiments

4.1. Experiment Overview

The following subsections present applications on the mapping of burned areas by using the proposed method. These applications include different areas and sensors. Further details are presented in Section 4.2.

For the sake of comparison, we assessed our methodology against two state-of-the-art approaches for burned area mapping: the well-established $\mathsf{\Delta }\mathrm{NBR}$ approach [59] for burn detection in multitemporal data series, referred to here as Multi-Temporal $\mathsf{\Delta }\mathrm{NBR}$ (MTDNBR), and the method introduced by Gholinejad and Khesali [26], termed here as GKM. Motivated by the minimum threshold for fire occurrence (i.e., the “low severity” category), as presented in Table 1, the MTDNBR and MRF+UFD methods assume $\mathsf{\Xi }=0.1$ as a parameter to detect burned areas. Regarding the GKM framework, a battery of preliminary tests showed that 0 and 0.2 are convenient threshold choices for the NDVI and NDWI indices. Posteriorly, the change-point procedure determines a third threshold that is taken to compute the $\mathsf{\Delta }$NBR index. More details are found in Gholinejad and Khesali [26].

The accuracies resulting from the analyzed methods were measured according to the DBFires product [60]. In summary, the DBFires database comprises a well-established governmental project that aims to detect fire focus and burned areas by integrating remotely sensed data and climatic services. According to this governmental cartographic product, burning areas are identified daily, reported as punctual occurrences, and then validated by taking real/media information. Further methodological details can be found in [61].

We employed Overall Accuracy (OA) [62], F1-score [63] and the kappa coefficient [64] as accuracy measures, as well as True/False Positive (TP; FP) and True/False Negative (TN; FN) rates [65] to evaluate and compare the performance of the examined methods. The Matthews Correlation Coefficient (MCC) [66] was also taken as a complementary measure to assess the obtained results w.r.t. the reference burned area dataset. Lastly, qualitative comparisons are carried out using the estimates from the “MCD64A1 MODIS Burned Area Monthly” [67] product.

The experiments were run on a computer with an Intel i7 processor (8-core, 3.5 GHz) and 16 GB of RAM running Ubuntu Linux version 20.04 as its operating system.

4.2. Study Areas and Data Description

In order to evaluate the proposed MRF+UFD method and compare it to the MTDNBR and GKM approaches, applications involving three study areas and different sensors were carried out. Figure 2 illustrates the study area locations. The first study area (Area 1) comprises a region of São Félix do Xingu city, state of Pará, Brazil. The second study area (Area 2) includes a region on the border between Brazil and Bolivia. The third study area (Area 3) contains a portion of Santo Antônio de Leverger city, state of Mato Grosso, Brazil. These regions have a recurrent history of forest fires between August and October.

Concerning Areas 1, 2 and 3, we took images obtained by the OLI (30 m spatial resolution of size $883\times 809$ pixels wide), MSI (10 m, $660\times 689$ pixels wide), and MODIS (500 m, $444\times 222$ pixels wide) sensors onboard the Landsat-8, Sentinel-2, and Terra satellites, respectively. The green, red, NIR, and SWIR wavelength bands were also used in order to properly run the analysis methods. Knowing the tendency of fire incidence in these regions, all images available from 1 August to 31 October 2020, limited to 50% cloud/shadow coverage over the study area, were considered in the analyses.

Table 2. Dates of the images considered in the experiments.

		August 2020	September 2020	October 2020	Sensor	Spatial Resolution
Area I	Days	13; 19	14	16	OLI	30 m
Area II		5; 10; 15; 25	4; 9; 14; 24	4; 9	MSI	10 m
Area III		4; 12; 20; 28	5; 13; 21; 29	7; 15; 23	MODIS	500 m

presents the dates of the images that meet these criteria.

Regarding the reference image (i.e., ${\mathcal{I}}^{\left(\mathrm{Ref}\right)}$, Section 3) required by the UFD+MRF method, we adopted the median image from all images acquired by the respective sensors between 1 June 2019 and 30 June 2020. A total of 9, 36, and 49 images were found in the considered period to compute the reference images for Areas 1, 2, and 3, respectively. Figure 3, Figure 4 and Figure 5 illustrate images of the study areas according to the reference image and the later instant.

Aiming at quantifying the accuracies of the methods when they are applied for mapping the burned areas over the study areas and period of analysis, the locations of fire events registered between 1 August and 31 October 2020 in the DBFires database [60] were taken as a benchmark for burned areas. Similarly, a 1 km-spaced point-wise grid covering regions not affected by fire was taken as ground-truth for non-burned areas. Figure 6 depicts the above-mentioned benchmark dataset.

As a reference for qualitative comparisons, we obtained estimates from “MCD64A1 MODIS Burned Area Monthly” [67] (500 m of spatial resolution), herein referred to as the MODIS Burned Areas (MBA) product.

Based on this product, positions/pixels in the study areas where the indication of fire shows confidence superior to 90% are taken as a “burn”; otherwise, they are considered “non-burn”. Despite the broad difference of spatial resolution between the MBA product and the OLI and MSI images, it should be stressed that this product allows for verifying the agreement between the mappings generated by the analyzed methods and the expected results. Consequently, the qualitative assessments are limited to verifying if the frequencies and spatial distribution of mapped burned areas are compatible with the expectations.

4.3. Results

Prior to applying the MTDNBR, GKM, and UFD+MRF methods, Figure 7, Figure 8 and Figure 9 show the behavior of the NBR index concerning the reference images at selected instants in each study area for data inspection purposes. With respect to Area 1 (Figure 7), it is possible to notice that, in comparison to the reference image, changes in NBR values occur progressively from August to October. In Area 2 (Figure 8), the most significant changes are only highlighted at the last instant (i.e., 9 October 2020). Conversely, the most relevant changes in Area 3 (Figure 9) occur in September and October.

Considering the evaluated study areas (Figure 2) and their corresponding sensors, all examined methods were applied to map the burned regions over the period of analysis (1 August to 31 October 2020—Table 2). The obtained results were quantitatively assessed in terms of OA, F1-score, MCC, and kappa coefficient by using the ground-truth reference dataset of burned and non-burned locations (Figure 6).

Table 3. Accuracy values summary.

Area	Method	OA	F1-Score	MCC	Kappa	Run-Time (s)
	MTDNBR	0.710	0.499	0.364	0.326	0.00248
1	GKM	0.862	0.635	0.550	0.551	1.27953
	UFD+MRF	0.885	0.655	0.598	0.588	180.44398
	MTDNBR	0.666	0.369	0.312	0.229	0.00278
2	GKM	0.924	0.532	0.561	0.500	0.67262
	UFD+MRF	0.953	0.741	0.745	0.721	66.09565
	MTDNBR	0.454	0.493	0.255	0.143	0.00172
3	GKM	0.861	0.746	0.650	0.651	0.16486
	UFD+MRF	0.874	0.750	0.671	0.667	10.47144

lists the computed accuracy values. The resulting burned area maps for each study area are presented in Figure 10, Figure 11 and Figure 12 and discussed in Section 4.4.

Regardless of the study area/adopted sensor, one can notice that the proposed methodology achieved better accuracy values in comparison to other state-of-the-art methods. While low accuracy values were assigned to the MTDNB method, more competitive performance was demonstrated by the GKM framework concerning Areas 1 and 3.

The run times for each examined method and study area are also included in Table 3. The ability to handle noisy data efficiently combined with the high-accuracy results produced by the proposed approach come with the price of computing the parameters of the Logit model, obtained by the descending gradient scheme, and the application of the ICM algorithm. Since the MTDNBR and GKM approaches only comprise thresholding-based operations, the assigned computational times for both methods are also compatible. The present UFD+MRF methodology took about 180 s, 66 s, and 10 s to processes the data sets for Areas 1, 2, and 3, respectively. Nonetheless, such a computational cost is relieved by more regular and behaved (i.e., less noisy) mappings, as shown in Figure 10, Figure 11 and Figure 12 and discussed in Section 4.4.

4.4. Discussion

Concerning the obtained results for Area 1 with OLI data, the MTDBNR method (Figure 10a) provides a map that is highly prone to false-positive labeling of burned areas in comparison to the GKM and UFD+MRF methods (Figure 10b,c). Conversely, the GKM and UFD+MRF methods deliver better discrimination between burned and non-burned areas. However, one can verify that our methodology outputs a more consistent and well-behaved mapping (i.e., less isolated pixels) in contrast to the GKM method. In terms of qualitative analysis, despite the distinct spatial resolution differences, the visual adherence of the mapping produced by the proposed approach over existing methods is evident when compared with the MBA product (Figure 10d).

Similar behavior can be observed concerning the qualitative analysis for Area 2, where 10 m spatial resolution data is used. The MTDNBR method is evidently assigned to frequent false-burn detection (Figure 11a). By inspecting the results from the GKM framework (Figure 11b) together with the ground-truth data from the MBA product (Figure 11d), one can find that the achieved mapping was unable to properly capture the burned segments. In contrast, our UFD+MRF approach (Figure 11c) delivers a fire mapping closer to the reference data.

Finally, focusing our inspections on the third study area, by using MODIS data with a spatial resolution (500 m) compatible with the reference MBA product (Figure 12d), the analyzed methods perform similarly to the previous study areas. While the MTDBNR method tends to overestimate the burning occurrences, the GKM method has frequent isolated regions/pixels (noisy aspect). Again, the output produced by the UFD+MRF approach (Figure 12c) is visually more adherent to the MBA product (Figure 12d).

It is worth stressing that the identified burned areas occurred in distinct instants of the analyzed period, as previously shown in Figure 7, Figure 8 and Figure 9. Hence, our proposal’s capability for dealing with multitemporal data in an accurate manner is also evident.

According to the obtained maps for Area 1, small burned areas were identified by the proposed approach (e.g., near coordinates 52.5W/6.7S or 65.2W/6.9S—Figure 10c) that also match the references made available in the DBFires database (Figure 6a). However, one may notice that these specific segments were not properly detected when checking the MBA product. One plausible reason behind this behavior is the different data resolutions used in each mapping (i.e., UFD+MRF/OLI—30 m; MBA/MODIS—500 m).

Similarly, regarding the burned area map obtained for Area 2, the input data spatial resolution (i.e., MS—10 m) leads to a consistent result in identifying isolated regions affected by fires (e.g., near 60.3W/15.1S and 60.5W/14.8S—Figure 11c). Indeed, these burned regions are confirmed by inspecting the DBFires reference samples (Figure 6b). On the other hand, the estimates from the MBA product indicate the fire incidence over the entire upper portion of the study area (Figure 11d), thus contradicting the DBFires records as well as the NBR variations over the evaluation period (Figure 8). Considering the accuracies in Table 3 and the above-discussed points, our approach can generate accurate burned-area estimates with better spatial resolution compared to the outputs generated by the MBA product.

Lastly, when using MODIS images as input data, the current methodology presents both similarities and differences compared to the MBA product for Area 3. For example, by taking as reference the reported locations in DBFire (Figure 6c), the lower-left (i.e., below coordinate 56.1W; 16.1S) and center-right (between coordinate 55.6–55.4W/15.9–16.0S) regions differ from the mapping produced by the MBA product. The same is not observed in the burned area mapping provided by the UFD+MRF method. A possible reason for such divergence, which also explains a few differences already discussed for Areas 1 and 2, comes from the minimum confidence level (90%) considered for the MBA product to detect burned areas. At the cost of increased inclusion errors and negative impacts on the representation consistency, identifying such regions through the MBA product can be achieved using lower confidence levels.

Both the quantitative and qualitative inspections presented above corroborate the scores outlined in Figure 13, where the proportions of true and false positives (burned areas) and negatives (non-burned areas) were computed from the ground-truth reference sets (Figure 6). The computed values reinforce that the proposed methodology is subject to inferior wrong decision proportions (i.e., false positive and negative errors) when compared against MTDNBR and GKM methods.

From the analyzed results, the proposed methodology turns out to be effective when mapping burning areas independently of the adopted multispectral sensors. Considering the operational utilization and the ground-truth data taken from the well-known MBA product, our data-driven approach ca be successfully used to generate mappings with higher spatial resolution. On the other hand, the resulting temporal estimates will also depend on the adopted sensor’s revisiting rate and cloud/shadow cover condition. However, when the main application interest lies in identifying burned areas in a given period of interest, the proposed methodology becomes more suitable.

Another aspect to observe when using our UFD+MRF approach is that it may be considered a ready-to-use core component of a monitoring system, designed to detect and quantify new burned areas as it is updated by newly acquired images. Finally, by using other spectral indices in a delta-type scheme (e.g., $\mathsf{\Delta }$NDVI and $\mathsf{\Delta }$NDWI [68,69]), the architecture described in this paper (Figure 1) can be easily customized for other application domains, as for example, deforestation and flooding mapping.

5. Conclusions

We proposed an unsupervised approach for identifying burned areas with multitemporal remote-sensing data. Spectral indices, logistic regression, and Markov random field modeling are the main concepts that support this proposal. The central idea lies in measuring the progressive deviations over time from a fixed reference image and then modeling these deviations as a probability function to map locations affected by fire in some instant of the time series. Three study cases, considering distinct regions and remote sensors, were carried out to assess the proposed method. Comparisons with two approaches found in the literature were included in the experiments.

The results showed a high adherence between the proposed method and a ground-truth reference data set. Qualitative comparisons with a coarse resolution product derived from the MODIS sensor for fire mapping purposes were also included in the analysis. While one concurrent method was unable to deliver consistent results, the other provided less regularized (noisy) mappings.

The proposed method’s robustness comes from a deviation measure computed over time, which is submitted to a probabilistic treatment to include the contextual information. On the other hand, the concurrent method (MTDNBR), as an alternative to mapping burned areas in temporal series with basis in the $\mathsf{\Delta }\mathrm{NBR}$, is highly subject to errors due to the use of fixed thresholds, as well as by the targets’ spectral variation over time.

Although the proposed method demonstrates greater computational run-time when compared to the simple threshold-based competitor, it is worth highlighting that it still comprises a light-cost algorithm.

In summary, we state as main contributions of this paper:

• A fully unsupervised learning formulation for discriminating and mapping burned areas using multispectral image time series;
• A conceptual formalization that can be conveniently customized for other remote sensing applications in addition to mapping fire-impacted areas;
• A comprehensive set of experiments and descriptive analyses using distinct sensors followed by quantitative and qualitative comparisons with existing methods from the specialized literature, thus advancing the field.

In future work, we plan to: (i) adapt and investigate the use of the proposed framework on other environmental issues, such as deforestation, flooding, oil spills, glaciers melting, etc.; (ii) analyze other probability modeling approaches (e.g., [70]); (iii) include a confidence measure concerning the resulting mappings; (iv) extend the proposed method beyond binary outputs; (v) adapt the proposed method to operate with radar images.