Comparative Analysis of Non-Negative Matrix Factorization in Fire Susceptibility Mapping: A Case Study of Semi-Mediterranean and Semi-Arid

1. Introduction

Forest fires’ rising frequency and severity have emerged globally as a critical issue fueled by natural factors and human activities. Extreme weather events, shifts in land use, and rapid urban development are crucial contributors to this growing problem, exacerbating the risk of wildfires [1,2,3,4]. In 2021, wildfires caused the global loss of approximately 7.02 million hectares of tree cover [5]. Climate change has exacerbated these incidents through rising temperatures, decreased precipitation, and more frequent droughts. In contrast, human activities such as deforestation and land-use changes have further contributed to this escalation. These effects are particularly severe in vulnerable regions like the Semi-Mediterranean (SM) and Semi-Arid (SA) forests of the Kurdo-Zagrosian mountains in western Iran and northern Iraq [6,7,8]. The consequences of forest fires extend beyond immediate environmental damage; they contribute to global carbon emissions and disrupt ecosystems, leading to soil degradation, erosion, loss of microfauna and flora, and the deterioration of water quality [9].

In response to these challenges, early identification of fire-prone areas has become increasingly crucial for effective risk mitigation. Forest Fire Susceptibility Mapping (FFSM) is vital in identifying fire-prone areas or areas at high risk of wildfires, enabling better ecological management and biodiversity conservation [10,11,12]. Such map-ping efforts are significant for regions like the Kurdo-Zagrosian mountains, where human activity and climate change are constantly threatened. By accurately mapping fire-prone zones, FFSMs help decision-makers prioritize resources, reduce vulnerability-ties, and implement preventive measures, ensuring long-term ecosystem health. These efforts rely on factors such as vegetation cover, temperature, humidity, rainfall, wind speed, proximity to roads and water bodies, elevation, slope, land use patterns, and the Topographic Wetness Index (TWI) [13,14,15,16].

Technological advances, particularly in combining Geographic Information Systems (GIS) with Remote Sensing (RS) data/techniques, have revolutionized the development of FFSMs, enabling more precise identification of fire-prone areas. Additionally, ensemble models, which combine the strengths of multiple algorithms, have shown superior accuracy in predicting fire susceptibility [17,18,19]. Studies have demonstrated that these models applied across regions like Gangwon-do in South Ko-rea and the Chaloos Rood watershed in Iran have yielded highly accurate predictions by integrating topography, climate, and human activity [20,21]. Beyond this, integrating geospatial and RS data has allowed researchers to account for various factors that influence wildfire susceptibility, such as temperature, land use, proximity to roads, slope, and vegetation cover [13,22]. Cloud-based platforms like Google Earth Engine (GEE) have facilitated the development of FFSMs by providing powerful processing capabilities for large-scale data analysis. Studies employing GEE, such as those conducted by [23], have achieved highly accurate fire susceptibility predictions for India, with Area Under the Curve (AUC) values exceeding 90%.

On the other hand, several studies have focused on classifying various land features and vegetation types for mapping fire-prone areas using only reflectance data from satellite sources, bypassing topographic features like slope, aspect, and distance from roads [4,24,25,26,27]. These models primarily utilize spectral indices from satellite imagery, such as the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Live Fuel Moisture Content (LFMC), and Zagros Gras Index (ZGI) derived from RS data [4,8,24,25,26,27,28,29]. These indices can capture vegetation’s health, moisture content, and stress, all critical factors in fire susceptibility.

Unsupervised techniques, such as Principal Component Analysis (PCA), K-means clustering, IsoData, and Non-negative Matrix Factorization (NMF), have shown significant potential in studying environmental studies [30,31,32]. PCA, K-means clustering, and IsoData are widely used unsupervised techniques that facilitate the analysis of complex datasets. PCA is a dimensionality reduction method that transforms correlated variables into a smaller set of uncorrelated variables, known as principal components, which capture the most variance in the data [33,34]. On the other hand, K-means and IsoData clustering partition the data into K distinct clusters based on feature similarity, making it effective for identifying patterns within large [35]. Collectively, these methods are instrumental in extracting meaningful information from RS data, supporting applications such as land cover classification and environmental monitoring [35,36]. Moreover, NMF is a dimensionality reduction method that breaks down data into non-negative components and their associated weights. This approach provides interpretable results, as the non-negativity constraint leads to distinct, additive parts, making it ideal for applications like image processing and spectral unmixing, where understanding individual, meaningful components is essential [37,38].

Satellite-based models and indices like NDVI, EVI, LFMC, and ZGI leverage spectral data to assess vegetation characteristics, relying on pre-existing knowledge of biomass and land surface properties while often excluding terrain factors such as slope or proximity to roads. These models evaluate vegetation’s health and stress through satellite band analysis [24,25,26,27,28,29].

