# Hybrid Topological Derivative-Gradient Based Methods for Nondestructive Testing

Authors: A. Carpio; M.-L. Rapún
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper is devoted to the reconstruction of objects buried in a medium and their material properties by hybrid topological derivative-gradient based methods. After illustrating the techniques in time-harmonic acoustic problems with different boundary conditions and in electrical impedance tomography problems with continuous Neumann conditions, we extend the hybrid method for a realistic model in tomography where the boundary conditions are
given at a discrete set of electrodes.


---

## Body

## 1. Introduction

In many practical aspects of life we face the need of reconstructing the geometry of objects buried in a medium and identifying their material properties. In medicine, precise detection systems are required to locate tumors, clots, or damaged tissues without harming the patient. Magnetic resonance imaging together with different kinds of tomography, X-rays, or ecography have become powerful diagnosis tools. Geophysicists have long relied on measuring scattered waves to obtain information on the structure of the ground: from locating water, oil, and gas reservoirs or mineral mines to identifying the nature of earth quake sources or fault structure. Assessing the structural integrity of bridges, planes, power plants, and buildings after natural disasters also requires sharp techniques for nondestructive testing to identify cracks or damaged areas.A classical strategy to obtain information about the inner structure of a medium is to illuminate it with some kind of radiation; see Figure1. The total wave is usually measured at a set of receptors distributed over a curve or surface Γmeas. From these measurements, we aim to reconstruct the inner geometry of the medium and their material coefficients. Objects are distinguished from the background because they differ in their material parameters. The nature of the incident radiation is chosen to enhance that contrast. Acoustic waves are used when the parameters of interest are the elastic constants. Electromagnetic scattering tracks the electric conductivity and permittivity. Thermal waves allow to distinguish objects by their thermal diffusivity and conductivity.Figure 1
Scheme of the inverse problem.There are two important mathematical problems involved in this process: the forward problem and the inverse problem. In the forward problem, the material properties of the background medium and the shape, size, location, and material parameters of all the inclusions are known. The goal is to compute the wave field at the receptors. It is usually a well-posed boundary value problem with a unique solution that depends continuously on the initial data. The equations governing the wave field will be Navier, Maxwell, or heat equations depending on the nature of the chosen incident radiation. In the inverse problem, measurementsumeas of the wave field at the receptors are known. The objective is to identify the geometry and nature of the objects located inside the background medium. Those objects are characterized with the fact that solving the forward problem using the known incident wave and placing those objects inside the medium, the resulting wave field evaluated at the receptors agrees with the measured data. This inverse scattering problem is nonlinear, and severely ill posed. It may not have a solution, the solution may not be unique, or it may not depend continuously on the data.The inverse problem can be reformulated in a more regular way as a constrained optimization problem. In practice, one does not need to know the exact geometry and material parameters of the objects for which the wave field evaluated at the receptors agrees exactly with the measured data. It is enough to find approximations making the error small enough. This motivates a weaker variational reformulation of the inverse problem: Find domainsΩ, and parameters ki minimizing
(1)J(ℛ∖Ω¯,ki)=12∫Γmeas‍|u-umeas|2dl,
where u is the solution of the forward problem with objects Ω immersed in the exterior media ℛ∖Ω¯ and parameters ki (i stands for “interior”). The design variables are Ω and ki. The forward problem is the constraint.Developing mathematical tools to solve this type of inverse problems is an active field, in need of even better techniques. Many beautiful mathematical theories have been proposed to produce reasonable reconstructions in a large variety of contexts from different standpoints, but the demand of more refined methods remains. Some of the first theories relied on linear approximations, such as the Kirchhoff or physical optics approximation [1] or the Born approximation [2]. Linear sampling methods became a step forward [3, 4]. These methods formulate an equation for a parametrization of the contour of the unknown domains in terms of its representative coefficients. Backpropagation principles [5] and modified gradient methods [6] have also been employed. More recently, different optimization strategies have been proposed to deform initial guesses of the contours of the objects in such a way that adequate cost functionals decrease. The techniques differ by the way the objects are parameterized and deformed. Classical shape deformation along vector fields was tried first [7–10]. This strategy has a drawback: it does not allow for topological changes in the objects. The number of boundaries (connected components, holes, etc.) has to be known from the beginning, and convergence of the method is unlikely without a priori information. Representing the objects as level sets of level set functions deformed following Hamilton-Jacobi equations, with velocities chosen to ensure the decrease of the cost functional, solves that drawback [11, 12]. Topological changes are allowed during the deformation process. However, level set methods need an external guess to start, and convergence may be slow unless the guess is good enough. Topological derivative techniques have arisen as a promising alternative. They have the advantage of providing good first guesses for the number, size, and location of the objects, serving also as a basis for iterative schemes that improve those initial approximations. Their potential to provide good first approximations was explored in [13] for exterior acoustic problems involving rigid bodies, in [14, 15] for three-dimensional transmission problems, and in [16] for waveguides with Dirichlet boundary conditions. Reference [17] discussed extensions to scattering by elastic waves. Reference [18] and references therein suggested possible applications in electromagnetism and tomography. A hybrid level set-topological derivative scheme was proposed in [19]. Iterative schemes entirely based on topological derivatives were introduced in [20] and later developed in [21–23]. Modified cost functionals allow to prove theoretical convergence results and study the effect of noise [24, 25].A vast amount of work in the field refers to time-harmonic incident waves, for which the boundary value problem governing the wave field becomes stationary. Less variety of methods are available when the incident waves are not time harmonic, for instance, for pulse-like excitations. For incident waves governed by time-dependent wave equations, see the classical work [26] and references therein. Methods for detection of inclusions near the surface by photothermal techniques are developed in [27].The previuosly mentioned references focus mostly on reconstructing the distribution of the scatterers. In practice, one also needs to predict their material parameters to know their nature. In medical applications, for instance, this is essential to be able to distinguish benign from malignant tumors. There is a vast amount of the literature about reconstructing the spatial variations of parameters dealing in particular with electrical impedance tomography; see the reviews [28, 29] and references therein. The complex case and the problem with discrete electrodes are much less studied. Schemes following discontinuities or objects with sharp interfaces usually consider a piecewise-constant conductivity [30–35]. Some of them seek for the objects without being able to precise the value of the conductivity [30, 36]. Others may require a large number of iterations [31–33, 35]. Level set techniques have been recently extended to track interfaces between regions of different conductivity [31, 35]. These methods are initialized resorting to topological derivatives in [32, 33]. For smooth conductivities, [37–39] describe promising methods.In elastic scattering, spatial dependence is unavoidable if the background matrix is heterogeneous or as a way to detect damage; see [40] for tests in functionally graded reference materials. Hybrid topological derivative-gradient based methods to reconstruct both the geometry of objects and their material parameters were introduced in [21, 22] for acoustic problems and in [23] for the electrical impedance tomography problem with continuous Neumann conditions.In this paper, we explore the potential of some specific optimization techniques to solve imaging problems. We focus on hybrid descent techniques that combine gradient methods to approximate the parameters and topological derivative based schemes to calculate the domains. We review some recent results and present new developments for tomography. The paper is organized as follows. Section2 exemplifies the hybrid schemes for time-harmonic incident waves. We focus on acoustic waves, though some types of polarized electromagnetic waves lead to inverse problems with exactly the same mathematical structure. Electrical impedance tomography is considered in Section 3, both the continuous and the discrete versions of the electrode model. We reformulate the problem with discrete electrodes as a constrained optimization problem. Expressions for the topological derivative with respect to the inclusions and the correctors of their parameters are given in terms of adjoint fields satisfying adequate variational equations.


## 2. Methods for Inverse Acoustic Problems

In this section, we deal with inverse acoustic scattering problems. We assume that the incident wave is time harmonic; that is,(2)Uinc(x,t)=Re[e-iωtuinc(x)],
where ω>0 is the frequency. In this case, after a sufficiently long time, the total wave field is also time harmonic, U(x,t)=Re[e-iωtu(x)], and the total amplitude u(x) is governed by a boundary value problem for the Helmholtz equation.Depending on the nature of the objects, different boundary conditions model the interaction between the scatterers, the medium, and the incident wave. For penetrable objects transmission conditions are imposed at the interfaces between the objects and the surrounding medium, whereas Dirichlet and Neumann conditions model sound-soft and sound-hard objects, respectively.We assume that the incident wave is a planar wave propagating in the directiond, uinc(x)=eikex·d with wave number ke>0. Other types of incident waves can be considered, for instance, in ℝ2, uinc(x)=(ı/4)H0(1)(ke|x-x0|) models a point source excitation at x0. Here, H0(1) is the Hankel function of the first kind and order zero (see [41]).Let us begin with describing the geometrical configuration of our inverse acoustic scattering problems (see Figure1). We have a finite number of bounded and simply connected open objects Ω1,…,Ωd having no pairwise contact. Their boundaries are assumed to be smooth, for instance 𝒞2. For simplicity, we will use the notation Ω=∪j=1dΩj. These objects are immersed in the exterior media ℝ2∖Ω¯. We focus on two-dimensional problems to illustrate the methods, though the theory extends to any dimension.Let us start by considering penetrable nonhomogeneous obstacles. To simplify, we assume that the mediumℝ2∖Ω¯ is homogeneous. Then, the total wave field
(3)u=uinc+uscinℝ2∖Ω¯,u=utrinΩ,
where usc and utr denote the scattered wave outside the obstacles ℝ2∖Ω¯ and the transmitted wave utr inside Ω, respectively, is the solution of
(4)Δu+ke2u=0inℝ2∖Ω¯Δu+ki2u=0inΩu-=u+,∂nu-=∂nu+on∂Ωlimr→∞r1/2(∂r(u-uinc)-ike(u-uinc))=0.
Here, ke is the wave number of the incident wave and ki(x)>κi>0 is related to the material properties of the inclusions Ω. We assume that ki is differentiable. When ki is constant inside each Ωj, we are dealing with homogeneous materials.The condition at infinity is the standard Sommerfeld radiation condition withr=|x|, which has to be satisfied uniformly in all directions. It implies that only outgoing waves are allowed. ∂r denotes radial derivatives. ∂n stands for outward normal derivatives, with n pointing outside Ω. u+ and u- denote the limits of u from the exterior and interior of Ω, respectively. The continuity condition for the normal derivative can be extended to αe∂nu-=αi∂nu+ when Δu is replaced by div(αe∇u) in ℝ2∖Ω¯ and by div(αi∇u) in Ω, as explored in [21].The optimization problem we want to solve is then as follows: FindΩ and ki(x) minimizing
(5)J(ℝ2∖Ω¯,ki)=12∫Γmeas‍|u(x)-umeas(x)|2dl,
where u(x) is the solution of the forward problem (4) with objects Ω and parameters ki. This cost functional is positive but vanishes for the true scatterers and the true parameters, which are characterized as a global minimum of the functional.A common strategy to approximate solutions of optimization problems is to resort to descent techniques. The idea is to find initial guesses for the parameterski,0 and the objects Ω0 and then generate sequences of approximations ki,m, Ωm along which the functional J decreases. In the next subsections, we explain how to generate those sequences in several stages. First, we assume that the domains Ω are known and present a gradient strategy to guess ki. Next, we will assume that the parameters ki are known and discuss how to approximate Ω by topological derivative techniques. Then, we will combine both methods performing a double iteration with respect to the domains and the parameters to approximate both the objects and their material constants. Once the method is presented for penetrable objects, we will discuss how it is modified to handle rigid or sound-soft obstacles.


### 2.1. Gradient Techniques to Identify Coefficients

Let us assume that the objectsΩ are known. The optimization problem becomes as follows: Find ki(x) minimizing
(6)J(ki)=12∫Γmeas‍|u(x)-umeas(x)|2dl,
where u(x) is the solution of the forward problem (4) with objects Ω and parameters ki.Fixed an initial guesski,0, we generate a sequence ki,m providing a decreasing sequence J(ki,m) by a gradient method [21]. At each stage, we fix the current guess ki,m and study the variations of the functional along directions v: J(η)=J(ki,m+ηv), where η>0 is a real parameter. The derivative with respect to η is given by
(7)dJ(η)dη=dJ(ki,m+ηv)dη=Re[∫Γmeas‍(uη(x)-umeas(x)¯)duη(x)dηdl],uη being the solution of the forward problem (4) with coefficient ki=ki,m+ηv inside the objects Ω. Computing duη/dη is avoided introducing an adjoint field p; see [21].Theorem 1.
The derivative with respect toη of the function J(η)=J(ki,m+ηv) defined by (6) evaluated at zero is given by
(8)dJ(η)dη|η=0=2Re[∫Ω‍vup],
where p is the solution of the adjoint problem
(9)Δp+ke2p=(umeas-u)¯δΓmeasinℝ2∖Ω¯Δp+ki2p=0inΩp-=p+,∂np-=∂np+on∂Ωlimr→∞r1/2(∂rp-ikep)=0
and u the solution of the forward problem (4) setting in both cases ki=ki,m.Choosingv in such a way that the derivative J′(0)<0, we ensure that J(ki,m+ηv)<J(ki,m) for η>0 small enough. For instance, we may set
(10)v(x)=-Re[u(x)p(x)].
The updated guess is then
(11)ki,m+1(x)=ki,m(x)-ηRe[u(x)p(x)].
This choice is reasonable if we are looking for space-dependent parameters. When Ω is composed of several disjoint domains Ωj, and ki is known to be almost piecewise constant, we may set
(12)v(x)=vj=-Re[∫Ωj‍up],x∈Ωj.
Technical details and proofs can be found in [21], where more general Helmholtz type constraints are considered:
(13)div(αe(x)∇u)+ke(x)2u=0inℝ2∖Ω¯,div(αi(x)∇u)+ki(x)2u=0inΩ,u-=u+,αe(x)∂nu-=αi(x)∂nu+on∂Ω,limr→∞r1/2(∂r(u-uinc)-iκe(u-uinc))=0.
The coefficients ke, ki, αe, αi are bounded from below by positive constants and depend smoothly on x, becoming constant at infinity. The constant κe at the Sommerfeld radiation condition is defined as κe=lim|x|→∞ke(x)/αe(x). It is the wave number of the incident wave. The interior parameter αi(x) can also be identified following the same technique.


