Validation Approach of the Tectonic Lineament Extraction Enhancement Using Sentinel 2A Images and Modified 3 × 3 Bidirectional Prewitt Filters. Case Study: Grombalia, Tunisia

In this study, filters are used to extract directional tectonic lineaments and, thus, to reveal the presence of faults or fractures in the satellite images. The extraction process consists first in eliminating other different types of lineaments (e.g. the lithological limit, ridgelines, hydrographic network, roads, etc.). Besides, a comparative and quantitative approach is applied to show that when the directional Prewitt filter has a zero setting of the opposite of the assumed direction, as is the case of the N-S direction, the East-West direction is also set to zero. This research work, based on a bidirectional filter (N-S; E-W), shows satisfactory results especially concerning the quantitative fluctuation of lineament directions by interval. This fluctuation is in perfect agree-ment with the lineaments provided by the digitizing of the different tectonic accidents directions extracted in a GIS environment from the geological map of Grombalia. In this research work, a quantitative approach was used to evaluate the result of the lineament extraction methodology based on one direction analysis and by an interval of directions. Indeed, the N45 lineament direction, well documented in the geological study of the region, was more clearly distinguished by applying the Prewitt filter than by using the Sobel filter. The result was validated by comparison with the results obtained by the digitization of the tectonic accidents mapped on the geological map.


Introduction
Tectonic lineament mapping facilitates the location of the appropriate site for mineral exploration [1] [2] [3]. It is also very effective in mineral prospection.
Indeed, this technique allows orienting prospection during the reconnaissance phase. In hydrogeology and hydrology, many recent works ( [4] [5] [6]) have shown that these fractures originated groundwater formation in fissured basement environments and constitute the ultimate zones of water flow. Thus, the investigation of this faults network during dam construction remains essential.
The main objective of this study is to improve the process of tectonic accidents lineaments extraction from optical satellite images, although this extraction could be carried out from other data sources (e.g. geomorphology and morphometric indices) as well as from other types of sensors (e.g. radar images) [7] [8].
This study uses Sentinel 2A multi-spectral optical images with a 10 m resolution and applies a two-step extraction approach.
The first step involves the tectonic accidents mapping based on a geological map constructed relying on surveys and field observations. The mapped accidents are digitized and classified in the form of an attribute database containing, essentially the orientation towards north by applying our algorithm to determine these directions in a GIS environment. The approach of quantifying the directions by unit or by interval is employed to compare with field reality. Indeed, it requires, first, digitizing the tectonic accidents mapped in a GIS environment using ArcGIS software.
It is very important to note that the access difficulties in several sectors of the studied region, because of natural obstacles (rivers, mountains, dense forests, and more often the large surface area of the territories), limits sometimes the identification and the mapping of tectonic accidents and make this approach costly and tedious [9]. The use of aerial photography or satellite images can overcome this uncertainty. Therefore, the objective of this approach is to improve the precision of the automatic lineaments extraction from satellite images in order to have a realistic and large-scale vision that reflects the reality of the tectonic history of the region under study.
The second step consists of testing several filters in an unsupervised manner for the automatic extraction of tectonic lineaments from SENTINEL 2A images.
To do so, our methodology is based on colour compositions, principal component analysis (PCA), mathematical morphology and band ratios, and Sobel and Prewitt directional filters in order to compare the results and to improve the Prewitt 3 × 3 filter selected by a bidirectional approach will process the resulting.

Study Area
The region under study is known as the plain of Grombalia (Figure 1), covering a total area of 598.5 Km 2 in North-Eastern Tunisia and characterized by a chaotic morphology.
This region constitutes a collapsed ditch linked to basement subsidence of more than 500 m and filled by Quaternary deposits [10] characterized by sandy substratum at the top of the series, which becomes marly, with gypsum and saliferous from a depth of 130 m [11].
A very dense fault network marks the region's tectonic history, which is the object of our study.

Data
In order to validate the tectonic lineaments extraction result obtained from satellite data based on a quantitative methodology, we used various image-processing techniques and multisource spatial data, as well as:  Digitized lineaments from the geological map of Grombalia [12] (Figure 2).
 Open source multi-spectral Sentinel 2A image of the study area Downloaded from https://www.usgs.gov/. The used bands (B02, B03, B04 and B08) have a resolution of 10 m and their characteristics are presented in (Table 1).

