Electrodeless electrode model for electrical impedance tomography
EELECTRODELESS ELECTRODE MODEL FOR ELECTRICALIMPEDANCE TOMOGRAPHY
J. DARD´E † , N. HYV ¨ONEN ‡ , T. KUUTELA ‡ , AND
T. VALKONEN § Abstract.
Electrical impedance tomography is an imaging modality for extracting informationon the interior structure of a physical body from boundary measurements of current and voltage.This work studies a new robust way of modeling the contact electrodes used for driving currentpatterns into the examined object and measuring the resulting voltages. The idea is to not define theelectrodes as strict geometric objects on the measurement boundary, but only to assume approximateknowledge about their whereabouts and let a boundary admittivity function determine the actuallocations of the current inputs. Such an approach enables reconstructing the boundary admittivity,i.e. the locations and strengths of the contacts, at the same time and with analogous methods as theinterior admittivity. The functionality of the new model is verified by two-dimensional numericalexperiments based on water tank data.
Key words. electrical impedance tomography, electrode models, extended electrodes, varyingcontact admittance, Bayesian inversion
AMS subject classifications.
1. Introduction. In electrical impedance tomography (EIT) contact electrodesare used for driving currents into an examined body and for measuring the resultingvoltages, with the aim of reconstructing the interior admittivity of the body basedon such boundary measurements. We refer to the review articles [3, 6, 33] for moreinformation on EIT, including theoretical results on the unique identifiability of theinterior admittivity assuming idealized boundary measurements. The goal of thiswork is to introduce a new way of modeling the electrodes in EIT based on a surfaceadmittivity function on the boundary of the imaged object.The most accurate model for electrode measurements of EIT is arguably the complete electrode model (CEM) that accounts for the shapes of the electrodes as wellas for the thin resistive layers at the electrode-object interfaces [7, 32]. In its basicform the CEM assumes a constant contact on each electrode, but it also generalizesstraightforwardly to the case of spatially varying contacts. Other useful electrodemodels include the shunt model [7] that can be considered as the limit of the CEM asthe contacts become perfect [12] and the point electrode model [13] that is the limit ofthe CEM for relative measurements as the electrodes converge to points on the objectboundary. A common property of all these models is that they treat the electrodesas subsets of the object boundary.In this work, we do not assume to know the precise locations of the electrodes butonly certain larger boundary sections, called extended electrodes (cf. [14]), that corre-spond to the available prior information on the whereabouts of the true electrodes. † Institut de Math´ematiques de Toulouse, UMR 5219, Universit´e de Toulouse, CNRS, UPS, F-31062 Toulouse Cedex 9, France ([email protected]). The work of JD was sup-ported by the Institut Fran¸cais de Finlande, the Embassy of France in Finland, the French Ministryof Higher Education, Research and Innovation, the Finnish Society of Sciences and Letters and theFinnish Academy of Science and Letters (2019 Maupertuis Programme). ‡ Aalto University, Department of Mathematics and Systems Analysis, P.O. Box 11100, FI-00076Aalto, Finland (nuutti.hyvonen@aalto.fi, topi.kuutela@aalto.fi). The work of NH and TK was sup-ported by the Academy of Finland (decision 312124). § Department of Mathematics and Statistics, University of Helsinki, Finland and ModeMat, Es-cuela Polit´ecnica Nacional, Quito, Ecuador (tuomo.valkonen@iki.fi). The work of TV was supportedby the Academy of Finland (decisions 314701 and 320022).1 a r X i v : . [ m a t h . NA ] F e b J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN
For each true electrode, there is exactly one extended electrode. The prior informa-tion on the positions of the true electrodes provided by the extended electrodes iscomplemented by allowing the contact admittance to vary spatially or even vanishon some parts of the extended electrodes. The spatial dependence of the contacts isencoded in the so-called boundary admittivity function. In particular, if the boundaryadmittivity is smooth over the whole object boundary, there are no singularities atthe edges of the electrodes (cf. [17]), which facilitates efficient numerical solution ofassociated forward problems.We consider two alternative parametrizations for the boundary admittivity inour two-dimensional numerical experiments: the first one is directly based on thepiecewise linear basis functions for the underlying finite element (FE) discretization,and the other corresponds to hat-shaped surface admittivities with variable positions,widths and heights on the extended electrodes (cf. [17]). We demonstrate that bothof these choices can be used for incorporating the estimation of the boundary admit-tivity in a certain iterative Bayesian output least squares reconstruction algorithm forEIT, leading to more accurate reconstructions of the main unknown, i.e. the interioradmittivity, than employing the standard CEM with inaccurate positions for the trueelectrodes. See [2, 4, 23] for more information on the artifacts caused by geometricmismodeling in EIT.Incorporating the estimation of the electrode positions in a reconstruction algo-rithm of EIT has previously been considered in [5, 9, 10, 11, 16]. These works arebased on computing (shape) derivatives of the measurements of EIT with respect tothe positions of electrodes with known shapes, not on modeling the uncertainty inthe measurement configuration via a boundary admittivity function. Compared to[5, 9, 10, 11, 16], an attractive aspect of our approach is that estimating the effectivelocations of the current inputs is conceptually similar to the reconstruction of theinterior admittivity, that is, one aims at reconstructing a spatially varying functiondefined on the employed FE mesh. For the sake of completeness, it should also bementioned that EIT with imprecisely known measurement setup has previously beenconsidered in [15, 21, 22, 27, 28] to name a few approaches, but none of these orother works on the topic have tackled the uncertainty in the electrode positions viaemploying a spatially varying boundary admittivity function. Moreover, it is wellknown that (time or frequency) difference imaging is a partial remedy for geometricmismodeling [1], but in this work we are solely working with absolute EIT measure-ments.This text is organized as follows. Section 2 introduces the new “electrodelesselectrode model” and tackles its Fr´echet differentiability with respect to the inte-rior and boundary admittivities. Section 3 describes our iterative Bayesian outputleast squares reconstruction algorithm, considers the two aforementioned alternativeparametrizations for the boundary admittivity function and presents a partial con-vergence result. In Section 4, the numerical tests based on experimental water tankdata are presented. Finally, Section 5 lists the concluding remarks.
2. New electrode model for EIT.
In this section, we first introduce our newelectrode model by generalizing the standard CEM and prove its unique solvability.Subsequently, the Fr´echet differentiability of the associated forward solution withrespect to the boundary and domain admittivities is considered.
Let the bounded Lipschitzdomain Ω ⊂ R n , n = 2 or 3, model the physical body that is under investigation byEIT. Its boundary ∂ Ω is partially covered by M ∈ N \{ } electrodes { e m } Mm =1 that are LECTRODELESS ELECTRODE MODEL ∂ Ω that they cover. However, aswe do not assume to know the precise locations of the actual electrodes, we introduceso-called extended electrodes { E m } Mm =1 that have the following properties:(i) For all m = 1 , . . . , M , e m ⊂ E m ⊂ ∂ Ω.(ii) The relatively open, connected extended electrodes are mutually disjoint: E m ∩ E l = ∅ for all m (cid:54) = l .In other words, { E m } Mm =1 characterize our prior information on the whereabouts of thetrue electrodes { e m } Mm =1 . Take note that we do not exclude the case E := ∪ E m = ∂ Ωat this stage.The contacts at the (extended) electrodes are characterized by the surface admit-tivity ζ : ∂ Ω → C . Since we must be able to feed current through all electrodes andcontact conductances with negative real parts are unphysical, the set of admissibleadmittivities is defined as Z := (cid:8) ζ ∈ L ∞ ( E ) (cid:12)(cid:12) Re( ζ ) ≥ ζ | E m ) (cid:54)≡ m = 1 , . . . , M (cid:9) ⊂ L ∞ ( ∂ Ω) , (2.1)where the conditions on ζ are to be understood in the topology of L ∞ ( E ). Hereand in what follows, Z and L ∞ ( E ) are always interpreted as subsets of L ∞ ( ∂ Ω)via zero continuation. On the other hand, the electromagnetic properties of Ω arecharacterized by an isotropic admittivity σ : Ω → C that belongs to S = (cid:8) σ ∈ L ∞ (Ω) (cid:12)(cid:12) ess inf(Re( σ )) > (cid:9) ⊂ L ∞ (Ω) . (2.2)In other words, the essential infimum of the real part of the conductivity inside Ω islarger than zero, meaning that current can flow anywhere within Ω.Suppose the net currents I m ∈ C , m = 1 , . . . , M , are driven through the cor-responding electrodes and the resulting constant electrode potentials U m ∈ C , m =1 , . . . , M , are measured. Due to conservation of charge, any physically reasonablecurrent pattern I = [ I , . . . , I M ] T belongs to the subspace C M (cid:5) := (cid:110) J ∈ C M (cid:12)(cid:12)(cid:12) M (cid:88) m =1 J m = 0 (cid:111) . To facilitate writing down the boundary condition for the CEM, U = [ U , . . . , U M ] T ∈ C M is identified with U = M (cid:88) m =1 U m χ m , (2.3)where χ m is the characteristic function of the m th extended electrode E m ⊂ ∂ Ω.Whether U refers to a piecewise constant function supported on E or a vector in C M should be clear from the context.We mimic the formulation of the CEM presented in [17]: The electromagneticpotential u inside Ω and the piecewise constant electrode potential U jointly satisfy ∇ · ( σ ∇ u ) = 0 in Ω ,ν · σ ∇ u = ζ ( U − u ) on ∂ Ω , (cid:90) E m ν · σ ∇ u d S = I m , m = 1 , . . . , M, (2.4) J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN interpreted in the weak sense and with ν ∈ L ∞ ( ∂ Ω , R n ) denoting the exterior unitnormal of ∂ Ω. Unlike in [17], the true electrodes do not appear in the forward problem(2.4) that is formulated using the extended electrodes. However, if the positions ofthe true electrodes were known, then one could return to the setting analyzed in [17]by simply defining ζ = 0 in E m \ e m for all m = 1 , . . . , M . The connection betweenthe original formulation of the CEM in [7, 32] and (2.4) is discussed in [17]; to put itshort, the contact impedances in the traditional formulation of the CEM are obtainedas z m := (1 /ζ ) | E m , m = 1 , . . . , M , if e m = E m and spatially varying contacts, whichmay have singularities, are allowed.We seek the solution of (2.4) in the quotient space H defined via H s := (cid:8) { ( v + c, V + c ) | c ∈ C } (cid:12)(cid:12) ( v, V ) ∈ H s (Ω) ⊕ C M (cid:9) , s ∈ R , with := [1 , . . . , T ∈ C M . In other words, all elements of H s (Ω) ⊕ C M that differby an additive constant are identified as an equivalence class — this corresponds tothe freedom in the choice of the common ground level of potential for the interior andelectrode potentials. In particular, when the second component of an element in H s isinterpreted as a piecewise constant function on the electrodes, the additive constantis also supported on E . The standard quotient norm for H s is defined as (cid:107) ( v, V ) (cid:107) H s := inf c ∈ C (cid:16) (cid:107) v − c (cid:107) H s (Ω) + (cid:107) V − c (cid:107) (cid:17) / , (2.5)where (cid:107) · (cid:107) denotes the Euclidean norm. We use the notation (cid:107) · (cid:107) C M / C for theassociated quotient norm of C M defined by dropping out the first term on the right-hand side of (2.5).As can easily be deduced based on material in [17, 32], the variational formulationof (2.4) amounts to finding ( u, U ) ∈ H such that B σ,ζ (cid:0) ( u, U ) , ( v, V ) (cid:1) = I · V for all ( v, V ) ∈ H , (2.6)where the sesquilinear form B σ,ζ : H × H → C is defined by B σ,ζ (cid:0) ( w, W ) , ( v, V ) (cid:1) = (cid:90) Ω σ ∇ w · ∇ v d x + (cid:90) ∂ Ω ζ ( W − w )( V − v ) d S, (2.7)with W, V ∈ C M identified with the corresponding piecewise constant functions. If σ ∈ S and ζ ∈ Z are admissible, the sesquilinear form B σ,ζ : H × H → C is boundedand coercive: (cid:12)(cid:12) B σ,ζ (cid:0) ( w, W ) , ( v, V ) (cid:1)(cid:12)(cid:12) ≤ C (cid:107) ( w, W ) (cid:107) H (cid:107) ( v, V ) (cid:107) H , (2.8)Re (cid:16) B σ,ζ (cid:0) ( v, V ) , ( v, V ) (cid:1)(cid:17) ≥ c (cid:107) ( v, V ) (cid:107) H , (2.9)where c = c ( σ, ζ, Ω , E ) > C = ( σ, ζ, Ω , E ) > w, W ) , ( v, V ) ∈H . This conclusion follows by repeating the argumentation in the proof of [17,Lemma 2.1], and it also leads to the following theorem that is essentially a restatementof [17, Theorem 2.2] with the extended electrodes playing here the role of the actualelectrodes in [17]. Theorem 2.1. If σ ∈ S and ζ ∈ Z , the problem (2.4) has a unique solution ( u, U ) ∈ H for any current pattern I ∈ C M (cid:5) . Moreover, (cid:107) ( u, U ) (cid:107) H ≤ C (cid:107) I (cid:107) , LECTRODELESS ELECTRODE MODEL where C = C (Ω , E, σ, ζ ) > is independent of I . If the extended electrodes are well separated, the regularity of u is dictated bythe smoothness of ζ ∈ Z ⊂ L ∞ (Ω) interpreted via zero continuation as a functionon the whole of ∂ Ω together with the regularity of the interior admittivity σ and theboundary ∂ Ω. In particular, smooth parametrizations for ζ and σ facilitate efficientnumerical solution of (2.4), say, by resorting to higher order FEs [17]. For simplicity,the following theorem is stated assuming ∂ Ω and σ are infinitely smooth; the effectthat irregularities in σ and ∂ Ω have on the regularity of u could be straightforwardlyanalyzed based on standard elliptic theory (see, e.g., [26]). Theorem 2.2.
Assume that σ ∈ S ∩ C ∞ (Ω) , ∂ Ω is of class C ∞ and E m ∩ E l = ∅ for m (cid:54) = l . If ζ ∈ Z ∩ H s ( ∂ Ω) for some s > ( n − / , then the solution to (2.6) satisfies (cid:107) ( u, U ) (cid:107) H s +3 / ≤ C (cid:107) I (cid:107) , (2.10) where C = C (Ω , E, σ, ζ, s ) > is independent of I ∈ C M (cid:5) .Proof . The claim is a consequence of standard elliptic regularity theory, and itcan be deduced by following the line of reasoning leading to [17, Theorem 2.5]. Remark 2.3.
The essential ingredient of Theorem 2.2 is that the smooth behaviorof ζ compensates to a certain extent for the singularities caused by the piecewiseconstant nature of U on the right-hand side of the second condition in (2.4) . Thisleads to the possibility to iterate a bootstrap argument that the Neumann trace of u has the same Sobolev regularity as its Dirichlet trace, which in turn leads to higherregularity for u itself up to a certain limit dictated by ζ . It seems intuitive thatone could use such an argument also in the case when the closures of the extendedelectrodes are allowed to touch if it is required that ζ goes smoothly enough to zero attheir common boundaries. However, as the case of touching (extended) electrodes wasnot considered in [17], we do not stress this matter any further to avoid rewriting thecorresponding proofs. Before consider-ing the differentiability of the second part of the solution to (2.4) with respect to ζ , itis natural to extend the set of admissible boundary conductivities so that it becomesan open subset of L ∞ ( E ). To this end, we simply define Z (cid:48) = (cid:8) ζ ∈ L ∞ ( E ) (cid:12)(cid:12) For any σ ∈ S , the condition (2.9) holds with c = c ( σ, ζ ) > (cid:9) , (2.11)which is, as usual, interpreted via zero continuation as a subset of L ∞ ( ∂ Ω). Sinceit is clear that (2.8) actually holds for any σ ∈ L ∞ (Ω) and ζ ∈ L ∞ ( E ), it followsimmediately from the Lax–Milgram lemma that Theorem 2.1 remains valid if Z isreplaced by Z (cid:48) in its assertion. Moreover, Z (cid:48) is indeed open as a subset of L ∞ ( E ),which is a byproduct of the following more general lemma. Lemma 2.4.
For any ( σ, ζ ) ∈ S × Z (cid:48) , there exists δ = δ ( σ, ζ, Ω , E ) > such that Re (cid:16) B σ + ω,ζ + η (cid:0) ( v, V ) , ( v, V ) (cid:1)(cid:17) ≥ c ( σ, ζ )2 (cid:107) ( v, V ) (cid:107) H ∀ ( v, V ) ∈ H for all ( ω, η ) ∈ L ∞ (Ω) × L ∞ ( E ) satisfying (cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) ≤ δ . In particular, Z (cid:48) is an open subset of L ∞ ( E ) . J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN
Proof . Consider arbitrary ( σ, ζ ) ∈ S × Z (cid:48) and ( ω, η ) ∈ L ∞ (Ω) × L ∞ ( E ). For any( v, V ) ∈ H ,Re (cid:16) B σ + ω,ζ + η (cid:0) ( v, V ) , ( v, V ) (cid:1)(cid:17) = Re (cid:16) B σ,ζ (cid:0) ( v, V ) , ( v, V ) (cid:1)(cid:17) + (cid:90) Ω Re( ω ) |∇ v | d x + (cid:90) ∂ Ω Re( η ) | V − v | d S ≥ (cid:0) c ( σ, ζ ) − (cid:107) ω (cid:107) L ∞ (Ω) − K (cid:107) η (cid:107) L ∞ ( E ) (cid:1) (cid:107) ( v, V ) (cid:107) H , where we used (2.11), K = K (Ω , E ) >
0, and the estimate for the integral over ∂ Ω fol-lows from the same reasoning as [17, Lemma 2.4]. Choosing δ = c ( σ, ζ ) / (4 max { , K } )completes the proof.In particular, the forward problem (2.4) remains uniquely solvable if ζ ∈ Z ⊂ Z (cid:48) is replaced by ζ + η , with a small enough 0 (cid:54) = η ∈ L ∞ ( E ), even if ζ + η is negative inits real part on some subsets of E (cf. (2.1)). Theelectrode potentials, i.e. the second part of the solution to (2.4), can obviously beconsidered as a function of three variables: U : (cid:40) ( σ, ζ, I ) (cid:55)→ U ( σ, ζ, I ) , S × Z (cid:48) × C M (cid:5) → C M (cid:5) , (2.12)where we have systematically selected the ground level of potential so that U hasalways zero mean. Since our choice for the (extended) set of admissible boundaryconductivities Z (cid:48) is nonstandard (cf. [18]), we need to verify some basic properties ofthe measurement map (2.12). Lemma 2.5.
The mapping U : S × Z (cid:48) × C M (cid:5) → C M (cid:5) is continuous.Proof . Let ( σ, ζ, I ) ∈ S × Z (cid:48) × C M (cid:5) be arbitrary and the first two componentsof ( ω, η, J ) ∈ L ∞ (Ω) × L ∞ ( E ) × C M (cid:5) small enough so that also ( σ + ω, ζ + η, I + J ) ∈ S × Z (cid:48) × C M (cid:5) ; see Lemma 2.4 and (2.2). Let ( u, U ) and ( u ω,η,J , U ω,η,J ) be theunique solutions of (2.4) for the parameter triplets ( σ, ζ, I ) and ( σ + ω, ζ + η, I + J ),respectively, and define ( w ω,η,J , W ω,η,J ) = ( u ω,η,J − u, U ω,η,J − U ). By subtractingthe corresponding variational formulations (2.7), one easily deduces that B σ + ω,ζ + η (cid:0) ( w ω,η,J , W ω,η,J ) , ( v, V ) (cid:1) = J · V − (cid:90) Ω ω ∇ u · ∇ v d x − (cid:90) ∂ Ω η ( U − u )( V − v ) d S (2.13)for all ( v, V ) ∈ H .Assume that (cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) ≤ δ , with δ = δ ( σ, ζ ) > v, V ) = ( w ω,η,J , W ω,η,J ) in (2.13) and employing Lemma 2.4 yields (cid:107) ( w ω,η,J ,W ω,η,J ) (cid:107) H ≤ c ( σ, ζ ) (cid:16) (cid:107) J (cid:107) (cid:107) W ω,η,J (cid:107) C M / C + (cid:90) Ω | ω ||∇ u · ∇ w ω,η,J | d x + (cid:90) ∂ Ω | η || U − u || W ω,η,J − w ω,η,J | d S (cid:17) ≤ C (cid:48) (cid:0) (cid:107) J (cid:107) + max {(cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) }(cid:107) ( u, U ) (cid:107) H (cid:1) (cid:107) ( w ω,η,J , W ω,η,J ) (cid:107) H , where C (cid:48) = C (cid:48) ( σ, ζ, Ω , E ) > ∂ Ω was handled by resorting tosimilar means as in the proof of Lemma 2.4. Dividing by (cid:107) ( w ω,η,J , W ω,η,J ) (cid:107) H and LECTRODELESS ELECTRODE MODEL (cid:107) W ω,η,J (cid:107) C M / C ≤ (cid:107) ( w ω,η,J , W ω,η,J ) (cid:107) H ≤ C (cid:48)(cid:48) (cid:0) (cid:107) J (cid:107) + max {(cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) }(cid:107) I (cid:107) (cid:1) , (2.14)where C (cid:48)(cid:48) = C (cid:48)(cid:48) ( σ, ζ, Ω , E ) >
0. Because the quotient norm (cid:107) · (cid:107) C M / C is equivalent to (cid:107) · (cid:107) on C M (cid:5) , the proof is complete.Let us next prove the differentiability of U : S × Z (cid:48) × C M (cid:5) → C M (cid:5) ; such a resultis of course well known for the standard formulation of the CEM (see, e.g., [18, 25]),but our extended definition for the admissible contact conductivities arguably merits(re)writing a formal proof. To this end, let us introduce an auxiliary variationalproblem of finding ( u (cid:48) , U (cid:48) ) = ( u (cid:48) ( σ, ζ, I ; ω, η, J ) , U (cid:48) ( σ, ζ, I ; ω, η, J )) ∈ H such that B σ,ζ (cid:0) ( u (cid:48) , U (cid:48) ) , ( v, V ) (cid:1) = J · V − (cid:90) Ω ω ∇ u ·∇ v d x − (cid:90) ∂ Ω η ( U − u )( V − v ) d S ∀ ( v, V ) ∈ H , (2.15)where σ ∈ S , ζ ∈ Z (cid:48) , J ∈ C M (cid:5) , ω ∈ L ∞ (Ω), η ∈ L ∞ ( E ) ⊂ L ∞ ( ∂ Ω), and ( u, U ) =( u ( σ, ζ, I ) , U ( σ, ζ, I )) ∈ H is the solution to (2.4). Since it is straightforward to checkthat the right-hand side of (2.15) satisfies (cf., e.g., [17]) (cid:12)(cid:12)(cid:12) J · V − (cid:90) Ω ω ∇ u · ∇ v d x − (cid:90) ∂ Ω η ( U − u )( V − v ) d S (cid:12)(cid:12)(cid:12) ≤ C (cid:0) (cid:107) J (cid:107) + max {(cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) }(cid:107) ( u, U ) (cid:107) H (cid:1) (cid:107) ( v, V ) (cid:107) H , it follows from the Lax–Milgram theorem that (2.15) has a unique solution that sat-isfies (cid:107) ( u (cid:48) , U (cid:48) ) (cid:107) H ≤ C (cid:0) (cid:107) J (cid:107) + max {(cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) }(cid:107) ( u, U ) (cid:107) H (cid:1) ≤ C (cid:48) (cid:0) (cid:107) J (cid:107) + max {(cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) }(cid:107) I (cid:107) (cid:1) , where C = C (Ω , E, σ, ζ ) and C (cid:48) = C (cid:48) (Ω , E, σ, ζ ) are positive constants. Moreover,the pair ( u (cid:48) , U (cid:48) ) depends linearly on ( J, ω, η ) as does the right-hand side of (2.15).
Theorem 2.6.
The mapping
S × Z (cid:48) × C M (cid:5) (cid:51) ( σ, ζ, I ) (cid:55)→ U ( σ, ζ, I ) ∈ C M (cid:5) is Fr´echet differentiable. To be more precise, for any ( σ, ζ, I ) ∈ S × Z (cid:48) × C M (cid:5) , (cid:13)(cid:13) U ( σ + ω, ζ + η, I + J ) − U ( σ, ζ, I ) − U (cid:48) ( σ, ζ, I ; ω, η, J ) (cid:13)(cid:13) = O (cid:0) max {(cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) , (cid:107) J (cid:107) } (cid:1) , ( ω, η, J ) ∈ L ∞ (Ω) × L ∞ ( E ) × C M (cid:5) , where the Fr´echet derivative U (cid:48) ( σ, ζ, I ; ω, η, J ) ∈ C M (cid:5) is the second part of the solutionto (2.15) . Moreover, this derivative can be assembled through the relation U (cid:48) ( σ, ζ, I ; ω, η, J ) · ˜ I = J · ˜ U − (cid:90) Ω ω ∇ u · ∇ ˜ u d x − (cid:90) ∂ Ω η ( U − u )( ˜ U − ˜ u ) d S, (2.16) where ( U, u ) , ( ˜ U , ˜ u ) ∈ H are the solutions for (2.4) at the common parameters ( σ, ζ ) ∈S × Z (cid:48) but for the current patterns I, ˜ I ∈ C M (cid:5) , respectively. J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN
Proof . Let us adopt the notation introduced in the proof of Lemma 2.5. Sub-tracting (2.15) from (2.13) results in B σ,ζ (cid:0) ( w ω,η,J − u (cid:48) ,W ω,η,J − U (cid:48) ) , ( v, V ) (cid:1) = − (cid:90) Ω ω ∇ w ω,η,J · ∇ v d x − (cid:90) ∂ Ω η ( W ω,η,J − w ω,η,J )( V − v ) d S for all ( v, V ) ∈ H . Setting ( v, V ) = ( w ω,η,J − u (cid:48) , W ω,η,J − U (cid:48) ) and employing (2.9)(cf. (2.11)), it is straightforward to deduce that (cid:107) ( w ω,η,J − u (cid:48) , W ω,η,J − U (cid:48) ) (cid:107) H ≤ C max {(cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) } (cid:107) ( w ω,η,J , W ω,η,J ) (cid:107) H × (cid:107) ( w ω,η,J − u (cid:48) , W ω,η,J − U (cid:48) ) (cid:107) H , where C = C ( σ, ζ, Ω , E ) >
0. After dividing by (cid:107) ( w ω,η,J − u (cid:48) , W ω,η,J − U (cid:48) ) (cid:107) H , theclaim on the Fr´echet differentiability follows from the second inequality in (2.14) underthe assumption that (cid:107) ω (cid:107) L ∞ (Ω) , (cid:107) η (cid:107) L ∞ ( E ) < δ , with δ = δ ( σ, ζ ) > U , ˜ u ) of the solution to the original problem (2.6) for ˜ I ∈ C M (cid:5) as the test functionpair in (2.15) and subsequently equating the left-hand sides of (2.6) and (2.15). Remark 2.7.
By the chain rule for Banach spaces (see, e.g., [8, Theorem 7.1-3]), the representation formula (2.16) also extends to an arbitrary parametrization of ζ as ∂ θ U ( σ, ζ ( θ ) , I ; ϑ ) · ˜ I = − (cid:90) ∂ Ω ∂ θ ζ ( θ ; ϑ ) ( U − u )( ˜ U − ˜ u ) d S, θ ∈ D, ϑ ∈ B , assuming the mapping B ⊃ D (cid:51) θ (cid:55)→ ζ ( θ ) ∈ Z (cid:48) ⊂ L ∞ ( E ) is Fr´echet differentiable,with B being a Banach space and D its open subset. A similar conclusion also appliesto parametrizations of the domain conductivity.
3. Bayesian reconstruction algorithm.
In this section we briefly describeour numerical forward solver and formulate a Gauss–Newton type reconstruction al-gorithm for simultaneous reconstruction of the isotropic domain admittivity σ andthe boundary admittivity ζ . Since both σ and ζ are modeled as real-valued in ournumerical experiments, they will be referred to as the domain and surface conductiv-ities in what follows. Similar algorithms have been previously employed for relatedproblems in, e.g., [16, 34].To facilitate a reasonable numerical implementation, we assume that the domainΩ ⊂ R is simply connected and sufficiently smooth, and that it can be easily tri-angulated. Furthermore, we assume that the boundary subsections corresponding tothe extended electrodes { E m } Mm =1 ⊂ ∂ Ω are well separated and connected. However,we still only assume that e m ⊂ E m for each m = 1 , . . . , M , i.e., that each extendedelectrode covers the corresponding true electrode, but it does not usually hold that e m = E m . Finally, we assume that the extended electrodes are enumerated andparametrized along the boundary curve ∂ Ω.Approximate solutions for (2.4) are computed using a finite element method (FEM). In our FEM solver, the domain Ω is discretized into a triangle mesh, denotedby Ω h , on which the domain and boundary potentials are solved as linear combina-tions of piecewise linear basis functions. The domain conductivity σ ∈ L ∞ + (Ω h , R )is discretized with respect to the same piecewise linear basis as the solution. In-spired by the observations in [18], our algorithm actually aims at reconstructing the LECTRODELESS ELECTRODE MODEL κ := log σ ∈ L ∞ (Ω h , R ) instead of the actualconductivity. The derivatives of the forward solutions with respect κ can be com-puted by resorting to (2.16) and the chain rule for Banach spaces; see Remark 2.7and [18, Eq. (9)].The piecewise linear boundary curve ∂ Ω h naturally also determines the edges onwhich the extended electrodes are defined. In order to avoid current leaking directlybetween electrodes, we require at least one boundary node (and the correspondingbasis function) in-between any adjacent extended electrodes on ∂ Ω h . The contactconductivity ζ is parametrized on the piecewise linear boundary curve of the mesh.We formulate two different parametrizations for ζ in Section 3.1, denoted generallyby a finite-dimensional real map θ (cid:55)→ ζ ( θ ).In order to simplify the notation, we define concatenated vectors for voltage mea-surements, current injections and forward model solutions, respectively, as V = (cid:2) ( V (1) ) (cid:62) , . . . , ( V ( M − ) (cid:62) (cid:3) (cid:62) ∈ R M ( M − , I = (cid:2) ( I (1) ) (cid:62) , . . . , ( I ( M − ) (cid:62) (cid:3) (cid:62) ∈ R M ( M − and U ( κ, θ, I ) = (cid:2) ( U (exp( κ ) , ζ ( θ ) , I (1) )) (cid:62) , . . . , ( U (exp( κ ) , ζ ( θ ) , I ( M − )) (cid:62) (cid:3) (cid:62) ∈ R M ( M − , where the measured (noisy) electrode voltages { V ( i ) } M − i =1 ⊂ R M (cid:5) correspond to thecurrent patterns { I ( i ) } M − i =1 ⊂ R M (cid:5) , and { U (exp( κ ) , ζ ( θ ) , I ( i ) ) } M − i =1 are the electrodepotential components of the solutions to the forward problem (2.4) for given domainlog-conductivity κ , parametrized boundary conductivity ζ ( θ ) and current pattern I ( i ) .The ‘current basis’ for R M − (cid:5) is defined via I ( i ) = c ( e − e i +1 ) with a physicallyreasonable c ∈ R and e , . . . , e m ∈ R M denoting the standard basis vectors, thatis, the first electrode feeds the current into the body and it is guided out in turnsthrough the other M − U ( κ, θ, I ) on I to shorten the notation.We formulate the reconstruction problem as minimization of a Tikhonov-typefunctional: arg min κ,θ (cid:16) (cid:107)U ( κ, θ ) − V(cid:107) − + (cid:107) κ − κ µ (cid:107) − κ + (cid:107) θ − θ µ (cid:107) − θ (cid:17) , (3.1)where the notation (cid:107) x (cid:107) A := x (cid:62) Ax is used. We have included in (3.1) prior informa-tion motivated by the Bayesian inversion paradigm: the symmetric positive definitecovariance matrices Γ noise , Γ κ and Γ θ correspond to, respectively, the measurementnoise, the domain log-conductivity and the parametrization of the boundary conduc-tivity, while κ µ and θ µ are the expected values for the to-be-estimated variables. Tobe more precise, (3.1) is equivalent to finding a maximum a posteriori estimate for theparameters defining the domain and boundary conductivities assuming that they andthe additive measurement noise follow a priori the Gaussian distributions N ( κ µ , Γ κ ), N ( θ µ , Γ θ ) and N (0 , Γ noise ), respectively; see, e.g., [20]. As the covariance matricesare assumed positive definite, the minimization functional is bounded from below.For the piecewise linear parametrization of the domain conductivity, we alwaysuse a Gaussian smoothness prior with the covariance structure(Γ κ ) ij = γ κ exp (cid:18) − | x i − x j | λ κ (cid:19) , (3.2) However, the conductivity σ = exp κ itself is visualized in the presented reconstructions tofacilitate comparison with previous works. J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN where x i , x j ∈ R are nodal positions in the FE mesh of Ω h , γ κ is the pointwisevariance, and λ κ is the correlation length controlling the expected spatial smoothnessof the reconstructed admittivity. The employed forms for Γ θ are considered in thefollowing subsection where our parametrizations for the boundary conductivity aredescribed.We aim at solving (3.1) iteratively using the Gauss–Newton algorithm [29]. Let J ( κ, θ ) = (cid:2) ∂ κ U ( κ, θ ) , ∂ θ U ( κ, θ ) (cid:3) be the Jacobian matrix of U ( κ, θ ) with respect tothe discretized domain log-conductivity κ and the parameters defining the boundaryconductivity θ . The Jacobian matrix can be assembled using the representation for-mula for the Fr´echet derivative (2.16). See Section 3.1 below for information on theemployed parametrization θ (cid:55)→ ζ ( θ ) as well as on using (2.16) in connection to thoseparametrizations.It is straightforward to see that applying the basic form of the Gauss–Newtonalgorithm to (3.1) corresponds to iteratively linearizing U ( κ, θ ) in the first term ofthe to-be-minimized functional in (3.1) around the current estimate for the solutionand subsequently defining the next estimate to be the minimizer of the resultingquadratic Tikhonov functional. To this end, observe that the aforementioned lin-earization around a given point τ := ( κ, θ ) leads via slight abuse of notation to theleast squares problemarg min ∆ τ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:32)(cid:34) J ( τ )I (cid:35) ∆ τ − (cid:34) U ( τ ) − V τ − τ µ (cid:35)(cid:33)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , (3.3)where ∆ τ denotes the additive change in τ , τ µ = ( κ µ , θ µ ) is the expected value forthe total parameter vector and I is an identity matrix of the appropriate size. Fur-thermore, L is a Cholesky factor for the block diagonal concatenation of the inversesof the covariance matrices Γ noise , Γ κ and Γ θ , i.e. L (cid:62) L = diag(Γ − , Γ − κ , Γ − θ ).Accompanying the basic Gauss–Newton algorithm with line search to deduce anoptimal step size for each update, we finally end up with the following algorithm: Algorithm 1:
Simultaneous reconstruction of κ and θ Set iteration count to j = 0. Choose an initial guess τ = ( κ , θ ). while Stopping criteria not met do Compute J j = J ( τ j ) Determine the step direction d j as the least squares solution of L (cid:34) J j I (cid:35) d j = L (cid:34) U ( τ j ) − V τ j − τ µ (cid:35) . Using a backtracking line search, seek a step length t ∈ [0 ,
1] such that τ j + td j obtains sufficient decrease on the original Tikhonov functionalin (3.1); compare Theorem 3.2. Update τ j +1 = τ j + t j d j . Update j = j + 1. LECTRODELESS ELECTRODE MODEL ζ . In the standard CEM the parametrization of the boundary conductivity can be inter-preted as a piecewise constant function defined on the boundary of the domain ζ CEM ( x ) = M (cid:88) m =1 θ m ˜ χ m ( x ) , θ m ∈ R + , x ∈ ∂ Ω , where ˜ χ m is the characteristic function of the true electrode e m on ∂ Ω. Althoughthis model is known to be precise enough to match measurement accuracy [7, 32], ithas some often ignored inconveniences. First, if the exact locations of the electrodecontacts are not known, significant reconstruction errors can be expected [2, 4, 9, 23].Second, the piecewise constant nature of the contact conductivity (or resistivity) hasa negative effect on the smoothness of the solution to (2.4), which in turn decreasesthe highest attainable convergence rate for, say, higher order FEMs when they areapplied to numerically solving (2.4) or evaluating (2.16) [17]. Finally, there seems tobe no other (physical) justification for a constant boundary conductivity over eachelectrode except the Occam’s razor: the simplest model is the most attractive one.In reality, there is no obvious reason to expect equally good contact across the wholeelectrode, especially if the contact is established by simply pressing a metal electrodetip against skin in some medical application of EIT.Our first novel parametrization for the boundary conductivity corresponds tosimply using the piecewise linear FE basis restricted to the boundary. This optionovercomes all issues listed above. However, as it carries numerous degrees of freedom,it also requires a strong prior or regularization. We call this option the piecewiselinear (PL) parametrization ζ PL : θ (cid:55)→ ζ PL ( θ ) , R N → H ( ∪ E m ) , (3.4)where N is the number of finite element nodes on ∪ E m . To assure nonnegativityof the boundary conductivity without having to introduce constraints in the opti-mization problem (3.1), the parameter vector θ ∈ R N defines the nodal values as ζ = θ , . . . , ζ N = θ N . The boundary conductivity on the rest of the extended elec-trodes is then simply the linear interpolant of the neighboring nodal values. We tacitlyassign the boundary nodes on ∂ Ω \ ∪ E m with zero boundary conductivities, thus ex-tending ζ PL to the whole boundary ∂ Ω as an element of H ( ∂ Ω). Consult [17] forinformation on the effect of such a parametrization on the convergence of FEMs forthe CEM.The PL parametrization is regularized assuming a similar Gaussian smoothnessprior as for the domain conductivity. To this end, let J m denote the index set for thenodal indices on the closure of the m th extended electrode. A preliminary covariancematrix for the parameters defining the nodal values of the boundary conductivity astheir squares is then formed via (cid:0) ˜Γ θ (cid:1) ij = γ θ exp (cid:16) − d ∂ Ω ( x i ,x j ) λ θ (cid:17) if i, j ∈ J m for some m ∈ { , . . . , M } , , (3.5)where d ∂ Ω ( · , · ) measures the Euclidean distance along the boundary curve, γ θ is thecomponentwise variance and λ θ is the correlation length along the boundary. Theexpected value θ µ is simply set to zero, as in our numerical experiments its value didnot significantly affect the reconstructions. The zero-mean Gaussian density defined2 J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN by these parameters is finally conditioned by fixing the contact conductivities at theend points of the electrodes to zero; the corresponding covariance matrix Γ θ for theparameters defining the nodal values on the interior of the extended electrodes canbe obtained by forming a suitable Schur complement based on ˜Γ θ (see, e.g., [20]).Our second parametrization is a modified version of the “hat-CEM” introducedin [17]. However, instead of letting the support of a hat-shaped boundary conductivityto cover a whole (extended) electrode, we allow its location and width to vary on eachextended electrode. More precisely, we define the parametric hat (PH) boundaryconductivity electrode-wise through( ζ PH ) m : (cid:40) ( h m , l m , w m ) (cid:55)→ ( ζ PH ) m ( h m , l m , w m ) , (cid:8) ( h, l, w ) ∈ R | l − w ≥ l + w ≤ (cid:9) =: Z (cid:48) PH → H ( E m ) , (3.6)with( ζ PH ) m ( h m , l m , w m )( t ) = h m ( − l m + w m +2 t ) w m if t ∈ ( l m − w m , l m ) , h m (2 l m + w m − t ) w m if t ∈ [ l m , l m + w m ) , t takes values in the interval (0 ,
1) as a normalized parametrization for thecurve segment E m . Observe that the parameters l m , w m and h m have clear physicalinterpretations; they correspond, respectively, to the (normalized) center point ofthe electrode contact, the (normalized) width of the contact and the net contactconductance. However, as with the standard CEM, there is no reason to actuallyexpect the boundary conductivity distribution to follow the shape of a hat function.The main motivation for the parametrization (3.7) is that it is arguably the simplest‘movable’ H ( ∂ Ω) parametrization for the boundary conductivity.Combining all extended electrodes and again setting the values of the boundaryconductivity on ∂ Ω \ ∪ E m to be zero, the parametrization (3.7) naturally extends tothe whole domain boundary as ζ PH : θ = ( h , l , w ) (cid:55)→ ζ PH ( h , l , w ) = ζ PH ( θ ) , ( Z (cid:48) PH ) M → H ( ∂ Ω) . As any reasonable finite element mesh has more than three edges per an (extended)electrode, the PH parametrization leads to a significantly lower number of degreesof freedom than the PL parametrization. For the sake of completeness, we haveincluded the representation formulas for computing the derivatives of the electrodemeasurements with respect to h , l and w , needed for constructing the Jacobian matrixwith respect to the parametrization of the boundary conductivity, in Appendix A;cf. Remark 2.7.We found the standard Tikhonov regularization with a simple diagonal weight ma-trix adequate for regularizing the PH boundary conductivity parametrization. Thatis, Γ θ = Γ h Γ l Γ w = γ h I γ l I γ w I , ∈ R M × M , (3.8)where γ h , γ l , γ w > Remark 3.1.
For the PH parametrization, we modify the line search step inAlgorithm 1 to ensure the contact conductivity parameters stay within Z (cid:48) PH and thus LECTRODELESS ELECTRODE MODEL define hat-shaped functions supported on the closures of the extended electrodes. Fora total parameter vector τ j + td j under consideration, denote by ( h m,j , l m,j , w m,j ) thecorresponding parameters defining the contact conductivity on the m th extended elec-trode according to (3.7) . Before evaluating the Tikhonov functional (3.1) at τ j + td j ,the parameter triplet ( h m,j , l m,j , w m,j ) is clamped inside Z (cid:48) PH for each m = 1 , . . . , M via the following steps (cf. (3.6) ): (i) If w m,j > , then redefine w m,j = 1 . (ii) If l m,j − w m,j < , then redefine l m,j = w m,j . (iii) If l m,j + w m,j > , then redefine l m,j = 1 − w m,j . In particular, this modification means the partial convergence proofpresented in the following section does not cover the case of the PH parametrization. For the PL parametrization of theboundary conductivity, a limited form of convergence of the Gauss–Newton methodof Algorithm 1 can be deduced from basic results [29]. The PH parametrization ismore involved due to the clamping of the conductivity parameters to the admissibledomain; see Remark 3.1. If a computationally expensive weighted projection wasperformed instead of simple clamping, and instead of performing a line search, localgrowth conditions were satisfied, the results of [30] could be used. Alternatively, ifwe were to add a proximal relaxation term to the method, the results of [19] could beused, which allow for such constraints and even nonsmooth regularization.
Theorem 3.2.
Assume that the PL parametrization (3.4) is employed for theboundary conductivity and suppose Γ κ , Γ θ ≤ γ I for some γ > . Write F for theobjective of (3.1) . Let { τ j } ∞ j =0 be generated by Algorithm 1 with the step lengths t j satisfying for given parameters < c < c < the Wolfe conditions F ( τ j + t j d j ) ≤ F ( τ j ) + c t j F (cid:48) ( τ j ) d j (sufficient decrease) and F (cid:48) ( τ j + t j d j ) d j ≥ c F (cid:48) ( τ j ) d j (curvature) . Then F (cid:48) ( τ j ) → , the set of accumulation points ˆ τ of { τ j } j ∈ N is non-empty, and all ofthem are critical, i.e., F (cid:48) (ˆ τ ) = 0 . Moreover, there always exist subintervals of [0 , ∞ ) on which t j satisfies the Wolfe conditions.Proof . It follows from Theorem 2.5 and the solution formula (2.16) that U iscontinuously differentiable as a function of ( σ, ζ ) with I fixed. The same is then truefor the finite dimensional function U ( κ, θ ) appearing in (3.1) due to the smoothnessof the functions s (cid:55)→ exp( s ) and s (cid:55)→ s ; see (3.4). It follows that F is continuouslydifferentiable. That the Wolfe conditions can be satisfied now follows from [29, Lemma3.1] and the continuous differentiability of F .For any real parameter r , we denote lev r F := { τ | F ( τ ) ≤ r } the r -sublevelset of F . Since F is coercive due to the regularization terms in (3.1) and the boundΓ κ , Γ θ ≤ γ I, the sublevel set lev F ( τ ) is bounded. Due to the finite-dimensionalityof the domain, and the continuous differentiability of U , it follows using a compactcovering argument that J = ( U σ , U θ ) is Lipschitz within N for any bounded open set N ⊃ lev F ( τ ) F . Due to the sufficient decrease in the Wolfe conditions as well as d j being a direction of descend for F , i.e., F (cid:48) ( τ j ) d j ≤
0, see the proof of [29, Theorem10.1], we also have the monotonicity F ( τ j +1 ) ≤ F ( τ j ). Consequently { τ j } j ∈ N ⊂ lev F ( τ ) .To prove convergence, we set R ( τ ) := L (cid:18) U ( τ ) − V τ − τ µ (cid:19) . Then F ( τ ) = (cid:107) R ( τ ) (cid:107) whereas the least squares problem (3.3) can be written more4 J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN compactly as min ∆ τ (cid:107) R ( τ ) + R (cid:48) ( τ )∆ τ (cid:107) , where R (cid:48) ( τ ) = L (cid:18) J ( τ )I (cid:19) . That F (cid:48) ( τ j ) → R (cid:48) is Lipschitz in an open bounded neighborhood N of lev F ( τ ) F ,b) R (cid:48) ( τ ) (cid:62) R (cid:48) ( τ ) ≥ ˜ γ I for some ˜ γ > L (cid:62) L = diag(Γ − , Γ − κ , Γ − θ ), we have L (cid:62) L ≥ (cid:104) Γ − γ − I (cid:105) . It follows that R (cid:48) ( τ ) (cid:62) R (cid:48) ( τ ) = J ( τ ) (cid:62) Γ − J ( τ ) + γ − I, so that b) holds with ˜ γ = γ − . On theother hand, condition a) holds if J is Lipschitz within N , which we have also shown.Thus F (cid:48) ( τ j ) → { j k } k ∈ N ⊂ N such that τ j k converges to some ˆ τ follows from finite-dimensionality and the boundedness of lev F ( τ ) F . Since, for anysuch subsequence, F (cid:48) ( τ j k ) → F (cid:48) is continuous by the continuity and continuousdifferentiability of U , necessarily F (cid:48) (ˆ τ ) = 0.
4. Numerical examples.
In this section we present some experimental resultsthat support the use of our new model in connection with real world EIT measure-ments. To put it short, including the estimation of a boundary contact function asa part of a reconstruction algorithm for EIT can result in qualitatively better re-construction of the primary unknown, i.e. the domain conductivity, if it is expectedthat the electrode positions may have been mismodeled. Similar conclusions havepreviously been reached in [5, 9, 10, 11, 16] by resorting to shape derivatives of EITmeasurements with respect to the electrode locations. Having (slightly) incompleteinformation on the electrode positions is an important case in practice as the contactlocations are exactly known only in trivial experimental settings such as water tanks.Be that as it may, we resort to water tank experiments in what follows.We consider the same data as in [17] . The measurements were performed usingthe Kuopio impedance tomography (KIT4) device [24] on a thorax-shaped tank havinga circumference of 106 cm; see Figure 4.1(b). The tank was equipped with 16 near-equidistant 2 cm wide rectangular electrodes that reached from the bottom of the tankto the water surface at the height of 5 cm. The measurements were conducted usinglow-frequency common ground current patterns with amplitudes of approximately1.2 mA. The phase information in the potential measurements is ignored and theiramplitudes are treated as if they had resulted from measurements with direct cur-rent. Before forming any reconstructions, the current amplitudes (and the respectiveelectrode potentials) are scaled to 1.0 mA to avoid uneven weighting between currentpatterns. As the measurement setup is homogeneous in the vertical direction, it canbe modeled by a two-dimensional version of the forward problem (2.4). For furtherdetails on the measurement setup, see [17]. The source code of our implementationhas been made available on Zenodo at doi:10.5281/zenodo.4475524In all experiments, the two-dimensional computational domain is triangulated asa FE mesh consisting of 2725 nodes, being densest near the (extended) electrodesand sparsest near the center of the domain. The triangulations are generated using aslightly modified version of Triangle [31] . Electrical impedance dataset of thorax shaped tank doi:10.5281/zenodo.4475417 Triangle with disable attribute interpolation switch by Kuutela 2020doi:10.5281/zenodo.4472025LECTRODELESS ELECTRODE MODEL γ κ
10 V − standard deviation for the log-conductivity λ κ · − m correlation length for the log-conductivity γ θ · S/V standard deviation in the PL parametrization λ θ · − m correlation length for the PL parametrization γ h S/(m V) standard deviation for h in the PH parametrization γ l . V − standard deviation for l in the PH parametrization γ w V − standard deviation for w in the PH parametrizationTable 4.1: The parameters employed in (3.2), (3.5) and (3.8) throughout the numericalexperiments.In order to provide a reasonably fair comparison between different approaches, weuse the same values for the constants appearing in the prior covariances (3.2), (3.5)and (3.8) in all tests. The choices are based on a grid search and manual verificationsacross several test cases. The employed parameter values are listed in Table 4.1. Takenote that we have made the naive choice of Γ noise = I ∈ R in our numerical tests,which means that the standard deviations in Table 4.1 are in fact the ratios betweenthe standard deviations for the listed parameters and that for the additive zero-meannoise that is assumed to have the same variance and be mutually independent overall potential measurements. Moreover, the weighted norm (cid:107) · (cid:107) Γ − appearing in(3.1) simply becomes the Euclidean norm (cid:107) · (cid:107) . When computing reconstructionsbased on the traditional CEM with constant contacts (over the actual electrodes), noregularization with respect to the contact conductances is needed or used.We consider three general geometries for the interplay between the true and ex-tended electrodes. First, we present results for (almost) exactly known electrodepositions, that is, the extended electrodes E m coincide with e m for m = 1 , . . . , M up to the available information on the precise positions of the latter on the interiorsurface of the water tank. In particular, the available information on the electrodepositions is used to its full extent when forming the reconstruction with all threeparametrization, i.e. the standard CEM as well as the PL and PH parametrizationsfor the boundary conductivity in the new electrodeless model. In the second andthird cases, the extended electrodes are approximately 12 mm and 22 mm wider, re-spectively, than the physical electrodes, with the corresponding extensions divided(uniformly) randomly between the two endpoints of the electrodes. When comput-ing reconstructions based on the standard CEM in these cases, we assume that theelectrode widths are known but intentionally mismodel the setup by placing eachcomputational electrode in the middle of its extended counterpart. This amounts onaverage to approximately 3 mm and 6 mm and at most 6 mm and 11 mm of displace-ment, respectively, for the computational electrodes compared to the true physicalsetup when the standard CEM is used for forming the reconstructions. In our first experiment, we con-sider a “trivial” case where the tank is filled with mere tap water without any embed-ded inhomogeneities. The domain conductivity is parametrized by a single variablewith no regularization, i.e. no penalty term for the conductivity is included in theminimized Tikhonov functional. We are particularly interested in the final residual (cid:107)U ( κ opt , θ opt ) − V(cid:107) (4.1)6 J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN
Residuals Domain conductivitiesExact 12 mm 22 mm Exact 12 mm 22 mmCEM 0.0795 0.150 0.301 0.0228 0.0229 0.0231PL 0.0450 0.0349 0.0971 0.0228 0.0227 0.0228PH 0.249 0.0382 0.0369 0.0227 0.0227 0.0227Table 4.2: Residuals (4.1) and the estimated conductivity levels (S/m) for a tankwithout embedded inhomogeneities.and the resulting optimal constant conductivity level σ opt = exp κ opt for the con-sidered models, i.e. the standard CEM and the electrodeless model with the PLand PH parametrizations for the boundary conductivity. Our secondary interest lieswith the logarithms of the estimated net electrode conductances log (cid:82) E m ζ ( θ opt ) d s , m = 1 , . . . , M , where θ opt denotes the optimal parameters for the boundary conduc-tivity produced by Algorithm 1 for the considered model and ζ = ζ CEM , ζ = ζ PL or ζ = ζ PH .As mentioned in Section 3.1, the expected value for the contact conductancein the PL parametrization is the zero function, whereas those for the parametersappearing in the PH parametrization are zero for the net conductance h , the center ofthe extended electrode for the center of the hat function l , and the relative length ofthe true electrode compared to the employed extended electrode for the hat width w .The initial guesses for the contact conductivities in the PL and PH parametrizations,as well as for the standard CEM, are chosen so that the net conductance over eachelectrode is 0.001 S. For the PL model, the initial nodal values are the same over thewhole extended electrode, whereas for the PH model, the initial hat is at the centerof the corresponding extended electrode, and it has the same width as the physicalelectrodes. In each case, Algorithm 1 is run until obvious convergence, or up to 50iterations.The residuals and the estimated domain conductivity levels for the three modelsand the three geometric settings for the interplay between the electrodes and the ex-tended electrodes are listed in Table 4.2, whereas the logarithms of the estimated netelectrode conductances and their variances over the electrodes are given in Table 4.3.For comparison, the conductivity levels estimated for the same measurement setupin [17] were 0.022722 S/m for the standard CEM and 0.022720 S/m for a hat shapedcontact conductivity model with known electrode positions and no regularization withrespect to the height of the hat in the inversion algorithm. The corresponding recon-structed contact conductance levels at the electrodes are shown in [17, Fig. 7].When the electrode positions are modeled (almost) exactly, the PH parametriza-tion results in a larger residual (4.1) than the standard CEM and the PL model. Theexact reason for this discrepancy is unknown, but it is presumably related to the priormodel associated to the PH parametrization. All estimated domain conductivities arewell aligned with the values from [17].When the extended electrodes are 12 mm wider than the true ones, i.e. the com-putational electrodes for the standard CEM are misplaced on average by 3.0 mm, theresiduals corresponding to the PL and, especially, the PH parametrization decreasewhereas the residual for the standard CEM increases. All estimates for the domainconductivity are still within one per cent of the values reported in [17]. It seems thatthe slight widening of the extended electrodes gives the PL and PH parametrizationsmore freedom to match the measurement data, whereas the mismodeling of the true LECTRODELESS ELECTRODE MODEL (cid:82) E m ζ d s Standard deviationsExact 12 mm 22 mm Exact 12 mm 22 mmCEM -3.65 -3.82 -4.19 0.481 0.363 0.246PL -2.61 -4.47 -4.83 0.518 0.402 0.312PH -3.39 -3.34 -3.40 0.179 0.272 0.335Table 4.3: Means and standard deviations for the logarithms of the reconstructed netelectrode conductances (S) for a tank without embedded inhomogeneities.electrodes slightly impairs the performance of the standard CEM in this regard.Increasing the difference in widths between the extended and true electrodes to22 mm, and thus the average misplacement of the computational CEM electrodes to5.6 mm, essentially highlights the observations regarding the residuals (4.1) in theprevious case: employing the standard CEM leads to a residual that is significantlylarger than those for the PL and PH parametrizations. Both PH and PL models leadto domain conductivity estimates that are still in the vicinity of those listed in [17],while the estimate produced by the CEM is somewhat inaccurate with the algorithmalso having some trouble converging.The means and standard deviations for the logarithms of the reconstructed netelectrode conductances are listed in Table 4.3. The net conductances are on averagelargest for the PH parametrization, which is once again in line with the observationsin [17]. In particular, the net integrals of the contact conductivity over the electrodesis obviously not the only property of the contacts that affects the (simulated) electrodemeasurements, but the shape of the contact conductivity over individual electrodesalso plays a role. Mismodeling of the electrode positions seems to decrease the esti-mated net conductances for the standard CEM and the PL parametrization. One ofthe electrodes, i.e. number 15, seems to have a significantly better contact than theothers, which does not have any clear explanation, but the observation is anywaysaligned with the material in [17].
We next consider recon-structing a nontrivial conductivity distribution inside the same water tank. Twoinclusions are placed in the tank: an insulating plastic cylinder and a rectangularmetal pipe as shown in Figure 4.1(b). They are both homogeneous in the verticaldirection and break the water surface, making the measurement configuration onceagain modelable by a two-dimensional version of (2.4). This time the conductivityfield is discretized using the same 2725 piecewise linear FE basis functions as the onesused for solving (2.4) in the two-dimensional domain. Encouraged by the experimentswith the empty water tank, we choose σ = σ µ = 0 .
02 S/m as both the homogeneousinitial guess and the expected value for the domain conductivity. The initial guessesand expected values for the PL and PH contact conductivity parametrizations are asin the previous experiment. In each reconstruction process, we let Algorithm 1 rununtil clear convergence or up to 50 iterations. The models for the considered extendedelectrodes and the resulting errors in the placement of the computational electrodesfor the standard CEM are the same as for the experiments with the empty tank.The reconstructions for the three forward models, as well as for the three con-figurations for the (computational) electrodes and extended electrodes, are shownin Figure 4.1. An immediate observation is that the quality of the reconstructionsproduced by the PL and PH parametrizations for the domain conductivity does notsuffer noticeably from the mismodeling of the electrode positions; the freedom in these8
J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN(a) Reconstructions(b) Target
Fig. 4.1: Domain conductivity reconstructions (S/m) for a water tank with an in-sulating (top left) and a well conducting (bottom right) inclusions. The columnscorrespond to the standard CEM and the new electrodeless model with the PL andPH parametrizations. The rows are associated with different levels of error in themodel for the electrode positions. The darkest extended electrode corresponds to m = 1 and the others are numbered in the positive directions.models for fine-tuning the contact conductivity shapes along the extended electrodesseems to provide enough flexibility for coping with the geometric modeling error. Onthe other hand, the reconstructions corresponding to the standard CEM significantlydeteriorate as the mean amount of misplacement in the (computational) electrodesincreases from 0 to 5.6 mm. To be more precise, conductivity artifacts start to appearespecially near the domain boundary.Table 4.4 shows the final values for the different quadratic terms appearing in theminimized Tikhonov functional (3.1) for the three forward models and the extendedelectrodes that are 22 mm wider than the true ones, corresponding on average toa misplacement of 5.6 mm for the computational electrodes in the standard CEM.(Recall that there is no penalization with respect to the constant contact conductancesin (3.1) for the standard CEM, which explains the missing value in Table 4.4.) It isobvious that with the considered prior model for the domain log-conductivity, the LECTRODELESS ELECTRODE MODEL (cid:107)U ( κ, θ ) − V(cid:107) (cid:107) κ − κ µ (cid:107) Γ − κ (cid:107) θ − θ µ (cid:107) Γ − θ CEM 0.255 0.0318 —PL 0.0182 0.0187 0.00864PH 0.0186 0.0195 0.0128Table 4.4: Final values of the three terms in (3.1) for the standard CEM and thePL and PH parametrizations for the boundary conductivity. There is no penalizationwith respect to the constant contact conductivities in the standard CEM.standard CEM with the erroneous electrode positions is not able to bring the datafit term (cid:107)U ( κ, θ ) − V(cid:107) down to the same level as the more flexible PL and PHparametrizations. A lower value for this residual could be achieved by relaxing theprior for κ = log σ , but this would inevitably lead to even more severe artifacts in theassociated reconstruction of σ . A natural question is whetherour proposed methods actually employ the flexibility in the boundary conductiv-ity parametrizations to reconstruct the electrode positions or do the reconstructionsmerely fit to the measurements with no clear physical interpretation. To this end,Figure 4.2 visualizes the reconstructed boundary conductivities for the experimentswith the empty water tank in Section 4.1 (left column) and with the tank enclosingthe two cylindrical inclusions in Section 4.2 (right column). All images in Figure 4.2correspond to extended electrodes that are 22 mm wider than the physical ones. Ac-cording to a qualitative assessment of the rows in Figure 4.2, the reconstructions of thecontact conductivity associated to the PH parametrization are similar for the emptytank and the variable conductivity case on each electrode, although the equivalencecannot be described as perfect. A similar conclusion does not hold on all electrodes forthe PL parametrization; for the empty tank, the contact is divided on some extendedelectrodes into two humps located near the end points, which must be considered abad approximation of the underlying physical reality.In particular, according to Figure 4.2, Algorithm 1 seems to usually move thecenter of mass for the PH contact conductivity parametrization toward the centerof the actual physical electrode. To quantify such a claim, let s m ∈ R denote themidpoint of the m th physical electrode in an arclength parametrization of ∂ Ω, and letˆ s m ∈ R be the center of mass for some reconstructed boundary conductivity on E m with respect to that same arclength parametrization. We define the correspondingmean error in the localization of the electrode contacts aserr E = 1 M M (cid:88) m =1 | s m − ˆ s m | . (4.2)Table 4.5 lists these mean errors in the reconstructed electrode positions for the exper-iments with 22 mm electrode extensions in Sections 4.1 and 4.2, with err E = 5 .
57 mmfor the standard CEM corresponding to the random variation in the positions of the(extended) electrodes in the computational model. The numbers in Table 4.5 confirmthat the contact centers indeed move closer to the midpoints of the physical targetelectrodes for the PH parametrization, but for the PL parametrization the results donot indicate a capability to locate the true electrodes. In particular, it seems that theunphysical division of some contacts into two parts leads to a worse performance thanthe trivial guess (cf. the CEM) for the PL parametrization in the case of the emptytank.0
J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN -1 E l ec t r od e Homogeneous domain -1 E l ec t r od e -1 E l ec t r od e -1 E l ec t r od e -1 E l ec t r od e -1 E l ec t r od e -1 E l ec t r od e -1 E l ec t r od e -1 Domain with inclusions -1 -1 -1 -1 -1 -1 -1 CEM PL PH
Fig. 4.2: Contact conductivity reconstructions on the odd electrodes for the experi-ments in Section 4.1 (left) and Section 4.2 (right) with 22 mm electrode extensions.The horizontal axis represents arclength parametrization in the counter-clockwise di-rection and in millimeters for the extended electrode, with the ticks indicating nodesin the FE mesh. The vertical axis is logarithmic in S/m . The physical electrodeposition is indicated by a gray background and the red line depicts the reconstructedpiecewise constant contact conductivity for the standard CEM at the (randomly)erroneous position. The dark and light blue curves are the reconstructed contact con-ductivities for the PL and PH parametrizations, respectively. Each inverted trianglemarks the center of the contact conductivity mass for a particular reconstruction,while the black upright triangle is the center of the true electrode.The general features of the numerical results documented in this (and the previ-ous) sections were consistently reproduced with different randomizations of the elec-trode extensions and for almost all experimental data for the considered water tankat our disposal. That is, when Algorithm 1 is applied to the PL and PH parametriza-tions, the reconstructed boundary conductivities tend to compensate for the mismod-eled electrode positions regardless of the domain conductance. However, it is unclearwhether this extends to more complicated geometric settings and/or for other regular-ization or statistical inversion strategies for tackling the inverse problem of practical LECTRODELESS ELECTRODE MODEL
5. Concluding remarks.
This work introduced a new approach to modelingelectrodes via introducing a spatially varying boundary admittivity function in theframework of the CEM of EIT under only generic information on the electrode lo-cations. Two parametrizations for the boundary admittivity were considered: onedirectly employing the underlying piecewise linear FE basis and another correspond-ing to a single hat-shaped admittivity with varying width, height and location oneach extended electrode. Our numerical experiments based on experimental watertank data demonstrated that reconstructing the boundary admittivity function as apart of an iterative Bayesian output least squares algorithm leads to far better recon-structions of the domain admittivity than simply ignoring the incomplete informationon the electrode positions.
Appendix A. Derivatives for the PH parametrization.
The derivatives ofelectrode measurements with respect to the free parameters in the PH electrode con-ductance parametrization (3.7) can be straightforwardly deduced using Remark 2.7: ∂ h m U ( σ, ζ PH ( h , l , w ) , I ) | E m | · ˜ I = − (cid:90) l m l m − w m / − l m + w m + 2 t ) w m ( U − u )( ˜ U − ˜ u ) d t − (cid:90) l m + w m / l m l m + w m − t ) w m ( U − u )( ˜ U − ˜ u ) d t,∂ l m U ( σ, ζ PH ( h , l , w ) , I ) | E m | · ˜ I = (cid:90) l m l m − w m / h m w m ( U − u )( ˜ U − ˜ u ) d t − (cid:90) l m + w m / l m h m w m ( U − u )( ˜ U − ˜ u ) d t,∂ w m U ( σ, ζ PH ( h , l , w ) , I ) | E m | · ˜ I = (cid:90) l m l m − w m / h m ( − l m + w + 4 t ) w m ( U − u )( ˜ U − ˜ u ) d t + (cid:90) l m + w m / l m h m (4 l m + w − t ) w m ( U − u )( ˜ U − ˜ u ) d t, where the integrals correspond to a normalized parametrization of the m th extendedelectrode curve segment. REFERENCES[1]
D. C. Barber and B. H. Brown , Applied potential tomography , J. Phys. E: Sci. Instrum., 17(1984), pp. 723–733. J. DARD´E, N. HYV ¨ONEN, T. KUUTELA, AND T. VALKONEN[2]
D. C. Barber and B. H. Brown , Errors in reconstruction of resistivity images using a linearreconstruction technique , Clin. Phys. Physiol. Meas., 9 (1988), pp. 101–104.[3]
L. Borcea , Electrical impedance tomography , Inverse problems, 18 (2002), pp. R99–R136.[4]
W. Breckon and M. Pidcock , Data errors and reconstruction algorithms in electricalimpedance tomography , Clin. Phys. Physiol. Meas., 9 (1988), pp. 105–109.[5]
V. Candiani, A. Hannukainen, and N. Hyv¨onen , Computational framework for applying elec-trical impedance tomography to head imaging , SIAM J. Sci. Comput., 41 (2019), pp. B1034–B1060.[6]
M. Cheney, D. Isaacson, and J. Newell , Electrical impedance tomography , SIAM Rev., 41(1999), pp. 85–101.[7]
K.-S. Cheng, D. Isaacson, J. S. Newell, and D. G. Gisser , Electrode models for electriccurrent computed tomography , IEEE Trans. Biomed. Eng., 36 (1989), pp. 918–924.[8]
P. G. Ciarlet , Linear and Nonlinear Functional Analysis with Applications , Society for In-dustrial and Applied Mathematics, 2013.[9]
J. Dard´e, H. Hakula, N. Hyv¨onen, and S. Staboulis , Fine-tuning electrode information inelectrical impedance tomography , Inverse Probl. Imag., 6 (2012), pp. 399–421.[10]
J. Dard´e, N. Hyv¨onen, A. Sepp¨anen, and S. Staboulis , Simultaneous reconstruction of outerboundary shape and admittivity distribution in electrical impedance tomography , SIAM J.Imaging Sci., 6 (2013), pp. 176–198.[11]
J. Dard´e, N. Hyv¨onen, A. Sepp¨anen, and S. Staboulis , Simultaneous recovery of admittivityand body shape in electrical impedance tomography: An experimental evaluation , InverseProblems, 29 (2013), p. 085004.[12]
J. Dard´e and S. Staboulis , Electrode modelling: The effect of contact impedance , ESAIM:Math. Model. Num., 50 (2016), pp. 415–431.[13]
M. Hanke, B. Harrach, and N. Hyv¨onen , Justification of point electrode models in electricalimpedance tomography , Math. Models Methods Appl. Sci., 21 (2011), pp. 1395–1413.[14]
N. Hyv¨onen , Approximating idealized boundary data of electric impedance tomography byelectrode measurements , Math. Models Methods Appl. Sci., 19 (2009), pp. 1185–1202.[15]
N. Hyv¨onen, V. Kaarnioja, L. Mustonen, and S. Staboulis , Polynomial collocation forhandling an inaccurately known measurement configuration in electrical impedance to-mography , SIAM J. Appl. Math., 77 (2017), pp. 202–223.[16]
N. Hyv¨onen, H. Majander, and S. Staboulis , Compensation for geometric modeling er-rors by positioning of electrodes in electrical impedance tomography , Inverse Problems, 33(2017), p. 035006.[17]
N. Hyv¨onen and L. Mustonen , Smoothened electrode model , SIAM J. Appl. Math., 77 (2017),pp. 2250–2271.[18] ,
Generalized linearization techniques in electrical impedance tomography , Numer. Math.,140 (2018), pp. 95–120.[19]
J. Jauhiainen, P. Kuusela, A. Sepp¨anen, and T. Valkonen , Relaxed Gauss–Newton methodswith applications to electrical impedance tomography , SIAM Journal on Imaging Sciences,13 (2020), pp. 1415–1445.[20]
J. Kaipio and E. Somersalo , Statistical and Computational Inverse Problems , Springer, 2005.[21]
V. Kolehmainen, M. Lassas, and P. Ola , Inverse conductivity problem with an imperfectlyknown boundary , SIAM J. Appl. Math., 66 (2005), pp. 365–383.[22] ,
The inverse conductivity problem with an imperfectly known boundary in three dimen-sions , SIAM J. Appl. Math., 67 (2007), pp. 1440–1452.[23]
V. Kolehmainen, M. Vauhkonen, P. A. Karjalainen, and J. P. Kaipio , Assessment oferrors in static electrical impedance tomography with adjacent and trigonometric currentpatterns , Physiol. Meas., 18 (1997), pp. 289–303.[24]
J. Kourunen, T. Savolainen, A. Lehikoinen, M. Vauhkonen, and L. M. Heikkinen , Suit-ability of a PXI platform for an electrical impedance tomography system , Meas. Sci. Tech-nol., 20 (2009), p. 015503.[25]
A. Lechleiter and A. Rieder , Newton regularizations for impedance tomography: A numer-ical study , Inverse Problems, 22 (2006), pp. 1967–1987.[26]
J. Neˇcas , Direct methods in the theory of elliptic equations , Springer Monographs in Mathe-matics, Springer, Heidelberg, 2012. Translated from the 1967 French original by GerardTronel and Alois Kufner, Editorial coordination and preface by ˇS´arka Neˇcasov´a and acontribution by Christian G. Simader.[27]
A. Nissinen, V. Kolehmainen, and J. P. Kaipio , Compensation of modelling errors due tounknown domain boundary in electrical impedance tomography , IEEE Trans. Med. Imag.,30 (2011), pp. 231–242.[28]
A. Nissinen, V. Kolehmainen, and J. P. Kaipio , Reconstruction of domain boundary and
LECTRODELESS ELECTRODE MODEL conductivity in electrical impedance tomography using the approximation error approach ,Int. J. Uncertain. Quantif., 1 (2011), pp. 203–222.[29] J. Nocedal and S. J. Wright , Numerical Optimization , Springer, 1999.[30]
S. Salzo and S. Villa , Convergence analysis of a proximal Gauss–Newton method , Compu-tational Optimization and Applications, 53 (2012), pp. 557–589.[31]
J. R. Shewchuk , Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Trian-gulator , in Applied Computational Geometry: Towards Geometric Engineering, M. C. Linand D. Manocha, eds., vol. 1148 of Lecture Notes in Computer Science, Springer-Verlag,May 1996, pp. 203–222. From the First ACM Workshop on Applied Computational Ge-ometry.[32]
E. Somersalo, M. Cheney, and D. Isaacson , Existence and uniqueness for electrode modelsfor electric current computed tomography , SIAM J. Appl. Math., 52 (1992), pp. 1023–1040.[33]
G. Uhlmann , Electrical impedance tomography and Calder´on’s problem , Inverse Problems, 25(2009), p. 123011.[34]