An inversion model applied to DC soundings interpretation 

C . D E L G I U D I C E * 

Received on S e p t e m b e r 10, 1982 

A B S T R A C T 

T h e inversion technique is used in DC soundings intepretation to determine the 
t h i c k n e s s e s a n d the t r u e resistivities of the layers, starting from the field apparent 
r e s i s t i v i t i e s . 

T h e r e l a t i o n s h i p between the predicted a p p a r e n t resistivities and the earth 
p a r a m e t e r s is not linear. S t a r t i n g from this relationship a methodology is described 
t o o b t a i n a wellposed system of M linear equations. This system permits to 
c a l c u l a t e , by m e a n s of a n iterative procedure, the e a r t h p a r a m e t e r s that minimize 
t h e d i f f e r e n c i e s (error) between the field and the predicted a p p a r e n t resistivities. 

T h r e e different iterative procedures are described. Practical examples have 
s h o w n t h a t all the iterative procedures are reliable and give comparable results in 
t e r m s of m i n i m u m e r r o r reached and CPU time. 

R I A S S U N T O 

La t e c n i c a dell'inversione è u s a t a nella interpretazione di sondaggi elettrici 
v e r t i c a l i p e r d e t e r m i n a r e gli spessori e le resistività vere degli elettrostrati 
( p a r a m e t r i del terreno) p a r t e n d o dalle resistività a p p a r e n t i . La relazione tra le 
r e s i s t i v i t à a p p a r e n t i ed i p a r a m e t r i del terreno non è lineare. 

* Presently employed by GEOMATH - Via Cavour, 43 - PISA (Italy). 



1 6 6 C. D E L G I U D I C E 

P a r t e n d o da q u e s t a relazione viene descritta u n a metodologia già nota per 
o t t e n e r e un s i s t e m a ben posto di M equazioni lineari. 

Q u e s t o s i s t e m a p e r m e t t e di calcolare, per mezzo di u n a procedura iterativa, i 
p a r a m e t r i del t e r r e n o che m i n i m i z z a n o la differenza tra le resistività a p p a r e n t i di 
c a m p o e q u e l l e calcolate. Vengono inoltre descritte tre diverse procedure iterative. 
Gli e s e m p i p r a t i c i h a n n o m o s t r a t o che t u t t e le procedure iterative utilizzate sono 
a t t e n d i b i l i e p a r a g o n a b i l i tra loro per q u a n t o riguarda i risultati del procedimento 
di m i n i m i z z a z i o n e e i tempi di calcolo utilizzati. 

1 . S T A T M E N T O F P R O B L E M 

T h e a p p a r e n t resistivities in ohm- m versus the AB/2 distances 
in M r e p r e s e n t the field d a t a . In t e r m s of the voltage V,, currents I, 
a n d g e o m e t r i c factor K, (i = 1, N) where N is the n u m b e r of 
s a m p l e s , t h e a p p a r e n t resistivity is: 

K, 
AV, 

I, 

We w a n t to find an e a r t h model, consisting of a distribution of 
t r u e resistivities a n d thicknesses t h a t minimize the error between 
t h e field a p p a r e n t resistivities and the predicted ones. 

T h e observed a p p a r e n t resistivities are, in vectorial form: 

Pa, 

P a2 

If w e call P' the resistivities predicted by a forward earth model P 

c o n s i s t i n g of the t r u e resistivities a n d the thicknesses of the ^ ^ ^ 

l a y e r s 



AN I N V E R S I O N M O D E L APPLIED, ECC. 1 6 7 

T h e a b o v e values are functions of the p a r t i c u l a r earth model P, 
( / = 1, M) w h e r e generally M < N. 

P i " 
P2 

T h e n P' = f(P) [1] 

T h e r e l a t i o n s h i p [1] between P and P' is not linear. By expanding 
t h e e q u a t i o n [1] in a Taylor series and by keeping only the linear 
t e r m s , w e o b t a i n : 

- - - - 8«P) _ 

Pi = f ( P + A P) = f (P) + - y ^ - AP 

w h e r e AP is t h e model i m p r o v e m e n t . 

W e n e e d the P' to fit the observed d a t a : 

P., = P' 

O r : 

ö f ( P ) _ 
Po = f ( P ) + " T T " AP 



1 6 8 C. D E L G I U D I C E 

C a l l i n g d = pu — f (P), this is: 

a f ( P ) -
d = AP 

dP 

Or: 

d = [A) AP [2] 

w h e r e [A] is t h e sensitivity matr ix a n d AP is the vector p a r a m e -
t e r c o r r e c t i o n (Jupp, Vozof 1975, Lanczos 1961, M a r q u a r d t 1963). 
W r i t i n g o u t the matrices of equation [2] we have: 

d , 
8 f, (P) 

¿ P , 
d 2 • 

• 

• 
• 
• 

• 
• 
• 

• 
• 

d f.v (P) 
0 P , 

8 f , ( P ) 

d P „ 

