Inverting for arbitrarily shaped targets using parameterized radial basis functions

Introduction¶

Inversion of geophysical data is a crucial step in mineral exploration for identifying potential deposits. Mineral deposits are rarely simple; they often exhibit curvy, irregular, and non-smooth geometries. Accurately representing these complex shapes is particularly important prior to drilling, but remains a challenge for conventional inversion techniques.

Traditional least squares regularized inversion techniques struggle with shapes that have sharp edges or small physical property contrasts. An alternative is parametric shape-based inversion, where shapes such as rectangular prisms Belliveau & Haber, 2023 or ellipsoids McMillan et al., 2015 are deformed using a set of associated parameters to recover the target. Parametric methods offer the significant advantage of reducing the total parameter space of the inverse problem. These methods, however, have notable drawbacks: (1) they require prior knowledge and initialization of the expected number of targets; (2) they perform poorly if the target deviates from the assumed parameterized shape.

Methods based on parametric basis functions have been shown to be able to address these issues Aghasi et al., 2011Kadu et al., 2017Ozsar et al., 2025. Instead of selecting a predefined parameterized shape, approaches based on Radial Basis Functions (RBFs) offer the benefits of parametric inversion without imposing strict shape constraints.

In this work, we adapt the RBF method originally described by Aghasi et al., 2011. We parameterize the model space with a user-defined number of RBFs, each associated with scalar weights that serve as the inversion parameters. The RBFs form a level-set function that is binarized using a chosen indicator function, which defines the width of the boundary between the shape and background model helping to delineate the target shape. The level-set is updated during the inversion to find a model that fits the data. We validate this method using two models. The first features two distinct targets: a rectangular block and a sphere, while the second contains an anomaly with an arbitrary shape. Both are imaged using a 2D DC Resistivity (DCR) survey. The anomaly is treated as a level-set with fixed conductivity. While we demonstrate the method in 2D DCR, it is not limited to this setting. The approach is both method- and dimension-agnostic and has been successfully applied to 3D problems and other geophysical data types, including potential field and seismic data Kadu et al., 2017.

Level-Set Approach to Inversion¶

In geophysical inverse problems, a common goal is to determine the location, shape, and physical properties of a target of interest in a heterogeneous background. Consider a general inverse problem where the forward model is:

with $\mathbf{m}$ representing the discretized model, $\mathbf{d}$ the observed data, and $F$ the forward operator that maps the model to the data space. Here, $\eta$ is the noise in the data, usually assumed to be Gaussian.

For level-set methods, we define the shape and location of the target as the $c$ level-set of some function ${\phi(\mathbf{x})} = c$ where $\mathbf{x}$ denotes the spatial coordinates in model space and $c$ is a scalar (usually taken to be zero). In cases where one or more anomalous bodies in present in a background, the overall model space can then be discretized as

where $\mathbf{m_0(x)}$ is the background model and $\mathbf{m_p(x)}$ describes the physical properties of the anomalous bodies; in the simplest case, $\mathbf{m_0(x)}$ and and $\mathbf{m_p(x)}$ are scalars. Here, $\mathbf{m_p}$ is multiplied by a level-set function $\phi(\mathbf{x},\mathbf{p})$, where $\mathbf{p}$ is a vector of parameters determined during inversion. The level-set is refined by an indicator function, $H$, which defines the boundary of the level-set domain $\Omega$ (and hence the boundary of the target). A general form of the indicator function is a step function:

For parametric level-set methods, the level-set function $\phi(\mathbf{x})$ is represented in terms of a set of parameters $\mathbf{p}$. In geometric parameterizations (e.g., ellipses or prisms), $\mathbf{p}$ might include the location of the shape’s centre, strike and dip extents, rotation angles, etc. Belliveau & Haber, 2023. In contrast, this work eschews geometric parameterization of the level-set function in favour of Radial Basis Functions (RBFs), given that any smooth level-set function can be approximated by a linear combination of a set of sufficiently smooth RBFs Aghasi et al., 2011.

