# The Calculation of High-Order Vertical Derivative in Gravity Field by Tikhonov Regularization Iterative Method

Authors: Wei Du; Yangyang Zhang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In mathematics, statistics, and computer science, particularly in the fields of machine learning and inverse problems, regularization is a process of introducing additional information in order to solve an ill-posed problem or to prevent overfitting. The Tikhonov regularization method is widely used to solve complex problems in engineering. The vertical derivative of gravity can highlight the local anomalies and separate the horizontal superimposed abnormal bodies. The higher the order of the vertical derivative is, the stronger the resolution is. However, it is generally considered that the conversion of the high-order vertical derivative is unstable. In this paper, based on Tikhonov regularization for solving the high-order vertical derivatives of gravity field and combining with the iterative method for successive approximation, the Tikhonov regularization method for calculating the vertical high-order derivative in gravity field is proposed. The recurrence formula of Tikhonov regularization iterative method is obtained. Through the analysis of the filtering characteristics of this method, the high-order vertical derivative of gravity field calculated by this method is stable. Model tests and practical data processing also show that the method is of important theoretical significance and practical value.


---

## Body

## 1. Introduction

Regularization is a process of introducing additional information in order to solve an ill-posed problem or to prevent overfitting, and Tikhonov regularization has been invented independently in many different contexts. It became widely known for its application to integral equations from the works of Tikhonov [1–3] and Phillips [4]. Some authors use the term Tikhonov-Phillips regularization. The finite-dimensional case was expounded by Hoerl who took a statistical approach, and Foster interpreted this method as a Wiener-Kolmogorov (Kriging) filter [5, 6]. Following Hoerl, it is known in the statistical literature as ridge regression.The vertical derivative is of important physical significance in the data processing of gravity exploration. The vertical derivative is used to classify the superimposed anomalies generated by the anomaly sources in different depths and different sizes. The calculation of vertical derivatives is an important part in other processing methods of gravity exploration, and the results directly affect the accuracy of calculations with these methods.The calculation of vertical derivative in gravity field can be divided into spatial-domain method and frequency-domain method. The spatial-domain methods include the finite element method and spline function method [7]. The frequency-domain methods include fast Fourier transform operator (FFT), Wiener filtering method [8], regularization method, new regularization method, and wavenumber domain iteration method. Other methods include the integrated second vertical derivative (ISVD) method [9], which combines the spatial domain with the frequency domain to calculate the vertical derivatives. ISVD method has become an important method to calculate the vertical derivatives stably, and it is widely used in data processing of gravity field currently.Because of the disadvantages such as the slow computing speed of the methods in the spatial domain, the high-order vertical derivatives of the gravity field are usually calculated in the frequency domain. However, the derivative response factor in the frequency domain has a high-frequency amplification effect. The higher the order of the derivative, the greater the sensitivity of the data to noise, making the obtained result unstable. In order to suppress the instability of the calculation results of high-order vertical derivatives, the Tikhonov regularization iteration method is proposed to obtain the vertical derivatives of gravity field. The physical meaning of the method is clear, the calculation is simple, and the result is stable. Through the model test and the application of the actual data, it is proved that the method has stability and advantages in calculating the vertical derivatives.


## 2. Method

### 2.1. The Tikhonov Regularization Iteration Method

In many applications such as spectroscopy [10], seismology [11], and medical imaging [12], the data are collected by convolving noise signals with detectors. The linear model of this process leads to the first class of integral equations:(1)∫01ks,txtde=y0s+es,where y0s+es is the measured signal and y0s is the real signal, es is the uncertain noise, the kernel function ks,t is the response function, and xt is the solution of the equation. Since the measured signal usually has only s finite values, the linear discrete form of the continuous formula (1) is(2)Kx=y0+e=y,where K is a matrix of dimension m×n, and it is assumed that m > n. Except for a few deconvolution problems, when a small change in data causes a serious deviation between the approximate solution and the true solution, this is an ill-posed problem. This is reflected in an ill-posed condition of matrix K of the discrete model, and it increases with the increase of the dimension. Therefore, if we try to solve formula (2) directly, we will get the solution vector seriously disturbed by noise.Therefore, some regularization operation is needed to reduce the influence of noise. The well-known regularization method is Tikhonov regularization, which is to choose a solution to solve the minimization problem:(3)minxKx−y22+λ2x22.In the above formula, the regularization parameterα mainly controls the relative size between the solution norm x22 and the residual norm Kx−y22. This method depends on how many regularization parameters of filter weight are introduced in the regularization calculation. Therefore, the key to this method is to find a regularization parameter that can reduce enough noise without losing too much information in the calculated solution.


