The Bayesian Formulation of EIT

The Bayesian Formulation of EIT: Analysis and Algorithms

Abstract.

We provide a rigorous Bayesian formulation of the EIT problem in an infinite dimensional setting, leading to well-posedness in the Hellinger metric with respect to the data. We focus particularly on the reconstruction of binary fields where the interface between different media is the primary unknown. We consider three different prior models – log-Gaussian, star-shaped and level set. Numerical simulations based on the implementation of MCMC are performed, illustrating the advantages and disadvantages of each type of prior in the reconstruction, in the case where the true conductivity is a binary field, and exhibiting the properties of the resulting posterior distribution.

Key words and phrases:
inverse problems, ill-posed problems, electrical impedance tomography, Bayesian regularisation, geometric priors, level set methodology.
1991 Mathematics Subject Classification:
Primary: 62G05, 65N21; Secondary: 92C55.

Matthew M. Dunlop

Computing & Mathematical Sciences,

California Institute of Technology,

Andrew M. Stuart

Computing & Mathematical Sciences,

California Institute of Technology,

1. Introduction

1.1. Background

Electrical Impedance Tomography (EIT) is an imaging technique in which the conductivity of a body is inferred from electrode measurements on its surface. Examples include medical imaging, where the methodology is used to non-invasively detect abnormal tissue within a patient, and subsurface imaging where material properties of the subsurface are determined from surface (or occasional interior) measurements of the electrical response; the methodology is often referred to as electrical resistance tomography – ERT – in this context and discussed in more detail below. The concept of EIT appears as early as the late 1970’s [15] and ERT the 1930’s [28].