### 2.2. Topological Derivative Techniques to Reconstruct Domains

We have seen how to generate descent directions for the parameterski fixing the objects Ω. Let us assume now that we know the parameters ki. The optimization problem becomes as follows: Find Ω minimizing
(14)J(ℝ2∖Ω¯)=12∫Γmeas‍|u(x)-umeas(x)|2dl,
where u(x) is the solution of the forward problem (4) with objects Ω and parameters ki. To generate descent directions for the domains, we must rely on some kind of derivative with respect to the domains. We choose a specific type of derivative for shape functionals: the topological derivative. It has the advantage of providing a first guess for the objects too.The topological derivative of a shape functionalJ(ℛ) at the point x∈ℛ is defined as [42]
(15)DT(x,ℛ)=limɛ→0J(ℛ∖Bɛ(x)¯)-J(ℛ)h(ɛ),
where Bɛ(x) is a ball centered at x with small radius ɛ>0; see Figure 2. The scalar function h(ɛ)>0 is chosen in such a way that the limit (15) is finite and does not vanish. For instance, for the transmission or Neumann conditions in two dimensions, we take h(ɛ)=Vol(Bɛ(x))=πɛ2. For Dirichlet conditions in 2D, we may choose h(ɛ)=-2π/log(keɛ).Figure 2
Geometry for the computation of topological derivatives.The topological derivative is a scalar function defined for everyx∈ℛ which measures the sensitivity of the functional to removing from ℛ a small object located at x. In regions where DT(x,ℛ) takes large negative values, we have a large probability to find an object. This remark is based on the following expansion:
(16)J(ℛ∖Bɛ(x)¯)=J(ℛ)+h(ɛ)DT(x,ℛ)+o(h(ɛ)),J(ℛ∖Bɛ(x)¯)=J(ℛ)h(ɛ)51..DT(x,ℛ)ɛ⟶0.
Whenever DT(x,ℛ)<0, J(ℛ∖Bɛ(x)¯)<J(ℛ) for ɛ>0 small. When J is our functional (14), this suggests a strategy to find first guesses for the objects and then improve them iteratively. Our algorithm to recover objects when the parameters are known is as follows:(i)
Compute the topological derivative whenΩ=∅ and choose an initial guess Ω0 for the objects as
(17)Ω0:={x∈ℛ∣DT(x,ℛ)≤-C0},
where C0 is a positive constant.(ii)
Check thatJ(ℛ∖Ω¯0)<J(ℛ). Otherwise increase the constant C0.(iii)
Form=1:mmax
(a)
Compute the topological derivative whenΩ=Ωm-1 and select
(18)Ωm=Ωm-1∪{x∈ℛ∖Ω¯m-1∣DT(x,ℛ∖Ω¯m-1)≤-Cm}
for decreasing thresholds Cm>0.(b)
Check thatJ(ℛ∖Ω¯m)<J(ℛ∖Ω¯m-1). Otherwise increase the constant Cm.(c)
Check the stopping criteria: the iteration stops if either meas(Ωm∖Ωm-1) is small enough, or J(ℛ∖Ω¯m) is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle |umeasɛ-um|<τɛ is satisfied, where ɛ describes measurement errors and τ>1 (in our numerical experiments τ=1.2). Here, um is the solution of the forward problem when the objects are Ωm.Computing topological derivatives at grids of points at each stage using (15) is computationally very expensive. Explicit formulas that require only the knowledge of forward and adjoint fields are used instead. For our particular Helmholtz type constraint (4), the explicit expression is [14, 21].Theorem 2.
The topological derivative ofJ(ℝ2∖Ω¯) defined by (14) is given by
(19)DT(x,ℝ2∖Ω¯)=Re[(ki2(x)-ke2)u(x)p(x)],DT(xℝ2∖Ω¯)=e[(ki2(x).-ke2)],x∈ℝ2∖Ω¯,
where u and p are solutions of the forward and adjoint problems (4) and (9), respectively.Notice that whenΩ and ki are the true objects and parameters for which the data are measured, the adjoint field p vanishes since the forward field agrees with the measured data at the receptors. Therefore, the topological derivative vanishes at the exact solution, which is a global minimum of the cost functional, as explained before.In real reconstruction tests, one usually has access to measurements for several incoming directionsdj, j=1,…,M. The shape cost functional (14) is then replaced with:
(20)J(ℝ2∖Ω¯)=12∑j=1M‍∫Γmeas|u(x;dj)-umeas(x;dj)|2dlx.
Here, u(x;dj) is the solution to the forward problem (4) with incident field uinc(x;dj)=eıkex·dj and umeas(x;dj) is the measured total field on Γmeas. The topological derivatives for this functional are obtained adding up the topological derivatives for each single incident wave.We test the algorithm for a geometrical configuration with two objects. The interior parameter in them iski=1 whereas ke=4. We have computed the topological derivative for data given at 14 sampling points uniformly distributed on the circle centered at (0,0) with radio 3, for 10 incident planar waves with uniformly distributed angles in [0,2π). The geometrical setting for our experiments is represented in Figure 3. Our synthetic data were generated by solving the corresponding forward problems by boundary element techniques using a finer mesh and adding a one percent relative random noise at each observation point.Figure 3
Geometrical description of the numerical experiments.The values of the topological derivative whenΩ=∅ at the sampling region [-1,1]×[-1,1] are represented in Figure 4(a). The largest negative values are attained inside each object. The first computation of the topological derivative provides the correct number of objects and their location, although they look like balls. Our initial guess Ω0 is represented in Figure 4(b). We compute now the topological derivative when Ω=Ω0, obtaining the color map in Figure 4(b). New regions close to the current objects appear where the topological derivative attains large negative values. Those regions should belong to the scatterers. Our objects are bigger and they are not balls. Figures 4(c)–4(j) represent the subsequent iterations. The values of the cost functional versus the number of iterations are given in Figure 4(k). After few iterations, our reconstructions are satisfactory.Reconstruction of objects whenki=1 is known using the monotone algorithm. (a) Topological derivative when Ω=∅. (b)–(j) Approximate domains Ωm superimposed to the topological derivative for Ω=Ωm, for m=0,…,8. (k) Cost functional versus the number of iterations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)The quality of the reconstructions as well as the number of iterations depends on the values of the constantsCm appearing in (18). Guidelines for their selection are given in [21]. When Cm are too big, the number of iterations required is high. A sequence of smaller Cm provides faster reconstructions, but we are prone to include spurious regions in Ωm that cannot be removed in the next iterations. We illustrate this fact in Figure 5. Figure 5(a) is exactly the same as Figure 4(a) and represents the values of the topological derivative when Ω=∅. Choosing a smaller constant C0, our initial guess Ω0, represented in Figure 5(b), has two objects that are bigger than the ones in Figure 4(b). The object on the right includes too many spurious points. Computing the topological derivative when Ω=Ω0, we detect a small region at the top of the left object that is included in the next iteration (see Figure 5(b)). No improvement for the rightmost object is pointed out. The reconstruction is given in Figure 5(c). In the next iteration (see Figure 5(d)), only the object on the left is slightly modified and the algorithm stopped because the difference between the consecutive iterations in the value of the functional was negligible, as observed in Figure 5(e). The quality of the reconstruction of the object on the left is reasonable, but, comparing with the previous example, the reconstruction of the object on the right is not acceptable.Reconstruction of objects whenki=1 is known, with the monotone algorithm (18). (a) Topological derivative when Ω=∅. (b)–(d) Approximate domains Ωm superimposed to the topological derivative for Ω=Ωm. (e) Evolution of the cost functional.
(a)
(b)
(c)
(d)
(e)The monotone algorithm cannot remove spurious regions in the rightmost object. The topological derivative is only defined at points not belonging to the current approximation. However, observing the values of the topological derivative in Figures5(b)–5(d), we realize that there are regions close to object where it takes large positive values. Computing the topological derivative when Ω=∅ (see Figure 5(a)), we did not see that effect. This indicates that spurious regions are included and that we should use a nonmonotone method if we want to correct them.The classical notion of topological derivative we used previously is defined only inℛ∖Ω¯. Therefore, any time we update our approximations Ωm, we have a criterion to add points. However, if at some stage we include in our approximation points that do not belong to the objects, we cannot remove them. This artifact is eliminated observing that for our particular type of functional we can extend the definition of topological derivative to points inside the objects. The topological derivative of a functional J(ℛ∖Ω¯) was extended to the points inside Ω in [21] as
(21)DT(x,ℛ∖Ω¯)=limɛ→0J(ℛ∖(Ω∖Bɛ(x)¯))-J(ℛ∖Ω¯)h(ɛ),DT(x,ℛ∖Ω¯)=limɛ→0J(ℛ∖(Ω∖Bɛ(x)¯))21.-Jh(ɛ)x∈Ω.
From this definition, the expansion below follows:
(22)J(ℛ∖(Ω∖Bɛ(x)¯))=J(ℛ∖Ω¯)+h(ɛ)DT(x,ℛ∖Ω¯)+o(h(ɛ)),=J(ℛ∖Ω¯)+h(ɛ)DT(x,ℛ∖1Ω¯).ɛ⟶0.
Since h(ɛ) is a positive function, if DT(x,ℛ∖Ω¯) is large and negative, the functional decreases when we remove the ball Bɛ(x) from Ω. This suggests a strategy to remove points from Ω. Our algorithm can generate a nonmonotone sequence of objects by substituting the definition of Ωm in the previous method (see (18)) with
(23)Ωm=Ωm-1∪{x∈ℛ∖Ω¯m∣DT(x,ℛ∖Ω¯m-1)≤-Cm}∖{x∈Ωm-1∣DT(x,ℛ∖Ω¯m-1)≤-cm}¯,
with decreasing sequences Cm, cm>0.Using the extended definition (21), we find the following expression for the topological derivative of the functional (14) of our transmission problem (see [21]).Theorem 3.
The extended topological derivative (21) of J(ℝ2∖Ω¯) defined by (14) is given by
(24)DT(x,ℝ2∖Ω¯)=Re[(ke2-ki2(x))u(x)p(x)],DT(Ω¯)=Re[(ke2-ki2(x))u(x)1p(x)]x∈Ω,
where u and p are solutions of the forward and adjoint problems (4) and (9), respectively.Notice that the expressions (19) and (24) are identical except for a minus sign. Moreover, since the solutions to the problems (4) and (9) are continuous functions in ℝ2, the expression
(25)Re[(ki2(x)-ke2)u(x)p(x)]
defines a continuous function in ℝ2. This motivated in [21] to modify the extended definition of the topological derivative inside Ω of a functional J(ℛ∖Ω¯) as
(26)DT(x,ℛ∖Ω¯)=limɛ→0J(ℛ∖Ω¯)-J(ℛ∖(Ω∖Bɛ(x)¯))h(ɛ),DT(x,ℛ∖Ω¯).=limɛ→0J(ℛ∖Ω¯)-J(ℛ∖(..5))h(ɛ),x∈Ω,
which is exactly the same as (21) multiplied by -1. Using this notion,
(27)J(ℛ∖(Ω∖Bɛ(x)¯))=J(ℛ∖Ω¯)-h(ɛ)DT(x,ℛ∖Ω¯)+o(h(ɛ)),=J(ℛ∖Ω¯)-5.h(ɛ)DT(x,ℛ∖Ω¯).ɛ⟶0,
and whenever DT(x,ℛ∖Ω¯) is large and positive, the functional decreases when removing the ball Bɛ(x) from Ω. This definition also allows to define a non-monotone scheme by substituting the definition of Ωm in (23) by
(28)Ωm=Ωm-1∪{x∈ℛ∖Ω¯m∣DT(x,ℛ∖Ω¯m-1)≤-Cm}∖{x∈Ωm-1∣DT(x,ℛ∖Ω¯m-1)≥cm}¯,
for decreasing thresholds Cm, cm>0.With this definition, the topological derivative of the functional for the transmission problem is as follows.Theorem 4.
The extended topological derivative (26) of J(ℝ2∖Ω¯) defined by (14) is given by
(29)DT(x,ℝ2∖Ω¯)=Re[(ki2(x)-ke2)u(x)p(x)],x∈Ω,
where u and p are solutions of the forward and adjoint problems (4) and (9), respectively.There is a missprint in [21]. Definition (21) was given, but in fact (26) was used. The topological derivative inside Ω was computed as in (29) and in the numerical examples the implemented formula is (29). The modified algorithm described by (28) was used to update the domains.Let us revisit now the example of Figure5. The topological derivative for Ω=∅ provided the initial guess Ω0 represented in Figure 6(b) by a dashed line. The values of the topological derivative inside Ω in Figure 6(b) are calculated using the definition (21), that is, formula (24). The object on the right includes spurious regions (the points with large negative values inside Ω0) that will be removed in the next iterations. Some points at the top of the object on the left should be included in that object. If we use the alternative definition (26), the topological derivative is a continuous function (see Figure 6(c)). A region inside the rightmost object where large positive values are attained should be removed. A region outside the leftmost object where it takes large negative values has to be included. The second definition provides a color map that is easier to interpret. Therefore, from now on we choose the definition (26) and update the domains using (28). The final reconstruction at the 11th iteration is represented in Figure 6(d) by a dashed red line. For comparison, we have included in Figure 6(e) the final reconstruction obtained for the same problem in Figure 4 with the monotone method (for a well-calibrated sequence of constants Cm) after nine iterations. Our non-monotone method is not very sensitive to the initial guess.Comparison of the alternative definitions of the topological derivative inside an object and of the monotone and nonmonotone methods. (a) Topological derivative whenΩ=∅. (b) Initial guess Ω0 (dashed lines) superimposed to the topological derivative for Ω=Ω0 using the extended definition (21) inside Ω. (c) Counterpart of (b) using the extended definition (26). (d) Final reconstruction at the 11th iteration using the non-monotone method (dashed red line). (e) Final reconstruction at the 9th iteration using the monotone method (dashed red line). (f) Evolution of the cost functional applying the non-monotone strategy.
(a)
(b)
(c)
(d)
(e)
(f)A wide gallery of reconstructions using the monotone and the non-monotone methods can be found in [21, 43], illustrating the ability of the method to reconstruct objects with different parameters ke, ki (simulating high and low frequencies), problems with different values of ki inside each object, scatterers of different sizes, different levels of noise in the data, reconstructions when data are given on a limited region of a circle, source point excitations instead of planar incident waves, and so forth.


