Modified Occam's Inversion

ARTICLE IN PRESS

Computers & Geosciences 31 (2005) 319–328 www.elsevier.com/locate/cageo

An efficient 1D OCCAM’S inversion algorithm using analytically computed first- and second-order derivatives for DC resistivity soundings$ Nimisha Vedanti, Ravi P. Srivastava, John Sagode1, V.P. Dimri National Geophysical Research Institute, Council of Scientific and Industrial Research, Post Bag No. 724, Uppal Road, Hyderabad 500 007, India

Abstract An efficient algorithm has been developed for 1D resistivity inversion problem using both first- and second-order derivatives, which are computed analytically. The second-order derivative matrix, which is not used in the OCCAM’s inversion, has been incorporated into the algorithm employing analytical expressions. Computation of complicated second-order derivatives in each iteration is circumvented by a new algorithm. These modifications result in stable convergence of the OCCAM’s inversion and in general, better misfit can be achieved specially for smoothing parameter, mo1: The modified inversion algorithm, coded in MATLAB was tested using two synthetic Schlumberger resistivity sounding examples. Its application has been illustrated with field data from south India. r 2004 Elsevier Ltd. All rights reserved. Keywords: Hessian matrix; Jacobian matrix; Schlumberger sounding; Newton’s method; Smoothing parameter

1. Introduction The OCCAM’s inversion algorithm was first introduced by Constable et al. (1987) to find the smoothest model that fits the magnetotelluric (MT) and Schlumberger geoelectric sounding data. The method gained popularity in inversion studies and was applied to many investigations (LaBrecque et al., 1996; Siripunvaraporn and Egbert, 1996; Qian et al., 1997). In this scheme a highly nonlinear problem is formulated in a linear fashion, which obviates the computation of second-

$

Code on server at http://www.iamg.org/CGEditor/index.htm.

Corresponding author. Tel.: +91 40 234 34 600;

fax: +91 40 234 34651. E-mail address: [email protected] (V.P. Dimri). 1 CSIR-TWAS Fellow.

order derivatives that carry useful curvature information of the objective function. In this paper the 1D OCCAM’s algorithm has been improved by inclusion of second-order derivative matrix known as Hessian that is computed analytically. This leads to a quadratic equation approximation of the objective function. The modified algorithm has been tested on synthetic and real field resistivity sounding data. It is found that the modified algorithm is more stable and convergent than OCCAM’s inversion. The computation of second-order derivatives in Schlumberger resistivity sounding involves cumbersome piece of algebra and therefore these derivatives are computed numerically using finite difference schemes. This introduces many unacceptable errors and requires more computational time, which results in inaccurate curvature information that decides the step of descent where as computation of the derivatives analytically

0098-3004/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2004.10.015

ARTICLE IN PRESS N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

320

solves the problem. The analytic approach based on recursive formulation provides faster computation of the derivatives.

Using the identity rxfWGðxÞÞT W DYðxÞg ¼ ðWHðxÞÞT W DYðxÞ  ðWGðxÞÞT WGðxÞ the Qk becomes

2. Formulation of the problem

Qk ¼ r2 xU The nonlinear resistivity inversion relates observed data and model parameters by equation Y ¼ gðxÞ þ n;

(1)

where Y ¼ ðy1 ; . . . ; yN ; Þ is a vector representing observations at different half-electrode separations in Schlumberger sounding, N is number of half-electrode separations, gðxÞ ¼ ðg1 ðxÞ; . . . ; gN ðxÞÞ represents predicted data at different half-electrode separations, and n is the measurement noise. Various schemes to treat this non-linear problem are described in detail by Dimri (1992). Terminology used here closely follows Chernoguz (1995) with some differences arising due to Constable et al. (1987). Following Constable et al. (1987) the inverse problem is posed as a constrained optimization problem, set forth to minimize misfit X ¼ jjW Y  WgðxÞjj; subject to the constraint that roughness R ¼ jj@xjj is also minimized. This can be converted to an unconstrained problem by the use of Lagrange parameter m as follows: 1 1 U ¼ jj@xjj2 þ fðW DYðxÞÞT ðW DYðxÞ  w2 Þg; 2 2m where @ is (1987) as 0 0 B 1 B B B 0 @¼B B : B B @ : 0

(2)

N N matrix defined by Constable et al. 1

:::

:::

0

1 ::: 1 1 : : : :

::: ::: : :

