Deidda.pdf 549 ANNALS OF GEOPHYSICS, VOL. 46, N. 3, June 2003 Inversion of electrical conductivity data with Tikhonov regularization approach: some considerations Gian Piero Deidda (1), Ernesto Bonomi (2) and Cristina Manzi (2) (1) Dipartimento di Ingegneria del Territorio, Università di Cagliari, Italy (2) Centro di Ricerca e Sviluppo di Studi Superiori della Sardegna, Cagliari, Italy Abstract 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 Key words inverse problems – Tikhonov regulariza- tion – projected conjugate gradient – high-frequency electromagnetics 1. 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). How- ever, in many situations the quantities we wish to determine are different from the ones we are able to measure. If the measured data de- pends, 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 best- known 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 con- sidered a simple least squares problem for the estimation of the electrical conductivity profile from high-frequency electromagnetic induction 550 Gian Piero Deidda, Ernesto Bonomi and Cristina Manzi 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, 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 conduct- ivity 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 con- ductivity profile from synthetic apparent con- ductivity curves solving a least squares prob- lem 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. 2. EM38 electrical conductivity data 2.1. 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 ¥ 103 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 mag- netic 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 (2.1) where µ π0 74 10= ⋅ − henry/m is the magnetic permeability of free space. The EM38 measures the quadrature com- ponent 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 N s sB = =δ µ ωσ0 2 Fig. 1. Instrument configurations: horizontal (left) and vertical dipoles (right). s s 551 Inversion of electrical conductivity data with Tikhonov regularization approach: some considerations the assumption of a homogeneous medium and NB <<1, the secondary magnetic field is a simple function of w, s and s (2.2) where w = 2p f (McNeill, 1980). It is worth- while mentioning that the field ratio, in the ap- proximation provided by eq. (2.2), is inde- pendent 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. 2.2. 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 (2.6) where s (z) is the conductivity at depth z. The sensitivity function fV(z) is described by eq. (2.4). Similarly, for the horizontal orientation, the apparent conductivity σa H is written as follows: (2.7) with fH(z) given by eq. (2.5). Collecting meas- urements of σa V and σa H recorded at different ele- vation, h1, h2, …, hN, 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. H H ss p ≅ ωµ σ0 2 4 σ ωµ a s ps H H =       4 0 2 φV z z z ( )= +( ) 4 4 12 3 2/ φ H z z z ( )= − +( ) 2 4 4 12 1 2/ σ φ σaV Vh z h z dz( )= +( ) ( ) ∞ ∫0 σ φ σa H Hh z h z dz( )= +( ) ( ) ∞ ∫0 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. 552 Gian Piero Deidda, Ernesto Bonomi and Cristina Manzi 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: mj = m0 for j = 1, 2, …, M. Let d T(s) = [sV a (h1), s V a (h2), ..., s V a (hN), s Ha(h1), s H a(h2), ..., s H a(hN)] 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). 3. Data inversion 3.1. The least squares problem From the (2N ¥ M)-linear system (2.9) the electrical conductivity profile can now be esti- Fig. 3. Stratified subsurface model Fig. 4a,b. a) Synthetic ground conductivity profile; b) apparent conductivity data computed from the ground conductivity profile in (a). a b d d d σ( )=       V H d Kσ σ( )= . K V H =       V z h dzi j V i z z j j , = +( ) − ∫ φ 1 H z h dzi j H i z z j j , = +( ) − ∫ φ 1 553 Inversion of electrical conductivity data with Tikhonov regularization approach: some considerations 3.2. 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 K dσ − 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 sn quantifies the regularity of s = sa, 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 = sa , 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 1≤ <κ κ( ˆ )A A( ). This behaviour is illustrated by the curves displayed in fig. 5. The simplest regularizing operator is L0 = I, where I denotes the identity matrix. Another form of control may be obtained through the implementation of L2, the second derivative operator. L2 enforces the smoothness of the conductivity profile while L0 controls its fluc- tuations. As illustrated in fig. 5, L0 provides a better conditioning to matrix Â. mated from the available experimental infor- mation d by solving the least squares prob- lem 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 = KTK and b = KTd. 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) 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 lM and l1 are, respectively the largest and smallest eigenvalues of A, plays the promi- nent role as a measure of the difficulty of computing s = A-1b 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.  = A + L Lα nT n εε ΚΚs s d( )= − 2 A bσ = A V V H H= +T T b V d H d H= +T V T . εα α αs Ks d L s( )= − + ≥ 2 2 0n , .  bσα = 554 Gian Piero Deidda, Ernesto Bonomi and Cristina Manzi As already mentioned, the larger the a the more important the diagonal dominance of  is; but increasing a causes the perturbed solution sa 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 cer- tainly the most used. Since the minimum of ε a is a linear combination of two terms, K dσα − and Lnσα , 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, sa tends to the solution s of the ori- ginal least squares problem. This implies the se- quence of points K d Lσ σα α−( ), 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 L0 (fig. 6a) and L2 (fig. 6b). Note that, as expected for large values of a , the norm Lnσα 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 corre- sponding to the most balanced perturbed solution of the inverse problem. Fig. 5. Effect of matrix regularization on the condition number: high values of a lower the condition number. Fig. 6a,b. L-curve obtained from data d in fig. 4b with regularization imposed by: a) L0, and b) L2 (16-byte arithmetic). b a Fig. 7. Inverse problem solutions for the input data in fig. 4b. Inversions with no regularization and with L0 regularization give the best solutions while inversion with L2 regularization gives a meaningful solution. 555 Inversion of electrical conductivity data with Tikhonov regularization approach: some considerations 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 L0, and a = 4 . 10-2 for L2. 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 K d Lσ σα α−( ), 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 sa which is physically meaningless (fig. 7). 3.3. The projected conjugate gradient algorithm When solving the least square problem (3.5), the solution representing the conductivity profile must satisfy the physical requirement sa 0. Under this condition, the conjugate gradient algorithm is no longer applicable because, even if we start inside the feasible set S = Œ{ }s R sM 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 16- byte 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 L2 regularization are spread over a large range, in the L0 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 ma- trix A. It is worthwhile mentioning that a care- ful 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. Fig. 8a,b. a) Eigenvalues of Â. The largest eigenvalue of each matrix, placed around 50, is not displayed. b) Con- vergence history of the conjugate gradient algorithm with different regularizations. 556 Gian Piero Deidda, Ernesto Bonomi and Cristina Manzi Although these values are the most suitable, the discrepancy between the resulting solutions, with respect to that corresponding to a = 0, is practi- cally unnoticeable. 4. 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 non- orthogonal to the error (Hansen, 1998). This signifies that conjugate gradient yields results close to the truncated singular value decom- position solution in which the truncation para- meter (an integer at which the singular values are deemed to be negligible) equals the num- ber of iterations for the convergence of the con- jugate 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 sa, 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 least- squares problem, S N= = ∈ℜ{ }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, Im A A = for all ( )= ∈ℜ{ }y v y v M, orthogonal (see Appendix) to the kernel of A defined as Ker A A( )= ={ }z z 0 . Note that each solution can be decomposed as follows: x = x0+ l z with x0 ŒS a particular solution, z ∈ ( )Ker A and l a real number. Hence S represents a translation of the kernel of the system. The conjugate gradient algorithm, starting from an initial guess in the subspace Im(A), will necessarily converge to the unique solution x A0 ∈ ( )S I Im . 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 = Ax - b , 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 so- lution x of the system is reached at the inter- section between the Im(A) and the set of all solutions S, x = AS 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 sub- space 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. 5. Conclusions In this article the simple inversion problem based on McNeill’s linear response model of Fig. 9. A schematic view of the evolution of the so- lution in the admissible subspace Im(A). 557 Inversion of electrical conductivity data with Tikhonov regularization approach: some considerations 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 algo- rithm 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 spite of the good performance of the pres- ent methodology on synthetic examples here presented, its development certainly needs additional computational experiments adopt- ing 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 νν∈ ( )Ker A and w A∈ ( )Im . Being Im(A) spanned by eigenvectors associated to non-zero eigenvalues, we can write the vector w as a linear combination of eigenvectors: w wi=∑α i ; if li are the 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). =     = ( )     =∑ ∑νν ΑΑ ΑΑννT T i i i T i i i α λ α λ w w 0 νν νν νν νν νν ΑΑ⋅ = = ( )=     =     =∑ ∑ ∑w w w w w T T i i T i i i i T i i iα α λ λ α λ REFERENCES BARKER, R.D. (1990): Investigation of groundwater salinity by geophysical methods, in Geotechnical and Environmental Geophysics, edited by S.H. WARD, Vol. 2, 201-211. BERTSEKAS, D.P. and J.N. TSITSIKLIS (1989): Parallel and Distributed Computation: Numerical Methods (Prentice Hall), pp. 718. BIRGIN, E.G., J.M. MARTINEZ and M. RAYDAN (1999): Non- monotone spectral projected gradient methods on convex set, SIAM J. Optim., 10, 1196-1211. BORCHERS, B., T. URAM and J.M.H. HENDRICKX (1997): Tikhonov regularization for determination of depth profiles of electrical conductivity using non-invasive electromagnetic induction measurements, Soil Sci. Soc. Am. J., 61, 1004-1009. CAMERON, D.R., E. DEJONG, D.W.L. READ and M. OOR- STEVELD (1981): Mapping salinity using resistivity and electromagnetic techniques, Can. J. Soil Sci., 61, 67-78. COOK, P.G., M.W. HUGHES, G.R. WALKER and G.B. ALLISON (1989): The calibration of frequency-domain electromagnetic induction meters and their possible use in recharge studies, J. Hydrol., 107, 251-265. 558 Gian Piero Deidda, Ernesto Bonomi and Cristina Manzi HANSEN, P.C. (1992): Analysis of discrete ill-posed problems by means of the L-curve, SIAM Rev., 34, 561-580. HANSEN, P.C. (1998): Rank-Defi cient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion, SIAM Monographs on Mathematical Modeling and Computation 4, pp. 263. HENDRICKX, J.M.H., B. BAERENDS, Z.I. RAZA, M. SADIQ and M.A. CHAUDRY (1992): Soil salinity assessment by electromagnetic induction of irrigated land, Soil Sci. Soc. Am. J., 56, 1933-1941. HILGENDORF, A. (1997): Linear and nonlinear models for inversion of electrical conductivity profiles in filed coil from EM38 measurements, M.Sc. Thesis, Institute of Mining and Technology, New Mexico. KACHANOSKI, R.G., E.G. GREGORICH and I.J. VANWESEN- BEECK (1988): Estimating spatial variations of soil water content using non-contacting electromagnetic inductive methods, Can. J. Soil Sci., 68, 715-722. KACHANOSKI, R.G., E. DEJONG and I.J. VANWESENBEECK (1990): Field scale patterns of soil water storage from non-contacting measurements of bulk electrical conductivity, Can. J. Soil Sci., 70, 537-541. MCNEILL, J.D. (1980): Electromagnetic terrain conductivity measurement at low induction numbers, Tech Note TN-6, Geonics Ltd., Ontario, Canada. NOBLE, B. (1969): Applied Linear Algebra (Prentice-Hall, Inc., Englewood Cliffs, New Jersey), pp. 523. RHOADES, J.D., P.J. SHOUSE, W.J. ALVES, N.A. MANTEGHI and S.M. LESCH (1990): Determining soil salinity from soil electrical conductivity using different models and estimates, Soil Sci. Soc. Am. J., 54, 46-54. SAUNDERS, W.R. and S.A. COX (1987): Use of an elec- tromagnetic induction technique in subsurface hydrocarbon investigation, in Proceedings of the fi rst National Outdoor Action Conferences on Aquifer Restoration, Ground Water Monitoring and Geophysical Methods, 585-599. SCANLON, B.R., J.G. PAINE and R.S. GOLDSMITH (1999): Evaluation of electromagnetic induction as a recon- naissance technique to characterize unsatured flow in an arid setting, Groundwater, 37, 296-304. SHEETS, K.R. and J.M.H. HENDRICKX (1995): Non-invasive soil water content measurement using electromagnetic induction, Water Resour. Res., 31, 2401-2409. VA L E N T I N E, R.M. and T. KW A D E R (1985): Terrain conductivity as a tool for delineating hydrocarbon plumes in a shallow aquifer- A case study, in Proceedings of the National Water Well Association Conference on Surface and Borehole Geophysical Methods in Ground Investigations, 52-63. VOIEVODINE, V. (1976): Algèbre Linéaire, Editions Mir, Moscou. WABHA, G. (1990): Spline Models for Observation Data, SIAM, Philadelphia, vol. 59.