### 2.3. Hybrid Approach to Reconstruct Objects and Identify Their Material Parameters

Combining the two strategies to approximate objects and parameters discussed in the previous two sections, we obtain a reconstruction scheme [21, 22]. (i)
Initialization
(a)
Choose an initial guesski,0 for the parameters ki, for example, a perturbation of the background parameter ke.(b)
Choose an initial guessΩ0 for the objects Ω, using the topological derivative: set ki=ki,0 and compute the topological derivative with object Ω=∅:
(30)DT(x,ℝ2)=Re[(ki,02-ke2)u(x)p(x)],whereu and p are forward and adjoint fields with object Ω=∅ and parameter ki,0. Then,
(31)Ω0={x∈ℝ2∣DT(x,ℝ2,ki,0)≤-C0},whereC0 is a certain large positive threshold.(ii)
Iteration
(a)
Updateki,m(x) using the gradient technique:
(32)ki,m+1(x)=ki,m(x)-ηRe[u(x)p(x)],whereu, p solve forward and adjoint problems with objects Ωm and parameter ki,m(x).(b)
Check thatJ(ℝ2∖Ω¯m,ki,m+1)<J(ℝ2∖Ω¯m,ki,m). Otherwise, reduce η.(c)
UpdateΩm computing the topological derivative with objects Ωm and parameter ki,m+1(x). We may resort to a monotone
(33)Ωm+1=Ωm∪{x∈ℝ2∖Ω¯m∣DT(x,ℝ2∖Ω¯m,ki,m+1)≤-Cm+1},or nonmonotone strategy(34)Ωm+1=Ωm∪{x∈ℝ2∖Ω¯m∣DT(x,ℝ2∖Ω¯m,ki,m+1)=Ωm∪55x∈ℝ2∖Ω¯m∣DT(x,ℝ2∖Ω¯m,ki,m+1)≤-Cm+1}∖{x∈Ωm∣DT(x,ℝ2∖Ω¯m,ki,m+1)≥cm+1}¯,with decreasing thresholdsCm+1, cm+1. The second one is more convenient if we wish to be able to remove spurious points at any stage, for instance, when reconstructing objects with holes.(d)
Check thatJ(ℝ2∖Ω¯m+1,ki,m+1)<J(ℝ2∖Ω¯m,ki,m+1). Otherwise, increase the thresholds cm+1, Cm+1.(e)
Stop when meas(Ωm∖Ωm-1) is small enough, or J(ℝ2∖Ω¯m) is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle |umeasɛ-u|<τɛ is satisfied.The iterative procedure described perviously generates sequences of domainsΩm and parameters ki,m along which the functional decreases. The magnitude of the topological derivative is also observed to decrease. Theoretically, convergence to the true solution is not guaranteed, since our cost functional might have local minima different from the global minimum we are looking for. This situation can be improved or worsened by the amount of incident waves we use, the number of receptors, and how we distribute them. In the tests we have performed, the thresholds were calibrated so that the scheme worked quite well in a few iterations. See [21, 22] for details about calibration, numerical schemes and extensions to constraints of the form (13).The forward and adjoint problems have explicit solutions whenke is constant and we are computing the first guess Ω0 on the whole space. They become
(35)Δu+ke2u=0inℝ2limr→∞r1/2(∂r(u-uinc)-ike(u-uinc))=0,(36)Δp+ke2p=(umeas-u)¯δΓmeasinℝ2limr→∞r1/2(∂rp-ikep)=0.
The solution to (35) is the incident wave and the solution to (36) can be computed using the fundamental solution of the Helmholtz equation. The forward and adjoint problems (4) and (9) are solved by either boundary element or finite element methods depending on whether ki is piecewise constant or vary smoothly with space [11, 13, 21, 44–46].In practice, updating the domains is more expensive than updating the parameters. Any time we change the guess of the objects, we must mesh again the new domains to be able to compute forward and adjoint fields. Therefore, we usually correct the parameterski,m several times for a fixed approximate domain Ωm before proceeding to update Ωm, to reduce the computational cost.We have tested the hybrid topological derivative-gradient based method considering the configuration with two objects of the previous examples, keeping the same observation points and incident waves (see Figure3). We also take the same values of the parameters: ke=4 is assumed to be known and ki=1 is assumed to be unknown. We computed the topological derivative when Ω=∅ choosing ki,0=3 as initial guess for ki. Notice that since the forward and adjoint fields when Ω=∅ do not depend on ki, the topological derivative is the same as in Figure 4(a) except for the scale (we have now in (19) the multiplying factor ((ki,0)2-ke2) instead of (ki2-ke2)). Therefore, our initial guess Ω0, represented in Figure 7(a), is the same as in Figure 4(b). We fix now Ω=Ω0 and iterate five times to update the value of ki without modifying the scatterers. Then, we fix the value of ki and compute the topological derivative to update the domains and so on. The reconstruction at the ninth iteration with respect to the domain is shown in Figure 7(b). The values of ki through the iterative procedure are represented in Figure 7(c). Two identical values of the parameter in the plot mark each iteration to improve the domains (ki is not updated). Our final approximation is 1.23 while the true value is 1. Some reconstructions with two or three objects with different constant parameters can be found in [21, 43].Reconstruction of objects with unknownki. (a) First guess of the objects. (b) Final reconstruction of the objects after nine updates of the domain (c) Values of ki through the iterative method.
(a)
(b)
(c)In the next example we test the hybrid algorithm for the reconstruction of an objectΩ with space-dependent coefficient ki. The scatterer is the ball centered at (-0.25,0.3) and radius 0.35 and ki is the function represented in Figure 8(c). It grows radially from 1 to 2 inside Ω. The observation points and the incident waves are the same as in all the previous examples. To start the algorithm, we took as initial guess the constant function ki,0=3. The topological derivative for Ω=∅ provided the initial guess Ω0 represented by a red dashed line in Figure 8(a). We applied the hybrid method alternating one iteration with respect to the domain with five steps of the gradient method. The final reconstructions of Ω and ki are shown in Figures 8(b) and 8(d), respectively. The parameter ki is overestimated with an average error about 0.2, but we clearly recognize a function that grows radially. The location, shape, and size of the object are reconstructed with accuracy. The interested reader can find in [22] some numerical experiments including the reconstruction of two homogeneous objects immersed in an heterogeneous material and of two heterogeneous scatterers buried in an homogeneous media.Reconstruction an object with space-dependent parameterki. (a) True object (blue solid line) and initial guess (red dashed line). (b) True object (blue solid line) and its final reconstruction (red dashed line) at the ninth iteration with respect to the domain. (c) True function ki. (d) Final reconstruction of ki.
(a)
(b)
(c)
(d)


### 2.4. Sound-Hard and Sound-Soft Objects

The methods presented in the previous sections apply when the objects are penetrable by the incident radiation. Part of the incident wave is transmitted inside the object and transmission boundary conditions are imposed at the interface between the objects and the surrounding medium.When the scatterers are sound-hard (rigid), no transmitted wave is generated. The incident radiation is scattered and a Neumann boundary condition is imposed at the surface of the objects. The constraint for the cost functional (14) becomes
(37)Δu+ke2u=0,inℝ2∖Ω¯∂nu=0,on∂Ωlimr→∞r1/2(∂r(u-uinc)-ike(u-uinc))=0.
Topological derivative methods to find first guesses of the number, shape, and location of rigid objects were developed in [13, 44]. However, to the best of the authors’ knowledge, iterative methods entirely based on the computation of topological derivatives for the reconstruction of sound-hard objects had not been previously used.For sound-soft scatterers, a Dirichlet boundary condition is imposed at the interface between the objects and the background medium:(38)Δu+ke2u=0,inℝ2∖Ω¯u=0,on∂Ωlimr→∞r1/2(∂r(u-uinc)-ike(u-uinc))=0.
Topological derivative methods were first employed to reconstruct sound-soft objects in [16, 20, 47].For rigid and sound-soft objects, we cannot predict their material parameters, since they do not appear in the constraints. The cost functionals only depend on the geometry of the inclusionsΩ. The topological derivative based iterative methods sketched in the previous sections can be extended to these two types of objects, replacing the explicit expressions of the topological derivatives with other adequate formulas, which involve modified forward and adjoint fields. The forward and adjoint problems are formally similar to the transmission ones, suppressing the equation for the transmitted wave inside Ω and imposing either Neumann or Dirichlet boundary conditions on ∂Ω.The new expressions for the topological derivatives are:(i)
For the Neumann problem as follows(39)DT(x,ℝ2∖Ω¯)=Re[2∇u(x)∇p(x)-ke2u(x)p(x)],
where the forward field is a solution of (37) and the adjoint field is a solution of
(40)Δp+ke2p=(umeas-u)¯δΓmeas,inℝ2∖Ω¯,∂np=0,on∂Ω,limr→∞r1/2(∂rp-ikep)=0.
This formula was obtained in [21], correcting the previous paper [13] where the multiplying factor 2 in (39) was missed.(ii)
For the Dirichlet problem (see [47])
(41)DT(x,ℝ2∖Ω¯)=Re[u(x)p(x)],
where the forward field is a solution of (38) and the adjoint field is a solution of
(42)Δp+ke2p=(umeas-u)¯δΓmeas,inℝ2∖Ω¯,p=0,on∂Ω,limr→∞r1/2(∂rp-ikep)=0.We test the performance of our iterative topological derivative method for rigid and sound-soft objects in Figures9 and 10, respectively. The geometrical setting, sketched in Figure 3, is the same as in the previous examples for the transmission problem with constant coefficients. We also keep the same observation points and incident waves. The (exterior) wave number is ke=4 too.Reconstruction of rigid objects. (a) Topological derivative whenΩ=∅. (b) First guess of the objects (dashed red line). (c) Final reconstruction of the objects after 9 steps (dashed red line).
(a)
(b)
(c)Reconstruction of sound-soft objects. (a) Topological derivative whenΩ=∅. (b) Initial guess Ω0 superimposed to the topological derivative for Ω=Ω0. (c)-(d) Approximate domains Ωm superimposed to the topological derivative for Ω=Ωm, for m=1 and m=2. (e) Final reconstruction of the objects after 13 steps (dashed red line).
(a)
(b)
(c)
(d)
(e)In Figure9, we illustrate the behavior of the method when the objects are sound-hard. Computing the topological derivative when Ω=∅ (Figure 9(a)), we detect the correct position and approximate size of the two objects. The initial guess is represented in Figure 9(b). After nine steps of the iterative procedure, the reconstruction is rather satisfactory as shown in Figure 9(c).When applying the method for sound-soft objects (see Figure10(a)) the topological derivative when Ω=∅ attains the largest negative values in a region inside the object on the right. This means that in this case the topological derivative ignores the object on the left and the first approximation consists of one object, the ball Ω0, represented in Figure 10(b). In the same plot we also represent the values of the topological derivative when Ω=Ω0. The object on the right is bigger. The object on the left is ignored again. In the next iteration (see Figures 10(c)-10(d)), we find the leftmost object. After thirteen steps both objects are reconstructed with accuracy, see Figure 10(e). A number of theoretical results [47–49] suggest that good reconstructions of sound-soft objects can be achieved with just one or very few incident waves. Numerically, this was tested in [47].Some inverse scattering problems involving polarized electromagnetic waves (TM or TE) can be reformulated as constrained optimization problems similar to the ones studied in Section2. The constraints are Helmholtz type equations, eventually with complex wave numbers [18]. The methods described in the previous section apply, with a slight variation of the formulas for topological derivatives about handling constraints involving complex wave numbers [27].


## 2.1. Gradient Techniques to Identify Coefficients

Let us assume that the objectsΩ are known. The optimization problem becomes as follows: Find ki(x) minimizing
(6)J(ki)=12∫Γmeas‍|u(x)-umeas(x)|2dl,
where u(x) is the solution of the forward problem (4) with objects Ω and parameters ki.Fixed an initial guesski,0, we generate a sequence ki,m providing a decreasing sequence J(ki,m) by a gradient method [21]. At each stage, we fix the current guess ki,m and study the variations of the functional along directions v: J(η)=J(ki,m+ηv), where η>0 is a real parameter. The derivative with respect to η is given by
(7)dJ(η)dη=dJ(ki,m+ηv)dη=Re[∫Γmeas‍(uη(x)-umeas(x)¯)duη(x)dηdl],uη being the solution of the forward problem (4) with coefficient ki=ki,m+ηv inside the objects Ω. Computing duη/dη is avoided introducing an adjoint field p; see [21].Theorem 1.
The derivative with respect toη of the function J(η)=J(ki,m+ηv) defined by (6) evaluated at zero is given by
(8)dJ(η)dη|η=0=2Re[∫Ω‍vup],
where p is the solution of the adjoint problem
(9)Δp+ke2p=(umeas-u)¯δΓmeasinℝ2∖Ω¯Δp+ki2p=0inΩp-=p+,∂np-=∂np+on∂Ωlimr→∞r1/2(∂rp-ikep)=0
and u the solution of the forward problem (4) setting in both cases ki=ki,m.Choosingv in such a way that the derivative J′(0)<0, we ensure that J(ki,m+ηv)<J(ki,m) for η>0 small enough. For instance, we may set
(10)v(x)=-Re[u(x)p(x)].
The updated guess is then
(11)ki,m+1(x)=ki,m(x)-ηRe[u(x)p(x)].
This choice is reasonable if we are looking for space-dependent parameters. When Ω is composed of several disjoint domains Ωj, and ki is known to be almost piecewise constant, we may set
(12)v(x)=vj=-Re[∫Ωj‍up],x∈Ωj.
Technical details and proofs can be found in [21], where more general Helmholtz type constraints are considered:
(13)div(αe(x)∇u)+ke(x)2u=0inℝ2∖Ω¯,div(αi(x)∇u)+ki(x)2u=0inΩ,u-=u+,αe(x)∂nu-=αi(x)∂nu+on∂Ω,limr→∞r1/2(∂r(u-uinc)-iκe(u-uinc))=0.
The coefficients ke, ki, αe, αi are bounded from below by positive constants and depend smoothly on x, becoming constant at infinity. The constant κe at the Sommerfeld radiation condition is defined as κe=lim|x|→∞ke(x)/αe(x). It is the wave number of the incident wave. The interior parameter αi(x) can also be identified following the same technique.


