Unsupervised machine learning and depth clusters of Euler deconvolution of magnetic data: a new approach to imaging geological structures

Abstract

We present a novel approach that determines the location and dip of geologic structures by clustering Euler deconvolution depth solutions using Density-Based Spatial Clustering Applications with Noise (DBSCAN). This method and workflow rely on the association of changes in the location and relationships between Euler depth clusters and cluster boundaries with changes in rock susceptibility. We applied our method to global magnetic and high-resolution aeromagnetic datasets over Phanerozoic-Precambrian zone-bounding faults in west and central Victoria. The architecture of these structures at different scales from this imaging technique is comparable to interpreted 2D seismic reflection data. The results from the global magnetic data resolved the architecture of these structures below 5 km, while the aeromagnetic data used were limited to structural information of faults above 2 km depth. Therefore, this method shows the structural relationship of the west-dipping Avoca Fault that soles into the east-dipping Moyston Fault at a depth of ∼22 km in central Victoria and at a shallower depth of ∼15 km southward beneath the Quaternary basaltic rocks of the Newer Volcanic Province. In the vicinity of the Heathcote Zone, the method resolves the location, dip, and overprinting relationship between faults and extrusive rocks, such as the relationship between the Heathcote and Mount William Faults and the granitic Cobaw Batholith. We show how combining magnetic data at various scales can track faults from the near-surface to deeper roots while avoiding possible over-interpretation. We demonstrate how to optimise the DBSCAN parameters and a sensitivity analysis of how to determine clusters and cluster boundaries that are geologically relevant in the absence of geological constraints. Our technique provides an effective and rapid tool for imaging structures and can supplement complex and expensive imaging techniques to resolve the architecture of structures in complex geologic terrains.

KEYWORDS:

• Fault architecture
• magnetic data
• unsupervised machine learning
• Euler deconvolution

Introduction

Mapping geological structures below the Earth’s surface has consistently posed challenges for geologists and geophysicists. Traditional methodologies, such as drilling and seismic surveys, often involve high costs and significant time investment and carry the potential for environmental disruption. Advancements in geophysical technology now allow us to survey and image these deep-seated geological formations in a much more efficient and non-invasive manner, revolutionising our approach to exploring and understanding Earth’s subsurface architectures.

The upper crust of numerous continents has been effectively imaged utilising potential field data, primarily through gravity and magnetic datasets (Blaikie et al. 2017; Fadel et al. 2015; Ganguli and Pal 2023; McLean et al. 2010; Pedreira et al. 2007; Sreejith et al. 2013; Williams et al. 2009). Collecting and processing potential field data is relatively affordable compared to other methods like seismic reflection transects, which offer detailed crust images but lack wide coverage. The signal in potential field datasets corresponds to the 3D arrangement of subsurface rocks, making it ideal for imaging geological structures (Aitken et al. 2009, 2013; Armit et al. 2014; Blaikie et al. 2014; Jessell et al. 1993; Le Pape et al. 2017). However, inherent ambiguity in the data often necessitates auxiliary a priori information, such as petrophysical or geological information, for refinement.

Machine learning (ML) has recently been applied to resolve ambiguities in interpreting geophysical responses. Recent advances in the Geoscientific application of ML applications include the determination of earthquake positions (DeVries et al. 2018; Mousavi and Langston 2016), earth system analysis (Reichstein et al. 2019), and seismic studies (Hu et al. 2019; Huang et al. 2006; Jia and Ma 2017; Lei et al. 2019; Li et al. 2018; Liang et al. 2014; Malfante et al. 2018; Mandelli et al. 2018). Although ML approaches aren’t yet broadly applied to gravity and magnetic analysis (Bergen et al. 2019; Gupta et al., 2022 ; Yu and Ma 2021), they hold promise for elucidating geological structures at varying scales. This potential stems from the ability to convert such data into the frequency domain so that the wavelengths of gravity and magnetic anomalies correspond to their source bodies. This capability underpins numerous spectral-based techniques that determine the locations and depths of sources linked to potential field anomalies (Chukwu et al. 2018; Ganguli and Pal 2023; Kumar et al. 2021; Phillips 2001; Spector and Grant 1970; Teknik and Ghods 2017).

The global magnetic data compilation is improved within southeast Australia with addition of air-based, ground and ocean data (Maus et al. 2009). These global data compilations capture longer wavelengths from the data, providing images of primary, deep crustal geological structures. Global data compilations are less effective at imaging near-surface geological structures associated with shorter wavelengths. However, when combined with an analysis of high-resolution potential field data that primarily captures shorter wavelength anomalies from the surface to shallow crustal levels allows for tracking near-surface structures to their deeper crustal sources (Guimarães et al. 2014). This permits the examination of the magnetic responses across a broad spectrum of wavelengths over an extensive area.

The Euler deconvolution technique is a widely adopted mathematical approach employed to delineate geological structures on multiple scales (Ravat 1996; Reid et al. 2002; Stavrev and Reid 2009). This technique is applied to gravity and magnetic data, to determine the position, depth, and geometry of underlying geological structures causing perturbations in the Earth’s gravitational and magnetic fields. Euler deconvolution is based on Euler’s homogeneity equation (Thompson 1982). It estimates source location and depth in its cartesian form, which can approximate the dip of geologic structures such as faults and folds (Amaral Mota et al. 2020; De Castro et al. 2014). This is due to its ability to determine subtle changes in magnetisation and density without prior knowledge of rock susceptibility and density (1). (1) $\begin{array}{r}(x-x_0)\frac{\mathrm{\delta T}}{\mathrm{\delta x}}+(y-y_0)\frac{\mathrm{\delta T}}{\mathrm{\delta y}}+(z-z_0)\frac{\mathrm{\delta T}}{\mathrm{\delta z}}+\mathrm{SI}(F-R)=0\end{array}$(1) The Euler method uses a deconvolution window size of x − x₀ to y − y₀ to search the gridded potential field data with local (F), regional components (R), and order of homogeneity (SI, known as the structural index), for reliable depth solutions (z) (Reid et al. 2002). The value of SI depends on the geometry of the source. For point or sphere magnetic source, SI = 3, for thick sheet edge, SI = 1, and SI = 0 for contacts (Reid and Thurston 2014). The Euler deconvolution method was designed to estimate the location and depths of gravity and magnetic anomalies associated with simple geological structures (Fairhead and Williams 2006; FitzGerald et al. 2004; Mikhailov et al. 2003), homogenous magnetic field (Fedi et al. 2015 ) and is well-suited to high-resolution potential field data (Khalil 2016; Mohammadzadeh 2015; Robert and Cooper 2015). The method is less sensitive to noise compared to other techniques, making it more robust in the presence of measurement errors or interference (Reid et al. 2002; Stavrev and Reid 2009). The method has been used successfully to define geologic contacts and faults (Aziz et al. 2013; Cooper and Manzi 2015; Khalil 2016; Mohammadzadeh 2015; Oruç and Selim 2011; Ravi Kumar et al. 2020; Wang et al. 2017), locate igneous intrusions (Shaole et al. 2020; Wu et al. 2020), map suture zones (Bournas et al. 2003; Milano et al. 2016), active seismogenic faults (Mazabraud et al. 2005; Minelli et al. 2016), and identify structures in and beneath sedimentary basins (Anand et al. 2009; De Castro et al. 2014; Fang et al. 2015; Ferraccioli et al. 2009; Sridhar et al. 2017).

