Euler inversion: locating sources of potential-field data through inversion of Euler’s homogeneity equation

SUMMARY

Locating the sources of observed disturbances in potential-field data is a challenging problem due to the non-unique nature of the inverse problem. The Euler deconvolution method was created to solve this issue, particularly for idealized sources (such as spheres and planar vertical dykes). Euler deconvolution has become widely used in potential-field methods due, in large part, to its low-computational cost and ease of implementation into software. However, it is widely known that Euler deconvolution suffers from some shortcomings: (1) non-uniqueness of the solution with respect to the depth of the source and the structural index (a parameter that represents the idealized shape of the source); (2) sensitivity to short-wavelength noise in the data derivatives which are used as inputs for the method. Here, we present a new method called Euler inversion which is a reformulation of the inverse problem of Euler’s homogeneity equation as an implicit mathematical model rather than a parametric one. Euler inversion is a constrained, nonlinear inverse problem capable of estimating both the model parameters (location of the source and constant base level) and the predicted data (potential field and its derivatives). We show that Euler inversion is less sensitive than Euler deconvolution to short-wavelength noise and to the presence of interfering sources in the data window. By also estimating the predicted data, Euler inversion is also able to estimate the best integer structural index to be used for inversion. Our results show that the estimated structural index minimizes the data misfit and coincides with those of the simulated sources. Furthermore, most matrices involved in the method are either sparse or diagonal, making Euler inversion computationally efficient. Tests on synthetic data and a real aeromagnetic data set from Rio de Janeiro, Brazil, demonstrate the effectiveness of Euler inversion to delineate sources with variable geometries and correctly estimate their depths.

1 INTRODUCTION

Estimating the depths of the sources of measured anomalies is a common challenge in potential-field geophysics. One of the most widely used techniques for providing depth estimates is Euler deconvolution (Thompson 1982; Reid et al. 1990). Its widespread adoption is due, in large part, to its low algorithmic complexity and fast computation times, both of which are orders of magnitude smaller than solutions from 3-D inverse problems. As a result, Euler deconvolution is widely available in both commercial and open-source software (Uieda et al. 2013, 2014). Unfortunately, this popularity has also led to abuses of the method, as reported in Reid & Thurston (2014) and Reid et al. (2014).

Euler deconvolution is a method that assumes potential-field data are generated by idealized sources, such as dikes, dipoles or pipes. The geometry of these sources is characterized by the structural index, a parameter that must be an integer to retain physical significance (Stavrev & Reid 2007; Reid & Thurston 2014). The technique involves performing a least-squares inversion of Euler’s homogeneity equation multiple times, in a moving window scheme. Each inversion estimates the base level, a constant shift in the data, and also the coordinates of a single idealized source potentially present within the study area.

It is well known that Euler deconvolution suffers from some limitations, of which we highlight:

• Separation of reliable and spurious solutions: The moving window scheme adopted in Euler deconvolution generates many estimated positions which are considered spurious and must be removed. Most of the spurious solutions happen when the moving window either lacks significant potential-field anomalies or only contains a truncated anomaly. FitzGerald et al. (2004) and Melo & Barbosa (2020) provide overviews of the many existing methods that have been developed to remove spurious solutions.

• Sensitivity to high-frequency noise: Random noise in the data is usually of high- frequency, which gets amplified in the derivative calculations. Since the field derivatives are used in the Jacobian matrix of the least-squares inversions, errors in the derivatives will have a large impact on the solution. Pašteka et al. (2009), Saleh & Pašteka (2012) and Florio et al. (2014) recommend using regularized derivatives or other smoothing techniques to reduce the noise amplification and obtain more reliable solutions. This is also why Euler deconvolution variants that rely on higher-order derivatives, like tilt-Euler deconvolution (Salem et al. 2007; Huang et al. 2019) and the combined analytic signal and Euler method (AN-EUL) (Salem & Ravat 2003), present a larger dispersion of estimated positions and are more sensitive to noise in general. Methods like finite-difference Euler deconvolution (Gerovska et al. 2005) and ratio-Euler deconvolution (Huang et al. 2022) were specifically developed to avoid the use of higher-order derivatives because of this noise-sensitivity issue.

• Correlation of the estimated depth and the structural index: Silva et al. (2001) demonstrated that the estimated depth from Euler deconvolution is directly correlated with the structural index used. The higher the structural index, the larger the estimated depth. This makes it very important to know the best integer structural index for the type of source being interpreted. Some Euler deconvolution variants have been developed that are able to estimate the structural index (e.g. Salem & Ravat 2003; Silva & Barbosa 2003; Gerovska et al. 2005; Salem et al. 2007; Melo et al. 2013; Florio & Fedi 2013; Florio et al. 2014; Melo & Barbosa 2018). However, most of them estimate real-valued structural indices instead of integers, are sensitive to noise, and tend to underestimate the structural index under realistic noise and signal overlap scenarios. Another approach is that of Mushayandebvu et al. (2004), who exploit the ill-conditioning of the Jacobian matrix of Euler deconvolution to detect the presence of 2-D sources (structural index of one for magnetic data and 0 for gravity data) in a data window and correctly estimate their position and strike.

Euler deconvolution and its variants are also known to struggle with models that have two or more contact points, like steps which have a top and a bottom. To solve this problem, methods like MaGSoundFDST method of Gerovska et al. (2010) were developed based on the similarity transform. MaGSoundFDST, in particular, is able to estimate structural index, source locations and the number of sources, hence side-stepping the problem of spurious solutions altogether.

We aim to tackle some of these issues by reformulating the inverse problem of solving Euler’s homogeneity equation. The issue of noise sensitivity can be traced back to the presence of data derivatives in the Jacobian matrix, which generally contain larger amounts of noise than the original potential field. We propose formulating the inverse problem as a non-linear optimization with Euler’s equation as a constraint. This is similar to ‘total least-squares’ in statistics (Van Huffel & Vandewalle 1991) and ‘combined adjustment’ in geodesy (Vaníček & Krakiwsky 1986). Another advantage of this new formulation is the ability to calculate predicted data for the potential field and its three derivatives, which is impossible in Euler deconvolution and all of its variants. We call our new method ‘Euler inversion’.

2 METHODOLOGY

Starting with Thompson (1982) and Reid et al. (1990), Euler’s equation has been used to estimate the source positions of gravity and magnetic data. In this section, we will review the solution of Euler’s equation for the source location $(x_o,y_o,z_o)$ by Euler deconvolution (Reid et al. 1990) and then present a new method, called Euler inversion, for solving Euler’s equation using total least-squares.

