Hybrid Topological Derivative-Gradient Based Methods for Nondestructive Testing

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.

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 Figure 1. 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, measurements 𝑢_meas 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 𝑘_𝑖 minimizing 𝐽ℛ⧵Ω,𝑘_𝑖=12_Γmeas𝑢−𝑢_meas²𝑑𝑙, (1) where 𝑢 is the solution of the forward problem with objects Ω immersed in the exterior media ℛ⧵Ω and parameters 𝑘_𝑖 (𝑖 stands for “interior”). The design variables are Ω and 𝑘_𝑖. 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. Section 2 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, 𝑈_inc(𝐱,𝑡)=Re𝑒^{−𝑖𝜔𝑡}𝑢_inc(𝐱), (2) where 𝜔>0 is the frequency. In this case, after a sufficiently long time, the total wave field is also time harmonic, 𝑈(𝐱,𝑡)=Re[𝑒^{−𝑖𝜔𝑡}𝑢(𝐱)], and the total amplitude 𝑢(𝐱) 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 direction 𝐝, 𝑢_inc(𝐱)=𝑒^{𝑖𝑘𝑒𝐱⋅𝐝} with wave number 𝑘_𝑒>0. Other types of incident waves can be considered, for instance, in ℝ², 𝑢_inc(𝐱)=(𝚤/4)𝐻₀⁽¹⁾(𝑘_𝑒|𝐱−𝐱₀|) models a point source excitation at 𝐱₀. Here, 𝐻₀⁽¹⁾ 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 Figure 1). We have a finite number of bounded and simply connected open objects Ω¹,…,Ω^𝑑 having no pairwise contact. Their boundaries are assumed to be smooth, for instance 𝒞². For simplicity, we will use the notation Ω=∪^𝑑_𝑗=1Ω^𝑗. These objects are immersed in the exterior media ℝ²⧵Ω. 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 ℝ²⧵Ω is homogeneous. Then, the total wave field 𝑢=𝑢_inc+𝑢_scinℝ²⧵Ω,𝑢=𝑢_trinΩ, (3) where 𝑢_sc and 𝑢_tr denote the scattered wave outside the obstacles ℝ²⧵Ω and the transmitted wave 𝑢_tr inside Ω, respectively, is the solution of Δ𝑢+𝑘²_𝑒𝑢=0inℝ²⧵ΩΔ𝑢+𝑘²_𝑖𝑢𝑢=0inΩ₋=𝑢⁺,𝜕_𝐧𝑢⁻=𝜕_𝐧𝑢⁺on𝜕Ωlim_𝑟→∞𝑟^1/2𝜕_𝑟𝑢−𝑢_inc−𝑖𝑘_𝑒𝑢−𝑢_inc=0. (4) Here, 𝑘_𝑒 is the wave number of the incident wave and 𝑘_𝑖(𝐱)>𝜅_𝑖>0 is related to the material properties of the inclusions Ω. We assume that 𝑘_𝑖 is differentiable. When 𝑘_𝑖 is constant inside each Ω^𝑗, we are dealing with homogeneous materials.

The condition at infinity is the standard Sommerfeld radiation condition with 𝑟=|𝐱|, which has to be satisfied uniformly in all directions. It implies that only outgoing waves are allowed. 𝜕_𝑟 denotes radial derivatives. 𝜕_𝐧 stands for outward normal derivatives, with 𝐧 pointing outside Ω. 𝑢⁺ and 𝑢⁻ denote the limits of 𝑢 from the exterior and interior of Ω, respectively. The continuity condition for the normal derivative can be extended to 𝛼_𝑒𝜕_𝐧𝑢⁻=𝛼_𝑖𝜕_𝐧𝑢⁺ when Δ𝑢 is replaced by div(𝛼_𝑒∇𝑢) in ℝ²⧵Ω and by div(𝛼_𝑖∇𝑢) in Ω, as explored in [21].

The optimization problem we want to solve is then as follows: Find Ω and 𝑘_𝑖(𝐱) minimizing 𝐽ℝ²⧵Ω,𝑘_𝑖=12_Γmeas𝑢(𝐱)−𝑢_meas(𝐱)²𝑑𝑙, (5) where 𝑢(𝐱) is the solution of the forward problem (4) with objects Ω and parameters 𝑘_𝑖. 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 parameters 𝑘_𝑖,0 and the objects Ω₀ and then generate sequences of approximations 𝑘_𝑖,𝑚, Ω_𝑚 along which the functional 𝐽 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 𝑘_𝑖. Next, we will assume that the parameters 𝑘_𝑖 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 𝑘_𝑖(𝐱) minimizing 𝐽𝑘_𝑖=12_Γmeas𝑢(𝐱)−𝑢_meas(𝐱)²𝑑𝑙, (6) where 𝑢(𝐱) is the solution of the forward problem (4) with objects Ω and parameters 𝑘_𝑖.

Fixed an initial guess 𝑘_𝑖,0, we generate a sequence 𝑘_𝑖,𝑚 providing a decreasing sequence 𝐽(𝑘_𝑖,𝑚) by a gradient method [21]. At each stage, we fix the current guess 𝑘_𝑖,𝑚 and study the variations of the functional along directions 𝑣: 𝐽(𝜂)=𝐽(𝑘_𝑖,𝑚+𝜂𝑣), where 𝜂>0 is a real parameter. The derivative with respect to 𝜂 is given by 𝑑𝐽𝜂=𝑘𝑑𝜂𝑑𝐽_𝑖,𝑚+𝜂𝑣𝑑𝜂=Re_Γmeas𝑢_𝜂(𝐱)−𝑢_meas(𝐱)𝑑𝑢_𝜂(𝐱),𝑑𝜂𝑑𝑙 (7)𝑢_𝜂 being the solution of the forward problem (4) with coefficient 𝑘_𝑖=𝑘_𝑖,𝑚+𝜂𝑣 inside the objects Ω. Computing 𝑑𝑢_𝜂/𝑑𝜂 is avoided introducing an adjoint field 𝑝; see [21].

Theorem 1. The derivative with respect to 𝜂 of the function 𝐽(𝜂)=𝐽(𝑘_𝑖,𝑚+𝜂𝑣) defined by (6) evaluated at zero is given by 𝑑𝐽𝜂𝑑𝜂_𝜂=0=2Re_Ω𝑣𝑢𝑝, (8) where 𝑝 is the solution of the adjoint problem Δ𝑝+𝑘²_𝑒𝑝=𝑢_{𝑚𝑒𝑎𝑠}−𝑢𝛿_{Γ𝑚𝑒𝑎𝑠}𝑖𝑛ℝ²⧵ΩΔ𝑝+𝑘²_𝑖𝑝𝑝=0𝑖𝑛Ω₋=𝑝⁺,𝜕_𝐧𝑝⁻=𝜕_𝐧𝑝⁺𝑜𝑛𝜕Ωlim_𝑟→∞𝑟^1/2𝜕_𝑟𝑝−𝑖𝑘_𝑒𝑝=0 (9) and 𝑢 the solution of the forward problem (4) setting in both cases 𝑘_𝑖=𝑘_𝑖,𝑚.

