Multifrequency Topological Derivative Approach to Inverse Scattering Problems in Attenuating Media

Abstract

Detecting objects hidden in a medium is an inverse problem. Given data recorded at detectors when sources emit waves that interact with the medium, we aim to find objects that would generate similar data in the presence of the same waves. In opposition, the associated forward problem describes the evolution of the waves in the presence of known objects. This gives a symmetry relation: if the true objects (the desired solution of the inverse problem) were considered for solving the forward problem, then the recorded data should be returned. In this paper, we develop a topological derivative-based multifrequency iterative algorithm to reconstruct objects buried in attenuating media with limited aperture data. We demonstrate the method on half-space configurations, which can be related to problems set in the whole space by symmetry. One-step implementations of the algorithm provide initial approximations, which are improved in a few iterations. We can locate object components of sizes smaller than the main components, or buried deeper inside. However, attenuation effects hinder object detection depending on the size and depth for fixed ranges of frequencies.

1. Introduction

Most inverse algorithms for shape reconstruction using acoustic data assume that the attenuation coefficient is negligible. However, in some applications, attenuation effects cannot be discarded. For example, in medical studies of biological tissues, the attenuation coefficient provides important information for diagnosis. It often differs noticeably from normal tissue to damaged tissue, while other acoustical parameters such as the density or the wave speed might vary only slightly [1]. Similarly, sound attenuation in subbottom tomography for the seabed plays an important role (see [2] and references therein). Attenuation effects are also important in photoacoustic imaging [3] and in helioseismology models [4], to mention a few. Seismic waves typically decrease in amplitude due to spherical spreading as well as to loss mechanisms in the rocks, due to mechanical or other causes.

Attenuation complicates the detection of buried scatterers to a great extent. Furthermore, in many applications, data can only be acquired in a limited part of the sounded region, as in seismic tomography or marine acoustics, where taking measurements around the sounded region is impossible. A graphical illustration of such situations is depicted in Figure 1.

A wide variety of algorithms have been proposed to address inverse scattering problems, that is, inverse problems in which the goal is to detect localized objects: linear approximations such as the Kirchhoff or physical optics approximation [5] and the Born approximation [6], linear sampling techniques [7,8], factorization methods [9], backpropagation principles [10], modified gradient methods [11], or single shot methods [12] to quote a few. General shape optimization strategies [13] have inspired many approaches, based, in particular, on level-set methods [14,15] or topological derivative techniques [16,17] (see [18,19] for reviews of recent applications). Here, we show that multifrequency topological derivative-based methods [20,21,22,23] have the potential of providing reconstructions of buried objects in attenuation regimes, within some size and depth limits. We propose a tool based on heuristic ideas, as most of the methods that are based on topological derivatives. For rigorous mathematical justifications of this kind of approach in some specific situations, we refer to [24,25,26].

The paper is organized as follows. Section 2 formulates the forward and inverse scattering problems for penetrable scatterers when both the exterior medium and the objects have nonzero attenuation coefficients. The problem is set in a half-plane (space), but with a symmetry argument, it can be reformulated as a problem in the whole plane (space). Section 3 explains the one-step topological derivative-based method with fixed frequencies, while Section 4 introduces the multifrequency approach. Results are further improved in Section 5, implementing an iterative topological derivative-based multifrequency scheme. Section 6 explains how to extend our method for the detection of sound-soft and sound-hard objects immersed in an attenuating medium. Finally, Section 7 summarizes our conclusions.

2. Forward and Inverse Problems

2.1. Forward Problem

When the incident wave is time harmonic, the resulting attenuated wave field is a time-harmonic solution $U(\mathbf{x},t)=\mathrm{Re}\left[e^{-\mathit{ı}\omega t}u\left(\mathbf{x}\right)\right]$ (where $\omega >0$ is the excitation frequency), of the damped wave equation

$$\Delta U-\frac{1}{c^2}{U}_{tt}-a{U}_t=F.$$

Therefore, the amplitude u is a solution to

$$\Delta u+\left(\frac{{\omega }^2}{c^2}+\mathit{ı}a\omega \right)u=f,$$

where $c>0$ is the wave speed and $a\ge 0$ is the attenuation coefficient. Here the function F (and its time-harmonic amplitude f in the case of time-harmonic excitations) depends on the selected type of incident wave. In the case of time-harmonic planar excitations, $f=0$, while a point source located at a point ${\mathbf{x}}_0$ is modeled by the Dirac delta function at ${\mathbf{x}}_0$, $f={\delta }_{{\mathbf{x}}_0}$.

In the setup depicted in Figure 1, we assume that the components of $\Omega$ are penetrable obstacles defined by the wave speed $c_i>0$, the attenuation coefficient $a_i\ge 0$, and the density ${\rho }_i>0$. The exterior medium is identified with the half-plane (space) ${\mathbb{R}}_{-}^d=\{(x_1,\dots ,x_d)\in {\mathbb{R}}^d,x_d<0\}$, where $d=2$ or $d=3$, and described by the exterior parameters $c_e,a_e,{\rho }_e$. The upper half-plane (space) is assumed to be filled by air, which will be modeled at the upper boundary $\mathrm{\Pi }:=\{(x_1,\dots ,x_d),x_d=0\}$ by the homogeneous Neumann condition

$${\partial }_{x_d}u=0,\mathrm{on} \mathrm{\Pi }.$$

(1)

To simplify notations, we introduce the positive ratio

$$\gamma =\frac{{\rho }_i{c}_i^2}{{\rho }_e{c}_e^2},$$

(2)

and the complex wave numbers

$$k_e\left(\omega \right)=\sqrt{\frac{{\omega }^2}{c_e^2}+\mathit{ı}{a}_e\omega },k_i\left(\omega \right)=\sqrt{\frac{{\omega }^2}{c_i^2}+\mathit{ı}{a}_i\omega },$$

(3)

where the square roots are selected with non-negative imaginary parts.

Incident fields will be generated at points ${\mathbf{x}}_0\in \mathrm{\Pi }$ and mathematically described by the fundamental solution of the Helmholtz equation with wave number $k_e\left(\omega \right)$:

$$u_{{\mathbf{x}}_0,\omega }^{inc}\left(\mathbf{x}\right)=G(k_e\left(\omega \right),\mathbf{x}-{\mathbf{x}}_0)$$

(4)

where

