1 Introduction

Identifying conductivity in electrical impedance tomography

with total variation regularization

Michael HinzeEmail: michael.hinze@uni-hamburg.de,  Barbara.Kaltenbacher@aau.at,  quyen.tran@uni-hamburg.de Barbara KaltenbacherM. Hinze gratefully acknowledges support of the Lothar Collatz Center for Computing in Science at the University of HamburgB. Kaltenbacher gratefully acknowledges support of the Austrian Wissenschaftsfonds through grant FWF I2271 entitled ”Solving inverse problems without forward operators”T.N.T. Quyen gratefully acknowledges support of the Alexander von Humboldt-Foundation Tran Nhan Tam QuyenAuthor to whom any correspondence should be addressed

University of Hamburg, Bundesstraße 55, D-20146 Hamburg, Germany
Alpen-Adria-Universität Klagenfurt, Universitätsstraße 65-67, A-9020 Klagenfurt, Austria

Abstract: In this paper we investigate the problem of identifying the conductivity in electrical impedance tomography from one boundary measurement. A variational method with total variation regularization is here proposed to tackle this problem. We discretize the PDE as well as the conductivity with piecewise linear, continuous finite elements. We prove the stability and convergence of this technique. For the numerical solution we propose a projected Armijo algorithm. Finally, a numerical experiment is presented to illustrate our theoretical results.

Key words and phrases: Conductivity identification, electrical impedance tomography, total variation regularization, finite element method, Neumann problem, Dirichlet problem, ill-posed problems.

AMS Subject Classifications: 65N21; 65N12; 35J25; 35R30

1 Introduction

Let be an open bounded connected domain in with polygonal boundary and be given. We consider the following elliptic boundary value problem

 −∇⋅(q∇Φ) =f~{}in~{}Ω, (1.1) q∇Φ⋅→n =j†~{}on~{}∂Ω~{}and~{} (1.2) Φ =g†~{}on~{}∂Ω, (1.3)

where is the unit outward normal on .

The system (1.1)–(1.3) is overdetermined, i.e. if the Neumann and Dirichlet boundary conditions and the conductivity

 q∈Q:={q∈L∞(Ω) ∣∣ q–≤q(x)≤¯¯¯q~{}a.e. in~{}Ω} (1.4)

are given, then there may be no satisfying this system. Here and are some given positive constants.

In this paper we assume that the system is consistent and our aim is to identify the conductivity and the electric potential in the system (1.1)–(1.3) from current and voltage i.e., Neumann and Dirichlet measurements at the boundary of the exact satisfying

 ∥∥jδ−j†∥∥H−1/2(∂Ω)+∥∥gδ−g†∥∥H1/2(∂Ω)≤δ~{}with~{% }δ>0.

Note that using the topology for the data is natural from the point of view of solution theory for elliptic PDEs but unrealistic with regard to practical measurements. We will comment in this issue in Remark 2.2 below.

For the purpose of conductivity identification — a problem which is very well known in literature and practice as electrical impedance tomography EIT, see below for some references — we simultaneously consider the Neumann problem

 −∇⋅(q∇u)=f~{}in~{}Ω~{}and~{}q∇u⋅→n=jδ~{}on~{}∂Ω (1.5)

and the Dirichlet problem

 −∇⋅(q∇v)=f~{}in~{}Ω~{}and~{}v=gδ~{}on~{}∂Ω (1.6)

and respectively denote by , the unique weak solutions of the problems (1.5), (1.6), which depend nonlinearly on , where is normalized with vanishing mean on the boundary. We adopt the variational approach of Kohn and Vogelius in [30, 31, 32] to the identification problem. In fact, for estimating the conductivity from the observation of the exact data , we use the functional

 Jδ(q):=∫Ωq∇(Nqjδ−Dqgδ)⋅∇(Nqjδ−Dqgδ)dx.

For simplicity of exposition we restrict ourselves to the case of just one Neumann-Dirichlet pair, while the approach described here can be easily extended to multiple measurements , see also Example 5.3 below. It is well-known that such a finite number of boundary data in general only allows to identify conductivities taking finitely many different values in the domain , see, e.g., [2].

Indeed, we are interested in estimating such piecewise constant conductivities and therefore use total variation regularization, i.e., we consider the minimization problem

 minq∈QadJδ(q)+ρ∫Ω|∇q|, (1.7)

where is the admissible set of the sought conductivities, is the space of all functions with bounded total variation (see §2.1 for its definition) and is the regularization parameter, and consider a minimizer of (1.7) as reconstruction.

