Inversion of electrical conductivity data with Tikhonov regularization approach: some considerations

Electromagnetic induction measurements, which are generally used to determine lateral variations of apparent electrical conductivity, can provide quantitative estimates of the subsurface conductivity at different depths. Quantitative inference about the Earth’s interior from experimental data is, however, an ill-posed problem. Using the generalised McNeill’s theory for the EM38 ground conductivity meter, we generated synthetic apparent conductivity curves (input data vector) simulating measurements at different heights above the soil surface. The electrical conductivity profile (the Earth model) was then estimated solving a least squares problem with Tikhonov regularization optimised with a projected conjugate gradient algorithm. Although the Tikhonov approach improves the conditioning of the resulting linear system, profile reconstruction can be surprisingly far from the desired true one. On the contrary, the projected conjugate gradient provided the best solution without any explicit regularization (a = 0) of the objective function of the least squares problem. Also, if the initial guess belongs to the image of the system matrix, Im(A), we found that it provides a unique solution in the same subspace Im(A). Mailing address: Dr. Gian Piero Deidda, Dipartimento di Ingegneria del Territorio, Sezione di Geologia Applicata e Geofisica Applicata, Università di Cagliari, P.zza d’Armi 16, 09123 Cagliari, Italy; e-mail: gpdeidda@unica.it