Unsupervised methods such as PCA, K-means clustering, and IsoData enable automatic detection of patterns in data, which is widely applied in environmental studies and landcover mapping for dimensionality reduction and clustering [39,40]. NMF, though primarily employed for hyperspectral data reduction and unmixing [41], offers unique benefits through interpretable non-negative components [38,41] yet remains less explored in mapping fire susceptibility, highlighting the potential for further research. NMF is particularly suited for post-fire analysis due to its non-negativity constraint, which aligns with the physical properties of spectral reflectance. Its ability to decompose mixed pixels—common in burned landscapes—into interpretable, additive components makes it effective for burn severity mapping [42]. Unlike supervised classifiers, NMF does not require labeled training data, making it advantageous for rapid post-fire assessment.

Besides proposing a modified NMF method, sparse and endmember-independent NMF (SEI-NMF), this study aims to perform a comparative analysis of the latest and three other already developed NMF methods (Basic or Standard NMF, L1 Sparsity, and L1/2 Sparsity) in FSSM regarding other advanced techniques such as PCA, K-means, and IsoData, to evaluate how these unsupervised methods can predict fire-prone in SM and SA area of Kurdo-Zagrosian region, using only satellite-derived data. By extracting features from these images, the study aims to create new satellite-based layers that enhance fire susceptibility understanding, moving beyond the traditional reliance on the NDVI for vegetation representation. Including the ZGI factor accounts for phenological characteristics and enriches the dataset with ecological insights, providing a more comprehensive approach to fire risk assessment.

2. Materials and Methods

2.1. Study Area

The study focuses on the SA and SM climate regions, particularly the forests of Marivan and Sarvabad in Kurdistan Province, western Iran, which have faced numerous wildfires over the past decades (see Figure 1) [8].

Fires typically begin in late May, driven by the accumulation of dry grass species, and persist until the onset of autumn rains. These forests are located between 750 and 1,700 meters above sea level (a.s.l.), neighboring different shrub species [43]. Grass species, highly flammable and widespread, further increase fire vulnerability, especially during the dry summer months, when most of these species dry out, contributing to fire susceptibility [8]. These forests are located near the Iraq border, between longitudes 45°57′50″E and 46°46′41″E and latitudes 35°1′1″N and 35°49′51″N [8]. Marivan sits at an average elevation of 1,287 meters within the northern Zagros mountain range, featuring mountainous terrains and valleys [8]. The region’s SM climate brings cold winters and hot summers, with annual precipitation of approximately 991 mm. The vegetation, predominantly Brant’s oak (Quercus brantii) and Crataegus pontica and Pistacia atlantica, covers dense and semi-dense forests, rangelands, and pastures [44,45]. These areas have experienced numerous fires in recent decades, highlighting the urgent need for studies to address the issue [8].

2.2. Data Sources

The data used in this study are presented in Table 1. Sentinel-2 multispectral imagery, spanning 2020 to 2023, was obtained from the European Space Agency’s (ESA) Copernicus archive via GEE [46]. The elevation data was derived from the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) at 30-meter resolution, later resampled to 10 meters using GEE to meet the resolution of the Sentinel-2 data to avoid scaling biases [47]. All datasets were processed within the Universal Transverse Mercator (UTM) coordinate system, Zone 38N (EPSG:32638). Additionally, the ZGI was calculated from the Sentinel-2 imagery [2].

2.3. C. Sparse and Endmember-Independent Non-Negative Matrix Factorization (SEI-NMF)

Besides applying three scenarios for the already developed NMF methods, this study introduces the sparse and endmember-independent non-negative matrix factorization (SEI-NMF) models. Three variants of NMF were applied: standard NMF, L1-NMF, and L1/2-NMF. The L1 sparsity constraint enforces a more localized, part-based decomposition suitable for isolating fire-affected spectral signatures [48]. At the same time, the L1/2 variant promotes greater sparsity and robustness to noise, which is particularly valuable in hyperspectral fire scene analysis [49]. Essential foundational methods are discussed to contextualize SEI-NMF’s approach. The Linear Mixture Model (LMM) is a widely used model in RS that decomposes pixel spectra into physical components, or "endmembers," representing distinct surface types [37]. NMF, closely aligned with LMM, similarly decomposes multispectral data into non-negative components, making it particularly effective for identifying additive spectral features. SEI-NMF extends NMF by adding sparsity and endmember independence, improving its utility in detecting surface materials and mapping fire-prone regions using remote sensing data [37,38].

