Parallel Algorithms for Solving Inverse Gravimetry Problems: Application for Earth’s Crust Density Models Creation

Abstract

The paper describes a method of gravity data inversion, which is based on parallel algorithms. The choice of the density model of the initial approximation and the set on which the solution is sought guarantees the stability of the algorithms. We offer a new upward and downward continuation algorithm for separating the effects of shallow and deep sources. Using separated field of layers, the density distribution is restored in a form of 3D grid. We use the iterative parallel algorithms for the downward continuation and restoration of the density values (by solving the inverse linear gravity problem). The algorithms are based on the ideas of local minimization; they do not require a nonlinear minimization; they are easier to implement and have better stability. We also suggest an optimization of the gravity field calculation, which speeds up the inversion. A practical example of interpretation is presented for the gravity data of the Urals region, Russia.

1. Introduction

Gravity data inversion is the main tool for obtaining the Earth’s crust density model. Our main goal of the study is to construct a stable high-resolution method of three-dimensional (3D) gravity data inversion. The simplest approach of finding 3D density values is by forward gravity modeling. We can vary the initial model in an interactive mode to reduce the gravity field residuals. However, such an approach cannot resolve the non-uniqueness and instability of inverse problem solving. Our method of inversion follows the procedure of Cordell and Henderson [1]. The density values are calculated using the observed gravity field by iterative regularized algorithm, without the time-consuming trial-and-error procedure.

The interpretation process involves solving the forward and inverse gravity problems. The inverse gravity problems are ill-posed since their solution is non-unique and unstable depending on the initial data [2], hence why one should define the set of correctness for the solution and select a reasonable initial condition. Density is one of the most important parameters for geophysical modelling. Density models are usually created using complete information (for example, seismic data) to select the initial model. We apply the iterative scheme for transforming the model slightly at each step to create the geological meaningful density model [3]. To extend the layer outside the study volume, we suggest selecting piecewise-constant density function ${\rho }_0\left(z\right)$, which depends only on depth.

We propose a method for constructing density models with flat topography based on inversion of gravity data with a complete Bouguer correction. In this article we describe the creation of a density model method for layered medium. The process of constructing density models based on gravity field anomalies proposed here contains five major steps: the fast forward modeling algorithm; the initial 3D density model creation; the upward and downward continuation (the extraction of the gravity fields of layers); the choice of the density as multiplicative function; and the stable adaptive algorithm for the inverse problem solving. The horizontal layered models may well correspond to reality in some areas but deviate from the actual geology in others; this is a limitation of the method.

For layer-by-layer separation of the gravitational field by depth, we use a filter based on the analytical continuation of the vertical field component upward and downward with Lavrentiev’s regularization. Several researchers have been engaged in the downward continuation of potential fields. In [4], the authors apply Tikhonov’s regularization to stabilize the process of solving this problem. They present realization of the method in the form of a Matlab-based program. The optimum regularization parameter value is selected as a local minimum of constructed L_p-norms functions-in the majority of cases. They demonstrate excellent stabilizing properties of this method on several synthetic models and one real-world example from high-definition magnetometry.

In [5], the authors propose using the combination of upward continuation and horizontal derivative to accomplish the downward continuation of potential field data. The proposed method was demonstrated on synthetic potential field data. The method was also applied to real potential field data, and the results show that the proposed method accomplishes the downward continuation of the real data, stably.

An alternative method of downward continuation is proposed in [6]. The authors use a combination of Taylor series expansion and upward continuation for computing vertical derivatives. This method has been tested on the gravitational anomaly of infinite horizontal cylinder in both cases with and without random noise. The vertical derivative method is applied to calculate the downward continuation according to the Taylor series expansion method. The method was also tested on both complex synthetic models and real data.

In [7], the authors formulate a functional model for a spectral downward continuation of selected gravitational field quantities to an irregular topographic surface. This functional model is further generalized to allow for transformation between different types of gravitational field quantities. In particular, spectral weights are derived for estimation of the disturbing potential or disturbing/anomalous gravity at the Earth’s surface by combining the first-, second- and third-order radial gradients of the disturbing potential (disturbing gradients). The combined spectral estimator is applied to simulated satellite disturbing gradients polluted by a realistic Gaussian noise. The spectral combination method requires no matrix inversion.

Another example of applying regularization to inversion of potential field data is presented in [8]. Authors have introduced a new processing and modeling technique of magnetic data conducted over mines or near-surface geophysical targets for accurate and precise determination of location and depth. The technique is based on the application of the Kaczmarz regularization method to the ill-posed magnetic inverse problem. The advantage of this method is the optimum transformation of regularized normal equations to an equivalent augmented regularized normal system of equations. The method is applied to an unexploded ordnance test site in the United Kingdom.

A fast algorithm for a forward gravity problem is absolutely essential for efficiently solving the inverse problem. We build the initial density model using deep seismic sounding (DSS) data along two-dimensional seismic profiles. The volume between profiles is filled with interpolated values of density on a regular 3D grid [3]. Each cell of the grid is a rectangular prism with constant density. The analytical form for the gravity field of the rectangular prism is well-known [9,10]. The finite element method (FEM) in forward modelling is used in this study. Mostafa [11] and Couder-Castañeda et al. [12] are relevant studies that consider relatively small grid sizes: the models consist of 10,000–15,000 prisms in each example. For our study, we required a technique that could handle more than a billion elements in a reasonable time. Therefore, we cannot rely on parallelization alone and it is necessary to modify the algorithm itself. We note that Dubey and Tiwari [13], Li and Oldenburg [14] consider elements of non-equal size, and affirm that FEM can also be applied in iterative field inversion.

The correctness of the 3D inversion depends on the technique of mass continuation outside the study area. One cannot simply consider that density is zero everywhere outside. In this paper, the authors use background density, calculated by averaging the density values of every horizontal layer of the interpolated model. This approach is similar to that presented by Cai and Wang [15].

The paper is organized as follows: we present new mathematical algorithms for downward continuation of the gravity field and the 3D inverse problem, and then both algorithms are applied to the gravity anomalies of the Middle Urals region, Russia.

2. Methods

2.1. Layer-by-Layer Separation of the Local Field Component

In [16], a technique was proposed for separating the sources of the gravitational field by depth. Examples of the implementation of the technique are given in and [17]. In this paper, we describe an implementation of this technique based on a new and more efficient numerical algorithm that enables parallel computing on graphics accelerators.

In order to separate the observed field anomalies by depth, we use the original method of elevated transformations [16,18] with the modification of the computational part, proposed by the authors of this paper, based on the method of local corrections [17,19,20]. We use this method for downward continuation and the density refinement. We assume that the original field $g\left(x,y,0\right)$ is specified on the plane z = 0.

The general scheme of the method for isolating the effects of the sources in a layer down to the certain depth (concerning the Earth’s surface) z = −H consists of three stages.

1.

We construct the analytical continuation of the field upward to the height z = H: $g\left(x,y,0\right)\stackrel{up\left(H\right)}{\Rightarrow }g\left(x,y,H\right)$, while assuming that the influence of the local near-surface sources (down to the depth z = −H) is significantly weakened or eliminated at all.

2.

In order to hide the influence of the local sources located in the horizontal layer between the plain z = 0 and the plain z = −H, the field $g\left(x,y,H\right)$ is recalculated upwards is then analytically continued downward to the depth z = −H: $g\left(x,y,H\right)\stackrel{down\left(2H,\kappa \right)}{\Rightarrow }u\left(x,y,-H|(-\infty ,-H]\right)$. Moreover, since the problem is ill-posed, it is necessary to use methods with regularization, hence κ in the formula is a regularization parameter. The resulting field is calculated on the plane z = −H from sources located below the boundary z = −H.

3.

At the last step, we recalculate the field $u\left(x,y,-H|(-\infty ,-H]\right)$ upward again to the level of the plain z = 0: $u\left(x,y,-H|(-\infty ,-H]\right)\stackrel{up\left(H\right)}{\Rightarrow }u\left(x,y,0|(-\infty ,-H]\right)$. The resulting field is calculated on the plane z = −H from sources located below the boundary z = 0. Further, we subtract this field from the original one obtaining the field for the layer $z\in (-H,0]$: $u\left(x,y,0|(-H,0]\right)=g\left(x,y,0\right)-u\left(x,y,0|(-\infty ,-H]\right)$.

If we want to obtain the field $u\left(x,y,0|(-H_2,-H_1]\right)$ on the plane z = 0 from the sources located in the depth interval $z\in (-H_2,-H_1]$, then we must perform three stages of the method for the heights H₁ and H₂, and then take the difference between the results: $u\left(x,y,0|(-H_2,-H_1]\right)=u\left(x,y,0|(-\infty ,-H_1]\right)-u\left(x,y,0|(-\infty ,-H_2]\right)$.

Thus, if the observed field g (x, y, 0) is given on the plane z = 0, and there is also an ordered set of depths ${\left(-H_i\right)}_{i=0}^L$, $0=H_0<H_1<\cdots <H_L$, then the problem of dividing g(x, y, 0) into components ${\left(u\left(x,y,0|(-H_{i+1},-H_i]\right)\right)}_{i=0}^{L-1}$ on the plane z = 0, corresponding to fields from sources located in depth intervals ${\left((-H_{i+1},-H_i]\right)}_{i=0}^{L-1}$, can be solved as follows. For every depth $-H_i$ independently, the three steps of the recalculation method described above are performed, and the field $u\left(x,y,0|(-\infty ,-H_i]\right)$ is obtained. In this case, at the second stage, the regularization parameter ${\kappa }_i$ is used, $u\left(x,y,0|(-\infty ,-H_0]\right)=g\left(x,y,0\right)$ Further, the difference for the successive depths is calculated $u\left(x,y,0|(-H_{i+1},-H_i]\right)=u\left(x,y,0|(-\infty ,-H_i]\right)-u\left(x,y,0|(-\infty ,-H_{i+1}]\right)$. It turns out that $g\left(x,y,0\right)=u\left(x,y,0|(-\infty ,-H_L]\right)+\sum _{i=0}^{L-1}u\left(x,y,0|(-H_{i+1},-H_i]\right)$.

2.2. Definition of the Operation Up(·)

Suppose the gravitating masses are located in the layer below the horizontal plane at the depth z. On this plane, we denote the gravitational field by $u\left(·,·,z\right)$ and consider it as the boundary function for the Dirichlet problem for the Laplace equation over a semi-infinite domain. From the values $u\left(·,·,\zeta \right)$ on the boundary, the solution $u\left(·,·,\zeta \right)$ reconstructs the harmonic field function everywhere above z. So, for the upper half-space $\zeta \ge z$ the solution $u\left(·,·,\zeta \right)$ of the problem can be written in terms of the Poisson integral [21]:

