Comparison of Multi-Resolution Optical Landsat-8, Sentinel-2 and Radar Sentinel-1 Data for Automatic Lineament Extraction: A Case Study of Alichur Area, SE Pamir

: Lineament mapping, which is an important part of any structural geological investigation, is made more efﬁcient and easier by the availability of optical as well as radar remote sensing data, such as Landsat and Sentinel with medium and high spatial resolutions. However, the results from these multi-resolution data vary due to their difference in spatial resolution and sensitivity to soil occupation. The accuracy and quality of extracted lineaments depend strongly on the spatial resolution of the imagery. Therefore, the aim of this study was to compare the optical Landsat-8, Sentinel-2A, and radar Sentinel-1A satellite data for automatic lineament extraction. The framework of automatic approach includes deﬁning the optimal parameters for automatic lineament extraction with a combination of edge detection and line-linking algorithms and determining suitable bands from optical data suited for lineament mapping in the study area. For the result validation, the extracted lineaments are compared against the manually obtained lineaments through the application of directional ﬁltering and edge enhancement as well as to the lineaments digitized from the existing geological maps of the study area. In addition, a digital elevation model (DEM) has been utilized for an accuracy assessment followed by the ﬁeld veriﬁcation. The obtained results show that the best correlation between automatically extracted lineaments, manual interpretation, and the preexisting lineament map is achieved from the radar Sentinel-1A images. The tests indicate that the radar data used in this study, with 5872 and 5865 lineaments extracted from VH and VV polarizations respectively, is more efﬁcient for structural lineament mapping than the Landsat-8 and Sentinel-2A optical imagery, from which 2338 and 4745 lineaments were extracted respectively. from the directional ﬁlter image as well as existing faults in the study area. This study demonstrates that the S1A radar images not only are able to provide a high number of lineaments but also can be used for edge enhancement and visual interpretation using directional ﬁltering techniques. These abilities demonstrate advantages of the radar S1A data over optical L8 and S2A images for structural geological studies. Consequently, the framework used in the present study has shown a great efﬁciency in automated lineament extraction using S1A SAR images and can be applied in surrounding areas.


Geographic and Geological Settings of the Study Area
The Pamir plateau, located in the junction zone between the Tibetan-Himalayan and the Tien-Shan orogenic belts, is the western tip of the India-Asia collision zone [28]. It is located between the Tien-Shan Remote Sens. 2019, 11, 778 4 of 29 in the north, the Tarim basin in the east, the Tajik basin in the west, the Kun-Lun mountains and Tibetan plateau in the southeast, and the Karakorum in the south. Today's Pamir plateau formed through the amalgamation (stitching) of several land masses (terranes) that detached from the Gondwanaland (now African continent) drifted northward and accreted to the southern margin of Eurasia during the Mesozoic Era [29][30][31], which subsequently were deformed by the Cenozoic collision between the Indian and Asian plates. These terranes are now marked by the tectonic zones that comprise the Pamir and are named as Northern Pamir, Central Pamir, and Southern Pamir [32]. The terranes of the Pamir are separated from each other and from the surrounding tectonic units by several suture zones (shown by black lines in Figure 1a), which mark the oceanic basins that once existed between these tectonic units and closed as a result of their amalgamation [30]. The general structure of the Pamir plateau is bent northward, which reflects the deformation related to the indentation of the Indian plate following the India-Asia collision, the onset of which is thought to be at about 50 Ma [33]. The deformation of the pre-Cenozoic geological structures within the Pamir plateau by the Cenozoic India-Asia collision is still active and is being caused by the subduction of two continental plates-the Asian and Indian beneath the Pamir from the north and south respectively-thus forming a unique region with an intracontinental subduction [2,28]. This subduction generates many strong earthquakes, some of which are devastating [34,35].
The Southern Pamir, which includes the study area of the present research, extends from Sarez lake in the north to the Afghan-Tajik boundary in the south and from Shahdara Range in the west to the Sarykol Range in the east (Figure 1a). It is often divided into Southwestern Pamir and Southeastern Pamir [28,29,31]. The SW Pamir is characterized by Precambrian metamorphic rocks and the SE Pamir shows a thick Permian to Cenozoic sedimentary succession, stacked into a polyphase Mesozoic-Cenozoic fold and thrust belt. The SE Pamir represents the upper unmetamorphosed crust while the SW Pamir represents the part of the crust that was metamorphosed on the middle and lower crustal levels followed by its exhumation in the Cenozoic [36]. The crystalline rocks and sedimentary deposits of the Southern Pamir are intruded by a series of Mesozoic-Cenozoic igneous plutons, which cover about 30-40% of the surface geological outcrops.
The Cenozoic deformation of the Pamir plateau makes studying the pre-Cenozoic geological evolution of the region and the assessment of the role of the pre-Cenozoic processes in the formation of the plateau difficult, hence making the models of the India-Asia collision less accurate. Therefore, it is important to distinguish the pre-Cenozoic and Cenozoic structures from each other. This can be done by mapping and analyzing the structural lineaments of the plateau with the use of remotely sensed data and image processing techniques, which need to be developed.
For this study, an area in the northern part of SE Pamir was chosen, which is less affected by the Cenozoic deformation compared to the SW Pamir, hence preserving both the pre-Cenozoic and Cenozoic deformation patterns. This area is located approximately 60 km southeast of Lake Sarez (Figure 1a), near Alichur village. It is situated between the following coordinates: Latitudes 38 • 3 41 N and 37 • 47 25 N, Longitudes 73 • 14 18 E and 73 • 36 16 E. According to the values extracted from the digital elevation model, the research area rises from 3733 m up to 5793 m above sea level [37]. Geologically, the study area comprises sedimentary rocks of the terrigenous and carbonaceous composition of the Permian to Jurassic age intruded by a polyphase pluton of dioritic to granitic composition (Figure 1b,c). The study area is characterized by a complex deformation pattern with different stress directions within pre-Jurassic and Jurassic sediments. The Permian to Triassic sediments were mostly folded in roughly the NW-SE and NE-SW striking syncline and anticlines, whereas the Jurassic fold belt in the southern part of the study area has mostly E-W-oriented folds (Figure 1b). Figure 13g,h demonstrates these directions on two rose diagrams that show the strike and dip azimuths respectively. These two deformation patterns could be related to two different deformation events that were previously described as being related to the Cimmerian and post-Cimmerian orogenies in Pamir, which mark the collision of the Southern Pamir with the Central Pamir and the deformation related to the subduction of Neotethys beneath the Eurasian margin during the Jurassic and Cretaceous [29,31]. The faults in the study area also have strike directions consistent with the fold belts, with the Northeast striking faults within the Permo-Triassic and E-W trending faults within the Jurassic deposits. The granite massifs seem to be crosscutting the faults, which might indicate that they intruded after the faults were formed. However, the absence of the extension of the faults into the granite massifs could be due to difficulties in mapping these faults within unlayered homogeneous granites.
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 29 which might indicate that they intruded after the faults were formed. However, the absence of the extension of the faults into the granite massifs could be due to difficulties in mapping these faults within unlayered homogeneous granites.   Voskonyanc (1966); and (c) a geological map of the study area mapped with the use of remote sensing techniques [37].

