Finite Element Approximation and Numerical Analysis of Three-dimensional Electrical Impedance Tomography

Electrical impedance tomography is solved by solving an inverse problem of elliptic equation, and a new numerical method or a new technique is argued to consider finite element (such as normal element and mixed element) in this paper on three dimensional region. Introducing different perturbations to boundary restrictions and using different spacial steps, the authors obtain numerical solutions and give comparison with exact solutions. Numerical data show that numerical solution can approximate exact solution well as spacial step taken small and the approximation of Neumann boundary condition is more stable than that of Dirichlet case. For Newton iterations on finite element method, a large-scaled system of massive linear equations is solved in each iteration, thus the computation is quite expensive. So two techniques are argued in the first half of this paper. Firstly, the invariance property of quasi-element stiffness matrix is used in the iterations and a type of special current model is introduced. Then the minimum number of direct problems solved is considered. Later a local conservative numerical approximation, low order mixed element (block-centered method) is presented in the latter part and the positive semi-definiteness and the existence of its solution are proved. Computational formula of error functional Jacobi matrix is derived and the least direct problems in each iteration are solved by using the symmetry of algorithm and a special current basis. This method has been applied successfully in actual numerical simulation of three-dimensional electrical impedance tomography.


Introduction
Electrical Impedance Tomography (EIT) is a new biomedicine imagining technique.Firstly some electrodes are attached on the skin and safe current is injected.Secondly, the values of potentials are measured and the distribution is obtained.Thirdly, the data are input in a computer and the distribution of electrical resistivity inside the body is illustrated by the monochrome graph after calculations and transformations.Since electrical resistivity is assigned by different values at different places, so the graph shows where and how the inner organs are.This technique is safe and economic, so EIT can be generalized in numerical simulation of more actual problems such as geophysical exploration, hydrogeology detection, dam body detection and underwater target detection.Therefore, this method is important and arouses more attention on theoretical research and its applications (Kirsch, 1996;Cheney, 1999;Webster, 1990;Tamburrino, 2002;Breckon, 1987;Du, 1997).
Electrical impedance tomography is essentially described by an inverse problem of elliptic partial differential equation, where u(x) denotes the potential, ρ > 0 is the electrical resistivity, ψ(x) is the boundary electric density and φ(x) denotes the boundary potential.
Given the electric resistivity ρ(x) and the boundary density ψ(x), then the potential distribution u(x) is obtained by (1a) and (1b).This is called Neumann boundary problem.Dirichlet boundary problem is interpreted by the fact that u(x) is computed by (1a) and (1c) as ρ(x) and φ(x) are given.The above two types are direct problems while EIT is an inverse problem.The functions ρ(x) and u(x) are unknown and ρ is an objective function as only two boundary value conditions ψ(x) and φ(x) are given.If ρ(x) is known, the potential distribution u(x) is determined by a proper boundary current density ψ(x) only different by a constant, then the boundary potential u(x)| ∂Ω = φ(x) is obtained.Therefore, the electrical resistivity distribution ρ(x) can define a mapping on the boundary current density set Ψ into the boundary potential set Φ(x): F ρ : ψ(x) → φ(x).
EIT shows how the electric resistivity ρ(x) is obtained under the mapping F ρ = F 0 satisfying It is hard to get the exact relation of ρ(x) and F ρ (direct problem), and the boundary current density ψ(x) is measured hardly.So this problem is solved by numerical approximations.
At present many numerical methods are discussed in solving this problem, and Newton iteration of normal finite method is a powerful tool.While Newton iteration gives rise to large-scaled computation especially for three-dimensional case.The method of finite element has high order of accuracy and has strong suitability in computational geometric regions.But the computation of three numerical integrals in generating coefficient matrices is generally quite expensive, where coefficient matrix is reformulated in each iteration to compute the gradient vector and Hesse matrix of objective function.Thus, in this paper we introduce the quasi-element stiffness matrix combined with normal finite element, because quasi-element has little-amounted and small-sized matrix fixed at each iteration.The element matrix has simple relation with coefficient matrix and the numerical integral computation cost is shortened.Another difficult of quasi-element stiffness matrix is to compute the gradient vector and Hesse matrix and to solve large-scaled linear equations (direct problem), so we introduce a special current forcing model.This uses sufficiently the potentials at all the nodes (solutions of direct problem), derives simply a Jacobi matrix of error vectors and finishes the computations of gradient vector and Hesse matrix.In the latter part the authors present a local conservative, low order mixed element scheme (block-centered finite volume element) (Russell, 1983;Weiser, 1988;Mishev, 1998;Lazarov, 1996;Larsson, 2003), and prove the scheme's positive semi-definiteness and the solution's existence.Computational formula is derived for Jacobi matrix of error functional and the number of direct problems is minimized at each iteration by the symmetry and special current basis vectors.Two different EIT models of continuous model and electrode model are considered.The model's correctness and the algorithm's reliability and feasibility are testified by exact solution simulation for continuous model, then actual algorithm is concluded for electrode model.
Numerical simulation is discussed in this paper for the inverse problem.On normal finite element method some discussion of iteration algorithm and numerical simulation is given in §2, and the algorithm's feasibility and the model's correctness are testified by comparing exact solutions on three dimensional domain (0, 1) 3 with numerical simulation data of different perturbations.In §3, considering actual simulation of EIT problem, the authors present quasi-element stiffness matrix technique to decrease numerical integral computation and discuss a handling method for Neumann boundary condition.
The projection F ρ is disretized and transformed into a mapping of currents and potentials at different electrodes, and a new current model and computational algorithm of Jacobi matrix are given.EIT problem on a cylinder is simulated numerically.In §4, a numerical method of low order mixed element (or named block-centered finite volume element) and its approximations are argued.This algorithm has the nature of mass conservation.In §5, image reconstruction and simulation experiments of mixed element method are considered.