3 t \ , ( P ) 

3 P M 

A P , 
A P 2 

A P. 

[3] 

T h i s is a system of N equations in M unknowns with N > M. The 
a b o v e s y s t e m of e q u a t i o n s is then over-determined and generally 
ill-posed in the sense that small changes in the d a t a lead to large 
c h a n g e s in the solutions. To solve this r e d u n d a n t system of 
e q u a t i o n s , we ca n apply the method of least squares. 

C a l l i n g e = d —[A] AP we have: 

N • 
S « 

; - i 
e t e 



AN I N V E R S I O N M O D E L A P P L I E D , ECC. 1 6 9 

w h e r e eT is the transpose of e or: 

e2 = [ d —[A] AP]T • [d — [A] AP] 

F o r e2 to be a m i n i m u m ; its derivative with respect to AP must 
be z e r o . 

T h e n d i f f e r e n t i a t i n g the above equation we have: 

[A]T [A] AP = [A]t d [4] 

We h a v e t r a n s f o r m e d the old system [2] of N equations in M 
u n k n o w n s i n t o a n M x M system [4], From the above expression 
[4] w e h a v e : ' 

AP = [ A T A ] " 1 [ A ] T d [5] 

2 . I T E R A T I V E P R O C E D U R E S 

In t h e i t e r a t i v e procedure we do not use the algorithm [5] but 
its m o d i f i c a t i o n m a d e by M a r q u a r d t (Marquardt 1963): 

AP = [AT A -I- k 2 1]" 1 [A]T d [6] 

w h e r e k2 is called « M a r q u a r d t p a r a m e t e r » and I is the unit 
m a t r i x . 

T h e s y s t e m of e q u a t i o n s [6] is well-posed. 
T h e a l g o r i t h m [6] has the advantage over [5] that the region of 

c o n v e r g e n c e is g r e a t e r a n d the a m p l i t u d e of the p a r a m e t e r correc-
t i o n A P i s s m a l l e r . 

T h e g e n e r a l expression for the iterative procedure is: 

pimi _ p (m-1) ^ p (m) [7] 



1 7 0 C. D E L G I U D I C E 

w h e r e m i n d i c a t e s the iteration n u m b e r a n d : 

Ap(m, = [ (A T A)<m_1) + ( k 2 I ) (AT d ) , m - " [8] 

I n e a c h i t e r a t i o n , we j^an use i values of k2 t h a t give i 
c o r r e s p o n d i n g values of A P. We then select the k2 t h a t gives the 
m i n i n u m e r r o r b e t w e e n the a p p a r e n t resistivities and the pre-
d i c t e d ones. Then: 

( m + i) = P ( m - 0 + AP (m + 0 [9] 

To a s c e r t a i n t h a t a m i n i m u m h a s been found so far, we have 
t o t e s t n e i g h b o u r i n g points, increasing the n u m b e r of k 2 p a r a m e -
t e r s a n d decreasing the range of variability of k 2 . Thus the 
e x p r e s s i o n [9] becomes: 

p (m + i + n) _ p (m - i - n) ^ p (m + i + n) [ 1 0 ] 

w h e r e n is t h e n u m b e r of times t h a t the "neighbourhood t e s t " is 
p e r f o r m e d . In expression [10] AP<", + , + ") is equal to: 

^ p (m + i + n) _ j^A T A ) <m ~ ' ~ + ( k 2 I) <m+ ,' + n) ] " ' (AT (J) (m- n) 

[11] 

F i g u r e 1 shows the chart of a c o m p u t e r p r o g r a m based on this 
f i r s t i t e r a t i v e a p p r o a c h . 

In o r d e r to find a faster procedure two other approaches have 
b e e n e x a m i n e d . 

In the second a p p r o a c h , the equation [10] has been modified: 

p (m + i + n) _ p (m — i +«) ^ p ( m + i + n) [ 1 2 ] 

H e n c e in the " n e i g h b o u r h o o d t e s t " we s u b s t i t u t e for the vector 
P t h e c u r r e n t value of it t h a t gives the m i n i m u m error Pa — P' 
w h e r e for A P we use the s a m e expression [11]. 



AN I N V E R S I O N M O D E L APPLIED, ECC. 1 7 1 

Figure 1 - First iterative a p p r o a c h flow chart 

In t h e t h i r d a p p r o a c h , we use the expression [12] also modify-
i n g t h e AP expression t h a t becomes: 

A P ( m + i + H) = [ ( A T A ) < T i + n) 4- ( k
2 J)(m + i + n) f j j T ^ ( m - i + n) ] - ' (AT d) 

[ 1 3 ] 



1 7 2 C. D E L G I U D I C E 

3 . P R A C T I C A L A P P L I C A T I O N S 

T w o sets of ten electrical soundings have been run utilizing 
t h e a b o v e described three approaches; no problem of convergence 
h a s b e e n e n c o u n t e r e d (fig. 2). 

