A Spectral-Domain Pseudo-Inverse Method for True 3D Gravity Inversion

1 Introduction

• Surface gravity data must be used to recover subsurface density in 3D, despite underdetermination, instability, and poor depth resolution from measurements on one plane.The gravity field decays with depth, attenuating signals from deeper sources at the surface.Source: §1 · pp. 1–2
• Existing spatial-domain iterative methods can be computationally expensive and sensitive to initial models and regularization parameters, while frequency-domain interface methods do not directly extend to volumetric density inversion.The introduction cites Tikhonov regularization, conjugate gradient techniques, and Parker–Oldenburg inversion.Source: §1 · p. 2
• The prior pseudo-inverse theory constructs a stable generalized inverse in the transform domain for unitary-diagonalizable forward operators, with regularization balancing resolution and stability.The construction has bounded stability and converges to the Moore–Penrose generalized inverse as α approaches zero from above.Source: §1 · p. 2
• This work extends that theory to 3D gravity inversion, where the Poisson-equation Green’s function must be unitarily diagonalizable and the input data must form a complete 3D volume.Surface observations alone do not obviously meet the volume-data requirement.Source: §1 · p. 2
• The framework uses upward continuation to construct a 3D data volume, then expresses the forward operator spectrally so inversion requires a forward 3D Fourier transform, spectral filtering, and an inverse transform.The paper presents this as a deterministic closed-form calculation and a unified framework for gravity and magnetic potential-field inversion, with different spectral responses.Source: §1 · pp. 2–3

2.1 Motivation

• The source-free upper half-space lets Laplace’s equation extend surface gravity data into a 3D volume, supplying input on a grid for the pseudo-inverse method.Given surface Dirichlet data, the continuation is unique and stable.Source: §2.1 · p. 3

2.2 Frequency-Domain Analytical Continuation Formula

• In the frequency domain, each Fourier component decays independently during upward continuation, with higher wavenumbers attenuating faster than lower ones.Source: §2.2 · p. 3

2.3 Mathematical Description of Data-Space Completion

• Continued levels and the measured surface form a 3D data volume that matches the boundary data and is uniquely determined without adjustable parameters.The upper-half-space spectrum is treated continuously; the inversion grid's vertical extent and spacing are set by target depth range and desired resolution.Source: §2.3 · p. 4

2.4 Methodological Significance

• Upward continuation both extends the data domain to enable a 3D Fourier transform and enriches spatial variation in a field satisfying the Poisson equation.Source: §2.4 · p. 4

2.5 Spectral-Domain Analytical Expression for Half-Space Data

• The half-space spectrum factors into the surface field’s horizontal spectrum and a universal vertical propagation kernel, whose pole informs the inversion contour.The analytical expression avoids numerical 3D-transform sampling, truncation, and leakage errors and can be used directly in the inversion formula.Source: §2.5 · p. 5

3.1 Unitary Diagonalization of the 3D Poisson Equation

• The 3D discrete Fourier transform unitarily diagonalizes the gravity forward operator, matching the structure required by the spectral-domain pseudo-inverse theory.The proposition identifies its spectral response as the one given in equation (12).Source: §3.1 · p. 5

3.2 Construction and Application of the Spectral-Domain Pseudo-Inverse Filter

• Each density Fourier mode is obtained independently by scalar-filtering the corresponding gravity-anomaly mode.This is the decoupling produced by Fourier diagonalization.Source: §3.2 · p. 6

3.3 Implementation of the 3D Inverse Transform

• For subsurface recovery (z > 0), the inverse kz transform closes the contour in the lower half-plane around its poles, and numerical integration on a pole-free contour is described as stable.The density distribution for z ≤ 0 is not considered in the inversion result.Source: §3.3 · pp. 6–7

3.4 Stability Analysis

• The filter has bounded magnitude at most 1/(2√α), preventing unbounded noise amplification, and tends to the exact inverse on nonzero spectral components as α → 0+.The zero component remains zero.Source: §3.4 · p. 8

3.5 Treatment of the Zero-Frequency Singularity

• At k = 0, the density spectrum is undefined (0/0), reflecting that the absolute average density produces no observable gravity anomaly.For kz = 0 with K² ≠ 0, the result vanishes, consistent with horizontal infinite slabs generating no gravity anomaly.Source: §3.5 · p. 8

4.1 Test Model

• The test model is a homogeneous sphere with center depth zc = 5000 m, radius R = 1000 m, and density anomaly ρ0 = 500 kg/m³.Its total mass is M = ρ0 · 4πR³/3.Source: §4.1 · p. 8

4.2 Forward Gravity Anomaly on the Surface and in the Air

• The sphere’s external field is equivalent to that of a point mass at its center, and its vertical gravity anomaly is specified for the surface and air.The total mass is used in the forward-field expression.Source: §4.1 · p. 8, §4.2 · p. 9

4.3 Inversion

• The inverse-transform derivation applies a 2D inverse Fourier transform to the intermediate result and omits the δ(z) term for z > 0.The resulting density expression is given for the subsurface.Source: §4.3 · p. 9

4.4 Analysis of Results

• Using half-space data, the inversion locates the anomaly, removes Green’s-function spreading, recovers the source singularity, and recovers total mass by integration.As α → 0, the recovered density tends to a Dirac delta at the source; total mass is recovered by integrating over all space.Source: §4.4 · p. 10
• For finite α > 0, regularization smooths the divergence into a finite-width distribution while preserving total mass; smaller α gives sharper results and more accurate internal density values.The mass conservation holds for any α > 0.Source: §4.4 · p. 10

5.1 Applicability and Limitations

• Finite surface patches cause spectral leakage and boundary artifacts, especially at low wavenumbers; zero-padding, smooth tapering, and avoiding boundaries can mitigate them.The inversion region should be restricted to a sub-area away from boundaries.Source: §5.1 · p. 10

