Microsoft Word - 001.docx


 CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 66, 2018 

A publication of 

 

The Italian Association 
of Chemical Engineering 

Online at www.aidic.it/cet 

Guest Editors: Songying Zhao, Yougang Sun, Ye Zhou 
Copyright © 2018, AIDIC Servizi S.r.l. 

ISBN 978-88-95608-63-1; ISSN 2283-9216 

Finite Difference Modeling for Resistivity Tomography Based 

on Reconstruction Algorithms 

Peng Liua,b, Jianhua Yuea,* 

aCollege of Resource & Geosciences, China University of Mining & Technology, Xuzhou 221008, China 
bYankuang Donghua Construction Co., Ltd, Zoucheng 273500, China 

15052020148@163.com 

The finite difference forward modeling solves the point source arbitrary 2D geoelectrical section resistivity model 

herein. A system of linear equations is finally formed. Two reconstruction algorithms, i.e. conjugate gradient 

method and Cholesky algorithm, are introduced into the numerical analysis to derive the theoretical formula. 

With uniform media and layered media chosen as the physical model for forward modeling, the numerical 

simulation results that low and high resistance objects correspond to different reconstruction algorithms in two 

media models are analyzed. The findings show that the Cholesky reconstruction algorithm works well in the 

numerical simulation of sounding low resistance object, while there is a little difference between the two 

reconstruction algorithms for high resistance object. But also, in the even semi-infinite medium, when the point 

source lies at the center of the surface mesh, the relative error Rms to the potential value and the theoretical 

value of the unit point source geoelectric field is obtained by the Cholesky algorithm, and the maximum error 

occurred when it converges for 33 iterations is more than 3.0%. As polar distance increases, relative error also 

has no tendency to amplify. Compared to the conjugate gradient method, the Cholesky algorithm not only greatly 

requires less in memory but also greatly increases the computing speed. The study provides the clues to future 

finite difference numerical simulation and interpretation of resistivity tomography. 

1. Introduction 

Resistivity Tomography is such that the medium resistivity distribution inside the detection area is studied by 

observing the potential or potential difference in different directions around detection area using the potential 

field excited by DC power source (LaBrecque, 1992). Since the potential field is vulnerable to the rock lithology 

and structure, the projection data space is generally smaller than that of image data needed to be inverted. In 

this sense, it often has multiple solutions. It is required for people to perform numerical simulations first, make 

a survey on how the geological structure affects the observation of electric fields, and find out what the 

characteristics of electric field propagation seem under geological conditions, in order to lay the foundation for 

the subsequent study and interpretation on inversion approaches. Numerical simulation plays an important role 

in processing electric field data, including finite element method, the integral equation method, the Newton-

Raphson method and so on. Shima et al. (1987) proposed the resistivity back-projection technology and 

thereafter Yorkey et al., calculated the potential distribution generated from a 2D structure to a 3D electrode 

using the finite difference method. Due to too heavy computation burden in the finite element method in forward 

modeling, Barker (1992), Lietal (1992) introduced the finite difference method into the forward modeling.  

In the study of resistivity tomography, the forward and inverse modeling using simulated annealing and genetic 

algorithm gets hot, while traditional and classical numerical methods has gradually been out of fashion (such as 

the conjugate gradient method and Cholesky decomposition). There are even less involvement on the forward 

modeling of the finite difference methods (Li and Oldenburg, 1994; Loke and Barker, 1995; Murai, 1985; Noel 

and Xu, 1991). 

This paper builds a forward model to carry out a case study on the relationship between conjugate gradient 

method, Cholesky decomposition algorithm and anomaly space distribution feature (geometric, burial features) 

in finite difference numerical simulation, and perform forward numerical simulation of the model. What the 

relationship of the forward field curve with the spatial position of the deeply buried object seems to be known. 

                                

 
 

 

 
   

                                                  
DOI: 10.3303/CET1866171 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Please cite this article as: Liu P., Yue J., 2018, Finite difference modeling for resistivity tomography based on reconstruction algorithms, 
Chemical Engineering Transactions, 66, 1021-1026  DOI:10.3303/CET1866171   

1021



The results provide the clues to interpretation on application conditions and observation data for different 

geological models with the conjugate gradient method and Cholesky decomposition algorithm in finite difference. 

2. Mathematical principle 

Resistivity tomography is highly nonlinear. Discretization of imaging area will cause certain errors in the imaging 

process when using numerical solutions. In order to capture more real and reliable imaging, it is required to 

depend on forward modeling or other apriori information to get actual imaging. Thus, we should build a 

mathematical model for the engineering problem, i.e. the differential integral equation with relevant boundary 

and initial conditions, before making analysis. 

To solve these mathematical models, there are methods roughly composed of two types: analytical and 

numerical methods. It is well known that the analytical method has some limitations. When the boundary 

conditions of a system are complex, it is often impossible to find analytical or approximate solutions to the 

problem. Therefore, for most engineering problems, the numerical solution gets more practical. In the case of 

dense mesh, the finite difference method can quickly identify the electrical anomaly features of complex 

geoelectrical sections and calculate electrical anomalies under complex terrain conditions. With this method, 

how the terrain affects electrical anomalies can easily come to light. 

Underground electric field built by resistivity method satisfies the differential equation (Price, 1979; Qian et al., 

1995; Sasaki, 1989, 1992): 

∇ ∙ (𝜎 ∙ ∇𝑈) = −𝐼𝛿(𝑥 − 𝑥0)𝛿(𝑦 − 𝑦0)𝛿(𝑧 − 𝑧0)                                                                                                    (1) 

Where, σ is the electrical conductivity of the rock; U is the potential at any point on the ground or in the medium; 

I is the supply current intensity; (𝑥0, 𝑦0, 𝑧0) is the position coordinate of the power supply point; δ is the Dirac 

function; σ and U are the functions of (x, y, z). 

For 2D geoelectrical profiles, there is 
𝜕

𝜕𝑦
[𝜎(𝑥，𝑦，𝑧)] = 0, then (1) 

∇ ∙ [𝜎(𝑥，𝑧)∇𝑈(𝑥，𝑦，𝑧)] = −𝐼𝛿(𝑥 − 𝑥0)𝛿(𝑦 − 𝑦0)𝛿(𝑧 − 𝑧0)                                                                             (2) 

After using Fourier transform 

ϕ(𝑥，𝑘，𝑧) = ∫ 𝑈(𝑥，𝑦，𝑧) cos(𝑘 ∙ 𝑦)𝑑𝑦
∞

0
                                                                                                         (3) 

Transform the potential U in the (x, y, z) into the potential φ in the (x, k, z). When 𝑦0 = 0, it derives from (2): 

∇ ∙ [𝜎(𝑥，𝑦)∇𝜙(𝑥，𝑘，𝑧)] − 𝑘2 𝜎(𝑥，𝑧)𝜙(𝑥，𝑘，𝑧) = −
1

2
𝛿(𝑥 − 𝑥0)𝛿(𝑧 − 𝑧0)                                                    (4) 

Where: σ is the conductivity; k is the space wave number; δ is the Dirac function. 

After the discretization and variational process, the potential distribution generated by the 3D electrode in the 

2D structure can be obtained. At last, the potential in the wave number domain is transformed to the spatial 

domain, and the forward simulation calculation is achieved. 

3. Reconstruction algorithms 

As the key to forward modeling, the reconstruction algorithm matters the imaging speed and quality, and it is 

also the key and aporia of tomography. 

After discretizing and varying the formula (4), we obtain: 

U(𝑥，𝑦，𝑧) =
2

𝜋
∫ 𝑣(𝑥，𝜆，𝑧) cos(𝜆𝑦)𝑑𝜆

∞

0
                                                                                                          (5) 

Using formula (5), the spatial domain potential field distribution of the stable current field can be available. 

It is not easy to solve the formula (5) directly. First, the area is discretized into a 3D mesh of 𝑁𝑥 × 𝑁𝑦 × 𝑁𝑧. The 

node number is i = 1,2, ⋯ , 𝑁𝑥  in the x direction, j = 1,2, ⋯ , 𝑁𝑦 in y direction, and k = 1,2, ⋯ , 𝑁𝑧 in z direction, set 