This study employs a range of notation to represent mathematical entities. Lowercase letters stand for scalars, while boldface lowercase letters and uppercase letters stand for vectors and matrices, respectively. For any matrix $M$, its $i$-th column and $j$-th row are denoted by $M_i$ and $M_{\left(j\right)}$, and $M_{ij}$ represents its (i, j)-element. The trace of $M$ is represented by $Tr\left(M\right)$, and the transposed matrix of $M$ is denoted by $M^T$. The Frobenius norm of a matrix $M\in {\mathbb{R}}^{d\times n}$ is described as Equation 1 which provides a single scalar value that measures the overall magnitude of the matrix and is commonly used for comparing matrices or for error minimization in optimization problems [38].

$${\left|\right|M\left|\right|}_F=\sqrt{{\sum }_i^n{\sum }_j^d{M}_{ji}^2}=\sqrt{Tr\left(M^{T}M\right)}=\sqrt{Tr\left(M{M}^T\right)}$$

(1)

To clarify the following steps, we will first discuss the concepts of LMM, NMF, and the sparsity regularizer. Subsequently, we will highlight the modifications that lead to the proposed NMF method (SEI-NMF).

2.3.1. Linear Mixture Model (LMM)

Given an observation matrix $X=[x_1,x_2,\dots ,x_n]\in {\mathbb{R}}^{d\times n}$ with $d$ bands and $n$ pixels, where each band $i$ corresponds to a particular wavelength ${\lambda }_i$, the LMM can be expressed as Equation 2.

$$X=WH+E$$

(2)

where $W=[w_1,w_2,\dots ,w_r]\in {\mathbb{R}}^{d\times r}$ and $H=[h_1,h_2,\dots ,h_n]\in {\mathbb{R}}^{r\times n}$ represent the endmember matrix and the abundance matrix, respectively, $r$ denotes the number of endmembers, and $E\in {\mathbb{R}}^{d\times n}$ is the additive noise. In contrast to the number of pixels or the bands, the number of endmembers existing in a multispectral image is often much smaller, i.e., $r<min(d,n)$.

Let $i$ stand for the location of a pixel in the image after putting all the pixels in order. Then, according to LMM, a single-pixel $x_i$ can be approximately represented as a linear combination of endmember vectors (Equation 3).

$$x_i\approx \sum _{p=1}^r{w}_p{h}_{pi}$$

(3)

In Equation 3, the value of $h_{pi}$ measures how much the endmember $w_p$ contributes to the pixel $x_i$.

2.3.2. Non-negative Matrix Factorization (NMF_Basic)

Since NMF can lead to the interpretable part-based representation of high-dimensional data [48], it owns a wide range of applications, such as clustering [49], feature selection [53], recommender system [54], community detection [55], and matrix completion [56]. Let $X$ be the given data matrix, represented as $X=[x_1,x_2,\dots ,x_n]\in {\mathbb{R}}^{d\times n}$, where $d$ is the number of features and $n$ is the number of samples. Each column vector $x_i$ denotes a non-negative data sample with $d$ dimensions. NMF aims to discover two non-negative matrices $W\in {\mathbb{R}}^{d\times r}$ and $H\in {\mathbb{R}}^{r\times n}$, which can accurately reconstruct the data matrix as $X\approx WH$, according to Equation 4.

$$\underset{W,H}{\min }\sum _{i=1}^n\mathsf{\Xi }\left(x_i,W{h}_i\right),s.t.W,H\ge 0$$

(4)

where, according to the measure function Ξ, each sample x_i​ can be reconstructed as a linear combination of the vector bases in $W$, using coefficients given by the vector $h_i$​. The basic NMF model utilizes the square error distance to quantify the difference between $X$ and $WH$. Its objective function is defined as Equation 5.

$$\underset{W,H}{\min }{\left|\right|X-WH\left|\right|}_F^2=\sum _{i=1}^n{\left|\right|{x_i-Wh}_i\left|\right|}^2,s.t.W,H\ge 0$$

(5)

where $\left|\right|{x_i-Wh}_i\left|\right|$ represents the reconstructed error of the i-th sample and ${\left|\right|.\left|\right|}_F$ refers to the Frobenius norm. Although Equation 5 is convex in $W$ and $H$ independently, it loses convexity when both variables are considered simultaneously. Therefore, obtaining the globally optimal solution is impractical, and a locally optimal solution can only be obtained by using optimization methods. It is worth noting that most NMF algorithms are iterative and leverage the fact that NMF can be simplified to a convex non-negative least square problem (NNLS) when either $W$ or $H$ is fixed. Specifically, one of the two factors are held constant throughout each iteration. At the same time, the other is updated to decrease the objective function. The most well-known and widely used optimization can be conducted using the following multiplicative update rules (Equations 6 and 7) [57].

