Topological Imaging Methods for the Iterative Detection of Multiple Impedance Obstacles

In this paper, we investigate shape inversion algorithms based on the computation of iterated topological derivatives for the detection of multiple particles coated by a complex surface impedance in two- and three-dimensional acoustic media. New closed-form formulae for the topological derivative of the misfit functional are derived when an approximate set of unknown particles has already been recovered. Proofs rely on the computation of shape derivatives followed by the topological asymptotic analysis of a boundary integral equation formulation of the forward and adjoint problems. The relevance of the theoretical results is illustrated by various 2D and 3D experiments using monochromatic imaging algorithms either fully or partially based on topological derivatives.


Introduction
In the continuation of our recent paper [35], we investigate the performance of the topological derivative (TD) for developing full shape inversion methods from noisy measurements when an impedance condition is imposed on the boundary of the multiple unknown particles. We already derived in [35] the topological gradient formulae for the 2D and 3D time-harmonic Helmholtz equation in the free space. In that work, we presented a large gallery of numerical experiments illustrating the performance of the instantaneous method when multi-frequency noisy data are available. In the single-frequency situation, we observed that first stage approximations were not always satisfactory, in the sense that any level set of the topological indicator function do not provide in many cases the true number and location of the defects or a reasonable shape identification. The aim of this paper is to improve such approximation by iteratively computing new TDs.
Full imaging algorithms of thin coated objects using very reduced data are of high importance in many fields of applied physics such as nondestructive testing, biomedical and radar imaging [8][9][10].
The inverse shape problem is commonly reformulated as the minimization problem of the least-squares functional, which accounts for the difference between the measured and reconstructed data. The topological gradient of this misfit functional measures the variations of such functional in the presence of an infinitesimal scatterer at each point of the region of interest. While a general abstract definition of the first-order TD has been set for any small shape domain [22,41,45], we restrict our study to small balls. Former numerical comparison in the literature of the TD formulae depending on various choices of shape perturbations did not demonstrate significant differences (see [26]). The topological gradient can also be reinterpreted as a limit of shape derivatives [40] for which an adjoint problem is usually introduced. In the absence of scatterers, the topological asymptotic analysis can then be performed using 2D Fourier and 3D Mie series expansions of the solution to the forward and adjoint problems [33,35]. However, extending the results in the presence of an already set of approximate shapes is not immediate. We adopt here a new strategy based on a boundary integral equation formulation of multiple scattering problems, first proposed by the authors in the electromagnetic case [34]. This approach is well suited to obstacle scattering frameworks. Indeed, using the 2D Fourier and 3D Mie series expansions of the potential operators and their traces, one can establish the asymptotic behaviors of the perturbed solution to the forward and adjoint problems far from and on the boundary of the infinitesimal particles. First-order TD formulae of classical misfit functionals are deduced readily.
From a computational point of view, the indicator function formulae are relatively easy to implement with either FEM or BEM solvers and a few parallel computations in the three-dimensional case. As a consequence, various TD-based iterative methods have been investigated in the literature. Iterative methods completely based on the computation of TDs were first introduced for optimal design of structures [22,25], and then adapted for inverse scattering problems for the Helmholtz equation in 2D in [13][14][15][16] and 3D in [12,46]. The TD can also be used to generate initial guesses in shape derivatives-based methods such as the level set method [20,21]. A new trend consists in crossing TDs with parametric derivatives through converging Gauss-Newton iterations [30]. The latter algorithm has been recently applied to the challenging imaging problem of a few electromagnetic particles from two-dimensional holographic measurements [11]. New unknown particles are automatically recovered by TDs when Gauss-Newton iterations stagnate. To our best knowledge, the performance of iterated TD-based methods for the impedance boundary condition case has not been analyzed yet. The numerical results presented here are proper to the Helmholtz equation case and by no means can be obtained from the Maxwell equation case with an impedance boundary condition [33].
The paper is organized as follows: In Sect. 2, we present an equivalent boundary integral reformulation of the direct scattering problems for multiple complex 2D and 3D impedance obstacles. The associated inverse obstacle problem is also briefly outlined. The topological derivative is introduced in Sect. 3 and algorithms based on iterated TDs are described. Using shape derivatives tools [28] and spectral expansions of the boundary integral operators for ball-shaped domains, we derive closed-form formulae of the TD of the leastsquares cost functional. We show numerically in Sect. 4 that iterating the single-frequency TD indicator provides a reasonable reconstruction under high impedance contrast, even when only partial aperture data is available. Moreover, combined with parametric derivatives through converging Gauss-Newton methods, it automatically and accurately distinguishes various particles with a lower complex impedance contrast. Finally, we outline concluding remarks and discuss possible research lines in Sect. 5. The paper ends with a short appendix where we collect some relations and asymptotic properties of Bessel functions.