$$G(k,\mathbf{z})=\left\{\begin{array}{cc}{\displaystyle \frac{i}{4}{H}_0^{\left(1\right)}\left(k\right|\mathbf{z}\left|\right),}\hfill & \mathrm{i}\mathrm{f} d=2,\hfill \\ {\displaystyle \frac{e^{\mathit{ı}k\left|\mathbf{z}\right|}}{4\pi \left|\mathbf{z}\right|},}\hfill & \mathrm{i}\mathrm{f} d=3,\hfill \end{array}\right.$$

(5)

$H_0^{\left(1\right)}$ being the Hankel function of the first kind and order zero [27]. Notice that since ${\mathbf{x}}_0\in \mathrm{\Pi }$, it follows that $u_{{\mathbf{x}}_0,\omega }^{inc}$ satisfies the condition (1). This incident field generates a scattered wave $u_{{\mathbf{x}}_0,\omega }^{scat}$ in ${\mathbb{R}}_{-}^d\backslash \overline{\mathrm{\Omega }}$ and a transmitted wave $u_{{\mathbf{x}}_0,\omega }^{tr}$ in $\mathrm{\Omega }$. The total field

$$u_{{\mathbf{x}}_0,\omega }:=\left\{\begin{array}{cc}{u}_{{\mathbf{x}}_0,\omega }^{inc}+u_{{\mathbf{x}}_0,\omega }^{scat},\hfill & \mathrm{in}{ \mathbb{R}}_{-}^d\backslash \overline{\mathrm{\Omega }},\hfill \\ u_{{\mathbf{x}}_0,\omega }^{tr},\hfill & \mathrm{in} \mathrm{\Omega },\hfill \end{array}\right.$$

is a solution of

$$\left\{\begin{array}{c}\Delta u_{{\mathbf{x}}_0,\omega }+k_e{\left(\omega \right)}^2{u}_{{\mathbf{x}}_0,\omega }={\delta }_{{\mathbf{x}}_0},\mathrm{i}\mathrm{n} {\mathbb{R}}_{-}^d\backslash \overline{\mathrm{\Omega }},\hfill \\ \Delta u_{{\mathbf{x}}_0,\omega }+k_i{\left(\omega \right)}^2{u}_{{\mathbf{x}}_0,\omega }=0,\mathrm{i}\mathrm{n} \mathrm{\Omega },\hfill \\ u_{{\mathbf{x}}_0,\omega }^{-}=u_{{\mathbf{x}}_0,\omega }^{+},\gamma {\partial }_{\mathbf{n}}{u}_{{\mathbf{x}}_0,\omega }^{-}={\partial }_{\mathbf{n}}{u}_{{\mathbf{x}}_0,\omega }^{+},\mathrm{o}\mathrm{n} \partial \mathrm{\Omega },\hfill \\ {\partial }_{x_d}{u}_{{\mathbf{x}}_0,\omega }=0\mathrm{o}\mathrm{n} \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r(u_{{\mathbf{x}}_0,\omega }-u_{{\mathbf{x}}_0,\omega }^{inc})-\mathit{ı}{k}_e\left(\omega \right)(u_{{\mathbf{x}}_0,\omega }-u_{{\mathbf{x}}_0,\omega }^{inc}))=0.\hfill \end{array}\right.$$

(6)

The superscripts + and − denote limit values from outside and inside the objects, respectively. Here, ${\partial }_{\mathbf{n}}$ stands for normal derivatives, $\mathbf{n}$ being the normal unit vector pointing inside the objects, whereas ${\partial }_r$ denotes the radial derivative. For our particular type of complex wave-numbers, problem (6) has a unique solution [28] and it can be solved efficiently by means of boundary element techniques [29,30]. In Reference [31] several direct, indirect, and mixed formulations are analyzed and compared for a thermal problem, which has exactly the same structure as (6) but where wave numbers are related to the time-harmonic heat equation and take the form $k\left(\omega \right)=\sqrt{\mathit{ı}b\omega }$ with $b>0$ instead of those in (3). However, for the wave numbers of the form (3), such formulations are also valid. Results in [28,31] are established by using a symmetry argument, where the problem in the half-plane (space) is transformed into an equivalent one in ${\mathbb{R}}^d$: the condition ${\partial }_{x_d}u=0$ on $\mathrm{\Pi }$ disappears and the objects are $\tilde{\mathrm{\Omega }}=\mathrm{\Omega }\cup {\mathrm{\Omega }}_{\mathrm{sym}}$, where $\mathrm{\Omega }\subset {\mathbb{R}}_{-}^d$ are the original objects and ${\mathrm{\Omega }}_{\mathrm{sym}}$ are their symmetrical (with respect to $\mathrm{\Pi }$) ones, namely, ${\mathrm{\Omega }}_{\mathrm{sym}}=\{(x_1,\dots ,x_d)\in {\mathbb{R}}^d,(x_1,\dots ,-x_d)\in \mathrm{\Omega }\}$. Boundary element methods proposed in [28,31] are based on potentials and operators defined on the symmetric boundary $\partial \tilde{\mathrm{\Omega }}=\partial \mathrm{\Omega }\cup \partial {\mathrm{\Omega }}_{\mathrm{sym}}$. In the forthcoming examples for the two-dimensional case, we opt for a fully discrete version based on trigonometric polynomials of the direct formulation described in Section 8 in [31] (which indeed is a generalization of the original method proposed in [32]), where the unknowns are the interior and the exterior Cauchy data of the solution of (6) at the boundary of the objects. Alternatively, we could have used the symmetric formulation proposed in [33].

For simplicity, in the initial simulations, we fix $\gamma =1$, $c_e=1$ and $c_i=4$ in the configuration sketched in Figure 2. We will comment on the case $\gamma \ne 1$ in the final simulations. Figure 3 and Figure 4 illustrate the effect of the attenuation coefficients $a_e$ and $a_i$ on the incident and scattered fields for two different frequencies $\omega$. For a fixed value of the interior attenuation parameter $a_i$, we observe that the intensity of the scattered waves diminishes as the exterior parameter $a_e$ increases and its maxima are smoothed, while the support of the incident wave becomes smaller. On the other hand, for a fixed value of $a_e$, the effect of increasing $a_i$ does not have a noticeable qualitative impact on the scattered wave, and of course, the incident wave does not change since it is independent of this parameter. As we increase the frequency $\omega$, the scattered wave becomes more oscillatory. Notice that the magnitude of the scattered wave depends on the frequency, which will have to be taken into account when combining multifrequency data for shape reconstruction algorithms.

2.2. Inverse Problem

Given objects $\mathrm{\Omega }$ and the incident wave $u_{{\mathbf{x}}_0,\omega }^{inc}$, the solution of (6) evaluated at the receptors ${\mathbf{z}}_r$, $r=1,\dots ,N_{receivers},$ is expected to agree with the measured data $u_{{\mathbf{x}}_0,\omega }^{meas,r}$ except for measurement errors and noise. When we have measurements of the total wave field at a given frequency $\omega$ for different source point locations ${\mathbf{x}}_s$, $s=1,\dots ,N_{emitters}$, a constrained optimization reformulation of the inverse problem seeks minimizers $\mathrm{\Omega }$ of the shape functional

$$\begin{array}{c}\hfill J_{\omega }({\mathbb{R}}_{-}^d\backslash \overline{\mathrm{\Omega }})=\frac{1}{2}\sum _{s=1}^{N_{emitters}}\sum _{r=1}^{N_{receivers}}{|u_{{\mathbf{x}}_s,\omega }^{\mathrm{\Omega }}\left({\mathbf{z}}_r\right)-u_{{\mathbf{x}}_s,\omega }^{meas,r}|}^2,\end{array}$$

(7)

where $u_{{\mathbf{x}}_s,\omega }^{\mathrm{\Omega }}$ denotes the solution of (6) with objects $\mathrm{\Omega }$ for known parameters $c_i,c_e,a_i,a_e$ and $\gamma$ when the incident wave is $u_{{\mathbf{x}}_s,\omega }^{inc}$ (defined in (4)).

For our two-dimensional experiments, data will be generated synthetically by solving problem (6) for the true defects $\mathrm{\Omega }$ via the direct formulation described in Section 8 in [31] to obtain the solution at the receivers, $u_{{\mathbf{x}}_s,\omega }^{\mathrm{\Omega }}\left({\mathbf{z}}_r\right)$, $r=1,\dots ,N_{receivers}$. Then, the associated scattered wave $u_{{\mathbf{x}}_s,\omega }^{scat,r}=u_{{\mathbf{x}}_s,\omega }^{\mathrm{\Omega }}\left({\mathbf{z}}_r\right)-u_{{\mathbf{x}}_s,\omega }^{inc}\left({\mathbf{z}}_r\right)$ is corrupted by adding a random number of size $\delta$ to generate corrupted scattered values ${\tilde{u}}_{{\mathbf{x}}_s,\omega }^{scat,r}$ satisfying

$$\frac{\sqrt{{\sum }_{r=1}^{N_{receivers}}{|{\tilde{u}}_{{\mathbf{x}}_s,\omega }^{scat,r}-u_{{\mathbf{x}}_s,\omega }^{scat,r}|}^2}}{\sqrt{{\sum }_{r=1}^{N_{receivers}}{\left|u_{{\mathbf{x}}_s,\omega }^{scat,r}\right|}^2}}\le \delta .$$

(8)

Finally, the synthetic measured data are defined by

$$u_{{\mathbf{x}}_s,\omega }^{meas,r}=u_{{\mathbf{x}}_s,\omega }^{inc}\left({\mathbf{z}}_r\right)+{\tilde{u}}_{{\mathbf{x}}_s,\omega }^{scat,r}.$$

(9)

In the forthcoming experiments, we set a one percent relative noise, namely, $\delta =0.01$. However, we have already checked by numerical tests that results are rather insensitive to $\delta$ and results are almost the same if we consider a ten percent relative noise (i.e., for $\delta =0.1$). To further avoid inverse crimes, we use a different boundary integral method to solve the direct problems of the form (6) that will appear during the iterative algorithm. Being more precise, we use the indirect formulation described in Section 3 in [31] (see also [28]).

In the next section, we explain how to generate first guesses of the objects by topological derivative methods.

3. Topological Derivative-Based Approach

The concept of the topological derivative was first introduced in the 1990s in [13,34,35] as a tool for optimal shape design. The topological derivative of a shape functional $J\left(\mathcal{R}\right)=J\left(u_{\mathcal{R}}\right)$, where $\mathcal{R}\subset {\mathbb{R}}^d$ and $u_{\mathcal{R}}$ is the solution of a given boundary value problem defined in $\mathcal{R}$, is a scalar field that measures the sensitivity of such functional to infinitesimal perturbations at each point $\mathbf{x}\in \mathcal{R}$. It provides expansions of the form

$$\begin{array}{c}\hfill J(\mathcal{R}\backslash \overline{B_{\epsilon }\left(\mathbf{x}\right)})=J\left(\mathcal{R}\right)+f\left(\epsilon \right)D_T(\mathcal{R},\mathbf{x})+o\left(f\left(\epsilon \right)\right),\end{array}$$

(10)

where $B_{\epsilon }\left(\mathbf{x}\right)$ are balls of radius $\epsilon$ centered at $\mathbf{x}$ and the positive function f satisfies $f\left(\epsilon \right)\to 0$ as $\epsilon \to 0$. The field $D_T(\mathcal{R},\mathbf{x})$ is the topological derivative of the shape functional at the point $\mathbf{x}\in \mathcal{R}$. The function f depends on the boundary conditions enforced in the boundary value problem. In our case, $\mathcal{R}={\mathbb{R}}_{-}^d$ and f will be the measure of the ball, that is, $f\left(\epsilon \right)=\pi {\epsilon }^2$ when $d=2$ and $f\left(\epsilon \right)=\frac{4}{3}\pi {\epsilon }^3$ when $d=3$. Whenever $D_T(\mathcal{R},\mathbf{x})$ takes large negative values, the expansion (10) implies that the shape functional decreases by placing an infinitesimal scatterer at $\mathbf{x}$. This idea will be exploited in the definition of our reconstruction algorithm.

The general definition of the topological derivative allows for any kind of infinitesimal shape, not necessarily a ball. However, to obtain the associated function f and the derivation of closed-form formulas in this case is more involved, while it is not clear that it would produce improvements in practice, as observed in [36], where a comparison between results for infinitesimal balls and ellipsoids at different orientations is presented for sound-soft and sound-hard objects immersed in a non-attenuating media.

A closed-form expression of the topological derivative of a shape functional similar to (7), involving complex wave numbers, but related to time-harmonic thermal waves in a half-plane is established in [37]. In an acoustic setup in ${\mathbb{R}}^3$ without attenuation, it was established in [36]. In a similar way, we obtain the following result:

Theorem 1.

For any $\mathbf{x}\in {\mathbb{R}}_{-}^d$, the topological derivative of the shape functional (7) is

$$\begin{array}{ccc}\hfill D_T^{\omega }({\mathbb{R}}_{-}^d,\mathbf{x})& =& \sum _{s=1}^{N_{emitters}}Re[\frac{d(1-\gamma )}{d-1+\gamma }\nabla u_{{\mathbf{x}}_s,\omega }^{\varnothing }\left(\mathbf{x}\right)\nabla \overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}\left(\mathbf{x}\right)\hfill \\ & & +\left(\gamma k_i{\left(\omega \right)}^2-k_e{\left(\omega \right)}^2\right)u_{{\mathbf{x}}_s,\omega }^{\varnothing }\left(\mathbf{x}\right)\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}\left(\mathbf{x}\right)]\hfill \end{array}$$