There are some difficulties with the Euler deconvolution technique. For example, it requires correct structural index selection, representing the field’s decay rate with distance from the source (Barbosa et al. 1999; de Melo and Barbosa 2017; Reid and Thurston 2014). Like all potential field analyses, depth solutions have inherent ambiguity, and the Euler deconvolution can produce multiple solutions, leading to several possible locations and depths for subsurface structures, making it difficult to pinpoint the exact position of geological features. The method also assumes that the subsurface geological structures have homogeneous petrophysical properties, which is geologically unreasonable.

In this paper, we present a novel approach that uses Euler deconvolution to determine the depth of magnetic sources, and Density-Based Clustering to identify clusters of irregular shapes and densities to determine the locations and dips of faults over large distances and depths. We examined the performance of this method in imaging crustal faults using global and high-resolution magnetic datasets at different spatial scales in southeast Australia. Results from a 2 arc minute global magnetic grid and aeromagnetic data of 40 m grid spacing were analysed and compared against published interpretations of crustal-scale 2D-seismic and geologic data of western and central Victoria, in southeast Australia. Using Euler results in plan view, we demonstrate how this method can be used alongside other potential geophysical algorithms to constrain the geological interpretation of magnetic data at various scales.

Geological setting

From the Early Cambrian to the Middle Triassic, southeast Australia experienced numerous cycles of basin formation and orogenesis in a convergent tectonic setting along the eastern margin of Gondwana (Cawood 2005; Glen and Cooper 2021; Gray et al. 2006; Gray and Foster 2004; Phillips and Offler 2011). Tectonism involved lithospheric extension and tectonic mode switching, leading to transient episodes of regional shortening and crustal thickening (Foster and Gray 2000; Gray and Foster 2004). These rocks are now preserved in the Delamerian and Lachlan orogens (Miller et al. 2002; VandenBerg 1978; VandenBerg et al. 2000).

The Delamerian and Lachlan orogens preserve prolonged sedimentation linked to Neoproterozoic to Cambrian break-up of Rodinia (Merdith et al. 2019, 2021), west-dipping subduction initiation (Cawood 2005; Gray and Foster 2004), back-arc basin development (Betts et al. 2002; Gray and Foster 2004), and transpressional deformation attributed to the ca 464–442sMa Benambran Orogeny and ca 400 Ma Bindian Orogeny (Huston et al. 2016; Phillips et al. 2012; Wilson et al. 2020). Foden et al. (2006) suggested that the Late Cambrian to Early Ordovician Delamerian Orogeny developed in response to west-dipping subduction along the eastern margin of Gondwana.

The early Ordovician to early Carboniferous (ca 485–340 Ma) Lachlan Orogen (Cawood 2005; Cayley 2011; Kemp et al. 2002) is a subduction-accretionary system that includes deep-marine sedimentary rocks, MORB-type tholeiitic volcanic and oceanic boninitic and arc-related volcanic rocks, and an entrained Proterozoic microcontinent (VanDieland: Cayley 2011a; Cayley et al. 2002; Moore et al. 2016; Moresi et al. 2014). Ordovician and Silurian arc-related volcanic rocks, continental fragments, and shallower marine sedimentary rocks are mainly exposed in the east, while deep-marine sedimentary rocks are primarily mapped in the west. The western Lachlan Orogen is divided into the Cambrian metasedimentary rocks of the Stawell Zone, Cambrian – Ordovician metasedimentary and mafic volcanic rocks of the Bendigo Zone, and Ordovician to Upper Devonian metasedimentary rocks of the Melbourne Zone (Fergusson et al. 2013) that overlies the VanDieland Proterozoic basement (Cayley 2011a; Moore et al. 2016).

Within the Delamerian and Lachlan orogens are several prominent north-northeast to north-striking crustal-scale faults. Their geometries have been constrained by geological mapping (Edwards et al. 1998; 2001; Morand et al. 1995, 2003; VandenBerg 1978; VandenBerg et al. 2000), potential field analysis (Cayley et al. 2002; Gibson et al. 2013) and interpretation of regional seismic transects (Cayley 2011b; Korsch et al. 2002; Willman et al. 2010). The key crustal faults delineating the limits of significant tectonic elements within the Delamerian and Lachlan orogens are discussed in the following sections. Our focus is on those faults that partition Proterozoic and Paleozoic rocks, especially the NNE-SSW trending, Late Cambrian Heathcote Fault Zone (Figure 1), where we showcase how our approach effectively delineates the structure of these faults.

Figure 1. Area covered by global magnetic data showing a geologic map of structural zones Delamerian and Lachlan Orogens with faults and major zone bounding faults within central and western Victoria (data source: Geoscience Australia and modified in ArcMap after VandenBerg et al. 2000). Faults around the northern segment of the Heathcote Fault Zone were modified after Edward et al. 2001. The insert Q is where further analysis is performed using aeromagnetic data. Location of profile BB′ across the Newer Volcanic province where the edges of the basaltic extrusives rocks are represented by light grey. The location of previously interpreted seismic lines (1, 2 and 3) and cross-sections of Euler depths (AA′ and BB′) are shown.

Yarramyljup fault

The Yarramyljup Fault is a regional steeply west-dipping thrust fault separating the Glenelg Zone in the west and the Grampians-Stavely Zone in the east (Morand et al. 2003). Cambrian mafic and ultramafic igneous rocks of the Glenelg Zone are exposed in the Yarramyljup Fault’s hanging wall within the Glenelg River Metamorphic Complex. The Cambrian Mount Stavely Volcanic Complex is exposed in the footwall (Morand et al. 2003). The fault’s magnetic fabric indicates an NNW-SSE orientation (Figure 1). Seismic data suggest it soles into a Proterozoic to Cambrian passive and syn-rift sequences of the Adelaide Rift Complex (Preiss 2000) that sit along a detachment that floors the Glenelg and Grampians-Stavely Zones (Cayley et al. 2011).

Magnetic data reveal the Yarramyljup Fault extends undercover to the south, where it becomes segmented and overlain by Devonian rocks of the Rocklands Volcanic Group, late Jurassic mafic to intermediate volcanic rocks of the Casterton Formation and the Quaternary basaltic rocks of the Newer Volcanic Province (VandenBerg et al. 2000; Welsh et al. 2011).

Moyston fault

