Validation Approach of the Tectonic Lineament Extraction Enhancement Using Sentinel 2A Images and Modified 3 × 3 Bidirectional Prewitt Filters. Case Study: Grombalia, Tunisia ()
1. 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.
2. Study Area
The region under study is known as the plain of Grombalia (Figure 1), covering a total area of 598.5 Km2 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.
3. 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).
Table 1. Characteristics of sentinel 2A band (https://gdla.org/drivers/raster/sentinel2.html).
Figure 2. Digitized tectonic accident lineaments from geological map of Grombalia region (sheet 29, scale 1:50,000, (National Office of Mines “ONM”, 1971).
4. 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 reliability 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).
Figure 3. Validation of the tectonic lineaments extraction from Sentinel 2A images by applying the modified 3 × 3 dual-directional filters methodology.
4.1. 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 ArcMapsoftware) 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) according to their directions towards the north.
Thus, the directions of the lineaments by class were generated (Figure 5), providing a quantitative tool for the analysis of the lineaments.
This applied transformation (Figure 4) allowed obtaining 810 rectilinear poly-lines from 451 curvilinear and rectilinear poly-lines.
To obtain the vectors file, the three following steps are applied [15] :
• Calculation of the coordinates of the extremities of each lineament.
Figure 4. Transformation of the digitized lineaments into straight polylines in ArcMap.
Figure 5. The principle of lineament orientation calculation [15].
• Geometrical determination of the orientation angle.
• Determination of the direction of the lineaments.
4.1.1. 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.
4.1.2. 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.
4.1.3. 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).
Figure 6. Map of lineaments with degree angle value of each direction.
Figure 7. Rose diagram for the lineaments directions.
Table 2. Attributes table of the tectonic accident lineaments extracted from the structural map of Grombalia.
Table 3. Number of lineaments by class of directions according to north.
4.1.4. Tectonic Accident Lineament Mapand Analysis
A quantitative and spatial analysis based on the coordinates and directions of the lineaments of the tectonic accidents was carried out (Figure 6).
First, a map of the lineaments, indexed by the value of each direction, was constructed. It shows around 810 fractures of variable lengths and is oriented in almost all directions.
The directional rosette, as well as the directional histogram, reveal that the classes of directions ]30, 45], ]45, 60] are the major directions of the lineament lattice considered in our previous works.
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):
(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).
Table 4. Spatial analysis of tectonic lineaments.
4.2. Tectonic Lineaments Extraction by Remote Sensing
The existing approaches applied to extract geological lineaments by remote sensing can be classified into three main categories:
• Manual extraction by photo-interpretation, digital image processing and texture analysis [6], image fusion [18] and spatial filtering ( [4] [5] [6] [19] [20] );
• Semi-automatic extraction ( [21] [22] );
• Automatic extraction ( [23] [24] ) based on morphology [25], cellular neural networks, Hough transform [26] and thresholding.
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].
4.2.1. 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.
Figure 9. Image of the Principal Component Analysis (PCA) of Sentinel 2A with 3 Blue (B02) Red (04) and Near Infrared (08) bands.
Table 5. Optimal Index Factor (OIF) per band combination.
The algorithms used to compute the IFM are generally executed according to the following formula in GRASS GIS manual (https://grass.osgeo.org/):
(2)
where:
: is the sum of the standard deviations (standard deviations) of the combinations of the 3 bands k1, k2, k3;
: denotes the sum of the absolute values of the correlation coefficient.
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.
4.2.2. 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.
Several authors, such as ( [30] [31] ), used 5 × 5 windows on Landsat-TM images with a ground resolution of 30 m × 30 m for structural study (Table 6).
Other researchers ( [9] [10] [11] [32] ) retained 7 × 7 windows (Table 7).
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.
4.2.3. 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
Table 10. Comparison of different calibrations of the Line PCI module parameters.
(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.
5. 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.
Figure 10. Synthetic map of lineaments obtained by 3 × 3 filters.
Table 11. Calibration parameters selected for the Line PCI module.
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).
Table 1 shows that the directions between N30 and N60 represent the major fault directions in the study area. Indeed, the number of lineaments in this range is 477 out of a total of 810, i.e. 58.9% of lineaments.
The result of the automatic lineament extraction (Table 10) shows that the number of lineaments in the same range is 1137 among 4786 lineaments extracted from the 3 × 3 filters; 2179 among 9745 lineaments extracted from 5 × 5 filters; and 2,700 of the 13,000 lineaments extracted from the 7 × 7 filters, representing a percentage of 23.76%, 22.36% and 20.77%, respectively.
Figure 11. Rose diagram obtained from 3 × 3 dual-directional adopted filters.
Table 12. Lineaments extracted from satellite image.
The average length of the lineaments varies between 195 m and 198 m while that of the digitized faults is about 629 m. However, the linear density of the lineaments obtained from the 3 × 3 filters is 1566 m/Km2, and it represents the closest value to the digitized lineaments, which is about 851 m/Km2.
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 (E-W-N and S) or (NO, NW, SE, NE and SW). The provided result prove the efficiency of our work in highlighting the major tectonic lineament directions for both directions family: (N30-45) and (N45-60).
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 adopted in this interval are, respectively, 951 and 651 of a total of 3233 and 2571, representing 29.42% and 25.32% of the total number of lineaments.
The average length of segments is 358.04 m, whereas the linear density of the lineaments obtained by using the Prewitt 3 × 3 filters is about 1934 m/Km2.
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/Km2 (Table 13).
Table 13. Spatial analysis of tectonic lineaments.
6. 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 statistical comparison demonstrated that the bidirectional Prewitt 3 × 3 filter provided the best result. This filter is most efficient in showing the major directions in the study area (30˚, 45˚) and (45˚, 60˚), compared to the Sobel and directional Prewitt filters. However, this modified Prewitt filter had reliability of about 50%, in comparison to the cartographic data.
This methodological approach based on quantitative study confirms that the large number and average length of segments obtained by automatic extraction from Sentinel 2A are relatively exaggerated. Hence, it was necessary to further examine the reliability of the automatic lineament extraction with a restricted classification of minor directions in order to highlight the source of estimation of the allowed excess number of lineaments.
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.