(11)

where $u_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $w_{{\mathbf{x}}_s,\omega }^{\varnothing }$ solve the forward and adjoint problems with $\mathrm{\Omega }=\varnothing$:

$$\left\{\begin{array}{c}\Delta u_{{\mathbf{x}}_s,\omega }^{\varnothing }+k_e{\left(\omega \right)}^2{u}_{{\mathbf{x}}_s,\omega }^{\varnothing }={\delta }_{{\mathbf{x}}_s},in {\mathbb{R}}_{-}^d,\hfill \\ {\partial }_{x_d}{u}_{{\mathbf{x}}_s,\omega }^{\varnothing }=0on \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r(u_{{\mathbf{x}}_s,\omega }^{\varnothing }-u_{{\mathbf{x}}_s,\omega }^{inc})-\mathit{ı}{k}_e\left(\omega \right)(u_{{\mathbf{x}}_s,\omega }^{\varnothing }-u_{{\mathbf{x}}_s,\omega }^{inc}))=0,\hfill \end{array}\right.$$

(12)

and

$$\left\{\begin{array}{c}\Delta \overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}+k_e{\left(\omega \right)}^2\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}=\sum _{r=1}^{N_{receivers}}\left(\overline{u_{{\mathbf{x}}_s,\omega }^{meas,r}-u_{{\mathbf{x}}_s,\omega }^{\varnothing }}\right){\delta }_{{\mathbf{z}}_r},in {\mathbb{R}}_{-}^d,\hfill \\ {\partial }_{x_d}\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}=0on \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}-\mathit{ı}{k}_e\left(\omega \right)\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }})=0.\hfill \end{array}\right.$$