In t h e first set, a c o m m o n t r u e resistivity value has been used 
a s a n initial e a r t h model for all the layers of each VES obtaining 
a s a f i r s t solution of the equation [1] a straightline (fig. 3). 

T h e a v e r a g e error, between the field resistivities and the 
c o m p u t e d ones, h a s been of 113.91%. In the second set, the initial 
e a r t h m o d e l s h a v e been obtained by a very rapid analysis of the 

ITERETION NUMBER 

F i g u r e 2 - E r r o r ' s trend relative to the three iterative a p p r o a c h e s 



AN I N V E R S I O N M O D E L APPLIED, ECC. 1 7 3 

V E S N ° 6 

F i r s t set 

+ f i r s t s o l u t i o n 

. f i e l d 

+ + + + + + + + + + + + + + + + + 

10 15 SO 3 0 

A B / 2 ( m ) 

100 150 2 0 0 3 0 0 4 0 0 

F i g u r e 3 - First solution relative to an initial model apparteining 
to the first set 

V E S N° 2 

• S e c o n d s e t 
• 

+ + • + + f i r s t s o l u t i o n 
+ + . f i e l d 

+ • + 
• 

' . + 

+ • + • 
+ 

+ 4- + + • 
r 1 1 1 1 1 1 1 1 i I 

I 15 2 3 5 7 10 IC 2 0 3 0 5 0 1 0 0 150 

A ß / 2 ( m ) 

F i g u r e 4 - First solution relative to an initial model apparteining 
to the second set 

f i e l d c u r v e , giving to the t r u e resistivities an a p p r o x i m a t e 
a s y m p t o t i c value a n d deducing the thicknesses from the abscissa 
of t h e inflection point (fig. 4). 



1 7 4 C . D E L G I U D I C E 

A • 
o 

Iterative First set Second set 
A p p r o a c h of V E S of V E S 

1 O O 

2 • • 
El 
A 3 A A 

0.9 1.0 1.1 1.2 1.3 1 4 
t 
A P P R O A C H 1 

F i g u r e 5 - Final e r r o r versus CPU t i m e 

T h e a v e r a g e error, between the field resistivities and the 
c o m p u t e d ones, h a s been of 27.90%. Figure 5, referred to the first 
a n d s e c o n d set of VES, shows the final e r r o r versus the CPU time 
r e l a t i v e to the first iterative a p p r o a c h . For both sets of VES the 
r e s u l t s o b t a i n e d utilizing the three different iterative approaches 
a r e c o m p a r a b l e r e g a r d i n g to the reached m i n i m u m error (differen-
ce in e r r o r 0.2% - 0.3%). 

R e g a r d i n g the CUP t i m e the first iterative a p p r o a c h is faster 
t h a n t h e o t h e r t w o w h e n the initial e r r o r is high (more than 100%) 
w h i l e t h e t h r e e iterative approaches are c o m p a r a b l e when the 
i n i t i a l e r r o r is small (20% - 30%). Results obtained by using the 
t h r e e i t e r a t i v e p r o c e d u r e s applied to the above practical examples 
h a v e s h o w n t h a t all the procedures are reliable and comparable 
a m o n g t h e m s e l v e s regards to the m i n i m u m error reached and the 
CPU t i m e . 



AN I N V E R S I O N M O D E L A P P L I E D , ECC. 1 7 5 

A C K N O W L E D G E M E N T S 

T h e w o r k w a s been carried out within a fellowship (N. 216849). 
f i n a n c e d by t h e italian «Consiglio Nazionale delle Ricerche» 
(CNR) a n d " N o r t h Atlantic Treaty Organization" (NATO) at the 
C o l o r a d o School of Mines-Department of Geophysics. 

T h e a u t h o r is p a r t i c u l a r l y grateful to Prof. G.V. Keller, head of 
t h e G e o p h y s i c a l d e p a r t m e n t and to the Prof. C.H. Stoyer. 



1 7 6 C . D E L G I U D I C E 

R E F E R E N C E 

JOANSEN H . K . , 1977 - A man/computer interpretation system for resistivity 
soundings over a horizontally stratified earth. " G e o p h y s . P r o s p e c t i n g " , 
25, p p . 667-691. 

J U P P D . , V O Z O F K . , 1 9 7 5 - Stable iterative method for the inversion of geophysical 
data. " G e o p h y s . J . F . Astron. Soc.", 42, pp. 9 5 7 - 9 7 6 . 

L A N C Z O S C., 1 9 6 1 - Linear differential operators. D . Van Nostrand Company 
L i m i t e d . 

M A R Q U A N R D T D . , 1 9 6 3 - An algorithm for least-square estimation of non-
linear parameters. " J . Soc. Ind. Appl. Math.", II, 2, J u n e . 

P A R K E R R., 1977 - Understanding inverse Theory. "Ann. Rev. E a r t h Planet. S." 5, 
p p . 35-64. 

V A R G A R., 1 9 6 2 - Matrix iterative analysis. Prentice - Hall, Inc.