Application of the normalized largest eigenvalue of structure tensor in the interpretation of potential field tensor data

Obtaining horizontal edges and the buried depths of geological bodies, using potential field tensor data directly is an outstanding question. The largest eigenvalue of the structure tensor is one of the commonly used edge detectors for delineating the horizontal edges without depth information of the potential field tensor data. In this study, we presented a normalized largest eigenvalue of structure tensor method based on the normalized downward continuation (NDC) to invert the source location parameters without any priori information. To improve the stability and accuracy of the NDC calculation, the Chebyshev–Pade´ approximation downward continuation method was introduced to obtain the potential field data on different depth levels. The new approach was tested on various models data with and without noise, which validated that it can simultaneously obtain the horizontal edges and the buried depths of the geological bodies. The satisfactory results demonstrated that the normalized largest eigenvalue of structure tensor can describe the locations of geological sources and decrease the noise interference magnified by the downward continuation. Finally, the method was applied to the gravity data over the Humble salt dome in USA, and the near-bottom magnetic data over the Southwest Indian Ridge. The results show a good correspondence to the results of previous work.


Introduction
The source parameters estimation of potential field tensor data is an important part for geological interpretation, which can provide the horizontal position and the buried depth of the geological bodies as an a priori information for potential field inversion (Li and Oldenburg 1998). The normalized full gradient (NFG) method proposed by the Russian Geophysical School Pašteka 2009, 2019) can be directly used to obtain (2011) applied different normalization factors on analytical signals modulus and the downward continued potential field itself, and called this generalized method as normalized downward continuation (NDC). Based on the generalized NDC method, some authors have applied it on the local wavenumber (Ma et al. 2014), total horizontal derivatives , directional analytic signals (Zhou 2015) and directional total horizontal derivatives (Zhou et al. 2017) to obtain the horizontal edges and the buried depth simultaneously.
The structure tensor is one of the image processing techniques and presents a local orientation in an n-dimensional space (Weickert 1999a, b). The largest eigenvalues of the structure tensor of the potential field data have been used to delineate the edges of the geological bodies (Sertcelik and Kafadar 2012;Oruç et al. 2013;Yuan et al. 2014). In contrast to traditional derivativebased edge detection methods, the structure tensor has a property of which can reduce noise in the data while enhancing discontinuity boundaries, due to a Gaussian envelop within it. Based on the advantages of the largest eigenvalue of the structure tensor, in this paper, we apply the NDC on the largest eigenvalue of the structure tensor of the downward continued field to estimate the depth and the horizontal edges. In order to increase the stability of the NDC process, the Chebyshev-Pade´ approximation downward continuation method developed by Zhou et al. (2018) was brought to compute each downward continuation level field. The synthetic models with and without noise and the real measured potential field data are utilized to validate the effectiveness of this new method.

NDC of largest eigenvalue of structure tensor
Structure tensor of potential field tensor data Potential field gradient tensor data are the second-order space derivatives of potential field f (gravitational field or magnetic field) in the three orthogonal directions x , y and z . The potential field gradient tensor matrix is The original structure tensor consists of a Gaussian envelope and horizontal gradient tensor of the potential field data (Sertcelik and Kafadar, 2012). The expression of the original structure tensor matrix is x 2 δx 2 + y 2 δy 2 , δx and δy are the standard deviations of Gaussian envelope in x and y directions. The homogeneous characteristic equation for 2D tensor M is The largest eigenvalue of matrix M is Sertcelik and Kafadar (2012) pointed that the largest eigenvalue can locate the edges of geological bodies. However, the balancing ability of equalizing the edge signal amplitude of large and small anomalies is weak. The larger standard deviation δx and δy can enhance the balancing ability, but reduce the resolution of the identified edges.

