Automatic Pattern Recognition of Tectonic Lineaments in Seaﬂoor Morphology to Contribute in the Structural Analysis of Potentially Hydrocarbon-Rich Areas

: This work presents novel pattern recognition techniques applied on bathymetric data from two large areas in Eastern Mediterranean. Our objectives are as follows: (a) to demonstrate the e ﬃ ciency of this methodology, (b) to highlight the quick and accurate detection of both hydrocarbon related tectonic lineaments and salt structures a ﬀ ecting seaﬂoor morphology, and (c) to reveal new structural data in areas poised for hydrocarbon exploration. In our work, we ﬁrst apply a multiple ﬁltering and sequential skeletonization scheme inspired by the hysterisis thresholding technique. In a second stage, we categorize each linear and curvilinear segment on the seaﬂoor skeleton (medial axis) based on the strength of detection as well as the length, direction, and spatial distribution. Finally, we compare the seaﬂoor skeleton with ground truth data. As shown in this paper, the automatic extraction of the bathymetric skeleton allows the interpretation of the most prominent seaﬂoor morphological features. We focus on the competent tracing of tectonic lineaments, as well as the e ﬀ ective distinction between seaﬂoor features associated with shallow evaporite movements and those related to intense tectonic activity. The proposed scheme has low computational demand and decreases the cost of the marine research because it facilitates the selection of targets prior to data acquisition.


Introduction
Submarine geomorphology studies the landforms (i.e., relief) and processes (tectonic, sedimentary, oceanographic, and biological) in the submarine realm, many of which comprise renewable and non-renewable resources in many maritime countries. Resources of importance to such countries include unique ecosystems, fishery resources, freshwater, aggregates, minerals, ocean-driven energy, and hydrocarbons [1]. In fact, ocean basins around the world host an extraordinary number of landforms (mid-ocean ridges, seamount/knoll/guyots, volcanic islands, rift basins, trenches, abyssal plains) where such resources are usually located. In addition, continental shelves, comprising the underwater portion of the continent in relatively shallow waters, are rich in marine life, minerals, and hydrocarbons. Marine geomorphometry deals with the characterization of the seafloor using geomorphometric techniques in digital bathymetric data, that is, terrain attributes such as slope, aspect, (a) marine geomorphology/geomorphometrics and further structural studies; (b) the accurate tracing of tectonic lineaments on the seafloor in large marine areas; (c) the detection of morphostructures on the seafloor related to shallow evaporite movement known as halokinesis.
In order to address the previous aims, bathymetric skeletons (medial axes from bathymetric data) are initially computed based on recently developed multiple filtering and skeletonization algorithms [6,16,18]. Then, the linear and curvilinear elements of the bathymetric skeleton are categorized based on criteria such as the strength of detection, their length, direction, and spatial distribution. Finally, we compare the automatically calculated bathymetric skeleton with ground truth data. Both case studies (the Southern Cretan offshore and the Levantine Basin, Figure 1) present bathymetric skeletons that are indicative of active tectonics. In addition, the bathymetric pattern in the Levantine Basin, southwest of the Eratosthenes Seamount, mostly presents morphological characteristics related to shallow salt movements ( Figure 1).
Evaporites are water-soluble, layered, crystalline sedimentary rocks formed from surface or near surface brines in areas where high solar evaporation rates prevailed in the geological past [22]. Evaporation of seawater initially produces calcium carbonate (CaCO 3 ) followed by calcium sulfate in the form of gypsum (CaSO 4 ·2H 2 O) or anhydrite (CaSO 4 ). The following mineral to form is halite (NaCl), commonly known as rock salt. Halite presents a low density (2160 kg/m 3 ) and is mechanically weak, particularly under heavy overburden loading and high temperatures derived from its progressive burial. Thus, salt can flow like a visco-elastic fluid, even in cases where geologically strain rates are rapid [23]. Differential loading induced by gravitational forces, variable thermal conditions, or salt boundary displacement is the primary force inducing salt movement and usually results in vertical and lateral spreading or gliding along an inclined substratum [23]. Buried evaporite often manages to reach shallow depths below the seafloor, giving rise to specific morphological structures and further lineaments, as shown later in this work. Generally, lineaments are the surface expression of geological features such as fractures; faults; folds; geological contacts; rivers; streams; and other physiographic structures showing different origin, size, age, and depth [24][25][26][27]. They are important in recognizing geological structures on Earth's surface related to the control and distribution of natural and mineral resources, geohazards, geothermal energy, and earthquakes. Tectonic lineaments relate to tectonic processes such as faulting and folding and are used in the production of structural maps.

