Lithological mapping and automatic lineament extraction using Aster and Gdem data in the Imini-Ounilla-Asfalou district, South High Atlas of Marrakech, Morocco

Lithological and lineament mapping using remote sensing is a fundamental step in various geological studies, as it forms the basis for the interpretation and validation of the results obtained. There were two objectives for this study, applied in the Imini-Ounilla-Asfalou district, South High Atlas of Marrakech region: first, lithological mapping by satellite image processing techniques such as ASTER L1B (hight spectral and spatial resolution), namely Principal Component Analysis (PCA), Minimum Noise Fraction (MNF), as well as the application of three types of supervised classification, namely Spectral Angle Mapper (SAM), Maximum Likelihood (ML) and Minimum Distance (MD), on the visible/near-infrared (VNIR) and short-wave infrared (SWIR) spectral bands of our ASTER image; second, an analysis of the distribution of lineaments by automatic extraction using a Global Digital Elevation Model (GDEM) and the PC1 image derived from the PCA transformation applied to the satellite image. The best results are highlighted by the delineation of new facies in relation to the existing map; after confirmation in the field, all of these facies, which include Eocene, Triassic and Jurassic formations, are represented on the new map. The results of lineaments showed that each of them systematically shows a similarity in terms of concentration and orientation, with four preferential oriented systems: NE-SW, E-W, NNE-SSW and NW-SE. The lineaments mainly follow those of the major fault zones, with high concentrations in the northeast and southwest parts of the study area.


