Abstract
In computational mathematics, the iterative method is a mathematical procedure. This method uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the n-th approximation is derived from the previous ones. The iterative method is widely used to solve complex problems in engineering. In this paper, the iterative method is applied to inverse the subsurface interface with the gravity anomaly. First, the classical Parker-Oldenburg interface inversion formula was introduced and analogized to the downward continuation formula. Then, combined with the regular-integral downward continuation method, the iterative inversion formula of the gravity interface is derived. The iterative mode of the improved method suppresses high-frequency signals effectively. At the same time, there is no need to perform forward calculations in the iterative process. The model test shows that the proposed method can accurately calculate the depth of the interface. Finally, the proposed interface inversion method is applied to the Qinghai-Tibet Plateau, and the calculated Moho interface provides some geophysical data support for the geological interpretation of the area in the future.
1. Introduction
The iteration method is a mathematical process to solve the problem by finding a sequence of approximate solutions from an initial value. In contrast, direct methods attempt to solve the problem by a finite sequence of operations. In general, a direct method is always preferred if possible. However, when we encounter nonlinear equations or linear problems involving many variables, iterative methods are often the only choice. Some researchers [1] use iterative methods to solve a considerable number of mathematical engineering problems. Liu et al. [2] presented a comparative study of direct and iterative inversion methods to reconstruct the shear modulus distribution of linearly elastic solids. Mei et al. [3] solved the inverse problem in elasticity with an iterative method.
It is one of the classical problems in engineering geophysics to know the geometry of a density interface. The gravity anomaly is a comprehensive reflection of the undulations of various physical interfaces underground and the unevenness of lithology. In general, the nonuniformity of underground rock density is often related to the variation in some geological structures or the distribution of some minerals. Therefore, the study of crustal density interfaces using gravity data is of great practical significance for understanding the internal structure of the crust, studying regional geology, and detecting metal and nonmetal minerals. The inversion of the Moho surface and density basement is a typical application [4–7].
There are many different techniques to compute the geometry of a density interface related to a known gravity anomaly [8, 9]. The Parker-Oldenburg method is one of the most popular interface inversion methods. Parker [10] introduced the fast Fourier transform (FFT) in the forward calculation of the potential field interface. Oldenburg [11] proposed an iterative inversion method based on Parker’s interface forward method. Many improved algorithms [12–14] for the Parker-Oldenburg method have been presented in recent decades. To suppress the high-frequency amplification factor in the Parker-Oldenburg inversion, Reamer and Ferguson [15] proposed a linear variable-density interface inversion method. Guspi [16, 17] improved the Taylor series high-order term coefficients in Oldenburg inversion and proposed a noniterative and nonlinear gravity anomaly inversion method that improved the calculation speed. Fedi [18] proposed a method in the frequency domain, which has a fast calculation speed and a small amount of calculation. However, some of these methods [19–21] may make the inversion density interface too smooth and deviate from the real underground interface fluctuation due to the removal of useful high-frequency information. In some methods, the inversion process is still divergent due to the improper selection of filter factors. These methods do not fundamentally meet the accuracy and convergence required by density interface inversion. This paper presents an improved iterative algorithm for classical Oldenburg inversion. The divergent term in Parker-Oldenburg’s algorithm is equivalent to the downward continuation factor in the Fourier domain [20]. We derive an iterative inversion formula to calculate the density interface based on the Parker-Oldenburg interface inversion method and regular-integral downward continuation method. The algorithm does not need to perform forward calculation in the iterative process and can suppress high-frequency signals well, which provides better results.
2. Method
2.1. Parker-Oldenburg Forward Principle
A substance interface S underground is shown in Figure 1. It is assumed that the residual density ρ in the upper part of the interface S is 0. The gravitational field produced by the volume element dv at Q (ξ, η, ζ) under the S plane to the observation point P (x, y, 0) is . The z axis is vertical downward:

and is the gravitational constant.
The average depth of interface S is z0, and the function K is expanded by the Taylor series at ζ = z0:
is the n-th derivative of ζ, and using the above formula, we integrate K with ζ from (z0 + Δh) to z0, where Δh is the distance between the S interface and z0. Δh above is negative and below is positive. Δh is a function of .
Assuming that is only a function of , the gravity anomaly produced by all the residual materials below S to point P is
According to the convolution theorem, let and be constants, and the spectrum expression of is obtained:
; u and are the wavenumbers in the x and y directions, respectively. Formula (5) is the interface forward formula of the gravity anomaly in the frequency domain.
Parker-Oldenburg interface inversion and steps:
Using formula (5), we obtained
Formula (6) is the interface inversion calculation formula of the frequency domain.
The iterative inversion process is as follows:(1)We ignore the summation term at the right end of formula (6), and is obtained by Fourier transform of the observed gravity anomaly. is substituted into the formula, and by inverse Fourier transform, we calculate the initial value Δh0 of the depth Δh:(2)We bring the initial depth value to the forward formula (5) and obtain the initial spectrum . The first-order spectral difference is obtained by making a difference with the spectrum of the measured gravity anomaly in the first step, and the first-order spectral difference is obtained as follows:(3)The first-order spectral difference is introduced into inversion formula (6), and then the first-order residual value of Δh is calculated according to step 1:(4) is the first-order approximate value of the undulating depth Δh. Then, the second-order spectral difference and the second-order residual value of Δh are calculated according to step 2 and step 3, respectively. The n-order approximate value of can be obtained by repeating steps 2 and 3.
Until the n-order spectrum difference is less than the specified iteration cutoff constant or the number of interations n reaches the specified iteration number, is a small positive number. The flowchart is shown as Figure 2.

2.2. Regular-Integral Downward Continuation Method
The formulas of the upward and downward continuation in wavenumber domain are
and are the results of upward and downward continuation in the frequency domain, and z0 is the height of continuation.
Zeng et al. [22] proposed the wavenumber domain regular-integral iterative method for downward continuation. The expression of the regular-integral iterative method in the wavenumber domain is
, , and is a regular parameter. This method improves the accuracy, stability, and convergence speed of the algorithm.
2.3. Improved Parker-Oldenburg Method Based on the Regular-Integral Downward Continuation Method
The exponential amplification factor in the interface inversion expression (formula (6)) is the same as the downward continuation factor in formula (12). The characteristics of high-pass filtering will slow down the inversion convergence or make the convergence diverge.
Analogous to the upward continuation formula (11), formula (5) can also be written in another form:
Then, we have
Its form is similar to the downward continuation formula (); we make the following analogy: ‘∼’ is the analogical symbol.
By analogy with the downward continuation formula (13), the regular-integral downward continuation method is suitable for the Parker-Oldenburg interface inversion method. The following gravity inversion formula is obtained:
Ignoring the term where n is greater than 2 in at the right end of formula (17), we obtain the following formula:
And the formula of interface fluctuation is
When the number of iterations m at the right end of formula (19) tends to be infinite, the following formula is obtained:
The above formula shows that when , formula (20) is equal to formula (15), and the spectrum of the m-th iteration tends to the spectrum in the original inversion formula, indicating that the iteration method converges.
For the selection of regular parameters, we use the L curve method.
2.4. The Selection of Regularization Parameters
There are two ways to select regularization parameters: a priori and a posteriori. It is based on whether the noise level of the original data needs to be estimated in advance. It is difficult to give the noise level of the original data in advance. Therefore, the posteriori method is more practical. The most commonly used criteria are the generalized cross-check (GCV) criterion and the L curve criterion. This paper mainly uses the L curve criterion method to select Tikhonov regularization parameters.
2.4.1. The generalized Cross-Check (GCV) [23, 24]
The basic idea of cross-checking is as follows: if any point yi of the measurement data is removed, the selected regular parameter should be able to predict the change caused by the removed item. Although ordinary cross-checking depends on the specific ordering of the data, generalized cross-checking is invariant to the orthogonal transformation of the data vector y. The generalized cross function to be minimized in this method is defined aswhere is an arbitrary matrix mapping y to the solution , and trace represents the sum of the principal diagonal elements in the matrix.
Although GCV can solve many problems, it is difficult to find a good regularization parameter in some cases. One problem mentioned in the related literature [25] is that the GCV function may get a very flat minimum value, it is difficult to determine the minimum numerically. Another problem is that GCV sometimes mistakes noise for useful signals. GCV is quite effective for the nonuniformity of square error and non-Gaussian error. However, if the errors are highly correlated, the method may not obtain satisfactory results.
2.4.2. L Curve Criterion [26–28]
Taking log-log as the scale, (, ) forms a monotonously decreasing curve, as shown in Figure 3(a). Since the shape of this curve is similar to the letter ‘L’, it is called the L curve.

(a)