Normal Finite Element Numerical Method and Simulation of Inverse Problem
An inverse problem is solved usually by the method of iterations.Firstly Neumann problem is solved by using normal finite element method.Then an ordinary least squares problem is given by combing its numerical solution and Dirichlet boundary condition.

Normal Finite Element Method of Neumann Boundary Problem
Neumann boundary problem of (1a) and ( 1b) is changed into a Galerkin variation as follows: to find u ∈ H 1 (Ω) such that where Ω is partitioned by a series of elements {K m } M m=1 whose vertexes make up a grid Ω (5) The quasi-element stiffness matrix is defined as follows.
Definition 1 Ā(m) = (ā (m) i j ) p×p is made up of ā(m) i j ordered by the node labels of K m , then Ā(m) is called the quasi-element stiffness matrix of the element K m .
The element stiffness matrix is defined by A (m) = (a (m)  i j ) p×p , where and it is easy to find the relation of A (m) and Ā(m) Noting that and letting [A (m) ] and [ Ā(m) ] of order n denote the expanded formulations of the matrices A (m) and Ā(m) of order p, then we have the following relation Similar to Neumann problem, the solutions of finite element equation (4) differ by constants, that is to say the solution exists solely by assigning the value of a node, u i 0 = ψ(x i 0 ).The equation (4) turns into an equivalent symmetric and positive definite system where A d is obtained by deleting the (i 0 , j) entries of the i 0 th row and the (i, i 0 ) entries of the i 0 th column from the matrix A, and the vectors u d , b d are given by deleting the i 0 th component from u, b.Given ρ h (x), the equation ( 11) is solved by U denotes the vector of potentials at all the boundary nodes regardless the reference node x i 0 , whose components are included in u d , then The (i, j) entries of matrix P are assigned either by the number 0 or by the number 1, corresponding to the reference nodes and the other nodes.

Iterations of Inverse Problem
Let U 0 denote the vector of the values of φ(x) at the boundary nodes except x i 0 , and ρ = (ρ 1 , ρ 2 , • • • , ρ M ) T is defined by the values of ρ h (x) at all the elements.Let r = U − U 0 , then we solve inverse problem by solving the following nonlinear optimization problem, where The gradient vector of objective function f (ρ) is given by where denotes the Jacobi matrix of r.
we find that it is adequate to only compute Using the technique in (Cheney, 1999), and noting that A −1 d A d = I, we get For A d is the remaining matrix of A deleted the i 0 th row and the i 0 th column, we get the expression of ∂A d ∂ρ m similarly from From ( 10), Substituting ( 18), ( 19) and ( 20) into ( 17), and noting that (12), we have Hesse matrix of objective function f (ρ) is defined by ≈ J T J, and Levenberg-Marquardt iteration is formulated as follows for solving (14), where I is an identity matrix.