### 2.2. The Selection of Tikhonov 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 to obtain the high-order vertical derivatives of the gravity field. Therefore, the a posteriori method is more practical. The most commonly used are the generalized cross-validation (GCV) criterion and theL curve criterion. This paper mainly uses the L curve criterion method to select Tikhonov regularization parameters.


#### 2.2.1. The Generalized Cross-Validation (GCV) [13, 14]

The basic idea of cross-validation 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-validation depends on the specific ordering of the data, generalized cross-validation is invariant to the orthogonal transformation of the data vector y. The generalized cross function to be minimized in this method is defined as(4)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 [15] is that the GCV function can have a very flat minimum value, so the minimum value itself may be difficult to locate by a number. Another problem is that GCV sometimes mistakes noise for useful signals. GCV is quite effective for nonuniformity of square error and non-Gaussian error. However, if the errors are highly correlated, the method may not get satisfactory results.


#### 2.2.2.L Curve Criterion [16–18]

Taking log-log as the scale, (Kxα−y2, xα2) forms a monotonously decreasing curve, as shown in Figure 1(a). Since the shape of this curve is similar to the letter “L,” it is called the L curve.Figure 1
(a) Schematic diagram ofL curve. (b) Curve of curvature function.
(a)(b)In the vertical part ofL curve, the regularization parameter and Kxα−y2 are small, and the regularized solution is in good agreement with the measured signal data. But xα2 is more sensitive to the change of regularization parameter, and the vertical part belongs to the underregularization state. In the horizontal part of L curve, the regularization parameter is relatively large, and the regularization error is dominant. With the increase of α, Kxα−y2 increases correspondingly, but xα2 almost does not change, so the horizontal part belongs to the overregularization state. Therefore, in order to balance underregularization and overregularization, the regularization parameter is selected at the corner of L curve uα (the angle between vertical part and horizontal part). Usually, people choose the point with the greatest curvature kα (in Figure 1(b)) on the L curve as the corner of the L curve.(5)uα=logKxα−y2,vα=logxα.Then the curvature function ofL curve with α as the parameter is(6)kα=u′v″−u″v′u′2+v′23/2,where ′ represents the derivation of α. Through the parametric expression of L curve, that is, the exact expressions of functions uα and vα, the maximum of curvature function can be directly calculated, and then the corresponding regularization parameters can be obtained.


### 2.3. The Basic Principle of Tikhonov Regularization Iteration Method in Frequency Domain

Many geophysical estimation problems are mathematically ill posed because they operate with insufficient data [19]. Regularization is a technique to make the estimation problem well posed by adding indirect constraints to the estimation model [20, 21]. The regularization method was first proposed by Tikhonov [22]. It has become an indispensable part of the inverse problem theory and has been widely used in geophysical problems [23–26].The relationship between gravity fieldTx,y and its vertical derivative Dmx,y (m represents the order) in the wavenumber domain is(7)Dmu,v=2πu2+v2mTu,v.Dmu,v and Tu,v are the spectra of the Fourier transform of Dmx,y and Tx,y. u and v are wavenumbers in the x and y directions, respectively. φm=2πu2+v2m is the vertical derivative operator. Due to the noise amplification characteristic of φm, the result Dmu,v is unstable, and the higher the derivative order is, the more unstable Dmu,v will be.For the instability problem of formula (7), Tikhonov regularization is a widely used method. Solve a minimized regularization functional; from formula (3),(8)minφ−m⋅Dmu,v−Tu,v22+αDmu,v22.2 represents the L2 norm. α is the regular parameter, used to balance the instability and smoothness. The regular operator in the article is selected according to the form of the derivative factor of the gravity field. The solution of the above minimization problem is(9)Dmαu,v=φm11+αφ2mTu,v,where Dmαu,v represents the conversion result of the regular vertical derivative. According to formula (9), Tikhonov regular vertical derivative conversion operator can be divided into two parts: conventional vertical derivative conversion operator part φm and Tikhonov regular low-pass filtering function part 1/1+αφ2m. It is equivalent to multiplying the conventional vertical derivative conversion operator by a regular low-pass filter. Then the approximation of Dmu,v (first-order approximation spectrum) after regular low-pass filtering can be expressed as(10)Dm1u,v=φm11+αφ2mTu,v.SubstituteDm1u,v into formula (7) to obtain an iterative approximate spectrum of T:(11)T1u,v=11+αφ2mTu,v.Then,Tu,v and T1u,v correct Dm1u,v to obtain the second iteration approximate spectrum of Dmu,v:(12)Dm2u,v=Dm1u,v+φm11+αφ2mTu,v−T1u,v.Repeat the above steps to obtain thenth-iteration approximate spectrum of Dmu,v:(13)Dmnu,v=Dmn−1u,v+φm11+αφ2mTu,v−Tn−1u,v.It can be seen from formula (7) that(14)11+αφ2mDmn−1u,v=φm1+αφ2mTn−1.Substitute it into formula (13):(15)Dmnu,v=φm1+αφ2mTu,v+1−11+αφ2mDmn−1u,v.It is easy to obtain the iterated general formula by mathematical induction:(16)Dmnu,v=1+1−11+αφ2m+⋯+1−11+αφ2mn−1φm1+αφ2mTu,v.Obviously,1+1−1/1+αφ2m+⋯+1−1/1+αφ2mn−1 is an equivalency sequence, with a common ratio of 0≤1−1/1+αφ2m<1; then,(17)Dmnu,v=1−1−11+αφ2mnφmTu,v.When the time of the iteration is infinite,(18)Dmu,v=limn⟶∞Dmnu,v=φmTu,v.The above formula shows that the Tikhonov regularization iteration method proposed in this paper can converge to the theoretical solution for high-order derivative conversion. When the order spectral difference is less than the specified iteration cut-off constantε or the number of the iterations reaches the specified iteration number, the iteration stops, where ε is a small positive number. The choice of regularization weight is to use the curvature of the L curve.


