One-step iterative reconstruction of conductivity inclusion via the concept of topological derivative

# One-step iterative reconstruction of conductivity inclusion via the concept of topological derivative

## Abstract

We consider an inverse problem of location identification of small conductivity inhomogeneity inside a conductor via boundary measurements which occurs in the EIT (Electrical Impedance Tomography). For this purpose, we derive topological derivative by applying the asymptotic formula for steady state voltage potentials in the existence of conductivity inclusion of small diameter. Using this derivative, we design only one-step iterative location search algorithm of small conductivity inhomogeneity completely embedded in the homogeneous domain by solving an adjoint problem. Numerical experiments presented for showing the feasibility of proposed algorithm.

###### keywords:
Topological derivative, Inverse conductivity problems, Asymptotic formula, Numerical experiments
1

## 1 Introduction

Electrical Impedance Tomography (EIT), an technique for imaging distribution of the conductivity distribution of conducting objects from surface electrical measurements, is an interesting and important problem which is arising in physics, medical science, material engineering, and so on, all domains highly related with human life. Related works can be found in (5); (6); (13); (15); (18); (19); (37) and reference therein. However, due to the ill-posedness and inherent non-linearity of problem, it remains a challenging research area. In recent research, various remarkable imaging techniques proposed, many of them based on the Newton-type iteration method, refer to (1); (18); (36) and references therein. For a successful application of such method, a good initial guess that is close to the unknown target is essentially required. If one proceed without a good initial guess, issue such as non-convergence, the occurrence of several minima, and large computational costs may arise. Moreover, iteration method often requires suitable regularization terms that highly depend on the specific problem at hand.

For this purpose, various non-iterative methods are also developed as an alternatives, e.g., a real-time algorithm for finding location of conductivity inhomogeneity (5); (6); (19); (37), simple pole method (20), and end-point location search algorithm of thin conductivity inhomogeneities (2); (3); (21).

The main purpose of this paper is to develop a fast, non-iterative imaging algorithm of small conductivity inclusion embedded in a homogeneous domain from boundary measurement by adopting famous topological derivative concept (14); (16); (39). By applying asymptotic expansion formula in the existence of conductivity inclusion of small diameter (5); (6); (9), we can easily derive desired derivative and design a fast (only requires one-step iteration procedure by solving an adjoint problem) imaging algorithm. Although, the result via topological derivative does not guarantee complete shape of unknown target, it will be a good initial guess of Newton-type iteration method.

This paper organized as follows. In section 2, we briefly review basic mathematical model of direct conductivity problem and introduce the asymptotic expansion formula due to the existence of small inclusion. In section 3, we apply asymptotic expansion formula in order to rigorously and easily derive the topological derivative and develop a one-step iterative imaging algorithm. In the following section 4 numerical experiments are shown. Section 6 contains a brief conclusion.

## 2 Mathematical survey on direct conductivity problems and asymptotic expansion formula

Let be a smooth, bounded domain that represents a homogeneous medium. We assume that this medium contains a conductivity inclusion with small diameter represented as

 D=z+εB,

where is some fixed bounded domain containing the origin that completely embedded in . Throughout this paper, we denote and are two-dimensional vectors and assume that does not touch the boundary , i.e., there exists positive constant such that

 dist(D,∂Ω)>h. (1)