We start with Euler’s homogeneity equation

in which f is a homogeneous function (in this case, a potential-field), ${\partial }_{\alpha }$ is the derivative operator in the $\alpha$ direction, $(x,y,z)$ are the coordinates of the observation point, $(x_o,y_o,z_o)$ are the coordinates of the field source, b is the base level representing a constant shift in the field and $\eta$ is the structural index, which is related to the nature of the source and how its potential-field values decay with distance (Ruddock et al. 1966; Reid & Thurston 2014). The coordinate system is defined with x pointing eastwards, y pointing northwards and z pointing upwards. Eq. (1) relates the coordinates of the source with the potential field and its gradient observed at the point $(x,y,z)$⁠.

Given N observations points in which we have measured f and its gradient (for a total $4N$ data), we can define the system of N equations and 4 unknowns

Both Euler deconvolution and Euler inversion aim to solve the equation system above to estimate the parameter vector

2.1 Euler deconvolution

Euler deconvolution starts by rearranging eq. (2) to place the parameters on the left-hand side and all other terms on the right-hand side. This is an attempt to form a parametric model which results in the equation system

which can be written in matrix form as

in which $\mathbf{A}$ is the Jacobian matrix of Euler’s equation (eq. 1) concerning the parameters (eq. 3) and $\mathbf{c}$ is a pseudo-data vector.

The solution proposed by Thompson (1982) and Reid et al. (1990) is a least-squares estimate of $\mathbf{p}$

The covariance matrix of the parameters ${\mathbf{C}}_p$ is obtained through standard error propagation assuming that the only variable with uncertainty is the pseudo-data vector $\mathbf{c}$

in which ${\hat{\sigma }}_0^2=\Vert \mathbf{c}-\mathbf{A}\mathbf{p}{\Vert }^2/(N-4)$ is the reduced chi-squared statistic and an estimate of the variance factor of $\mathbf{c}$⁠.

The solution in eq. (6) above is valid only if the contents of the Jacobian matrix $\mathbf{A}$ are assumed to be error-free. As can be seen from eq. (5), the Jacobian contains the derivatives of f, which are often computed numerically by finite-differences or Fourier transforms and are known to amplify the high-frequency random noise in the data. This presents a problem, particularly for the estimation of $z_o$⁠, which has been widely explored in the literature (Silva et al. 2001; Pašteka et al. 2009; Florio et al. 2014; Melo & Barbosa 2020).

2.2 Euler inversion: formulation

Euler inversion starts by assigning the potential field f to a $N\times 1$ vector

We can then define a $4N\times 1$  data vector which contains all of the values of f and its gradient

in which ${\mathrm{\nabla }}_{\alpha }$ is the gradient operator in the $\alpha$ dimension.

Next, we formulate the $N\times 4$ equation system from Euler’s equation (eq. 2) as a nonlinear function of both parameters and data

which is known in geodesy as an implicit mathematical model (Vaníček & Krakiwsky 1986).

We then wish to solve the following constrained optimization problem with nonlinear equality constraints to estimate both the parameters and the predicted data simultaneously

in which ${\mathbf{d}}^o$ is the observed data vector which contains all of the $4N$ observations of f and its gradient, $\mathbf{d}$ is the predicted data vector from eq. (9), and $\mathbf{W}$ is a $4N\times 4N$ diagonal weight matrix. The first N terms of the diagonal of $\mathbf{W}$ are the weights for the potential field observations and the following $3N$ terms are the weights of x-, y-, and z-derivatives of the potential field, in order. In practice, the weight matrix is usually diagonal because covariances of observations are seldom available. The weights can be used to reduce the importance of each datum in the fitting process, allowing for mitigation of outliers or data components that have higher levels of noise.

The constrained problem in eq. (11) can be transformed into an unconstrained problem by using the Lagrangian

in which $\mathit{\lambda }$ is an $N\times 1$ vector of Lagrange multipliers. The non-linear Lagrangian is minimized through Newton’s method (Aster et al. 2018). We start with initial estimates ${\mathbf{p}}_0$ and ${\mathbf{d}}_0$ and then iteratively apply corrections ${\mathbf{\Delta }\mathbf{p}}_k$ and ${\mathbf{\Delta }\mathbf{d}}_k$ until convergence is achieved. To calculate the corrections, we introduce a new variable $\mathbf{u}=[{\mathbf{d}}^T{\mathit{\lambda }}^T{\mathbf{p}}^T{]}^T$⁠, expand the Lagrangian $\mathcal{L}(\mathbf{u})$ (eq. 12) in a Taylor series around point ${\mathbf{u}}_k$⁠, and disregard terms of order higher than two

in which $\mathrm{\nabla }$ is the gradient operator and ${\mathbf{H}}_k$ is the Hessian matrix of $\mathcal{L}$ evaluated at ${\mathbf{u}}_k$⁠. eq. (13) is a quadratic function of ${\mathbf{\Delta }\mathbf{u}}_k$ and we can obtain its minimum by taking its gradient and equating it to the null vector

The equation above is the system of normal equations, which can also be written in terms of $\mathbf{p}$⁠, $\mathit{\lambda }$ and $\mathbf{d}$

in which ${\mathrm{\nabla }}_{\alpha }$ is the gradient operator with respect to variable $\alpha$ and ${\mathbf{H}}_k^{\alpha \beta }$ is the Hessian matrix of $\mathcal{L}$ with respect to variables $\alpha$ and $\beta$⁠, evaluated at ${\mathbf{u}}_k$⁠. Since the order of derivation can be swapped in the Hessian matrices and the Hessian of $\mathcal{L}$ is symmetric, the above equation can be simplified to

Now, we must derive the three gradient vectors and six Hessian matrices in eq. (16). We will start with the gradient vectors.

in which ${\mathbf{e}}_k=\mathbf{e}({\mathbf{p}}_k,{\mathbf{d}}_k)$ (eq. 10), ${\mathbf{A}}_k$ is the $N\times 4$  parameter Jacobian matrix of Euler’s equation (eq. 5) evaluated on $({\mathbf{p}}_k,{\mathbf{d}}_k)$⁠, and ${\mathbf{B}}_k$ is the $N\times 4N$  data Jacobian of Euler’s equation, also evaluated on $({\mathbf{p}}_k,{\mathbf{d}}_k)$⁠. The data Jacobian ${\mathbf{B}}_k$ contains the first derivatives of Euler’s equation (eq. 1) with respect to the data vector $\mathbf{d}$ (eq. 9). It is composed of four diagonal matrices

The diagonal elements of each of the four matrices are