(13)

Both the forward and the adjoint problems can be solved in closed-form (recall that since both emitters and receivers are located on the upper boundary $\mathrm{\Pi }$, the fundamental solution $G(k_e\left(\omega \right),\mathbf{x}-{\mathbf{x}}_{*})$ satisfies the boundary condition (1) at $\mathbf{x}\in \mathrm{\Pi }$ when ${\mathbf{x}}_{*}$ is either a receiver or an emitter):

$$\begin{array}{ccc}\hfill u_{{\mathbf{x}}_s,\omega }^{\varnothing }\left(\mathbf{x}\right)& =& u_{{\mathbf{x}}_s,\omega }^{inc}\left(\mathbf{x}\right)=G(k_e\left(\omega \right),\mathbf{x}-{\mathbf{x}}_s),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill \overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}\left(\mathbf{x}\right)& =& \sum _{r=1}^{N_{receivers}}G(k_e\left(\omega \right),\mathbf{x}-{\mathbf{x}}_r)\left(\overline{u_{{\mathbf{x}}_s,\omega }^{\varnothing }\left({\mathbf{x}}_r\right)-u_{{\mathbf{x}}_s,\omega }^{meas,r}}\right).\hfill \end{array}$$

(15)

Formula (11) can be directly evaluated in an observation region ${\mathcal{R}}_{obs}$ when the recorded data $u_{{\mathbf{x}}_s,\omega }^{meas,r}$, the incident waves $u_{{\mathbf{x}}_s,\omega }^{inc}$ and the material properties are known. In case the properties of the objects $c_i,{\rho }_i,a_i$ are not known, we can implement iterative procedures to calculate both shapes and material properties following [38]. We will explain how to iterate with respect to the objects using several frequencies in Section 5.

A simpler situation where $\mathcal{R}={\mathbb{R}}^d$ in the absence of attenuation effects was studied in [36,39]. The formula for the topological derivative is analogous to that in Theorem 1, but the wave numbers $k_e\left(\omega \right)$ and $k_i\left(\omega \right)$ are real. The forward and adjoint problems are similar, but there is no boundary condition on $\mathrm{\Pi }$ and the forward solution satisfies $\Delta u_{{\mathbf{x}}_s,\omega }^{\varnothing }+k_e{\left(\omega \right)}^2{u}_{{\mathbf{x}}_s,\omega }^{\varnothing }={\delta }_{{\mathbf{x}}_s}$ in ${\mathbb{R}}^d$, while the adjoint solution solves $\Delta \overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}+k_e{\left(\omega \right)}^2\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}={\sum }_{r=1}^{N_{receivers}}\left(\overline{u_{{\mathbf{x}}_s,\omega }^{meas,r}-u_{{\mathbf{x}}_s,\omega }^{\varnothing }}\right){\delta }_{{\mathbf{z}}_r}$, in ${\mathbb{R}}^d$. In this case, when emitters and receivers are distributed around the observation region, the largest negative values of the topological derivative clearly locate the objects for different frequencies, as it happens in Figure 5. The physical parameters selected for this simulation are $c_e=1$, $c_i=4$ and $\gamma =1$. Notice that the shape and orientation of the two objects are better described when we increase the frequency, in spite of oscillations.

Let us turn now to the original problem set in ${\mathbb{R}}_{-}^d$ for the limited aperture configurations under study, see Figure 1. Locating buried scatterers with limited-aperture data is an issue of practical interest that has been the object of intensive research [36,40]. Placing emitters and receivers on the upper boundary, and keeping the same physical parameters $a_e=a_i=0$, $c_e=1$, $c_i=4$ and $\gamma =1$, the information provided by topological derivatives is less clear, as illustrated in Figure 6. Now, we observe stripes alternating positive and negative values of the topological derivative. Nevertheless, the largest negative values are attained around the two true objects, though resolution for the object placed farthest from the emitters worsens, especially as the frequency increases. In comparison with the experiments presented in Figure 5, we observe that while frequencies $\omega \ge 8$ are the best for shape identification when emitters and receivers are located around the objects, in our setting with emitters and receivers on $\mathrm{\Pi }$, the use of frequencies higher than $\omega =4$ promotes the appearance of spurious regions. Only frequencies below $\omega =8$ will be considered in the forthcoming experiments. If we increase the number of emitters and/or receivers but we keep their location on $\mathrm{\Pi }$, reconstructions do not improve.

In the presence of attenuation, we lose even more information, see Figure 7 where $a_e=a_i=0.5$. The topological derivative still detects the closest object, though its orientation is uncertain. We fail to detect the distant object, while a number of spurious minima appear. Depending on the frequency, we may end up with completely wrong reconstructions. For instance, for $\omega =4$, we conjecture the presence of a close object at about $y=-0.5$, whereas for $\omega =6$, we seem to have three defects. If we increase the interior attenuation coefficient to $a_i=1.5$ and keep $a_e=0.5$, the topological derivatives are qualitatively equal to those presented in Figure 7 and therefore plots are omitted for brevity.

Further increasing the exterior attenuation coefficient to $a_e=1.5$ while keeping the interior one unchanged $a_i=0.5$ (or increasing it to $a_i=1.5$), reconstructions worsen. As shown in Figure 8, only the lowest frequency can track the closest object. The true objects are usually placed in regions where the topological derivative takes negative values, but the largest negative values are attained in variable regions as the frequency changes. Furthermore, there are regions where large negative values for one frequency become large positive values for another (compare panels (c) and (f), for instance). These remarks suggest the possibility of combining information from different frequencies to enhance the correct contributions provided by all of them and to cancel the variable, spurious details, which we analyze in the next section.

4. Multifrequency Method

To incorporate the contributions from several frequencies ${\omega }_1,\dots ,{\omega }_{N_{freq}}$, we consider a new shape functional defined as a linear combination of monochromatic functionals (7), namely, a function of the form [22]:

$$J({\mathbb{R}}_{-}^d\backslash \overline{\mathrm{\Omega }})=\frac{1}{2}\sum _{\ell =1}^{N_{freq}}{P}_{{\omega }_{\ell }}{J}_{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{\mathrm{\Omega }})=\frac{1}{2}\sum _{\ell =1}^{N_{frec}}\sum _{s=1}^{N_{emitters}}\sum _{r=1}^{N_{receivers}}{P}_{{\omega }_{\ell }}{|u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{\mathrm{\Omega }}\left({\mathbf{z}}_r\right)-u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{meas,r}|}^2,$$