For each let and be corresponding approximations of and in the finite dimensional space of piecewise linear, continuous finite elements and denote a minimizer of the discrete regularized problem corresponding to (1.7), i.e. of the following minimization problem

 (1.8)

with and being a positive functional of the mesh size satisfying .

In Section 4 we will show the stability of approximations for fixed positive . Furthermore as and with an appropriate a priori regularization parameter choice , there exists a subsequence of converging in the -norm to a total variation-minimizing solution defined by

 q†∈argmin{q∈Qad | Nqj†=Dqg†}∫Ω|∇q|.

In particular, if is uniquely defined, then this convergence holds for the whole sequence . The corresponding state sequences and converge in the -norm to solving the system (1.1)-(1.3). Finally, for the numerical solution of the discrete regularized problem (1.8), in Section 5 we employ a projected Armijo algorithm. Numerical results show the efficiency of the proposed method and illustrate our theoretical findings.

We conclude this introduction with a selection of references from the vast literature on EIT, which has evolved to a highly relevant imaging and diagnostics tool in industrial and medical applications and has attracted great attention of many scientists in the last few decades.

To this end, for any fixed we define the Neumann-to-Dirichlet map , by

 j↦Λqj:=γNqj,

where is the Dirichlet trace operator. Calderón in 1980 posed the question whether an unknown conductivity distribution inside a domain can be determined from an infinite number of boundary observations, i.e. from the Neumann-to-Dirichlet map :

 p,q∈Q⊂L∞(Ω)~{}with~{}Λp=Λq⇒p=q ? (1.9)

Calderón did not answer his question (1.9); however, in [15] he proved that the problem linearized at constant conductivities has a unique solution. In dimensions three and higher Sylvester and Uhlmann [41] proved the unique identifiability of a -smooth conductivity. Päivärinta el al. [37] and Brown and Torres [12] established uniqueness in the inverse conductivity problem for -smooth conductivities with and , respectively. In the two dimensional setting, Nachman [34] and Brown and Uhlmann [13] proved uniqueness results for conductivities which are in with and with , respectively. Finally, in 2006 the question (1.9) has been answered to be positive by Astala and Päivärinta [3] in dimension two. For surveys on the subject, we refer the reader to [10, 17, 20, 33, 43] and the references therein.

Although there exists a large number of papers on the numerical solution of the inverse problems of EIT, among these also papers considering the Kohn-Vogelius functional (see, e.g., [28, 29]) and total variation regularization (see, e.g., [21, 36]), we have not yet found investigations on the discretization error in a combination of both functionals for the fully nonlinear setting, a fact which motivated the research presented in this paper.

Throughout the paper we use the standard notion of Sobolev spaces , , , etc from, for example, [1]. If not stated otherwise we write instead of .

2 Problem setting and preliminaries

2.1 Notations

Let us denote by

 γ:H1(Ω)→H1/2(∂Ω)

the continuous Dirichlet trace operator while

 γ−1:H1/2(∂Ω)→H1(Ω)

is the continuous right inverse operator of , i.e. for all . With (with a slight abuse of notation) in (1.1) being given, let us denote

 cf:=(f,1),

where the expression denotes the value of the functional at . We also denote

where the notation stands for the value of the functional at . Similarly, we denote

while is the closed subspace of consisting of all functions with zero mean on the boundary, i.e.

 H1⋄(Ω):={u∈H1(Ω) ∣∣ ∫∂Ωγu=0}.

Let us denote by the positive constant appearing in the Poincaré-Friedrichs inequality (see, for example, [38])

 CΩ⋄∫Ωφ2≤∫Ω|∇φ|2~{}for all~{}φ∈H1⋄(Ω). (2.1)

Then for all defined by (1.4), the coercivity condition

 ∥φ∥2H1(Ω)≤1+CΩ⋄CΩ⋄∫Ω|∇φ|2≤1+CΩ⋄CΩ⋄q–∫Ωq∇φ⋅∇φ (2.2)

holds for all . Furthermore, since , the inequality (2.2) remains valid for all .

Finally, for the sake of completeness we briefly introduce the space of functions with bounded total variation; for more details one may consult [4, 24]. A scalar function is said to be of bounded total variation if

 TV(q):=∫Ω|∇q|:=sup{∫Ωqdiv Ξ ∣∣ Ξ∈C1c(Ω)d, |Ξ(x)|∞≤1, x∈Ω}<∞.

Here denotes the -norm on defined by and the space of continuously differentiable functions with compact support in . The space of all functions in with bounded total variation is denoted by

 BV(Ω)={q∈L1(Ω) ∣∣ ∫Ω|∇q|<∞}