The Hessian matrices are calculated using a Gauss–Newton approximation disregarding second-order derivatives. The six independent Hessians are given by

Substituting the gradients (eq. 17) and Hessians (eq. 20) into the system of normal equations of Newton’s method (eq. 16) we arrive at

Since the data weight matrix $\mathbf{W}$ is diagonal and invertible, we can use the following identity to eliminate one equation from the equation system above (Wells & Krakiwsky 1971)

Applying the identity above to eq. (21) with $\mathbf{g}={\mathbf{\Delta }\mathbf{d}}_k$ and $\mathbf{h}={\left[\mathbf{\Delta }{\mathit{\lambda }}_k^T{\mathbf{\Delta }\mathbf{p}}_k^T\right]}^T$ leads to

in which ${\mathbf{Q}}_k={\mathbf{B}}_k{\mathbf{W}}^{-1}{\mathbf{B}}_k^T$ and ${\mathbf{r}}_k=[{\mathbf{d}}^o-{\mathbf{d}}_k]$ is the residual vector. Applying the identity in eq. (22) once more to the equation system above leads to a solution for the parameter correction vector

We can obtain an expression for $\mathbf{\Delta }{\mathit{\lambda }}_k$ as a function of ${\mathbf{\Delta }\mathbf{p}}_k$ from eq. (23), which results in

Finally, we can substitute the expression above into the first equation of the system of normal equations (eq. 21) to obtain the data correction as a function of ${\mathbf{\Delta }\mathbf{p}}_k$

It is worth mentioning that the Lagrange multipliers ${\mathit{\lambda }}_k$ and their corrections $\mathbf{\Delta }{\mathit{\lambda }}_k$ are not explicitly calculated during the optimization. They are a mathematical tool for enforcing the equality constraints and have no evident interpretation.

The covariance matrix of $\mathbf{p}$ is used to rank and filter solutions during the moving window procedure. It can be estimated by propagating uncertainties from the observed data ${\mathbf{d}}^o$ to the parameter correction vector (eq. 24) and, hence, to the parameter vector (Wells & Krakiwsky 1971). The propagation requires the observed data covariance matrix ${\mathbf{C}}_d$⁠, which can be approximated by ${\mathbf{C}}_d={\hat{\sigma }}_0^2{\mathbf{W}}^{-1}$⁠. Recalling that matrix $\mathbf{Q}$ is diagonal, the parameter covariance matrix is estimated at the last iteration of the Gauss–Newton method (iteration L) as

in which ${\hat{\sigma }}_0^2=\Vert {\mathbf{d}}^o-{\mathbf{d}}_L{\Vert }^2/(4N-4)$ is the reduced chi-squared statistic of the Euler inversion and an estimate of the variance factor of the observed data ${\mathbf{d}}^o$⁠.

2.3 Euler inversion: practical implementation

2.3.1 Initial estimates and convergence

Unlike a traditional Gauss–Newton inversion of a parametric model, the Euler inversion procedure estimates corrections to both the parameter vector $\mathbf{p}$ and the predicted data vector $\mathbf{d}$ at each iteration. Hence, the optimization requires initial values for both the parameters and the predicted data. The initial value of the parameters is taken as the solution of traditional Euler deconvolution ${\mathbf{p}}_0={\left[{\mathbf{A}}^T\mathbf{A}\right]}^{-1}{\mathbf{A}}^T\mathbf{c}$ (eq. 6). The initial value for the predicted data should be close to the observed data. We found that in practice a reasonably fast convergence is achieved by assigning ${\mathbf{d}}_0=0.9{\mathbf{d}}^o$⁠.

Convergence of the solution cannot be directly evaluated by the value of the Lagrangian (eq. 12) because values $\mathit{\lambda }$ are not calculated. Instead, we specify a merit function  $\mathcal{M}$ which combines the data misfit as well as the adherence to the constraints

in which $\sqrt{{\mathbf{r}}_k^T\mathbf{W}{\mathbf{r}}_k}$ is the weighted root-mean-squared error (WRMSE) and $\nu$ is a trade-off parameter that balances fitting the data and strict adherence to the constraints. In practice, we have found that a value of $\nu =0.1$ works well in all of our synthetic data tests and our field data application. The merit function is evaluated at every iteration. The nonlinear optimization stops when a given maximum number of iterations is reached, the merit function increases, or when the change in its value drops below a given threshold.

An outline of the entire Euler inversion procedure is given in Algorithm 1. Note that eqs (24) and (26) for calculating ${\mathbf{\Delta }\mathbf{p}}_k$ and ${\mathbf{\Delta }\mathbf{d}}_k$ do not depend on ${\mathit{\lambda }}_k$ or $\mathbf{\Delta }{\mathit{\lambda }}_k$⁠. Thus, eq. (25) does not need to be calculated in practice.

2.3.2 Structural index estimation

An advantage of Euler inversion over Euler deconvolution is its ability to obtain predicted values of the potential field and its gradient. In Section 3.2, we demonstrate that the weighted root-mean-squared error

of the predicted data at the L-th iteration ${\mathbf{d}}_L$ appears to be smallest when the correct structural index $\eta$ is used. Given this observation, we can estimate the optimal value of $\eta$ by running the Euler inversion in a given data window for different values of $\eta$ and choosing the one that produces the smallest WRMSE. This procedure is summarized in Algorithm 2. In all of our synthetic data tests, we used the range ${\eta }_{\min }=0$ and ${\eta }_{\max }=3$⁠.

2.3.3 Moving window procedure

For cases with multiple sources in a given data set, we adopt a moving window procedure similar to the classic Euler deconvolution. We divide the data region into M overlapping windows. For each window, we run Algorithm 2 to obtain an estimate of the parameters and the structural index $\eta$⁠. This procedure leads to spurious solutions, much like standard Euler deconvolution, in cases where there are no sources inside windows or when sources are heavily truncated. To filter out spurious sources, we rank the solutions for each structural index separately by the variance of the $z_o$ estimate, which can be obtained by selecting the corresponding element from the covariance matrix ${\mathbf{C}}_p$ (eq. 27), and keep only a given percentage of those with the smallest variance. This procedure is summarized in Algorithm 3.

2.3.4 Choice of data weights

The data weight matrix $\mathbf{W}$ (eq. 11) can be used to assign different weights to different observations. This matrix should be diagonal, with each diagonal element corresponding to the weights of the corresponding datum. Weights should be normalized to the range $[0,1]$⁠. Assigning a weight smaller than 1 to a datum will result in a larger residual for that datum at the end of the inversion procedure. Hence, weights can be used to mitigate outliers in the data by assigning a small weight to them. Different weights can also be assigned to the potential field and its derivatives. Assigning a weight smaller than 1 to the derivatives will cause the inversion to prioritize fitting the observed potential field instead of the, often noisier, derivatives.