We first define a set of RBFs, $\psi(\mathbf{r})$, centred on a grid with a spacing ($\mathbf{x_{rbf}}$) that is sparser than the simulation grid (denoted by $\mathbf{x}$). Here, $\mathbf{r} = \mathbf{x} - \mathbf{x_{rbf}}$. Then, we represent the level-set function as a weighted sum of these RBFs:

where $\alpha_j$ is the scaling parameter for the $j^{th}$ RBF, and $\beta$ is the dilation parameter that determines its width (Figure 1).

There are various choices of RBF, such as Gaussian, thin-plate splines, multiquadratic, compactly-supported RBFs, etc. In shape-reconstruction inverse problems, Gaussian and Wendland functions are commonly used. In this study, we employ the $C^6$ Wendland function as it is shown to improve the sparsity of the system matrices, reduce the computational cost, and recover sharper edgesAghasi et al., 2011. The $C^6$ Wendland function is defined as

where $(1-r)_+$ denotes that the term is non-zero only when $r<1$.

In general, one can choose to parameterize the level-set function with $\mathbf{p} = [\boldsymbol{\alpha}, \beta, \mathbf{x}_{rbf}]$. For even greater flexibility, further parameters (such as stretching or skewing of each basis) can be introduced Ozsar et al., 2025. However, in this study, we only treat $\boldsymbol{\alpha}$ as the inversion parameters, keeping $\beta$ and $\mathbf{x}_{rbf}$ fixed. This choice reduces the dimensionality and improves numerical stability, since multiple pairs of $[\alpha, \beta]$ can yield the same $\phi(\mathbf{x}, \mathbf{p})$. A positive value of $\alpha_j$ produces a positive “bump” in the level-set function, with higher values increasing the sharpness—thus aiding in resolving sharp corners or edges of a target.

We fix the dilation parameter as $\beta = 1/\Delta_{rbf}$, where $\Delta_{rbf}$ is the spacing between RBF centers (typically between 2 and 4 times the simulation mesh spacing). The RBF centers $\mathbf{x}_{rbf}$ are defined on an equally spaced grid and remain fixed during the inversion.

Another choice is the selection of the Heaviside function. Since an ideal Heaviside function (a step function) is non-differentiable, a smooth sigmoid approximation is preferred. A common choice is a sigmoid function. We use a regularized ${C^2}$ Heaviside function following Aghasi et al. (2011):

Figure 2 illustrates this Heaviside function. Here, $\epsilon$ controls the steepness of the transition; choosing an appropriate $\epsilon$ is crucial, as an overly steep or flat transition can lead to numerical instabilities (e.g., vanishing gradients).

We update $\epsilon$ at each iteration and define it as a fraction of the dynamic range of the level-set function Kadu et al., 2017:

We include $\gamma$ as an additional inversion parameter.

Thus, the inverse problem is formulated as a minimization of a misfit function with respect to the inversion parameters $\mathbf{p} = [\boldsymbol{\alpha}, \gamma]$:

where $\lambda$ is a weighting factor for the regularization term $R$ that can be applied to the parameters $\mathbf{p}$.

We then invert the data by solving the minimization problem in equation (8) using a Gauss-Newton optimization method. The necessary derivatives are computed automatically via the autograd functionality in PyTorch Ansel et al., 2024.

Examples¶

We test our inversion approach using a synthetic 2D DC resistivity survey. Two model scenarios are examined: one containing compact targets (a block and a circular target) in a halfspace, and one where we aim to estimate depth to a basement layer. The data are generated for a dipole-dipole configuration with electrodes spaced every 50 m along $x = [-1000, 1000]$ on a flat topography. All simulations and inversions are performed using the SimPEG package Heagy et al., 2017Cockett et al., 2015.

Example 1: Multiple Compact Targets¶

We consider the model shown in Figure 3a consisting of a 400 m $\times$ 200 m rectangular block centred at [250, -200], and a circular target with a 150 m radius centred at [-600, -200] in a homogeneous halfspace. This model is chosen to evaluate the method’s ability to resolve multiple targets with different geometries. The simulation mesh uses a 40 $\times$ 40 cell grid with 50 m spacing, resulting in 800 active cells.

