Subscribe Now Subscribe Today
Research Article

A Damped Least-Squares Inversion Program for the Interpretation of Schlumberger Sounding Curves

Yunus Levent Ekinci and Alper Demirci
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail

An iterative inversion program for 1D interpretation of Schlumberger sounding curves is presented in this research. The program, written in MATLAB v.5.3 (R11), optimizes the fit between the observed and the calculated apparent resistivity data. The procedure is initialized by user-selected initial model parameters and employs the damped least-squares solution with the Singular Value Decomposition (SVD). The applicability of the present inversion program was tested on three synthetic models which symbolize ideal resistivity distributions of geological relevance. Noise was added to two models as a further concession to reality. Moreover, the program was used to interpret real data from Turkey and the results obtained were compared with published results. In all cases, applications indicated that the model parameters obtained are in good agreement with the real ones and the published results. The program is suitable for shallow small scale investigations if the lateral changes in the study area are expected to be smooth.

Related Articles in ASCI
Similar Articles in this Journal
Search in Google Scholar
View Citation
Report Citation

  How to cite this article:

Yunus Levent Ekinci and Alper Demirci, 2008. A Damped Least-Squares Inversion Program for the Interpretation of Schlumberger Sounding Curves. Journal of Applied Sciences, 8: 4070-4078.

DOI: 10.3923/jas.2008.4070.4078



The Direct Current (DC) resistivity sounding method is one of the oldest geophysical techniques which is intensively used for the evaluation of the deep and shallow structure of the subsurface. By introducing the electrical current directly into the ground through a pair of electrodes (current electrodes), the difference of the resulting voltage can be measured between the other pair of electrodes (potential electrodes). The apparent resistivity of the subsurface which is the reciprocal of conductivity can be calculated by this way. The method aims to determine the resistivity variation with depth. The depth of the penetration depends on the distance between the current electrodes. Increasing the depth of the penetration can be carried out by enlarging the distance between the current electrodes from a small distance in the beginning to larger distances at the end of the array.

The investigations performed by using DC resistivity sounding method can yield satisfactory results in the areas where the ground is horizontally layered with a little lateral variation. The method has been applied to a wide variety of geophysical exploration problems so far. The most successful and beneficial surveys have been carried out in hydrogeological, geothermal, environmental and engineering applications (e.g., Batayneh et al., 2001; Mota et al., 2004; Agnesi et al., 2005; Özurlan and Şahin, 2006; Asfahani and Radwan, 2007; Park et al., 2007).