## 2.2. Topological Derivative Techniques to Reconstruct Domains

We have seen how to generate descent directions for the parameterski fixing the objects Ω. Let us assume now that we know the parameters ki. The optimization problem becomes as follows: Find Ω minimizing
(14)J(ℝ2∖Ω¯)=12∫Γmeas‍|u(x)-umeas(x)|2dl,
where u(x) is the solution of the forward problem (4) with objects Ω and parameters ki. To generate descent directions for the domains, we must rely on some kind of derivative with respect to the domains. We choose a specific type of derivative for shape functionals: the topological derivative. It has the advantage of providing a first guess for the objects too.The topological derivative of a shape functionalJ(ℛ) at the point x∈ℛ is defined as [42]
(15)DT(x,ℛ)=limɛ→0J(ℛ∖Bɛ(x)¯)-J(ℛ)h(ɛ),
where Bɛ(x) is a ball centered at x with small radius ɛ>0; see Figure 2. The scalar function h(ɛ)>0 is chosen in such a way that the limit (15) is finite and does not vanish. For instance, for the transmission or Neumann conditions in two dimensions, we take h(ɛ)=Vol(Bɛ(x))=πɛ2. For Dirichlet conditions in 2D, we may choose h(ɛ)=-2π/log(keɛ).Figure 2
Geometry for the computation of topological derivatives.The topological derivative is a scalar function defined for everyx∈ℛ which measures the sensitivity of the functional to removing from ℛ a small object located at x. In regions where DT(x,ℛ) takes large negative values, we have a large probability to find an object. This remark is based on the following expansion:
(16)J(ℛ∖Bɛ(x)¯)=J(ℛ)+h(ɛ)DT(x,ℛ)+o(h(ɛ)),J(ℛ∖Bɛ(x)¯)=J(ℛ)h(ɛ)51..DT(x,ℛ)ɛ⟶0.
Whenever DT(x,ℛ)<0, J(ℛ∖Bɛ(x)¯)<J(ℛ) for ɛ>0 small. When J is our functional (14), this suggests a strategy to find first guesses for the objects and then improve them iteratively. Our algorithm to recover objects when the parameters are known is as follows:(i)
Compute the topological derivative whenΩ=∅ and choose an initial guess Ω0 for the objects as
(17)Ω0:={x∈ℛ∣DT(x,ℛ)≤-C0},
where C0 is a positive constant.(ii)
Check thatJ(ℛ∖Ω¯0)<J(ℛ). Otherwise increase the constant C0.(iii)
Form=1:mmax
(a)
Compute the topological derivative whenΩ=Ωm-1 and select
(18)Ωm=Ωm-1∪{x∈ℛ∖Ω¯m-1∣DT(x,ℛ∖Ω¯m-1)≤-Cm}
for decreasing thresholds Cm>0.(b)
Check thatJ(ℛ∖Ω¯m)<J(ℛ∖Ω¯m-1). Otherwise increase the constant Cm.(c)
Check the stopping criteria: the iteration stops if either meas(Ωm∖Ωm-1) is small enough, or J(ℛ∖Ω¯m) is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle |umeasɛ-um|<τɛ is satisfied, where ɛ describes measurement errors and τ>1 (in our numerical experiments τ=1.2). Here, um is the solution of the forward problem when the objects are Ωm.Computing topological derivatives at grids of points at each stage using (15) is computationally very expensive. Explicit formulas that require only the knowledge of forward and adjoint fields are used instead. For our particular Helmholtz type constraint (4), the explicit expression is [14, 21].Theorem 2.
The topological derivative ofJ(ℝ2∖Ω¯) defined by (14) is given by
(19)DT(x,ℝ2∖Ω¯)=Re[(ki2(x)-ke2)u(x)p(x)],DT(xℝ2∖Ω¯)=e[(ki2(x).-ke2)],x∈ℝ2∖Ω¯,
where u and p are solutions of the forward and adjoint problems (4) and (9), respectively.Notice that whenΩ and ki are the true objects and parameters for which the data are measured, the adjoint field p vanishes since the forward field agrees with the measured data at the receptors. Therefore, the topological derivative vanishes at the exact solution, which is a global minimum of the cost functional, as explained before.In real reconstruction tests, one usually has access to measurements for several incoming directionsdj, j=1,…,M. The shape cost functional (14) is then replaced with:
(20)J(ℝ2∖Ω¯)=12∑j=1M‍∫Γmeas|u(x;dj)-umeas(x;dj)|2dlx.
Here, u(x;dj) is the solution to the forward problem (4) with incident field uinc(x;dj)=eıkex·dj and umeas(x;dj) is the measured total field on Γmeas. The topological derivatives for this functional are obtained adding up the topological derivatives for each single incident wave.We test the algorithm for a geometrical configuration with two objects. The interior parameter in them iski=1 whereas ke=4. We have computed the topological derivative for data given at 14 sampling points uniformly distributed on the circle centered at (0,0) with radio 3, for 10 incident planar waves with uniformly distributed angles in [0,2π). The geometrical setting for our experiments is represented in Figure 3. Our synthetic data were generated by solving the corresponding forward problems by boundary element techniques using a finer mesh and adding a one percent relative random noise at each observation point.Figure 3
Geometrical description of the numerical experiments.The values of the topological derivative whenΩ=∅ at the sampling region [-1,1]×[-1,1] are represented in Figure 4(a). The largest negative values are attained inside each object. The first computation of the topological derivative provides the correct number of objects and their location, although they look like balls. Our initial guess Ω0 is represented in Figure 4(b). We compute now the topological derivative when Ω=Ω0, obtaining the color map in Figure 4(b). New regions close to the current objects appear where the topological derivative attains large negative values. Those regions should belong to the scatterers. Our objects are bigger and they are not balls. Figures 4(c)–4(j) represent the subsequent iterations. The values of the cost functional versus the number of iterations are given in Figure 4(k). After few iterations, our reconstructions are satisfactory.Reconstruction of objects whenki=1 is known using the monotone algorithm. (a) Topological derivative when Ω=∅. (b)–(j) Approximate domains Ωm superimposed to the topological derivative for Ω=Ωm, for m=0,…,8. (k) Cost functional versus the number of iterations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)The quality of the reconstructions as well as the number of iterations depends on the values of the constantsCm appearing in (18). Guidelines for their selection are given in [21]. When Cm are too big, the number of iterations required is high. A sequence of smaller Cm provides faster reconstructions, but we are prone to include spurious regions in Ωm that cannot be removed in the next iterations. We illustrate this fact in Figure 5. Figure 5(a) is exactly the same as Figure 4(a) and represents the values of the topological derivative when Ω=∅. Choosing a smaller constant C0, our initial guess Ω0, represented in Figure 5(b), has two objects that are bigger than the ones in Figure 4(b). The object on the right includes too many spurious points. Computing the topological derivative when Ω=Ω0, we detect a small region at the top of the left object that is included in the next iteration (see Figure 5(b)). No improvement for the rightmost object is pointed out. The reconstruction is given in Figure 5(c). In the next iteration (see Figure 5(d)), only the object on the left is slightly modified and the algorithm stopped because the difference between the consecutive iterations in the value of the functional was negligible, as observed in Figure 5(e). The quality of the reconstruction of the object on the left is reasonable, but, comparing with the previous example, the reconstruction of the object on the right is not acceptable.Reconstruction of objects whenki=1 is known, with the monotone algorithm (18). (a) Topological derivative when Ω=∅. (b)–(d) Approximate domains Ωm superimposed to the topological derivative for Ω=Ωm. (e) Evolution of the cost functional.
(a)
(b)
(c)
(d)
(e)The monotone algorithm cannot remove spurious regions in the rightmost object. The topological derivative is only defined at points not belonging to the current approximation. However, observing the values of the topological derivative in Figures5(b)–5(d), we realize that there are regions close to object where it takes large positive values. Computing the topological derivative when Ω=∅ (see Figure 5(a)), we did not see that effect. This indicates that spurious regions are included and that we should use a nonmonotone method if we want to correct them.The classical notion of topological derivative we used previously is defined only inℛ∖Ω¯. Therefore, any time we update our approximations Ωm, we have a criterion to add points. However, if at some stage we include in our approximation points that do not belong to the objects, we cannot remove them. This artifact is eliminated observing that for our particular type of functional we can extend the definition of topological derivative to points inside the objects. The topological derivative of a functional J(ℛ∖Ω¯) was extended to the points inside Ω in [21] as
(21)DT(x,ℛ∖Ω¯)=limɛ→0J(ℛ∖(Ω∖Bɛ(x)¯))-J(ℛ∖Ω¯)h(ɛ),DT(x,ℛ∖Ω¯)=limɛ→0J(ℛ∖(Ω∖Bɛ(x)¯))21.-Jh(ɛ)x∈Ω.
From this definition, the expansion below follows:
(22)J(ℛ∖(Ω∖Bɛ(x)¯))=J(ℛ∖Ω¯)+h(ɛ)DT(x,ℛ∖Ω¯)+o(h(ɛ)),=J(ℛ∖Ω¯)+h(ɛ)DT(x,ℛ∖1Ω¯).ɛ⟶0.
Since h(ɛ) is a positive function, if DT(x,ℛ∖Ω¯) is large and negative, the functional decreases when we remove the ball Bɛ(x) from Ω. This suggests a strategy to remove points from Ω. Our algorithm can generate a nonmonotone sequence of objects by substituting the definition of Ωm in the previous method (see (18)) with
(23)Ωm=Ωm-1∪{x∈ℛ∖Ω¯m∣DT(x,ℛ∖Ω¯m-1)≤-Cm}∖{x∈Ωm-1∣DT(x,ℛ∖Ω¯m-1)≤-cm}¯,
with decreasing sequences Cm, cm>0.Using the extended definition (21), we find the following expression for the topological derivative of the functional (14) of our transmission problem (see [21]).Theorem 3.
The extended topological derivative (21) of J(ℝ2∖Ω¯) defined by (14) is given by
(24)DT(x,ℝ2∖Ω¯)=Re[(ke2-ki2(x))u(x)p(x)],DT(Ω¯)=Re[(ke2-ki2(x))u(x)1p(x)]x∈Ω,
where u and p are solutions of the forward and adjoint problems (4) and (9), respectively.Notice that the expressions (19) and (24) are identical except for a minus sign. Moreover, since the solutions to the problems (4) and (9) are continuous functions in ℝ2, the expression
(25)Re[(ki2(x)-ke2)u(x)p(x)]
defines a continuous function in ℝ2. This motivated in [21] to modify the extended definition of the topological derivative inside Ω of a functional J(ℛ∖Ω¯) as
(26)DT(x,ℛ∖Ω¯)=limɛ→0J(ℛ∖Ω¯)-J(ℛ∖(Ω∖Bɛ(x)¯))h(ɛ),DT(x,ℛ∖Ω¯).=limɛ→0J(ℛ∖Ω¯)-J(ℛ∖(..5))h(ɛ),x∈Ω,
which is exactly the same as (21) multiplied by -1. Using this notion,
(27)J(ℛ∖(Ω∖Bɛ(x)¯))=J(ℛ∖Ω¯)-h(ɛ)DT(x,ℛ∖Ω¯)+o(h(ɛ)),=J(ℛ∖Ω¯)-5.h(ɛ)DT(x,ℛ∖Ω¯).ɛ⟶0,
and whenever DT(x,ℛ∖Ω¯) is large and positive, the functional decreases when removing the ball Bɛ(x) from Ω. This definition also allows to define a non-monotone scheme by substituting the definition of Ωm in (23) by
(28)Ωm=Ωm-1∪{x∈ℛ∖Ω¯m∣DT(x,ℛ∖Ω¯m-1)≤-Cm}∖{x∈Ωm-1∣DT(x,ℛ∖Ω¯m-1)≥cm}¯,
for decreasing thresholds Cm, cm>0.With this definition, the topological derivative of the functional for the transmission problem is as follows.Theorem 4.
The extended topological derivative (26) of J(ℝ2∖Ω¯) defined by (14) is given by
(29)DT(x,ℝ2∖Ω¯)=Re[(ki2(x)-ke2)u(x)p(x)],x∈Ω,
where u and p are solutions of the forward and adjoint problems (4) and (9), respectively.There is a missprint in [21]. Definition (21) was given, but in fact (26) was used. The topological derivative inside Ω was computed as in (29) and in the numerical examples the implemented formula is (29). The modified algorithm described by (28) was used to update the domains.Let us revisit now the example of Figure5. The topological derivative for Ω=∅ provided the initial guess Ω0 represented in Figure 6(b) by a dashed line. The values of the topological derivative inside Ω in Figure 6(b) are calculated using the definition (21), that is, formula (24). The object on the right includes spurious regions (the points with large negative values inside Ω0) that will be removed in the next iterations. Some points at the top of the object on the left should be included in that object. If we use the alternative definition (26), the topological derivative is a continuous function (see Figure 6(c)). A region inside the rightmost object where large positive values are attained should be removed. A region outside the leftmost object where it takes large negative values has to be included. The second definition provides a color map that is easier to interpret. Therefore, from now on we choose the definition (26) and update the domains using (28). The final reconstruction at the 11th iteration is represented in Figure 6(d) by a dashed red line. For comparison, we have included in Figure 6(e) the final reconstruction obtained for the same problem in Figure 4 with the monotone method (for a well-calibrated sequence of constants Cm) after nine iterations. Our non-monotone method is not very sensitive to the initial guess.Comparison of the alternative definitions of the topological derivative inside an object and of the monotone and nonmonotone methods. (a) Topological derivative whenΩ=∅. (b) Initial guess Ω0 (dashed lines) superimposed to the topological derivative for Ω=Ω0 using the extended definition (21) inside Ω. (c) Counterpart of (b) using the extended definition (26). (d) Final reconstruction at the 11th iteration using the non-monotone method (dashed red line). (e) Final reconstruction at the 9th iteration using the monotone method (dashed red line). (f) Evolution of the cost functional applying the non-monotone strategy.
(a)
(b)
(c)
(d)
(e)
(f)A wide gallery of reconstructions using the monotone and the non-monotone methods can be found in [21, 43], illustrating the ability of the method to reconstruct objects with different parameters ke, ki (simulating high and low frequencies), problems with different values of ki inside each object, scatterers of different sizes, different levels of noise in the data, reconstructions when data are given on a limited region of a circle, source point excitations instead of planar incident waves, and so forth.


