Improved Lithological Map of Large Complex Semi-Arid Regions Using Spectral and Textural Datasets within Google Earth Engine and Fused Machine Learning Multi-Classifiers

Abstract

In this era of free and open-access satellite and spatial data, modern innovations in cloud computing and machine-learning algorithms (MLAs) are transforming how Earth-observation (EO) datasets are utilized for geological mapping. This study aims to exploit the potentialities of the Google Earth Engine (GEE) cloud platform using powerful MLAs. The proposed method is implemented in three steps: (1) Based on GEE and Sentinel 2A imagery (spectral and textural features), that cover 1283 km² area, a variety of lithological maps are generated using five supervised classifiers (random forest (RF), support vector machine (SVM), classification and regression tree (CART), minimum distance (MD), naïve Bayes (NB)); (2) the accuracy assessments for each class are performed, by estimating overall accuracy (OA) and kappa coefficient (K) for each classifier; (3) finally, the fusion of classification maps is performed using Dempster–Shafer Theory (DST) for mapping lithological units of the northern part of the complex Paleozoic massif of Rehamna, a large semi-arid region located in the SW of the western Moroccan Meseta. The results were quantitatively compared with existing geological maps, enhanced color composite and validated by field survey investigation. In comparison of individual classifiers, the SVM yields better accuracy of nearly 88%, which was 12% higher than the RF MLA; otherwise, the parametric MLAs produce the weakest lithological maps among other classifiers, with a lower OA of approximately 67%, 54% and 52% for CART, MD and NB, respectively. Noticeably, the highest OA value of 96% is achieved for the proposed approach. Therefore, we conclude that this method allows geoscientists to update previous geological maps and rapidly produce more precise lithological maps, especially for hard-to-reach regions.

1. Introduction

One of the most challenging geological applications using remote-sensing-based satellite data is mapping lithological features, especially for large and geologically complex areas that require high-precision lithological maps [1,2]. Furthermore, the limited availability of the spatial and spectral quality of open access image data (e.g., Landsat-5 TM, Landsat-7 ETM+ (enhanced thematic mapper plus), Landsat-8 OLI (operational land imager), ASTER and Sentinel-2) is widely utilized to extract lithological [3,4,5], mineral [6,7] and structural information [8,9,10]. The majority of studies investigating the potential of remote sensing for geological mapping [4,6,11,12,13] have been conducted over relatively small geographical areas using individual machine learning algorithms (MLAs). Therefore, an innovative approach, based on several factors, is implemented in this research to overcome the aforementioned issues. First, we choose well-discriminating spectral and textural information; second, we perform an efficient approach based on fused multi-classifiers to optimize accuracy; third, we exploit the high computational processing capacities of Google Earth Engine (GEE) to process free-access remote sensing (RS) geospatial big data and use MLAs over a large area.

Sentinel-2 multispectral data was used due to its high potential for geological mapping [14]. A handful of studies have investigated the potential of combining spatio-spectral bands for lithological classification over wide geologically complex regions [15,16,17]. Furthermore, different rock types with similar spectral properties can lead to spectral overlap and misclassification [18,19,20]. Hence, many classification approaches have been proposed in the literature to improve classification accuracy by introducing textural characteristics, which can be considered the most discriminating property between different lithological units when the spatial resolution of the sensor is enhanced [16,21,22].

There are several texture analysis approaches, such as the grayscale co-occurrence matrix (GLCM) [23,24], fractal analysis [25], discrete wavelet transform [26], Laplacian filter [27], Markov random field [28], granulometric analysis [29,30]. The most commonly used in geological applications is the grayscale co-occurrence matrix (GLCM) [31,32,33]. It is a second-order statistical technique that measures the relative occurrence of reference pixels compared to neighboring pixels and incorporates spatial information in the form of index images that quantify this texture. Second-order statistics have an advantage over first-order statistics since they include spatial variability information. The three principal computational parameters of GLCM are displacement, neighborhood size and quantization. According to previous researches, it is preferable to analyze the textural characteristics of high-resolution imagery with a small neighborhood size [16].

Haralick et al. [32] proposed a total of 14 second-order statistical measures that may be obtained using the GLCM probability matrix. Among them, the following four statistical metrics are frequently used: contrast, ASM, entropy and correlation [28,34]. Entropy is the most complex statistical parameter because it includes an effect related to image heterogeneity [35]. The quantification of the texture of the image facies is performed by calculating texture indices in the form of additional layers that are combined with the spectral layers in the different bands of the multispectral image. Generally, the integration of textural information within the classification algorithms enhances the lithological mapping when fewer (lower) correlated indices are used [16,22,33]. The contrast index and entropy are the lowest correlated and most discriminating indices for the texture of geological bedrock [36]. Applying these indices enhances the accuracy of classification algorithms in lithological mapping, although the amount of improvement depends on the type of rock [11]. However, the integration of spatial information in image classification approaches exposes these algorithms to a new challenge related to the processing of this multi-dimensional data (spectral, texture, geometry and shape) for large-scale geological mapping. Other issues that make lithological mapping more challenging are the limited robustness of the classification methods offered by both free and commercial software.

