INTEGRAL FORMULATION OF THE COMPLETE ELECTRODE MODEL OF ELECTRICAL IMPEDANCE TOMOGRAPHY

,


1.
Introduction. Assume the conductivity distribution σ of a bounded domain Ω ⊂ R 3 is of interest. To infer σ, L electrodes are mounted on the boundary ∂Ω of the domain. Through these electrodes, direct current is either injected into or extracted out of the region. All electrodes are assumed to be perfect conductors. A thin resistive layer exists between each electrode and the boundary of the domain. Assume the amount of current through each electrode, the conductivity distribution within the region and all the thin resistive layers are known. The forward problem of electrical impedance tomography (EIT) is to compute the distribution of potential and current density within the region, and the voltages on each electrode. This paper proposes a new model for the forward problem of EIT described above. The model is based on the minimum energy principle [11,Section 9.6]. The new formulation characterizes the resistive layer underneath each electrode with the conductivity distribution within it; while the complete electrode model (CEM) uses the surface impedance. From the new formulation, we derive both the so-called shunt model and the CEM of EIT. We also propose a new method to compute numerically the forward problem of EIT. The new forward solver is formulated in terms of current. Its performance is tested and compared with the traditional finite element method (FEM), on a 2D EIT model. This paper is organized as follow. The rest of this section is a review of the standard CEM [1] and its weak form [9]. In Section 2.1, we introduce the integral formulation of EIT. We prove that it has a unique solution within appropriate function spaces in Section 2.2, and characterize its solution in Section 2.3. In Section 2.4, we show how the shunt model and CEM of EIT can be derived from the integral formulation. In Section 3, we propose a new numerical method for solving the forward problem of EIT, based on the integral formulation. In Section 4, we test the performance of the new forward solver on a 2D EIT model and compare it with that of the traditional FEM. Conclusions are made in the final section.
According to [9], the CEM for the EIT described above is Here u is the potential distribution over a bounded domain Ω with piecewise smooth boundary ∂Ω. The n is the unit out normal to ∂Ω. For i = 1, 2, ..., L, I i is the amount of current extracted out of the region through the i-th electrode; U i is the voltage at the i-th electrode; e i is the part of boundary ∂Ω of Ω that is covered by the i-th electrode; z i is contact impedance for the i-th electrode. The contact impedances are allowed to be spatially dependent [5,10]. Particularly, the choice for the shape of the contact conductance function affects the smoothness of the forward solution, and also the speed of convergence for any higher order numerical methods [6].
Following [9], we assume σ is continuously differentiable in Ω up to its boundary; the σ and {z i } L i=1 are all positive-valued; σ is bounded from above. Let U = (U 1 , U 2 , ..., U L ) T . The weak solution of boundary value problem 1 in The weak solution exists and is unique inḢ [9]. It is also the unique solution of (6 Here ς i is the conductivity within Ω i , i = 1, 2, ..., L. Like σ, they are assumed to be positive-valued and bounded from above. The J(j) is the rate of total Joule heat consumed over Ω .
For i = 1, 2, ..., L, e i is the surface of the i-th layer of contact impedance. They are assumed to have smooth boundary, and are relatively open with regard to ∂Ω . They also satisfy The n is the unit out normal to ∂Ω , and the j · n is the normal component of j on ∂Ω . According to [4,Theorem 2.5], the linear continuous mapping (11) γ n : j ∈ H(div; Ω ) −→ j · n ∈ H −1/2 (∂Ω ) is defined by a continuous extension on the mapping The continuous extension works because C ∞ (Ω ) 3 is dense in H(div; Ω ) with respect to the norm Theorem 2.4], and H(div; Ω ) is a Hilbert space for the norm 13. By its definition 9, H s is a closed subspace of H(div; Ω ), hence also a Hilbert space with respect to the norm 13. The definition 9 indicates that currents follow through the region Ω only through the electrodes, and the total amount of current flows into the region is the same as that of current flowing out of the region.
In definition 7c, χ e i is the characteristic function for the surface e i and (14) j · n , χ e i ∂Ω def = ∂Ω j · n χ e i ds = e i j · n ds, i = 1, 2, ..., L.
The constraint (15) j · n , χ e 1 ∂Ω − I 1 , j · n , χ e 2 ∂Ω − I 2 , . . . , j · n , χ e L ∂Ω − I L = 0 R L corresponds to that current of amount I i is extracted out of the region through the i-th electrode, i = 1, 2, ..., L. We have In constraint 7b with definition 7c, the corresponds to that no source of current exists within the region. The equation 18 is well-defined because of the following lemma.

