diff --git a/doc/api/index.rst b/doc/api/index.rst index 8a6770f96..c3a1c20fe 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -101,8 +101,8 @@ Isostatic Moho isostatic_moho_airy -Source estimation ------------------ +Source position estimation +-------------------------- .. autosummary:: :toctree: generated/ diff --git a/doc/references.rst b/doc/references.rst index 781f2f843..e8f496db9 100644 --- a/doc/references.rst +++ b/doc/references.rst @@ -15,8 +15,11 @@ References .. [Nagy2002] Nagy, D., Papp, G. & Benedek, J.(2002). Corrections to "The gravitational potential and its derivatives for the prism". Journal of Geodesy. doi:`10.1007/s00190-002-0264-7 `__ .. [Oliveira2021] Oliveira Jr, Vanderlei C. and Uieda, Leonardo and Barbosa, Valeria C. F.. Sketch of three coordinate systems: Geocentric Cartesian, Geocentric Geodetic, and Topocentric Cartesian. figshare. doi: `10.6084/m9.figshare.15044241.v1 `__ .. [Reid1990] Reid, A. B., Allsop, J. M., Granser, H., Millett, A. J., & Somerton, I. W. (1990). Magnetic interpretation in three dimensions using Euler deconvolution. GEOPHYSICS. doi:`10.1190/1.1442774 `__ +.. [Reid2014] Reid, A. B., J. Ebbing, and S. J. Webb (2014), Avoidable Euler Errors – the use and abuse of Euler deconvolution applied to potential fields, Geophysical Prospecting, doi:`10.1111/1365-2478.12119 `__. +.. [ReidThurston2014] Reid, A., and J. Thurston (2014), The structural index in gravity and magnetic interpretation: Errors, uses, and abuses, GEOPHYSICS, 79(4), J61-J66, doi:`10.1190/geo2013-0235.1 `__. .. [Soler2019] Soler, S. R., Pesce, A., Gimenez, M. E., & Uieda, L. (2019). Gravitational field calculation in spherical coordinates using variable densities in depth, Geophysical Journal International. doi: `10.1093/gji/ggz277 `__ .. [Soler2021] Soler, S. R. and Uieda, L. (2021). Gradient-boosted equivalent sources, Geophysical Journal International. doi:`10.1093/gji/ggab297 `__ +.. [Uieda2014] Uieda, L., V. C. Oliveira Jr., and V. C. F. Barbosa (2014), Geophysical tutorial: Euler deconvolution of potential-field data, The Leading Edge, 33(4), 448-450, doi:`10.1190/tle33040448.1 `__. .. [Uieda2015] Uieda, Leonardo (2015). A tesserioid (spherical prism) in a geocentric coordinate system with a local-North-oriented coordinate system. figshare. Figure. doi: `10.6084/m9.figshare.1495525.v1 `_ .. [Vajda2004] Vajda, P., Vaníček, P., Novák, P. and Meurers, B. (2004). On evaluation of Newton integrals in geodetic coordinates: Exact formulation and spherical approximation. Contributions to Geophysics and Geodesy, 34(4), 289-314. .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press. diff --git a/harmonica/_euler_deconvolution.py b/harmonica/_euler_deconvolution.py index 92adda96b..ad3d5562a 100644 --- a/harmonica/_euler_deconvolution.py +++ b/harmonica/_euler_deconvolution.py @@ -4,8 +4,12 @@ # # This code is part of the Fatiando a Terra project (https://www.fatiando.org) # +""" +Classes for Euler Deconvolution of potential field data +""" import numpy as np import scipy as sp +import verde.base as vdb class EulerDeconvolution: @@ -18,6 +22,17 @@ class EulerDeconvolution: Euler's homogeneity equation. **Assumes a single data window** and provides a single estimate. + .. hint:: + + Please read the paper [Reid2014]_ to avoid doing **horrible things** + with Euler deconvolution. [Uieda2014]_ offer a practical tutorial using + `legacy Fatiando a Terra `__ code to show + some common misinterpretations. + + .. note:: + + Does not yet support structural index 0. + Parameters ---------- structural_index : int @@ -30,15 +45,52 @@ class EulerDeconvolution: Attributes ---------- - location_ : numpy.ndarray + 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. + + Notes + ----- + + Works on any potential field that satisfies Euler's homogeneity equation + (like gravity, magnetic, and their gradients caused by **simple sources**): + + .. math:: + + (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 :math:`f_i` is the given potential field observation at point + :math:`(e_i, n_i, u_i)`, :math:`b` is the base level (a constant shift of + the field, like a regional field), :math:`\eta` is the structural index, + and :math:`(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 deconvolution estimates :math:`(e_0, n_0, u_0)` and :math:`b` + given a potential field and its easting, northing, and upward derivatives + and the structural index. However, **this 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 + ===================================== ======== ========= - References - ---------- - [Reid1990]_ """ def __init__(self, structural_index): @@ -46,26 +98,32 @@ def __init__(self, structural_index): # The estimated parameters. Start them with None self.location_ = None self.base_level_ = None + self.covariance_ = None - def fit(self, coordinates, field, east_deriv, north_deriv, up_deriv): + def fit(self, coordinates, data): """ 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 east, north, and upward directions. + in easting, northing, and upward directions. + + .. tip:: + + Data does not need to be gridded for this to work. Parameters ---------- - coordinates : tuple of 1d-arrays - Arrays with the coordinates of each data point, in the order of - (x, y, z), representing easting, nothing, and upward directions, - respectively. - field : 1d-array - Field measurements at each data point. - east_deriv, north_deriv, up_deriv : 1d-array - Partial derivatives of the field with respect to east, north, and - upward directions, respectively. + coordinates : tuple 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. + data : tuple 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. Returns ------- @@ -73,19 +131,31 @@ def fit(self, coordinates, field, east_deriv, north_deriv, up_deriv): The instance itself, updated with the estimated `location_` and `base_level_`. """ - n_data = field.shape[0] - matrix = np.empty((n_data, 4)) - matrix[:, 0] = east_deriv - matrix[:, 1] = north_deriv - matrix[:, 2] = up_deriv - matrix[:, 3] = self.structural_index - data_vector = ( - coordinates[0] * east_deriv - + coordinates[1] * north_deriv - + coordinates[2] * up_deriv + coordinates, data, _ = vdb.check_fit_input(coordinates, data, weights=None) + field, east_deriv, north_deriv, up_deriv = vdb.n_1d_arrays(data, 4) + easting, northing, upward = vdb.n_1d_arrays(coordinates, 3) + n_data = field.size + jacobian = np.empty((n_data, 4)) + jacobian[:, 0] = east_deriv + jacobian[:, 1] = north_deriv + jacobian[:, 2] = up_deriv + jacobian[:, 3] = self.structural_index + pseudo_data = ( + easting * east_deriv + + northing * north_deriv + + upward * up_deriv + self.structural_index * field ) - estimate = sp.linalg.solve(matrix.T @ matrix, matrix.T @ data_vector) - + hessian = jacobian.T @ jacobian + # Invert the Hessian instead of solving the system because this is a + # 4x4 system and it won't cost much more, plus we need the inverse + # anyway to estimate the covariance matrix (used as a filtering + # criterion in windowed implementations) + hessian_inv = sp.linalg.inv(hessian) + estimate = hessian_inv @ jacobian.T @ pseudo_data + pseudo_residuals = pseudo_data - jacobian @ estimate + chi_squared = np.sum(pseudo_residuals**2) / (n_data - 4) + self.covariance_ = chi_squared * hessian_inv self.location_ = estimate[:3] self.base_level_ = estimate[3] + return self diff --git a/harmonica/tests/test_euler_deconvolution.py b/harmonica/tests/test_euler_deconvolution.py index 069d39f06..1338a23a2 100644 --- a/harmonica/tests/test_euler_deconvolution.py +++ b/harmonica/tests/test_euler_deconvolution.py @@ -49,10 +49,7 @@ def test_euler_with_numeric_derivatives(): coordinates = (grid_table.easting, grid_table.northing, grid_table.upward) euler.fit( (grid_table.easting, grid_table.northing, grid_table.upward), - grid_table.tfa, - grid_table.d_east, - grid_table.d_north, - grid_table.d_up, + (grid_table.tfa, grid_table.d_east, grid_table.d_north, grid_table.d_up), ) npt.assert_allclose(euler.location_, dipole_coordinates, atol=1.0e-3, rtol=1.0e-3) @@ -84,11 +81,8 @@ def test_euler_with_analytic_derivatives(): euler = EulerDeconvolution(structural_index=2) euler.fit( - (coordinates[0].ravel(), coordinates[1].ravel(), coordinates[2].ravel()), - gz.ravel(), - gze.ravel(), - gzn.ravel(), - gzz.ravel(), + (coordinates[0], coordinates[1], coordinates[2]), + (gz, gze, gzn, gzz), ) npt.assert_allclose(euler.location_, masses_coordinates, atol=1.0e-3, rtol=1.0e-3)