Materials and Methods
Medium-and high-resolution bathymetric data, from the open-file database EMODNet (European Marine Observatory and Data Network [28]) or compiled from various other sources [29][30][31][32], were primarily used in our research. Medium-resolution bathymetry derived from EMODNet comprises a grid/cell size of 0.25 arc minutes, while high-resolution bathymetry has a grid/cell size of 250 m.

Enhancement of the Seafloor Texture
Firstly, slope and aspect images as well their derivatives are computed. Then, the slope and aspect images are efficiently combined for the computation of the enhancement image F of seafloor texture (e.g., lineaments). The enhancement image F is computed as proposed in [16] and tested using onshore and offshore data [6,[35][36][37][38][39][40]. F(p) = (S 2 (p) × SS(p) × SA(p)) 1/4 is calculated based on the slope S(p), the aspect A(p), the first derivative of the slope SS (p), and the first derivative of the aspect SA(p) at point p (S(p)) of the topographic surface.

Multiple Filtering
In the enhancement image F, a step filter G(a,w) is applied, proposed in [34], to decrease noise and emphasize the linear and curvilinear structures of orientation a and width w.
The filter G(a, w) was constrained to be zero mean. In addition, the filter energy is normalized to one under any orientation and width (filter parameters), so that the responses of different angles and widths are comparable, respectively.
In order to enhance the curvilinear structures under any orientation and width, the image I m is computed by getting the maximum of the corresponding pixel values of I g (a, w).

Pixel Labeling
The preliminary goal of skeletonization is to classify I m pixels into three classes, C 1 , C 2 , and C 3 , with label numbers 1, 2, and 3, respectively: C 1 : The pixels corresponding to curvilinear structures. C 2 : The pixels uncertain to correspond to curvilinear structures. We decide for them in the last step, based on the other classes' classification (C 1 and C 3 ). C 3 : The pixels that do not correspond to curvilinear structures. The definition of the above classes is based on the hysterisis thresholding technique, which has been used on edge detection problem [41]. According to hysterisis thresholding, two thresholds, T l (low) and T h (high), are used to initially classify the pixels into three classes (C 1 , C 2 , and C 3 ). The advantage of this thresholding type is the exclusion of some connected point groups [42]. In the proposed scheme, the thresholds T l and T h are automatically estimated based on the median value of I m values (Med).
• T l is given by the mean value of I m values that reveal a value lower than Med. • T h is given by the mean value of I m values that reveal a value higher than Med.
Let B i be the image of initial pixel classification into classes C 1 , C 2 , and C 3 and I m (p) and mv(p) to denote the value of image I m on pixel p and the median value of the nine pixel-neighborhood of pixel p in I m , respectively.

1.
If I m (p) ≥ T h and I m (p) > mv(p), p is classified to C 1 , as its value is very high compared with the image (I m (p) ≥ T h ) and with its neighbourhood (I m (p) > mv(p)).

2.
Else, if I m (p) ≥ T h or (I m (p) > T l and I m (p) > mv(p)), p is classified to C 2 . This is true if the pixel value is high compared with the image, but it is not high enough compared with its neighbourhood, or reversely.

3.
Otherwise, p is classified to C 3 .
Finally, a region growing-based method is applied on pixels of C 2 class to provide the final pixel labeling into classes C 1 and C 3 . Let B f be the image of final pixel classification into classes C 1 and C 3 .
According to the method, the pixels of C 2 class are classified to C 1 if they are connected to a pixel of C 1 ; otherwise, they are classified to C 3 class.
Image skeletonization (thin curvilinear structure detection) (B t ) is provided, if we change the second rule of classification to class C 2 , removing the case of I m (p) ≥ T h .

Automatically Calculated Bathymetric Skeleton
Image skeletonization aims to enhance shape analysis and data interpretation. Figures 3 and 4 illustrate the results of the multiple filtering and skeletonization for both the Southern Cretan offshore and the area on, east, and southwest of the Eratosthenes Seamount, using the bathymetric data as input. The output of the multiple filtering approach for the Southern Cretan offshore and the area on, east, and southwest of the Eratosthenes Seamount is illustrated in Figures 3a and 4a, respectively. The strong and weak detections are depicted with yellow and blue colors, respectively. Finally, Figures 3b and 4b indicate the outcome of the thin curvilinear structure detection (skeletonization) using colour lines, projected on the original bathymetric relief. The red lines correspond to strongly detected linear and curvilinear elements, while the blue lines correspond to weak ones. The texture of the skeleton and changes in specific segments of the images are shown in Figures 3b and 4b. The differences in the skeleton texture possibly relate to (a) sedimentation rates, (b) the degree of geodynamic activity in a certain area, (c) fluid and mass transfers, (d) slope instability processes, and (e) the activity of gravity (gliding) tectonics. Because of the differences in the skeleton texture, we divided both study areas in sub-regions (A, B, C and A', B', C'), respectively, based on the following criteria: • The strength of the detection using the color scale shown in Figures 3a and 4a ranging from 0 (blue) to 550 (yellow) with no units; • The length of lineaments (short, medium, long); • The direction of lineaments; • The spatial distribution of the lineaments (sparse, medium, dense).  More specifically, the sub-regions A, B, C and A', B', C' (Figure 3a,b and Figure 4a,b) show the following characteristics taking into account the above criteria: • Sub-region A: Located in the Southern Cretan offshore (Ptolemy and Pliny Trenches including Gavdos high), characterized by medium to long in length and medium to sparse distributed lineaments generally trending NE-SW, NW-SE, and E-W, which present high detection strength (Figure 3a,b). • Sub-region B: Located in the Southern Cretan offshore (Strabo Trench and the region west of Gavdos high) characterized by short to medium in length and medium to dense distributed lineaments trending to all directions (Figure 3a,b). The lineaments in the southeastern part of sub-region B present high detection strength, while those in the southwestern part of sub-region B present low detection strength (Figure 3a).

•
Sub-region C: Located in the Southern Cretan offshore (Mediterranean Ridge) characterized by medium to long in length and medium to sparse distributed lineaments. Lineaments, close to the limit with sub-regions B and C, present similar orientation with this limit (Figure 3a,b). In addition, the lineaments near to the limit with sub-regions B and C present high detection strength (Figure 3a). The lineaments away from the limit with sub-regions B and C are generally short in length, medium to sparse distributed, and show no prevailing orientation (Figure 3a,b). Sub-region C': Located east and southeast of the Eratosthenes Seamount characterized by generally weak detected, short in length and medium to sparse distributed lineaments that generally trend E-W to NW-SE (Figure 4a,b).

Comparison with Ground Truth Data
In Figure 5a,b, we present the comparison of the automatically calculated bathymetric skeleton with ground truth data for two selected areas. A significant part of the linear and curvilinear elements of the bathymetric skeleton is expected to correspond to tectonic lineaments, as both areas are strongly influenced by intense geodynamic activity related to different mechanisms. The seafloor topography of the Southern Cretan Margin is strongly associated with transtensional movements in the upper 10-15 km of the crust, while transpression prevails below these depths [43]. Small-scale domes due to Messinian evaporitic intrusions may locally deform the recent stratigraphic successions (Pliocene-Quaternary) of the Southern Cretan Margin. Specifically, Figure 5a presents the automatically detected linear and curvilinear elements for the wide area of Gavdos Rise in South Crete. The red arrows show the most prominent automatically detected tectonic lineaments related to already reported geological faults [43][44][45]. Figure 6a,b present in more detail how we implemented the validation of the automatically detected tectonic lineaments. This figure resulted from the overlay of the geological faults (in red) [43][44][45] with the automatically detected skeleton. According to our previous work [16], the proposed method with linear patterns selection algorithm (LPSA) yields approximately 75% of the automatic detected lineaments to coincide either in location or in direction with the geological faults from ground truth data. The methodology in this work yields approximately 85% of the geological faults (thick red lines) to coincide in both location (less than 0.5 arc-minutes) and direction (less than 10 degrees) with the corresponding part of the bathymetric skeleton, meaning that the recall of the proposed method is about 85% (Figure 6b). Thus, there is a strong correlation between the automatically detected tectonic lineaments from the bathymetric skeleton ( Figure 6) and those reported from previous studies [43][44][45]. In addition, a significant amount of information concerning the presence of additional tectonic lineaments is shown in Figure 5a, possibly related to primary or secondary fracture systems that should be investigated in the future. The deformation of the area southwest of the Eratosthenes Seamount is mostly driven by salt tectonics related to the presence of underlying thick Messinian evaporites triggering processes of regional gravity spreading and gliding [46]. Specifically, the authors of [46], based on seismic and swath bathymetric data for this area, reported the presence of basins anchored within pre-Messinian sequences. This is because of massive salt removal and displacement of the salt/sediment system to areas with smoother slopes. Two examples of such areas, indicated by green line, are shown in Figure 5b. These areas (Figure 5b) are limited by distinctive, linear or curvilinear, continuous, strongly detected lineaments that probably correspond to large-scale geological faults [46,47]. In addition, these areas enclose curved morphological structures, whose distribution and morphology refer to salt flow. This argument is based on the correlation of the lineaments in these two areas (Figure 5b) with the seafloor morphology reported by [6,47,48]. In addition, the authors of [48] reported that, for the Central Red Sea, strongly influenced by salt flow, the morphology of the sedimentary cover on top of the moving salt often includes downslope ridges and troughs, along-slope ridges, areas of rough topography, as well as step like morphology at the flow fronts. The formation of brines in several deeps is the result of the dissolution of the Miocene evaporites upper part, and this can explain the irregular seafloor topography (vertical relief) close to the flow fronts [48]. While the origin of the downslope and along-slope ridges observed on the proposed salt flows cannot be definitely solved, their striking morphological similarity, where present, with features observed on ice glaciers or related to submarine mass movement indicates that sediment flow actually takes place near the Red Sea central trough [48]. In this work, we used the findings of [46][47][48] to interpret and further demonstrate how the automatically detected skeleton can assist the detection of shallow salt movements strongly related to potentially hydrocarbon-rich areas such as the wide area southwest of the Eratosthenes Seamount.

Seafloor Deformation due to Halokinesis and Tectonism
The term "halokinesis" is used for autonomous salt movements due to heterogeneous density stratification in the crust [49,50]. The Mediterranean Basin is considered as the largest intracontinental deep-water basin bordered by Europe and Afro-Arabia, being one the world's largest salt-bearing super giants (Figure 1). Two main salt geological intervals occur in the Mediterranean Basin; Triassic-Lower Jurassic salt and Upper Miocene (Messinian) salt (see Figures 8.5,8.7,and 8.8 in [8]). Examples of seismic reflection (i.e., the most reliable ground truth method for imaging earth's crust), showing overburden and related seafloor deformation due to halokinesis, are shown in Figure 7. Figure 7a,b present migrated seismic sections from line ION-7 crossing the Ionian Basin [51], corresponding to both intrusions of Messinian and Triassic evaporites. Figure 7a images the seafloor and underlying strata in the Ionian Abyssal Plain to reveal that Messinian salt deforms the overlying Pliocene-Quaternary strata. This kind of morphology is the so-called 'cobblestone' morphology [52,53]. The term 'cobblestone' morphology is used when recognizing small-scale, topographic bathymetric features with short wavelengths, a pattern primarily associated with small scale, diapiric movements. Figure 7b presents the seafloor deformation due to diapiric movements of Triassic evaporites in the hanging-wall anticlines of pre-existing thrusts [51,54,55]. Seafloor deformation due to the diapiric movements of Triassic evaporites is clearly characterized by larger scale topographic highs compared with the deformation due to the diapiric movements of the Messinian salts. Figure 7c shows an example from Southern Crete (Ptolemy Trough), where the seafloor is only slightly deformed as a result of halokinesis [43]. In Figure 7d, very moderate seafloor deformation in the form of small-scale ridges and troughs is visible in the NE-SW trending profile located SW of the Eratosthenes Seamount [56]. At this point, we point out that the deformation of the Pliocene-Pleistocene sequence overlying the Messinian salt in salt-rich basins around the Eratosthenes Seamount in the Eastern Mediterranean has been confirmed by abundant studies over many decades [46][47][48]. The bathymetric pattern related to salt presents distinct differences when compared with the bathymetric pattern related to intense tectonic activity. In [6], a comparison of the bathymetric patterns on and north of the Eratosthenes Seamount in Eastern Mediterranean proved that the present seafloor morphology north of Eratosthenes Seamount relates to shallow salt movements, creating a chaotic pattern. On the contrary, the bathymetric pattern on the Eratosthenes Seamount relates to the action of E-W trending geological faults. Going one step forward in this work, we apply multiple filtering and skeletonization to (a) correlate the automatically detected pattern of Gavdos Rise in Southern Crete with ground truth data (Figures 5a and 6) and (b) detect patterns related to shallow salt movements southwest of the Eratosthenes Seamount (Figure 5b). Taking into account [6] and the present work, we report the following: • Well-distinguished groups of lineaments trending to specific directions (Figures 5a and 6) characterize the bathymetric pattern related to areas mostly affected by tectonism. The detection strength of the lineaments for these groups is generally high (red and yellow colors in Figures 3a and 6); Considerable work on the automatic extraction of geological lineaments has been published in the last years. The authors of [19] developed an algorithm that is based on tensor voting coupled Hough transform, aiming to extract geological lineaments from DEM and Landsat 8 OLI for the Loess areas in Baoji, China. Later, the authors of [78] proved that the extraction of lineaments is affected by the differences in the vertical accuracy. Specifically, the number of extracted lineaments as well as their length and density increase when the accuracy of the DEM increases. Recently, the authors of [79] managed to extract geological lineaments for a mine area in China based on wavelet edge detection and tracking by hillshade. However, all three previously mentioned works used land DEMs in the experimental phase supported by an abundance of ground truth data to validate their results. On the contrary, the present work used marine DEMs of considerably large areas in Eastern Mediterranean with limited ground truth data and not of very high resolution. However, multiple filtering and skeletonization seem to be effective for the automatic extraction of tectonic lineaments and salt-related morphostructures in bathymetric data even when bathymetric data are not of very high accuracy.

The Contribution of Multiple Filtering and Skeletonization in Marine Geology and Geophysics
Multiple filtering and skeletonization present a wide application on different problems of computer vision and pattern recognition with significant real-world applications in remote sensing [19,34,[78][79][80][81], security [82], medicine [83], and computer graphics [84], especially for the identification and modelling of structures. Thus, these techniques, as previously shown, can significantly assist the interpretation of bathymetric and geophysical images to recognize seafloor and further subsurface structures. Specifically, multiple filtering and skeletonization can be applied in the following ways: • prior to marine mapping at the stage where already available digitized data are reprocessed and evaluated, aiming to accurately specify the research targets and reduce the time and financial cost of the entire marine research; • in digital bathymetric data acquired through remote sensing techniques or marine surveying.
The combination of multiple filtering and skeletonization techniques to identify tectonic lineaments and patterns related to evaporite movements with seafloor expression presents several advantages:

•
Both are fast and simple methods with linear computational complexity concerning the image area, yielding high accuracy detections, as shown in Figures 3-6 of Section 3.

•
Multiple filtering is scale and orientation invariant [78].

•
The skeleton accurately reflects the shape of the original image, preserving the topological properties (homotopy) and symmetry of the object. It holds that the stability of these methods is connected to the data quality (Figures 3-5).

•
Automatic identification does not require any ground truth data or training process as the popular deep learning (neural networks)-based methods [85], which cannot be applied in most of the cases, owing to unavailable or very limited ground truth data.

•
The detection strength, the length, the direction of lineaments, as well as their spatial distribution are robust criteria for their sorting and further geological interpretation.
The only drawback of the proposed scheme is some possible false detections after the skeletonization process, as skeletonization is sensitive to minor boundary deformations and image noise. Some of them can be automatically removed via skeleton pruning [86], if they concern skeleton parts of well-detected medial axes. In addition, segments with unpredictable and irregular orientation can be automatically detected as outliers, if we compare their orientation with the orientation of the other detected medial axes in a neighborhood. This part can be also included in our future work.

Conclusions
This work demonstrates a fast and effective methodology of bathymetric analysis based on image multiple filtering and skeletonization. The experimental results with bathymetric data from the Southern Cretan offshore and the Levantine Basin (on, east, and southwest of the Eratosthenes Seamount) and comparison with ground truth data indicate the reliable, time-effective, and cost-effective performance of the proposed scheme. As a conclusion:

•
The skeleton, provided by the application of the present scheme, preserves the shape and symmetry of the original lineaments. In addition, the tracing of the lineaments is accurate, taking into account their location, direction, and spatial extend. Obviously, the accuracy of lineaments automatic extraction depends on DEM's accuracy.

•
Multiple filtering and skeletonization are independent in training process, so it works even in cases of limited ground truth data.

•
Robust criteria, such as the strength of the detection, the length, the direction, the spatial distribution, and/or density of the lineaments, are valuable for the interpretation of the bathymetric patterns. Consequently, either seafloor morphology related with intense tectonic activity or evaporite movements with seafloor expression can be recognized and further interpreted to assist the structural analysis of potentially hydrocarbon-rich areas.
Ongoing work targets (a) the automatic erasing of false detections and the combination of this method with deep learning when more ground data are available, and (b) the prediction of possible natural gas or oil deposits in combination with related datasets.
Funding: This research received no external funding.