harmonica.EulerInversion

class harmonica.EulerInversion(*, structural_index=(1, 2, 3), max_iterations=20, tol=0.1, euler_misfit_balance=0.1)

Estimate source location, base level, and SI using Euler Inversion.

Implements Euler Inversion [Uieda2025] to estimate subsurface source location from potential field data and their directional derivatives. Also estimates any constant shifts or biases of the data (called the base level), as well as the structural index (SI; a parameter reflecting the source geometry; see below). The approach employs a non-linear total-least-squares approach to solve the inverse problem of Euler’s homogeneity equation.

Assumes a single data window and provides a single estimate.

Hint

Euler Inversion is much more stable than Euler Deconvolution. It’s less sensitive to noise in the field derivatives and to interfering sources within the data window. It can also estimate integer-valued structural indices (SI).

Note

Does not yet support structural index 0.

Parameters:

structural_indexint or sequence, optional

Defines the nature of the source of the potential field data. Should be an integer value, usually between 1 and 3 (but can vary depending on the nature of the field). If a sequence of values is passed, the inversion will try each one and pick the value that provides the smallest misfit to the observed data (default). It’s the degree of the field’s rate of change with distance from the source, influencing the decay rate of the field and the formulation of Euler’s homogeneity equation. Correlated with the depth estimate, so larger structural index will lead to larger estimated depths. Choose based on known source geometry (see table below) or allow the inversion to estimate the optimal value.

max_iterationsint, optional

The maximum number of iterations allowed in the non-linear Gauss-Newton inversion. If the value is too small, there is a risk of exiting the inversion without the solution converging to the minimum of the goal function. Larger values won’t necessarily lead to longer computation times since the inversion will stop if convergence is reached.

tolfloat, optional

The tolerance in decimal percentage that is needed to continue the iterations. If the change in the merit function (see below) is less than tol times the current merit function value, the iterations will be terminated. Use smaller values to allow for longer inversions.

euler_misfit_balancefloat, optional

The trade-off parameter \(\nu\) between fitting the data and obeying Euler’s homogeneity equation (see below).

Notes

Works on any potential field that satisfies Euler’s homogeneity equation (like gravity, magnetic, and their gradients caused by simple sources):

\[(e_i - e_0)\dfrac{\partial f_i}{\partial e} + (n_i - n_0)\dfrac{\partial f_i}{\partial n} + (u_i - u_0)\dfrac{\partial f_i}{\partial u} = \eta (b - f_i),\]

in which \(f_i\) is the given potential field observation at point \((e_i, n_i, u_i)\), \(b\) is the base level (a constant shift of the field, like a regional field), \(\eta\) is the structural index, and \((e_0, n_0, u_0)\) are the coordinates of a point on the source (for a sphere, this is the center point).

The Euler Inversion estimates \((e_0, n_0, u_0)\) and \(b\) given a potential field and its easting, northing, and upward derivatives. If the structural index is not given, it can estimate an integer valued \(\eta\) by running the inversion multiple times and choosing the \(\eta\) that produces the best fit to the data. This is a big advantage of the Euler Inversion approach over the Deconvolution since the latter is unable to calculate predicted data and thus cannot evaluate true data misfit.

The convergence of the solution is measured through a merit function

\[\mathcal{M}(\mathbf{p}, \mathbf{d}) = \sqrt{\mathbf{r}^T\mathbf{W}\mathbf{r}} + \nu\sqrt{\mathbf{e}^T\mathbf{e}}\]

in which \(\mathbf{p}\) is the parameter vector, \(\mathbf{d}\) is the predicted data vector, \(\mathbf{W}\) is the weight matrix, \(\mathbf{r}\) is the residual vector, \(\mathbf{e}\) is the evaluation of Euler’s equation using the current data and parameters, and \(\nu\) is trade-off parameter that balances fitting the data with obeying Euler’s equation.

As with Euler Deconvolution, Euler Inversion still assumes that the sources are ideal (see the table below). We recommend reading [ReidThurston2014] for a discussion on what the structural index means and what it does not mean. After [ReidThurston2014], values of the structural index (SI) can be:

Source type	SI (Mag)	SI (Grav)
Point, sphere	3	2
Line, cylinder, thin bed fault	2	1
Thin sheet edge, thin sill, thin dyke	1	0

Attributes:

location_1d-array

Estimated (easting, northing, upward) coordinates of the source after model fitting.

base_level_float

Estimated base level constant of the anomaly after model fitting.

covariance_2d-array

The 4 x 4 estimated covariance matrix of the solution. Parameters are in the order: easting, northing, upward, base level. This is not an uncertainty of the position but a rough estimate of their variance with regard to the data.

structural_index_int

The estimated structural index.

Methods

fit(coordinates, data, *[, weights])

Fit the model using potential field measurements and their derivatives.

EulerInversion.fit(coordinates, data, *, weights=(1, 0.1, 0.1, 0.025))

Fit the model using potential field measurements and their derivatives.

Solves Euler’s homogeneity equation to estimate the source location and base level by utilizing field values and their spatial derivatives in easting, northing, and upward directions. Constructs an implicit mathematical model and estimates both the parameters and the predicted data (field and its derivatives). Will also estimate the structural index if a single value was not provided.

Tip

Data does not need to be gridded for this to work.

Parameters:

coordinatestuple of arrays

Tuple of 3 with the coordinates of each data point. Should be in the following order: (easting, northing, upward). Arrays can be n-dimensional but must all have the same shape.

datatuple of arrays

Tuple of 4 arrays with the observed data in the following order: (potential_field, derivative_easting, derivative_northing, derivative_upward). Arrays can be n-dimensional but must all have the same shape as the coordinates. Derivatives must be in data units over coordinates units, for example nT/m or mGal/m.

weightstuple, list, 1d-array, optional

Weights assigned to each of the four data types (field and its derivatives) in the inversion. Reducing the weights of the derivatives helps reduce the influence of random noise in the results. By default, weights are 1 for the field, 0.1 for its eastward and northward derivatives, and 0.025 for its upward derivative. The upward derivative has a smaller weight because it usually contains more errors arising from FFT-based processing. If uncertainties are available, the weights can be 1 / uncertainty but should be normalized so that the largest weight is 1.

Returns:

self

The instance itself, updated with the estimated location_, base_level_, and structural_index_.

---