Enhancement of the balanced total horizontal derivative of gravity data using the power law approach

Abstract

Mapping the boundaries of subsurface geological structures constitutes a fundamental aspect of interpreting gravity data. Various methods relying on gravity derivatives have been introduced to delineate these boundaries. However, these approaches are not without their limitations, often leading to the generation of false edges or resulting in divergent outcomes. This article introduces a novel method that utilizes the power law approach and balanced total horizontal derivative for mapping gravity data edges. To assess the robustness of the proposed technique, we conducted tests on model examples and applied it to a real case study in northwest Vietnam. The results demonstrate that our method can accurately and clearly map edges in comparison to existing methodologies, while also preventing the generation of spurious edges in the output maps.

Keywords:

• Edge detection
• gravity data
• power transform
• balanced total horizontal derivative

1. Introduction

Gravity surveys play a crucial role in revealing subsurface structures (Prasad et al. 2018; Saibi et al. 2019; Ekinci et al. 2023). Gravity anomaly maps provide insights into the spatial distribution of buried density structures at varying depths, each exhibiting distinct geometric and physical characteristics, making the detection of corresponding geological bodies challenging (Saibi et al. 2012a, 2012b; Oksum et al. 2019, 2021; Aprina et al. 2024). Various enhancement methods are commonly employed to highlight specific features, such as density source boundaries (Ekinci and Yiğitbaş 2015; Yuan et al. 2016; Narayan et al. 2017, 2021; Kafadar 2022; Sahoo et al. 2022a, 2022b; Altınoğlu 2023). These techniques typically rely on vertical or horizontal derivatives or their combinations (Ekinci 2017; Nasuti et al. 2019; Eldosouky et al. 2022a, 2022b). Identifying source boundaries often involves utilizing maximum, minimum, or zero values in filtered anomalies (Li et al. 2023; Kamto et al. 2023a, 2023b).

Historically, methods like the vertical gradient (Evjen 1936), total horizontal derivative (Cordell and Grauch 1985), and analytic signal (Roest et al. 1992) have been employed for this purpose. Beiki (2010) introduced an edge detector based on directional analytic signals to mitigate interference effects from neighboring structures. However, these methods face challenges in accurately determining the edges of deep bodies (Ghomsi et al. 2022; Eldosouky et al. 2022a).

To address the delineation of shallow and deep bodies, several normalization methods have been proposed (Pham et al. 2020, 2022; Prasad et al. 2022a, 2022b). The tilt angle normalization technique was first introduced by Miller and Singh (1994), utilizing the total horizontal derivative to normalize the vertical derivative. Subsequent methods, such as those presented by Verduzco et al. (2004), Wijns et al. (2005), and Oruç (2011) introduced variations and improvements to this approach. Nasuti and Nasuti (2018) proposed an improved version of the tilt angle, demonstrating its effectiveness in enhancing sharp signals over the source edges. Zareie and Moghadam (2019) introduced the total directional theta method based on potential field gradient tensor data, which is more effective than the traditional methods in detecting the edges of complex geological structures. Al-Bahadily et al. (2023) introduced the inverse tilt angle of second order gradients, and applied it to highlight reactivated regions and faults of the Iraq Southern Desert. Some authors also introduced other formats of normalization methods (Li 2013; Ma 2013; Yuan and Yu 2014; Ma et al. 2016; Chen et al. 2017; Alvandi et al. 2023; Ghiasi et al. 2023; Pham 2023; Pham et al. 2023).

Despite the benefits of normalization methods in balancing gravity anomalies from sources at various depths, some tend to produce false edges, leading to inaccurate geological interpretations (Pham and Prasad, 2023 2024). Ma et al. (2015) addressed this by introducing the balanced total horizontal derivative (BTHD), capable of outlining large and small amplitude edges simultaneously. In fact, this technique has already been used for interpreting magnetic anomalies by Ferreira et al. (2013). In recent years, many studies have applied the BTHD to identify buried geological structures from potential field anomalies. Weihermann et al. (2018) applied the BTHD method to interpret aeromagnetic data of the Paranaguá area (Brazil). Sahoo and Pal (2019) used the BTHD to identify fault zones in the Indian Ocean from gravity data. Zanella et al. (2020) applied the BTHD to map structures of the Campo Alegre area (Brazil) from aeromagnetic data. Reis et al. (2020) used this method to identify structures in Western Sao Francisco (Brazil) from aeromagnetic anomalies. Safani et al. (2023) determined the structural configuration of the Sengkang Basin of Indonesia using the BTHD of gravity data. While the BTHD avoids false information, it has been reported to produce divergent edges (Nasuti and Nasuti 2018; Alvandi and Ardestani 2023). To increase the sharpness of the edges, some techniques based on third-order have been developed (Zhang et al. 2015; Nasuti et al. 2019; Ibraheem et al. 2023), although the use of the high-order derivatives is more sensitive to high-frequency noise (Oliveira and Pham 2022).

