Full-tensor gravity gradient eigenvector analysis for locating complex geological source positions

We develop an eigenvector analysis method to locate the centroids of geological structures of full tensor gravity 10 gradient (GGT) data. Although the boundary detection method for Bouguer gravity and CGGT (curvature gravity gradient tensor) has been widely discussed and applied, the source location method for GGT data remains an area of active research. In this paper, we first discuss the theoretical basis and physical meaning of the eigenvector analysis on GGT data, and then a new source centriod location method is derived. The proposed method uses eigenvector analysis to extract the source centroid information. The interference of multiple and overlapping sources and the parameter identification related with the 15 multiple scales of the GGT eigenvector analysis are presented in the theoretical and experimental sections. Finally, the proposed method is applied to synthetic and field data.


Introduction
A gravity survey can reflect the response of density contrasts in the subsurface, such as high-density mineral deposit, or lowdensity oil deposit, with respect to the host or country rock. The gravity gradient anomaly as defined here is an anomalous 20 response with respect to the background at the particular scale and magnitude of interest. For example, in geotechnical applications, a mass deficit due to a mine would result in a short-wavelength negative anomaly with respect to background geological signal of larger wavelength. Conversely, a mass excess due to an ore body would produce a positive anomaly imprinted on the background tectonic/geoid-based signals. Location of source positions is critical in the interpretation of gravity gradient data, which has been widely used in mineral exploration and resource surveys (e.g. Dransfield 1994;25 Mikhailov et al., 2007). Various edge detection methods for Bouguer gravity anomaly data, such as the derivative filter, analytic signal and others, have been proposed to delineate the outline of a geologic target and provide vital information for Nonlin. Processes Geophys. Discuss., doi:10.5194/npg-2016-75, 2017 Manuscript under review for journal Nonlin. Processes Geophys. Published: 13 January 2017 c Author(s) 2017. CC-BY 3.0 License. data interpretation (Gordell 1989;Wijns et. al., 2005;Li 2006;Cooper and Cowan 2006;Foks 2013;Ma and Li 2012;Phillips 2015;).
In recent years, gravity gradient tensor (GGT) measurement devices and methods have been widely researched (e.g. Pedersen and Rasmussen, 1990;Zhdanov et al., 2004;Fedi, et al., 2005;Dransfield, 2010;Beiki et al., 2011). Various platforms, such as airborne, ground, and satellites have been utilized for GGT data observation. Although, there are many 5 boundary delineation methods that perform well in the user coordinate system for potential fields, recent research (Li, 2015) shows that rotation of the coordinate system can reveal additional information in GGT data, especially when choosing the new coordinate system via eigenvector analysis. Pedersen and Rasmussen (1990) introduced eigenvector analysis on the gradient tensor of potential field data and derived the relationship between the eigenvalues and eigenvectors of the potential field source parameters. Beiki and Pedersen (2010) combined Euler deconvolution and the eigenvector analysis to estimate 10 the depth of the geological body. Sertcelik and Kafadar (2012) suggested that the eigenvalues of CGGT (curvature gravity gradient tensor) could be used to detect the edges and corners of a subsurface source, while Oruc et al., (2013) applied the CGGT eigenvalue analysis on a complex field data interpretation application. Zuo and Hu (2015) proposed the eigenvalue analysis method for GGT data incorporating signal normalization. In addition to these eigenvalue methods, Cevallos (2016) proved that using tidal tensor theory to analyze the equipotential surface can, to some extent, reach an equivalent result to the 15 eigenvalue analysis method.
Although the utility of eigenvalues of GGT data for boundary detection performs well on source boundary detection as discussed in previous papers, the corresponding eigenvector analysis for locating the source position has not been presented.
Identifying the source position has important utility not only in geologic applications, but also in environmental, time-lapse, and engineering applications as well. In this paper, we apply eigenvector analysis on GGT data to delineate the source 20 positions and boundaries. The relationship between the proposed eigenvector analysis method and the source parameters is developed in detail. Unlike previous Bouguer gravity anomaly boundary location methods that mainly indicate edges and CGGT methods that detect the corners, the proposed method is designed to locate the centroid of sources while delineating the boundary of them, tying edges to sources, functioning even when sources overlap. In the experimental section, both synthetic and field data sets are used to illustrate the effectiveness of the method. 25

