Mapping Fire Susceptibility in the Brazilian Amazon Forests Using Multitemporal Remote Sensing and Time-Varying Unsupervised Anomaly Detection

Abstract

The economic and environmental impacts of wildfires have leveraged the development of new technologies to prevent and reduce the occurrence of these devastating events. Indeed, identifying and mapping fire-susceptible areas arise as critical tasks, not only to pave the way for rapid responses to attenuate the fire spreading, but also to support emergency evacuation plans for the families affected by fire-related tragedies. Aiming at simultaneously mapping and measuring the risk of fires in the forest areas of Brazil’s Amazon, in this paper we combine multitemporal remote sensing, derivative spectral indices, and anomaly detection into a fully unsupervised methodology. We focus our analysis on recent forest fire events that occurred in the Brazilian Amazon by exploring multitemporal images acquired by both Landsat-8 Operational Land Imager and Modis sensors. We experimentally confirm that the current methodology is capable of predicting fire outbreaks immediately at posterior instants, which attests to the operational performance and applicability of our approach to preventing and mitigating the impact of fires in Brazilian forest regions.

1. Introduction

Worldwide, wildfires comprise a phenomenon of paramount importance due to their severe economic and environmental consequences [1]. It is well-known that the origin of a fire may be natural, accidental, or even criminal [2]. However, in recent decades, climate changes and the intense human activity have substantially contributed to increasing the occurrence and severity of fire-related incidents on terrestrial biomes and ecosystems [3]. Furthermore, since fires are often used to renew pasture and cultivation areas [4], such dangerous strategies can trigger fire-caused accidents, leading to widespread fires such as those that devastated the Brazilian Amazon and Pantanal biomes in 2021 [5].

Fire in tropical forests are potentially hazardous for forest conservation and regeneration, which compromise plant species richness and vegetation biomass and structure [6].

In addition, according to the Intergovernmental Panel on Climate Change (IPCC), fire-occasioned events are considered the primary source of greenhouse gas emissions, thus deliberately contributing to global warming [7]. Moreover, as stated by Brando et al. [8] and Pan et al. [9], environmental issues such as temperature increasing, precipitation reduction, and other weather-related changes in vegetation areas are responsible for the emergence of less dense humid forest environments, which are more susceptible to fires.

Although fire monitoring represents a critically important activity, mapping the risk of fire events is a very challenging task in practice [10,11] since the success in predicting future fire outbreaks is limited by the lack of understanding of what triggers and controls the dynamics of such events [8]. Aiming at circumventing this issue, one can apply remote sensing (RS) tools, which include a vast repository of multitemporal data and robust computer vision apparatus. Indeed, new RS-based technologies have been proposed in recent years, allowing for the acquisition of reliable data over large monitored areas with high temporal frequency and affordable operating costs [12]. The increase in quantity and quality of RS data, as provided by a wide variety of sensors, favors the development of new and more robust data-driven methods for tracking, monitoring, and predicting fires [13,14].

Among the approaches devoted to detecting forest fires, the assessment of fire susceptibility is a particularly effective strategy, since it allows the forecasting of potential events as well as mapping areas of risk [15]. For example, the combination of varied RS data with the potentiality of unsupervised learning (UL) has proven to be consistent, not only for predicting wildfires [16,17,18], but also for other extreme weather incidents such as flood [19,20] and landslide [21] events. Moreover, unsupervised learning classification enables the automated extraction of knowledge from the Earth’s surface while mapping spatial domains. The works carried out by Dickson et al. [22], Kamalakannan et al. [23], and Hong et al. [15] are good representatives of fire susceptibility risk models that integrate the UL paradigm and large collections of RS images. The unification of UL-based techniques and multiple RS data sources is also advocated by Pourghasemi et al. [24], where the authors state the possibility of linking fire-prone areas with an incendiary event, including forest burning.

The problem of mapping fire vulnerability areas also appears in the context of anomaly detection (AD), where a classifier is trained to model and categorize nonconforming patterns that considerably deviate from the well-behaved image instances [25,26]. For instance, Xie et al. [27] merge geostationary earth observation and polar-orbiting data into a contextual AD model that relies on diurnal temperature adjustment schemes and classic image processing filtering such as Kalman and Otsu algorithms to detect active fire outbreaks. More recently, Coca and Datcu [28] discussed the use of AD for postfire assessment, estimating the fire-affected vegetation through multispectral imagery, while Saad [29] applied neural networks (NN) to forecast fire ignition areas in Western Africa by taking meteorological features and vegetation indices to solve a supervised learning classification problem. Such a learning paradigm was also adopted by Mohammed et al. [30], where a supervised support vector machine (SVM) and other classic machine learning methods are compared for early forest fire detection using geodata. Similarly, Gholamnia et al. [31] provide a comparative study between NN, random forest, decision trees, and SVM calibrated with data collected from a moderate resolution imaging spectrometer thermal anomalies product to estimate fire zones in Iran.