which is a Banach space with the norm

 ∥q∥BV(Ω):=∥q∥L1(Ω)+∫Ω|∇q|.

Furthermore, if is an open bounded set with Lipschitz boundary, then .

2.2 Neumann operator, Dirichlet operator and Neumann-to-Dirichlet map

2.2.1 Neumann operator

We consider the following Neumann problem

 −∇⋅(q∇u)=f~{}in~{}Ω~{}and~{}q∇u⋅→n=j~{}on~{}∂Ω. (2.3)

By the coercivity condition (2.2) and the Riesz representation theorem, we conclude that for each and there exists a unique weak solution of the problem (2.3) in the sense that and satisfies the identity

 ∫Ωq∇u⋅∇φ=⟨j,γφ⟩+(f,φ) (2.4)

for all . By the imposed compatibility condition , i.e.

 ⟨j,1⟩+(f,1)=0, (2.5)

and the fact that , equation (2.4) is satisfied for all . Furthermore, this solution satisfies the following estimate

 ∥u∥H1(Ω) ≤1+CΩ⋄CΩ⋄q–(∥γ∥L(H1(Ω),H1/2(∂Ω))∥j∥H−1/2(∂Ω)+∥f∥H−1(Ω)) ≤CN(∥j∥H−1/2(∂Ω)+∥f∥H−1(Ω)), (2.6)

where

 CN:=1+CΩ⋄CΩ⋄q–max(1,∥γ∥L(H1(Ω),H1/2(∂Ω))).

Then for any fixed we can define the Neumann operator

 N:Q→H1⋄(Ω)% ~{}with~{}q↦Nqj

which maps each to the unique weak solution of the problem (2.3).

Remark 2.1.

We note that the restriction instead of preserves the compatibility condition (2.5) for the pure Neumann problem. In case this condition fails, the strong form of the problem (2.3) has no solution. This is the reason why we require . However, its weak form, i.e. the variational equation (2.4), attains a unique solution independently of the compatibility condition. By working with the weak form only, all results in the present paper remain valid for .

2.2.2 Dirichlet operator

We now consider the following Dirichlet problem

 −∇⋅(q∇v)=f~{}in~{}Ω~{}and~{}v=g~{}on~{}∂Ω. (2.7)

For each and , by the coercivity condition (2.2), the problem (2.7) attains a unique weak solution in the sense that , and satisfies the identity

 ∫Ωq∇v⋅∇ψ=(f,ψ) (2.8)

for all . We can rewrite

 v=v0+G, (2.9)

where and is the unique solution to the following variational problem

 ∫Ωq∇v0⋅∇ψ=(f,ψ)−∫Ωq∇G⋅∇ψ

for all . Since

 ∥G∥H1(Ω)≤∥∥γ−1∥∥L(H1/2(∂Ω),H1(Ω))∥g∥H1/2(∂Ω),

we thus obtain the priori estimate

 ∥v∥H1(Ω) ≤∥v0∥H1(Ω)+∥G∥H1(Ω) ≤1+CΩ⋄CΩ⋄q–∥f∥H−1(Ω)+1+CΩ⋄CΩ⋄q–¯¯¯q∥G∥H1(Ω)+∥G∥H1(Ω) ≤CD(∥g∥H1/2(∂Ω)+∥f∥H−1(Ω)), (2.10)

where

The Dirichlet operator is for any fixed defined as

 D:Q→H1(Ω)~{}with~{}q↦Dqg

which maps each to the unique weak solution of the problem (2.7).

2.2.3 Neumann-to-Dirichlet map

For any fixed we can define the Neumann-to-Dirichlet map

 Λq:H−1/2−cf(∂Ω) →H1/2⋄(∂Ω) j ↦Λqj:=γNqj.

Since

 ∫Ωq∇Nqj⋅∇ψ=(f,ψ)

for all , in view of (2.8) we conclude that

 Λqj=g~{}if and only if~{}Nqj=Dqg.

2.3 Identification problem

The inverse problem is stated as follows.

Given , with , find .

In other words, the problem of interest is, given , and a single Neumann-Dirichlet pair , to find and such that the system (1.1)–(1.3) is satisfied in the weak sense.

2.4 Total variation regularization

Assume that is the measured data of the exact boundary values with

 ∥∥jδ−j†∥∥H−1/2(∂Ω)+∥∥gδ−g†∥∥H1/2(∂Ω)≤δ (2.11)

for some measurement error . Our problem is now to reconstruct the conductivity from this perturbed data . For this purpose we consider the cost functional

 Jδ(q):=∫Ωq∇(Nqjδ−Dqgδ)⋅∇(Nqjδ−Dqgδ), (2.12)