In this article, we propose an enhanced version of the BTHD method, named PBTHD, which employs the power law approach. This proposed method aims to map source boundaries with increased resolution and without false edges. Additionally, to provide more stable results for edge detection, we employ the vertical derivative calculated from the finite-difference approach of Oliveira and Pham (2022) instead of those from the standard wavenumber domain approach. We illustrate the application of the PBTHD method through model examples and a real case study from northwest Vietnam.

2. New method

The total horizontal derivative (THD) method is widely utilized for delineating the edges of gravity data, and is defined by the equation (Cordell and Grauch 1985): (1) $$\mathit{\text{THD}}=\sqrt{{\left(\frac{\partial F}{\partial x}\right)}^2+{\left(\frac{\partial F}{\partial y}\right)}^2},$$(1) where F represents the gravity anomaly.

Recognizing the THD’s limitations in detecting edges of deep density sources, Ma et al. (2015) introduced a balanced version, BTHD, capable of outlining both shallow and deep source edges simultaneously. The BTHD utilizes the arctangent of the ratio between derivatives of the THD to map the edges: (2) $$\mathit{\text{BTHD}}=\mathit{\text{atan}}\frac{\frac{\partial \mathit{\text{THD}}}{\partial z}}{\sqrt{{\left(\frac{\partial \mathit{\text{THD}}}{\partial x}\right)}^2+{\left(\frac{\partial \mathit{\text{THD}}}{\partial y}\right)}^2}}.$$(2)

The BTHD yields more reliable source edges than traditional methods, making it recommended for geological body extraction from gravity and magnetic data (Pham et al. 2021a, 2021b).

To enhance signal resolution in the edge map, we employ the power law approach derived from digital image processing (Wahab et al. 2018). Specifically, we apply the power transform to $\left(\mathrm{\text{BTHD}}+\frac{\pi }{2}\right)$ to obtain sharper signals over source edges: (3) $$\mathit{\text{PBTHD}}={\left(\mathit{\text{atan}}\frac{\frac{\partial \mathit{\text{THD}}}{\partial z}}{\sqrt{{\left(\frac{\partial \mathit{\text{THD}}}{\partial x}\right)}^2+{\left(\frac{\mathit{\text{THD}}}{\partial y}\right)}^2}}+\frac{\pi }{2}\right)}^n,$$(3) where n is a positive number defined by the interpreter.