Numerical Simulation and Data
In this subsection we give several experimental tests to show the successful applications.Take Ω = (0, 1) 3 , and let Γ i denote boundary surfaces, Γi .Each (0, 1) interval is divided into n 1 subintervals with equal step h = 1/n 1 , and Ω is divided into n 3 1 cubes.Basis node functions {N i (x)} are defined by three piecewise defined linear interpolation functions.Relative error function is defined by .

Experiment 1
Suppose that exact solution of (1a) is ρ(x) = u(x) = e x 1 +x 2 +x 3 , satisfying two boundary conditions (1b) and ( 1c), where To show how the initial data affect numerical results, we consider ψ(x) and φ(x) perturbed in the simulation by and give the perturbation effects in Table 2.
, and numerical data are illustrated in Table 3.The boundary conditions ψ(x) and φ(x) are perturbed by ( 23) and ( 24) (see Table 4).The functions ρ(x), u (1) (x), ψ (1) (x) are taken as ρ(x), u(x), ψ(x) of Experiment 2. u (2) (x) = x 1 x 2 x 3 is exact solution of (1a) corresponding to ρ(x), and Neumann boundary condition is defined by The objective function ( 14) is modified by ] T r (l) .Let r (l) = U (l) − U (l) 0 , where U (l) is defined by ( 13) and U (l)  0 is defined by the values at boundary nodes of u (l) (x).Gradient operator and Hesse matrix are formulated by ] T J (l) , where J (l) denotes Jacobi matrix of r (l) .Numerical data are shown in Table 5 and Table 6 after repeated simulations of Experiment 1 and Experiment 2. Notes.Since the problem has distinct semi-convergence, so numerical data in Table 5 and Table 6 denote optimal values of iteration computations.

Conclusions and Discussions
From numerical experiments, we conclude that numerical solution of inverse problem can approximate exact solution well as h approaches zero.However, the convergence rate is slow and the computational cost increases rapidly for threedimensional problems.The problem is ill-conditioned, while numerical solution of Neumann boundary condition is more stable than that of Dirichlet case.Comparing Experiment 2 with Experiment 3, we can develop the accuracy by considering boundary condition but don't increase supplemental work, because the computation cost on modifying the L-M parameter µ k decreases.The model ( 2) is testified correctly.For numerical simulation of inverse problems, at present there only has the primary regularization theory of Tikhonov and Lavnentiv (Borcea, 2003,Newman2000;Quarteroni, 2000;Stoer, 1993).Therefore, it is necessary to give a careful discussion on the convergence and stability analysis, and this is a major problem in theoretical research on inverse problem (Wang, 2007;Liu, 2005;Xiao, 2003).

Application of Normal Finite Element in Electrical Impedance Tomography
In actual imagining problem, some computation and measures should be solved as follows.

Quasi-element Stiffness Matrix Strategy
In the previous section, the iteration algorithm is formulated on finite element method, where massive three-dimensional numerical integrals are considered in expressing coefficient matrix, and this is most expensive.Furthermore, coefficient matrices are different in the iterations and give a large amount of additional computation.Thus we introduce quasi-element stiffness matrix and use constant function to approximate ρ(x).Quasi-element stiffness matrix only depends on the values of the node and basis function ( 5), and doesn't vary during the iteration.The square matrix of order p is described by (10), where p denotes the number of nodes in an element, and all the small-scaled matrices derived from a numerical integral calculation are stored and expanded into the coefficient matrix with ρ h (x).As shown in the previous numerical experiments, the domain is divided uniformly and all the quasi-element stiffness matrices have a same definition (5).So only an 8 × 8 square matrix is computed and stored, and this computation cost can be ignored.
Actual imaging region Ω is usually irregular, and it is partitioned into a series of rectangular subdomains and an irregular subdomain.Rectangular domain is divided equally into lots of small rectangles, and the element size determines a quasielement stiffness matrix only stored one time in the computation.Irregular domain is partitioned by isoparametric element method and all the quasi-element stiffness matrices are stored.For regular elements are more greatly than irregular elements, so the size and quality of stored quasi-element stiffness matrices are very small and their storage and computation are simple.
The quasi-element stiffness matrix of irregular element is determined as follows.In (Borcea, 2003), let K denote reference unit, and let node basis function N i (ξ) be defined by a three piecewise linear interpolation.Define B m is Jacobi matrix of isoparametric transformation of reference unit K to grid unit K m , and by ( 5), the quasi-element stiffness matrix Ā(m) is defined by

