Open AccessArticle

FFT-Based Probability Density Imaging of Euler Solutions

Shujin Cao

^1,2,3

Peng Chen

¹,

Guangyin Lu

^2,*

Zhiyuan Ma

¹,

Bo Yang

¹ and

Xinyue Chen

School of Earth Sciences and Spatial Information Engineering, Hunan University of Science and Technology, Xiangtan 411201, China

School of Geosciences and Info-Physics, Central South University, Changsha 410083, China

Institute of Geophysics & Geomatics, China University of Geosciences, Wuhan 430074, China

Author to whom correspondence should be addressed.

Entropy 2024, 26(6), 517; https://doi.org/10.3390/e26060517

Submission received: 5 February 2024 / Revised: 13 June 2024 / Accepted: 13 June 2024 / Published: 15 June 2024

(This article belongs to the Topic Interfacing Statistics, Machine Learning and Data Science from a Probabilistic Modelling Viewpoint)

Download

Browse Figures

Versions Notes

Abstract

When using traditional Euler deconvolution optimization strategies, it is difficult to distinguish between anomalies and their corresponding Euler tails (those solutions are often distributed outside the anomaly source, forming “tail”-shaped spurious solutions, i.e., misplaced Euler solutions, which must be removed or marked) with only the structural index. The nonparametric estimation method based on the normalized B-spline probability density (BSS) is used to separate the Euler solution clusters and mark different anomaly sources according to the similarity and density characteristics of the Euler solutions. For display purposes, the BSS needs to map the samples onto the estimation grid at the points where density will be estimated in order to obtain the probability density distribution. However, if the size of the samples or the estimation grid is too large, this process can lead to high levels of memory consumption and excessive computation times. To address this issue, a fast linear binning approximation algorithm is introduced in the BSS to speed up the computation process and save time. Subsequently, the sample data are quickly projected onto the estimation grid to facilitate the discrete convolution between the grid and the density function using a fast Fourier transform. A method involving multivariate B-spline probability density estimation based on the FFT (BSSFFT), in conjunction with fast linear binning appropriation, is proposed in this paper. The results of two random normal distributions show the correctness of the BSS and BSSFFT algorithms, which is verified via a comparison with the true probability density function (pdf) and Gaussian kernel smoothing estimation algorithms. Then, the Euler solutions of the two synthetic models are analyzed using the BSS and BSSFFT algorithms. The results are consistent with their theoretical values, which verify their correctness regarding Euler solutions. Finally, the BSSFFT is applied to Bishop 5X data, and the numerical results show that the comprehensive analysis of the 3D probability density distributions using the BSSFFT algorithm, derived from the Euler solution subset of

\{x_{0}, y_{0}, z_{0}\}

, can effectively separate and locate adjacent anomaly sources, demonstrating strong adaptability.

Keywords:

Bishop model; Euler deconvolution; B-spline density estimation; probability density imaging; fast linear binning approximation; fast Fourier transform

1. Introduction

Euler deconvolution is suitable for the semi-automatic estimation of source locations in large-scale potential field data [1,2,3,4,5]. It originates from the overdetermined system of Euler with a large condition number, which causes significant perturbations and often shows divergent trends in Euler solutions [6,7,8,9,10]. Therefore, many researchers use the standard deviation [6] or the truncation error [11] of the overdetermined system of Euler as criteria to filter out spurious solutions. As indicated by FitzGerald et al. [7], correct Euler solutions exhibit the clustering characteristic: “Euler deconvolution yields many similar or duplicate solutions, which may be tightly clustered in the vicinity of real sources.” This is frequently exploited by clustering algorithms to identify anomalous sources. Some additional constraint equations or conditions, such as the horizontal components of the gravity gradient tensor or the relationship between “worming” and Euler equations [12,13], have been proposed for the construction of mixed Euler deconvolution methods [14], seeking to obtain optimal solutions by determining the cluster relationships formed by the Euler solutions and eliminating any spurious solutions that do not belong to a cluster [6,7,10,15,16]. An adaptive fuzzy clustering algorithm has been introduced to address challenges such as the need for predetermined cluster numbers, local optimization in fuzzy clustering analysis, and classification uncertainty [17,18]. However, traditional cluster analysis methods struggle to find the best Euler solution cluster [6,10,14,15].

As a simple application of the clustering characteristic of Euler solutions, the histogram method divides the entire range of structural index values into a series of adjacent and non-overlapping equally spaced bins. It then calculates the number of values in each interval to determine the optimal structural index. This is practical for simple geological structures. However, for complex geological structures, or when there are multiple anomaly sources, it is unrealistic to simply use the histogram method to determine one or more of the “best” structural indices. The density histogram is employed to visualize each formation and to check the coherency of the inversion process [19]. However, it is difficult to efficiently divide a multi-dimensional space into non-coherent histogram bins while keeping the error rates low [20].

A probability density distribution is a continuous version of a histogram with densities, and it is widely used in geophysical inversions. FitzGerald et al. [7] pointed out that the Euler deconvolution yields many similar or duplicate solutions, which may be tightly clustered in the vicinity of real sources, i.e., the 14th strategy in the paper [7]. Based on this strategy, Cao et al. [21] proposed a non-parametric estimation method based on the normalized B-spline probability density (BSS), aiming to separate the Euler solution clusters to mark different anomaly sources according to the similarity and density characteristics of the Euler solutions [22,23]. B-spline curves are linear combinations of B-spline basis functions connected by node vectors [24]. The c + 1 order B-spline basis functions can all be obtained recursively from the lower-order basis functions [25]. Faenza et al. [26] developed a novel non-parametric multivariate model to depict the spatial and temporal distributions of large earthquakes. The model enables the intuitive examination of diverse hypotheses, including all forms of time dependency (e.g., seismic gaps, clustering, and Poisson assumptions) [27,28]. Mautz et al. [29] proposed the so-called B-spline wavelets to represent the spatial time series of the GPS-derived global ionosphere maps of the vertical total electron content of the Earth’s surface to the mean altitudes of GPS satellites over Japan. Herrmann and Hennenfent [30] proposed a recovery method based on non-parametric transformation, utilizing recently developed curvelet frames to compress seismic data [31].

To display and calculate the probability density value of the sample, it is generally necessary to construct an estimation grid, which consists of a set of equidistant or non-equidistant spatial points and generally needs to cover the entire input sample. However, when the data sample is large or the estimation grid is extensive, the BSS algorithm, using brute force computation [32,33], i.e., traversing each grid node over the data sample, becomes computationally intensive and consumes excessive amounts of memory. Wang et al. [34] proposed a flexible and high-precision algorithm to perform transformations of the potential field and gradient components in 3D and 2D, as well as derivative calculations in the spatial domain using cubic B-splines, applicable at any point within the computational domain. To reduce the computational intensity of kernel density estimation, Wand [32] proposed the binned kernel density estimation, which first uses the fast linear binning technique to map the sample data onto the predefined estimation grid and then obtains the grid counts of the sample data concerning the nodes of the estimation grid, to compute the approximate value of the kernel estimation on the estimation grid over the sample data using the fast Fourier transform (FFT) [35,36]. Based on Wand [32], Cao et al. [37] focused on the computational details of the sample projection in the estimation grid. Based on fast linear binning approximation and Fourier-based fast convolution, the multivariate kernel density derivative estimation (KDDE) was proposed to compute the probability values of Euler solutions derived from tensor gravity data using tensor Euler deconvolution. The algorithm is an extension of kernel density estimation [38], which is a mathematical process used to estimate the probability density function of a random variable. The principle and procedure of the fast linear binning approximation, Fourier-based convolution, and the KDDE algorithm for computing the probability values of Euler solutions were detailed. Then, the effects of the bandwidths and the estimation grid’s size on the computational efficiency of fast linear binning approximation and the KDDE algorithm were analyzed. The results obtained from synthetic models and field data showed that the multivariate KDDE, in conjunction with fast linear binning approximation and the FFT, demonstrated high computational efficiency to overcome this difficulty. Moreover, the resulting probability density isosurface successfully locates meaningful geological targets.

This paper proposes a B-spline probability density estimation method based on the FFT (BSSFFT) in conjunction with a fast linear binning approximation, aiming to rapidly map Euler solutions to the estimation grid. Using a fast Fourier transform, the discrete convolution of the estimation grid and the density function is achieved. Compared with the BSS, BSSFFT can effectively avoid the issue of brute force calculation. Simultaneously, it provides a more focused probability density estimation result for identifying various anomaly sources. Numerical tests and synthetic models demonstrate the reliability of the BSS and BSSFFT. Finally, the algorithm was applied to magnetic field data from the Bishop model, which was generated from topographic data from a portion of the volcanic tablelands area north of Bishop, California. The probability density isosurfaces of the synthetic models and the field data obtained by BSSFFT distinguish the clusters formed by the Euler solutions, effectively separating and locating adjacent anomaly sources.

2. Theory

2.1. Principle of Gravity Euler Deconvolution

Following the suggestions of Reid [39] and based on the 14th strategy proposed by FitzGerald et al. [7], the Euler deconvolution can be expressed as follows:

x_{0} \frac{\partial g_{z}}{\partial x} + y_{0} \frac{\partial g_{z}}{\partial y} + z_{0} \frac{\partial g_{z}}{\partial z} + N (B_{z} - g_{z}) = x \frac{\partial g_{z}}{\partial x} + y \frac{\partial g_{z}}{\partial y} + z \frac{\partial g_{z}}{\partial z}

(1)

where

α \in \{x, y, z\}

;

B_{z}

is the anomalous background field;

g_{z}

is the gravity field; and

\partial g_{z} / \partial β

is the partial derivative of

g_{z}

with respect to

β

(x, y, z)

are the coordinates of the observation point; and

(x_{0}, y_{0}, z_{0})

are the spatial positions of the sought anomaly source. When the term regarding

\frac{\partial g_{z}}{\partial y}

is removed, Equation (1) can be applied to the 2D profile data used to obtain

{x_{0}, z_{0}, N}

. The structural index N depends on the nature (type) of the geologic source, which makes it challenging to characterize multiple anomalous sources using a single structural index [40,41]. Due to the difficulty in determining the N of the geological target of interest [7,42], the traditional Euler deconvolution uses a series of predefined structural indices and yields many spurious solutions [10]. Many researchers have attempted to determine the best cluster via manual evaluation using a series of tentative N values, which significantly increases the interpreter’s workload, particularly for complex geological structures with multiple anomalous sources [6,10,43,44,45]. To overcome these problems, N will not be predicted in this paper [6,7,42].

Utilizing a square sliding window of size

w_{x}

w_{y}

=w to gradually traverse the entire survey area, and using Equation (1), each sliding window can form an equation system, as follows:

[\begin{matrix} {(\partial g_{z} / \partial x)}^{1} & {(\partial g_{z} / \partial y)}^{1} & {(\partial g_{z} / \partial z)}^{1} & {(g_{z} - B_{z})}^{1} \\ \dots & \dots & \dots & \dots \\ {(\partial g_{z} / \partial x)}^{j} & {(\partial g_{z} / \partial y)}^{j} & {(\partial g_{z} / \partial z)}^{j} & {(g_{z} - B_{z})}^{j} \\ \dots & \dots & \dots & \dots \\ {(\partial g_{z} / \partial x)}^{n_{w}} & {(\partial g_{z} / \partial y)}^{n_{w}} & {(\partial g_{z} / \partial z)}^{n_{w}} & {(g_{z} - B_{z})}^{n_{w}} \end{matrix}] [\begin{matrix} x_{0}^{i} \\ y_{0}^{i} \\ z_{0}^{i} \\ N^{i} \end{matrix}] = [\begin{matrix} b_{z}^{1} \\ \dots \\ b_{z}^{j} \\ \dots \\ b_{z}^{n_{w}} \end{matrix}]

(2)

where

b_{z}^{j} = {(x \partial g_{z} / \partial x)}^{j} + {(y \partial g_{z} / \partial y)}^{j} + {(z \partial g_{z} / \partial z)}^{j}

n_{w} = w_{x} \times w_{y}

; the Euler deconvolution solution set is denoted as

m = {m^{i}}

; and the i-th Euler deconvolution solution is denoted as

m^{i} = {[x_{0}^{i}, y_{0}^{i}, z_{0}^{i}, N_{0}^{i}]}^{T}

, representing the spatial location and structural indices of the sought anomaly source. The superscript T denotes the transpose operator.

2.2. Principle of BSS for Euler Solutions