Introduction
The goal of collecting geophysical data is to gain meaningful information about a given Earth property (for example, density, seismic velocity or conductivity of a geologic body).However, in many situations the quantities we wish to determine are different from the ones we are able to measure.If the measured data depends, in some way, on the quantities we want, then the data contains at least some informa-tion about those quantities.Thus, measured data can be used to predict the quantities we really want.Unfortunately, the prediction is not straightforward because an inverse problem must be solved.A typical feature of inverse problems is that they are ill-posed.A unique solution may not exist and small errors in the data may cause prohibitively large variations in the estimates of the quantity sought.To overcome these difficulties one has to regularize the original problem, that is the original problem has to be replaced by a nearby well-posed problem in order to obtain a stable solution.One of the bestknown and most used regularization methods is Tikhonov regularization.
The aim of this work is to illustrate how Tikhonov regularization may sometimes be the cause of misleading results, far from the expected solutions.For such a purpose we considered a simple least squares problem for the estimation of the electrical conductivity profile from high-frequency electromagnetic induction measurements.High-frequency electromagnetics have enjoyed increased popularity in near surface geophysics and nowadays are routinely used in environmental and hydrogeological studies.Shallow applications of electromagnetic methods have been carried out for salinity monitoring in agricultural land (Cameron et al., 1981;Rhoades et al., 1990;Hendrickx et al., 1992), detection of contaminants in soils and shallow aquifers (Valentine and Kwader, 1985;Saunders and Cox, 1987;Barker, 1990), soil water content measurements (Kachanoski et al., 1988(Kachanoski et al., , 1990;;Sheets and Hendrickx, 1995), and vadose zone characterisation (Cook et al., 1989;Scanlon et al., 1999).A popular method employed in these applications is the frequency-domain electromagnetic (FEM) induction technique (McNeill, 1980) which uses a ground conductivity meter to measure the apparent electrical conductivity of the subsurface.Although this technique is generally used for detection of lateral changes in the apparent electrical conductivity, it can also provide quantitative vertical variations: FEM measurements acquired at different heights above the soil surface can be used to predict the electrical conductivity at different depths.
Adopting the same approach as Borchers et al. (1997), we estimated the electrical conductivity profile from synthetic apparent conductivity curves solving a least squares problem with Tikhonov regularization optimised with a projected conjugate gradient algorithm.Although the Tikhonov approach improved the conditioning of the problem, the calculated solution was surprisingly far from the desired true one.In particular, we found that the choice of the point on the L-curve corresponding to a = 0 provided the best solution, that is no smoothing penalty term has to be added to the objective function.The conjugate gradient provided itself sufficient regularization to assure stable solution.

EM38 instrument
Various non-invasive devices can be used to acquire information from the subsurface through electromagnetic induction measurements.The instrument used in this study was the EM38, an electromagnetic induction sensor manufactured be Geonics Limited, Ontario, Canada.It consists of two coils on a lightweight bar 1 m long, which includes calibration controls and a digital readout of apparent electrical conductivity in mS/m.The EM38 instrument operates at a frequency of f = 14.6 kHz which corresponds to w = 91.7 ¥ 10 3 rad/s.The coil spacing s is equal to 1 m.The instrument can be held so that the two coils are either oriented horizontally or vertically with respect to the soil surface, as illustrated in fig. 1. Alternate current is sent through the transmitter coil; this generates a magnetic field H p that induces current to flow on the second (receiver) coil, which in turn generates a secondary magnetic field H s .
Defining the skin depth d as the depth at which the primary magnetic field has been attenuated to 1/e of its original strength, we can introduce the induction number N B , which is the ratio of the intercoil spacing s to the skin depth d.For a soil with uniform conductivity s, it can be shown that The EM38 measures the quadrature component of the ratio of the two magnetic fields.In general, the secondary magnetic field is a complicated function of the intercoil spacing s, the operating frequency f and the ground conductivity s.It can be shown that, under the assumption of a homogeneous medium and N B << 1, the secondary magnetic field is a simple function of w, s and s where w = 2p f (McNeill, 1980).It is worthwhile mentioning that the field ratio, in the approximation provided by eq.(2.2), is independent of the dipole orientation.
Since the ratio between the two magnetic fields is linearly proportional to the electrical conductivity of the soil, the conductivity is thus evaluated by measuring this ratio.Whatever the structure and composition of the medium under investigation, eq.(2.2) is used to define the apparent conductivity s a , the experimental information required when solving the inverse problem of electromagnetic sounding . (2.3) Under the assumption of homogeneity, and normalizing all spatial dimensions with respect to s, McNeill (1980) described the sensitivity f of the instrument to conductivity at depth z, for both vertical and horizontal modes (2.4) . (2.5) Figure 2 displays the sensitivity of the EM38 for both vertical and horizontal dipole configurations, versus the depth z measured in units of distance between the two coils.

The forward propagation model
Borchers et al. (1997) describe and discuss a more general linear model for the instrument response, which can be extrapolated from the model of McNeill (1980) under the following assumption: 1) The subsurface model represents a horizontally stratified medium in which the current flow is entirely horizontal.
2) The current flow at any point of the subsurface is independent of the current flow at any other point, since the magnetic coupling between all current loops is negligible.
With the two coils in vertical mode, assuming the instrument at a given height h above the soil, σ a V takes the form where s(z) is the conductivity at depth z.The sensitivity function f V (z) is described by eq.
(2.4).Similarly, for the horizontal orientation, the apparent conductivity σ a H is written as follows: (2.7) with f H (z) given by eq.(2.5).Collecting measurements of σ a V and σ a H recorded at different elevation, h 1 , h 2 , …, h N , above the soil surface, the two integral eqs.(2.6) and (2.7) provide the linear forward model to invert, from which the electrical conductivity profile s(z) can be estimated.Assuming a stratified medium model (fig.3) the subsurface is divided into M layers with specified thickness dz j , electrical conductivity s j and magnetic permeability m j equal, in this context, to that of the free space: T denote the vector gathering data relative to apparent conductivity measurements . (2.8) Using the instrument response model described by (2.6) and (2.7), the following system of linear equations establishes a correspondence between the subsurface conductivity profile and the apparent conductivity measurements: (2.9) The system matrix K is constructed as follows: (2.10) where the elements of V and H are, respectively (2.11) and (2.12) for i = 1, 2, …, N and j = 1, 2, …, M. Figure 4b shows an example of input data vector d(s) obtained from the model of the electrical conductivity profile s displayed in fig.4a, for N =11 and M = 30.Each value of the apparent conductivity was simulated for different elevations of the instrument above the soil, implementing the forward propagation model (2.9).

The least squares problem
From the (2N ¥ M)-linear system (2.9) the electrical conductivity profile can now be esti-