Assume that every materials are fully characterized by their electrical conductivity. Let and denote the conductivity of the domain and inclusion , respectively. By using this notation, we can introduce the following piecewise constant conductivity:

 ^σ(x)={σ0forx∈Ω∖¯¯¯¯DσDforx∈D. (2)

In this paper, for the sake, we assume that .

Let be the steady state voltage potential in the presence of the inclusion , that is, the unique solution to

 ∇⋅(^σ(x)∇u(x))=0forx∈Ω (3)

with the Neumann boundary condition

 σ0∂u∂ν(x)=g(x)% forx∈∂Ω (4)

and with the compatibility condition

 ∫∂Ωu(x)dS(x)=0.

Here denote the unit outer normal to at and the function represents the applied boundary current satisfies the compatibility condition

 ∫∂Ωg(x)dS(x)=0.

Let us denote be the background potential induced by the current in the domain without , that is, the unique solution to

 σ0Δu0(x)=0forx∈Ω

with the Neumann boundary condition

 σ0∂u0∂ν(x)=g(x)forx∈∂Ω

and normalization condition to restore uniqueness

 ∫∂Ωu0(x)dS(x)=0.

Then asymptotic expansion formula in the presence of can be written as follows:

###### Theorem 2.1 (See (6); (9)).

Let be a bounded, smooth domain. Then for sufficiently small , the asymptotic expansion formula due to the presence of small inclusion can be expressed in terms of the magnitude of diameter as:

 u(y)−u0(y)=−πε2∇u0(x)⋅M(x)⋅∇N(x,y)+O(ε3) (5)

for each on . Here the remaining term is uniform in , denotes the Neumann function for the domain

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩σ0ΔN(x,y)=−δ(x,y)inΩ\@@LTX@noalign\vskip6.0ptplus2.0ptminus2.0pt\omitσ0∂N(x,y)∂ν(y)=−1|∂Ω|,∫∂ΩN(y,z)dS(y)=0on∂Ω,

and is the symmetric matrix associated with the inclusion and the conductivities and . Specially, if is a ball with radius , is given by

 M(x)=2σD−σ0σD+σ0(1001). (6)

## 3 Topological derivative and its structure

At this stage, we derive the topological derivative in order to establish an imaging algorithm of small conductivity inclusion . For a more detailed discussion about the topological derivative, we recommend research articles (7); (10); (12); (14); (25); (28); (29); (30); (39).

The topological derivative measures the influence of creating a small ball (or crack, etc.) with small radius at a certain point inside the domain – we denote as the such domain. If we assume that the boundaries and are sufficiently smooth, by considering a cost functional , the topological derivative is defined as

 dTD(z)=limr→0+D(Ω∖Σ)−D(Ω)ρ(r), (7)

where as . Notice that the function is mainly defined by geometrical factors of the created shape (here, ). With this definition (7), we have the asymptotic expansion:

 D(Ω∖Σ)=D(Ω)+ρ(r)dTD(z)+o(ρ(r)).

Suppose that contains a thin inclusion , , , be given functions denotes the boundary condition on and is the solution to the problem in the presence of :

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∇⋅(^σ(x)∇u(l)D(x))=0inΩ\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus2.0ptminus2.0pt\omitσ0∂u(l)D∂ν(x)=g(l)(x)on∂Ω.

With this, construct as the solutions to the following problem in the absence of inclusion:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩σ0△u(l)0(x)=0%inΩ\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus2.0ptminus2.0pt\omitσ0∂u(l)0∂ν(x)=g(l)(x)on∂Ω.

Let us define following discrepancy function:

 D(Ω):=12||uD(x)−u0(x)||2L2(∂Ω)=12∫∂Ω|uD(x)−u0(x)|2dS(x). (8)

In order to derive the topological derivative , let us create inside the domain and denote be the solutions of the following problem in the presence of :

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∇⋅(σ(x)∇u(l)Σ(x))=0inΩ,\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus2.0ptminus2.0pt\omitσ0∂u(l)Σ∂ν(x)=g(l)(x)on∂Ω.

where a piecewise constant can be defined similarly with (2). With this, calculation of the topological derivative can be carried out as follows:

###### Theorem 3.2.

The topological derivative of the discrepancy function of (8) can be written as follows:

 D(Ω∖Σ)=D(Ω)+ρ(r)dTD(z)+o(r2), (9)

where

 ρ(r)=2(σ−σ0)σ+σ0πr2 (10)

and

 dTD(z)=−Re(L∑l=1∇v(l)0(z)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z)). (11)