(16)

where $u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{\mathrm{\Omega }}$ is the solution of (6) with frequency $\omega ={\omega }_{\ell }$ and incident wave $u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{inc}$ of the form (4). Here, $P_{{\omega }_{\ell }}>0$ are weights to be selected. By linearity, it follows that the topological derivative is a linear combination of the monochromatic topological derivatives with $\omega ={\omega }_{\ell }$ given in (11), i.e.,

$$\begin{array}{c}\hfill D_T({\mathbb{R}}_{-}^d,\mathbf{x})=\sum _{\ell =1}^{N_{frec}}{P}_{{\omega }_{\ell }}{D}_T^{{\omega }_{\ell }}({\mathbb{R}}_{-}^d,\mathbf{x}),\mathbf{x}\in {\mathbb{R}}_{-}^d.\end{array}$$

(17)

As in [22], we select the weights $P_{{\omega }_{\ell }}$ after computing each individual derivative, selecting

$$P_{{\omega }_{\ell }}=-\underset{\mathbf{y}\in {\mathcal{R}}_{obs}}{min}{D}_T^{{\omega }_{\ell }}({\mathbb{R}}_{-}^d,\mathbf{y}),$$

(18)

where ${\mathcal{R}}_{obs}$ is the observation region, that is, the region where topological derivatives are evaluated. This choice balances the contribution of all frequencies since it guarantees that each term in the sum (17) satisfies ${min}_{\mathbf{x}\in {\mathcal{R}}_{obs}}{P}_{{\omega }_{\ell }}{D}_T^{{\omega }_{\ell }}({\mathbb{R}}_{-}^d,\mathbf{x})=-1$. Notice that depending on the selected frequency, $P_{{\omega }_{\ell }}$ has a different order of magnitude (compare panels (a) and (f) in Figure 6, Figure 7 and Figure 8).

Figure 9 illustrates the results obtained combining the six single frequency topological derivatives for the same ellipses ($c_i=4,a_i=0.5,\gamma =1$) when the exterior velocity is $c_e=1$ and $a_e=0$ (no exterior attenuation), $a_e=0.5$ (moderate attenuation) and $a_e=1.5$ (high attenuation). For the attenuation coefficients $a_e=0.5$ and $a_e=1.5$, we have combined the six topological derivatives displayed in Figure 7 and Figure 8, while for $a_e=0$ the individual ones are qualitatively identical to those in Figure 6, which corresponds to $a_i=0$ instead of $a_i=0.5$. When $a_e=0$, we clearly distinguish two elliptical defects, obtaining good approximations for their size, shape and orientation. For $a_e=0.5$, the object closest to emitters/receivers is reasonably reproduced, similarly to the case $a_e=0$, but the orientation slightly worsens. We conjecture the presence of another object in the second region where the multifrequency topological derivative becomes negative, and no spurious negative regions where no true objects are located appear. Increasing the attenuation in the medium, the information on the orientation of the closest object worsens and the other object cannot really be guessed.

For a quantitative visualization, Figure 10 illustrates the initial approximations ${\mathrm{\Omega }}_{\alpha }$ we would obtain plotting the topological derivative (17) for the example in Figure 9c (which corresponds to the most demanding case $a_e=1.5$) and selecting the points for which the topological derivative falls below a threshold:

$$\begin{array}{c}\hfill \mathrm{\Omega }\approx {\mathrm{\Omega }}_{\alpha }=\{\mathbf{x}\in {\mathcal{R}}_{obs},D_T({\mathbb{R}}_{-}^d,\mathbf{x})\le (1-\alpha )\underset{\mathbf{y}\in {\mathcal{R}}_{obs}}{min}{D}_T({\mathbb{R}}_{-}^d,\mathbf{y})\},\end{array}$$

(19)

where $0<\alpha <1$ is the selected threshold coefficient and ${\mathcal{R}}_{obs}=[-2.5,2.5]\times [-2.5,-0.2]$ is the sounded region. Obviously, the larger $\alpha$, the larger the area. However, notice that reconstructions are robust in terms of this parameter and only a rough calibration is needed. Notice also that ${\mathrm{\Omega }}_{\alpha }$ not only depends on $\alpha$ but also on the selected mesh for ${\mathcal{R}}_{obs}$ (both for the definition of each weight $P_{{\omega }_{\ell }}$, and for the final evaluation of the multifrequency topological derivative). In all our experiments, we have considered a grid of 120 points in the x-coordinate and 60 points in the y-coordinate. We have also checked that the method is very robust in connection with this issue. For a discussion about the selection of optimal meshes for the evaluation of indicator functions, we refer to [41].

These results provide the initial stage to implement iterative corrections, as we explain in Section 5.

5. Multifrequency Iterative Method

In this section, we adapt the iterative algorithm proposed by us in [37,39] for monochromatic measurements to the multifrequency case. It is based upon an expression of the topological derivative when we already have a non-empty guess ${\mathrm{\Omega }}_{app}$ for the scatterers.

Theorem 2.

Assume that ${\mathrm{\Omega }}_{app}$ is an approximation of the true defects Ω. Then, for any $\mathbf{x}\in {\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{app}}$ the topological derivative of the shape functional $J_{\omega }$ defined in (7) is given by

$$\begin{array}{ccc}\hfill D_T^{\omega }({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{app}},\mathbf{x})& =& \sum _{s=1}^{N_{emitters}}Re[\frac{d(1-\gamma )}{d-1+\gamma }\nabla u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}\left(\mathbf{x}\right)\nabla \overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\left(\mathbf{x}\right)\hfill \\ & & +\left(\gamma k_i{\left(\omega \right)}^2-k_e{\left(\omega \right)}^2\right)u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}\left(\mathbf{x}\right)\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\left(\mathbf{x}\right)]\hfill \end{array}$$

(20)

where $u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ and $w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ are solutions of forward and adjoint problems with objects ${\mathrm{\Omega }}_{app}$:

$$\left\{\begin{array}{c}\Delta u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}+k_e{\left(\omega \right)}^2{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}={\delta }_{{\mathbf{x}}_s},in {\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{app}},\hfill \\ \Delta u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}+k_i{\left(\omega \right)}^2{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}=0,in {\mathrm{\Omega }}_{app},\hfill \\ {\left(u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}\right)}^{-}={\left(u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}\right)}^{+},\gamma {\left({\partial }_{\mathbf{n}}{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}\right)}^{-}={\left({\partial }_{\mathbf{n}}{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}\right)}^{+},on \partial {\mathrm{\Omega }}_{app},\hfill \\ {\partial }_{x_d}{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}=0on \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r(u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}-u_{{\mathbf{x}}_s,\omega }^{inc})-\mathit{ı}{k}_e\left(\omega \right)(u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}-u_{{\mathbf{x}}_s,\omega }^{inc}))=0,\hfill \end{array}\right.$$

(21)

and