Eigenvector analysis of GGT
The gradient of the gravitational field, T, can be written in the form (Eq. (1)): Nonlin. Processes Geophys. Discuss., doi:10.5194/npg-2016-75, 2017 Manuscript under review for journal Nonlin. Processes Geophys. Published: 13 January 2017 c Author(s) 2017. CC-BY 3.0 License. y, z; j = x, y, z; ) denotes the elements of the symmetric gravity gradient tensor. We assumed that the gradient tensor is decomposed as: is the eigenvector matrix and Λ=diag[λ 1 , λ 2 , λ 3 ] is the eigenvalue matrix. Eigenvalue and eigenvector decomposition has been described in many textbooks; we use the notation as described in Anton and Rorres (2000). Pedersen and Rasmussen (1990) suggested the eigenvalue decomposition algorithm 5 for gradient tensor data, and expressed the relationship between the eigenvalue decomposition and coordinate system rotation. Dransfield (1994) proposed a method that rotates the coordinate system around the z-axis with a rotation transform, and Cevallos (2016) introduced the expression of the eigenvector obtained from 3-D tensor data as Eq.
(2) is theoretically derived from the classical matrix eigenvector decomposition method, and the numerical calculated 10 results are equivalent. Pedersen and Rasmussen (1990) illustrated that the eigenvector v 1 (corresponding to the largest eigenvalue λ 1 ) display a relationship between the observer and the source. They deduced the relationship based the classical gravity gradient forward equation with a 3D point mass model. Their deduced result shows that the eigenvector v 1 indicates the direction from the observer to the point source ( v 1 = (x ', y', z ') R , where R is the distance between the point source and observation point). Although the practical geological structure is much complicated than a mass point, it can be considered 15 as a combination of many point mass cells, and the eigenvector of these sets of cells contains the boundary information of the total anomalous source. Beiki (2010) suggested that the eigenvalue decomposition could be considered as an equivalent point source method for a causative gravimetric body. Li (2015) presented theoretical and practical aspects related to the eigenvectors and rotations of the gravity gradient tensor. Li (2015) suggested that the rotation of the coordinate system aligned with the source direction will lead to the diagonalization of a gravity gradient matrix (in the absence of noise, the 20 gradient matrix is symmetric positive-semidefinite). The magnitude of Φ xy , Φ xz and Φ yz will be much smaller than the Nonlin. Processes Geophys. Discuss., doi:10.5194/npg-2016-75, 2017 Manuscript under review for journal Nonlin. Processes Geophys.
diagonal terms after the transforms. He also discussed the necessity of performing such a spatially-dependent coordinate rotation to accurately calculate curvature.
For every observation point, we assume that the eigenvector v 1 is unique, and the diagonalization can be applied for local GGT data matrix effectively. Physically, the diagonalization of the GGT data matrix can reduce the value of components Φ xy , Φ xz and Φ yz . The eigenvector decomposition on every single observation data point can extract the 5 correct source centroid position direction, and this procedure will not be influenced by complex geological structure and noise.
In this paper, we propose the source location method according to the position vector of the equivalent source with eigenvector analysis. The position vector contains the rotated angle φ and the horizontal vector r . As Fig. 1 shows, angle φ rotates along the z-axis with eigenvector decomposition. Physically, the essence of the eigenvector analysis is to decompose 10 GGT data into two parts -the eigenvector dataset related with the source position, and the eigenvalue data set is related with the anomalous intensity of the source. The eigenvector in every observation position ( v 1 (x ', y ', z ')' ) is a unit vector that is directly points to a source centroid, while the eigenvalue data sets contain magnitude information of the source. For example, changing a point source density can result in a different eigenvalue, but the eigenvector will not be changed. In this paper, we utilize the eigenvector v 1 to measure the source centroid position from every observation point in observation plane. The 15 angle φ is used to measure the position direction along the z-axis. It is the complementary angle between the vector v 1 that points to the source centroid and the vertical z-axis. In other words, for every observation data point, the eigenvector transform can be considered as a rotation from z-axis to z'-axis. r is the horizontal vector of the equivalent source to the observation point in a rotated horizontal plane, while r is the horizontal magnitude of the eigenvector in the horizontal plane. The rotated angle φ can be identified according to the eigenvector v 1 . 20 As Fig.1 shows, if the buried source is immediately below the observation point, then the angle φ will reach the maximum value ( π / 2 radians) and r will equal zero. Thus, the rotated angle φ indicates the actual source direction from the observation point. We propose to utilize the angle φ and vector r to indicate the centroid and the profile of the equivalent sources.
The eigenvector v 1 can be calculated from the gravity tensor components directly (Eq. (2)); this eigenvector should be unit 25 vectors. However due to systematic and random errors, usually eigenvectors v 1 calculated in this manner are not unit vectors. Therefore, we apply a balance operator n(⋅) on v 1 to reduce it to a unit vector. According to Eq. (2) and above analysis, we locate the source centroid and the outlines by using tanφ and r , as in Eq. (3) Nonlin. Processes Geophys. Discuss., doi:10.5194/npg-2016Discuss., doi:10.5194/npg- -75, 2017 Manuscript under review for journal Nonlin. Processes Geophys. Published: 13 January 2017 c Author(s) 2017. CC-BY 3.0 License.
Briefly, tanφ can be also expressed as tanφ =v 1 . From the view of the eigenvector, if there exists only a single point source, r will equal 0 at the position of source centroid in the observation plane. Since v 1 is a unit vector, z' (Eq. (2)) will nearly equal 1 and tanφ → ∞ . As is known, in practical data processing and interpretation, the sources are inhomogeneously distributed in all directions, and multiple sources commonly overlap one another, so the modulus of vector 5 r when directly over the source will not exactly equal 0. tanφ will present a relative large value above the centroid of sources on observation plane. With the considering of negative and positive anomaly in data, we add the vertical gravity data Φ z in and define GTA (Gradient Tensor Angle) as expressed in Eq. (4): where Φ z = ∂U ∂z . Φ z is the spatial change rate of gravitational acceleration in the vertical direction. Comparing with 10 Bouguer gravity data, Φ z can present a more significant anomalous value at the observation position above a buried source, while the spectral power of Φ z contains more high frequencies anomalous information. Although Φ z commonly is not directly measured in GGT survey, it can be calculated by using the measured vertical gravity gradient ( Φ zz ) and regional ground gravity data (Dransfield, 2010).
Small sources, or inhomogeneous distribution of a source can present some small local tanφ peak values. ⋅ ⎡ ⎤ β β is design for 15 filtering out these local small values. The filter threshold parameter β can be automatically identified, as in Eq.(5) where According to Eq. (4), GTA value dramatically decreases with the distance increasing from the centroid position, and it will be close to zero value at the source boundary position. So with drawing the contours on the GTA map, it is possible to delineate the boundary of sources.
The value of tanφ will decay dramatically with an increasing |r|. In the multiple source condition, the GTA value of at every source centroid position will be influenced by others sources. It will not reach an infinitely large value, but it still relates with 5 the mess of the source strongly. To proving this point, we deduce the following analysis on the basis of Eq. (4) Assuming there is a point source s 1 and its horizontal direction vector is r s 1 ( r s 1 → 0 ), we add another similar point source According to Eq. (7), in the overlapping condition the tanφ of the two sources does not interfere with each other dramatically, the tanφ ratio of them still is equal to value m approximately. Both the local maximum values of weak 20 anomaly sources and strong anomaly sources can be easily located.

Synthetic Experiments
In this section, the performance of the proposed method is tested with synthetic data. The selection of the scalar parameter ( β ) is discussed in detail in the field example.
Considering the multiple-source scenario, we design a synthetic model that contains nine sources which are distributed in 5 different depths, as shown in Fig.2. The depth to the top of the nine sources ranges from 1km to 2.6km, each with equal depth extents of 0.2km. The density contrast of each source is 0.3 g/cm 3 . There are 200x200 observation points involved in this experiment, and both x-y direction lateral extent of this synthetic model are set as 100m. The data are simulated as a fulltensor gradiometer response observed at ground.
The vertical gravity data Φ z and the tensor component Φ zz are shown in Fig.3 with a 0.1km spaced observation grid. We 10 calculate GTA according to Eq. (4) with β =10 (Fig. 4). The value of β is identified according to Eq. (5).
As Fig.4 shows, although there are multiple sources and the anomalies are weak for the deep sources, the result of GTA is not obviously influenced by this complexity. The GTA value at the centroid of the 1 st (shallowest) source reaches a value of 0.6098, while the GTA value at the centroid of the 9 th (deepest) source is 0.1443. This GTA map shows all of the centroids of these sources with high precision. The boundaries of the anomaly sources, including the boundaries of the deeper buried 15 sources, are all delineated correctly in the contour maps.
To test the stability and robustness of the GTA map with noise, we add Gaussian noise with a standard deviation equal to 30% of the max magnitude of the tensor components to all of the gravity gradient components. The data added noise and the output GTA map are shown in Fig.5.
As Fig. 5 shows, there is obvious noise in the tensor data (Fig.5), but the source location results (Fig. 6(a)) are the same as 20 the GTA map produced by the clean data ( Fig. 4(a)). The noise interferences that show in the contour map ( Fig. 6 (b)) mainly come from the product of the vertical gravity data ( Φ z )(Eq. (4)), which are not amplified in the GTA map calculation.
In practice, the sources are distributed in various depths and may overlap each other. Here we design another synthetic model to test the performance of the GTA method in this condition. The deeper source (label A) is buried at a depth of 2km with depth extent 1km. Two shallow sources (label B and C) overlap on the top of source A with a depth extent of 0.2km. 25 The three sources are joined with each other, so these synthetic models can alternatively be considered as one whole source with an undulating upper surface. The 3D view, the plan view and the Φ zz data of the synthetic model are shown in Fig. 7.
The gravity gradient anomaly shown in Fig. 7(d) is simulated based on a 0.1 km regular grid, and the gradients of the synthetic model are derived according to the formulas that were proposed by Cooper (2006). We calculate the eigenvector v 1 , and then delineate the source centroids and outlines (Fig. 8) with the GGT data using Eq.(3) and Eq.(4). 30 Numerically, the values of tanφ are related to the source size--larger sources will show larger tanφ values. The tanφ value Nonlin. Processes Geophys. Discuss., doi:10.5194/npg-2016-75, 2017 Manuscript under review for journal Nonlin. Processes Geophys. clearly. This experiment also indicates that the value of tanφ is related with the horizontal area of the source, for the only model difference between the source B and the source C is their areas in horizontal plane.

Field data experiment
In the field data experiment, we apply the proposed method to a high-resolution airborne gravity gradient dataset over northeast Iowa and southeast Minnesota, U.S. (Drenth, et al., 2015). This GGT data is collected by FTG-003 Full Tensor 10 Gradiometer (FTG) system in 2013. The survey contains 94 east-west traverse-flying lines in 400 m apart. The field data were acquired with an 80m nominal flight height and subsequently terrain corrected (2400kg/m 3 ). The GGT data contain 481x481 observation points in the study region. The Φ z and Φ xy gravity gradient components and the Proterozoic geophysical interpretation map (Drenth et al. 2015) of the survey area are shown in Fig. 9 and Fig. 10.
In the center of the survey area, there exists an obvious low-density source (unit Ysp in Fig.10). Drenth et al. (2015) 15 interpreted the geophysical characteristics of it as a silicic pluton. They suggested that the large-amplitude gravity response of the Decorah complex has notable geophysical similarities to Keweenawan alkaline ring complexes. In this paper, we apply the GTA method on this data set with the various β parameter selections. The results are shown in Fig. 11.
In this experiment, the field data contain both positive and negative anomalies. As Fig. 11 shows, the GTA locates the centroid of positive and negative causative sources of varying scales simultaneously, while delineating the edges of them 20 clearly. With a smaller threshold value ( β <=80), Fig.11 (a)-(d) roughly outlines Proterozoic and part of Decorah complex sources. While using the larger threshold value ( β >160), the results indicate the centroid position and the boundaries of the primary sources in the survey area, as Fig. 11 (e)-(k) shows. The results show that there are 5671, 329, 82, 20, 4 peak values display for the parameter Beta set as 10, 40, 80,160 and 320, respectively. The number of local relative small peak value (e.g. 5671 for Beta>10) is much larger than the overall peak value (e.g. 4 for Beta>320). With the value of Beta increasing, the 25 centroid positions of larger sources start to appear in GTA map, while the small sources, or noise are filtered out gradually.
By increasing the value of β , the proposed eigenvector analysis method will become more insensitive to the presence of noise. Both of the lower density unit Ysp and the higher density gabbro of the Decorah complex can be well located in this procedure.
The parameter β is identified according to Eq.(5) to local the primary anomalous sources. We also display the procedure of 30 this strategy with this field data, as shown in Fig.12. Nonlin. Processes Geophys. Discuss., doi:10.5194/npg-2016-75, 2017 Manuscript under review for journal Nonlin. Processes Geophys. Published: 13 January 2017 c Author(s) 2017. CC-BY 3.0 License.
The index of the β parameter is incremented in powers of 2. As shown in Fig.12, the maximum curvature is reached at nearly 80, which is coincident with the experimental results in Fig.11(e) that shows the clear boundary range of the two primary sources. With a β value larger than 128, the number of detected sources is less than 10, and these points are mainly indicators of the position of the primary source centroids.
In addition, a lower bound for the β parameter can be identified via data resolution--if the approximate depth to sources of 5 interest is known, minimum resolvable feature width can be used to place a lower bound on the parameter. For example, in this experiment, the Proterozoic features of interest are at approximately 500 meters depth; therefore features smaller than approximately 750 meters to a kilometer are not resolvable. So the distance between the detected source centroids should not be less than 750 meters. According to this standard, the parameter β should be set equal or larger than 80 for outputting credible results. 10

Comparison with alternative boundary detection algorithms
We also compare the proposed algorithm with another GGT data boundary detection algorithm NGTE, and the usual Bouguer gravity data boundary detection algorithms on Φ z data, such as Tilt Angle (Miller & Singh, 1994), as Fig.13 shows.
As Fig. 13 (a)-(d) shows, although many of the methods commonly show good performance on the Bouguer gravity anomaly data processing, they are not very suitable for the GGT data. Fig.13 (f) shows the tanφ map that is defined in Eq.(3). For 15 visibility, the maximum value is limited to 12. This map shows that there is a huge negative-density source (Ysp) in the centroid and the positive-density sources comprising the Decorah complex in the southern direction. All of the centroids and the boundary details of these bodies are well-delineated by the proposed method.

Conclusions
We have introduced a new method to locate the centroids and boundaries of sources for full gravity gradient tensor data 20 across different scales. Based on the eigenvector decomposition of GGT data, we describe the theoretical basis and derive the procedure of the proposed method in detail. In the context of realistic source geometries, the interference of multiple sources has been specially discussed. We use two different synthetic models to demonstrate the efficacy and robustness of the method in the presence of multiple, overlapping sources with noise. These synthetic experiments also demonstrate the insensitivity of the method to the amplitude of the anomaly, working effectively even for weak anomalies. 25 We applied the technique to a field GGT dataset from northeast Lowa to test the proposed method, and two strategies to choose the scale parameter were derived and discussed in detail. The proposed method delineates source centroid positions and horizontal boundaries of sources in varying scales, for both positive and negative anomalies. Compared with conventional potential field data, our proposed method utilizes more dimensions and information that are contained in GGT data, and provides a clear map of the buried source bodies. This method is essentially based on an eigenvector analysis, so 30 individual small sources or random noise do not interfere significantly with the result. However, like other detection methods Nonlin. Processes Geophys. Discuss., doi:10.5194/npg-2016-75, 2017 Manuscript under review for journal Nonlin. Processes Geophys. Published: 13 January 2017 c Author(s) 2017. CC-BY 3.0 License.
it is also difficult to clearly detect individual small sources, which is an area of future research.
Ultimately, the proposed method is feasible and reliable to locate source positions with full gravity gradient tensor data, even in areas where the sources overlap. This can be used as an effective tool for geological interpretation, locating the positions of exploration wells, or in seeding 3D density inversion algorithms (e.g. Uieda and Barbosa,2012).