In practice, data weights can be calculated based on known data uncertainties or determined by trial and error. We have found that the following weights work well in most of the applications we have undertaken: 1 for the total-field anomaly, 0.1 for the east-derivative, 0.1 for the north-derivative and 0.025 for the upward-derivative. The decrease in weight for the derivatives helps mitigate the effect of high-frequency noise, which is amplified by the numerical derivation, on the inversion estimate.

3 RESULTS

In this section, we demonstrate the effectiveness and limitations of the Euler inversion method by applying it to a series of synthetic data sets and to real aeromagnetic data from Rio de Janeiro, Brazil. The applications are organized as follows:

• Method demonstration: This test uses a single data window and a single dipolar source. Its aim is to demonstrate the convergence of the Euler inversion method and its ability to correctly estimate the source position when the structural index is known, even in the presence of random noise.

• Effect of structural index choice: This test uses several different sources, each in a separate data window, and runs the Euler inversion method on each with different values of the structural index $\eta$⁠. Its aim is to determine the effect of the choice of $\eta$ on the estimated coordinates and the weighted root-mean-squared error (eq. 28).

• Effect of random noise: This test uses a single dipolar source and a single data window with data contaminated with increasing levels of pseudo-random noise. Its aim is to investigate the effect of random high-frequency noise on the Euler inversion estimated source coordinates, base level and structural index.

• Effect of interfering sources: This test uses models composed of two sources at increasing distances from each other. The test is performed for two dipolar sources and also two dykes. Its aim is to investigate the effect of interfering sources inside the data window on the Euler inversion, Euler deconvolution and finite-difference Euler deconvolution results. We use a single data window so that we can understand what happens at each step of the moving window procedure.

• Moving window procedure with multiple sources: This test combines several sources and uses the moving window procedure from Algorithm 3. Its aim is to show how the Euler inversion method behaves on a more complex data set and provide a comparison with standard Euler deconvolution and finite-difference Euler deconvolution.

• Aeromagnetic data from Rio de Janeiro, Brazil: This test applies the Euler inversion method to a real data set which contains multiple sources. Its aim is to demonstrate the effectiveness of the method on a real data set with realistic levels of noise, signal overlap and geometry of sources.

The Python source code used to produce the results presented here, as well as extra explanation of the models and procedures, can be found in the supplementary information at https://doi.org/10.6084/m9.figshare.26384140 (Uieda et al. 2024).

3.1 Method demonstration

The main goal of this synthetic data test is to demonstrate the general effectiveness of the Euler inversion method to estimate the position and base level of a single source. To this end, we created a model composed of a single dipole located at $(x_o=15000\mathrm{m},y_o=12000\mathrm{m},z_o=-3000\mathrm{m})$ with a dipole moment magnitude of $5\times {10}^{11}{\mathrm{A}\mathrm{m}}^{-1}$⁠, inclination of –30${}^{\circ }$⁠, and declination of 15${}^{\circ }$⁠. The reference field direction was the same as the dipole moment direction. The synthetic total-field magnetic anomaly data was calculated on a regular grid with point spacing of 300 m at a height of 800 m. To the data, we added a base level of 100 nT and pseudo-random Gaussian noise with 0 nT mean and 10 nT standard deviation. The eastwards and northwards derivatives of the total-field anomaly grid were calculated with a central-difference scheme. The upward derivative was calculated by Fast Fourier Transform (FFT). The synthetic anomaly and its three derivatives are shown in Figs 1(a)–(d).

The Euler inversion method described in Algorithm 1 was applied to the synthetic data. We chose a fixed structural index of $\eta =3$⁠, which is the correct index for a magnetic dipole. For data weights, we used 1 for the total-field anomaly, 0.1 for the east-derivative, 0.1 for the north-derivative and 0.025 for the upward-derivative. These weights were chosen to counteract the increased effect of noise on the derivatives, particularly the upward derivative which was calculated through FFT. Figs 1(e)–(h) show the inversion residuals after convergence was achieved (⁠$L=6$ iterations) for the total-field anomaly and its eastwards, northwards, and upwards derivatives, respectively. Also shown are the true source location, the initial source location and the predicted source location from Euler inversion. The initial estimate of the source location was $(x_o=14626\mathrm{m},y_o=11865\mathrm{m},z_o=-1553\mathrm{m})$ and the base level was $b=94\mathrm{nT}$⁠, which are the Euler deconvolution results. The final Euler inversion prediction of the source location was $(x_o=15045\mathrm{m},y_o=12028\mathrm{m},z_o=-2663\mathrm{m})$ and the estimated base level was $b=93\mathrm{nT}$⁠, which is an improvement on the estimated values by Euler deconvolution (Fig. 1i).

Fig. 1(i) shows the error in the estimated source coordinates and base level. We can see from the figure that the error in the $x_o$ (easting) and $y_o$ (northing) coordinates, as well as the base level, do not vary greatly from the initial solution at iteration 0. However, the error in the $z_o$ (upward) coordinate drops from over 1400 m to less than 400 m in two iterations. The merit function (eq. 28) also drops sharply in value by two iterations, as can be seen in Fig. 1(j), confirming the rapid convergence of the Euler inversion method.

3.2 Effect of structural index choice

In this synthetic data test, we created data sets using four different models: a dipole, a horizontal cylinder composed of a right-rectangular prism stretched in the southwards direction, a vertical pipe composed of a right-rectangular prism stretched in the downwards direction, and a vertical dyke composed of a right-rectangular prism stretched in the southwards, northwards and downwards directions. All models share the same true location of $(x_o=15000\mathrm{m},y_o=10000\mathrm{m},z_o=0\mathrm{m})$⁠, base level of 300 nT, and induced magnetization with inclination of 35${}^{\circ }$⁠. The data were generated on a regular grid with spacing of 300 m, height of 1000 m and contaminated with pseudo-random Gaussian noise with 0 nT mean and 15 nT standard deviation. Figs 2(a)–(d) show the synthetic noise-corrupted total-field anomaly data.

