# FM-BEM and topological derivative applied to acoustic inverse scattering

Marc Bonnet, Nicolas Nemitz

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marc Bonnet, Nicolas Nemitz. FM-BEM and topological derivative applied to acoustic inverse scattering. IABEM 2006 Conference, 2006, Graz, Austria. 10.1007/978-3-540-47533-0_8. hal-00128912

HAL Id: hal-00128912

https://hal.science/hal-00128912v1

Submitted on 17 Oct 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# FM-BEM and Topological Derivative Applied to Acoustic Inverse Scattering

Marc Bonnet and Nicolas Nemitz

Solid Mechanics Laboratory (UMR CNRS 7649), Ecole Polytechnique,
91128 Palaiseau, France
{bonnet,nemitz}@lms.polytechnique.fr

**Summary.** This study is set in the framework of inverse scattering of scalar (e.g. acoustic) waves. A qualitative probing technique based on the distribution of topological sensitivity of the cost functional associated with the inverse problem with respect to the nucleation of an infinitesimally-small hard obstacle is formulated. The sensitivity distribution is expressed as a bilinear formula involving the free field and an adjoint field associated with the cost function. These fields are computed by means of a boundary element formulation accelerated by the Fast Multipole method. A computationally fast approach for performing a global preliminary search based on the available overspecified boundary data is thus defined. Its usefulness is demonstrated through results of numerical experiments on the qualitative identification of a hard obstacle in a bounded acoustic domain, for configurations featuring $O(10^5)$ nodal unknowns and $O(10^6)$ sampling points.

## 1 Introduction

Defect identification problems are often solved by minimization of a cost function featuring the experimental data and (if available) prior information. Such cost functions are non-convex and exhibit local minima. Despite that fact, traditional iterative minimization or equation-solving methods are usually preferred to global search techniques such as evolutionary algorithms due to the prohibitive computational cost of solving large numbers of forward wave scattering problems. To perform optimally, gradient-based iterative algorithms are used in conjunction with shape sensitivity techniques see e.g. [15, 18, 20].

Still, the stand-alone use of gradient-type minimization for such purposes is not satisfactory for its success is strongly dependent on a reliable prior information about the geometry of the hidden object. This has prompted the development of ‘sampling’ or ‘probe’ non-iterative methods, which may be used in isolation or as a preliminary step for choosing adequate initial guesses to be used in subsequent standard optimization schemes. Such methods are surveyed in a recent review article [21] and include the linear sampling
method [3, 17], not pursued here, and the application of topological sensitiv-
ity, which is the subject of this article. The concept of topological sensitivity
consists in quantifying the sensitivity of a cost function with respect to the cre-
ation of a scatterer $B_{\varepsilon}(\mathbf{x}_{\text{s}})$ of small characteristic radius $\varepsilon$ and given location
$\mathbf{x}_{\text{s}}$, as a function of $\mathbf{x}_{\text{s}}$. It appeared first in [6, 22] in connection with topo-
logical optimization of mechanical structures, allowing to define algorithms
where “excess” material is iteratively removed until a satisfactory shape and
topology is reached [9]. More recently, other investigations have studied the
topological sensitivity as a sampling tool for inverse scattering problems, in
the context of identification of cavities in 3D semi-infinite and infinite elastic
media [13] and in elastic 3D bounded bodies [2], and of elastic inclusions [14]
(see also [8] for 2D elastostatics and [7] for 2D linear acoustics).

The distribution of topological sensitivity can be expressed in terms of a
bilinear formula involving the free field and an adjoint field associated with the
cost function. However, the computational cost of solving the forward and ad-
joint problems and evaluating the topological sensitivity distribution on a fine
sampling grid increases rapidly with the non-dimensional wavenumber. The
purpose of this article is to propose the topological sensitivity field computed
by means of the Fast Multipole BEM (FM-BEM) [11, 12, 19] as the basis of a
computationally fast tool for probing acoustic media for hidden hard obstacles
on the basis of overdetermined boundary data, within the model framework
of forward scattering problems governed by the scalar Helmholtz equation.
To that end, the FM-BEM is in particular applied to evaluate in a fast way
the integral representation formulae expressing the free and adjoint fields at
a large number of sampling points inside the medium.

This article is organized as follows. After some preliminaries concerning
the forward and inverse problems of interest (Section 2), the concept of topo-
logical sensitivity is presented in Section 3. The FM-BEM approach for the
scalar Helmholtz equation is then summarized in Section 4. Finally, results of
numerical experiments on qualitative scatterer identification using computed
distributions of topological sensitivity are presented in Section 5, for configu-
rations featuring $O(10^5)$ nodal BE unknowns and $O(10^6)$ sampling points.

## 2 Forward and Inverse Problems

This article is concerned with the identification of rigid obstacles embedded in
acoustic media. The generic acoustic scattering problem of interest is defined
as follows. Let $\Omega$ denote a three-dimensional open domain, either bounded or
unbounded, with a sufficiently regular boundary $S$ and filled with an acoustic
medium characterized by wave velocity $c$ and mass density $\rho$; this configura-
tion will be referred to as the reference (i.e. obstacle-free) medium. Let $B^*$
denote a rigid scatterer (or a set thereof) bounded by the closed surface $\Gamma^*$, so
that $\Omega^* = \Omega \setminus \bar{B}^*$ is the acoustic region surrounding the scatterer. Steady-state
excitations on $S$ with angular frequency $\omega$ generate an acoustic pressure field
$u^*$ in the acoustic domain $\Omega^*$, governed by the following set of field equations and boundary conditions (collectively denoted by $\mathcal{P}(B^*)$ for later reference):

$$ (\Delta + k^2)u^* = 0 \quad (\text{in } \Omega^*), \\ \mathcal{P}(B^*): \quad p[u^*] = p^D \quad (\text{on } S), \\ p[u^*] = 0 \quad (\text{on } \Gamma^*), \tag{1} $$

where $k = \omega/c$ is the wavenumber, $w \to p[w] \equiv w_{,n} = \nabla w \cdot n$ denotes the normal derivative operator, $n$ is the normal on $S \cup \Gamma^*$ outward to $\Omega^*$, and $p^D$ is the prescribed Neumann data over $S$ (other types of boundary conditions may be considered as well). It is assumed that $\omega$ is not an eigenfrequency of any of the boundary-value problems arising in the ensuing developments.

In the inverse scattering problem of interest, an unknown obstacle $B^{\text{true}}$, of boundary $\Gamma^{\text{true}}$, is to be identified. The corresponding exact acoustic field $u^{\text{true}}$ is then governed by problem $\mathcal{P}(B^{\text{true}})$. With reference to problem $\mathcal{P}(B^{\text{true}})$, supplementary information is needed for the identification of $B^{\text{true}}$. Here, measured values $u^{\text{obs}}$ of acoustic pressure are assumed to be available over the measurement surface $S^{\text{obs}} \subset S$. Ideally (i.e. assuming that the physics is exactly described by the chosen linear acoustics setting and that no measurement errors are present), $u^{\text{obs}}$ is the trace of $u^{\text{true}}$ on $S^{\text{obs}}$. The identification of $B^{\text{true}}$ may then be formulated in terms of the minimization of a cost function. Generic cost function of format

$$ \mathcal{J}(\Omega^*) = \int_{S^{\text{obs}}} \varphi(u_R^*(\xi), u_I^*(\xi), \xi) \, d\Gamma \tag{2} $$

are considered, where $u^*$ is the boundary trace of the solution to the forward problem $\mathcal{P}(B^*)$ for an assumed obstacle configuration $B^*$, the subscripts 'R' and 'I' being used to indicate the real and imaginary parts of a complex quantity (i.e. $w_R = \text{Re}(w)$ and $w_I = \text{Im}(w)$). For instance, the output least-squares cost function associated to measurement $u^{\text{obs}}$ on $S^{\text{obs}}$, commonly used for such purposes, corresponds to

$$ \varphi(w_R, w_I, \xi) = \frac{1}{2} |w(\xi) - u^{\text{obs}}(\xi)|^2, \tag{3} $$