∆𝑉𝑖,𝑗,𝑘 represents the volume element near the node (i, j, k), (5) integrates in ∆𝑉𝑖,𝑗,𝑘, then use the Gauss theorem 

∭ ∇ ∙ [𝜎∇𝜑]
 

∆𝑉𝑖,𝑗,𝑘
𝑑𝑣 = ∬ 𝜎∇𝜑 ∙ 𝑛𝑑𝑠

 

𝑠𝑖,𝑗,𝑘
                                                                                                                (6) 

Where, 𝑠𝑖,𝑗,𝑘  is the surface of ∆𝑉𝑖,𝑗,𝑘，then 

− ∬ 𝜎
𝜕𝜑

𝜕𝑛
𝑑𝑠

 

𝑠𝑖，𝑗，𝑘
=

𝐼

∆𝑉𝑖，𝑗，𝑘
∭ 𝛿(𝑥 − 𝑥0)𝛿(𝑦 − 𝑦0)𝛿(𝑧 − 𝑧0)𝑑𝑥𝑑𝑦𝑑𝑧

 

∆𝑉𝑖，𝑗，𝑘
= {

1，(𝑥0，𝑦0，𝑧0) ∈ ∆𝑉𝑖，𝑗，𝑘
0. (𝑥0，𝑦0，𝑧0) ∉ ∆𝑉𝑖，𝑗，𝑘

        (7) 

1022



Again, the center difference is used to calculate𝜕𝜑 /𝜕𝑛, for each node (i, j, k), the above formula can be written 

as the following difference equation: 

𝐶𝑡𝜑𝑖，𝑗，𝑘−1 + 𝐶𝑏0𝜑𝑖，𝑗，𝑘+1+𝐶𝑙 𝜑𝑖−1，𝑗，𝑘 +𝐶𝑟 𝜑𝑖+1，𝑗，𝑘 +𝐶𝑓 𝜑𝑖，𝑗−1，𝑘 +𝐶𝑏𝑎𝜑𝑖，𝑗+1，𝑘 +𝐶𝑝𝜑𝑖，𝑗，𝑘 =

{
1，(𝑥0，𝑦0，𝑧0) ∈ ∆𝑉𝑖，𝑗，𝑘
0. (𝑥0，𝑦0，𝑧0) ∉ ∆𝑉𝑖，𝑗，𝑘

                                                                                                                              (8) 

Where, 𝐶𝑡, 𝐶𝑏0, 𝐶𝑙 , 𝐶𝑟 , 𝐶𝑓 , 𝐶𝑏𝑎, 𝐶𝑝 are the connection coefficients around the node (i, j, k) and their own nodes, 

respectively, written in a matrix form and the equation set function (9) is obtained. 

Aφ=S                                                                                                                                                                  (9) 

4. Conjugate gradient method 

From minimized function F(x)=xTAx/2-bTx from (9), let r(0)=b-Ax(0), p(0)=r(0)，then 

 

   

( ) ( ) ( ) ( )

( 1) ( ) ( ) ( )

( 1) ( ) ( ) ( )

( ) ( 1) ( 1) ( ) ( )

( 1) ( 1) ( ) ( )

, / ( , ),

,

0,1, 2,

, / , ,

.

i i i i i

i i i i

i i i i

i i i i i

i i i i

a r r p Ap

x x a p

ir r a Ap

r r r r

p r p









 

 




  


  


 
                                                                                                   (10) 

Since A is a capacity matrix, it is a large sparse symmetric positive definite band matrix, as shown below: 

0

1

0

0

0

0

0

0

ba

r

t p b

t

C

C

A C C C

C

C

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
                                                                              (11) 

It can be seen that there are at most 7 non-zero elements in each line of the matrix A. S is the right-end term 

relevant to the supply current, and it only has values on the power supply node ρ, i.e. S=(0,…,0, Sp, 0,…,0)T, 

while Sp=I, so the whole process of the numerical decomposition is to solve the product of matrix A and a column 

vector. A is decomposed into two 2D arrays, i.e. a real array (storing non-zero elements of matrix A) and an 

integer array (storing the index of the positions of the elements in the real array in matrix A), from which the 

product of the matrix A and any column vector is obtained to further acquire the potential φ(x,y,z). 

5. Cholesky decomposition algorithm 

The Cholesky decomposition is  

A≈CCT                                                                                                                                                              (12) 

where C is the lower triangular matrix, obtained from the diagonal matrix. The diagonal element of D is defined 

by: 

2
/

jj jj jk kk

k j

d a a d


 
                                                                                                                                          (13) 

Where ajk is the element of matrix A. Obviously, the element of ajk=0 in the summation has no contribution. D 

can be simply obtained. C is determined by the formula 

1023



1/ 2
C UD


                                                                                                                                                         (14) 

Where U is the lower triangular matrix, its diagonal elements ujj=djj, for non-diagonal elements, k<j, then

jk jk
u a

. It is obvious that the decomposition factor C is sparse, and the off-diagonal element is also a 

symmetric positive definite band matrix. then (9) can be rewritten as: 

 
1

1 1T T
C A C C x C b


   

                                                                                                                                        (15) 

If (CCT)-1 is the inverse approximation of the matrix A, C-1A(CT)-1 will be an approximate unit matrix. Kershaw 

expounded the changes in the eigenvalue of the matrix C-1A(CT)-1in practices, most of which approach to 1.0 

and close to the unit matrix, which shows that the Cholesky decomposition of (9) is feasible. Substitute this 

matrix into (10), that is, r0=b-Ax0, p0=(CCT)-1r0, and then 

    

     
 

1

1

1

1 1

1 1

1

1 1

, / ,

, / ,

0,1, 2,

T

i i i i

i i i i

i i i i

T T

i i i i i

T

i i i i

a r CC r p AP

x x a p

r r a Ap

r CC r r CC r

p CC r p

i











 

 



 

 



 

  





  

                                                                                                                     (16) 

The iterative decomposition of (9) is performed quickly by (15), from which the product of matrix A and any 

column vector is also obtained and the potential φ(x,y,z) is solved. 

6. Forward modeling and analysis  

According to the established forward model, we obtain the comparison of the numerical simulation results from 

different reconstruction algorithms on the similar anomaly object at different resistivity. 

6.1 Numerical simulation of different reconstruction algorithms for high and low resistance object in 
uniform medium 

The test model is shown in Fig. 1(a). The dipole electrode sounding observation system is adopted with the 

electrode spacing of 5m, the model thickness of 33.7m, the surrounding rock resistivity of 100, the anomaly 

object ρ1 (red) has a resistivity of 1000 Ω·m, and the electrical property contrast is 10:1. The anomaly object ρ2 

(green) has a resistivity of 10 Ω·m and the electrical property contrast is 1:10. The model simulates high- and 

low-resistance anomaly objects in a uniform space. The forward model and coordinate system are shown in Fig. 

1(a): 

 

Figure 1: (a) Forward modeling of high- and low-resistance in uniform space  

 

Figure 1: (b) Forward simulation imaging based on conjugate gradient reconstruction algorithm  

1024



  

Figure 1: (c) Forward simulation imaging based on Cholesky reconstruction algorithm  

As shown in Figs. 1(b) and 1(c), it can be seen that imaging position and range of anomaly objects are subjected 

to change with the forward reconstruction algorithms. The conjugate gradient method has a large imaging range 

on low resistance object. In relation to the position of low resistance model, the imaging range exceeds a lot, 

while the Cholesky has a more accurate imaging position for the low resistance object, close to the size of the 

low resistance model. As compared to the imaging sizes of the high-resistance object, the two reconstruction 

algorithms have little difference. 

6.2 Numerical simulation of different reconstruction algorithms for low-resistance anomaly object in 
layered media  

The experimental model is shown in Fig. 2(a). A dipole electrode sounding observation system is used with the 

electrode spacing of 5m, the model thickness of 28.5m, the upper surrounding rock (green) resistivity of 200 

Ω·m, and the lower surrounding rock (yellow) resistivity of 500 Ω·m. The anomaly object ρ1 (blue) resistivity of 

1 Ω·m. Pass through two layers of media. This model simulates a low resistance anomaly object in the layered 

media space. The forward model and coordinate system are shown in Fig. 2(a): 

 

Figure 2: (a) Forward modeling for low resistance anomalous object in layered media 

 

Figure 2: (b) Forward modeling results of low resistance object based on conjugate gradient reconstruction 

algorithm 

 

Figure 2: (c) Forward modeling results of low resistance object based on Cholesky reconstruction algorithm 

As shown in Figs. 2(b) and 2(c), the imaging range and position of the anomaly objects in the layered medium 

are significantly different. In Fig. 2(b), the anomaly object imaging position is relatively close to the upper 

medium, while for anomaly objects in physical model, the imaging size seems smaller. In Fig. 2(c), the imaging 

1025



position of anomaly objects is relatively accurate, and the range is closer to that of the low resistance model. 

Comparing these two forward reconstruction algorithms, we can see that Cholesky's imaging position for the 

low-resistance object is relatively precise in the range closer to the size of the low-resistance model. 

In Fig. 2(c), the potential value calculated based on the Cholesky reconstruction algorithm, as prior information, 

is substituted into the model's inversion iteration, namely, the matrix equation is solved for the damped least  

7. Conclusion 

Finite difference method is used to explore the numerical simulation of low-resistance and high-resistance 

anomaly objects in homogeneous and layered media by using different reconstruction algorithms (conjugate 

gradient method and Cholesky). It is found that the Cholesky reconstruction algorithm works well on the numeral 

simulation for sounding low resistance object, while there is a little difference between two algorithms for high 

resistance object, such as air-raid shelters and holes. 

With the increase of the number of dissection meshes, the numerical simulation results based on the two 

reconstruction algorithms are more significantly different. The Cholesky has lower computational speed and 

requires less memory than the conjugate gradient method. The 3D space geoelectric simulation is more complex 

since there are many problems needed to be solved. The algorithm for the calculation program also has a high 

requirement on hardware. The work done hereof is only part of the whole. 

Reference 

Barker R., 1992, A simple algorithm for electrical imaging of the subsurface, First Break, 10(2), 53–62, DOI: 

10.3997/1365-2397.1992004 

LaBrecque D.J., Owen E., Dailey W., Tamirez A.L., 1992, Noise and Occam's inversion of resistivity tomography 

data, 62th SEG, Expanded Abstracts, 11(1), 1410, DOI: 10.1190/1.1822100 

Li Y.G., Oldenburg D.W., 1992, Approximate inverse mappings in DC resistivity problems, Geophys J Int, 109, 

343–362. 

Li Y.G., Oldenburg D.W., 1994, Inversion of 3-D DC resistivity data using an approximate inverse mapping, 

Geophys J Int,116(3), 527-537, DOI: 10.1111/j.1365-246X.1994.tb03277.x 

Loke M.H., Barker R.D., 1995, Least-squares deconvolution of apparent resistivity pseudosection, Geophysics, 

60(6), 1682-1690, DOI: 10.1190/1.1443900 

Murai T., 1985, Electrical impedance computed tomography based on a finite element model, IEEE Trans. 

Biomed Eng, BME-32(3), 654-712, DOI: 10.1109/TBME.1985.325526 

Noel M., Xu B.W., 1991, Arachaeological investgation by electrical resistivity tomography: a preliminary study, 

Geophys J Int, 107, 95-102. 

Price L.R., 1979, Electrical impedance computer tomography(ICT): a new CT imaging technique, Ibid, NS-26(2), 

2736-2739, DOI: 10.1109/TNS.1979.4330526 

Qian F.Y., Tanaka Y., Zhao Y.L., 1995, Resistivity tomography with surface electrode arrays and its application 

in ar-chaeology in Backwater region of Three Gorge Reservoir, Computerized Tomography Theory & 

Applications, 4(2), 40-50. 

Sasaki Y., 1989, Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data, 

Geophysics, 54(2), 254-262, DOI: 10.1190/1.1442649 

Sasaki Y., 1992, Resolution of resistivity tomography inferred from numerical simulation, Geophysical 

Prospecting, 40(4), 453-463, DOI: 10.1111/j.1365-2478.1992.tb00536.x 

Shima H., Sakayama T., 1987, Resistivity tomography: An approach to 2-D resistivity inverse problems, 57th 

SEG, Ex-panded Abstracts, 204–207. 

1026