Notation 1 For a bounded open domain
and H s (∂ ) the standard (local in the case of the exterior domain) complexvalued, Hilbert-Sobolev space of order s ∈ R defined on , R d \ and ∂ , respectively (with the convention H 0 = L 2 ).

Forward and Inverse Problems
In this section, we introduce first the forward impedance scattering problem and recall an equivalent boundary integral formulation proposed in [36]. Next, the inverse impedance scattering problem is presented.

Forward Problem
Let us consider a bounded domain ⊂ R d , for d = 2 or 3. We assume that consists of N ≥ 1 disjoint components with smooth simply connected boundaries, namely, The scattering of time-harmonic (exp(−iωt) time dependence) acoustic waves when each is an impedance object can be modeled by the following forward problem for the total field u [19]: for a given incident acoustic wave u inc solving the Helmholtz equation in the absence of any object, find the total field u satisfying where κ = ω/c is the wave number (ω > 0 is the frequency and c > 0 is the wave speed), and λ ∈ C is the impedance parameter, assumed to satisfy Here, n is the outward unit normal vector to ∂ (to each disjoint component of ). The last condition in (1) is the wellknown Sommerfeld radiation condition at infinity that has to be satisfied uniformly for every unitary direction x = x |x| . It is well known that under the assumption (2), problem (1) has a unique solution (see [36]). From the Sommerfeld radiation condition at infinity, we also have [19,Theorem 2.6.]: uniformly for all unitary directions x. The function u ∞ , which is defined on the unit sphere S d−1 of R d , is called the far-field pattern. Let us introduce now a boundary integral formulation of the forward problem (1) on which the derivation of our forthcoming TD formula will be based. To do that, we define for any densities ψ ∈ H − 1 2 (∂ ), ϕ ∈ H 1 2 (∂ ), and any point x ∈ R d \∂ the standard single and double-layer operators, S κ and D κ , respectively, as where G is the fundamental solution to the Helmholtz equation [39, Section 2.1] and n is the outward unit normal vector to ∂ . Then, by the Helmholtz integral representation formula in R d \ [18,Theorem 3.3], the solution u ∈ H 1 loc (R d \ ) to the impedance problem (1) can be written as: Taking in (6) the exterior Dirichlet and Neumann traces of u on ∂ for = 1, . . . , N , by the jump relations of the single and double-layer operators, one finds the identities where the boundary integral operators S κ , D κ , D κ and N κ are defined for any densities ψ ∈ H − 1 2 (∂ ), ϕ ∈ H 1 2 (∂ ) and any point x ∈ ∂ by We remark here that we are using the standard notation in the boundary integral community for the operator D κ , where the prime symbol is used to indicate that D κ and D κ are adjoint to each other with respect to the duality product in Theorem 1 [36,Theorems 3.4 and 3.5] Under the hypothesis (2), the solution u ∈ H 1 loc (R d \ ) of the impedance scattering problem (1) can be expressed as where the Dirichlet traces ϕ := u |∂ , = 1, . . . , N are the solution to the following system of boundary integral equations:

Inverse Problem
In this work, we aim at finding the number N , the location, and the approximate shapes of a set of unknown impedance obstacles = ∪ N =1 from measurements of either the near-field or far-field patterns resulting from the scattering of a finite number N inc of incident waves. For the numerical experiments, we will consider plane waves of the form u inc r (x) = e iκx·d r , r = 1, . . . , N inc , where d r ∈ S d−1 . If we denote by u r the total field solution to the forward impedance problem (1) for a given incident wave u inc r , then in case of near-field data, the inverse problem consists in recovering the scatterers that minimize the cost function 1 2 meas u r (y obs ) − u meas,r (y obs ) 2 dμ(y obs ), (11) where meas is the set of observation points and u meas,r (y obs ) are the measured data at the points y obs ∈ meas . The measure μ represents the Lebesgue measure on a surface of R d or, in practical applications, it rather represents a finite sum of Dirac measures. In the later case, the cost function reads where y s obs , for s = 1, . . . , N obs , are the observation points in R d . In case of far-field measurements, the cost function is where y s obs , for s = 1, . . . , N obs , are a finite number of observation points in S d−1 and u ∞ meas,r ( y s obs ) is the measured far-field at those points.

Topological Derivative-Based Methods
The minimization strategy will be based on the computation of the topological derivative (also known as the topological gradient or the topological sensitivity) of the functional J defined in (11) or (12), for near-field or far-field observations, respectively. In this section, we first provide the general definition of the topological derivative, to introduce then our iterative minimization strategy (Sect. 3.1). Closedform formulae, needed for the practical implementation of the algorithm, will be derived in Sect. 3.2.