Here, and denotes the real-part and complex conjugate of , respectively. Adjoint state is defined as the solution to

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩σ0△v(l)0(x)=0%inΩ\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus2.0ptminus2.0pt\omitσ0∂v(l)0∂ν(x)=u(l)D(x)−u(l)0(x)on∂Ω.
###### Proof.

Let us apply equation (5) to (9) then we can compute as follows:

 D(Ω∖Σ) =12L∑l=1∫∂Ω∣∣u(l)D(x)−u(l)Σ(x)∣∣2dS(x) =12L∑l=1∫∂Ω∣∣u(l)D(x)−u(l)0(x)∣∣2dS(x) =+L∑l=1∫∂Ω(u(l)D(x)−u(l)0(x))(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯u(l)0(x)−u(l)Σ(x))dS(x)+o(r2) =D(Ω)+DΣ(z)+o(r2).

where

 DΣ(z)=L∑l=1∫∂Ω(u(l)D(x)−u(l)0(x))(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯u(l)0(x)−u(l)Σ(x))dS(x).

By applying asymptotic expansion formula (5) and equations (2.1) and (3.2), can be written:

 DΣ(z)= L∑l=1∫∂Ω(u(l)0(x)−u(l)D(x))(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯u(l)Σ(x)−u(l)0(x))dS(x) = L∑l=1∫∂Ωσ0∂v(l)0∂ν(x)(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯u(l)Σ(x)−u(l)0(x))dS(x) = πr2L∑l=1∫∂Ωσ0∂v(l)0∂ν(x)(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z)⋅M(z)⋅∇N(x,z)+o(r2))dS(x) = πr2L∑l=1∫Ωσ0△v(l)0(x)(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z)⋅M(z)⋅∇N(x,z))dx πr2L∑l=1∫Ω∇v(l)0(x)⋅(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z)⋅M(z)⋅σ0△N(x,z))dx = Extra open brace or missing close brace = −πr2L∑l=1∇v(l)0(z)⋅M(z)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z).

for . Therefore, by (6),

 DΣ(z)=−2σ−σ0σ+σ0πr2L∑l=1∇v(l)0(z)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z). (12)

Finally, by taking real part of (12), we can obtain equations (10) and (11). ∎

## 4 Structure of topological derivative

From now on, we identify the structure of topological derivative (11) and discuss certain properties.

Since satisfies (3.2), following relations hold for ,

 v(l)0(z) =∫∂Ωσ0∂v(l)0(y)∂ν(y)N(z,y)dS(y) ∇zv(l)0(z) =∫∂Ωσ0∂v(l)0(y)∂ν(y)∇zN(z,y)dS(y).

With them, applying boundary condition of (3.2) and asymptotic expansion formula (5) yields

 ∇zv(l)0(z)=∫∂Ωσ0(u(l)D(y)−u(l)0(y))∇zN(z,y)dS(y)=πr2∫∂Ω(∇u(l)0(x)⋅M(x)⋅∇xN(x,y))∇zN(z,y)dS(y). (13)

Therefore, (11) can be written as follows:

 dTD(z) =πr2Re{[∫∂Ω(∇u(l)0(x)⋅M(x)⋅∇xN(x,y))∇zN(z,y)dS(y)]⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z)} =πr2Re{∇u(l)0(x)⋅M(x)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z)(∫∂Ω∇xN(x,y)∇zN(z,y)dS(y))}.

From the fact that in the two-dimensional space, Neumann function can be decomposed with the singular and regular parts:

 Extra open brace or missing close brace (14)

where for any and solves

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩ΔR(x,y)=0% forx∈Ω\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus2.0ptminus2.0pt\omit∂R(x,y)∂ν(x)=−1|∂Ω|+1πx−y|x−y|2⋅ν(x)forx∈∂Ω.

Particularly, when the domain is a ball with radius centered at origin, Neumann function is written by (see (5))

 N(x,y)=−12πln|x−y|−12πln∣∣∣R|x|x−|x|Ry∣∣∣+lnRπ.