Choosing 𝑣 in such a way that the derivative 𝐽^′(0)<0, we ensure that 𝐽(𝑘_𝑖,𝑚+𝜂𝑣)<𝐽(𝑘_𝑖,𝑚) for 𝜂>0 small enough. For instance, we may set 𝑣(𝐱)=−Re𝑢(𝐱)𝑝(𝐱). (10) The updated guess is then 𝑘_𝑖,𝑚+1(𝐱)=𝑘_𝑖,𝑚(𝐱)−𝜂Re𝑢(𝐱)𝑝(𝐱). (11) This choice is reasonable if we are looking for space-dependent parameters. When Ω is composed of several disjoint domains Ω^𝑗, and 𝑘_𝑖 is known to be almost piecewise constant, we may set 𝑣(𝐱)=𝑣_𝑗=−Re_Ω𝑗𝑢𝑝,𝐱∈Ω^𝑗. (12) Technical details and proofs can be found in [21], where more general Helmholtz type constraints are considered: div𝛼_𝑒(𝐱)∇𝑢+𝑘_𝑒(𝐱)²𝑢=0inℝ²⧵𝛼Ω,div_𝑖(𝐱)∇𝑢+𝑘_𝑖(𝐱)²𝑢𝑢=0inΩ,₋=𝑢⁺,𝛼_𝑒(𝐱)𝜕_𝐧𝑢⁻=𝛼_𝑖(𝐱)𝜕_𝐧𝑢⁺on𝜕Ω,lim_𝑟→∞𝑟^1/2𝜕_𝑟𝑢−𝑢_inc−𝑖𝜅_𝑒𝑢−𝑢_inc=0. (13) The coefficients 𝑘_𝑒, 𝑘_𝑖, 𝛼_𝑒, 𝛼_𝑖 are bounded from below by positive constants and depend smoothly on 𝐱, becoming constant at infinity. The constant 𝜅_𝑒 at the Sommerfeld radiation condition is defined as 𝜅_𝑒=lim_|𝐱|→∞𝑘_𝑒(𝐱)/𝛼_𝑒(𝐱). It is the wave number of the incident wave. The interior parameter 𝛼_𝑖(𝐱) 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 parameters 𝑘_𝑖 fixing the objects Ω. Let us assume now that we know the parameters 𝑘_𝑖. The optimization problem becomes as follows: Find Ω minimizing 𝐽ℝ²⧵Ω=12_Γmeas𝑢(𝐱)−𝑢_meas(𝐱)²𝑑𝑙, (14) where 𝑢(𝐱) is the solution of the forward problem (4) with objects Ω and parameters 𝑘_𝑖. 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 functional 𝐽(ℛ) at the point 𝐱∈ℛ is defined as [42] 𝐷_𝑇(𝐱,ℛ)=lim_𝜀→0𝐽ℛ⧵𝐵_𝜀(𝐱)−𝐽(ℛ)ℎ(𝜀), (15) where 𝐵_𝜀(𝐱) is a ball centered at 𝐱 with small radius 𝜀>0; see Figure 2. The scalar function ℎ(𝜀)>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 ℎ(𝜀)=Vol(𝐵_𝜀(𝐱))=𝜋𝜀². For Dirichlet conditions in 2D, we may choose ℎ(𝜀)=−2𝜋/log(𝑘_𝑒𝜀).

Figure 2 

Geometry for the computation of topological derivatives.

The topological derivative is a scalar function defined for every 𝐱∈ℛ which measures the sensitivity of the functional to removing from ℛ a small object located at 𝐱. In regions where 𝐷_𝑇(𝐱,ℛ) takes large negative values, we have a large probability to find an object. This remark is based on the following expansion: 𝐽ℛ⧵𝐵_𝜀(𝐱)=𝐽(ℛ)+ℎ(𝜀)𝐷_𝑇(𝐱,ℛ)+𝑜(ℎ(𝜀)),𝜀⟶0. (16) Whenever 𝐷_𝑇(𝐱,ℛ)<0, 𝐽(ℛ⧵𝐵_𝜀(𝐱))<𝐽(ℛ) for 𝜀>0 small. When 𝐽 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 Ω₀ for the objects as Ω₀∶=𝐱∈ℛ∣𝐷_𝑇(𝐱,ℛ)≤−𝐶₀, (17) where 𝐶₀ is a positive constant. (ii) Check that 𝐽(ℛ⧵Ω₀)<𝐽(ℛ). Otherwise increase the constant 𝐶₀. (iii) For 𝑚=1∶𝑚max(a)  Compute the topological derivative when Ω=Ω_𝑚−1 and select Ω_𝑚=Ω_𝑚−1∪𝐱∈ℛ⧵Ω_𝑚−1∣𝐷_𝑇𝐱,ℛ⧵Ω_𝑚−1≤−𝐶_𝑚 (18) for decreasing thresholds 𝐶_𝑚>0. (b) Check that 𝐽(ℛ⧵Ω_𝑚)<𝐽(ℛ⧵Ω_𝑚−1). Otherwise increase the constant 𝐶_𝑚. (c) Check the stopping criteria: the iteration stops if either meas (Ω_𝑚⧵Ω_𝑚−1) is small enough, or 𝐽(ℛ⧵Ω_𝑚) is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle |𝑢^𝜀_meas−𝑢_𝑚|<𝜏𝜀 is satisfied, where 𝜀 describes measurement errors and 𝜏>1 (in our numerical experiments 𝜏=1.2). Here, 𝑢_𝑚 is the solution of the forward problem when the objects are Ω_𝑚.

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 of 𝐽(ℝ²⧵Ω) defined by (14) is given by 𝐷_𝑇𝐱,ℝ²⧵Ω=Re𝑘²_𝑖(𝐱)−𝑘²_𝑒𝑢(𝐱)𝑝(𝐱),𝐱∈ℝ²⧵Ω, (19) where 𝑢 and 𝑝 are solutions of the forward and adjoint problems (4) and (9), respectively.

Notice that when Ω and 𝑘_𝑖 are the true objects and parameters for which the data are measured, the adjoint field 𝑝 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 directions 𝐝^𝑗, 𝑗=1,…,𝑀. The shape cost functional (14) is then replaced with: 𝐽ℝ²⧵Ω=12^𝑀_𝑗=1Γmeas𝑢𝐱;𝐝^𝑗−𝑢_meas𝐱;𝐝^𝑗2𝑑𝑙_𝐱. (20) Here, 𝑢(𝐱;𝐝^𝑗) is the solution to the forward problem (4) with incident field 𝑢_inc(𝐱;𝐝^𝑗)=𝑒^{𝚤𝑘𝑒𝐱⋅𝐝𝑗} and 𝑢_meas(𝐱;𝐝^𝑗) 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 is 𝑘_𝑖=1 whereas 𝑘_𝑒=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 Ω₀ is represented in Figure 4(b). We compute now the topological derivative when Ω=Ω₀, 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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 4 