Topological Derivative: Theory and Minimization Algorithm
The topological derivative (TD) of a given cost function J defined in a region R ⊂ R d measures the sensitivity of J to ball-like perturbations of infinitesimal radius at each point of R [41,45], namely, for each x ∈ R it provides an asymptotic expansion of the form where B ε (x) is an infinitesimal ball centered at x of radius ε > 0, and f > 0 is a monotonically increasing function such that f (ε) → 0 as ε → 0. In our previous paper [35] we proved that for the impedance problem in R = R d , the function f must behave as ε d−1 . In the forthcoming section we will prove that also for R = R d \ app , where app is a bounded domain (an approximated guess for ), f must behave as ε d−1 .
A more general definition of the TD is also possible by replacing the role of B ε (x) by a general infinitesimal domain containing the point x, but in practice, closed-form formulae become rather involved. For this reason, we will restrict our definition to balls.
In view of the expansion (13), intuitively one can conclude that placing infinitesimal objects at the points x ∈ R where the topological derivative D T (x; R) attains the most pronounced negative values, the cost functional should decrease. Therefore, one can consider the following set 0 as an initial guess for the unknown objects : (14) where 0 < C 0 < 1 is a tunable constant. Some partial justifications of this heuristic choice for some specific problems in acoustics can be found in [5,27].
We propose to update the initial guess 0 by subsequent TD computations: the region R \ 0 plays now the role of R in the expansion (13) to find the TD indicator function D T (x; R \ 0 ) for all x ∈ R \ 0 . Then, 0 can be updated by adding the set of points where the largest negative values of this derivative are attained, namely, the new approximated domain is defined as for some tunable constant 0 < C 1 < 1. The method can be straightforwardly iterated now. This method was firstly proposed in [14] for the 2D-transmission problem for the Helmholtz equation, and later used in [13] for the Dirichlet problem. It has also been used in other contexts, like for 2D photothermal imaging [15], the 2D electrical impedance problem [16] or for 3D electromagnetic scattering in [34].
Notice that the iterative algorithm generates a nested sequence of approximated domains, namely j ⊂ j+1 . Therefore, if at some step a spurious region is added, it cannot be removed in the subsequent iterations. For the transmission problem, the TD can be also computed inside the current approximations, being able to both include and remove regions from one iteration to the next one [14].
In the next section we give closed-form formulae of the TD for the near-field and far-field least-squares cost functions both for R = R d and for R = R d \ app , with app being a given approximation of the true domain .