First, we show an inverted model that is recovered using a regularized weighted least-squares method, starting from a halfspace model (Figure 3b). The inverted model is able to image the targets in approximately the correct locations but the boundaries of the shapes are smooth and difficult to interpret. For the parametric RBF inversion, we initialize the model (Figure 3c) with 187 RBFs distributed uniformly every 150 m within the active region of the mesh. The $\mathbf{\alpha}$ parameters are drawn from a normal distribution with a standard deviation of 0.1. To reduce the influence of RBFs located below the survey’s sensitivity (below $z = -500$ m), these are assigned a large negative weight ($\alpha = -10$) to dampen their impact on the inversion. No explicit regularization is enforced (i.e., $\lambda=0$ in equation (8)) in the parametric RBF inversion, making this step necessary to prevent overfitting. We initialize the $\gamma$ parameter in equation 7 as 0.1, corresponding to 10% of the dynamic range of the initial level-set function. The conductivity values for the background and targets, $m_0$ and $m_p$ respectively, are set to match those of the true model. Absent a-priori information, the physical properties can also be included as additional parameters in the inversion.

Our inverted model (Figure 3d) successfully recovers the geometry of both targets. The inversion forms a parametric level-set function whose zero-level set delineates two distinct, compact targets. This result highlights the advantage of our approach: it can recover multiple shapes without prior information, unlike other parametric methods that require the number of targets to be specified at initialization.

Example 2: Depth to Basement¶

To further evaluate our method, we test it on a model that simulates depth-to-basement variations. The true model (Figure 4a) represents a two-layer Earth with a higher-conductivity lower layer (the basement) whose depth varies sinusoidally. At its shallowest the basement is at $z=-200$ m and further along the line the basement deepens up to $z=-500$ m. The simulation mesh uses a 40 $\times$ 40 cell grid with 50 m spacing, resulting in 800 active cells.

For the RBF inversion (Figure 4b), the RBFs are arranged similarly to the previous example but are spaced every 200 m, resulting in 128 parameters to be estimated in the inversion. The $\alpha$ parameters are again initialized from a normal distribution with a standard deviation of 0.1. Here, RBFs located above $z = -600$ m are assigned a large negative weight ($\alpha = -10$) to introduce a conductive layer in the initial model at a depth where the survey starts losing sensitivity. This constraint is necessary in the absence of explicit regularization.

The resulting inverted model (Figure 4c) effectively recovers the undulating structure of the basement layer, despite the limited sensitivity of the survey at these depths.

Conclusion¶

In this work, we have demonstrated that the RBF level-set approach offers a flexible framework for geophysical inversion. Our results highlight the method’s ability to recover arbitrarily shaped features, as illustrated by both the compact targets in a halfspace and the undulating basement in a two-layer model. Importantly, the approach does not require a priori knowledge of the number or geometry of targets, setting it apart from traditional parametric methods.

While our examples focused on 2D DC resistivity surveys, the method is inherently generalizable and can be applied to 3D problems as well as other geophysical techniques, including magnetics, gravity, electromagnetic, and seismic methods.

Our future work will focus on combining this method with hybrid inversion techniques where a weighted least squares inversion can be used to estimate any heterogeneity in the background physical properties Belliveau & Haber, 2023. This information can be incorporated into the model parameterization ($\mathbf{m_0}$ term in equation 2) or as a reference model for regularization. Another focus will be designing appropriate regularization in the RBF parametric model space. Early efforts show that adding regularization is helpful in introducing a priori knowledge, particularly outside the core sensitivity region of the surveys.

The inversion process still requires careful selection of bounds and an appropriate initial model to ensure numerical stability. In our experience, the RBF approach is more robust to the initial model as compared to geometric parameterization approaches. Additionally, resolving multiple bodies with different physical properties remains challenging without specialized modifications. Our ongoing research focuses on addressing this limitation by incorporating anisotropic RBFs with stretching and skewing Ozsar et al., 2025 and associating a physical property parameter with each RBF to recover multi-contrast targets within a single level-set function.