Despite the development of specific algorithms for the identification and quantification of wildfires, the RS literature still lacks robust and automatic techniques [32]. In fact, the need for eligible ground-truth data or a huge amount of annotated images for training is a challenge for most supervised learning methods, which include both machine [31] and deep learning [28] approaches. Another critical issue commonly found in most fire detection methods is related to their difficulty in simultaneously dealing with temporal data, illumination conditions, and the complexity of several spectral bands to properly work. Finally, as pointed out by Khosravi et al. [33] and Barmpoutis et al. [31], there is no definitive data-driven model that is best for certain hazards, since this strongly depends on the study area and abundance of data for that specific region.

Aiming at addressing most of the issues raised above, in this paper we propose a new methodology for mapping and quantifying fire-susceptible areas. More specifically, we design a fire mapping model as a time-varying anomaly detection problem, which is solved in terms of the unsupervised learning paradigm and the derivatives of spectral indices, computed from multitemporal remote sensing data. The proposed methodology is capable of both extracting and classifying fire-related features from different localities affected in a past period. In contrast to most fire detection and prediction systems, which usually depend on current weather information and do not consider features of local fire dynamics, our approach is independent of weather or other climatic variables, so that it only takes the local land cover temporal information to generate a fire-susceptibility inference model. We attest to the accuracy performance of the proposed methodology by investigating and assessing two real cases of forest fires in Brazil.

In summary, the main contributions of this study are:

• The proposal of a fully automatic methodology for both mapping and quantifying fire-susceptible areas which relies on unsupervised anomaly detection, spectral indices differences, and satellite image time series towards better detecting patterns in complex data by learning from examples automatically.
• The applicability assessment of two anomaly detection techniques as time-varying models to select the best-performing approach for the tasks of simultaneously classifying and quantifying fire-prone areas in Brazilian Amazon rainforest portions.
• The development of an entire unsupervised training approach that integrates multiple sources of freely available satellite imagery and does not require any labeled data to generate a suitable fire detection model.
• A comparative and statistical significance analysis for each implemented method regarding areas assigned as fire against areas of true fire for two real events of wildfires in the Brazilian Amazon.

2. Theoretical Aspects and Background

2.1. Anomaly Detection as a Classification Problem

In an elementary point-of-view, a classifier is a function $F:\mathcal{X}\to \mathcal{Y}$ that assigns an element $\mathbf{x}$ from the attribute space $\mathcal{X}$ to a specific class in $\Omega = \left\{{\omega }_1,{\omega }_2,\dots ,{\omega }_c\right\}$ by setting a class indicator in a subset of natural numbers $\mathcal{Y} = \left\{1,2,\dots ,c\right\}$. Under these conditions, if $\mathbf{x}$ $\in \mathcal{X}$ and $y\in \mathcal{Y}$, $y = F\left(\mathbf{x}\right)$ indicates that $\mathbf{x}$ belongs to the class ${\omega }_y$.

The image classification task comprises the application of F on each pixel of a given image $\mathcal{I}$. More specifically, the image $\mathcal{I}$ is defined on a support $\mathcal{S}\subset {\mathbb{N}}^2$, where each pixel $s\in \mathcal{S}$ is assigned to a vector $\mathbf{x}$. Reciprocally, $\mathcal{I}\left(s\right) = \mathbf{x}$ may be used to denote the assignment between s and the vector $\mathbf{x}$, giving rise to the formula $\mathcal{C}\left(s\right) = {\omega }_y$ to express that $F\left(\mathbf{x}\right) = y$.

Basically, most image classification methods stand for different ways of defining $F:\mathcal{X}\to \mathcal{Y}$ and applying it on $\mathcal{I}$ [34]. Supervised and unsupervised are the most usual learning paradigms adopted to model F. In the supervised approach, formulating F demands taking a set of training samples whose class indicator is known in advance. On the other hand, in the unsupervised case, the classifier does not depend on prelabeled training data so that the method accounts for automatically identifying patterns, clusters, and particular relationships over the data.