Traditional Euler deconvolution optimization strategies are not able to differentiate optimal solutions from pseudo-traces (misplaced solutions, which must be removed or marked) among Euler solutions. Therefore, the BSS algorithm, which is based on normalized B-splines, is adopted to separate each cluster based on the similarity and clustering degree of the Euler solutions in order to identify various sources of anomalies [7]. To introduce the BSS algorithm [21], let

X^{1}, \dots, X^{i}, \dots, X^{n}

be an independent and identically distributed random sample of a d-dimensional probability density function (pdf) f. The sample ith can be represented as

X^{i} = {[X_{1}^{i}, \dots, X_{j}^{i}, \dots, X_{d}^{i}]}^{T}

. Specifically,

X^{i}

can denote different combinations of Euler solutions,

x_{0}^{i}

y_{0}^{i}

z_{0}^{i}

N^{i}

, such as

x_{0}^{i}

{[x_{0}^{i}, y_{0}^{i}]}^{T}

{[x_{0}^{i}, y_{0}^{i}, z_{0}^{i}]}^{T}

{[x_{0}^{i}, y_{0}^{i}, z_{0}^{i}, N^{i}]}^{T}

1 \leq d \leq 4

Projecting samples

X^{i}

onto an estimation grid

χ = {[χ_{1}, \dots, χ_{j}, \dots, χ_{d}]}^{T}

involves a grid designed as a unified structure to accurately cover the input samples, with a size of

λ = {[λ_{1}, \dots, λ_{j}, \dots, λ_{d}]}^{T}

, a bandwidth of

h = {[h_{1}, \dots, h_{j}, \dots, h_{d}]}^{T}

, and upper and lower bounds of

\{a_{j}\}

\{b_{j}\}

. The relationship between these variables is as follows:

h_{j} = (b_{j} - a_{j}) / (λ_{j} - 1)

(3)

On an estimate point

\bar{x} = [{\bar{x}}_{1}, \dots, {\bar{x}}_{j}, \dots, {\bar{x}}_{d}]

of the estimation grid

χ

, for a sample

X^{i}

, the estimation of the probability density of the d-dimensional B-spline can be expressed using c uniform B-spline functions

B_{l_{j}}^{c} (\bar{x})

, denoted as

\hat{f}

, given by the following Gehringer and Redner [46]:

\hat{f} (\bar{x}) = \sum_{l_{1} = 1}^{λ_{1}} \dots \sum_{l_{d} = 1}^{λ_{d}} \frac{a_{l_{1} \dots l_{j} \dots l_{d}}}{h_{1} \dots h_{j} \dots h_{d}} B_{l_{1}}^{c} ({\bar{x}}_{1}) \dots B_{l_{j}}^{c} ({\bar{x}}_{j}) \dots B_{l_{d}}^{c} ({\bar{x}}_{d})

(4)

Where the subscript

l_{j}

specifies a position along the jth dimension of the estimation grid

χ

, and the coefficient

a_{l} = a_{l_{1} \dots l_{j} \dots l_{d}}

in the equation can be expressed by sample

X^{i}

as follows:

a_{l_{1} \dots l_{j} \dots l_{d}} = \frac{1}{n} \sum_{i = 1}^{n} B_{l_{1}}^{c} (X_{1}^{i}) \dots B_{l_{j}}^{c} (X_{j}^{i}) \dots B_{l_{d}}^{c} (X_{d}^{i})

(5)

where, for sample

X^{i}

, based on the normalized cth-order uniform B-spline

B^{c}

and the jth-dimensional bandwidth

h_{j}

, the B-spline basis functions

B_{l_{j}}^{c} ({\bar{x}}_{j})

can be redefined as follows:

B_{l_{j}}^{c} ({\bar{x}}_{j}) \equiv B_{l_{j}}^{c} (X_{j}^{i}) = B^{c} (\frac{X_{j}^{i}}{h_{j}} - ⌊\frac{X_{j}^{i}}{h_{j}}⌋)

(6)

where

B^{c}

is a normalized

c^{t h}

-order uniform B-spline.

Theorem 1.

The uniform B-splines,

\{B_{k}^{c} (\bar{x})\}

, form a partition of unity [47], i.e.,

\sum_{k = - \infty}^{\infty} B_{k}^{c} (\bar{x}) = \sum_{k = j - c + 1}^{j} B_{k}^{c} (\bar{x}) = 1

(7)

for all

\bar{x} \in [{\bar{x}}_{j}, {\bar{x}}_{j + 1}]

where

c =

1, 2, and 3,

B^{c} (\bar{x})

in Equation (6) can be expressed as follows:

B^{1} (\bar{x}) = \{\begin{matrix} \bar{x}, & 0 \leq \bar{x} < 1 \\ (2 - \bar{x}), & 1 \leq \bar{x} \leq 2 \end{matrix}

(8)

B^{2} (\bar{x}) = \{\begin{matrix} {\bar{x}}^{2} / 2, & 0 \leq \bar{x} < 1 \\ (- 2 {\bar{x}}^{2} + 6 \bar{x} - 3) / 2, & 1 \leq \bar{x} < 2 \\ {(3 - \bar{x})}^{2} / 2, & 2 \leq \bar{x} \leq 3 \end{matrix}

(9)

B^{3} (\bar{x}) = \{\begin{matrix} {\bar{x}}^{3} / 6, & 0 \leq \bar{x} < 1 \\ (- 3 {\bar{x}}^{3} + 12 {\bar{x}}^{2} - 12 \bar{x} + 4) / 6, & 1 \leq \bar{x} < 2 \\ (3 {\bar{x}}^{3} - 24 {\bar{x}}^{2} + 60 \bar{x} - 44) / 6, & 2 \leq \bar{x} < 3 \\ (- {\bar{x}}^{3} + 12 {\bar{x}}^{2} - 48 \bar{x} + 64) / 6, & 3 \leq \bar{x} \leq 4 \end{matrix}

(10)

2.3. Principle of BSSFFT

Through Equation (4), the probability density values of each Euler solution can be calculated on the estimation grid

χ

. To overcome the issue of brute force computations [37,48], an improved BSS method is proposed based on fast linear binned approximation [32]. Equation (4) can be reconstructed as follows:

\hat{f} (χ_{j}) = n^{- 1} \sum_{l_{1} = 1}^{λ_{1}} \dots \sum_{l_{d} = 1}^{λ_{d}} B_{K}^{c} (χ_{j}^{l_{j}} - χ_{l}^{l_{j}}) C_{l}

(11)

where

C_{l}

is the grid count obtained via approximate fast binning at point

χ_{l}

of the estimation grid,

1 \leq l_{j} \leq λ_{j}

To compute

C_{l}

in Equation (11), we define it in the j-dimensional direction with edge

[χ_{j}^{l_{j}}, χ_{j}^{l_{j + 1}}]

and denote it as

χ_{j}^{l_{j}, l_{j} + 1}

. Assume that

X^{i}

is precisely located in a hyper-rectangle (also called a bin) V formed by the edges

χ_{1}^{l_{1}, l_{1} + 1}, \dots, χ_{j}^{l_{j}, l_{j} + 1}, \dots χ_{d}^{l_{d}, l_{d} + 1}

, with

2^{d}

nodes and

2^{d}

faces [49]. Then, the projection of

X_{j}^{i}

χ_{j}^{l_{j}, l_{j} + 1}

can be used to represent weights

u_{j}^{0}

and

u_{j}^{1} = 1 - u_{j}^{0}

at both ends of edge

χ_{j}^{l_{j}, l_{j} + 1}

, as follows:

u_{j}^{0} = \frac{X_{j}^{i}}{h_{j}} - ⌊ \frac{X_{j}^{i}}{h_{j}} ⌋ \equiv \frac{X_{j}^{i} - χ_{j}^{l_{j}}}{h_{j}}

(12)

where

⌊ ⌋

represents the floor operation.

In other words,

u_{j}^{0}

and

u_{j}^{1}

are equivalent to the natural coordinate values of

X^{i}

induced by the two endpoints of edge

χ_{j}^{l_{j}, l_{j} + 1}

. Taking a vertex v of V as an example, according to Rao [50] and Chacón and Duong [49], the corresponding grid count

{\bar{C}}_{l}

of sample

X^{i}

is equivalent to the multiplication of the natural coordinate values (

u_{j} \in \{u_{j}^{0}, u_{j}^{1}\}

) induced by

X^{i}

with respect to the adjacent edges of the vertex of V, which means that

{\bar{C}}_{l}

is equal to the natural coordinate values of v with respect to

X^{i}

{\bar{C}}_{l}

can be expressed as follows:

{\bar{C}}_{l} = \prod_{j = 1}^{d} u_{j}

(13)

As shown in Figure 1, suppose that a trivariate point

X^{i}

is located exactly in a bin surrounded by 8 points

χ_{j_{1} j_{2} j_{3}}

, ⋯,

χ_{(j_{1} + 1) (j_{2} + 1) (j_{3} + 1)}

. To obtain the count

{\bar{C}}_{l}

of the point with respect to its 8 neighboring vertices, we project sample

X^{i}

onto the 8 edges and 6 faces of the bin, and the points on the edges

v_{1}, v_{2}, \dots, v_{8}

and the points on the faces

s_{1}, s_{2}, \dots, s_{6}

are obtained. Then, three lines,

s_{1} s_{3}, s_{2} s_{4}

, and

s_{5} s_{6}

, are formed to pass through sample

X^{i}

. Using planes parallel to the coordinate axes and passing through these lines, the bin is divided into sub-bins

A_{1}, \dots, A_{8}

from top to bottom in a clockwise direction. Taking the point

χ_{j_{1} (j_{2} + 1) j_{3}}

as an example, according to Equation (12), the lines

s_{1} s_{3}, s_{2} s_{4}

, and

s_{5} s_{6}

are each divided into two segments and labeled

u_{1}^{0}, u_{1}^{1}, u_{2}^{0}, u_{2}^{1}, u_{3}^{0}

, and

u_{3}^{1}

, respectively. Thus, according to Equation (13), the grid count of

χ_{j_{1} (j_{2} + 1) j_{3}}

with respect to

X^{i}

is as follows:

{\bar{C}}_{l} = \prod_{j = 1}^{d} u_{j} = u_{1}^{1} u_{2}^{1} u_{3}^{1} = A_{6}

(14)

Thus, when a sample

X^{i}

falls within a bin, for a node v of the bin, its corresponding grid count with respect to the sample is determined by the length (in 1D), area (in 2D), or volume (in 3D) that the diagonal node from this node v encloses with the sample, as shown in Figure 1. For the calculation of the count

{\bar{C}}_{l}

at an estimated point for univariate and bivariate

X^{i}

, please see Cao et al. [37]. For each estimation point

χ_{l}

, its grid count

C_{l}

is the sum of the grid counts

{\bar{C}}_{l}

of all samples pertaining to it.

According to the algorithm of Wand [32], Equation (11) can be rewritten in the form of convolution, as follows:

{\hat{f}}_{l} = \sum_{k_{1} = - λ_{1} + 1}^{λ_{1} - 1} \dots \sum_{k_{j} = - λ_{j} + 1}^{λ_{j} - 1} \dots \sum_{k_{d} = - λ_{d} + 1}^{λ_{d} - 1} C_{l - k} r_{k}

(15)

where

r_{k} = \frac{1}{n} B_{k}^{c} (h_{1} k_{1}, \dots, h_{j} k_{j}, \dots, h_{d} k_{d})

(16)

where

k_{j}

is the value of f along the j dimension of the estimation grid, and the convolution of C and r in Equation (15) can be computed using the discrete convolution theorem. When matrices C and r differ in size, techniques such as zero-padding and wrap-around ordering (see Chapter 13 of Teukolsky et al. [51]) are used to ensure that C and r match in size, leveraging the special structured matrix properties discussed in [51,52]. For ease of understanding, only the 2D scenario, as expressed in Equations (17) and (18), is provided here.

C = [\begin{matrix} C_{1, 1} & C_{1, 2} & \dots & C_{1, λ_{2}} \\ ⋮ & ⋮ & ⋱ & ⋮ & 0 \\ C_{λ_{1}, 1} & C_{λ_{1}, 2} & \dots & C_{λ_{1}, λ_{2}} & \dots \\ 0 & ⋮ & 0 \end{matrix}]

(17)

r = [\begin{matrix} r_{0, 0} & r_{0, 1} & \dots & r_{0, λ_{2}} & r_{0, λ_{2}} & \dots & r_{0, 1} \\ r_{1, 0} & r_{1, 1} & \dots & r_{1, λ_{2}} & r_{1, λ_{2}} & \dots & r_{1, 1} \\ ⋮ & ⋮ & ⋱ & ⋮ & 0 & ⋮ & ⋱ & ⋮ \\ r_{λ_{1}, 0} & r_{λ_{1}, 1} & \dots & r_{λ_{1}, λ_{2}} & \dots & r_{λ_{1}, λ_{2}} & \dots & r_{λ_{1}, 1} \\ 0 & ⋮ & 0 & ⋮ & 0 \\ r_{λ_{1}, 0} & r_{λ_{1}, 1} & \dots & r_{λ_{1}, λ_{2}} & \dots & r_{λ_{1}, λ_{2}} & \dots & r_{λ_{1}, 1} \\ ⋮ & ⋮ & ⋱ & ⋮ & 0 & ⋮ & ⋱ & ⋮ \\ r_{1, 0} & r_{1, 1} & \dots & r_{1, λ_{2}} & r_{1, λ_{2}} & \dots & r_{1, 1} \end{matrix}]

(18)

Then, C and r will have the same dimensions as

λ_{1} \times, \dots, \times λ_{j} \times, \dots, \times λ_{d}

. According to the discrete convolution theorem [32], Equation (19) can be re-expressed as follows:

S = F^{- 1} (F (C) F (r))

(19)

where F and

F^{- 1}

are denoted as the Fourier transform and the Fourier inverse transform.

The sought density estimate

\hat{f}

can be obtained by dividing the submatrix of S in Equation (19) by the product of

λ_{j}, j \in \{1, \dots, d\}

\hat{f} = \prod_{j = 1}^{d} λ_{j}^{- 1} S [1 : λ_{1}, \dots, 1 : λ_{d}]

(20)

where for the 2D case, S

[a : b, c : d]

represents the subset of rows from a to b and the submatrix of columns from c to d for matrix S.

The BSSFFT technique first converts the samples into grid counts at the estimation grid points using fast linear binning approximation via the FFT. Referring to the works of Wand [32], Raykar et al. [33], the computational complexity of brute force computation is

O (n^{2})

, and the computational complexity of the algorithm described in this paper is only

O (\prod_{i = 1}^{d} λ_{i} log (\prod_{i = 1}^{d} λ_{i}))

. Cao et al. [37] provided a detailed analysis of the computational performance of fast linear binning approximation via the FFT for the KDDE algorithm, which is in line with the expectations regarding computational complexity described in [33]. Moreover, interested readers can refer to the relevant open-source codes, such as KS [48], KDEpy [53], and Vikas [54], to build their own applications. All the tests in this paper were carried out on a server equipped with an Intel(R) Xeon(R) Gold 5117 CPU and 256 GB of memory.

3. Algorithm Verification

A 1D random dataset with a sample size of

n = 3000

was constructed using three normal distributions with a mean of

3^{- 1} φ (0, 0.25) + 3^{- 1} φ (5, 0.5) + 3^{- 1} φ (11, 0.75)

. The Gaussian kernel smoothing estimation algorithm (KS) [55,56] is introduced as a comparison to validate the reliability of the BSS and BSSFFT algorithms.

In Figure 2, the BSS (Curve 2), KS (Curve 3), and BSSFFT (Curves 4–6) methods all exhibit three probability peaks at the same positions, with peak coordinates of 0, 5, and 11, consistent with the theoretical parameters of three normal distributions. Under the same bandwidth conditions, compared to the true pdf (Curve 1), the KS estimate is too smooth and fails to provide indicative information; the BSS estimate approximates the true pdf, demonstrating the correctness of the BSS algorithm; and the BSSFFT results differ slightly from the true pdf, but the difference is insignificant. It arises from the approximation process of obtaining the grid counts on the estimation grid, which utilizes the fast linear binning approximation method with the samples. In Curves 4–6, the results were obtained using the default bandwidth (Curve 4), which closely approximates the true probability density function (pdf). Using a large bandwidth results in an estimate that is larger than the true pdf, due to the small number of estimation grid points. This result is affected by probability density normalization, which inflates the probability values at the estimation points. Conversely, with a small bandwidth—that is, a large number of estimation grid points under the same input sampling conditions—the probability density normalization causes the values at the estimation points to be deflated. By examining the differences in the estimation results with a small bandwidth (Curve 6), it is evident that the zero points in Curve 7, which are derived from the differences in Curve 6, align precisely with the theoretical positions of normal distributions in the artificially generated sample data. This demonstrates that the proposed BSSFFT algorithm effectively reflects the distribution of the true pdf.

Furthermore, a random 2D dataset with a sample size of n= 30,000 is generated using three normal distributions with 10,000 data points. The specific random normal distribution functions are as follows:

f (x) = \frac{1}{3} \sum_{i = 1}^{3} φ (μ_{i}, σ_{i})

(21)

The parameters of each normal distribution of Equation (21) are as follows:

\begin{matrix} μ_{1} = (\begin{matrix} 0 \\ 0 \end{matrix}), μ_{2} = (\begin{matrix} 5 \\ - 5 \end{matrix}), μ_{3} = (\begin{matrix} - 5 \\ 5 \end{matrix}) \\ σ_{1} = (\begin{matrix} 0.75 & - 0.5 \\ - 0.5 & 0.75 \end{matrix}), σ_{i} = (\begin{matrix} 0.75 & 0.5 \\ 0.5 & 0.75 \end{matrix}), i \in \{2, 3\} \end{matrix}

(22)

As shown in Figure 3, the 2D BSS results all show three probability density peaks whose horizontal coordinates are

(0, 0)

(5, - 5)

, and

(- 5, 5)

, which are consistent with the settings of Equation (22). The results of both the BSS and BSSFFT methods are not significantly different from the theoretical values; however, the results of BSSFFT are smoother than those of the BSS method due to the use of fast linear binning approximation. The reliability of the BSS and BSSFFT algorithms can be verified by analyzing the results of the random 1D and 2D normal distribution data.

4. Model Experiment

4.1. Single Cubic Model

A homogeneous cube has analytical solutions and its theoretical structural index is 2 [21,57]. Therefore, a single cube model with a center

(- 1000, - 2000, 1500)

, a size of 2000 m × 2000 m × 2000 m, and a residual density of 0.36 g/cm³ is constructed to further test the adaptability of the B-spline density estimation to Euler solutions and the reliability of the algorithm. The observation height is 50 m; the range of the survey area is x: 0 m∼10,000 m, y: 0 m∼10,000 m; and the size of the survey grid is 100 m × 100 m. The forward data are contaminated by Gaussian random noise with a zero mean and a standard deviation of

σ (\bar{p})

σ (\bar{p}) = (\max (d_{o b s}) - \min (d_{o b s})) \times \bar{p}

(23)

where the min and max functions return the minimum and maximum values of the input data.

Euler deconvolution with a sliding window of 15

\times 15

grid points is used to process the contaminated gravity data (

\bar{p} = 3 %

; see Figure 4) to obtain the Euler solutions, which contain structural indices and spatial positions

(x_{0}, y_{0}, z_{0})

. In this paper, all Euler solutions that are outside the moving window width or the observed measurement data,

N < 0

N > 3

, are eliminated. A total of 4112 Euler solutions remain after filtering. Figure 5 shows the scatter plots of the Euler solutions. The Euler deconvolution yields very pure solutions when there is no noise in the data (see Figure 5a); in contrast, as the Euler equations are easily affected by noise, the Euler solutions tend to diverge when there is noise in the data (see Figure 5b). For this reason, we use an empirical formula proposed by Thompson [4] to cull out spurious solutions [4,6]:

σ_{r} = N \frac{σ}{z_{0}}

(24)

where

z_{0}

is the estimated depth and

σ_{r}

is the standard deviation

z_{0}

of the solution

σ

(i.e., m) relative to the estimated depth. Most of the spurious solutions are eliminated when

σ \leq 10^{- 3}

, as shown in Figure 5c,d.

Regarding the computation time for the analysis of one-dimensional subsets

\{x_{0}\}

\{y_{0}\}

\{z_{0}\}

, and N, the BSS algorithm requires 0.1083, 0.0817, 0.0788, and 0.1011 s, and the BSSFFT algorithm requires 0.05667, 0.0338, 0.0329, and 0.0334 s, respectively. Figure 6 shows the probability curves derived from subsets of Euler solutions using the 1D BSS and BSSFFT methods with the same bandwidth conditions. It is clear from the curves in Figure 6 that both 1D BSS and BSSFFT show multiple peaks in the estimation results for subsets

{x_{0}}

{y_{0}}

{z_{0}}

, and

{N}

. When analyzing Figure 5 and Figure 6 simultaneously, we note that the rightmost peaks of curves 1–2 and 5–6 in Figure 6 are correlated with the x/z coordinates of the anomaly source center, while curves 3–4 show multiple peaks, but it is not possible to determine which peak corresponds to the y coordinate of the anomaly source center, and due to the excessive number of spurious solutions, the peaks correlated with the theoretical value of N are not efficiently extracted. Although the abscissas of some peaks in the curves of Figure 6 are partially consistent with the model settings, the estimation results of the 1D subset all exhibit multiple peaks. This increases the difficulty in manually selecting peaks and is likely to cause ambiguity. The 4D set

{x_{0}, y_{0}, z_{0}, N}

of the single cube model is now arranged into six groups of 2D subsets:

{x_{0}, y_{0}}

{x_{0}, z_{0}}

{y_{0}, z_{0}}

{x_{0}, N_{0}}

{y_{0}, N_{0}}

{z_{0}, N_{0}}

. A total of 12 results are obtained using the 2D BSS and BSSFFT for all Euler solution subsets.

Figure 7 shows the probability density distributions obtained through the 2D BSS and BSSFFT analysis of 2D subsets for a single cubic model. Regarding the computation time for the analysis of the 2D subsets

{x_{0}, y_{0}}

{x_{0}, z_{0}}

{y_{0}, z_{0}}

{x_{0}, N}

{y_{0}, N}

, and

{z_{0}, N}

, the BSS algorithm requires 0.5212, 0.4752, 0.5231, 0.4862, 0.5186, and 0.5102 s, while the BSSFFT algorithm requires 0.0423, 0.0397, 0.0425, 0.0420, 0.0436, and 0.0413 s, respectively. The values of

x_{0}

y_{0}

z_{0}

, and N corresponding to the main probability density peaks of each subplot are 7000 m, 4000 m, 2000 m, and

\approx 1.8

, respectively, due to the presence of a large number of spurious solutions near

z = 0

and

N = 0

in the shallow part of Figure 5. This leads to the appearance of peaks with small tails, which correspond to misplaced solutions, in each subplot, with the exception of Figure 7a,g. For example, see Figure 7c–f,h–i. By comparing the probability densities of the subplots in Figure 7a,g,c,i,e,k with those in Figure 7f,i, it can be seen that BSSFFT effectively eliminates spurious solutions and determines the anomalous source centers by converting the samples into grid counts through fast linear binning approximation. Compared to the 1D results, the 2D results of BSS and BSSFFT show multiple peaks for both the depth and structural index estimation of the model, but only one peak is dominant, providing a basis for the determination of the N and spatial locations of the anomalies. When analyzing the presence of multiple peaks, the equivalence effect [58] is considered, i.e., two adjacent cubes can also produce such anomalies, but there is some ambiguity [37]. When analyzing subsets, the probability density trailing phenomenon is more obvious. Furthermore, the graphical interpretation of the 5D results of the 4D probability density estimation is too complex. To solve these problems, four 3D subsets,

\{x_{0}, y_{0}, z_{0}\}

\{x_{0}, y_{0}, N\}

\{x_{0}, z_{0}, N\}

, and

\{y_{0}, z_{0}, N\}

, are derived from the entire Euler solution set

\{x_{0}, y_{0}, z_{0}, N\}

and analyzed using the 3D BSS and BSSFFT.

Figure 8 and Figure 9 show the Euler solution subsets’ 3D BSS and BSSFFT results, respectively. Through the comparative analysis of these two figures, it can be observed that in the 3D probability density estimation results, with the exception of a subset

\{x_{0}, y_{0}, z_{0}\}

, which has a single probability density peak, the results of all subsets show multiple peaks due to spurious solutions, with the most prominent peak corresponding to the theoretical settings of the model. Compared to the BSS result, it can be seen that the BSSFFT method gives more compact results, with only one probability density peak, and is unaffected by the spurious solutions. Since the probability density values corresponding to spurious solutions are very small, as shown in Figure 9, it is convenient to identify the influence of these spurious solutions. By comparing the 1D, 2D, and 3D probability density results, it can be concluded that the probability density peaks/contours/isosurfaces obtained by BSSFFT can be used to analyze Euler solution subsets, avoid the influence of spurious solutions, and effectively determine anomaly sources. Regarding the computation time for the analysis of the 3D subsets

\{x_{0}, y_{0}, z_{0}\}

\{x_{0}, y_{0}, N\}

\{x_{0}, z_{0}, N\}

, and

\{y_{0}, z_{0}, N\}

, the BSS algorithm requires 29.1144, 26.2609, 23.9917, and 32.2006 s, and the BSSFFT algorithm requires 3.4662, 2.6576, 2.3155, and 3.5115 s, respectively. Therefore, BSSFFT is considered to be superior to the BSS algorithm to some extent.

4.2. Synthetic Model

To further verify the correctness of the proposed BSSFFT algorithm for multiple anomalies, based on a single cube model setup, a synthetic model is studied. It is composed of two cubes, both 1000 m × 1000 m × 1000 m in size, with center-of-mass coordinates of (1500, −1500, 2000) and (1500, −1500, 2000). Euler deconvolution with a sliding window of

15 \times 15

is used to process gravity data contaminated by

3 %

Gaussian noise based on Equation (23) in order to obtain Euler solutions (see Figure 10 and Figure 11).

As shown in Figure 12, by applying 1D BSS and BSSFFT to the Euler solutions derived from the synthetic model, probability density curves are obtained for each subset of the Euler deconvolution solutions. The 1D BSS and BSSFFT results provide accurate estimates for each subset; this can be seen by comparing the horizontal coordinates of each density peak with the theoretical value of the synthetic model, as shown in Table 1. However, similar to the results in Figure 6, the presence of multiple peaks makes it difficult to effectively use 1D probability density curves to analyze Euler solutions in practical applications.

Figure 13 shows the probability density distribution of each subset obtained from the 2D Euler solutions of the combined model using BSS and BSSFFT. By comparing the models’ positions and their structural indices, the 2D BSSFFT results can be used to determine the spatial positions of the anomaly sources and their boundaries, showing good focusing performance compared to the 2D BSS method. However, there is some ambiguity due to the model setup, i.e., both the cubic burial depths and structural indices are the same, resulting in only one peak in Figure 13f, as shown in Table 2. For this reason, 3D Euler solution subsets analyzed with the synthetic model using only the 3D BSSFFT are considered.

Figure 14 shows the 3D BSSFFT density estimation results for the subset of the synthetic model. Comparing the estimation results of the 1D, 2D, and 3D BSS, and BSSFFT, it can be observed that the BSSFFT method produces more compact results, with only one density peak unaffected by the instabilities of the Euler homogeneity equation. The results are remarkably stable, with very few instances of probability density tails [7] that correspond to spurious solutions.

From Table 1, Table 2 and Table 3, it can be seen that the results of the probability density estimation for the two cubes are slightly shifted from the theoretical values. This phenomenon is visualized in Figure 14, particularly in Figure 14a. Through a comparison with Figure 9, it can be hypothesized that the mutual interference of adjacent anomalous sources causes this.

A comprehensive analysis of the above results indicates that the BSSFFT can effectively localize anomaly sources. Therefore, in our subsequent research, only the BSSFFT method proposed in this paper will be used to analyze Euler solutions.

4.3. Sensitivity of 3D BSSFFT to Separations

To evaluate the performance of the proposed algorithms in distinguishing neighboring anomaly sources, three models presented are used in the work by Cao et al. [37], each consisting of two cubes of the same size (side length of 2.0 (km)) and different densities (2.36 g/cm³ and 1.27 g/cm³ for the left and right cubes, respectively), at separation L = 4.0 (km), 2.5 (km), and 1.0 (km), respectively. These values of separation L allow the two sources to approach each other simultaneously in the x and y directions. Table 4 lists the two cubes’ parameters. A subset of

\{x_{0}, y_{0}, z_{0}\}

, obtained from FTG data and contaminated by

3 %

Gaussian noise based on Equation (23), is used in this 3D BSSFFT study.

As the separation L decreases, it becomes more difficult to distinguish between two adjacent cubes using the Euler solutions in Figure 15d–f, even with Thompson’s empirical formula [4], as shown in Figure 15g–i. Moreover, as seen in Figure 15j–l, similar to the fuzzy C-means clustering algorithm (FCM) [37], the k-means clustering algorithm (K-means) depends on the number of predefined clusters; thus, it can effectively discriminate between adjacent anomaly sources but cannot eliminate spurious solutions. On the other hand, some spurious solutions, such as the green Euler solutions in Figure 15m–o, can be found by using the density-based spatial clustering of applications with noise (DBSCAN) algorithm, when the distance between nearby anomaly sources is relatively large. Regarding the computation time for different values of L (4.0 (km), 2.5 (km), and 1.0 (km)), the K-means algorithm requires 0.1223, 0.1927, and 0.2839 s, and the DBSCAN algorithm requires 69.2537, 64.1225, and 72.1126 s, respectively.

Figure 16 illustrates the separation of the anomaly sources by using 3D BSSFFT to estimate the subset

\{x_{0}, y_{0}, z_{0}\}

. Unlike other techniques used to separate adjacent anomaly sources, the two probability density peaks directly distinguish two sources at different separations L, as shown in Figure 15.

4.4. Sensitivity of 3D BSSFFT to Gaussian Noise

If the signal-to-noise ratio of the input data is low, the Euler solutions will tend to diverge. Many strategies/methods have been proposed to eliminate spurious solutions and determine the optimal solutions [4,7,43]. For example, clustering methods cannot effectively distinguish spurious solutions from shallow sources and neighboring bodies [15,59]. Therefore, a model containing two cubes with a separation L of 2.5 (km) (see Figure 16) is selected, and the corresponding Euler solutions are derived from the corresponding forward data, which are contaminated with

\bar{p} =

0%, 4%, and 8% noise based on Equation (23). On this basis, the sensitivity of the BSSFFT method to different levels of Gaussian noise is analyzed.

As shown in Figure 17, the contours in the gravity map become increasingly distorted as the percentage of noise increases, and the number of Euler solutions with low SIs increases dramatically. However, this is not obvious in Figure 17b,e,h or Figure 17c,g,k, due to the fact that the previously drawn solutions overlap with the later ones, and backward solutions are obscured by forward ones. Compared to the results in Cao et al. [37], as we use a series of sliding windows of different sizes to generate the Euler solutions, the estimated positions derived from the 3D probability density distribution are very close to the theoretical positions of the synthetic model, despite the percentage of noise being equal to 8%. However, the estimated positions are all far from the center of the anomaly sources due to the influence of neighboring anomaly sources.

4.5. Sensitivity of 2D Euler Deconvolution to Gaussian Noise

The two cubes described in Table 4 were placed on the y-axis to construct four models to verify the correctness of the 2D Euler deconvolution algorithm, as shown in Figure 18 and Figure 19.

The Euler solutions are obtained by means of Euler deconvolution in conjunction with the gravity analytic solution of the rectangular, as shown in Figure 19. When the data are free of noise, in Figure 19a,d,g,j, the Euler deconvolution yields “pure” Euler solutions for a single cube model; for two cube models, as L decreases, it is clearly seen that most of the Euler solutions are in the tails and have deviated from the anomaly source. When the data are contaminated with noise (

p = 4 %

), the 2D Euler deconvolution can no longer accurately obtain the optimal solutions, and only part of the Euler solution is inside the anomalous source; however, when the two anomalous sources are close to each other (L = 1.0 (km)), the 2D Euler deconvolution has detected a deeper anomalous source, but the Euler solution is farther away from its center.

Compared with the results of 2D Euler deconvolution (see Figure 19), the results of 2D BSSFFT can locate the anomaly source more clearly; compared with 3D BSSFFT’s results (see Figure 16 and Figure 17d,h,l), the estimation results of 2D BSSFFT may deviate from the center of the anomaly source when the noise level is high or when adjacent anomaly sources are close to each other, as shown in Figure 20.

5. Bishop Model

To improve the algorithm’s credibility, we introduce a widely used model, i.e., the Bishop model. It is a 3D synthetic basement model [60], which was first proposed by Fairhead et al. [61]. The terrain data used are derived from a part of the digital elevation model in the volcanic highlands of northern Bishop, California, covering an area of about

10.5

km ×

10.5

km, scaled and shifted in the depth direction, i.e., has been upscaled by a factor of 30 in the x, y, and z dimensions and then subtracted from a fixed value, resulting in the surface layer being translated below sea level [62,63]. Then, a 3D basement model with a horizontal grid size of 500 m × 500 m is created using the terrain data [63]. Because there are two large offset faults in this area [64], the elongated deep tectonic low pressure and other minor faults are formed from west to east and from north to south, and the top structure of the magnetic base layer similar to the basement is further formed [64,65].

Based on Bishop’s 3D basement model, Fairhead et al. [61] used the Northwest Geophysical GM-SYS 3D software (version 4.10) to calculate and obtain the gravity and magnetic field forward results, thus constructing the Bishop 5X dataset [66,67]. Because it is a shared complex gravity/magnetic anomaly dataset, it is welcomed and studied by many researchers in the field of geophysics [60,61,65,66,67,68,69,70].

The Bishop 5X dataset contains only gravity anomalies derived from strata and complex magnetic anomalies. Therefore, to obtain the corresponding pseudo-gravity anomaly, Ekinci and Yigitbas [71] used pseudo-gravity transformation to convert magnetic data into pseudo-gravity data with gravity data characteristics to assist geological interpretation [71,72,73]. Pseudo-gravity transformation [72] used Fourier technology [74] as a “bridge” to integrate the grid data of total magnetic field intensity to obtain pseudo-gravity grid data [71,75]. The pseudo-gravity data can reduce the advantages of the shallow magnetic source and enhance the magnetic anomaly related to the deep magnetic source. Because pseudo-gravity data represent intuitive depictions of the underground magnetic field, they are not directly related to gravity anomalies or density. Essentially, they pertain to magnetic anomalies, which can only be explained by the distribution of magnetization [76,77]. Therefore, the data are rarely valued and applied in practical exploration problems [72,78]. Thus, we apply the Euler deconvolution Equations (1) and (2) to the magnetic data of the Bishop model. Considering that the data point interval of the Bishop 5X dataset is much smaller than the magnetic sources’ size, we use 1000 m × 1000 m to re-grid the Bishop 5X magnetic grid data.

Considering that it is difficult to obtain the background field information of an input field data by Euler deconvolution, we use correlation analysis to determine the background field [37,79,80,81] and obtain the corresponding residual magnetic anomaly field (i.e, the residual field between the magnetic anomaly and the background field, as shown in Figure 21c,d) for the magnetic anomaly of Bishop 5X (the magnetic inclination

i = 45^{\circ}

and

i = 90^{\circ}

), as shown in Figure 21. To distinguish between gravity data and magnetic data in this paper, the marker g of Equations (1) and (2) is changed to T for magnetic data. Then, the residual magnetic anomaly (also denoted as

T_{z}

in this paper) is converted into corresponding first-order derivatives (

T_{z x}

T_{z y}

, and

T_{z z}

), as shown in Figure 22 and Figure 23. It can be seen that there is no data distortion or oscillation in each subgraph.

Here, it is assumed that the magnetic sources’ spatial distribution information is unknown. For this reason, we use Equations (1) and (2) in 2D and 3D space with different-sized sliding windows (

w_{x} = 10 \sim 35

) to obtain the corresponding Euler solutions. For the 2D application, a red survey line spanning three sources of anomalies is selected, and

T_{z}

and its corresponding derivatives are obtained by Fourier transform [82], as shown in Figure 22a and Figure 23a. By comparing and analyzing Figure 24 and Figure 25, it can be seen that both Euler solutions show a relatively scattered trend, and only the optimal solutions may be seen to show a linear aggregation trend. However, under the interference of spurious solutions, it is difficult to determine optimal solutions and use them to locate the source of anomalies. In contrast, the BSSFFT’s results are effective in identifying the anomaly source with comparable depth resolution. Compared to Figure 25c, there is a probability density peak interference in the shallow part of Figure 24c, which we hypothesize is related to the interference caused by the magnetic inclination.

In the 3D case, for

i = 45^{\circ}

and

90^{\circ}

, the computation time required for the Euler deconvolution to traverse the entire survey area is 572.25 s and 601.39 s, respectively. A total of 237,292 and 286,370 Euler solutions are obtained (filtered by

N < 0.25

N > 3.0

), as shown in Figure 26 and Figure 27. A large number of Euler solutions with small N values are scattered in the shallow part, which are likely to be spurious solutions. The structural index N increases continuously as depth z increases, and when

z > 15, 000

m, Euler solutions (

N > 2.5

) are somewhat distorted. In Figure 26b,d and Figure 27b,d, with the help of geological boundary lines, some Euler solutions can mark magnetic sources. However, these Euler solutions cannot be distinguished from other scattered Euler solutions.

For this reason, we use the BSSFFT algorithm with an estimated grid of

100 \times 100 \times 50

to perform probability density analysis on these Euler solutions in Figure 26 and Figure 27. For

i = 45^{\circ}

and

90^{\circ}

, the calculation times are

3.0502

and

3.1829

s, respectively.

The different depth probability density slices of the Euler solutions for the Bishop 5X magnetic anomaly with

i = 45^{\circ}

and

90^{\circ}

are shown in Figure 28 and Figure 29, showing multiple density peaks in these subfigures. Peak A consists of several sub-peaks surrounding adjacent anomalous sources, which, combined with stratigraphic information, suggests that peak A may be associated with adjacent faults and strata; peak B extends from

z = 4500

m down to 13,500 m, but fails to accurately indicate the corresponding magnetic source; peak D extends from

z = 4500

m down to

z = 7500

m, but occurs only at the southwestern end of the corresponding anomalous source and affects a limited depth range. In contrast to the other peaks associated with the corresponding anomalous sources, peak C extends from

z = 1500

m down to

z = 7500

m, pinpointing the source at

z = 4500 \sim 6000

m (at the maximum probability density value).

The correlation between the probability density peaks in Figure 29, the geological structure boundary, and the magnetic anomaly sources is similar to that in Figure 28. However, peaks B and E are different. In Figure 29, peak B extends downward from

z = 3000

m to

z = 7500

m, showing two separate peaks, similar to the two peaks of the cylinder model affected by the equivalent effect described in reference Cao et al. [37]. Compared to Figure 29, Figure 28 shows a large number of very small amplitude probability dense peaks that fail to dominate any subgraph. It can be concluded that these peaks are due to the boundary effect of the survey area and the influence of magnetic inclination.

Figure 30 and Figure 31 show that the corresponding probability density isosurfaces of the magnetic anomaly data obtained by the proposed method can effectively determine the anomaly source’s spatial distributions while avoiding the influence of magnetic inclination (see Figure 29). Referring to Cao et al. [37], more focused results will be obtained if smaller bandwidths or larger estimation grids are used.

6. Discussion

In practice, violations of the Euler homogeneity conditions (e.g., inappropriate data grid spacing), low signal-to-noise ratios (either natural or FFT-induced) [7,9,42,83,84], and source overlay have resulted in many spurious solutions, leading to a sparsely dispersed distribution of the entire Euler solution set, as opposed to a dense distribution [9,83,84,85]. When interpreting Euler solutions, the abundance of solutions generated by Euler deconvolution leads to rendering issues where previously drawn solutions overlap with later ones, and backward solutions are obscured by forward ones. Therefore, the end goal of the processing and interpretation of Euler solutions is not merely to reject spurious solutions [39,42,44,84]. To overcome the above problems, we propose the BSS and BSSFFT methods based on probability density estimates, which are used to apply Euler solutions for more intuitive geological interpretation.

The apparent N may appear to be different due to numerical issues, but the true value is always constant [37,40]. Therefore, in practical applications, Euler deconvolution uses a series of predefined structural indices, which leads to the presence of many solutions in the Euler deconvolution process [6,10]. Since it is challenging to determine N for the geological target of interest [41,42], this work adopts the Euler deconvolution of gravity/magnetic data with an unprescribed structural index to obtain the Euler solutions. Furthermore, as fixed-size sliding windows cannot effectively identify anomaly sources of varying sizes, the Euler deconvolution of gravity data in conjunction with differently sized sliding windows is applied to the field data.

For a given sample, the bandwidth is inversely proportional to the grid size. The bandwidth determines the level of smoothness in the resulting probability density curve, reflecting the proportion of observed data points contributing to the formation of the curve. A larger bandwidth results in a flatter curve due to a lesser contribution of individual data points to the final shape; conversely, a smaller bandwidth leads to a steeper curve because of a greater contribution of the data points. Therefore, the quality of the BSS’s final probability density distribution is influenced by the bandwidth parameter or the grid size. If the sample size is not considered, the memory consumption of the BSS and BSSFFT algorithms depends solely on the size of the estimation grid and the peak memory usage of the n-dimensional FFT. Therefore, the memory consumption of BSSFFT is higher than that of BSS.

Increasing the grid size significantly reduces the computational efficiency of BSS, particularly for large data samples. Furthermore, combining sample data with high-dimensional sparse data, especially in scenarios involving multiple-density clusters, complicates the efficient identification of the “best” bandwidth [86].

In this study, to effectively overcome these challenges associated with increased computational efficiency, BSSFFT is implemented using fast linear binning approximation and FFT technology, which significantly improve the computational efficiency in the estimation of BSS [37]. Moreover, experimenting with grids of different sizes or selecting isosurfaces from different formations could be more effective in identifying interesting geological structures.

For adjacent anomaly sources, the results obtained using FCM, DBSCAN, and probability density were compared and analyzed, and 3D probability density contours/isosurfaces could effectively separate adjacent anomaly sources [37]. Due to the fast linear binning approximation and FFT techniques, we can apply different grid sizes and/or select different formations of isosurfaces to identify interesting geological structures more efficiently. For 3D applications, the 1D or 2D subsets produce too many probability density distribution curves or images, making simultaneous interpretation difficult. The 5D results

\{x_{0}, y_{0}, z_{0}, N, p - v a l u e\}

are obtained using BSSFFT for Euler solutions but are challenging to demonstrate.

In this work, we use BSSFFT to analyze the subsets

\{x_{0}, z_{0}\}

(2D case) and

\{x_{0}, y_{0}, z_{0}\}

(3D case) obtained by the conventional Euler deconvolution with some success. By analyzing the 4D probability density distribution, we can obtain information with a good depth and lateral resolution for the analysis of the geological structures in the study area.

7. Conclusions

This study addresses the optimization of the computational efficiency problem found in a previous work. Inspired by the work of Wand [32], fast linear binning approximation with FFT technology is applied, which significantly improves the efficiency in probability density estimation. The current study uses fast linear binning approximation and FFT technology to provide a quick and efficient BSSFFT algorithm. Then, by analyzing the computational accuracy of KS, BSS, and BSSFFT, it is concluded that BSSFFT has high computational accuracy. Supported by the fast linear binning approximation technique, BSSFFT demonstrates strong focusing capabilities.

The procedure of the BSSFFT method for representing geological bodies using Euler solutions, as proposed in this paper, bears similarities to the KDDE algorithm process. It involves selecting moving windows of varying or fixed sizes and using Equation (2) to traverse grid data for Euler solutions. An estimation grid is constructed with a predefined size

λ

or bandwidth h. The subset

\{x_{0}, y_{0}, z_{0}\}

is projected onto this grid using a fast linear binning approximation to produce a grid count C. Using Equation (15) to perform the kernel function evaluation, the process constructs zero-padded versions of the grid count and the kernel evaluation. Fast convolution between these elements via the FFT yields the probability density distribution of the subset

\{x_{0}, y_{0}, z_{0}\}

The tests with Bishop 5X data and synthetic models show that the BSSFFT approach presented in this paper, in conjunction with subset

\{x_{0}, y_{0}, z_{0}\}

derived from the Euler deconvolution, achieves meaningful geological results that are independent of pre-existing geological data.

Author Contributions

Conceptualization, S.C. and G.L.; data curation, S.C., Z.M. and B.Y.; formal analysis, S.C.; funding acquisition, S.C. and G.L.; investigation, S.C., B.Y. and X.C.; methodology, S.C.; project administration, X.C. and G.L.; resources, G.L. and X.C.; software, S.C., P.C., Z.M. and B.Y.; supervision, G.L. and X.C.; validation, S.C., P.C. and Z.M.; visualization, S.C.; writing—original draft preparation, S.C., P.C. and Z.M.; writing—review and editing, S.C. and P.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China under grant 41704138 and grant 41974148, in part by the Hunan Provincial Science & Technology Department of China under grant 2017JJ3069, in part by the Project of the Doctoral Foundation of Hunan University of Science and Technology under Grant E51651, and in part by the Hunan Provincial Key Laboratory of Share Gas Resource Exploitation under grant E21722.

Data Availability Statement

Data are contained within the article.

Acknowledgments

We thank Gerovska and Araúzo-Bravo [6] for releasing the source code of the traditional Euler deconvolution with an unprescribed structural index. We would like to thank the three anonymous reviewers for their helpful comments and suggestions, which allowed us to further improve the quality of the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

KDDE	kernel density derivative estimation
BSS	B-spline probability density estimation method
BSSFFT	B-spline probability density estimation method based on fast Fourier transform
KS	Gaussian kernel smoothing estimation
FFT	fast Fourier transform
pdf	probability density function
DBSCAN	density-based spatial clustering of applications with noise
K-means	k-means clustering algorithm
FCM	fuzzy C-means clustering algorithm

References

Smellie, D. Elementary approximations in aeromagnetic interpretation. Geophysics 1956, 21, 1021–1040. [Google Scholar] [CrossRef]
Hood, P. Gradient measurements in aeromagnetic surveying. Geophysics 1965, 30, 891–902. [Google Scholar] [CrossRef]
Choudhury, S.; Amaravadi, R. The Direct Approach to Magnetic Interpretation and Its Practical Application; discussion and reply. Geophysics 1972, 37, 181–182. [Google Scholar] [CrossRef]
Thompson, D. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data. Geophysics 1982, 47, 31–37. [Google Scholar] [CrossRef]
Reid, A.; Allsop, J.; Granser, H.; Millett, A.; Somerton, I. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics 1990, 55, 80–91. [Google Scholar] [CrossRef]
Gerovska, D.; Araúzo-Bravo, M.J. Automatic interpretation of magnetic data based on Euler deconvolution with unprescribed structural index. Comput. Geosci. 2003, 29, 949–960. [Google Scholar] [CrossRef]
FitzGerald, D.; Reid, A.; McInerney, P. New discrimination techniques for Euler deconvolution. Comput. Geosci. 2004, 30, 461–469. [Google Scholar] [CrossRef]
Huang, L.; Zhang, H.; Li, C.F.; Feng, J. Ratio-Euler deconvolution and its applications. Geophys. Prospect. 2022, 70, 1016–1032. [Google Scholar] [CrossRef]
Farrelly, B. What is Wrong with Euler Deconvolution? In Proceedings of the 59th EAGE Conference & Exhibition, Geneva, Switzerland, 26–30 May 1997; European Association of Geoscientists & Engineers: Utrecht, The Netherlands, 1997; pp. 1–2. [Google Scholar] [CrossRef]
Melo, F.F.; Barbosa, V.C. Reliable Euler deconvolution estimates throughout the vertical derivatives of the total-field anomaly. Comput. Geosci. 2020, 138, 104436. [Google Scholar] [CrossRef]
Beiki, M. TSVD analysis of Euler deconvolution to improve estimating magnetic source parameters: An example from the Åsele area, Sweden. J. Appl. Geophys. 2013, 90, 82–91. [Google Scholar] [CrossRef]
FitzGerald, D.; Milligan, P. Defining a deep fault network for Australia, using 3D “worming”. ASEG Ext. Abstr. 2013, 2013, 1–4. [Google Scholar] [CrossRef]
Agarwal, B.; Srivastava, S. Analyses of self-potential anomalies by conventional and extended Euler deconvolution techniques. Comput. Geosci. 2009, 35, 2231–2238. [Google Scholar] [CrossRef]
Keating, P.; Pilkington, M. Euler deconvolution of the analytic signal and its application to magnetic interpretation. Geophys. Prospect. 2004, 52, 165–182. [Google Scholar] [CrossRef]
Mikhailov, V.; Galdeano, A.; Diament, M.; Gvishiani, A.; Agayan, S.; Bogoutdinov, S.; Graeva, E.; Sailhac, P. Application of artificial intelligence for Euler solutions clustering. Geophysics 2003, 68, 168–180. [Google Scholar] [CrossRef]
Goussev, S.A.; Peirce, J.W. Magnetic basement: Gravity-guided magnetic source depth analysis and interpretation. Geophys. Prospect. 2010, 58, 321–334. [Google Scholar] [CrossRef]
Cao, S.; Zhu, Z.; Lu, G. Gravity tensor Euler Deconvolution solutions based on adaptive fuzzy cluster analysis. J. Central South Univ. 2012, 43, 1033–1039. (In Chinese) [Google Scholar]
Zhang, H.; Wang, Q.; Shi, W.; Hao, M. A Novel Adaptive Fuzzy Local Information C -Means Clustering Algorithm for Remotely Sensed Imagery Classification. IEEE Trans. Geosci. Remote Sens. 2017, 55, 5057–5068. [Google Scholar] [CrossRef]
Husson, E.; Guillén, A.; Séranne, M.; Courrioux, G.; Couëffé, R. 3D Geological modelling and gravity inversion of a structurally complex carbonate area: Application for karstified massif localization. Basin Res. 2018, 30, 766–782. [Google Scholar] [CrossRef]
Lee, J.H.; Kim, D.H.; Chung, C.W. Multi-dimensional selectivity estimation using compressed histogram information. SIGMOD Rec. 1999, 28, 205–214. [Google Scholar] [CrossRef]
Cao, S.; Deng, Y.; Yang, B.; Lu, G.; Zhu, Z.; Chen, P.; Xie, J.; Chen, X. 3D Probability Density Imaging of Euler Solutions using Gravity Data: A Case Study of Mount Milligan. Acta Geophys. 2024, 1, 1–21. [Google Scholar] [CrossRef]
Chen, J.; Chen, H.; Chen, Q.; Song, X.; Wang, H. Vessel sailing route extraction and analysis from satellite-based AIS data using density clustering and probability algorithms. Ocean Eng. 2023, 280, 114627. [Google Scholar] [CrossRef]
Peng, J.; Wang, X.; Zhao, H.; Dong, Z. Single-sample unmixing and parametric end-member modelling of grain-size distributions with transformed probability density functions and their performance comparison using aeolian sediments. Sediment. Geol. 2023, 445, 106328. [Google Scholar] [CrossRef]
Gu, J.; Zhang, J.; Chen, L.; Cai, Z. An isogeometric BEM using PB-spline for 3-D linear elasticity problem. Eng. Anal. Bound. Elem. 2015, 56, 154–161. [Google Scholar] [CrossRef]
Eilers, P.H.; Marx, B.D. Flexible smoothing with B-splines and penalties. Stat. Sci. 1996, 11, 89–121. [Google Scholar] [CrossRef]
Faenza, L.; Marzocchi, W.; Boschi, E. A non-parametric hazard model to characterize the spatio-temporal occurrence of large earthquakes; an application to the Italian catalogue. Geophys. J. Int. 2003, 155, 521–531. [Google Scholar] [CrossRef]
Liao, J.; Liu, H.; Li, W.; Guo, Z.; Wang, L.; Wang, H.; Peng, S.; Hursthouse, A.S. 3-D Butterworth Filtering for 3-D High-density Onshore Seismic Field Data. J. Environ. Eng. Geophys. 2018, 23, 223–233. [Google Scholar] [CrossRef]
Xiao, W.; Ma, H.; Zhou, L.; Li, H. Adaptive Fuzzy Fixed-Time Formation-Containment Control for Euler-Lagrange Systems. IEEE Trans. Fuzzy Syst. 2023, 31, 3700–3709. [Google Scholar] [CrossRef]
Mautz, R.; Ping, J.; Heki, K.; Schaffrin, B.; Shum, C.; Potts, L. Efficient spatial and temporal representations of global ionosphere maps over Japan using B-spline wavelets. J. Geod. 2005, 78, 662–667. [Google Scholar] [CrossRef]
Herrmann, F.J.; Hennenfent, G. Non-parametric seismic data recovery with curvelet frames. Geophys. J. Int. 2008, 173, 233–248. [Google Scholar] [CrossRef]
Li, C.; Wen, X.; Liu, X.; Zu, S. Simultaneous Seismic Data Interpolation and Denoising Based on Nonsubsampled Contourlet Transform Integrating With Two-Step Iterative Log Thresholding Algorithm. IEEE Trans. Geosci. Remote Sens. 2022, 60, 1–10. [Google Scholar] [CrossRef]
Wand, M. Fast computation of multivariate kernel estimators. J. Comput. Graph. Stat. 1994, 3, 433–445. [Google Scholar] [CrossRef]
Raykar, V.C.; Duraiswami, R.; Zhao, L.H. Fast Computation of Kernel Estimators. J. Comput. Graph. Stat. 2010, 19, 205–220. [Google Scholar] [CrossRef]
Wang, B.; Krebes, E.S.; Ravat, D. High-precision potential-field and gradient-component transformations and derivative computations using cubic B-splines. Geophysics 2008, 73, I35–I42. [Google Scholar] [CrossRef]
Togbenou, K.; Xiang, H.Y.; Li, Y.; Chen, N. Improved Spectral Representation Method for the Simulation of Stochastic Wind Velocity Field Based on FFT Algorithm and Polynomial Decomposition. J. Eng. Mech. 2018, 144, 04017171. [Google Scholar] [CrossRef]
Fang, B.; Chen, S.; Dong, Z. Density Distillation for Fast Nonparametric Density Estimation. IEEE Trans. Neural Netw. Learn. Syst. 2022, 34, 9424–9438. [Google Scholar] [CrossRef] [PubMed]
Cao, S.; Deng, Y.; Yang, B.; Lu, G.; Hu, X.; Mao, Y.; Hu, S.; Zhu, Z. Kernel Density Derivative Estimation of Euler Solutions. Appl. Sci. 2023, 13, 1784. [Google Scholar] [CrossRef]
Harfouche, L.; Zougab, N.; Adjabi, S. Multivariate generalised gamma kernel density estimators and application to non-negative data. Int. J. Comput. Sci. Math. 2020, 11, 137–157. [Google Scholar] [CrossRef]
Reid, A.B. Euler deconvolution: Past, present and future—A review. SEG Tech. Prog. Expand. Abstr. 1995, 03, 272–273. [Google Scholar] [CrossRef]
Ravat, D. Analysis of the Euler method and its applicability in environmental magnetic investigations. J. Environ. Eng. Geophys. 1996, 1, 229–238. [Google Scholar] [CrossRef]
Melo, F.F.; Barbosa, V.C. Correct structural index in Euler deconvolution via base-level estimates. Geophysics 2018, 83, J87–J98. [Google Scholar] [CrossRef]
Reid, A.B.; Thurston, J.B. The structural index in gravity and magnetic interpretation: Errors, uses, and abuses. Geophysics 2014, 79, J61–J66. [Google Scholar] [CrossRef]
Barbosa, V.C.; Silva, J.B.; Medeiros, W.E. Stability analysis and improvement of structural index estimation in Euler deconvolution. Geophysics 1999, 64, 48–60. [Google Scholar] [CrossRef]
Ugalde, H.; Morris, W.A. Cluster analysis of Euler deconvolution solutions: New filtering techniques and geologic strike determination. Geophysics 2010, 75, L61–L70. [Google Scholar] [CrossRef]
Pan, Q.; Liu, D.; Feng, S.; Feng, M.; Fang, H. Euler deconvolution of the analytic signals of the gravity gradient tensor for the horizontal pipeline of finite length by horizontal cylinder calculation. J. Geophys. Eng. 2017, 14, 316–330. [Google Scholar] [CrossRef]
Gehringer, K.R.; Redner, R.A. Nonparametric probability density estimation using normalized b-splines. Commun. Stat. Simul. C 1992, 21, 849–878. [Google Scholar] [CrossRef]
Schumaker, L. Spline Functions: Basic Theory; Cambridge University Press: Cambridge, UK, 2007; pp. 35–42. [Google Scholar] [CrossRef]
Gramacki, A.; Gramacki, J. FFT-based fast computation of multivariate kernel density estimators with unconstrained bandwidth matrices. J. Comput. Graph. Stat. 2017, 26, 459–462. [Google Scholar] [CrossRef]
Chacón, J.E.; Duong, T. Multivariate Kernel Smoothing and Its Applications; CRC Press: Boca Raton, FL, USA, 2018; pp. 34–37. [Google Scholar] [CrossRef]
Rao, S. Interpolation Models. In The Finite Element Method in Engineering, 6th ed.; Butterworth-Heinemann: Oxford, UK, 2017; Section 3; pp. 81–127. [Google Scholar] [CrossRef]
Teukolsky, S.A.; Flannery, B.P.; Press, W.H.; Vetterling, W.T. Numerical Recipes in C: The Art of Scientific Computing; Cambridge University Press: Cambridge, UK, 1992; pp. 100–135. [Google Scholar] [CrossRef]
Arndt, J. Matters Computational: Ideas, Algorithms, Source Code; Springer Science and Business Media: Berlin/Heidelberg, Germany, 2010; pp. 153–170. [Google Scholar] [CrossRef]
Odland, T. tommyod/KDEpy: Kernel density estimation in python. Zenodo 2018, 8, 45–48. [Google Scholar] [CrossRef]
Raykar, V.C.; Duraiswami, R. Fast optimal bandwidth selection for kernel density estimation. In Proceedings of the 2006 SIAM International Conference on Data Mining, Bethesda, MD, USA, 20–22 April 2006; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2006; pp. 524–528. [Google Scholar] [CrossRef]
Duong, T. ks: Kernel density estimation and kernel discriminant analysis for multivariate data in R. J. Stat. Softw. 2007, 21, 1–16. [Google Scholar] [CrossRef]
Botev, Z.I.; Grotowski, J.F.; Kroese, D.P. Kernel density estimation via diffusion. Anal. Stat. 2010, 38, 2916–2957. [Google Scholar] [CrossRef]
Xu, J.C.; Li, W.Y.; Liu, Y.X. Application of euler deconvolution method in airborne gravity exploration. Prog. Phys. 2016, 31, 390–395. (In Chinese) [Google Scholar]
Li, Y.; Oldenburg, D.W. 3-D inversion of gravity data. Geophysics 1998, 63, 109–119. [Google Scholar] [CrossRef]
Gvishiani, A.D.; Mikhailov, V.O.; Agayan, S.M.; Bogoutdinov, S.R.; Graeva, E.M.; Diament, M.; Galdeano, A. Artificial intelligence algorithms for magnetic anomaly clustering. Izvestiya Phys. Solid Earth 2002, 38, 545–559. [Google Scholar]
Williams, S.E.; Fairhead, J.D.; Flanagan, G. Comparison of grid Euler deconvolution with and without 2D constraints using a realistic 3D magnetic basement model. Geophysics 2005, 70, L13–L21. [Google Scholar] [CrossRef]
Fairhead, J.D.; Williams, S.E.; Flanagan, G. Testing Magnetic Local Wavenumber Depth Estimation Methods using a Complex 3D Test Model. SEG Tech. Prog. Expand. Abstr. 2004, 10, 742–745. [Google Scholar] [CrossRef]
Florio, G. ITRESC: A fast and efficient method to recover the basement morphology from potential fields data. SEG Tech. Prog. Expand. Abstr. 2018, 08, 1415–1419. [Google Scholar] [CrossRef]
Florio, G. Mapping the Depth to Basement by Iterative Rescaling of Gravity or Magnetic Data. J. Geophys. Res. Solid Earth 2018, 123, 9101–9120. [Google Scholar] [CrossRef]
Salem, A.; Green, C.; Campbell, S.; Fairhead, J. A Practical Approach to 3D Inversion of Pseudo-gravity for Depth to Basement Mapping—A Test Using the Bishop Model. In Proceedings of the 74th EAGE Conference and Exhibition incorporating EUROPEC 2012, Copenhagen, Denmark, 4–7 June 2012; p. P338. [Google Scholar] [CrossRef]
Reid, A. Hybrid Euler magnetic basement depth estimation: Bishop 3D tests. SEG Tech. Prog. Expand. Abstr. 2005, 24, 6–11. [Google Scholar] [CrossRef]
Gerovska, D.; Araúzo-Bravo, M.J.; Whaler, K.A.; Stavrev, P.; Reid, A.B. Three-dimensional interpretation of magnetic and gravity anomalies using the finite-difference similarity transform. Geophysics 2010, 75, 1JA–Z98. [Google Scholar] [CrossRef]
Dwivedi, D.D.; Chamoli, A. Source Edge Detection of Potential Field Data Using Wavelet Decomposition. Pure Appl. Geophys. 2021, 178, 919–938. [Google Scholar] [CrossRef]
Salem, A.; Williams, S.E.; Fairhead, D.J.; Smith, R.S.; Ravat, D. Interpretation of magnetic data using tilt-angle derivatives. Geophysics 2008, 73, 1–10. [Google Scholar] [CrossRef]
Zhou, W.; Nan, Z.Y.; Li, J. Self-Constrained Euler Deconvolution Using Potential Field Data of Different Altitudes. Pure Appl. Geophys. 2016, 173, 2073–2085. [Google Scholar] [CrossRef]
Li, X. Terracing gravity and magnetic data using edge-preserving smoothing filters. Geophysics 2016, 81, G37–G43. [Google Scholar] [CrossRef]
Ekinci, Y.L.; Yigitbas, E. A geophysical approach to the igneous rocks in the Biga Peninsula (NW Turkey) based on airborne magnetic anomalies: Geological implications. Geodin. Acta 2012, 25, 267–285. [Google Scholar] [CrossRef]
Baranov, V. A new method for interpretation of aeromagnetic maps: Pseudo-gravimetric anomalies. Geophysics 1957, 22, 359–382. [Google Scholar] [CrossRef]
Oni, O.A.; Aizebeokhai, A.P. Aeromagnetic data processing using MATLAB. IOP Conf. Ser. Earth Environ. Sci. 2022, 993, 012017. [Google Scholar] [CrossRef]
Hildenbrand, T.G. FFTFIL; a Filtering Program Based on Two-Dimensional Fourier Analysis of Geophysical Data; U.S. Geological Survey: Reston, VA, USA, 1983; pp. 23–37. [CrossRef]
Grauch, V. Limitations of determining density or magnetic boundaries from the horizontal gradient of gravity or pseudogravity data. Geophysics 1987, 52, 118–121. [Google Scholar] [CrossRef]
Bott, M.H.P.; Smith, R.A.L.; Stacey, R.A. Estimation of the direction of magnetization of a body causing a magnetic anomaly using a pseudo-gravity transformation. Geophysics 1966, 31, 803–811. [Google Scholar] [CrossRef]
Salem, A.; Green, C.M.; Cheyney, S.; Fairhead, J.D.; Aboud, E.; Campbell, S. Mapping the depth to magnetic basement using inversion of pseudogravity: Application to the Bishop model and the Stord Basin, northern North Sea. Interpretation 2014, 2, 1M–T127. [Google Scholar] [CrossRef]
Pratt, D.A.; Shi, Z. An improved pseudo-gravity magnetic transform technique for investigation of deep magnetic source rocks. ASEG Ext. Abstr. 2004, 2004, 1–4. [Google Scholar] [CrossRef]
Zeng, H.; Xu, D.; Tan, H. A model study for estimating optimum upward-continuation height for gravity separation with application to a Bouguer gravity anomaly over a mineral deposit, Jilin province, northeast China. Geophysics 2007, 72, I45–I50. [Google Scholar] [CrossRef]
Montaj, G. The Core Software Platform for Working with Large Volume Gravity and Magnetic Spatial Data; Geosoft Inc.: Toronto, ON, Canada, 2008. [Google Scholar]
Setiadi, I.; Marjiyono; Nainggolan, T.B. Gravity data analysis based on optimum upward continuation filter and 3D inverse modelling (Case study at sedimentary basin in volcanic region Malang and its surrounding area, East Java). IOP Conf. Ser. Earth Environ. Sci. 2021, 873, 012008. [Google Scholar] [CrossRef]
Mickus, K.L.; Hinojosa, J.H. The complete gravity gradient tensor derived from the vertical component of gravity: A Fourier transform technique. J. Appl. Geophys. 2001, 46, 159–174. [Google Scholar] [CrossRef]
Zhang, C.; Mushayandebvu, M.F.; Reid, A.B.; Fairhead, J.D.; Odegard, M.E. Euler deconvolution of gravity tensor gradient data. Geophysics 2000, 65, 512–520. [Google Scholar] [CrossRef]
Reid, A.B.; Ebbing, J.; Webb, S.J. Avoidable Euler errors—The use and abuse of Euler deconvolution applied to potential fields. Geophys. Prospect. 2014, 62, 1162–1168. [Google Scholar] [CrossRef]
Pašteka, R.; Richter, F.; Karcol, R.; Brazda, K.; Hajach, M. Regularized derivatives of potential fields and their role in semi-automated interpretation methods. Geophys. Prospect. 2009, 57, 507–516. [Google Scholar] [CrossRef]
Duong, T.; Cowling, A.; Koch, I.; Wand, M.P. Feature significance for multivariate kernel density estimation. Comput. Stat. Data Anal. 2008, 52, 4225–4242. [Google Scholar] [CrossRef]

Figure 1. Linear binning counts: a trivariate datum

χ^{i}

is converted into the counts assigned to its eight nearest grid points. Following Rao’s work [50] and Chacón and Duong’s work [49], their respective counts are equal to their natural coordinate values.

Figure 1. Linear binning counts: a trivariate datum

χ^{i}

Figure 2. The 1D verification results. BSS: B-spline probability density estimation method; KS: Gaussian kernel smoothing estimation method; BSSFFT: B-spline probability density estimation method based on fast Fourier transform.

Figure 3. The 2D verification results: (a) the 2D random dataset,

n = 3000

; (b) the BSS result; (c) the BSSFFT result; (d) the true probability density function (pdf); (e) the relative error between the BSS result and true pdf; (f) the relative error between the BSSFFT result and true pdf.

Figure 3. The 2D verification results: (a) the 2D random dataset,

n = 3000

Figure 4.

g_{z}

of the cubic model (a) without noise and (b) with

3 %

Gaussian noise. The residual density of the cube is 0.36 g/cm³.

Figure 4.

g_{z}

of the cubic model (a) without noise and (b) with

3 %

Gaussian noise. The residual density of the cube is 0.36 g/cm³.

Figure 5. Scatter plots of Euler solutions. (a)

g_{z}

without noise and (b) with

3 %

Gaussian noise, Euler solutions filtered by

σ_{r} \leq 1 \times 10^{- 3}

[4] in different views: (c)

(- 37.5, 30)

and (d)

(90, 0)

Figure 5. Scatter plots of Euler solutions. (a)

g_{z}

without noise and (b) with

3 %

Gaussian noise, Euler solutions filtered by

σ_{r} \leq 1 \times 10^{- 3}

[4] in different views: (c)

(- 37.5, 30)

and (d)

(90, 0)

Figure 6. Probability density curves obtained using 1D BSS and BSSFFT with

λ = 100

for subsets of the Euler solutions derived from noise-corrupted data. (a) subsets

{x_{0}}

{y_{0}}

, and

{z_{0}}

; (b) subset

{N}

Figure 6. Probability density curves obtained using 1D BSS and BSSFFT with

λ = 100

for subsets of the Euler solutions derived from noise-corrupted data. (a) subsets

{x_{0}}

{y_{0}}

, and

{z_{0}}

; (b) subset

{N}

Figure 7. Probability density distribution obtained using 2D BSS and BSSFFT with

λ = [100, 100]

for subsets

{x_{0}, y_{0}}

{x_{0}, z_{0}}

{y_{0}, z_{0}}

{x_{0}, N}

{y_{0}, N}

, and

{z_{0}, N}

of the cubic model.

Figure 7. Probability density distribution obtained using 2D BSS and BSSFFT with

λ = [100, 100]

for subsets

{x_{0}, y_{0}}

{x_{0}, z_{0}}

{y_{0}, z_{0}}

{x_{0}, N}

{y_{0}, N}

, and

{z_{0}, N}

of the cubic model.

Figure 8. Probability density isosurface obtained using 3D BSS with

λ = [100, 100, 100]

for Euler solution subsets: (a)

\{x_{0}, y_{0}, z_{0}\}

; (b)

\{x_{0}, y_{0}, N\}

; (c)

\{x_{0}, z_{0}, N\}

; (d)

\{y_{0}, z_{0}, N\}

Figure 8. Probability density isosurface obtained using 3D BSS with

λ = [100, 100, 100]

for Euler solution subsets: (a)

\{x_{0}, y_{0}, z_{0}\}

; (b)

\{x_{0}, y_{0}, N\}

; (c)

\{x_{0}, z_{0}, N\}

; (d)

\{y_{0}, z_{0}, N\}

Figure 9. Probability density isosurface obtained using 3D BSSFFT with

λ = [100, 100, 100]

for subsets (a)

\{x_{0}, y_{0}, z_{0}\}

; (b)

\{x_{0}, y_{0}, N\}

; (c)

\{x_{0}, z_{0}, N\}

; (d)

\{y_{0}, z_{0}, N\}

Figure 9. Probability density isosurface obtained using 3D BSSFFT with

λ = [100, 100, 100]

for subsets (a)

\{x_{0}, y_{0}, z_{0}\}

; (b)

\{x_{0}, y_{0}, N\}

; (c)

\{x_{0}, z_{0}, N\}

; (d)

\{y_{0}, z_{0}, N\}

Figure 10.

g_{z}

of the synthetic model (a) without noise and (b) with

3 %

Gaussian noise. The density values of left and right anomalous sources are 0.3 g/cm³ and 0.7 g/cm³, respectively.

Figure 10.

g_{z}

of the synthetic model (a) without noise and (b) with

3 %

Gaussian noise. The density values of left and right anomalous sources are 0.3 g/cm³ and 0.7 g/cm³, respectively.

Figure 11. Scatter plots of Euler solutions. The red box and black box are representative of the abnormal source. (a)

g_{z}

without noise and (b) with

3 %

Gaussian noise, Euler solutions filtered by

σ_{r} \leq 1 \times 10^{- 3}

in different views: (c)

(- 37.5, 30)

and (d)

(90, 0)

Figure 11. Scatter plots of Euler solutions. The red box and black box are representative of the abnormal source. (a)

g_{z}

without noise and (b) with

3 %

Gaussian noise, Euler solutions filtered by

σ_{r} \leq 1 \times 10^{- 3}

in different views: (c)

(- 37.5, 30)

and (d)

(90, 0)

Figure 12. Probability density curves obtained using 1D BSS with

λ = 100

for subsets of the combination models. (a) subsets

{x_{0}}

{y_{0}}

, and

{z_{0}}

; (b) subset

{N}

Figure 12. Probability density curves obtained using 1D BSS with

λ = 100

for subsets of the combination models. (a) subsets

{x_{0}}

{y_{0}}

, and

{z_{0}}

; (b) subset

{N}

Figure 13. Probability density distribution obtained using 2D BSS and BSSFFT with

λ = [100, 100]

for subsets

{y_{0}, x_{0}}

{z_{0}, x_{0}}

{N, x_{0}}

{z_{0}, y_{0}}

{N, y_{0}}

, and

{N, z_{0}}

of the combination model.

Figure 13. Probability density distribution obtained using 2D BSS and BSSFFT with

λ = [100, 100]

for subsets

{y_{0}, x_{0}}

{z_{0}, x_{0}}

{N, x_{0}}

{z_{0}, y_{0}}

{N, y_{0}}

, and

{N, z_{0}}

of the combination model.

Figure 14. Probability density isosurface obtained using 3D BSSFFT with

λ = [100, 100, 100]

for subsets (a)

\{x_{0}, y_{0}, z_{0}\}

, (b)

\{x_{0}, y_{0}, N\}

, (c)

\{x_{0}, z_{0}, N\}

, and (d)

\{y_{0}, z_{0}, N\}

of the combination model. The red box and black box are representative of the abnormal source.

Figure 14. Probability density isosurface obtained using 3D BSSFFT with

λ = [100, 100, 100]

for subsets (a)

\{x_{0}, y_{0}, z_{0}\}

, (b)

\{x_{0}, y_{0}, N\}

, (c)

\{x_{0}, z_{0}, N\}

, and (d)

\{y_{0}, z_{0}, N\}

of the combination model. The red box and black box are representative of the abnormal source.

Figure 15. Clustering methods used to separate the adjacent clusters formed by Euler solutions. The red box and blue box are representative of the abnormal source. (a–c) Contaminated gravity; (d–f) Euler solutions obtained using

w_{x} \equiv w_{y} = 3 \sim 15

; (g–i) Euler solutions after rejection by

σ_{r} \leq 1 \times 10^{- 3}

; (j–l) clusters obtained via the K-means clustering algorithm (K-means), where the number of clusters is predetermined by 2; and (m–o) three clusters automatically obtained using the density-based spatial clustering of applications with noise (DBSCAN) algorithm, setting the number of targets in their neighborhood to 1% of the total number of samples. There are 32,362, 42,237, and 48,266 solutions after filtering in (g–i), respectively.

w_{x} \equiv w_{y} = 3 \sim 15

; (g–i) Euler solutions after rejection by

σ_{r} \leq 1 \times 10^{- 3}

Figure 16. Probability density isosurfaces obtained using 3D BSSFFT for subset

x_{0}, y_{0}, z_{0}

with varied separations L: (a) 4.0 (km); (b) 2.5 (km); and (c) 1.0 (km);

λ = [100, 100, 100]

. The red box and blue box are representative of the abnormal source.

Figure 16. Probability density isosurfaces obtained using 3D BSSFFT for subset

x_{0}, y_{0}, z_{0}

with varied separations L: (a) 4.0 (km); (b) 2.5 (km); and (c) 1.0 (km);

λ = [100, 100, 100]

. The red box and blue box are representative of the abnormal source.

Figure 17. Illustration of the 3D BSSFFT method’s sensitivity to different levels of Gaussian noise. The red box and blue box are representative of the abnormal source. The top, middle, and bottom rows correspond to

\bar{p} =

0%, 4%, and 8%, respectively; the columns from left to right correspond to noise-corrupted gravity, Euler solutions, Euler solutions filtered by

σ_{r} \leq 1 \times 10^{- 3}

, and the probability density distributions obtained from the BSSFFT with

λ = [100, 100, 100]

, respectively.

w_{x} \equiv w_{y} = 3 \sim 15

. There are 43,258, 42,552, and 48,923 solutions after filtering in (c,g,k), respectively.

\bar{p} =

0%, 4%, and 8%, respectively; the columns from left to right correspond to noise-corrupted gravity, Euler solutions, Euler solutions filtered by

σ_{r} \leq 1 \times 10^{- 3}

, and the probability density distributions obtained from the BSSFFT with

λ = [100, 100, 100]

, respectively.

w_{x} \equiv w_{y} = 3 \sim 15

. There are 43,258, 42,552, and 48,923 solutions after filtering in (c,g,k), respectively.

Figure 18. The profiles

g_{z}

of single and two anomalous sources at different noise levels (

p = 0, 4

), respectively.

Figure 18. The profiles

g_{z}

of single and two anomalous sources at different noise levels (

p = 0, 4

), respectively.

Figure 19. Illustration of the 2D Euler deconvolution sensitivity to different levels of Gaussian noise. The red box and green box are representative of the abnormal source. The left (a,c,e,g) and right (b,d,f,h) columns correspond to

\bar{p} =

0% and 4%, respectively; each row from top to bottom corresponds to the single cube and two cubes with L = 4.0 (km), 2.5 (km), and 1.0 (km), respectively. The density of left and right anomalous sources is 2.36 g/cm³ and 1.27 g/cm³, respectively.

w_{x} = 10 \sim 35

\bar{p} =

w_{x} = 10 \sim 35

Figure 20. Probability density distribution obtained using 2D BSSFFT with

λ = [100, 100]

for subsets

{x_{0}, z_{0}}

of the single cube and two cubes with L = 4.0 (km), 2.5 (km), and 1.0 (km), respectively. The red box and green box are representative of the abnormal source. The left (a,c,e,g) and right (b,d,f,h) columns correspond to

\bar{p} =

0% and 4%, respectively.

Figure 20. Probability density distribution obtained using 2D BSSFFT with

λ = [100, 100]

for subsets

{x_{0}, z_{0}}

\bar{p} =

0% and 4%, respectively.

Figure 21. The Bishop 5X dataset, (a) topographic map of Basement, (b) magnetic susceptibility, (c) total field magnetic anomaly,

i = 45^{\circ}

, and (d) total field magnetic anomaly,

i = 90^{\circ}

Figure 21. The Bishop 5X dataset, (a) topographic map of Basement, (b) magnetic susceptibility, (c) total field magnetic anomaly,

i = 45^{\circ}

, and (d) total field magnetic anomaly,

i = 90^{\circ}

Figure 22. The components (a)

T_{z}

and first-order derivatives (b)

T_{z x}

, (c)

T_{z y}

, and (d)

T_{z z}

of the magnetic anomaly (magnetic inclination

45^{\circ}

) of Bishop 5X data. The red line is the survey line applying the 2D Euler deconvolution.

Figure 22. The components (a)

T_{z}

and first-order derivatives (b)

T_{z x}

, (c)

T_{z y}

, and (d)

T_{z z}

of the magnetic anomaly (magnetic inclination

45^{\circ}

) of Bishop 5X data. The red line is the survey line applying the 2D Euler deconvolution.

Figure 23. The components (a)

T_{z}

and first-order derivatives (b)

T_{z x}

, (c)

T_{z y}

, and (d)

T_{z z}

of the magnetic anomaly (magnetic inclination

90^{\circ}

) of Bishop 5X data. The red line is the survey line applying the 2D Euler deconvolution.

Figure 23. The components (a)

T_{z}

and first-order derivatives (b)

T_{z x}

, (c)

T_{z y}

, and (d)

T_{z z}

of the magnetic anomaly (magnetic inclination

90^{\circ}

) of Bishop 5X data. The red line is the survey line applying the 2D Euler deconvolution.

Figure 24. Illustration of the 2D BSSFFT method for Bishop 5X profile data (

i = 45^{\circ}

). (a)

T_{z}

, (b) Euler solutions, (c) BSSFFT’s result. The computation time of the 2D Euler deconvolution and BSSFFT is 12.19 and 0.0253 s, respectively.

n = 7220, λ = [100, 100]

Figure 24. Illustration of the 2D BSSFFT method for Bishop 5X profile data (

i = 45^{\circ}

). (a)

T_{z}

, (b) Euler solutions, (c) BSSFFT’s result. The computation time of the 2D Euler deconvolution and BSSFFT is 12.19 and 0.0253 s, respectively.

n = 7220, λ = [100, 100]

Figure 25. Illustration of the 2D BSSFFT method for Bishop 5X profile data (

i = 90^{\circ}

). (a)

T_{z}

, (b) Euler solutions, (c) BSSFFT’s result. The computation times of the 2D Euler deconvolution and BSSFFT are 11.25 and 0.0372 s, respectively.

n = 7019, λ = [100, 100]

Figure 25. Illustration of the 2D BSSFFT method for Bishop 5X profile data (

i = 90^{\circ}

). (a)

T_{z}

, (b) Euler solutions, (c) BSSFFT’s result. The computation times of the 2D Euler deconvolution and BSSFFT are 11.25 and 0.0372 s, respectively.

n = 7019, λ = [100, 100]

Figure 26. The Euler solutions of Bishop 5X (magnetic inclination

45^{\circ}

Figure 26. The Euler solutions of Bishop 5X (magnetic inclination

45^{\circ}

Figure 27. The Euler solutions of Bishop 5X (magnetic inclination

90^{\circ}

Figure 27. The Euler solutions of Bishop 5X (magnetic inclination

90^{\circ}

Figure 28. Probability density slices of the Euler solution at various depths: (a) 1500, (b) 3000, (c) 4500, (d) 6000, (e) 7500, (f) 9000, (g) 10,500, (h) 12,000, and (i) 13,500 (magnetic inclination

i = 45^{\circ}

). A–D are denoted as probability density peaks.

i = 45^{\circ}

). A–D are denoted as probability density peaks.

Figure 29. Probability density slices of the Euler solution corresponding to (a) 1500, (b) 3000, (c) 4500, (d) 6000, (e) 7500, (f) 9000, (g) 10,500, (h) 12,000, and (i) 13,500 (magnetic inclination

i = 90^{\circ}

). A–D are denoted as probability density peaks.

i = 90^{\circ}

). A–D are denoted as probability density peaks.

Figure 30. Probability density isosurfaces derived from Bishop 5X magnetic data (

i = 45^{\circ}

) in (a) the perspective view and (b) the plane view (from bottom to top). A–D are denoted as probability density peaks.

Figure 30. Probability density isosurfaces derived from Bishop 5X magnetic data (

i = 45^{\circ}

) in (a) the perspective view and (b) the plane view (from bottom to top). A–D are denoted as probability density peaks.

Figure 31. Probability density isosurfaces derived from Bishop 5X magnetic data (

i = 90^{\circ}

) in (a) the perspective view and (b) the plane view (from bottom to top). A–D are denoted as probability density peaks.

Figure 31. Probability density isosurfaces derived from Bishop 5X magnetic data (

i = 90^{\circ}

) in (a) the perspective view and (b) the plane view (from bottom to top). A–D are denoted as probability density peaks.

Table 1. The 1D probability density estimation results.

	Cube 1 ( $- 1500, - 1500, 2000$ )				Cube 2 ( $1500, 1500, 2000$ )
1D Subset	$x_{0}$	$y_{0}$	$z_{0}$	$N$	$x_{0}$	$y_{0}$	$z_{0}$	$N$
Theoretical values	−1500	−1500	2000	2	1500	1500	2000	2
BSS results	−1632	−1731	1051	1.33	1634	1631	1051	1.33
BSSFFT’s results	−1731	−1731	1051	1.40	1634	1668	1052	1.40

Table 2. The 2D probability density estimation results.

2D Subset	$(x_{0}, y_{0})$	$(z_{0}, x_{0})$	$(N, x_{0})$	$(z_{0}, y_{0})$	$(N, y_{0})$	$(N, z_{0})$
BSS results for Cube 1	$(1724, - 1591)$	$(1494, 1784)$	$(1.33, 1754)$	$(1515, - 1644)$	$(1.32, - 1664)$	$(1.3, 1275)$ *
BSS results for Cube 2	$(- 1584, 1892)$	$(1041, - 1644)$	$(1.34, - 1674)$	$(1041, 1724)$	$(1.35, 1655)$	$(1.3, 1275)$ *
BSSFFT results for Cube 1	$(1605, - 1642)$	$(1598, 1605)$	$(1.58, 1575)$	$(1639, - 1703)$	$(1.58, - 1644)$	$(1.50, 1580)$
BSSFFT results for Cube 2	$(- 1614, 1690)$	$(1103, - 1733)$	$(1.62, - 1644)$	$(1103, 1665)$	$(1.60, 1605)$	$(1.58, 1129)$

* Unable to separate adjacent anomalies.

Table 3. The 3D BSSFFT results.

	$\{x_{0}, y_{0}, z_{0}\}$	$\{x_{0}, y_{0}, N\}$	$\{x_{0}, z_{0}, N\}$	$\{y_{0}, z_{0}, N\}$
Cube 1	$(1472, - 1901, 2200)$	$(1386, - 1860, 1.98)$	$(1371, 2200, 1.89)$	$(- 1874, 2280, 2.01)$
Cube 2	$(- 1860, 1404, 1600)$	$(- 1860, 1324, 2.04)$	$(- 1963, 1640, 2.01)$	$(1329, 1680, 2.10)$

Table 4. The geometric parameters of the three synthetic models with various separations.

Separation L	Centroid of Left Cube	Centroid of Right Cube
4000	$(- 4000, 4000, 3000)$	$(4000, - 4000, 2000)$
2500	$(- 2500, 2500, 3000)$	$(2500, - 2500, 2000)$
1000	$(- 1000, 1000, 3000)$	$(1000, - 1000, 2000)$

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Cao, S.; Chen, P.; Lu, G.; Ma, Z.; Yang, B.; Chen, X. FFT-Based Probability Density Imaging of Euler Solutions. Entropy 2024, 26, 517. https://doi.org/10.3390/e26060517

AMA Style

Cao S, Chen P, Lu G, Ma Z, Yang B, Chen X. FFT-Based Probability Density Imaging of Euler Solutions. Entropy. 2024; 26(6):517. https://doi.org/10.3390/e26060517

Chicago/Turabian Style

Cao, Shujin, Peng Chen, Guangyin Lu, Zhiyuan Ma, Bo Yang, and Xinyue Chen. 2024. "FFT-Based Probability Density Imaging of Euler Solutions" Entropy 26, no. 6: 517. https://doi.org/10.3390/e26060517

APA Style

Cao, S., Chen, P., Lu, G., Ma, Z., Yang, B., & Chen, X. (2024). FFT-Based Probability Density Imaging of Euler Solutions. Entropy, 26(6), 517. https://doi.org/10.3390/e26060517

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

FFT-Based Probability Density Imaging of Euler Solutions

Abstract

1. Introduction

2. Theory

2.1. Principle of Gravity Euler Deconvolution

2.2. Principle of BSS for Euler Solutions

2.3. Principle of BSSFFT

3. Algorithm Verification

4. Model Experiment

4.1. Single Cubic Model

4.2. Synthetic Model

4.3. Sensitivity of 3D BSSFFT to Separations

4.4. Sensitivity of 3D BSSFFT to Gaussian Noise

4.5. Sensitivity of 2D Euler Deconvolution to Gaussian Noise

5. Bishop Model

6. Discussion

7. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

Abbreviations

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI