{'fraction_non_alphanumeric': 0.03481802654053344, 'fraction_numerical': 0.008146104322690842, 'mean_word_length': 4.687173100871731, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FDEMtools: a MATLAB package for FDEM data inversion

Electromagnetic induction surveys are among the most popular techniques for non-destructive investigation of soil properties in order to detect the presence of either ground inhomogeneities or of particular substances. This work introduces a MATLAB package for the inversion of electromagnetic data collected by a ground conductivity meter. Based on a nonlinear forward model used to describe the interaction between an electromagnetic field and the soil, the software reconstructs either the electrical conductivity or the magnetic permeability of the soil with respect to depth, by a regularized damped Gauss-Newton method. The regularization part of the algorithm is based on a low-rank approximation of the Jacobian of the nonlinear model. Both the relaxation parameter and the regularization parameter are chosen by automatic procedures. The package allows the user to experiment with synthetic data sets and different regularization strategies, in order to compare them and draw conclusions.

1. Introduction. Electromagnetic induction (EMI) techniques are often used for nondestructive investigation of soil properties which are affected by electromagnetic features of the subsurface layers, namely the electrical conductivity σ and the magnetic permeability µ. Knowing such parameters allows one to ascertain the presence of particular substances, which is essential in many important applications: hydrological characterizations [3,5], hazardous waste studies [22], archaeological surveys [21,26], precision agriculture [12,34], unexploded ordnance detection [20], etc.
A ground conductivity meter is the typical measuring instrument for frequency domain electromagnetic (FDEM) induction techniques. It is composed by two coils (a transmitter and a receiver) placed at a fixed distance. The dipoles may be aligned vertically or horizontally with respect to the ground level. An alternating sinusoidal current in the transmitter produces a primary magnetic field H P , which induces small eddy currents in the subsurface. These currents produce, in turn, a secondary magnetic field H S , which is sensed by the receiver. The ratio H S /H P of the secondary to the primary magnetic fields is measured by the device, providing information about the amplitude and the phase of the signal. The real part, or the inphase component, of the measured signal is mainly affected by the magnetic permeability of the subsoil; the imaginary part, also called the out-of-phase or quadrature component, mainly by the electrical conductivity.
In this work, we present a MATLAB package for the numerical inversion of a nonlinear model which describes the interaction between an electromagnetic field and the soil. The computation is performed by a Gauss-Newton method and regularized by means of either the truncated singular value decomposition (TSVD), or the generalized truncated singular value decomposition (TGSVD). The computation of the forward model, as well as the analytical expression of its Jacobian matrix, is performed by some functions in the package. Either the electrical conductivity or the magnetic permeability of the soil with respect to depth can be reconstructed, if an assumption can be made on the behavior of the other quantity.
The software development started during the research work which lead to [9]. In this paper, assuming the permeability of the soil to be known, the aim was to reconstruct the electrical conductivity taking as input only the quadrature component of the measurements. Here, the loop-loop device used a single scanning frequency, and multiple measurements were obtained by varying the height of the instrument above the ground and the orientation of the two coils. An inversion method was proposed, based on a low-rank approximation of the Jacobian of the forward model. Also, the analytical expression of the Jacobian with respect to the electrical conductivity was computed for the first time, and it was compared to its finite difference approximation for what regards accuracy and computing time.
The algorithm was extended to deal with multiple-frequency data sets in [11]. This approach was motivated by the availability of devices which can take simultaneous readings at multiple scanning frequencies. In the same paper, experiments were performed taking as input either the in-phase or the quadrature components of the signal, in order to investigate the possibility of extracting further information from the available data.
Later, the paper [6] focused on the identification of the magnetic permeability. The main result was to obtain stable analytical formulas for the computation of the Jacobian of the forward model with respect to the variation of the magnetic permeability. The paper numerically investigated the conditioning of the problem, as well as the performance of the analytical expression of the Jacobian in comparison to its finite difference approximation. Then, under the working assumption that the conductivity was known in advance, the reliability of the inversion algorithm was verified, and the information content of data sets obtained by varying either the scanning frequency of the device or its height above the ground was investigated. The generated data sets, contaminated by a noise level compatible with real-world applications, allowed the authors to analyze the behavior of the algorithm in a controlled setting.
Suitably extended version of the software were used in [7] to investigate the effect of a particular sparsity enhancing regularization technique, and in [3,12] to process specific experimental data sets. The algorithm was finally adapted in [8] to the inversion of the whole complex signal, rather than just one of its components. The resulting software, that is essentially the one presented in this paper, was applied to an experimental data set collected at the Molentargius Saline Regional Park, Sardinia, Italy [8].
The above mentioned references focus on a nonlinear model described, e.g., in [33], and further studied in [18]. We remark that when the soil conductivity takes relatively small values (below 0.5 S/m) a linear model can be used to predict the measurements produced by a particular device [23]. A numerical method for this model was initially proposed in [4]. The same model has recently been studied from the theoretical point of view in [10], where an optimized solution method was proposed too.
The present paper is structured as follows. In Section 2 we give a brief description of the nonlinear forward model and we introduce the inversion algorithm, together with the regularization procedure. Section 3 presents the FDEMtools package as well as its graphical user interface (GUI), describing the installation process and how to use the software by means of some numerical examples. For extensive numerical experiments, we refer to the above mentioned papers. The final Section 4 summarizes the content of the paper.
2.1. The nonlinear forward model. A forward model which predicts the data measured by a FDEM induction device, when the distribution of the conductivity and the magnetic permeability in the subsoil is known, has been described in [17].
The model assumes the soil to be layered, so that both the electrical conductivity and the magnetic permeability are piecewise constant functions. Each subsoil layer, of thickness d k (m), is characterized by an electrical conductivity σ k (S/m) and a magnetic permeability µ k (H/m), for k = 1, . . . , n; see [6]. The thickness of the deepest layer d n is considered infinite. The two coils of the measuring device are at height h above the ground, their distance is ρ.
Let u k (λ) = λ 2 + iσ k µ k ω, where ω is the angular frequency of the instrument, that is, 2π times the frequency in Hertz. The variable λ ranges from zero to infinity, it measures the ratio between the depth below the ground surface and the inter-coil distance ρ. If we denote the characteristic admittance in the k-th layer by N k (λ) = u k (λ)/(iµ k ω), then it is shown in [32] that the surface admittance Y k (λ) at the top of the same layer verifies the recursion for k = n − 1, . . . , 1. For k = n, the characteristic admittance and the surface admittance are assumed to coincide, that is, the value Y n (λ) = N n (λ) is used to initialize the recursion. Both the characteristic and the surface admittances are functions of the frequency ω via the functions u k (λ). Let us define the reflection factor as where Y 1 (λ) is computed by the recursion (2.1). The ratio of the secondary to the primary field for the vertical and horizontal orientation of the coils, respectively, are given by the Hankel transform, where J 0 , J 1 are first kind Bessel functions of order 0 and 1, respectively.
We remark here that the functions in (2.3) take complex values and that the measuring devices, in general, return both the real (in-phase) and the imaginary (quadrature) components of the fields ratio. The quadrature component, suitably scaled, can be interpreted as an apparent conductivity, while the in-phase component is related to the magnetic permeability of the soil.
Let us consider the complex residual vector as a function of the conductivities σ k and the permeabilities µ k , k = 1, . . . , n; the vector function M returns the readings predicted by the model (2.3) in the same order the values b ν tij were arranged in the vector b. To process the data set, we solve the following nonlinear least-squares problem by the Gauss-Newton method. At each step of the iterative algorithm we minimize the 2norm of a linear approximation of the residual r(σ k+1 , µ k+1 ) = r(σ k + s k , µ k + t k ), namely, for i = 1, . . . , m and k = 1, . . . , n.
The residual function r being complex-valued, we solve problem (2.5) by stacking the real and imaginary part of the residual as follows [8] The positive parameter β allows for balancing the contributions of the in-phase and the quadrature component. It is in general set to one, unless it is modified by the user on the basis of available a priori information, e.g., on the soil composition or on the noise level. So, we replace (2.5) by with J k = J(σ k , µ k ), and the iterative method becomes where J † k is the Moore-Penrose pseudoinverse of J k and α k is a damping parameter which ensures the convergence. It is determined by coupling the Armijo-Goldstein principle [2] to the positivity constraint (σ k+1 , µ k+1 ) > 0 (componentwise) [6,9].
The analytical expression of the Jacobian matrices with respect to the electrical conductivity and the magnetic permeability were computed in [9] and [6], respectively. In fact, the same papers show that the analytical Jacobian is more accurate and faster to compute, than resorting to a finite difference approximation.

Regularization.
It is well known that the minimization problem (2.4) is extremely ill-conditioned. In particular, it has been observed in [6,9] that the Jacobian matrix J has a large condition number virtually for all values of the variables σ and µ in the solution domain. A common remedy to address ill-conditioning consists of resorting to regularization, i.e., replacing the linearized least-squares problem (2.6) by a nearby problem, whose solution is less sensitive to the propagation of the errors in the data.
A regularization method which particularly suits the problem, given the size of the matrices involved, is the truncated singular value decomposition (TSVD). The best rank approximation ( ≤ κ = rank( J k )) to the Jacobian matrix, according to the Euclidean norm, is easily obtained by the SVD decomposition. This factorization allows us to replace the ill-conditioned Jacobian J k by a well-conditioned low-rank matrix; see [15].
When some kind of a priori information for the problem is available, e.g., the solution is a smooth function, it is useful to introduce a regularization matrix L ∈ R p×2n (p ≤ 2n), whose null space approximately contains the sought solution [29]. Under the assumption N ( J k ) ∩ N (L) = {0}, problem (2.6) is replaced by Very common choices for L are the discretization of the first or second derivative operators. Let the generalized singular value decomposition (GSVD) [13] of the matrix pair ( J k , L) be where U and V are matrices with orthogonal columns u i and v i , respectively, Z is a nonsingular matrix with columns z i , and Σ J , Σ L are diagonal matrices with diagonal entries c i and s i . Under the assumption that m = 2m ρ m h m ω < 2n, common in the generality of cases, the truncated GSVD (TGSVD) solution q ( ) (see [15] for details) can be written as where κ = rank( J k ), = 1, . . . , κ + p − 2n is the regularization parameter, and r k = r(σ k , µ k ).
The resulting regularized damped Gauss-Newton method reads with fixed and α k determined at each step as explained in Section 2.2. We iterate until The solution at convergence is denoted by (σ ( ) , µ ( ) ). The choice of the regularization parameter is crucial. In real applications, experimental data are affected by noise. We express the data vector in the residual function (2.4) by b = b + e, where b contains the exact data and e is the noise vector. If the noise is Gaussian and an accurate estimate of e is available, we determine by the discrepancy principle [15,25]. When such an estimate is not at hand, we adopt heuristic methods, such as the L-curve, in the implementation presented in [14], or an hybrid method based on the quasi-optimality criterion, described in [24,29]; see [19,27,29] for a review of similar methods.
Of course, the choice of the regularization matrix L has strong effects on the results, as it incorporates the available a priori information regarding the solution. Consistently, whenever we have information about a blocky behavior of the solution, in order to maximize the spatial resolution of the result, we consider an alternative regularization term called the minimum gradient support (MGS), which has been successfully applied in various geophysical settings [31,35,36]. This approach consists of substituting the term Lq 2 in (2.7) by the nonlinear stabilizing term [8] (2.8) where L is a regularization matrix. The nonlinear regularization term S τ (q) favors the sparsity of the solution and the reconstruction of blocky features. Indeed, it can be shown [28,30] that, when τ approaches 0, this approach minimizes the number of components where the vector Lq is nonzero. Therefore, if L is chosen to be the discretization of the first derivative D 1 , the stabilizer (2.8) selects the solution update corresponding to minimal nonvanishing spatial variation; see [8] for further comments and for a description of our implementation.
3. The FDEMtools software package. The package consists of a set of MATLAB routines which implements the algorithms for the inversion of FDEM data sketched in the previous sections. In addition to the analysis and solution routines, the package includes 1D and 2D test problems, showing how to process both synthetic and experimental data sets, and investigate different device configurations and numerical strategies, in order to compare them and draw conclusions.
The software is available from Netlib (http://www.netlib.org/numeralgo/) as the na53 package. It requires P. C. Hansen's Regularization Tools package [16] to be installed on the computer, and its directory to be added to MATLAB search path. FDEMtools is distributed as an archive file; by uncompressing it, a new directory will be created, containing the functions of the toolbox as well as a user manual. This directory must be added to MATLAB search path in order to be able to use the software from other directories. More information on the installation procedure can be found in the README.txt file, contained in the toolbox directory. All the routines are documented via the usual MATLAB help function. Table 3.1 groups the MATLAB routines by subject area, giving a short description of their purpose. The first group "Forward Model Routines" includes the functions for computing the forward model, that is, the model prediction for a given conductivity and permeability distribution. The user can optionally compute also the Jacobian matrix, by its analytical expression. The section "Computational Routines" describes the codes for the inversion algorithm, and the last section "Auxiliary Routines" lists some routines needed to complete the whole process, that are unlikely to be called directly by the user. See the file Contents.m for details.
For example, the command computes the function M V in (2.3) (corresponding to the "vertical" orientation of the device), as well as its Jacobian with respect to the σ variables, under the assumption that the conductivity and the permeability in each layer are defined, respectively, by the vectors sigma and mu, and the layers thickness by the vector d. The variables h, R, and omega define the height of the device above the ground, the inter-coil distance ρ, and the angular frequency ω, respectively; the last three variables may be vectors, in case of multiple measurements.
The computed predictions are returned in a vector, according to a particular ordering which is described in the documentation; see help hratio. The reflection factor (2.2) is computed by the function reflfact, while the computation of the Hankel transform in (2.3) is performed by a quadrature formula described in [1], which is constructed by two routines from [17]. The "Forward Model" functions are totally independent of the rest of the package, so they could be used to test the performance of other inversion algorithms.
The computational routines which are intended to be called directly by the user are emsolvenlsig and emsolvenlmu. The first one reconstructs the conductivity distribution with respect to depth, taking as input the device measurements, and assuming the permeability distribution is known. The second function does the same for the magnetic permeability.
The following example shows how to call the inversion function in order to perform data After initializing some variables, the model function is chosen. In this case, it is the complex signal in its entirety considering both orientations for the device, i.e., both components of (2.3). The permeability mu is assumed to be known and the conductivity sigma is determined by the regularized algorithm, assuming the data set is contained in the vector y. Some options are set before the call and encoded in the opts variable: the forward model, the regularization matrix (in this case the discrete approximation of the second derivative), the orientation of the coils (both vertical and horizontal), and the request to display intermediate results.
In the case of multiple measurements, taken at different locations, each data vector should be stored in a column of y. In this case the inversion routine processes a column at a time, and returns a matrix containing in each column the computed distribution of the electrical conductivity.
The package provides the driver program as an example of how to perform the actual computations. The first part of the driver defines the parameters which control the construction of the synthetic data set, the device setting, the soil discretization, the chosen model, and the tuning of the algorithm. Then, the synthetic data set is constructed, it is contaminated by Gaussian noise, and the inversion routine is called.
The driver program provided with the package generates a synthetic data set corre-  The graph on the left shows the solution selected by the "corner" criterion; the one on the right the data fitting between the synthetic data set and the model prediction. This two graphs can be updated by clicking with the left mouse button on one of the solutions displayed in Figure 3 sponding to a device operating with a single pair of coils, at 6 operating frequencies, and a single height above the ground. Both the vertical and the horizontal orientations are considered so that there are m = 12 data points, but since the complex signal is processed the number of data values to fit is doubled. The subsoil is discretized in n = 35 equal layers, up to a depth of 3.5m. We assume a constant magnetic permeability µ 0 (the value in the empty space) and reconstruct the distribution of the electrical conductivity.
During the computation, information about the data set and the intermediate results are displayed in the MATLAB main window and in some figures, which we reproduce here.  Figure 3.2 displays the solution selected by the "corner" criterion, and the one on the right the fitting between the input data set and the model prediction corresponding to the computed solu-tion. We remark that actual calls to driver may produce different results, because different realizations of the noise vector propagate differently in the computation.
A second test program driver2D is provided, which performs the processing of an experimental data set composed by 6 measurements. An additional figure displays the computed solutions in a 2D pseudo-color plot; see Figure 3.3.   Table 3.2 reports the options that can be set to tune the functionality of the package. All the options have a default value; see the driver program for an example of their use, as well as the help pages of the computational routines.
Further information on the use of the package, as well as on the data structure adopted to read experimental data sets, is provided in the User Manual, included in the package main directory. While for processing large experimental data sets it is probably preferable to use a script if > 0 show information on identification of reg. parameter rowcol size of grid in the iteration graphical window xlim x-limit for the graphs in the iteration window ylim y-limit for the graphs in the iteration window truesol exact solution, if available (to compute errors) resembling the driver program, the toolbox includes a MATLAB graphical user interface (GUI) in order to simplify experimenting with the functions of the package. It is started by issuing the command FDEMgui in the main MATLAB window; see Figure 3.4. It is composed of a set of input panels that lets the user generate or load data sets, choose different approaches and parameters for the inversion algorithm and, finally, visualize the computed results. The panels, which are extensively described in the User Manual, are the following: • Physical quantity to be inverted, • Data to be inverted, • Device configuration, • Data management, • Synthetic dataset, • Discretization, • Noise, • Inversion options, • Regularization. The GUI is completed by three buttons that allow the user to start the computation, interrupt it, in case something goes wrong, and save the computed solution to a .mat file.