The apparent resistivity data do not represent the true resistivity distribution of the subsurface. Thus, to determine the true resistivity distribution of the subsurface, inversion of the observed apparent resistivity data must be carried out. DC resistivity inverse problem was investigated for the first time in the early 1930s. From that time until the late 1980s, the methodology of field survey and character of data from the measurements remained unchanged. From the late 1980s and early 1990s until today, there has been considerable improvement in data collection and interpretation (Muiuane and Pedersen, 1999). In the last decade, depending on high technological development of computer controlled multi-electrode survey systems and the inversion software packages, much of the investigations carried out on DC resistivity method have mainly focused on 2D and 3D data acquisition and interpretation techniques (e.g., Nassir et al., 2000; Dahlin, 2001; Bauer et al., 2006; Drahor et al., 2007; Ekinci et al., 2008). However, 1D techniques can be used to interpret the data and may yield good approximation of the earth`s subsurface if the investigation is carried out in areas with gently sloping or flat lying ground and where the resistivity and thickness variations between the adjacent layers are smooth.

In the frame of this study we present a MATLAB-based inversion program that uses damped least-squares solution for the interpretation of Schlumberger sounding curves. The reliability of the present program was tested on both noisy and noise-free synthetic data. Additionally, a field survey was also conducted in Gökçeada-Turkey to scrutinize the efficiency of the program on real data. We compared the results with our earlier study which indicates almost horizontally layered resistivity variations obtained by the interpretation of 2D resistivity inverse model sections (Ekinci et al., 2008).


Schlumberger electrode configuration is very simple to use in the survey area and hence is the most popular DC resistivity sounding array. The relationship between the apparent resistivity data and the layer parameters is expressed by an integral equation. The very well-known expression of the computation of a forward response for Schlumberger electrode configuration (ρa) over an earth model consisting of homogeneous and isotropic layers is related to the kernel function through the Hankel integral and can be formulated as follows (Koefoed, 1970):


where, s is half the current electrode spacing in Schlumberger electrode configuration, J1 denotes the first-order Bessel function of the first kind and λ denotes the integral variable. The resistivity transform function, T(λ), is given by the recurrence relationship as follows:


where, n denotes the number of layers, ρi and hi are the resistivity and thickness of the ith layer, respectively. Equation 2 is used to compute the forward model response for DC resistivity sounding data.

Non-linear least-squares inversion scheme iteratively updates the model parameters in each step with the use of a correction vector which is the solution of a set of normal equations. Inversion of geoelectrical data is an ill-posed problem. Therefore, small changes in the data may lead large changes in the model. This and the ensuing sub-optimality restrict the initial model to being in the near vicinity of the true model. The problem may be reduced by introducing damping into the system of equations (Roy, 1999). This parameter is added to the principal diagonal of ATA which helps to increase the DC level of eigenvalues such that none of its eigenvalues can become zero (Raju, 2003). The damped least-squares solution is given by the Eq. 3 as follows:


where, Δp is the parameter correction vector, Δd is the data difference vector, A is the Jacobian matrix. I is the identity matrix and the term ε is called damping factor (Levenberg, 1944; Marquardt, 1963).

Singular Value Decomposition (SVD) is well-known technique used in many areas of applied sciences including the earth sciences. It can be easily applied to small scale geophysical problems. It is mathematically robust and numerically stable and also provides other vital information on the state of the model and data thus enabling model resolution and covariance studies (Meju, 1994). However it must be noted that the use of SVD is not logical for large scale problems. The large scale problems can be solved either explicitly or implicitly using iterative methods like conjugate gradients. In this study, because the examples used are small scale problems, Eq. 3 was solved using SVD in the inversion scheme.

An nxn or nxp matrix A can be factored into a product of three other matrices as follows:


where, for n data and p parameters, U (nxp) and V (pxp) are, respectively the data space and parameter space eigenvectors and S is a pxp diagonal matrix containing at most r non-zero eigenvalues of A, with r≤p. These diagonal entries in S (λ1, λ2, ..., λp) are called singular values of A. This factorization is known as SVD of A (Meju, 1994). If SVD is used in damped least-squares solution, Eq. 5 is obtained:


By adding the damping factor to the diagonal elements we obtain:


The inverse of the Eq. 6 is given by as follows:


Substituting the Eq. 7 in 5 we obtain:


and the parameter correction vector can be expressed as:


Equation 9 provides damped least-squares solution with SVD. For automated inversion, the common practice is to set the damping factor first to a large positive value thus taking the advantage of the good initial convergence properties of the steepest descent method and thereafter the damping factor is multiplied by a factor less that unity after each iteration so that the least-squares method predominates near the solution (Meju, 1994). A way of determining the damping factor has been developed by Arnason and Hersir (1988) and given by as follows:


where, L is the test number performed for the damping factor at any iteration, λ is the parameter eigenvalue and the term of Δx is given by:


where, xr-1 is the misfit value obtained previous iteration and xr is the misfit value of the current iteration. In this study, Eq. 10 and 11 were used to set the damping factor in each iteration.


The inversion program presented in this paper is composed of three parts. One is the main program and the others are the functions.

The main program VES1dinv reads the observed apparent resistivity data and AB/2 spacings from the data file, uses user-selected initial model parameters (number of layers, resistivity and the thicknesses of these layers) and inverts the apparent resistivity data by using damped least-squares solution with SVD. The program also uses the basic plotting commands for the illustration of the inversion results.
The function Jacobian calculates the partial derivatives of parameters numerically by using a user-selected difference interval value and forms the Jacobian matrix.
The function VES1dmod calculates the apparent resistivity data from the given values (resistivity, thickness and AB/2) for Schlumberger electrode configuration by using the filter coefficients (Nyman and Landisman, 1977).

The full codes of the main program and the functions are given in the Appendix A-C, respectively.


Synthetic models with ideal resistivity distributions of geological relevance were used in the applications. Three different examples consisting of three- or four-layered earth models were preferred. In general, as a similar explanation mentioned by Muiuane and Pedersen (1999), the subsurface models where resistivity increases or decreases with depth are not difficult to resolve. On the other hand, if a resistor layer is embedded between two conductor layers or a conductor layer is embedded between two resistor layers convergence is more difficult and slow in comparison to first type. We used these types of models for the numerical examples. The synthetic data were generated at a predetermined set of 20 AB/2 values for all models by using the forward modeling function of the program. The 5% noise level is reasonable for sites with high resistivity contrasts and higher electrode contact resistances (Dahlin, 1993). Thus, as a further concession to reality, Model 2 and 3 were contaminated with ±5% random noise depending on the magnitude of each data point.

First example demonstrates a noise-free three-layered earth model with decreasing resistivity with depth. Second example deals with a three-layered earth model consisting of a resistor layer between two conductors. A four-layered, decreasing resistivity with depth model was used in the last example. We used homogeneous half spaces for the initial models. The inversion procedure was controlled by two misfit functions in the iterations as given in Appendix A. When the predetermined error limit was reached, the inversion procedure was stopped automatically. The convergences were calculated as an rms error function between the observed and calculated apparent resistivity values at the end of the process as follows:


where, obs is observed, cal is calculated apparent resistivity values and ND is the number of the apparent resistivity data. In all cases, although the initial guesses were far removed from the true ones, the model parameters were resolved and the program yielded satisfactory results. The inversion results are illustrated as best-fit curves and resistivity-depth distributions for Model 1, 2 and 3 in Fig. 1-3, respectively.

Fig. 1: Model 1, (a) observed and calculated apparent resistivity curves and (b) resistivity-depth distribution of the initial and final models

Fig. 2: Model 2, (a) observed and calculated apparent resistivity curves and (b) resistivity-depth distribution of the initial and final models

Fig. 3: Model 3, (a) observed and calculated apparent resistivity curves and (b) resistivity-depth distribution of the initial and final models

Table 1: Observed and calculated apparent resistivity data of the models

Table 2: True and obtained model parameters of the models

Observed and calculated apparent resistivity values and obtained model parameters after inversion procedure are given in Table 1-2, respectively.


To test the validity of the present inversion program on real data, a Schlumberger sounding curve was interpreted. The survey area, named Kaleköy is located in the northeastern part of Gökçeada (Fig. 4). The area has been previously studied (Ekinci et al., 2008). The seawater and the freshwater saturated layer were determined by the interpretation of 2D electrical resistivity tomography survey. The resistivity model section shows almost horizontally layered resistivity variations (Ekinci et al., 2008). Therefore, we expected that Schlumberger sounding should provide similar results.

The geological settings of Gökçeada are represented by a thick sedimentary sequence, volcanic rock and alluvium. Regarding of the survey site, western and eastern parts of it in which Mezardere Formation outcrops are hilly and the formation consists of sandstone-shale sequence. The plain is also represented by alluvium (Fig. 4). Hill slope allows rapid surface runoff wherefore infiltration occurs easily into alluvium. The existence of water saturation in this area is suggested due to both hydrogeological conditions (shallow water level, rapid surface runoff) and the high P-wave velocity in the alluvium (Ekinci et al., 2008). The site is also appropriate to investigate the seawater existence. Thus, we applied a Schlumberger sounding in order to test the efficiency of the developed program on the determination of seawater existence and its depth and thickness. The measurements were performed with Iris Syscal-R1+ resistivity meter. The maximum half-current electrode separation was 120 m. In order to improve the quality of the data, the resistivity meter performed a stack of a minimum of three for each data point.

Fig. 4: The location map of site and the Schlumberger sounding (not to scale)

Fig. 5: (a) observed and calculated apparent resistivity curves and (b) resistivity-depth distribution of the initial and final models for the field example

Table 3: The obtained model parameters of the real data

When the standard deviation between the measurements reached 3%, a maximum of five measurements were performed for each data point during the investigation. The standard deviations were generally below 1%.

We used a four-layered earth model for the initial guess. The inversion converged after the 65 th iteration with an rms error of 0.77. The fit between the observed and calculated apparent resistivity data are shown in Fig. 5 and obtained model parameters are given in Table 3. The sounding curve represents decreasing resistivity with depth but ends in a resistive basal. According to the geological evaluation of the inversion results, the top layer represents the consolidated alluvium of ~2.5 m overlying a water saturated zone, which extends to a depth of ~13.5 m. This water saturated zone exists above the very conductive layer (< 2.5 Ωm). The source of the freshwater is thought to be the groundwater recharge process, which occurs directly from rainfall and from rapid surface runoff wherefore infiltration occurs easily into alluvium due to the hill slopes. The very conductive layer which extends to a depth of ~23 m corresponds to seawater. Additionally, the increase of resistivity values below that layer may trace the existence of Mezardere formation. This may be interpreted to be decreasing porosity in response to increasing compaction and this compaction may indicate Mezardere formation represented by sandstone-shale sequence. These results are in close agreement with those obtained by Ekinci et al. (2008).


A MATLAB-based 1D inversion program for the interpretation of the Schlumberger sounding curves was presented with synthetic and real data examples in this paper. The presented inversion scheme is computationally simple, fast and easy to carry out. The program uses damped least-squares solution with SVD and runs iteratively by updating the model parameters.

The program was tested on small scale synthetically produced noise-free and noisy apparent resistivity data sets which have different curve types. Although, the initial guesses were far removed from the true ones, the present inversion program produced satisfactory results. Additionally, the results obtained by the inversion of the real data are in close agreement with the results of an earlier study. The existence of seawater and its depth and thickness were easily determined.

The program is suitable to use by the interpreter for small scale problems where the field conditions do not allow gathering enough data for multi-dimensional modeling scheme. It has been concluded that the presented program can provide reliable results where the resistivity and thickness variations between the adjacent layers are smooth. This inversion program is developed for Schlumberger electrode configuration but due to the fact that it is an open source code it can be easily modified for the other electrode configurations by changing the forward modeling function or extended to the needs of the other users.


The authors would like to express their gratitude to Dr. E.U. Ulugergerli for his useful suggestions. Thanks to Dr. T. Bekler for his permission to use the data set. They also thank to C. Ertekin, scientist, for his valuable information about the geological and hydrogeological properties of the surveyed area. The following undergraduate students are thanked for their aid in the acquisition of the field data: H. Şahlan, M. Özben, S. Kutluk, S. Güneş and O. Zop. The authors are especially grateful to Dr. M.A. Riahi for his critiques that have helped us to improve the paper and also thanks to three anonymous reviewers for their valuable comments.


% VES1dinv %
% 1D Inversion of Schlumberger Sounding Data %
% By Yunus Levent Ekinci and Alper Demirci %
close all;
clear all;
format long;
load ab2.dat;
x = ab2(:,1);
roa = ab2(:,2);
r = [ρ1 ρ2 ρn];
t = [t1 tn-1];
m = [r t];
rinitial = r;
tinitial = t;
lr = length(r);
lt = length(t);
kr = 10e-10;
iteration = 1;
maxiteration = 500;
dfit = 1;
while iteration<maxiteration
r = m(1:lr);
t = m(1+lr:lr+lt);
for i = 1:length(x)
s = ab2(i);
[g] = VES1dmod (r,t,s);
roa1(i,:) = g;
e1 = [log(roa)-log(roa1)];
dd = e1;
misfit1 = e1`*e1;
if misfit1<kr
axis([1 1000 1 1000])
xlabel(`AB/2 (m)`);
ylabel(`Apparent Resistivity (Ohm-m)`);
[A] = jacobian(ab2,x,r,t,lr,lt,roa,roa1);
[U S V] = svd(A,0);
ss = length(S);
say = 1;
k = 0;
while say<ss
diagS = diag(S);
beta = S(say)*(dfit^(1/say));
if beta<10e-5
beta = 0.001*say;
for i4 = 1:ss
SS(i4,i4) = S(i4,i4)/(S(i4,i4)^2+beta);
dmg = V*SS*U`*dd;
mg = exp(log(m)+dmg`);
r = mg(1:lr);
t = mg(1+lr:lr+lt);
for i5 = 1:length(x)
s = ab2(i5);
[g] = VES1dmod (r,t,s);
roa4(i5,:) = g;
e2 = [log(roa)-log(roa4)];
misfit2 = e2`*e2;
if misfit2>misfit1
(`Beta control`)
say = say+1;
k = k+1;
if k = = ss-1
iteration = maxiteration;
say = ss+1;
say = ss+1;
m = mg;
dfit = (misfit1-misfit2)/misfit1;
iteration = iteration+1;
a = iteration;
if dfit<kr
iteration = maxiteration;
say = say+1;
axis([1 1000 1 1000])
xlabel(`AB/2 (m)`);
ylabel(`Apparent Resistivity (Ohm-m)`);
observed = roa;
calculated = roa4;
h = legend(`obs`,`cal`,1);
format bank;
rr = [0,r];
tt = [0,cumsum(t),max(t)*10];
rrr = [0,rinitial];
ttt = [0,cumsum(tinitial),max(t)*10];
hold on;
ylim([0 50]);
xlim([0 1000]);
xlabel(`Resistivity (Ohm-m)`);
ylabel(`Depth (m)`);
rms = norm(roa4-roa) /sqrt(length(roa));


% Function Jacobian %
function[A] = jacobian(ab2,x,r,t,lr,lt,roa,roa1)
par = 0.1;
r2 = r;
for i2 = 1:lr
r2(i2) = (r(i2)*par)+r(i2);
for ii = 1:length(x)
s = ab2(ii);
[g] = VES1dmod (r2,t,s);
roa2(ii,:) = g;
A1(:,i2) = [(roa2-roa1)/(r(i2)*par)]*r(i2)./roa;
r2 = r;
t2 = t;
for i3 = 1:lt
t2(i3) = (t(i3)*par)+t(i3);
for ii = 1:length(x)
s = ab2(ii);
[g] = VES1dmod (r,t2,s);
roa3(ii,:) = g;
A2(:,i3) = [(roa3-roa1)/(t(i3)*par)]*t(i3)./roa;
t2 = t;
A = [A1 A2];


% Function VES1dmod %
function[g] = VES1dmod(r,t,s)
q = 13;
f = 10;
m = 4.438;
x = 0;
e = exp(0.5*log(10)/m);
h = 2*q-2;
u = s*exp(-f*log(10)/m-x);
l = length(r);
n = 1;
for i = 1:n+h
w = l;
v = r(l);
while w>1
w = w-1;
aa = tanh(t(w)/u);
v = (v+r(w)*aa)/(1+v*aa/r(w));
a(i) = v;
u = u*e;
i = 1;
g = 105*a(i)-262*a(i+2)+416*a(i+4)-746*a(i+6)+1605*a(i+8);
g = g-4390*a(i+10)+13396*a(i+12)-27841*a(i+14);
g = g+16448*a(i+16)+8183*a(i+18)+2525*a(i+20);
g = (g+336*a(i+22)+225*a(i+24))/10000;
1:  Agnesi, V., M. Camarda, C. Conoscenti, C. Di Maggio, I.S. Diliberto, P. Madonia and E. Rotigliano, 2005. A multidisciplinary approach to the evaluation of the mechanism that triggered the Cerda landslide (Sicily, Italy). Geomorphology, 65: 101-116.
CrossRef  |  Direct Link  |  

2:  Asfahani, J. and Y. Radwan, 2007. Tectonic evolution and hydrogeological characteristic of the Khanaser Valley, Northern Syria, derived from the interpretation of vertical electrical sounding. P. Applied Geophys., 164: 2291-2311.
CrossRef  |  Direct Link  |  

3:  Batayneh, A.T., A.S. Al Zoubi and A.A. Abueladas, 2001. Geophysical investigations for the location of a proposed dam in Al Bishriyya (Al Aritayn) area, Northeast Badia of Jordan. Environ. Geol., 40: 918-922.
CrossRef  |  Direct Link  |  

4:  Bauer, P., R. Supper, S. Zimmermann and W. Kinzelbach, 2006. Geoelectrical imaging of groundwater salinization in the Okavango Delta, Botswana. J. Applied Geophys., 60: 126-141.
CrossRef  |  Direct Link  |  

5:  Dahlin, T., 1993. On the automation of 2D resistivity surveying for engineering and environmental applications. Ph.D Thesis, Lund Institut of Technology.

6:  Dahlin, T., 2001. The development of DC resistivity imaging techniques. Comput. Geosci., 27: 1019-1029.
CrossRef  |  Direct Link  |  

7:  Drahor, M.G., G. Gokturkler, M.A. Berge, T.O. Kutulmq and N. Tuna, 2007. 3D resistivity imaging from an archaeological site in South-Western Anatolia: A case study. Near Surf. Geophys., 5: 195-201.
Direct Link  |  

8:  Ekinci, Y.L., A. Demirci and C. Ertekin, 2008. Delineation of the seawater-freshwater interface from the coastal alluvium of Kalekoy-Gokceada, NW Turkey. J. Applied Sci., 8: 1977-1981.
CrossRef  |  Direct Link  |  

9:  Koefoed, O., 1970. A fast method for determining the layer distribution from the raised kernel function in geoelectrical soundings. Geophys. Prospect., 18: 564-570.
CrossRef  |  Direct Link  |  

10:  Levenberg, K., 1944. A method for the solution of certain non-linear problems in least-squares. Q. Applied Math., 2: 164-168.
Direct Link  |  

11:  Marquardt, D.W., 1963. An algorithm for least-squares estimation of non linear parameters. J. Soc. Ind. Applied Math., 11: 431-441.
Direct Link  |  

12:  Nyman, D.C. and M. Landisman, 1977. VES dipole-dipole filter coefficients. Geophysics, 42: 1037-1044.
CrossRef  |  Direct Link  |  

13:  Meju, M.A., 1994. Geophysical Data Analysis: Understanding Inverse Problem Theory and Practice. Society of Exploration Geophysics Course Notes Series, No. 6, 1st Edn., SEG Publishers, Tulsa, Oklahoma, ISBN: 1-56080-027-5, pp: 296.

14:  Mota, R., F.A.M. Santos, A. Mateus, F.Q. Marques, M.A. Goncalves, J. Figueiras and H. Amaral, 2004. Granite fracturing and incipient pollution beneath a recent landfill facility as detected by geoelectrical surveys. J. Applied Geophys., 57: 11-22.
CrossRef  |  Direct Link  |  

15:  Muiuane, E.A. and L.B. Pedersen, 1999. Automatic 1D interpretation of DC resistivity sounding data. J. Applied Geophys., 42: 35-45.
CrossRef  |  Direct Link  |  

16:  Nassir, S.S.A., M.H. Loke, C.Y. Lee and M.N.M. Nawawi, 2000. Salt-water intrusion mapping by geoelectrical imaging surveys. Geophys. Prospect., 48: 647-661.
CrossRef  |  Direct Link  |  

17:  Ozurlan, G. and M.H. Sahin, 2006. Integrated geophysical investigations in the Hisar geothermal field, Demirci, Western Turkey. Geothermics, 35: 110-122.
CrossRef  |  Direct Link  |  

18:  Park, Y.H., S.J. Doh and S.T. Yun, 2007. Geoelectric resistivity sounding of riverside alluvial aquifer in an agricultural area at Buyeo, Geum River watershed, Korea: An application to groundwater contamination study. Environ. Geol., 53: 849-859.
CrossRef  |  Direct Link  |  

19:  Raju, D.C.V., 2003. LIMAT: A computer program for least-squares inversion of magnetic anomalies over long tabular bodies. Comput. Geosci., 29: 91-98.
CrossRef  |  Direct Link  |  

20:  Roy, I.G., 1999. An efficient non-linear least-squares 1D inversion scheme for resistivity and IP sounding data. Geophys. Prospect., 47: 527-550.
CrossRef  |  Direct Link  |  

21:  Arnason, K. and G.P. Hersir, 1988. One dimensional inversion of Schlumberger resistivity soundings: Computer Program, Description and User’s Guide: The United Nations University, Geothermal Training, Report 8, pp: 59.

©  2021 Science Alert. All Rights Reserved