We ran the Euler inversion method on each data grid four times, each time changing the structural index between zero, one, two and three. Fig. 2(e) shows the upwards coordinate $z_o$ estimated for each of the four models as a function of the structural index $\eta$⁠. The Euler inversion estimated $z_o$ correlates with $\eta$⁠, with larger values of the structural index leading to deeper source estimates. Values closest to the true $z_o=0\mathrm{m}$ are achieved when the correct structural index is used (⁠$\eta =1$ for the dyke, $\eta =2$ for the cylinder and pipe, and $\eta =3$ for the dipole). Fig. 2(f) shows the WRMSE (eq. 29) at the final iteration of the Euler inversion method for all four models as a function of structural index. The WRMSE is a measure of goodness-of-fit between the predicted total-field anomaly and its three derivatives and their observed counterparts. The WRMSE is minimum for all four models when the correct structural index is used.

3.3 Effect of random noise

We conducted another experiment to determine the effect of random high-frequency noise on the Euler inversion estimates. To this end, we created synthetic data from a dipole model located at $(x_o=15000\mathrm{m},y_o=11000\mathrm{m},z_o=-5000\mathrm{m})$ and with a dipole moment magnitude of $2\times {10}^{12}{\mathrm{A}\mathrm{m}}^{-1}$⁠, inclination of –30${}^{\circ }$⁠, and declination of 15${}^{\circ }$⁠. The total-field anomaly data were generated on a regular grid with a spacing of 500 m and a constant height of 800 m. The reference field direction was the same as the dipole moment direction. A base level of 100 nT was added to the data. We generated different data sets by adding pseudo-random Gaussian noise with 0 nT mean and standard deviations varying from 0 nT to 40 nT with a step of 0.2  nT. Figs 3(a)–(d) show the synthetic data for noise levels 0, 10, 25 and 40 nT, while Figs 3(e)–(h) show the upwards derivative calculated from the total-field anomaly through FFT.

On each data set, we ran Euler deconvolution (eq. 6), Euler inversion with unit weights, and Euler inversion with weights 1 for the total-field anomaly, 0.1 for the eastwards derivative, 0.1 for the northwards derivative and 0.025 for the upwards derivative. Both Euler inversion runs used the structural index estimation procedure (Algorithm 2 with ${\eta }_{\min }=0$ and ${\eta }_{\max }=3$⁠). Figs 2(i)–(l) show the error in the estimated easting, northing and upwards coordinates as well as the base level for each of the methods as a function of noise level. The error in each of three coordinates raises sharply with noise level for Euler deconvolution, particularly for the upwards $z_o$ coordinate. The unweighted Euler inversion results vary less regularly but the present errors are just as large as Euler deconvolution for the upwards coordinate. However, the weighted Euler inversion presented overall smaller errors and a slower growth curve for the upwards coordinate error than the other two methods. The base level error is nearly constant at approximately 10 nT for Euler deconvolution and the weighted Euler inversion, but varies to as much as 40 nT for the unweighted Euler inversion.

Fig. 3(m) shows the estimated structural index $\eta$ for the weighted and unweighted Euler inversion as a function of noise level. The unweighted Euler inversion estimated the wrong structural index $\eta =2$ from approximately noise level 7 nT and $\eta =1$ from approximately noise level 20 nT. These jumps in the estimated structural index appear to correlate with jumps in the base level and $z_o$ coordinate errors. The weighted Euler inversion was able to estimate the correct structural index (⁠$\eta =3$⁠) for all noise levels tested.

3.4 Effect of interfering sources

Another common issue encountered during the application of Euler-based methods is the presence of interfering sources within the data window. To test this effect on Euler inversion, we create two different scenarios, one with two dipoles and another with two dykes. In both scenarios, we created several synthetic total-field anomaly data sets by varying the distance between the two sources. No noise was added to these synthetic data in order to isolate the effect of the interfering source from the effect of random noise. We added to all data sets a base level of 100 nT. On each data set, we ran Euler deconvolution (eq. 6) with the correct structural index for the source (⁠$\eta =3$ for the dipoles and $\eta =1$ for the dykes), the finite-difference Euler deconvolution method of Gerovska et al. (2005), and Euler inversion with the structural index estimation (Algorithm 2 with ${\eta }_{\min }=0$ and ${\eta }_{\max }=3$⁠) and data weights of 1 for the total-field anomaly, 0.1 for the eastwards derivative, 0.1 for the northwards derivative and 0.025 for the upwards derivative.

The dipole models contain a main dipole located at $(x_o=7000\mathrm{m},y_o=4000\mathrm{m},z_o=-3000\mathrm{m})$ with a dipole moment amplitude of $5\times {10}^{11}{\mathrm{A}\mathrm{m}}^{-1}$⁠, inclination of –30${}^{\circ }$ and declination of –10${}^{\circ }$⁠. The interfering dipole was located at $y_o=5000\mathrm{m}$ and $z_o=-1500\mathrm{m}$⁠, with the $x_o$ coordinate varying from –1000 m to 5000 m. The reference field direction was the same as the dipole moment direction. The total-field anomaly data were generated on regular grids with a spacing of 200 m and at a constant height of 400 m. Figs 4(a)–(d) show the total-field anomaly of four out of the 31 models.

The error in the estimated eastwards, northwards and upwards coordinates are shown in Figs 4(e)–(g). For the eastwards and northwards coordinates, the three methods are mostly compatible, with Euler inversion being slightly closer to the true source for most distances between sources. For the upward coordinate, Euler inversion and Euler deconvolution are roughly equivalent and both have smaller errors than finite-difference Euler deconvolution for all but the largest distances. For the structural index estimates (Fig. 4h), finite-difference Euler deconvolution underestimates $\eta$ for all but the largest distances, while Euler inversion estimates the correct index of $\eta =3$ for all distances.

The dyke models contain a main dyke at the east with a centre point at $(x_o=7000\mathrm{m},y_o=4500\mathrm{m})$ and a top at $z_o=0\mathrm{m}$ with a dipole moment amplitude of $2\times {10}^1{\mathrm{A}\mathrm{m}}^{-1}$⁠, inclination of –30${}^{\circ }$ and declination of 20${}^{\circ }$⁠. The interfering dyke has a top at $z_o=300\mathrm{m}$⁠, with the $x_o$ coordinate varying from –2000 to 6000 m. The reference field direction was the same as the dipole moment direction. The total-field anomaly data were generated on regular grids with a spacing of 150 m and at a constant height of 400 m. Figs 5(a)–(d) show the total-field anomaly of four out of the 33 models.