$$W\leftarrow W\odot {\displaystyle \frac{X{H}^T}{WH{H}^T}}$$

(6)

$$H\leftarrow H\odot {\displaystyle \frac{W^{T}X}{W^{T}WH}}$$

(7)

where the operator ${(.)}^T$ represents the transpose of the matrix, and $\odot$ and / denote the elementwise multiplication and division, respectively. Note that the rules denoted by Equations 6 and 7 will not change the non-negativity of $W$ and $H$, provided that the initial $W$ and $H$ are non-negative.

2.3.3. Sparsity regularizer (NMF_L1 and L1/2 Sparsity)

In NMF, the objective function is non-convex, which means it can have many local minimum points. As a result, solutions can vary, making the outcome less stable or unique. To address this, additional constraints are often added to standard NMF. Since multispectral data naturally exhibit sparse abundance patterns—where only a few components significantly contribute to each pixel—a sparsity constraint is applied to the objective function. This approach, introduced in studies like [48] and [58], helps improve stability by focusing the solution on more realistic, sparse data representations. Equation 8 shows the objective function.

$$\underset{W,H}{\min }{\left|\left|X-WH\right|\right|}_F^2+\gamma g\left(H\right)s.t.W,H\ge 0$$

(8)

where $\gamma \ge 0$ is the parameter used to control the contribution of the sparsity measure function $g(·)$ of the matrix $H$, which is regarded as the regularization term. In this study, we introduce two kinds of sparsity regularizers: $L_1$ and $L_{1/2}$. The corresponding L₁-NMF is given as Equations 9 and 10.

$$\underset{W,H}{\min }{\left|\left|X-WH\right|\right|}_F^2+\gamma {\left|\left|H\right|\right|}_{1}s.t.W,H\ge 0$$

(9)

The $L_{1/2}$-NMF is then written as:

$$\underset{W,H}{\min }{\left|\left|X-WH\right|\right|}_F^2+\gamma {\left|\left|H\right|\right|}_{1/2}s.t.W,H\ge 0$$

(10)

where ${\left|\left|H\right|\right|}_{1/2}=\sum _{i,j}^{d,n}{\left(H_{ij}\right)}^{1/2}$ and $H_{ij}$ is the abundance fraction for the $i$th endmember at the $j$th pixel in the multispectral data.

2.3.4. Proposed Method (SEI-NMF)

To develop the proposed SEI-NMF method, the L1-sparsity NMF model was extended by incorporating specific constraints. This enhancement allows the model to capture sparse, and independent features better. The solution space of the NMF model is vast, and the identification of endmembers plays a critical role in fire susceptibility research. Incorporating the characteristics of endmembers as prior knowledge into the NMF model can significantly enhance its performance. This approach enables the identification of more accurate endmembers, leading to improved results. Since different components of the multispectral data contribute with specific proportions, it is essential to ensure their independence. To achieve this, the autocorrelation matrix can be adapted as a constraint, where the independence of these components is reflected in a diagonal autocorrelation matrix. That is, the off-diagonal elements of its autocorrelation matrix should be as close to 0 as possible. Therefore, the NMF model with endmember independence constraint is as Equation 11.

$$\underset{W,H}{\min }{\left|\right|X-WH\left|\right|}_F^2+\alpha \left({\left|\right|W^{T}W\left|\right|}_1-{\left|\left|W\right|\right|}_F^2\right)s.t.W,H\ge 0$$

(11)

where α is the parameter to balance the data fidelity and endmember independence term. The second term refers to the sum of the off-diagonal elements of the autocorrelation matrix for endmembers, i.e., the difference between the sum of all the components (the first sub-term) and the sum of the diagonal elements (the second sub-term). The purpose of the second term in Equation 10 is to make the endmembers independent of each other as much as possible; that is, the correlation between different endmembers should be as small as possible.

Equation 12 describe the ultimate loss function proposed to our SEI-NMF.

$$\underset{W,H}{\min }{\left|\right|X-WH\left|\right|}_F^2+\gamma {\left|\left|H\right|\right|}_q+\alpha \left({\left|\right|W^{T}W\left|\right|}_1-{\left|\left|W\right|\right|}_F^2\right)s.t.W,H\ge 0$$

(12)

where q∈ {1⁄2,1} determines sparse regularization term.

2.3.5. Optimization

Similar to the NMF case, the cost functions (Equation 12) are not convex concerning $W$ and $H$ together. To this end, a standard method is to iteratively optimize the objective function by minimizing one variable while keeping the other variable fixed. We then perform two types of updates. First, we update $H$ while $W$ is fix. Let $\Phi \ge 0$ and $\Psi \ge 0$ be the corresponding Lagrange multipliers. Consider the Lagrange $\mathcal{L}$ for q=1⁄2 as Equation 13.