A very influential mathematical formulation of the inverse problem associated with EIT dates back to 1980, due to Calderón. He formulated an abstract version of the problem, in which the objective is recovery of the coefficient of a divergence form elliptic PDE from knowledge of its Neumann-to-Dirichlet or Dirichlet-to-Neumann operator. Specifically, in the Dirichlet-to-Neumann case, if and is given, let solve

 {−∇⋅(σ∇u)=0 in Du=g on ∂D.

The problem of interest is then to ask does the mapping given by

 g↦σ∂u∂ν

determine the coefficient in ? Physically, corresponds to boundary voltage measurements, and corresponds to the current density on the boundary. Much study has been carried out on this problem – some significant results, in the case where all conductivities are in and , concern uniqueness [42], reconstruction [35], stability [2] and partial data [22]. Details of these results are summarised in [38].

In 1996, Nachman proved global uniqueness and provided a reconstruction procedure for the case , involving the use of a scattering transform and solving a D-bar problem [36]. The D-bar equation involved is a differential equation of the form , where denotes the conjugate of the complex derivative and depends on the scattering transform. A regularised D-bar approach, involving the truncation of the scattering transform, was provided in [23, 24], enabling the recovery of features of discontinuous permeabilities. The regularised D-bar approach is also used in [25], for the case when the data is not of infinite precision. Other work in the area includes joint inference of the shape of the domain and conductivity [26].

For applications, a more physically appropriate model for EIT was provided in 1992 in [40]. This model, referred to as the complete electrode model (CEM), replaces complete boundary potential measurements with measurements of constant potential along electrodes on the boundary, subject to contact impedances. The authors show that predictions from this model agree with experimental measurements up to the measurement precision of 0.1%. For this model they also prove existence and uniqueness of the associated electric potential. It is this model that we shall consider in this paper, and it is outlined in section 2.

When using the CEM, there is a limitation on the number of measurements that can be taken to provide additional information. This is due to the fact that there are only a finite number of electrodes through which current can be injected and voltages read, combined with the linear relationship between current and voltage. The data is therefore finite dimensional in the inverse problem, as distinct from the Calderón problem where knowledge of an infinite dimensional operator is assumed. As a consequence, reconstruction using the CEM often makes use of Tikhonov regularisation; such a regularisation can also be used to account for noise on the data. The paper [13] analyses numerical convergence when an or TV penalty term is used, with a finite element discretisation of the problem. We will adopt a Bayesian approach to regularisation, and this is discussed below.

A closely related problem to EIT is Electrical Resistivity Tomography (ERT), which concerns subsurface inference from surface potential measurements, see for example [28] which discussed the problem as early as 1933. Physically the main difference between EIT and ERT is that alternating current is typically used for the former, and direct current for the latter. Additionally, due to the scale of ERT, it is a reasonable approximation to model the electrodes as points, rather than using the CEM. Another difference between the two is that the relative contrast between the conductivities of different media are typically higher in subsurface applications than medical applications, which permits the approximation of the problem by a network of resistors in some cases [37]. Nonetheless, the Bayesian theory and MCMC methodology introduced here will be useful for the ERT problem as well as the EIT problem.

Statistical, in particular Bayesian, approaches to EIT inversion have previously been studied, for example in [18, 19, 31]. In [18], the authors prove certain regularity of the forward map associated with the CEM, formulate the Bayesian inverse problem in terms of the discretised model, and investigate the effect of different priors on reconstruction and behaviour of the posterior. The paper [31] focuses on Whittle-Matèrn priors, using EIT and ERT as examples for numerical simulation. The paper [19] presents a regularised version of the inverse problem, which admits a Bayesian interpretation.

The Bayesian approach to inverse problems, especially in infinite dimensions, is a relatively new technique. Two approaches are typically taken: discretise first and use finite dimensional Bayesian techniques, or apply the Bayesian methodology directly on function space before discretising. The former approach is outlined in [20]. The latter approach for linear problems has been studied in [12, 32, 33, 34]. More recently, this approach has been applied to nonlinear problems [10, 29, 30, 41]. It is this approach that we will be taking in this paper.

1.2. Our Contribution

The main contributions in our paper are as follows:

1. This is the first paper to give a rigorous Bayesian formulation of EIT in infinite dimensions.

2. This setting leads to proof that the posterior measure is Lipschitz in the data, with respect to the Hellinger metric, for all three priors studied; further stability properties of the posterior with respect to perturbations, such as numerical approximation of the forward model, may be proved similarly.

3. We employ a variety of prior models, based on the assumption that the underlying conductivity we wish to recover is binary. We initially look at log-Gaussian priors, before focusing on priors which enforce the binary field property. These binary field priors include both single star-shaped inclusions parametrised by their centre and by a radial function [8], and arbitrary geometric interfaces between the two conductivity values defined via level set functions [17].

4. Numerical results using state of the art MCMC demonstrate the importance of the prior choice for accurate reconstruction in the severely underdetermined inverse problems arising in EIT.

1.3. Organisation of the Paper

In section 2 we describe the forward map associated with the EIT problem, and prove relevant regularity properties. In section 3 we formulate the inverse problem rigorously and describe our three prior models. We then prove existence and well-posedness of the posterior distribution for each of these choices of prior. In section 4 we present results of numerical MCMC simulations to investigate the effect of the choice of prior on the recovery of certain binary conductivity fields. We conclude in section 5.

2. The Forward Model

In subsection 2.1 we describe the complete electrode model for EIT as given in [40]. In subsection 2.2 we give the weak formulation of this model, stating assumptions required for the quoted existence and uniqueness result. Then in subsection 2.3 we define the forward map in terms of this model, and prove that this map is continuous with respect to both uniform convergence and convergence in measure.

2.1. Problem Statement

Let , , be a bounded open set with smooth boundary representing a body, with conductivity . A number of electrodes are attached to the surface of the body. We treat these as subsets of the boundary , and assume that they have contact impedances . A current stimulation pattern is applied to the electrodes. Then the electric potential within the body and boundary voltages on are modelled by the following partial differential equation.

 (1) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−∇⋅(σ(x)∇v(x))=0x∈D∫elσ∂v∂ndS=Ill=1,…,Lσ(x)∂v∂n(x)=0x∈∂D∖⋃Ll=1elv(x)+zlσ(x)∂v∂n(x)=Vlx∈el,l=1,…,L

This model was first proposed in [40]; a derivation can be found therein. Note that the inputs to this forward model are the conductivity , input current and contact impedances . The solution comprises the function and the vector of voltages. Also note that solutions to this equation are only defined up to addition of a constant: if solves the equation, then so does for any . This is because it is necessary to choose a reference ground voltage.

2.2. Weak Formulation

We first define the space in which the solution to equation (1) will live. Using the notation of [40] we set

 H =H1(D)⊕RL, ∥(v,V)∥2H =∥v∥2H1(D)+∥V∥2ℓ2 =∥v∥2L2(D)+∥∇v∥2L2(D)+∥V∥2ℓ2.

Since solutions are only defined up to addition of a constant, we define the quotient space by

 ˙H =H/R, ∥(v,V)∥˙H =infc∈R∥(v−c,V−c)∥H.

We will often use the notation for brevity. It is more convenient to equip with an equivalent norm, as stated in the following lemma from [40]:

Lemma 2.1.

Define by

 ∥v′∥2∗=∥∇v∥2L2(D)+L∑l=1∫el|v(x)−Vl|2dS.

Then and are equivalent.

We can now state the weak formulation of the problem as derived in [40]. For this let .

Proposition 2.2.

Let and be defined by

 B(v′,w′;σ) =∫Dσ∇v⋅∇wdx+L∑l=11zl∫el(v−Vl)(w−Wl)dS, r(w′) =L∑l=1IlWl.

Then if is a strong solution to the problem (1), it satisfies

 (2) B(v′,w′;σ)=r(w′)for all% w′∈˙H.

We will use the weak formulation (2) to define our forward map for the complete electrode model (1). In order to guarantee a solution to this problem, we make the following assumptions.

Assumptions 2.3.

The conductivity , the contact impedances and the current stimulation pattern satisfy

1. ;

2. ;

3. .

Under these assumptions, existence of a unique solution to (2) is proved in [40] and stated here for convenience:

Proposition 2.4.

Let Assumptions 2.3 hold, then there is a unique solving (2). We may, without loss of generality, choose the element of the equivalence class of solutions to be that which satisfies

 (3) L∑l=1Vl=0.
Remark 2.5.

Assumptions 2.3 (i) and (ii) are to ensure coercivity and boundedness of . Assumption 2.3 (iii) is necessary for continuity of , and physically may be thought of as a conservation of charge condition. Choosing a solution from the equivalence class corresponds to choosing a reference ground voltage.

2.3. Continuity of the Forward Map

In what follows we will restrict to the set of admissible conductivities, which is defined as follows.

Definition 2.6.

A conductivity field is said to be admissible if

1. there exists , open disjoint subsets of for which ;

2. ; and

3. there exist such that for all .

The set of all such conductivities will be denoted .

Note that any will satisfy Assumptions 2.3(i). Assume that the current stimulation pattern and contact impedances are known and satisfy Assumptions 2.3(ii)-(iii). Then we may define the solution map to be the unique solution to (2) satisfying (3). The above existence and uniqueness result tells us that this map is well-defined.

In [18] it is shown that is Fréchet differentiable when we equip with the supremum norm. Though this is a strong result, this choice of norm is not appropriate for all of the conductivities that we will be considering. We hence establish the following continuity result.

Proposition 2.7.

Fix a current stimulation pattern and contact impedances satisfying Assumptions 2.3. Define the solution map as above. Let and let be such that either

1. converges to uniformly; or

2. converges to in (Lebesgue) measure, and there exist such that for all and , .

Then .

Proof.

Define the maps and as in Lemma 2.2, but on rather than . Then denoting and , we have for all ,

 B(v′ε,w′;σε)=r(w′),B(v′,w′;σ)=r(w′).

It follows that

 0 =B(v′ε,w′;σε)−B(v′,w′;σε)+B(v′,w′;σε)−B(v′,w′;σ) =∫Dσε∇(vε−v)⋅∇wdx+L∑l=11zl∫el((vε−v)−(Vεl−Vl))(w−Wl)dS +∫D(σε−σ)∇v⋅∇wdx.

Letting , we see that

 ∫Dσε|∇(vε−v)|2dx+L∑l=11zl∫el((vε−v) −(Vεl−Vl))2dS ≤∫D|σε−σ||∇v⋅∇(vε−v)|dx.

In both cases (i) and (ii), we have that is bounded uniformly below by a positive constant. Hence for small enough , the left hand side above may be bounded below by . We then have by Cauchy-Schwarz

 ∥v′ε−v′∥2∗ ≤C∫D|σε−σ||∇v⋅∇(vε−v)|dx ≤C(∫D|σε−σ|2|∇v|2dx)1/2⋅∥∇(vε−v)∥L2 (4) ≤C(∫D|σε−σ|2|∇v|2dx)1/2⋅∥v′ε−v′∥∗ (5) ≤C∥σε−σ∥∞∥∇v∥L2∥v′ε−v′∥∗.

If uniformly, we deduce from (5) that and the result follows. If in measure, then since , it follows that the integrand in (4) tends to zero in measure, see for example Corollary 2.2.6 in [6]. Since is assumed to be uniformly bounded, the integrand is dominated by a scalar multiple of the integrable function . We claim that this implies that the integrand tends to zero in . Suppose not, and denote the integrand . Then there exists and a subsequence such that for all . This subsequence still converges to zero in measure, and so admits a further subsequence that converges to zero almost surely. An application of the dominated convergence theorem leads to a contradiction, hence we deduce that tends to zero in and the result follows. ∎

Denote the projection , . The following lemma shows that the above result still holds if we replace by .

Corollary 2.8.

Let the assumptions of Proposition 2.7 hold. Then

 |Π∘M(σε)−Π∘M(σ)|ℓ2→0.
Proof.

We show that there exists such that for all with , . By the equivalence of and , Lemma 2.1, we have

 ∥(v,V)∥∗≥Cinfc∈R(∥v−c∥H1+|V−c|ℓ2)≥Cinfc∈R|V−c|ℓ2

The infimum on the right-hand side is attained at

 c=1LL∑l=1Vl=0.

Then by Proposition 2.7, we have

 0≤|Π∘M(σε)−Π∘M(σ)|ℓ2≤∥M(σε)−M(σ)∥∗→0

3. The Inverse Problem

We are interested in the inverse problem of determining the conductivity field from measurements of the voltages on the boundary, for a variety of input currents on the boundary. To this end we introduce the following version of Ohm’s law. Observe that the mapping , taking the current stimulation pattern to the solution of (2), is linear. Then given a conductivity field , there exists a resistivity matrix such that the boundary voltage measurements arising from the solution of the forward model are related to via

 V(σ)=R(σ)I

By applying several different current stimulation patterns we should be able to infer more about the conductivity . Note however that since the mapping is linear, only linearly independent stimulation patterns will provide more information111If there is noise on the measurements, additional linearly dependent observations can be made to effectively reduce the noise level on the original measurements. We can assume that this has been done and scale the noise appropriately.. Since we have the conservation of charge condition on , there are at most linearly independent patterns we can use.

Assume that linearly independent current patterns , , are applied, and noisy measurements of are made:

 yj=V(j)+ηj,ηj∼N(0,Γ0) iid.

We have

 yj=Gj(σ)+ηj

where . Concatenating these observations, we write

 y=G(σ)+η, η∼N(0,Γ)

where . The inverse problem is then to recover the conductivity field from the data . This problem is highly ill-posed: the data is finite dimensional, yet we wish to recover a function which, typically, lies in an infinite dimensional space. We take a Bayesian approach by placing a prior distribution on . The choice of prior may have significant effect on the resulting posterior distribution, and different choices of prior may be more appropriate depending upon the prior knowledge of the particular experimental set-up under consideration.

In subsection 3.1 we outline three different families of prior models, and show the appropriate regularity of the forward maps arising from them. In subsection 3.2 we describe the likelihood and posterior distribution formally, before rigorously proving that the posterior distribution exists and is Lipschitz with respect to the data in the Hellinger metric.

3.1. Choices of Prior

In this section we consider three priors, labelled by , defined by functions which map draws from prior measures on the Banach spaces to the space of conductivities . Our prior conductivity distributions will then be the pushfoward of the prior measures by these maps . We describe these maps, and establish continuity properties of them needed for the study of the posterior later. In what follows, the space of continuous functions on will be equipped with the supremum norm and the corresponding Borel -algebra.

3.1.1. Log-Gaussian prior

We first consider the simple case that the coefficient is given by the exponential of a continuous function. Let be defined by . Then it is easily seen that does indeed map into . Furthermore, since is bounded, if and is a sequence such that , then .

In this case, we will take our prior measure on to be a non-degenerate Gaussian measure on . Note that the push forward of a Gaussian measure by is a log-Gaussian measure.

Example 3.1.

Consider the case . Suppose that is drawn from a Gaussian measure . Typical samples from are shown below222Given a measure on and a measurable map between measurable spaces, denotes the pushforward of by , i.e. the measure on given by for all .. The covariance is chosen such that the samples almost surely have regularity for all , where from left to right respectively. Here the samples are generated on and then restricted to , for computational simplicity.

3.1.2. Star-shaped prior

We now consider star-shaped inclusions, that is, inclusions parametrised by their centre and a radial function. These were studied in two-dimensions in the paper [8] to parametrise domains for a Bayesian inverse shape scattering problem. In [8] the authors prove well-posedness of the inverse problem in an infinite dimensional setting through the use of shape derivatives and Riesz-Fredholm theory.

Let , and . Let be the continuous function representing the mapping from Cartesian to angular polar coordinates. Define the mapping by

 A(r,x0)={x∈D∣∣|x−x0|≤r(h(x−x0))}

where is the space of continuous periodic functions on . Then describes the set of points in which lie within the closed surfaceparametrised in polar coordinates centred at by

 Γ(Θ)=(Θ,r(Θ)),Θ∈Rd−1.

In two dimensions, we have and the mapping is given by

 h(x,y)=atan2(y,x)≡2arctan(y√x2+y2+x)

where atan2 is the two-parameter inverse tangent function.

In three dimensions, we have and the mapping is given by

 h(x,y,z)=(atan2(y,x),arccot(z√x2+y2)).

Similar expressions for exist in higher dimensions, though for applications we are only interested in the case .

Define now the map by

 F2(r,x0) =u+\mathds1A(r,x0)+u−\mathds1D∖A(r,x0) =(u+−u−)\mathds1A(r,x0)+u−,

where are the scalar conductivity values. Again it can easily be seen that does indeed map into . We claim that this map is continuous in the following sense:

Proposition 3.2.

Define the map as above. Let and let be Lipschitz continuous.

1. Suppose that is a sequence of functions such that. Then in measure333A sequence of functions , , is said to converge in measure to a function if for all , . Here denotes the Lebesgue measure of a set ..

2. Suppose that is a sequence of points such that . Then in measure.

3. Let , be as above. Then in measure.

Proof.

In order to show that a sequence of functions converges to in measure, it suffices to show that there exists a sequence of sets with such that . Then for each we have

 |{x∈D||fε(x)−f(x)|>δ}|≤|{x∈D||fε(x)−f(x)|≠0}|≤|Zε|→0.
1. Fix the centre . Denote . Let and let be a sequence of functions such that . Then there exists such that . By definition we then have

 r(x)−γ(ε)≤rε(x)≤r(x)+γ(ε)for all x∈D and ε>0.

It follows that we have the inclusions

 A(r−γ(ε))⊆A(rε)⊆A(r+γ(ε)), A(r−γ(ε))⊆A(r)⊆A(r+γ(ε)).

Let denote the symmetric difference. We deduce that

 A(rε)ΔA(r)⊆A(r+γ(ε))∖A(r−γ(ε)).

Now the right-hand side is given by

 A(r+γ(ε)) ∖A(r−γ(ε)) ={x∈D∣∣r(h(x−x0))−γ(ε)<|x−x0|≤r(h(x−x0))+γ(ε)}.

As , this set decreases to the boundary set

 ∂A(r)={x∈D∣∣|x−x0|=r(h(x−x0))}.

Since the graph of a continuous function has Lebesgue measure zero, we deduce that . It follows that

 limε→0|A(rε)ΔA(r)|=0.

To conclude, note that

 |F2(rε,x0)−F2(r,x0)|≤|u+−u−||\mathds1A(rε)−\mathds1A(rε)|=C\mathds1A(rε)ΔA(r).
2. Let be Lipschitz continuous. Denote . Let be a sequence of points such that . Note that we may write

 A(xε0) ={x∈D||x−xε0|≤r(h(x−xε0))} ={x∈Rd||x−xε0|≤r(h(x−xε0))}∩D =((xε0−x0)+{x∈Rd||x−x0|≤r(h(x−x0))})∩D =:((xε0−x0)+A(x0)∗)∩D.

By the distributivity of intersection over symmetric difference, we then have that

 A(xε0)ΔA(x0) =[((xε0−x0)+A(x0)∗)∩D]Δ[A(x0)∗∩D] =[((xε0−x0)+A(x0)∗)ΔA(x0)∗]∩D ⊆((xε0−x0)+A(x0)∗)ΔA(x0)∗.

Therefore, using Theorem 1 from [39], we see that

 |A(xε0)ΔA(x0)| ≤|((xε0−x0)+A(x0)∗)ΔA(x0)∗| ≤|xε0−x0|Hd−1(∂A(x0)∗)

where is the -dimensional Hausdorff measure. Since we assume that is Lipschitz, the surface area of the boundary of is finite, and so it follows that

 limε→0|A(xε0)ΔA(x0)|=0.

As before, we conclude by noting that

 |F2(r,xε0)−F2(r,x0)|≤|u+−u−||\mathds1A(xε0)−\mathds1A(x0)|=C\mathds1A(xε0)ΔA(x0).
3. We have that

 |F2(rε,xε0)−F2(r,x0)| ≤|F2(rε,xε0)−F2(r,xε0)|+|F2(r,xε0)−F2(r,x0)| ≤C(\mathds1A(rε,xε0)ΔA(r,xε0)+\mathds1A(r,xε0)ΔA(r,x0)) ≤C\mathds1[A(rε,xε0)ΔA(r,xε0)]∪[A(r,xε0)ΔA(r,x0)].

Now note that

 |A(rε,y0)ΔA(r,y0)|≤|A(rε,y0)∗ΔA(r,y0)∗|.

The right hand-side is independent of by translation invariance of the Lebesgue measure. By the same argument as part (i) we conclude that it tends to zero. We then have that

 |[A(rε,xε0)ΔA(r,xε0)] ∪[A(r,xε0)ΔA(r,x0)]| ≤|A(rε,xε0)ΔA(r,xε0)|+|A(r,xε0)ΔA(r,x0)| ≤supy0∈D|A(rε,y0)ΔA(r,y0)|+|A(r,xε0)ΔA(r,x0)|

which tends to zero by the discussion above and part (ii).

Remark 3.3.

Above we assumed that was Lipschitz continuous. This assumption is only used in the proof of part (ii) of the proposition. If the centre of the star-shaped region is known, this assumption may then be dropped to allow for rougher boundaries.

We need to choose a prior measure on . We equip with the norm and corresponding Borel -algebra. We assume that and are independent under the prior so that we may factor where is a measure on and is a measure on . We will assume that is such that for all balls .

Example 3.4.

Consider the case . Suppose that is drawn from a log-Gaussian measure on , and is drawn from . Note that . Typical samples from are shown below. The covariance of is chosen such that the samples almost surely have regularity for all , where from left to right respectively.

3.1.3. Level set prior

We finally consider the case where the inclusions can be described by a single level set function, as in [17]. Let and fix constants . Given , define by

 Di={x∈D|ci−1≤u(x)

so that and for , . Define also the level sets

 D0i=¯¯¯¯¯Di∩¯¯¯¯¯Di+1={x∈D|u(x)=ci},i=1,…,n−1.

Now given strictly positive functions , we define the map by

 F3(u)=n∑i=1fi\mathds1Di.

Since each is continuous and strictly positive on a compact set , they are uniformly bounded above and below by positive constants, and so does indeed map into .

In this paper we are primarily concerned with the case of binary fields, and constant above, however the theory in proved in the general case. We have the following result regarding continuity of this map, by the same arguments as in [17].

Proposition 3.5.

Define the map as above. Let be such that for . Suppose that is an approximating sequence of functions so that . Then in measure.

Proof.

Denote by and the sets as defined above associated with the approximating functions . We can write

 F3(uε)−F3(u)=n∑i=1n∑j=1