The minimization of such cost functions can be performed using many methods, all of which are iterative and need repeated evaluations of $J(\Omega^*)$. Traditional gradient-based minimization may converge within a moderate number of evaluations of $\mathcal{J}(\Omega^*)$ if the trial surface $\Gamma^*$ can be described in terms of a few geometrical parameters, but reach a local minimum which depends on the choice of initial guess. Global search techniques, e.g. evolutionary algorithms [16] or sampling methods based on the Metropolis algorithm [25], perform a global search (i.e. identify absolute and/or multiple minima), but at the cost of very large numbers of cost functions evaluations. In this article, the topological sensitivity is proposed as a tool for performing a qualitative global search at a computational cost which is far below that entailed by a true
global optimization technique. The results may then (for example) be used as
initial guesses in subsequent minimization-based inversion procedures.

3 Topological Sensitivity of the Cost Function

3.1 Notations

Let $B_{\varepsilon}(\mathbf{x}_s) = \mathbf{x}_s + \varepsilon\mathcal{B}$, where $\mathcal{B} \subset \mathbb{R}^3$ is a fixed bounded open set with boundary $\mathcal{S}$ and volume $|\mathcal{B}|$ containing the origin, define the region of space occupied by a hard obstacle of (small) size $\varepsilon > 0$ containing a fixed sampling point $\mathbf{x}_s$. It is convenient to introduce the scaled position vector $\tilde{\mathbf{\xi}}$ defined by

$$
\xi = x_s + \varepsilon \tilde{\xi} \quad (\xi \in B_{\varepsilon}, \tilde{\xi} \in \mathcal{B}) \qquad (4)
$$

In particular, this mapping recasts integrals over $B_{\varepsilon}$ and $\Gamma_{\varepsilon}$ into integrals over $\mathcal{B}$ and $\mathcal{S}$, respectively, and transforms the differential volume and area elements according to

$$
dV_{\xi} = \varepsilon^3 d\bar{V}_{\xi} \quad (\xi \in B_{\varepsilon}, \tilde{\xi} \in \mathcal{B}), \quad d\Gamma_{\xi} = \varepsilon^2 d\bar{\Gamma}_{\xi} \quad (\xi \in \Gamma_{\varepsilon}, \tilde{\xi} \in \mathcal{S}) \quad (5)
$$

Without loss of generality, $\mathbf{x}_s$ can be chosen as the center of $B_\varepsilon$, i.e. such that

$$
\int_{B_{\varepsilon}} (\boldsymbol{x} - \boldsymbol{x}_{\text{s}}) \, dV_{\varepsilon} = \mathbf{0}, \quad \text{i.e.} \quad \int_{\mathcal{B}} \tilde{\boldsymbol{\xi}} \, d\bar{V}_{\tilde{\varepsilon}} = \mathbf{0}. \qquad (6)
$$

Let $u^* = u^\varepsilon(\xi; \mathbf{x}_s)$ denote the solution to the scattering problem $\mathcal{P}(B_\varepsilon(\mathbf{x}_s))$ defined by (1), where $\Omega^* = \Omega_\varepsilon(\mathbf{x}_s) = \Omega \setminus \bar{B}_\varepsilon(\mathbf{x}_s)$ and $\bar{B}_\varepsilon(\mathbf{x}_s)$ is the closure of $B_\varepsilon(\mathbf{x}_s)$. Further, let $J(\varepsilon; \mathbf{x}_s)$ be defined by

$$
J(\varepsilon; \mathbf{x}_s) = \mathcal{J}(\Omega_\varepsilon(\mathbf{x}_s)) = \int_{S_{\text{obs}}} \varphi(u_\text{R}^\varepsilon(\xi), u_\text{I}^\varepsilon(\xi), \xi) \, d\Gamma, \quad (7)
$$

For convenience, explicit references to $\mathbf{x}_s$ will often be omitted in the sequel,
e.g. by writing $J(\varepsilon)$ or $u^\varepsilon(\xi)$ instead of $J(\varepsilon; \mathbf{x}_s)$ or $u^\varepsilon(\xi; \mathbf{x}_s)$.

The evaluation of $J(\varepsilon)$ entails solving for $u^\varepsilon$ the forward problem $\mathcal{P}(B_\varepsilon)$.
It is convenient, and customary, to decompose $u^\varepsilon$ according to

$$
u^{\epsilon} = u + v^{\epsilon}, \tag{8}
$$

where $u$, the *free field* defined as the response of the obstacle-free (reference)
medium $\Omega$ due to the given excitation $p^D$, solves

$$
\begin{align}
(\Delta + k^2)u &= 0 \quad (\text{in } \Omega), \tag{9} \\
p[u] &= p^D \quad (\text{on } S), \nonumber
\end{align}
$$

while $v^\epsilon$, the *scattered field*, solves

$$
\begin{align}
(\Delta + k^2)v^\epsilon &= 0 && (\text{in } \Omega_\epsilon), \nonumber \\
p[v^\epsilon] &= 0 && (\text{on } S), \tag{10} \\
p[v^\epsilon] &= -p[u] && (\text{on } \Gamma_\epsilon). \nonumber
\end{align}
$$
## 3.2 Expansion of $J(\varepsilon)$

To establish the topological sensitivity of $J(\varepsilon)$, one starts with the expansion

$$J(\varepsilon) = J(0) + \int_{S^{\text{obs}}} \text{Re}[\varphi_{,u} v^{\varepsilon}] d\Gamma + o(|v^{\varepsilon}|_{S^{\text{obs}}}), \quad (11)$$

with

$$\varphi_{,u} \equiv \left( \frac{\partial \varphi}{\partial u_{\text{R}}^{\varepsilon}} - i \frac{\partial \varphi}{\partial u_{\text{I}}^{\varepsilon}} \right) \Big|_{u^{\varepsilon}=u} \quad (12)$$

Let the adjoint field $\hat{u}$ be defined by

$$\begin{align}
(\Delta + k^2)\hat{u} &= 0 && (\text{in } \Omega), \nonumber \\
p[\hat{u}] &= \varphi_N && (\text{on } S^{\text{obs}}), \tag{13} \\
p[\hat{u}] &= 0 && (\text{on } S \setminus S^{\text{obs}}). \nonumber
\end{align}$$