The error in the estimated eastwards, northwards and upwards coordinates are shown in Figs 5(e)–(g). The error for the eastwards and northwards coordinates was calculated with respect to the centre point of the eastern dyke. For the eastwards and northwards coordinates, the three methods are mostly compatible. When the two dykes intersect, all three methods estimate a horizontal position at the intersection point. When they do not intersect, the estimates for all three methods fall outside of the data window. This is a well known issue for Euler deconvolution methods because the Hessian matrix ${\mathbf{A}}^T\mathbf{A}$ (eq. 6) is ill-conditioned for 2-D sources (Mushayandebvu et al. 2004). For the upward coordinate, Euler inversion and Euler deconvolution are roughly equivalent and approach zero at the largest distances. Both methods have smaller errors than finite-difference Euler deconvolution for all distances tested. For the structural index estimates (Fig. 5h), finite-difference Euler deconvolution estimates incorrect $\eta$ at most distances, presenting no clear relationship between distance and $\eta$ estimate. Euler inversion estimates the correct index of $\eta =1$ for all distances.

3.5 Moving window procedure with multiple sources

To simulate a more realistic data set, we created a model composed of 10 sources combining dipoles at various locations and depths and vertical dykes at various orientations. All sources had induced magnetization in the direction of the regional field with a inclination of –30${}^{\circ }$ and declination of –20${}^{\circ }$⁠. The total-field anomaly of the model was calculated on a regular grid with a spacing of 500 m and at a constant height of 1000 m. We added to the data a base level of 1000 nT, pseudo-random Gaussian noise with 0 nT and 50 nT standard deviation, and a regional field composed of a first-degree polynomial with angular coefficients of 0.02 nT m${}^{-1}$ in the eastwards and –0.03 nT m${}^{-1}$ in the northwards directions. The noise-corrupted total-field anomaly data are shown in Fig. 6(a).

To the data set, we applied the moving window Euler inversion method (Algorithm 3 and Algorithm 2 with ${\eta }_{\min }=0$ and ${\eta }_{\max }=3$⁠), the finite-difference Euler deconvolution method of Gerovska et al. (2005), and standard Euler deconvolution (using structural indices 1, 2 and 3). Euler inversion was performed with data weights of 1 for the total-field anomaly, 0.1 for the eastwards derivative, 0.1 for the northwards derivative and 0.025 for the upwards derivative. All three methods used the same moving window procedure described in Algorithm 3 for the sake of comparison. The windows had a size of 10 000 m and were moved by 5000 m at a time. The ratio of estimates kept to form the final solution was $\gamma =0.3$ for Euler deconvolution, $\gamma =0.35$ for finite-difference Euler deconvolution, and $\gamma =0.25$ for Euler inversion.

Figs 6(b)–(f) show the estimated source positions and structural indices for finite-difference Euler deconvolution, Euler inversion and Euler deconvolution with structural indices 1, 2 and 3, respectively. The finite-difference method estimates a non-integer structural index, as a result Fig. 6(b) shows the closest integer value to the actual estimated $\eta$⁠. The finite-difference Euler deconvolution method underestimates the structural indices of all sources and, therefore, also underestimates their depths. The finite-difference method solutions are also more scattered than their Euler deconvolution and Euler inversion counterparts. The Euler deconvolution results are closer to the correct depths when the correct structural index is used. They present larger dispersion than Euler inversion in areas where the signals of multiple sources overlap. With the exception of the deeper dykes in the northwest and southeast and the small dipole with $z_o=-500\mathrm{m}$⁠, Euler inversion is able to estimate the correct structural index for most sources. The upward coordinate estimates for Euler inversion are also closer than Euler deconvolution to their true values when the correct structural index was estimated. Euler inversion notably estimates an incorrect $\eta$ and $z_o$ for smaller sources when there is a large amount of interference in the anomalies and for dykes that are deeper and produce a smoother signal.

It is also notable that there are solutions that outline the simulated dykes for all three methods. This seems to be in contradiction of the results presented in Section 3.4 and the theoretical proof in Mushayandebvu et al. (2004), which show that the $x_o$ and $y_o$ coordinates cannot be estimated for 2-D sources. However, these results are widely known in practice where Euler deconvolution-based methods are routinely used to map dykes and lineaments. We believe that this is the effect of high-frequency noise. The derivatives, which appear in the Jacobian matrix (eq. 5), will contain this high-frequency noise as well which is not 2-D in nature. This in turn causes the Hessian matrix (eq. 6) to not be singular in practice.

3.6 Aeromagnetic data from Rio de Janeiro

The geology of Rio de Janeiro state (Southeastern Brazil) consists primarily of high-grade metamorphic rocks and granitoid magmatism related to the Ribeira Belt (Heilbron et al. 2020). Fig. 7(a) shows a simplified geologic map of the area, which was modified from Heilbron et al. (2016) and Dantas et al. (2017). The Ribeira Belt is traditionally interpreted as a thrust belt formed by diachronous collisions mainly between the São Francisco and Congo paleocontinents (Trouw et al. 2000; Heilbron et al. 2008) or by an intracontinental orogeny (e.g. Meira et al. 2015, 2019), during the Brasiliano orogeny. This process culminated in an orogen-parallel, steep strike-slip shear system (Egydio-Silva et al. 2005), which deformed the Paleoproterozoic basement rocks and reworked the Meso- to Neoproterozoic metasedimentary units (for example, the Italva and São Fidelis groups) and syn-orogenic granitoid plutons (for example, the Rio Negro complex) which formed during the orogeny (Heilbron & Machado 2003; Heilbron et al. 2020). These tectonic events imprinted a distinct NE–ENE-trending structural pattern onto these rocks.

The late Neoproterozoic to Cambrian period witnessed post-orogenic magmatism (e.g. Valeriano et al. 2011), marking the final stages of the West Gondwana amalgamation. After this, the region remained tectonically quiescent until the Lower Cretaceous, when reactivation occurred with the emplacement of the NE-trending Serra do Mar mafic dyke swarm, preceding the break-up of West Gondwana and the opening of the South Atlantic Ocean (Almeida et al. 2013). Lastly, thermal anomalies in the region during the Upper Cretaceous to Paleocene period led to the emplacement of alkaline complexes and dykes (Thompson et al. 1998). The geological complexity of the Ribeira Belt, marked by the interplay of diverse tectonic regimes and magmatic events (Fig. 7a), makes the Rio de Janeiro region an ideal test case for Euler inversion.