Closed-Form Topological Derivative Formulae
In absence of any initial guess for , we set R = R d in the definition (13). The following results, corresponding to near-field and far-field data, were proved in [35]: Theorem 3.4] For the cost function J defined in (11), we have the expansion, for x ∈ R d , and the topological derivative is given by Theorem 3 [35,Theorem 3.5] For the far-field cost function J defined in (12) we have the expansion, for x ∈ R d , where f is given by (15) and the topological derivative is where The remaining of the section is devoted to the derivation of closed-form formulae when R = R d \ app . To do that, we follow the strategy proposed in [40] (and used in our previous paper [35]), to obtain the TD as a limit of shape derivatives. We will start with the near-field case. For notational simplicity, we assume that N inc = 1 in the definition of the functional (11) and omit the subindex r . The extension to N inc > 1 is straightforward by linearity. Theorem 4 Let app be a given approximation of the true domain . Let x ∈ R d \ app . Then, the topological derivative at x of the cost function (11) (for N inc = 1) can be obtained as (13), and where u ε solves the forward problem (1) with = app ∪ B ε (x) and w ε is the complex-conjugate of the solution w ε to the associated adjoint problem, namely, it solves Proof As stated in [40], the TD of a given functional J defined in a region R ⊂ R d can be computed at any x ∈ R as follows: take ε > 0 small enough satisfying B ε (x) ⊂ R, and define a vector field V such that n ε being the outward normal vector on ∂ B ε (x), and such that it vanishes outside a narrow neighborhood of ∂ B ε (x). Consider then the deformations ϕ τ for τ > 0 defined as Then, the TD of J at x ∈ R is given by where f is the derivative of the function f appearing in the definition of the topological derivative (13). Following now the same steps as in the beginning of Sect. 3.2 in [35], we obtain that where u ε is the solution to (1) with = app ∪ B ε (x) anḋ u ε solves the following problem [28, Theorem 2.3]: The proof now is a simple adaptation to the proof of Theorem 3.1 in [35], In order to derive a closed-form formula for the TD, we should study now the asymptotic behavior of the functions u ε and w ε to derive the asymptotic behavior of g ε , which in turn, determines how f must behave (we will see that it must be f (ε) = O(ε d−2 )), and therefore how f (ε) must behave. This will be carried out by using Fourier and Mie series expansions (for the two-and three-dimensional cases, respectively), combined with the boundary integral formulation of problem (1) presented in Theorem 1.
Let us briefly recall that for any x ∈ R d , a general incident field u inc satisfying u inc + κ 2 u inc = 0 in R d can be written for any z ∈ R d as follows [17,Eqs. (1) and (2)]: where the complex-valued sequence (α (1) n, j ) only depends on x, and for ε > 0 and ξ ∈ S 2 . The functions j n and Y n, j are the spherical Bessel functions of the first kind [ where the complex-valued sequence (α (1) n ) only depends on x, and u (1) n (κ, ε(cos θ, sin θ)) = J |n| (κε)e inθ .
Here, J n is the Bessel function of the first kind and order n ( [1, Section 9.1]).

=1
, and then, by Theorem 1 the field u ε can be expressed in terms of its Dirichlet traces on the boundaries ∂ app and ∂ B ε (x).
Let us introduce the function τ ε : ∂ B ε (x) → S 2 defined by τ ε (y) = (y−x)/ε and consider the family of shifted spherical harmonics Y ε n, j = Y n, j • τ ε for n ≥ 0 and | j| ≤ n, that are eigenfunctions of the acoustic boundary integral operators. To study the asymptotic behavior of u ε when ε → 0, we expand its Dirichlet trace on ∂ B ε (x) in terms of Y ε n, j : for any Let us define now the single S ε κ and double D ε κ layer operators related to the boundary ∂ B ε (x) (defined as in (4)). We want to study first the asymptotic behavior of (S ε κ Y ε n, j )(y) and To do that, we remark that for any y, z ∈ R 3 such that |y − x| > |z − x|, we have [19,Theorem 2.11] In the special case z ∈ ∂ B ε (x) and y ∈ ∂ app ⊂ R 3 \B ε (x), then Using now that and the expansions (47) collected in the appendix, we obtain the following formulae for y ∈ ∂ app : where j n is the derivative of the spherical Bessel function j n . Let us denote now by S ε κ , D ε κ , D ε κ and N ε κ the boundary integral operators related to the boundary ∂ B ε (x), namely, the counterpart operators to those defined in (9). They are mean values of the interior and exterior Dirichlet and Neumann traces on ∂ B ε (x) of the layer potentials S ε κ and D ε κ . Since the shifted spherical harmonics are eigenfunctions of these operators, we get by (47) and (48) the following formulae: being the derivative of the spherical Bessel function h (1) n .
(ii) We denote by I op , I op|∂ B ε (x) and I op|∂ app the invertible boundary integral operators associated with equation (10) for = app ∪ B ε (x), = B ε (x) and = app , respectively. Then, in view of (28), we have and where the asterisk * in the first equality refers to the series of Dirichlet and Neumann traces of the potential operators defined by integrals over app and evaluated on ∂ B ε (x) that are located in the system (10) between brackets. We only emphasize that the matrix operators are upper triangular, since the unspecified operators are irrelevant for us and the structure of the inverse operators is known to be identical. The integral operator I op|∂ B ε (x) is reduced to an invertible diagonal operator whose eigenvalues are up to O(1) for n = 0 and up to O(ε −1 ) for n = 0 when ε → 0. From (22) and (47), we can write for any where A (1) n, j (x) can be expressed in terms of α (1) n, j (x) and do not depend on ε. It follows, by inverting the system, that the Dirichlet trace of u ε on ∂ B ε (x) admits the same asymptotic behavior as the incident field. It means that the leading term is only given by n = 0 and we have for Now, equation (7) on the boundary ∂ B ε (x) can be written as We have S ε κ Y ε 0,0 = O(ε) and We derive then that Taking the limit ε → 0 in the above equation we obtain the assertion (25). To obtain (26), we apply now the tangential gradient operator to the boundary integral equation (7). The first term (for n = 0) vanishes, then the leading terms are given by n = 1.
we get Taking the limit ε → 0 in the above equation we obtain the assertion (26). Let us consider now the case d = 2. To show (24)- (26) we just need to adapt the previous arguments. First of all, for any z ∈ R 2 we define the angle θ(z) such that z = x+|z−x|e iθ(z) . Then, for y = x + ε(cos θ(y), sin θ(y)) ∈ ∂ B ε (x) we define the shifted cylindrical function E ε n (y) = e inθ(y) . On the other hand, for any y, z ∈ R 2 such that |y − x| > |z − x| we have Then, if z ∈ ∂ B ε (x) and y ∈ ∂ app ⊂ R 2 \ B ε (x) (then |y − x| > |z − x|), it follows that We likewise obtain the expression for the double-layer potential for y ∈ ∂ app : where J |n| is the derivative of the Bessel function J |n| . Using now (44) and (46), we find the asymptotic expansions for J |n| (kε) and J |n| (kε), from where we deduce for y ∈ ∂ app : Using now that the boundary integral operators are mean values of the interior and exterior Dirichlet and Neumann traces on ∂ B ε (x), we deduce that E ε n are eigenfunctions of these operators and furthermore, they behave as follows: |n| being the derivative of the Hankel function H (1) |n| . The remaining of the proof for deriving (24)-(26) is a straightforward adaptation of its counterpart in 3D: (24) is proved as in (i). For (25) and (26), we reason as in (ii), noticing that for any z ∈ ∂ B ε (x) we can write the two-dimensional version of (29): and deduce that u ε behaves as u ε (z) = A (1). We obtain then that (30) is also valid in 2D, and since S ε (32), and therefore (25) holds. To show (26), we apply the tangential gradient operator to the boundary integral equation (7) and use that S ε κ E ε n = O(ε) and D ε κ E ε n = O(ε) for |n| = 1. From here, we get and taking the limit ε → 0 we find (26).
Finally, to derive the results for w ε (for d = 2 and d = 3), we introduce the solution w inc ε to the adjoint problem which is explicitly given by w inc ε (z) = meas G(κ, z−y obs )(u ε (y obs ) − u meas (y obs ))dμ(y obs ).
It suffices now to decompose w ε = w scat ε + w inc ε , and notice that w scat Using the corresponding Mie (if d = 3) or Fourier (if d = 2) series expansion of the fundamental solution (5) we can write, when |z − x| < |y obs − x| (and in particular when where a (1) n (x, ε) and a (1) n, j (x, ε) do not depend on z. We also decompose w 0 = w scat 0 + w inc 0 , with u ε replaced by u 0 to define w inc 0 as in (35). Then, from (24) we find that lim ε→0 w inc ε (y obs ) = w inc 0 (y obs ).
The counterpart results to (25) and (26) follow now by considering the boundary integral formulation of problem (36) and repeating the arguments above replacing the role of u ε by w ε .
Finally, we can derive the formula for the TD when an approximation app has been already set: (11), we have the expansion

Theorem 5 For the cost function J defined in
where f is the function defined in (15), and the topological derivative is given by where u r ,0 is the solution to (1) with = app for the incident field u inc r , and w r ,0 solves (27) with u 0 = u r ,0 and u meas = u meas,r .
Proof Let u r ,ε be the solution to (1) with = app ∪ B ε (x) for the incident wave u inc r , and w r ,ε the solution to (21) for u ε = u r ,ε and u meas = u meas,r . Then, by Proposition 1 we have Then, by Theorem 4 we deduce now that f (ε) must behave as ε d−2 , and by choosing f namely, defining f as in (15), we find formula (37).
Finally, for the far-field data case we have the following result.

Theorem 6
For the far-field cost function J defined in (12) we have the expansion

where f is given by (15) and the topological derivative is
where u r ,0 is the solution to (1) for = app for the incident field u inc r , and w r ,0 can be computed as w r ,0 = w scat r ,0 + w inc with δ(d) as in (19), and w scat r ,0 solving Proof In [33] we proved that the TD for far-field data can be obtained as the limit when R → ∞ of the near-field topological derivative associated with the measurement curve/surface ∂ B R (0), where B R (0) is the ball centered at the origin 0 and radius R. Using the asymptotic behavior (3) of the scattered fields, and formula (45) for the two-dimensional case, the result follows readily (see also [35,Theorem 3.5] for the counterpart result for app = ∅).

Numerical Experiments
In this section, we provide some numerical examples to illustrate the performance of the proposed iterative method. For the sake of brevity, we just concentrate in a few examples to visualize specific features of the method, without providing a wide variety of simulations varying all the involved parameters in the scattering problem (namely, value of the wave number, of the impedance parameter, number of defects, their location and shape, number of incident waves and aperture angles, number of observation points and their location, nearfield or far-field data, et cetera). For a very ample gallery of experiments dealing with these parameters for the one-step method, we refer to our previous paper [35]. The remaining of the section is as follows. First we comment about the generation of synthetic data and the implementation of the TD formulae. Then, numerical experiments for the two-and three-dimensional situations will be shown.

Synthetic Data
Synthetic data both in 2D and 3D were generated by solving the forward problems for the true defects. True objects will be described by smooth parameterizations. Being more precise: (i) In the two-dimensional case, each connected component will be described by a smooth regular 1-periodic function x : R → ∂ , satisfying |x (t)| = 0 for all t ∈ R and x (t) = x (s) for s − t / ∈ Z. (ii) In the three-dimensional case, each connected component will be described by a smooth spherical parameterization x : S 2 → ∂ , such that every object is described by its center c and a spherical parameterization q , namely The forward problem is then solved by considering the equivalent system of boundary integral equations given in Theorem 1 fully discretized by Galerkin approaches: (i) In the two-dimensional case, we use the fully discrete version of the spectral method described in [19,Section 3.5] (see [38] for further details), which has superalgebraic convergence order, as tested in [43,44] for a related transmission problem where the same single-and double-layer boundary integral operators are involved. (ii) In the three-dimensional case, we implemented the spectral algorithm proposed in [24,31], which enjoys also superalgebraic convergence order.
To both avoid inverse crimes and to test robustness with respect to noise, in all the experiments we add a relative random noise to the numerically generated data described above. Specifically, once we obtain numerically the scat- a random noise of size δ is added to generate the synthetic measured data defined as u δ meas = u inc + u δ scat , where u δ scat − u scat 2 / u scat 2 ≤ δ. Unless otherwise stated, we set δ = 0.05.
To further avoid inverse crimes, most of the true defects will not be star-shaped ones, meaning that they will not belong to the class of objects where they will be approximated, as explained in the forthcoming subsection.

Implementation of the Algorithm
As stated in (14), the initial guess for is computed by selecting the points where the topological derivative D T ( · , R d ) (in case of near-field data) or D ∞ T ( · , R d ) (in case of far-field data) attains its largest negative values. For simplicity, we only consider here the near-field case, but the results for farfield data are analogous. In practice, rather than in the whole space R d , we evaluate it in a big enough bounded region R d obs ⊂ R d where the scatterers are assumed to be located. In our experiments, we set R d obs = [−3, 3] d . We discretize this region in a structured grid with (M + 1) d points, considering the discrete region Then, the initial guess is defined as where 0 < C 0 < 1 is a tunable parameter: the higher C 0 the more points contains 0 .
The evaluation of D T (x, R d ) for x ∈ R d M is very simple. By Theorem 2, the formula is (16), and we just need to evaluate the incident wave u inc r (x) = e ikx·d r , and the adjoint field w inc r (x) at the points x ∈ R d M using formula (17). The computational cost is therefore almost negligible since we do not have to solve any forward or adjoint problem numerically, we just need to evaluate functions that are known in closed form.
As explained in [33], the condition can be collected in a binary matrix of size (M+1) d , which can be processed in matlab with the command bwconncomp to obtain the number of connected components of 0 as well as the points belonging to each component. Then, each connected component is approximated by a star-shaped object, namely, we look for a radial parametrization of the boundary of each genus-zero component. Details about the implementation of the star-shaped approximation can be found in [14] for the two-dimensional case, and in [33] for the threedimensional one. A wide gallery of numerical examples for the one-step method, namely, when only 0 is computed following the procedure above, can be found in [35]. Now, 0 is iteratively updated by considering, for j = 1, 2, . . . , the approximations (40) As stated in Theorem 5, the evaluation of M requires now to solve a couple of problems for each incident wave u inc r : we need to solve numerically problems (1) and (27) Then, the same boundary integral equation method as for the forward problem can be applied. Indeed, both systems share the same matrix which has only to be computed once. Furthermore, this matrix is independent of the considered incident wave. Therefore, we only have to solve 2N inc linear systems of equations, all of them sharing the same matrix. In principle, the constant C j in (40) is equal to C j−1 . However, at each iteration we check if the cost function decreases.
, then j is accepted as new approximated domain. Otherwise, we replace the value of C j by C j / √ 2 and compute again j by (40). We stop the iterations if either: (i) A maximum number of iterations is reached. (ii) The method stagnates, namely, if the difference between consecutive objects produces almost no effect in the cost function, i.e., if where u j r is the solution to the forward problem (1) for the incident wave u r inc and for the scatterers j , while u δ meas,r (y s obs ) and u δ scat,r (y s obs ) are the total and scattered noisy data measured at the observation point y s obs .

Some Test Cases
We consider first a two-dimensional example with two impedance obstacles: a ball and a rounded triangle, with impedance parameter λ = 10 + i (see Fig. 1a). Synthetic measured data is generated as explained in Sect. 4.1 at 64 observation points uniformly distributed on the circle of radius 10 (marked by blue crosses in Fig. 1a) for 35 plane incident waves u inc r (x) = e iκx·d r with κ = 1 and directions d r = (cos(2π r /35), sin(2π r /35)), for r = 1, . . . , 35 (marked in Fig. 1a by black arrows). A relative random noise of size δ = 0.05 is added.  Fig. 1b, where we observe that the largest negative values of the TD (dark-red colors in the plot) are attained in two disjoint regions, inside the true defects. The boundaries of the true defects are represented by a black line for an easier visualization. By selecting the value C 0 = 0.05 to define the initial guess 0 by formula (39), we obtain the two objects represented in white in Fig. 1c. Both resemble ellipsoid-like defects, being the one corresponding to the true ball the biggest in size. The values of the TD computed for these objects are also represented in the same plot. We realize now that the other object, corresponding to the rounded triangle should be increased in size, since the largest negative values of the TD are mainly concentrated inside that object. During the next iterations, the size of this object increases while the other one remains almost unchanged, as can be seen in Fig. 1d, corresponding to the 7-th iteration. At the iteration 14, both objects seem to be comparable in size (see Fig. 1e). The algorithm stopped at the 27-th iteration, providing the reconstruction shown in Fig. 1g. The values of the cost functional versus the number of iterations are represented in Fig. 1h.
When considering very small values of the constant C 0 in (39), the initial guess is expected to underestimate the size of the scatterers, and even to disregard the presence of some objects. The number of iterations in this case is expected to increase. However, on the other hand, since the size of the objects is expected to increase less from one iteration to the next one, the method may be in this case less prone to include spurious regions. To illustrate this, we select in the previous example C 0 = 0.01. Then, as shown in Fig. 2c, the initial guess in this case has only one defect with a very small area. However, after a few iterations, the second object is detected, as can be seen in Fig. 2d, e. The final reconstruction, shown in Fig. 2g is qualitatively equal to the one obtained for C 0 = 0.05. Comparing the decay of the cost function in both situations (compare Figs. 1h and 2h), we see that now the decay is less steep during the first iterations, while at the 9-th iteration a sudden decrease occurs due to the identification of the second scatterer. This example shows, therefore, that although a priori our method seems to highly depend on the selected constant and to need a careful calibration, in practice, it is very robust in terms of it and the final reconstruction is almost insensitive to this parameter (at least for a wide range of C 0 ). Unless otherwise stated, we will keep C 0 = 0.01 to be more conservative.
We illustrate in Fig. 3 the performance of the method when considering objects of different sizes. By reducing the size of the ball, as expected, the largest negative values of the TD are attained inside the object bigger in size, and the initial guess has only one defect when considering the value C 0 = 0.01. However, after a few iterations, the method is able to capture the presence of the second object and the final approximation at the iteration 25 is satisfactory for both objects. Our results in Figs. 5 and 6 evidence that during the first iterations only the objects located closer to the observa-tions points are identified, but during the iterative procedure it is possible to capture and satisfactorily approximate the objects far from them. The reconstructions are slightly worse than in the full-aperture case, specially for the object located further from the observation points and for the parts of both objects in the opposite part to them. We have checked that the same happens when considering the counterpart examples with incident directions corresponding to a half-circle and full-aperture observation points. If we consider limited aperture for both incident directions and observation points while reducing also their numbers (in Fig. 6, we have N inc = 18, N obs = 32, while in the original Fig. 2, these numbers are N inc = 35, N obs = 64), the quality of the reconstructions worsens, as can be seen in Fig. 6. As expected, the reconstruction for the object closer to the observation points is sharper than for the other one, which underestimates its size. The obtained results are in accordance to limited aperture experiments obtained for other model problems by either using topological derivative-based methods (see [2,14,15,34,42]) or by other techniques like the (generalized) linear sampling, factorization, or MUSIC algorithms (see [3,4,6,37]) In all the previous examples with full-aperture configurations, the correct number of defects, as well as their shape are accurately reconstructed. However, we have observed that depending on the wave number and/or the impedance parameter values, reconstructions could be not precise due  to the stagnation of the method before reaching the discrepancy principle stopping criteria. This is illustrated in the next example, where the geometry is the same as in Figs. 1 and 2, but the impedance parameter is now λ = 1 + i instead of λ = 10 + i. The algorithm stopped at the 9th iteration, since it was not able to further increase the size of the current objects (compare Fig. 7d and e). It is clear also, by observing the values of the cost functional in Fig. 7f that the value of the cost function is still rather large in comparison with the values reached in the previous examples, in which the algorithm stopped due to the discrepancy principle. Therefore, we conclude that our method not always provides accurate reconstructions. In cases where the method stagnates, we propose to use our approximation for some other iterative method to further minimize the cost functional. For instance, it could be used in combination with parametric derivatives through iteratively regularized Gauss-Newton methods (IRGNM) [11,30,32]. Roughly speaking, the method consists in computing first the TD to have an initial guess for the IRGNM. Standard IRGNM computes updates now of the parameterization of the shapes of the objects. One of the main drawbacks of standard IRGNM is that the number of defects has to be known, since the initial guess must have the correct number of defects. However, TDs can be combined with IRGNM in a more sophisticated way, by improving the shape of the initial guesses as in standard IRGNM, but when stagnation of the cost functional occurs, then the Gauss-Newton iterations automatically stop and a new TD computation is performed for the current approximate domains (using therefore formula (37)) to identify new possible objects that were left out at the first TD iteration. In case new objects appear, then the Gauss-Newton iterations are restarted. For more details about the combined method, we refer to our recent paper [11].
In the next two examples, we illustrate the performance of the TD/IRGNM combination. Once the indicator function provides the location c * δ and the boundary parametrization q * δ of the boundary * of a new particle, the latter is optimized by inverting the operator formulation where F r denotes the boundary to near-field operator that maps a parametrization of the boundary onto the scattered field u scat r of the solution to the problem (1)-(3) for the incident wave u inc r . We solve (42) in a set of admissible parameterizations V which forms an open subset of a Hilbert space X = H s ( ), for some s > 1.
Then, the iterations of the IRGN method can be computed by solving the nonlinear least-squares problem [30] q k+1 Here, the operator F is Fréchet differentiable and its firstorder Fréchet derivative is given by the solution of a new impedance exterior problem [28, Theorem 2.5] that can be solved with indirect boundary integral equation methods too, namely the adjoint form of (10). We consider a 3D configuration with two obstacles of comparable sizes: a sphere and a rounded tetrahedron (see Fig. 8a) with impedance parameter λ = 1 + i, namely, again with a low impedance contrast. Mimicking the previous twodimensional examples, we consider 35 incident directions for the wave number κ = 2 and 64 observation points on the sphere of radius 10. A noise level of a 5% is added to the synthetic data. By computing the TD in the observation cube [−3, 3] 3 we obtain the color maps in Fig. 8b, c, d, which represent the values of the TD at some parallel planes to the Cartesian planes x = 0, y = 0 and z = 0, respectively. It is clear that the largest negative values concentrate inside both defects. If we select the constant C 0 = 0.1 to define our initial guess, we get the sphere-like defects shown in Fig. 8e, which roughly have the same size and shape. The iterative TD-based algorithm only improves the size approximation, but it does not provide satisfactory results (omitted for the sake of brevity) and stagnates before being able to distinguish the shapes. However, using the initial guess in Fig. 8e for the IRGNM we obtain an accurate approximation after only 15 IRGNM iterations. Some intermediate and the final reconstructions are shown in Fig. 8f, g, h, i. We observe that in the first iterations, the method captures the approximate volumes, while in the subsequent iterations, shapes are better identified. The values of the cost function versus the number of iterations are presented in Fig. 8j.
Finally, in our last example, we reduce the size of the sphere (see Fig. 9a) while keeping the remaining parameters unchanged, and repeat the experiment. The values of the TD at some parallel planes to the Cartesian planes x = 0, y = 0 and z = 0 are represented in Fig. 9b, c, d. Due to the small size of the ball in comparison with the size of the tetrahedron, the largest negative values of the TD concentrate now inside the tetrahedron. For a wide range of constants C 0 , the initial guess only consists in one object. Keeping the same value as in the previous experiment, C 0 = 0.1, we find the initial guess shown in Fig. 9e. In this case, again, the iterative TD-based algorithm does not converge to a satisfactory set of objects (figure omitted for brevity). The application of the IRGNM by using this initial guess somehow improves the reconstruction of the tetrahedron, but the method stagnates without reaching the discrepancy principle stopping criteria (see the intermediate reconstructions in Fig. 9f, g, h and the values of the cost function in the first iterations, collected in Fig. 9p). Then, the combined TD/IRGNM automatically performs a new TD computation where as can be observed in Fig. 9i, j, k the second object is clearly detected. In this figure, the largest negative values concentrate inside the small spherical object, while no new regions close to the current approximate domain (in cyan color) appear. The IRGNM iterations restart now with an approximation that has the correct number of defects which are reasonably well approximated in size. After a few iterations, we observe an improvement in the shape identification (panel (n)), and at the 16th iteration, both defects are accurately recovered, as can be seen in panel 9(o).

Conclusions
We have derived a closed-form formula for the topological derivative (TD) of the least-squares misfit functional in presence of a given set of coated impedance scatterers, which mimics the one associated with the absence of such defects (already found by us in our previous paper [35]), but whose derivation is different and requires the use of more involved mathematical tools.
This ready-to-implement formula is then used to propose an iterative method to improve the initial guesses obtained by a TD computation in the free space. Our numerical experiments evidence that in case of medium/high complex impedance contrast, our iterative method fully based on sequential TD computations can highly improve onestep approximations, not only in finding the correct number of objects, but also in improving their shape identification (even in some demanding cases dealing with limited aperture data). However, when the complex impedance contrast is low, our method could stagnate without finding a satisfactory reconstruction. The formula derived in the present paper can be used then in hybrid algorithms not only based on TD computations. We have illustrated this fact by considering a combined TD-Gauss-Newton method proposed in [11], finding accurate approximations in demanding cases with low impedance contrast and objects of different sizes. Some other hybrid methods could be used in combination with our formula, as for example TD-level set approaches [7,23,29].
Conceptually, we can extend the underlying idea of this paper to more complicated elliptic systems with high-order generalized impedance boundary conditions and to steady (visco)-elasticity equations. This will be the subject of forthcoming works.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.