### 2.4. The Convergence and Filtering Properties of Tikhonov Regularization Iterative Method

Figures2 and 3 are the filtering characteristics of Tikhonov regularized iterative operators for the first vertical derivative and second vertical derivative, respectively. Figures 2(a) and 3(a) are the filter response characteristic curves of Tikhonov regularized iteration operator of the first vertical derivative and the second vertical derivative with different numbers of iterations, respectively.Figure 2
Filtering characteristics of Tikhonov regularized iterative operator for first vertical derivative. (a)α = 0.2. (b) The number of iterations is 3.
(a)(b)Figure 3
Filtering characteristics of Tikhonov regularized iterative operator for second vertical derivative. (a)α = 0.01. (b) The number of iterations is 3.
(a)(b)ω in the x-coordinate is the radial circular frequency. The value of regularization parameter of first vertical derivative α is 0.2, and α of the second vertical derivative is 0.01. By comparison and the filtering characteristics of the theoretical vertical derivative operators φm and Tikhonov regularization vertical derivative operators, it can be obtained that the theoretical vertical derivative factor rises sharply with the increase of frequency, and the higher the derivative order is, the stronger the amplification of the high-frequency components will be. Tikhonov regularization operator approximates the theoretical vertical derivative operator in low-frequency band, which ensures the accuracy of derivative conversion of useful signals in low frequency, but can suppress the noise in the high frequency. As the number of iterations increases, Tikhonov regularized iterative operators gradually approach the theoretical vertical derivative operators.In order to analyze the influence of the value of regularization parameterα on the filtering characteristics of Tikhonov regularization operator, Figures 2(b) and 3(b) are the filtering response characteristics of the Tikhonov regularized iterative operator of the first and second vertical derivatives when the number of iterations is 3, respectively. It can be seen that when the number of iterations is constant, with the increase of the order, the suppression of Tikhonov regularized iterative operator in high-frequency components is gradually enhanced, which ensures the stability of the high-order derivative’s calculation. The smaller the value of regularization parameter α, the weaker the suppression effect of Tikhonov regularization operator in high frequency. When the value is 0, it is the theoretical derivative operator. When the value of regularization parameter α is too large, it will suppress the low frequency, which will affect the filtering effect.It can be seen that the Tikhonov regularization iteration operator proposed in this paper can accurately approximate the theoretical derivative operator in the low frequency and effectively suppress the high-frequency component. Therefore, this method has the advantages of amplitude maintenance and stability.


## 2.1. The Tikhonov Regularization Iteration Method

In many applications such as spectroscopy [10], seismology [11], and medical imaging [12], the data are collected by convolving noise signals with detectors. The linear model of this process leads to the first class of integral equations:(1)∫01ks,txtde=y0s+es,where y0s+es is the measured signal and y0s is the real signal, es is the uncertain noise, the kernel function ks,t is the response function, and xt is the solution of the equation. Since the measured signal usually has only s finite values, the linear discrete form of the continuous formula (1) is(2)Kx=y0+e=y,where K is a matrix of dimension m×n, and it is assumed that m > n. Except for a few deconvolution problems, when a small change in data causes a serious deviation between the approximate solution and the true solution, this is an ill-posed problem. This is reflected in an ill-posed condition of matrix K of the discrete model, and it increases with the increase of the dimension. Therefore, if we try to solve formula (2) directly, we will get the solution vector seriously disturbed by noise.Therefore, some regularization operation is needed to reduce the influence of noise. The well-known regularization method is Tikhonov regularization, which is to choose a solution to solve the minimization problem:(3)minxKx−y22+λ2x22.In the above formula, the regularization parameterα mainly controls the relative size between the solution norm x22 and the residual norm Kx−y22. This method depends on how many regularization parameters of filter weight are introduced in the regularization calculation. Therefore, the key to this method is to find a regularization parameter that can reduce enough noise without losing too much information in the calculated solution.


