Inversion of multiconfiguration complex EMI data with minimum gradient support regularization: A case study

TL;DR

This paper presents a novel inversion algorithm for multiconfiguration complex EMI data that incorporates a minimum gradient support regularization to enhance sparsity, improving subsurface electrical conductivity imaging.

Contribution

It introduces a sparsity-promoting regularization method with a minimum gradient support stabilizer within a Gauss-Newton inversion framework for complex EMI data.

Findings

Enhanced resolution of electrical conductivity distributions.

Effective depth estimation of the investigation area.

Validated results on synthetic and real data sets.

Abstract

Frequency-domain electromagnetic instruments allow the collection of data in different configurations, that is, varying the intercoil spacing, the frequency, and the height above the ground. Their handy size makes these tools very practical for near-surface characterization in many fields of applications, for example, precision agriculture, pollution assessments, and shallow geological investigations. To this end, the inversion of either the real (in-phase) or the imaginary (quadrature) component of the signal has already been studied. Furthermore, in many situations, a regularization scheme retrieving smooth solutions is blindly applied, without taking into account the prior available knowledge. The present work discusses an algorithm for the inversion of the complex signal in its entirety, as well as a regularization method that promotes the sparsity of the reconstructed electrical conductivity distribution. This regularization strategy incorporates a minimum gradient support stabilizer into a truncated generalized singular value decomposition scheme. The results of the implementation of this sparsity-enhancing regularization at each step of a damped Gauss-Newton inversion algorithm (based on a nonlinear forward model) are compared with the solutions obtained via a standard smooth stabilizer. An approach for estimating the depth of investigation, that is, the maximum depth that can be investigated by a chosen instrument configuration in a particular experimental setting is also discussed. The effectiveness and limitations of the whole inversion algorithm are demonstrated on synthetic and real data sets.

Equations

$$Y_k(\lambda )=N_k(\lambda )\frac{Y_{k+1}(\lambda )+N_k(\lambda )\tanh (d_k{u}_k(\lambda ))}{N_k(\lambda )+Y_{k+1}(\lambda )\tanh (d_k{u}_k(\lambda ))},$$

$$M_{\nu }(\mathbf{\sigma },\mathbf{\mu };h,\omega ,\rho )=-{\rho }^{3-\nu }{\int }_0^{\infty }{\lambda }^{2-\nu }{e}^{-2h\lambda }{R}_{\omega ,0}(\lambda )J_{\nu }(\rho \lambda )\mathrm{d}\lambda ,$$

$$R_{\omega ,0}(\lambda )=\frac{N_0(\lambda )-Y_1(\lambda )}{N_0(\lambda )+Y_1(\lambda )},$$

$${\mathcal{H}}_{\nu }[f](\rho )={\int }_0^{\infty }f(\lambda )J_{\nu }(\rho \lambda )\lambda \mathrm{d}\lambda ,$$

$$M_{\nu }(\mathbf{\sigma },\mathbf{\mu };h,\omega ,\rho )=-{\rho }^{3-\nu }{\mathcal{H}}_{\nu }[{\lambda }^{1-\nu }{e}^{-2h\lambda }{R}_{\omega ,0}(\lambda )](\rho ),\nu =0,1.$$

