IGUG: A MATLAB package for $3$D inversion of gravity data using graph theory

The open source MATLAB package IGUG for 3D inversion of gravity data is presented. It is based on methodology that was introduced by Bijani et al (2015), in which a homogeneous subsurface body is modeled by an ensemble of simple point masses. The model parameters are the Cartesian coordinates of the point masses and their total mass. Associating the point masses to the vertices of a weighted full graph with weights computed by the Euclidean pairwise distances separating vertices, Kruskal’s algorithm is used to find the minimum spanning tree of the graph. Stabilization is achieved using an equidistance function that restricts the spatial distribution of the masses and favors a homogeneous subsurface structure. A regularization parameter $\lambda$ is introduced to balance the two terms of the objective function, and reasonable physically-relevant bound constraints are imposed on the model parameters. Then the bound constrained objective function is solved using a genetic algorithm. A new diagnostic approach is presented for determining a suitable choice for $\lambda$, requiring a limited number of solutions for a small set of $\lambda$ without using the L-curve as suggested by Bijani et al. Simulations for synthetic examples demonstrate the efficiency and effectiveness of the implementation of the algorithm. It is verified that the constraints on the model parameters are not restrictive. Included in the package is the script GMD.m which is used for generating synthetic data and for putting measurement data in the format required for the inversion implemented by IGUG.m The script Diagnostic_Results.m is included for analyzing and visualizing the results. The software can be used to verify the simulations and the analysis of real data that is presented here, The real data set uses gravity data from the Mobrun ore body, north east of Noranda, Quebec, Canada.

💡 Research Summary

**
The paper presents IGUG, an open‑source MATLAB toolbox for three‑dimensional (3‑D) inversion of gravity data that builds on the graph‑theoretic approach introduced by Bijani et al. (2015). In this framework the subsurface is represented by an ensemble of simple point masses. The model parameters consist of the Cartesian coordinates of each point mass and the total mass of the ensemble. By treating the point masses as vertices of a fully connected weighted graph—where the edge weights are the Euclidean distances between vertices—Kruskal’s algorithm is employed to compute the Minimum Spanning Tree (MST). An “equidistance” regularization functional Θ(p) is defined as the variance of the MST edge lengths; minimizing Θ forces the edge lengths to be as equal as possible, thereby encouraging a homogeneous spatial distribution of the point masses and yielding a skeletal representation of the subsurface body.

The inversion objective function combines a data‑misfit term Φ(q) and the equidistance regularizer, weighted by a scalar regularization parameter λ: Γ(q)=Φ(q)+λ Θ(p).
Φ(q) is a weighted least‑squares misfit between observed gravity values g_obs and the forward model predictions g_z(q). The weighting matrix W_d incorporates prior knowledge of measurement noise (σ_i) by assigning each datum a weight proportional to σ_i⁻¹. This ensures that noisy observations have reduced influence on the solution.

Because Γ(q) is highly non‑linear and may possess many local minima, the authors adopt a Genetic Algorithm (GA) for global optimization. An initial population of candidate solutions is generated randomly within physically realistic bounds for coordinates and total mass. Each individual is evaluated using Γ, and standard GA operators—selection, crossover, and mutation—are applied iteratively. Constraints on the model parameters are enforced throughout the GA to respect geological prior information. The algorithm terminates when either the misfit falls below a statistical noise threshold (Φ≤N+√(2N), where N is the number of observations) or a maximum number of generations K_max is reached. The best individual across all generations is taken as the final estimate q_opt.

A key contribution of the paper is a new diagnostic procedure for selecting an appropriate λ without resorting to the traditional L‑curve. The authors run the GA for a modest set of λ values, record the convergence curves of Φ and Θ, and perform linear regression on the Θ‑versus‑Φ relationship. The regression’s coefficient of determination (R²) and residual statistics provide an objective measure of how well a given λ balances data fidelity and regularization. This approach dramatically reduces the computational burden associated with exhaustive λ sweeps and offers a statistically grounded criterion for λ selection.

The software package comprises three main scripts:

1. GMD.m – Generates synthetic point‑mass models and corresponding gravity data for a user‑specified survey geometry and noise level.
2. IGUG.m – Implements the inversion workflow, including graph construction, MST computation, evaluation of Φ and Θ, and the GA optimizer.
3. Diagnostic_Results.m – Analyzes the λ‑diagnostic runs, produces regression plots, and visualizes the recovered point‑mass skeleton.

The authors validate IGUG on two synthetic scenarios and on real gravity data from the Mobrun ore body in northeastern Noranda, Quebec, Canada. In the synthetic tests, IGUG accurately recovers the true locations, depths, and total masses of both a single spherical body and a more complex multi‑body configuration, outperforming conventional linear inversion which tends to smear the anomalies. For the Mobrun dataset, the recovered point‑mass skeleton aligns with known geological structures and previously published mass estimates, demonstrating the method’s applicability to real‑world problems.

Overall, IGUG integrates graph‑theoretic regularization, a robust global optimizer, and an automated λ‑selection scheme into a user‑friendly MATLAB toolbox. The equidistance regularizer effectively mitigates the non‑uniqueness and noise sensitivity inherent in gravity inversion, while the GA avoids dependence on initial guesses and local minima. By providing open‑source code and detailed documentation, the authors enable researchers and practitioners in geophysics, mineral exploration, and environmental geology to apply a modern, efficient, and interpretable 3‑D gravity inversion technique.