With this decomposition, we consider the following value

 ∫∂Ω ∇xN(x,y)∇zN(z,y)dS(y) = ∫∂Ω(−12πx−y|x−y|+∇xR(x,y))(−12πz−y|z−y|+∇zR(z,y))dS(y) = 14π2∫∂Ω(x−y|x−y|⋅z−y|z−y|)dS(y):=T1−12π∫∂Ω(x−y|x−y|⋅∇zR(z,y))dS(y):=T2 −12π∫∂Ω(z−y|z−y|⋅∇xR(x,y))dS(y):=T3+∫∂Ω∇xR(x,y)⋅∇zR(z,y)dS(y):=T4.

Note that since , there is no blowup of so that there exists a constant such that

 |∇xR(x,y)|≤C.

Hence applying Hölder’s inequality yields

 |T4| Extra open brace or missing close brace ≤C2length(∂Ω),

where means the length of . Moreover, since and , . Hence

 |T2| Missing \left or extra \right ≤Clength(∂Ω).

Since and , we must consider the singularity in and . For this purpose, we generate a ball of small radius centered at and separate into and , refer to Figure 1. Then

 |T3| =∣∣ ∣∣∫∂ΩS(z−y|z−y|⋅∇xR(x,y))dS(y)+∫∂ΩR(z−y|z−y|⋅∇xR(x,y))dS(y)∣∣ ∣∣ ≤limδ→0+(∫∂ΩS∣∣∣z−y|z−y|⋅∇xR(x,y)∣∣∣dS(y)+∫∂ΩR∣∣∣z−y|z−y|⋅∇xR(x,y)∣∣∣dS(y)) ≤limδ→0+(Cδδlength% (∂ΩS)+Clength(∂ΩR))=Clength% (∂Ω)

and

 |T1| =∣∣ ∣∣∫∂Ω(x−y|x−y|⋅z−y|z−y|)dS(y)∣∣ ∣∣≤∫∂Ω∣∣∣x−y|x−y|∣∣∣∣∣∣z−y|z−y|∣∣∣dS(y) ≤limδ→0+(∫∂ΩS∣∣∣z−y|z−y|∣∣∣dS(y)+∫∂ΩR∣∣∣z−y|z−y|∣∣∣dS(y)) ≤limδ→0+(δδlength(∂ΩS)+length(∂ΩR))=length(∂Ω).

Therefore, we can conclude that since

 ∣∣∣∫∂Ω∇xN(x,y)∇zN(z,y)dS(y)∣∣∣≤|T1|+|T2|+|T3|+|T4| ≤14π2length(∂Ω)+C2πlength(∂Ω)+C2πlength(∂Ω)+C2length(∂Ω)<∞,

there is no blowup. Therefore, the structure of normalized topological derivative (11) will be of the form:

 dTD(z)max|dTD(z)|≈L∑l=1∇u(l)0(x)⋅M(x)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z).

### 4.1 Limitation of topological derivative to EIT problem

In practice, most of EIT systems generally apply constant current sources in several direction so that the boundary condition (4) can be written with for a constant vector , ,

 g(l)(x)=σ0∂u∂ν(x)=al⋅ν(x)forx∈∂Ω. (15)

Then since background potential is

 u(l)0(x)=al⋅x−1|∂Ω|∫∂Ωal⋅ν(y)dS(y),

normalized topological derivative (11) becomes

 dTD(z)max|dTD(z)|≈∇u(l)0(x)⋅M(x)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)0(z)=al⋅al=4σD−σ0σD+σ0.

Therefore, (11) does not offers any information of . This is a reason that why topological derivative concept cannot be applied to the EIT problem.

### 4.2 Alternative boundary condition and corresponding topological derivative

This is motivated from the original idea in (15). We would like to mention (6); (8) for its application. Suppose that there is only one inhomogeneity exists in . Let us consider the following boundary conditions

 g(l)1(x) =ik(θl+iθ⊥l)⋅ν(x)exp(ik(θl+iθ⊥l)⋅x) g(l)2(x) =ik(θl−iθ⊥l)⋅ν(x)exp(ik(θl−iθ⊥l)⋅x).