## 2.3. Hybrid Approach to Reconstruct Objects and Identify Their Material Parameters

Combining the two strategies to approximate objects and parameters discussed in the previous two sections, we obtain a reconstruction scheme [21, 22]. (i)
Initialization
(a)
Choose an initial guesski,0 for the parameters ki, for example, a perturbation of the background parameter ke.(b)
Choose an initial guessΩ0 for the objects Ω, using the topological derivative: set ki=ki,0 and compute the topological derivative with object Ω=∅:
(30)DT(x,ℝ2)=Re[(ki,02-ke2)u(x)p(x)],whereu and p are forward and adjoint fields with object Ω=∅ and parameter ki,0. Then,
(31)Ω0={x∈ℝ2∣DT(x,ℝ2,ki,0)≤-C0},whereC0 is a certain large positive threshold.(ii)
Iteration
(a)
Updateki,m(x) using the gradient technique:
(32)ki,m+1(x)=ki,m(x)-ηRe[u(x)p(x)],whereu, p solve forward and adjoint problems with objects Ωm and parameter ki,m(x).(b)
Check thatJ(ℝ2∖Ω¯m,ki,m+1)<J(ℝ2∖Ω¯m,ki,m). Otherwise, reduce η.(c)
UpdateΩm computing the topological derivative with objects Ωm and parameter ki,m+1(x). We may resort to a monotone
(33)Ωm+1=Ωm∪{x∈ℝ2∖Ω¯m∣DT(x,ℝ2∖Ω¯m,ki,m+1)≤-Cm+1},or nonmonotone strategy(34)Ωm+1=Ωm∪{x∈ℝ2∖Ω¯m∣DT(x,ℝ2∖Ω¯m,ki,m+1)=Ωm∪55x∈ℝ2∖Ω¯m∣DT(x,ℝ2∖Ω¯m,ki,m+1)≤-Cm+1}∖{x∈Ωm∣DT(x,ℝ2∖Ω¯m,ki,m+1)≥cm+1}¯,with decreasing thresholdsCm+1, cm+1. The second one is more convenient if we wish to be able to remove spurious points at any stage, for instance, when reconstructing objects with holes.(d)
Check thatJ(ℝ2∖Ω¯m+1,ki,m+1)<J(ℝ2∖Ω¯m,ki,m+1). Otherwise, increase the thresholds cm+1, Cm+1.(e)
Stop when meas(Ωm∖Ωm-1) is small enough, or J(ℝ2∖Ω¯m) is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle |umeasɛ-u|<τɛ is satisfied.The iterative procedure described perviously generates sequences of domainsΩm and parameters ki,m along which the functional decreases. The magnitude of the topological derivative is also observed to decrease. Theoretically, convergence to the true solution is not guaranteed, since our cost functional might have local minima different from the global minimum we are looking for. This situation can be improved or worsened by the amount of incident waves we use, the number of receptors, and how we distribute them. In the tests we have performed, the thresholds were calibrated so that the scheme worked quite well in a few iterations. See [21, 22] for details about calibration, numerical schemes and extensions to constraints of the form (13).The forward and adjoint problems have explicit solutions whenke is constant and we are computing the first guess Ω0 on the whole space. They become
(35)Δu+ke2u=0inℝ2limr→∞r1/2(∂r(u-uinc)-ike(u-uinc))=0,(36)Δp+ke2p=(umeas-u)¯δΓmeasinℝ2limr→∞r1/2(∂rp-ikep)=0.
The solution to (35) is the incident wave and the solution to (36) can be computed using the fundamental solution of the Helmholtz equation. The forward and adjoint problems (4) and (9) are solved by either boundary element or finite element methods depending on whether ki is piecewise constant or vary smoothly with space [11, 13, 21, 44–46].In practice, updating the domains is more expensive than updating the parameters. Any time we change the guess of the objects, we must mesh again the new domains to be able to compute forward and adjoint fields. Therefore, we usually correct the parameterski,m several times for a fixed approximate domain Ωm before proceeding to update Ωm, to reduce the computational cost.We have tested the hybrid topological derivative-gradient based method considering the configuration with two objects of the previous examples, keeping the same observation points and incident waves (see Figure3). We also take the same values of the parameters: ke=4 is assumed to be known and ki=1 is assumed to be unknown. We computed the topological derivative when Ω=∅ choosing ki,0=3 as initial guess for ki. Notice that since the forward and adjoint fields when Ω=∅ do not depend on ki, the topological derivative is the same as in Figure 4(a) except for the scale (we have now in (19) the multiplying factor ((ki,0)2-ke2) instead of (ki2-ke2)). Therefore, our initial guess Ω0, represented in Figure 7(a), is the same as in Figure 4(b). We fix now Ω=Ω0 and iterate five times to update the value of ki without modifying the scatterers. Then, we fix the value of ki and compute the topological derivative to update the domains and so on. The reconstruction at the ninth iteration with respect to the domain is shown in Figure 7(b). The values of ki through the iterative procedure are represented in Figure 7(c). Two identical values of the parameter in the plot mark each iteration to improve the domains (ki is not updated). Our final approximation is 1.23 while the true value is 1. Some reconstructions with two or three objects with different constant parameters can be found in [21, 43].Reconstruction of objects with unknownki. (a) First guess of the objects. (b) Final reconstruction of the objects after nine updates of the domain (c) Values of ki through the iterative method.
(a)
(b)
(c)In the next example we test the hybrid algorithm for the reconstruction of an objectΩ with space-dependent coefficient ki. The scatterer is the ball centered at (-0.25,0.3) and radius 0.35 and ki is the function represented in Figure 8(c). It grows radially from 1 to 2 inside Ω. The observation points and the incident waves are the same as in all the previous examples. To start the algorithm, we took as initial guess the constant function ki,0=3. The topological derivative for Ω=∅ provided the initial guess Ω0 represented by a red dashed line in Figure 8(a). We applied the hybrid method alternating one iteration with respect to the domain with five steps of the gradient method. The final reconstructions of Ω and ki are shown in Figures 8(b) and 8(d), respectively. The parameter ki is overestimated with an average error about 0.2, but we clearly recognize a function that grows radially. The location, shape, and size of the object are reconstructed with accuracy. The interested reader can find in [22] some numerical experiments including the reconstruction of two homogeneous objects immersed in an heterogeneous material and of two heterogeneous scatterers buried in an homogeneous media.Reconstruction an object with space-dependent parameterki. (a) True object (blue solid line) and initial guess (red dashed line). (b) True object (blue solid line) and its final reconstruction (red dashed line) at the ninth iteration with respect to the domain. (c) True function ki. (d) Final reconstruction of ki.
(a)
(b)
(c)
(d)


## 2.4. Sound-Hard and Sound-Soft Objects

The methods presented in the previous sections apply when the objects are penetrable by the incident radiation. Part of the incident wave is transmitted inside the object and transmission boundary conditions are imposed at the interface between the objects and the surrounding medium.When the scatterers are sound-hard (rigid), no transmitted wave is generated. The incident radiation is scattered and a Neumann boundary condition is imposed at the surface of the objects. The constraint for the cost functional (14) becomes
(37)Δu+ke2u=0,inℝ2∖Ω¯∂nu=0,on∂Ωlimr→∞r1/2(∂r(u-uinc)-ike(u-uinc))=0.
Topological derivative methods to find first guesses of the number, shape, and location of rigid objects were developed in [13, 44]. However, to the best of the authors’ knowledge, iterative methods entirely based on the computation of topological derivatives for the reconstruction of sound-hard objects had not been previously used.For sound-soft scatterers, a Dirichlet boundary condition is imposed at the interface between the objects and the background medium:(38)Δu+ke2u=0,inℝ2∖Ω¯u=0,on∂Ωlimr→∞r1/2(∂r(u-uinc)-ike(u-uinc))=0.
Topological derivative methods were first employed to reconstruct sound-soft objects in [16, 20, 47].For rigid and sound-soft objects, we cannot predict their material parameters, since they do not appear in the constraints. The cost functionals only depend on the geometry of the inclusionsΩ. The topological derivative based iterative methods sketched in the previous sections can be extended to these two types of objects, replacing the explicit expressions of the topological derivatives with other adequate formulas, which involve modified forward and adjoint fields. The forward and adjoint problems are formally similar to the transmission ones, suppressing the equation for the transmitted wave inside Ω and imposing either Neumann or Dirichlet boundary conditions on ∂Ω.The new expressions for the topological derivatives are:(i)
For the Neumann problem as follows(39)DT(x,ℝ2∖Ω¯)=Re[2∇u(x)∇p(x)-ke2u(x)p(x)],
where the forward field is a solution of (37) and the adjoint field is a solution of
(40)Δp+ke2p=(umeas-u)¯δΓmeas,inℝ2∖Ω¯,∂np=0,on∂Ω,limr→∞r1/2(∂rp-ikep)=0.
This formula was obtained in [21], correcting the previous paper [13] where the multiplying factor 2 in (39) was missed.(ii)
For the Dirichlet problem (see [47])
(41)DT(x,ℝ2∖Ω¯)=Re[u(x)p(x)],
where the forward field is a solution of (38) and the adjoint field is a solution of
(42)Δp+ke2p=(umeas-u)¯δΓmeas,inℝ2∖Ω¯,p=0,on∂Ω,limr→∞r1/2(∂rp-ikep)=0.We test the performance of our iterative topological derivative method for rigid and sound-soft objects in Figures9 and 10, respectively. The geometrical setting, sketched in Figure 3, is the same as in the previous examples for the transmission problem with constant coefficients. We also keep the same observation points and incident waves. The (exterior) wave number is ke=4 too.Reconstruction of rigid objects. (a) Topological derivative whenΩ=∅. (b) First guess of the objects (dashed red line). (c) Final reconstruction of the objects after 9 steps (dashed red line).
(a)
(b)
(c)Reconstruction of sound-soft objects. (a) Topological derivative whenΩ=∅. (b) Initial guess Ω0 superimposed to the topological derivative for Ω=Ω0. (c)-(d) Approximate domains Ωm superimposed to the topological derivative for Ω=Ωm, for m=1 and m=2. (e) Final reconstruction of the objects after 13 steps (dashed red line).
(a)
(b)
(c)
(d)
(e)In Figure9, we illustrate the behavior of the method when the objects are sound-hard. Computing the topological derivative when Ω=∅ (Figure 9(a)), we detect the correct position and approximate size of the two objects. The initial guess is represented in Figure 9(b). After nine steps of the iterative procedure, the reconstruction is rather satisfactory as shown in Figure 9(c).When applying the method for sound-soft objects (see Figure10(a)) the topological derivative when Ω=∅ attains the largest negative values in a region inside the object on the right. This means that in this case the topological derivative ignores the object on the left and the first approximation consists of one object, the ball Ω0, represented in Figure 10(b). In the same plot we also represent the values of the topological derivative when Ω=Ω0. The object on the right is bigger. The object on the left is ignored again. In the next iteration (see Figures 10(c)-10(d)), we find the leftmost object. After thirteen steps both objects are reconstructed with accuracy, see Figure 10(e). A number of theoretical results [47–49] suggest that good reconstructions of sound-soft objects can be achieved with just one or very few incident waves. Numerically, this was tested in [47].Some inverse scattering problems involving polarized electromagnetic waves (TM or TE) can be reformulated as constrained optimization problems similar to the ones studied in Section2. The constraints are Helmholtz type equations, eventually with complex wave numbers [18]. The methods described in the previous section apply, with a slight variation of the formulas for topological derivatives about handling constraints involving complex wave numbers [27].


## 3. Methods for Impedance Tomography