NDC of the largest eigenvalue
The normalized full gradient (NFG) method (Elysseieva and Pasteka 2009;Zeng et al. 2002) is the normalization of the analytic signal modulus at different downward continuation levels. Fedi and Florio (2011) applied the normalization on the analytic signal and on the downward continuation potential field itself, called this generalized method as normalized downward continuation (NDC). The expression of NDC applied to the largest eigenvalue of the structure tensor matrix can be expressed as: where x, y, z is the edge detector of the structure tensor at point (x, y, z); N (z) is the normalization function.
Here, arithmetic mean, median and geometric mean are used as the normalization function, shown as: where M is the number of all calculation points. (3) (4) During the processing, there are many issues that needed to be overcome. The downward continuation plays the key role in the NDC method. The accuracy of the NDC method is determined by the stability of downward continuation method directly. Therefore, it is necessary to use a stable downward continuation in the calculation process. Many new stable algorithms have been introduced to implement the downward continuation method (Fedi and Florio, 2002;Cooper, 2004;Ma et al., 2013;Zeng et al., 2013Zeng et al., , 2014Zhang et al., 2013;Zhou et al., 2018). Zhou et al. (2018) has compared the errors between the function exp(x) and its different approximation functions, including Taylor series, Chebyshev approximation, Pade´ approximation, and Chebyshev-Pade´ approximation. Comparison of results indicate that downward continuation based on Chebyshev-Pade´ approximation can obtain a more precise result. Therefore, in this study, the Chebyshev-Pade´ approximation downward continuation method was chosen for the NDC of the largest eigenvalue of the structure tensor.
In the frequency domain, the downward continuation can be denoted as: where the downward continuation operator e �hω γ can be approximated by the (three and two) terms of Chebyshev-Pade´ as follows (Zhou et al. 2018): Therefore, the expression of the (three and two) terms of Chebyshev-Pade´ approximation downward continuation is where f h and f are the gravity and magnetic data at two observation heights separated by a vertical distance h , respectively. f h and f denotes the Fourier transform of f h and f , ω x and ω y are the wavenumbers x-and y-direction, and ω γ = ω 2 x + ω 2 y is the radial wavenumber. h is a positive number for the downward continuation distance and vice versa for the upward continuation. Zhou et al. (2018) have pointed out that Chebyshev-Pade´ approximation can give a more precise approximation in (9) the low frequency and it can suppress the high frequency. During the calculation of downward continuation, the short wavelength (high frequency) components can be suppressed by Chebyshev-Pade´ approximation. Therefore, we can obtain a stability and noise reduced potential field gradient data at different depths by using the Chebyshev-Pade´ downward continuation. The residual noise of the potential field gradient data at different depths can be further removed by the Gaussian envelop when using the NDC of the largest eigenvalue of structure tensor to estimate the geology source. Therefore, the calculation procedures of the NDC of the largest eigenvalue of structure tensor are: 1. Calculate the potential field gradient data from the observed potential field data, or obtain from the real measurement. 2. Set a maximum depth H (the depth that we want to know the maximum range of geologic source) for downward continuation, and divide the underground into n layers (i = 1, 2, 3, …, n). 3. alculate potential field gradient data at the different depth levels by using the Chebyshev-Pade´ approximation downward continuation method.
4. Calculate the NDC of the largest eigenvalue of structure tensor of potential field data according to Eqs. 2, 3, 4, 5.

Application to synthetic model data
Model 1: prism model In this section, we construct a prism gravity model with the center depth of 10 m to validate the effectiveness of the new proposed method. The thickness is 5 m and the contrasted density is 1 g/cm 3 . The 3D view of the synthetic model and the synthetic gravity anomaly and gravity gradient components f zx and f zy are shown in Fig. 1. The Chebyshev-Pade´ approximation downward continuation method is utilized to calculate different depth-level fields with the interval of 1 m, and the maximum depth is set as 25 m. In each depth level, arithmetic mean, median value and geometric mean are used as the normalization function to normalize the largest eigenvalue of the structure tensor. In order to compare the results conveniently, a cross section (the white line in Fig. 1b) is intercepted from the 3D result to display the character of the estimated depth results. Figure 2a-c shows the cross-sectional results normalized by arithmetic mean, median and geometric mean with δx = δy = 0.01 . All of their maximum values (the white dot) can display the edge position and the center depth of the prism at 10 m. Comparing these results, the maximum values normalized by the median display the center depth of the prism more clearly, shown in Fig. 2b. Figure 2d shows the 3D view of the median normalized result at different depths. We can see that at the depth of 10 m, the median normalized result can display the horizontal edge of the prism clearly. It demonstrates that we can obtain the source edge position and the buried depth information simultaneously by our new proposed method.