The application of artificial intelligence algorithms has been widely used in the field of geological applications [4,5,13]. Although many studies have attempted to compare different algorithms to determine the best classification approach for geological mapping, no consensus on the validity of the available MLAs has yet been achieved. Each classification algorithm can have very good performance for the discrimination of a lithology but may have some disadvantages when compared to other lithologies [37]. In other words, each MLA has its own pros and cons that depend on the geological complexity of the study area, the data quality of available RS imagery and the spatiospectral characteristics of each lithological unit. Therefore, the fusion approach based on multiple classifiers that merge into a single map may yield an improved result [38]. Due to its potential to increase accuracy, this technique has been gaining prospects over the last decade and is being widely used for several mapping applications, including land cover mapping [37,38,39,40], wildfire susceptibility mapping [41,42], water edge detection [43] and change detection [44], with limited application in lithological mapping. Notably, the classifier fusion method performed better when the accuracy of each classifier was more than 50%. According to theoretical and practical investigations, an optimum fusion outcome is obtained when individual classifiers are both accurate and contradictory to some classes. The Dempster–Shafer fusion (DST) [45] is one of the most frequently used fusion methods in the geospatial field; the information acquired from several sources is represented by the degree of belief/mass function, then merged/aggregated using the DST combination rule [46]. Therefore, the D-S theory is a multi-source data fusion approach for collecting more reliable data [47,48]. It is an accurate fusion technique that uses belief uncertainty intervals to provide the beliefs of hypotheses based on evidence from multiple observations [46]. The algorithm uses the reasoning, weight and probability of evidence contained in the dataset [46,47,48]. This method (DST) has been evaluated in many applications, including automatic pattern recognition [49], forensics [50] and fingerprint verification [51]. However, its application in the field of lithological mapping is non-existent [52].