Existence and uniqueness.
Theorem 2.2. The constrained minimization problem 7 has a unique solution.
Proof. First, we prove that is a closed convex subset of H s . Assume j 1 , j 2 ∈ H ss . Then for any α ∈ (0, 1), So αj 1 + (1 − α)j 2 ∈ H ss . Therefore H ss is a convex set. Assume {j n } ∈ H ss → j * , with respect to the norm 13. Then j * ∈ H s because H s is closed with respect to the norm 13. Since G is a linear continuous map, G (j * ) = 0 R L ⊗(L 2 (Ω )/R) . Therefore, j * ∈ H ss , and H ss is a closed subset of H s . Second, assume j n , j * ∈ H s . We have Since σ, ς i > 0, i = 1, 2, ..., L, the inequality 26 indicates that if j n → j * in H ss with respect to the norm 13, then Because of these, the generalized Weierstrass Existence Theorem e.g. [12,Theorem 2D] concludes that the constrained minimization problem 7 has a unique solution in H ss .

Solution characterization.
Theorem 2.3. The solution j of the constrained minimization problem 7 satisfies Here J (j) and G (j) is the Fréchet derivative of J and G at j ∈ H s , respectively. The R L ⊗ L 2 (Ω )/R * is the dual space of R L ⊗ L 2 (Ω )/R . The 1/2 in equation 28 ensures that multiplier λ can be interpreted as potential as seen later on.
Proof. The conclusion follows from the generalized Lagrange multiplier rule e.g. [12, p. 270] as long as The claim 30 is because with 1/C 1 being the minimum conductivity over Ω . The identity 33 is due to The claim 31 is because Finally, the claim 32 is because of identity 19.
2.4. Relation to traditional models. From the integral formulation 7, we derive both the shunt model and the CEM 1 of EIT. Under the framework of integral formulation, the CEM can be considered as a special example of the shunt model. To see these, we first prove the following theorem which implies that the Lagrange multiplier in condition 28 consists of potential distribution over the region and voltages on the electrodes.
for some u ∈ H 1 (Ω )/R and U i ∈ R, i = 1, 2, ..., L. Here Proof. First, we interpret the λ in condition 29 as Then condition 29 is equivalent to Second, we claim the u satisfying the constraint 41 should be in H 1 (Ω )/R. This is because the constraint 41 leads to 3 , the right-hand side of identity 42 is a linear functional on h which is considered as an element of L 2 (Ω ) 3 . So there exists g ∈ L 2 (Ω ) 3 such that The identity 43 implies that g = ∇u. So ∇u ∈ L 2 (Ω ) 3 and u ∈ H 1 (Ω )/R. Third, since u ∈ H 1 (Ω )/R, the Green's formula gives With definition 9 and the definition of Ω , the identity 44 becomes With identities 33 and 45, the constraint 41 is equivalent to Finally, the constraint 46 is equivalent to the constraint 38.
In terms of u in Theorem 2.4, all the constraints on the current density in equilibrium are The constraints 47 can be considered as a shunt model of EIT over the region Ω . The constraints 47a and 47b are due to the constraint 38a and the constraint 7b. The constraint 47c is because of constraint 38a and the definition 9 on H s . The constraint 47d is constraint 38b.
With Theorem 2.4, the CEM 1 of EIT can be derived from the integral formulation 7 as follow. The equation 1a is from equation (47a). The equation 1b is obtained from constraints 38a,7b and definition 9 as follows: (48) − ei σ∇u · n ds = ei j · n ds = e i j · n ds − ∂Ωi j · n ds = I i i = 1, 2, ..., L.
The equation 1c is because of constraint 47c and condition 8. Finally, the equation 1d is obtained when , on e i , i = 1, 2, ..., L.
Here Γ in the line integral in definition 49 is a path within Ω i that goes from the point x on the surface e i to an arbitrary point on the surface e i . The definition 49 indicates z i , i = 1, 2, ..., L is approximately proportional to the resistance of the i-th layer Ω i when Ω i is extremely thin. The integral formulation 7 is also closely related to the weak form 6 of CEM. With condition 38, the Lagrange function for the constrained minimization problem 7 is actually the negative of the minimization functional in problem 6 i.e.
This is because and 3. Algorithm. The integral formulation 7 suggests we solve the forward problem of EIT based on the current density j. To begin with, all the constraints on the solution j of the constrained problem 7 are summarized as follows: The equations 54a and 54b are from the constraint 7b. The equation 54c is because of the definition 9 of H s . The equation 54d is equivalent to equation 38a according to [4,Theorem 2.9]. Given conditions 38a, the equation 54e is equivalent to equation 38b i.e. u is constant on each e i , i = 1, 2, ..., L [4, Remark 2.3]. According to Theorem 2.4, the conditions 54d and 54e characterize the distribution of current density that minimizes the total Joule heat in equilibrium.
A new numerical method for solving the forward problem of EIT can be obtained by discretizing directly the constraints 54 on the current density j. Assume Ω is a polygon and partitioned by an admissible mesh (as defined e.g. in [2]) of M triangles. The current density is constant within each element. Let I i,j be the amount of current flowing out of the i-th triangle through its j-th edge, i = 1, 2, ..., M ; j = 1, 2, 3. In terms of I ij , the constraints 54 suggest I i l,s ,j = I l , l = 1, 2, ..., L, (55b) I ki,kj = 0, k = 1, 2, ..., n I , (55d) Here equation 55a corresponds to equation 54a. It says that the total amount of current flowing out of each triangle is zero.
The equation 55b corresponds to equation 54b. In equation 55b, the j -th edge of the i l,s -th triangle is the s-th line segment of the l-th electrode which consists of m l line segments.
The equation 55c corresponds to equation 54c. In equation 55c, the j b -th edge of the i b -th triangle is on the boundary where no electrode is covered.
The equation 55d corresponds to equation 54d. In equation 55d, the k-th node of all the n I nodes that are within the mesh, is one of the two vertices of the k j -th edge of the k i -th triangle. Summation is over all the edges that share the k-th internal node.
The equation 55e corresponds to equation 54e. In equation 55e, the i e -th triangle has an edge that belongs to some electrode. I ie,j1 and I ie,j2 are the amount of current flowing out of this triangle through its two other edges which do not belong to any electrodes. L j1 and L j2 are the projections of the j 1 -th edge and the j 2 -th edge of this triangle onto the edge that longs to the electrode.
Finally, the condition 55f is for current continuity. In equation 55f, the j-th edge of the i-th triangle is the j -th edge of the i -th triangle. The condition 55f should hold true for every edge within the mesh.
3.1. Implementations. The equations 55 suggest Algorithm 1 for solving the forward problem of EIT. In this algorithm, the system of linear equations 56 is over- 3. With {I i,j }, find the current density j which is assumed to be constant within each element. 4. Find the electric field over the region, given j and the conductivity distribution over the region. 5. Given the electric field, find the potential difference between any two vertices of each element. 6. Reconstruct the potential distribution from the potential difference between the two vertices of every element.
determined. This is because K has 3 · M columns and (57) Here E B is the number of edges on the boundary of the mesh; L is the number of electrodes; V I and E I are the number of nodes and edges within the mesh, respectively. The row number in 57 is larger than 3 · M because In step 6, a system of linear equations for the potentials at all the nodal points could be formed in the first place. Grounding is also required in order for the uniqueness of the solution.
4. Simulation results. We test the current-based forward solver on a 2D EIT model depicted in Figure 2. In this figure, two identical electrodes are mounted symmetrically on the boundary of a region. The angle subtended by each electrode is π/3. One ampere current is injected into and extracted out of the disk through the bottom and top electrode respectively. For simplicity, we specified that the extended region Ω , as defined in Section 2, is a disk of radius one meter. We also let the conductivity ς i , i = 1, 2, ..., L within the thin layer underneath each electrode be the same as the conductivity σ nearby within the region but outside the thin layer. The Figure 3 shows the underlying mesh for our numerical simulations. The mesh was generated by 'distmesh' [7] and consists of 4423 triangles.
We solved the forward problem of this EIT with both the traditional FEM solver and the new current-based solver. In FEM, linear triangle elements were employed.
When the region Ω has conductivity one Siemens per meter (S/m) uniformly, the potential distributions estimated by the current-based solver and the FEM solver (potential-based) were shown and compared in Figure 4. As seen they are  extremely close to each other, expect that the voltages on the electrode estimated by the current-based solver is a bit rough. As regard to the voltages on electrodes, it is which is obtained with conformal techniques e.g. [8](see also Appendix A), the current-based solver estimate slightly more accurate for this quantity than the traditional solver with linear elements, based on the same mesh. This, however, does not indicate that the new solver outperforms the standard FEM solver. Because the standard solver can use higher order elements to speed up its convergence for the CEM and get higher accuracy for the FEM solution [6].
When the region Ω has a nonuniform conductivity distribution as shown in Figure 5, the distributions of current density estimated by the two forward solvers were shown in Figure 6.  As seen, the distribution of current density estimated by the two methods were also close to each other. In both graphs, most of the currents flow around the two low conductivity region as expected.
In the simulation, both forward solvers assume that the potential distribution and the current density are constant within each element. The two solvers differ in the continuity constraints across the boundary of each element. The current-based solver requires the amount of current flowing across the neighboring element to be the same (i.e. condition 55f). This inevitably leaves the potential be discontinuous at each nodal point before the least square method is applied (i.e. step 6 in Algorithm 1). On the contrary, the traditional potential-based solver normally insists the potential distribution be continuous over the region, leaving the amount of current flowing across the neighboring elements be different.

5.
Conclusions. We found that EIT can be modeled based on the minimum energy principle. It results in a constrained minimization problem in terms of the distribution of current density. We proved that the new formulation has a unique solution within appropriate function spaces. We also derived both the shunt model and CEM of EIT from the new formulation, by characterizing its solution with the generalized Lagrange multiplier method. These derivations are consistent with the fact that the shunt model is the limit of the CEM in the special case when the contact resistances are zero [3]. The integral formulation also suggests computing the forward problem of EIT in terms of current density. We proposed a current-based forward solver of EIT. The new solver was shown to give similar results as those of the traditional FEM solver, for solving the forward problem of a 2D EIT model.
Further work is needed to study the consistency of the current-based solver, to improve its accuracy when it is generalized to solve the 3D forward problem. It is also important to test the new solver's accuracy with experimental data. The standard FEM solver usually can use higher order elements to improve the accuracy for its solution. So future work should explore the possibility for introducing 'higher order' numerical methods for the new current-based solver. In the η-plane, the set-up of EIT becomes that current of amount I is injected into the rectangular region ABCD through the side AB and extracted out of the region through side CD. Current flows vertically within the region and does not flow through either side AD or BC. The difference of potential on the side AB and CD is (63) U AB − U CD = cI 2bσ .