The problem of inverting measured gravity data for large regions is of a great importance for planetary structure studies. Unfortunately, the usual methods of local gravity field inversion do not scale up well. There are three primary factors that start to play significant role: topography or terrain surface with large height differences, spherical geometry of the planet, and high computational complexity. In our previous work we were separately considering each of those problems in detail. In this paper however, we will address those issues simultaneously, offering a complete and computationally effective method of recovering spherical density model of Earth's crust with the upper topography layer. The method utilizes a closed form expression for the discretized model's gravity field which allows for great accuracy and speed without enforcing restrictions on model geometry or gravity field data grid. Inversion process is based on the conjugate gradient method. An example of inversion for a synthetic regional model is presented.
spherical density model, terrain density model, gravity field inversion, gravimetry
Computationally Effective Gravity Inversion Allows for High-Resolution Regional Density Modeling of Earth's Crust with the Inclusion of the Topography Layer
P. S. Martyshko∗,1, D. D. Byzov1, and A. I. Chernoskutov1
1Bulashevich Institute of Geophysics, Ural branch of Russian Academy of Sciences, Yekaterinburg, Russia
Received 20 December 2021; accepted 7 April 2022; published 18 April 2022.
The problem of inverting measured gravity data for large regions is of a great importance for planetary structure studies. Unfortunately, the usual methods of local gravity field inversion do not scale up well. There are three primary factors that start to play significant role: topography or terrain surface with large height differences, spherical geometry of the planet, and high computational complexity. In our previous work we were separately considering each of those problems in detail. In this paper however, we will address those issues simultaneously, offering a complete and computationally effective method of recovering spherical density model of Earth’s crust with the upper topography layer. The method utilizes a closed form expression for the discretized model’s gravity field which allows for great accuracy and speed without enforcing restrictions on model geometry or gravity field data grid. Inversion process is based on the conjugate gradient method. An example of inversion for a synthetic regional model is presented.
Keywords: spherical density model, terrain density model, gravity field inversion, gravimetry
Introduction
Numerical solution of inverse gravimetry problems is one of the tools of interpreting gravimetric data. During preprocessing of the measured field, various methods are used to recalculate the field from the measurement surface (usually, the topography) to a certain reference surface (plane, surface of the Earth ellipsoid). In order to perform this step without loss of accuracy, it would require to have complete data on the density distribution in the region between the measurement surface and the reference surface. In practical applications however, it is common to use speculative or heuristic data on the medium density, since this parameter is unknown and, in fact, it is the final goal of solving the initial problem. To overcome this contradiction and obtain a more accurate result (which is especially important when studying large areas in high resolution), it is necessary to use gravimetric measurements on the topography, without recalculation to an auxiliary reference surface.
Direct calculation of the gravitational field from an object of complex shape, such as a section of the Earth’s crust, which bounded by the surface of the topography above and by the surface of an ellipsoid below, is a computationally difficult task. In , we proposed a computationally efficient method for solving the direct problem, the computational complexity of which does not depend on the geometric shape of the gravitating object (only on its discretization size). The method is a variant finite element method (FEM), based on the approximation of the elements of the partition with polyhedrons. This step is crucial, as it and the subsequent calculation of the field with a closed-form expression, without using the formulas of numerical integration. This leads to high accuracy and computational speed of the method.
In this paper, we will use this method as a basis for constructing a solution to the linear inverse problem of gravimetry. Reconstruction of the density distribution function is carried out for a layered model of the Earth’s crust including the upper topography layer and taking into account Earth’s sphericity. As a model of the initial approximation, we will use previously obtained solution (for a concrete study case), but in a simpler framework in which the whole model volume was represented by a cuboid with flat boundaries (we will refer to that framework as “flat”). It is assumed that this solution already includes all the available priori data, which will enable us to resolve the non-uniqueness of the solution.
Considering similar approaches to solving the inverse problem, it is necessary to note paper , the authors of which describe in detail two examples (synthetic and real) of density model reconstruction from the observed gravitational field, also taking into account the topography. The results of calculations using various weight functions that determine the nature of the separation of features by depth are presented. A non-uniform mesh is used, which increases its step with the depth. That allows obtaining good resolution in the topography layers. It is highlighted, that the gravitational effect of the topography surface significantly affects the measured (or calculated) field. Unfortunately, the paper does not elaborate the method used for solving the direct and inverse problems, optimization techniques, nor provide software/hardware details. The described examples have resolution of less than partition elements, which is not enough for construction of models of regional scale in adequate resolution. Due to the small area of the region under consideration, the authors did not need to take into account the spherical shape of the planet. The grid nodes of the measured field are located on a uniform grid in a plane (the data obtained by the airborne gravimetry), and not on the surface of the topography. In our article we will try to address these issues in more detail.
In , the authors carry out gravity modeling for the structural shells of the moon, including the topography. Calculations are performed in spherical coordinates using the Gauss-Legendre quadrature method to approximate the gravitational integral. High computational complexity of the considered method did not allow for model reconstruction of high detail; however, the work presents estimations of the average structural density in crust, mantle, and core. The computational optimizations in our method allows us to overcome this limitation.
Direct Gravity Calculation for a Spherical Density Model With Topography Layer
Let us define spherical density model in 3D space up to depth with upper topography layer. Models “upper” boundary
(Earth-air interface) is represented by the topography surface (heightmap) relative to the surface S of an ellipsoid of revolution (for example, the WGS84 reference ellipsoid) along its external normal. All points located at a distance of no more than
, along the inner normal to
, and, no more than
, along the outer normal are included in the model, where
is the elevation of
over
,
is latitude,
is longitude associated with the ellipsoid. In the described region
, the density distribution is defined as
, where
.
The vertical component of the gravitational field gradient induced by the region
at the outer point
is defined by the integral:
where
is the gravitational constant,
is the element of the volume of integration,
is the outward normal to
in the orthogonal projection of the point
onto
,
and
are the radius vectors of the points
and
, respectively.
Let us choose some discretization of . Now we can approximate the volume of an element of the partition
by a polyhedron
, obtaining an approximated model, to which we can apply the previously proposed algorithm for solving the direct gravity problem. Here we will reproduce the basic steps of the algorithm.
Let the density of each be uniform and equal to
. The field
for the model
at the point
is calculated through the sum of the field of each discretization element:
where
is the field of the element
with unit density at the point
up to the coefficient
.
It is apparent that the integral [eq:e1] for the field cannot be expressed in a closed form. It is also problematic to calculate the integral numerically with cubature formulas, since the boundaries of
can have complex description. That would require first or second order formulas with a large number of nodes in order to achieve acceptable accuracy. Therefore, we will calculate the integral [eq:e1] not for
, but for the approximating polyhedron
. We denote the set of faces
as
, then,
Next, we convert the volume integral to the surface integral, applying the divergence theorem to [eq:e2]:
The surface integral can be split to a sum of integrals for each of the faces of . Note, that external normal
is constant in every point of integration for a face:
where
is external normal to the face
.
As a final step we need to find a closed form expression for the integral over a triangle. In fact, this integral is the gravity potential of the triangle plate with unit density up to
. Let’s denote
are radius vectors of the triangle vertices
is the point of field calculation;
;
;
is the normal to the triangle with the length equal to its double surface area;
is the unit normal to the triangle;
;
is the (signed) distance from the point
to the triangle surface. Then,
Thus, we’ve devised an algorithm for calculating of the described model without using the formulas of approximate integration. The accuracy of the method depends only on the quality of the approximation of the elements
by the polyhedrons
.
In we’ve presented an example of calculating the gravitational field on the topography surface for a practical model, taking into account the Earth’s sphericity. The introduction of the complex-shaped surfaces (such as the topography) to the model does not affect the performance of calculations. This is one of the defining qualities of the method.
Inverse Linear Problem for the Gravity Field
Having in our disposal a computationally efficient method for solving the direct gravity problem for the described model, we can now proceed to consider the inverse problem. Finding solution to the inverse problem in an iterative process during which multiple computation of the direct problem will be required. This fact necessitates a high-performance computation method for [eq:e1].
Here we will consider a three-dimensional linear inverse problem, which consists in restoring the density distribution function in the form of
for the volume of the model under consideration from the input field
measured (or calculated) on the topography. Due to the fact that the problem does not have a unique solution and is unstable, it becomes necessary to impose additional constrains on the solution, as well as apply regularization techniques. When considering inversion for practical models, it is important to ensure that the set of feasible solutions remains geologically meaningful. For this purpose, various methods have been proposed that allow to incorporate additional a priori data in the solution, e.g. .
It is our opinion, during applying these steps, it is reasonable to introduce a “flat” analogue for the target density model, when the partitioning elements of the model are represented as right rectangular prisms. The gravity field of the simplified model should defer from the original no more than of order of 5–10%. In this setting, the direct and, as a consequence, the inverse problem has a much more computationally efficient way of acquiring numerical solution . This circumstance is extremely important, since, in order to obtain the desired, geologically justified morphology, finding multiple alternative solutions of the inverse problem may prove to be necessary. When a satisfactory solution in the “flat” framework is found, it can be used as an initial approximation for the original inverse problem (taking into account the topography and sphericity). This step can be performed with transformation of the “flat” geometry into spherical . Thus, we will consider the problem of refinement of the already existing model, which includes all the necessary a priori data, to fit the target field.
Let us denote the set of all from the domain
as the vector
, and introduce discretization for
and denote the set of its values at the grid nodes by the vector
. Then, the direct gravity problem can be written as
, where
is the discretized operator of the direct problem.
We will seek a solution to the inverse problem by minimizing the functional that takes into account the field residual and the deviation from the initial approximation model (the least-squares problem with regularization): That minimization problem can be reduced to a system of linear equations:
The system can be efficiently (in terms of computational resources) solved using the conjugate gradient method. The iterative process is described as follows :
Here, is iteration count,
,
. The initial value of
is set to zero, since the initial problem targets for minimum deviation from the initial approximation (which field is subtracted from the measured field, resulting in
). Stop condition is
for two consecutive iterations.
Gravity Inversion: a Synthetic Example
Initial Data Generation
To construct an example of solving the inverse problem, we will calculate the “measured” field by solving the direct problem for a synthetic density model. For the area with geographical coordinates 60–68
N, 48
–72
E the model has a layered structure along the inner normal to the surface of the Earth’s ellipsoid WGS84, the layer thickness of 1 km (along the normal) and the lateral discretization is
(
), the number of layers is 81. The inner space of the model has density values (from top to bottom) 2.5, 2.9 and 3.3 g/cm
with two horizontal curvilinear boundaries (sedimentary cover and Moho discontinuity) with asymptotes of
km and
km.
A synthetic example is constructed: for a field on topography from a model of a three-layer medium (two boundaries with a density jump), a linear inverse problem is solved (determining the density values in a selected volume). The purpose of constructing this example is to demonstrate the convergence of the method and performance for large amounts of input data.
The topography model is built on the top surface of the Earth ellipsoid using the heightmap of the designated area. Each point of the heightmap is represented with a vertical mass column of the corresponding height. The column density is set in proportion to its height in the range 1.8–3.2 g/cm, where the first value corresponds to zero height, the second – to the maximum height (1.4 km for the given area). This distribution was chosen rather arbitrarily, but, as noted in , density variations in the topography layer do not significantly affect the resulting gravitational field; the main contribution is made by the shape of the air-mass interface itself. When calculating the field of the model, the average density of each horizontal and the topography layer is subtracted from that layer density distribution. The calculation of the field is carried out at the points on the heightmap (the topography). shows the calculated field, the topography and layered medium of the model (excess density is shown in each 1 km layer). All images are pictured in the transverse Mercator (Gauss-Kruger) projection.
For the described numerical experiments, the heightmap of the topography and geoid from the ICGEM online resource were used.
Recovering Density Distribution
In order to construct a solution using the described method, it is necessary to determine a set of coefficients providing additional information on the distribution of the depth features. For this example, the coefficients were selected subjectively through a series of experiments, with a goal of minimizing the maximum deviation from the model of the initial approximation (zero), while maintaining the presence of anomalous masses in all layers of the model. The final expression is a linear mapping
(z is in km). The result of the solution is shown in . Density values of the topography layer are in the range [
] (5-percentile
, 95-percentile
); in the first top layer under the topography the range is [
; 0.24] (5-percentile
, 95-percentile
). For the subsequent layers these parameters rapidly decrease with depth and do not exceed in absolute value those that given for the upper layers. The final relative error in the field was 2.6%, the value of the stop condition has reached 0.006.
It is necessary to note that in the example described, restoring the surface of the boundaries shown in was not one of the goals. Indeed, when solving the inverse problem in the described form, it is impossible to unambiguously recover vertical density discontinuities; for this purpose, one should consider a structural inverse problem . In our case, the vertical discontinuities will be roughly reflected in the initial approximation model (in which the density of each layer is constant and equal to the average value of the layer density of the original model). shows the result of solving the “refinement” problem of the initial approximation model.
Conclusion
The proposed method for solving the linear inverse gravity problem makes it possible to consider models of complex geometry, bounded by the topography surface, while taking into account sphericity of the planet. The algorithm developed by the authors makes it possible to calculate the magnitude of the gravitational field on the topography surface without the impact on the computational efficiency: the method does not rely on the grid regularity of the density model and the calculated field. These circumstances allow the use of gravity measurements taken on the surface of the topography, directly, avoiding significant errors associated with the recalculation of the field to the reference surface. The recalculation step can be completely excluded from the interpretation process.
When calculating the direct gravity problem (which is needed in any iterative process of solving the inverse problem), we’ve developed software based on the program GRAFEN (also developed by the authors). The computation time for the example model was 2 min. per iteration of the inverse problem algorithm using five AMD Radeon VII GPUs.
Acknowledgments.
The research was supported by the Russian Fund for Basic Researches (project 20-05-00230 A).
-
Corresponding author: pmart3@mail.ru↩︎
1. Chernoskutov, A. I., and D. D. Byzov, GRAFEN v0.1 - Gravity Field Ellipsoidal Density Model Numerical Computations for CUDA-Enabled Distributed Systems, https://github.com/charlespwd/project-title, 2019.
2. Ince, E. S., F. Barthelmes, S. Reißland, K. Elger, C. Forste, F. Flechtner, and H. Schuh, ICGEM - 15 Years of Successful Collection and Distribution of Global Gravitational Models, Associated Services, and Future Plans, Earth System Science Data, 11(2), 647-674, doihttps://doi.org/10.5194/essd-11-647-2019, 2019.
3. Ladovskii, I. V., P. S. Martyshko, D. D. Byzov, and V. V. Kolmogorova, On Selecting the Excess Density in Gravity Modeling of Inhomogeneous Media, Izvestiya, Physics of the Solid Earth, 53(1), 130-139, doihttps://doi.org/10.1134/S1069351316060057, 2017.
4. Martyshko, P. S., I. V. Ladovskii, and A. G. Tsidaev, Construction of Regional Geophysical Models Based on the Joint Interpretation of Gravitaty and Seismic Data, Izvestiya, Physics of the Solid Earth, 46(11), 931-942, doihttps://doi.org/10.1134/S1069351310110030, 2010.
5. Martyshko, P. S., I. V. Ladovskii, D. D. Byzov, and A. G. Tsidaev, Gravity Data Inversion with Method of Local Corrections for Finite Elements Models, Geosciences, 8(10), doihttps://doi.org/10.3390/geosciences8100373, 2018a.
6. Martyshko, P. S., I. V. Ladovskij, D. D. Byzov, and A. I. Chernoskutov, On Solving the Forward Problem of Gravimetry in Curvilinear and Cartesian Coordinates: Krasovskii’s Ellipsoid and Plane Modeling, Izvestiya, Physics of the Solid Earth, 54(4), 565-573, doihttps://doi.org/10.1134/S1069351318040079, 2018b.
7. Martyshko, P. S., D. D. Byzov, and A. I. Chernoskutov, Interpretation of Gravity Data Measured by Topography, Doklady Earth Sciences, 495(2), 914-917, doihttps://doi.org/10.1134/S1028334X20120077, 2020.
8. Pedersen, L. B., J. Kamm, and M. Bastani, A Priori Models and Inversion of Gravity Gradient Data in Hilly Terrain, Geophysical Prospecting, 68(3), 1072-1085, doihttps://doi.org/10.1111/1365-2478.12897, 2020.
9. Potts, L. V., and R. R. von Frese, Comprehensive Mass Modeling of the Moon From Spectrally Correlated Free-Air and Terrain Gravity Data, Journal of Geophysical Research: Planets, 108(E4), doi:https://doi.org/10.1029/2000JE001440, 2003.
10. van der Vorst, H. A., Iterative Krylov Methods for Large Linear Systems, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, doihttps://doi.org/10.1017/CBO9780511615115, 2003.