The GEE (https://earthengine.google.com; Casablanca, Morocco; accessed on 5 September 2022) [53] is a free cloud computing platform that enables users to access and manage petabytes of satellite data on a large scale using JavaScript and Google’s cloud [54]. GEE uses Google’s computing capabilities to decrease uptime and provide a repository for script storage and exchange, allowing extensive user collaboration with minimal cost and equipment [55]. It provides a variety of packages for image collection, analysis, processing, classification and export to ensure that users are no longer solely dependent on expensive commercial software [55,56,57]. Several researches use GEE to take advantage of its massive data catalogue to examine dynamic processes over long time-series data and to generate large-scale thematic classifications for a variety of applications, including LU/LC mapping [58], cropland classifications [59], forest habitats mapping [60], surface water detection [61], urban and rural settlement [62], mine mapping [63], natural hazard mapping and snow and shoreline detection [64,65]. However, its application in geological mapping remains very limited [54]. This study provides an improved approach that is developed first in the GEE environment, collecting data and testing various variable combinations to enhance each classification step. Secondly, the DST fusion method was used to optimize lithological mapping and subsequently improve accuracy by considering the advantages of five MLAs classified outputs.

The novelty of this research lies in the proposed methodology for analyzing the additional value of integrating spatio-spectral bands and evaluating the performance of several MLAs available in the GEE platform for lithological classification over large geological complexes regions. Furthermore, in order to optimize and improve the classification accuracy of the exposed rock units, the fusion of multiple classifiers based on the DST method has been employed.

The study described in this paper has the following three main objectives:

• Evaluating the potential of the GEE platform and spatio-spectral bands of Sentinel-2 data to classify the lithological units exposed in a large complex semi-arid region within the GEE code editor environment;
• Evaluating and assessing the performance of different MLAs (RF, SVM, CART, MD and NB) in terms of classification accuracy of each class;
• Optimizing and enhancing lithological mapping accuracy for all the classes using the DST fusion approach.

The main goal of this research was to develop a scheme that could be adapted and deployed to produce an accurate lithological map of large geologically complex regions in a more operational and automatic approach.

2. Study Area and Materials

2.1. Study Area

The northern part of the complex Paleozoic massif of Rehamna is chosen as a geological test site to verify the proposed approach in lithological mapping. It is characterized by good bedrock exposure because it is a dry and semi-arid region, which is situated in the SW of the western Moroccan Meseta between latitudes 32°31′N and 32°47′N and longitudes 7°31′W and 8°03′07.1″W, with a total area of approximately 1283 km² (Figure 1). Geologically, the northern part of the Paleozoic massif of Rehamna is divided into the following two main subareas: the Paleozoic basement, which was produced during the Hercynian orogeny and eroded by preserving the paleoreliefs, and the tabular Meso-Cenozoic cover, which is slightly deformed by small normal faults [66].

The Hercynian basement formations were formed from the Cambrian to the Permian ages [66]. The Cambrian formations of the Imfout syncline are composed of siliciclastic and pyroclastic [67]. The Ordovician is also represented by coarse- to fine-grained siliciclastic terrains, particularly pelites and sandstones, with an alternation of large quartzitic bars displaying the important reliefs of Oued Kibane [67]. Clay and pelitic rocks from the Silurian are poorly exposed in the study area. The lower and middle Devonian are composed of pelites with intercalation of thin, gray bioclastic and massive limestones [68]. Otherwise, the formations of the upper Devonian are also pelitic but intercalated with quartzitic bars that are occasionally capped by a carbonate–sandstone formation that has been dated as Strunian, which is exposed in Mechra Ben Abbou and Foum el Mejez [69,70]. The Carboniferous (Upper Visean) comprises carbonate rocks and the alternation of limestone layers with pelites, in addition to pelite formation with intercalation of basic magmatic rocks, which are exposed as sills and dikes in Gada Jenabia [71]. The erosion products of the Hercynian edifice chain represent a large outcrop of the continental red deposits in the Paleozoic massif of Rehamna that comprise conglomeratic, sandstone and pelitic materials [66].

The epicontinental red cover base sediments are composed of conglomerates, sandstones and pelites. However, the underlying series comprises limestone and marl marine formations. The quaternary terraces with alluvium (conglomerate and sand) have been deposited along the Oum Er-Rbia River [66].

2.2. Data Sources and Preprocessing

In general, the preparation of the base dataset is a crucial step for any lithological classification scheme [7]. Spectral and textural features of Sentinel-2A (S2A) level-2A (L2A) optical data were used in this study. The S2A is a wide-swath, high-resolution and free-access multi-spectral instrument (MSI). It collects remote sensing datasets in 13 spectral bands situated in the visible, near infrared (VNIR), red edge (Re) and shortwave infrared (SWIR) spectral domains, with spatial resolution ranging from 10 to 60 m [72,73] (

Table 1. Spectral and spatial characteristics of Sentinel 2A MSI.

Band Number	Spectral Characteristic	Central Wavelength (nm)	Spatial Resolution (m)
B1	Coastal aerosol	443	60
B2	Blue (B)	490	10
B3	Green (G)	560	10
B4	Red (R)	665	10
B5	Vegetation red edge 1 (Re1)	705	20
B6	Vegetation red edge 2 (Re2)	740	20
B7	Vegetation red edge 3 (Re3)	783	20
B8	Near infrared (NIR)	842	10
B8a	Near infrared narrow (NIRn)	865	20
B9	Water vapor	945	60
B10	Shortwave infrared Cirrus	1380	60
B11	Shortwave infrared 1 (SWIR1)	1910	20
B12	Shortwave infrared 2 (SWIR2)	2190	20

).

It is worth noting that S2A datasets used in this study were filtered by the region of interest (ROI), the period of acquisition (23 September 2020), and the highest percentage of cloud cover (10). The S2A cloud cover masking is generated using the QA60 band (60 m), which comprises dense clouds and cirrus [74]. This band is given in GEE along with the S2 surface reflectance data, as GEE reduces most of the time-consuming preprocessing steps needed in traditional computing software [54].

In the next step, satellite image spectral enhancement techniques were performed by running various methods, such as minimum noise fraction (MNF) and principal component analysis (PCA).

3. Methodology

The investigation started with the integration of image spectral enhancement and spatial texture features for better discrimination of lithological units. Subsequently, a multi-classification scheme is fused for improved lithological mapping of a large, complex semi-arid area. The general processing chain of the proposed approach is displayed in Figure 2.

3.1. Data Enhancement and Processing

Various algorithms and techniques, such as principal component analysis (PCA), minimum noise fraction analysis (MNFA) and textural feature indices, were applied to Sentinel 2-A imagery in order to improve the discrimination of the lithological units exposed in the study area. The details of the processes are described below.

3.1.1. Spectral Enhancement Techniques

• Principal Component Analysis (PCA)

PCA [75] is a multivariate statistical and feature reduction process [5,76,77,78], to emphasize and highlight spectral information related to geological features by eliminating the irradiance effects of all MSI bands. It can be applied to multispectral datasets by transforming a number of correlated and high dimensional spectral bands into a set of linearly uncorrelated lower dimension output principal components (PCs) [75,79,80]. This transformation technique is performed by computing the covariance and correlation matrix and defining new orthogonal axes data, by removing redundancy in the dataset, filtering noise in the latest PCs and thereby generating enhanced contrast images represented by a false color composite (FCC) [81]. Two image components have been selected in this study (PC2 and PC5) from the application of standard PCA transformation to the VNIR and SWIR bands of S2A datasets.

2.

Minimum Noise Fraction Analysis (MNFA)

The second technique that we used is the minimum noise fraction analysis (MNFA), which is used to derive noise-free PCs, which show higher spectral contrast and showcase earth surface materials [82]. The MNF transformation involves two successive PCA rotations; the first one is also called noise-whitening and is based on the application of noise covariance matrix to decorrelate and rescale the noise in the datasets [83,84,85,86,87,88]; the second step is a standard PCA transform derived from the noise whitened data [89]. Depending on the MNFA band data, the first three components contained the largest amount of geological information and allows delineation of corresponding lithological boundaries of the study area (Figure 3).

3.1.2. Textural Feature Processing

The spatial features of a geological object are its shape, size and texture [21,90,91]. The latter is the simplest feature to extract from spatial datasets since it does not require a segmentation step [90]. When the texture is well quantified and exploited, it represents the spatial pattern of the gray levels of pixels and improves the discrimination of lithological units [21,90,91,92]. In geological mapping, the following several textural analysis approaches have been used: gray level co-occurrence matrix (GLCM) [23,24], fractal analysis [25], discrete wavelet transformation [26], Laplace filters [27], Markov random fields [28], granulometric analysis [29,30] or mathematical morphology analysis [93,94]. In addition to its availability on the GEE platform, the GLCM matrices are the most popular parameters for the quantification of textural information. It is a statistical computation approach of second-order histograms that uses a lot of processing techniques and produces a massive data set with many features [94]. Hence, the use of the GEE platform is preferred, which allows the extraction of GLCM matrices rapidly and easily for a large study area. Average, homogeneity, variance, entropy, contrast, dissimilarity and correlation are the texture measurements extracted from the GLCM in GEE. These parameters are frequently redundant when derived from multispectral imagery [93]. Therefore, PCA was used to select the richest image in textural information with the highest variance. In this study, only the first principal component (PC) was used to calculate the texture data since it accounts for 98.67% of total variations.

The GLCM defines a square matrix whose size N_g × N_g, with N_g is equal to the largest gray level appearing in the image. It shows that each image pixel’s spatial distribution contains textural information and is represented in Equation (1) as follows:

$$P_{n,m}\left(d_{\theta },N_g\right)=\left\{\left(\left(x,y\right),\left(x+a,y+b\right)/ I\left(x,y\right)=n;I\left(x+a,y+b\right)=m\right)\right\}$$

(1)

where the parameter d is the distance between each pair of pixels in the image, θ is the angle (generally 0°, 45°, 90° or 135°) and N_g is the number of gray levels considered in the image, n, m with n, m = 1, 2……, N_g are the gray levels of the two pixels separated by the displacement vector d of coordinates (a, b) and (x, y) are the spatial coordinates of the pixel having gray level n and (x + a, y + b) of the pixel with gray level m.

According to previous studies, it is preferable to analyze the textural properties of open access images with medium to coarse resolutions with a small moving window size [31]. Since the GLCM results depend on the values of N_g and d_θ, in this study the rotation invariant was ensured by using N_g = 256, d = 1 and the average value of the GLCM obtained for 0°, 45°, 90° and 135°, with 3 × 3 moving window covering the neighboring pixels. Two texture indices, contrast (CON) and entropy (ENT), were used due to their high decorrelation as defined in Equations (2) and (3).

$$\mathrm{CON}={\displaystyle \sum _{i,j=1}^{N_g}}{\left(i-j\right)}^2{P}_{ij}$$

(2)

$$\mathrm{ENT}=-{\displaystyle \sum _{i,j=1}^{N_g}}{P}_{ij}log{P}_{ij}$$

(3)

3.2. Machine Learning-Based Techniques

3.2.1. Random Forest (RF)

Random forest (RF), developed by Leo Breiman in 2001 [95], is a supervised non-parametric machine learning classifier that is based on bagging technique [96]. RF produces multiple trees by randomly splitting a predefined number of variables to divide at each node of the decision trees. It predicts classes using a majority vote based on the partitioning of data from numerous decision trees by implementing the Gini Index to choose the best split threshold of input variables [97,98]. This classifier requires the selection of the input variables for each decision tree (ntree) and the number of possible characteristic parameters (mtry) that can be randomly sampled for splitting at each node of the trees in the forest [95]. Therefore, among current remote sensing classification algorithms, RF has become the most extensively used in lithological mapping [5,99,100].

3.2.2. Support Vector Machine (SVM)

The support vector machine (SVM) is a popular non-parametric MLA based on the statistical learning theory proposed by Vapnik in 1995 [101]. The main aim of the SVM algorithm is to determine the optimal hyperplane (decision boundary) separating various input training samples (support vectors) from each class [102]. The parameters that need to be adjusted and optimized when using the SVM for classifying remote sensing images are the kernel function type (polynomial, sigmoid, linear and radial basis function), gamma kernel function (GKF) and the penalty parameter [103]. Studies have shown that the SVM has advantages in terms of small samples and high dimensionality [104]. However, the fast performance and accurate results make SVM one of the widely used lithological mapping classifiers [4,105].

3.2.3. Classification and Regression Tree (CART)

Classification and regression tree (CART) is a binary non-parametric decision tree classifier. It is based on a hierarchical decision tree (DT) framework, by recursively splitting parent nodes into child nodes until they reach a pre-defined threshold of non-splitable terminal nodes (also known as a “leaf node”) [106]. CART implements the information gain ratio (Gini impurity index) [107] to determine which input features will provide the optimum splitting at each node [106]. The higher Gini index value, the smaller the purity of the dataset. The great disadvantage of this technique is its extreme sensitivity to the training dataset. Even little changes in the training data can result in entirely different tree structures. However, this algorithm can create complex trees to solve complicated issues with large datasets. This classifier has been mildly used in lithological mapping [108].

3.2.4. Minimum Distance (MD)

The minimum distance (MD) classification is a classic supervised machine learning method, widely used in lithological mapping and mineral exploration [109,110,111]. It is based on computing the Euclidean distance between each image pixel vector and the mean vector for each class, in a multi-dimensional space using the mean vectors of each class [112].

3.2.5. Naïve Bayes (NB)

Naïve Bayes (NB) [113,114] is a statistical learning algorithm that uses the Bayes theorem to predict a class by estimating the highest posterior probability for each class under the naïve assumption of the inputs for each class that are conditionally independent [115]. NB has been found to be competitive against many different classification methods in a variety of practical tasks [116,117]. NB minimizes the issues of class discrimination to find class conditional marginal densities, which indicates the probability that a given sample belongs to one of the target classes [118].

3.3. Lithological Mapping Based on Dempster–Shafer Fusion

Dempster–Shafer theory (DST), also known as classical Dempster–Shafer evidence theory, is a well-defined fusion technique that uses the belief confusion matrix to combine different classified features into a new class [46]. It is a formal framework for defining and reasoning with uncertain information and which makes it possible to represent both imprecision and uncertainty using the mass function or evidence mass (m), the belief function or credibility (Cr) and the disbelief function or plausibility (pl) [119]. The masses of evidence are defined on a set called the framework of discernment Ω. This set contains the N hypotheses $H{y}_i$ necessary for the complete description of the situation, as demonstrated in Equation (4) as follows [120]:

$$\Omega =\left\{H{y}_1,H{y}_2,\dots \dots , H{y}_N\right\}$$

(4)

In this study, the hypotheses that describe the situation are the classes S_i to which the pixels of the image are assigned as demonstrated in Equation (5) as follows:

$$\Omega =\left\{{\mathrm{S}}_1,{\mathrm{S}}_2,\dots {\mathrm{S}}_j,\dots ,{ \mathrm{S}}_N\right\},$$

(5)

where, N the total number of classes.

The reasoning concerns the set E of parts S of Ω given in Equation (6) as follows:

$$\mathrm{E}=\left\{\mathrm{S}/\mathrm{S}\subset \mathsf{\Omega }\right\}=\left\{\varnothing ,\left\{{\mathrm{S}}_1\right\}, \left\{{\mathrm{S}}_2\right\},\dots \left\{{\mathrm{S}}_N\right\}, \mathsf{\Omega }\right\}$$

(6)

where, Card (E) = 2Ω = 2N.

The probability assignment is an evidence mass function m defined from 2N to [0, 1] that satisfies two requirements in the Equation (7) as follows:

$$\left\{\begin{array}{c}m\left(\mathrm{A}\right)\ge 0, \forall \mathrm{A}\subset \Omega \\ {\displaystyle \sum _{\mathrm{A}\in \mathsf{\Omega }}}\mathrm{m}\left(\mathrm{A}\right)=1\\ m(\varnothing )=0\end{array}\right.$$

(7)

The measure of credibility as well as that of plausibility are deduced from the mass of evidence by the following Equations (8) and (9):

$$Cr\left(\mathrm{S}\right)={\displaystyle \sum }_{\mathrm{A}\subseteq \mathrm{S}, \mathrm{A}\ne \varnothing }m\left(\mathrm{A}\right) ,$$

(8)

$$pl\left(\mathrm{S}\right)={\displaystyle \sum }_{\mathrm{S}\cap \mathrm{A}\ne \varnothing }m\left(\mathrm{A}\right)=1-m\left(\mathrm{A}\right)$$

(9)

with $Cr\left(\mathrm{S}\right)\le pl\left(\mathrm{S}\right)$.

Where, for every S $\subset$ Ω, Cr(S) corresponds to the amount of information that is contained entirely within the considered subset of S by m. It contains all the knowledge crediting the veracity of this subset. pl(S) corresponds to the amount of information not discretizing S, i.e., all the information contained in the sub-assemblies having an intersection with S. These two functions are considered as the lower and higher probability [121]. Figure 4 shows graphically the relationships between Cr, Pl and uncertainty.

Parallel combination of classifiers is essentially based on the errors of individual classifiers, which are usually represented in the confusion matrix for each classifier j given in Equation (10) as follows:

$$P^j=\left(\begin{array}{ccc}{n}_{11}& \cdots & n_{1N}\\ ⋮& n_{ij}& ⋮\\ n_{N1}& \cdots & n_{NN}\end{array}\right)$$

(10)

where, N represents the number of classes. The learning class S_i is represented by row i while the class determined by classifier j is represented by column j.

The confusion matrix for each classifier is subsequently represented. When classifier j selects class ${\mathrm{S}}_k^j$, the masses of evidence will be determined by the following Equations (11) and (12):

$$m_j\left({\mathrm{S}}_k^j\right)=\frac{n_{kk}}{{\sum }_{i=1}^N{n}_{ki}}$$

(11)

$$m_j\left({\overline{\mathrm{S}}}_k^j\right)=1-\frac{n_{kk}}{{\sum }_{i=1}^N{n}_{ki}}$$

(12)

Therefore, the recognition and the confusion probability of classifier j are defined in Equations (13) and (14) as ${\tau }_R^j$ and ${\tau }_C^j$, respectively.

$${\tau }_R^j=\frac{{\sum }_{i=1}^N{n}_{ii}}{{\sum }_{i,j=1}^N{n}_{ij}}$$

(13)

$${\tau }_C^j=1-{\tau }_R^j$$

(14)

The following Equations (15) and (16) based on the classical Dempster’s rule are used to make it more flexible in terms of combining information from n classifiers represented by their mass functions m_i at the same time [122,123,124]. The crucial point for merging is to add the probabilities of all the sets.

$$m\left(\mathrm{S}\right)=m_1\left({\mathrm{S}}_1\right)\oplus m_2\left({\mathrm{S}}_2\right)\oplus \dots \oplus m_j\left({\mathrm{S}}_j\right)+\dots \oplus m_n\left({\mathrm{S}}_n\right)$$

(15)

$$m\left(\mathrm{S}\right)=\left\{\begin{array}{c}\frac{1}{1-K}{\displaystyle \sum _{{\mathrm{S}}_k^1\cap {\mathrm{S}}_k^2\dots \cap {\mathrm{S}}_k^N=\mathrm{S}}}{m}_1\left({\mathrm{S}}_k^1\right).m_2\left({\mathrm{S}}_k^2\right)\dots m_j\left({\mathrm{S}}_k^j\right)\dots m_n\left({\mathrm{S}}_k^N\right)\\ 0 \mathrm{S}=\varnothing \end{array}\right.\mathrm{S}\ne \varnothing$$

(16)

where, K is defined in Equation (17),

$$K={\displaystyle \sum }_{{\mathrm{S}}_k^1\cap {\mathrm{S}}_k^2\dots \cap {\mathrm{S}}_k^N=\mathrm{S}}{m}_1\left({\mathrm{S}}_k^1\right).m_2\left({\mathrm{S}}_k^2\right)\dots m_n\left({\mathrm{S}}_k^N\right)$$

(17)

where, n is the number of masses of evidence or number of classifiers, N is the number of classes or hypotheses, k is the kth hypothesis, K is termed as degree of conflict in the range of [0, 1] to represent the conflict level among multiple mass of evidence, and m(S) is the orthogonal sum of n masses of evidence.

In this paper, for assigning each pixel in the fusing classifiers’ decision rules, we select the class label for which the belief function is maximal, as shown in Equation (18) as follows:

$$\forall x\in { \mathrm{S}}_w, m\left({\mathrm{S}}_w\right)\left(x\right)=max\left\{m\left({\mathrm{S}}_k\right)\left(x\right), 1\le k\le N\right\}$$

(18)

where, S_k (k = 1, 2, …, N) represents N safety levels, S_w is the result of the classification fusion of n classifiers.

The belief function is calculated by means of the Dempster–Shafer combination of Masses of Belief and indicates the belief that each input classification map presents for each label value. Moreover, the masses of belief are based on the input confusion matrices of each classification map, either by using the PRECISION or RECALL rates, or the OVERALL ACCURACY, or the KAPPA coefficient. Thus, each input classification map needs to be associated with its corresponding input confusion matrix file for the Dempster–Shafer fusion.

4. Results

In this study, Sentinel 2A imagery was used because its bands are more adequate for lithological discrimination than other satellite images due to its high spatial and spectral resolution in the VNIR to SWIR range [3]. Preprocessing, spectral subset and spatial resampling were performed on S2A remote sensing imagery using the GEE cloud platform, before attempting lithological classification using spectral and textural features. This was carried out to achieve reflectance data with bands in the same spectral range (VNIR-SWIR) and resample all image bands into a single spatial resolution (10 m). The areas covered by waterbodies, vegetation and terrace alluvium were mapped since the study area is a dry and semi-arid region. The lithological maps produced from each approach were assessed using the overall accuracy, the average accuracy of each class and the Kappa coefficient, derived from the confusion matrix. It should be noted that GEE enables us to solve the challenges concerning the data availability, data storage, data preprocessing and free computing resources, which facilitates the computationally cumbersome feature space used in this research (14 bands of spectral and textural information), in addition to the geological complexity of the study area that is characterized by multiple lithological classes (14 classes), which is time-consuming in local computing mode.

4.1. Selection of Training and Validation Samples

The selection of accurate training and validation samples is one of the most important stages in the lithological mapping process [125]. The geological maps, in situ field verification, data sources, Google Earth and enhanced color composite imagery of the study area allowed us to improve the selection of training and validation polygons. Since the classification approach employed was pixel-based, training and validation samples matched the defined lithological units to which they belonged, hence avoiding mixed pixel issues. A visual approach was used on the reference datasets to collect 480 polygons that were manually selected, with the same number of samples for each category (approximately 30 polygons for each class) distributed over the study area for all the classes in the GEE map area environment, where the reference samples were divided into two sets of 75% training and 25% validation samples.

Figure 5 shows the variation of an average of all training samples derived from the study area for Sentinel 2A bands. It illustrates that vegetation and water features have different spectral properties than lithological units.

4.2. Classification Schemes

Machine learning algorithms, including RF, SVM, CART, MD and NB, were used to classify spectral and textural features of Sentinel 2A, within the GEE cloud computing platform, in the northern part of the complex Paleozoic massif of Rehamna, to identify lithological units. All models were well trained by optimal parameters. The five classification algorithms and the proposed method results for the rocky outcrops are shown in Figure 6.

The classification results obtained using parametric classifiers, particularly CART (Figure 6d), MD (Figure 6e) and NB (Figure 6f), provide more noise than other methods. However, the lithological units comprised of siltstone were well-mapped using parametric algorithms. According to the pre-existing geological maps [66], the lithological units extracted using RF (Figure 6b) and SVM (Figure 6c), non-parametric MLAs, are more precise in detecting bedrock units (pelites, quartzite) and sedimentary rocks (limestone, sandstone, dolomite), respectively. However, the proposed fusion method using DST had high consistency in terms of all rock types exposed in the study area, as shown in Figure 6g.

4.3. Accuracy Assessments of Classified Maps

The accuracy estimates of the lithological mapping results of the ML classifiers available in GEE and the proposed method based on DST are displayed in Figure 7. Among the results obtained using different MLAs, it can be observed that the average value of OAs of non-parametric classifiers, particularly SVM and RF, were higher than those of parametric MLAs, namely, CART, MD and NB. Otherwise, the proposed approach, based on DST, effectively maximizes the classification accuracy by registering values of and for OA and K coefficient, respectively, which indicates that this method is more accurate in recognizing and classifying the lithological units in remote sensing images.

5. Discussion

Unlike the majority of previous studies that adopted individual MLAs using third-party software for lithological mapping of small geographical areas [3,4,13,126,127,128,129], we exploit the capabilities of the GEE cloud computing platform to perform the MLAs procedure using the spatio-spectral data from Sentinel 2. The spatial information was used as a pixel local variability in band measurements (entropy and contrast) contained in the texture index images extracted from the GLCMs. Then, we transferred it to a local computer to train the DST fusion technique by integrating all the MLAs results to produce an optimized lithological map. This is a point of novelty in our study.

The use of the GEE environment offers a variety of open-access geospatial big data and various algorithms (preprocessed spatial models and classifiers), which allow us to create and manage multiple workflows, in order to produce detailed large-scale maps in a cost- and time-efficient manner relative to traditional computing software [54,56]. In this platform, we used five supervised machine learning classification algorithms (RF, SVM, CART, MD, NB) to construct a lithological map in the northern part of the complex Paleozoic massif of Rehamna, based on the spatio-spectral characteristics of Sentinel 2A without needing to purchase or download any software or datasets.

According to the statistical results summarized in Figure 7, the highest accuracies were produced using nonparametric MLAs such as SVM and RF, which were completely in line with the previous studies in this field [4,13,126], compared to parametric MLAs, namely, CART, MD and NB, as they require normally distributed data and statistical parameters are directly inferred from training datasets [125,126]. In contrast, no assumptions or statistical parameters are needed for non-parametric MLAs [130,131]. The classification findings produced in this investigation were compared with pre-existing geological maps [66] and enhanced color composite images to assess the resulting maps acquired in this study. Figure 8 shows a visual comparison between the individual classification results from our research in three complex magnified regions, particularly Machraa ben Aabou, Foum el M’jez and Oued Quibane. It can be seen in Figure 8 that basement rocks, especially pelite (Pt) (Figure 8i,h) and quartzite (Qz) (Figure 8h), were well identified using the RF algorithm. Otherwise, the SVM MLA outperformed in distinguishing different sedimentary cover rocks, namely, clay and sandstone (ClS) (Figure 8j), limestone and dolomite (LD), as well as marl–limestone (ML) (Figure 8k,l). The worst lithological maps were produced when performing classification with parametric algorithms, as illustrated in Figure 8. The mapping results of most rock units exposed in the study area using CART, MD and NB are largely misclassified into surrounding lithological units. Although the three methods achieved low classification accuracy, there were still some lithological units well identified using those algorithms. For example, the sandstone and siltstone (SSt) were well distinguished (Figure 8m,n) and achieved the highest accuracy (Figure 8a) using the CART classifier, in addition to the MD algorithm that yielded good performances in identifying the argillite and siltstone (ASt) as well as the schist and quartzite sandstone (ScQS) (Figure 7a), besides NB shows a better capability for discriminating the psammite and siltstone (PSt) (Figure 7a and Figure 8s). This further illustrated the robustness of the proposed approach in improving classification accuracy in lithological mapping by fusing multiple classifiers. Thus, adding more MLAs increases the probability that the accuracy of the proposed approach will be enhanced and higher than the highest accuracy attained by a single algorithm.

The spectral profile from S2 bands (Figure 5) show that the rocks composed of argillites (ASt) in addition to those composed of silicates, particularly the magmatic rocks (Mg and RMo), are more distinct compared to other lithological units’ spectra; therefore, all the MLAs discriminate them correctly. Otherwise, a lesser separability in the spectral profile is shown by sedimentary rocks composed of the following (

Table 2. Lithological and mineralogical description of the rock units in the study area.

Era	Lithological Units	Mineral Characteristics
Quaternary	Alluvium and terraces	Clay, silt, sand, gravel or similar unconsolidated detrital materials
Meso-Cenozoic	Conglomerate	Silica, calcite or iron oxide
	Sandstone	Quartz sand, feldspar and sometimes silt and clay
	Pelite	Clay minerals (silica, kaolinite, alumina), quartz, feldspar and micas
	Limestone	Calcite and dolomite
	Marl	Calcium carbonate, clay and silt
Paleozoic	Magmatic rocks	Silicates (quartz, feldspars, feldspathoids, colored minerals containing iron and magnesium)
	Carbonate rocks (limestone)	Calcite, dolomite
	Pelite and quartzite	Quartz, feldspar and micas
	Clay and pelitic rocks	Clay minerals (silica, kaolinite, alumina), quartz, feldspar and micas
	Siliclastic and pyroclastic rocks (psammite, sandstone, siltstone)	Quartz and feldspar

): sandy rocks (PSt and SSt); clay (CCl and ClS); calcite and dolomite (LD); in addition to the basement rocks (Pt and Qz), their mineralogical composition is predominantly quartz, which leads to a confused discrimination in individual MLAs, especially parametric classifiers that were affected by the heterogeneity of these classes. Thus, lithological mapping using MLAs depends not only on the data quality and robustness of the classifiers but is also proportional to the nature of lithological units, their mineral compositions, their spectral separability, weatherability and the magnitude of bedrock outcrops relative to other lithological units.

Besides visual comparison with pre-existing geological data, extensive fieldwork has been performed to update and validate the results of the lithological map produced using the proposed approach. Three fieldwork missions were carried out in March 2022, by investigating the mapped area.

The geological map (Figure 9a) and the enhanced color composite (Figure 9c) were used along with the field investigation to validate the proposed approach (Figure 9b) in the region of Machraa ben Aabou. The folds of psammite and siltstone (PSt) and clay and sandstone (ClS) were tuned by the in-situ field verification shown in Figure 9d. Those folds are not mapped in the pre-existing geological map; however, they are observed in the color composite image and well identified in the proposed method as shown in the blue circle in Figure 9a–c. For the cover units, the geological map shows only one rock unit type, which is clay and sandstone (ClS); otherwise, the color composite image and field verification (Figure 9e) confirmed the presence of two cover-type lithologies in the resulting map, which are marl–limestone (ML) and clay and sandstone (ClS), as shown with a yellow circle in Figure 9a–c.

Furthermore, the geological map (Figure 10c), the enhanced color composite (Figure 10e) and the field survey performed along the complex region of Foum el M’jez (shown in Figure 10) were used to validate the optimized lithological map results using the DST approach. The verification locations show considerable compatibility with the predicted lithologies.

Figure 11 shows the in-situ field verification of the basement outcrops of the complex region of Oued Quibane. Additionally, through the field survey, the geological map (Figure 11a) and the enhanced color composite (Figure 11c) have been provided to validate the lithological map produced using the proposed approach. Most of the lithological units exposed in the study area are somewhat consistent and accurate; otherwise, the quartzite (Qz) exposed on surface lithology has been replaced by Quaternary/terraces (QT), as marked by red circles in Figure 11a–c, as the highest overall accuracy of the Qz lithological units achieved using RF (80%) is lower than the highest accuracy of the QT (89%) obtained using SVM MLA (Figure 7a).

Summarily, lithological classification of a large complex semi-arid region, such as the northern part of the Paleozoic of Rehamna in the Moroccan Meseta, is impossible to process by commonly used traditional classification software packages. Therefore, this study exploits the capabilities of cloud-based platforms jointly with local computer software to optimize lithological mapping. Based on spatio-spectral information of S2A imagery, the results reveal that the proposed approach improves the classification performance compared to individual classifiers.

6. Conclusions

Due to the complexity and limited data quality in free satellite-borne remote sensing imagery as well as the restricted performance of single classification algorithms, this study demonstrates a fused multi-classifier (SVM, RF, MD, NB) approach for optimizing lithological mapping using the DST method applied to spectral and textural information of Sentinel 2A imagery. The free cloud-based computing platform, GEE, has been proposed in order to access and process the datasets in a vast part of the study area. This process occurred in the north of the complex Paleozoic massif of Rehamna. The results demonstrate that the classification accuracy is mostly affected by the types of lithological units and their mineral compositions, variations in reflectance spectra, non-linear spectral mixing and sensitivity to weathering, in addition to the classification methods. The resulting lithological maps were compared visually and assessed statistically using the confusion matrix. The proposed DST classified map shows a high accuracy over lithological maps resulting from the other classifiers, which indicates that this approach is more robust and effective in increasing the accuracy in the lithological mapping of complex semi-arid regions due to the complementary contribution of MLAs. In future works, the lithological mapping will be based on a deep learning algorithm, which will be developed within GEE, as cloud computing services are making deep learning more accessible and easier to manage a large-scale computing capacity as well as training algorithms on a distributed architecture, which can be immensely helpful in a geological investigation.