Reconstruction of objects when 𝑘𝑖=1 is known using the monotone algorithm. (a) Topological derivative when Ω=∅. (b)–(j) Approximate domains Ω𝑚 superimposed to the topological derivative for Ω=Ω𝑚, for 𝑚=0,…,8. (k) Cost functional versus the number of iterations.

The quality of the reconstructions as well as the number of iterations depends on the values of the constants 𝐶_𝑚 appearing in (18). Guidelines for their selection are given in [21]. When 𝐶_𝑚 are too big, the number of iterations required is high. A sequence of smaller 𝐶_𝑚 provides faster reconstructions, but we are prone to include spurious regions in Ω_𝑚 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 𝐶₀, our initial guess Ω₀, 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 Ω=Ω₀, 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.

(a)

(b)

(c)

(d)

(e)

Figure 5 

Reconstruction of objects when 𝑘𝑖=1 is known, with the monotone algorithm (18). (a) Topological derivative when Ω=∅. (b)–(d) Approximate domains Ω𝑚 superimposed to the topological derivative for Ω=Ω𝑚. (e) Evolution of the cost functional.

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 Figures 5(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 Ω_𝑚, 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 𝐽(ℛ⧵Ω) was extended to the points inside Ω in [21] as 𝐷_𝑇𝐱,ℛ⧵Ω=lim_𝜀→0𝐽ℛ⧵Ω⧵𝐵_𝜀(𝐱)−𝐽ℛ⧵Ωℎ(𝜀),𝐱∈Ω. (21) From this definition, the expansion below follows: 𝐽ℛ⧵Ω⧵𝐵_𝜀(𝐱)=𝐽ℛ⧵Ω+ℎ(𝜀)𝐷_𝑇𝐱,ℛ⧵Ω+𝑜(ℎ(𝜀)),𝜀⟶0. (22) Since ℎ(𝜀) is a positive function, if 𝐷_𝑇(𝐱,ℛ⧵Ω) is large and negative, the functional decreases when we remove the ball 𝐵_𝜀(𝐱) from Ω. This suggests a strategy to remove points from Ω. Our algorithm can generate a nonmonotone sequence of objects by substituting the definition of Ω_𝑚 in the previous method (see (18)) with Ω_𝑚=Ω_𝑚−1∪𝐱∈ℛ⧵Ω_𝑚∣𝐷_𝑇𝐱,ℛ⧵Ω_𝑚−1≤−𝐶_𝑚⧵𝐱∈Ω_𝑚−1∣𝐷_𝑇𝐱,ℛ⧵Ω_𝑚−1≤−𝑐_𝑚, (23) with decreasing sequences 𝐶_𝑚, 𝑐_𝑚>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 𝐽(ℝ²⧵Ω) defined by (14) is given by 𝐷_𝑇𝐱,ℝ²⧵Ω=Re𝑘²_𝑒−𝑘²_𝑖(𝐱)𝑢(𝐱)𝑝(𝐱),𝐱∈Ω, (24) where 𝑢 and 𝑝 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 ℝ², the expression Re𝑘²_𝑖(𝐱)−𝑘²_𝑒𝑢(𝐱)𝑝(𝐱) (25) defines a continuous function in ℝ². This motivated in [21] to modify the extended definition of the topological derivative inside Ω of a functional 𝐽(ℛ⧵Ω) as 𝐷_𝑇𝐱,ℛ⧵Ω=lim_𝜀→0𝐽ℛ⧵Ω−𝐽ℛ⧵Ω⧵𝐵_𝜀(𝐱)ℎ(𝜀),𝐱∈Ω, (26) which is exactly the same as (21) multiplied by −1. Using this notion, 𝐽ℛ⧵Ω⧵𝐵_𝜀(𝐱)=𝐽ℛ⧵Ω−ℎ(𝜀)𝐷_𝑇𝐱,ℛ⧵Ω+𝑜(ℎ(𝜀)),𝜀⟶0, (27) and whenever 𝐷_𝑇(𝐱,ℛ⧵Ω) is large and positive, the functional decreases when removing the ball 𝐵_𝜀(𝐱) from Ω. This definition also allows to define a non-monotone scheme by substituting the definition of Ω_𝑚 in (23) by Ω_𝑚=Ω_𝑚−1∪𝐱∈ℛ⧵Ω_𝑚∣𝐷_𝑇𝐱,ℛ⧵Ω_𝑚−1≤−𝐶_𝑚⧵𝐱∈Ω_𝑚−1∣𝐷_𝑇𝐱,ℛ⧵Ω_𝑚−1≥𝑐_𝑚, (28) for decreasing thresholds 𝐶_𝑚, 𝑐_𝑚>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 𝐽(ℝ²⧵Ω) defined by (14) is given by 𝐷_𝑇𝐱,ℝ²⧵Ω=Re𝑘²_𝑖(𝐱)−𝑘²_𝑒𝑢(𝐱)𝑝(𝐱),𝐱∈Ω, (29) where 𝑢 and 𝑝 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 Figure 5. The topological derivative for Ω=∅ provided the initial guess Ω₀ 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 Ω₀) 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 𝐶_𝑚) after nine iterations. Our non-monotone method is not very sensitive to the initial guess.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

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 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 𝑘_𝑒, 𝑘_𝑖 (simulating high and low frequencies), problems with different values of 𝑘_𝑖 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 guess 𝑘_𝑖,0 for the parameters 𝑘_𝑖, for example, a perturbation of the background parameter 𝑘_𝑒. (b) Choose an initial guess Ω₀ for the objects Ω, using the topological derivative: set 𝑘_𝑖=𝑘_𝑖,0 and compute the topological derivative with object Ω=∅: 𝐷_𝑇𝐱,ℝ²=Re𝑘²_𝑖,0−𝑘²_𝑒𝑢(𝐱)𝑝(𝐱), (30) where 𝑢 and 𝑝 are forward and adjoint fields with object Ω=∅ and parameter 𝑘_𝑖,0. Then, Ω₀=𝐱∈ℝ²∣𝐷_𝑇𝐱,ℝ²,𝑘_𝑖,0≤−𝐶₀, (31) where 𝐶₀ is a certain large positive threshold. (ii) Iteration(a) Update 𝑘_𝑖,𝑚(𝐱) using the gradient technique: 𝑘_𝑖,𝑚+1(𝐱)=𝑘_𝑖,𝑚(𝐱)−𝜂Re𝑢(𝐱)𝑝(𝐱), (32) where 𝑢, 𝑝 solve forward and adjoint problems with objects Ω_𝑚 and parameter 𝑘_𝑖,𝑚(𝐱). (b) Check that 𝐽(ℝ²⧵Ω_𝑚,𝑘_𝑖,𝑚+1)<𝐽(ℝ²⧵Ω_𝑚,𝑘_𝑖,𝑚). Otherwise, reduce 𝜂.(c) Update Ω_𝑚 computing the topological derivative with objects Ω_𝑚 and parameter 𝑘_𝑖,𝑚+1(𝐱). We may resort to a monotone Ω_𝑚+1=Ω_𝑚∪𝐱∈ℝ²⧵Ω_𝑚∣𝐷_𝑇𝐱,ℝ²⧵Ω_𝑚,𝑘_𝑖,𝑚+1≤−𝐶_𝑚+1, (33) or nonmonotone strategy Ω_𝑚+1=Ω_𝑚∪𝐱∈ℝ²⧵Ω_𝑚∣𝐷_𝑇𝐱,ℝ²⧵Ω_𝑚,𝑘_𝑖,𝑚+1≤−𝐶_𝑚+1⧵𝐱∈Ω_𝑚∣𝐷_𝑇𝐱,ℝ²⧵Ω_𝑚,𝑘_𝑖,𝑚+1≥𝑐_𝑚+1, (34) with decreasing thresholds 𝐶_𝑚+1, 𝑐_𝑚+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 that 𝐽(ℝ²⧵Ω_𝑚+1,𝑘_𝑖,𝑚+1)<𝐽(ℝ²⧵Ω_𝑚,𝑘_𝑖,𝑚+1). Otherwise, increase the thresholds 𝑐_𝑚+1, 𝐶_𝑚+1. (e) Stop when meas (Ω_𝑚⧵Ω_𝑚−1) is small enough, or 𝐽(ℝ²⧵Ω_𝑚) is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle |𝑢^𝜀_meas−𝑢|<𝜏𝜀 is satisfied.