In this section, we deal with another imaging technique that uses time-harmonic electromagnetic waves: impedance imaging.A derivation of the impedance imaging problem from Maxwell’s equations is given in [29]. Inside the body ℛ, the electric potential u(x) satisfies
(43)∇·γ(x,ω)∇u=0,
where γ is the admittivity, γ=σ+ıωɛ, σ>0 is the electric conductivity, ɛ is the electric permittivity, and ω is the angular frequency of the applied current. When ω=0, we are left with an inverse conductivity problem.Currentsj are applied through electrodes located on the boundary ∂ℛ of the body, where the resulting voltage is measured. In real experiments, only the values of the current and the voltages at a “discrete” set of electrodes are known. Realistic models supplement (43) with a “discrete” Neumann boundary condition:
(44)γ∂nu=I,on∂ℛ,
where I satisfies
(45)I=0,on∂ℛ∖∪k=1Ne¯k,∫ek‍Idl=Ik,k=1,…,N,∑k=1N‍Ik=0,Ik is the current sent to the kth electrode and ek denotes the part of ∂ℛ that corresponds to the kth electrode. We assume that each ek is an open subset of ∂ℛ with positive surface measure and that dist(ek,ej)>0 for k≠j. ℛ and Ω have smooth boundaries. γe is continuous up to the boundary, γi is continuous, and Umeas and I define H1/2(∂ℛ) functions.At the electrodes, we measure the voltagesumeas,k. Dirichlet conditions
(46)u=umeas,kattheelectrodesek,k=1,…,N,
have been shown not to be realistic [29] and are commonly replaced with the constraints
(47)u+zkγ∂nu=umeas,k,attheelectrodesek,k=1,…,N,
where zk is the effective contact impedance or surface impedance and
(48)∑k=1N‍umeas,k=0.
As shown in [50], problem (43)-(44) admits the variational formulation: Find (u,U)∈H1(ℛ)⊕ℝ0N such that
(49)b((u,U),(w,W))=ℓ(w,W)∀(w,W)∈H1(ℛ)⊕ℝ0N,
where
(50)b((u,U),(w,W)):=∫ℛ‍γ∇u∇w¯dz+∑k=1N‍1zk∫ek‍(u-Uk)(w¯-W¯k)dl,ℓ(w,W):=∑k=1N‍IkW¯k.
The subscript 0 refers to vectors of vanishing mean; that is, ∑k=1N‍Uk=0. Existence and uniqueness are guaranteed by the Lax-Milgram lemma for bounded coefficients and smooth domains. A discretization scheme by FEM has been proposed in [51]; see also references therein.A simpler model imposes a “continuous” Neumann boundary condition:(51)γ∂nu=j,on∂ℛ,
and measures the voltage on the boundary
(52)u=umeas,on∂ℛ.Here,n denotes the outward unit normal and j stands for the outward pointing normal component of the current density generated by the electrodes on the surface. This Neumann problem (43)–(51) admits solutions when the current satisfies the compatibility condition ∫∂ℛ‍jdl=0, which turns to be the law of conservation of charge. To select a unique potential u, we impose ∫∂ℛ‍udl=0. The measured voltage also satisfies the condition ∫∂ℛ‍umeasdl=0.The goal is to reconstruct the structure of the admittivityγ inside ℛ from measurements at the boundary. Two slightly different settings may be considered. If ℛ contains a number of inclusions Ωj, the admittivity γ is a piecewise smooth function in ℛ with discontinuities at the boundaries of the inclusions. The admittivity outside the objects γe is known. The admittivity inside γi is unknown and must be reconstructed, together with the geometry of Ω, also unknown. Otherwise, one may consider admittivities γ to vary smoothly over ℛ, being only known at its boundary. For brevity, we assume here that the background contains a number of inclusions. The second setting was studied in [23]. First, we illustrate the idea applying hybrid topological derivative-gradient based techniques to the simplified model (43)–(51) following [23]. Reasonable reconstructions of the geometry of the objects are obtained in a few iterations, though the values of the admittivity inside inclusions might be better estimated replacing the gradient approach with other techniques [37–39]. Then, we extend the theory to the realistic tomography problem with discrete electrodes, much less studied.


### 3.1. Tomography with the Simplified Model

The “continuous” impedance tomography problem can be reformulated as an optimization problem. We assume thatℛ contains a number of inclusions Ωj and the admittivity γ is a piecewise function in ℛ with discontinuities at the boundaries of the inclusions. For simplicity, we take ℛ to be a bounded set of ℝ2 with smooth boundary. We set Ω=∪j=1dΩj with Ωj simply connected open bounded sets with smooth boundaries without pairwise contact; see Figure 11. The admittivity γe in the background medium ℛ∖Ω¯ is known and γi is unknown inside Ω. The functional to be optimized becomes
(53)J(ℛ∖Ω¯,γi)=12∫∂ℛ‍|u-umeas|2dl,
where u is a solution of
(54)∇·γe∇u=0inℛ∖Ω¯,∇·γi∇u=0inΩ,u--u+=0on∂Ω,γi∂nu--γe∂nu+=0on∂Ω,γe∂nu=jon∂ℛ.
The unit normal n points outside ℛ∖Ω¯ but inside Ω. u- and u+ denote the limit values of u on ∂Ω from outside and inside Ω, respectively; see Figure 11.Figure 11
Geometry for the impedance tomography problem.As observed in Section2, to apply hybrid topological derivative-gradient based methods to this functional, we need two ingredients: simple expressions for its topological derivatives and formulas for variations with respect to the admittivity. Formulas for topological derivatives of this functional are given in [23, 32, 33].Theorem 5.
The topological derivative of the cost functional (53) is
(55)DT(x,ℛ∖Ω¯,γi)=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)∇u(x)∇p¯(x)],=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)u]x∈ℛ∖Ω¯,(56)DT(x,ℛ∖Ω¯,γi)=Re[2γi(x)(γi(x)-γe(x))γi(x)+γe(x)∇u(x)∇p¯(x)],x∈Ω,
where u is the solution of the forward problem (54) and p is a solution of the adjoint problem
(57)∇·γ¯e∇p=0inℛ∖Ω¯,∇·γ¯i∇p=0inΩ,p--p+=0on∂Ω,γ¯i∂np--γ¯e∂np+=0on∂Ω,γ¯e∂np=umeas-uon∂ℛ.In [23] there is a misprint in the definition of the topological derivative inside Ω. It is written in (21) but there is a missing minus sign. The formula considered to obtain (56) was (26).When we fixΩ and perturb guesses of the admittivity γi,m in a certain direction, the derivatives of the resulting functional are given below. See [23] for the proof. Related differential calculus and formulas are presented in [28, 37, 52].Theorem 6.
The derivative with respect toη of the function J(η):=J(ℛ∖Ω¯,γi,m+ηv) with J defined in (53) is
(58)dJ(γi,m+ηv)dη|η=0=Re[∫Ω‍v∇u∇p¯dz],
where u and p are solutions of the forward (54) and adjoint (57) problems with objects Ω and admittivity γi,m.Choosing(59)v(x)=-∇u¯(x)∇p(x),x∈Ω,
we ensure J(γi,m+ηv)<J(γim) if η>0 is small enough. We have defined v only in Ω, but formula (59) holds for all x∈ℛ. When γi is piecewise constant, we expect γi,m to be almost piecewise constant. In this case, another choice for the corrector v is a piecewise constant approximation:
(60)v(x)=vj=-Re[∫Ωj‍∇u¯∇pdz],x∈Ωj,j=1,…,d.The reconstruction algorithm for the impedance tomography problem is similar to the algorithm described in Section2.3 with γi, γe playing the role of ki, ke, and using (55)-(56) for the topological derivatives, (59) or (60) for the parameter correctors. We test this method in the next experiment where ℛ is the unit ball. Data are given at 30 uniformly distributed electrodes on ∂ℛ (represented by crosses in Figure 12) for the ten current patterns
(61)jℓ(t)=cos(ℓt),jℓ+5=sin(ℓt),ℓ=1,…,5,jℓ(t)=cos(ℓt),.jℓ+5=sin(ℓt),t∈[0,2π).
We consider complex admittivities of the form γ=σ+ıωɛ with ω=1. Inside the two objects represented by solid blue lines in Figures 12(a)-12(b), the unknown admittivity is γi=8+2ı and γe=1+ı is assumed to be known. We started the hybrid algorithm with γi,0=5+3ı. The first guess for the domains was obtained observing the topological derivative with Ω=∅ and γi=γi,0. It has only one object, as shown in Figure 12(a). We applied the algorithm alternating one iteration with respect to the domain with five iterations with respect to the parameter. Although we initialized the method with an incorrect number of objects, after 17 iterations with respect to the domain we obtained a reasonable reconstruction of the two defects; see Figure 12(b). When the algorithm stops, the true admittivity γi=8+2ı is approximated by 7.63+2.30ı. The values of the conductivity σi (the real part of γi) and of the permittivity ɛi (the imaginary part of γi) throughout the iterative procedure are represented in Figures 12(c) and 12(d), respectively. We refer to [23] for an example of the reconstruction of an object with unknown spatial dependent admittivity. The interested reader can also find in [23] a gallery of reconstructions mainly focused on the conductivity problem (with γ∈ℝ).Reconstruction of objects by impedance tomography. (a) Initial guess. (b) Final reconstruction of the objects after 17 updates with respect to the domain. (c) Convergence of the conductivity. (d) Convergence of the permittivity.
(a)
(b)
(c)
(d)


### 3.2.  Tomography with Discrete Electrodes