$${\mathbf{\sigma }}^{\ast }=\arg \underset{\mathbf{\sigma }\in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert \mathbf{r}(\mathbf{\sigma }){\Vert }^2,$$

$$\mathbf{r}({\mathbf{\sigma }}_{k+1})\simeq \mathbf{r}({\mathbf{\sigma }}_k)+J_k{\mathbf{q}}_k,$$

$$\widetilde{\mathbf{r}}(\boldsymbol{\sigma})=\left[\begin{array}[]{c}\operatorname{Re}(\mathbf{r}(\boldsymbol{\sigma}))\\ \operatorname{Im}(\mathbf{r}(\boldsymbol{\sigma}))\end{array}\right]\in\mathbb{R}^{2m},\qquad\widetilde{J}(\boldsymbol{\sigma})=\left[\begin{array}[]{c}\operatorname{Re}(J(\boldsymbol{\sigma}))\\ \operatorname{Im}(J(\boldsymbol{\sigma}))\end{array}\right]\in\mathbb{R}^{2m\times n},$$

$$\underset{\mathbf{q}\in {\mathbb{R}}^n}{\min }\Vert \tilde{\mathbf{r}}({\mathbf{\sigma }}_k)+{\tilde{J}}_k\mathbf{q}\Vert .$$

$${\mathbf{\sigma }}_{k+1}={\mathbf{\sigma }}_k+{\alpha }_k{\mathbf{q}}_k,$$

$$\underset{\mathbf{q}\in \mathcal{S}}{\min }\Vert L\mathbf{q}{\Vert }^2,\mathcal{S}=\{\mathbf{q}\in {\mathbb{R}}^n:\mathbf{q}=\arg \min \Vert {\tilde{J}}_k\mathbf{q}+{\tilde{\mathbf{r}}}_k\Vert \},$$

$$\widetilde{J}_{k}=U\left[\begin{array}[]{ll}\Sigma&\quad 0\\ 0&\quad I_{n-p}\end{array}\right]Z^{-1},\qquad L=V\left[\begin{array}[]{ll}M&\quad 0\end{array}\right]Z^{-1},$$

$${\mathbf{q}}_k^{(\ell )}=\sum _{i=p-\ell +1}^p\frac{{\mathbf{u}}_i^T{\tilde{\mathbf{r}}}_k}{{\gamma }_i}{\mathbf{z}}_i+\sum _{i=p+1}^n({\mathbf{u}}_i^T{\tilde{\mathbf{r}}}_k){\mathbf{z}}_i,$$

$$S_{\tau }(\mathbf{q})=\sum _{r=1}^p\frac{{\left(\frac{(L\mathbf{q}{)}_r}{\tau q_r}\right)}^2}{{\left(\frac{(L\mathbf{q}{)}_r}{\tau q_r}\right)}^2+{ϵ}^2},$$

$${\left(\frac{(L\mathbf{q}{)}_k}{\tau q_k}\right)}^2<1,$$

$$S_{\tau }(\mathbf{q})\approx \Vert D_{\tau ,k}L\mathbf{q}{\Vert }^2,$$

$$(D_{\tau ,k}{)}_{i,i}=\frac{1}{\tau ({\mathbf{q}}_{k-1}{)}_r}{\left[{\left(\frac{(L{\mathbf{q}}_{k-1}{)}_r}{\tau ({\mathbf{q}}_{k-1}{)}_r}\right)}^2+{ϵ}^2\right]}^{-\frac{1}{2}}.$$

$${\tilde{\mathbf{b}}}_{\delta }\approx \tilde{\mathbf{M}}(\mathbf{\sigma })+\tilde{J}(\mathbf{\sigma })\mathbf{\delta },$$

$$\delta \tilde{\mathbf{b}}={\tilde{\mathbf{b}}}_{\delta }-\tilde{\mathbf{b}}\approx \tilde{J}(\mathbf{\sigma })\mathbf{\delta }.$$

$$\Vert \delta \tilde{\mathbf{b}}{\Vert }^2=\sum _{i=1}^{2m}(\delta {\tilde{b}}_i{)}^2=\sum _{i=1}^{2m}{\left(\tilde{J}(\mathbf{\sigma })\mathbf{\delta }\right)}_i^2,$$

$$\Vert \delta \tilde{\mathbf{b}}{\Vert }^2={ϵ}^2\sum _{i=1}^{2m}({\tilde{J}}_{i,r}{)}^2.$$

$${\Sigma }_r=\frac{\Vert \delta \tilde{\mathbf{b}}{\Vert }^2}{{ϵ}^2}=\Vert \tilde{J}{\mathbf{e}}_r{\Vert }^2,$$

$${\tilde{\mathbf{b}}}_{\delta }=\tilde{\mathbf{b}}+\frac{\delta \Vert \tilde{\mathbf{b}}\Vert }{\sqrt{m}}\mathbf{w}.$$

$${\text{SNR}}_{\delta }=10{\log }_{10}\frac{\Vert \tilde{\mathbf{b}}{\Vert }^2}{\Vert \tilde{\mathbf{b}}-{\tilde{\mathbf{b}}}_{\delta }{\Vert }^2}=60\text{dB}.$$

$${\sigma }_b=\frac{{\sigma }_w}{F},$$

Full text

1 Introduction

Frequency-domain electromagnetic induction (EMI) methods have been used extensively for near-surface characterization [23, 32, 30, 29, 25, 45, 4]. The typical measuring device is composed of two electric coils (the transmitter and the receiver) that are separated by a fixed distance and placed at a known height above the ground. The two coil axes are generally aligned either vertically or horizontally with respect to the surface of the soil. The transmitting coil generates a primary electromagnetic field $(H_{\mathrm{P}})$, which induces eddy currents in the ground, generating a secondary field $(H_{\mathrm{S}})$. The amplitude and phase components of both fields are then sensed by the receiving coil. The device stores the ratio between the secondary and the primary fields as a complex number.

Initially, raw EMI measurements were used directly for fast mapping of the electrical conductivity at specific depths, with no time spent on the inversion. Recent devices are endowed with multiple receivers (multicoil) or use alternating currents at different frequencies as probe signals (multifrequency). Because of their availability, and with the development of efficient inversion algorithms and powerful computers, EMI data are increasingly used for reliable (pseudo) three- and four-dimensional quantitative assessment of the spatial and temporal variability of the electrical conductivity in the subsurface [4, 11]. These data are usually collected with both ground-based and airborne systems [26], and they have begun to be used not only to infer the soil conductivity but also to determine its magnetic permeability [14, 10, 7]. The ratio between the secondary and the primary electromagnetic fields provides information about both the amplitude and the phase of the signal. The real part (in-phase component) is affected mainly by the magnetic permeability of the soil, while the imaginary part (out-of-phase or quadrature component) is affected mainly by its electrical conductivity. The in-phase or the quadrature component of the signal is typically inverted to reconstruct either the electrical conductivity or the magnetic permeability of the soil [9, 7, 10].

In general, an EMI survey consists of many soundings; in the case of airborne acquisition, for example, there can be hundreds of thousands. These soundings, measured with multiconfiguration devices at each specific location, are usually inverted separately and are only a posteriori stitched together in a (pseudo) two-/three-dimensional fashion. This is still a common practice even though inversion schemes based on two-/three-dimensional forward modeling are becoming available and practical for use. However, the advantages of truly two-/three-dimensional inversion with respect to one-dimensional approaches are still debatable [39]. Sometimes, in order to enforce lateral continuity between the one-dimensional inversion results, the one-dimensional approaches have been extended to incorporate spatial constraints connecting the model parameters from adjacent models [38].

As in many other fields of application, regularization is usually performed by imposing smooth constraints. However, this approach is not always consistent with the true nature of the system under investigation, as sharp interfaces might be present, for example. In these situations, a stabilizer selecting the smoothest solution among all the possible solutions the data, can produce a misleading solution, whereas compatible with the data can produce a misleading outcome, whereas a regularizing term promoting blocky solutions would definitely be more coherent with the expectations about the target. For these reasons, over the years, several approaches have been implemented to retrieve model solutions characterized by sharp boundaries. A particularly promising strategy is based on so-called minimum gradient support (MGS) stabilizers [46]. This type of stabilizer has been applied to several kinds of data and has been implemented in diverse inversion frameworks, for example, inversion of travel-time measurements [47, 40], electrical resistivity tomography [13], and spatially constrained reconstruction of time-domain electromagnetic data [24, 41, 42]. A preliminary application to frequency-domain EMI data was performed by [8]. The MGS stabilizer is a function of a focusing parameter, which influences the sparsity of the final reconstruction. Assigning a small value to this variable promotes the presence of blocky features in the solution, while a large value produces smooth results.

In this work, attention is focused on the inversion of complex-valued frequency-domain EMI data collected with different configurations by extending a numerical algorithm discussed by [9, 10, 7]. The new results are compared to those obtained by inverting the quadrature component of the signal. Additionally, the implementation of an MGS-like regularization technique is studied, coupled with the truncated generalized singular value decomposition (TGSVD) within a Gauss–Newton algorithm. For a better interpretation of the reconstructed conductivity, a possible strategy for assesing the depth of investigation (DOI) is also presented and used.

The paper is structured in six sections. Section 2 introduces the nonlinear forward modeling approach. In Sect. 3, an inversion algorithm based on the damped Gauss–Newton method coupled with a TGSVD regularization scheme, which can process the whole complex signal, is described. A minimum gradient support stabilizer and a procedure to incorporate it into the above algorithm is discussed in Sect. 4. To better evaluate the performance of the investigated inversion strategies, an approach for estimating the DOI is presented in Sect. 5, leading to the numerical experiments on synthetic and real data sets reported in Sect. 6. Section 7 concludes the paper, summarizing its content.

2 The nonlinear forward model

A forward model for predicting the EM response of the subsoil was discussed by [43]. This approach is based on Maxwell’s equations and takes into account the layered symmetry of the problem. The soil is assumed to have an $n$-layered structure below ground level ($z_{1}=0$). Each horizontal layer, of thickness $d_{k}$, ranges from depth $z_{k}$ to $z_{k+1}$, $k=1,\ldots,n-1$; the deepest layer, starting at $z_{n}$, is considered to have infinite thickness $d_{n}$; see Fig. 1. The $k$th layer is characterized by an electrical conductivity $\sigma_{k}$ and a magnetic permeability $\mu_{k}$.

The two coils of the EMI measuring device, separated by a distance $\rho$ and operating at frequency $f$ in Hz, are located at height $h$ above the ground with their axes oriented either vertically or horizontally with respect to the ground surface. Let $\omega=2\pi f$ be the angular frequency of the device, and let $\lambda$, ranging from zero to infinity, denote the depth below the ground, normalized by the intercoil distance $\rho$. As discussed by [9], if $u_{k}(\lambda)=\sqrt{\lambda^{2}+i\sigma_{k}\mu_{k}\omega}$ and $N_{k}(\lambda)=u_{k}(\lambda)/(i\mu_{k}\omega)$ are the propagation constant and the characteristic admittance in the $k$th layer, respectively, then the surface admittance $Y_{k}(\lambda)$ at the top of the layer satisfies the recursion equation

[TABLE]

for $k=n-1,n-2,\ldots,1$. The recursion relationship in Eq. (1) is initiated by setting $Y_{n}(\lambda)=N_{n}(\lambda)$ for the deepest layer. It is worth remarking that both the characteristic and surface admittances depend on the frequency through the functions $u_{k}$.

The ratio between the secondary and primary fields for the vertical ($\nu=0$) and horizontal ($\nu=1$) orientations is given by the expression

[TABLE]

where $\boldsymbol{\sigma}=(\sigma_{1},\ldots,\sigma_{n})^{\mathrm{T}}$ and $\boldsymbol{\mu}=(\mu_{1},\ldots,\mu_{n})^{\mathrm{T}}$ represent the conductivity and permeability vectors, respectively, and $J_{s}(\lambda)$ denotes the first-kind Bessel function of order $s$ [1, Sect. 4.5]. The reflection factor

[TABLE]

can be calculated by setting $N_{0}(\lambda)=\lambda/(i\mu_{0}\omega)$ and computing $Y_{1}(\lambda)$ via the recursion in Eq. (1), where $\mu_{0}$ represents the magnetic permeability of free space. The reader should note that the integrand function in Eq. (2) depends on the angular frequency $\omega$, as well as on the vectors $\boldsymbol{\sigma}$ and $\boldsymbol{\mu}$, through the functions $N_{0}(\lambda)$ and $Y_{1}(\lambda)$, which define the reflection factor in Eq. (3).

The complex-valued functions $M_{0}$ and $M_{1}$ can be expressed in a more compact form in terms of the Hankel transform [1, Sect. 4.11]

[TABLE]

as follows

[TABLE]

In general, EMI devices record both the real (in-phase) and imaginary (quadrature) parts of the field ratio.

3 The inversion scheme

To investigate different depths and be able to infer both the electrical conductivity and the magnetic permeability profiles for each measurement location, it is necessary to record EMI data in different configurations. Therefore, the measurements can be acquired with different intercoil distances, operating frequencies, and heights. To further increase the information content in the data, arbitrary combinations of those configurations can be utilized. Hence, by indicating with $m_{\omega}$, $m_{h}$, and $m_{\rho}$, respectively, the number of used frequencies, heights, and intercoil distances, the total number of data measurements, $b_{tij}^{\nu}$ (with $t=1,\ldots,m_{\rho}$, $i=1,\ldots,m_{h}$, $j=1,\ldots,m_{\omega}$, and $\nu=0,1$) available at each sounding location is $m=2m_{\rho}m_{h}m_{\omega}$. Of course, the ultimate goal is to retrieve an estimate of the electrical conductivity vector $\boldsymbol{\sigma}$ and the magnetic permeability vector $\boldsymbol{\mu}$ that produce the best approximation $M_{\nu}(\boldsymbol{\sigma},\boldsymbol{\mu})\approx b_{tij}^{\nu}$ of the observations.

In the following, it is assumed that the contribution of the permeability distribution to the overall response is negligible (i.e., $\mu_{k}=\mu_{0}$ for $k=1,\ldots,n$), so that the measurements are considered to be sensitive merely to the conductivity values. However, in principle, the regularization approach discussed here can also be extended to include the inversion for the $\boldsymbol{\mu}$ components. This would require fixing an estimate for the conductivity and determining the permeability from the data [7], or computing both quantities by considering the readings defined in Eq. (2) as functions of $2n$ variables $\sigma_{k}$ and $\mu_{k}$, for $k=1,\ldots,n$.

To retrieve the conductivity values $\sigma_{k}$ ($k=1,\ldots,n$) associated with the best data approximation, frequency-domain observations $b_{tij}^{\nu}$ can be rearranged in a data vector $\mathbf{b}\in\mathbb{C}^{m}$; the same is true for the corresponding calculated responses $M_{\nu}$, which can be represented as a vector $\mathbf{M}(\boldsymbol{\sigma})\in\mathbb{C}^{m}$. Disregarding, for the moment, the ill-posedness of the problem, the best approximation $\boldsymbol{\sigma}^{\ast}$ can be found by minimizing the Euclidean norm of the residual $\mathbf{r}(\boldsymbol{\sigma})$, that is,

[TABLE]

where $\mathbf{r}(\boldsymbol{\sigma})=\mathbf{b}-\mathbf{M}(\boldsymbol{\sigma})$ takes complex values.

The adopted inversion scheme is based on the Gauss–Newton method, consisting of the iterative minimization of the norm of a linear approximation of the residual. Hence, assuming the Fréchet differentiability of $\mathbf{r}(\boldsymbol{\sigma})$

[TABLE]

where $\boldsymbol{\sigma}_{k}$ is the current approximation, and $J_{k}=J(\boldsymbol{\sigma}_{k})\in\mathbb{C}^{m\times n}$ is the Jacobian of $\mathbf{r}(\boldsymbol{\sigma})=(r_{1}(\boldsymbol{\sigma}),\ldots,r_{m}(\boldsymbol{\sigma}))^{\mathrm{T}}$, defined by $\left[J(\boldsymbol{\sigma})\right]_{ij}=\frac{\partial r_{i}(\boldsymbol{\sigma})}{\partial\sigma_{j}}$, with $i=1,\ldots,m$ and $j=1,\ldots,n$.

To determine the step length $\mathbf{q}_{k}$, as is customary, the real and the imaginary parts of the arrays involved in the computation are stacked,

[TABLE]

and the following linear least squares problem is solved:

[TABLE]

The vectors $\widetilde{\mathbf{M}}(\boldsymbol{\sigma}),\widetilde{\mathbf{b}}\in\mathbb{R}^{2m}$ are similarly defined. In fact, this approach shows that inverting the full complex signal doubles the number of available data measurements.

The analytical expression of the Jacobian was derived by [9, 7]. In the same papers, it was proven that such an expression is both more accurate and faster to compute than its finite difference approximation.

To ensure the convergence, while also enforcing the positivity of the solution, the Gauss–Newton scheme can be implemented by incorporating a damping factor. The iterative method becomes

[TABLE]

where the step size $\alpha_{k}$ is determined according to the Armijo–Goldstein principle [3], with the additional constraint that the solution must be positive ($\boldsymbol{\sigma}_{k+1}>0$) at every iteration. This choice of $\alpha_{k}$ ensures the convergence of the iterative method, provided that $\boldsymbol{\sigma}_{k}$ is not a critical point, as well as the physical meaningfulness of the solution.

The inversion of frequency-domain EMI measurements is known to be ill-posed [46], so that the linearized problem in Eq. (5) is severely ill-conditioned for each value of $k$. A strategy to tackle the ill-posedness and find a unique and stable solution consists of including available physical information in the inversion process via regularization.

A way to incorporate such a priori information in the process is to couple the original least squares problem, expressed by Eq. (5), with an additional term, leading to the new minimization problem

[TABLE]

where $L$ is a suitable regularization matrix, which defines the $L$-weighted minimum norm least squares solution [3]. The lower the value of $\|L\mathbf{q}\|$ at the selected model, the better the matching between the solution and the a priori information. By far, the most commonly used regularization matrices favor solutions that are smoothly varying (either spatially or with respect to a reference model). In such cases, $L$ is often chosen to be the identity matrix or a discrete approximation of the first or second spatial derivative.

To cope with the ill-conditioning of the problem, if $L$ is the identity matrix, the minimum norm solution of Eq. (5) at each iteration of the Gauss–Newton method can be computed by the truncated singular value decomposition (TSVD) of the Jacobian $\widetilde{J}_{k}$ [17]. If $L\in\mathbb{R}^{p\times n}$, with $p\leq n$, is different from the identity matrix, then, assuming the intersection of the null spaces of $\widetilde{J}_{k}$ and $L$ to be trivial, Eq. (7) can be solved by means of the TGSVD. SVD (TGSVD).

The following discussion will be limited to the case $2m\geq n\geq p$, as the situation characterized by $2m<n$ can be treated in a similar manner. In this case, the generalized singular value decomposition (GSVD)

[TABLE]

where $U\in\mathbb{R}^{2m\times n}$ and $V\in\mathbb{R}^{p\times p}$ have orthonormal columns, $Z\in\mathbb{R}^{n\times n}$ is nonsingular, and $\Sigma=\mathop{\operator@font diag}\nolimits[\gamma_{1},\ldots,\gamma_{p}]$, $M=\mathop{\operator@font diag}\nolimits[\xi_{1},\ldots,\xi_{p}]$ are diagonal matrices with nonnegative entries, normalized so that $\gamma_{i}^{2}+\xi_{i}^{2}=1$, for $i=1,\ldots,p$.

The TGSVD solution of Eq. (7), with parameter $\ell=0,1,\ldots,p$, is then defined as

[TABLE]

in which $\mathbf{u}_{i}$ and $\mathbf{z}_{i}$ are the columns of $U$ and $Z$, respectively. Removing the first $\ell$ terms in the first summation of Eq. (9) eliminates the contribution associated with the smallest $\gamma_{i}$. This leads to an approximated solution that is more stable, so $\ell$ acts as a regularization parameter. For an implementation of the above discussed inversion algorithm, see [6].

At each step of the Gauss–Newton iteration, the regularized minimizer of Eq.  (4) is found by solving Eq. (7) through the TGSVD defined in Eq. (9) for a fixed value of the regularization parameter $\ell$. Thus, the solution at convergence $\boldsymbol{\sigma}^{(\ell)}$ depends on the specific choice of $\ell$. If a reliable estimation of the noise level in the data is available, the regularization parameter can be chosen by means of the discrepancy principle, which requires that the data fitting must match the noise level in the data. In contrast, other heuristic strategies can be adopted. One of the most frequently used approaches is the L-curve criterion [17], based on the reasonable assumption that the most appropriate choice for the regularization parameter is the one that guarantees the optimal trade-off between the best data fitting and the most appropriate stabilization. A comparison of different strategies for estimating the regularization parameter was presented by [34]. Clearly, the inverse problem can also be approached in a probabilistic framework; in this case, the solution consists of a posterior probability distribution that naturally provides an estimation of the uncertainty of the reconstruction [21, 36]. The empirical Bayes method presented by [15] supplies a method for estimating the regularization parameter, in addition to the overall model uncertainties.

The forward model $\mathbf{M}(\boldsymbol{\sigma})$, described in Sect. 2, is strongly nonlinear and nonconvex, so its inversion is rather sensitive to the starting solution $\boldsymbol{\sigma}_{0}$ used to initialize the iterative method defined by Eq. (6). In our experience, when the noise in the data is normally distributed and relatively small, as in the numerical experiments on synthetic data of Sects. 6.1 and 6.2, any reasonable choice of $\boldsymbol{\sigma}_{0}$ converges in general to a solution that may not be the best possible but still maintains physical significance. In contrast, when the noise type is consistent with real-world applications (see Sect. 6.3), an accurate choice of $\boldsymbol{\sigma}_{0}$ becomes essential for obtaining meaningful results. In this paper, the simple procedure of repeating the computation with a few different constant starting models was adopted, selecting the solution which produced the minimal residual at convergence. In the future, we plan to investigate the application of global optimization techniques [19] in order to reduce the importance of a priori information for choosing the initial solution. Such global strategies incur a high computational cost, but they are gaining increasing popularity [15] because high-performance parallel computers are now commonly available.

4 MGS regularization

Both the estimation of the regularization parameter and the choice of the stabilizing term, which incorporates the available a priori information on the solution, play an essential role in the accuracy of the final result. When the solution is known (or assumed) to be smooth, a common choice for $L$ is the discrete approximation of either the first or second spatial derivative of the conductivity distribution. Following the same rationale, to maximize the spatial resolution of the result, whenever the solution is expected to exhibit a blocky structure, a stabilizer promoting the sparsity of the computed solution and the retrieval of sharp interfaces should be used instead.

An example of such stabilizers is the minimum gradient support (MGS) approach [33, 46, 37]. It consists of replacing the term $\|L\mathbf{q}\|^{2}$ in Eq. (7) with

[TABLE]

where $L$ is a regularization matrix, while $\tau$ and $\epsilon$ are free parameters. As can be immediately observed, Eq. (10) depends only upon the product $\tau\epsilon$, so in the following, $\epsilon=1$ is fixed, and only $\tau$ varies.

[37] introduced a generalized stabilizing term that reproduces, for particular values of two parameters, the $L_{2}$ and $L_{1}$ norms, the MGS stabilizer, and others. The authors showed that for small values of $\tau$ Eq. (10) approximates an approach proposed by [22] that minimizes the pseudonorm $\|L\mathbf{q}\|_{0}$, that is, the number of nonzero entries in the vector $L\mathbf{q}$; see also [33], [40], and [44]. Therefore, the nonlinear regularization term $S_{\tau}(\mathbf{q})$ favors the sparsity of the solution and the reconstruction of blocky features. If $L$ is chosen to be the discretization of the first derivative $D_{1}$, the stabilizer introduced in Eq. (10) selects the solution update corresponding to minimal nonvanishing spatial variation. This is the origin of the “minimum gradient” descriptor. Its clear advantage is that it can mitigate the smearing and blurring effects of the more standard smooth regularization strategies.

The parameter $\tau$ determines how each term in Eq. (10) affects the overall value. In particular, as discussed by [41], the model updates with

[TABLE]

are weakly penalized, as the corresponding term in Eq. (10) is small, while updates with derivatives larger than the threshold defined by $\tau q_{k}$ may give a contribution close to 1. Thus, the MGS stabilizer penalizes the occurrences of variations larger than the threshold $\tau q_{k}$, rather than the magnitude of the variations itself. This, in turn, favors spatially sparse updates. The threshold defining when an update is to be considered large enough to be penalized is dynamically chosen, via the parameter $\tau$, as a fraction of the actual conductivity update $q_{k}$. In conclusion, the MGS stabilizer allows for reconstruction of sharp features while maintaining the smoothing effect of the regularization $L$ for small variations in the conductivity updates.

In general, applying the nonlinear regularizing term in Eq. (10) to a linear least squares problem requires a larger computational effort than the standard first/second derivative approach. In this case, the least squares problem is nonlinear itself, so Eq. (10) can be treated by the main iterative algorithm: it is linearized at each step of the Gauss–Newton method by evaluating the terms in the denominator at the previous iteration $\mathbf{q}_{k-1}$. At each step, Eq. (5) is solved by Eq. (7), replacing $\|L\mathbf{q}\|^{2}$ by the approximation

[TABLE]

where $D_{\tau,k}$ is the diagonal matrix with elements

[TABLE]

In the numerical simulation described in Sect. 6, the regularization matrix $L$ is always $D_{1}$.

Every time the forward model is linear, MGS regularization leads to a convex problem; this was proved, for example, by [33]. For nonlinear forward problems, such as the one discussed in the present study, the further nonlinearity introduced by the MGS stabilizer emphasizes the nonconvex nature of the data fitting problem; see the discussion at the end of Sect. 6.3. As already noted, the nonconvexity issue could be approached through global optimization algorithms [19], but approaches that employ available prior information for the starting model selection are still of some practical interest for their efficiency.

5 Depth of investigation

The depth of investigation (DOI) usually refers to the depth below which data collected at the surface are not sensitive to the physical properties of the subsurface. In short, the DOI provides an estimation of the maximum depth that can be investigated from the surface, given a specific device (in a specific configuration) and the physical properties of the subsoil. Without a DOI assessment, it is difficult to judge wether the reconstruction at depth is produced by the data or is merely an effect of the specific choice of the starting model and/or the inversion strategy.

One way to assess the DOI can be based on the skin depth calculation, function of the frequency and the medium conductivity [27]. Alternative methods rely on the study of the variability of the solution as a function of the starting model. For example, [28] discussed the effectiveness of inverting the data with very different initial half-space conductivities and subsequently comparing the results to determine the depth threshold up to which the reconstruction was influenced by the data.

Similar to the strategy by [5], the approach proposed here is based on the integrated sensitivity matrix, as discussed by [46]. Hence, in the following, the DOI is defined as the depth where, for each individual sounding, the integrated sensitivity values drop below a certain threshold. With the aim of studying the sensitivity of the data vector $\widetilde{\mathbf{b}}=\widetilde{\mathbf{M}}(\boldsymbol{\sigma})$ to a perturbation vector $\boldsymbol{\delta}$, the perturbed data $\widetilde{\mathbf{b}}_{\delta}=\widetilde{\mathbf{M}}(\boldsymbol{\sigma}+\boldsymbol{\delta})$ is taken into account. The linearized version of the problem produces the approximation

[TABLE]

which implies

[TABLE]

Then,

[TABLE]

where $\delta\widetilde{b}_{i}$ denotes the $i$th component of $\delta\widetilde{\mathbf{b}}$, $i=1,\ldots,2m$.

Now, assuming $\boldsymbol{\delta}=\varepsilon\mathbf{e}_{r}$, where $\mathbf{e}_{r}\in\mathbb{R}^{n}$ has zero entries except $(\mathbf{e}_{r})_{r}=1$, and denoting by $\widetilde{J}_{i,r}$ the $(i,r)$ entry of the Jacobian $\widetilde{J}(\boldsymbol{\sigma})$, the norm of the perturbation takes the form

[TABLE]

Then, the integrated sensitivity of the data is defined by

[TABLE]

where $\widetilde{J}\mathbf{e}_{r}$ denotes the $r$th column of the Jacobian matrix. This measure represents the relative sensitivity of the data vector to a perturbation in the conductivity of the ground layer at depth $z_{r}$.

When $\Sigma_{r}$ decreases significantly with respect to $\Sigma_{1}$, that is, when $\Sigma_{r}<\eta\Sigma_{1}$ for a fixed tolerance $\eta$, the recovered conductivity for the $r$th layer is not strictly related to the data or, thus, to the physical properties of the subsoil. Then, the depth $z_{r}$, at which the reduction $\Sigma_{r}<\eta\Sigma_{1}$ occurs is where the DOI is set. Evidently, there is some degree of arbitrariness in the choice of the threshold $\eta$ for the decrease in $\Sigma_{r}$.

6 Numerical experiments

Numerical experiments were run on a Xeon Gold 6136 computer running the Debian GNU/Linux operating system, using a MATLAB software package implementing the algorithms described in this paper [6]. The software is available at the web page http://bugs.unica.it/cana/software/ as the FDEMtools package.

In the numerical tests illustrated in this section, the electrical conductivity is determined starting from synthetic and experimental data sets under the assumption that the magnetic permeability can be approximated by that of empty space. The results obtained by processing the quadrature component of the signal will be compared with those derived from the complex signal in its entirety. Additionally, the MGS sparsity promoting strategy will be compared with traditional smooth stabilizers.

6.1 One-dimensional synthetic data

A synthetic data set is generated by representing the conductivity as a function of depth by the following test functions

• •

Gaussian profile: $\sigma_{1}(z)=e^{-(z-1.2)^{2}}$,

• •

Step profile: $\sigma_{2}(z)=\begin{cases}0.2,&\quad z<1,\\ 1,&\quad z\in[1,2],\\ 0.2,&\quad z>2.\end{cases}$

Assuming the magnetic permeability to be that of free space ($\mu=\mu_{0}$) and the subsoil to be divided into 60 layers ($n=60$) between $z=0$ m and $z=3.5$ m, the forward model described in Sect. 2 is applied in conjunction with a chosen device configuration to reproduce the instrument readings. Since the experimental data studied in Sect. 6.3 was recorded by a CMD Explorer ($\rho=1.48,2.82,4.49$ m; $f=10$ kHz), the instrument readings are constructed according to such configuration, assuming the measurements were acquired at heights $h=0.9,1.8$ m. This leads to six readings for each coil orientation ($m_{h}=2$, $m_{\rho}=3$, $m_{\omega}=1$).

To simulate experimental errors, given a vector $\mathbf{w}$ with normally distributed entries having zero mean and unit variance, the perturbed data vector $\widetilde{\mathbf{b}}_{\delta}$ is determined from the exact data $\widetilde{\mathbf{b}}$ by the formula

[TABLE]

This implies that $\|\widetilde{\mathbf{b}}-\widetilde{\mathbf{b}}_{\delta}\|\approx\delta\|\widetilde{\mathbf{b}}\|$. In the computed example, $\delta=10^{-3}$. The equivalent signal-to-noise ratio (in decibels) is

[TABLE]

This noise level is unrealistic in real-world applications, in which the experimental error may be non-Gaussian and highly correlated. Here, the aim is to test the performance of the inversion algorithm in an ideal situation.

For all numerical experiments, the regularization parameter $\ell$ (see Eq. (9)) is chosen by applying the discrepancy principle, as the noise is Gaussian and its level is exactly known.

In Fig. 2, the results obtained by the inversion of the complex signal are compared with those obtained by inverting only the quadrature component. In this experiment, the smooth test profile $\sigma_{1}(z)$ and the regularization term $L=D_{2}$, the discretization of the second derivative, are adopted. The graphs in the top row show the reconstruction of the conductivity when both orientations of the coils are used, that is, the data set is composed of 12 readings. The top-left graph represents the solution obtained by inverting the complex data, while the top-right image reports the reconstruction resulting from inverting just the quadrature component of the signal. It is clear that the inversion of the complex signal provides better results.

The graphs in the bottom row of Fig. 2 show the results for the same settings but processing data only for the vertical orientation. The reconstructions are very similar to those in the top row, showing that repeating the data acquisition with two orientations of the coils does not necessarily produce sensibly better results, especially if complex measurements are processed.

To investigate the performance of the algorithm in the presence of strong noise in the data, the above computation was repeated raising the Gaussian noise level to $\delta=0.2$, that is, $20\%$ of the signal, corresponding to $\text{SNR}_{\delta}=14\,\text{dB}$. This noise level has been considered realistic in urban sites [20], but it is much larger than that encountered in average real-world standards [12]. The remaining parameters are kept unchanged, i.e., $L=D_{2}$, the discrepancy principle is used to select the regularization parameter $\ell$, and the 12 readings correspond to both orientations of the coils. The results, displayed in Fig. 3, show that processing the complex signal leads to reasonably localizing the maximum conductivity. In contrast, the quadrature component alone does not allow computation of a meaningful reconstruction. This example underlines the importance of processing the whole complex data set when experimental measurements are considered.

Figure 4 displays the results concerning the second synthetic example, namely, the reconstruction of the discontinuous test profile $\sigma_{2}(z)$ for the electrical conductivity. The same CMD Explorer configuration as before is considered, but the data are generated only for the vertical orientation of the coils; the noise level is $\delta=10^{-3}$. The graphs in the top row illustrate the performance of the smooth regularizing matrix $L=D_{1}$ for both the complex signal and the quadrature component. The bottom row displays the same results for the nonlinear regularizing term $S_{\tau}(\mathbf{q})$, after setting $L=D_{1}$ in Eq. (10). The results in Fig. 4 show that the MGS stabilizer is able to approximate the presence of sharp boundaries with good accuracy in the model function. Again, processing the complex signal produces more accurate results.

6.2 Pseudo two-dimensional synthetic data

The example described in this section concerns the reconstruction of a series of one-dimensional models (more precisely, 50 soundings along a 10 m straight line) characterized by an abrupt change in conductivity (from 0.5 S/m to 2 S/m) occurring at an increasing depth. At the top of Fig. 5, the one-dimensional models are depicted side by side in a pseudo two-dimensional fashion. This facilitates comparison and assessment of the effectiveness of the methods as the depth of the conductivity transition varies. The synthetic data simulate an acquisition performed by CMD Explorer ($\rho=1.48,2.82,4.49$ m, $f=10$ kHz), with two orientations of the coils and two measurement heights $h=0.9,1.8$ m. The data values are finally perturbed by uncorrelated Gaussian noise with standard deviation $\delta=10^{-3}$. To simulate an experimental setting, in which often no information on the noise level is available, the regularization parameter is estimated in each one-dimensional inversion by the L-curve criterion.

The left-hand side graphs of Fig. 5 show the smooth inversion results corresponding to $L=I,D_{1},D_{2}$, obtained with a 60-layer parameterization up to a depth of 3.5 m, with layers of constant thicknesses and a homogeneous 0.5 S/m starting model. The right-hand side graphs of Fig. 5 correspond to sharp MGS inversions with three different values of the focusing parameter $\tau=10^{-1},10^{-2},10^{-4}$.

From these results, it is evident that the smooth inversion for $L=D_{1},D_{2}$ produces acceptable results, but with an excess of smoothness. Indeed, it correctly retrieves the transition between the upper resistivity layer and the lower conductive background, but the transition is not well identified in space. It is worth mentioning that the data for each sounding location are generated independently and that during the inversion, no lateral constraints are imposed, so the inversion proves to be quite stable.

Not surprisingly [41, 13], the MGS stabilizer with a large $\tau$ produces results very similar to the smooth results. Decreasing the value of the focusing parameter $\tau$ in Eq. (10) corresponds to penalizing the number of small vertical relative variations in the conductivity updates, as $\tau$ defines the variability range, allowing the derivative update to be considered “relevant” for the MGS stabilizer summation. The sparsity-enhancing effects of the MGS stabilizer are particularly effective when $\tau=10^{-2}$, where the discontinuity in the solution is more clearly identified. When further reducing the focusing parameter, for example, for $\tau=10^{-4}$, the reconstructions start to exhibit unrealistic blocky features. This is even more clear in Fig. 6, where the reconstruction of a single sounding (the 30th column of the two-dimensional synthetic model of Fig. 5) is reported, comparing the one-dimensional smooth reconstruction corresponding to $L=D_{1}$ (top left) with that of the MGS stabilizer for $\tau=10^{-1}$ (top right), $\tau=10^{-2}$ (bottom left), and $\tau=10^{-4}$ (bottom right).

A drawback of the MGS inversion process is the instability of the reconstructed solution. Close one-dimensional reconstructions, that is, close columns of the “two-dimensional” solution, are sometimes very different. This may be due to the concurrent effects of an incorrect estimation of the regularization parameter and the nonconvexity of the nonlinear regularized objective function, leading to solutions becoming trapped in inferior local optima in some spatial locations. The reconstruction may be forced to be more regular by imposing lateral constraints, for example, migrating additional pieces of information from adjacent models [41]. This will be the subject of future work.

6.3 Real Survey

The proposed algorithm was tested on an experimental data set collected with a multiconfiguration EMI device at the Molentargius Saline Regional Nature Park, located east of Cagliari in southern Sardinia, Italy, and displayed in Fig. 7. At this site, [16] investigated the flow dynamics associated with freshwater injection in a hypersaline aquifer through hydrogeophysical monitoring and modeling, using five 20 m deep boreholes (Fig. 7b, c). The park is a wetland characterized by the presence of very salty groundwater, with salinity levels as high as 3 times the NaCl concentration of seawater due to the long-term legacy of infiltration of hypersaline solutions from nearby salt pans (Fig. 7a) dating back to Roman times. This site appears to be ideal for testing the MGS regularized inversion procedure, as the very high electrical conductivity of the aquifer makes the unsaturated/fully saturated soil interface a sharp electrical conductivity interface.

Prior to the freshwater injection experiment, laboratory petrophysical measurements and different surface, in-hole, and cross-hole electrical resistivity surveys were carried out to characterize the background of unsaturated/saturated sedimentary succession dominated by sands. Figure 8 shows some results of these preliminary investigations, which were used as a reference against which to compare the reliability of the inversion results.

Groundwater conductivity ($\sigma_{w}$) logs recorded in boreholes (see Fig. 8a) enable two zones to be discriminated, with a transitional 2 m thick layer in between: (1) from the water table at 5.2 m depth to a depth of 7.5 m, the water electrical conductivity is approximately 2 S/m, and (2) below 9.5 m depth, the water electrical conductivity reaches 18.5 S/m.

Using a Terrameter SAS Log (ABEM Instrument) with an electrode separation of 64 inches, a long normal resistivity log was carried out in borehole #1 to estimate the bulk conductivity of the fully saturated soil. Figure 8b shows bulk conductivities $\sigma_{b}$ calibrated with the values (red dots) obtained from Archie’s empirical relationship [2]

[TABLE]

where $\sigma_{w}$ is the groundwater conductivity and $F=\phi^{-m}$ is the formation factor, which is a function of the porosity $\phi$ and the cementation factor $m$. The specific values for $F$ were measured from soil samples from borehole #1 in the laboratory and are shown in Fig. 8c, together with the values of the porosities of the soil samples (in brackets). These measured bulk conductivities are apparent conductivities and are representative of a cylindrical volume with a radius of $\sim 1.5$ m around the borehole. They reach values of up to 4 S/m, but they could be overestimated due to the very high conductivity of the water present in the borehole, which acts as a preferential path for the current.

Figure 8d shows the cross-hole conductivity image resulting from the inversion of apparent electrical conductivities measured with a bipole–bipole electrode configuration (one current and one potential electrode placed in each borehole). Black diamonds denote the position of the electrodes, and the blue line shows the groundwater table at 5.2 m below the ground surface. Above the water table, the electrical conductivity is low, ranging between 1 and 10 mS/m, while in the saturated zone it is very high, and vertical changes due to the layering of lithologies are not visible. A gradual change to lower conductivities can be seen only in the upper part, just below the water table. This is consistent with the water conductivity (Fig. 8a) and bulk conductivity (Fig. 8b) logs. The conductivity reaches its highest value below 9.5 m depth, even if it is slightly smaller than the highest bulk conductivity, and so it is probably underestimated. This is also an expected feature for the lack of resolution associated with measurements with large electrode spacings.

EMI data were collected along a 200 m straight-line path (Fig. 7b) with a topographic elevation varying from 1.6 m at the southeastern end to 5.7 m at the northwestern end, using CMD Explorer (Gf-Instruments). This system operates at a frequency of 10 kHz and has one transmitter coil paired with three coplanar receiver coils at 1.48, 2.82, and 4.49 m from the transmitter, allowing three simultaneous measurements of the apparent soil electrical conductivity using vertical (VCP, vertical coplanar) or horizontal (HCP, horizontal coplanar) dipole configurations. Two surveys were carried out along the same profile. Data were recorded in continuous mode, with a 0.5 s time step, and the system was carried at a height of 0.9 m above the ground, first using the HCP and then the VCP dipole configurations. Measurement locations (UTM coordinates) were logged using a Trimble differential GPS receiver able to ensure submeter accuracy. Before merging the HCP and VCP data sets, prior to the inversion, they were spatially resampled at 0.5 m intervals from a common starting point to ensure the same number of equally spaced measurement points. This allowed us to set up a data set consisting of a series of 400 geometric depth soundings with six complex (quadrature and in-phase components) CMD Explorer responses each (Fig. 9), suitable for imaging the water table and recovering (by inversion) the recover (by inversion) the soil electrical conductivities along the surveyed profile. At the southeastern end of the survey line, both quadrature and in-phase responses show higher values than those recorded along the remaining part, as they were recorded with sensors (transmitter and receiver coils) closest to the water table. The data set is available at the web page http://bugs.unica.it/cana/datasets/.

The complex response recorded at each sounding point was inverted individually to infer the electrical conductivity depth profile using the smooth inversion scheme described in Sect. 3, with the regularization terms $L=D_{1}$, $L=D_{2}$, and the MGS regularization described in Sect. 4 with focusing parameter $\tau=10^{-4}$, and $L=D_{1}$. For all one-dimensional inversions, the same homogeneous $0.07$ S/m starting model was used, and the regularization parameter was estimated by the L-curve criterion. The discrepancy principle might be used for choosing the regularization parameter, if a reliable estimation of the noise level were available. This approach was not pursued in these experiments because our experience suggests that the noise in EM data is seldom equally distributed with respect to varying the device configuration; see, for example, [9, Fig. 10].

Obtaining a noise estimate, even if not essential to perform the computation, could be useful to better characterize the experimental setting. Our impression is that the available data set, displayed in Fig. 9, is not sufficient to obtain a trustable noise estimate since repeated measurements in the same geographical location are missing. In principle, the linear regularization procedure introduced by [18] and [31] could be adapted to this task. The method is based on comparing two regularization techniques, for example, TGSVD and Tikhonov, to select the regularization parameter. This result is then used to estimate the noise level in the data. We are currently working on the extension of this method to nonlinear regularization.

The resulting one-dimensional models, with 100 layers to a depth of 10 m below the ground surface ($d_{k}=0.1\,\text{m}$), are stitched together and plotted as a pseudo two-dimensional section in Fig. 10. In each section, the DOI is also plotted to indicate the maximum depth at which the recovered conductivity is still related to data and not a numerical artifact. The DOI is represented by the black curves close to the lower boundary of each section and was estimated by setting the threshold value at $\eta=10^{-2}$ (see Sect. 5), which was the value that produced results consistent with the findings of the borehole investigations.

All three solutions satisfactorily capture the overall picture, although each of them appears better or worse than the others in certain respects. They clearly retrieve the unsaturated/saturated soil interface at approximately 0 m elevation; in the same way, in the southeastern part of the section, they show the same conductive anomaly due to the saltwater intrusion from the nearby third evaporation pan of the old saltworks (Fig. 7).

In the smooth solutions (Fig. 10a and 10b), the water table interface is more or less recognizable, but it is not easy to resolve its exact depth since it does not appear as a sharp interface, as it should actually be in this case. This undesirable effect, due to the imposed vertical smoothness constraints, is less noticeable for the solution obtained with the second derivative regularization term than for the other solutions. When analyzing the sections in the area between boreholes #1 and #4, indicated by magenta vertical lines, it is clear that this solution obtains a better fit to the tomography in Fig. 8d. The solution obtained with $L=D_{2}$ is also better for the absolute values of electrical conductivity, which are generally underestimated when compared with those obtained from the measurements in the boreholes. Finally, note that in both smoothed solutions, the electrical conductivities vary gradually in the lateral direction, although they have been obtained by inverting data, sounding by sounding, without any lateral constraint.

Compared with the previous results, the MGS reconstruction (Fig. 10c) is less blurred and more reliable in retrieving the sharp water table interface along the whole section. Electrical conductivities are generally consistent with those expected on the basis of the results of past surveys. In particular, to the right of distance 93 m, the one-dimensional inverted models show electrical conductivity profiles in very good agreement with those of the cross-hole tomography. In the MGS solution, however, the electrical conductivities do not vary gradually in the lateral direction, and they show sharp lateral changes that do not correspond to real features of the subsoil under investigation; see, for example, the changes indicated by the arrows in Fig. 10c. The reconstruction is particularly erratic in the south-eastern part of the section, where the nonlinearity of the forward model is amplified by the high conductivity due to the closeness to the old saltworks; see Fig. 7.

This is again an illustration of the strong dependence of the reconstruction on the initialization of the iterative method. In these experiments, the same starting model was adopted for all the data columns. The approach was successful for the first two regularization matrices, while the MGS stabilizer would require a more accurate initialization. This drawback, which was already highlighted at the end of Sect. 4, could be overcome by adopting global optimization techniques. Another possibility is to impose a correlation either between the data to be inverted corresponding to neighboring points or between the obtained one-dimensional inverse models.

Implementing lateral constraints in the reconstruction, for example, an approach based on total variation [35] would couple all the one-dimensional inversion problems into a large two-dimensional problem, requiring a suitable solution algorithm. Indeed, each linearized step of the Gauss–Newton method would lead to a least squares problem (see Eq. (5)) too large to be solved by TGSVD. Another possible approach is based on MGS regularization, followed by [41] for another kind of electromagnetic data.

7 Conclusions

Obtaining relevant solutions to inverse problems requires processing meaningful data by an effective regularization technique. In the case of EMI data, taking advantage of both the in-phase and the quadrature component of the available measurements enriches the data information content, allowing for the computation of more accurate solutions. Proper regularization consists of the formalization of a priori information via a stabilizing term. Thus, smoothing schemes might not always provide the most adequate solution. Whenever sharp interfaces are expected, it may be wiser to use regularization terms promoting the sparsity of the retrieved model. In this manuscript, a Gauss–Newton algorithm regularized by a TGSVD approach, initially designed for either real (in-phase component) or imaginary (quadrature component) data inversion, was extended to process complex measurements and to accommodate an MGS stabilizer. The performance of the new algorithm was tested on both synthetic and experimental data sets and compared with alternative approaches.

Synthetic examples over both one-dimensional and pseudo two-dimensional discontinuous conductivity profiles show that the new algorithm can provide better detail in the reconstruction than other algorithms. The MGS solution is not always preferable to the smooth one; it depends on the expectations/assumptions about the target. Nevertheless, it is also true that the focusing parameter can be selected such that the model maintains a certain degree of smoothness.

The enhanced information stemming from complex data values always improved the quality of the results. The one-dimensional inversion models produced by the complex-valued experimental data set were able to provide pseudo two-dimensional earth models, consistent with the findings of in-hole and cross-hole electrical conductivity investigations, and reliable down to the DOI.

In summary, the new one-dimensional inversion algorithm can be effectively applied to retrieve smooth and sharp electrical conductivity interfaces in hydrogeological, soil, and environmental investigations. The instability remains its chief drawback, which may be overcome by adopting global optimization techniques by developing a more reliable strategy for the regularization parameter selection and by imposing appropriate lateral constraints. This, as well as the application of Bayesian uncertainty quantification ideas, will be the subject of future work.

Acknowledgements

The authors wish to thank the reviewers for their comments, which led to improvements in the presentation. This research was supported in part by the Fondazione di Sardegna 2017 research project “Algorithms for Approximation with Applications [Acube]”, the INdAM-GNCS research project “Metodi numerici per problemi mal posti”, the INdAM-GNCS research project “Discretizzazione di misure, approssimazione di operatori integrali ed applicazioni”, the Regione Autonoma della Sardegna research project “Algorithms and Models for Imaging Science [AMIS]” (RASSR57257, intervento finanziato con risorse FSC 2014-2020 - Patto per lo Sviluppo della Regione Sardegna), the Visiting Scientist Programme 2015/2016 (University of Cagliari), and an RAS/FBS Grant (Grant No. F71/17000190002).