The iterative procedure described perviously generates sequences of domains Ω_𝑚 and parameters 𝑘_𝑖,𝑚 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 when 𝑘_𝑒 is constant and we are computing the first guess Ω₀ on the whole space. They become Δ𝑢+𝑘²_𝑒𝑢=0inℝ²lim_𝑟→∞𝑟^1/2𝜕_𝑟𝑢−𝑢_inc−𝑖𝑘_𝑒𝑢−𝑢_inc=0,Δ𝑝+𝑘²_𝑒𝑝=𝑢_meas−𝑢𝛿_Γmeasinℝ²lim_𝑟→∞𝑟^1/2𝜕_𝑟𝑝−𝑖𝑘_𝑒𝑝=0. (35)(36) 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 𝑘_𝑖 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 parameters 𝑘_𝑖,𝑚 several times for a fixed approximate domain Ω_𝑚 before proceeding to update Ω_𝑚, 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 Figure 3). We also take the same values of the parameters: 𝑘_𝑒=4 is assumed to be known and 𝑘_𝑖=1 is assumed to be unknown. We computed the topological derivative when Ω=∅ choosing 𝑘_𝑖,0=3 as initial guess for 𝑘_𝑖. Notice that since the forward and adjoint fields when Ω=∅ do not depend on 𝑘_𝑖, the topological derivative is the same as in Figure 4(a) except for the scale (we have now in (19) the multiplying factor ((𝑘_𝑖,0)²−𝑘²_𝑒) instead of (𝑘²_𝑖−𝑘²_𝑒)). Therefore, our initial guess Ω₀, represented in Figure 7(a), is the same as in Figure 4(b). We fix now Ω=Ω₀ and iterate five times to update the value of 𝑘_𝑖 without modifying the scatterers. Then, we fix the value of 𝑘_𝑖 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 𝑘_𝑖 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 (𝑘_𝑖 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].

(a)

(b)

(c)

Figure 7 

Reconstruction of objects with unknown 𝑘𝑖. (a) First guess of the objects. (b) Final reconstruction of the objects after nine updates of the domain (c) Values of 𝑘𝑖 through the iterative method.

In the next example we test the hybrid algorithm for the reconstruction of an object Ω with space-dependent coefficient 𝑘_𝑖. The scatterer is the ball centered at (−0.25,0.3) and radius 0.35 and 𝑘_𝑖 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 𝑘_𝑖,0=3. The topological derivative for Ω=∅ provided the initial guess Ω₀ 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 𝑘_𝑖 are shown in Figures 8(b) and 8(d), respectively. The parameter 𝑘_𝑖 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.

(a)

(b)

(c)

(d)

Figure 8 

Reconstruction an object with space-dependent parameter 𝑘𝑖. (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 𝑘𝑖. (d) Final reconstruction of 𝑘𝑖.

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 Δ𝑢+𝑘²_𝑒𝑢=0,inℝ²⧵Ω𝜕_𝐧𝑢=0,on𝜕Ωlim_𝑟→∞𝑟^1/2𝜕_𝑟𝑢−𝑢_inc−𝑖𝑘_𝑒𝑢−𝑢_inc=0. (37) 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: Δ𝑢+𝑘²_𝑒𝑢=0,inℝ²⧵Ω𝑢=0,on𝜕Ωlim_𝑟→∞𝑟^1/2𝜕_𝑟𝑢−𝑢_inc−𝑖𝑘_𝑒𝑢−𝑢_inc=0. (38) 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 𝐷_𝑇𝐱,ℝ²⧵Ω=Re2∇𝑢(𝐱)∇𝑝(𝐱)−𝑘²_𝑒𝑢(𝐱)𝑝(𝐱), (39) where the forward field is a solution of (37) and the adjoint field is a solution of Δ𝑝+𝑘²_𝑒𝑝=𝑢_meas−𝑢𝛿_Γmeas,inℝ²⧵𝜕Ω,_𝐧𝑝=0,on𝜕Ω,lim_𝑟→∞𝑟^1/2𝜕_𝑟𝑝−𝑖𝑘_𝑒𝑝=0. (40) 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]) 𝐷_𝑇𝐱,ℝ²⧵Ω=Re𝑢(𝐱)𝑝(𝐱), (41) where the forward field is a solution of (38) and the adjoint field is a solution of Δ𝑝+𝑘²_𝑒𝑝=𝑢_meas−𝑢𝛿_Γmeas,inℝ²⧵Ω,𝑝=0,on𝜕Ω,lim_𝑟→∞𝑟^1/2𝜕_𝑟𝑝−𝑖𝑘_𝑒𝑝=0. (42)

We test the performance of our iterative topological derivative method for rigid and sound-soft objects in Figures 9 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 𝑘_𝑒=4 too.

(a)

(b)

(c)

Figure 9 

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)

(d)

(e)

Figure 10 

Reconstruction of sound-soft objects. (a) Topological derivative when Ω=∅. (b) Initial guess Ω0 superimposed to the topological derivative for Ω=Ω0. (c)-(d) Approximate domains Ω𝑚 superimposed to the topological derivative for Ω=Ω𝑚, for 𝑚=1 and 𝑚=2. (e) Final reconstruction of the objects after 13 steps (dashed red line).