Problem & Methodology
In this work, the result of the tectonic lineaments extraction from the geological map is used as a reference to estimate, qualitatively and quantitatively, the relia-bility of the tectonic lineaments extraction results obtained from the Sentinel 2A image. We also present the number of the tectonic lineaments directions per class and evaluate the lineament extraction procedure by applying different types of directional filters on Sentinel 2A satellite images. Therefore, a methodological approach was implemented to ameliorate the results of lineament extraction from high-resolution satellite images.
In fact, this approach was used to improve the reliability of the automatic lineament extraction, achieve the reliability of field observations, which are relatively costly in terms of time and processing ( Figure 3).

Digitization of Fault Lineaments and Quantitative Analysis
Before digitizing the tectonic accident lineaments of the entire geological map ( Figure 2), geo-referencing using Datum WGS84 was used. In fact, this datum is generally employed in order to superimpose all the lineaments obtained either from the map or extracted automatically from satellite images and, thus, to compare them quantitatively. Obviously, this phase is very sensitive due to the discrepancy between the satellite and geological lineaments. Some researchers proved that this difference could result from geocoding [13] or the inaccuracy in the processed images after the re-sampling process [14]. In addition, geological maps show more inaccuracies due to scale changes, compared to the automatic extracted lineaments. In the present study, the scale of the map from which the lineaments were obtained is 1:50,000. It can result in errors in X and Y of the order of 5 m.
These errors were added to those caused by the conversions of coordinates from one geodetic system to another. It is worth noting that, in this study, there is a transition from Clarke 1880 Lambert to WGS 84 UTM 32 North Tunisia.
After scanning the faults in the map, curvilinear lines were obtained.
To extract automatically faults, the network of tectonic accidents lineaments was considered as a vectorial object.
A script (introduced in the ArcMap software) that provides rectilinear poly-lines from the digitized curvilinear lineaments ( Figure 4) was applied in order to quantify these poly-lines by 15˚ interval (Table 3)  To obtain the vectors file, the three following steps are applied [15]:  Calculation of the coordinates of the extremities of each lineament.   [15].
 Geometrical determination of the orientation angle.  Determination of the direction of the lineaments.

Calculation of Coordinates of Each Lineament Ends
In the lineaments attribute table, the coordinates of the start (X1, Y1) and end (X2, Y2) points were calculated using a script tool in ArcGIS software.

Geometrical Determination of the Orientation Angle
Once the coordinates were obtained, the orientation angle of the lineament was computed by applying the following formula [16]: The experimental result corresponds to the value of the angle θ expressed in radian. Then, this value was converted into degrees to indicate the orientations towards to North, as follows: Angle in degrees = (Angle in radian*180)/II.

Determination of the Lineaments Direction
The direction values α of the faults must be between 0˚ and 180˚ ( Figure 6).
The north should be considered as the origin of the measurements and the positive direction is clockwise.
If α < 0, then the value of the direction will be equal to the absolute value of the additional angle, as revealed in (Table 2).
In the quantitative analysis performed to identify the major directions ( Figure  6), a table showing the number of lineaments per directions intervals was presented. A direction interval of 15˚ was chosen, as demonstrated in (Table 3).
To highlight the studied directions, a directional rosette was drawn from the quantitative table of lineaments ( Figure 7). Journal of Geographic Information System

Tectonic Accident Lineament Map and Analysis
A quantitative and spatial analysis based on the coordinates and directions of the lineaments of the tectonic accidents was carried out ( Figure 6). The linear density is a very important characteristic that should be further examined [17]. It is determined by the ratio of the total lineament length in meters to the area of the studied zone in square kilometres. Apart from a few places that seem to be significantly less dense, the density map of lineaments of Grombalia shows an almost homogeneous density (Figure 8).
It is worth noting that several rectilinear and curvilinear features emerged from the lineament map.
After transforming all entities into a rectilinear form by applying a script from Arc Toolbox, the length of each lineament was calculated by applying the following formula of Vladimir Levenshtein (1): Using GIS environment, the result of lineaments extraction from geological map allowed the quantitative, the spatial and the statistical analysis (Table 4) necessary to validate the lineaments extraction results obtained from satellite images.
The lineaments lengths range from 62,879 m to 2,666,319 m, with an average of 628,727 m (Table 4).  Contrary to the manual and semi-automatic approaches, the interpreter does not influence automatic extraction. In fact, it rather depends on the application performance and the data provided in the examined image [27].