## 2.2. The Selection of Tikhonov 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 to obtain the high-order vertical derivatives of the gravity field. Therefore, the a posteriori method is more practical. The most commonly used are the generalized cross-validation (GCV) criterion and theL curve criterion. This paper mainly uses the L curve criterion method to select Tikhonov regularization parameters.


### 2.2.1. The Generalized Cross-Validation (GCV) [13, 14]

The basic idea of cross-validation 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-validation depends on the specific ordering of the data, generalized cross-validation is invariant to the orthogonal transformation of the data vector y. The generalized cross function to be minimized in this method is defined as(4)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 [15] is that the GCV function can have a very flat minimum value, so the minimum value itself may be difficult to locate by a number. Another problem is that GCV sometimes mistakes noise for useful signals. GCV is quite effective for nonuniformity of square error and non-Gaussian error. However, if the errors are highly correlated, the method may not get satisfactory results.


### 2.2.2.L Curve Criterion [16–18]

Taking log-log as the scale, (Kxα−y2, xα2) forms a monotonously decreasing curve, as shown in Figure 1(a). Since the shape of this curve is similar to the letter “L,” it is called the L curve.Figure 1
(a) Schematic diagram ofL curve. (b) Curve of curvature function.
(a)(b)In the vertical part ofL curve, the regularization parameter and Kxα−y2 are small, and the regularized solution is in good agreement with the measured signal data. But xα2 is more sensitive to the change of regularization parameter, and the vertical part belongs to the underregularization state. In the horizontal part of L curve, the regularization parameter is relatively large, and the regularization error is dominant. With the increase of α, Kxα−y2 increases correspondingly, but xα2 almost does not change, so the horizontal part belongs to the overregularization state. Therefore, in order to balance underregularization and overregularization, the regularization parameter is selected at the corner of L curve uα (the angle between vertical part and horizontal part). Usually, people choose the point with the greatest curvature kα (in Figure 1(b)) on the L curve as the corner of the L curve.(5)uα=logKxα−y2,vα=logxα.Then the curvature function ofL curve with α as the parameter is(6)kα=u′v″−u″v′u′2+v′23/2,where ′ represents the derivation of α. Through the parametric expression of L curve, that is, the exact expressions of functions uα and vα, the maximum of curvature function can be directly calculated, and then the corresponding regularization parameters can be obtained.


## 2.2.1. The Generalized Cross-Validation (GCV) [13, 14]

The basic idea of cross-validation 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-validation depends on the specific ordering of the data, generalized cross-validation is invariant to the orthogonal transformation of the data vector y. The generalized cross function to be minimized in this method is defined as(4)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 [15] is that the GCV function can have a very flat minimum value, so the minimum value itself may be difficult to locate by a number. Another problem is that GCV sometimes mistakes noise for useful signals. GCV is quite effective for nonuniformity of square error and non-Gaussian error. However, if the errors are highly correlated, the method may not get satisfactory results.


## 2.2.2.L Curve Criterion [16–18]

Taking log-log as the scale, (Kxα−y2, xα2) forms a monotonously decreasing curve, as shown in Figure 1(a). Since the shape of this curve is similar to the letter “L,” it is called the L curve.Figure 1
(a) Schematic diagram ofL curve. (b) Curve of curvature function.
(a)(b)In the vertical part ofL curve, the regularization parameter and Kxα−y2 are small, and the regularized solution is in good agreement with the measured signal data. But xα2 is more sensitive to the change of regularization parameter, and the vertical part belongs to the underregularization state. In the horizontal part of L curve, the regularization parameter is relatively large, and the regularization error is dominant. With the increase of α, Kxα−y2 increases correspondingly, but xα2 almost does not change, so the horizontal part belongs to the overregularization state. Therefore, in order to balance underregularization and overregularization, the regularization parameter is selected at the corner of L curve uα (the angle between vertical part and horizontal part). Usually, people choose the point with the greatest curvature kα (in Figure 1(b)) on the L curve as the corner of the L curve.(5)uα=logKxα−y2,vα=logxα.Then the curvature function ofL curve with α as the parameter is(6)kα=u′v″−u″v′u′2+v′23/2,where ′ represents the derivation of α. Through the parametric expression of L curve, that is, the exact expressions of functions uα and vα, the maximum of curvature function can be directly calculated, and then the corresponding regularization parameters can be obtained.