The Moyston Fault is a major east-dipping crustal boundary that separates sub-greenschist facies Cambrian rocks of the Grampians-Stavely Zone (Delamerian orogen) in its footwall to the west from amphibolite grade Cambrian mafic volcanic and turbidite rocks of the Moornambool Metamorphic Complex in the Stawell Zone (Lachlan orogen) in the hanging wall to the east (Cayley and Taylor 2001; Squire et al. 2006). Seismic reflection data suggest this zone of high-strain rocks extends to depths greater than 30 km (Willman et al. 2010), where it bounds second-order west-dipping antithetic reverse faults that truncate ca 500–490 Ma Cambrian volcanic rocks (Miller et al. 2006) and metaturbidites of the Stawell Zone (Cayley et al. 2011; VandenBerg et al. 2000).

Avoca fault

The west-dipping Avoca Fault generally separates northwest-striking Cambrian metaturbidites in the Stawell Zone from north-striking Ordovician metaturbidites in the Bendigo Zone (Gray and Willman 1991; Offler et al. 1998; Wilson et al. 1992). The fault is defined by a two km-wide zone of poly-deformed Cambrian metaturbidites with small entrained lenses of volcanic rocks (Morand et al. 1995). Rocks in the hanging wall record a more complex poly-deformation evolution (Phillips et al. 2002) compared with a single regional deformation event in the footwall (Gray and Foster 1988; Gray and Foster 2004; Morand et al. 1995), which is inferred to have occurred around 455–440 Ma (Leader et al. 2013; Phillips et al. 2002); based on ⁴⁰Ar-³⁹Ar dating of metamorphic mica within cleavages and during the Benambran Orogeny (Glen et al. 2007). Seismic data suggest that the Stawell and Bendigo zones appear as a 160–250 km-wide V-shaped geometry and have been interpreted as imbricates of back-arc inversion and oceanic crust (Cayley and Taylor 2001; Crawford et al. 2003; Miller et al. 2005; Spaggiari et al. 2003; Willman et al. 2010) (Figure 2).

Figure 2. Structural architecture along cross-section AA′ which covers the Stawell, Bendigo and Melbourne Zones. Modified from previous interpretations of seismic lines 1, 2 and 3 covered during 2006 seismic, shows major lithologic units and faults across these structural zones (Cayley et al. 2011; Korsch et al. 2002).

Heathcote fault zone

The Heathcote Fault Zone (Figure 2) is a major west-dipping crustal boundary with a distinct seismic character that extends to the Moho. It separates the Palaeozoic rocks of the Bendigo Zone in its hanging wall from the Proterozoic rocks of the VanDieland microcontinent and its overlying Ordovician to Devonian sedimentary succession of the Melbourne Zone (Cayley et al. 2011; Moore et al. 2016; Moresi et al. 2014; Willman et al. 2010). The fault zone is defined in the west by the Heathcote Fault and to the east by the Mount William Fault (Edward et al. 2001; Edwards et al. 1998). Movement along the fault resulted in the thrusting of Early Cambrian mafic igneous rocks along the Heathcote Fault Zone during the Benambran Orogeny, with another episode of reactivation during the Middle Devonian Tabberabberan orogeny (Edwards et al. 1998). Blueschist “knockers” along the Mount William Fault indicate proximity to a plate edge (Spaggiari et al. 2003).

Method

General processing workflow

Our novel approach for determining the dip and strike of geologic structures combines Euler deconvolution with unsupervised machine learning (see Figure 3 for workflow). We initially precondition the data by removing the geomagnetic reference field in preparation for preprocessing. Subsequently, the data resolution is evaluated to determine deconvolution parameters.

Figure 3. The flowchart illustrates the stages involved in applying our method. The pathway in obtaining optimised clusters is shorter where there are geological or geophysical constraints.

Unsupervised machine learning algorithms learn to recognise patterns in data without being explicitly told what those patterns are. Unlike supervised learning, where the algorithm is trained using labelled data, unsupervised learning is used when the data is unstructured, or we do not know what we are looking for. This study uses Density-Based Clustering Application with Noise (DBSCAN), a machine-learning algorithm for clustering data points based on their spatial density. DBSCAN is an unsupervised machine learning algorithm that relies on the reachability and connectivity of data points to establish a cluster. Data points are reachable if they fall within a given minimum distance (or radius of the neighbourhood) and are connected if all possible paths connecting two points are within the minimum radius (Figure 5a) (Georgoulas et al. 2013; Marques and Orger 2019; Tran et al. 2013). DBSCAN requires the neighbourhood radius, ε, and the minimum number of points, minPts,as inputs (Ester et al. 1996; Tran et al. 2013).

Other clustering approaches such as partitioning, e.g. K-means and K-center (Gupta et al. 2022; MacQueen 1967), hierarchical, e.g. Chameleon and HDBSCAN (Campello et al. 2013; Johnson 1967) and grid-based, e.g. STING and CLIQUE (Amini et al. 2011; Duan et al. 2012; Yu et al. 2015) perform better in datasets with clusters of linear, spherical or convex shapes with similar data densities (Zhou et al. 2017). DBSCAN is a density-based clustering algorithm that can identify clusters of arbitrary shapes and variable densities and detect outliers in the dataset (Georgoulas et al. 2013; Marques and Orger 2019; Tran et al. 2013), making it suitable for geoscientific applications. An assumption made by applying this method to the Euler depth solution is that distinct clusters only exist around distinct magnetised lithologic units of a minimum thickness equivalent to the neighbourhood radius ε. Hence, the clusters’ shape is associated with the geometry of geologic structures, while the vertical variations or jumps in clusters are distinct fault planes or fold limbs and hinges. The neighbourhood radius ε can be determined from sources such as outcrop exposures, drill-hole data or other geophysical techniques that provide information about the thickness of a stratigraphic unit or dyke, fault planes, fault core or damaged zone.

Oasis Montaj was used to pre-process the data and execute Euler deconvolution. Then, cross-sections of Euler depth solutions were clustered using the density-based clustering algorithm with noise, DBSCAN, implemented using a bespoke Python code.

Potential field datasets

To demonstrate the effectiveness of our approach in determining the location and dip of geological structures across multiple scales and different geologic terrains, we employ two magnetic datasets of different scales in this study. The Earth Magnetic Anomaly Grid version 2 (EMAG2) (Maus et al. 2009), which provides global coverage. We focus on the southern Tasmanides and the onshore region of Victoria. The EMAG2 dataset consists of ship track and airborne magnetic data from various worldwide sources (Maus et al. 2009), including data donated by Australia and New Zealand to the National Geophysical Data Center (NGDC), significantly improving the resolution (from 2 to 3 arc minute, ∼3 km) of magnetic anomalies in our study region (Milligan and Franklin 2004). The resulting data is 4 km upward continued from the geoid, preserving longer wavelengths of more than 300 km. The processing of this data is detailed in Maus et al. (2009). The upward continuation and grid cell size of EMAG2 effectively suppress noisy short-wavelength signals associated with near-surface anomalies and enhances long-wavelength signals linked to deep geologic structures.