(b)
In the vertical part of the L curve, the regularization parameter and are small, and the regularized solution is in good agreement with the measured signal data. However, is more sensitive to the change in the regularization parameter, and the vertical part belongs to the underregularization state. In the horizontal part of the L curve, the regularization parameter is relatively large, and the regularization error is dominant. With the increase in , increases correspondingly. While almost does not change, the horizontal part belongs to the overregularization state. Therefore, to balance under regularization and overregularization, the regularization parameter is selected at the corner of the L curve (the angle between the vertical part and horizontal part). Usually, people choose the point with the greatest curvature (in Figure 3(b)) on the L curve as the corner of the L curve.
Then, the curvature function of the L curve with as the parameter is represents the derivation of . Through the parametric expression of the L curve, that is, the exact expression of functions and , the maximum curvature function can be directly calculated, and the corresponding regularization parameters can be obtained.
3. Examples
To verify the accuracy of the inversion method and show the inversion effect clearly, a two-dimensional -shaped density interface is established for the experiment. The -shaped density interface is a kind of common density interface, such as the ocean trench, half graben, and Moho surface beneath a subduction zone. It is of great significance to describe the density interface shape by gravity data for regional structure research, oil and gas exploration, and physical oceanography. The location of the two-dimensional model with simple morphology is shown in Figure 4(a). The density interface is composed of a -shaped depression, and the density difference between the upper and lower interfaces is 0.3 kg/cm3. The gravity anomaly generated by this interface is shown in Figure 4(b). is 0.06, and the number of iterations is 50. Then, we obtain the inverse interface through the proposed method in Figure 5. We can see that the shape of the inversed interface is close to the real model.

(a)

(b)

We have also designed a gravity interface with a density contrast of 0.3 kg/cm3 on both sides, and the range of the model is 50 50 km. It is composed of two bulges and a depression, as shown in Figure 6. Figure 6(a) is the interface depth map, and Figure 6(b) is the gravity anomaly caused by the interface. The minimum value of the gravity anomaly is −0.473 mGal, and the maximum value is 1.947 mGal. The shallow depths of the two points (positions A and B in Figure 6) are 1.8829 km and 1.7992 km, respectively, and the coordinates are (21 km, 35 km) and (35 km, 25.5 km); the largest depression (position C in Figure 6) is 2.0774 km, and its coordinates are (19 km, 14 km). Figures 7(a) and 7(b) show the L curve and the curve of the curvature function. The calculated regular parameter is 0.4824, substituted into formula (19), and the iteration time is 10. The result of the calculated interface inversion depth is shown in Figure 8(b). Figure 8(a) is the result of the classic Parker-Oldenburg method. Table 1 is a comparison table of the inversion results of the two methods.

(a)

(b)

(a)

(b)

(a)

(b)
Comparing the inversion results of Figure 8, both the Parker-Oldenburg method and the method proposed in this paper can invert the approximate depth of the interface, and the morphology is similar to the real interface. From Table 1, for the inversion depth of point A, the improved method is closer to the theoretical value than the result of the Parker-Oldenburg method. The coordinate position (21 km, 35 km) calculated by the improved method is the same as the coordinates of the real position. The coordinate result calculated by Parker-Oldenburg method is (21 km, 34.5 km). For point B, the coordinate positions calculated by the two methods are the same, but the proposed method calculates that the depth is closer to the true depth than the Parker method. The same situation is for point C. It can be considered that the inversion results of the improved Parker-Oldenburg method are closer to the depth of the real interface than the results of the classical Parker-Oldenburg inversion method.
Since there is a certain amount of noise in the measured data, we add Gaussian white noise with an average value of 0 mGal and a standard deviation of 0.1 mGal to the gravity anomaly, which is shown in Figure 9. The regular parameter calculated from the L curve method in Figure 10 is 1.2174. The inversion result after 15 iterations is shown in Figure 11(b), and Figure 11(a) is the calculation result of the classical Parker-Oldenburg inversion method. Table 2 is the comparison of the inversion results in the case of adding noise. It can be seen from the table that the improved Parker-Oldenburg iterative inversion method is closer to the true depth of the model for point A and point B. Although the calculation results for the depth of point C are not as good as the classical Parker-Oldenburg inversion method, the root mean square error is 0.0015 km, less than 0.0027 km calculated by the classical Parker-Oldenburg inversion method.



(a)

(b)
The theoretical model tests preliminarily determine that the Parker-Oldenburg method based on the wave-domain regular-integral iterative downward continue method can quickly invert the anomaly interface. The method could preserve the true shape of the interface and has a strong ability against noise. This shows that this method can solve the problem of interface morphology in geophysical engineering, and it can be applied to actual data processing. Next, we will proceed with the actual data processing with the iterative method.
4. Real Data
4.1. Geological Background
The collision and convergence between the Indian plate and the Asian plate in the Cenozoic created the present main structure of the Qinghai-Tibet Plateau. The formation and uplift mechanism of the Qinghai-Tibet Plateau is a very complex problem [29]. The plateau is not a homogeneous as a whole but is composed of several blocks with different development histories [30]. From south to north (Figure 12), the Himalayan, Lhasa, Qiangtang, Songpan-Ganzi, and Kunlun-Qaidam terranes were formed, which were separated by a series of nearly E-W trending suture zones [31–33]. Geophysical exploration has supported this view. It is speculated that the different terranes in the Qinghai-Tibet Plateau were uplifted unevenly in the process of amalgamation. The process of terrane amalgamation also leads to variations in crustal thickness between different terranes [34]. However, due to the limitation of the accuracy of the exploration and the combination process of terranes, the variation of the crustal thickness between different terranes and the shape of the Moho surface has not been accurately revealed.