$$\begin{array}{l}\mathcal{L}={\left|\left|X-WH\right|\right|}_F^2+\gamma {\left|\left|H\right|\right|}_{1/2}+\alpha \left({\left|\right|W^{T}W\left|\right|}_1-{\left|\left|W\right|\right|}_F^2\right)+Tr\left(\Phi W^T\right)+Tr(\Psi H^T)\\ =Tr\left(X^{T}X-2X^{T}WH+H^T{W}^{T}WH\right)+\gamma {\left|\left|H\right|\right|}_{1/2}+\alpha \left({\left|\right|W^{T}W\left|\right|}_1-Tr\left(W^{T}W\right)\right)+Tr\left(\Phi W^T\right)+Tr(\Psi H^T)\end{array}$$

(13)

Then, taking the partial derivative concerning $H$ on both sides leads to Equation 14.

$${\displaystyle \frac{\partial \mathcal{L}}{\partial H}}=-2W^{T}X+2W^{T}WH+{\displaystyle \frac{1}{2}}\gamma H^{-{\displaystyle \frac{1}{2}}}+\Psi$$

(14)

According to the KKT conditions $\Psi \odot H=0$, so it follows to Equation 15.

$$(-2W^{T}X+2W^{T}WH+{\displaystyle \frac{1}{2}}\gamma H^{-{\displaystyle \frac{1}{2}}})\odot H=0$$

(15)

Therefore, the updating rule can be acquired as Equation 16.

$$H\leftarrow H\odot {\displaystyle \frac{W^{T}X}{W^{T}WH+\frac{\gamma }{4}{H}^{-\frac{1}{2}}}}$$

(16)

In the same way, the updating rule for the matrix $H$, based on (9) is represented by Equation 17.

$$H\leftarrow H\odot {\displaystyle \frac{W^{T}X}{W^{T}WH+\gamma }}$$

(17)

Second, we update $W$ while $H$ is fix, taking the partial derivative concerning $W$ on both sides of $\mathcal{L}$ leads to Equation 18.

$${\displaystyle \frac{\partial \mathcal{L}}{\partial W}}=-2X{H}^T+2WH{H}^T+2\alpha \left(W{I}_1-W\right)+\Phi$$

(18)

where $I_1$ represents a matrix in which all elements are equal to 1. Following the similar steps, the update rule for $W$ is as Equation 19.

$$W\leftarrow W\odot {\displaystyle \frac{X{H}^T}{WH{H}^T+\alpha (W{I}_1-W)}}$$

(19)

Figure 2 displays the methodologic framework briefly describing the steps to define the proposed method. As can be seen, the proposed method results from applying new constraints on L1 Sparsity NMF.

2.4. The Number of Components (Endmembers) and Iteration

In the context of dimensionality reduction, mainly when working with large datasets, it is crucial to identify the optimal number of dimensions that can capture the essential features of the data while minimizing complexity. Domain knowledge is often used to estimate the likely number of clusters or components before applying algorithms like K-means [59]. Accordingly, experimenting with different dimensionalities (3, 4, and 5), it was found that using three dimensions yielded the optimal results [59]. This suggests that three dimensions are sufficient to represent the data effectively, allowing for a more manageable and interpretable biomass analysis without significant information loss and less time-consuming.

Additionally, iterations play a crucial role in converging iteration-based algorithms, such as NMF. Therefore, various iterations should be considered to determine the optimum iteration for each tested method. This finding will indicate that the models can reach a stable solution within this number of iterations, balancing computational efficiency and model performance.

2.5. Labelling and Validation

PCA, K-means, and NMF methods create components based on the inherent structure in the data. These components are derived from patterns and relationships in the data and do not rely on arbitrary thresholds for classification or pre-defined labeled data. Instead, the methods aim to find the underlying variation or similarity in the data, grouping values or creating components that explain this structure.

Wildfire prediction is inherently challenging due to the complex interactions between natural factors, such as vegetation and climate, and human activities, which often act as triggers. Since human actions rather than natural conditions cause most wildfires, validating predictions becomes more complicated. However, by combining spatial consistency and consistency across models [60], we can build a solid foundation for assessing and labeling fire-prone areas despite these challenges.