Model 2: prism model with noise
To capture the effect of noise in the new proposed method, we add 1% Gaussian noise to the synthetic gravity gradient anomalies in Fig. 1c, d. The noisy gravity gradients f zx and f zy are shown in Fig. 3. Figure 4 shows the cross section estimated depth results by the NDC of the largest eigenvalue of the structure tensor normalized by arithmetic mean, median and geometric mean with different standard deviations δx and δy . We can see that all the results can display the depth information with the maximum values. As we have discussed in the theoretical part, the effect of noise can be suppressed by Chebyshev-Padé downward continuation. However, there are still some residual noise in the calculated potential field gradient data. The large standard deviation δx and δy can remove the residual noise disturbance more efficiently, and make the normalized results more clearly. However, with the increasing of the standard deviation δx and δy , the resolution of the edge information of the prism have been reduced and may bring an error to the result. From Fig. 4i, when standard deviation δx = 10 and δy = 10 , the estimated depth by geometric mean is 11 m, with an error of 10%. Therefore, a suitable standard deviation δx and δy can balance denoising ability and the resolution, and obtain a better result. However, the choice of a suitable standard deviation δx and δy depend on the noise level of the original measured potential data. A large noise level need choosing a large standard deviation, vice versa. In this model, we can see that Fig. 4d, e with the standard deviation δx = δy = 1 get a better result than others, which can reduce the noise disturbance and retain the resolution of the identified horizontal edges. Figure 5 displays the horizontal and vertical profiles of the distribution of the largest eigenvalues across, and along the edge of the model 1. From Fig. 5a, we can see that the normalized results by arithmetic mean, median and geometric mean all can demonstrate the horizontal edges model well. However, from Fig. 5b, we can see that the extent of the maximum values of the result normalized by median are corresponding to the vertical edges of the causative body. Therefore, we consider the normalized result by median can get a better result normalization factors arithmetic mean and geometric mean.

Model 3: complex model with two prisms
In order to test the application ability, a complex gravity model was constructed which contains two vertical-sided prisms with the center depths are 10 km and 15 km. The thickness of both prisms is 5 km. The contrasted density of the prisms is 1 g/cm 3 and − 2 g/ cm 3 , respectively. The 3D view of the synthetic model and the synthetic gravity anomaly and gravity gradient anomalies f zx and f zy are shown in Fig. 6. The downward continuation is utilized to calculate different depthlevel fields with an interval of 0.5 km, and the maximum depth is set as 25 km. Here, based on the results of model 1 and model 2, only the median was chosen as the normalization function to normalize the largest eigenvalue of the structure tensor in each depth level on this model. Figure 7c shows the cross section (white line in Fig. 6b) estimated depth results, with the maximum values of the normalized results by median are 10 and 15 km for prism 1 and prism 2, respectively. Figure 7d shows the 3D view of the median value normalized result at different depths. We can see that at the depth of 10 and 15 km, the median value normalized result can demonstrate the horizontal edge of the prism 1 and prism 2 clearly. Figure 7e displays the identified horizontal edges of the two prisms in plan view, which was extracted from Fig. 7d by the method proposed by Blakley and Simpson (1986), which compared with the eight nearest neighbors of the identified horizontal edge grid in four directions (along the row, column, and both diagonals) to see if a maximum is present. We can see that the newly proposed method can obtain the source edge position and the buried depth information simultaneously.
We should reduce the magnetic anomaly to the pole before applying the new method to the magnetic anomaly to detect the horizontal edges and estimate the depth, as the horizontal derivatives of the magnetic anomaly are sensitive to magnetization direction ).