5.2 Spatial Resolution

• The method has uniform depth resolution theoretically, but lower signal-to-noise ratios for deep sources reduce their practical resolution.Regularization suppresses weak spectral components along with high-frequency noise, preferentially affecting deeper sources.Source: §5.2 · pp. 10–11

5.3 Background Density Invisibility

• Because the DC component is unobservable, inversion recovers perturbation density rather than absolute density; absolute values require independent prior information.Examples of such information include borehole or seismic constraints.Source: §5.3 · p. 11

5.4 Relation to Classical Methods

• The method is numerically equivalent to zeroth-order Tikhonov regularization, but derives from the operator’s unitary diagonalizable structure rather than an optimization principle.Its filter resembles Wiener filtering, but α is a deterministic regularization parameter, not a statistical signal-to-noise ratio parameter.Source: §5.4 · p. 11

5.5 Importance and Generality of Equation (7)

• Equation (7) separates the half-space spectrum into a horizontal spectrum and a universal vertical propagation kernel, and applies to arbitrary surface gravity data.The derivation depends on Laplace’s equation in the air layer and data occupying a half-space, not on a particular density model.Source: §5.5 · p. 11

5.6 Pole Structure of the Half-Space Fourier Transform and Inversion Stability

• The data-domain pole kz = −ia sets the contour closure direction, while regularization poles shift with α and govern the stability–resolution trade-off.For α > 0, shifting poles away from the real axis supports stable contour integration; as α → 0, the exact inverse recovers singular behavior at the source.Source: §5.6 · p. 12

5.7 Generality of the Method

• The framework also applies to magnetic inversion and, more generally, to linear inverse problems with translation-invariant forward operators diagonalizable by the Fourier transform.For magnetic inversion, the spectral response replaces the gravity response.Source: §5.7 · p. 13

5.8 Choice of α

• The parameter α controls the resolution–stability trade-off and can be selected using L-curve or generalized cross-validation criteria.Because the filter is analytic, its effect can be assessed directly in the wavenumber domain.Source: §5.8 · p. 13

Data-space completion theory

• Upward continuation uniquely extends surface data into a physically self-consistent 3D volume, whose half-space spectrum separates horizontal data from a vertical propagation kernel.This provides the input framework for the subsequent spectral-domain inversion.Source: §6 · p. 13

Spectral-domain pseudo-inverse inversion theory

• The 3D Poisson forward operator is unitarily diagonalized with spectral response λ(k) = −i4πGkz/K2, and its regularized inverse filter has bounded stability and converges to the Moore–Penrose inverse.The filter’s stability bound is ∥qα∥∞ ≤ 1/(2√α).Source: §6 · p. 13

Closed-form inversion

• The closed-form inversion reduces to one forward 3D Fourier transform, one spectral scaling, and one inverse transform.A contour-integral numerical scheme handles the kz inverse transform; the inversion chain is analytical in theory and stable in numerical implementation.Source: §6 · p. 13

Theoretical positioning

• The method is numerically equivalent to zeroth-order Tikhonov regularization but derives from structural decomposition rather than an optimization principle.The conclusions describe this approach as unifying other potential-field inversions.Source: §6 · p. 14

Numerical validation

• Homogeneous-sphere validation recovers the anomaly’s 3D location, removes Green’s function spreading, recovers source singularity, and reproduces total mass by integration.These results are presented as analytical validation of the method’s quantitative capability.Source: §6 · p. 14

Scope and future work

• The results are theoretically self-consistent under ideal continuous-spectrum data; real-data applications, numerical implementation, and robustness analysis remain future work.The conclusions also identify image processing of high-resolution density reconstructions as a possible route to subsurface structure images.Source: §6 · p. 14

A.1 Term-by-term transformation

• For z > 0, the contour is closed in the lower half-plane, where the data factor contributes a pole at kz = −ia and the regularization denominator contributes additional poles.The regularization poles are determined by the equation associated with the denominator.Source: §A.1 · pp. 14–15
• The kernel Kα(kx, ky, z) is obtained by summing residues at kz = −ia, −ip, and −iq.Source: §A.1 · p. 15
• For finite α, the shifted poles make the regularized kernel smooth the δ-singularity into a finite-width distribution.Source: §A.1 · p. 15

A.2 Degenerate Limit: α →0+

• The appendix evaluates the kernel in the limit α = 0 using two standard integrals.Source: §A.2 · p. 15
• In the limit, Kα→0 takes the form −1/(8πG) [a sgn(z) + 2δ(z)].Source: §A.2 · p. 15
• As α → 0+, the inversion recovers the spatial δ-singularity expected for a point source.Source: §A.2 · p. 15

• The 2D inverse Fourier transform over kx and ky is applied to (B1).Source: Appendix B · p. 16

B.1 Term-by-term transformation

• For z > 0, the second term vanishes, and linearity allows the inverse transform to be handled term by term.Source: §B.1 · p. 16

B.2 Computation of F−1{ae−zca}

• Using polar coordinates for horizontal wavenumbers and integrating over angle reduces the inverse transform to an expression involving J0.Source: §B.2 · p. 16

B.3 Use of a standard integral formula

• The derivation differentiates a standard integral with respect to zc and substitutes the resulting derivative into the transform expression.Source: §B.3 · pp. 16–17

B.4 Substitution into (B3)

• The evaluated integral is substituted into (B3) to obtain the displayed expression.Source: §B.3 · p. 17, §B.4 · p. 17

B.5 Conversion to distance-variable form

• Using distance variable zc − z gives the density a Green’s-function form, and the result agrees with direct inversion of the sphere forward formula (28).Source: §B.5 · p. 17