To assign appropriate labels to the output, geographic information categorizes the results into classes such as High, Average, or Low fire susceptibility. In the unsupervised methods, the labeling process occurs after the methods are applied based on the interpretation of the results concerning geographical and field data. In contrast, supervised models rely on pre-labeled data, where the labeling is completed before applying the method. This distinction highlights the exploratory nature of unsupervised approaches, which derive insights from the data, requiring subsequent interpretation to assign meaningful categories. In this case, the fire-prone regions identified by the model were cross-referenced with real-world fire-prone characteristics, such as dense vegetation zones or urban interfaces where human activities are likely to ignite fires. For instance, fire-prone areas are often closely related to specific geographic features such as slope, elevation, vegetation types, and proximity to water bodies [61,62]. When these geographic factors are used to label classification results, they can help distinguish between areas with higher and lower fire risk, creating meaningful classes (e.g., high, average, and low fire risk). This approach provides interpretive validation, ensuring that high-risk areas align with regions typically associated with wildfire susceptibility.

To further evaluate the effectiveness of our method, the proposed SEI-NMF technique was compared with multiple classification approaches, including K-means, PCA, IsoData, and three variations of previously developed NMF methods.

This comparison was aimed at assessing the consistency of fire-prone area classification. Regions consistently identified as high-risk across these models are likely reliable indicators of fire-prone areas. By identifying stable patterns across different unsupervised techniques, we validate the robustness of our findings. Moreover, adjusting the number of components or classes may change the classification outcomes, but the results remain interpretable.

The validation process was conducted through spatial overlap analysis, comparing the areas classified as high, average, and low fire susceptibility with burned areas recorded after 2020. To verify and assess the effectiveness of the applied models in identifying fire-prone areas, burned area data derived from Sentinel-2 satellite imagery from 2021 to 2023 were employed. These burned areas were calculated based on the Normalized Burn Ratio (NBR) alongside statistical records from the Kurdistan Province Natural Resources Administration in Iran. Equation 20 represents the NBR index calculated using Sentinel-2 imagery, utilizing Band 8 (Near-Infrared (NIR)) and Band 12 (Short-Wave Infrared (SWIR)), as these bands are susceptible to burn-related changes in vegetation [63,64].

$$\mathrm{N}\mathrm{B}\mathrm{R}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}{NIR+SWIR}}$$

(20)

Figure 3 presents the methodological framework used.

3. Results

The results are presented in two sets: first, the maps generated using the applied methods, and second, the statistical results, provided numerically and in chart format for clear visualization and comparison. The results compare model performance in two scenarios, with and without elevation data.

Figure 4 represents the convergence chart of each NMF scenario used. Convergence in NMF methods refers to the process by which the algorithm iteratively improves its solution until it reaches a stable point where further iterations do not significantly change the results. This property is essential for ensuring the reliability and robustness of NMF, especially in unsupervised tasks like feature extraction and dimensionality reduction. However, convergence can be sensitive to initialization and the chosen number of components, highlighting the need for careful model setup. In practice, methods like NMF-L1 Sparsity, which enforce sparsity constraints, tend to improve convergence by promoting more interpretable and meaningful features for real-world applications, such as identifying fire-prone zones, compared to the Basic NMF method.

In the first scenario, Sentinel-2 data, ZGI, and elevation were used to assess fire susceptibility (Figure 5). In contrast, only Sentinel-2 data and ZGI were applied in the second scenario, making it solely based on multispectral information (Figure 6).

By testing these two scenarios, the study examines whether including elevation improves model performance or disrupts clustering when combined with Sentinel-2 data. DEM, a static variable representing topographical elevation, is recognized in many fire susceptibility models, especially supervised methods. However, DEM differs fundamentally from Sentinel-2 reflectance data, which captures dynamic surface characteristics like vegetation health, moisture, and land use. These differences can influence clustering in unsupervised methods, potentially altering cluster composition. Comparing both scenarios highlights how terrain and surface data interact, affecting the adaptability of unsupervised clustering in fire susceptibility assessments.

Figure 5 displays the results of applying discussed unsupervised methods, including Basic NMF (Figure 5a), L1-Sparsity (Figure 5b), L1/2-Sparsity (Figure 5c), SEI-NMF (Figure 5d), PCA (Figure 5e), K-means (Figure 5f), and IsoData (Figure 5g) on Sentinel-2 satellite data, ZGI, and elevation data.

The same way, Figure 6 presents the resulted maps of applying Basic NMF (Figure 6a), L1-Sparsity (Figure 6b), L1/2-Sparsity (Figure 6c), SEI-NMF (Figure 6d), PCA (Figure 6e), K-means (Figure 6f), and IsoData (Figure 6g) on Sentinel-2 satellite data, and ZGI. The dark areas represent the high-fire susceptible areas, while the orange and blue areas represent the average (moderate) and low-fire susceptible areas, respectively. Furthermore, the red polygons represent the post-2020 burned areas, and the yellow polygons represent the water bodies (Zrebar Lake and Garan Lake). Marivan town, as an urban area, is displayed as a black polygon within the maps.

Figure 7a shows the statistical distribution of each class for the 2020 maps, focusing specifically on areas that experienced fires from 2021 to 2023 (post-2020 fires) regarding the first scenario within elevation included. Figure 7b represents the second scenario without the elevation factor.

4. Discussion

In applying NMF, PCA, K-means, and IsoData clustering to Sentinel-2 data, we observed variability across different runs, mainly due to these algorithms’ random initialization processes and non-convex nature. Without a fixed random seed, these methods may converge to different local minima, producing slightly different decom-positions or cluster assignments [48,65]. This is particularly evident in NMF, where initial values significantly influence results. While setting a fixed random seed could improve consistency, some variability remains inherent due to the complexity and dimensionality of RS data [48,66]. Such considerations are crucial for ensuring reproducibility and robustness in analyses using high-dimensional satellite imagery.

Regarding the first scenario encompassing Sentinel 2, ZGI, and elevation, the results demonstrate that the NMF methods offer superior predictive performance compared to PCA, K-means, and IsoData. Specifically, the L1/2-Sparsity NMF model provided the best results [48] and identified 72.2% of the burned areas as high-risk, with an overall combined classification of 95.1% when high and average risk areas were considered. The L1-Sparsity, Basic NMF, and SEI-NMF models also performed robustly, capturing 60.8%, 55.3%, and 58.8 % as high-risk and 78.2%, 77.7%, and 78.4% as high plus average risk zones, respectively. This consistency across NMF variants suggests a heightened sensitivity to critical fire-prone features, which may be attributed to NMF’s decomposition of data into parts-based representations, helping isolate areas with characteristics strongly associated with fire susceptibility [38,48].

In contrast, while yielding 79.6% in the high and average categories, PCA classified only 41.4% of burned areas as high-risk, underscoring its lower precision in identifying zones with the highest fire susceptibility. K-means and IsoData further demonstrated limitations, with high-risk classifications of 31.3% and 33.9%, respectively (Figure 6a), and struggled to identify non-flammable features like water bodies accurately (Figure 5f-5g, and 6f-6g). It can be due to the high sensitivity of K-means and IsoData to outliers, high dimensionality, and noise [35].

On the other hand, when only Sentinel-2 bands and the ZGI index were included, the results show that the SEI-NMF method outperformed all other models in predicting fire-prone areas, capturing approximately 72% of burned areas as high-risk and nearly 90% as high and average risk combined. The L1/2-Sparsity NMF model followed closely, identifying 70.4% as high-risk and 99% as high-plus-average areas. Basic NMF and L1-Sparsity models achieved around 57% and 64% in high-risk classifications, respectively (Figure 7b). When DEM data is included, its inherent characteristics may introduce noise, which can disrupt accurate spectral discrimination compared to using Sentinel data alone, which has distinct and consistent spectral properties. DEM data may dilute spectral signals by emphasizing elevation-driven features that are less relevant to the land cover targets (endmembers), reducing the algorithm’s capacity to focus on the critical spectral characteristics of endmembers. By relying solely on Sentinel-2 data, the modified NMF can better isolate and analyze the pure spectral signals, free from the noise introduced by mixed-data characteristics of DEM. This characteristic can also be generalized to other data like slope, aspect, and distance from roads, which have been widely used to detect fire-prone areas, especially supervised methods.

PCA and K-means results remained largely unchanged from previous analyses (where elevation was included), with K-means and PCA showing limited ability to differentiate high-risk areas. Additionally, the NMF models showed a solid ability to recognize non-fire-prone regions such as water bodies, mountainous areas with minimal vegetation cover, and urban zones, which typically exhibit the lowest susceptibility to wildfires (yellow polygons in Figure 4a- 4d and 5a-5c). These regions were consistently classified as low-risk areas, highlighting NMF’s effectiveness in distinguishing between fire-prone and non-fire-prone landscapes. PCA similarly identified these low-susceptibility areas accurately, although with slightly less precision than the NMF models (Figure 4e and 5e). This ability to detect and categorize non-flammable zones with high accuracy further validates the NMF methods’ applicability in comprehensive fire-risk mapping, ensuring non-vulnerable areas are correctly excluded from high-risk classifications. K-means and IsoData failed to detect water bodies or urban areas as low fire-prone areas. If the K-means and IsoData methods are adjusted such that water bodies are classified as low fire-prone areas, it could lead to inconsistencies with ground observations and the performance of other methods. Typically, water bodies, due to their lack of vegetation and low flammability, are expected to be classified as low-risk areas for wildfires. However, when not correctly distinguishing water bodies, these methods may misclassify them. This would conflict with observations and other models’ classifications, which align water bodies with areas of minimal fire susceptibility. Adjusting the classification could disrupt the model’s consistency, as real-world fire-prone zones (such as areas with dense vegetation) would be misrepresented, undermining the overall validity of the classification. This could lead to a less accurate representation of fire-prone areas and affect resource allocation and fire management strategies based on these classifications.

