3-D Kernel Density Imaging based on the Euler Deconvolution of Tensor Gravity Data

: Traditional discrimination techniques for Euler deconvolution use only the color spectrums of structural indexes, without considering the spatial distribution characteristics and inherent relationships among the Euler solutions to separate adjacent causative sources. In the present study, a new approach was developed for discriminating uncorrelated Euler solutions from coherent solutions based on the focusing levels indicated by the probability density distributions generated using multivariate kernel density estimations (KDE). A novel multiple coverage technique was proposed by using a series of different sized moving windows over gridded gravity data, which formed tight clusters of Euler solutions for different sized causative sources. The results of the probability density distributions were obtained using a 3-D KDE method for the Euler solution subsets {x, y, z} of synthetical models, and real data from a survey conducted in British Columbia (Canada) which had successfully established more credible and meaningful geological models when compared with three other subsets.


Introduction
The crucial features of gravitational and magnetic fields is that they enhance and maintain the edges of geologic bodies. With the development of gravity gradient measurement methods and their corresponding interpretation techniques, many research findings regarding persistent edges have been published in recent years which have utilized such potential data sources as edge recognition [1]; eigenvector analyses [2]; and focusing inversion techniques [3]. Full tensor gravity gradient (FTG) technology is a multi-component gravimetry technique for measuring the different components of gravity gradient tensor [4]. In recent years, full tensor measurements have been widely carried out, and large amounts of high-precision tensor gravity data have been obtained. It has been found that the three-dimensional inversion of full tensor gravity data can achieve good results. In addition, full tensor gravity gradiometric surveys tend to have much higher bandwidths than traditional gravity surveys, and are able to deliver images with higher resolution [5].
However, FTG data impose higher requirements on their corresponding technical interpretations due to the fact that they often generate more complex anomalies for given sources [6].
Euler deconvolution processes based on Euler's homogeneity equation use gridded potential field data (for example, total fields and gradient fields) with series of square moving windows for the purpose of determining the source locations and depths of potential field anomalies for assumed structure indexes. Therefore, such methods are suitable for the analysis and interpretation of large-area potential field data [7,8]. The aforementioned methods have strong adaptability and flexibility, as well as simple calculation processes. The convergence rates of the tensor Euler deconvolution of full tensor gravity data have been observed to be better than those of traditional Euler deconvolution [9].
A linear system of equations is generally composed for each moving window, and solved using singular value decomposition or an equivalent technique in order to obtain such Euler solutions as the coordinates x, y, z, and the structure index N. The regional B and the  (assume x z    for a 3-D application) are constants which normally vanish, with the exception of the cases where N = 0 [10,11]. However, in practice, traditional Euler deconvolution has been found to easily produce spurious solutions when a system of equations is over-determined. Therefore, it is often difficult to determine the best solutions and discriminate the spurious solutions within an entire Euler solution dataset for geological bodies. In order to overcome the aforementioned difficulties, several methods have been proposed based on various posterior estimation strategies for the purpose of filtering out spurious solutions [12]. In previous studies, Gerovska et al. (2003) [13] filtered spurious Euler solutions using the properties of differential similarity transformation, which was based on the error estimations of singular point coordinates and structure indices. In another related study, Yao et al. (2004) [14] proposed three new criteria: Horizontal gradient filtering criterion; distance constraint criterion; and convergence constraint criterion, which could be utilized to filter spurious Euler solutions for the interpretations of magnetic and gravitational anomalies.
Kernel density estimations are used to reject solutions that are sparsely located and which are not parts of dense clusters. These estimations provide solid foundations for density-based clustering algorithm construction. Ugalde and Morris (2010) [8] applied clustering and kernel density distribution techniques to Euler solutions in order to obtain coherent filtering solutions from uncorrelated solutions. However, a drawback was encountered in that the aforementioned algorithm typically over-relied on the number of manually determined abnormal sources. The clustering analyses of Euler solutions are used to classify Euler solution sets according to similarity, and then classify the sets with the greatest similarity into one category, which occupy local regions of featured spaces. The aggregation centers of each local area play representative roles and indicate the spatial information of the corresponding geological bodies. Beiki (2013) [15] confirmed that the dimensionality of magnetic anomalies is the ratio between the smallest and intermediate eigenvalues acquired. Then, the removal of spurious Euler solutions can be completed using truncated singular value decomposition for the quasi-2D magnetic sources. On the other hand, the interpretation processes of the Euler solutions regarding the locations and depths of sources tend to utilize color levels or different sized circles. However, such methods have been found to have difficulty discriminating among different focusing degrees, such as the dense distributions (for example, the centroid of geological bodies) and sparse distributions (for example, the gaps among geological bodies) of the Euler solutions.
In this research investigation, Euler solutions were obtained using tensor Euler deconvolution of FTG data, which avoided the calculation of the potential field derivatives for coordinates x and y, and depth z, and did not predict the structure index N. Then, using a comparative analysis process, the probability density curves, images, and iso-surfaces were obtained using different dimensions of KDE for the combinations of the Euler solution datasets, and 1D and 2D KDE producing too many probability density distribution curves or images, making it difficult to analyze. Furthermore, illustration of the results of the probability density distributions, which were 5D data obtained by quaternion KDE for the entire Euler solution dataset, were also found to be difficult to achieve. Consequently, this study utilized 3D KDE to determine the focusing levels for finding the probability density peaks and removing the spurious Euler solutions. It was found that the proposed algorithm was able to successfully establish a greater amount of credible and meaningful geological bodies when compared with other discrimination techniques for Euler deconvolution.

