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.
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:where is the measured signal and is the real signal, is the uncertain noise, the kernel function is the response function, and 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) iswhere 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:
In the above formula, the regularization parameter mainly controls the relative size between the solution norm and the residual norm . 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 the L 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 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-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 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 [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, (, ) 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.
(a)
(b)
In the vertical part of L curve, the regularization parameter and are small, and the regularized solution is in good agreement with the measured signal data. But 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 , increases correspondingly, but 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 (the angle between vertical part and horizontal part). Usually, people choose the point with the greatest curvature (in Figure 1(b)) on the L curve as the corner of the L curve.
Then the curvature function of L curve with as the parameter iswhere ′ represents the derivation of . Through the parametric expression of L curve, that is, the exact expressions of functions and , 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 field and its vertical derivative (m represents the order) in the wavenumber domain is
and are the spectra of the Fourier transform of and . u and are wavenumbers in the x and y directions, respectively. is the vertical derivative operator. Due to the noise amplification characteristic of , the result is unstable, and the higher the derivative order is, the more unstable will be.
For the instability problem of formula (7), Tikhonov regularization is a widely used method. Solve a minimized regularization functional; from formula (3),
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 iswhere 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 and Tikhonov regular low-pass filtering function part . It is equivalent to multiplying the conventional vertical derivative conversion operator by a regular low-pass filter. Then the approximation of (first-order approximation spectrum) after regular low-pass filtering can be expressed as
Substitute into formula (7) to obtain an iterative approximate spectrum of T:
Then, and correct to obtain the second iteration approximate spectrum of :
Repeat the above steps to obtain the nth-iteration approximate spectrum of :
It can be seen from formula (7) that
Substitute it into formula (13):
It is easy to obtain the iterated general formula by mathematical induction:
Obviously, is an equivalency sequence, with a common ratio of ; then,
When the time of the iteration is infinite,
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
Figures 2 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.
(a)
(b)
(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 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 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.
(a)
(b)
(c)
(a)
(b)
(a)
(b)
From Figures 5 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 Figures 4(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 Table 1. We calculate the second vertical derivative of the gravity anomaly 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).
(a)
(b)
Figure 10 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.
(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. Figure 10(b) shows the second vertical derivative calculated by the conventional FFT method after the upward continuation.
Figure 10(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.
It can be seen from Figure 10 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 Table 2 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 Table 2 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 Figure 10, 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. Figure 11 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.
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 Figure 12. 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.
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 Figure 13 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.
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 Figure 14. 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.
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.
Data Availability
The data are downloaded from International Gravimetric Bureau (https://bgi.obs-mip.fr/data-products/grids-and-models/wgm2012-global-model/).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge the International Gravimetric Bureau for providing the gravity data. This work was supported by the National Natural Science Foundation of China (41904129)