We then conducted reduction-to-pole (RTP) (Figure 4) to remove the anomaly asymmetry caused by the inclined Earth’s magnetic field and to locate the anomalies directly above their causative bodies. We initially calculated the RTP on the EMAG2 grid (Figure 4) using an inclination of −68.73° and a declination of 11.73°. Although the resulting data is less detailed than high-resolution aeromagnetic or ground gravity data, it imaged first-order zone-bounding structures (Figure 4).

Figure 4. The wavelength components of a 4 × 4° grid of the Reduced to Pole of (a) Earth Magnetic Anomaly Grid (EMAG2; Maus et al. 2009) and (b) high-resolution aeromagnetic data (HRAM; Geoscience Australia, 2016). (c) Power spectrum plot of (a) and (b) shows three major spectral zones with a strong correlation of both data at long to medium wavelength zones (see Zone A and B).

Also, a relatively high-resolution magnetic data from the Australian National Magnetic Grid over an area covered by a major fault zone in central Victoria, the Heathcote Fault Zone (Poudjom et al. 2019). These data were re-levelled to a baseline wavelength based on the Australia Wide Airborne Geophysical Survey and had an approximate grid cell size of 40 m (Foss et al. 2021). An RTP was calculated using a geomagnetic inclination of −67° and a declination of 11⁰ (Figure 4) All data were georeferenced to MGA zone 54.

Depth resolution of EMAG2 in southeast Australia

In order to determine the acceptable Euler solutions, we compared the depth resolution of the EMAG2 data relative to the high-resolution aeromagnetic (HRAM) grid (Sreejith et al. 2013; Srinivasa and Radhakrishna 2017). By examining the spectral content of 4 × 4° RTP grids of EMAG2 and HRAM (Figure 4c) at different crustal levels. This can be done because the depth, d, and wavelength, ƛ, could be related by d∼2ƛ (Abraham et al. 2015; Chukwu et al. 2018; Ganguli et al. 2022). Figure 4c shows two significant crustal boundaries separating zones A, B and C. Zone A shows that EMAG2 and HRAM strongly correlate at long-wavelength magnetic anomalies (ƛ > 50 km) associated with the lower and middle crusts. In Zone B, magnetic anomalies with wavelengths of (2.5 km > ƛ > 50 km) are linked to mid- and upper crustal anomalies. However, in zone B, EMAG2 and HRAM have a strong-to-medium linear correlation. Here, EMAG2, with lower resolution, plots below HRAM. The relative lower power spectrum in EMAG2 indicates that HRAM contains more spectra content (or information).

Zone C has short wavelength (ƛ < 2.5 km) responses and a vague relationship between EMAG2 and HRAM that reflects upper crustal magnetic anomalies. Due to the scattering of EMAG2 spectra in Zone C, it is evident that its resolution is lower than that of HRAM. This lower resolution means it may not resolve geological structures with a wavelength of less than 2.5 km. Because of this limitation, we have chosen to use the upper limit of Zone C, 2.5 km, as a cut-off wavelength that EMAG2 can resolve geologic structures (see Figure 5).

Figure 5. (a) Principle and optimisation of parameters used for DBSCAN required to generate clusters. The K-NN graph shows the number of nearest neighbourhood points (NN) across geologically constrained epsilon (ε) is sorted from farthest away from the ε in Group A to closest in Group D. Cross-section of raw Euler depth results from EMAG2 along profile AA′ (b) and BB′ (c) provides no structural /lithologic information. Location of the profiles is shown in Figure 1.

Euler deconvolution calculation

To determine the strike and depth of magnetic bodies that illustrate the subsurface structure of the region, the deconvolution window size must be more than half of the expected depth to magnetic source body and large enough to produce reliable and well-clustered depth results (Barbosa et al. 1999; Reid et al. 2014; Whitehead 2012). A deconvolution window size of 40 km was used for the EMAG2 to accurately determine source depth of >20 and 1 km for the aeromagnetic data to determine source depths > 0.5 km. A maximum error in depth locations ($\mathrm{\delta xy}$) of 40% was maintained for all Euler depth solutions (Barbosa et al. 1999; Whitehead 2012).

Integer values were used for the structural index to avoid spurious depth solutions SI (Reid et al. 2014; Reid and Thurston 2014). Non-integer SI values in Euler deconvolution applications imply the degree of homogeneity of the anomalous potential field source changes and are prone to misleading depth solutions (Stavrev and Reid 2007). We have chosen 0 and 1 because they are associated with contact and dyke models linked mainly to magnetic contrast from geologic structures such as kilometre-scale lithologic boundaries, faults and shear zones. Only Euler deconvolution using SI = 1 was done for the aeromagnetic dataset, which was required for machine learning as it represents the geometry of geologic structures of interest in this study.

Theory and concept of the application of DBSCAN to Euler depth solutions from EMAG2

DBSCAN, or Density-Based Spatial Clustering of Applications with Noise, is a widely used unsupervised machine learning algorithm. It groups together points in a dataset based on their proximity to one another, following a density criterion. Notably, DBSCAN excels at handling datasets with noise or outliers, as it can detect and disregard these points during the clustering process.

We used DBSCAN to identify relevant clusters within zones of valid and spurious depth solutions to determine the architecture of faults and lithologic boundaries (see Figure 3). We applied our method to cross-sections of Euler depth solutions derived from EMAG2 across major zone-bounding faults in southeast Australia. Furthermore, we chose an area where we could test the method, both with and without constraints, against a well-understood 2D reflection seismic survey.

Initially, we clustered Euler depth solutions across profile A-A′ along the previously interpreted deep 2D seismic survey 06GA-V1 to 06GA-V3 across central Victoria (Figures 1 and 2) (Cayley et al. 2011). Rock exposures and seismic interpretation suggest the Cambrian mafic volcanic units are thinnest at the western edge of the seismic profile and range from ∼2 to 3 km. Based on this constraint, we used a neighbourhood radius, ε, of 2500 ± 200 m (Figure 5). This constraint limits the number of infinite clusters and allows clustering to be varied by changes in the minimum number of points (K-NN) alone.

We then used constrained clustering parameters derived from profile A-A′ on profile B-B which is across the Newer Volcanic Province (Figure 1). Euler solutions from SI = 1 are used for this stage since they were likely to yield solutions related to two-dimensional structures like faults and dykes (Reid et al. 1990; Reid and Thurston 2014). Clusters are expected to change drastically with slight changes in the number of clustering points. Hence, we performed clustering with predetermined ε over a wide range of neighbourhood points (minNumPoints and maxNumPoints) with a significant jump in their midpoint when plotted in ak-nearest neighbourhood (k-NN) distance graph (Figure 5a).

This approach allowed us to determine first-degree clusters. First-degree clusters have cluster boundaries that persist for several changes within the minNumPoints and maxNumPoints. A step-change in the number of points may be necessary for a limited constraint or high uncertainty in the neighbourhood radius, ε. Subsurface structures interpreted as first-degree clusters have the highest degree of confidence. In contrast, we expect more cluster changes from other possible NumPoints and therefore their resulting cluster boundaries are regarded as less reliable. The range of minNumPoints and maxNumPoints that captures a jump and fits an average of our chosen neighbourhood radius for our Euler depth solutions is between 28 and 36 points (see Figure 5a).

Application of machine learning to Euler depth solutions from the HRAM

To test our approach at the outcrop scale, clustering using unconstrained parameters that are automatically derived from the datasets based on the average distances between data points (Figure 6), for a cross-section of Euler depth results extracted from HRAM across the Heathcote Fault Zone in central Victoria. The parameters required for applying DBSCAN, such as the minimum number of points and the neighbourhood radius, ε, were not geologically constrained but obtained by searching for a noticeable jump in the K-NN graph (Figure 6). A significant jump is observed in a minimum number of points within the ranges of 12 and 19 (Figure 5a). All clusters were analysed and compared to the geology of the fault zone.

Figure 6. Unconstrained parameters used for clustering depth results obtained from HRAM. (a) The K-NN graph shows a noticeable jump over an ε of 800 in a window between 12 and 19-the minimum number of points derived from the raw cross-section of Euler depth solutions shown in (b).

Results

Interpretation of Euler deconvolution (ED) depth results: EMAG2 and HRAM

The results of the Euler depth solutions obtained from the EMAG2 and HRAM datasets provides insight into the geology of fault networks at multiple scales. By examining the depth variations across the two magnetic datasets at different scales, we identify areas where the locations of clusters of depth solutions deviate from major geologic faults.

Euler depth results from EMAG2

In map view, Euler results of structural index SI = 0 yield depth range from 4 to 25 km and cluster along strike of previously known extensive zone bounding faults in southeast Australia (see Figures 1 & 7). Similarly, Euler depth results from SI = 1 in map view cluster along strike and on major faults separating the structural zones, providing two-dimensional geometric information for these faults. Euler depth results from SI = 1 (Figure 7c) are more numerous than depth solutions from SI = 0 (Figure 7a) and can be used to determine the strike and dip of the faults, providing additional information to identify these faults at depth (see Figure 7c and d).

Figure 7. Plan view of Euler depth solutions derived from EMAG2 on the topographic map and their correlation along major zone bounding faults. (a) Euler depth solutions from SI = 0 and their interpretation in (b). (c) Euler depth solutions from SI = 1 and its corresponding interpretation in (d). Euler depth solutions in both cases appear to cluster more along the hangingwall of the major zone-bounding faults.

Interpretation of ED depth results from SI = 0

Euler solutions from SI = 0 provide almost no values within areas of strong magnetic contrast. However, they provide depth values correlating with major faults in areas with little interfering magnetic contrast, such as along the west-dipping Avoca Fault and the northern segment of the east-dipping Moyston Fault (see Figure 7a and b). To the north of the Newer Volcanics Province, linear NNW-SSE depth clusters are aligned with a short segment of the west-dipping Yarramyljup Fault but do not give depth results across the Newer Volcanic Province (Figure 7b and d). It can be difficult to establish the dip of the faults using depth clusters in a plan view.

Interpretation Euler deconvolution depth results from SI = 1

Shallower depth solutions (from <5100 to 17,500 m) from SI = 1 are present west of the Mount William Fault. Depths greater than 17,500 m east of the Mount William Fault correspond to the thicker and older crustal unit beneath the Melbourne Zone (Cayley et al. 2002; Edwards et al. 1998) (Figure 7d). The transition of deeper to shallower depths from east to west continues southward and connects with the Bambra Fault with both shallow and deeper depths (Figure 7d). Deeper clusters along the Selwyn Fault suggest that this fault lies within the relatively thick crust of VanDieland (Berry et al. 2005; Moore et al. 2015).

Depth clusters in the Newer Volcanic Province generally show the orientation of depth clusters from SI = 1 that trend east–west. These Euler depth solutions show the orientations of the young basaltic rocks with much higher magnetic susceptibilities than the basement rocks. Thus, within this zone, the orientation of the depth solutions does not reflect the NNW, NNE and N-S trends of the major faults.

Euler depth results from HRAM

The results of using the high-resolution aeromagnetic data for Euler deconvolution with SI = 1 show clear linear clusters of depth solutions around several faults close to the Heathcote Fault Zone (Figure 8). These clusters are aligned with zones of abrupt changes in the magnetic signatures around the major faults. The depth solutions obtained through Euler deconvolution range from 100 meters to 2 kilometres and correlate with different geophysical responses associated with the various geological structures around the fault zone (Slater and Haydon 1999).

Figure 8. Correlation of Euler depth solutions generated from RTP of the aeromagnetic data using SI = 1. (a) shows the data with no annotation, while (b) shows the locations of the major faults in the vicinity of the Heathcote Fault Zone (locations in Figure 1). Concentric clusters of relatively deep depth values correlate with mappable late Devonian granites G1, G2, G3 and G4 (see Edwards et al. 2001). Profile C-C′ shows where the cross-section of the Euler depth from the aeromagnetic data is located (Figure 9).

Figure 9. DBSCAN clustering of cross-section AA′ using ε = 2500 m and individual K-NN values and grouped according to their mirror distance from the inflexion point on the K-NN graph in Figure 5a. (a) End-member cluster assemblages farthest and (D) is the closest to the inflexion point or optimal cluster assemblages. The distinct cluster assemblages are identified by their colours. The n on each cluster assemblage for individual K-NN values indicates the number of clusters or cluster labels. All depth solutions are projected to GDA94/MGA 54.

The NNE and NNW trending linear clusters runs along the major faults, including faults within the Heathcote Fault Zone (Figure 8). In plan view, linear clustered depth solutions align within the hangingwall of the Heathcote and the Mount William Faults. Linear clusters of depth values from 200-750 m are aligned with the hanging wall of the west-dipping Mount William Fault, the eastern-most fault in the zone that separates Ordovician to Devonian sedimentary rocks of the Melbourne Zone in the east from Cambrian to Ordovician turbidite rocks of the Bendigo Zone to the west (Edwards et al. 1998; 2001). Linear depth clusters are also aligned along the west-dipping Heathcote Fault, the westernmost fault segment in the fault zone in this area (Figures 2 and 8b). Other NNW and NNE linear depth clusters lie along the west-dipping Meadow Valley Fault and west-dipping Fosterville Fault to the west of the Heathcote Fault Zone (Figure 8b).