Regularization in the sense of Tikhonov
A way to enhance the stability of the inverse problem is Tikhonov regularization, a method which allows a form of optimal tuning on the sensitivity of the solution to input data errors.This is obtained by the trade-off between the residual norm and some desirable property resulting from the action of a discrete differential operator L n on the profile of s.Here n denotes the order of the differentiation.
A perturbed solution of the inverse problem is computed by solving the least squares problem associated to the new functional (3.5)The norm L s n quantifies the regularity of s = s a , the electrical conductivity profile that minimises (3.5).Clearly, for a = 0, ε 0 corresponds to the function ε , eq. (3.1), and thus the optimal con- ductivity is the solution of the system of eq.(3.2).In the general case of a > 0, the minimum of ε a is reached for the conductivity profile s = s a , solution of the linear system (3.6)where Â, a symmetric, positive definite matrix, has the form (3.7) and A is defined, as before, by eq.(3.3).
Since L n is always a diagonally dominant matrix, the perturbing term in (3.7) becomes absolutely crucial to improve the conditioning of Â with respect to that of A (Bertsekas and Tsitsiklis, 1989).Hence, for a > 0, we necessarily have This behaviour is illustrated by the curves displayed in fig. 5.
The simplest regularizing operator is L 0 = I, where I denotes the identity matrix.Another form of control may be obtained through the implementation of L 2 , the second derivative operator.L 2 enforces the smoothness of the conductivity profile while L 0 controls its fluctuations.As illustrated in fig.5, L 0 provides a better conditioning to matrix Â.Since the instrument response s a depends on the cumulative effect of all sudden changes of the subsurface conductivity profile s, the measured data field is weakly sensitive to the perturbations of the medium conductivity.Conversely, this physical property is mathematically translated into a strong sensitivity of the inverse problem solution s with respect to the perturbation of the apparent conductivity s a .
The condition number, defined as κ A ( ) = λ λ = M 1 where l M and l 1 are, respectively the largest and smallest eigenvalues of A, plays the prominent role as a measure of the difficulty of computing s = A -1 b in face of data uncertainty and roundoff errors.A classical result for non-singular operators (Voievodine, 1976) states that, for large values of κ A ( ), the system solution s might be highly perturbed even in the case of weak perturbations of both A and b, or one of them.In such a situation, the problem is said to be ill-conditioned.
Matrix A, defined in (3.3), may be extremely ill-conditioned with the condition number κ A ( ) of the order of hundred of thousands.Obviously, from the point of view of stability, the condition number should be as close as possible to one.
As already mentioned, the larger the a the more important the diagonal dominance of Â is; but increasing a causes the perturbed solution s a to deviate further from the exact solution s.

It is then important to achieve an acceptable balance between stability and accuracy of the solution by tuning carefully the regularization parameter a.
There are several heuristic ways to proceed in order to select a (Wabha, 1990;Hansen, 1992;Hilgendorf, 1997), but the criterion described below, based on the L-curve construction, is certainly the most used.Since the minimum of ε a is a linear combination of two terms, K d σ α − and L n σ α , the idea behind this criterion is to display one term as a function of the other for different values of the parameter a.The resulting plot is called the L-curve.
According to the Tikhonov theory, for a going to zero, s a tends to the solution s of the original least squares problem.This implies the sequence of points , n moves along a trajectory, denoted as L-curve, presenting a limit point.This point is indicated with a cross in both examples of fig.6a,b, where the regularization is imposed by L 0 (fig.6a) and L 2 (fig.6b).Note that, as expected for large values of a , the norm L n σ α decreases while the residual K d σ α − increases.
In the spirit of the L-curve criterion, the most suitable value of the regularizing parameter a is determined by selecting one intermediate point on the corner of the L-curve (Hansen, 1992).Such a point, indicated with a circle in fig.6a,b, is supposed to provide, in terms of accuracy and regularity, the value of the parameter corresponding to the most balanced perturbed solution of the inverse problem.   Figure 7 shows three solutions of the system (3.6).Two of them, marked with circles in fig.6a,b, correspond to a = 10 -4 for L 0 , and a = 4 .10 -2 for L 2 .The solution profile obtained with no regularization (fig.7) gives rise to the limit point framed by a square symbol on the L-curves.Note that on the , n plane the point corresponding to the true conductivity profile may be far away above the point selected by the L-curve criterion.Although this last point represents a compromise between accuracy and regularity of the perturbed solution, the resulting a may lead to a conductivity profile s a which is physically meaningless (fig.7).