Gravity anomaly of Humble salt dome
To demonstrate the real application effect, the new method is applied to the residual Bouguer gravity anomaly caused by the Humble salt dome near Houston, Texas (Nettleton 1976), which was digitized by Shaw and Agarwal (Shaw and Agrawal 1997). Figure 8a shows the Bouguer anomaly where the grid interval is 100 m. The horizontal gravity gradients f zx and f zy are calculated numerically from gravity data using the method proposed by Minkus and Hinojosa (2001), shown in Fig. 8b, c. The interval of the depth level is 0.3 km, and the maximum depth is set as 9 km. Figure 9 shows the depth estimated results of NDC of the largest eigenvalue of structure tensor by median. The estimated center depth of the salt dome is 4.2 km by the maximum value of the    Fig. 9b. Besides, several other researchers have used different methods to interpret this anomaly. Table 1 is the comparison of the new method with previous salt dome depth estimation results, which proves the depth information obtained through the new method is reasonable. Figure 9c shows the horizontal slices of the NDC results normalized by the median at different depths, and Fig. 9d is the horizontal slice at the buried depth level 4.2 km in Fig. 9c, which indicates the horizontal extent of the dome in E-W direction is about 4 km while in N-S direction is about 2.5 km.

Near-bottom magnetic anomaly of Longqi hydrothermal sulfides
Magnetic surveys have been established as an effective method to explore hydrothermal sulfides and study the internal structure of sulfide deposits (Tivey et al. 1993(Tivey et al. , 1996Szitkar et al. 2014a, b;Szitkar and Dyment 2015). We use the new method to a near-bottom magnetic data over the first observed active high-temperature hydrothermal vent field at the ultraslow spreading Southwest Indian Ridge (SWIR). The data were collected by the autonomous underwater vehicle ABE (Autonomous Benthic Explorer) on board R/V DaYangYiHao with a fixed driving depth 2625 m in Feb.-Mar. 2007(Zhu et al. 2010;    Fig. 8a is the cross section for Fig. 9 Tao et al. 2017). The measured line space is 250 m. The measured magnetic anomaly reveals a strong correlation with the topography. Zhu et al. (2010) pointed out that the entire local topography high is an excess magnetic volume. Figure 10 displays the measured topography, magnetic anomaly and the calculated two magnetic gradient components f zx and f zy . From Fig. 10a, the magnetic anomaly deviates from the local topography, which indicates that the root of this excess magnetic volume is tilt. Here, our new method was used to estimate the depth of the magmatic volume. In the downward continuation, the interval of the depth level is 0.05 km, and the maximum depth is set as 2 km. Figure 11 shows the depth estimated results of NDC of largest eigenvalue of structure tensor by the median value. We can see that the estimated center depth is between the 0.55 and 0.80 km by the maximum value of the normalized results by median below the measurement plane. When removing the depth of seawater, the estimated depth of the magmatic volume is between 0.404 and 0.583 km below the  seafloor. The estimated depth shows that magmatic volume is tilted, consistent with the interpretation result in Fig. 10a.

Conclusions
We have presented a source parameters estimation method based on the NDC of the largest eigenvalue of the structure tensor of the potential field tensor data, which can simultaneously locate the horizontal edges and the buried depth of the geological sources. The maximum value in 3D space can locate the position. Due to the use of a Gaussian envelop filter in the structure tensor, the new proposed method can reduce the noise effect that is magnified by the downward continuation. The method has been tested by synthetic data with and without noise, and satisfactory results are obtained. Finally, the method was applied to the real measured gravity data over the Humble salt dome in USA and the near-bottom magnetic data over the Southwest Indian Ridge successfully. The results show a good correspondence to the previous work results.