## 2.3. The Basic Principle of Tikhonov Regularization Iteration Method in Frequency Domain

Many geophysical estimation problems are mathematically ill posed because they operate with insufficient data [19]. Regularization is a technique to make the estimation problem well posed by adding indirect constraints to the estimation model [20, 21]. The regularization method was first proposed by Tikhonov [22]. It has become an indispensable part of the inverse problem theory and has been widely used in geophysical problems [23–26].The relationship between gravity fieldTx,y and its vertical derivative Dmx,y (m represents the order) in the wavenumber domain is(7)Dmu,v=2πu2+v2mTu,v.Dmu,v and Tu,v are the spectra of the Fourier transform of Dmx,y and Tx,y. u and v are wavenumbers in the x and y directions, respectively. φm=2πu2+v2m is the vertical derivative operator. Due to the noise amplification characteristic of φm, the result Dmu,v is unstable, and the higher the derivative order is, the more unstable Dmu,v will be.For the instability problem of formula (7), Tikhonov regularization is a widely used method. Solve a minimized regularization functional; from formula (3),(8)minφ−m⋅Dmu,v−Tu,v22+αDmu,v22.2 represents the L2 norm. α is the regular parameter, used to balance the instability and smoothness. The regular operator in the article is selected according to the form of the derivative factor of the gravity field. The solution of the above minimization problem is(9)Dmαu,v=φm11+αφ2mTu,v,where Dmαu,v represents the conversion result of the regular vertical derivative. According to formula (9), Tikhonov regular vertical derivative conversion operator can be divided into two parts: conventional vertical derivative conversion operator part φm and Tikhonov regular low-pass filtering function part 1/1+αφ2m. It is equivalent to multiplying the conventional vertical derivative conversion operator by a regular low-pass filter. Then the approximation of Dmu,v (first-order approximation spectrum) after regular low-pass filtering can be expressed as(10)Dm1u,v=φm11+αφ2mTu,v.SubstituteDm1u,v into formula (7) to obtain an iterative approximate spectrum of T:(11)T1u,v=11+αφ2mTu,v.Then,Tu,v and T1u,v correct Dm1u,v to obtain the second iteration approximate spectrum of Dmu,v:(12)Dm2u,v=Dm1u,v+φm11+αφ2mTu,v−T1u,v.Repeat the above steps to obtain thenth-iteration approximate spectrum of Dmu,v:(13)Dmnu,v=Dmn−1u,v+φm11+αφ2mTu,v−Tn−1u,v.It can be seen from formula (7) that(14)11+αφ2mDmn−1u,v=φm1+αφ2mTn−1.Substitute it into formula (13):(15)Dmnu,v=φm1+αφ2mTu,v+1−11+αφ2mDmn−1u,v.It is easy to obtain the iterated general formula by mathematical induction:(16)Dmnu,v=1+1−11+αφ2m+⋯+1−11+αφ2mn−1φm1+αφ2mTu,v.Obviously,1+1−1/1+αφ2m+⋯+1−1/1+αφ2mn−1 is an equivalency sequence, with a common ratio of 0≤1−1/1+αφ2m<1; then,(17)Dmnu,v=1−1−11+αφ2mnφmTu,v.When the time of the iteration is infinite,(18)Dmu,v=limn⟶∞Dmnu,v=φmTu,v.The above formula shows that the Tikhonov regularization iteration method proposed in this paper can converge to the theoretical solution for high-order derivative conversion. When the order spectral difference is less than the specified iteration cut-off constantε or the number of the iterations reaches the specified iteration number, the iteration stops, where ε is a small positive number. The choice of regularization weight is to use the curvature of the L curve.


## 2.4. The Convergence and Filtering Properties of Tikhonov Regularization Iterative Method