Concentric clusters with relatively deeper solutions at their centres lie around major granite intrusions. Depth range from 750 to 2000 m southwest and southeast of the Heathcote Fault Zone (G1 and G2 in Figure 8b). Another concentric cluster lies between linear along the Meadow Valley Fault and the Heathcote Fault Zone (G3 in Figure 8b). All these concentric depth clusters are located within the vicinity of the east–west oriented late Devonian granite known as the Cobaw Batholith, which is bounded by linear Euler depth clusters associated with major faults near the Heathcote Fault Zone (Clemens et al. 2016; VandenBerg et al. 2000).

Analysis of geologically constrained cluster assemblages

In Figure 5b and c, unclustered cross-sections of Euler depth solutions are shown along profiles A-A′ and B-B′. The top of the Euler deconvolution cross-section along profile A-A′ represents an undulating surface minimum depth of approximately 5 km. The lower part of the cross-section is less coherent at a depth of 15 km at about 200 km from the eastern edge of the traverse (Figure 5b). The unclustered transverse (Figures 5b and c) provides no structural or lithologic information as these key geological boundaries cannot be identified from the section.

To perform the sensitivity to the number of clusters and cluster boundaries due to step changes in the minimum number of points (minPts) and epsilon distance (ε), of 2500 ± 200 m (see Figure 5a). We categorized the cluster assemblages based on the proximity of their minPts to the inflection point on the K-NN graph (see Figure 5a), beginning with the end-members (Figure 9a) and proceeding to the assemblages nearest to the inflection point (Figure 9d).

End-member clusters across profile AA′

The end-member clusters that are the least optimised are those produced using the 28-NN and 36-NN (Figure 9a). The consistent cluster boundaries in these least optimised end-member clusters provide a means to determine cluster boundaries that will persist in all cluster assemblages. These consistent cluster boundaries are regarded as first-level cluster boundaries and are most likely to be associated with geological structures.

The end-member group have the highest variability in the number of clusters, as the 28-NN and 36-NN produced 6 and 12 clusters, respectively, yet share cluster boundaries with similar geometries (Figure 9a). Both cluster assemblages show cluster boundaries that runs from the western edge at a depth of 10 km and dips east down to 30 km at easting 675,000 m. Another cluster boundary in both cluster assemblages dips west from easting 860,000–790,000 m (Figure 9a). The boundaries in both end-member cluster assemblages are few and are first-level cluster boundaries. The first-level cluster boundaries are likely to persist in all cluster assemblages as we converge towards the inflexion point and so are considered the boundaries that are most likely to reflect faults or lithological boundaries.

Significant differences between the cluster assemblages are the number of clusters and the geometry of the resulting cluster boundaries. The 36-NN cluster assemblage has the highest number of clusters associated with the top, resulting in cluster boundaries with varying geometries. Within the 36-NN cluster assemblage, at least six major cluster boundaries extends to the top. Two clusters exist at easting of 830,000 m from depth of 5–15 km forming a west-dipping cluster boundary. A cluster boundary at easting of 670,000 m that dips west extends from depths of 10 km to depth of 20 km. In contrast, the cluster assemblage derived with 28-NN has fewer clusters and boundaries within the traverse. A cluster boundary exists within 28-NN cluster assemblage extending from the western fringe of profile A-A′ and terminating at the easting of 860,000 m (Figure 9a).

Other clusters across profile A-A′

Different cluster parameters impact properties of resulting cluster boundaries, allowing us to identify how geological constraints result in the most optimal parameter values that produce the most meaningful and useful clustering. We analyse other cluster assemblage pairs (Figure 9b–d) with corresponding parameters as they approach the inflexion point for clustering done with useful geological constraints. Variations in the number and geometry of other clusters reflect how these the geometry and number of cluster boundaries change their K-NN values approaches the inflexion point. Generally, the number of clusters tends to increase relative to the number of neighbourhood points at lower k-NN values before the inflexion point. The number of clusters becomes constant at higher numbers of neighbourhood points, usually after the inflexion point.

For clusters before the inflexion point, such as cluster assemblages, 29-NN and 30-NN produced the same number, but 31-NN has a higher number (12) of clusters. There is a large increase in the number of clusters produced after the inflexion point, as 33-NN and 34-NN give 13 clusters. Clusters continue to increase until a large drop at cluster assemblage 35-NN, which gives 11 clusters. This drastic fall in the number of clusters, after the inflexion point, suggests over-optimisation and the likely production of unreasonable clusters. The domain where the quantity of clusters stays the same or marginally increases as K-NN values increase past the inflexion point is deemed the domain with optimised clusters (e.g. cluster assemblages of 32-NN, 33-NN and 34-NN). It’s important to note the changes in the number of clusters directly affects the resulting geometry and number cluster boundaries formed within each cluster assemblage. Eliminating cluster assemblages based on their proximity to the inflexion point can be used to identify optimised clusters in absence of constraints.

The K-NN distance graph (Figure 5a) shows that the clustering pairs 30-NN and 32-NN, and 31-NN and 33-NN, are the closest to the inflexion point. This means that these pairs are most likely to produce optimised clusters. They demonstrate cluster boundaries with geometries and depths to those from previous clusters, including first-order clusters. The 31-NN and 33-NN pairs are the most optimised clusters and reflect the geological complexity across the traverse. While other clustering assemblages show only two major west-dipping clusters within eastings of 800,000 and 820,000 m, the 31-NN and 33-NN pairs have more varied cluster boundaries at different depths. Only these clustering pairs show similar west-dipping boundaries at 675,000 m.

Optimised parameters: clusters across B-B′

We apply clustering parameters of 33 points and an epsilon (ε) of 2500 m from the A-A traverse to Euler depths in cross-section B-B′ (Figure 10). Section B-B′ differs from section A-A′ because the orientations of the Euler depth solutions from the highly magnetic basaltic cover rocks rather than the basement geology. Nevertheless, 10 distinct cluster assemblages were generated at varying depths.

Figure 10. Euler depth solutions across profile B-B′ clustered according to optimised parameters derived from profile A-A′. Also, note the short path to achieve geologically relevant clusters along B-B′ (see Figure 3).

These cluster assemblages show generally similar geometries to those found in the optimised clusters for the A-A′ traverse. Notably, there is a similar east-dipping cluster boundary at the western edge starting at 5 km depth and extending to a depth of 30 km at easting 750,000 m (Figure 10). This boundary demarcates some minor west-dipping cluster boundaries with the least number of clustering points.

A major west-dipping boundary at the western edge extends from a shallow iso-surface to depths of ∼30 km along the traverse at easting 775,000 m. Below this boundary, depth solutions are absent at greater depths. This major west-dipping cluster boundary also bounds most of the narrow east-dipping boundaries in the eastern half of the traverse. Another steep west-dipping cluster boundary exists at easting 725,000 m with shallow iso-surface. This boundary terminates at the major east-dipping boundary at a depth of ∼15 km.