The urban area, such as Marivan town (shown as a black polygon inside the maps), is also recognized as low fire-prone. The performance of the methods used within the two scenarios differed, with and without elevation data. When elevation was considered, almost all methods failed to detect it as a low fire-prone area, while excluding elevation, practically all NMF methods and PCA recognized it accurately, and even K-means and IsoData partly recognized it. However, the NMF-L1/2 method recognizes the town as partly average and low fire-prone. The results also showed that 150 iterations were generally enough to achieve stable NMF outcomes (Figure 4a-4d). However, testing different iteration sizes or using predefined, reliable initial values could further improve model performance.

NMF has shown great promise in various fields, particularly in dimensionality reduction and uncovering hidden patterns within complex data. Despite its proven capabilities, NMF has not yet been applied in the context of wildfire detection and burn severity mapping [41,48]. This gap presents a significant opportunity, as NMF’s ability to handle heterogeneous data could enhance wildfire research, especially when dealing with the varied and complex features found in wildfire datasets. Some methods, like Segmented Non-negative Matrix Factorization (SNMF), have outperformed PCA and NMF in hyperspectral image analysis [64], indicating a significant advantage in data processing. Yet, the potential of NMF in wildfire applications remains untapped despite its successes in other fields [64,67].

In contrast, other techniques like PCA have been utilized in this domain. Zheng et al. [68] demonstrated the effectiveness of PCA with Polarimetric Synthetic Aperture Radar (PolSAR) data and random forest algorithm, achieving promising correlations by integrating vegetation structure. Similarly, PCA has been used in multiple studies, such as those by Henry and Epting et al., to differentiate burned areas using Landsat data [69,70].

Additionally, hybrid methods like K-means clustering combined with neural networks and the Forest Fire Hybrid Detection Model (FFHDM) have demonstrated remarkable success in fire detection. These techniques utilize clustering and advanced classifiers to improve detection efficiency [71,72].

By integrating NMF into wildfire research, especially with other advanced techniques like PCA and machine learning models, researchers could unlock new potential for more accurate and efficient wildfire detection and burn severity mapping. While advanced deep learning methods could potentially improve performance, they often require large labeled datasets and significant computational resources, which may not be feasible in time-sensitive post-fire scenarios. Our choice of interpretable, lightweight methods such as NMF and K-means enables more practical deployment in operational wildfire monitoring. Future work may investigate hybrid approaches combining deep feature extraction with unsupervised clustering.

5. Conclusions

This study underscores the critical need for effective fire susceptibility mapping in SM and SA regions, which face increasing wildfire threats. By proposing a modified NMF method and comparing it with other unsupervised techniques such as PCA, K-means, and IsoData, the research provides valuable insights into the potential of satellite-derived data for wildfire risk prediction. The results demonstrate that NMF, especially the L1/2 sparsity method, performs best in detecting high-risk fire areas when using Sentinel-2 bands, ZGI, and elevation data, achieving up to 72% high-risk areas. When elevation data was not considered, the proposed NMF method performed similarly well, identifying 72% of burned areas as high-risk and nearly 90% as high and average combined.

However, challenges remain in improving the accuracy of these unsupervised models. Variability in results due to random initialization and the non-convex nature of these algorithms suggests the need for more controlled experimentation. To address these challenges, future research could explore integrating NMF results into supervised machine learning models, where they could serve as features to enhance classification accuracy. Additionally, using predefined initial values and testing different iteration sizes may improve model stability and reduce outcome variability.

Further improvements could include refining the classification of non-fire-prone areas, such as water bodies and urban zones, which K-means and IsoData struggled to classify correctly. A more robust approach would involve integrating these unsupervised methods with supervised classification techniques, where ground truth data or expert labels could provide the final validation layer. By doing so, building more reliable and contextually accurate fire risk maps would be possible. Moreover, combining these satellite-based models with environmental factors like weather patterns, human activities, and land cover changes could offer a more holistic view of fire susceptibility.

In conclusion, this study provides a valuable foundation for future wildfire management strategies in SM and SA regions. This framework is potentially applicable on other areas with similar weather regimes, such as Portugal, Spain, and Greece. Future work should focus on improving the consistency of unsupervised methods and exploring their integration into supervised machine learning methods, which could result in more accurate, reliable, and actionable fire-prone area predictions.