Symmetry-Based Reformulation of the Gravity Sensitivity Matrix for Efficient Regional-Scale Gravity Imaging

Regional-scale gravity imaging is limited by the size of the sensitivity matrix: a direct three-dimensional formulation of the gravity operator requires on the order of 1013 coefficients, making explicit storage infeasible on standard hardware. We exploit two symmetry properties of the gravitational sensitivity kernel: translational invariance, giving the operator a Block Toeplitz with Toeplitz Blocks (BTTB) structure, and central symmetry, halving the number of unique kernel evaluations, to reformulate the Linear Back-Projection (LBP) adjoint operator as depth-wise two-dimensional convolutions. With finite-support truncation at cutoff radius Rmax, memory requirements drop from terabyte-scale dense arrays to compact per-layer kernels. We implement and benchmark two equivalent backends: an OpenMP spatial-domain stencil and an FFT convolution via BTTB circulant embedding. For a production-scale model (701×401×512 cells, Rmax=150), the stencil completed the LBP computation in 7187 s on one thread and 354 s on 32 threads (S32=20.29), while the FFT backend required 13.8 s; both agree to within an invariant checksum (3.3464×1010). Applied to Bouguer gravity data over the Chicxulub impact crater, northern Yucatán Peninsula, the raw LBP amplitudes were calibrated to physical density contrasts through a two-parameter rescaling solved in closed form as a box-constrained convex quadratic program and cross-validated against an independent simulated-annealing search, reducing the RMS misfit from 27.3 to 5.3 mGal (an 81% reduction). The resulting model shows density patterns consistent with the crater’s known structural organization. As a single-pass adjoint method, LBP does not resolve the non-uniqueness of the gravity inverse problem, and the model is interpreted for regional structure rather than for pointwise density values; the formulation extends directly to iterative regularized inversion.