Neumann Boundary Condition Argument
Boundary current density ψ(x) (Neumann boundary condition) is immeasurable, so it is necessary to be properly handled.Boundary potential φ(x) (Dirichlet boundary condition) can be measured by U 0 only at the electrodes (partial boundary nodes).
Electrodes, whose total quantity is n 0 , are supposed to locate at boundary nodes (see Fig. 1), and each electrode only covers a node.The curved face covered by the kth electrode is denoted by γ k , where the node is x j k and the measurable current intensity is ψdS .Let I represent the current vector of I k .γ k , we make integrals on both sides of (1a) on Ω and use Gauss formula to get It is obviously seen that the input equals the output, and permissible current vector can be expanded to an n 0 −1 dimensional space.The right-hand side term of (4) is treated by where k i is the electrode number corresponding to the node x i , That is to say that , where b is an n dimensional vector expanded from the n 0 dimensional vector I according the node numbers.Therefore, ( 4) is changed into The right expression is given by the measured current directly, that is to say that the function ψ(x) and numerical integrals are not involved.
Similar analysis as in §2 is considered, Q d is the remaining matrix of Q deleted the i 0 th row.P is an (n 0 − 1) × (n 0 − 1) matrix made up of the numbers 0 and 1, and it can extract the potentials at the electrodes except the k i 0 electrode of u d .U denotes an n 0 − 1 dimensional vector consisting of the electrode potentials except the reference node.