We used aeromagnetic data from the state of Rio de Janeiro which are distributed by the Serviço Geológico do Brasil (https://geosgb.sgb.gov.br). The data were collected in two phases: Subarea 1 was surveyed between 1978 March 25 and May 27, using an Islander aircraft (PT-KRP), while Subarea 2 was surveyed between 1978 April 6 and July 19, using a Bandeirante aircraft (PT-GKJ), both funded by the Brazilian government. As shown in Fig. 7(b), the survey followed a pattern of north–south flight lines spaced approximately 1 km apart, with east–west tie lines. Data were recorded at 100-m intervals using a Geometrics G-803 magnetometer. Some of the notable features of the data are the NE–SW linear features (interpreted here as dykes), which coincide with known dyke outcrops, and complex dipolar anomalies which coincide with some of the post-collisional magmatism and alkaline intrusions. A subset of 50 882 data points were used in our analysis.

The data were not interpolated on a regular grid to avoid any smoothing effects that the interpolation might have on the linear features. This could result in an over-estimation of their depth, as discussed in Section 3.5. Instead, we used the gradient-boosted equivalent sources method of Soler & Uieda (2021) to fit a model to the observed line data. We then used the model to make predictions of the three spatial derivatives at the original measurement locations by a central-difference method with a coordinate shift of 1 m. Further details about the data processing can be found in the source code archive that accompanies this article https://doi.org/10.6084/m9.figshare.26384140 (Uieda et al. 2024).

We performed the moving-window Euler inversion (Algorithm 3 and Algorithm 2 with ${\eta }_{\min }=1$ and ${\eta }_{\max }=3$⁠) on the observed total-field anomaly line data using windows of size of 12 000 m which were moved 2400 m at a time. The proportion of solutions kept was $\gamma =0.15$⁠. The inversion was performed with data weights of 1 for the total-field anomaly, 0.1 for the eastwards derivative, 0.1 for the northwards derivative and 0.05 for the upwards derivative. To aid in the geological interpretation of the results, we converted the estimated upwards source coordinates $z_o$ to depths below the surface of the Earth. We did so by subtracting the estimated $z_o$ from the interpolated topographic height of the Shuttle Radar Topography Mission (NASA JPL 2013). The estimated positions and structural indices are shown in Fig. 8. We also performed Euler deconvolution with structural indices one, two and three, as well as finite-difference Euler deconvolution on the same data set using the same window size and window step for the sake of consistency. These results are shown in Fig. 9.

The Euler inversion estimated source positions shown in Fig. 8 highlight the NE–SW lineaments as well as some of the more dipolar anomalies. The lineaments are estimated with a mix of $\eta =1$⁠, $\eta =2$ and $\eta =3$⁠. The southernmost lineament is mostly estimated with $\eta =1$ and depths suggesting that it does not outcrop in its southernmost parts (depths of 400 to 600 m), which is consistent with the geologic information in Fig. 7(a). The southernmost part of this lineament, in particular, has an estimated $\eta =3$⁠, which is known to happen for deeper dykes in our synthetic data tests (Section 3.5). Conversely, the northernmost part of the lineament has a larger prevalence of $\eta =1$ with shallower depths which coincide with a known dyke outcrop. Other known dyke outcrops coincide with estimated sources with $\eta =1$⁠, however their depths range from 100 to 300 m. This may be caused by an excess of smoothing in the vertical derivative or effects of noise in the estimated coordinates. The lineaments in the northwestern part of the region are also highlighted by estimated sources. However, their structural indices are a mix of $\eta =2$ and $\eta =3$⁠, suggesting deeper sources. This is inline with the geologic information, which includes no outcrops of linear structures in the area.

The dipolar anomalies are associated with post-collisional and alkaline intrusions, many of which are also cut by known outcropping dykes or have known dykes with magnetic signals that significantly overlap with the dipolar anomalies. The Euler inversion estimated structural indices for them range from $\eta =2$ to $\eta =3$⁠. We have highlighted four dipolar anomalies, marked as A, B, C and D in Fig. 8, to aid in our discussion.

• Anomaly A: Has a reversed polarity and linear feature to its north that is not associated with any known dyke outcrop. The linear feature is highlighted by Euler inversion estimates with $\eta =1$ and depth of 300–400 m, which can be interpreted as a non-outcropping dyke. The dipolar anomaly itself has Euler inversion solutions with $\eta =3$ and depth of 1000–2000 m. The solutions in the centre of the anomaly present a shallower depth than the solutions to the north and south of the anomaly centre. From the results on synthetic data in Section 3.5, we can interpret the depth range to be caused by the moving window procedure and the effect of interfering sources. The depth to the centre of the anomaly source is likely close to 1000 m.

• Anomaly B: The dipolar anomaly is likely associated with a non-outcropping portion of the post-collisional magmatism. The anomaly is cut by several NE–SW linear features, some of which overlap with known dyke outcrops. The linear feature to the north is associated with Euler inversion results with $\eta =1$ and depths ranging from 300 to 600 m, suggesting a non-outcropping dyke. At the centre of the anomaly are Euler inversion estimates with $\eta =3$ and depth estimate of approximately 1400 m. The Euler inversion solutions surrounding these central solutions are likely caused by interference from other sources.

• Anomaly C: A dipolar anomaly associated with an outcropping portion of the post-collisional magmatism. There is a known outcropping dyke to the south of the anomaly, which is associated with Euler inversion estimates with $\eta =1$ and depths ranging from 500 to 1000 m. These depth estimates are likely overestimated because of the interference of the dipolar anomaly. The main anomaly has Euler inversion solutions with $\eta =2$ and $\eta =3$ and depths varying from 1400 to 1800 m. There is no clear indication of which of these estimates is more reliable.

• Anomaly D: A small dipolar anomaly associated with an outcropping alkaline intrusion. The Euler inversion estimates have $\eta =3$ and depths 1700–2000 m. There are known outcropping dykes around the main intrusion but they have no discernible magnetic anomalies and no Euler inversion solutions associated with them.

Overall, the Euler inversion solutions in Fig. 8 are consistent with the known geology in Fig. 7(a). The main linear features are mostly associated with Euler inversion estimates with $\eta =1$ and shallow depths, particularly where known dyke outcrops are located. Deeper linear features are estimated with $\eta =2$ and $\eta =3$⁠, which is consistent with the synthetic data results (Section 3.5). The dipolar anomalies have consistent Euler inversion estimates with $\eta =3$ when they are well isolated from interfering sources. Otherwise, they are estimated with a mix of structural indices and depths, as was demonstrated in Section 3.5.

When compared to the Euler deconvolution and finite-difference Euler deconvolution results (Fig. 9), the Euler inversion results are less dispersed and better delinear the linear features present in the data. The structural index results from finite-difference Euler deconvolution (Fig. 9d) are underestimated, with most values of $\eta$ being less than one, which is not in accordance with the geology of the area.

4 CONCLUSION

Euler deconvolution is a widely used method for locating the sources of potential-field data. It has its limitations in real-world scenarios due to its dependence on the chosen value of the structural index $\eta$ and its sensitivity to high-frequency noise and signal overlap from interfering sources. We have developed a new method to solve Euler’s homogeneity equation for the source position, base level and integer structural index, which we call Euler inversion. Unlike Euler deconvolution, Euler inversion is also able to estimate the predicted field and its spatial derivatives, as well as assign different weights to each type of data. Our method can be applied to gridded and non-gridded data, which can be useful to limit the effects of smoothing from interpolation in the final results when the original data have large spacing between flight lines. The Euler inversion algorithm is computationally efficient because most of the large matrices involved in the computations are diagonal or block-diagonal. We found that, in practice, the computation time of Euler inversion and Euler deconvolution are on the same order of magnitude.

Tests on synthetic data show that Euler inversion outperforms Euler deconvolution and finite-difference Euler deconvolution (a variant that estimates $\eta$ but does not rely on second-order derivatives) in terms of robustness to random noise and interfering sources inside the data window. Our tests also show that the estimated $z_o$ coordinate is correlated with the structural index, as is the case for Euler deconvolution. We have also found that the data misfit from Euler inversion is minimal when the integer structural index used is equal or close to the true one for idealized sources. This led us to develop an algorithm for estimating the best integer structural index based on the data misfit. A test on complex synthetic data from a model of dykes and dipoles with overlapping signals shows that Euler inversion is able to estimate the structural index and position of the sources within expected error bounds when the signal overlap is not larger than the data window. For deeper dykes in particular, Euler inversion was not able to estimate the correct $\eta =1$⁠, leading to an overestimation of the depths.

We applied Euler inversion to an aeromagnetic dataset from Rio de Janeiro, Brazil, to analyse its performance under real-world scenarios. Euler inversion was able to locate the NE–SW linear features in the data with an $\eta =1$ which are associated with known dyke outcrops. For the deeper linear features, Euler inversion was not able to estimate the correct $\eta =1$⁠. Some of the dipolar anomalies present in the data were picked out with $\eta =3$⁠, while the sources with a large signal overlap with other features provided a mix of $\eta =2$ and $\eta =3$⁠. These results are consistent with the synthetic data tests and show the benefits and limitations of the proposed method.

Euler inversion outperforms Euler deconvolution and finite-difference Euler deconvolution in most cases. Its reduced sensitivity to noise and interfering sources, in particular, may prove beneficial for magnetic microscopy studies, in which high-frequency noise and interference from multiple dipolar sources are a significant hurdle (Souza-Junior et al. 2024). However, it still suffers from some of the same limitations. While Euler inversion is less sensitive to signal overlap, it still fails to correctly estimate the position and structural index when the overlap is large. The windowing procedure still generates a large amount of spurious solutions which need to be filtered out. This could be improved with techniques like the source detection method proposed by Castro et al. (2020), for example. Euler inversion can also be coupled with other inverse problems by following our methodology to add Euler’s equation as a non-linear constraint. This could help with issues of non-uniqueness and stability in traditional 3-D inverse problems in potential-field methods.

As is the case with other Euler deconvolution-based methods, Euler inversion also suffers from instability when sources are 2-D (Mushayandebvu et al. 2004) and a lack of support for sources that are defined by multiple points, for example, steps which have a top and bottom (Gerovska et al. 2010). The issue of instability was evident in our synthetic data tests which simulated dykes and did not contaminate the data with random noise. However, both in the synthetic data tests with noise and the real data application, this did not appear to be a significant issue. In both cases, Euler inversion was better at outlining the 2-D sources than Euler deconvolution. Nonetheless, it would be worthwhile to investigate this issue further and explore the use of regularization and an adaptation of the method of Mushayandebvu et al. (2004) to the Euler inversion mathematical formulation. A thorough comparison of Euler inversion with the similarity transform-based methods, like Gerovska et al. (2010), would also be worth pursuing. Euler inversion could complement such methods, which often require gridded data, in cases where data are of poor quality or flight-line spacing is large.

ACKNOWLEDGMENTS

We are indebted to the developers and maintainers of the open-source software without which this work would not have been possible. We thank Dr Valéria C. F. Barbosa for many insightful discussions over the years which helped shape our research. We are grateful to editor Dr Kosuke Heki and reviewers Dr Alan B. Reid and Dr Saulo Oliveira for their constructive comments. LU would like to thank Prof. Spiros Pagiatakis for being an incredible instructor and teaching him the mathematics which formed the foundations of this work during his undergraduate exchange at York University. LU was supported in part by start-up grant PRPI 22.1.09345.01.2 from Universidade de São Paulo. GFS-J was supported by scholarship 2021/08379-5 from the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP). The opinions, hypotheses and conclusions or recommendations expressed in this material are the responsibility of the authors and do not necessarily reflect the views of FAPESP.