Figures2 and 3 are the filtering characteristics of Tikhonov regularized iterative operators for the first vertical derivative and second vertical derivative, respectively. Figures 2(a) and 3(a) are the filter response characteristic curves of Tikhonov regularized iteration operator of the first vertical derivative and the second vertical derivative with different numbers of iterations, respectively.Figure 2
Filtering characteristics of Tikhonov regularized iterative operator for first vertical derivative. (a)α = 0.2. (b) The number of iterations is 3.
(a)(b)Figure 3
Filtering characteristics of Tikhonov regularized iterative operator for second vertical derivative. (a)α = 0.01. (b) The number of iterations is 3.
(a)(b)ω in the x-coordinate is the radial circular frequency. The value of regularization parameter of first vertical derivative α is 0.2, and α of the second vertical derivative is 0.01. By comparison and the filtering characteristics of the theoretical vertical derivative operators φm and Tikhonov regularization vertical derivative operators, it can be obtained that the theoretical vertical derivative factor rises sharply with the increase of frequency, and the higher the derivative order is, the stronger the amplification of the high-frequency components will be. Tikhonov regularization operator approximates the theoretical vertical derivative operator in low-frequency band, which ensures the accuracy of derivative conversion of useful signals in low frequency, but can suppress the noise in the high frequency. As the number of iterations increases, Tikhonov regularized iterative operators gradually approach the theoretical vertical derivative operators.In order to analyze the influence of the value of regularization parameterα on the filtering characteristics of Tikhonov regularization operator, Figures 2(b) and 3(b) are the filtering response characteristics of the Tikhonov regularized iterative operator of the first and second vertical derivatives when the number of iterations is 3, respectively. It can be seen that when the number of iterations is constant, with the increase of the order, the suppression of Tikhonov regularized iterative operator in high-frequency components is gradually enhanced, which ensures the stability of the high-order derivative’s calculation. The smaller the value of regularization parameter α, the weaker the suppression effect of Tikhonov regularization operator in high frequency. When the value is 0, it is the theoretical derivative operator. When the value of regularization parameter α is too large, it will suppress the low frequency, which will affect the filtering effect.It can be seen that the Tikhonov regularization iteration operator proposed in this paper can accurately approximate the theoretical derivative operator in the low frequency and effectively suppress the high-frequency component. Therefore, this method has the advantages of amplitude maintenance and stability.


## 3. Theoretical Model Test

In this paper, the gravity anomaly generated by a 2D vertical cylinder is used for the numerical tests. The model body is 2.0 km wide, and its centre position is 5.0 km. The top depth is 2.0 km, and the bottom depth is 3.0 km. The residual density is 1.0 g/cm3. Its location and the gravity anomaly Δg are shown in Figure 4(a). Figures 4(b) and 4(c) are the theoretical second vertical derivative and third vertical derivative, respectively. The random noise with 0.1% gravity anomaly amplitude is added. Figures 5 and 6 are the second and the third vertical derivatives calculated by the conventional FFT method and Tikhonov regularization iteration method in the case of noise. Figures 7 and 8 are regularization parameters calculated by the curvature function of L curve. The regularization parameters of the second and the third vertical derivatives are 114 and 26.8, respectively.Figure 4
Gravity anomalies generated by the model and their related theoretical derivatives. (a) The gravity anomaly and model location. (b) The second vertical derivative. (c) The third vertical derivative.
(a)(b)(c)Figure 5
The calculation of the vertical second derivative. (a) Conventional FFT. (b) Tikhonov regularization iterative method.
(a)(b)Figure 6
The calculation of the vertical third derivative. (a) Conventional FFT (b) Tikhonov regularization iterative method.
(a)(b)Figure 7
Curvature ofL-curve of the second vertical second derivative.Figure 8
Curvature ofL-curve of the third vertical second derivative.From Figures5 and 6, due to the interference of noise, the results of the vertical derivative calculated by conventional FFT transform are highly volatile, and the higher the order of the derivative is, the worse the stability is. It is impossible to identify the useful anomalies, which is due to the amplification effect of theoretical derivative operator in high-frequency components. The Tikhonov regularization iterative method is stable in calculating the second and third vertical derivatives and has a strong ability to suppress noise. This is due to the fact that, because of the filtering characteristics of the regularized iterative method, it can better approximate the theoretical factor in the low-frequency band and can suppress noise better in the high-frequency band.By comparing the second and third vertical derivatives of the theory in Figures4(b) and 4(c), it can be seen that, with the increase of the derivative’s order, the influence of high-frequency amplification becomes greater. From Figures 5(b) and 6(b), the derivatives calculated by Tikhonov regularization iteration method are of the same order of magnitude as the theoretical vertical derivatives (Figures 4(b) and 4(c)). There is no obvious oscillation in the third vertical derivative, which shows the accuracy and stability of Tikhonov regularization iteration method. At the same time, it can be seen that, with the increase of the calculated derivative order, the corresponding effect between the zero-point position and the boundary of the geological body becomes better and better.In order to test the application effect of Tikhonov regularization iteration method on plane data, the combined cuboid models are adopted in this paper. The model parameters are shown in Table1. We calculate the second vertical derivative of the gravity anomaly Δg data of the model. In order to illustrate the noise suppression effect of this method, add Gaussian white noise with a mean value of 0 mGal and a standard deviation of 0.03 mGal to the theoretical gravity anomaly of the model. The gravity anomaly and the noise-added gravity anomaly are shown in Figures 9(a) and 9(b).Table 1
The model parameters in Figure9.
Model numberTop/bottom (km)x direction boundary (km)y direction boundary (km)Residual density (103 kg/m3)A0.5/1.35.0/7.02.0/8.00.3B0.3/0.52.0/4.02.0/4.00.1C0.3/0.43.0/4.07.0/8.00.1D0.2/0.38.0/9.05.0/6.00.1Figure 9
The plane position and gravity anomaly contour map of the theoretical model. (a) Gravity anomaly of theoretical model. (b) Gravity anomaly after adding random noise to Figure9(a).
(a)(b)Figure10 shows the results of three methods for calculating the second vertical derivative of the combined models with noise (Figure 9(b)), where Figure 10(a) is the contour map of the theoretical second vertical derivative in the case of no noise.Figure 10
The second vertical derivative results. (a) The theoretical second vertical derivative. (b) The conventional FFT method after upward extension. (c) The ISVD method. (d) The Tikhonov regularization iteration method.
(a)(b)(c)(d)The upward continuation method can suppress the noise. Generally, in order to ensure the stability of the calculation results, the second vertical derivative is calculated after the upward continuation of the gravity data. Figure10(b) shows the second vertical derivative calculated by the conventional FFT method after the upward continuation.Figure10(b) is the second vertical derivative calculated by the ISVD method. Figure 10(d) shows the results by the method proposed in the paper, and the regularization parameter is 0.001. Table 2 shows the comparison between the results calculated by the three methods and the theoretical second vertical derivative results.Table 2
Comparison of the three methods.
MethodMaximum (g.u./km2)Minimum (g.u./km2)Average value (g.u./km2)Root mean square error (g.u./km2)Theoretical value55.00−29.760.31Conventional FFT42.95−20.310.294.50ISVD54.88−30.820.041.62Method in this paper57.61−31.770.331.29It can be seen from Figure10 that the three methods can calculate the second vertical derivative of the model stably. Among them, the zero-value line of the second vertical derivative calculated by conventional FFT can roughly identify the boundary of the model body. However, if the upward continuation continues, the noise will be suppressed more strongly, resulting in the reduction of high-frequency signal, and the boundary morphology of small-scale shallow geological bodies cannot be identified.It can be seen from Table2 that the maximum and minimum of ISVD are the closest to the theoretical second vertical derivative among the three methods, and the recognition effect of ISVD on the shallow model body is also good. Since the error of the method is cumulative in calculating the vertical derivative of each order, the error increases with the increase of the derivative order. As can be seen from Figure 10(c), the noise is relatively large, so denoising must be carried out before calculating the vertical derivative with this method.It can be seen from Table2 that the root mean square error of the method proposed in this paper is the smallest, and the zero-value line of the second vertical derivative coincides well with the boundary of the model body. This method has a strong ability of noise suppression and the proper selection of regularization parameters can balance the noise and useful signals.Through the comparison of the methods in Figure10, it can be concluded that the method based on Tikhonov regularization iteration method to calculate the vertical derivative of gravity data has a strong ability to suppress noise and the results are stable. The zero-value line of the second-order vertical derivative can better correspond to the boundary of the model body, so this method can be applied to the actual data processing.