Current Model and Simplified Calculation of Jacobi Matrix and Hesse Matrix
(34) shows that the relation of current vector and boundary potential vector (a discretization of F ρ ) is linear.Different boundary currents I (l) , l = 1, 2, • • • , L are imposed on Ω, where {I (l) } at least includes a basis.Measuring the potential at the boundary electrodes U (l) 0 , then we give an objective function where r (l) = U (l) − U (l) 0 and U (l) is defined by (34), L ≥ n 0 − 1.Similarly, where J (l) denotes the Jacobi matrix of r (l) .Adopting the expression of (31), then we reformulate ( 33) and (34) as follows where [I] d is an n − 1 dimensional vector by deleting the i 0 th element from [I].
Jacobi matrix, gradient vector and Hesse matrix are computed by ( 21), and the following equation system is solved The quantity is determined by the number of partitioned elements and the measurement times, and there at least n×(n 0 −1) systems are solved at each iteration for finite element problem.So its computation is greatly complicated.A technique of (Du, 1997) is introduced and improved.Noting that A d and (A −1 d ) T = (A T d ) −1 = A −1 d are symmetric and positive definite, we give another expression of ( 21), P T is an (n − 1) × (n 0 − 1) matrix, thus the calculation of A −1 d P T turns into solving (n 0 − 1) systems of equations.Moreover, the current model is arranged as follows.A unit current is input at the reference electrode (the i 0 th node) and is output at the others.The reference electrode is supposed to be labeled by the node n, and the current model is equivalent to the following current basis, Collecting the vectors as W = ([I (1) ], [I (2) ], • • • , [I (n 0 −1) ]), and deleting the i 0 th row, we have This shows the entries −1 are deleted from W, and W d = P T .Then, d computed previously is the solution of ( 32), and it is the potential vector of the lth time imposing current I (l) d except the reference node i 0 .Then ( 39) is equivalent to Therefore, it is adequate to solve n 0 − 1 finite element problems (systems of equations) at each iteration, and the computations of Jacobi matrix, gradient vector and Hesse matrix end out of the complicated formula (38).
A numerical experiment follows on a cylinder domain Ω = {x : The function ψ(x) has no analytical expression, and I (l) is known.Using this method above we can get numerical solutions of u(x) and boundary potential U (l) to replace measurement solution U (l)  0 .Initial approximation is taken by ρ m ≡ 1, and the electrical resistivity ρ is solved by L-M method.It aims to testify how the approximation ρ h (x) approaches the preestablished function ρ(x).The domain Ω is divided into 10 element strips from the top to the bottom, 300 units in each strip and the total number is 3000.The domain is separated by 11 node layers, 331 nodes in each layer, and the total number is 3641.In the simulation 261 electrodes are arranged and numerical data are shown in Table 7.For convenience, we take the section x 3 = 0.55 as an illustration, and the comparisons of true data and numerical simulations are shown in Fig. 2-Fig.5.All the numerical data are obtained by Matlab software in our personal computer.

Low Order Mixed Element Method
A local conservative low order mixed element method, or named block-centered finite volume method, is discussed in this section (Russell, 1983,Weiser1988;Mishev, 1998;Lazarov, 1996;Larsson, 2003).The domain Ω is enclosed by several Vol. 8, No. 4; 2016 faces parallel to the coordinate planes.The space R 3 is partitioned by lots of rectangular solids with the length, the weight and the height equal to (h 1 , h 2 , h 3 ), whose centers x i = (x i,1 , x i,2 , x i,3 ) form a grid ω = {x i } n i=1 , and the center of boundary block is supposed to lie in the boundary ∂Ω (see Fig. 6 and Fig. 7).The subset of cuboid and Ω is called control volume, in symbol V i .It holds obviously where γ i j is a rectangle and its area is denoted by m(γ i j ).Let ω = ω ∩ Ω, ∂ω = ω ∩ ∂Ω.Introduce the notations, Making integral on both sides of (1a) on any subset V ⊂ Ω, and using Gauss formula, we get the electric quantity balance equation, Taking V = V i , and letting then, we have As ∂V i ∩∂Ω ψ(x)dS = 0. Therefore, we get the algorithm of finite volume element, where Let then, we have the following statement Theorem 1 Suppose that the grid ω is connected and the coefficient matrix of ( 47) is denoted by A = (a i j ) n×n , then A is symmetric and positive semi-definite, the rank r(A) = n − 1, and any equation of ( 47) can be derived from the other n − 1 equations.
Proof: From (50) it follows clearly that A is symmetric.
y(x) and z(x) are supposed to be grid functions on ω, then by the symmetry of A we have Take z = y in the above expression, we have Then A is positive semi-definite.
Considering L h y = 0 and the above expression together, we have y j = y i , j ∈ Σ(i), x i ∈ ω, then by the connected property of ω we have y i ≡ C. Then the kernel space of A is one-dimensional, therefore r(A) = n − 1.

Making summation on
Making integrals on both sides of (1a) on Ω, and using Gauss formula and (1b), we have That is to say that any equation of ( 47) can be expressed by the other n − 1 equations.
From the proof of Theorem 1, we can get the following corollary.
Corollary 1: Any two solutions of ( 47) are different by a constant.
Given ρ(x), it is next to solve the equation (47

and express (47) in matrices
By Theorem 1 and Corollary 1, the potential at a node x i 0 (electrode) is taken as the reference potential u h,i 0 = 0, then the solution of (47) exists and is unique.Deleting the i 0 th equation (as a redundant equation), we have where the matrix A d is obtained by deleting the entries of the i 0 th row and the i 0 th column and it is symmetric and positive definite.u d and b d are obtained from u and b similarly.Its solution is U denotes a potential vector of all the boundary nodes except the reference node, and U is made up of the partial components of u d .Therefore, where the entries of P are assigned by the number 0 or by the number 1, and the values correspond to the potentials of the nodes except the reference node.

Mixed Element Iteration Algorithm of Inverse Problem
Let U (l) 0 denote the vector of the values of ϕ (l) at the boundary nodes except x i 0 , l = 1, 2, • • • , L. U is defined by ( 55), an approximation to U (l)  0 by using finite volume element, and ρ = (ρ 1 , ρ 2 , • • • , ρ n ) T is defined by the values of ρ h (x) at all the grids.Let r (l) = U (l) − U (l)  0 , then we solve inverse problem by solving the following nonlinear optimization problem, where Here we adopt the method of Levenberg-Marquardt, and we have to formulate the gradient vector and the Hesse matrix of objective function.

The gradient vector of
For convenience to discuss, we omit all the superscript (l).By the definition (55), it holds and it is adequate to only compute ∂ρ m = 0, then we have Then it follows from ( 58) and ( 54), From the definition of A d , we get that ∂A d ∂ρ m is derived from ∂A ∂ρ m by deleting the i 0 th row and the i 0 th column, Let ∂A ∂ρ m = (α m i j ) n×n , then Hesse matrix of objective function f (ρ) is defined by where The Levenberg-Marquardt iteration is formulated as follows to solve ( 56) where I is an identity matrix.

Mixed Finite Element Numerical Approximation of Inverse Problem
In this subsection we give several experimental tests to show the feasibility.Take Ω = (0, 0.8) × (0, 1) × (0, 1.2), and define discrete norms on ω as follows |v(x)| 0,h = ( ∑ where m(V i ) denotes the volume of V i , and relative error is
, and exact solution of (1a) and Neumann boundary condition are given as follows, Initial approximations are taken by ρ Notes.Since the problem is distinctly semi-convergent, so numerical data in Table 9 and Table 10 denote optimal values of iteration computations.

Conclusions and Discussions
From numerical experiments, we conclude that numerical solution of inverse problem can approximate exact solution well as h approaches zero.Meanwhile, the convergence rate is slow and the computational cost increases rapidly for three-dimensional problems.The problem is ill-posed, while numerical solution of Neumann boundary condition is more stable than that of Dirichlet case.For direct problems, low order mixed element (block-centered finite volume element) is argued as a powerful tool in numerical simulation of oil reservoir (Russell, 1995;Jones, 1995;Rui, 2012;Rui, 2013;Rui, 20122;Yuan, 2016;Yuan, 2015), which originates from the principle of mass conservation in physical science and is expressed in a discrete formulation.It has many advantages such as the simplicity and the small-scaled calculation of the scheme.The coefficient matrix is symmetric and positive definite and the problem is solved easily.Its numerical solution is convergent of second order accuracy.Therefore, the convergence and stability analysis of inverse problem should be paid more attention.

Mixed Element Simulation Experiment of Image Reconstruction
In actual image reconstruction, Neumann boundary condition ψ(x) denotes the boundary current density and is immeasurable.What we can measure is the current intensity through the electrodes, and it is just handled by finite volume element.The electrodes only cover partial boundary nodes, so the boundary potentials (boundary condition) we measure denote the potential values of partial boundary nodes.Generally, the computational region is not a cuboid.The simulation experiment of image reconstruction is discussed and the simple algorithm of Jacobi matrix is shown.
Later we discuss the relation of current vector J and the right-hand side terms of (47) or (52).By (49), we have for For x i ∈ ∂ω, there is not any electrode at x i , ∂V i ∩ ∂Ω ⊂ Γ 0 , then we have

Discussions
In this paper the highlights are concluded as follows.
(I) Electrical impedance tomography is a new technique in simulating biological medicine imagining technology.So this research has important theoretical and applicable values.
(II) Quasi-element stiffness matrix is applied in the iterations and has invariance property, so numerical integral computation of coefficient matrix is shortened greatly.
(III) A special current forcing model is put forward to compute gradient vector and Hesse matrix by using enough data.It is not necessary to consider large-scale linear equations.
(IV) Based on standard finite element method, a conservative low-order mixed finite element is established and its semipositive definiteness and existence are testified.
(V) From lots numerical simulations and data, we conclude that this method is a powerful, efficient and effective tool to solve some science problems such as geophysical exploration and hydrogeology.
and the discretized linear equations are formulated by A(ρ h )u = b, (4) where A = (a i j ) n×n and b = (b 1 , b 2 , • • • , b n ) T .Let N i (x) be the node basis function on i = 1, 2, • • • , n, and let the number of nodes of each element be denoted by p,

Figure
Figure 2. True distribution Figure 3. EIT simulation without perturbation

Figure
Figure 10.Actual electrical resistivity distribution Figure 11.Numerical distribution without perturbation

Table 1 .
Without perturbation and numerical data are shown in Table1.

Table 7 .
Simulation results on a cylinder

Table 9 .
Perturbation effects (ω = 10) and numerical data are shown in Table10after a similar process of the above experiment.

Table 11 .
Simulation data on stepped region