DATA AVAILABILITY

The Python source code and data that were used to produce all results and figures presented here are available at https://github.com/compgeolab/euler-inversion and https://doi.org/10.6084/m9.figshare.26384140 (Uieda et al. 2024) under the CC-BY license and the MIT license. This study made use of the following open-source scientific software: matplotlib (Hunter 2007) and PyGMT (Tian et al. 2024) for generating figures and maps, Numpy (Harris et al. 2020) and Scipy (Virtanen et al. 2020) for linear algebra, Pandas for manipulating tabular data (McKinney 2010; The pandas development team 2024), GeoPandas for reading and plotting shapefiles (den Bossche et al. 2024), pyproj for data projection (Snow et al. 2024), xarray (Hoyer & Hamman 2017) for working with gridded data, Verde (Uieda 2018) for moving windows and interpolation, and Harmonica (Fatiando a Terra Project et al. 2023) for potential-field data processing and modelling. The aeromagnetic and geologic data are available from Serviço Geológico do Brasil (https://geosgb.sgb.gov.br) under a CC-BY-NC license. The magnetic data are part of survey 1038 ‘Projeto Aerogeofísico São Paulo–Rio de Janeiro’. Both are also available in our source code and data archive (Uieda et al. 2024).

AUTHOR CONTRIBUTIONS

Leonardo Uieda (Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Visualization, Writing—original draft), Gelson Ferreira Souza-Junior (Data curation, Formal analysis, Resources, Software, Visualization, Writing—original draft, Writing—review & editing), India Uppal (Data curation, Formal analysis, Investigation, Software, Writing— review & editing) and Vanderlei Coelho Oliveira Jr. (Conceptualization, Methodology, Writing—review & editing).

REFERENCES