$$u\left(\xi ,\eta ,\zeta \right)=\frac{\zeta -z}{2\pi }\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}\frac{u\left(x,y,z\right)dxdy}{{\left({\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+{\left(z-\zeta \right)}^2\right)}^{\frac{3}{2}}}.$$

(1)

The operation of the upward recalculation $u\left(·,·,z\right)\stackrel{up\left(H\right)}{\Rightarrow }u\left(·,·,\zeta \right)$ is a direct calculation of the integral in the Formula (1). The operation of the downward recalculation $u\left(·,·,\zeta \right)\stackrel{down\left(H,\kappa \right)}{\Rightarrow }u\left(·,·,z\right)$ is the solution of the Fredholm integral equation of the first kind (1), where the values of the field $u\left(·,·,\zeta \right)$ are assumed to be given, and the values of$u\left(·,·,z\right)$ under the integral sign are to solve for. In the general case, finding a solution to this equation is an ill-posed problem requiring regularization, which is, again, reflected in the notation: κ is the regularization parameter.

The implementation of the $up\left(H\right)$ operation for the upward recalculation that accounts for the difference in heights $=\zeta -z\ge 0;$ $up\left(H\right)$, described in detail here, will be used later in the paper. We assume that the field u(x,y,z) is defined everywhere on the z plane as a piecewise constant function:

$$u\left(x,y,z\right)=\{\begin{array}{c}{u}_{i,j}\left(z\right),\left(x,y\right)\in [x_i;x_{i+1})\times [y_j;y_{j+1}) \\ u_a,\left(x,y\right)\notin [x_{min};x_{max})\times [y_{min};y_{max})\end{array},$$

(2)

where $x_{min}=x_0<x_1<\cdots <x_{M_x}=x_{max}$, $y_{min}=y_0<y_1<\cdots <y_{M_y}=y_{max}$, $u_a$ can be considered as the horizontal asymptote of the field. In this case, the integral (1) is obtained analytically in the following form:

$$\begin{array}{c}u\left(\xi ,\eta ,\zeta \right)=u_a+\hfill \\ \hspace{1em}\hspace{1em}+\frac{1}{2\pi }{\displaystyle \sum _{i=0}^{M_x-1}}{\displaystyle \sum _{j=0}^{M_y-1}}\left(u_{i,j}\left(z\right)-u_a\right)\left(\mathrm{arctg}\left(\frac{\left(x-\xi \right)\left(y-\eta \right)}{HR}\right)\right)|\begin{array}{c}{x}_{i+1}\\ \\ x=x_i\end{array}|\begin{array}{c}{y}_{j+1}\\ \\ y=y_j\end{array}, \hfill \end{array}$$

(3)

where $R=\sqrt{{\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+H^2}$.

The field asymptote does not change during recalculation. Formula (3) allows for calculating the values $u\left(\xi ,\eta ,\zeta \right)$ of the field recalculated to the height ζ at the point. For the downward recalculation $down\left(·,·\right)$ we must approximate the upward recalculated field $u\left(\xi ,\eta ,\zeta \right)$ by a piecewise constant function at the same intervals as the field $u\left(x,y,z\right)$. To do this, we need to average $u\left(\xi ,\eta ,\zeta \right)$ over the said intervals:

$$u\left(\xi ,\eta ,\zeta \right)\approx \{\begin{array}{c}{u}_{n,m}\left(\zeta \right),\left(\xi ,\eta \right)\in [x_n;x_{n+1})\times [y_m;y_{m+1})\\ u_a,\left(\xi ,\eta \right)\notin [x_{min};x_{max})\times [y_{min};y_{max}) \end{array},$$

(4)

$$\begin{gathered} u_{n,m}\left(\zeta \right)=\frac{1}{\left(x_{n+1}-x_n\right)\left(y_{m+1}-y_m\right)}\underset{x_n}{\overset{x_{n+1}}{\int }}\underset{y_m}{\overset{y_{m+1}}{\int }}u\left(\xi ,\eta ,\zeta \right)d\eta d\xi = \\ =u_a+\frac{1}{2\pi \left(x_{n+1}-x_n\right)\left(y_{m+1}-y_m\right)}\times \\ \times \sum _{i=0}^{M_x-1}\sum _{j=0}^{M_y-1}\left(u_{i,j}\left(z\right)-u_a\right)\nu \left(x-\xi ,y-\eta ,H\right)|\begin{array}{c}{x}_{n+1}\\ \\ \xi =x_n\end{array}|\begin{array}{c}{y}_{m+1}\\ \\ \eta =y_m\end{array}|\begin{array}{c}{x}_{i+1}\\ \\ x=x_i\end{array}|\begin{array}{c}{y}_{j+1}\\ \\ y=y_j\end{array}, \end{gathered}$$

(5)

where $\nu \left(x,y,z\right)=xy·\mathrm{arctg}\left(\frac{xy}{zR}\right)-\frac{xz}{2}\ln \left(\frac{R-x}{R+x}\right)-\frac{yz}{2}\ln \left(\frac{R-y}{R+y}\right)-zR$, $R=\sqrt{x^2+y^2+z^2}$.

Note that in the sequential calculation $u_{n,m}\left(\zeta \right)$ for all possible indices $0\le n\le M_x-1$ and $0\le m\le M_y-1$ according to the Formula (5), we would have to compute the value of the ν function $4M_x^2{M}_y^2$ times. However, it is possible to reduce this number if we assume that the field sampling grid is uniform in x and y, then $\forall 0\le i\le M_x-1, 0\le j\le M_y-1: x_{i+1}-x_i=\Delta x, y_{j+1}-y_j=\Delta y$ and we can rewrite (5) as:

$$u_{n,m}\left(\zeta \right)=u_a+\frac{1}{2\pi \Delta x\Delta y}\left(T_{n+1,m+1}-T_{n+1,m}-T_{n, m+1}+T_{n,m}\right),$$

(6)

$$T_{n,m}={\displaystyle \sum }_{i=0}^{M_x-1}{\displaystyle \sum }_{j=0}^{M_y-1}\left(u_{i,j}\left(z\right)-u_a\right)\nu \left(x-x_n,y-y_m,H\right)|\begin{array}{c}{x}_{i+1}\\ \\ x=x_i\end{array}|\begin{array}{c}{y}_{j+1}\\ \\ y=y_j\end{array}.$$

(7)

Applying the same optimization technique for calculating the convolution of functions on a uniform grid, which was used in [20], for a fast algorithm for solving the forward gravity problem, we can write (7) in the form:

$$T_{n,m}={\displaystyle \sum }_{i=-n}^{M_x-n}{\displaystyle \sum }_{j=-m}^{M_y-m}{\overline{u}}_{i+n,j+m}\left(z\right){\nu }_{i,j}\left(H\right),$$

(8)

where ${\overline{u}}_{n,m}\left(z\right)=u_{n,m}\left(z\right)-u_{n,m-1}\left(z\right)-u_{n-1,m}\left(z\right)+u_{n-1,m-1}\left(z\right)$, suppose that $u_{n,m}\left(z\right)=u_a$, if $n=-1\vee n=M_x\vee m=-1\vee m=M_y$; ${\nu }_{i,j}\left(H\right)=\nu \left(i\Delta x,j\Delta y,H\right)$.

Thus, $up\left(H\right)$ recalculation consists of following steps.

• The preprocessing stage. This stage is done once and is suitable for any initial field, depends only on the difference in the heights H of the recalculation and the grid discretization steps Δx, Δy. It is based on the fact that the elements ${\nu }_{i,j}\left(H\right)$ do not depend on the input amplitudes $u_{n,m}\left(z\right)$ of the field, only on H, Δx, and Δy. The stage consists of calculating $\left(2M_x+1\right)\left(2M_y+1\right)$ elements of the set ${\left({\nu }_{i,j}\left(H\right)\right)}_{\left(i,j\right)=\left(-M_x,-M_y\right)}^{\left(M_x,M_y\right)}$. Note that the function $\nu \left(x,y,z\right)$ is even for x and y and odd for z, hence ${\nu }_{i,j}\left(H\right)={\nu }_{-i,j}\left(H\right)={\nu }_{i,-j}\left(H\right)={\nu }_{-i,-j}\left(H\right)$, so in fact, one only needs to calculate $\left(M_x+1\right)\left(M_y+1\right)$ values of $\nu \left(x,y,z\right)$ for different arguments, which is an order less than $4M_x^2{M}_y^2$ values when calculating “directly” as per (5).
• Setup initial data: set ${\left(u_{n,m}\left(z\right)\right)}_{\left(n,m\right)=\left(0,0\right)}^{\left(M_x-1,M_y-1\right)}$ from $M_x{M}_y$ elements (field amplitudes in piecewise constant representation).

1.

Calculation of $\left(M_x+1\right)\left(M_y+1\right)$ elements of the set ${\left({\overline{u}}_{n,m}\left(z\right)\right)}_{\left(n,m\right)=\left(0,0\right)}^{\left(M_x,M_y\right)}$.

2.

Calculation of $\left(M_x+1\right)\left(M_y+1\right)$ elements of the set ${\left(T_{n,m}\right)}_{\left(n,m\right)=\left(0,0\right)}^{\left(M_x,M_y\right)}$ (8).

3.

Calculation of $M_x{M}_y$ elements of the output set ${\left(u_{n,m}\left(\zeta \right)\right)}_{\left(n,m\right)=\left(0,0\right)}^{\left(M_x-1,M_y-1\right)}$ (6).

As one can see, the preprocessing step, which is lengthy in terms of the computation time, is not directly included in the recalculation up(H) at all (steps 1–3). This fact makes it possible to significantly reduce the calculation time of the downward recalculation operation down(H,$\kappa$), which is defined later.

Similar to (8), one can write the “fast” version of (3) for the upward recalculation for a set of points of the form $\left({\xi }_0+n\Delta x,{\eta }_0+m\Delta y,\zeta \right)$:

$$\begin{gathered} u\left({\xi }_0+n\Delta x,{\eta }_0+m\Delta y,\zeta \right)=u_a+ \\ +\frac{1}{2\pi }{\displaystyle \sum }_{i=-n}^{M_x-n}{\displaystyle \sum }_{j=-m}^{M_y-m}{\overline{u}}_{i+n,j+m}\left(z\right)\mathrm{arctg}\left(\frac{\left(i\Delta x-{\xi }_0\right)\left(j\Delta y-{\eta }_0\right)}{HR}\right), \end{gathered}$$

(9)

where $R=\sqrt{{\left(i\Delta x-{\xi }_0\right)}^2+{\left(j\Delta y-{\eta }_0\right)}^2+H^2}$. Moreover, if ${\xi }_0$ and ${\eta }_0$ are fixed, $0\le n\le N_x-1$,$0\le m\le N_y-1$, then the value of the arctangent must be computed for $\left(N_x+M_x\right)\left(N_y+M_y\right)$ different arguments. If, under the same conditions, we use Formula (3), then we would need $4N_x{N}_y{M}_x{M}_y$ of arctangent computations. As it becomes evident, among these $4N_x{N}_y{M}_x{M}_y$ values there would be no more than $\left(N_x+M_x\right)\left(N_y+M_y\right)$ unique ones, which was not completely evident from the expression (3).

2.3. Description of the Operation Down(·,·)

The implementation of the down(H, κ) operation for downward recalculation for the height difference $H=\zeta -z\ge 0$, with a regularization parameter κ, is described in detail hereAs mentioned above, it is necessary to solve the Fredholm integral equation of the first kind (1) for the unknown function $u\left(x,y,z\right)$, where z is fixed.

We will use the same piecewise constant function representation of the field as its parameterization as for the operation up(H), we seek $u\left(x,y,z\right)$ in the form of (2) and assume that the field $u\left(\xi ,\eta ,\zeta \right)$ on the left side of (1) is given by its mean values (4) within the indicated boundaries. We seek $u\left(x,y,z\right)$ using the modified method of local corrections, which was described as a solution of linear inverse gravimetry problem in [20]. Its stages for solving Equation (1) with the chosen parameterization are reproduced as follows. The initial data at the iteration “zero” (the iteration number is denoted by the superscript in parentheses): the set ${\left(\delta u_{n,m}^{\left(0\right)}\left(\zeta \right)\right)}_{\left(n,m\right)=\left(0,0\right)}^{\left(M_x-1,M_y-1\right)}$ of the “residual” field amplitudes at a height ζ (we take $\delta u_{n,m}^{\left(0\right)}\left(\zeta \right)=u_{n,m}\left(\zeta \right)$); the set ${\left(u_{n,m}^{\left(0\right)}\left(z\right)\right)}_{\left(n,m\right)=\left(0,0\right)}^{\left(M_x-1,M_y-1\right)}$ of the field amplitudes at height z (we take $u_{n,m}^{\left(0\right)}\left(z\right)=0$). The definition of the “residual” field is $\delta u_{n,m}^{\left(\theta \right)}\left(\zeta \right)$ at the iteration θ: assuming that $u^{\left(\theta \right)}\left(z\right)\stackrel{up\left(H\right)}{\Rightarrow }{u}^{\left(\theta \right)}\left(\zeta \right)$, then $\delta u_{n,m}^{\left(\theta \right)}\left(\zeta \right)=u_{n,m}\left(\zeta \right)-u_{n,m}^{\left(\theta \right)}\left(\zeta \right)$. At each iteration the horizontal asymptote for all fields appearing in the algorithm remains equal $u_a$. Notice that it is not the absolute values of the fields that are important, but their deviations from the asymptote. To explain the used notation: when we write $u_{n,m}^{\left(\theta \right)}\left(z\right)$, $u_{n,m}^{\left(\theta \right)}\left(\zeta \right)$, $\delta u_{n,m}^{\left(\theta \right)}\left(\zeta \right)$, $w_{n,m}^{\left(\theta \right)}$, $S_{n,m}$ we mean the average value of these functions in the region $[x_n;x_{n+1})\times [y_m;y_{m+1})$, when we write $u^{\left(\theta \right)}\left(z\right)$, $u^{\left(\theta \right)}\left(\zeta \right)$, $\delta u^{\left(\theta \right)}\left(\zeta \right)$, $w^{\left(\theta \right)}$, $S$ (without indices) we mean piecewise constant functions depending on x and y in the domain ${ℝ}^2$. For the method below, it is essential that the region in x and y, on which we will vary the required function $u\left(x,y,z\right)$, is $[x_{min};x_{max})\times [y_{min};y_{max})$ and coincides with the region of the actual definition of the known field $u\left(x,y,\zeta \right)$ at the height ζ. Moreover, the values of all functions $u^{\left(\theta \right)}\left(z\right)$, $u^{\left(\theta \right)}\left(\zeta \right)$, $\delta u^{\left(\theta \right)}\left(\zeta \right)$, $w^{\left(\theta \right)}$, $S$ appearing in the method that are outside of the region $[x_{min};x_{max})\times [y_{min};y_{max})$ are set to the value of the horizontal asymptote $u_a$. For example, the results of up(H) recalculation outside the specified region are replaced with $u_a$.

So, the stages of the modified method of local corrections for solving (1) with the chosen piecewise constant parameterization are as follows:

1.

Calculate operation up(·) for $\delta u^{\left(\theta -1\right)}\left(\zeta \right)\stackrel{up\left(H\right)}{\Rightarrow }{w}^{\left(\theta \right)}$. On a side note, according to the physical meaning of the up(H) recalculation, the field $w^{\left(\theta \right)}$ is defined in the plane with the applicate ζ + H. But in what follows, we only need $w^{\left(\theta \right)}$ as a formal function of x and y, its physical meaning is of no matter.

2.

Calculate ${\alpha }^{\left(\theta \right)}$ and ${\beta }^{\left(\theta \right)}$:

$$\begin{array}{c}{\alpha }^{\left(\theta \right)}=\frac{1}{Q^{\left(\theta \right)}}(❬S-u_a,S-u_a❭❬\delta u^{\left(\theta -1\right)}\left(\zeta \right)-u_a,w^{\left(\theta \right)}-u_a❭\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-❬S-u_a,w^{\left(\theta \right)}-u_a❭❬\delta u^{\left(\theta -1\right)}\left(\zeta \right)-u_a,S-u_a❭),\end{array}$$

$$\begin{array}{c}{\beta }^{\left(\theta \right)}=\frac{1}{Q^{\left(\theta \right)}}(❬w^{\left(\theta \right)}-u_a,w^{\left(\theta \right)}-u_a❭❬\delta u^{\left(\theta -1\right)}\left(\zeta \right)-u_a,S-u_a❭\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-❬S-u_a,w^{\left(\theta \right)}-u_a❭❬\delta u^{\left(\theta -1\right)}\left(\zeta \right)-u_a,w^{\left(\theta \right)}-u_a❭),\end{array}$$

$$Q^{\left(\theta \right)}=❬w^{\left(\theta \right)}-u_a,w^{\left(\theta \right)}-u_a❭❬S-u_a,S-u_a❭-❬S-u_a,w^{\left(\theta \right)}-u_a❭{}^2,$$

where $❬u,w❭={\int }_{x_{min}}^{x_{max}}{\int }_{y_{min}}^{y_{max}}u\left(x,y\right)w\left(x,y\right)dxdy$ for the selected parameterization will be written in the form $❬u,w❭=\Delta x\Delta y{\sum }_{n=0}^{M_x-1}{\sum }_{m=0}^{M_y-1}{u}_{n,m}{w}_{n,m}$; S is upward recalculation for the difference in heights H of the unit field which is equal to $u_a+1$ if $\left(x,y\right)\in [x_{min};x_{max})\times [y_{min};y_{max})$ or $u_a$ if $\left(x,y\right)\notin [x_{min};x_{max})\times [y_{min};y_{max})$

$$S\left(x,y\right)=u_a+\mathrm{arctg}\left(\frac{\left(x-\xi \right)\left(y-\eta \right)}{H\sqrt{{\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+H^2}}\right)|\begin{array}{c}{x}_{max} \\ \\ \xi =x_{min}\end{array}|\begin{array}{c}{y}_{max} \\ \\ \eta =y_{min}\end{array},$$

$$S_{n,m}=u_a+\frac{1}{2\pi \Delta x\Delta y}\nu \left(x-\xi ,y-\eta ,H\right)|\begin{array}{c}{x}_{max} \\ \\ \xi =x_{min}\end{array}|\begin{array}{c}{y}_{max} \\ \\ \eta =y_{min}\end{array}|\begin{array}{c}{x}_{n+1} \\ \\ x=x_n\end{array}|\begin{array}{c}{y}_{m+1} \\ \\ y=y_m\end{array}.$$

3.

Calculate $u_{n,m}^{\left(\theta \right)}\left(z\right)=u_{n,m}^{\left(\theta -1\right)}\left(z\right)+{\alpha }^{\left(\theta \right)}\left(\delta u_{n,m}^{\left(\theta -1\right)}\left(\zeta \right)-u_a\right)+{\beta }^{\left(\theta \right)}$.

4.

As we can see, $u^{\left(\theta \right)}\left(z\right)\stackrel{up\left(H\right)}{\Rightarrow }{u}^{\left(\theta \right)}\left(\zeta \right)=u^{\left(\theta -1\right)}\left(\zeta \right)+{\alpha }^{\left(\theta \right)}\left(w^{\left(\theta \right)}-u_a\right)+{\beta }^{\left(\theta \right)}\left(S-u_a\right),$ hence we calculate $\delta u_{n,m}^{\left(\theta \right)}\left(\zeta \right)=u_{n,m}\left(\zeta \right)-u_{n,m}^{\left(\theta \right)}\left(\zeta \right)=\delta u_{n,m}^{\left(\theta -1\right)}\left(\zeta \right)+u_{n,m}^{\left(\theta -1\right)}\left(\zeta \right)-u_{n,m}^{\left(\theta \right)}\left(\zeta \right)=\delta u_{n,m}^{\left(\theta -1\right)}\left(\zeta \right)-{\alpha }^{\left(\theta \right)}\left(w_{n,m}^{\left(\theta \right)}-u_a\right)-{\beta }^{\left(\theta \right)}\left(S_{n,m}-u_a\right)$.

5.

Evaluate termination criteria: $\Vert \delta u^{\left(\theta \right)}\left(\zeta \right)-u_a\Vert =\sqrt{\delta u^{\left(\theta \right)}\left(\zeta \right)-u_a,\delta u^{\left(\theta \right)}\left(\zeta \right)-u_a}<\epsilon$ (reaching the required accuracy ε of the field approximation).

If the expressions (6, 8) are used in the algorithm above at stage 1 to recalculate upwards, then the approximated value $u_{n,m}\left(\zeta \right)$ will denote the average of the field over the regions $[x_n;x_{n+1})\times [y_m;y_{m+1})$, which should be obtained as a result of upward recalculation of the required field $u\left(z\right)$. If, at the first stage, we use the expression (9) with the condition $\left({\xi }_0,{\eta }_0\right)\in [x_0;x_1)\times [y_0;y_1)$, then the values $u_{n,m}\left(\zeta \right)$ will denote the values of the recalculated field at the specific points of the form $\left({\xi }_0+n\Delta x,{\eta }_0+m\Delta y,\zeta \right)$. Both options can be used, but it is worth noting that the iterative process will be more stable if we use the average over the region instead of the value at just one point in that region. Also, in the second case $\delta u_{n,m}^{\left(\theta -1\right)}\left(\zeta \right)$is the value of $\delta u^{\left(\theta -1\right)}\left(\zeta \right)$ at one specific point in the region $[x_n;x_{n+1})\times [y_m;y_{m+1})$. For the operation up(H) in the chosen piecewise constant parameterization, it is required to “continue” $\delta u_{n,m}^{\left(\theta -1\right)}\left(\zeta \right)$ over the entire region, which introduces more error than the “fair” averaging. In addition, the formula for the dot product $\langle u,w\rangle$ given in step 2 yields better result with averaging.

As our experiments show the above method converges even without classical regularization (in both versions: with averaging and “pointwise”), but, perhaps, will not reach zero error ε. However, for the purposes described in the next section, we introduce a formal regularization. Since the kernel of the integral (1) is symmetric and positive definite, we can apply Lavrentiev’s regularization [22]. The regularized equation is:

$$u\left(\xi ,\eta ,\zeta \right)=\kappa u\left(\xi ,\eta ,z\right)+\frac{\zeta -z}{2\pi }\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\frac{u\left(x,y,z\right)dxdy}{{\left({\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+{\left(z-\zeta \right)}^2\right)}^{\frac{3}{2}}},$$

(10)

where $\zeta$ and $z$ are fixed, $\kappa >0$ is regularization parameter. This equation has exactly one solution in $L_2$ for any left side function in $L_2$. The solution depends continuously on $\kappa$ [22]. The Equation (10) can be solved by the modified method of local corrections proposed earlier for Equation (1) with piecewise constant parameterization, but requires a modified method of calculating $w^{\left(\theta \right)}$ on the first step: instead of $\delta u^{\left(\theta -1\right)}\left(\zeta \right)\stackrel{up\left(H\right)}{\Rightarrow }{w}^{\left(\theta \right)}$ should be $\delta u^{\left(\theta -1\right)}\left(\zeta \right)\stackrel{up\left(H\right)}{\Rightarrow }{\widehat{w}}^{\left(\theta \right)}$, $w^{\left(\theta \right)}={\widehat{w}}^{\left(\theta \right)}+\kappa \delta u^{\left(\theta -1\right)}\left(\zeta \right)$.

2.4. Applying the Recalculations to Recover the Density Model

We consider the proposed method for separating the gravitational effect $u\left(x,y,0|(-\infty ,-H]\right)$ of sources from the total observed field $g\left(x,y,0\right)$ in the half-space below certain depth $z=-H<0$: $g\left(x,y,0\right)$ $\stackrel{up\left(H\right)}{\Rightarrow }$ $g\left(x,y,H\right)$ $\stackrel{down\left(2H,\kappa \right)}{\Rightarrow }$ $u\left(x,y,-H|(-\infty ,-H]\right)$ $\stackrel{up\left(H\right)}{\Rightarrow }$ $u\left(x,y,0|(-\infty ,-H]\right)$. Let $g\left(x,y,0\right)$ be specified exactly by its piecewise constant representation. If the operations up(H) and down(H,0) were performed precisely, in the analytical form, then the result would be $u\left(x,y,0|(-\infty ,-H]\right)\equiv g\left(x,y,0\right)$. In the proposed recalculations, however, there are inaccuracies: the averaging of the exact results of up(H) in the cells $[x_n;x_{n+1})\times [y_m;y_{m+1})$and equating the values up(H) and down(H,0) to horizontal asymptote $u_a$ outside the target region $[x_{min};x_{max})\times [y_{min};y_{max})$.

The authors have carried out a large number of numerical experiments with the data of the real measured fields. In the worst cases of “high-frequency” fields and large H (of the order of 100 km), the residual $\frac{\Vert g\left(x,y,0\right)-u\left(x,y,0|(-\infty ,-H]\right)\Vert }{\Vert g\left(x,y,0\right)\Vert }$ was not more than 10% and not more than 1% for H < 40 km. So, if we strive for the greatest accuracy when implementing the recalculation scheme using the analytic continuation of harmonic functions, then the separation of the fields will not work. Actually, for the purposes of separation, the formal regularization was introduced into the operation down(H,κ).

The result $u\left(x,y,0|(-\infty ,-H]\right)$ of the three-stage recalculation scheme continuously depends on the regularization parameter κ. We separate the observed field $g\left(x,y,0\right)$ along L horizontal layers with depth intervals ${\left((-H_{i+1},-H_i]\right)}_{i=0}^{L-1}$ and consider $u\left(x,y,0|(-H_{i+1},-H_i]\right)=u\left(x,y,0|(-\infty ,-H_i]\right)-u\left(x,y,0|(-\infty ,-H_{i+1}]\right)$ is the field of this layer, while the fields $u\left(x,y,0|(-\infty ,-H_i]\right)$ are obtained using the regularization parameters ${\kappa }_i$, with the exception of $u\left(x,y,0|(-\infty ,-H_0]\right)=g\left(x,y,0\right)$, $H_0=0$, ${\kappa }_0=0$. If we assume that the entire observed field is generated by the masses in the layer $(-H_L,0]$, then it is necessary to ensure $u\left(x,y,0|(-\infty ,-H_L]\right)\approx u_a$, for this ${\kappa }_L$ should large enough. To fulfill the condition of “continuous” joining of the separated fields $u\left(x,y,0|(-H_{i+1},-H_i]\right)$ of adjacent layers, it is required to choose ${\kappa }_i$ in ascending order between ${\kappa }_0=0$ and ${\kappa }_L$ without abrupt transitions from ${\kappa }_{i-1}$ to ${\kappa }_i$. Increasing ${\kappa }_i$ with depth produces a sequence of separated fields in accordance with the physical meaning: the deeper layers are more “smoothed”. Of course, there are infinitely many options for “continuously” increasing sequence of ${\kappa }_i$, and different options can give significantly different recalculations and, accordingly, different density models. This is precisely where the non-uniqueness of the solution to the linear inverse gravimetry problem is manifested. Here, the choice falls on the shoulders of the interpreter and is subjective. One helpful fact is that small changes in the sequence of ${\kappa }_i$ will lead to small changes in the separated fields and the resulting density model.

Thus, κ is used in recalculations not for regularization really, but as a continuous smoothing factor, which cannot be replaced by, for example, several passes through the observed field with a discrete Gaussian filter in the “time” domain, because even one such path for the upper layers is too many. Frequency domain filters can also be used for the purpose of separating fields, but it seems they are less justified in terms of binding to the specific depths.

2.5. Inverse Gravity Problem for the Layered Medium Model

Calculations of the 3D density distribution $\rho \left(x,y,z\right)$ in an inhomogeneous region D based on field values $g\left(\xi ,\eta ,\zeta \right)$ are implemented by inversion of integral operator [21] on the right hand side of (11), where g is a known function, γ is the gravitational constant.

$$g\left(\xi ,\eta ,\zeta \right)=\gamma \underset{D}{\overset{}{{\displaystyle \int }}}\frac{\left(z-\zeta \right)\rho \left(x,y,z\right)dxdydz}{{\left({\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+{\left(z-\zeta \right)}^2\right)}^{\frac{3}{2}}},$$

(11)

Mathematically, this is an ill-posed problem and its solution depends drastically on small variations in the initial field data g. However, if we select a density class with only lateral density variations, the determination of the density distribution in the horizontal layer will be stable [23]. Let $D$ be the rectangular prism $D=[x_{min};x_{max})\times [y_{min};y_{max})\times [z_{min};z_{max})$ filled with mass and $g\left(\xi ,\eta ,\zeta \right)$ be the observed gravity field separated for the layer $[z_{min};z_{max})$ using described above recalculation technique. Weseek the density in the form of $\rho \left(x,y,z\right)={\rho }^{\left(0\right)}\left(x,y,z\right)+{\rho }_0\left(z\right)\mathsf{\Phi }\left(x,y\right)$, where ${\rho }^{\left(0\right)}\left(x,y,z\right)$ is the density of initial model and it is assumed to be known. Further, we assume ${\rho }_0\left(z\right)$ may be approximated by some kind of the initial model analysis.

The 2D function $\mathsf{\Phi }\left(x,y\right)$ is defined based on the Fredholm integral equation:

$$\delta g^{\left(0\right)}\left(\xi ,\eta ,\zeta \right)=\gamma \underset{x_{min}}{\overset{x_{max}}{\int }}\underset{y_{min}}{\overset{y_{max}}{\int }}\mathsf{\Phi }\left(x,y\right)K\left(x,y,\xi ,\eta ,\zeta \right)dydx,$$

(12)

where

$$K\left(x,y,\xi ,\eta ,\zeta \right)=\underset{z_{min}}{\overset{z_{max}}{\int }}\frac{\left(z-\zeta \right){\rho }_0\left(z\right)}{{\left({\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+{\left(z-\zeta \right)}^2\right)}^{3/2}}dz,$$

$$\delta g^{\left(0\right)}\left(\xi ,\eta ,\zeta \right)=g\left(\xi ,\eta ,\zeta \right)-U^{\left(0\right)}\left(\xi ,\eta ,\zeta \right),$$

$$U^{\left(0\right)}\left(\xi ,\eta ,\zeta \right)=\gamma \underset{D}{\overset{}{\int }}\frac{\left(z-\zeta \right){\rho }^{\left(0\right)}\left(x,y,z\right)}{{\left({\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2+{\left(z-\zeta \right)}^2\right)}^{3/2}}dxdydz,$$

where $U^{\left(0\right)}\left(\xi ,\eta ,\zeta \right)$ is the field of the initial model.

The numerical solution of the Equation (12) is based on the local corrections method [20]. This approach allows one to solve the problem without usage of non-linear minimization and, thus, reduce the calculation time. Iterative algorithm for the refinement values of $\mathsf{\Phi }\left(x,y\right)$ minimizes the discrepancy $\Vert \delta g\Vert$ of the observed $g\left(\xi ,\eta ,\zeta \right)$ and model $U\left(\xi ,\eta ,\zeta \right)$ fields. First, the field of the initial model $U^{\left(0\right)}\left(\xi ,\eta ,\zeta \right)$ is calculated. Thus, the residual field is $\delta g^{\left(0\right)}=g-U^{\left(0\right)}$. We assume that ${\mathsf{\Phi }}^{\left(0\right)}\equiv 0$, $\zeta =const$, $\zeta \notin \left(z_{min},z_{max}\right)$, $\left(\xi ,\eta \right)\in \left[x_{min},x_{max}\right]\times \left[y_{min},y_{max}\right]$.

The following steps are repeated in a loop (for the iterations θ ≥ 1).

1.

Calculate

$$\delta U^{\left(\theta \right)}\left(\xi ,\eta \right)=\gamma \underset{x_{min}}{\overset{x_{max}}{\int }}\underset{y_{min}}{\overset{y_{max}}{\int }}\delta g^{\left(\theta -1\right)}\left(x,y,\zeta \right)\mathrm{K}\left(x,y,\xi ,\eta ,\zeta \right)dydx.$$

2.

Determine ${\alpha }^{\left(\theta \right)}$ and ${\beta }^{\left(\theta \right)}$:

$${\alpha }^{\left(\theta \right)}=\frac{1}{Q^{\left(\theta \right)}}(❬S,S❭❬\delta g^{\left(\theta -1\right)},\delta U^{\left(\theta \right)}❭-❬S,\delta U^{\left(\theta \right)}❭❬\delta g^{\left(\theta -1\right)},S❭),$$

$${\beta }^{\left(\theta \right)}=\frac{1}{Q^{\left(\theta \right)}}(❬\delta U^{\left(\theta \right)},\delta U^{\left(\theta \right)}❭❬\delta g^{\left(\theta -1\right)},S❭-❬S,\delta U^{\left(\theta \right)}❭❬\delta g^{\left(\theta -1\right)},\delta U^{\left(\theta \right)}❭),$$

$$Q^{\left(\theta \right)}=❬\delta U^{\left(\theta \right)},\delta U^{\left(\theta \right)}❭❬S,S❭-❬S,\delta U^{\left(\theta \right)}❭{}^2,$$

$$S\left(\xi ,\eta \right)=\gamma \underset{x_{min}}{\overset{x_{max}}{\int }}\underset{y_{min}}{\overset{y_{max}}{\int }}\mathrm{K}\left(x,y,\xi ,\eta ,\zeta \right)dydx.$$

3.

Calculate ${\mathsf{\Phi }}^{\left(\theta \right)}={\mathsf{\Phi }}^{\left(\theta -1\right)}+{\alpha }^{\left(\theta \right)}\delta g^{\left(\theta -1\right)}+{\beta }^{\left(\theta \right)}.$

4

Calculate $\delta g^{\left(\theta \right)}=\delta g^{\left(\theta -1\right)}-{\alpha }^{\left(\theta \right)}\delta U^{\left(\theta \right)}-{\beta }^{\left(\theta \right)}S.$

5.

Evaluate termination criteria: $\Vert \delta g^{\left(\theta \right)}\Vert =\sqrt{\delta g^{\left(\theta \right)},\delta g^{\left(\theta \right)}}<\epsilon$ (achieving the required accuracy ε of the field approximation).

Once the loop is complete, the excess density distribution ${\rho }_0\left(z\right){\mathsf{\Phi }}^{\left(\theta \right)}\left(x,y\right)$ approximates (up to a constant) the difference $\delta g^{\left(0\right)}$ between the observed and initial model fields. After summing the initial model with obtained distribution, we obtain the density model with the field that differs from the observed one by the amount of error of $\Vert \delta g^{\left(\theta \right)}\Vert$.

The described algorithm relies on the solution of the forward problem only (step 1). The lateral density is calculated using the difference between the observed field and the model field at each iteration step. For this reason, no classical regularization (in Tikhonov sense) is required. One can regulate convergence of the iterative process by coefficients ${\alpha }^{\left(\theta \right)}$ and ${\beta }^{\left(\theta \right)}$. Thus, the problem of calculation of corrective additive could be solved independently for all layers. Therefore, the uniqueness of the solution to the equation of the inverse problem for the lateral density [23] and the whole problem is guaranteed.

3. Experimental Results. Three-Dimensional Density Model of Middle-Urals Region

In this section, we demonstrate the described method on the case study of Middle-Urals, Russia, and neighboring regions. The territory is located at 56–60° N and 54–66° E, the model thickness is 80 km down from the Earth’s surface level. The given gravity field with a complete Bouguer correction (the combined global gravity field model XGM2019e_2159 [24]) is presented in Figure 1. The position of DSS profiles is drawn as black and white lines. The maximum depth of the profile sections is 80 km [25].

The initial density model frame is formed with the 2D density cuts along the regional profiles (Figure 2). The data between profiles are interpolated values (layer by layer) [3]. The mean value of density ${\rho }_0\left(z\right)$ (Figure 3) is calculated for each depth of the initial model (Figure 4). This 1D density distribution is continued outside the study volume and the excess density is calculated relative to it. As we seek 3D density distribution as a product of this function, ${\rho }_0\left(z\right)$ and some 2-dimensional corrective addition $\mathsf{\Phi }\left(x,y\right)$, the inverse problem for lateral density in a flat layer will be sufficiently stable [23,26].

The described interpretation algorithm of the gravity data was applied to the difference $\delta g^{\left(0\right)}$ of the observed field and the initial model field (Figure 5).

In order to evaluate the acceleration of calculations using the graphics accelerator (GPU) in comparison with the central processing unit (CPU), a series of field upward recalculations was carried out with different sampling grid resolutions. It should be noted that in the calculations, the discretization parameters of the grids of the recalculated and resulting fields were the same. Discretization parameters and computation times are shown in

Table 1. Dependence of the computation time on the discretization parameters of the field evaluation grid.

Nx × Ny	N = Nx × Ny, 106 Points	CPU Calculation Time, min	GPU Calculation Time, min
200 × 200	0.04	7.3 × 10−4	5 × 10−5
500 × 500	0.25	0.072	4.9 × 10−3
1000 × 1000	1	0.95	0.065
2000 × 2000	4	16.5	0.86
3000 × 3000	9	99.6	4.4
4000 × 4000	16	354.6	13.9
5000 × 5000	25	958.7	34.5

. The calculations were carried out with 24 cores (48 threads) of the AMD EPYC 7451 processor with a clock frequency of 2.9 GHz and one GPU AMD Radeon VII (Vega 20). Figure 6 and Figure 7 show the corresponding graphs of performance.

The Figure 8 and Figure 9 show the results of the layer-by-layer separation of the field $\delta g^{\left(0\right)}$.

A density model that satisfies the observed data was constructed using the described method. The initial ${\rho }_0\left(z\right)$ distribution (Figure 3) was used as the input data. Resulting density model (Figure 10) consists of prismatic elements; resolution of this model is 1236 × 1314 × 80. Figure 11 shows the change in the intervals of densities between the initial and the fitted model.

The stability of the proposed method is demonstrated here. To the observed field (Figure 5-(1)), at a random 33% of all points, we add Gaussian noise with zero mean and standard deviation of 6.72 mGal (33% of the standard deviation of the field). Figure 12 shows the result.

We calculate the difference between the noised field and the initial model field, shown in Figure 5-(2) Then, we apply the above-described recalculation filter to this difference for the layer-by-layer separation of the field up to 80 km in depth with a step of 1 km. Leaving the set ${\left\{{\kappa }_i\right\}}_{i=0}^{80}$ of smoothing parameters the same as in the example without noise, the separation result for four thick layers is shown in Figure 13. “Cuboid of the separated fields” for 1 km thick layers is shown in Figure 14.

As can be seen, most of the noise has been filtered out. This is not surprising, because the basis of our filter is Lavrentiev’s regularization. The remnants of the noise “penetrated” to a depth of about 7 km. It should be noted that if we chose a set ${\left\{{\kappa }_i\right\}}_{i=0}^{80}$ of smoothing parameters specifically for this case, then the noise could be completely neutralized. This, of course, would also be facilitated by the fact that our original observed field does not have a “high-frequency” component. From a practical point of view, the “noisy” (Figure 15) and “clean” (Figure 10) models coincide.

4. Conclusions

In this study, we offered a method of gravity field interpretation, which uses complete a priori information. The algorithm’s stability is ensured by the selection of the initial model and solving of the problem on the correctness set. The new algorithm is suggested as a preliminary pre-processing of gravity observations that allows the extraction of the gravitational effect of the layers. The modified method of local corrections is developed for the downward continuation of the gravity field. A new method for reconstruction of density distribution using gravity data is presented. Both algorithms are applied to the real gravity data for the Ural region, Russia.

The density model of the initial approximation and the fast original algorithms for solving the gravity problem with high resolution (using parallel computations) allows for the calculation of large-scale density models. Based on the results of the numerical modeling, it is possible to construct a volumetric (gradient) model of the layer-by-layer density distribution in the inhomogeneous layer and, within the framework of the obtained solution, restore the zones of the local inhomogeneity. We have created such 3D model for the study area. Volumetric density models are of considerable practical importance and can be used to justify the location of the prospecting and exploration work in the vicinity of promising areas.