$$\left\{\begin{array}{c}\Delta \overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}+k_e{\left(\omega \right)}^2\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=\sum _{r=1}^{N_{receivers}}\left(\overline{u_{{\mathbf{x}}_s,\omega }^{meas,r}-u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\right){\delta }_{{\mathbf{z}}_r},in {\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{app}},\hfill \\ \Delta \overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}+k_i{\left(\omega \right)}^2\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=0,in {\mathrm{\Omega }}_{app},\hfill \\ {\left(\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\right)}^{-}={\left(\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\right)}^{+},\gamma {\left({\partial }_{\mathbf{n}}\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\right)}^{-}={\left({\partial }_{\mathbf{n}}\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\right)}^{+},on \partial {\mathrm{\Omega }}_{app},\hfill \\ {\partial }_{x_d}\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=0on \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}-\mathit{ı}{k}_e\left(\omega \right)\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}})=0.\hfill \end{array}\right.$$

(22)

Notice that Formulas (11) and (20) are exactly the same, but now $u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ and $w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ play the role of their counterparts $u_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $w_{{\mathbf{x}}_s,\omega }^{\varnothing }$ in Theorem 1. In contrast to the former case, problems (21) and (22) cannot be solved in closed form except for very simple geometries (i.e., when ${\mathrm{\Omega }}_{app}$ is a ball or a square), and they must be solved numerically. This will be the general situation for the initial guesses defined by (19). As already mentioned, we will use the indirect boundary integral formulation proposed in [28].

The iterative method we propose to approximate the true objects from knowledge of the physical parameters $c_e,a_e,c_i,a_i$ and $\gamma$ from measured data $u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{meas,r}$ for $r=1,\dots ,N_{receivers}$, $s=1,\dots ,N_{emitters}$, $\ell =1,\dots ,N_{freq}$ corresponding to the excitation of the medium via incident waves $u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{inc}$, is a multifrequency version of the one proposed by us in [37,39]. The algorithm proceeds in the following steps:

• Initialization. We initialize the method by computing the multifrequency topological derivative $D_T({\mathbb{R}}_{-}^d,\mathbf{x})$ defined in (17) from the individual ones (as described in Theorem 1) with weights $P_{{\omega }_{\ell }}^0$ defined by (18). Then, we select a value $0<{\alpha }_0<1$ and define the set

  $$\begin{array}{c}\hfill {\mathrm{\Omega }}_0=\{\mathbf{x}\in {\mathcal{R}}_{obs},D_T({\mathbb{R}}_{-}^d,\mathbf{x})\le (1-{\alpha }_0)\underset{\mathbf{y}\in {\mathcal{R}}_{obs}}{min}{D}_T({\mathbb{R}}_{-}^d,\mathbf{y})\}.\end{array}$$

  (23)
  For our numerical experiments, we will set in principle ${\alpha }_0=0.1$. For numerical purposes, the set ${\mathrm{\Omega }}_0$ needs now to be decomposed in its connected components, which is performed by using the Matlab function bwconncomp as explained in [42], and then each connected component is approximated by a star-shaped parameterization as detailed in [39].
• Iteration. For each step $j=1,\dots ,J_{max}$ we proceed as follows. Given the current approximation ${\mathrm{\Omega }}_{j-1}$,

  -

  We compute the monocromatic topological derivatives $D_T^{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}},\mathbf{x})$ for $\mathbf{x}\in {\mathcal{R}}_{obs}\backslash \overline{{\mathrm{\Omega }}_{j-1}}$ making use of Theorem 2 with ${\mathrm{\Omega }}_{app}={\mathrm{\Omega }}_{j-1}$, to define the new multifrequency topological derivative:

  $$D_T({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}},\mathbf{x})=\sum _{\ell =1}^{N_{frec}}{P}_{{\omega }_{\ell }}^j{D}_T^{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}},\mathbf{x}),\mathbf{x}\in {\mathcal{R}}_{obs}\backslash \overline{{\mathrm{\Omega }}_{j-1}},$$

  with weights $P_{{\omega }_{\ell }}^j=-{min}_{\mathbf{y}\in {\mathcal{R}}_{obs}\backslash \overline{{\mathrm{\Omega }}_{j-1}}}{D}_T^{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}},\mathbf{y})$. Notice that now the corresponding forward and adjoint problems involve the objects ${\mathrm{\Omega }}_{j-1}$ and therefore they have to be numerically solved.

  -

  We update the current set ${\mathrm{\Omega }}_{j-1}$ by adding to ${\mathrm{\Omega }}_{j-1}$ the points where the new topological derivative attains the largest negative values:

  $$\begin{array}{c}\hfill {\mathrm{\Omega }}_j={\mathrm{\Omega }}_{j-1}\cup \{\mathbf{x}\in {\mathcal{R}}_{obs}\backslash \overline{{\mathrm{\Omega }}_{j-1}},D_T({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}},\mathbf{x})\le \\ \hfill (1-{\alpha }_j)\underset{\mathbf{y}\in {\mathcal{R}}_{obs}\backslash \overline{{\mathrm{\Omega }}_{j-1}}}{min}{D}_T({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}},\mathbf{y})\}.\end{array}$$

  (24)

  In principle, ${\alpha }_j={\alpha }_{j-1}$. We calculate the number of connected components of ${\mathrm{\Omega }}_j$ and obtain star-shaped approximations of each of them, which form the final ${\mathrm{\Omega }}_j$.

  -

  We check if each monofrequency shape functional decreases. If $J_{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_j})<J_{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}})$ for $\ell =1,\dots ,N_{frec}$, then ${\mathrm{\Omega }}_j$ is accepted. Otherwise, we replace ${\alpha }_j$ by ${\alpha }_j/\sqrt{2}$ and compute again ${\mathrm{\Omega }}_j$ as in (24). The algorithm stops if either:

  *

  A maximum number of iterations $J_{max}$ is reached.

  *

  The measure of two consecutive approximations is negligible, namely if

  $$|\mathrm{meas}\left({\mathrm{\Omega }}_j\right)-\mathrm{meas}\left({\mathrm{\Omega }}_{j-1}\right)|<0.01\mathrm{meas}\left({\mathrm{\Omega }}_{j-1}\right).$$

  *

  Updating objects produces almost no variation for all the monofrequency functionals, i.e., if

  $$J_{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}})-J_{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_j})<0.01J_{{\omega }_{\ell }}({\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{j-1}}),\forall \ell =1,\dots ,N_{freq}.$$

  *

  The discrepancy principle for at least one of the frequencies is reached, that is, for any $\ell =1,\dots ,N_{freq}$

  $$\frac{\sqrt{{\sum }_{s=1}^{N_{emitters}}{\sum }_{r=1}^{N_{receivers}}{|u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{{\mathrm{\Omega }}_j}\left({\mathbf{z}}_r\right)-u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{meas,r}|}^2}}{\sqrt{{\sum }_{s=1}^{N_{emitters}}{\sum }_{r=1}^{N_{receivers}}{|u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{{\mathrm{\Omega }}_j}\left({\mathbf{z}}_r\right)-u_{{\mathbf{x}}_s,{\omega }_{\ell }}^{inc}\left({\mathbf{z}}_r\right)|}^2}}\le 1.2\delta ,$$

  where $\delta$ is the noise level (see the end of Section 2.2).