Tensor Euler deconvolution
Full tensor gravity data is composed of the derivatives along the tangent vectors of gravity anomalies. The first derivative, second derivative, and the extension of a gravity field also satisfy the Euler equation. The tensor Euler deconvolution equation is composed of the traditional Euler Equation (1) [16] and two similar equations (2) and (3) [17,18].
  0 0 0 y y y y y y y y y where a B is the abnormal background of gravity; g  is the gradient gravity; g  represents the component of FTG,   , , , x y z    ; and a B is very small. Therefore, a B can be removed from Eqs. (1) to (3), and N will not be predicted [19]. Next, square moving windows over gridded gravity data were used to determine the locations of the causative sources. Each of the windows led to a solution involving the location of the causative sources {x, y, z} and structure index N using singular value decomposition or equivalent techniques.
As a result, Eqs. (1) to (3) were rewritten in matrix form as follows:

Multivariate KDE of the Euler solution datasets
In the present study, the linear system (4)  This study observed that in the available statistical data, KDE was considered to be a non-parametric approach for calculating the probability density distributions of data instances. For example, when compared with parametric approaches, non-parametric approaches make full use of all the data points and provide several advantages, such as asymptotic unbiasedness, square consistency, and uniform convergence. At the present time, KDE is a very popular tool for analyzing low-dimensional data instances [23,24].
Therefore, by letting provide a sample of independent and identically distributed random variables with a common density f , the kernel density estimator of f  can be defined as follows: where Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 October 2020 doi:10.20944/preprints202010.0059.v1 Then, for the purpose of convenience, where function K is the 1-D kernel function, 0 h  , which is the bandwidth. In the formula, the multivariate kernel function , d h K is the inner product of the multiple 1-D kernel functions [24], and d indicates the dimension of the dataset X .
Theoretically speaking, any function can be used as a kernel function. However, for the convenience and rationality of density function estimations, this study used the Gaussian Kernel as follows: where Therefore, in accordance with above definition, the extended form of the d-dimensional multivariate kernel density estimator f  of f can be written as follows: Therefore, by taking the partial derivative of Eq. (10) with respect to f , the yield kernel density derivative estimator will be as follows: where  is the column vector of the d partial first derivatives, and This study used the multivariate KDE to analyze the various combinations, such as   x , , , x y z , of the Euler solution datasets for the purpose of avoiding the involvement of the above-mentioned discrimination techniques [25,26].  among all the cubes, had displayed relative stability. Therefore, by comprehensively considering the aforementioned curves, this study concluded that the selection of the sizes of the moving windows was difficult for the small causative bodies. However, the multiple coverage technique for Euler deconvolution had effectively formed tight clusters.