where and is the unique weak solutions of the problems (2.3) and (2.7), respectively, with in (2.3) and in (2.7) being replaced by and . Furthermore, to estimate the possibly discontinuous conductivity, we here use the total variation regularization (cf., e.g., [14, 21, 22]), i.e., we consider the minimization problem

 minq∈QadΥρ,δ(q):=minq∈Qad(Jδ(q)+ρ∫Ω|∇q|),

where

 Qad:=Q∩BV(Ω)

is the admissible set of the sought conductivities.

Remark 2.2.

The noise model (2.11) is to some extent an idealized one, since in practice, measurement precision might be different for the current and the voltage , and, more importantly, it will be first of all be given with respect to some norm (e.g., corresponding to normally and to uniformly distributed noise) rather than in . While the Neumann data part is unproblematic, by continuity of the embedding of in for , we can obtain an version of the originally Dirichlet data e.g. by Tikhonov regularization (cf. [22] and the references therein) as follows. For simplicity, we restrict ourselves to the Hilbert space case and assume that we have measurements such that

 ∥~gδg−g†∥L2(∂Ω)≤δg

Tikhonov regularization applied to the embedding operator amounts to finding a minimizer of

 ming∈H1/2(∂Ω)∥Kg−~gδg∥2L2(∂Ω)+α∥g∥2H1/2(∂Ω)

where we use

 ∥g∥H1/2(∂Ω):=∥γ−1g∥H1(Ω)=(∫Ω(|∇γ−1g|2+|γ−1g|2)dx)1/2

as a norm on . The first order optimality conditions for this quadratic minimization problem yield

 ∫∂Ωϕ(gδgα−~gδg)ds+α∫Ω(∇γ−1gδgα⋅∇γ−1ϕ+γ−1gδgαγ−1ϕ)dx=0 for % all ϕ∈H1/2(∂Ω),

which is equivalent to

 ∫∂Ωγφ(γw−~gδg)ds+α∫Ω(∇w⋅∇φ+wφ)dx=0 for all φ∈H1(Ω),