Data Description
The satellite dataset covering the study area includes a scene of Landsat-8 OLI (L8), a mosaic of two scenes of Sentinel-2A (S2A), and a scene of Sentinel-1A (S1A) Interferometric Wide Mode (IW) [26,38,39]. In this study, a terrain-corrected L8 scene (path/raw = 151/034) covering the study area, which comprises 11 bands (coastal aerosol, blue, green, red, near infrared, two shortwave infrareds, panchromatic, cirrus, and two thermal infrareds) and is free of cloud, was downloaded from the USGS Earth Explorer website. The visible and near-infrared (VNIR) and short-wave infrared (SWIR) (Bands 1-7) channels of the L8 data with a spatial resolution of 30 m were utilized. The acquisition date of the provided L8 data is 15 August 2017. Also, the S2A images with the tile/granule (100 × 100 km) numbers of T43SCB and T43SCC acquired on 14 August 2017 were downloaded. These images possess all the spatial and spectral characteristics corresponding to an S2A Multispectral Instrument (MSI) level 1C (L1C) product. The spectral channels include four bands of 10 m (2, 3, 4, and 8), six bands of 20 m (5, 6, 7, 8a, 11, and 12), and three bands of 60 m spatial resolution (1, 9, and 10). A radar image of S1A Synthetic Aperture Radar (SAR) acquired with the IW swath mode in ascending orbit and dual polarization (Vertical transmit Vertical receive (VV) and Vertical transmit Horizontal receive (VH)) on 21 Aug. 2017 was used in addition to the optical images. The utilized S2A MSI as an L1C product (geometrically corrected top-of-atmosphere reflectance) and S1A SAR C-band (5.405 GHz) imaging and a level-1 high-resolution Ground Range Detected (GRD) were downloaded at free of costs from the Copernicus Open Access Hub website [40]. In addition, a SRTM (Shuttle Radar Topography Mission) 3-arc second elevation data product was used for the orthorectification of radar image as well as for a shaded relief analysis.

Preprocessing and Image Enhancement
The preprocessing stage in digital image processing is an important step, which needs to be done in order to reduce the atmospheric effects that exist in optical remotely sensed imagery [41]. Therefore, this stage was performed on both the L8 and S2A multispectral data to enhance the quality of the raw imagery utilized in this study. Similarly, the radar image should also be preprocessed by speckle filtering and terrain correction (orthorectification) in order to reduce the existing speckle noise and to geometrically correct it so that it reflects the real Earth's surface as close as possible. The abstracted graph of the framework followed in the present study, including the image preprocessing and processing stages, is shown in Figure 2.

Data Description
The satellite dataset covering the study area includes a scene of Landsat-8 OLI (L8), a mosaic of two scenes of Sentinel-2A (S2A), and a scene of Sentinel-1A (S1A) Interferometric Wide Mode (IW) [26,38,39]. In this study, a terrain-corrected L8 scene (path/raw = 151/034) covering the study area, which comprises 11 bands (coastal aerosol, blue, green, red, near infrared, two shortwave infrareds, panchromatic, cirrus, and two thermal infrareds) and is free of cloud, was downloaded from the USGS Earth Explorer website. The visible and near-infrared (VNIR) and short-wave infrared (SWIR) (Bands 1-7) channels of the L8 data with a spatial resolution of 30 m were utilized. The acquisition date of the provided L8 data is 15 August 2017. Also, the S2A images with the tile/granule (100 × 100 km) numbers of T43SCB and T43SCC acquired on 14 August 2017 were downloaded. These images possess all the spatial and spectral characteristics corresponding to an S2A Multispectral Instrument (MSI) level 1C (L1C) product. The spectral channels include four bands of 10 m (2, 3, 4, and 8), six bands of 20 m (5, 6, 7, 8a, 11, and 12), and three bands of 60 m spatial resolution (1, 9, and 10). A radar image of S1A Synthetic Aperture Radar (SAR) acquired with the IW swath mode in ascending orbit and dual polarization (Vertical transmit Vertical receive (VV) and Vertical transmit Horizontal receive (VH)) on 21 Aug. 2017 was used in addition to the optical images. The utilized S2A MSI as an L1C product (geometrically corrected top-of-atmosphere reflectance) and S1A SAR C-band (5.405 GHz) imaging and a level-1 high-resolution Ground Range Detected (GRD) were downloaded at free of costs from the Copernicus Open Access Hub website [40]. In addition, a SRTM (Shuttle Radar Topography Mission) 3-arc second elevation data product was used for the orthorectification of radar image as well as for a shaded relief analysis.

Preprocessing and Image Enhancement
The preprocessing stage in digital image processing is an important step, which needs to be done in order to reduce the atmospheric effects that exist in optical remotely sensed imagery [41]. Therefore, this stage was performed on both the L8 and S2A multispectral data to enhance the quality of the raw imagery utilized in this study. Similarly, the radar image should also be preprocessed by speckle filtering and terrain correction (orthorectification) in order to reduce the existing speckle noise and to geometrically correct it so that it reflects the real Earth's surface as close as possible. The abstracted graph of the framework followed in the present study, including the image preprocessing and processing stages, is shown in Figure 2.