4.2. Using Gravity Data to Invert the Moho Depth
Moho depth (crustal thickness) is one of the most important parameters to characterize the structure of the Earth. The crust of the Qinghai-Tibet Plateau is the thickest in the world. The location of the study area is 74–105° E, 26–40° N, which is the main part of the Qinghai-Tibet Plateau. Geophysicists have studied the deep structure of the Qinghai-Tibet Plateau for more than 60 years. One of the important studies is to explore the Moho depth of the Tibetan Plateau. It is mainly revealed by seismic methods, magnetotelluric methods, gravity, and magnetic exploration methods. Although the seismic method can well describe the shape of the Moho, the seismic profile data are limited, so it is impossible to fully describe the fluctuation of the Moho surface in the whole Qinghai-Tibetan Plateau. Therefore, the depth of the Moho should be inversed by the gravity anomaly. First, the gravity anomaly [35] is shown in Figure 13. It can be seen from the figure that the gravity anomaly in the central part of the study area is low, and the amplitude of the gravity anomaly varies from −550 mGal to 200 mGal. Such a large negative value reflects the existence of an extremely thick crust in this area. The variation in the gravity gradient at the suture zone is also relatively large, which is related to the drastic change in crustal thickness in this area.

The seismic profile data show that the crustal thickness of the Qinghai-Tibet Plateau changes greatly, and there are some faults on the Moho surface near the plate suture zone. From the Moho inversion results in Figure 14, it can be seen that the Moho inversed by this method is consistent with the seismic data. It can be seen from the figure that the Moho depth is between 50 and 76 km, the depth on the southern and northern margins of the plateau is relatively shallow, deepening to the interior of the plateau, and there are different degrees of dislocation on both sides of the main suture zone. For the whole plateau, the thickness of the plateau's crust seems to be a dustpan with thin sides and thick inside. The results show that the crustal mantle structure of the Qinghai-Tibet Plateau is obviously heterogeneous, with a thick crust and a thin lithosphere. The thickness of the crust is approximately 75 km in the interior of the Qinghai-Tibet Plateau and 54–64 km in the periphery.

To analyze the north–south and east–west Moho distribution characteristics of the Qinghai-Tibet Plateau, the north–south profile AB and the east–west profile CD are intercepted with blue lines in Figure 14. The Moho depth of north-south profile AB is as follows.
From Figure 15, it can be seen that the Moho depths of the Lhasa terrane, Qiangtang terrane, and Kunlun-Qaidam terrane in the central Qinghai-Tibet Plateau are the deepest [36, 37], ranging from 70 to 75 km, and the Moho surface near the suture zone has been faulted. The Moho depth on both sides of the profile gradually uplifted, and in the south, the depth was shallower than that in the north, reaching 54 km.

The Moho depth of the east–west profile CD is shown in Figure 16.

It can be seen from the above figure that the depth of the Moho depth reaches its lowest point at 79°E, the Moho depth continues to 74 km eastward until 85°E, and then it rises to the shallowest depth of 64 km. It can be inferred that the depth of the Moho in the east–west direction of the Qinghai-Tibet Plateau starts at 79°E, gradually rises eastward and rapidly rises to 65 km westward [38, 39].
5. Conclusion
In this paper, the downward continuation of the regularized integral iteration method is used to improve the Parker-Oldenburg interface inversion method of the gravity interface. At the same time, the reliability of the inversion interface is proven by model tests. Finally, the actual data are processed, and the Moho depth of the Qinghai-Tibet Plateau is inversed. We draw the following conclusions:(1)Through the derivation of the iterative mode, the final iterative scheme can be obtained. The scheme is simple and easy to operate, and the inversion results can be controlled by the number of iterations.(2)The interface inversion formula is compared with the downward continuation formula, and the regularization coefficient is used to suppress the noise in the inversion process. This method reduces the influence of noise and makes the inversion interface smoother.(3)The model tests verified the feasibility of the method. In the real data application, the inversion of the Moho surface in the Qinghai-Tibet Plateau is consistent with previous studies. It is proven that the iterative Parker-Oldenburg interface inversion method based on regular integrals is an effective method to calculate the gravity interface fluctuation, which can solve actual engineering geological problems.
Data Availability
The data are downloaded from https://topex.ucsd.edu.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors thank the International Gravimetric Bureau for providing the gravity data. This research was supported by the National Natural Science Foundation of China (41904129) and the National Key R&D Program of China (2020YFE0201300).