The projected conjugate gradient algorithm
When solving the least square problem (3.5), the solution representing the conductivity profile must satisfy the physical requirement s a 0. Under this condition, the conjugate gradient algorithm is no longer applicable because, even if we start inside the feasible set S = OE { } s R s M 0 , an update may take solutions outside that set (non-physical solutions).The projected conjugate gradient generalized the conjugate gradient algorithm to the case where there are constraints (Bertsekas and Tsitsiklis, 1989;Birgin et al., 1999) assuring the attainment of the physical requirement of the non-negativity of the solution.
All the computations run for the analysis illustrated in figs.6a,b and 7 are performed in 16byte arithmetic.With this level of accuracy, the choice of the point on the L-curve corresponding to a = 0 provides the best solution, in the sense of proximity to the point representing the desired true solution.As a consequence, in this framework the regularization of system (3.2) is simply not necessary.
As illustrated in fig.5, the problem (3.2) is ill-conditioned.However, the computation of the eigenvalues of Â, all strictly positive, shows an intriguing result concerning their aggregation close to zero as displayed in fig.8a.While the eigenvalues relative to the L 2 regularization are spread over a large range, in the L 0 case and in the case with no regularization, namely Â = A, there are only three distinct eigenvalues and the remaining group is practically coincident near the smallest one.This property is of fundamental importance to obtain a fast convergence of the conjugate gradient algorithm.In particular, it can be shown that the method converges faster if most of the eigenvalues of the system matrix Â are clustered in a small interval and the remaining eigenvalues lie to the right of the interval (Bertsekas and Tsitsiklis, 1989).This result is eloquently illustrated by the three convergence curves presented in fig.8b.It also shows the excellent performance of the case a = 0, in spite of the ill-conditioning of matrix A. It is worthwhile mentioning that a careful construction of the L-curve may present values of a > 0 providing a better solution of the least squares problem, in the sense of proximity to the desired true solution, as shown in the example of fig.6a where a = 10 -12 .Although these values are the most suitable, the discrepancy between the resulting solutions, with respect to that corresponding to a = 0, is practically unnoticeable.

Discussion
We have seen that clustered eigenvalues allow faster convergence of the conjugate gradient.In fact, the conjugate gradient converges after a number of iterations equal to the number of distinct eigenvalues whose eigenvector is nonorthogonal to the error (Hansen, 1998).This signifies that conjugate gradient yields results close to the truncated singular value decomposition solution in which the truncation parameter (an integer at which the singular values are deemed to be negligible) equals the number of iterations for the convergence of the conjugate gradient itself.Therefore, the conjugate gradient can also be seen as a regularization tool operating on the original least-squares problem.It is important to note, however, that this is not actually the case of Tikhonov regularization since it is only used as a trick to improve the conditioning of the original problem, providing a solution s a , approximating s.
Another aspect that needs discussion refers to the uniqueness of the solution.The synthetic example here presented belongs to the class of the underdetermined system, the number of unknowns, M, being larger than the number of equations, N.This kind of problem does not admit a unique solution (Noble, 1969).However, the strategy implemented with the conjugate gradient algorithm is able to select, in the set of all possible solutions of the leastsquares problem, x Ax b b , with , a unique solution which is in a suitable subspace of unknowns of dimension N, where N is the rank of the system matrix A. More precisely this subspace is the image of the matrix A, The conjugate gradient algorithm, starting from an initial guess in the subspace Im(A), will necessarily converge to the unique solution