In Figure 9, 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 Figure 10(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 Ω₀, represented in Figure 10(b). In the same plot we also represent the values of the topological derivative when Ω=Ω₀. 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 Section 2. 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 𝑢(𝐱) satisfies ∇⋅𝛾(𝐱,𝜔)∇𝑢=0, (43) 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.

Currents 𝑗 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: 𝛾𝜕_𝐧𝑢=𝐼,on𝜕ℛ, (44) where 𝐼 satisfies 𝐼=0,on𝜕ℛ⧵∪^𝑁_𝑘=1𝑒_𝑘,_𝑒𝑘𝐼𝑑𝑙=𝐼_𝑘,𝑘=1,…,𝑁,_𝑁𝑘=1𝐼_𝑘=0, (45)𝐼_𝑘 is the current sent to the 𝑘th electrode and 𝑒_𝑘 denotes the part of 𝜕ℛ that corresponds to the 𝑘th electrode. We assume that each 𝑒_𝑘 is an open subset of 𝜕ℛ with positive surface measure and that dist(𝑒_𝑘,𝑒_𝑗)>0 for 𝑘≠𝑗. ℛ and Ω have smooth boundaries. 𝛾_𝑒 is continuous up to the boundary, 𝛾_𝑖 is continuous, and 𝑈_meas and 𝐼 define 𝐻^1/2(𝜕ℛ) functions.

At the electrodes, we measure the voltages 𝑢_meas,𝑘. Dirichlet conditions 𝑢=𝑢_meas,𝑘attheelectrodes𝑒_𝑘,𝑘=1,…,𝑁, (46) have been shown not to be realistic [29] and are commonly replaced with the constraints 𝑢+𝑧_𝑘𝛾𝜕_𝐧𝑢=𝑢_meas,𝑘,attheelectrodes𝑒_𝑘,𝑘=1,…,𝑁, (47) where 𝑧_𝑘 is the effective contact impedance or surface impedance and ^𝑁_𝑘=1𝑢_meas,𝑘=0. (48) As shown in [50], problem (43)-(44) admits the variational formulation: Find (𝑢,𝐔)∈𝐻¹(ℛ)⊕ℝ^𝑁₀ such that 𝑏((𝑢,𝐔),(𝑤,𝐖))=ℓ(𝑤,𝐖)∀(𝑤,𝐖)∈𝐻¹(ℛ)⊕ℝ^𝑁₀, (49) where 𝑏((𝑢,𝐔),(𝑤,𝐖))∶=_ℛ𝛾∇𝑢∇𝑤𝑑𝐳+^𝑁_𝑘=11𝑧_𝑘𝑒𝑘𝑢−𝑈_𝑘𝑤−𝑊_𝑘ℓ()𝑑𝑙,𝑤,𝐖∶=^𝑁_𝑘=1𝐼_𝑘𝑊_𝑘. (50) The subscript 0 refers to vectors of vanishing mean; that is, ∑^𝑁_𝑘=1𝑈_𝑘=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: 𝛾𝜕_𝐧𝑢=𝑗,on𝜕ℛ, (51) and measures the voltage on the boundary 𝑢=𝑢_meas,on𝜕ℛ. (52)

Here, 𝐧 denotes the outward unit normal and 𝑗 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 ∫_𝜕ℛ𝑗𝑑𝑙=0, which turns to be the law of conservation of charge. To select a unique potential 𝑢, we impose ∫_𝜕ℛ𝑢𝑑𝑙=0. The measured voltage also satisfies the condition ∫_𝜕ℛ𝑢_meas𝑑𝑙=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 Ω^𝑗, the admittivity 𝛾 is a piecewise smooth function in ℛ with discontinuities at the boundaries of the inclusions. The admittivity outside the objects 𝛾_𝑒 is known. The admittivity inside 𝛾_𝑖 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 Ω^𝑗 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 ℝ² with smooth boundary. We set Ω=∪^𝑑_𝑗=1Ω^𝑗 with Ω^𝑗 simply connected open bounded sets with smooth boundaries without pairwise contact; see Figure 11. The admittivity 𝛾_𝑒 in the background medium ℛ⧵Ω is known and 𝛾_𝑖 is unknown inside Ω. The functional to be optimized becomes 𝐽ℛ⧵Ω,𝛾_𝑖=12_𝜕ℛ𝑢−𝑢_meas²𝑑𝑙, (53) where 𝑢 is a solution of ∇⋅𝛾_𝑒∇𝑢=0inℛ⧵Ω,∇⋅𝛾_𝑖𝑢∇𝑢=0inΩ,₋−𝑢⁺𝛾=0on𝜕Ω,_𝑖𝜕_𝐧𝑢⁻−𝛾_𝑒𝜕_𝐧𝑢⁺𝛾=0on𝜕Ω,_𝑒𝜕_𝐧𝑢=𝑗on𝜕ℛ. (54) The unit normal 𝐧 points outside ℛ⧵Ω but inside Ω. 𝑢⁻ and 𝑢⁺ denote the limit values of 𝑢 on 𝜕Ω from outside and inside Ω, respectively; see Figure 11.

Figure 11 

Geometry for the impedance tomography problem.

As observed in Section 2, 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 𝐷_𝑇𝐱,ℛ⧵Ω,𝛾_𝑖=Re2𝛾_𝑒(𝐱)𝛾_𝑖(𝐱)−𝛾_𝑒(𝐱)𝛾_𝑒(𝐱)+𝛾_𝑖(𝐱)∇𝑢(𝐱)∇𝑝(𝐱),𝐱∈ℛ⧵𝐷Ω,_𝑇𝐱,ℛ⧵Ω,𝛾_𝑖=Re2𝛾_𝑖(𝐱)𝛾_𝑖(𝐱)−𝛾_𝑒(𝐱)𝛾_𝑖(𝐱)+𝛾_𝑒(𝐱)∇𝑢(𝐱)∇𝑝(𝐱),𝐱∈Ω, (55)(56) where 𝑢 is the solution of the forward problem (54) and 𝑝 is a solution of the adjoint problem ∇⋅𝛾_𝑒∇𝑝=0𝑖𝑛ℛ⧵Ω,∇⋅𝛾_𝑖𝑝∇𝑝=0𝑖𝑛Ω,₋−𝑝⁺=0𝑜𝑛𝜕Ω,𝛾_𝑖𝜕_𝐧𝑝⁻−𝛾_𝑒𝜕_𝐧𝑝⁺=0𝑜𝑛𝜕Ω,𝛾_𝑒𝜕_𝐧𝑝=𝑢_meas−𝑢𝑜𝑛𝜕ℛ. (57)

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 𝛾_𝑖,𝑚 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 𝐽(𝜂)∶=𝐽(ℛ⧵Ω,𝛾_𝑖,𝑚+𝜂𝑣) with 𝐽 defined in (53) is 𝑑𝐽𝛾_𝑖,𝑚+𝜂𝑣𝑑𝜂_𝜂=0=Re_Ω𝑣∇𝑢∇𝑝𝑑𝐳, (58) where 𝑢 and 𝑝 are solutions of the forward (54) and adjoint (57) problems with objects Ω and admittivity 𝛾_𝑖,𝑚.

Choosing 𝑣(𝐱)=−∇𝑢(𝐱)∇𝑝(𝐱),𝐱∈Ω, (59) we ensure 𝐽(𝛾_𝑖,𝑚+𝜂𝑣)<𝐽(𝛾_𝑖𝑚) if 𝜂>0 is small enough. We have defined 𝑣 only in Ω, but formula (59) holds for all 𝐱∈ℛ. When 𝛾_𝑖 is piecewise constant, we expect 𝛾_𝑖,𝑚 to be almost piecewise constant. In this case, another choice for the corrector 𝑣 is a piecewise constant approximation: 𝑣(𝐱)=𝑣_𝑗=−Re_Ω𝑗∇𝑢∇𝑝𝑑𝐳,𝐱∈Ω^𝑗,𝑗=1,…,𝑑. (60)

The reconstruction algorithm for the impedance tomography problem is similar to the algorithm described in Section 2.3 with 𝛾_𝑖, 𝛾_𝑒 playing the role of 𝑘_𝑖, 𝑘_𝑒, 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 𝑗_ℓ(𝑡)=cos(ℓ𝑡),𝑗_ℓ+5()[).=sinℓ𝑡,ℓ=1,…,5,𝑡∈0,2𝜋 (61) 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 𝛾_𝑖=8+2𝚤 and 𝛾_𝑒=1+𝚤 is assumed to be known. We started the hybrid algorithm with 𝛾_𝑖,0=5+3𝚤. The first guess for the domains was obtained observing the topological derivative with Ω=∅ and 𝛾_𝑖=𝛾_𝑖,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 𝛾_𝑖=8+2𝚤 is approximated by 7.63+2.30𝚤. The values of the conductivity 𝜎_𝑖 (the real part of 𝛾_𝑖) and of the permittivity 𝜀_𝑖 (the imaginary part of 𝛾_𝑖) 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 𝛾∈ℝ).