Landsat-8 (L8) Dataset
The geometric correction of the acquired L8 data was not necessary to be performed since it is already geometrically corrected by the United States Geological Survey (USGS) in the Universal  The geometric correction of the acquired L8 data was not necessary to be performed since it is already geometrically corrected by the United States Geological Survey (USGS) in the Universal Transfer Mercator (UTM) projection with the WGS-1984 map datum, Zone 43N. The further preprocessing of the raw L8 data, including radiometric calibration, atmospheric correction, and image normalization, was essential to performed for the dataset to be prepared for further processing. Radiometric calibration and the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction model were applied using the Environment for Visualizing Images (ENVI) software version 5.1 as it was described in the previous study [37]. Image normalization, then, was performed to this dataset in the same environment, and further sub-setting was done to extract the targeted study area. The L8 VNIR and SWIR channels with a 30 m pixel size were fused with the panchromatic channel, having a 15 m pixel size to enhance the spatial characteristics of those with lower spatial resolution. The automatic image fusion was performed using the pansharp module with no generation of a refined pansharpened image option, which is recommended for digital classification purposes [42]. The module used for image fusion is implemented in PCI Geomatica software. This image fusion technique is known as pansharpening and can be applied to the real data types (8,16, and 32 bit). Different image fusion techniques, such as Intensity, Hue, and Saturation (IHS) are also widely used [43,44]. However, the image fusion technique used in this study tends to produce superior sharpening results, while preserving the spectral characteristics of the original images [42].