The horizontal gradients in Equation (3)(3) $$\mathit{\text{PBTHD}}={\left(\mathit{\text{atan}}\frac{\frac{\partial \mathit{\text{THD}}}{\partial z}}{\sqrt{{\left(\frac{\partial \mathit{\text{THD}}}{\partial x}\right)}^2+{\left(\frac{\mathit{\text{THD}}}{\partial y}\right)}^2}}+\frac{\pi }{2}\right)}^n,$$(3) are computed in the spatial domain using finite differences (Blakely 1995; Ekinci et al. 2020), while the vertical derivative is computed using the β-VDR method (Oliveira and Pham 2022) to attenuate noise: (4) $$\frac{\partial \mathit{\text{THD}}}{\partial z}=\frac{p_1\mathit{\text{THD}}\left(h_1\right)+p_2\mathit{\text{THD}}\left(h_2\right)+p_3\mathit{\text{THD}}\left(h_3\right)+p_4\mathit{\text{THD}}\left(h_4\right)+p_5\mathit{\text{THD}}\left(h_5\right)}{\Delta h},$$(4) where p₁, …, p₅ are given by: (5) $$\{\begin{array}{l}{p}_1=\left(2{\beta }^3+15{\beta }^2+35\beta +25\right)/12,\\ p_2=\left(-8{\beta }^3-54{\beta }^2-104\beta -48\right)/12,\\ \begin{array}{c}{p}_3=\left(12{\beta }^3+72{\beta }^2+114\beta +36\right)/12,\\ \begin{array}{l}{p}_4=\left(-8{\beta }^3-42{\beta }^2-56\beta -16\right)/12,\\ p_5=\left(2{\beta }^3+9{\beta }^2+11\beta +3\right)/12,\end{array}\end{array}\end{array}$$(5) and $\mathit{\text{THD}}(h_i)$ represents the upward-continued data to $h_i=z_0-\beta \Delta h-\left(i-1\right)\Delta h$ with $z_0$ as the observed plane height, $\beta$ = 30 and $\Delta h=\frac{1}{10}$ of grid spacing (Oliveira and Pham 2022).

The primary objective of the proposed method is to identify locations of sudden density variations. The maximum values in the PBTHD map are employed to delineate source boundaries similar to the BTHD filter. The main advantage of the PBTHD lies in its ability to produce boundaries with a high sharpness. By utilizing the ratio of THD gradients, the PBTHD offers a clearer mapping of deep body boundaries compared to the THD.

3. Edge detectors of analysis

To assess the efficacy of the proposed method, we conducted a comparative analysis with several standard and recent techniques. The standard techniques include the total horizontal derivative (THD), analytic signal (AS), tilt angle (TA), total horizontal derivative of the tilt angle (THDTA), and theta map (TM). Additionally, we compared our method with recent techniques, such as the balanced total horizontal derivative (BTHD), improved normalized horizontal derivative (ITDX), and total horizontal derivative of NTilt (THDNTilt).

The analytic signal (AS) is designed to accentuate source boundaries by leveraging peaks, and is defined as follows (Roest et al. 1992): (6) $$\mathit{\text{AS}}=\sqrt{{\left(\frac{\partial F}{\partial x}\right)}^2+{\left(\frac{\partial F}{\partial y}\right)}^2+{\left(\frac{\partial F}{\partial z}\right)}^2}.$$(6)

The tilt angle (TA) represents the first normalized technique employing zero values to extract source boundaries (Miller and Singh 1994): (7) $$\mathit{\text{TA}}=\mathit{\text{atan}}\frac{\frac{\partial F}{\partial z}}{\sqrt{{\left(\frac{\partial F}{\partial x}\right)}^2+{\left(\frac{\partial F}{\partial y}\right)}^2}}.$$(7)

The total horizontal derivative of the tilt angle (THDTA) enhances the display of bodies by providing maxima over the boundaries (Verduzco et al. 2004): (8) $$\mathit{\text{THDTA}}=\sqrt{{\left(\frac{\partial \mathit{\text{TDR}}}{\partial x}\right)}^2+{\left(\frac{\partial \mathit{\text{TDR}}}{\partial y}\right)}^2}.$$(8)

The theta map (TM) normalizes the THD by the AS, yielding minima over the edges (Wijns et al. 2005): (10) $$\mathit{\text{TM}}=\mathit{\text{acos}}\frac{\sqrt{{\left(\frac{\partial F}{\partial x}\right)}^2+{\left(\frac{\partial F}{\partial y}\right)}^2}}{\sqrt{{\left(\frac{\partial F}{\partial x}\right)}^2+{\left(\frac{\partial F}{\partial y}\right)}^2+{\left(\frac{\partial F}{\partial z}\right)}^2}}.$$(10)

The improved normalized horizontal derivative (ITDX) relies on the second derivatives of gravity data to detect edges (Ma et al. 2016): (11) $$\mathit{\text{ITDX}}=\mathit{\text{atan}}\frac{\sqrt{{\left(\frac{{\partial }^{2}F}{\partial z\partial x}\right)}^2+{\left(\frac{{\partial }^{2}F}{\partial z\partial y}\right)}^2}}{\left|\frac{{\partial }^{2}F}{{\partial \mathrm{z}}^2}\right|}.$$(11)

The second vertical derivative in Equation (11)(11) $$\mathit{\text{ITDX}}=\mathit{\text{atan}}\frac{\sqrt{{\left(\frac{{\partial }^{2}F}{\partial z\partial x}\right)}^2+{\left(\frac{{\partial }^{2}F}{\partial z\partial y}\right)}^2}}{\left|\frac{{\partial }^{2}F}{{\partial \mathrm{z}}^2}\right|}.$$(11) is computed using the Laplace equation to mitigate noise effects: (12) $$\frac{{\partial }^{2}F}{{\partial \mathrm{z}}^2}=-\frac{{\partial }^{2}F}{{\partial \mathrm{x}}^2}-\frac{{\partial }^{2}F}{{\partial \mathrm{y}}^2}.$$(12)

The total horizontal derivative of the NTilt (THDNTilt) normalizes the second vertical derivative using the enhanced analytic signal (Nasuti and Nasuti 2018): (13) $$\mathit{\text{THDNTilt}}=\sqrt{{\left(\frac{\partial \mathit{\text{NTilt}}}{\partial x}\right)}^2+{\left(\frac{\partial \mathit{\text{NTilt}}}{\partial y}\right)}^2},$$(13) where NTilt is defined as (Nasuti and Nasuti 2018): (14) $$\mathit{\text{NTilt}}=\mathit{\text{atan}}\left(k^2\frac{\frac{{\partial }^{2}F}{{\partial z}^2}}{\sqrt{{\left(\frac{\partial {\mathit{\text{AS}}}_2}{\partial x}\right)}^2+{\left(\frac{\partial {\mathit{\text{AS}}}_2}{\partial y}\right)}^2}}\right),$$(14)

In Equation (14)(14) $$\mathit{\text{NTilt}}=\mathit{\text{atan}}\left(k^2\frac{\frac{{\partial }^{2}F}{{\partial z}^2}}{\sqrt{{\left(\frac{\partial {\mathit{\text{AS}}}_2}{\partial x}\right)}^2+{\left(\frac{\partial {\mathit{\text{AS}}}_2}{\partial y}\right)}^2}}\right),$$(14) , AS_n and k are given by: (15) $${\mathit{\text{AS}}}_2=\sqrt{{\left(\frac{\partial F_{\mathit{\text{zz}}}}{\partial x}\right)}^2+{\left(\frac{\partial F_{\mathit{\text{zz}}}}{\partial y}\right)}^2+{\left(\frac{\partial F_{\mathit{\text{zz}}}}{\partial z}\right)}^2},$$(15) (16) $$k=\frac{P}{\sqrt{{\Delta x}^2+{\Delta y}^2}}$$(16) with $F_{\mathit{\text{zz}}}=\frac{{\partial }^{2}F}{{\partial z}^2},$ P is regional gravity value, $\Delta x$ and $\Delta y$ are the grid spacings along the x and y directions, respectively.

4. Synthetic data experiments

The presented filter’s robustness is demonstrated using gravity datasets with and without noise. A comparison with other filters is also conducted. The synthetic model comprises two overlapped prismatic bodies (G1 and G2) and three prismatic bodies (G3, G4, and G5), each of the same size (Figure 1a). Geometrical and density information for the sources is detailed in Table 1. The anomaly of the model is computed on a 201 × 201 grid with 1 km spacing. Four cases are considered: (1) all sources have positive density contrasts, (2) the model contains both positive and negative density contrasts, (3) 10% Gaussian noise is added to gravity data in the second case, (4) noisy gravity data is upward continued to a height of 1 km.

Figure 1. (a) Perspective view of the gravity model, (b) gravity anomaly generated by sources with positive density contrasts, (c) gravity anomaly generated by sources with both positive and negative density contrasts. The dashed lines represent the actual borders of the sources. The white line denotes the profile AB over the model (see Figure 3).

Table 1. Parameters of the sources.

Parameters	G1	G2	G3	G4	G5
Center coordinates (km; km)	50; 100	50; 100	144; 155	144; 100	144; 45
Width (km)	10	60	30	30	30
Length (km)	60	90	70	70	70
Top depth (km)	2	3	3	6	9
Bottom depth (km)	3	6	6	9	12
Density contrast (g/cm^3)	0.1	0.2	±0.3	0.3	±0.3

In the first case, the synthetic gravity anomaly generated by the bodies is presented in Figure 1b. Figure 2a and b show the edges determined by the THD and AS methods, respectively. Both methods exhibit dominance by strong amplitude responses from the shallow body G3, while weak amplitude responses from the overlapped body G1 and other bodies are blurred. The TA in Figure 2c balances different anomaly amplitudes but lacks sharp responses at the boundaries. The THDTA in Figure 2d successfully determines the edges of the bodies G2, G3, and G4 but is less effective for the deepest source G5, introducing false edges around source G1. The edges obtained from the TM filter (Figure 2e) indicate the source edges’ locations, but it fails to show the edges of the body G1 and presents the detected edges of the deepest source G5 as larger than their actual size. Applying the BTHD filter (Figure 2f) detects all edges without false edges but results in diffuse edges. Figure 2g and h show the edges obtained from applying the ITDX and THDNTilt methods, respectively. While these methods provide higher resolution edge maps compared to the THD, AS, TA, THDTA, TM, and BTHD, they introduce false boundaries around source G1. For the presented method (PBTHD), Figure 2i–o show the edges with varying parameter values (n = 2, 3, 4, 6, 8, 12, and 20). Increasing n enhances the resolution of edges in the PBTHD maps, but the edges of bodies G1 and G5 become faint.

Figure 2. Enhanced maps of data in Figure 1b (a) THD, (b) as, (c) TA, (d) THDTA, (e) TM, (f) PBTHD, (g) ITDX, (h) THGNTilt, (i) PBTHD with n = 2, (j) PBTHD with n = 3, (k) PBTHD with n = 4, (l) PBTHD with n = 6, (m) PBTHD with n = 8, (n) PBTHD with n = 12, (o) PBTHD with n = 20. The dashed lines on the maps indicate the actual borders of the sources.

To visually assess the boundary determination results with varying n values, we analyze the horizontal profile AB (Figure 1b). Figure 3a and b display the model and gravity anomaly along this profile. Results of the BTHD and PBTHD with n = 2, 3, 4, 6, 8, 12, and 20 are shown in Figure 3c–j. The PBTHD profiles exhibit sharper signals over the edges than BTHD, although the transformed anomaly amplitude over body G1 decreases with increasing n. The PBTHD yields weak amplitude responses over body G1 for n > 4; hence, we use n = 3 for other examples.

Figure 3. (a) The sources along the profile AB in Figure 1b (b) Gravity data along the profile, (c) PBTHD, (d) PBTHD with n = 2, (e) PBTHD with n = 3, (f) PBTHD with n = 4, (g) PBTHD with n = 6, (h) PBTHD with n = 8, (i) PBTHD with n = 12, (j) PBTHD with n = 20.

Examining the model with positive and negative density contrasts, the second scenario involves density contrasts of −0.3 and −0.3 g/cm³ for the sources G3 and G5, respectively. The synthetic gravity anomaly of this model is displayed in Figure 1c. Results from the THD and AS methods are shown in Figure 4a and b, respectively. Both methods detect the edges of the shallow source G3 but fail to yield clear edges for the overlapped body G1 and other sources. Figure 4d shows the TA of gravity data, equalizing different amplitudes but not clearly distinguishing edges. The THDTA in Figure 4d clearly highlights the edges of the sources G2 and G3, with faint edges for the sources G4 and G5, introducing false information around the source G1. Figure 4e and f present edges obtained from the TM and BTHD methods, respectively. While the BTHD detects all edges, the TM fails to show the edges of the source G1 and presents false boundaries around the sources G3 and G5. Both the TM and BTHD bring edges with low resolution. Figure 4g and h depict edges obtained from the ITDX and THDNTilt methods, respectively. The THDNTilt provides results with very high resolution, but both the ITDX and THDNTilt introduce false boundaries around the sources G1, G3, and G5. Figure 4i depicts edges determined from the PBTHD method, demonstrating less dependence on structure depth and more convergent edges than the BTHD. The PBTHD provides high-resolution edges without generating spurious edges.

Figure 4. Enhanced maps of data in Figure 1c (a) THD, (b) as, (c) TA, (d) THDTA, (e) TM, (f) BTHD, (g) ITDX, (h) THGNTilt, (i) PBTHD with n = 3. The dashed lines indicate the actual borders of the sources.

To test the stability of PBTHD with respect to data contamination, 10% Gaussian noise is added to gravity data in Figure 1c. Noise-contaminated gravity data is shown in Figure 5a. Figure 6a and b present results from the THD and AS methods, which remain less effective in mapping edges. Figure 6c and d show results from the TA and THDTA methods. The THDTA is more sensitive than the TA, producing faint edges. Figure 6e and f present edges obtained from the TM and BTHD methods, respectively. While the BTHD provides clearer results than the TA, THDTA, and TM, the edges are diffuse. Figure 6g–i show the edges obtained from the ITDX, THDNTilt, and PBTHD methods. The ITDX and THDNTilt are more sensitive to noise, with the PBTHD demonstrating less sensitivity to noise than normalized techniques. The PBTHD remains effective in mapping all edges without introducing false information.

Figure 5. (a) Synthetic gravity data in Figure 4c corrupted with random noise, (b) Noisy synthetic gravity data after upward continuation of 1 km. The dashed lines indicate the actual borders of the sources.

Figure 6. Enhanced maps of data in Figure 5a (a) THD, (b) as, (c) TA, (d) THDTA, (e) TM, (f) BTHD, (g) ITDX, (h) THGNTilt, (i) PBTHD with n = 3. The dashed lines indicate the actual borders of the sources.

In the noisy case, the edges obtained from the PBTHD are accurate, but upward continuation is advised before applying this method to noise-contaminated data, as recommended by other authors (e.g. Ekinci et al. 2013; Nasuti et al. 2019; Zareie and Moghadam 2019; Alvandi and Ardestani 2023). Figure 5b shows noise-contaminated data after 1 km upward continuation. Outputs obtained from the THD, AS, TA, THDTA, TM, BTHD, ITDX, THDNTilt, and PBTHD are shown in Figure 7a–i, respectively. The THD and AS methods remain less effective in detecting the edges of deep sources (G1, G2, G4, and G5). In this case, the ITDX and THDNTilt are still noisy. While the TM and BTHD provide results with low resolution or additional edges, the PBTHD demonstrates better performance, providing high-resolution edges without introducing false boundaries.

Figure 7. Enhanced maps of data in Figure 5b (a) THD, (b) as, (c) TA, (d) THDTA, (e) TM, (f) BTHD, (g) ITDX, (h) THGNTilt, (i) PBTHD with n = 3. The dashed lines indicate the actual borders of the sources.

5. Application to real data

We applied the proposed method to gravity datasets from northwest Vietnam (Figure 8). The geology of northwest Vietnam includes two major units: the Indochina and Sino-Vietnam composite terranes (Tran et al. 2022). The studied area, documented as one of the most seismically active regions in Indochina (Zuchiewicz et al. 2004), experiences its most robust seismic activity along northwest-southeast trending faults (Zuchiewicz et al. 2004). Many big earthquakes occurred in this area (Figure 8) (Duong et al. 2013; Pham et al. 2023). The geological map of the study region and its surroundings is illustrated in Figure 8 (Koszowska et al. 2007). The area encompasses five northwest-southeast-oriented main structures: Late Palaeozoic-Mesozoic rift, Tay Bac fold system, Volcano-tectonic rift trough, Mesozoic superimposed rift, Trung Son fold system (Figure 8). The gravity data employed in this study is based on free-air anomalies (version 28.1) from the satellite gravity model of Sandwell et al. (2014), which can be downloaded from https://topex.ucsd.edu/cgi-bin/get_data.cgi. Bouguer gravity anomaly data were computed using the Parker method with a reference crustal density of 2670 kg/m³ (Pham 2020), as depicted in Figure 9. This dataset ranges from −217 to 34 mGal, with a mean value of 101 mGal, and exhibits a dominant northwest-southeast anomaly trend.

Figure 8. Geological map of the study area and its surroundings.

Figure 9. Bouguer gravity anomaly map of the study area. The black lines show surface geological boundaries.

In the actual implementation, we applied the THD, AS, TA, THDTA, TM, BTHD, ITDX, THDNTilt, and PBTHD to the Bouguer data of the study area. The results obtained from these applications are presented in Figure 10. The THD, illustrated in Figure 10a, is characterized by high amplitude responses primarily in the northern part of the study area. Although the THD map indicates northwest-southeast trending structures, many of these structures appear diffuse and faint. The AS results, shown in Figure 10b, mirror the THD with dominant high amplitude responses in the northern part of the area. However, the AS lacks sharp signals and often presents single bell-shaped anomalies over the bodies, as evident in Figure 10b. The TA technique, depicted in Figure 10c, effectively balances anomalies with different amplitudes but does not distinctly delineate structural boundaries. The THDTA result (Figure 10e) exhibits less effectiveness in providing a clear structural map for the study area, dominated by high amplitude isolated anomalies. The edges computed from applying the TM (Figure 10e) show structural boundaries in the southern part more clearly than the THD, AS, and THDTA. However, some adjacent boundaries are connected, complicating the structural interpretation. The BTHD, illustrated in Figure 10f, provides a balanced edge map without additional edges, facilitating geological interpretations compared to other methods, albeit with lower resolution. Figure 10g and h show the edges obtained from applying the ITDX and THDNTilt, respectively. Both methods present boundaries with high resolution, especially the THDNTilt. However, some boundaries in these maps are connected, potentially leading to incorrect structural interpretations. Additionally, both methods introduce additional edges, complicating real edge estimation. Figure 10i displays the edges obtained from applying the presented PBTHD method to gravity data in Figure 9. The PBTHD method effectively identifies various crucial geological boundaries, showing sharp signals for source boundaries and providing a balanced structural image for the study area. Our model examples demonstrate that this method outlines all edges, including additional boundaries.

Figure 10. Enhanced maps of Bouguer anomaly data in Figure 9, (a) THD, (b) as, (c) TA, (d) THDTA, (e) TM, (f) BTHD, (g) ITDX, (h) THGNTilt, (i) PBTHD with n = 3.

Given that the PBTHD is more effective than other methods in mapping density structure edges, we extracted the structural boundaries (black lines) from the PBTHD map. We compared these boundaries with faults (red lines) reported by Zuchiewicz et al. (2004) (Figure 11a). It is evident from Figure 11a that PBTHD aligns with northwest-southeast trending structures, consistent with the main structures’ trend in the area. The PBTHD edges coincide with some faults with northwest-southeast trends in the northeastern, southwestern, and central regions. Rose diagrams of detected boundaries (Figure 11b) and known faults (Figure 11c) reveal a major trend in the northwest-southeast direction, indicating a significant correlation between gravity boundaries and known faults. Notably, the presented method has uncovered some new structures not identified on the geological map.

Figure 11. (a) Structural boundaries (black lines) extracted by the PBTHD and the faults (red lines) reported by Zuchiewicz et al. (2004), (b) Rose diagram of the structural boundaries, (c) Rose diagram of the faults.

6. Conclusions

A novel technique based on the power law approach and balanced total horizontal derivative has been introduced to delineate the edges of gravity sources. The effectiveness of this approach was compared with several existing techniques, including the THD, AS, TA, THDTA, TM, BTHD, ITDX, and THDNTilt, using both noise-free and noisy synthetic gravity data, as well as real gravity data from northwest Vietnam. In synthetic studies, our proposed method successfully outlines all edges with varying amplitude anomalies. Importantly, it avoids generating any additional edges when the model incorporates both negative and positive density contrasts. Furthermore, our presented method demonstrates superior precision and clarity in recognizing edges, achieving high-resolution results. The results from noisy synthetic gravity data show that the presented method is less sensitive to noise than other second and third order derivative-based methods. Examining Bouguer gravity anomaly data as an illustration, our findings yield a more interpretable image compared to other methods. The results obtained from the presented method indicate that the most prominent structural lineaments in northwest Vietnam are in the northwest-southeast direction. These results align with known faults in the study area, providing further validation of the method’s efficacy.

Acknowledgement

Deep thanks and gratitude to the Researchers Supporting Project number (RSP2024R351), King Saud University, Riyadh, Saudi Arabia, for funding this research article. This research was funded by the research project QG.23.64 of Vietnam National University, Hanoi.

Data availability statement

Data are available on request from the authors.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Deep thanks and gratitude to the Researchers Supporting Project number (RSP2024R351), King Saud University, Riyadh, Saudi Arabia, for funding this research article. This research was funded by the research project QG.23.64 of Vietnam National University, Hanoi.