(a)

(b)

(c)

(d)

Figure 12 

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.

3.2. Tomography with Discrete Electrodes

The impedance tomography problem with discrete electrodes (43)-(44) can be reformulated as an optimization problem: Minimize 𝐽ℛ⧵Ω,𝛾_𝑖=12^𝑁_{𝑘=1𝑒𝑘}𝑢+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑢−𝑢_meas,𝑘²𝑑𝑙, (62) where 𝑢 is a solution of ∇⋅𝛾_𝑒∇𝑢=0inℛ⧵Ω,∇⋅𝛾_𝑖𝑢∇𝑢=0inΩ,₋−𝑢⁺𝛾=0on𝜕Ω,_𝑒𝜕_𝐧𝑢⁻−𝛾_𝑖𝜕_𝐧𝑢⁺𝛾=0on𝜕Ω,_𝑒𝜕_𝐧𝑢=0on𝜕ℛ⧵∪^𝑁_𝑗=1𝑒_j,_𝑒𝑘𝛾_𝑒𝜕_𝐧𝑢𝑑𝑙=𝐼_𝑘𝑘=1,…,𝑁. (63)

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 𝛾_𝑒 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 𝛾_𝑖 optimizing with respect to 𝛾_𝑖 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 𝐽(𝛾_𝑖)=𝐽(ℛ⧵Ω,𝛾_𝑖) when 𝑢 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 𝛾_𝑖,0 is usually available. Otherwise, we select a first approximation by perturbing 𝛾_𝑒∶𝛾_𝑖,0=𝛾_𝑒+𝜖. We may improve this first approximation iteratively by a gradient method. Given a guess 𝛾_𝑖,𝑚, we generate a descent direction 𝑣 to improve this guess studying the function of a real variable 𝐽(𝜂)=𝐽(ℛ⧵Ω,𝛾_𝑖,𝑚+𝜂𝑣). As in the acoustic setting and in the “continuous” model, 𝜂>0 and 𝑣 are selected in such a way that 𝐽^′(0)<0. An explicit expression for this derivative provides the following formula for the descent directions 𝑣.

Theorem 7. The derivative with respect to 𝜂 of the function 𝐽(𝜂)=𝐽(ℛ⧵Ω,𝛾_𝑖,𝑚+𝜂𝑣), with 𝐽 defined in (62) is 𝑑𝐽𝜂𝑑𝜂_𝜂=0=Re_Ω𝑣∇𝑢∇𝑝𝑑𝐳, (64) where 𝑢 is a solution of the forward problem (63) with parameter 𝛾_𝑖=𝛾_𝑖,𝑚 in Ω and 𝑝 is a solution of the adjoint problem (65) given below 𝑝,𝐏∈𝐻¹(ℛ)⊕ℝ^𝑁₀,𝑏(𝑤,𝐖),𝑝,𝐏=ℓ_𝑎(𝑤,𝐖),∀(𝑤,𝐖)∈𝐻¹(ℛ)⊕ℝ^𝑁₀, (65)(66) where 𝑏Ω;(𝑤,𝐖),𝑝,𝐏∶=_𝑅⧵Ω𝛾_𝑒∇𝑤∇𝑝𝑑𝐳+_Ω𝛾_𝑖∇𝑤∇+𝑝𝑑𝐳_𝑁𝑘=11𝑧_𝑘𝑒𝑘𝑤−𝑊_𝑘𝑝−𝑃_𝑘𝑑𝑙, (67)ℓ_𝑎(𝑤,𝐖)∶=−^𝑁_{𝑘=1𝑒𝑘}𝑢+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑢−𝑢_{𝑚𝑒𝑎𝑠,𝑘}𝑤+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑤𝑑𝑙, (68) with 𝛾_𝑖=𝛾_𝑖,𝑚 in Ω. Choosing 𝑣(𝐱)=−Re∇𝑢(𝐱)∇𝑝(𝐱),𝐱∈Ω (69) and 𝜂>0 small enough, one ensures 𝐽(𝛾_𝑖,𝑚+𝜂𝑣)<𝐽(𝛾_𝑖,𝑚).

