# The Iterative Inversion Method of the Gravity Interface Based on the Regular-Integral Downward Continuation Method

Authors: Heyu Wu; Wei Du; Yangyang Zhang
Journal: Mathematical Problems in Engineering
 (2022)

---

## 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 then-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.


---

## Body

## 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 interfaceS 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 dΔg. The z axis is vertical downward:(1)dΔg=GKρdv.Figure 1
The schematic diagram of interfaceS.K=ζ/x−ξ2+y−η2+z−ζ23/2 and G is the gravitational constant.The average depth of interfaceS is z0, and the function K is expanded by the Taylor series at ζ = z0:(2)Kx−ξ,y−η,ζ=∑n=0∞Knx−ξ,y−η,z0n!ζ−z0n.Knx−ξ,y−η,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 ζ=z0 is negative and below ζ=z0 is positive. Δh is a function of ξ,η.(3)∫h+ΔhhKx−ξ,y−η,ζdζ=∑n=0∞∫z0+Δhz0Knx−ξ,y−η,z0n!ζ−z0ndζ=∑n=0∞Knx−ξ,y−η,z0n+1!ζ−z0n+1|z0+Δhz0=−∑n=1∞Kn−1x−ξ,y−η,z0n!Δhnξ,η.Assuming thatρ is only a function of ξ,η, the gravity anomaly produced by all the residual materials below S to point P is(4)Δgx,y,0=∫−∞∞∫−∞∞Gρξ,ηdξdη∫h+ΔhhKx−ξ,y−η,z0dζ=−G∫−∞∞∫−∞∞∑n=1∞Kn−1x−ξ,y−η,z0n!Δhnξ,ηρξ,ηdξdη.According to the convolution theorem, letρξ,η and Mξ,η be constants, and the spectrum expression of Δg is obtained:(5)FΔg=−2πGe−z0kρ∑n=1∞−kn−1n!FΔhn.k=u2+v2; u and v 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(6)FΔh=−ez0k2πGρFΔg−∑n=2∞kn−1n!Fhn.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 FΔg is obtained by Fourier transform of the observed gravity anomaly. FΔg is substituted into the formula, and by inverse Fourier transform, we calculate the initial value Δh0 of the depth Δh:(7)Δh0=F−1−ez0k2πGρΔg.(2)
We bring the initial depth value to the forward formula (5) and obtain the initial spectrum FΔg. The first-order spectral difference is obtained by making a difference with the spectrum FΔg0 of the measured gravity anomaly in the first step, and the first-order spectral difference is obtained as follows:(8)F∂Δg1=FΔg−FΔg0.(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:(9)Δ∂h1=F−1−ez0k2πGρF∂Δg1.(4)
Δh1=Δh0+Δ∂h1 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 Δhn of Δh can be obtained by repeating steps 2 and 3.(10)Δhn=Δh0+∑i=1nΔ∂hn.Until the n-order spectrum differenceF∂Δgn 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.Figure 2
Flowchart of the iterative inversion process.


### 2.2. Regular-Integral Downward Continuation Method

The formulas of the upward and downward continuation in wavenumber domain are(11)FΔgup=FΔge−z0k,(12)FΔgdown=FΔgez0k.FΔgup and FΔgdown 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(13)FΔgdown1=RλFΔg,FΔgdownn+1=FΔgdownn+RλFΔg−ΦFΔgdownn.Rλ=Φ/Φ2+λ, Φ=ez0k, 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 factorez0k in the interface inversion expression (formula (6)) is the same as the downward continuation factor in formula (12). The characteristics of high-pass filtering ez0k 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:(14)FΔg=F−2πGρ∑n=1∞−kn−1n!Δhne−z0k.Then, we have(15)F−2πGρ∑n=1∞−kn−1n!Δhn=FΔgez0k.Its form is similar to the downward continuation formula (FΔgdown=FΔgez0k); we make the following analogy: ‘∼’ is the analogical symbol.(16)−2πGρ∑n=1∞−kn−1n!Δhn∼Δgdown.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:(17)F−2πGρ∑n=1∞−kn−1n!Δhn1=RλFΔg,F−2πGρ∑n=1∞−kn−1n!Δhnm+1=F−2πGρ∑n=1∞−kn−1n!Δhnm+RλFΔg−ΦF−2πGρ∑n=1∞−kn−1n!Δhnm.Ignoring the term wheren is greater than 2 in ∑n=1∞−kn−1/n!Δhn at the right end of formula (17), we obtain the following formula:(18)F−2πGρΔh1=RλFΔg,F−2πGρΔhm+1=F−2πGρΔhm+RλFΔg−ΦF−2πGρΔhm,=1−RλΦm+1−RλΦm−1+⋯+1−RλΦ+1F−2πGρΔh1=1−RλΦm+1−RλΦm−1+⋯+1−RλΦ+1RλFΔg=Φ−11−1−RλΦm+1FΔg.And the formula of interface fluctuation is(19)Δhm+1=−12πGρF−1Φ−11−1−RλΦm+1FΔg.When the number of iterationsm at the right end of formula (19) tends to be infinite, the following formula is obtained:(20)F−2πGρ∑n=1∞−kn−1n!Δhnm+1=limm⟶∞Φ−11−1−RλmFΔg=FΔgez0k.The above formula shows that whenm⟶∞, 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 theL 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 theL 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 pointyi 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 as(21)Gλ=Kxα−y22traceI−KKαI2,where KαI is an arbitrary matrix mapping y to the solution xα, 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, (Kxα−y2, xα2) 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.Figure 3
Schematic diagram of the L curve: (a) Schematic diagram of theL curve. (b) Curve of the curvature function.
(a)(b)In the vertical part of theL curve, the regularization parameter and Kxα−y2 are small, and the regularized solution is in good agreement with the measured signal data. However, xα2 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 α, Kxα−y2 increases correspondingly. While xα2 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 uα (the angle between the vertical part and horizontal part). Usually, people choose the point with the greatest curvature kα (in Figure 3(b)) on the L curve as the corner of the L curve.(22)uα=logKxα−y2,vα=logxα.Then, the curvature function of theL curve with α as the parameter is(23)kα=u′v′′u′′v′u′2+v′23/2,′ represents the derivation of α. Through the parametric expression of the L curve, that is, the exact expression of functions uα and vα, the maximum curvature function can be directly calculated, and the corresponding regularization parameters can be obtained.


## 2.1. Parker-Oldenburg Forward Principle

A substance interfaceS 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 dΔg. The z axis is vertical downward:(1)dΔg=GKρdv.Figure 1
The schematic diagram of interfaceS.K=ζ/x−ξ2+y−η2+z−ζ23/2 and G is the gravitational constant.The average depth of interfaceS is z0, and the function K is expanded by the Taylor series at ζ = z0:(2)Kx−ξ,y−η,ζ=∑n=0∞Knx−ξ,y−η,z0n!ζ−z0n.Knx−ξ,y−η,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 ζ=z0 is negative and below ζ=z0 is positive. Δh is a function of ξ,η.(3)∫h+ΔhhKx−ξ,y−η,ζdζ=∑n=0∞∫z0+Δhz0Knx−ξ,y−η,z0n!ζ−z0ndζ=∑n=0∞Knx−ξ,y−η,z0n+1!ζ−z0n+1|z0+Δhz0=−∑n=1∞Kn−1x−ξ,y−η,z0n!Δhnξ,η.Assuming thatρ is only a function of ξ,η, the gravity anomaly produced by all the residual materials below S to point P is(4)Δgx,y,0=∫−∞∞∫−∞∞Gρξ,ηdξdη∫h+ΔhhKx−ξ,y−η,z0dζ=−G∫−∞∞∫−∞∞∑n=1∞Kn−1x−ξ,y−η,z0n!Δhnξ,ηρξ,ηdξdη.According to the convolution theorem, letρξ,η and Mξ,η be constants, and the spectrum expression of Δg is obtained:(5)FΔg=−2πGe−z0kρ∑n=1∞−kn−1n!FΔhn.k=u2+v2; u and v 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(6)FΔh=−ez0k2πGρFΔg−∑n=2∞kn−1n!Fhn.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 FΔg is obtained by Fourier transform of the observed gravity anomaly. FΔg is substituted into the formula, and by inverse Fourier transform, we calculate the initial value Δh0 of the depth Δh:(7)Δh0=F−1−ez0k2πGρΔg.(2)
We bring the initial depth value to the forward formula (5) and obtain the initial spectrum FΔg. The first-order spectral difference is obtained by making a difference with the spectrum FΔg0 of the measured gravity anomaly in the first step, and the first-order spectral difference is obtained as follows:(8)F∂Δg1=FΔg−FΔg0.(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:(9)Δ∂h1=F−1−ez0k2πGρF∂Δg1.(4)
Δh1=Δh0+Δ∂h1 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 Δhn of Δh can be obtained by repeating steps 2 and 3.(10)Δhn=Δh0+∑i=1nΔ∂hn.Until the n-order spectrum differenceF∂Δgn 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.Figure 2
Flowchart of the iterative inversion process.


## 2.2. Regular-Integral Downward Continuation Method

The formulas of the upward and downward continuation in wavenumber domain are(11)FΔgup=FΔge−z0k,(12)FΔgdown=FΔgez0k.FΔgup and FΔgdown 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(13)FΔgdown1=RλFΔg,FΔgdownn+1=FΔgdownn+RλFΔg−ΦFΔgdownn.Rλ=Φ/Φ2+λ, Φ=ez0k, 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 factorez0k in the interface inversion expression (formula (6)) is the same as the downward continuation factor in formula (12). The characteristics of high-pass filtering ez0k 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:(14)FΔg=F−2πGρ∑n=1∞−kn−1n!Δhne−z0k.Then, we have(15)F−2πGρ∑n=1∞−kn−1n!Δhn=FΔgez0k.Its form is similar to the downward continuation formula (FΔgdown=FΔgez0k); we make the following analogy: ‘∼’ is the analogical symbol.(16)−2πGρ∑n=1∞−kn−1n!Δhn∼Δgdown.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:(17)F−2πGρ∑n=1∞−kn−1n!Δhn1=RλFΔg,F−2πGρ∑n=1∞−kn−1n!Δhnm+1=F−2πGρ∑n=1∞−kn−1n!Δhnm+RλFΔg−ΦF−2πGρ∑n=1∞−kn−1n!Δhnm.Ignoring the term wheren is greater than 2 in ∑n=1∞−kn−1/n!Δhn at the right end of formula (17), we obtain the following formula:(18)F−2πGρΔh1=RλFΔg,F−2πGρΔhm+1=F−2πGρΔhm+RλFΔg−ΦF−2πGρΔhm,=1−RλΦm+1−RλΦm−1+⋯+1−RλΦ+1F−2πGρΔh1=1−RλΦm+1−RλΦm−1+⋯+1−RλΦ+1RλFΔg=Φ−11−1−RλΦm+1FΔg.And the formula of interface fluctuation is(19)Δhm+1=−12πGρF−1Φ−11−1−RλΦm+1FΔg.When the number of iterationsm at the right end of formula (19) tends to be infinite, the following formula is obtained:(20)F−2πGρ∑n=1∞−kn−1n!Δhnm+1=limm⟶∞Φ−11−1−RλmFΔg=FΔgez0k.The above formula shows that whenm⟶∞, 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 theL 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 theL 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 pointyi 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 as(21)Gλ=Kxα−y22traceI−KKαI2,where KαI is an arbitrary matrix mapping y to the solution xα, 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, (Kxα−y2, xα2) 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.Figure 3
Schematic diagram of the L curve: (a) Schematic diagram of theL curve. (b) Curve of the curvature function.
(a)(b)In the vertical part of theL curve, the regularization parameter and Kxα−y2 are small, and the regularized solution is in good agreement with the measured signal data. However, xα2 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 α, Kxα−y2 increases correspondingly. While xα2 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 uα (the angle between the vertical part and horizontal part). Usually, people choose the point with the greatest curvature kα (in Figure 3(b)) on the L curve as the corner of the L curve.(22)uα=logKxα−y2,vα=logxα.Then, the curvature function of theL curve with α as the parameter is(23)kα=u′v′′u′′v′u′2+v′23/2,′ represents the derivation of α. Through the parametric expression of the L curve, that is, the exact expression of functions uα and vα, the maximum curvature function can be directly calculated, and the corresponding regularization parameters can be obtained.


## 2.4.1. The generalized Cross-Check (GCV) [23, 24]

The basic idea of cross-checking is as follows: if any pointyi 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 as(21)Gλ=Kxα−y22traceI−KKαI2,where KαI is an arbitrary matrix mapping y to the solution xα, 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, (Kxα−y2, xα2) 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.Figure 3
Schematic diagram of the L curve: (a) Schematic diagram of theL curve. (b) Curve of the curvature function.
(a)(b)In the vertical part of theL curve, the regularization parameter and Kxα−y2 are small, and the regularized solution is in good agreement with the measured signal data. However, xα2 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 α, Kxα−y2 increases correspondingly. While xα2 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 uα (the angle between the vertical part and horizontal part). Usually, people choose the point with the greatest curvature kα (in Figure 3(b)) on the L curve as the corner of the L curve.(22)uα=logKxα−y2,vα=logxα.Then, the curvature function of theL curve with α as the parameter is(23)kα=u′v′′u′′v′u′2+v′23/2,′ represents the derivation of α. Through the parametric expression of the L curve, that is, the exact expression of functions uα and vα, 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-dimensionalV-shaped density interface is established for the experiment. The V-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 V-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.Figure 4
The model position: (a) the depth of theV-shaped interface, (b) the gravity anomaly generated by Figure 4(a).
(a)(b)Figure 5
The comparison of calculation depth and model depth.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.Figure 6
(a) The depth of the model and (b) the gravity anomaly.
(a)(b)Figure 7
(a) L curve, (b) curve of curvature function.
(a)(b)Figure 8
The inversion results: (a) Parker-Oldenburg inversion calculation results, (b) the results of the improved method in this paper.
(a)(b)Table 1
The comparison of inversion results in the model.
Bulge A (km)Bulge B (km)Depression C (km)Root mean square error (km)CoordinateDepthCoordinateDepthCoordinateDepthModel(21, 35)1.8829(35, 25.5)1.7992(19, 14)2.0774Parker-Oldenburg method(21, 34.5)1.8240(35, 25.5)1.7891(19, 14)2.07310.0027Improved method(21, 35)1.8279(35, 25.5)1.7956(19, 14)2.07510.0015Comparing the inversion results of Figure8, 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 Figure9. 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.Figure 9
Gravity anomaly of the model after adding noise.Figure 10
Curve of curvature function.Figure 11
The inversion results: (a) Parker-Oldenburg inversion calculation results; (b) the results of the improved method in this paper.
(a)(b)Table 2
The comparison of inversion results in the model with noise.
Bulge A (km)Bulge B (km)Depression C (km)Root mean square error (km)CoordinateDepthCoordinateDepthCoordinateDepthModel(21, 35)1.8829(35, 25.5)1.7992(19, 14)2.0774Parker-Oldenburg method(21, 34.5)1.8236(35, 25.5)1.7894(19, 14)2.07340.0027Improved method(21, 35)1.8306(35, 25.5)1.7993(19, 14)2.07230.0015The 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.Figure 12
Research area. IYS : Indus-Yarlung Zangbo suture; BNS : Bangong-Nujiang suture; JS : Jinsha River suture; AKMS : Ayimaqin-Kunlun-Mutztagh suture; SQS : South Qilian suture; KDS : Kudi suture; STDS : South Tibet Detachment System; MCT : Main Central Thrust.


### 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.Figure 13
Gravity anomaly in the study 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 Figure14, 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.Figure 14
Moho depth in study area.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 Figure14. The Moho depth of north-south profile AB is as follows.From Figure15, 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.Figure 15
Moho depth of AB in Figure14.The Moho depth of the east–west profile CD is shown in Figure16.Figure 16
Moho depth of CD in Figure14.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].


## 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.Figure 12
Research area. IYS : Indus-Yarlung Zangbo suture; BNS : Bangong-Nujiang suture; JS : Jinsha River suture; AKMS : Ayimaqin-Kunlun-Mutztagh suture; SQS : South Qilian suture; KDS : Kudi suture; STDS : South Tibet Detachment System; MCT : Main Central Thrust.


## 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.Figure 13
Gravity anomaly in the study 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 Figure14, 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.Figure 14
Moho depth in study area.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 Figure14. The Moho depth of north-south profile AB is as follows.From Figure15, 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.Figure 15
Moho depth of AB in Figure14.The Moho depth of the east–west profile CD is shown in Figure16.Figure 16
Moho depth of CD in Figure14.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.



---

*Source: 1836384-2022-02-18.xml*