0C C C 0C C: :C C C :A

:::

1

1

0

0

W is weighting matrix, w is acceptable misfit value and, m is a Lagrange parameter used to optimize the constrained functional ‘U’ (Smith, 1974) and DYðxÞ ¼ W Y  WgðxÞ: If we expand the functional in Taylor’s series at x ¼ xk (say) we get 1 Uðxk þ d; m; YÞ ¼ Uðxk ; m; YÞ þ J Tk d þ dT Qk d; 2 where 1 J k ¼ rxU ¼ @T @x  ðWGðxÞÞT W DYðxÞ m and 1 Qk ¼ r2 xU ¼ @T @  rx fðWGðxÞÞT W DYðxÞg: m

1 ¼ @T @  fðWGðxÞÞT WGðxÞ m  ðWHðxÞÞT W DYðxÞg where GðxÞ is Jacobian of gðxÞ and HðxÞ is Hessian of gðxÞ: If we define X ðWHÞT W DY ¼ WH j W DYj as q; j

where in H j is Hessian of gðxÞ evaluated at jth data point and DYj ¼ yj  gj ðxÞ; q is the nonlinear part of the Hessian, then minimization of the functional (2) using Newton’s method for ith iteration step di yields  1 di ¼  @T @ þ m1 fWGðxÞT WGðxÞ  qg 

@T @x  m1 WGðxÞT W DYðxÞ x¼xi :

ð3Þ

Thus xiþ1 ¼ xi þ di forms the iterative basis for the optimization of functional (2). Eq. (3) gives generalized OCCAM’s correction steps. By setting q as a null matrix, the equation gives the model correction of the popular OCCAM’s inversion algorithm. The OCCAM’s optimization in Eq. (3) can be viewed as two subalgorithms, where primary optimizes the functional U for different values of x and m and secondary optimizes only misfit function. The difficulty may arise when the primary suggests corrections in the direction of decreasing U and secondary moves in search of decreasing misfit without regarding the roughness. Thus due to lack of curvature information in the OCCAM’s inversion, secondary algorithm becomes blind in the direction of true minimum of w2 ; when the requirement of the primary algorithm to reduce U is overpowering. Another problem in OCCAM’s inversion is the choice of Lagrange’s parameter m. If we take mo1 for minimization of functional used in standard OCCAM’s inversion, the algorithm tends to be blind to minimize misfit function in the absence of the curvature information. Hence, there is need to incorporate curvature information in terms of Hessian matrix. If we take mp1 in Eq. (3) the nonlinear part carrying curvature information, will contribute to the convergence. In our work we include the curvature information in the OCCAM’s model correction steps. From Newton–Gauss method we have, xiþ1 ¼ xi þ ai di ;

(4)

ARTICLE IN PRESS N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