Proof. The derivative with respect to 𝜂 is 𝑑𝐽𝛾𝑑𝜂^𝑖,𝑚+𝜂𝑣_𝜂=0=Re^𝑁_{𝑘=1𝑒𝑘}𝑢+𝑧_𝑘𝛾𝜕_𝐧𝑢−𝑢_meas,𝑘𝑢̇+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑢̇𝑑𝑙, (70) where 𝑢̇=(𝑑/𝑑𝜂)𝑢_𝜂|_𝜂=0 and 𝑢_𝜂 is a solution of (63) with 𝛾_𝑖=𝛾_𝑖,𝑚+𝜂𝑣; that is, 𝑢_𝜂,𝐔_𝜂∈𝐻¹(ℛ)⊕ℝ^𝑁₀𝑏𝜂;𝑢_𝜂,𝐔_𝜂,()(),∀()𝑤,𝐖=ℓ𝑤,𝐖𝑤,𝐖∈𝐻¹(ℛ)⊕ℝ^𝑁₀, (71) where 𝑏𝜂;(𝑢,𝐔),(𝑤,𝐖)∶=_ℛ⧵Ω𝛾_𝑒∇𝑢∇𝑤𝑑𝐳+_Ω𝛾_𝑖,𝑚+𝜂𝑣∇𝑢∇+𝑤𝑑𝐳_𝑁𝑘=11𝑧_𝑘𝑒𝑘𝑢−𝑈_𝑘𝑤−𝑊_𝑘𝑑𝑙, (72) and ℓ is defined in (50). Notice that when 𝜂=0 and 𝛾_𝑖=𝛾_𝑖,𝑚, the sesquilinear forms (67) and (72) are the same.
We eliminate 𝑢̇ by introducing the adjoint problem. Let us set ℒ(),𝜂𝑏(),,𝜂;𝑢,𝐔𝑝,𝐏∶=𝐽+Re𝜂;𝑢,𝐔𝑝,𝐏−ℓ𝑝,𝐏 (73) for (𝑢,𝐔),(𝑝,𝐏)∈𝐻¹(ℛ)⊕ℝ^𝑁₀. Then, choosing (𝑢,𝐔)=(𝑢_𝜂,𝐔_𝜂), it follows that 𝐽𝜂=ℒ𝜂;𝑢_𝜂,𝐔_𝜂,𝑝,𝐏,∀𝑝,𝐏∈𝐻¹(ℛ)⊕ℝ^𝑁₀ (74) and therefore 𝑑𝐽𝜂𝑑𝜂_𝜂=0=𝑑ℒ𝑢𝑑𝜂𝜂;_𝜂,𝐔_𝜂,𝑝,𝐏_𝜂=0𝑑=Re𝑏(),𝑑𝜂𝜂;𝑢,𝐔𝑝,𝐏_𝜂=0+Re𝑏0;𝑢̇,𝐔̇,𝑝,𝐏+Re^𝑁_{𝑘=1𝑒𝑘}𝑢+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑢−𝑢_meas,𝑘×𝑢̇+𝑧_𝑘𝛾𝜕_𝐧𝑢̇𝑑𝑙, (75) for all (𝑝,𝐏)∈𝐻¹(ℛ)⊕ℝ^𝑁₀. Selecting as (𝑝,𝐏) a solution of (65), we obtain the desired formula 𝑑𝐽𝜂𝑑𝜂_𝜂=0𝑑=Re𝑏(),𝑑𝜂𝜂;𝑢,𝐔𝑝,𝐏_𝜂=0=Re_Ω𝑣∇𝑢∇𝑝𝑑𝐳. (76)

Notice that if 𝑤∈𝐻¹(ℛ), then 𝑤∈𝐿²(𝜕ℛ)∩𝐻^1/2(𝜕ℛ) and 𝜕𝑤/𝜕𝑥_𝑖∈𝐻^−1/2(𝜕ℛ). Thus, the integrals appearing in the definition (68) of ℓ_𝑎 are in fact duality products.