Analysis of unconstrained cluster

Using constrained parameters allows us to achieve optimal clustering that quickly outlines the geological complexity, as there are an infinite number of potential clusters and cluster boundaries that could be formed. We assess the effectiveness of our method when the clustering parameters are unconstrained but are instead identified by a notable increase across an average epsilon within a range of K-NN values. Here, the epsilon is generated directly from the average distance of each Euler depth value from the other. A significant jump is observed in the K-NN window ranging from minPts of 12–19.

In general, an increase in the number of clusters aligns with an increase in K-NN values, except for outliers. As the number of nearest-neighbour distance assemblages increases, the number of clusters also increases. The number of minPts increases relative to clusters and the resulting cluster boundaries. However, an outlier appears at a minPts of 13 and 19, which shows an unusually high 27 clusters (Figure 11b and h). The outlier is a cluster assemblage that does not follow the general trend or the expected number of clusters. An outlier assemblage can also be identified when the number of clusters before the inflexion point is higher than that after the inflexion point. However, the outlier cluster assemblage should show the general characteristics of clusters and cluster boundaries identified in cluster assemblages derived from an adjacent number of neighbourhood points on the K-NN graph.

Figure 11. DBSCAN clustering of Euler depth solutions derived from HRAM across profile C-C′, the Heathcote Fault Zone, for each cluster assemblage using ε = 800 m and K-NN values, from 12 to 19, as shown in Figure 6. Clustering shown in (a) and (g) are end-member cluster assemblages. Cluste assemblages (b) and (h) are regarded as outliers. (d) and (e) are regarded as the optimised cluster assemblages. All depth solutions are projected to GDA94/MGA 54.

Most cluster boundaries in the end-member assemblages produced complex clusters with varied depths and geometries (Figure 11a and g). Consistent cluster boundaries across all cluster assemblages, arising from clusters with similar geometries, locations, and depths. These are considered first-level cluster boundaries and reflect the geology at depth. A notable difference between the cluster boundaries of end-member cluster assemblages is that minPts = 12 (Figure 11a) shows a single cluster between easting 830,000 and 850,000 m while minPts = 19 shows multiple clusters (Figure 11 h). As a result, no cluster boundary can be identified in minPts = 12 cluster assemblage while minPts = 19 cluster assemblage shows several cluster boundaries from easting 830,000 and 850,000 m. In contrast, similar or near-identical clusters occur using minPts of 17 and 18 despite the difference in the generated number of clusters (Figure 11f and g).

Cluster assemblages closest to the inflexion point typically show a maximum positive difference in the number of clusters. The two cluster assemblages below and above the inflexion point are minPts = 15, which gives 25 clusters and minPts = 16 with 27 clusters (Figure 11d and e). Though these clusters show differences in their number of clusters, they are the most optimised as they are closest to the inflexion point. Generally, they share similarities in the geometry of major cluster boundary are considered cluster assemblages that best reflect the complexity along the traverse.

Implications of cluster assemblages

Optimum clusters across A-A′ and B-B′

We initially compared the optimised cluster assemblage on profile A-A′ with the near-parallel previously interpreted 2D reflection seismic survey 06GA-V1 to 06GA-V3 and then used the knowledge gained to examine other optimised Euler depth cross-sections (Figures 2, 9d & 10). For traverse A-A′, the optimised cluster assemblages were produced using a neighbourhood radius epsilon (ε) of 2500 ± 200 m and 33-NN (Figure 12).

Figure 12. Optimised clustered Euler depth solutions of the cross-section along profile A-A′. Cluster boundaries strongly correlate with the locations of the major zone-bounding faults at depth. The first-order cluster boundaries are in red and, other cluster boundaries, in black.

Figure 13. Crustal architecture beneath the Newer Volcanic Province derived from interpreted clustered Euler depth solutions along the B-B′ traverse. (a)Uninterpreted cluster assemblage (b) Interpreted cluster assemblages.

Comparing the optimised cluster assemblage with the pre-interpreted seismic line reveals that the cluster boundaries align with major boundaries along the seismic line associated with Proterozoic to Paleozoic faults (Cayley et al. 2011; Willman et al. 2010). Our interpretation of these clusters is also constrained by the depth solutions clustering along these faults in plan view (Figure 7). The geometries and locations of the major faults interpreted in previous studies correspond with most of the first-level cluster boundaries in the end-member and other cluster assemblages. Furthermore, cluster boundaries retain the same cross-cutting relationships and preserve the structures’ geological history.

The steeply east-dipping cluster boundary that extends from the western edge to a depth of 30 km at easting 675,000 m correlates with the Moyston Fault (see Figures 2 and 11). Other first-level cluster boundaries, such as the west-dipping cluster boundaries at 790,000 and 860,000 mE, correlate with the Whitelaw and Mount William Faults. One significant west-dipping boundary at 700,000 mE terminates at a depth of ∼22 km on the east-dipping cluster boundary associated with the Moyston Fault. This cluster boundary aligns with the Avoca Fault. Other cluster boundaries correlate with less significant faults, such as west-dipping cluster boundaries at easting of 695,000, 710,000, and 760,000 m, which align with imbricated west-dipping faults identified in the seismic section (Cayley et al. 2011; Willman et al. 2010).

Using the optimised cluster parameters from profile A-A′ for the cluster assemblage on B-B′, we observed a strong correlation of geometry, location, and overprinting relationships of cluster boundaries and geologic structures despite the masking effect of highly magnetic basalts of the Newer Volcanic Province (Figure 13). In the cluster assemblage along profile B-B′, the east-dipping boundary associated with the Moyston Fault truncates the west-dipping boundary associated with the Avoca Fault at a depth of ∼15 km, which is a shallower depth of than that observed in profile A-A′. The cluster boundary that coincides with the Mount William Fault continues to deeper levels in traverse B-B′ compared with profile A-A′. In profile B-B′ there are several minor east and west-dipping cluster boundaries that are not evident in profile A-A′, suggesting a subtly different structural pattern in the south. The distance between the Moyston and the Mount William Fault is greatly reduced in profile B-B′, suggesting they converge to the south (see Rawling et al. 2011).

Figure 14. Geologic interpretation of optimised clustering of Euler depth solutions generated from higher resolution airborne magnetic data along profile C-C′ (Figure 8). Depth clusters denoted as G1, G2 and G4 are located beneath concentric clusters in plan view (see Figure 8) associated with exposed granitic rocks.

Optimum cluster across the Heathcote Fault Zone

We apply our method across the Heathcote Fault Zone using unconstrained clustering parameters by allowing clustering parameters to be automatically determined based on a significant jump of K-NN values over an average epsilon distance (Figure 8). For optimised cluster assemblages, we discuss the geologic significance and how cluster assemblages provide insight into the geometry of faults at depth. We correlate cluster interpretations with the location of depth clusters in plan view.