## 4. Real Data

### 4.1. Geological Background

In order to test the effect of Tikhonov regularization iteration method on the processing of actual data, the second derivative of the gravity data of Qinghai-Tibet Plateau was calculated. Figure11 shows the location of the study area and its structural divisions. The Qinghai-Tibet Plateau is composed of many terranes, and each terrane presents a shape that is long from east to west and narrow from north to south. The lithostratigraphy and organisms also have some East-West continuity. It can be seen from Figure 11 that the suture zone is the contact zone of each block, that is, the boundary of each block. As mentioned earlier, the zero-value line of the vertical derivative can identify the structural boundary of the geological body. Using the method proposed in this paper, the second vertical derivative of gravity data in the study area is calculated, and the zero-value line is used to detect the fault structure in the study area.Figure 11
Research area and the tectonic division. 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. The Calculation of the Second Vertical Derivative Using Tikhonov Regularization Iteration Method

The location of the study area is 74–105° E, 26–40° N, which is the main part of the Qinghai-Tibet Plateau. The gravity anomaly results in the study area are shown in Figure12. It can be seen that the gravity value in the middle of the study area is the lowest and in the surrounding edge the value is high. There is an obvious anomalous gradient zone at the low gravity anomaly boundary, which corresponds to the structures in Figure 11.Figure 12
Gravity anomaly in the study area.In this paper, the second vertical derivative of gravity anomaly in the study area is calculated using Tikhonov regularization iteration method. The regularization coefficient is 0.1 and the number of iterations is 5. The results are shown in Figure13 and it can be seen that the position of the zero-value line of the second vertical derivative is in good agreement with the position of the gradient zone around the low density of the original gravity anomaly.Figure 13
The second vertical derivative of the study area.In order to make a clearer comparison, the zero-value line of the second vertical derivative is extracted, and the structural map is superimposed with it, as shown in Figure14. As can be seen from Figure 14, each suture zone can also be well reflected from the zero-value line. Although Jinsha River suture in the middle part of the study area is not continuously reflected, it can be outlined by scattered zero-value lines. For the faults, Altyn fault in the north and Jiali fault in the south can be identified from the zero-value line of the second vertical derivative. All of the above shows that the vertical derivative of gravity calculated by Tikhonov regularization iteration method is reliable and has good application effect on the actual data.Figure 14
Zero-value line of the second vertical derivative anomaly of the study area.