Then, the reciprocity identity (i.e. third Green's formula) applied to the states $\hat{u}$ and $v^\varepsilon$ over the domain $\Omega_\varepsilon$ leads, by virtue of the boundary conditions in (10) and (13), to the identity

$$\int_{S^{\text{obs}}} \varphi_{,u} v^{\varepsilon} \, d\Gamma + \int_{\Gamma_{\varepsilon}} p[\hat{u}] v^{\varepsilon} \, d\Gamma + \int_{\Gamma_{\varepsilon}} \hat{u} p[u] \, d\Gamma = 0 \quad (14)$$

As a result, the integral in the r.h.s. of (11) is converted into integrals over the vanishing cavity surface. Besides, since both $u$ and $\hat{u}$ are also defined inside $B_\varepsilon$, the last integral in (14) can be recast into a domain integral over $B_\varepsilon$ by means of the divergence formula. Expansion (11) then takes the form

$$\begin{align}
J(\varepsilon) = J(0) + \operatorname{Re} \left\{ \int_{B_{\varepsilon}} \left[ \nabla u \cdot \nabla \hat{u} - k^2 u \hat{u} \right] dV - \int_{\Gamma_{\varepsilon}} v^{\varepsilon} p[\hat{u}] d\Gamma \right\} \nonumber \\
+ o(|v^{\varepsilon}|_{S^{\text{obs}}}), \tag{15}
\end{align}$$

The first integral in (15) features a density whose definition does not depend on $\varepsilon$, and its expansion about $\varepsilon = 0$ can therefore be obtained by simply using the scaled coordinates (4), (5) and expanding $[\nabla u \cdot \nabla \hat{u} - k^2 u \hat{u}] (x_s + \varepsilon \bar{\xi})$. In contrast, the second integral of (15) features the scattered field $v^\varepsilon$, which depends on $\varepsilon$. Its asymptotic behaviour must then be obtained from that of $v^\varepsilon$ on $\Gamma_\varepsilon$ (taking into account the fact that $\Gamma_\varepsilon$ also depends on $\varepsilon$). This step is based on exploiting an integral equation reformulation of equations (10).

## 3.3 Governing Integral Equation Formulation for the Scattered Field

The governing problem (10) for the scattered field $v^\varepsilon = v^\varepsilon(\cdot; x_s)$ can be recast as the boundary integral equation [1, 4]:

$$\frac{1}{2} v^{\varepsilon}(\boldsymbol{x}) + \int_{\Gamma_{\varepsilon}} \mathcal{H}(\boldsymbol{x}, \boldsymbol{\xi}; k) v^{\varepsilon}(\boldsymbol{\xi}) \, d\Gamma_{\varepsilon} = - \int_{\Gamma_{\varepsilon}} \mathcal{G}(\boldsymbol{x}, \boldsymbol{\xi}; k) p[u](\boldsymbol{\xi}) \, d\Gamma_{\varepsilon} \quad (\boldsymbol{x} \in \Gamma_{\varepsilon}), \quad (16)$$
where the Green's function $G(x, \xi, k)$ is defined by

$$ (\Delta_{\xi} + k^2)G(x, \xi, k) + \delta(\xi - x) = 0 \quad (\xi \in \Omega), \quad \mathcal{H}(x, \xi, k) = 0 \quad (\xi \in S), \quad (17) $$

and $\mathcal{H}(x, \xi; k) = \nabla G(x, \xi, k) \cdot n(\xi)$ is the normal derivative of $G(x, \xi, k)$ (the nabla symbol $\nabla$, when used in front of such kernel, conventionally indicates a gradient with respect to the second argument $\xi$). Moreover, the free and adjoint fields have the explicit expressions

$$ u(x) = \int_S G(x, \xi; k) p^D(\xi) \, d\Gamma_{\xi}, \quad \hat{u}(x) = \int_S G(x, \xi; k) \varphi_{,u}(\xi) \, d\Gamma_{\xi}. \quad (18) $$

It is convenient for the present purposes to split $(\mathcal{G}, \mathcal{H})$ according to:

$$ \begin{aligned} \mathcal{G}(x, \xi; k) &= G(x, \xi; k) + G_C(x, \xi; k) \\ \mathcal{H}(x, \xi; k) &= H(x, \xi; k) + H_C(x, \xi; k), \end{aligned} \quad (19) $$

where $(G, H)$ is the well-known singular free-space fundamental solution for the Helmholtz equation, given by

$$ G(x, \xi; k) = \frac{1}{4\pi r} e^{ikr}, \quad H(x, \xi; k) = [r \cdot n(\xi)] \frac{\mathrm{ir}r - 1}{4\pi r^3} e^{ikr}, \quad (20) $$

with $r = \xi - x$ and $r = |\xi - x| = |r|$, and the complementary part $(G_C, H_C)$ is not singular at $\xi = x$. On using decomposition (19), equation (16) becomes

$$ \begin{aligned} & \frac{1}{2} f(x) + \int_{\Gamma_{\varepsilon}} H(x, \xi; k) v^{\varepsilon}(\xi) \, d\Gamma_{\varepsilon} + \int_{\Gamma_{\varepsilon}} H_C(x, \xi; k) v^{\varepsilon}(\xi) \, d\Gamma_{\varepsilon} \\ & = - \int_{\Gamma_{\varepsilon}} G(x, \xi; k) p[u](\xi) \, d\Gamma_{\varepsilon} - \int_{\Gamma_{\varepsilon}} G_C(x, \xi; k) p[u](\xi) \, d\Gamma_{\varepsilon} \quad (x \in \Gamma_{\varepsilon}). \end{aligned} \quad (21) $$

## 3.4 Leading Asymptotic Contribution to the Scattered Field

To study the asymptotic behaviour of integral equation (21) as $\varepsilon \to 0$, it is useful to introduce further scaled geometric quantities:

$$ x = \varepsilon \bar{x}, \quad r = \varepsilon \bar{r}, \quad r = \varepsilon \bar{r} \quad (x, \xi \in \Gamma_{\varepsilon}; \bar{x}, \bar{\xi} \in S) \quad (22) $$

in addition to definition (4) of $\bar{\xi}$. The leading contributions as $\varepsilon \to 0$ to the fundamental kernels featured in equation (21) are

$$ \begin{aligned} G(x, \xi; k) &= \varepsilon^{-1} G(\bar{x}, \bar{\xi}) + O(1) \\ H(x, \xi; k) &= \varepsilon^{-2} H(\bar{x}, \bar{\xi}) + O(1) \end{aligned} \quad (x, \xi \in \Gamma_{\varepsilon}), \quad (23) $$

for the singular kernels $(G, H)$ defined by (20), where

$$ G(\bar{x}, \bar{\xi}) = \frac{1}{4\pi \bar{r}}, \quad H(\bar{x}, \bar{\xi}) = -\frac{\bar{r} \cdot n(\bar{\xi})}{4\pi \bar{r}^3} \quad (24) $$
are the fundamental kernels for the Laplace equation, and

$$
\begin{align}
G_C(\mathbf{x}, \boldsymbol{\xi}; k) &= G_C(\mathbf{x}_s, \mathbf{x}_s; k) + O(\varepsilon) \tag{25} \\
H_C(\mathbf{x}, \boldsymbol{\xi}; k) &= H_C(\mathbf{x}_s, \mathbf{x}_s; k) + O(\varepsilon) \quad (\mathbf{x}, \boldsymbol{\xi} \in \Gamma_{\varepsilon}). \nonumber
\end{align}
$$

for the nonsingular kernels ($G_C, H_C$).

On performing the coordinate transformations (4), (22) and using esti-
mates (23), (25) together with (5b), one finds that

$$
\begin{equation}
\begin{aligned}
\int_{\Gamma_{\varepsilon}} H_C(\mathbf{x}, \boldsymbol{\xi}; k) v^{\varepsilon}(\boldsymbol{\xi}) \, d\Gamma_{\varepsilon} &= O(\varepsilon^2 |v^{\varepsilon}|) \\
\int_{\Gamma_{\varepsilon}} G_C(\mathbf{x}, \boldsymbol{\xi}; k) p[u](\boldsymbol{\xi}) \, d\Gamma_{\varepsilon} &= O(\varepsilon^3).
\end{aligned}
\tag{26}
\end{equation}
$$

The limiting form of integral equation (21) as $\varepsilon \to 0$, retaining only the leading
contributions, is hence found to be

$$
\bar{\mathcal{L}}_S [v^\varepsilon(\mathbf{x}_S + \varepsilon \bar{\boldsymbol{\xi}})](\bar{\mathbf{x}}) = -\varepsilon \nabla u(\mathbf{x}_S) \cdot \text{int}_S \frac{1}{4\pi\bar{r}} \mathbf{n}(\bar{\boldsymbol{\xi}}) \, d\bar{\Gamma}_{\varepsilon} + o(\varepsilon) \quad (27)
$$

where $\bar{\mathcal{L}}_S$, defined by

$$
[\bar{\mathcal{L}}_S f](\bar{x}) = \frac{1}{2} f(\bar{x}) + \int_S H(\bar{x}, \bar{\xi}) f(\bar{\xi}) \, d\bar{\Gamma}_{\bar{\xi}} \quad (\bar{x} \in S), \quad (28)
$$

is in fact the governing integral operator associated with exterior Neumann
problems for the Laplace equation in the normalized domain $\mathbb{R}^3 \setminus \bar{\mathcal{B}}$. Equa-
tion (27) indicates that the scattered field is of order $O(\varepsilon)$ on $\Gamma_{\varepsilon}$:

$$
v^{\varepsilon}(\boldsymbol{\xi}) = v^{\varepsilon}(\mathbf{x}_s + \varepsilon \bar{\boldsymbol{\xi}}) = \varepsilon \nabla u(\mathbf{x}_s) \cdot \mathcal{V}(\bar{\boldsymbol{\xi}}) + o(\varepsilon) \quad (\boldsymbol{\xi} \in \Gamma_{\varepsilon}, \bar{\boldsymbol{\xi}} \in \mathcal{S}) \quad (29)
$$

where the vector function $\mathcal{V}(\bar{\xi})$ solves the exterior Laplace problem

$$
[\bar{\mathcal{L}}_S \mathcal{V}](\bar{x}) = -\nabla u(\bar{x}_s) \cdot \int_S \frac{1}{4\pi\bar{r}} \mathbf{n}(\bar{\xi}) \, d\bar{\Gamma}_{\bar{\xi}} \quad (30)
$$

i.e.

$$
\begin{align}
\Delta \mathcal{V} &= 0 && (\bar{\xi} \in \mathbb{R}^3 \setminus \bar{\mathcal{B}}), \nonumber \\
\nabla_{\bar{\xi}} \mathcal{V} \cdot \mathbf{n} &= -\mathbf{n} && (\bar{\xi} \in \mathcal{S}), \tag{31} \\
|\bar{\xi}|^2 \mathcal{V} &= O(1) && (|\bar{\xi}| \to +\infty) \nonumber
\end{align}
$$

It is important to note that $\mathcal{V}(\bar{\xi})$ does not depend on the sampling point $\mathbf{x}_s$,
and hence needs to be computed only once.

3.5 Topological Derivative

On substituting (29) into (15) and taking (5) into account, one finally arrives
at the following expansion of $J(\varepsilon)$:
$$J(\varepsilon) = J(0) + \varepsilon^3 \mathcal{T}(\mathbf{x}_s) + o(\varepsilon^3) \quad (32)$$

where $\mathcal{T}(\mathbf{x}_s)$, the topological derivative of $J(\varepsilon)$, is given in terms of the free and adjoint fields by

$$\mathcal{T}(\mathbf{x}_s) = \mathrm{Re}\{\nabla \hat{u} \cdot \mathcal{A}(S) \cdot \nabla u - |\mathcal{B}| k^2 \hat{u} u\}(\mathbf{x}_s) \quad (33)$$

and with the second-order tensor $\mathcal{A}(S)$ defined by

$$\begin{aligned} \mathcal{A}(S) &= |\mathcal{B}| I - \int_S [\mathcal{V}(\bar{\xi}) \otimes \mathbf{n}(\bar{\xi})] \, d\bar{\Gamma}_{\bar{\xi}} && (34) \\ &= |\mathcal{B}| I - \int_S [\mathbf{n}(\bar{\xi}) \otimes \mathcal{V}(\bar{\xi})] \, d\bar{\Gamma}_{\bar{\xi}} \end{aligned}$$

where the second equality (i.e. the fact that $\mathcal{A}(S)$ is symmetric) easily stems from the third Green's identity and the definition (31) of $\mathcal{V}$.

For arbitrary surfaces $S$ which are sufficiently regular for integral equation (30) to be mathematically meaningful (this includes surfaces with edges and corners, e.g. box-shaped scatterers, but precludes infinitely-thin screens), the vector density $\mathcal{V}$ may be found by e.g. solving numerically three sets of BEM equations for exterior Laplace problems, a computationally modest task.

For the simplest case of a rigid spherical obstacle (where $\mathcal{B}$ is the unit sphere, $|\mathcal{B}| = 4\pi/3$, and on which $\mathbf{n}(\bar{\xi}) = -\bar{\xi}$) one easily finds by analytical means that $\mathcal{V}(\bar{\xi}) = \bar{\xi}/2$. Then, (34) is readily found to be given by

$$\mathcal{A}(S) = 2\pi I. \quad (35)$$

## 3.6 Practical Computation of Topological Sensitivity

The developments thus far are based on the Green's function $\mathcal{G}$ defined by (17), and lead to almost explicit formulae for $\mathcal{T}(\mathbf{x}_s)$, their only non-explicit component being the auxiliary density $\mathcal{V}$, which must be computed numerically except for simple shapes of the trial scatterer $\mathcal{B}$.

In practice, this explicit character is retained only for geometries $\Omega$ such that the corresponding Green's function is known analytically. Such cases correspond with geometrically simple configurations, e.g. the acoustic half-space. For configurations where the Green's function is not available, the free and adjoint fields and the nonsingular kernels may be sought as solutions of boundary integral equations [1, 4]. The free and adjoint fields, defined by (9) and (13), satisfy the well-known integral identities

$$c(\mathbf{x})u(\mathbf{x}) + \int_S H(\mathbf{x}, \boldsymbol{\xi}; k)u(\boldsymbol{\xi}) \, d\Gamma_{\boldsymbol{\xi}} = \int_S G(\mathbf{x}, \boldsymbol{\xi}; k)p^D(\boldsymbol{\xi}) \, d\Gamma_{\boldsymbol{\xi}} \quad (36)$$

$$c(\mathbf{x})\hat{u}(\mathbf{x}) + \int_S H(\mathbf{x}, \boldsymbol{\xi}; k)\hat{u}(\boldsymbol{\xi}) \, d\Gamma_{\boldsymbol{\xi}} = - \int_{S^{\text{obs}}} G(\mathbf{x}, \boldsymbol{\xi}; k)\varphi_{\text{N}}(\boldsymbol{\xi}) \, d\Gamma_{\boldsymbol{\xi}} \quad (37)$$
which provide integral equations for $x \in S$ (with $c(x) = 1/2$, except at points $x$ where $S$ is only piecewise smooth, such as edges or corners, for which $c(x)$ is also known) and integral representation formulae for $x \in \Omega$ (with $c(x) = 1$). Differentiation of the latter under the integral sign provide integral representations of $\nabla u(x)$ and $\nabla u(x_s)$ in $\Omega$.

Hence, evaluation of the topological sensitivity field as given by (33) entails the solution of integral equations (36) and (37) for $u$ and $\hat{u}$ on $S$, followed by an evaluation of $u(x_s)$, $\nabla u(x_s)$ and $\hat{u}(x_s)$, $\nabla u(x_s)$ by means of (36) and (37) used as integral representations. The first step may involve large numbers of unknowns if the diameter of $\Omega$ spans more than a few wavelengths. Moreover, the present objective being to comprehensively explore a 3-D region of space for hidden scatterers by examining the distribution of $T(x_s)$, formula (33) is to be evaluated at a large number of sampling points $x_s$. Both steps can then be considerably accelerated using the FM-BEM.

# 4 The Fast Multipole Method for Helmholtz Equation

## 4.1 BEM Discretization

Equations (36) and (37) are in this article solved by means of the simplest BEM discretization, which employs flat triangular boundary elements with straight edges and piecewise-linear $C^0$ interpolation of $u$ and $\hat{u}$ (other choices, e.g. $C^0$ quadratic interpolation, would of course have been possible). All numerical results presented hereinafter have been obtained on that basis. The primary unknowns are the values of $u$ or $\hat{u}$ at the mesh nodes, i.e. at all the element vertices. Equations (36) and (37) are collocated at all mesh nodes. All singular element integrals associated with the kernel $H(x, \xi; k)$ are zero because of the piecewise-flat geometry representation, which of course simplifies the implementation. Denoting by $N$ the total number of nodes (and hence of unknowns), this procedure gives rise to the linear systems of equations

$$[A]\{u\} = \{b\} \qquad (38)$$

$$[A]\{\hat{u}\} = \{\hat{b}\} \qquad (39)$$

where the $N$-vectors $\{u\}$ and $\{\hat{u}\}$ collect all nodal values of $u$ and $\hat{u}$. The discussion to follow will focus on solving system (38), the adjoint system (39) being of course solved in exactly the same way with $\{b\}$ replaced by $\{\hat{b}\}$.

As the problem size $N$ grows, direct solvers become impractical or infeasible with respect to both computing time and storage, mainly due to the fully-populated nature of the BEM matrix $[A]$, and iterative solvers are used instead. Since $[A]$ is a non-symmetric, invertible matrix (except when $k$ is an eigenvalue for $\Omega$ and homogeneous Neumann BCs), the iterative solution technique most often used for systems such as (38) is the generalized minimal residual (GMRES) algorithm (see e.g. [10]), which is applicable to general
invertible square matrices. Such algorithms are based on matrix-vector evaluations, and therefore do not require actual assembly and storage of [A]. The GMRES algorithm requires repeated evaluations of the residual

$$ \{b\} - [A]\{u\} \qquad (40) $$

where {$u$} is a given solution candidate, which is updated at each GMRES iteration. Hence, one needs to compute (discretized versions of) the double-layer and single-layer potentials featured in the left-hand and right-hand sides, respectively, of (36) and (37) for known densities.

Traditional BE methods lead to a $O(N^2)$ computational cost for each residual evaluation, because element integrals computed for a collocation point cannot be reused for another collocation point. By adopting the Fast Multipole boundary element method (FM-BEM), each residual can be computed within a $O(N\log N)$ time. The implementation used here, concisely described in the remainder of this section, follows Darve [5] and Sylvand [23, 24].

## 4.2 Expansion of the Fundamental Solution

The starting point for the FM-BEM is the following representation of the full-space fundamental solution (20):

$$ G(\mathbf{x}, \boldsymbol{\xi}; k) = \lim_{p \to \infty} \frac{ik}{4\pi} \int_{\hat{S}} e^{-ik(\hat{\mathbf{s}} \cdot \hat{\mathbf{x}})} T_p(\hat{\mathbf{s}}, \mathbf{r}_0) e^{ik(\hat{\mathbf{s}} \cdot \tilde{\boldsymbol{\xi}})} \, d\Gamma_{\hat{s}} \qquad (41) $$

where $\hat{S} = \{\hat{\mathbf{s}}, |\hat{\mathbf{s}}| = 1\}$ is the unit sphere, the position vector $\mathbf{r} = \boldsymbol{\xi} - \mathbf{x}$ has been decomposed as

$$ \mathbf{r} = (\boldsymbol{\xi}_0 - \mathbf{x}_0) + (\boldsymbol{\xi} - \boldsymbol{\xi}_0) - (\mathbf{x} - \mathbf{x}_0) = \mathbf{r}_0 + \tilde{\boldsymbol{\xi}} - \tilde{\mathbf{x}} \qquad (42) $$

in terms of two poles $\mathbf{x}_0$ and $\boldsymbol{\xi}_0$, and with the transfer function $T_p(\hat{\mathbf{s}}, \mathbf{r}_0)$ defined by

$$ T_p(\hat{\mathbf{s}}, \mathbf{r}_0) = \frac{(-i)^n}{4\pi} \sum_{n=0}^{p} (-1)^n (2n + 1) h_n^{(1)}(kr_0) P_n(\hat{\mathbf{s}} \cdot \hat{\mathbf{r}}_0) \qquad (43) $$

In (43), the $h_n^{(1)}$ and $P_n$ are respectively the spherical Hankel functions of first kind and the Legendre polynomials, and $r_0 \equiv |\boldsymbol{\xi}_0 - \mathbf{x}_0|$. Moreover, for any vector $\mathbf{z} \in \mathbb{R}^3$, a hat symbol indicates the corresponding unit vector, i.e. $\hat{\mathbf{z}} = \mathbf{z}/|\mathbf{z}|$. Representation (41) holds under the condition

$$ |r_0| > |\tilde{\boldsymbol{\xi}} - \tilde{\mathbf{x}}| \qquad (44) $$

The poles $\mathbf{x}_0$ and $\boldsymbol{\xi}_0$ are actually meant to be chosen close to $\boldsymbol{\xi}$ and $\mathbf{x}$, respectively, so as to satisfy the stronger condition

$$ |\boldsymbol{\xi} - \mathbf{x}_0| > |\mathbf{x} - \mathbf{x}_0| \quad \text{and} \quad |\mathbf{x} - \boldsymbol{\xi}_0| > |\boldsymbol{\xi} - \boldsymbol{\xi}_0| \qquad (45) $$
The representation of $H(\boldsymbol{x}, \boldsymbol{\xi}; k)$ obtained by differentiating (41) is then

$$H(\boldsymbol{x}, \boldsymbol{\xi}; k) = \lim_{p \to \infty} \frac{-k^2}{4\pi} \int_{\hat{S}} e^{-ik(\hat{s} \cdot \hat{x})} T_p(\hat{s}, \boldsymbol{r}_0) [\hat{s} \cdot \boldsymbol{n}(\boldsymbol{\xi})] e^{ik(\hat{s} \cdot \hat{\xi})} \, d\Gamma_{\hat{s}} \quad (46)$$

In practice, representations (41) and (46) are approximated by (i) using the transfer function $T_p(\hat{s}, \boldsymbol{r}_0)$ with a finite index $p$ and (ii) replacing the integral over the unit sphere with a quadrature rule with $Q$ points $\hat{s}_q$ and weights $w_q$. So, one replaces (41) and (46) with the approximations

$$G(\boldsymbol{x}, \boldsymbol{\xi}; k) \approx \frac{ik}{4\pi} \sum_{q=1}^{Q} w_q e^{-ik(\hat{s}_q \cdot \hat{x})} T_p(\hat{s}_q, \boldsymbol{r}_0) e^{ik(\hat{s}_q \cdot \hat{\xi})} \quad (47)$$

$$H(\boldsymbol{x}, \boldsymbol{\xi}; k) \approx \frac{ik}{4\pi} \sum_{q=1}^{Q} w_q e^{-ik(\hat{s}_q \cdot \hat{x})} T_p(\hat{s}_q, \boldsymbol{r}_0) [\hat{s}_q \cdot \boldsymbol{n}(\boldsymbol{\xi})] e^{ik(\hat{s}_q \cdot \hat{\xi})} \quad (48)$$

The choice of points $\hat{s}_q$ and weights $w_q$, and their number $Q$, depends on the truncation order $p$ used in the transfer function $T_p$. On parameterizing unit vectors $\hat{s} \in \hat{S}$ using spherical angular coordinates $(\theta, \phi)$, a commonly used choice [23, 5] consists in using the $Q(p) = (p+1)(2p+1)$ quadrature points of the form $\hat{s}_q = (\theta_a, \phi_b)$, where $\cos\theta_a$ ($1 \le a \le p+1$) are the abscissae for the $(p+1)$-point Gauss-Legendre 1-D quadrature rule over $[-1, 1]$, $\phi_b = 2\pi b/(2p+1)$ ($1 \le b \le 2p+1$) are uniformly-spaced abscissae on $[0, 2\pi]$. The associated weights are $w_q = 2\pi w_a^\theta/(2p+1)$, where $w_a^\theta$ are the Gauss-Legendre weights for the $(p+1)$-point rule.

Now, let $S_x$ and $S_\xi$ denote two disjoint portions of $S$, and let the poles $\boldsymbol{x}_0$ and $\boldsymbol{\xi}_0$ be chosen close to $S_x$ and $S_\xi$, respectively, in such a way that (45) holds for any $\boldsymbol{x} \in S_x$ and $\boldsymbol{\xi} \in S_\xi$. Consider the computation of

$$I(\boldsymbol{x}) = \int_{S_{\xi}} G(\boldsymbol{x}, \boldsymbol{\xi}; k) v(\boldsymbol{\xi}) \, d\Gamma_{\xi} \quad (\boldsymbol{x} \in S_x)$$

for a given density $v(\boldsymbol{\xi})$, which is a typical contribution to the evaluation of the residual (40). On substituting (47), one obtains

$$I(\boldsymbol{x}) \approx \frac{ik}{4\pi} \sum_{q=1}^{Q(p)} w_q e^{-ik(\hat{s}_q \cdot \hat{x})} T_p(\hat{s}_q, \boldsymbol{r}_0) \int_{S_{\xi}} e^{ik(\hat{s}_q \cdot \hat{\xi})} v(\boldsymbol{\xi}) \, d\Gamma_{\xi} \quad (\boldsymbol{x} \in S_x) \quad (49)$$

So, the same integral over $S_\xi$ can be re-used for all collocation points $\boldsymbol{x} \in S_x$. Moreover, for a chosen portion $S_\xi$, this is true for any portion $S_x$ and associated pole $\boldsymbol{x}_0$ such that condition (45) holds. Computations of the form (49) can be decomposed into three stages: (i) compute for each quadrature point of $\hat{S}$ the multipole moment $\mathcal{R}(\hat{s}_q; \boldsymbol{\xi}_0)$:

$$\mathcal{R}(\hat{s}_q; \boldsymbol{\xi}_0) = \int_{S_{\xi}} e^{ik(\hat{s}_q \cdot \hat{\xi})} v(\boldsymbol{\xi}) \, d\Gamma_{\xi} \quad (50)$$
(ii) multiply $\mathcal{R}(\hat{s}_q; \xi_0)$ by the transfer function $T_p(\hat{s}_q, r_0)$, to obtain local expansion coefficients $\mathcal{L}(\hat{s}_q; x_0)$ at $x_0$:

$$ \mathcal{L}(\hat{s}_q; x_0) = T_p(\hat{s}_q, \xi_0 - x_0)\mathcal{R}(\hat{s}_q; \xi_0) \quad (51) $$

(iii) for all $x \in S_x$, transfer $\mathcal{L}(\hat{s}_q; x_0)$ locally from $x_0$ to $x$ and perform the numerical quadrature over $\hat{S}$, to obtain (an approximation of) $I(x)$:

$$ I(x) \approx \frac{ik}{4\pi} \sum_{q=1}^{Q(p)} w_q e^{-ik(\hat{s}_q \cdot \hat{x})} \mathcal{L}(\hat{s}_q; x_0) \quad (52) $$

The one-level fast multipole method consists in partitioning the spatial region containing $S$ into cubic cells of identical sizes whose vertices lie on a regular cubic grid. Each pair $(S_x, S_\xi)$ is such that $S_x = S \cap C_x$ and $S_\xi = S \cap C_\xi$, where $(C_x, C_\xi)$ is any pair of disjoint cubic cells. The poles $x_0$ and $\xi_0$ are the respective cell centroids. The one-level FM-BEM has a complexity of $O(N^{3/2})$ per iteration for equations of type (38), which is of course better than the $O(N^2)$ complexity of traditional BEM, but not optimal. Further acceleration is provided by using the multi-level fast multipole method, where the size of clusters $S_x, S_\xi$ depends of their distance.

## 4.3 Multilevel FM-BEM Algorithm

To exploit optimally the acceleration afforded by (49), a hierarchical oct-tree structure of elements is introduced. For that purpose, a cube containing the boundary $S$, called ‘level-0 cell’, is divided into eight cubes (level-1 cells), each of which is divided in the same fashion, and so on. A level-$\ell$ cell is divided into level-$(\ell+1)$ cells unless it contains less than a preset (relatively small) number $E$ of boundary elements (such cells are termed leaves). A noteworthy feature of the FM-BEM applied to Helmholtz-type equations is that to achieve the same accuracy in approximations (47), (48) at all levels, the truncation parameter $p$ is level-dependent. A often-used formula [5] for the adjustment of $p$ is of the form

$$ p(\ell) = \sqrt{3kd}(\ell) + C\mathrm{Log}_{10}(\sqrt{3kd}(\ell) + \pi) \quad (53) $$

where $d(\ell)$ is the size of a cubic $\ell$-level cell (so, $d(\ell+1) = 2d(\ell)$ and $p(\ell+1)$ is roughly $2p(\ell)$) and $C$ is a constant. As a result, the set of quadrature points $\hat{s}_q$ on $\hat{S}$ is also level-dependent: each level necessitates a distinct set $\hat{s}_q^\ell$ of quadrature points and associated weights. For the present implementation, values of $C$ such that $2 \le C \le 8$ were found to provide a acceptable compromise between accuracy and cost, and $C = 4$ was actually used.

The FM-BEM algorithm then consists of:

* An upward pass where multipole moments (50) are first computed for the lowest-level cells and then recursively aggregated by moving upward in the
tree until level 2 (for which there are 4×4×4 cells) is reached. Letting $C(\xi_0)$ denote the set of children of a given $\ell$-level cell $C(\xi_0)$, i.e. of $(\ell+1)$-level cells $C(\xi'_0)$ contained in $C(\xi_0)$, this operation relies on the identity

$$ \mathcal{R}(\hat{s}_q^{\ell+1}; \xi_0) = \sum_{C(y-0) \in \mathcal{C}(x_0)} e^{ik(\hat{s}_q^{\ell+1} \cdot (\xi'_0 - \xi_0))} \mathcal{R}(\hat{s}_q^{\ell+1}; \xi'_0) \quad (54) $$

for shifting the origin from the center $\xi'_0$ of a level-$(\ell+1)$ cell to the center $\xi_0$ of a level-$\ell$ cell in the contribution of a given cell.

Then, it is necessary to interpolate, i.e. compute $\mathcal{R}(\hat{s}_q^{\ell}; \xi_0)$ at the quadrature points $\hat{s}_q^{\ell}$ from the previously determined values $\mathcal{R}(\hat{s}_q^{\ell+1}; \xi_0)$. The procedure used follows [23, 5] and is not detailed here.

* A downward pass where the coefficients of local expansions are first computed at level $\ell = 2$ and then evaluated at lower-level cells by tracing the tree structure downwards. This operation relies on the identity

$$ \mathcal{L}(\hat{s}_q^{\ell}; x'_0) = e^{-ik(\hat{s}_q^{\ell+1} \cdot (x_0 - x'_0))} \mathcal{L}(\hat{s}_q^{\ell}; x_0) \quad (55) $$

for shifting the origin from the center $x_0$ of a level-$(\ell)$ cell to the center $x'_0$ of a level-$(\ell+1)$ cell. This operation is not performed at the root level, i.e. when $\ell+1 = 2$.

Then, the contributions of all level-$(\ell+1)$ cells belonging to the interaction list $\mathcal{I}(x_0)$ of the level-$\ell+1$ cell $C(x_0)$ (i.e. all such cells which are non-adjacent to $C(x_0)$ but whose father is adjacent to the father of $C(x_0)$) are aggregated:

$$ \sum_{C(y-0) \in \mathcal{I}(x_0)} T_p(\hat{s}_q^{\ell}, y-0-x_0) \mathcal{R}(\hat{s}_q^{\ell}; y-0) $$

and the result is added to $\mathcal{L}(\hat{s}_q^{\ell}; x'_0)$ given by (55)

Then, the values $\mathcal{L}(\hat{s}_q^{\ell}; x'_0)$ are converted to values at the quadrature points $\hat{s}_q^{\ell+1}$ by a 'reverse interpolation', or 'anterpolation', procedure.

* When the lowest level is reached, all quadratures of the form (52) are finally performed, where $x_0$ is the centers of a leaf cell, thus evaluating all far-field contributions to the residual at all collocation points.

Moreover, for all leaf cells $C(x_0)$ and all collocation points $x \in C(x_0)$, the near-field contributions are computed by evaluating using traditional integration methods the element contributions for all elements located in $C(x_0)$ and all cells of same level adjacent to $C(x_0)$.

The computation of integral representation formulae for $u(x_s), \nabla u(x_s)$ and $\hat{u}(x_s), \nabla u(x_s)$ at all chosen sampling points follows the same approach, with collocation points $x \in S$ replaced with sampling points $x_s \in \Omega$. For sampling points lying in leaf cells not adjacent to any same-level cell intersecting $S$, the integral representations result from far-field interactions (i.e. fast-multipole contributions) only. Besides, all multipole moments used in this step are those corresponding to the solution of (38) or (39), i.e. those evaluated at the last iteration of the GMRES solution algorithm
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** CPU time for one evaluation of residual (40) as a function of the number $N$ of nodal unknowns, with matrix-vector products computed using classical BEM techniques (lozenges, theoretical complexity $O(N^2)$) or the FM-BEM (circles, theoretical complexity $O(N\log N)$).</figcaption></figure>

## 4.4 Numerical Verification of Theoretical Complexity

To check that the theoretical complexity of the FM-BEM is achieved and illustrate the computational advantage brought by the FM-BEM over the conventional BEM, the simple situation of a spherical acoustic domain subjected on its surface to a uniform normal velocity, is considered. BEM solutions for this problem have been computed for a sequence of meshes with decreasing element size. For each mesh, the prescribed frequency is selected so that the number of nodes per wavelength is (approximately) the same for all meshes. Figure 1 shows the CPU time used for one evaluation of residual (40), using either classical BEM integration techniques (with theoretical complexity $O(N^2)$) or the FM-BEM (with theoretical complexity $N\log N$), as a function of the number $N$ of nodal unknowns. Functions of the form $C_1N^2$ and $C_2N\log N$ closest to the actual recorded values of CPU($N$) are also plotted on the same graph. The theoretical complexity for both the classical BEM and the FM-BEM are very well verified by the actual CPU times. The FM-BEM, as expected, performs much better for large problems.

# 5 Preliminary Identification via Topological Sensitivity: Numerical Examples

To illustrate the approach described in Sections 3 and 4, the following configuration has been considered: the bounded acoustic domain is the cube defined by $\Omega(L) = \{|\xi_i| \le L (i = 1, 2, 3)\}$, with $L = 8a$ or $L = 16a$ in terms of a
<figure><img src="image_4.png" /><figcaption>Fig. 2. Pattern of excitation surfaces $S_q$ on each face of external boundary $S$.</figcaption></figure>

reference length $a$. A simulated testing configuration is based on 30 experiments, each of which consists in applying a uniform excitation over a small region $S_q$ ($1 \le q \le 30$) of $S$ and recording the acoustic pressure over the whole boundary $S$ (i.e. at all BE mesh nodes), so that $S^{\text{obs}} = S$. The acoustic excitation is such that the wavelength is $\lambda = 3a$. Each of the six faces of the cubical domain $\Omega$ supports five excitation surfaces $S_q$, each a disk of radius $a$, arranged as depicted on Fig. 2. The cost function for the inverse problem is defined by

$$ \mathcal{J}(\Omega^*) = \frac{1}{2} \sum_{q=1}^{30} \int_S |u^* - u_q^{\text{obs}}|^2 d\Gamma_{\xi} \quad (56) $$

where $u_q^{\text{obs}}$ is the data obtained for the $q$-th applied excitation, with $u_q^{\text{obs}} = u^{\text{true}}$ in the absence of data noise. The centroid $\mathbf{x}^{\text{true}}$ of true scatterer $B^{\text{true}}$ to be identified is located at $\mathbf{x}^{\text{true}} = (2a, 3a, 2a)$.

To facilitate the graphical interpretation, a thresholded variant $\hat{T}(\mathbf{x}_s)$ of $T(\mathbf{x}_s)$ is introduced according to

$$ \hat{T}(\mathbf{x}_s) = \begin{cases} T(\mathbf{x}_s), & T \le C T_{\min} \\ 0, & T > C T_{\min} \end{cases} \quad \text{with } C = 0.25. \quad (57) $$

The BE meshes used for computing the free field $u$, the adjoint field $\hat{u}$ and the simulated data $u_q^{\text{true}}$ are made of three-noded flat triangular elements, arranged in a regular mesh with approximately 15 nodes per wavelength. For the purposes of computing the simulated error-free data $u^{\text{true}}$ for each synthetic experiment, BE meshes of $\Gamma^{\text{true}}$ have been set up as well. Table 1 indicates the numbers of nodes and elements supported by the BE meshes.

<table><caption>Table 1. Number of element and DOFs supported by the BE meshes.</caption><thead><tr><th rowspan="2">Cube size</th><th colspan="2">Cube</th><th colspan="2">Obstacle</th><th colspan="2">Total</th></tr><tr><th>Elements</th><th>nodes</th><th>Elements</th><th>nodes</th><th>Elements</th><th>nodes</th></tr></thead><tbody><tr><td>2L = 16a</td><td>76800</td><td>38402</td><td>336</td><td>170</td><td>77136</td><td>38572</td></tr><tr><td>2L = 32a</td><td>307200</td><td>153602</td><td>336</td><td>170</td><td>307536</td><td>153772</td></tr></tbody></table>
First, the identification of one spherical scatterer is considered, for three cases with increasing scatterer radii 0.2a, 0.4a and 0.8a. The field $\mathcal{T}$ has been computed for each case on the basis of error-free synthetic data, over a sampling grid $\mathcal{S}$ made of 100 × 100 × 100 sampling points located on the vertices of a regular cubic grid, centered at the origin and with grid spacing $\Delta x_s = 16a/101$, uniformly filling the whole acoustic domain bounded by $\mathcal{S}$. Figures 3, 4 and 5 show, for each scatterer configuration considered in turn, the distribution of the thresholded topological sensitivity $\hat{T}$ defined by (57) in the three coordinate planes containing the true scatterer centroid $\mathbf{x}^{\text{true}}$. In all cases, $\hat{T}$ (and hence $\mathcal{T}$) is seen to attain its lowest values in zones corresponding to, or close to, the actual true scatterer location. It should however be mentioned that $\mathcal{T}$ has been observed to also achieve low values in zones located near edges or corners of $\mathcal{S}$, where there is no scatterer. To emphasize and illustrate this remark, Fig. 6 shows the iso-surfaces of the field $\mathcal{T}$ corresponding to $\mathcal{T} = 0.55\mathcal{T}_{\min}$, computed on the entire search grid $\mathcal{S}$ (Fig. 6a) and on truncated search grids such that $\{|\xi_i| \le 6.5a \ (i = 1, 2, 3)\}$ (Fig. 6b) and $\{|\xi_i| \le 5a \ (i = 1, 2, 3)\}$ (Fig. 6c), i.e. in which sampling points situated at distances less than 1.5a and 3a, respectively, from $\mathcal{S}$ have been taken out of $\mathcal{S}$. Figure 6c shows that, at least in the central region $\{|\xi_i| \le 5a \ (i = 1, 2, 3)\}$ of the acoustic domain, low negative values of $\mathcal{T}$ occur only in a small region which is consistent with the actual location of the scatterer.

Then, the effect of data errors is considered for the spherical scatterer of radius $0.8a$, by using synthetic data $u^{\text{obs}}$ in the form $u^{\text{obs}} = u^{\text{true}}(1 + \eta)$, where $\eta$ are random numbers with zero mean and uniform distribution over the interval $[-0.1, 0.1]$. Figure 7 shows the distribution of $\hat{T}$ in the three coordinate planes containing $\mathbf{x}^{\text{true}}$. A comparison between Figs. 5 and 7 reveals that the distribution of $\hat{T}$ is only marginally affected by the data noise.

The examples shown thus far illustrate the capability of $\mathcal{T}$, here defined on the basis of an asymptotic analysis involving vanishing spherical obstacles, to identify the location of spherical obstacles of finite extent. Now, the identification of a non-spherical, box-shaped obstacle whose sides are aligned along the coordinate axes and of finite size $0.8a \times 0.8a \times 1.6a$ and whose centroid is still $\mathbf{x}^{\text{true}} = (2a, 3a, 2a)$, is considered. Figure 8 shows the distribution of $\hat{T}$ in the three coordinate planes containing $\mathbf{x}^{\text{true}}$. The true obstacle is again satisfactorily located.

A last example illustrates the case of a larger acoustic domain $\Omega(16a)$, instead of $\Omega(8a)$ considered up to this point, with the same wavelength $\lambda = 3a$ as before. The ‘true’ scatterer (again a sphere of radius $0.8a$) is still located at $\mathbf{x}^{\text{true}} = (2a, 3a, 2a)$, and hence is located at a larger distance (expressed in wavelengths) from the measurement surfaces. The sampling grid $\mathcal{S}$ is now made of 150 × 150 × 150 regularly-spaced sampling points, with a grid spacing now of $\Delta x_s = 32a/151$. Figure 9 shows the distribution of $\hat{T}$ in the three coordinate planes containing $\mathbf{x}^{\text{true}}$. Presumably as a result of greater remoteness (and hence lower identifiability) of the scatterer, these distributions show, in addition to the correct one, secondary spatial zones where the presence of a
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Identification of spherical hard scatterer of radius 0.2a: distribution of thresholded topological sensitivity $\tilde{T}(\mathbf{x}_s)$ for sampling points $\mathbf{x}_s$ in coordinate planes $\xi_1 = x_1^{\text{true}}$ (**a**), $x_2 = x_2^{\text{true}}$ (**b**) and $x_3 = x_3^{\text{true}}$ (**c**).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Identification of spherical hard scatterer of radius 0.4a: distribution of thresholded topological sensitivity $\tilde{T}(\mathbf{x}_s)$ for sampling points $\mathbf{x}_s$ in coordinate planes $\xi_1 = x_1^{\text{true}}$ (**a**), $x_2 = x_2^{\text{true}}$ (**b**) and $x_3 = x_3^{\text{true}}$ (**c**).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Identification of spherical hard scatterer of radius 0.8a: distribution of thresholded topological sensitivity $\mathcal{T}(\mathbf{x}_s)$ for sampling points $\mathbf{x}_s$ in coordinate planes $\xi_1 = x_1^{\text{true}}$ (**a**), $x_2 = x_2^{\text{true}}$ (**b**) and $x_3 = x_3^{\text{true}}$ (**c**).</figcaption></figure>
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Identification of spherical hard scatterer of radius 0.4a: iso-surfaces of $T(\mathbf{x}_s)$ for $T = 0.55T_{\min}$, computed on the entire search grid $S$ (a) and on truncated search grids such that $\{|\xi_i| \le 6.5a \ (i = 1, 2, 3)\}$ (b) and $\{|\xi_i| \le 5a \ (i = 1, 2, 3)\}$ (c). Values of $T(\mathbf{x}_s)$ lower than the iso-value are inside the iso-surface.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 7.** Identification of spherical hard scatterer of radius $0.4a$, synthetic data with 10% noise: distribution of thresholded topological sensitivity $\hat{T}(\mathbf{x}_s)$ for sampling points $\mathbf{x}_s$ in coordinate planes $\xi_1 = x_1^{\text{true}}$ (**a**), $x_2 = x_2^{\text{true}}$ (**b**) and $x_3 = x_3^{\text{true}}$ (**c**).</figcaption></figure>
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 8.** Identification of box-shaped scatterer: distribution of thresholded topological sensitivity $\bar{T}(\mathbf{x}_s)$ for sampling points $\mathbf{x}_s$ in coordinate planes $\xi_1 = x_1^{\text{true}}$ (**a**), $x_2 = x_2^{\text{true}}$ (**b**) and $x_3 = x_3^{\text{true}}$ (**c**).</figcaption></figure>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 9.** Identification of spherical hard scatterer of radius 0.4a in large domain: distribution of thresholded topological sensitivity $\bar{T}(\mathbf{x}_s)$ for sampling points $\mathbf{x}_s$ in coordinate planes $\xi_1 = x_1^{\text{true}}$ (a), $x_2 = x_2^{\text{true}}$ (b) and $x_3 = x_3^{\text{true}}$ (c).</figcaption></figure>
small scatterer is also consistent with the data. Still, the lowest values of $\hat{T}$ furnish a reasonable indication of the true obstacle location, as seen on the 3-D plots of iso-surfaces $\mathcal{T} = 0.55\mathcal{T}_{\min}$ and $\mathcal{T} = 0.7\mathcal{T}_{\min}$ of Fig. 10 for the truncated grid defined by $\{|\xi_i| \le 14a (i = 1, 2, 3)\}$. Again, one notes that values of $\mathcal{T}$ close to its minimum $\mathcal{T}_{\min}$ occur only in the vicinity of the correct obstacle location (excluding, as before, regions close to the external surface $S$).

<figure><img src="image_20.png" /><figcaption>**Fig. 10.** Identification of spherical hard scatterer of radius $0.8a$ embedded in large domain $\Omega(16a)$: iso-surfaces of (a) $\mathcal{T}(\mathbf{x}_s)$ for $\mathcal{T} = 0.55\mathcal{T}_{\min}$ and (b) $\mathcal{T} = 0.7\mathcal{T}_{\min}$ computed on the truncated search grid such that $\{|\xi_i| \le 14a (i = 1, 2, 3)\}$. Values of $\mathcal{T}(\mathbf{x}_s)$ lower than the iso-value are inside the iso-surface.</figcaption></figure>
Typical CPU times and GMRES iteration counts are provided in Table 2. It is interesting to observe that the overall CPU times for $\Omega(16a)$, which involves roughly 4 times as many nodal unknowns as $\Omega(8a)$, are about 5 times higher than those for $\Omega(8a)$ (while the expected ratio for traditional BEM would be $4^3 = 64$), and that the GMRES iteration counts are only fractionally higher for $\Omega(16a)$. All computations have been performed on a Linux PC computer with one 3 GHz processor.

<table><caption>**Table 2.** CPU times and (in parentheses) GMRES iteration count for computing the true, free and adjoint solutions on the boundary, and CPU times for computing the topological sensitivity over the whole sampling grid $S$.</caption><thead><tr><th></th><th>u<sup>true</sup> on S &cup; I<sup>true</sup></th><th>u on S</th><th>u&#770; on S</th><th>I on S</th></tr></thead><tbody><tr><td>2L = 16a</td><td>1444s (435)</td><td>969s (282)</td><td>1163s (342)</td><td>852s</td></tr><tr><td>2L = 32a</td><td>6461s (439)</td><td>5615s (388)</td><td>6818s (476)</td><td>1860s</td></tr></tbody></table>

## 6 Conclusion

In this article, a computationally fast qualitative technique for probing acoustic media for hidden hard obstacles on the basis of overdetermined boundary data, based on the computation via the FM-BEM of the distribution of topological sensitivity of the cost functional associated with the inverse problem, has been presented. Its usefulness has been demonstrated through results of numerical experiments on the qualitative identification of a hard obstacle in a bounded acoustic domain, for configurations featuring $O(10^5)$ nodal unknowns and $O(10^6)$ sampling points, resulting in overall computing times of a few hours on a 3 GHz PC computer. There is ample scope for increasing these computational sizes. Besides, the proposed approach can be developed for many other physical models, e.g. elastodynamics or electromagnetic waves.

## References

1. M. Bonnet: Boundary Integral Equations Methods for Solids and Fluids. John Wiley and Sons, 1999.

2. M. Bonnet, B. B. Guzina: Sounding of finite solid bodies by way of topological derivative. Internat. J. Numer. Methods Engrg. 61 (2004) 2344–2373.

3. D. Colton, A. Kirsch: A simple method for solving inverse scattering problems in the resonance region. Inverse Problems 12 (1996) 383–393.

4. D. Colton, R. Kress: Integral Equation Method in Scattering Theory. John Wiley and Sons, 1983.

5. E. Darve: The fast multipole method: numerical implementation. J. Comp. Phys. 160 (2000) 195–240.
6. H. A. Eschenauer, V. V. Kobelev, A. Schumacher: Bubble method for topology and shape optimization of structures. Structural Optimization 8 (1994) 42–51.

7. G. R. Feijóo: A new method in inverse scattering based on the topological derivative. Inverse Problems 20 (2004) 1819–1840.

8. R. Gallego, G. Rus: Identification of cracks and cavities using the topological sensitivity boundsary integral equation. Comp. Mech. 33 (2004) 154–163.

9. S. Garreau, P. Guillaume, M. Masmoudi: The topological asymptotic for PDE systems: the elasticity case. SIAM J. Control Opt. 39 (2001) 1756–1778.

10. A. Greenbaum: Iterative methods for solving linear systems. SIAM, Philadelphia, 1997.

11. L. Greengard, V. Rokhlin: A fast algorithm for particle simulations. J. Comp. Phys. 73 (1987) 325–348.

12. N. A. Gumerov, R. Duraiswami: Fast multipole methods for the Helmholtz equation in three dimensions. Elsevier Series in Electromagnetism, 2005.

13. B. B. Guzina, M. Bonnet: Topological derivative for the inverse scattering of elastic waves. Quart. J. Mech. Appl. Math. 57 (2004) 161–179.

14. B. B. Guzina, I. Chikichev: From imaging to material identification: a generalized concept of topological sensitivity. In: Proceedings of ECCOMAS 2004 (P. Neittanmäki, et al. eds.).

15. A. Kirsch: The domain derivative and two applications in inverse scattering theory. Inverse Problems 9 (1993) 81–96.

16. Z. Michalewicz, D. B. Fogel: How to solve it: modern heuristics. Springer, 2004.

17. S. Nintcheu Fata, B. B. Guzina: A linear sampling method for near-field inverse problems in elastodynamics. Inverse Problems 20 (2004) 713–736.

18. S. Nintcheu Fata, B. B. Guzina, M. Bonnet: A computational basis for elasto-dynamic cavity identification in a semi-infinite solid. Comp. Mech. 32 (2003) 370–380.

19. N. Nishimura: Fast multipole accelerated boundary integral equation methods. Appl. Mech. Rev. 55 (2002) 299–324.

20. R. Potthast: On the convergence of a new Newton-type method in inverse scattering. Inverse Problems 17 (2001) 1419–1434.

21. R. Potthast: A survey on sampling and probe methods for inverse problems. Inverse Problems 22 (2006) R1–R47.

22. A. Schumacher: Topologieoptimierung von Bauteilstrukturen unter Verwendung von Lochpositionierungskriterien. Ph.D. thesis, Universität Siegen, 1995.

23. G. Sylvand: La méthode multipôle rapide en électromagnétisme: performances, parallélisation, applications. Ph.D. thesis, Ecole Nationale des Ponts et Chaussées, Noisy le Grand, 2002.

24. G. Sylvand: Complex industrial computations in electromagnetism using the fast multipole method. In: Mathematical and numerical methods aspects of wave propagation (G. Cohen, E. Heikkola, P. Joly, P. Neittaanmäki eds.), pp. 657–662, Springer, 2003.

25. A. Tarantola: Inverse problem theory and methods for model parameter estimation. SIAM, 2005.