Then it is easy to observe that

 u(l)1(x)=exp(ik(θl+iθ⊥l)⋅x)andu(l)2(x)=exp(ik(θl−iθ⊥l)⋅x).

are satisfy (3) when the boundary conditions are and , respectively. Here, is a positive real number, is an arbitrary vector on the two-dimensional unit circle and is orthogonal to with . In this paper, we set .

We introduce an alternative topological derivatives with respect to and as

 DA(z;k):=L∑l=1(∇v(l)1(z)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)1(z))(∇v(l)2(z)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)2(z)), (16)

where and satisfy (3.2) with respect to and , respectively. Then, if is sufficiently large enough, we can obtain

 DA(z;k)≈ L∑l=1(∇u(l)1(x)⋅M(x)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)1(z))(∇u(l)2(x)⋅M(x)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)2(z)) = ∫S1(4kσD−σ0σD+σ0exp(ik(θ+iθ⊥)⋅x)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯exp(ik(θ+iθ⊥)⋅z)) ×(4kσD−σ0σD+σ0exp(ik(θ−iθ⊥)⋅x)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯exp(ik(θ−iθ⊥)⋅z))dθ = (4kσD−σ0σD+σ0)2∫S1exp(2ikθ⋅(x−z))dθ=16k2π(σD−σ0σD+σ0)2J0(2k|x−z|).

Here, following identity used (see (17); (24)); for sufficiently large ,

 L∑l=1exp(ikθl⋅x)≈∫S1exp(ikθ⋅x)dθ=2πJ0(kx),

where is the Bessel function of order zero and of the first kind. Therefore, normalized topological derivative will be of the form

 dTDA(z;k)=DA(z;k)max|DA(z;k)|=J0(2k|x−z|).

This gives some certain properties of summarized as follows.

1. Since reaches its maximum value at , will plot its maximum value at so that we can identify the location of .

2. Resolution of image is highly depends on the value of and . Based on the property of (see Figure 2), one can obtain a result with high resolution if is sufficiently large. In contrast, if the value of is small, one cannot identify the location of , refer to Figure 3.

3. Due to the assumption of (1), if is (nearly) touching the boundary , it is very hard to conclude that yields a good image.

4. Asymptotic expansion formula (5) holds for small inhomogeneity in theory. Therefore, when is an extended target, further analysis of is required.

## 5 Numerical experiments and related discussions

In this section, a numerical result is shown for showing the effectiveness of (16). For this, we take the domain to be the unit circle and we insert one inhomogeneity in the shape of ball with the radius , the location , and the conductivity . incident directions are applied and every forward problems (3), (3), (3.2) are solved via traditional Finite Element Method (FEM) in order to avoid the inverse crime. After the generation of boundary measurement data, a noise is added as follows

 u(l)D,noise(y)=(1+ξ{rand1(−1,1)+irand2(−1,1)}u(l)D(y)),

where and are arbitrary real values between and . Throughout this paper, we take . The values and , , of (16) are evaluated by the Matlab command pdegrad included in the Partial Differential Equation Toolbox.

From Figure 3, we can observe that the location of is clearly identified when is sufficiently large enough, but on the other hand we cannot identity the location when the value of is small.

We apply (16) for finding locations of multiple inhomogeneities. For this, we add another one inhomogeneity in the shape of ball with the radius , the location , and the conductivity . Figure 4 shows the maps of with various values of . Unfortunately, in contrast with the imaging of single inhomogeneity, we cannot identify and .

In order to reveal the reason, we reconsider (16) in the existence of two-different inhomogeneities and . In this case, by a simple calculation, (16) becomes

 DA(z;k):= L∑l=1(∇v(l)1(z)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)1(z))(∇v(l)2(z)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)2(z)) = L∑l=1{2∑m=1(∇u(l)1(xm)⋅M(xm)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇u(l)1(z