## 4.1. Geological Background

In order to test the effect of Tikhonov regularization iteration method on the processing of actual data, the second derivative of the gravity data of Qinghai-Tibet Plateau was calculated. Figure11 shows the location of the study area and its structural divisions. The Qinghai-Tibet Plateau is composed of many terranes, and each terrane presents a shape that is long from east to west and narrow from north to south. The lithostratigraphy and organisms also have some East-West continuity. It can be seen from Figure 11 that the suture zone is the contact zone of each block, that is, the boundary of each block. As mentioned earlier, the zero-value line of the vertical derivative can identify the structural boundary of the geological body. Using the method proposed in this paper, the second vertical derivative of gravity data in the study area is calculated, and the zero-value line is used to detect the fault structure in the study area.Figure 11
Research area and the tectonic division. 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. The Calculation of the Second Vertical Derivative Using Tikhonov Regularization Iteration Method

The location of the study area is 74–105° E, 26–40° N, which is the main part of the Qinghai-Tibet Plateau. The gravity anomaly results in the study area are shown in Figure12. It can be seen that the gravity value in the middle of the study area is the lowest and in the surrounding edge the value is high. There is an obvious anomalous gradient zone at the low gravity anomaly boundary, which corresponds to the structures in Figure 11.Figure 12
Gravity anomaly in the study area.In this paper, the second vertical derivative of gravity anomaly in the study area is calculated using Tikhonov regularization iteration method. The regularization coefficient is 0.1 and the number of iterations is 5. The results are shown in Figure13 and it can be seen that the position of the zero-value line of the second vertical derivative is in good agreement with the position of the gradient zone around the low density of the original gravity anomaly.Figure 13
The second vertical derivative of the study area.In order to make a clearer comparison, the zero-value line of the second vertical derivative is extracted, and the structural map is superimposed with it, as shown in Figure14. As can be seen from Figure 14, each suture zone can also be well reflected from the zero-value line. Although Jinsha River suture in the middle part of the study area is not continuously reflected, it can be outlined by scattered zero-value lines. For the faults, Altyn fault in the north and Jiali fault in the south can be identified from the zero-value line of the second vertical derivative. All of the above shows that the vertical derivative of gravity calculated by Tikhonov regularization iteration method is reliable and has good application effect on the actual data.Figure 14
Zero-value line of the second vertical derivative anomaly of the study area.


## 5. Conclusion

In order to solve the problem that the calculation of gravity high-order vertical derivative is unstable in frequency domain and sensitive to noise, a new method based on Tikhonov regularization iteration method is proposed to calculate the high-order vertical derivative of gravity field. The reliability of vertical derivative calculation is verified by model test. Finally, the structures of the study area are obtained according to the zero-value line of the second vertical derivative of Qinghai-Tibet Plateau. We draw the following conclusions:(1)
The calculation of gravity vertical derivative is often carried out in the frequency domain, but the vertical derivative operator has the amplification effect in the high-frequency domain, which leads to instability and sensitivity to noise. The higher the derivative order is, the more unstable the calculation result is.(2)
Through the analysis of the convergence and filtering characteristics of Tikhonov regularization iterative method, it can be seen that the method can accurately approximate the theoretical derivative operator in the low-frequency band and it can effectively suppress the high-frequency component. Therefore, the method can suppress the noise, and the calculated vertical derivative results are stable and close to the theoretical vertical derivative.(3)
In the algorithm, the smaller the regularization parameter is, the smaller the suppression effect on high-frequency components is, and the calculation is greatly affected by noise. The larger the regularization parameter is, the stronger the suppression effect is in the high-frequency components, and the low-frequency components will also be suppressed. Therefore, choosing the proper regularization parameter is the advance of the method to calculate the vertical derivative accurately.(4)
As an effective method to calculate the high-order vertical derivative of gravity, the Tikhonov regularization iteration method is used to calculate the gravity second vertical derivative of Qinghai-Tibet Plateau. The calculated results are in good agreement with the geological data.



---

*Source: 8818552-2021-05-08.xml*