orthogonal (see Appendix) to the kernel of A defined as
. As a matter of fact, the iterative procedure of the conjugate gradient performs only transformations by multiplying matrix A to a residual vector.As the gradient, defined by g = Axb , is necessarily in Im(A), each iteration of the conjugate gradient will keep the search for the solution within the subspace Im(A).
Obviously, b = Ax is in Im(A).
As schematically depicted in fig.9, the solution x of the system is reached at the intersection between the Im(A) and the set of all solutions S, x = A S I Im( ).It is the orthogonal projection of all solutions on the image of the matrix of the system Ax = b.This solution is unique and it is the best possible within the subspace Im(A) being at the intersection with the subset S of all possible solutions of the least-squares problem.This analysis does not contemplate constraints on the solution.The introduction of projections in the conjugate gradient confines the search within a convex set of Im(A) and results in highly non-linear iterates converging to the best possible profile x in the sense of proximity to that of the unconstrained problem x. the EM38 ground conductivity meter has been formulated and developed.Computer experiments on synthetic data provide credibility to results and conclusions presented in this work.Computations have shown how Tikhonov regularization and in particular the L-curve criterion can be the cause of misleading conductivity profiles far from the desired true ones.The analysis of the system matrices to invert their structure and eigenvalues shows that, although the original system problem is ill-conditioned, the conjugate gradient algorithm is a robust method for the solution of the electrical conductivity data inversion without any additional regularization of the least squares problem.The projection strategy, implemented together with the conjugate gradient, enforces the positivity of the solution and provides the best possible profile in the sense of proximity to that of the unconstrained problem.

In this article the simple inversion problem based on McNeill's linear response model of
In spite of the good performance of the present methodology on synthetic examples here presented, its development certainly needs additional computational experiments adopting more adequate models of two loops EM soundings over a stratified half-space.For a more complete validation, it will be necessary to undertake a more detailed analysis on real data acquired in sites where the stratification of the near subsurface and the relative electric properties are well known.

Appendix.
To observe the orthogonality of the subspaces Im(A) and ker (A), let us consider the scalar product v .w with ν ν ∈ ( ) /m is the magnetic permeability of free space.

Fig. 2 .
Fig. 2. Instrument sensitivity curves for the two dipole configurations.Horizontal mode configuration is more sensitive to contributions from materials at the very near subsurface while the vertical mode configuration better discriminates contributions at lower depth, with a maximum value at about 0.4 times the distance between the coils.

Fig. 3 .
Fig. 3. Stratified subsurface model mated from the available experimental information d by solving the least squares problem associated with the functional (3.1) at the moment without any requirement, physical or mathematical, on the solution profile.The minimum of the function ε (s) is reached for an electrical conductivity profile s = s, solution of the following system: (3.2) where A = K T K and b = K T d.From eq. (2.10) it follows that: (3.3)where A, an (M ¥ M)-matrix, is symmetric and positive definite.The right-hand side of (3.2), see eq. (2.8), takes the form (3.4)

Fig. 5 .
Fig. 5. Effect of matrix regularization on the condition number: high values of a lower the condition number.

Fig. 7 .
Fig. 7. Inverse problem solutions for the input data in fig.4b.Inversions with no regularization and with L 0 regularization give the best solutions while inversion with L 2 regularization gives a meaningful solution.

Fig
Fig. 8a,b.a) Eigenvalues of Â.The largest eigenvalue of each matrix, placed around 50, is not displayed.b) Convergence history of the conjugate gradient algorithm with different regularizations.
that each solution can be decomposed as follows: x = x 0 +lz with x 0 OES a particular solution, z ∈ ( )Ker A and l a real number.Hence S represents a translation of the kernel of the system.

Fig. 9 .
Fig. 9.A schematic view of the evolution of the solution in the admissible subspace Im(A).
Im(A)  spanned by eigenvectors associated to non-zero eigenvalues, we can write the vector w as a linear combination of eigenvectors:w w i = ∑ α i ; if l i arethe eigenvalues respectively associated to the eigenvectors w i (that is Aw i =l i w i ), and in the hypothesis of symmetric matrix A, for each v and for each w we have This means that Im(A) is orthogonal to Ker(A).