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

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

Problem Statement
The paper addresses the fundamental challenges inherent in 3D gravity inversion: underdetermination, instability, and poor depth resolution. Traditional surface observations provide data only on a 2D plane ($z=0$), lacking vertical wavenumber information required to resolve the 3D subsurface density distribution $\rho(x,y,z)$. This leads to an insufficient number of independent constraints. Furthermore, the gravity field decays rapidly with depth, causing signals from deep sources to be severely attenuated, rendering the inverse problem ill-conditioned and highly sensitive to noise. Existing approaches, such as spatial-domain iterative regularization (e.g., Tikhonov), are computationally expensive and sensitive to initial models, while frequency-domain methods are typically limited to interface reconstruction rather than full volumetric inversion.

Methodology
The author proposes a "true 3D" inversion method based on the spectral-domain pseudo-inverse theory for unitary diagonalizable systems. The methodology proceeds in two distinct theoretical levels:

1. Data-Space Completion via Analytical Continuation:
   To satisfy the requirement of the spectral-domain pseudo-inverse for 3D input data, the paper exploits the fact that the gravity anomaly $g_z$ satisfies Laplace's equation ($\nabla^2 g_z = 0$) in the source-free air layer ($z < 0$). By applying a 2D Fourier transform to the surface data $g_z(x,y,0)$, the author derives an analytical upward continuation formula. This allows the construction of a complete, physically self-consistent 3D data volume above the observation plane. Crucially, the paper derives a general analytical expression for the half-space spectrum:
   $$\tilde{g}^{obs}_z(k_x, k_y, k_z) = \tilde{g}_z(k_x, k_y, 0) \frac{1}{a - i k_z}$$
   where $a = \sqrt{k_x^2 + k_y^2}$. This expression separates the horizontal spectrum from a universal vertical propagation kernel, providing the necessary vertical wavenumber sampling for the 3D Fourier transform without numerical interpolation errors.

2. Spectral-Domain Pseudo-Inverse Inversion:
   The forward operator for 3D gravity is shown to be unitarily diagonalizable in the Fourier domain. The relationship between the density spectrum $\tilde{\rho}$ and the gravity spectrum $\tilde{g}_z$ is expressed as:
   $$\tilde{g}_z(\mathbf{k}) = \lambda(\mathbf{k}) \tilde{\rho}(\mathbf{k}), \quad \text{where } \lambda(\mathbf{k}) = \frac{-i 4\pi G k_z}{K^2}$$
   Here, $K^2 = k_x^2 + k_y^2 + k_z^2$. Based on the author's prior theory for unitary diagonalizable systems, a stable generalized inverse filter is constructed:
   $$q_\alpha(\lambda) = \frac{\bar{\lambda}}{|\lambda|^2 + \alpha}$$
   where $\alpha > 0$ is a regularization parameter. This filter is applied to the constructed 3D data spectrum to obtain the density spectrum. The final spatial density distribution is recovered via an inverse 3D Fourier transform. Due to the pole structure of the data spectrum (specifically the pole at $k_z = -ia$), the $k_z$ inverse transform is implemented via a contour integral in the complex plane, closed in the lower half-plane for $z > 0$.

Key Contributions

• Data-Space Completion: The paper demonstrates that upward continuation via Laplace's equation can rigorously construct a 3D data volume from 2D surface data, thereby mitigating the underdetermination of the problem and providing the vertical wavenumber sampling required for 3D spectral analysis.
• Closed-Form Inversion: The method proves that the 3D gravity forward operator is exactly unitarily diagonalizable. Consequently, the inversion is reduced to a deterministic closed-form calculation consisting of one forward 3D Fourier transform, a spectral scaling step using the derived filter, and one inverse transform, eliminating the need for iterative optimization.
• Stability and Consistency: The constructed inverse filter is proven to have bounded stability ($\|q_\alpha\|_\infty \leq 1/(2\sqrt{\alpha})$) and consistency (converging to the Moore–Penrose generalized inverse as $\alpha \to 0^+$).
• Analytical Framework: The derivation of the half-space spectrum (Equation 7) provides a unified input framework that separates horizontal and vertical components, applicable to arbitrary surface gravity data regardless of the subsurface model complexity.

Results
The method was validated using a homogeneous sphere model. The results indicate that:

• The inversion correctly recovers the 3D location of the anomaly.
• As the regularization parameter $\alpha \to 0^+$, the method successfully removes the Green's function spreading effect, recovering the singular behavior (Dirac delta function) at the source point.
• For finite $\alpha$, the solution is a smoothed distribution, and the integral over the entire space strictly conserves the total mass of the source.
• The method recovers the perturbation density (relative to an unknown background) with a zero mean; absolute densities require independent prior information. Quantitative recovery of the density distribution depends on the choice of $\alpha$: smaller $\alpha$ yields sharper results and more accurate internal values, approaching the divergent limit at the source.

Significance and Claims
The paper claims to establish a theoretically self-consistent, analytical framework for "true 3D" gravity inversion. Its significance lies in shifting the paradigm from iterative optimization to a structural decomposition approach. By proving that the problem can be solved via a single spectral filtering operation, the method offers a unified framework for stable inversion of convolution-type problems. The author notes that while the method is numerically equivalent to zeroth-order Tikhonov regularization, it differs conceptually by deriving directly from the unitary diagonalizable structure of the Poisson equation rather than an optimization principle. The paper concludes that while the method achieves ideal results under continuous-spectrum assumptions, future work will focus on applications to real data, numerical implementation details, and robustness analysis.