Sentinel-2A (S2A) Dataset
The S2A images that are used correspond to the L1C products processed by the Payload Data Ground Segment (PDGS). The provided Top-Of-Atmosphere (TOA) reflectance in cartographic geometry (UTM/WGS84) coded in JPEG2000 format data products are already processed for radiometric and geometric corrections, including orthorectification and multispectral and multi-temporal registration on a global reference system with sub-pixel accuracy [38,45]. These L1C images have usually been processed to Level-2A (L2A) data by integrating the Sen2Cor processor in the Sentinel Application Platform (SNAP) [46][47][48][49], which consists of a rich set of visualization, analysis, and processing tools for the exploitation of MSI data from the Copernicus Sentinel-2 mission. However, the L1C data was alternatively scaled to surface reflectance [50,51], using the Dark Object Subtraction (DOS) method [52][53][54][55].
In the present work, the acquired S2A L1C dataset was atmospherically corrected to the surface reflectance, using Sen2Cor atmospheric correction processor version 2.4.0 [56,57]. Sen2cor processor performs a preprocessing of the L1C TOA data product and applies a scene classification with an atmospheric, terrain, and cirrus correction and a subsequent conversion into an ortho-image L2A Bottom of Atmosphere (BOA) reflectance image data [58,59]. As it has been concluded, the higher the spatial resolution, the higher the quality of the lineament map [20]; therefore, a pixel size of 10 m for the resulting level 2A product was chosen in the correction procedure. The output BOA reflectance dataset contained bands from 2 to 7, 8, 8a, 11, and 12. The coastal aerosol-band 1 (60 m), water vapor-band 9 (60 m), and Cirrus-band 10 (60 m) were excluded by the atmospheric correction processor. The next step was to merge the resultant corrected S2A images. The mosaicking of two or more S2A L2A multi-pixel size products with different spatial resolutions in SNAP requires the same spatial resolution for all corresponding channels. Therefore, the SWIR, Red Edge, and near-infrared (NIR) bands with a 20 m spatial resolution were first resampled using the nearest neighbor resampling method to correspond to the VNIR 10 m spatial resolution. The principle of using resampling methods is described by Tempfli, Kerle, Huurneman and Janssen [41]. Further, the datasets were mosaicked in SNAP software version 5.0.0 using Mosaic Processor. Later, the merged in a single image dataset was spatially resized in order to reduce the image size to cover only the targeted study area. Bands 5, 6, and 7 that have been utilized for the detection of vegetation red edge [49] were excluded from the dataset, leaving the VNIR and SWIR channels. For the Interferometry Synthetic Aperture Radar (InSAR) data, the multi-looking of the S1A level-1 GRD data products was not necessary, since according to the slant ranges and azimuthal directions (5 × 1 respectively), these imageries are already multi-looked and have a homogenous number of looks which is equal to four. Moreover, these data are projected to the ground range based on the WGS84 coordinate system-the Earth ellipsoid model by the European Space Agency (ESA). The incidence angle within the provided swath ranges from 31.1 • to 46.5 • . The Sentinel-1 Toolbox (S1TBX) implemented in SNAP software version 5.0.0, developed by the ESA, was used to preprocess the S1A data [60].
In the preprocessing stage, firstly, speckle noise reduction filtering was applied to the S1A SAR data, allowing the preservation of the structure in the image data by filtering the homogeneous surfaces and preserving the edges, which is essential in the field of structural geology. The speckle noises in SAR images are generated by the coherent interference of waves backscattered from the rough surface of the Earth and complicate the image interpretation problem [61]. Therefore, the Refined Lee filter [62], widely used in SAR image speckle effect reduction for better backscatter analysis and visual interpretation, [49,[63][64][65] was applied.
Next, for the image to geometrically represent the real Earth's surface as close as possible and to be ready for processing and interpretation, a geocoding and orthorectification of the SAR data should be performed. To this end, terrain correction, using the Range-Doppler Terrain Correction method, was performed. To do this, the terrain correction operator of the S1TBX SAR Processing module auto-downloaded a Shuttle Radar Topography Mission (SRTM) image with a spatial resolution of 3 arc seconds (90 m), which is recommended for radar data. The resultant output from the preprocessed data contains the amplitude (A) and intensity (I) layers in both the VV and VH polarizations respectively. These two components (A and I) contain information about the amplitude and the phase of the detected signal and are stored in different layers. The (A) layer stands for the position of a point in time on a waveform cycle of a transmitted microwave signal by the sensor, while the (I) layer stands for the measured backscattered signal intensities by the sensor.
The two components derived in the previous step were combined together using the so-called ratioing technique. The ratioing is done by dividing the digital number (DN) values of each pixel in one band by the DN values of another band [66]. This technique enhances the spectral differences between the bands and reduces the shadow effects caused by topography [37,67]. Although the ratioing technique has been used mainly for optical imageries, it was suggested that the same technique can be implemented on the amplitude and intensity bands of the radar data also [41,68]. The slopes in the study area that are facing away from the sensor cannot be illuminated by the radar beam. Therefore, there is no energy that can be backscattered to the sensor, and those regions remain dark in the image [41]. Thus, the ratioing of these two components gave significant and detailed information, with a maximum contrast of lineaments over the study area. Moreover, this combination allowed us to perform the directional convolution filtering on the radar data for visual interpretation and manual lineament extraction.

The LINE (Lineament Extraction) Module
Automatic lineaments extraction is performed on the basis of two main processing steps. The first step is the detection of edges that give information on areas of abrupt changes in the values of neighboring pixels, while the second step is the lines detection [52,69,70]. These steps are often carried out using the LINE module of the PCI Geomatica software, which is a widely used module for automatic lineament extraction [23,71]. There are two categories of parameters applied through this module [20,69] which are listed below and govern the edge detection and curve extraction stages of the automatic lineament extraction procedure: The maximum distance between two vectors to be linked [69].
In the first stage mentioned above, the Canny filter is applied to produce a contour strength image. This edge detection algorithm has been proven to be a better performer than the other edge detection operators and is well-known as a good automatic edge detection technique [12,72,73]. This step is done using the radius provided by the RADI parameter. Further, the gradient is computed from the filtered image. The gradient has been thresholded with the value provided by the GTHR parameter in the second stage to obtain a binary image. In the third stage, the curves are extracted from the binary edge image. According to the value assigned in the LTHR parameter, any curve with a number of pixels smaller than the value provided in this parameter is not taken into consideration for further processing. The extracted curves are, then, converted to a vector form by fitting line segments to them. After that, polylines are generated if they respect the tolerance value defined in the FTHR parameter. Finally, two polylines are linked to form a lineament if they satisfy the criteria specified in ATHR and DTHR parameters [74].

Directional Filtering
Directional filtering is a spatial domain filtering technique and edge enhancement filter that selectively enhances image features having specific direction gradients [75][76][77]. Filtering consists of changing the value of a pixel based on that of its neighbors [78][79][80]. This technique is very useful in radar and optical image processing for edge enhancement, visual interpretation, and manual lineament extraction of geo-structural features using local spatial frequency. This is a widely used technique in geological applications to highlight geological structures and lineaments [8,21,70,81]. The objective aimed in the application of filters by studies was to find the best way for lineament identification and the enhancement of geo-structural information in the images corresponding to structural discontinuities [19,82,83]. Directional filtering is designed to highlight the specific characteristics of an image based on their frequency related to the texture [84]. The algorithm of directional filtering produces the first difference of the image input in the vertical, horizontal, and diagonal directions. As a result, many additional edges of diverse orientations are enhanced. Edge enhancement is performed by convolving the original data with a weighted mask or kernel [8,85,86].

Automatic Lineament Extraction
The first step of automatic lineament extraction was to perform an experimental test using different value combinations for each parameter of the LINE module. The test was performed over the optical and radar datasets to select the optimal parameter settings that produced suitable output information of the lineaments extracted. The application of these parameters over the utilized data was carried out with respect to the literature and validated by visual interpretation. The default values, the proposed values, and the values defined for the LINE module parameters in this study are presented in Table 2. The second step was to select the optimal bands from the two optical datasets (L8 and S2A) that should be used for lineament extraction [20]. Each of the individual bands from these datasets was tested for automatic lineament extraction based on the ability to identify the linear segments and features over the study area using the defined optimal values of parameters. Band 5 (NIR, 1.55-1.75 µm) and band 8 (NIR, 835.1µm) from the L8 and S2A respectively were selected to be used as optimal inputs for lineament extraction. The final prepared images (L8-band 5, S2A-band 8, S1A-VH, and S1A-VV) of the study area were used as the input data for LINE module to produce the final lineament maps.

Manual Lineament Extraction
In order to evaluate the automatically extracted lineaments by the LINE module, a reference lineament map of the study area is required. The reference map for the evaluation is created based on the manual extraction of lineaments [18,87,88]. Some advantages and disadvantages of manual and automatic lineament extraction are discussed in previous studies [1,9]. However, the main advantage of manual extraction is that it is easy to distinguish the type of lineaments, whether the linear features are the tectonic origin or man-made.
Further, all seven bands of L8 and a same number of bands (VNIR and SWIR) of S2A datasets were compared in terms of the contrast and definition of geological structural features. As a result of the visual inspection, band-6 from the L8 and band-11 from the S2A datasets, which record the information at the central wavelengths of 1.6085 mm and 1.6137 mm, respectively, corresponding to the SWIR-1 channel, were selected for this stage of the study. These channels show better contrast on the geological structures compared to the other channels. The selected channels from the two used optical datasets are equivalent to band 5 of Landsat-7 which is widely used for directional filtering in numerous studies [17,81,89,90].
After that, the directional convolution filtering was applied on band 6 of L8 (Figure 3a), band 11 of S2A (Figure 3b), and VH and VV polarizations of S1A (Figure 3c,d respectively) imagery to enhance the geological structural information and to highlight the main lineament directions in the study area. The task was carried out using ENVI software version 5.1. Four principle directional filters: N-S, E-W, NE-SW, and NW-SE with a 7 × 7 kernel size were applied (Table 3). The angles of the directional filter were assigned as follows: N-S by 0 • , NE-SW by 45 • , E-W by 90 • , and NW-SE by 135 • . The kernel size selected in this study was suggested as an optimal kernel size typically used in edge enhancement [91]. The image add-back value was assigned as 60%, which specifies the percentage of the original image that is included in the output image. This part of the original image preserves the spatial context and is typically done to sharpen an image. Later, the filtered images from the N-S, NE-SW, and E-W of optical and radar data were displayed in an RGB color space (Figure 4). The results show that the composed RGB image from the filtered images in the four directions of VH polarization of S1A ( Figure 5) displayed a good contrast on the geological structural features such as folds and faults (Figure 4c), while the composed RGB images of optical sensors mostly highlight the lithological boundaries and with less contrast. Therefore, the RGB image composed of the filtered VH polarization (Figure 6) was selected for further manual lineament extraction. Finally, after the filtering operation and the selection of a most  Later, the filtered images from the N-S, NE-SW, and E-W of optical and radar data were displayed in an RGB color space (Figure 4). The results show that the composed RGB image from the filtered images in the four directions of VH polarization of S1A ( Figure 5) displayed a good contrast on the geological structural features such as folds and faults (Figure 4c), while the composed RGB images of optical sensors mostly highlight the lithological boundaries and with less contrast. Therefore, the RGB image composed of the filtered VH polarization (Figure 6) was selected for further manual lineament extraction. Finally, after the filtering operation and the selection of a most appropriate color composite image from the filtered images, the lineaments were defined and delineated manually with ArcMap software (version 10.3).        (90)) of S1A VH polarization used for manual lineament extraction.

Shaded Relief
To verify the results, the lineaments obtained automatically using the LINE module and manually from the filtered images are spatially correlated against the shaded relief thematic map [52]. Shaded relief thematic maps are a visually pleasing representation of the terrain [92] consisting of gray values stored in raster images [92]. It assumes an illumination of the raster surface given at a defined azimuth and altitude of the sun. The boundaries between the shaded and unshaded areas may indicate the presence of lineaments [22,93]. In addition, the analysis of the shading map extracted from the Digital Elevation Model (DEM) contributes to the characterization of lineaments. Therefore, a shaded relief thematic map was derived from DEM and enhanced by applying a minimummaximum stretch along a color ramp. After comparing the sun's different relative position angles, the azimuth of 0° (north), which indicates the sun's relative position along the horizon, was chosen because it can expose the shaded and unshaded areas better than other values. The altitude angle of 45°, indicating the sun's angle of elevation above the horizon, was chosen. The task was performed using the spatial analyst hillshade tool implemented in the ArcMap software.

Field
A field verification survey was conducted using a Garmin® eTrex® 20 Global Positioning System (GPS) to provide accurate locations for the structural features in the study area. Several sites of geological structural features in the study area were chosen for verification that was extracted from remote sensing data through automatic and manual techniques. A field view of the geomorphology and structural elements photographed are shown in Figure 7. In addition, the extracted lineaments were compared with the faults existing in the geological map of the study area (1:200,000 scale) [94,95].

Shaded Relief
To verify the results, the lineaments obtained automatically using the LINE module and manually from the filtered images are spatially correlated against the shaded relief thematic map [52]. Shaded relief thematic maps are a visually pleasing representation of the terrain [92] consisting of gray values stored in raster images [92]. It assumes an illumination of the raster surface given at a defined azimuth and altitude of the sun. The boundaries between the shaded and unshaded areas may indicate the presence of lineaments [22,93]. In addition, the analysis of the shading map extracted from the Digital Elevation Model (DEM) contributes to the characterization of lineaments. Therefore, a shaded relief thematic map was derived from DEM and enhanced by applying a minimum-maximum stretch along a color ramp. After comparing the sun's different relative position angles, the azimuth of 0 • (north), which indicates the sun's relative position along the horizon, was chosen because it can expose the shaded and unshaded areas better than other values. The altitude angle of 45 • , indicating the sun's angle of elevation above the horizon, was chosen. The task was performed using the spatial analyst hillshade tool implemented in the ArcMap software.

Field
A field verification survey was conducted using a Garmin ® eTrex ® 20 Global Positioning System (GPS) to provide accurate locations for the structural features in the study area. Several sites of geological structural features in the study area were chosen for verification that was extracted from remote sensing data through automatic and manual techniques. A field view of the geomorphology and structural elements photographed are shown in Figure 7. In addition, the extracted lineaments were compared with the faults existing in the geological map of the study area (1:200,000 scale) [94,95].

Results and Discussion
Image preprocessing techniques such as radiometric calibration and atmospheric correction for Landsat-8 (L8) and Sentinel-2A (S2A) satellites, pan-sharpening for L8 and mosaicing for S2A, speckle filtering and terrain correction for Sentinel-1A (S1A) satellite data as well as sub-setting applied on these images allowed us to preliminary enhance and prepare them for further analysis. For the L8 and S2A images, the comparative analysis of the used bands from these datasets in terms of the identification ability of linear features showed that the spatially enhanced band-5 (NIR) of L8 and band-8 (NIR) of S2A are the optimal bands which could be used for the auto-extraction of lineaments. This task was performed using defined optimal parameter values ( Table 2). The number of automatically extracted lineaments from different bands of the L8 and S2A images is shown in Figure  8. The NIR bands of these datasets give higher lineaments number compared with visible and shortwave infrared (SWIR) bands. According to Figure 8a, the biggest number of lineaments (2338)

Results and Discussion
Image preprocessing techniques such as radiometric calibration and atmospheric correction for Landsat-8 (L8) and Sentinel-2A (S2A) satellites, pan-sharpening for L8 and mosaicing for S2A, speckle filtering and terrain correction for Sentinel-1A (S1A) satellite data as well as sub-setting applied on these images allowed us to preliminary enhance and prepare them for further analysis. For the L8 and S2A images, the comparative analysis of the used bands from these datasets in terms of the identification ability of linear features showed that the spatially enhanced band-5 (NIR) of L8 and band-8 (NIR) of S2A are the optimal bands which could be used for the auto-extraction of lineaments. This task was performed using defined optimal parameter values ( Table 2). The number of automatically extracted lineaments from different bands of the L8 and S2A images is shown in Figure 8. The NIR bands of these datasets give higher lineaments number compared with visible and shortwave infrared (SWIR) bands. According to Figure 8a, the biggest number of lineaments (2338) is extracted using band-5 (NIR) of the L8 dataset. Similarly, it can be seen that the highest amount (4745) of extracted lineaments is obtained from band-8 (NIR) of S2A (Figure 8b). Therefore, these bands were considered the optimal bands for further auto-extraction of lineaments besides the two VH and VV polarization images of the radar S1A. is extracted using band-5 (NIR) of the L8 dataset. Similarly, it can be seen that the highest amount (4745) of extracted lineaments is obtained from band-8 (NIR) of S2A (Figure 8b). Therefore, these bands were considered the optimal bands for further auto-extraction of lineaments besides the two VH and VV polarization images of the radar S1A. Further, the steps of automatic lineament extraction were applied over band 5, band 8, and VH and VV polarizations of the L8, S2A, and S1A images respectively, using the defined optimal parameters. The resulting extracted lineaments are then overlain to the classified geological map of the study area as shown in Figure 9a-d. The purpose of overlaying of structural information with the classified geologic map was to analyze their distribution with regards to the different existing lithological units. According to the obtained results, one can see that the lineaments extracted from the S1A VH and VV polarization images are denser and have a high concentration compared to those resulting from the L8 and S2A satellite images. For the L8 and S2A satellites, the concentration of lineaments occurs mainly in the lithological units of the study area. On the other hand, the main concentration of lineaments obtained from the VH and VV images occurs in the shadow areas, steep slopes, dikes, and structural linear features.
The statistical information presented in Figure 10a-d shows that 2338 and 4745 geological lineaments were identified and extracted from the L8 and S2A images respectively, while the S1A VH and VV polarization images provide the higher number of lineaments as 5872 and 5865 respectively. According to these statistics, the values of the lineament length corresponding to the L8 range from 300 m to a maximum of 2360 m (Figure 10a). The values show a range between 200 m and around 2600 m for the S2A image (Figure 10b), while the values for the S1A VH (Figure 10c) vary between 200 m and around 2100 m and the values vary between 200 and around 2550 for the VV polarization images (Figure 10d). The maximum length for almost all four datasets is the same except the S1A VH polarization image which the maximum length is lower from those of the other images, especially the optical images having different natures. This difference shows the usefulness of the VH radar image for the automatic extraction of geological lineaments using the proposed framework.
The most abundant lineament lengths are between 300 and 400 m (1050) for the L8, 200 and 300 m (3704) for the S2A, 200 and 250 m (3131) for the S1A VH, as well as 200 and 300 m (4348) for the S1A VV that occupy 45%, 78%, 53%, and 74% of the total number of automatically extracted lineaments respectively. Comparing these numbers, it can be noted that the lineaments identified from the S1A VH polarization image are mostly characterized by smaller lengths with regards to those extracted from the other images used in this study. It, again, shows the effectiveness of the S1A Further, the steps of automatic lineament extraction were applied over band 5, band 8, and VH and VV polarizations of the L8, S2A, and S1A images respectively, using the defined optimal parameters. The resulting extracted lineaments are then overlain to the classified geological map of the study area as shown in Figure 9a-d. The purpose of overlaying of structural information with the classified geologic map was to analyze their distribution with regards to the different existing lithological units. According to the obtained results, one can see that the lineaments extracted from the S1A VH and VV polarization images are denser and have a high concentration compared to those resulting from the L8 and S2A satellite images. For the L8 and S2A satellites, the concentration of lineaments occurs mainly in the lithological units of the study area. On the other hand, the main concentration of lineaments obtained from the VH and VV images occurs in the shadow areas, steep slopes, dikes, and structural linear features.
The statistical information presented in Figure 10a-d shows that 2338 and 4745 geological lineaments were identified and extracted from the L8 and S2A images respectively, while the S1A VH and VV polarization images provide the higher number of lineaments as 5872 and 5865 respectively. According to these statistics, the values of the lineament length corresponding to the L8 range from 300 m to a maximum of 2360 m (Figure 10a). The values show a range between 200 m and around 2600 m for the S2A image (Figure 10b), while the values for the S1A VH (Figure 10c) vary between 200 m and around 2100 m and the values vary between 200 and around 2550 for the VV polarization images (Figure 10d). The maximum length for almost all four datasets is the same except the S1A VH polarization image which the maximum length is lower from those of the other images, especially the optical images having different natures. This difference shows the usefulness of the VH radar image for the automatic extraction of geological lineaments using the proposed framework.
VH polarization image, which is independent of the soil properties compared to the optical images [52].   The most abundant lineament lengths are between 300 and 400 m (1050) for the L8, 200 and 300 m (3704) for the S2A, 200 and 250 m (3131) for the S1A VH, as well as 200 and 300 m (4348) for the S1A VV that occupy 45%, 78%, 53%, and 74% of the total number of automatically extracted lineaments respectively. Comparing these numbers, it can be noted that the lineaments identified from the S1A VH polarization image are mostly characterized by smaller lengths with regards to those extracted from the other images used in this study. It, again, shows the effectiveness of the S1A VH polarization image, which is independent of the soil properties compared to the optical images [52].
For a detailed mapping of the linear structures in the study area and an evaluation of those provided through an automatic approach, directional filters were applied on the L8 band-6, S2A band-11, S1A VH, and VV polarization images. Based on the results from various color combinations, the first three principal directions (N-S, NE-SW, and E-W) were found to be more suitable for edge enhancement than the other images used in this study. Geological structures are more recognizable after a directional filtering in the VH polarization image. Therefore, the RGB color composite image was displayed to screen assigning the filtering directions derived from the VH polarization image as follows: N-S (Figure 5a) as red, NE-SW (Figure 5b) as green, and E-W (Figure 5c) as blue ( Figure 6). This color composite image was used to extract the geological linear features by visual photointerpretation. [8,81]. As a result, the lineament map (Figure 9e) provided by the directional filtering technique shows a significantly lower number of lineaments compared to those obtained automatically (Figures 9 and 10). However, the lineaments in this map mostly indicate the geological structural features rather than the boundaries of lithological units.
As shown in Figure 10e, in total, 1485 geological lineaments were visually identified and extracted from the filtered color composite image. The values of lineament length show a range between around 200 m and a maximum of around 3900 m. The most dominating lengths of extracted lineaments are between 300 and 500 m (476) which is 32% of the total number of extracted lineaments.
Additionally, by overlaying our automatically and manually extracted lineaments on a shaded relief thematic map, the correlation between the extracted lineaments and shading illumination was analyzed ( Figure 11). The results of this analysis with a visual interpretation show that most lineaments extracted through visual interpretation from the image obtained by applying the directional filter are perfectly located between the shaded and unshaded areas (Figure 11e). It means that the obtained map can be used as a suitable reference map for the validation of those obtained by an automatic method from the radar and optical imagery. Figure 11a,b shows that the lineaments extracted from optical (L8 and S2A) data are mostly located in the boundaries between the lithological units compared to those extracted from the S1A images, which show a much less correlation with these areas (Figure 9a,b). The lineaments extracted from the S1A VH and VV polarization images (Figure 11c,d respectively) are mainly located in the slope and shading and even in the areas with no value classified as a shading area having slops and structural features, thus, indicating the high sensitivity of S1A radar data to geomorphology compared to optical L8 and S2A data.