The impedance tomography problem with discrete electrodes (43)-(44) can be reformulated as an optimization problem: Minimize
(62)J(ℛ∖Ω¯,γi)=12∑k=1N‍∫ek‍|u+zkγe∂nu-umeas,k|2dl,
where u is a solution of
(63)∇·γe∇u=0inℛ∖Ω¯,∇·γi∇u=0inΩ,u--u+=0on∂Ω,γe∂nu--γi∂nu+=0on∂Ω,γe∂nu=0on∂ℛ∖∪j=1Ne¯j,∫ek‍γe∂nudl=Ikk=1,…,N.As before, we assume thatℛ contains inclusions Ω of a different material. For simplicity, we work in the two dimensional setting, but the method remains valid in 3D with straightforward modifications. The admittivity of the background medium γe is known. To determine its spatial variation inside the inclusions we must optimize (62) with respect to the interior domains and their admittivity. As shown before, even if the admittivity is unknown, we may locate the regions where it undergoes noticeable changes by computing the topological derivative of the shape functional in the whole domain. Once a first guess for Ω is found, we may approximate the admittivity γi optimizing with respect to γi for Ω fixed by a gradient technique. The guesses for Ω can be updated computing new topological derivatives. Below, we give expressions for the correctors to be used in the gradient method and for the topological derivatives.Let us fixΩ and seek to minimize J(γi)=J(ℛ∖Ω¯,γi) when u is a solution of the forward problem with discrete electrodes (63). When the nature of the inclusions we look for is known, an initial approximate value γi,0 is usually available. Otherwise, we select a first approximation by perturbing γe:γi,0=γe+ϵ. We may improve this first approximation iteratively by a gradient method. Given a guess γi,m, we generate a descent direction v to improve this guess studying the function of a real variable J(η)=J(ℛ∖Ω¯,γi,m+ηv). As in the acoustic setting and in the “continuous” model, η>0 and v are selected in such a way that J′(0)<0. An explicit expression for this derivative provides the following formula for the descent directions v.Theorem 7.
The derivative with respect toη of the function J(η)=J(ℛ∖Ω¯,γi,m+ηv), with J defined in (62) is
(64)dJ(η)dη|η=0=Re[∫Ω‍v∇u∇p¯dz],
where u is a solution of the forward problem (63) with parameter γi=γi,m in Ω and p is a solution of the adjoint problem (65) given below
(65)(p,P)∈H1(ℛ)⊕ℝ0N,(66)b((w,W),(p,P))=ℓa(w,W),∀(w,W)∈H1(ℛ)⊕ℝ0N,
where
(67)b(Ω;(w,W),(p,P)):=∫R∖Ω¯‍γe∇w∇p¯dz+∫Ω‍γi∇w∇p¯dz+∑k=1N‍1zk∫ek‍(w-Wk)(p¯-P¯k)dl,(68)ℓa(w,W):=-∑k=1N‍∫ek‍(u+zkγe∂nu-umeas,k¯)(w+zkγe∂nw)dl,
with γi=γi,m in Ω. Choosing
(69)v(x)=-Re[∇u¯(x)∇p(x)],x∈Ω
and η>0 small enough, one ensures J(γi,m+ηv)<J(γi,m).Proof.
The derivative with respect toη is
(70)ddηJ(γi,m+ηv)|η=0=Re[∑k=1N‍∫ek‍(u+zkγ∂nu-umeas,k¯)(u˙+zkγe∂nu˙)dl],
where u˙=(d/dη)uη|η=0 and uη is a solution of (63) with γi=γi,m+ηv; that is,
(71)(uη,Uη)∈H1(ℛ)⊕ℝ0Nb(η;(uη,Uη),(w,W))=ℓ(w,W),=ℓ(w,W)∀(w,W)∈H1(ℛ)⊕ℝ0N,
where
(72)b(η;(u,U),(w,W)):=∫ℛ∖Ω¯‍γe∇u∇w¯dz+∫Ω‍(γi,m+ηv)∇u∇w¯dz+∑k=1N‍1zk∫ek‍(u-Uk)(w¯-W¯k)dl,
and ℓ is defined in (50). Notice that when η=0 and γi=γi,m, the sesquilinear forms (67) and (72) are the same.
We eliminateu˙ by introducing the adjoint problem. Let us set
(73)ℒ(η;(u,U),(p,P)):=J(η)+Re[b(η;(u,U),(p,P))-ℓ(p,P)],
for (u,U),(p,P)∈H1(ℛ)⊕ℝ0N. Then, choosing (u,U)=(uη,Uη), it follows that
(74)J(η)=ℒ(η;(uη,Uη),(p,P)),∀(p,P)∈H1(ℛ)⊕ℝ0N
and therefore
(75)ddηJ(η)|η=0=ddηℒ(η;(uη,Uη),(p,P))|η=0=Re[ddηb(η;(u,U),(p,P))|η=0]+Re[b(0;(u˙,U˙),(p,P))]+Re[∑k=1N‍∫ek(u+zkγe∂nu-umeas,k¯)×(u˙+zkγ∂nu˙)dl{∑k=1N‍∫ek(u+zkγe∂nu-umeas,k¯)}],
for all (p,P)∈H1(ℛ)⊕ℝ0N. Selecting as (p,P) a solution of (65), we obtain the desired formula
(76)ddηJ(η)|η=0=Re[ddηb(η;(u,U),(p,P))]|η=0=Re∫Ω‍v∇u∇p¯dz.Notice that ifw∈H1(ℛ), then w∈L2(∂ℛ)∩H1/2(∂ℛ) and ∂w/∂xi∈H-1/2(∂ℛ). Thus, the integrals appearing in the definition (68) of ℓa are in fact duality products.We have definedv only in Ω, but formula (69) holds for all x∈ℛ. When γi is piecewise constant, we replace v defined in (69) with piecewise constant approximations inside each object Ωj,
(77)v(x)=vj=-Re[∫Ωj‍∇u¯∇pdz],x∈Ωj,j=1,…,d.
An initial approximation of the objects buried in the medium as well as their successive corrections is computed using the topological derivative. For the tomography problem with discrete electrodes, the explicit expression for the topological derivative in terms of forward and adjoint fields is given in the theorem below. Combining the expressions for the approximation of the admittivity γi with the topological derivative computation to approximate the domains Ω, we can define a hybrid reconstruction scheme for the discrete electrode model as in the previous sections. The implementation of the numerical codes is out of the scope of this paper and will be the subject of our future work.Theorem 8.
The topological derivative of the cost functionalJ(ℛ∖Ω¯;γi) defined in (62) is given by
(78)DT(x,ℛ∖Ω¯,γi)=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)∇u(x)∇p¯(x)],55=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)]x∈ℛ∖Ω¯,DT(x,ℛ∖Ω¯,γi)=Re[2γi(x)(γi(x)-γe(x))γi(x)+γe(x)∇u(x)∇p¯(x)],5=Re[2γi(x)(γi(x)-γe(x))γi(x)+γe(x)∇.(x)]x∈Ω,
where u and p solve the forward and adjoint problems (63) and (65) with coefficient γi in Ω.Proof.
We adapt the strategy used in [13, 21] for inverse scattering problems. The topological derivative of a shape functional can be computed as a limit of shape derivatives [13]. Consider the family of deformations depending on the real parameter τ,
(79)φτ(x):=x+τV(x),x∈ℝ2,
where the vector field V is defined as
(80)V(z)=-n(z),z∈∂Bɛ(x),
and extended to ℝ2 in such a way that it vanishes away from a narrow neighborhood of ∂Bɛ(x). Then,
(81)DT(x,ℛ)=limɛ→01h′(ɛ)ddτJ(φτ(ℛ∖Bɛ(x)¯))|τ=0,x∈ℛ,
where h′(ɛ) is the derivative of the function h(ɛ) that appears in the definition of the topological derivative (15). For the electrical impedance tomography problem in 2D, we can take h(ɛ)=πɛ2 as we will shortly see.
The proof is organized in three steps. We assume first thatγe and γi are constant and compute the topological derivative when Ω=∅. Afterwards, we adapt the proof to smoothly varying coefficients. Finally we consider the case Ω≠∅.
Step
1 (computation of the topological derivative whenΩ=∅andγeandγiare constant parameters). Let us compute first
(82)ddτJ(φτ(ℛ∖Bɛ(x)¯))|τ=0.
Since the field V decays fast enough away from ∂Bɛ(x), it follows that φτ(∂ℛ)=∂ℛ and φτ(ℛ∖Bɛ(x)¯)=ℛ∖Ω¯τ, where Ωτ:=φτ(Bɛ(x)). The cost functional in the transformed domains is then
(83)J(ℛ∖Ω¯τ)=12∑k=1N‍∫ek‍|uτ+zkγe∂nuτ-umeas,k|2dl,
where uτ solves
(84)(uτ,Uτ)∈H1(ℛ)⊕R0N,b(Ωτ;(uτ,Uτ),(w,W))=ℓ(w,W),∀(w,W)∈H1(ℛ)⊕R0N,
with b defined in (67) and ℓ defined in (50). Notice that when τ=0, Ωτ=Bɛ(x), and the solution of (84) is the solution (uɛ,Uɛ) to (63) with Ω=Bɛ(x).
The derivative of (83) with respect to τ involves the derivative u˙ of uτ with respect to τ. To avoid its computation we introduce the Lagrangian functional
(85)ℒ(Ωτ;(u,U),(p,P)):=J(ℛ∖Ω¯τ)+Re[b(Ωτ;(u,U),(p,P))-ℓ(p,P)].
Then, for all (p,P)∈H1(ℛ)⊕R0N,
(86)ddτJ(ℛ∖Ω¯τ)|τ=0=ddτℒ(Ωτ;(uτ,Uτ),(p,P))|τ=0=Re[ddτb(Ωτ;(uɛ,Uɛ),(p,P))|τ=0]+Re[b(Bɛ(x);(u˙,U˙),(p,P))]+Re[∑k=1N‍∫ek(uɛ+zkγe∂nuɛ-umeas,k¯)×(u˙+zkγe∂nu˙)dl{∑k=1N‍∫ek(uɛ+zkγe∂nuɛ-umeas,k¯)}].
If (pɛ,Pɛ) satisfies (65) with object Ω=Bɛ(x) and parameter γi, then
(87)ddτJ(ℛ∖Ω¯τ)|τ=0=Re[ddτb(Ωτ;(uɛ,Uɛ),(pɛ,Pɛ))|τ=0]=Re[∫ℛ∖Bɛ(x)¯‍γe(∇uɛ·∇p¯ɛ)∇·Vdz-∫ℛ∖Bɛ(x)¯‍γe((∇V+∇V⊤)∇uɛ)·∇p¯ɛdz+∫Bɛ(x)‍γi(∇uɛ·∇p¯ɛ)∇·Vdz-∫Bɛ(x)‍γi((∇V+∇V⊤)∇uɛ)·∇p¯ɛdz]=Re[∫∂Bɛ(x)‍(γi(γiγe-1)∂nuɛ-∂np¯ɛ-+(γi-γe)∂tuɛ-∂tp¯ɛ-{(γiγe-1)})dl{∫∂Bɛ(x)‍}].
To obtain the last equality in (87), we have integrated by parts and used that uɛ and pɛ solve (63) and (65) with Ω=Bɛ(x), respectively, and that V vanishes on ∂ℛ.
By (81), the topological derivative at x can be computed now as
(88)DT(x,ℛ)=limɛ→01h′(ɛ)Re[∫∂Bɛ(x)‍(γi(γiγe-1)∂nuɛ-∂np¯ɛ-+(γi-γe)∂tuɛ-∂tp¯ɛ-{(γiγe-1)})dl{∫∂Bɛ(x)‍}].
The asymptotic behavior of uɛ and pɛ at ∂Bɛ(x) is obtained expressing them as corrections of u0 and p0, the solutions to (63) and (65) with Ω=∅,
(89)uɛ(z)=u0(z)χℛ∖Bɛ(x)¯(z)+vɛ(z),pɛ(z)=p0(z)χℛ∖Bɛ(x)¯(z)+qɛ(z),
and expanding the remainders in powers of ɛ. Reasoning as in [23] we find that
(90)uɛ(z)⟶u0(x),∇uɛ-(z)⟶2γeγe+γi∇u0(x),uɛ(z)⟶u0(x),55.∇uɛ-(z)⟶asɛ⟶0,
uniformly when |z-x|=ɛ. For pɛ there is a slight technical difference because pɛ is defined by the variational problem (65). However, this is due to the source term at ∂ℛ, that is included in p0. Once that source term is eliminated, the remainder qɛ(z) satisfies equations similar to the equations for vɛ and we find
(91)pɛ(z)⟶p0(x),∇pɛ-(z)⟶2γeγe+γi∇p0(x),pɛ(z)⟶p0(x),∇pɛ-(z)5.0⟶asɛ⟶0,
uniformly when |z-x|=ɛ.
Replacing the limits (90) and (91) in (88) and choosing h(ɛ)=πɛ2, formula (78) with Ω=∅ follows.
Step
2 (variable parameters). The previous derivation holds when the parameters are constant. When they vary smoothly, we expand them about x and repeat the previous computations. Notice that due to the choice of V, only the values of γi and γe in a narrow band about ∂Bɛ(x) are relevant for the computation of the shape derivative. For the topological derivative, we set γi(z)=γi(x)+(z-x)·∇γi(x)+O(ɛ2) and an analogous expression for γe. Then, we check that the limits (90)-(91) also hold with γi=γi(x) and γe=γe(x). For continuous γe and γi, we use a regularizing approximation, see [21] for details.
Step
3 (topological derivative with inclusionsΩ). We have to redo the previous steps taking as objects Ω∖Bɛ(x)¯ for x∈Ω, and Ω∪Bɛ(x) when x∈ℛ∖Ω¯ depending on whether we are computing the derivative inside or outside Ω; see [21] for details.


## 3.1. Tomography with the Simplified Model

The “continuous” impedance tomography problem can be reformulated as an optimization problem. We assume thatℛ contains a number of inclusions Ωj and the admittivity γ is a piecewise function in ℛ with discontinuities at the boundaries of the inclusions. For simplicity, we take ℛ to be a bounded set of ℝ2 with smooth boundary. We set Ω=∪j=1dΩj with Ωj simply connected open bounded sets with smooth boundaries without pairwise contact; see Figure 11. The admittivity γe in the background medium ℛ∖Ω¯ is known and γi is unknown inside Ω. The functional to be optimized becomes
(53)J(ℛ∖Ω¯,γi)=12∫∂ℛ‍|u-umeas|2dl,
where u is a solution of
(54)∇·γe∇u=0inℛ∖Ω¯,∇·γi∇u=0inΩ,u--u+=0on∂Ω,γi∂nu--γe∂nu+=0on∂Ω,γe∂nu=jon∂ℛ.
The unit normal n points outside ℛ∖Ω¯ but inside Ω. u- and u+ denote the limit values of u on ∂Ω from outside and inside Ω, respectively; see Figure 11.Figure 11
Geometry for the impedance tomography problem.As observed in Section2, to apply hybrid topological derivative-gradient based methods to this functional, we need two ingredients: simple expressions for its topological derivatives and formulas for variations with respect to the admittivity. Formulas for topological derivatives of this functional are given in [23, 32, 33].Theorem 5.
The topological derivative of the cost functional (53) is
(55)DT(x,ℛ∖Ω¯,γi)=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)∇u(x)∇p¯(x)],=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)u]x∈ℛ∖Ω¯,(56)DT(x,ℛ∖Ω¯,γi)=Re[2γi(x)(γi(x)-γe(x))γi(x)+γe(x)∇u(x)∇p¯(x)],x∈Ω,
where u is the solution of the forward problem (54) and p is a solution of the adjoint problem
(57)∇·γ¯e∇p=0inℛ∖Ω¯,∇·γ¯i∇p=0inΩ,p--p+=0on∂Ω,γ¯i∂np--γ¯e∂np+=0on∂Ω,γ¯e∂np=umeas-uon∂ℛ.In [23] there is a misprint in the definition of the topological derivative inside Ω. It is written in (21) but there is a missing minus sign. The formula considered to obtain (56) was (26).When we fixΩ and perturb guesses of the admittivity γi,m in a certain direction, the derivatives of the resulting functional are given below. See [23] for the proof. Related differential calculus and formulas are presented in [28, 37, 52].Theorem 6.
The derivative with respect toη of the function J(η):=J(ℛ∖Ω¯,γi,m+ηv) with J defined in (53) is
(58)dJ(γi,m+ηv)dη|η=0=Re[∫Ω‍v∇u∇p¯dz],
where u and p are solutions of the forward (54) and adjoint (57) problems with objects Ω and admittivity γi,m.Choosing(59)v(x)=-∇u¯(x)∇p(x),x∈Ω,
we ensure J(γi,m+ηv)<J(γim) if η>0 is small enough. We have defined v only in Ω, but formula (59) holds for all x∈ℛ. When γi is piecewise constant, we expect γi,m to be almost piecewise constant. In this case, another choice for the corrector v is a piecewise constant approximation:
(60)v(x)=vj=-Re[∫Ωj‍∇u¯∇pdz],x∈Ωj,j=1,…,d.The reconstruction algorithm for the impedance tomography problem is similar to the algorithm described in Section2.3 with γi, γe playing the role of ki, ke, and using (55)-(56) for the topological derivatives, (59) or (60) for the parameter correctors. We test this method in the next experiment where ℛ is the unit ball. Data are given at 30 uniformly distributed electrodes on ∂ℛ (represented by crosses in Figure 12) for the ten current patterns
(61)jℓ(t)=cos(ℓt),jℓ+5=sin(ℓt),ℓ=1,…,5,jℓ(t)=cos(ℓt),.jℓ+5=sin(ℓt),t∈[0,2π).
We consider complex admittivities of the form γ=σ+ıωɛ with ω=1. Inside the two objects represented by solid blue lines in Figures 12(a)-12(b), the unknown admittivity is γi=8+2ı and γe=1+ı is assumed to be known. We started the hybrid algorithm with γi,0=5+3ı. The first guess for the domains was obtained observing the topological derivative with Ω=∅ and γi=γi,0. It has only one object, as shown in Figure 12(a). We applied the algorithm alternating one iteration with respect to the domain with five iterations with respect to the parameter. Although we initialized the method with an incorrect number of objects, after 17 iterations with respect to the domain we obtained a reasonable reconstruction of the two defects; see Figure 12(b). When the algorithm stops, the true admittivity γi=8+2ı is approximated by 7.63+2.30ı. The values of the conductivity σi (the real part of γi) and of the permittivity ɛi (the imaginary part of γi) throughout the iterative procedure are represented in Figures 12(c) and 12(d), respectively. We refer to [23] for an example of the reconstruction of an object with unknown spatial dependent admittivity. The interested reader can also find in [23] a gallery of reconstructions mainly focused on the conductivity problem (with γ∈ℝ).Reconstruction of objects by impedance tomography. (a) Initial guess. (b) Final reconstruction of the objects after 17 updates with respect to the domain. (c) Convergence of the conductivity. (d) Convergence of the permittivity.
(a)
(b)
(c)
(d)


## 3.2.  Tomography with Discrete Electrodes