The anomaly detection problem can be understood as a particular application of unsupervised classification which aims at detecting events of rare occurrences or incidents that may conflict with a set of observations [35]. One-class support vector machine [36] and isolation forest [37] are effective representatives of anomaly detection techniques. These techniques can be appropriately fit and applied to extensive datasets with acceptable computational run time.

Such approaches have been successfully applied to detect bank frauds, verify intruders into security systems, and support medical examinations [38]. In addition, anomaly detection techniques have been used as very useful tools to cope with environmental monitoring issues [26,39,40].

Next, we briefly describe both aforementioned anomaly detection techniques, which are used to model our fire mapping methodology.

2.1.1. One-Class Support Vector Machine

Support vector machine (SVM) is a very popular classification method in remote sensing applications. Its solid mathematical formulation, simple algorithmic architecture, and high generalization capability are a few of the attractive characteristics of this well-established method [41,42].

From the original conception of SVM, different variants have been proposed, for example, Laplacian [43], transductive [44], context sensitive [41,45] and “one-class” SVMs [36]. Particularly, the one-class SVM (OC-SVM) relies on the problem of quantile estimation for anomaly detection.

Conceptually, starting from a set of observations $\mathcal{I}$, the OC-SVM provides an unsupervised model that distinguishes the attribute vectors $\mathbf{x}$ as part of a set of nonanomalous elements according to a probability $\nu$ of false positive occurrences. Formally, we can write $F:\mathcal{I}\subset \mathcal{X}\to \left\{+1,-1\right\}$, where the output $+1$ implies that the data input is in $\mathcal{I}$ and $-1$ otherwise. The definitive classifier, F, is given as follows:

$$F\left(\mathbf{x}\right) = \mathrm{sgn}\left(\sum _{i = 1}^m{\alpha }_{i}K(\mathbf{x},{\mathbf{x}}_i) - b\right) ,$$

(1)

where $b = {\sum }_{j = 1}^m{\alpha }_{j}K({\mathbf{x}}_i,{\mathbf{x}}_j)$ to some ${\mathbf{x}}_i\in \mathcal{I}$, and $K(·, ·)$ is a kernel function.

The coefficients ${\alpha }_i$, $i = 1,\dots ,$ and m, are obtained by solving the following optimization problem:

$$\begin{array}{c}\underset{{\alpha }_1,\dots ,{\alpha }_m}{\min }\sum _{i,j = 1}^m{\alpha }_i{\alpha }_{j}K({\mathbf{x}}_i,{\mathbf{x}}_j)\hfill \\ s.t.\left\{\begin{array}{c}{\alpha }_i\in [0,\frac{1}{\nu m}]\hfill \\ \sum _{i = 1}^m{\alpha }_i = 1\hfill \end{array}\right.\hfill \end{array}$$

(2)

It is worth noting that the OC-SVM is parameterized by $\nu \in [0,1]$ as well as by other parameters related to the kernel function used. Additional details on kernel functions are discussed in [46].

2.1.2. Isolation Forest

Isolation forest (IF) is one of the most recent methods applied to detect anomalies. Unlike other anomaly detection methods, IF relies on the fact that it does not depend on distance measurements or density-based models to select anomalous features in the image [47]. More precisely, IF focuses on looking for isolated anomalies rather than profiling regular patterns (i.e., nonanomalies). The anomaly characterization is performed through a set of binary decision trees, denominated “isolation trees” (IT), which are components of ensemble-type classifiers, the so-called “isolation forest”.

Two main properties are the basis of the IF classification [35]: (i) anomalous instances are minority present in a dataset; (ii) the data gather attribute values very different from those observed in regular instances. In addition, the IF method can be tuned by varying three major parameters: the number $n_{it}$ of IT in the ensemble; the subsampling size $m_{sub}$ of elements randomly drawn from the analyzed dataset; and the maximum depth d allowed for the IT growth.

Formally speaking, let $\mathcal{I}$ be a dataset whose elements $\mathbf{x}$ are defined on an n-dimensional space. We denote $\tilde{\mathcal{I}}$ as the set resulting from a subsampling performed on $\mathcal{I}$ with $m_{sub}$ elements. From $\tilde{\mathcal{I}}$, a binary decision tree is constructed where each node T, and its associated data, can be subdivided into the descendants $T_l$ and $T_r$. Such division is accomplished by randomly choosing an attribute $q\in \{1,\dots ,n\}$ and a value p, which may assume $m_{sub}$ possible values. The division is then accomplished when the dataset associated with the node is greater than one; there are different values w.r.t. of the selected attribute, and the binary tree has not yet reached the maximum depth d.

We now assume that the process leading to the construction of an IT is replicated $n_{it}$ times under the condition of different $\tilde{\mathcal{I}}$, with distinct q and p selected for its construction. Therefore, as mentioned, the set of IT allows the creation of the definitive IF model.

Once the IF is properly built, a given vector $\mathbf{x}$ is classified as anomaly according to a root-to-leaf path length $h\left(\mathbf{x}\right)$ taken over each IT in the IF. Based on this concept, the following function is used to measure the anomaly level of $\mathbf{x}$:

$${\displaystyle R(\mathbf{x},m_{sub}) = 2^{-\frac{E\left(h\right(\mathbf{x}\left)\right)}{c\left(n\right)}},}$$

(3)

where $c\left(n\right) = 2ln(n-1)+\zeta -\frac{2(n-1)}{n}$; $\zeta = 0.577215664$, representing an approximation for the Euler’s constant.

Therefore, whereas $R(\mathbf{x},m_{sub})\to 1$, $\mathbf{x}$ is characterized as an anomaly; otherwise, if $R(\mathbf{x},m_{sub})<0.5$, the pattern $\mathbf{x}$ is considered regular. Finally, if $R(\mathbf{x},m_{sub})\approx 0.5$, the whole set is classified as regular.

2.2. Spectral Indexes and Burn Detection

The use of spectral indices is a well-established, robust strategy to obtain discriminative features and useful information for a desirable set of targets [48]. These include as targets the wildfire areas captured by a remotely sensed scene, so that the indices allow the assessment of burn severity of these areas while still making use of low computational resources [49,50]. Some spectral indices traditionally used in remote sensing applications are: the normalized difference vegetation index (NDVI) [51], the normalized difference water index (NDWI) [52], and the normalized burn ratio (NBR) [53].

In general, spectral indices such as NDVI and NBR combine information from visible to shortwave-infrared spectral bands (which are sensitive to variations in color, soil composition, moisture, vegetation chlorophyll, etc.) in order to capture different characteristics of soil and vegetation areas affected by fire [54]. Moreover, it is worth mentioning that due to its low-reflectance behavior, the NBR index may be unable to distinguish water bodies and burned areas. In order to circumvent this issue, the NDWI becomes a useful index to quantify areas not associated with fires.

In more mathematical terms, let $\mathcal{I}\left(s\right) = \mathbf{x}$ be an attribute vector positioned at pixel $s\in \mathcal{S}$, where $\mathbf{x}$ is formed by the reflectance intensities at green, red, near-infrared, and shortwave-infrared wavelengths bands, $x_{\mathrm{Green}}$, $x_{\mathrm{Red}}$, $x_{\mathrm{NIR}}$, and $x_{\mathrm{SWIR}}$, respectively. ${\mathcal{I}}_{\mathrm{NDVI}}$, ${\mathcal{I}}_{\mathrm{NDWI}}$, and ${\mathcal{I}}_{\mathrm{NBR}}$ establish the representation of image $\mathcal{I}$ in terms of NDVI, NDWI, and NBR spectral indices:

$$\begin{array}{c}\hfill {\mathcal{I}}_{\mathrm{NDVI}} = \frac{x_{\mathrm{NIR}}-x_{\mathrm{Red}}}{x_{\mathrm{NIR}}+x_{\mathrm{Red}}},{\mathcal{I}}_{\mathrm{NDWI}} = \frac{x_{\mathrm{Green}}-x_{\mathrm{NIR}}}{x_{\mathrm{Green}}+x_{\mathrm{NIR}}},{\mathcal{I}}_{\mathrm{NBR}} = \frac{x_{\mathrm{NIR}}-x_{\mathrm{SWIR}}}{x_{\mathrm{NIR}}+x_{\mathrm{SWIR}}}.\end{array}$$

Notice that the above-described indices can be used as time-series images, as the difference (derivative) between two images taken at different instants may allow identifying and mapping specific events. As stated by Sobrino et al. [55], most fire mapping methods that make use of remote sensing data utilize postfire or prefire images. For example, the “delta-NBR” model ($\mathsf{\Delta }\mathrm{NBR}$), introduced by Key and Benson [53], was designed to numerically assess the fire severity on vegetation areas by discriminating fire events as a function of the values obtained by:

$$\mathsf{\Delta }\mathrm{NBR} = {\mathrm{NBR}}^{\left(\mathrm{pre}\right)} - {\mathrm{NBR}}^{\left(\mathrm{pos}\right)},$$

(4)

where ${\mathrm{NBR}}^{\left(\mathrm{pre}\right)}$ and ${\mathrm{NBR}}^{\left(\mathrm{pos}\right)}$ account for the NBR index computed at two distinct instants, before and after the fire event. In our approach, burning signs are determined by values above $0.1$ [53]. Moreover, a high value of derivative $\mathsf{\Delta }\mathrm{NBR}$ implies an elevated level of severity for burning.

3. Fire Susceptibility Mapping

3.1. Computational Methodology

Now, we describe the main steps of the proposed methodology for mapping fire-prone areas by combining multitemporal remote sensing resources and unsupervised anomaly detection methods. Figure 1 shows an overview of the main stages of our approach, i.e., (i) analysis configuration and data acquisition, (ii) design of the fire-susceptibility mapping model, and (iii) generation of output results for analysis and validation purposes.

3.1.1. Multitemporal Arrangement of Remote Sensing Data

We start by describing the sequence of images taken in our study as multitemporal data, i.e., a set of remotely sensed images that were captured during three different periods.

• Prefire period: Comprises an image series taken before the fire occurred. The goal here is to capture the central tendency at each position in the study area and then use this information to generate the NBR index as a benchmark before the presence of fire (i.e., ${\mathrm{NBR}}^{\left(\mathrm{pre}\right)}$) according to the $\mathsf{\Delta }\mathrm{NBR}$ model.
• Modeling period: Covers a time interval whose data instances are exploited to identify fire events and, subsequently, build a time-varying anomaly detection model which learns the behavior of the fires immediately before they spread.
• Analysis period: Consists of the test period, where our trained anomaly detection model is applied to classify the fire-susceptible areas.

3.1.2. Spectral Mapping, $\mathsf{\Delta }\mathrm{NBR}$, and Modeling Dataset

Once the set of input images is collected, they are used to automatically build a training database, so that the unsupervised fire identification models are adjusted and generated (see Figure 2).

First, considering the image times series composed by “prefire period” data, the respective median image is computed, followed by the computation of the corresponding NBR index. As a result, the reference ${\mathrm{NBR}}^{\left(\mathrm{pre}\right)}$ is then generated. After that, concerning each instant of the “analysis period”, the NDVI, NDWI, and NBR indices are computed, where the latter is taken as ${\mathrm{NBR}}^{\left(\mathrm{pos}\right)}$, which is used to calculate the derivative spectral index $\mathsf{\Delta }\mathrm{NBR}$.

In our approach, if there is more than one image within the same time-lapse interval, we prioritize the most recent one. Moreover, areas affected by cloud/shadow occurrences are removed and recursively filled with the corresponding areas (i.e., not affected by cloud/shadow) taken from the immediately preceding instants.

The computation of $\mathsf{\Delta }\mathrm{NBR}$ is performed by taking the median of NBR, as calculated in the “prefire” period. In this process, the observed NBR values in the training period are taken as “postfire” (i.e., ${\mathrm{NBR}}^{\left(\mathrm{pos}\right)}$). From the $\mathsf{\Delta }\mathrm{NBR}$ computed at each instant, if it exceeds a prefixed threshold $\tau \ge 0.1$ (Preliminary tests show that $\tau \ge 0.1$ provides a suitable mapping of fire-affected locations.), then it is labeled as an event associated with fire occurrence. Aiming at minimizing false positives, the identified fire-affected areas are passed through the Modis burned area product [56]. This product, provides monthly maps with 500 $\mathrm{m}$ spatial resolution covering the burned areas and their respective confidence estimations. Areas that achieve $\mathsf{\Delta }\mathrm{NBR}>\tau$ and confidence above 50%, representing a burn area in the respective period, are “double-checked” as a fire-affected location. The NDVI, NDWI, and NBR spectral indices are recorded at each fire-affected position but with respect to the immediately preceding instant.

In summary, we build a fire-prone database in a fully unsupervised fashion so that potential fire-related events are discriminated in terms of their spectral signatures and derivatives. Finally, the obtained database is used to train an anomaly detection model such as OC-SVM or IF.

3.1.3. Time-Varying Unsupervised Anomaly Detection

The anomaly detection models described in Section 2.1.1 and Section 2.1.2 are then applied on each image of the analysis period immediately after the computation of the spectral indices NDVI, NDWI, and NBR. Consequently, a collection of maps is obtained so that these mappings express the anomaly occurrences at each instant.

Notice that, in our fire mapping proposal, the identification of nonanomalous locations is interpreted as areas that demonstrate fire susceptibility. As a result, a high percentage of nonanomaly occurrences w.r.t. the number of instants in the “analysis period” at each position of the region of interest gives us the definitive fire-susceptibility map.

When the percentage values tend to 100% in a certain position, it means that such a location is similar to other fire-affected areas in all instants of the analysis period. It is important to emphasize that, once a certain AD model is built from fire-affected areas, a regular behavior (i.e., nonanomaly) is assigned to elements with areas similar to the observed fire-affected regions. Reversely, if the percentage values tend to zero, the assigned location rarely behaves like a fire-affected target, thus inferring a small susceptibility to fire.

A useful characteristic of our approach is that once the definitive model is trained, it is capable of learning the local environmental conditions, which include fire occurrences. Since this is a learned behavior, the changes occasioned by the fire incidence in a certain region make the image patterns—expressed in terms of spectral indices—less susceptible to being labeled as fire at later instants.

3.2. Data Sets, Computational Resources, and Parameter Tuning

In this section, we focus on covering a few computational aspects of our methodology as well as the data and tools used during our research.

To implement our computing prototypes, we use the Python 3.8 programming language [57]. Concerning data organization and processing, we use the Numpy [58] and Pandas [59] libraries. Finally, to train and build the OC-SMV and IF anomaly detection models, we run the Scikit-Learn [60] library.

In order to access the multitemporal remote sensing images, we adopt the Google Earth Engine application programming interface (API) [61], which is compatible with the Python language. This API favors access to high-performance computing resources for processing geospatial datasets, allowing the automation of the image search process for a given period and region of interest. The images used to define the prefire, training, and analysis periods are taken from the Operational Land Imager (OLI) sensor, onboard the Landsat-8 satellite, with 30 $\mathrm{m}$ of spatial resolution, 16 days of temporal resolution, and surface reflectance information from ultra-blue to shortwave-infrared wavelengths. Concerning the ancillary Modis burned area data, used to identify the fire events and build the training dataset, it comprises a 500 $\mathrm{m}$ resolution product containing the monthly burned area per pixel according to a confidence value [56].

The images selected to define the prefire, training, and analysis periods are limited to 50% of cloud/shadow occurrences with respect to the area comprised by the region of interest (as defined in the “analysis configuration”). This percentage level was established after a battery of preliminary tests. The cloud/shadow detection is determined by using the so-called image band “pixel_qa”. More details about this procedure are found in [62]. Regarding the parameter “time-lapse” (ℓ), once our implementation focuses on Landsat-8 images, we adopt $\ell = 15$ towards considering only one image per instant in the full time series.

As mentioned in Section 2.1, the OC-SVM and IF methods require parameters tuning. Faced with the high freedom degrees associated with the process of selecting appropriate parameters for both AD methods, we apply the well-established grid search [63] procedure to calibrate the definitive models. Basically, this procedure consists of exhaustively testing, over a defined search space, the best set of parameters that ensure higher performance. For the OC-SVM method, the search space that determines the tested settings is given by $\nu \in \{{10}^{-1},{10}^{-2},\dots ,{10}^{-7}\}$ and an RBF kernel represented by $\gamma \in \{{10}^{-1},{10}^{-2},\dots ,{10}^{-7}\}$. Regarding the IF parameters, we take $n_{it}\in \{50,100,200\}$, $m_{sub}\in \{\sqrt{dim\left(\mathcal{X}\right)},100\%,75\%,50\%\}$, and $d\in \{1,2,\dots ,30\}$. In addition, this procedure is replicated into the decision rules of the data for each possible configuration, according to a 10-fold cross-validation process.

Finally, we freely provide the codes and data used to run our experiments at https://github.com/rogerionegri/fsm, accessed on 12 May 2022.

4. Study Areas and Assessment Periods

We apply the proposed methodology on two study areas and distinct epochs. Figure 3 depicts the spatial location of such regions. The first study area (Area 1) comprises an area of the Brazilian Amazon, i.e., a portion of São Félix do Xingu city, State of Pará, Brazil. The second area (Area 2) includes a portion of Cáceres city, State of Mato-Grosso, Brazil, which is another area of legal Amazon. Both areas were the scene of recent burn events with critical consequences.

Multitemporal image series registered by the OLI sensor are used as input data to the current methodology. For the sake of data exhibition, Figure 4 shows the median image computed from 1 January 2017 to 31 December 2019 regarding each study area.

The experiments comprise three distinct epochs: 2018, 2019, and 2020.

Table 1. Epochs and respective periods, expressed in terms of the reference year Y.

Pre-Fire	Modeling	Analysis	Assessment
1 January (Y-3) to	1 June (Y-1) to	1 July Y to	1 September Y to
31 December (Y-1)	31 March Y	31 August Y	31 December Y
Epochs	I	II	III
Reference year (Y)	2018	2019	2020

summarizes the epochs (I, II, and III) and the respective prefire, modeling, and analysis periods expressed as a function of the reference year Y. In addition, the “assessment period” shown in Table 1 stands for the burned areas registered by the Modis burned area product over the period from 1 September to 31 December Y, which are adopted to assess the fire susceptibility mappings. Figure 5 and Figure 6 present such reference dataset for Areas 1 and 2 at distinct epochs.

5. Experiments and Results

We now evaluate the current methodology for two study areas covering real fire incidents that occurred in the Brazilian Amazon rainforest from 2018 to 2020, as presented in Section 4. The obtained fire-susceptibility maps are then assessed by taking the burn events identified by the Modis burned area product as a benchmark (Figure 5 and Figure 6). For the sake of representation and comparison, the fire-susceptibility values are grouped into four classes: low $[0,0.25[$, mid-low $[0.25,0.5[$, mid-high $[0.5,0.75[$, and high $[0.75,1]$.

The experiments are conducted taking both OC-SVM and IF methods as anomaly detection models, so that the respective outputs are compared in order to determine the best fire mapping model for the investigated areas of the Amazon rainforest.

Figure 7 and Figure 8 comprise the obtained fire maps for Areas 1 and 2 at epochs I, II, and III, where the IF anomaly detection model was taken as part of our full classification methodology. Similarly, Figure 9 and Figure 10 present the results when employing OC-SVM to detect fire-prone regions. One can notice that the IF method provides more regularized mappings, i.e., results carrying fewer “isolated” pixels, when compared against OC-SVM. Targets with no burning event history are allocated in the low-susceptibility class (i.e., $[0,0.25[$).

Notice that these targets appear as water bodies or low biomass portions, which means that such learned behavior reveals the consistency of the current approach.

The accuracy of the obtained fire-susceptibility maps can also be visually verified and assessed from the postfire events registered by the Modis burned area product (Figure 5 and Figure 6). As previously mentioned, only the burned areas are considered in our adherence analysis, so that the expected outputs should stand for high-susceptibility values.

The plots displayed in Figure 11 and Figure 12 present the observed frequencies for each susceptibility class when using the IF and OC-SVM methods, respectively. Both plotted results are quite similar to each other. As expected, the fire-affected areas (according to the adopted benchmark dataset—Figure 5 and Figure 6) are usually assigned to high-susceptibility levels (i.e., the class $[0.75,1.0]$).

Aiming at statically assessing the significance level of the results (i.e., mid-high and high—$[0.5,1.0]$) in comparison with the low-susceptibilities frequencies (i.e., low and mid-low $[0,0.5[$, a single-tailed (unilateral) hypothesis test for population proportion [64] was computed. Under a significance of 1%, the statistical tests reveal that the higher susceptibilities occur in proportions above 75%.

With the purpose of quantitatively comparing the susceptibility maps achieved by IF and OC-SVM at each study area and epochs, we compute $2\times 2$ contingency error matrices. These reference matrices assume that the pixels of fire-affected areas, according to the Modis burn product (burn areas—Figure 5 and Figure 6), should be assigned to mid-high or high classes on the susceptibility maps, hence leading to a True Positive (TP) label. Reversely, pixels assigned to low and mid-low classes on the susceptibility maps are expected to fall within the nonburnt portions of Modis burn product images, resulting in a True Negative (TN) classification. False positives (FP) and negatives (FN) occur when expected allocations of TP or TN do not occur, respectively. For the sake of simplicity, FP and FN values are summed to indicate a global error rate, herein denoted by “False”.

Figure 13 depicts the TP, TN, and False proportions for each method and epoch. Since the number of pixels assigned to low/mid-low may change on each susceptibility map, the mentioned proportions may also change. In general, one can check that the IF method delivered slightly more False errors than OC-SVM. Nonetheless, the TP and TN proportions are quite similar for both methods regarding the results for Areas 1 and 2 in 2018 and 2019. As a result, high susceptibility regions are correctly assigned to fire-affected areas, and low-susceptibility regions are not places where fire events occur. Concerning 2020, both IF and OCSVM variants reached small False detection ratios in the period with higher occurrences of fire events (Figure 5c and Figure 6c). Moreover, the IF model achieved higher TP ratios than OCSVM in both study areas.

Additionally, we compute the F1-Score [65,66] and kappa coefficient and its variance [67] for the obtained contingency error matrices.

Table 2. Accuracy values summary. Variance of kappa scores are multiplied by ${10}^4$.

Method	Y	Area 1			Area 2
		F1-Score	Kappa	Var. Kappa	F1-Score	Kappa	Var. Kappa
IF	2018	1.00	1.00	0	1.00	0.94	356.6
	2019	0.97	0.57	1.1	0.99	0.98	18.94
	2020	0.95	0.91	1.5	0.95	0.90	1.1
OC-SVM	2018	1.00	1.00	0	1.00	0.90	226.9
	2019	0.99	0.61	1.2	1.00	0.99	35.8
	2020	0.96	0.94	3.8	0.88	0.86	2.2

lists the quality metrics for both variants of our approach in all evaluation scenarios. From the tabulated kappa values for Area 1, one can verify that the OC-SVM achieves slightly better scores than IF. In terms of F1-Score, the trained models behave very similarly. Concerning Area 2, except for 2019, the IF model produced better performance. In a generalist view, it is valid to mention that the proposed methodology equipped with IF or OCSVM AD models is capable of delivering similar results.

Finally, in order to verify the statistical significance between IF and OC-SVM variants, unilateral hypothesis tests were performed to compare the kappa and variances scores. The corresponding p-values and decisions are presented in

Table 3. Comparison of kappa values assigned to IF and OC-SVM models according to a 5% significance level.

		2018	2019	2020
Area 1	p-value	0.5	0.003	0.077
	decision	non-significant	OC-SVM	non-significant
Area 2	p-value	0.442	0.437	0.013
	decision	nonsignificant	nonsignificant	IF

. Under a 5% level, the statistical analysis indicates that there are nonsignificant differences between the analyzed methods, except for Area 1/2019 and Area 2/2020, where OC-SVM and IF deliver the most accurate results, respectively.

6. Conclusions

In this paper, we proposed a new automatic methodology for mapping fire susceptibility in Amazon forest areas using multitemporal remote sensing images, derivative spectral indices, and unsupervised anomaly detection methods. We demonstrate the effectiveness and accuracy of the introduced approach by quantitatively as well as qualitatively analyzing two case studies covering areas heavily affected by fire in the legal Amazon forest in Brazil.

From the obtained results, we showed that the current methodology was capable of assigning a high fire-susceptibility level to the real fire-damaged areas by comparing the resulting maps with reference data representing the real fire occurrences in periods posterior to those adopted to infer the susceptibility maps. Concerning the effectiveness and consistency of IF and OC-SVM, taken as time-varying unsupervised anomaly detection techniques, we found that despite IF delivering more behaved susceptibility maps, both methods produced similar results in terms of accurateness and capability of mapping fire-prone areas.

Future perspectives for this research include: (i) adapting other anomaly detection models; (ii) taking other spectral indices as part of the anomaly detection modeling; (iii) applying the proposed methodology for other study areas; (iv) extending the implementation to other multispectral sensors such as Sentinel-2 MSI and Terra/Aqua Modis; (v) evaluating the proposed methodology through adaptations to other environmental issues, such as deforestation, floods, oil spills, melting glaciers, etc. Finally, we also intend to adapt our framework to support manually labeled data, i.e., building a semisupervised strategy. A feasible attempt could integrate labeled samples together with the conception of “cluster labeling” or even using methods such as the transduction support vector machine [68]. By doing so, we could build a training database composed of labeled data with high confidence, so that this database could be used to train anomaly detection models or even classifiers to discriminate areas susceptible to fire events in the “analysis period”.