Notice that since $P_{{\omega }_{\ell }}^{j-1}\ne P_{{\omega }_{\ell }}^j$, at each step we are considering a different shape functional, which is a different linear combination of the same monochromatic functionals. In our previous works, only monochromatic data were considered, and therefore, the same shape functional was preserved during the iterations. As stated, this iterative procedure only allows us to create objects and add points to existing objects. Modifying the iterative scheme as indicated in [43], we can also remove points or destroy existing spurious components.

Next, we will illustrate the performance of the iterative multifrequency algorithm for the two-dimensional case. Figure 11 shows how the iterative method works for the demanding situation considered in the previous section where $a_e=1.5$ and $a_i=0.5$. Plot (a) represents the values of the topological derivative $D_T({\mathbb{R}}_{-}^2,\mathbf{x})$ for $\mathbf{x}\in {\mathcal{R}}_{obs}=[-2.5,2.5]\times [-2.5,-0.2]$ computed as in Theorem 1. Choosing ${\alpha }_0=0.1$, we find the initial approximation ${\mathrm{\Omega }}_0$ shown in panel (b) in white. In the same plot, we represent the values of the topological derivative when ${\mathrm{\Omega }}_{app}={\mathrm{\Omega }}_0$ computed by applying Theorem 2. The subsequent plots, indexed by the value of j in the title, represent the current approximation ${\mathrm{\Omega }}_{j-1}$ and the values of the topological derivative when ${\mathrm{\Omega }}_{app}={\mathrm{\Omega }}_{j-1}$ (which are used to define ${\mathrm{\Omega }}_j$ by Formula (24)). We observe that although our initial guess contains only one component, the second one is detected when $j=3$ and ${\mathrm{\Omega }}_3$ has the correct number of components. In fact, the topological derivative for $j=2$ already suggests it, but since the value of ${\alpha }_j$ is rather conservative, the approximate domain ${\mathrm{\Omega }}_2$ has only one component. The last approximation to the closest object is quite good (see panel (f)), including size, shape and orientation. For the deepest one, the location is correct, though it appears smaller. We only see the top illuminated part.

Figure 12 considers a more challenging situation: the deepest object is smaller and deeper. However, our method is also able to detect it. As we move this object deeper, we reach a depth at which we do not see it with our choice of frequencies, attenuation and thresholds, as can be seen in Figure 13. Nevertheless, the other object is still reasonably described. Attenuation effects hinder the detection of objects, when they are deep and small enough.

In all previous examples, the physical parameters $c_e=1$, $c_i=4$ and $\gamma =1$ were fixed. Since in realistic situations that appear, for instance, in geophysics or in biological problems [44], the simplification $\gamma =1$ could be unrealistic, we have repeated the experiments in Figure 11 and Figure 12 when replacing $\gamma =1$ by $\gamma =1/2$ and $\gamma =2$. Figure 14 shows the results for the configuration with two ellipses, and Figure 15 corresponds to the configuration where the deepest object is a ball.

Results evidence that the effect of attenuation is similar in the cases $\gamma =1$, $\gamma =1/2$ and $\gamma =2$, being the case $\gamma =2$ the one that seems to be less favorable. In particular, in Figure 15d, we observe that the small defect is not detected when $\gamma =2$. Although not illustrated here for brevity, when considering the configuration in Figure 13, in all of the cases ($\gamma =1$, $\gamma =1/2$ and $\gamma =2$) the deepest object is missed. We want to emphasize that we have not performed a thorough parametric study and it might happen that for other selections of speeds and frequencies, the case $\gamma =2$ could lead to more accurate reconstructions than $\gamma =1/2$ or $\gamma =1$.

Finally, let us remark that topological derivative-based methods are very robust with respect to noise. This has been extensively tested by many authors for different inverse problems (see for instance [17,20]). Synthetically generated data are usually corrupted by random noise in a similar way as we do in (8) and (9). When real experimental data are processed, noise is inevitably present. The reconstructions obtained in [23,45] evidence the robustness of the method when dealing with experimental data. For our inverse scattering problem in an attenuating half space, we have checked that increasing the level of noise $\delta$ defined in (8) from $\delta =0.01$ to $\delta =0.1$ results are almost identical. A detailed study of the effect of noise on a topological derivative-based approximation is performed in [46] in a Bayesian framework, though it is exemplified in different emitter/receiver configurations.

6. Other Boundary Conditions

The ideas developed in this paper do not depend on the conditions on the boundary of the defects. The results can be extended straightforwardly to Dirichlet (sound-soft defects) or Neumann (sound-hard defects) conditions. For sound-soft defects embedded in an attenuating media occupying the half plane/space ${\mathbb{R}}_{-}^d$, the counterpart forward problem corresponding to (6) for an incident wave $u_{{\mathbf{x}}_s,\omega }^{inc}$ emitted from the source point ${\mathbf{x}}_s$ at a frequency $\omega$ is:

$$\left\{\begin{array}{c}\Delta u_{{\mathbf{x}}_s,\omega }+k_e{\left(\omega \right)}^2{u}_{{\mathbf{x}}_s,\omega }={\delta }_{{\mathbf{x}}_s},\mathrm{in} {\mathbb{R}}_{-}^d\backslash \overline{\mathrm{\Omega }},\hfill \\ u_{{\mathbf{x}}_s,\omega }=0,\mathrm{on} \partial \mathrm{\Omega },\hfill \\ {\partial }_{x_d}{u}_{{\mathbf{x}}_s,\omega }=0\mathrm{on} \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r(u_{{\mathbf{x}}_s,\omega }-u_{{\mathbf{x}}_s,\omega }^{inc})-\mathit{ı}{k}_e\left(\omega \right)(u_{{\mathbf{x}}_s,\omega }-u_{{\mathbf{x}}_s,\omega }^{inc}))=0,\hfill \end{array}\right.$$

(25)

where the complex wave number $k_e\left(\omega \right)$ is defined in terms of the exterior wave velocity $c_e$ and the attenuation coefficient $a_e$ as in (3). In this case, adapting the results in [47,48,49] it can be proved that the formula for the topological derivative in Theorem 1 should be replaced by:

$$D_T^{\omega }({\mathbb{R}}_{-}^d,\mathbf{x})=\sum _{s=1}^{N_{emitters}}\mathrm{Re}\left[u_{{\mathbf{x}}_s,\omega }^{\varnothing }\left(\mathbf{x}\right)\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}\left(\mathbf{x}\right)\right]$$

(26)

where $u_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}$ solve (12) and (13). In case an initial guess ${\mathrm{\Omega }}_{app}$ of $\mathrm{\Omega }$ is available, we obtain the counterpart formula to (26) with $u_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}$ replaced by $u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ and $\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}$, and the solutions to the associated Dirichlet forward and adjoint problems with objects ${\mathrm{\Omega }}_{app}$:

$$\left\{\begin{array}{c}\Delta u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}+k_e{\left(\omega \right)}^2{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}={\delta }_{{\mathbf{x}}_s},\mathrm{in} {\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{app}},\hfill \\ u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}=0,\mathrm{on} \partial {\mathrm{\Omega }}_{app},\hfill \\ {\partial }_{x_d}{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}=0\mathrm{on} \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r(u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}-u_{{\mathbf{x}}_s,\omega }^{inc})-\mathit{ı}{k}_e\left(\omega \right)(u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}-u_{{\mathbf{x}}_s,\omega }^{inc}))=0,\hfill \end{array}\right.$$

(27)

and

$$\left\{\begin{array}{c}\Delta \overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}+k_e{\left(\omega \right)}^2\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=\sum _{r=1}^{N_{receivers}}\left(\overline{u_{{\mathbf{x}}_s,\omega }^{meas,r}-u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}\right){\delta }_{{\mathbf{z}}_r},\mathrm{in} {\mathbb{R}}_{-}^d\backslash \overline{{\mathrm{\Omega }}_{app}},\hfill \\ \overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=0,\mathrm{on} \partial {\mathrm{\Omega }}_{app},\hfill \\ {\partial }_{x_d}\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=0\mathrm{on} \mathrm{\Pi },\hfill \\ \underset{r\to \infty }{lim}{r}^{(d-1)/2}({\partial }_r\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}-\mathit{ı}{k}_e\left(\omega \right)\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}})=0.\hfill \end{array}\right.$$

(28)

These formulas can be derived by adapting the results in [49].

For sound-hard objects, modeled by the analogous problem to (25) obtained by replacing in this problem the boundary condition $u_{{\mathbf{x}}_s,\omega }=0$ on $\partial \mathrm{\Omega }$ by ${\partial }_{\mathbf{n}}{u}_{{\mathbf{x}}_s,\omega }=0$ on $\partial \mathrm{\Omega }$, the counterpart formula to that in Theorem 1 can be obtained either by adapting the results in [17,39] or by formally taking the limit as $k\to 0$ in (11):

$$D_T^{\omega }({\mathbb{R}}_{-}^d,\mathbf{x})=\sum _{s=1}^{N_{emitters}}\mathrm{Re}\left[\frac{d}{d-1}\nabla u_{{\mathbf{x}}_s,\omega }^{\varnothing }\left(\mathbf{x}\right)\nabla \overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}\left(\mathbf{x}\right)-k_e{\left(\omega \right)}^2{u}_{{\mathbf{x}}_s,\omega }^{\varnothing }\left(\mathbf{x}\right)\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}\left(\mathbf{x}\right)\right],$$

(29)

where again, $u_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}$ solve (12) and (13). The counterpart formula to that in Theorem 2 looks exactly as (29), with $u_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}$ replaced by $u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ and $\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}$, the solutions to the associated Neumann forward and adjoint problems with objects ${\mathrm{\Omega }}_{app}$. The forward problem is (27) with the condition $u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}=0$ on $\partial {\mathrm{\Omega }}_{app}$ replaced by ${\partial }_{\mathbf{n}}{u}_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}=0$ on $\partial {\mathrm{\Omega }}_{app}$, while the adjoint problem is given by (28) when replacing $\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=0$ on $\partial {\mathrm{\Omega }}_{app}$ by the Neumann condition ${\partial }_{\mathbf{n}}\overline{w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}}=0$ on $\partial {\mathrm{\Omega }}_{app}$.

Notice that for the three types of defects (penetrable, sound-soft and sound-hard), the same forward and adjoint problems are involved in the topological derivative Formulas (11), (26) and (29). However, depending on the type of defect, the formulas are different linear combinations of the terms $u_{{\mathbf{x}}_s,\omega }^{\varnothing }{w}_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $\nabla u_{{\mathbf{x}}_s,\omega }^{\varnothing }·\nabla \overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}$. In the case of not having information about the nature of the scatterers (and even in the case when objects of different nature are embedded in the same attenuating media), we could make use of a related mathematical tool, the topological energy [50], which for our model problem would be defined by

$$E_T^{\omega }\left(\mathbf{x}\right)=\sum _{s=1}^{N_{emitters}}|u_{{\mathbf{x}}_s,\omega }^{\varnothing }{|}^2{\left|\overline{w_{{\mathbf{x}}_s,\omega }^{\varnothing }}\right|}^2,$$

(30)

where $u_{{\mathbf{x}}_s,\omega }^{\varnothing }$ and $w_{{\mathbf{x}}_s,\omega }^{\varnothing }$ solve (12) and (13). For the use of a multifrequency version of this indicator function for a problem in ${\mathbb{R}}^2$ where objects of different nature are immersed in a non-attenuating media, we refer to [51]. The extension of the iterative method to cope with objects of different nature is not immediate, since the formula would be the counterpart to (30), but the corresponding fields $u_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ and $w_{{\mathbf{x}}_s,\omega }^{{\mathrm{\Omega }}_{app}}$ depend on the boundary conditions imposed on the components of ${\mathrm{\Omega }}_{app}$. Therefore, at its current state, we could only generalize the use of iterated multifrequency topological energies to the situation where all objects have the same nature. A monofrequency iterative method based on topological energy computations was explored in [52] for a related problem in ${\mathbb{R}}^3$ with no attenuation and penetrable objects.

7. Conclusions

We have proposed a topological derivative-based multifrequency iterative algorithm to identify objects buried in attenuated media with limited aperture data. In principle, these methods use no a priori information other than the knowledge of the measured data, the incident waves, the employed frequencies and the properties of the background medium. To simplify, here we have assumed that we know the material properties of the scatterers, and we look for their location, size, orientation and shape. One-step implementations of the algorithm provide initial approximations, which are improved in a few iterations. Components of smaller size than the main components, or buried deeper inside, can be detected. However, attenuation effects may prevent detecting objects depending on their size and depth for fixed ranges of frequencies.

A key point in the new multifrequency method defined in this paper is the way to combine monochromatic data. We have defined a strategy based on weighting individual topological derivatives, obtaining therefore a method where the global shape functional changes from one step to the next one, and consequently, the discrepancy principle must be tested individually. It would be possible to define the weights in a different way to keep the same functional or to define a different way to combine single-frequency data in a more efficient way, exploring the alternatives proposed in [21,53,54,55,56]. Indeed, the weights could be considered part of the inverse problem. Optimal selection of weights could extract more information from the measured data, perhaps increasing resolution of small deep objects.

The formulas obtained in this paper could be combined with other strategies to define iterative methods. For instance, the reconstructions obtained by the first computation of the topological derivative in ${\mathbb{R}}_{-}^d$ could be used as initial guesses for other types of iterative methods, such as Newton-type algorithms [57], level set methods [14,15] or other topological derivative-based approaches [58,59,60,61,62]. Iterative methods that alternate topological derivative computations with regularized Gauss-Newton iterations [63] could also be explored.

When the material properties of the objects are unknown, we could generate initial approximations to their geometry by related topological energy techniques [45,50,51,52] and to their material parameters by hybrid gradient methods [38]. Uncertainty due to noise can be addressed by combining topological priors with Bayesian techniques [46].