where ai is the extra smoothing parameter which uses curvature information as provided by Hessian terms. Following Chernoguz (1995) we have ( j0 ð0Þj00 ð0Þ if 0oai p1; ai ¼ (5) 1 otherwise; where j0 ð0Þ is the first-order derivative and j00 ð0Þ is the second-order derivative of Uðxi þ ai di ; m; YÞ that are computed as j0 ð0Þ ¼ dTi ½@@xi  m1 WGðxi ÞT W DYðxi Þ and j00 ð0Þ ¼  dTi ½@@xi  m1 fWGðxi ÞT W DYðxi Þ  WHðxi ÞT W DYðxi Þdi g; where the computation of Hessian matrix is involved in each iteration. Generally, derivatives are calculated using finitedifference techniques. The two most commonly used techniques are forward difference, which is accurate to first order and central difference, which is accurate to second order. In general, errors in forward difference approximation become unacceptably large during the last few iterations and often in such cases, central difference scheme is used. Computation of partial derivatives using central difference requires twice evaluations of forward functional than forward difference (Constable et al., 1987), and takes more computational time. Hence it is not worthwhile to use central difference for computation of Jacobian and Hessian matrices. Errors associated with computation of Hessian terms have a substantial effect on inversion algorithm. The use of analytical expressions instead of finite difference solves all the above problems. In OCCAM’s inversion Constable et al. (1987) have computed Jacobian matrix for Schlumberger sounding using analytical expressions and have concluded that the Schlumberger resistivity inversion with finite-difference derivatives takes 15 times more computational time than the analytical derivatives. Hence in our algorithm, the Hessian terms are calculated analytically. For very large data sets computation of the Hessian in each iteration becomes very tedious and time consuming, hence we follow Chernoguz (1995) to minimize the computations and at the same time to preserve the curvature information using Hessian matrix in the extra smoothing parameter a. In our algorithm Eq. (3) can be used for first two iterations to get the two consecutive values of a that are less than 1, and then the next correction step can be deduced using a relation s

ai ¼ 1  expðti Þ;

(6)

321

where t and s are positive scalar constants (sX1) in ith iteration, given as t ¼  ln ð1  a1 Þ; s ¼ log2 ½ln ð1  a2 Þ= ln ð1  a1 Þ:

ð7Þ

Eq. (5) is used to find the two consecutive values of a, i.e. a1 and a2 that are less than 1. In subsequent iterations, Eq. (6) can be applied directly to find the values of extra smoothing parameter to be used in correction steps of Eq. (4). This avoids computation of the tedious Hessian matrix in each iteration. In the modified algorithm we choose initial model at random as in the case of global optimization techniques. The above-described inversion algorithm is coded in MATLAB and can run on any machine. The weighting matrix used in the code is taken as identity if uncertainty in the data is not known. Step wise description of the MATLAB code is given in Appendix A.

3. Analytical expressions of derivatives The forward calculation for Schlumberger apparent resistivity over a layered earth is given as  ra ¼

AB 2

2 Z1 T 1 ðlÞJ 1

 AB l l dl; 2

(8)

0

where J1 is the first-order Bessel’s function of the first kind, AB/2 is the half-electrode spacing, l is electrode parameter, T 1 is the resistivity transform (Koefoed, 1976) that can be calculated recursively as T i1 ¼

T i þ ri1 tanhðlti1 Þ ; 1 þ T i tanhðlti1 Þ=ri1

where ri and ti are the resistivity and thickness of ith layer, respectively. For the last layer T M ¼ rM ; where M is number of layers. Using fast Hankel transforms apparent resistivity can be written as (Ghosh, 1971): X ra ¼ T 1 ðlK Þf K ; (9) K

where f K are the filter coefficients and K is the number of coefficients. 3.1. First-order derivatives The first-order derivatives of Eq. (9) are qra X qT 1 ðlK Þ ¼ f K; qrj qrj K

(10)

ARTICLE IN PRESS N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

322

where qT j1 qT j qT 1 qT 1 qT 2 ¼  ; qrj qT 2 qT 3 qT j qrj

(10a)

Table 2 Synthetic data generated using input model given in (a) Table 1A and (b) Table 1B AB/2 (m)

qT j ¼ ½1  tanh2 ðtj lK Þ=C j ; qT jþ1

Apparent resistivity (O m)

Log 10 apparent resistivity (O m)

(a) 2.50 3.00 3.70 4.60 5.80 7.20 8.40 10.00 12.50 16.00 20.00 25.00 30.00 37.00 46.00 58.00 72.00 84.00 100.00 125.00 160.00 200.00 250.00 300.00

8.69 9.05 9.71 10.79 12.61 15.26 17.95 22.08 29.67 42.05 57.78 78.35 98.68 125.21 154.72 185.90 212.54 229.22 245.38 261.82 275.04 283.37 289.15 292.49

0.94 0.96 0.99 1.03 1.10 1.18 1.25 1.34 1.47 1.62 1.76 1.89 1.99 2.10 2.19 2.27 2.33 2.36 2.39 2.42 2.44 2.45 2.46 2.47

(b) 2.50 3.00 3.70 4.60 5.80 7.20 8.40 10.00 12.50 16.00 20.00 25.00 30.00 37.00 46.00 58.00 72.00 84.00 100.00 125.00 160.00 200.00 250.00 300.00

8.22 8.36 8.64 9.10 9.91 11.08 12.25 14.04 17.28 22.70 30.02 40.60 52.44 70.47 94.83 127.02 161.69 187.91 217.66 253.67 288.15 313.16 332.26 343.99

0.91 0.92 0.94 0.96 1.00 1.04 1.09 1.15 1.24 1.36 1.48 1.61 1.72 1.85 1.98 2.10 2.21 2.27 2.34 2.40 2.46 2.50 2.52 2.54

(10b)

qT j ¼ tanhðtj lK Þ½1 þ T 2jþ1 =r2j qrj þ 2T jþ1 tanhðtj lK Þ=rj =C j ;

ð10cÞ

and C j ¼ ½1 þ tanh2 ðtj lK ÞT jþ1 =rj 2 ;

(10d)

also qT M ¼ 1; qrM

(11)

since T M ¼ rM for the last layer. 3.2. Second-order derivatives The second-order derivatives of Eq. (9) are X q2 T 1 ðlK Þ qra ¼ f : qrj qri K qrj qri K

(12)

Table 1 Input synthetic models Input model Input model (O m) log10 (resistivity) (O m)

Inverted model Error (O m) predicted log 10 (resistivity) (O m)

(a) Model 1a 8.10 15.40 166.40 215.50 300.40

0.91 1.19 2.22 2.33 2.48

0.90 1.31 1.91 2.32 2.48

1.02 0.76 2.04 1.02 1.00

(b) Model 2b 8.00 15.60 60.30 112.20 210.30 375.80

0.90 1.19 1.78 2.05 2.32 2.58

0.93 1.25 1.63 1.97 2.29 2.57

0.93 0.87 1.41 1.20 1.07 1.02

Error (O m) ¼ 10(log10(input model) log10 (inverted model). a No. of layers is 5 and thickness of each layer is 2.5 m. b No. of layers is 6 and thickness of each layer is 3.5 m.

ARTICLE IN PRESS N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

323

Fig. 1. (A–E) Comparative convergence of both the algorithms for different values of m for synthetic data given in Table 2B. Solid lines denote the observed data, asterisks (*) denote the predicted values using the modified algorithm and circles (o) denote the predicted values obtained using OCCAM’s algorithm. The Starting model is a half-space of 105 O m.

For Hessian, we!need to evaluate q2 T 1 q qT 1 ¼ : qri qrj qri qrj

(13)

Hence the Hessian matrix for each electrode spacing is formed as: 2 q2 T 3 2 q2 T 1 q2 T 1 T1 1 : : : qrq qr qr1 qr2 qr1 qr3 qr2 1 M 6 21 7 6 q T1 q2 T 1 q2 T 1 q2 T 1 7 6 qr @r 7 : : : 2 qr2 qr3 qr2 qrM 7 qr2 6 2 1 6 : : : : : : : 7 6 7 6 7 (14) 6 : : : : : : : 7: 6 7 6 : 7 : : : : : : 6 7 6 7 6 : : : : : : : 7 4 2 5 2 2 q T1 q T1 : : : : qqrT2 1 qr qr qr qr M

1

M

2

M

The Hessian matrix is N M M but corresponding to each electrode spacing it is a 1 M M symmetric matrix. Thus we need only either upper triangular or lower triangular elements to evaluate the Hessian matrix. Diagonal elements (i ¼ j) are computed using Eq. (13) as ! q2 T 1 q qT 1 ¼ qrj qrj qr2j ! qT j1 qT j q qT 1 qT 2 qT 3 ... ; ð14aÞ ¼ qT j qrj qrj qT 2 qT 3 qT 4

Table 3 Results obtained using synthetic data shown in (a) Table 2Aa and (b) Table 2Bb m

RMS misfit (modified algorithm)

RMS misfit (OCCAM’s algorithm)

No. of iterations

(a) 10 1.5 1.0 0.85 0.5 0.25 0.015

0.2987 0.0869 0.0645 0.0575 0.0380 0.0257 0.0076

0.2678 0.0703 0.0499 0.0429 15.5581 36.2329 2.5276

7 5 5 5 4 6 8

(b) 1.5 1.0 0.75 0.5 0.25

0.0741 0.0529 0.0416 0.0297 0.0170

0.0553 1.0576 1.2609 1.8186 1.2405

5 5 5 5 6

a

Layer thickness is 2.5 m. Layer thickness is 3.5 m.

b

where we get the sum of i terms as     qT j2 2 qT j1 2 q2 T 1 q2 T 1 qT 2 2 qT 3 2 ¼ . . . qT 4 qT j1 qT j qr2j qT 22 qT 3

ARTICLE IN PRESS N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

324

qT j qrj



!2 þ

qT 1 qT 2

q2 T 2 qT 23 !2

!  qT j2 2 qT 3 2 ... qT 4 qT j1

 qT j1 2 qT j qT 1 qT 2 q2 T k

þ ... þ ... 2 qT j qrj qT 2 qT 3 qT kþ1  2  2 qT kþ1 qT j1 qT 1 qT 2

... þ ... þ qT kþ2 qT j qT 2 qT 3

qT j2 qT j1 q2 T j qT 3 ... qT 4 qT j1 qT j qr2j

ð14bÞ

and nondiagonal elements are represented as   q2 T 1 q2 T 1 qT 2 2 qT 3 2 qT i qT i qT j ¼ ... qri qrj qT 4 qT j qri qrj qT 22 qT 3 ! qT 1 q2 T 2 qT 3 2 qT 4 qT i qT j þ ... qT 2 qT 23 qT 4 qT 5 qri qrj

with the inverted models. Filter coefficients used in the algorithm (Das and Kumar, 1976) are given in Appendix B. First column of Table 1A and B shows the input models with constant layer thickness used to generate the synthetic data. This synthetic data shown in Table 2A and B is used to test the inversion algorithm. The inverted model is compared with the true model to demonstrate the accuracy and efficacy of the modified algorithm. Errors of the inverted parameters are shown in the fourth columns of Table 1, which ranges between 0.8 and 2.0 O m. The convergence of both the modified and original OCCAM’s inversion algorithms for different values of m has been shown in Fig. 1(A–E) for input synthetic model No. 2 (Table 2B). The results are shown in Table 3A and B. It clearly shows that for mo1 the modified algorithm works whereas OCCAM’s algorithm fails. For m41; the OCCAM’s algorithm weighted more on the smoothness terms. Hence

qT 1 qT 2 q2 T 3 qT i qT i qT j ... þ ... qT 2 qT 3 qT 24 qT j qri qrj  qT 1 qT 2 qT 3 qT i1 qT j q qT i ... ; þ qT 2 qT 3 qT 4 qT i qrj qri qT j þ

3.5 3

µ

ð14cÞ where in 

qT j q qrj qT jþ1

¼

q qT jþ1

qT j qrj

¼ 2 tanhðtj lK ÞT jþ1 ½1 þ T jþ1

! (14d)

RMS misfit

2.5 2 1.5

=

1.5

=

1.0

=

0.75

=

0.50

=

0.25

1



qT j tanhðtj lK Þ=rj  qT jþ1

 r2j C j

0.5

ð14eÞ

0 2

3

þ tanhðtj lk ÞT jþ1 g ;

1

2

(A)

q Tj ¼  2 tanhðtj lk Þf1  tanh2 ðtj lk ÞgT 2jþ1 =frj qr2j

3

4

5

250 µ

ð14fÞ

=

1.50

=

1.00

=

0.75

=

0.50

=

0.25

q2 T j ¼  2r2j tanhðtj lK Þ½1 qT 2jþ1 ð14gÞ

Cj is the same as in Jacobian computations. It is advantageous to compute all the terms of the Hessian and the adjoining products simultaneously to save computational time and memory.

RMS misfit

200

 tanh2 ðtj lk Þ=frj þ tan hðtj lk ÞT jþ1 g;

6

No. of Iterations

150

100

50

0

4. Applications 4.1. Synthetic examples To assess the performance of the modified inversion algorithm, two synthetic models have been compared

(B)

1

2

3

4

5

6

No. of Iterations

Fig. 2. (A) Plot of iterations vs. RMS misfit for different values of m obtained by modified algorithm using synthetic data as shown in Table 2B and (B) plot of iterations vs. RMS misfit for different values of m using obtained by OCCAM’s algorithm using synthetic data as shown in Table 2B.

ARTICLE IN PRESS N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

the modified algorithm can be seen as generalized OCCAM’s inversion. It is worth pointing here that smaller values of m should be preferred to give less weight to the smoothness/roughness part as described in Eq. (2). The bias of algorithm to generate smooth models should be avoided whenever it is not required. Variation of RMS misfit with iteration number has been shown in Fig. 2 for different values of m: The modified algorithm shows stable convergence, while erratic variation in RMS misfit of OCCAM’s algorithm has been observed for mo1: We observed a remarkable difference in convergence between modified inversion algorithm and existing OCCAM’s algorithm for mo1: Results obtained using two different filters is shown in Appendix C.

325

half-space of 105 O m, which was far from the observed one. The modified algorithm searches for the lowest misfit until it becomes constant with further iterations as shown in Fig. 4. This is the global optimization strategy of our inversion scheme, where the existing OCCAM’s inversion fails for mo1: Results obtained using the modified algorithm for different values of m, are shown in Table 4.

2.5

µ 2

= 1.00

RMS misfit

= 0.10

4.2. Field example

Log10 [Apparent Resistivity] [ohm-m]

Log10 [Apparent Resistivity] [ohm-m]

The developed algorithm is used to invert 1D Schlumberger resistivity sounding data. The data from Southern Granulitic Terrain (SGT) of India (111340 5400 N, 78130 1800 E) over a 10 km long profile is given in Appendix D. This data has been collected by Deep Resistivity Sounding (DRS) Group of NGRI, Hyderabad, India using the Scintrex make Deep Resistivity Equipment, TSQ4-10 KVA. The convergence of the modified algorithm has been shown in Fig. 3(A–F). We have assumed starting model as a

Iteration No. 1

= 0.05

1.5

= 0.0125 = 0.005

1

0.5

0

1

2

8

(A)

4

5

6

Fig. 4. Plot of iterations vs. RMS misfit for different values of m using field data of SGT.

Iteration No. 2

5

3

No. of Iterations

Iteration No. 3 4.5

(B)

(C)

4

4

6 3.5

3

4

2

2

0

4

2

Observed Predicted

1

3 2.5 2

0 4

1.5 0

2

Iteration No. 4

4

Iteration No. 5 4

(D)

(E)

4

3.5

3.5

3.5

3

3

3

2.5

2.5

2.5

2

2

2

1.5 2 log 10(AB/2) [m]

4

4

Iteration No. 6

(F)

1.5

1.5 0

2

0

0

4 2 log 10(AB/2) [m]

0

2 4 log 10(AB/2) [m]

Fig. 3. (A–F) Convergence of modified algorithm for m ¼ 0:0125 using field data of SGT as given in Appendix D. Solid lines denote observed data, asterisks (*) denote predicted values using modified algorithm. The starting model is a half-space of 105 O m.

ARTICLE IN PRESS 326

N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

Table 4 Results obtained by modified algorithm using field data over SGT Iteration No.

1 2 3 4 5 6

RMS misfit m ¼ 1:0

m ¼ 0:1

m ¼ 0:05

m ¼ 0:0125

m ¼ 0:005

2.3147 1.9339 0.493 0.1751 0.1682 0.168

2.3147 2.3385 0.1925 0.0977 0.0952 0.0951

2.3147 2.4075 0.1703 0.0851 0.0795 0.0792

2.3147 2.4804 0.2209 0.0777 0.0593 0.0584

2.3147 2.4967 0.3221 0.0986 0.0529 0.0494

5. Concluding remarks The modified OCCAM’s inversion algorithm is very efficient, robust and simple to be used for 1D resistivity inversion. The modified algorithm finds an alternate way to simplify the nonlinear resistivity inversion problem without linearizing it. The code computes the first- and second-order derivatives of Schlumberger resistivity analytically and in general the convergence is obtained within 5–6 iterations. We have demonstrated the efficacy of the algorithm with two synthetic examples. For mo1; the modified gives less RMS error than OCCAM’s with equal number of iterations. For m41 our modified algorithm yields similar results as OCCAM’s. Hence our inversion algorithm may be considered as generalized OCCAM’s inversion.

2. 3.

4.

5. 6. 7.

Acknowledgements Authors express heartfelt thanks towards late Professor P.S. Moharir for his invaluable suggestions to develop the MATLAB codes. We are thankful to Professor Peter Weidelt and reviewers of the manuscript for many constructive comments that helped us in improving quality of the work. We acknowledge Dr. S.B. Singh, Head, Deep Resistivity Sounding Group, NGRI Hyderabad for providing DC geoelectric sounding data of SGT. One of the authors John Sagode acknowledges CSIR, India and TWAS, Italy for the award of fellowship.

Appendix A. Stepwise code description of 1D nonlinear resistivity inversion algorithm 1. Input to the program: (a) digital resistivity filters (as per the choice of user);

8. 9.

10.

(b) AB/2 vs. apparent resistivity (in log10); (c) starting model. Predicted apparent resistivity values corresponding to each AB/2 are computed using starting model by the forward algorithm. First-order and second-order derivatives are computed analytically which is multiplied with the filter coefficients to obtain apparent resistivity values. First-order derivatives are arranged to form Jacobian and second-order derivative are arranged to form Hessian matrix. Predicted resistivity values are matched with the observed data, and the difference is computed. With the above-mentioned steps the functional to be minimized is ready. Search for m that minimizes the objective functional depending upon user’s choice (Golden section search may be used for this purpose). Iteration steps are computed using Eq. (4) and values of a are computed using Eq. (5). If in any two consecutive iteration different values of a are obtained such that these are o1, then Hessian computation is bypassed and the value of a obtained by consecutive a values are used in iteration steps, otherwise a ¼ 1: Above procedure is repeated until two consecutive misfit values o1 having difference of o0.001 is obtained.

Appendix B. Schlumberger coefficients as given by Das and Kumar (1976) are

Sampling interval Shift Total number of filter points

ln 10/6 +0.1343155 14

ARTICLE IN PRESS N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

Number of filter points with negative indices Number of filter points with positive indices

3

327

Appendix D. Field data acquired over Southern Granulite Terrain (SGT) of South India

10

Index number

Filter coeff.

3 2 1 0 1 2 3 4 5 6 7 8 9 10

0.0140 0.0282 0.0838 0.2427 0.6217 1.1877 0.3954 3.4531 2.7568 1.2075 0.4595 0.1975 0.1042 0.0359

Appendix C. Results obtained using different filters: m ¼ 0:012

AB/2 (m)

Apparent resistivity (O m)

Log 10 (apparent resistivity) (O m)

4 5 6 15 25 40 60 120 150 200 300 500 600 800 1000 1500 2500 5000

67 68 74 170 260 320 280 150 150 200 540 2400 4200 6300 6600 6400 7400 9950

1.8261 1.8325 1.8692 2.2304 2.415 2.5051 2.4472 2.1761 2.1761 2.301 2.7324 3.3802 3.6232 3.7993 3.8195 3.8062 3.8692 3.9978

Model Log 10 (resistivity) Layer thickness 6.0 m Input model

Computed model using Das Kumar filter

Computed model using Ghosh filter

1.8261 1.8325 1.8692 2.2304 2.4150 2.5051 2.4472 2.1761 2.1761 2.3010 2.7324 3.3802 3.6232 3.7993 3.8195 3.8062 3.8692 3.9978

1.8243 1.8354 1.9806 2.1132 2.2144 2.3003 2.3882 2.4883 2.6046 2.7369 2.8825 3.0380 3.1997 3.3642 3.5287 3.6908 3.8486 4.0006

1.8213 1.8275 1.9759 2.1112 2.2142 2.3017 2.3907 2.4916 2.6084 2.7408 2.8865 3.0418 3.2033 3.3675 3.5316 3.6933 3.8505 4.0018

RMS

0.0039

0.0040

Iteration

6

6

References Chernoguz, N.G., 1995. A smoothed Newton–Guass method with application to bearing only position location. IEEE Transactions on Signal Processing 43 (8), 2011–2013. Constable, S.C., Parker, R.L., Constable, C.G., 1987. OCCAM’s inversion: a practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 52 (3), 289–300. Dimri, V.P., 1992. Deconvolution and Inverse Theory. Elsevier Science Publishers, Amsterdam 230pp. Das, U.C., Kumar, R., 1976. Improved digital filters for computing Schlumberger and Wenner type curves, unpublished. Ghosh, D.P., 1971. Inverse filter coefficients for the computation of apparent resistivity standard curves for a horizontally layered earth. Geophysical Prospecting 29, 769–775. Koefoed, O., 1976. Error propagation and uncertainty in the interpretation of resistivity sounding data. Geophysical Prospecting 24, 31–48. LaBrecque, D.J., Miletto, M., William, D., Ramirez, A., Owen, E., 1996. The effects of noise on OCCAM’s inversion of resistivity tomography data. Geophysics 61 (2), 538–548. Qian, W., Jeffrey, G.T., Scott, H.J., Richard, L., Dennis, A., 1997. Inversion of airborne electromagnetic data using an OCCAM technique to resolve a variable number of layers. In: Proceedings of the Symposium on the Application of

ARTICLE IN PRESS 328

N. Vedanti et al. / Computers & Geosciences 31 (2005) 319–328

Geophysics to Environmental and Engineering Problems (SAGEEP), pp. 735–739. Siripunvaraporn, W., Egbert, G.D., 1996. An efficient data space OCCAM’s inversion for MT and MV data. EOS,

Transactions, American Geophysical Union 77(46) (Suppl.) p. 156. Smith, D.R., 1974. Variational Methods in Optimization. Prentice-Hall, Inc., Englewood, NJ.