for , i.e., the weak form of the Robin problem

 {−Δw+w=0~{}in~{}Ω,α∇w⋅→n+w=~gδg~{}on~{}∂Ω. (2.13)

Thus, according to well-known results from regularization theory (cf., e.g. [22]), the smoothed version (where weakly solves (2.13)) of converges to as tends to zero, provided the regularization parameter is chosen appropriately. The latter can, e.g., be done by the discrepancy principle, where is chosen such that

 ∥Kgδgα−~gδg∥2L2(∂Ω)=∫∂Ω|gδgα−~gδg|2dx∼δ2g.

We also wish to mention the complete electrode model cf., e.g., [40], which fully takes into account the fact that current and voltage are typically not measured pointwise on the whole boundary, but via a set of finitely many electrodes with finite geometric extensions as well as contact impedances.

2.5 Auxiliary results

Now we summarize some useful properties of the Neumann and Dirichlet operators. The proof of the following result is based on standard arguments and therefore omitted.

Lemma 2.3.

Let be fixed.

(i) The Neumann operator is continuously Fréchet differentiable on the set . For each the action of the Fréchet derivative in direction denoted by is the unique weak solution in to the Neumann problem

 −∇⋅(q∇ηN)=∇⋅(ξ∇Nqj)~{}in~{}Ω~{}and~{}q∇ηN⋅→n=−ξ∇Nqj⋅→n~{}on~{}∂Ω

in the sense that the identity

 ∫Ωq∇ηN⋅∇φ=−∫Ωξ∇Nqj⋅∇φ (2.14)

holds for all . Furthermore, the following estimate is fulfilled

 ∥ηN∥H1(Ω)≤(1+CΩ⋄)CNCΩ⋄q–(∥j∥H−1/2(∂Ω)+∥f∥H−1(Ω))∥ξ∥L∞(Ω). (2.15)

(ii) The Dirichlet operator is continuously Fréchet differentiable on the set . For each the action of the Fréchet derivative in direction denoted by is the unique weak solution in to the Dirichlet problem

 −∇⋅(q∇ηD)=∇⋅(ξ∇Dqg)~{}in~{}Ω~{}and~{}ηD=0~{}on~{}∂Ω

in the sense that it satisfies the equation

 ∫Ωq∇ηD⋅∇ψ=−∫Ωξ∇Dqg⋅∇ψ

for all . Furthermore, the following estimate is fulfilled

 ∥ηD∥H1(Ω)≤(1+CΩ⋄)CDCΩ⋄q–(∥g∥H1/2(∂Ω)+∥f∥H−1(Ω))∥ξ∥L∞(Ω).
Lemma 2.4.

If the sequence converges to in the -norm, then and for any fixed the sequence converges to in the -norm. Furthermore, there holds

 limn→∞Jδ(qn)=Jδ(q),

where the functional is defined in (2.12).

Proof.

Since converges to in the -norm, up to a subsequence we assume that it converges to a.e. in , which implies that . For all we infer from (2.4) that

 ∫Ωqn∇Nqnjδ⋅∇φ =⟨jδ,γφ⟩+(f,φ)=∫Ωq∇Nqjδ⋅∇φ

and so that

 ∫Ωqn∇(Nqnjδ−Nqjδ)⋅∇φ=∫Ω(q−qn)∇Nqjδ⋅∇φ. (2.16)

Taking , by (2.2), we get

 CΩ⋄q–1+CΩ⋄∥∥Nqnjδ−Nqjδ∥∥2H1(Ω) ≤∫Ωqn∇(Nqnjδ−Nqjδ)⋅∇(Nqnjδ−Nqjδ) =∫Ω(q−qn)∇Nqjδ⋅∇(Nqnjδ−Nqjδ) ≤(∫Ω|q−qn|2∣∣∇Nqjδ∣∣2)1/2(∫Ω∣∣∇(Nqnjδ−Nqjδ)∣∣2)1/2

and so that

 ∥∥Nqnjδ−Nqjδ∥∥H1(Ω)≤1+CΩ⋄CΩ⋄q–(∫Ω|q−qn|2∣∣∇Nqjδ∣∣2)1/2.

Hence, by the Lebesgue dominated convergence theorem, we deduce from the last inequality that

 limn→∞∥∥Nqnjδ−Nqjδ∥∥H1(Ω)=0. (2.17)

Similarly to (2.16), we also get

 ∫Ωqn∇(Dqngδ−Dqgδ)⋅∇ψ=∫Ω(q−qn)∇Dqgδ⋅∇ψ

for all . Since , taking in the last equation, we also obtain the limit

 limn→∞∥Dqngδ−Dqgδ∥H1(Ω)=0. (2.18)

Next, we rewrite the functional as follows

 Jδ(qn) =∫Ωqn∇Nqnjδ⋅∇Nqnjδ−2∫Ωqn∇Nqnjδ⋅∇Dqngδ+∫Ωqn∇Dqngδ⋅∇Dqngδ =⟨jδ,γNqnjδ⟩+(f,Nqnjδ)−2(⟨jδ,gδ⟩+(f,Dqngδ))+∫Ωqn∇Dqngδ⋅∇Dqngδ (2.19)

and, by (2.17)–(2.18), have that

 (2.20)

as tends to . We now consider the difference

 ∫Ωqn∇Dqngδ⋅∇Dqngδ =∫Ωqn∇(Dqngδ−Dqgδ)⋅∇(Dqngδ+Dqgδ)−∫Ω(q−qn)∇Dqgδ⋅∇Dqgδ

and note that

 ∫Ω(q−qn)∇Dqgδ⋅∇Dqgδ→0

as goes to , by the Lebesgue dominated convergence theorem. Furthermore, then applying the Cauchy-Schwarz inequality, we also get that

 ∣∣∣∫Ωqn ∇(Dqngδ−Dqgδ)⋅∇(Dqngδ+Dqgδ)∣∣∣ ≤¯¯¯q(∫Ω|∇(Dqngδ−Dqgδ)|2)1/2(∫Ω|∇(Dqngδ+Dqgδ)|2)1/2 ≤¯¯¯q∥Dqngδ−Dqgδ∥H1(Ω)(∥Dqngδ∥H1(Ω)+∥Dqgδ∥H1(Ω))→0

as approaches , here we used (2.2.2) and (2.18). We thus obtain that

 ∫Ωqn∇Dqngδ⋅∇Dqngδ→∫Ωq∇Dqgδ⋅∇Dqgδ (2.21)

as tends to . Then we deduce from (2.5)–(2.21) that

 limn→∞Jδ(qn) =Jδ(q),

which finishes the proof. ∎

Lemma 2.5 ([24]).

(i) Let be a bounded sequence in the -norm. Then a subsequence which is denoted by the same symbol and an element exist such that converges to in the -norm.

(ii) Let be a sequence in converging to in the -norm. Then and

 ∫Ω|∇q|≤liminfn→∞∫Ω|∇qn|. (2.22)

We mention that equality need not be achieved in (2.22