Introduction
The Moroccan High Atlas is an orogenic system considered to be the result of the recent Cenozoic evolutionary tectonic inversion of the Triassic and Jurassic rift systems [1].
It presents the result of a long and complex tectonosedimentary evolution linked to the Mesozoic continental rifting evolution of Pangaea, the subsequent opening of the Atlantic and north western Tethys Oceans (pre-orogenic period [2,3] and subsequently, the Cenozoic convergence between Africa and the European plates that led to a complete tectonic inversion of the striated basins (orogenic period [1,3].
The WSW-ENE trend chain is delimited by the North Atlas fault to the north and the South Atlas fault to the south, which are called the master faults of the striated basins active during the pre-orogenic period and then hosting inversions during the orogenic period. This intracontinental inversion chain is subdivided into three zones: The Western High Atlas, the Eastern High Atlas and the Central High Atlas, which is the field of study of this work. Specifically, the study area includes the districts of Imini, Asfalou and Telouet; it is a semi-arid zone south of the Central High Atlas covering an area covering a 60x60 km swath of the ASTER satellite image (Fig. 1).
However, the works of several authors have studied the power of these ASTER images for the geological mapping of different regions in the world. [4] in the Moroccan High Atlas.

Geological setting
In our study area, the adjacent massifs contain igneous (mainly granite) and metamorphic rocks of pre-Triassic age that were deformed during previous orogenesis. They are unconformably overlain by Cenozoic sedimentary rocks deposited in various sedimentary environments [5]. The region has undergone only slight west-to-east buckling. The brittle accidents affecting the cover are materialized by reverse faults with the following main directions: WNW-ESE (Bouazzer and Assaoud faults) and WSW-ENE (Bulgir, Tighermit and Timkit faults). Sometimes, there are abrupt folds [6] (Fig. 2).

Data and methodology
The multispectral ASTER data used were acquired on July 30, 2016 at 09:17 TU. These data are level 1B, corrected for radiometric stripping errors and geometric distortions (ERSDAC 2001). The scene is geo-referenced in the UTM 29 projection and for the WGS-84 ellipsoid. The ASTER sensor records the reflectance in three spectral bands of the VNIR, six of the SWIR and five of the TIR. Only the VNIR and SWIR data are used in our study.
The multispectral ASTER data set available to us is an important source of geological information. In order to extract and synthesize this information into a final geological map, we have set up a semi-automated methodology for processing multispectral data. In the first part of this methodology ( Figure. 3), we present the multispectral image processing techniques used for the lithological identification of the layers, namely PCA, MNF, ICA, IHS to RBG, band ratio and supervised classification. In the second part, we present the techniques used for the automatic lineament extraction, which was done using the line module algorithm of the PCI Geomatica software for the identification of faulted and fractured networks in the study area based on the DEM data and the PC1 band of the ASTER image. This stereoscopic capability makes ASTER ideal for geological and geomorphological interpretation (NASA ASTER 2004). ENVI (Environment for Visualizing Images) version 5.3, PC Geomatica 2016 and ArcGIS 10.5 software were used for ASTER image processing and GIS layer preparation and map layout, respectively.

Principal Component Analysis
After the pre-processing step, the remote sensing data set to be studied has nine ASTER spectral bands and 10 m of spatial resolution for the VNIR bands. Principal Component Analysis (PCA) is an efficient technique for accentuating a multispectral image for geological interpretation [7]. The nine ASTER spectral bands available are highly correlated (a correlation coefficient between 0.7 and 0.9).
A Principal Component Analysis (PCA) was performed on these bands, and the nine new components obtained were unrolled and ordered by decreasing variance rate ( Table 1). The visual study shows the colored compound of the new components 1, 2 and 3 in RGB (Fig. 4A).
The Cretaceous soils are generally marked by green tones, while the Jurassic is well colored by blue pixels. Additionally, Eocene and Quaternary soils appear from light red to pink in the southern and eastern parts of the region.

Transformation into the Maximum Noise Fraction (MNF)
The MNF transform has been extensively used in multispectral and hyperspectral data for feature extraction, noise whitening and spectral data reduction. PCA and the MNF transform have been used extensively in lithological and alteration mapping [8,9,10]. Rotational MNF transformations are used to determine the inherent dimensionality of image data, to separate out noise in the data and to reduce computational requirements for further processing [11].
The Maximum Noise Fraction transformation is applied to the same nine strips. The choice of the number of main components to keep is, as in PCA, a "delicate step left to the user's discretion". After a visual study of the new components, the first 11 new components appear to be noisefree, and therefore, they are retained for further study. The colored compound of MNF 1, 2 and 4 in RGB was chosen because of its potential in providing information about the discrimination of geological formations. The results of the MNF reports provided results in agreement with the reality in the field (Fig. 4B).
The report raised the geological contours of the Upper Cretaceous to the Eocene with well distinguished hues; additionally, the Ordovician, which is extended east-west, appears in light blue pixels in the southern sections of the two maps in Fig. 4.

Band Ratio
ASTER colored band ratio compounds (4/6, 4/7, 8/9), (4/5, 4/3, 3/4) in RGB were used in the present study and proved to be more or less effective in mineralogical and lithological discrimination. The reports show obvious contacts for some rock units: basaltic intrusions of the Triassic finish have a pale red to light green color (Fig. 5A), while Cretaceous clays, sandstones and gypsum barrens appear clearly in a dark yellow color (Fig. 5B). These ratios did not reveal the geological contours any more than the zones of alteration minerals.
A comparison with other previously published research on ASTER band ratios was made [12] this ratio is (6+9/7+8), (5+7/6), 2∕1] in RGB. As the ASTER ratio (6+9/7+8) is used to map the chlorite-carbonate ensembles, the ASTER ratio (5+7/6) is used to detect minerals such as smectite and illite, and the ASTER ratio (2/1) is used for the identification of ferric iron (Fe 3+ ). The color composition of these ratios has indeed allowed this author to better discriminate between existing lithological units. Clays have a green to yellow color, schists appear in dark blue, altered basalts are represented by a bright red color and rhyolites are represented by a dark red color only.

Supervised Classification
A series of supervised classifications were carried out and then analyzed. The purpose of this study was to map the geological formations in detail from the reduced remote sensing data set.
The supervised classification methods used are of two types: one is based on the spectral signatures of the sampled points (SAM) and the other uses portions of the image known as "Regions of Interest (ROIs)". In our case, Maximum Likelihood (ML) and Minimum Distance (MD) are used.

Classification Supervised by Maximum Likelihood
This classification assumes that the statistics for each class in each band were normally distributed and calculates the probability that a given pixel belongs to a specific class. Unless a probability threshold is set, all pixels were classified. Each pixel was assigned to the class with the highest probability (that is, the maximum likelihood). If the highest probability is smaller than a threshold that we specify, the pixel remains unclassified. ENVI implements a maximum likelihood classification by calculating the following discriminating functions for each pixel in the image [13] ( ) = ( ) − where i = the class; x = n-dimensional data (where n is the number of bands); p (ωi) = the probability that class ωi occurs in the image (assumed the same for all classes); |Σi| = the determinant of the covariance matrix of the data in class ωi; Σi-1 = its inverse matrix; and mi = the mean vector.
In the present study, maximum likelihoods were applied on the ASTER bands in the study area using ROIs, with an overall accuracy of 89.0 % (Figure. 6A).

Supervised Classification by Minimum Distance
Minimum Distance (MD): This technique uses the mean vectors of each endmember and calculates the Euclidean distance from each unknown pixel to the mean.
All pixels are classified in the nearest class unless a standard deviation or distance threshold is specified. This technique was applied on the ASTER bands (Fig. 6B).

Supervised Classification Using Spectral Signature Type and Spectral Angle Mapper (SAM)
The SAM method is a supervised classification approach developed [14]. It considers all spectral bands of the image in an "N-dimensional" spectral space and requires radiometrically and atmospherically calibrated and standardized image data [15,16]. It is based on a physical concept that measures the angular similarity between the spectrum of each pixel of the image and the reference spectra, known as prototype spectra or "endmembers". The latter can be measured directly in the field using a spectroradiometer, just as they can be extracted from the image [11,17,18] (Fig. 6C).
Comparing the two types of classification, it seems that lithological mapping using the region of interest (minimum distance and maximum likelihood) gave a better result in discriminating a maximum of geological facies, and statistically the classification showed an accuracy of 91.3%, 90.1% of the overall precision. On the other hand, even for the SAM classification, the kappa coefficient was calculated at 81.2%, which influences the accuracy of the classification.
Visually, the classification showed that there are zones on the geological map, such as the finitriase-triassic basalts west of the IMINI sector; they belong to the Triassic class, and a spectral signature specific to finitriassic basalts was underlined. Even the Jurassic formations of the northeast are classified in the Cretaceous category: that is why we carried out field missions to validate these results. Thanks to all of the techniques applied in our study area, followed by a field mission, we have improved and increased the precision of the geological map (lithology and lineaments). Consequently, some polygons that are considered as average Lias in the northern part of the map are verified and classified in the Triassic class; this also occurs in the southern and south-western parts. Several polygon confusions (terrains) have been redrawn and updated thanks to the strong potential of the GIS and the ASTER image that we used (Fig. 7).

Fig. 7.
Geological map of the region at 1/100000, with significant precision of some facies distinguished in discontinuous black frames.

Lineament Extraction
Automatic extraction was done using PCI Geomatica software via the "LINE MODULE" algorithm. The geological and topographic map was georeferenced and projected in the coordinate system WGS_1984 UTM Zone 29. Then, the overlay of the spatial dataset in a GIS (ArcGIS 10.6) is an operation that facilitates the validation of the lineaments.
After processing and extraction of the lineaments, followed by the removal of some non-geological lineaments (based on a topographic map covering the study area and also using google earth imagery), we obtained the maps of the major lineaments (Figs. 8A, B and C) with their analyses, namely the directional rosettes and the statistics of the frequency distribution of the lineament lengths (Fig. 8D, E and F) and their spatial densities (Figs. 9A, B and C).
The results obtained from the extraction of lineaments from the PC1 image derived from the satellite image and the Digital Field Module (DEM) make it clear that each map produced shows a strong similarity in terms of concentration and orientation with three preferential systems: NE-SW, E-W and ENE-WSW, with a high concentration in the central and northern parts of the study area.
This comparison of fracture orientation shows a good correlation with previous work in the Central High Atlas [19]. A major orientation of the ENE-WSW, WNW-ESE, NE-SW, and NW-SE lineaments was highlighted based on the extraction of lineaments from the processing of a Landsat ETM satellite image. [20] showed that the direction of the lineaments is from NE-SW to E-W, whereas, the research of Bentaher and co-workers [21] also showed that the majority of the lineaments have a NE-SW and E-W direction based on ASTER L1B image processing.

Conclusion
In this work, an ASTER L1B to GDEM 30 m satellite image was used for lithological mapping and lineament extraction. In view of the results obtained and with reference to the reality of the field, it is observable that the envisaged ASTER L1B image processing methods significantly highlight their potential for lithological mapping and lineament extraction.