The impedance tomography problem with discrete electrodes (43)-(44) can be reformulated as an optimization problem: Minimize
(62)J(ℛ∖Ω¯,γi)=12∑k=1N‍∫ek‍|u+zkγe∂nu-umeas,k|2dl,
where u is a solution of
(63)∇·γe∇u=0inℛ∖Ω¯,∇·γi∇u=0inΩ,u--u+=0on∂Ω,γe∂nu--γi∂nu+=0on∂Ω,γe∂nu=0on∂ℛ∖∪j=1Ne¯j,∫ek‍γe∂nudl=Ikk=1,…,N.As before, we assume thatℛ contains inclusions Ω of a different material. For simplicity, we work in the two dimensional setting, but the method remains valid in 3D with straightforward modifications. The admittivity of the background medium γe is known. To determine its spatial variation inside the inclusions we must optimize (62) with respect to the interior domains and their admittivity. As shown before, even if the admittivity is unknown, we may locate the regions where it undergoes noticeable changes by computing the topological derivative of the shape functional in the whole domain. Once a first guess for Ω is found, we may approximate the admittivity γi optimizing with respect to γi for Ω fixed by a gradient technique. The guesses for Ω can be updated computing new topological derivatives. Below, we give expressions for the correctors to be used in the gradient method and for the topological derivatives.Let us fixΩ and seek to minimize J(γi)=J(ℛ∖Ω¯,γi) when u is a solution of the forward problem with discrete electrodes (63). When the nature of the inclusions we look for is known, an initial approximate value γi,0 is usually available. Otherwise, we select a first approximation by perturbing γe:γi,0=γe+ϵ. We may improve this first approximation iteratively by a gradient method. Given a guess γi,m, we generate a descent direction v to improve this guess studying the function of a real variable J(η)=J(ℛ∖Ω¯,γi,m+ηv). As in the acoustic setting and in the “continuous” model, η>0 and v are selected in such a way that J′(0)<0. An explicit expression for this derivative provides the following formula for the descent directions v.Theorem 7.
The derivative with respect toη of the function J(η)=J(ℛ∖Ω¯,γi,m+ηv), with J defined in (62) is
(64)dJ(η)dη|η=0=Re[∫Ω‍v∇u∇p¯dz],
where u is a solution of the forward problem (63) with parameter γi=γi,m in Ω and p is a solution of the adjoint problem (65) given below
(65)(p,P)∈H1(ℛ)⊕ℝ0N,(66)b((w,W),(p,P))=ℓa(w,W),∀(w,W)∈H1(ℛ)⊕ℝ0N,
where
(67)b(Ω;(w,W),(p,P)):=∫R∖Ω¯‍γe∇w∇p¯dz+∫Ω‍γi∇w∇p¯dz+∑k=1N‍1zk∫ek‍(w-Wk)(p¯-P¯k)dl,(68)ℓa(w,W):=-∑k=1N‍∫ek‍(u+zkγe∂nu-umeas,k¯)(w+zkγe∂nw)dl,
with γi=γi,m in Ω. Choosing
(69)v(x)=-Re[∇u¯(x)∇p(x)],x∈Ω
and η>0 small enough, one ensures J(γi,m+ηv)<J(γi,m).Proof.
The derivative with respect toη is
(70)ddηJ(γi,m+ηv)|η=0=Re[∑k=1N‍∫ek‍(u+zkγ∂nu-umeas,k¯)(u˙+zkγe∂nu˙)dl],
where u˙=(d/dη)uη|η=0 and uη is a solution of (63) with γi=γi,m+ηv; that is,
(71)(uη,Uη)∈H1(ℛ)⊕ℝ0Nb(η;(uη,Uη),(w,W))=ℓ(w,W),=ℓ(w,W)∀(w,W)∈H1(ℛ)⊕ℝ0N,
where
(72)b(η;(u,U),(w,W)):=∫ℛ∖Ω¯‍γe∇u∇w¯dz+∫Ω‍(γi,m+ηv)∇u∇w¯dz+∑k=1N‍1zk∫ek‍(u-Uk)(w¯-W¯k)dl,
and ℓ is defined in (50). Notice that when η=0 and γi=γi,m, the sesquilinear forms (67) and (72) are the same.
We eliminateu˙ by introducing the adjoint problem. Let us set
(73)ℒ(η;(u,U),(p,P)):=J(η)+Re[b(η;(u,U),(p,P))-ℓ(p,P)],
for (u,U),(p,P)∈H1(ℛ)⊕ℝ0N. Then, choosing (u,U)=(uη,Uη), it follows that
(74)J(η)=ℒ(η;(uη,Uη),(p,P)),∀(p,P)∈H1(ℛ)⊕ℝ0N
and therefore
(75)ddηJ(η)|η=0=ddηℒ(η;(uη,Uη),(p,P))|η=0=Re[ddηb(η;(u,U),(p,P))|η=0]+Re[b(0;(u˙,U˙),(p,P))]+Re[∑k=1N‍∫ek(u+zkγe∂nu-umeas,k¯)×(u˙+zkγ∂nu˙)dl{∑k=1N‍∫ek(u+zkγe∂nu-umeas,k¯)}],
for all (p,P)∈H1(ℛ)⊕ℝ0N. Selecting as (p,P) a solution of (65), we obtain the desired formula
(76)ddηJ(η)|η=0=Re[ddηb(η;(u,U),(p,P))]|η=0=Re∫Ω‍v∇u∇p¯dz.Notice that ifw∈H1(ℛ), then w∈L2(∂ℛ)∩H1/2(∂ℛ) and ∂w/∂xi∈H-1/2(∂ℛ). Thus, the integrals appearing in the definition (68) of ℓa are in fact duality products.We have definedv only in Ω, but formula (69) holds for all x∈ℛ. When γi is piecewise constant, we replace v defined in (69) with piecewise constant approximations inside each object Ωj,
(77)v(x)=vj=-Re[∫Ωj‍∇u¯∇pdz],x∈Ωj,j=1,…,d.
An initial approximation of the objects buried in the medium as well as their successive corrections is computed using the topological derivative. For the tomography problem with discrete electrodes, the explicit expression for the topological derivative in terms of forward and adjoint fields is given in the theorem below. Combining the expressions for the approximation of the admittivity γi with the topological derivative computation to approximate the domains Ω, we can define a hybrid reconstruction scheme for the discrete electrode model as in the previous sections. The implementation of the numerical codes is out of the scope of this paper and will be the subject of our future work.Theorem 8.
The topological derivative of the cost functionalJ(ℛ∖Ω¯;γi) defined in (62) is given by
(78)DT(x,ℛ∖Ω¯,γi)=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)∇u(x)∇p¯(x)],55=Re[2γe(x)(γi(x)-γe(x))γe(x)+γi(x)]x∈ℛ∖Ω¯,DT(x,ℛ∖Ω¯,γi)=Re[2γi(x)(γi(x)-γe(x))γi(x)+γe(x)∇u(x)∇p¯(x)],5=Re[2γi(x)(γi(x)-γe(x))γi(x)+γe(x)∇.(x)]x∈Ω,
where u and p solve the forward and adjoint problems (63) and (65) with coefficient γi in Ω.Proof.
We adapt the strategy used in [13, 21] for inverse scattering problems. The topological derivative of a shape functional can be computed as a limit of shape derivatives [13]. Consider the family of deformations depending on the real parameter τ,
(79)φτ(x):=x+τV(x),x∈ℝ2,
where the vector field V is defined as
(80)V(z)=-n(z),z∈∂Bɛ(x),
and extended to ℝ2 in such a way that it vanishes away from a narrow neighborhood of ∂Bɛ(x). Then,
(81)DT(x,ℛ)=limɛ→01h′(ɛ)ddτJ(φτ(ℛ∖Bɛ(x)¯))|τ=0,x∈ℛ,
where h′(ɛ) is the derivative of the function h(ɛ) that appears in the definition of the topological derivative (15). For the electrical impedance tomography problem in 2D, we can take h(ɛ)=πɛ2 as we will shortly see.
The proof is organized in three steps. We assume first thatγe and γi are constant and compute the topological derivative when Ω=∅. Afterwards, we adapt the proof to smoothly varying coefficients. Finally we consider the case Ω≠∅.
Step
1 (computation of the topological derivative whenΩ=∅andγeandγiare constant parameters). Let us compute first
(82)ddτJ(φτ(ℛ∖Bɛ(x)¯))|τ=0.
Since the field V decays fast enough away from ∂Bɛ(x), it follows that φτ(∂ℛ)=∂ℛ and φτ(ℛ∖Bɛ(x)¯)=ℛ∖Ω¯τ, where Ωτ:=φτ(Bɛ(x)). The cost functional in the transformed domains is then
(83)J(ℛ∖Ω¯τ)=12∑k=1N‍∫ek‍|uτ+zkγe∂nuτ-umeas,k|2dl,
where uτ solves
(84)(uτ,Uτ)∈H1(ℛ)⊕R0N,b(Ωτ;(uτ,Uτ),(w,W))=ℓ(w,W),∀(w,W)∈H1(ℛ)⊕R0N,
with b defined in (67) and ℓ defined in (50). Notice that when τ=0, Ωτ=Bɛ(x), and the solution of (84) is the solution (uɛ,Uɛ) to (63) with Ω=Bɛ(x).
The derivative of (83) with respect to τ involves the derivative u˙ of uτ with respect to τ. To avoid its computation we introduce the Lagrangian functional
(85)ℒ(Ωτ;(u,U),(p,P)):=J(ℛ∖Ω¯τ)+Re[b(Ωτ;(u,U),(p,P))-ℓ(p,P)].
Then, for all (p,P)∈H1(ℛ)⊕R0N,
(86)ddτJ(ℛ∖Ω¯τ)|τ=0=ddτℒ(Ωτ;(uτ,Uτ),(p,P))|τ=0=Re[ddτb(Ωτ;(uɛ,Uɛ),(p,P))|τ=0]+Re[b(Bɛ(x);(u˙,U˙),(p,P))]+Re[∑k=1N‍∫ek(uɛ+zkγe∂nuɛ-umeas,k¯)×(u˙+zkγe∂nu˙)dl{∑k=1N‍∫ek(uɛ+zkγe∂nuɛ-umeas,k¯)}].
If (pɛ,Pɛ) satisfies (65) with object Ω=Bɛ(x) and parameter γi, then
(87)ddτJ(ℛ∖Ω¯τ)|τ=0=Re[ddτb(Ωτ;(uɛ,Uɛ),(pɛ,Pɛ))|τ=0]=Re[∫ℛ∖Bɛ(x)¯‍γe(∇uɛ·∇p¯ɛ)∇·Vdz-∫ℛ∖Bɛ(x)¯‍γe((∇V+∇V⊤)∇uɛ)·∇p¯ɛdz+∫Bɛ(x)‍γi(∇uɛ·∇p¯ɛ)∇·Vdz-∫Bɛ(x)‍γi((∇V+∇V⊤)∇uɛ)·∇p¯ɛdz]=Re[∫∂Bɛ(x)‍(γi(γiγe-1)∂nuɛ-∂np¯ɛ-+(γi-γe)∂tuɛ-∂tp¯ɛ-{(γiγe-1)})dl{∫∂Bɛ(x)‍}].
To obtain the last equality in (87), we have integrated by parts and used that uɛ and pɛ solve (63) and (65) with Ω=Bɛ(x), respectively, and that V vanishes on ∂ℛ.
By (81), the topological derivative at x can be computed now as
(88)DT(x,ℛ)=limɛ→01h′(ɛ)Re[∫∂Bɛ(x)‍(γi(γiγe-1)∂nuɛ-∂np¯ɛ-+(γi-γe)∂tuɛ-∂tp¯ɛ-{(γiγe-1)})dl{∫∂Bɛ(x)‍}].
The asymptotic behavior of uɛ and pɛ at ∂Bɛ(x) is obtained expressing them as corrections of u0 and p0, the solutions to (63) and (65) with Ω=∅,
(89)uɛ(z)=u0(z)χℛ∖Bɛ(x)¯(z)+vɛ(z),pɛ(z)=p0(z)χℛ∖Bɛ(x)¯(z)+qɛ(z),
and expanding the remainders in powers of ɛ. Reasoning as in [23] we find that
(90)uɛ(z)⟶u0(x),∇uɛ-(z)⟶2γeγe+γi∇u0(x),uɛ(z)⟶u0(x),55.∇uɛ-(z)⟶asɛ⟶0,
uniformly when |z-x|=ɛ. For pɛ there is a slight technical difference because pɛ is defined by the variational problem (65). However, this is due to the source term at ∂ℛ, that is included in p0. Once that source term is eliminated, the remainder qɛ(z) satisfies equations similar to the equations for vɛ and we find
(91)pɛ(z)⟶p0(x),∇pɛ-(z)⟶2γeγe+γi∇p0(x),pɛ(z)⟶p0(x),∇pɛ-(z)5.0⟶asɛ⟶0,
uniformly when |z-x|=ɛ.
Replacing the limits (90) and (91) in (88) and choosing h(ɛ)=πɛ2, formula (78) with Ω=∅ follows.
Step
2 (variable parameters). The previous derivation holds when the parameters are constant. When they vary smoothly, we expand them about x and repeat the previous computations. Notice that due to the choice of V, only the values of γi and γe in a narrow band about ∂Bɛ(x) are relevant for the computation of the shape derivative. For the topological derivative, we set γi(z)=γi(x)+(z-x)·∇γi(x)+O(ɛ2) and an analogous expression for γe. Then, we check that the limits (90)-(91) also hold with γi=γi(x) and γe=γe(x). For continuous γe and γi, we use a regularizing approximation, see [21] for details.
Step
3 (topological derivative with inclusionsΩ). We have to redo the previous steps taking as objects Ω∖Bɛ(x)¯ for x∈Ω, and Ω∪Bɛ(x) when x∈ℛ∖Ω¯ depending on whether we are computing the derivative inside or outside Ω; see [21] for details.



---

*Source: 816134-2013-09-02.xml*