Density
In lineament analysis and mapping studies, a density classification map is a parameter widely used for a correlation analysis with other geological information [4,20,73,96]. This classification map provides information about the concentration of lineaments per unit area [97]. The lineament density maps in this study are derived for each of the optical, radar, as well as manual interpretation lineament maps to analyze the dispersion pattern of the lineaments (Figure 12). The task was performed using the line density tool implemented in the ArcMap software. The higher density values in the maps are represented by a red color, while the lower values are in blue. Generally, most of the high lineament density values are located in the northern part of the images, which is mainly occupied by granites and granodiorites. These higher densities might reflect the slope angle under the lithological control. However, since granite is more compositionally and structurally homogenous than sedimentary strata that have a heterogeneous structure due to layering, one would expect smooth slopes in the areas of granite distribution and, hence, lower lineament densities. During our field campaign, this was also verified that granite slopes are smoother than the sedimentary strata. Therefore, we assume that the higher density of lineaments in the northern part of the study area could represent a higher level of deformation of rocks in this area. Moreover, the lower densities of the lineaments in the river valleys, which are the domain of quaternary deposits, suggests that deformation that is visible within the granite domain that is older than the quaternary deposits and, hence, covered by them. The analysis shows that the density maps of lineaments extracted through an automatic approach from the S1A radar images (Figure 12c,d) are highly correlated with the reference density map of lineament occurrences derived by visual interpretation (Figure 12e) rather than those from optical L8 and S2A data (Figure 12a,b respectively). The differences in lineament density patterns between optical and radar imageries is most likely related to the illumination effects, spatial resolution, and other factors that give an advantage to radar data over optical imageries.
These density classification maps can be more useful when they are combined with other information by overlaying them on the maps. Therefore, the correlation between the density maps of the extracted lineaments and the known existing faults in the study area is also investigated through digitizing and overlaying them on the derived density maps. The density classification maps were used to find the correlation between the concentration of lineaments and the distribution of existing faults in the study area. Overlaying the existing fault map of the study area on the density maps (Figure 12a-e), it can be seen that the density maps of lineaments extracted automatically from the S1A VH and VV polarization images respectively are highly correlated with the existing faults occurrences, which is also similar with the density map provided by the manual extraction from the directional filter image (Figure 12e). An area in the northern part of the density maps, which exhibits a high lineament density but do not contain any fault, corresponds to the area of distribution of granitic bodies (Figure 9), in which the faults are hard to document due to similar lithological compositions on both sides of the fault lines. This shows that the mapping of some structural elements within the big granitic bodies using image processing techniques is even more efficient than the ground-based surveys. Amongst other obtained maps, the lineament density map of the L8 image (Figure 12a) shows the lowest correlation with the existing fault occurrences. For the optical L8 and S2A satellites, the high-density values are concentrated mainly in the northern part of the study area (Figure 12a,b), whereas the results obtained from the S1A images showed a better correlation with the distribution of the existing faults in the southern part. The directions of the faults correspond with the directions of the high-density patterns in the southern part of these maps. However, the directional correlation of faults and density patterns can also be seen in the lineament density map from the optical S2A image (Figure 12b). This analysis also proves that the results obtained from the S1A radar data are the best correlated in relation to the faults with respect to the conclusion of the shading criteria described previously.