Principal Component Analysis (PCA)
For geological interpretation purposes, the Principal Component Analysis (PCA) is an effective technique employed to enhance a multi-spectral image [28]. The information in several bands is reduced to a smaller number of components generally representing up to 97% of the total variance of the original data set [29]. Sometimes the information in 5 or 6 bands is minimized by the PCA to only 3 components. This analysis allows the creation of colored compounds of the first three components that constitute an excellent visual interpretation product, which increases the contrast between the various objects on the ground.
In order to obtain the best color composition (Table 5), the optimal index factor technique was used in the experiments.
In addition, the Optimal Index Factor (OIF), (2) is a statistical value that provides the optimal combination of three bands out of all the bands in a satellite image from which a color composite can be obtained (Figure 9). This combination contains the largest amount of information (the maximum value of standard deviations), with the minimum duplication (lowest correlation among the pairs of bands).
To calculate the optimal index factors, a map list containing at least three raster maps was created and a correlation matrix or variance-covariance matrix was calculated for the map list, which provided the standard deviations and correlation coefficients applied to calculate the IFM. Journal of Geographic Information System  The algorithms used to compute the IFM are generally executed according to the following formula in GRASS GIS manual (https://grass.osgeo.org/): The calculation of the OIF was performed by the (ILWIS) software that integrated the GIS platform with remote sensing for vector and matrix processing the objects used in our study.
In order to improve the Sentinel 2A satellite image of the area under this study, Principal Component Analysis (PCA) was applied.

Application of the Different Filters (3 × 3) (5 × 5) (7 × 7)
Filtering aims at eliminating noise in data. In geology, scientists focus on finding discontinuities in the textures present in the images (e.g. the contours of relatively homogeneous zones that can reveal the presence of faults or fractures. Lineaments are enhanced by highlighting strong reflectance transitions in the image and the associated spatial high frequencies. Then, the numerical value of a pixel is modified according to its relationship with the values of neighbouring pixels obtained by spatial convolution filtering. The size of the window in pixels is directly proportional to the order of the magnitude of the spatial changes related to the lineaments to be detected.

Other researchers ([9] [10] [11] [32]) retained 7 × 7 windows (
The two tables presented above show that directional filters improve the perception of lineaments by causing an optical shadow effect on the image.
In this approach, these filters were applied in four main directions: N-S, E-W, NE-SW, SE-NW, in order to improve the perception of lineaments not exposed to the illumination source.
Only lineaments with a length greater than half of the convolution window were detected. Therefore, we can conclude that the choice of filter size depends on the length of the lineament to be extracted. In our study, we chose a convolution window of 3 by 3 to detect lineaments with a few meters length (Table 8).
With an angular range of 45˚ of the directional filters around the main enhancement direction, four drift images were first generated from four Sobel directional filters (NO NW, OW, SO SW and S). Subsequently, the modified 3 × 3 Prewitt filters were retained to detect the lineaments in all possible directions.
Given the size of the study area and the order of magnitude of the structures, a 3 × 3 window was used to detect the major lineaments and generate the finer

structures ([23] [24] [32]).
In order to enhance the linear characteristics, directional filters were employed, which improved the perception of the lineaments by causing an optical shadow effect on the image.
All previous works used Sobel's directional filters. However, Prewitt's filters are often compared to Sobel's filters to detect contours.
Thus, the Prewitt filters were applied in the four directions (N-S, E-W, NE-SW and NW-SE) knowing that the NE-SW and NW-SE directions are rarely used in this type of filter. Afterwards, the directional Prewitt filters (Table 9) were modified to obtain better results. Journal of Geographic Information System Table 6. The different 5 × 5 filters.     Table 9. The different adopted 3 × 3 filters.

Extraction of Lineaments with PCI Geomatica
PCI Geomatica's LINE module is generally utilized to extract linear shapes from the image and produce vector files, using 6 parameters to be identified by the user (Table 9 and Table 10). The used parameters are chosen by the photo-interpreter, which results in directed extraction [25].
Indeed, the parameters defining the discontinuities in the filtered images are complex and depend on the length, the angle and the gradient level thresholds of a pixel or a set of pixels considered as a single linear or curvilinear element (Table 10).
The extracted lineaments were processed in vector (.shp) format to conduct the statistical analyses.
The LINE algorithm consists of three steps:  Edge detection;  Thresholding;  Curve extraction.
In the first step, Canny's edge detection algorithm was applied to produce an edge intensity image. The application of this algorithm can be divided into three sub-steps. First, the input image was filtered with a Gaussian function whose radius was provided by the Filter Radius parameter. Second, the gradient was calculated from the filtered image. Finally, pixels whose gradient does not correspond to the local maximum were removed by setting the edge intensity to 0.
In the second step, the edge intensity image was thresholded to obtain a binary image. Each activated pixel of the binary image represents an edge element.
The threshold value was finally defined by the GTHR (Edge Gradient Threshold) parameter.
In the third step, the curves were extracted from the binary edge image. This step can be divided into several sub-stages. First, a thinning algorithm was applied to the binary edge image to produce the pixel skeleton curves. Subsequently, a sequence of pixels for each curve was extracted from the image. Further processing ignored curves containing a number of pixels less than the value of the LTHR (curve length threshold) parameter. After that, an extracted pixel curve was converted into vector form by fitting line segments to it. The resulting polyline approximates the original pixel curve where the maximum fitting error Journal of Geographic Information System (distance between the two) is specified by the FTHR (Line Fitting Threshold) parameter. Finally, the algorithm linked pairs of poly-lines that meet the following criteria:  The two end segments of the two poly-lines are apposite to each other and have similar orientation (the angle between the two segments is less than the value specified by ATHR);  The two end segments are close to each other (the distance between the end points is less than the value of DTHR); The final polylines were stored in a vector segment. Several values of Line PCI module calibration parameters have been used and we notice that the best result has been obtained by applying the following parameters (Table 11).
Four images were produced for each filter related to the 4 directions N-S, E-W, NE-SW and NW-SE. This set of resulting images is used as input for the automatic lineament extraction methods.
By superimposing the four lineament maps obtained in the four directions, a synthesis map of the lineaments was created ( Figure 10).
Automatic methods often generate a lot of insignificant lineaments as well as other types of lineaments, such as roads, tracks, outcrops, river systems, etc., these shapes have been identified and eliminated.

Results and Discussion
After applying the same script applied to the digitized lineaments, rectilinear poly-lines were obtained. Then, the directions of the lineaments by class were generated, which provided quantitative data that can be used for lineaments analysis.  Linking Distance Threshold 100 Table 12 indicates the number of lineaments per 15˚ direction interval.
Then, from the quantitative table of lineaments, a directional rose diagram was created for each type of filter ( Figure 11). This result shows clearly that the quantity of lineaments extracted from the 3 × 3, 5 × 5 or 7 × 7 filters exceeds the number of those obtained from the digitization method. Therefore, the automatic method allows the extraction of tectonic lineaments as well as other types of lineaments. The statistical analysis reveals that the best result was obtained by using the 3 × 3 filters. However, the results show some qualitative improvement in detecting and highlighting the major directions.
As well as the Sobel directional filters, the Prewitt filters were used in the four directions (N-S, E-W, NE-SW and NW-SE). Afterwards, the opposite direction of the filter was set to zero (Table 9), as it is the case for the bidirectional filter The result of this approach, based on the choice of filters to be adopted, shows that the directions between N30 and N60 represent the major directions family in the study area.
In accordance with digitized lineaments extracted from the geological map, the numbers of lineaments extracted from the Prewitt 3 × 3 filters and the filters The average length of segments obtained by using the 3 × 3 filters is 408.997 m, while the linear density of lineaments is 1507 m/Km 2 (Table 13).

Conclusions
In order to evaluate the different types of the used directional filters, Sobel and Prewitt or Prewitt improved bi-directional 3 × 3 filters specifically dedicated to the extraction of tectonic lineaments, this research work developed a methodological approach based essentially on quantitative comparative analysis.
After various tests, the 3 × 3 directional Prewitt filters were maintained, which allowed detecting two major directional intervals. Then, the directional Prewitt was transformed into a bidirectional Prewitt in order to improve lineaments extraction. This modification allowed the best result to extract tectonic lineaments.
The Despite the results with high benefits obtained by this approach, this study also shows that the automatic extraction of tectonic lineaments from satellite imagery is still deficient and requires improvements to the preliminary filters concerning the road network, the hydrographic network, ridgelines, edge of outcrops and alignment of vegetation.