The boundaries of the optimised cluster assemblage along profile C-C′ strongly correlate with the known structures in and around the Heathcote Fault Zone (Figures 8 and 14). A major west-dipping cluster boundary correlates with the Mount William Fault (Figure 14). This west-dipping boundary appears to truncate an east-dipping boundary west of Mount William Fault at a depth of 300 m. Along the Heathcote Fault (Figure 14), there is a west-dipping cluster boundary that terminates on an unnamed east-dipping boundary at 1 km. The unnamed east-dipping boundary continues deeper to 1.8 km, where it terminates on a less evident west-dipping boundary at easting 830,000 m that appears to continue deeper as marked by dashed lines.

NNE trending linear clusters in the map view along the Meadow Valley Fault (Figure 8), along the B-B′ traverse, correspond to a west-dipping cluster boundary extends to deeper depths (Figure 1 and 14). Less prominent west-dipping cluster boundaries at 830,000 mE suggests the presence of a west-dipping faults at depth of 2 km. At the western edge of the traverse is a west-dipping cluster boundary that lies along the Fosterville Fault (Figures 8 and 14).

Concentric clusters, evident in plan-view, align with known granitic rock exposures (Slater and Haydon 1999) and are bounded by major basement faults. East of the Mount William Fault, along profile C-C′, subcircular Euler depth cluster G1 (Figure 8) coincides with basaltic outcrops. Similar elliptical clusters (G2, G3 and G4) occur between the Heathcote Fault and the Fosterville Fault. The cluster boundaries associated with major faults extend deeper than the concentric clusters (Figure 14).

Discussion and conclusions

Effect of data resolution and method limitations

Data resolution is a critical factor in how combines unsupervised machine learning with Euler deconvolution solutions performs, both in determining the architecture of structures at depth and for other geophysical techniques. The method uses the capacity of the deconvolution process to distinguish subtle magnetic changes resulting from geological differences despite the homogeneous nature of the conventional Euler deconvolution equation used in this study. The number of clusters and resulting cluster boundaries are directly limited by the inherent trade-off that exist between data resolution and depth of investigation determined by the choice of the window size used during the deconvolution process.

For example, the global magnetic data grid with 2 arc-minute spacing and window size of 40 km can be used to image deep-crustal geologic structures, such as kilometre-wide fault damage zones and shear zones. In contrast, aeromagnetic data with a 40 m cell size gives depth results along individual faults within the Heathcote Fault Zone (Figures 11 & 14). Clustering of deconvolution depth results can determine the architecture of faults at metre scales and their dip at depths less than 100 m. However, this window size provides no structural information beyond 2 km due to Euler depth results’ limitations. The strong correlation of the DBSCAN clustering results with pre-interpreted 2D reflection seismic data (see Figure 12), supports the effectiveness our technique even in presences of unaccounted inhomogeneity in the magnetic field during the Euler deconvolution process. Nevertheless, DBSCAN clusters of depth solutions from the Euler deconvolution that considers the inhomogeneity in the magnetic field will yield more clusters that represent geological complexities (Fedi et al. 2015).

The method works best when DBSCAN’s parameters are constrained, particularly the minimum radius of the cluster (see Figures 12 and 13). Using previously known constraints can quickly generate optimised clusters representing geological complexities along the traverse. However, if the minimum radius is smaller than the data resolution, the clustering process yields clusters and cluster boundary artefacts that do not reflect geological structures. Also, while most cluster boundaries determined using optimum parameters correlate with known geologic structures, not all cluster boundaries represent structures, and over-interpretation is possible. Cluster boundaries should first be interpreted from coarser geological features, and only after that should finer-scale cluster boundaries from parameters nearer to the end-members on the K-NN graph be used.

Potential applications

Our new technique demonstrates extensive potential in a variety of applications. As a swift and reliable method, it helps resolve the dip and strike of both deep and shallow geologic structures. The method can help constrain well-known, more complex and expensive imaging techniques such as seismic refraction and reflection, potential field modelling, and geophysical inversion. It provides critical information on the first 500 m depth from high-resolution magnetic data aiding the interpretation of other near-surface techniques like electromagnetic and electrical tomographic methods (Boaga et al. 2020; Romero-Ruiz et al. 2018) and multichannel analysis of surface waves (Loo and Leong 2018; Pelton 2005) also deeper crustal depth imaging techniques such as seismic reflection (Peron-Pinvidic et al. 2022; Stephenson et al. 2021), refraction (Mooney 2015; Stixrude and Jeanloz 2015) and ambient noise (Bem et al. 2020; O’Donnell et al. 2023; Zeng et al. 2021).

Imaging of crustal structures under igneous or volcanic provinces has long been a challenge (Ernst et al. 2005). The strong attenuation properties of the overlying volcanic rocks make these structures difficult to identify in seismic sections (Grijalva et al. 2018; Holt et al. 2013; Song et al. 2018; Yule and Spandler 2022). The high magnetic susceptibilities of many volcanic rocks also mask the magnetic signatures of underlying structures and rocks (Blaikie et al. 2014; Lillis et al. 2009; Raimi et al. 2014). However, our technique can accurately determine the locations and dips of faults in a basaltic province in southeast Australia. This has implications for not only better understanding the basic structures in these areas but also in their relationship to the volcanic plumbing systems and their ascent pathways (Ernst et al. 2005; Guardo et al. 2022). These structures help in shaping continental margins that are often sites of mafic volcanic or thick sedimentary sequences (Anudu et al. 2014; Bladon et al. 2015; Gibson et al. 2013; Samsu et al. 2021) and potentially, their links with orebodies (Blundell et al. 2005; Boyce et al. 2014; Diakov et al. 2002; Graham et al. 2017; Ji et al. 2023; Pinotti et al. 2016).

The potential to extend our approach to resolve the architecture of geologic structures at depth in 3D space represents another important but challenging area that requires attention in the future. Having resolved the location and dip of geologic structures at various scales along profiles in 2D which can improve sub-surface models from the same data, the current approach should be able to resolve the architecture of structures in 3D space if the parameters of the DBSCAN algorithm are optimised. Clustering and/or segmentation of points in 3D or multi-space has largely lagged behind 2D machine learning algorithms but has huge promise in finding relationships in Euler depth solutions that are geologically meaningful (Ahmed and Chew 2020; Liu et al. 2023).

Acknowledgement

The authors thank Geoscience Australia and Geological Survey Victoria for providing and maintaining publicly available high-resolution magnetic data. This study was supported by the Australian Society of Exploration Geophysics Foundation Grant (RF21P01) and Monash-IITB scholarship awarded to Chukwu. We also thank the editor, Mark Lackie, the associate editor, Sam Matthews, and three additional anonymous reviewers for their constructive comments that have helped improved the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work was supported by Australian Society of Exploration Geophysics Foundation Grant [Grant Number RF21P01]; Monash-IITB scholarship awarded.