Orientations
The lineament directions extracted using an automated approach can be compared with the directions relating to the manually obtained lineaments and existing faults in the study area [74,96]. Therefore, in the next accuracy assessment section, the lineament directions from the optical and radar data are compared to the directional filter image as well as the existing faults to evaluate the potential of optical L8, S2A, and S1A radar images for automatic structural lineament extraction. This allows the identification of the most frequent directions of lineaments. The lineament directions reflect the result of tectonic processes that produced them. Therefore, they also can be correlated with the geological history of the area to verify the plausibility of the obtained results. Lineament orientation is analyzed by creating rose diagrams for each of the lineament maps including the existing faults and bedding measurements from a geological map, which present the number of lineaments dominating in a particular direction. An angular spacing of 15 • without using the frequency of lineament population was set ( Figure 13). As a result of the analysis, two main dominant directions are identified for each lineament set. These directions for optical data are WNW-ESE (W270-280N) and NNE-SSW (N0-10E) (Figure 13a,b). For the radar data, they are NW-SE (W310-320N) and NNE-SSW (N0-10E) (Figure 13c,d). The lineament set manually extracted from the directional filter reference image shows a trend at the WNW-ESE (W290-300N) and NNE-SSW (N0-10E) directions ( Figure 13e) and for the fault set from the geological map, it is WNW-ESE (W270-280N) and NW-SE (W330-340N) (Figure 13f). Figure 13g,h demonstrates the bulk W-E direction of the strike azimuths of strata and their NNE-SSW dipping directions, respectively, which are digitized from the geological map. The rivers sometime follow the geological structural elements but mostly are the result of erosion processes, which are not necessarily related to deformation. Therefore, it is important to evaluate the contribution of drainage-controlled extracted lineaments in the rose diagram directions of the lineaments. To this end, the drainage network of the study area was extracted from an SRTM digital elevation model using the automatic watershed delineation module of the ArcSWAT model [98,99] built for the ArcMap software. The directions of the drainage network are shown in Figure 13i. As one can see on this rose diagram, the rivers are distributed in four main directions with an angle of approx. 45 degrees between each of them, which is quite natural for rivers. Comparing these directions with the lineament directions extracted from the radar data, we see that the contribution of the drainage network on the directions of the extracted lineaments is small because the lineaments' directions have one dominant direction-NW-SE-and are not homogenously distributed in all directions.
It is interesting to note that one of the dominant directions of the faults (WNW-ESE) correlates well with one of the dominant directions of the lineaments automatically extracted from the optical imagery, while the other dominant fault direction (NW-SE) correlates with one of the lineament-dominant directions extracted from the S1A radar image. As we can see from the geological map (Figure 1b), the faults which are roughly E-W are distributed along the lithological boundaries, which are best extracted from the optical imagery because this type of imagery is more sensitive to compositional changes. This is also supported by the direction of strike azimuths, which is also roughly E-W (Figure 13g). The second dominant faults direction, which correlates with the lineaments extracted from the radar image, may indicate the main direction of the structural geological lineaments, such as faults which are not parallel to the lithological boundaries, and therefore, are not well-extracted from the optical images and best extracted from the radar image. This can also be observed on the geological map (Figure 1b), on which part of the faults striking to NW crosscut the lithological boundaries. The only dominant population of lineaments, the direction of which is the same in the optical and radar images (NNE-SSW), correlates well with the dipping direction of the geological strata and compares them with the topography and geological map; one can conclude that this direction reflects the direction of such features as watersheds and river valleys. Although, we should admit that the uneven distribution of the dipping measurements and the scarcity of the faults mapped does not allow us to excellently correlate the extracted lineaments with the geological features mapped. Therefore, a good correlation between the automatically and manually extracted lineaments is important to back up the verification of automatic lineament extraction results. The orientation analysis indicates that the major trends in the lineament maps obtained through automatic extraction from radar images have a strong correlation with the lineament map derived manually than those from optical imagery. Moreover, the lineaments extracted from the optical imagery are partly controlled by lithological units, which are not always consistent with the deformation pattern. Therefore, although the optical images generate a good number of extracted lineaments, which are parallel to the lithological boundaries, the S1A radar image is best for extracting structural lineaments.

Orientations
The lineament directions extracted using an automated approach can be compared with the directions relating to the manually obtained lineaments and existing faults in the study area [74,96]. Therefore, in the next accuracy assessment section, the lineament directions from the optical and radar data are compared to the directional filter image as well as the existing faults to evaluate the potential of optical L8, S2A, and S1A radar images for automatic structural lineament extraction. This allows the identification of the most frequent directions of lineaments. The lineament directions reflect the result of tectonic processes that produced them. Therefore, they also can be correlated with the geological history of the area to verify the plausibility of the obtained results. (c,d) S1A VH and VV polarizations respectively, and (e) directionally filtered images in comparison to those of (f) the faults digitized from a geologic map, (g,h) the bedding directions of strata-strike azimuths and dipping azimuths respectively, and (i) the drainage network orientations of the study area Lineament orientation is analyzed by creating rose diagrams for each of the lineament maps including the existing faults and bedding measurements from a geological map, which present the number of lineaments dominating in a particular direction. An angular spacing of 15° without using the frequency of lineament population was set ( Figure 13). As a result of the analysis, two main (c,d) S1A VH and VV polarizations respectively, and (e) directionally filtered images in comparison to those of (f) the faults digitized from a geologic map, (g,h) the bedding directions of strata-strike azimuths and dipping azimuths respectively, and (i) the drainage network orientations of the study area We chose several sites of extraction from the S1A satellite image lineaments to verify during our field verification campaign, some examples of which are shown in Figure 7.

Conclusions
The main objective of the study was to compare the performance of optical Landsat-8 (L8) and Sentinel-2A (S2A) in addition to the Sentinel-1A (S1A) radar data in structural lineament extraction. A comparison between the different bands of optical data is also performed to select the best bands for automatic lineament extraction. The results indicate that the best bands for automatic lineament extraction are the fused band-5 (NIR) of L8 and band-8 (NIR) of the S2A. The comparison analysis of automatic lineament detection has been successful with S1A radar images performing better at identifying and extracting more lineament number with a better accuracy compared to the L8 and S2A satellite data. The S1A image resulted in 5872 and 5865 extracted lineaments from the VH and VV components respectively, whereas from the S2A and optical images, 4745 and 2338 lineaments respectively were extracted. This indicates that the optical images are 20% (S2A) and 60% (L8) less effective compared to the S1A VH component in automatic lineament extraction. This success can be attributed to S1A VH and VV polarizations and to the fact that the amplitude and intensity components of this image are combined together. Edge enhancement through directional filtering has also shown its potential to be effective in lineament visualization and their extraction using S1A VH polarization radar image compared to optical L8 and S2A imagery.
The comparison of the gained results, including the shaded relief, density maps, and rose diagrams of the extracted lineaments compared to the directional filter image and existing faults, as well as the field survey showed that the S1A data give a high number of lineaments and better results in terms of accuracy in structural lineament detection than the optical sensors. The lineaments detected from the optical L8 and S2A data follow mainly the boundaries between the lithological units, while those extracted from S1A show a strong correlation with the abrupt changes in areas of slope and shading. In addition, the lineaments extracted automatically from S1A images converge very well with the localization and the orientation of those from the directional filter image as well as existing faults in the study area. This study demonstrates that the S1A radar images not only are able to provide a high number of lineaments but also can be used for edge enhancement and visual interpretation using directional filtering techniques. These abilities demonstrate advantages of the radar S1A data over optical L8 and S2A images for structural geological studies. Consequently, the framework used in the present study has shown a great efficiency in automated lineament extraction using S1A SAR images and can be applied in surrounding areas.