Verification of the validity of the KDE algorithm
On the premise of not rejecting any Euler deconvolution solutions, three simple models were designed in this study for the purpose of verifying the validity of the multivariate KDE algorithm. Analytical solutions were used to calculate the three simple models. The dimensions of the cubes were 1,000 m × 1,000 m × 1,000 m, with the centroid located at (-1,000, -2,000, and 1,500). The radii, lengths, and centroid of the horizontal cylinders were 1,000 m, 4,000 m, and (0, 0, 2,500), respectively. The sphere radius was 1,000 m, with its centroid was located at (2500, 2500, and 3,500). There were three bodies with the same density contrast (0.36 Then, the tensor Euler deconvolution algorithm and inverted FTG data contaminated by 3% noise were used for the x, y, z, and N with singular value decomposition. The probability values of a 1-D subset (for example, x, y, z, and N) of the Euler solutions were obtained using a 1-D KDE method, as shown in Fig. 3. The kernel function in the KDE method was a Gaussian function, and its density function was a probability density estimation function. Further details can be referenced from Bowman and Azzalini (1997). Fig. 3a shows the results of the coordinate x of probability density peaks, which were determined to be equal to -1,000 m, 0 m, and 2,600 m, for the cube, cylinder, and sphere, respectively. Those three values were found to be in accordance with the designed models. Fig. 3b shows that there was a very good agreement between the coordinate y of the probability density peaks and this study's designed models of the cube and sphere. It was found that when compared with other methods, there were two equal peak values of the probability density at the top and bottom ends of the cylinder (y = ± 4,000m). Fig. 3c details that there were two peaks in the probability density curve for each of the designed models. The first peak had corresponded to the shallow buried Euler solutions, which were considered to be spurious Euler solutions due to the instability of the Euler equation. The second peaks of three curves corresponded exactly to the theoretical values of the designed models. The same situation can also be observed in Fig. 3d.
The best locations and depths of the causative sources obtained from the optimum estimates had successfully verified the validity of the 1-D KDE algorithm. However, only single probability density peaks were found to correspond to the theoretical values. Therefore, in order to overcome these difficulties, this study used a 2-D KDE method to analyze the six As shown in Fig. 4, the probability density of the cube and the sphere had only one peak, and there were two peaks in the probability density image corresponding to the ends of the cylinder. In regard to the other combinations, due to the limited length of this report, only the relevant results were attached. The probability density of the Euler solutions along two poles of the cylinder were characterized by a scattering distribution. In accordance with the equivalence principle, that type of phenomenon could also be produced by two adjacent cubes, indicating that ambiguities existed.
The locations and structural index were easily interpreted in this study using 2-D KDE. However, the 2-D KDE for the different combinations of Euler solution datasets had the drawback of producing too many results for the probability density distributions, making the data difficult to analyze. Furthermore, the visualization of the results using the 5-D data obtained by quaternion KDE for the entirety of the Euler solution datasets was found to be very challenging. Consequently, this study introduced 3-D KDE for the Euler solution subsets, such as {x, y, z}, {x, y, N}, {x, z, N}, and {y, z, N}. Then, there were a total of twelve results for the probability density distributions for the previously mentioned three designed models. In addition, for ease of presentation, this study only illustrated the results of the probability density distributions of the cylinder. For all of the experimental tests of the 3-D KDE in the remainder of this study, the significance level  was equal to 0.05 for three levels of the transparency color of the 3-D contours of the probability density data. The colors indicated the focusing levels of the clustered dataset of the Euler solutions. Further details are available from the research presented by [25,27]. The bandwidth parameters of the 3-D KDE were observed to be more stable than those obtained with the 2-D KDEs, as shown in Table 2. Further details of the bandwidth parameters are available in the research presented by [25]. In this study's experimental tests of the cylinder, the iso-surfaces for the probability density distributions obtained using 3-D KDE for all the subsets, including the y component, were found to have two cluster sources along the Y-axis direction, as detailed in Figs. 5a, 5b, and 5d. The above-mentioned three results provided mutual corroboration. Moreover, those results were found to be in accordance with Figs. 3 and 4. Therefore, by comprehensively considering the aforementioned figures, it was concluded that the final solution for the cylinder tests was a single causative body with a structural index between 1.3 and 2.4.

Influences of the adjacent causative bodies on the 3-D KDE for the Euler solution subset {x, y, z}
It should be pointed out that the results of the three subsets {x, y, N}, {x, z, N}, and {y, z, N} were easily confused with the adjacent causative bodies or similar causative bodies, as shown in Fig. 5. The interpretations of the Euler solution datasets using 3-D KDE relied mainly on the subset {x, y, z}, while making the assistance of the other subsets subsidiary. Therefore, this study chose to illustrate the results for subset {x, y, z} in the remainder of this research report. Fig. 6a details the forward modeling of the FTG for the designed models. There were cubic bodies with densities of 0.36 × 103 kg/m3 and 1.27 × 103 kg/m3 for the left and right, respectively, in which the centroids were (-4,000, 4,000, 2,500) and (-4,000, 4,000, 2,500) for the left and right, respectively. The size of each body was 2,000 m × 2,000 m × 2,000 m. The rectangular survey grid was divided into 40,000 = 200 × 200 observation points located on the Earth's surface. The observation height was set as 25 m, and the sampling interval was 100 m in the x-and y-directions. Then, the two causative bodies were allowed to approach each other with L = 1,500 (moving step) in the x-and y-directions at the same time. In the current investigation, 3-D KDE was used to analyze the Euler solutions subset {x, y, z} for the FTG data contaminated by 3% noise. Figs. 6a to 6c show that the two geological bodies with different gaps had been directly discriminated using 3-D KDE, and the depths and probability density distributions were in accordance with the designed model.

Case study completed in British Columbia (Canada)
The province of British Columbia (Canada) has been thoroughly explored for copper and gold porphyry deposits. However, the region is covered by a thick layer of sand and gravel left behind by glaciers. In the aforementioned case study, airborne gravity surveys were effectively used to identify potential kimberlite targets and assist geophysicists understand the geology of the British Columbia area. The survey blocks consisted of two parallelograms with two pairs of parallel sides. The first parallelogram measured approximately 386 km in length and 120 km in width. The second was 120 km in length and 60 km width and appended to the north, as illustrated in Fig. 7a [28]. The measured heights ranged from 240 m to 2,850 m above the mean sea level, making it necessary to continue the irregular gravity anomaly upward to a horizontal plane. Fig. 8a details the Bouguer gravity anomaly. The gravity anomaly data were continued upward to 3,000 m and 6,000 m above the mean sea level in order to obtain a new gravity anomaly and the background field, respectively. Then, the background field from the new gravity anomaly was selected in order to obtain the local gravity anomaly. The FTG data were converted from the local gravity anomaly using a Fast Fourier Transform Method [29]. However, due to the space limitations of this report, only the Gzz component of the FTG was illustrated, as shown in Fig. 8b. In order to form tightly clustered datasets of the Euler solutions, the entire survey gravity map was covered with a series of different sized moving windows, which ranged in size from seven to 25. There were a total of 1,088,276 solutions generated using tensor Euler deconvolution for the FTG data before the solutions were filtered. The Euler solutions were illustrated with 2-D and 3-D color spectra, respectively, as shown in Fig. 9. It was only possible to roughly determine the number of causative sources, as shown in Fig. 9a. In particular, difficulties were encountered in discriminating the adjacent causative bodies (Fig. 9a), and also determining the depths of the causative bodies (Fig. 9b).  Fig. 9, the proposed algorithm had successfully established more credible and meaningful geological models, as shown in Fig. 11. This study was able to easily determine the locations and the depths of the causative sources, as well as the weights among all the causative sources from the probability density distributions.

Conclusions
This study developed a new approach to discriminating uncorrelated Euler solutions from coherent solutions. The proposed method was based on the focusing levels indicated by probability densities generated using multivariate KDE. A novel multiple coverage technique applied using a series of different sized moving windows which forced the causative sources to form tighter clusters. A number of conclusions could be drawn by calculating one synthetic model. This study's conclusions were as follows: (1) The clusters generated by Euler deconvolution had easily deviated from each other when influenced by the equivalence principle.
(2) It was observed that difficulties were encountered when choosing the most effective sizes for the moving windows for the smaller causative bodies.
(3) The multiple coverage technique for the Euler deconvolution had effectively formed tighter clusters.
In the present study, the comparative analysis probability density curves and images were obtained using various dimensions of KDE for the different combinations of Euler solution datasets. The conclusions were as follows: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 October 2020 doi:10.20944/preprints202010.0059.v1 (1) Multiple probability density peaks existed for one geological body in the results obtained using 1-D and 2-D KDE. Therefore, it was difficult to determine which of the probability density peaks corresponded to the true values.
(2) For the 2-D KDE, it was found that there were too many results for the probability density distributions, which made it difficult to analyze the entirety of Euler solution datasets.
(3) The visualization of the results with the 5-D data obtained by 4-D KDE for the entirety of the Euler solution datasets was found to be difficult to achieve.
(4) In regard to the 3-D KDE, the results of the three subsets {x, y, N}, {x, z, N}, and {y, z, N} were easily confused with the adjacent causative bodies or similar causative bodies. Therefore, the 3-D KDE for the Euler solutions subset {x, y, z} of the synthetical model, and the real data from a survey conducted in British Columbia (Canada), were used to successfully establish more credible and meaningful geological models, when compared with the other three subsets.
The approach presented in this study was considered to have the potential to also be applied to gravity, magnetic, and full tensor gradient magnetic data.