We have defined 𝑣 only in Ω, but formula (69) holds for all 𝐱∈ℛ. When 𝛾_𝑖 is piecewise constant, we replace 𝑣 defined in (69) with piecewise constant approximations inside each object Ω^𝑗, 𝑣(𝐱)=𝑣_𝑗=−Re_Ω𝑗∇𝑢∇𝑝𝑑𝐳,𝐱∈Ω^𝑗,𝑗=1,…,𝑑. (77) 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 𝛾_𝑖 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 functional 𝐽(ℛ⧵Ω;𝛾_𝑖) defined in (62) is given by 𝐷_𝑇𝐱,ℛ⧵Ω,𝛾_𝑖=Re2𝛾_𝑒(𝐱)𝛾_𝑖(𝐱)−𝛾_𝑒(𝐱)𝛾_𝑒(𝐱)+𝛾_𝑖(𝐱)∇𝑢(𝐱)∇𝑝(𝐱),𝐱∈ℛ⧵𝐷Ω,_𝑇𝐱,ℛ⧵Ω,𝛾_𝑖=Re2𝛾_𝑖(𝐱)𝛾_𝑖(𝐱)−𝛾_𝑒(𝐱)𝛾_𝑖(𝐱)+𝛾_𝑒(𝐱)∇𝑢(𝐱)∇𝑝(𝐱),𝐱∈Ω, (78) where 𝑢 and 𝑝 solve the forward and adjoint problems (63) and (65) with coefficient 𝛾_𝑖 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) where the vector field 𝐕 is defined as 𝐕(𝐳)=−𝐧(𝐳),𝐳∈𝜕𝐵_𝜀(𝐱), (80) and extended to ℝ² in such a way that it vanishes away from a narrow neighborhood of 𝜕𝐵_𝜀(𝐱). Then, 𝐷_𝑇(𝐱,ℛ)=lim_𝜀→01ℎ^′(𝜀)𝑑𝐽𝜑𝑑𝜏^𝜏ℛ⧵𝐵_𝜀(𝐱)_𝜏=0,𝐱∈ℛ, (81) where ℎ^′(𝜀) is the derivative of the function ℎ(𝜀) that appears in the definition of the topological derivative (15). For the electrical impedance tomography problem in 2D, we can take ℎ(𝜀)=𝜋𝜀² as we will shortly see.
The proof is organized in three steps. We assume first that 𝛾_𝑒 and 𝛾_𝑖 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 𝛾_𝑒 and 𝛾_𝑖 are constant parameters). Let us compute first 𝑑𝐽𝜑𝑑𝜏^𝜏ℛ⧵𝐵_𝜀(𝐱)_𝜏=0. (82) Since the field 𝐕 decays fast enough away from 𝜕𝐵_𝜀(𝐱), it follows that 𝜑_𝜏(𝜕ℛ)=𝜕ℛ and 𝜑_𝜏(ℛ⧵𝐵_𝜀(𝐱))=ℛ⧵Ω_𝜏, where Ω_𝜏∶=𝜑_𝜏(𝐵_𝜀(𝐱)). The cost functional in the transformed domains is then 𝐽ℛ⧵Ω_𝜏=12^𝑁_{𝑘=1𝑒𝑘}𝑢_𝜏+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑢_𝜏−𝑢_meas,𝑘²𝑑𝑙, (83) where 𝑢_𝜏 solves 𝑢_𝜏,𝐔_𝜏∈𝐻¹(ℛ)⊕𝑅^𝑁₀,𝑏Ω_𝜏;𝑢_𝜏,𝐔_𝜏,(()()𝑤,𝐖)=ℓ𝑤,𝐖,∀𝑤,𝐖∈𝐻¹(ℛ)⊕𝑅^𝑁₀, (84) with 𝑏 defined in (67) and ℓ defined in (50). Notice that when 𝜏=0, Ω_𝜏=𝐵_𝜀(𝐱), and the solution of (84) is the solution (𝑢^𝜀,𝐔^𝜀) to (63) with Ω=𝐵_𝜀(𝐱).
The derivative of (83) with respect to 𝜏 involves the derivative 𝑢̇ of 𝑢_𝜏 with respect to 𝜏. To avoid its computation we introduce the Lagrangian functional ℒΩ_𝜏;(𝑢,𝐔),𝑝,𝐏∶=𝐽ℛ⧵Ω_𝜏+Re𝑏Ω_𝜏;(𝑢,𝐔),𝑝,𝐏−ℓ𝑝,𝐏. (85) Then, for all (𝑝,𝐏)∈𝐻¹(ℛ)⊕𝑅^𝑁₀, 𝑑𝐽𝑑𝜏ℛ⧵Ω_𝜏𝜏=0=𝑑ℒΩ𝑑𝜏^𝜏;𝑢_𝜏,𝐔_𝜏,𝑝,𝐏_𝜏=0𝑑=Re𝑏Ω𝑑𝜏^𝜏;𝑢^𝜀,𝐔^𝜀,𝑝,𝐏_𝜏=0+Re𝑏𝐵_𝜀(𝐱);𝑢̇,𝐔̇,𝑝,𝐏+Re^𝑁_{𝑘=1𝑒𝑘}𝑢^𝜀+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑢^𝜀−𝑢_meas,𝑘×𝑢̇+𝑧_𝑘𝛾_𝑒𝜕_𝐧𝑢̇𝑑𝑙. (86) If (𝑝^𝜀,𝐏^𝜀) satisfies (65) with object Ω=𝐵_𝜀(𝐱) and parameter 𝛾_𝑖, then 𝑑𝐽𝑑𝜏ℛ⧵Ω_𝜏𝜏=0𝑑=Re𝑏Ω𝑑𝜏^𝜏;𝑢^𝜀,𝐔^𝜀,𝑝^𝜀,𝐏^𝜀_𝜏=0=Re_{ℛ⧵𝐵𝜀(𝐱)}𝛾_𝑒∇𝑢^𝜀⋅∇𝑝^𝜀−∇⋅𝐕𝑑𝐳_{ℛ⧵𝐵𝜀(𝐱)}𝛾_𝑒∇𝐕+∇𝐕^⊤∇𝑢^𝜀⋅∇𝑝^𝜀+𝑑𝐳_{𝐵𝜀(𝐱)}𝛾_𝑖∇𝑢^𝜀⋅∇𝑝^𝜀−∇⋅𝐕𝑑𝐳_{𝐵𝜀(𝐱)}𝛾_𝑖∇𝐕+∇𝐕^⊤∇𝑢^𝜀⋅∇𝑝^𝜀𝑑𝐳=Re_{𝜕𝐵𝜀(𝐱)}𝛾_𝑖𝛾_𝑖𝛾_𝑒−1𝜕_𝐧𝑢^𝜀−𝜕_𝐧𝑝^𝜀−+𝛾_𝑖−𝛾_𝑒𝜕_𝐭𝑢^𝜀−𝜕_𝐭𝑝^𝜀−𝑑𝑙. (87) To obtain the last equality in (87), we have integrated by parts and used that 𝑢^𝜀 and 𝑝^𝜀 solve (63) and (65) with Ω=𝐵_𝜀(𝐱), respectively, and that 𝐕 vanishes on 𝜕ℛ.
By (81), the topological derivative at 𝐱 can be computed now as 𝐷_𝑇(𝐱,ℛ)=lim_𝜀→01ℎ^′(𝜀)Re_{𝜕𝐵𝜀(𝐱)}𝛾_𝑖𝛾_𝑖𝛾_𝑒−1𝜕_𝐧𝑢^𝜀−𝜕_𝐧𝑝^𝜀−+𝛾_𝑖−𝛾_𝑒𝜕_𝐭𝑢^𝜀−𝜕_𝐭𝑝^𝜀−𝑑𝑙. (88) The asymptotic behavior of 𝑢^𝜀 and 𝑝^𝜀 at 𝜕𝐵_𝜀(𝐱) is obtained expressing them as corrections of 𝑢⁰ and 𝑝⁰, the solutions to (63) and (65) with Ω=∅, 𝑢^𝜀(𝐳)=𝑢⁰(𝐳)𝜒_{ℛ⧵𝐵𝜀(𝐱)}(𝐳)+𝑣^𝜀(𝐳),𝑝^𝜀(𝐳)=𝑝⁰(𝐳)𝜒_{ℛ⧵𝐵𝜀(𝐱)}(𝐳)+𝑞^𝜀(𝐳), (89) and expanding the remainders in powers of 𝜀. Reasoning as in [23] we find that 𝑢^𝜀(𝐳)⟶𝑢⁰(𝐱),∇𝑢^𝜀−(𝐳)⟶2𝛾_𝑒𝛾_𝑒+𝛾_𝑖∇𝑢⁰(𝐱),as𝜀⟶0, (90) uniformly when |𝐳−𝐱|=𝜀. For 𝑝^𝜀 there is a slight technical difference because 𝑝^𝜀 is defined by the variational problem (65). However, this is due to the source term at 𝜕ℛ, that is included in 𝑝⁰. Once that source term is eliminated, the remainder 𝑞^𝜀(𝐳) satisfies equations similar to the equations for 𝑣^𝜀 and we find 𝑝^𝜀(𝐳)⟶𝑝⁰(𝐱),∇𝑝^𝜀−(𝐳)⟶2𝛾_𝑒𝛾_𝑒+𝛾_𝑖∇𝑝⁰(𝐱),as𝜀⟶0, (91) uniformly when |𝐳−𝐱|=𝜀.
Replacing the limits (90) and (91) in (88) and choosing ℎ(𝜀)=𝜋𝜀², 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 𝐱 and repeat the previous computations. Notice that due to the choice of 𝐕, only the values of 𝛾_𝑖 and 𝛾_𝑒 in a narrow band about 𝜕𝐵_𝜀(𝐱) are relevant for the computation of the shape derivative. For the topological derivative, we set 𝛾_𝑖(𝐳)=𝛾_𝑖(𝐱)+(𝑧−𝐱)⋅∇𝛾_𝑖(𝐱)+𝑂(𝜀²) and an analogous expression for 𝛾_𝑒. Then, we check that the limits (90)-(91) also hold with 𝛾_𝑖=𝛾_𝑖(𝐱) and 𝛾_𝑒=𝛾_𝑒(𝐱). For continuous 𝛾_𝑒 and 𝛾_𝑖, 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 Ω⧵𝐵_𝜀(𝐱) for 𝐱∈Ω, and Ω∪𝐵_𝜀(𝐱) when 𝐱∈ℛ⧵Ω depending on whether we are computing the derivative inside or outside Ω; see [21] for details.

Acknowledgments

The authors are partially supported by the Spanish Government Research Project TRA2010-18054 and the Spanish Ministerio de Economia y Competitividad Grants no. FIS2011-28838-C02-02 and no. FIS2010-22438-E.