Inverse acoustic scattering by solid obstacles: topological sensitivity and its preliminary application

A to pological sensitivity approach is developed for the imaging of penetrable solid obstacles embedded in a fluid background medium by way of inverse acoustic scattering. To this end, an asymptotic form for the scattered field caused by a nucleating spherical solid obstacle in the reference fluid medium is derived within the boundary integral equation framework, where the required limiting behaviour of the acceleration field inside the perturbation is established based on solutions of two simplified fluid–solid interaction problems. With this result, the equivalence of acoustic scattering by fluid and solid scatterers in terms of asymptotic behaviour is observed and numerically validated. The direct and adjoint-field topological sensitivity expressions for transient and time-harmonic excitations are obtained accordingly. The utility of the proposed method as a tool for preliminary obstacle reconstruction and bulk modulus characterization is illustrated through numerical examples and an exploratory experimental study. On the basis of adjusted topological sensitivity formulas for fluid reference domains containing pre-existing solid inhomogeneities, an iterative 3D obstacle reconstruction approach is established and its performance with synthetic data is demonstrated.


Introduction
The inverse scattering problem of using either far-or near-field patterns of the scattered wave fields to identify obstacles embedded in an acoustic medium is a challenging subject, [1,2] with applications to a wide range of areas such as underwater sonar acoustics, seismic imaging and medical diagnosis. In the context of full waveform analysis, this class of problems is often affiliated with non-linear minimization and thus considerable computational cost due to repeated evaluation of the forward solution [e.g. 3,4]. As shown in a number of recent studies, however, sampling-based non-iterative techniques [5,6] that reconstruct obstacles from the distribution of an established indicator function provide a computationally effective preliminary imaging tool. The topological sensitivity method, [7,8] as one of its variants and the focus of this study, formulates the indicator function from a physical standpoint.
The concept of topological sensitivity (TS), since its inception in the context of structural shape optimization, [9,10] has been generalized and applied to deal with inverse scattering problems in acoustics, [7,[11][12][13][14][15][16][17] electromagnetism, [18][19][20][21][22][23] and elastodynamics. [24][25][26][27][28][29][30][31][32][33] In the reconstruction approach, the TS indicator function is defined as the sensitivity of a given cost functional (that would commonly be used as a platform for non-linear minimization) with respect to the nucleation of an infinitesimal obstacle at a prescribed sampling point in the reference (background) medium. Accordingly, the support of hidden obstacles is exposed through the spatial distribution of topological sensitivity, namely the regions where TS attains pronounced negative values. Typically, the topological sensitivity can be formulated explicitly in terms of either the germane Green's function or suitable wave fields computed for the reference (scatterer-free) domain, which makes this class of inverse solutions appealing from the computational point of view. With the help of a selected threshold value and minimization of TS with respect to obstacle material parameters, this method may also provide preliminary information about the size and material properties of the hidden obstacle. [28][29][30] Concerning the inverse scattering of acoustic waves, [11] introduced the idea of using the topological sensitivity to construct impenetrable obstacles of finite extent in a homogeneous background medium by way of time-harmonic waveforms. This concept was furthered in [7,15] to admit penetrable scatterers and heterogeneous background media. In [12] and [16], the foregoing analyses were extended to the time domain, where the feasibility of exposing (resp. impenetrable and penetrable) scatterers using acoustic wavelets of sufficient smoothness was demonstrated. In these topological sensitivity developments, the vanishing scatterer is assumed to be acoustic, i.e. of perfectly incompressible fluid with zero shear stiffness (Poisson's ratio = 1/2). However, in a number of relevant engineering applications such as underwater sonar imaging, the hidden scatterer is often a solid body, which is incompatible with the above assumption. For the involved inverse fluid-solid interaction scattering problem, some previous attempts using different methods have been proposed, such as the optimization method, linear sampling method or factorization method [e.g. [34][35][36][37][38]. In the context of topological sensitivity method, the case of a scattering solid obstacle in an acoustic (fluid) background medium has not been addressed yet.
To help mitigate the problem, the concept of topological sensitivity is herein particularized to quantify the perturbation of a given cost functional due to the creation of an infinitesimal solid obstacle in the reference acoustic medium. Through preliminary limiting analysis of the boundary integral representation for the scattered field, its asymptotic behaviour is affiliated with the acceleration field inside the vanishing obstacle, which is then derived from solutions of two simplified fluid-solid interaction problems involving spherical obstacles. Here, it is observed and numerically validated that the limiting behaviour of the scattered acoustic field caused by a penetrable solid obstacle is equivalent to that of an acoustic scatterer with the same mass density and bulk modulus, which in a certain sense supports the solid obstacle to acoustic obstacle simplification considered by engineers at a lower computational cost. The topological sensitivity formulas for transient and timeharmonic excitations are expressed explicitly in terms of the incident acoustic field and either (i) the Green's function or (ii) the adjoint field for the reference (scatterer-free) domain.
Through numerical examples and a preliminary experimental investigation, it is shown that the featured developments carry the potential of furnishing both geometric and material information of the hidden scatterers. To further improve the geometric reconstruction of hidden obstacles, an iterative 3D obstacle reconstruction approach is established based on topological sensitivity formulas adjusted for acoustic reference domains with known solid obstacles. Its performance with synthetic data is demonstrated.

Direct problem
With reference to Figure 1, consider the fluid layer = {ξ | − < ξ 3 < 0} ∈ R 3 in a reference Cartesian frame {O; ξ 1 , ξ 2 , ξ 3 }, housing a solid obstacle B with smooth boundary S = ∂ B. The background fluid, bounded by the top Dirichlet boundary T = {ξ | ξ 3 = 0} and bottom Neumann boundary B = {ξ | ξ 3 = − }, is characterized by constant mass density ρ and bulk modulus κ. The prescribed top and bottom boundary conditions are to represent free fluid surface and impenetrable bottom. The solid obstacle is assumed to be elastic, homogeneous and isotropic with mass density ρ * , bulk modulus κ * and shear modulus μ * . With the above definitions, − = \ B is the exterior region surrounding the obstacle and ∂ − = ∪ S with = B ∪ T denotes its boundary. Hereon, it is assumed that (i) the obstacle is acoustically illuminated via transient source distribution , and (iii) all featured acoustic and elastodynamic fields have quiescent past.

Free field
In the obstacle-free (reference) medium , the prescribed excitation E(ξ , t) gives rise to a pressure distribution p F (ξ , t), termed the free field, that satisfies the field equation and boundary conditions: Here p F ,n = n ·∇ p F , where n is the normal on oriented outward from , c = √ κ/ρ denotes the acoustic wave speed, andṗ F = ∂ p F /∂t indicates differentiation with respect to time. In the case of time-harmonic excitation E(ξ , t) = g(ξ )e −iωt with angular frequency ω, (1a) can be reduced to where the complex-valued P F (ξ , ω) denotes the time-harmonic free-field pressure and the time-dependent term e −iωt is omitted. Here, it is assumed that ω is not an eigenfrequency of any of the boundary-value problems appearing in the analysis.

Scattered field
Following the usual treatment of scattering problems, perturbation of the pressure field due to the presence of an obstacle B ⊂ can be conveniently quantified in terms of the scattered fieldp where p is the (total) pressure field in − generated by E(ξ , t). Assuming that the free field p F can be obtained separately beforehand, the forward scattering problem (i.e. the task of calculating p in − given E and B) can be conveniently written in terms ofp as follows: where the boundary conditions in (3e) represent continuity of surface traction and normal acceleration on the fluid-solid interface, see e.g. [39] for the counterpart of (3) in terms of time-harmonic scattering. Here, u is the displacement field inside the solid obstacle, n is the normal on ∂ − oriented towards the exterior of − , is the fourth-order elasticity tensor of the solid obstacle, where I j denotes the symmetric jth-order identity tensor, and t = −(C * : ∇u) · n is the traction vector on S = ∂ B. The radiation condition accompanying (1)-(3) will be introduced next in the form of boundary integral representation.

Boundary integral representation
To facilitate the ensuing developments, it is useful to recall the Riemann (time) convolution between two spatiotemporal fields, namely and to introduce the Green's function for the reference domainĜ(ξ , x, t), which satisfies (1) with where δ is the Dirac delta function whose dimension reflects that of its argument. Let G(ξ , x, t) and H (ξ , x, t) ≡ n·∇ ξ G(ξ , x, t) denote, respectively, the time-domain acoustic fundamental solution and its normal derivative for the free space R 3 due to a point source δ(ξ − x)δ(t). Using the method of images, [40,41] the Green's functionĜ(ξ , x, t) and its normal derivativeĤ (ξ , x, t) ≡ n·∇ ξĜ (ξ , x, t) for the reference fluid layer can be written as: where e 3 signifies the unit vector in the ξ 3 −direction, and the terms in the series account for the multiple reflections off the top and bottom boundaries of the layer, T and B . Due to the fact that the reference domain extends to infinity in the ξ 1 − and ξ 2 −directions, (1)-(3) are complemented by the requirement that there is no energy influx from infinity. In the context of infinite media, this requirement is commonly imposed in terms of the Sommerfeld radiation condition. [42] For problems involving waveguides such as that examined in this study, a relevant statement of the radiation condition can be written, following the work of [43] for elastic waves, as where q ∈ {p F ,p}, q ⊆ is the support of q, R = {ξ | ξ 2 1 + ξ 2 2 = R 2 , − < ξ 3 < 0} is a cylindrical surface of radius R and n is the normal on R oriented away from the origin.
On the basis of (1)-(6), the scattered pressure fieldp in the exterior of the obstacle B can be expressed in terms of the boundary distributions ofp andp ,n over ∂ B as:

Inverse problem and topological sensitivity
With reference to Figure 1, consider the inverse problem where a solid obstacle, B true , characterized by mass density ρ true , bulk modulus κ true and shear modulus μ true , is to be reconstructed and identified from the knowledge of: (i) the illuminating transient source distribution E(ξ , t) = g(ξ ) f (t), (ii) the background medium, including the acoustic parameters (ρ, κ), geometry and boundary conditions, and (iii) observations of the induced pressure field, p true , taken over the measurement surface obs ⊂ − . In the following analysis, the latter measurements are denoted by p obs , so that p obs (ξ , t) = p true (ξ , t), ξ ∈ obs , t ∈ [0, T ] under ideal measurement conditions.
To tackle the featured inverse problem, a cost functional quantifying the misfit between the sensory data p obs and their simulations p (computed for a trial obstacle B) is introduced, where − = \ B, m = (ρ * , κ * , μ * ) collects the material parameters of B, and ϕ ≥ 0 is a distance function, assumed to be differentiable with respect to its first argument. A commonly used example of distance function assumes the weighted least-squares format, namely where W ξ (ξ ) ≥ 0 and W t (t) ≥ 0 denote, respectively, the spatial and temporal weighting functions.

Topological sensitivity
To cater for the reconstruction of penetrable solid obstacles in acoustic media, the concept of topological sensitivity is hereon particularized to quantify the perturbation of the cost functional where is the cost functional computed for the free field p F and h(a)> 0, characterizing the leading term, is a function of the vanishing obstacle size. For completeness, it is further assumed that With reference to (2), letp a = p a − p F denote the scattered field caused by the infinitesimal obstacle (i.e. solution of the forward transmission problem (3) wherein B is replaced by B a ≡ B a (x o )). As a → 0, the scattered fieldp a is on physical grounds expected to vanish, i.e.
whereas the free field p F is independent of a. By virtue of this result, (2) written in terms of p a and the definition of the cost functional (8), J ( − a , m; E) can be expanded with respect top a as which, together with (10) and (11), results in the expression for topological sensitivity, In light of (15), the evaluation of topological sensitivity requires the limiting distribution of the scattered fieldp a over the measurement surface obs as a → 0. This problem is addressed next.

Scattered pressure field due to vanishing solid obstacle
To obtain a computable expression for topological sensitivity, consider the scattering problem (3) for a vanishing obstacle the sampling point is separated from the support of acoustic sources used to illuminate the obstacle), and B is a ball of unit radius centred at the origin of the reference Cartesian system. By virtue of the integral representation (7), the governing Equation (1a) and the fluid-solid interfacial conditions over S a = ∂ B a as in (3e), the scattered pressure field in − a = \ B a can be written as: where n is oriented towards the interior of B a , index "a" is used to indicate quantities that are dependent on the vanishing obstacle size, and the second convolution operator, indicated in boldface, signifies the "vector" convolution defined (using Einstein summation convention) as in terms of vectors g and h with respective Cartesian components g i and h i .

Limiting analysis
For any given x o ∈ \ E , x ∈ \ {x o } and ξ ∈ S a , one may expandp F (ξ , t), ∇ p F (ξ , t), G(ξ , x, t) and ∇ ξĜ (ξ , x, t) in Taylor series about ξ = x o as where ∇∇ p denotes the gradient of ∇ p (whereby it is a second-order tensor with compo- (17), the leading behaviour of (16) as a = |ξ − x o | → 0 can then be written as: On the basis of (3e), the fact that S a is a sphere (whereby ξ ∈ S a can be expressed as ξ = x o − a n), and the Taylor expansion of p F over S a , i.e.
where the operator ":" denotes double dot product [44] (e.g. the double dot product of two second-order tensors A and B is a scalar given by A i j B i j ), one finds that By virtue of (20) and the balance of linear momentum (3b), (18) reduces tõ Further evaluation ofp a from (21) requires the acceleration fieldü a (ξ , t) in the vanishing obstacle B a .
As shown in [45], the governing Equation (3b) for B = B a corresponds to a regularized boundary integral equation: where u k (ξ , x, t) is the elastodynamic Green's function for an infinite solid R 3 with elastic parameters (κ * , μ * ) and mass density ρ * , defined as the displacement at ξ ∈ R 3 due to a time-impulsive point force δ(ξ − x)δ(t)e k acting at x ∈ R 3 in the kth direction (k = 1, 2, 3), t k = −(C * : ∇u k ) · n is the associated traction, and − signifies integral in the Cauchy principle value sense.
To expose the limiting behaviour of (22) as a → 0, the following scaled coordinates, are introduced and the leading asymptotic behaviours of u a and t a are decomposed as where u 0 (t) denotes the rigid body motion of B a and S = ∂B is the surface of the unit ball B. As implied by (24), the elastodynamic fields u a and t a in the vanishing solid obstacle are decomposed into two parts. The first part denotes the regular term of O (1). It corresponds to a simplified fluid-solid interaction problem where a rigid sphere moves in an inviscid fluid. The second part signifies the infinitesimal term of O(a) as a goes to 0. With further asymptotic analysis, the second part turns out to be associated with another simplified fluid-solid interaction problem that will be found out later. As for the first simplified fluid-solid interaction problem, considering B a (a → 0) as an infinitesimal particle in the fluid, its rigid body accelerationü 0 can be calculated using the concept of virtual added mass in fluid mechanics. [46] For a sphere in an inviscid fluid, the added mass is equal to one half of the displaced mass. So that the rigid body acceleration of B a isü which makes use of the balance of linear momentum in the obstacle-free fluid, whereü F denotes the free field acceleration in the reference fluid medium. On the basis of (23) and the singular behaviour of the Stokes fundamental solution, [45] it can be shown [29] that whereū k is the elastostatic Kelvin's fundamental solution for an infinite solid R 3 with elastic parameters (κ * , μ * ) of B a ,t k = −(C * : ∇u k ) · n is the associated traction, and (24) and (27) into (22) yields which makes use of S a = a 2 S. On dropping the (underlined) higher order term, the integral Equation (28) turns out to be associated with an elastostatic transmission problem where a unit ball B with elastic parameters κ * and μ * is subject to a uniform pressure p F (x o , t) on its boundary. Therefore, the second part of the u a decomposition in (24) is associated with the simplified fluid-solid problem that describes the deformation of a sphere subject to a uniform pressure applied by the surrounding fluid. The solution to this problem can be computed as where c and ω are, respectively, the rigid body translation and rotation as permitted by the Neumann boundary condition of this problem.
With (24), (25) and (29) in place, the two volume integrals in expression (21) for the scattered fieldp a can be evaluated as Note that the rigid body accelerationü 0 (t) has a zero divergence with respect to ξ ∈ B a . From (30), the asymptotic representation (21) forp a finally reduces tõ In case of time-harmonic scattering, (31) can be simplified as whereĜ(ξ , x, ω) and ∇ ξĜ (ξ , x, ω) are the respective counterparts ofĜ(ξ , x, t) and ∇ ξĜ (ξ , x, t) for time-harmonic excitation. Owning to the linearity of the forward scattering problem with respect to the excitation, (31) and (32) are interchangeable via Fourier transform and its inversion, i.e.
where f (t) is the excitation wavelet.
Remark Here, it is important to be aware of that the equivalence between (31) and (32) is formal as long as the conditions for invertibility of Fourier transforms [47] are met, which is, however, often troublesome in wave equations.

Discussion and validation
From (31) and (32), one may note that the limiting behaviour of the scattered pressure field p a does not involve all the material parameters of the trial solid obstacle. Indeed, it depends on the densities and bulk moduli of the background and the trial solid obstacle. Nonetheless, it is important to observe that expressions (31) and (32) are identical to their counterparts for a vanishing acoustic obstacle with mass density ρ * and bulk modulus κ * . [7,16] This implies that in terms of leading-order behaviour, an acoustic source in the surrounding fluid cannot induce shear in a homogeneous isotropic solid obstacle.
Because of the equivalence of expressions (31) and (32), the time-harmonic expression (32) is validated for simplicity. First, the o(a 3 ) validity of (32) for configurations with both ρ * = ρ and κ * = κ is tested. Consider a solid ball B a of radius a, located at x o = (0.5d, −1.0d, −9.0d) in the full space R 3 characterized by mass density ρ and bulk modulus κ, where d is the reference length. The ball is elastic and isotropic with the same mass density ρ * = ρ and bulk modulus κ * = κ as the background. The Poisson's ratio ν * of the ball takes 0.25 (the shear modulus μ * = 3κ * (1 − 2ν * )/2(1 + ν * )). A time-harmonic excitation with frequency ω = 2c/d is applied at (−10.0d, −1.0d, −9.0d). The induced scattered pressure fieldP a is computed at a number of points x ∈ \B a using a boundary element method [48] tailored for fluid-solid coupling, [49] while the free field and Green's functions are calculated analytically. Figure 2 shows the variation of |P a |/a 3 with respect to the obstacle size a/d at a few field points. As a/d → 0, |P a |/a 3 goes to zero rapidly. A number of similar configurations have been tested and the same trend is observed, indicating that the scattered pressure fieldP a is of smaller order than a 3 , which thus verifies the o(a 3 ) term in expression (32). Next, the coefficients of the two O(a 3 ) terms in (32) are validated. For obstacles with either ρ * = ρ or κ * = κ, (32) can be rewritten as where Q 1 and Q 2 denote the coefficients related to κ * and ρ * , respectively. Note that as a → 0, the values of Q 1 and Q 2 do not depend on the location of the obstacle x o , the measurement point x and the excitation frequency ω. To facilitate the validation of (34), let denote the counterparts of Q 1 and Q 2 for small a, where the scattered fieldP a (x, ω) can be computed numerically. The correctness of (34) can be validated if for any Let the solid ball B a in the testing configuration take material parameters different from those of the background. Figure 3 shows the |Q a versus a/d curves at a field point x = x 1 for different sets of (κ * , ν * ) and (ρ * , ν * ). Note that ν * = 0.5 indicates that the scatterer is acoustic, i.e. having zero shear modulus. As can be seen from the display, both Q a 1 /Q 1 and Q a 2 /Q 2 tend to 1 as a goes to 0. A number of configurations with different x o , x and ω values have been tested and the same trend is observed, which thus validates the correctness of the coefficients. On comparing the last two curves in each plot, it can be seen that when the mass density ρ * and bulk modulus κ * of the obstacle are kept the same, the change of shear modulus μ * (Poisson's ratio ν * ) does affect the scattered field, however, its influence should be of higher order and therefore the coefficients are expressed in terms of the mass density and bulk modulus.

Formulas for topological sensitivity
By analogy to the developments in [7,16], two alternative expressions for the topological sensitivity can be established. On substituting (31) into (15) and taking h(a) = a 3 |B| = 4πa 3 /3, the "direct" topological sensitivity formula is expressed as: in terms of the Green's functionĜ for the reference domain .
In case that such Green's function is unavailable, (38) can be formulated alternatively using a so-called adjoint field. The adjoint field, p , is defined as the response of the reference domain due to virtual excitation E (ξ , t) = ∂ϕ/∂p( p F , p obs , ξ , T −t) applied over the measurement surface obs . Subject to applicable radiation condition, this adjoint field solves the boundary-value problem where n is the unit normal on obs and [[g]] = lim →0 g(ξ + n ) − g(ξ − n ). The boundary conditions (39b) and (39c) indicate that across obs , the pressure is continuous but its normal derivative is discontinuous. Such discontinuity is to be interpreted in the sense of a single-layer potential and treated as a distribution of acoustic sources over obs .
Here, it should be noted that the homogeneous boundary condition in (39d) is a deliberate choice to allow recognization in (38) the integral representation of the adjoint field. For the least-squares distance function in (9), the virtual excitation is given by the time-reversed data residual, i.e.
On the basis of (39), the adjoint field and its gradient can be expressed in terms of the Green's functionĜ as: By virtue of this result, change of integration variables and commutativity of the convolution operator, (38) can be rewritten in terms of the adjoint field as: which also makes use of the Green's function reciprocity, i.e.Ĝ (ξ , x, t) =Ĝ(x, ξ , t). Without the requirement of a Green's function, such adjoint-field topological sensitivity formula (41) can be applied to homogeneous reference domains of arbitrary shape and boundary condition, and in this case the adjoint field p can be computed numerically via a suitable (e.g. boundary element) solution scheme. An extension of these developments to heterogeneous reference domains may be established as those in [15,30]; such a generalization is, however, not addressed here.
For inverse scattering problems with time-harmonic excitations, (38) and its adjoint field counterpart (41) can be simplified as (42) and From (38) and (41)- (43), the computation of topological sensitivity does not involve the shear modulus μ * of the vanishing obstacle; therefore, let m = (ρ * , κ * ) from here on for simplicity.
Remark #1 As derived and validated in Section 3.2, the limiting behaviours of the transient and time-harmonic scattered pressure fieldsp a (x, t) andP a (x, ω) caused by an infinitesimal solid obstacle in a fluid background medium are identical to their counterparts for a vanishing acoustic obstacle with the same mass density and bulk modulus, and therefore the resulting direct and adjoint field topological sensitivity formulas (38) and (41)-(43) for the inverse acoustic scattering of solid obstacles are also the same as their counterparts for the acoustic/acoustic interaction, see e.g. [7] and [16] for the time-harmonic and timedependent cases, respectively, which again implies that the shear wave in the vanishing obstacle plays at most a second-order role.
Remark #2 As for the elastic/elastic interaction, the topological sensitivity formulas for time-harmonic and transient excitations can be found in e.g. [28] and [29]. For comparison, the direct time-harmonic TS formula [28] is given here, where f denotes the time-harmonic force distribution, u F and σ F are the free-field displacement vector and stress tensor, u obs denotes the displacement measurements, e k is the unit vector in the kth coordinate direction,û k (x o , ξ , ω) andσ k (x o , ξ , ω) constitute the timeharmonic elastodynamic Green's function for the reference elastic medium by denoting the respective displacement vector and stress tensor at x o ∈ due to a unit point force acting at ξ ∈ in the kth direction, and A is a fourth-order tensor depending on the shape B and material parameters m of the infinitesimal trial obstacle. The corresponding TS formula (42) for solid obstacle in fluid background takes a similar form as the expression (44) for elastic/elastic scattering. Actually, the case of solid obstacle in fluid background can be considered as a particular case of elastic/elastic scattering where the shear modulus of the reference medium is zero. However, we did not find a way to particularize the results for elastic/elastic scattering to this particular case.

Numerical examples
In previous studies [7,[28][29][30] involving acoustic scattering by acoustic scatterers or elastic scattering by solid scatterers, the topological sensitivity has been shown to be an effective tool for preliminary obstacle reconstruction and material characterization through an assembly of sampling points with pronounced negative TS values and a pointwise identification of "optimal" obstacle properties that minimize the TS at each sampling location. In what follows, these possibilities in terms of acoustic scattering of solid obstacles are explored using the example of a fluid layer containing penetrable solid obstacles.

Testing configuration
With reference to Figure 4, consider the synthetic testing configuration where a solid obstacle B true is embedded in a homogeneous fluid layer = {ξ | −20d < ξ 3 < 0} with Dirichlet top surface T = {ξ | ξ 3 = 0} and Neumann bottom surface B = {ξ | ξ 3 = −20d}, where d is the reference length. The obstacle is a ball of radius 0.5d centred at (0.5d, −1.0d, −9.0d). The background acoustic medium is characterized by mass density ρ and bulk modulus κ. The solid obstacle B true is elastic, homogeneous and isotropic with mass density ρ true = ρ/2, bulk modulus κ true = 2κ and Poisson's ratio ν true = 0.25. Such material parameters are chosen to mimic an obstacle without acoustic impedance (defined as ρc) contrast that is usually difficult for conventional reflection-based methods to detect.
To expose the hidden obstacle, the problem domain is sequentially illuminated by M = 25 point sources equally spaced on a source grid src = {ξ | ξ 1 = −10.0d, ξ 2 = −8d : 4d : 8d, ξ 3 = −18d : 4d : −2d} and for each source location x j , j = 1, 2, . . . , M the induced pressure p obs is monitored by all the receivers on a measurement grid obs = {ξ | ξ 1 = 10.0d, ξ 2 = −10d : 2d : 10d, ξ 3 = −18d : 2d : −2d}, as shown in Figure 4. For such multi-source illumination, the topological sensitivity can be computed as where T j is the single-source topological sensitivity, obtained by applying corresponding topological sensitivity formulas (38) or (41)- (43) to the testing configuration with only the jth point source E j . Again, the synthetic observations p obs are computed with the aid of a boundary element method [48] tailored for fluid-solid coupling, [49] while the free field p F and the adjoint field p are calculated using the Green's function (5), truncated at |n| ≤ N = 9. Note that in this case the topological sensitivity can also be computed directly using the "direct" formulas (38) and (42) since the Green's function is available. Furthermore, it should be noted that the topological sensitivity formulas are independent of the numerical implementation used to compute the wave fields required.

Time-harmonic excitation
For the time-harmonic example, the synthetic measurements are generated at an excitation frequency where λ is the acoustic wavelength in the background medium. Here, the selected frequency falls into the so-called resonance region, characterized by wavelengths comparable to or larger than the characteristic dimension of the scatterer. [50] Assuming that m = m true = (ρ true , κ true ) for the vanishing obstacle, the topological sensitivity values T (x o , m true ; E) over the three cross-sections passing through the centre of the obstacle (Figure 4) are calculated using (43). Note that a uniform spatial weighting function W ξ (ξ ) = 1, ξ ∈ obs is used, i.e. all the measurement points are equally weighted. The spacing between sampling points is 0.4d. Figure 5 plots the sectional distributions of T (x o , m true ; E), where the location and approximate geometry of the obstacle (white circle in the plot) are indicated by pronounced negative TS values. From Figure 5, it can also be seen that the reconstructed obstacle is elongated in ξ 1 -direction due to limited aperture.
In practice, experimental data always involves noise. To test the sensitivity of the TS imaging method to measurement noise, random noise is added to the synthetic pressure data p obs and in this case p obs = p true + p noise . As an illustration, the topological sensitivity distributions for the same testing configuration with uniformly distributed random measurement errors in the range of ±5%, ±10% and ±20% are computed, respectively. Figure 6 displays the ξ 1 = 0.5d and ξ 2 = −1.0d TS distributions with different levels of  measurement noise (the quality of TS distribution at ξ 3 = −9.0d is in between and thus not displayed). On comparing these diagrams to their no-noise counterparts in Figure 5, it is seen that no significant image deterioration is observed for 5% measurement noise, the TS distributions with 10% noise are somewhat distorted but still acceptable, whereas the TS results with 20% noise are noticeably distorted especially for the cross-section ξ 2 = −1.0d. From Figure 6, it can also been seen that the TS distribution at ξ 2 = −1.0d is more sensitive to noise due to limited aperture. As discussed in Section 3.2.2, the shear modulus μ true of the solid obstacle will affect the scattered acoustic field, i.e. p obs depends on μ true , and thus the topological sensitivity depends on μ true . To investigate the influence of μ true , the topological sensitivity for obstacles with the same mass density ρ true = ρ/2 and bulk modulus κ true = 2κ but different Poisson's ratio ν true (shear modulus μ true ) is computed. The sectional distributions of TS for obstacles with different shear moduli are compared in Figure 7. The TS values are very close and the TS distributions are almost the same, which indicates that the shear modulus of the solid obstacle has a very small influence on the topological sensitivity when its mass density and bulk modulus are kept fixed.
Remark From another point of view, the pressure fields p obs computed for the fluid-solid problems with 0 < ν true < 0.5 may be considered as noisy observation data of p true ν=0.5 (solution of the acoustic problem) if their differences are less than 5%. Here, the mean relative errors when comparing the solutions for ν true = 0.05 − 0.45 with an interval of 0.05 with that of the acoustic case p true ν=0.5 are computed. It turns out that the mean relative errors are between 0.3% and 2% although the maximum relative error can reach 35%. Different sets of ρ true /ρ, κ true /κ parameters for selected values of ν true = 0.5 are also tested and the mean relative errors between p true ν =0.5 and p true ν=0.5 are less than 5%. These results also support the solid obstacle to acoustic obstacle simplification.

Material characterization
In a number of TS studies, [7,[28][29][30] the topological sensitivity has also been employed as a preliminary material characterization tool, where the topological sensitivity value at each sampling point is minimized with respect to the material parameters of the trial obstacle and the trial material parameters generating the most negative TS value are considered as an estimation of the relative (compared to the background) material parameters of the true obstacle. Here, the above T -distributions are generated assuming full prior knowledge about the material properties m true of the hidden obstacle, an information that is often not available. To explore the potential of topological sensitivity towards material characterization, the TS at the centre of the obstacle x c = (0.5d, −1.0d, −9.0d) is evaluated with varying trial material property κ * or ρ * while the other parameter takes a "predefined" fixed value. Figure 8 plots the variations of T (x c , m; E) with respect to κ * /κ and ρ * /ρ for two material configurations, (2κ, 0.5ρ) and (0.5κ, 0.5ρ). The Poisson's ratio of the obstacle  Owing to the revealing dependence of T (x o , m; E) on κ * as shown in Figure 8, a preliminary material characterization approach can be established through a pointwise minimization of T (x o , m; E) with respect to the trial bulk modulus κ * . For each sampling point, the optimized bulk modulus κ opt is given by where κ is the search interval for κ * and ρ * takes predefined value ρ pre . Consequently, the TS value corresponding to κ opt is denoted by With such an approach, simultaneous geometric and material characterization of hidden obstacles can be accomplished by: (1) Geometrically reconstruct the obstacle using an isosurface computed with respect to the global minimum, where 0 < η < 1 is a suitable threshold and is the sampled region; and (2) Estimate the bulk modulus of the obstacle using the point-optimal bulk modulus computed by As an illustration, the topological sensitivity for the example in Section 4.2 is minimized as κ * varies in the interval κ = {κ * | 0.6κ < κ * < 15κ} and ρ * conveniently takes the background mass density ρ as the predefined value. The sectional distributions of optimized topological sensitivity T (x o , m opt ; E) and normalized point-optimal bulk modulus κ est /κ are displayed in Figure 9, where the obstacle is reconstructed by the 80%-level isosurface (dark line in the plot). Here, the true obstacle is indicated by a red circle for better contrast. For the purpose of obstacle reconstruction, the TS distributions in Figure 9 obtained using variable m = m opt are clearly comparable to companion images with m = m true in Figure 5, regardless of the apparent differences between them. Furthermore, the obstacle is identified correctly from the distributions of κ est /κ as "stiffer" than the background medium, i.e. having the bulk modulus κ true > κ. In practice, such preliminary geometric and material information could be used as the initial "guess" for gradient-based, non-linear minimization approaches (e.g. [51]) or form the basis of iterative obstacle reconstruction schemes (e.g. [15,52]).

Transient excitation
To take advantage of multi-frequency illumination, a raised cosine wavelet [6] given by is used as the input signal, where t 0 denotes the offset of the peak load. On taking α = 0.85, f 0 = 1.0, and t 0 = 3π/c, the time-domain wavelet and its frequency spectrum are plotted in Figure 10. The temporal synthetic data p obs (ξ , t), ξ ∈ obs are still simulated using the frequency-domain boundary element code. [49] Time-harmonic simulations are performed over the frequency range ωd/c = 0 − 16 with a step of 0.08. To avoid zero frequency, zero is replaced by a small value 10 −5 . The resulting frequency-domain data are then multiplied by the frequency spectrum in Figure 10(b) and inverse Fourier-transformed to give the temporal signal.
With reference to (41), (45) and the material optimization method established above, the sectional distributions of optimized topological sensitivity T (x o , m opt ; E) on the three cross-sections are displayed in Figure 11. Note that uniform spatial and temporal weighting functions W ξ (ξ ) = 1, ξ ∈ obs and W t (t) = 1, t ∈ [0, T ] are used. The predefined mass density takes ρ pre = ρ and the search interval for the trial bulk modulus is the same as in the time-harmonic example. Similar to the results in Figure 9, the regions where T (x o , m opt ; E) takes pronounced negative values point towards the true obstacle and the obstacle is identified correctly as "stiffer" than the background medium because κ est /κ = 15 > 1 inside the reconstructed obstacle. In terms of locating the obstacle, transient excitation gives much clearer TS distributions.
To examine the influence of ρ pre and seek a better estimation of κ true , the true mass density of the obstacle ρ true is taken as the predefined mass density for TS computation. However, the results are very similar to their counterparts shown in Figure 11 (and thus not displayed here). In the reconstructed obstacle, the TS goes to negative as κ * /κ increases and does not reach a most negative value at κ * = κ true = 2κ, and thus cannot provide a better estimation other than "stiffer" or "softer".

Multiple obstacles
Next, the applicability of the foregoing developments to more complex situations is examined by a testing problem involving two solid obstacles with different material properties. The reference fluid layer and source-observation setup are the same as in previous examples. The first obstacle is located at (−0.5d, 1.4d, −8.2d) and has the same geometry and material properties as the obstacle in Figure 4. The second obstacle is a "softer" ellipsoidal solid object characterized by mass density 2ρ, bulk modulus 0.5κ and Poisson's ratio 0.25. It is centred at (−0.5d, −1.8d, −11.4d) and has semi-axes (0.5d, 0.6d, 0.4d) aligned with the reference Cartesian frame. The obstacles are illuminated by transient excitations in the form of raised cosine ( Figure 10). On taking the same predefined mass density and search interval for the trial bulk modulus, Figure 12 plots the sectional distributions of optimized topological sensitivity T (x o , m opt ; E) and normalized point-optimal bulk modulus κ est /κ on three selected crosssections, where both obstacles are detected and correctly identified as one "stiffer" and the other "softer".

Iterative 3D obstacle reconstruction
In a few recent studies [15,[52][53][54] dealing with acoustic and electrical impedance tomography problems, the topological sensitivity has been effectively employed to iteratively update the geometric approximations of hidden/buried obstacles. By analogy to these developments, an iterative 3D obstacle reconstruction approach based on the topological sensitivity established for acoustic scattering of solid obstacles is established and verified in the following analysis.

Iterative algorithm
To facilitate the iterative application, the adjoint-field topological sensitivity formulations (41) and (43) are conveniently adjusted to acoustic reference domains containing preexisting solid inhomogeneities/obstacles. In this case, both the free and adjoint fields in (41) can be redefined with respect to the reference domain with known solid obstacles where u i and t i denote the displacement and surface traction for the ith obstacle, respectively. Similarly, the free and adjoint fields in (43) can be redefined by the counterparts of (51) and (52) with time-harmonic excitations.
Remark Note that the topological sensitivity formulas (41) and (43) with free and adjoint fields computed using (51) and (52) or their time-harmonic counterparts are under the assumption that the sampling points are outside existing solid obstacles. In [15], the topological sensitivity formula for sampling points inside existing acoustic scatterers has been derived. However, the situation here is different in that the existing obstacles are solids. For a sampling point x o inside an existing solid obstacle i , the free and adjoint fields defined by (51) and (52) would be in terms of displacements u F i , u i and their associated stresses. Thus, the TS formulas (41) and (43) obtained for a vanishing obstacle in a fluid background and expressed in terms of pressure fields p F and p would be inapplicable. In this case, the infinitesimal trial solid defect B a should be introduced into a solid reference body surrounded by an acoustic medium. With some generalization, the topological sensitivity formulas established for inverse elastic scattering such as those in [28][29][30] may be applicable to this situation, which is however beyond the scope of this study.
With topological sensitivity formulas for reference domains with known solid obstacles in place, the obstacles reconstructed at a certain step can be considered as pre-existing obstacles for the computation of topological sensitivity at the next step, which thus allows the implementation of an iterative procedure. To focus on the performance of topological sensitivity for obstacle reconstruction, here the material properties (mass density, bulk modulus and shear modulus) of hidden obstacles are assumed to be known beforehand. The 3D distribution of topological sensitivity, calculated using proper TS formulas, is used to iteratively update the geometry of hidden obstacles. The procedures, similar to those in [15,[52][53][54], are described in the following.
• Find the initial reconstructions of hidden obstacles 1 i using the distribution of topological sensitivity as illustrated in Sections 4.3 and 4.4. At this step, a tight threshold value η between 0.95 and 1.0 is applied to avoid overestimation of the obstacle and the condition that • Solve the forward (free field) and adjoint problems with

Numerical experiments
The numerical example presented in Section 4.2 is used for a first illustration. The initial 3D reconstruction based on TS images displayed in Figure 5 is shown in Figure 13(a), where the reconstructed obstacle and its projections on the three coordinate planes are in red and the projections of the true obstacle are represented with blue solid lines. The ball-shaped obstacle is reconstructed as an ellipsoid, elongated in the ξ 1 −direction. On taking this ellipsoid as a pre-existing obstacle, the topological sensitivity in the exterior of the ellipsoid is computed using the adjoint-field TS formula with redefined free and adjoint fields. The sectional distributions of T (x o , m true ; E) and the 3D reconstruction at the second iteration are displayed in Figures 14 and 13(b), respectively, where the reconstructed obstacle is expanded in the ξ 2 − and ξ 3 −directions. The 3D reconstruction results for the third and fourth iterations are shown in Figure 13(c) and (d). The location and geometry of the 3D reconstruction have stopped to improve although they are still not very close to those of the true obstacle. The reconstructed obstacle is considerably shifted and deformed in the ξ 1 −direction, due to limited aperture in this dimension.
To improve the 3D reconstruction, the testing configuration shown in Figure 4 is supplemented by an opposite set of source and receiver grids, located at ξ 1 = 10d and −10d, respectively. The topological sensitivity values, computed using both sets of source and receiver grids, are added together for 3D obstacle reconstruction. The obtained sequence  of 3D approximations is displayed in Figure 15, where the location and geometry of the reconstructed obstacle are significantly improved after only seven iterations. The reconstructed obstacle can be considered as an ellipsoid centred at (0.51d, −0.99d, −8.94d), with semi-axes (0.51d, 0.51d, 0.47d). These estimations are very close to their true values, i.e. centre (0.5d, −1.0d, −9.0d) and radius 0.5d. Considering the limited aperture of the source and receiver grids, the 3D obstacle reconstruction results are satisfying.

Preliminary experimental application
Notwithstanding the progress made on the theoretical formulation of the topological sensitivity method, there are only a few number of studies that demonstrate its usefulness through experimental data. [55][56][57] In what follows, the effectiveness of the topological sensitivity formulas obtained above is examined using preliminary experimental data.

Experimental setup
As shown in Figure 16, respectively, and thus its bulk modulus is κ true = 7.10 GPa. Two RESON hydrophones with model numbers TC4033 and TC4013 are used as the transducer and receiver, respectively. A metal frame is designed to hold the hydrophones and adjust their positions. An Agilent 33220A waveform generator is used to generate the input signal for the transducer and the received signal is recorded using a Sigma S60-4 oscilloscope. The hydrophones used here are sensitive in the frequency range of 5-170 kHz. In this experiment, only one receiver is used since receivers also act as obstacles in water.
To obtain meaningful and high-quality data, special attentions are paid to the following factors. (1) Signal noise. To reduce random signal noise, a number of measurements of exactly the same experiment are averaged as the measured output and the recovery period for repeated tests should be long enough to guarantee a quiescent past. (2) Point source and receiver. To satisfy the point source and receiver assumption made in the numerical calculation, the two hydrophones (although specified as omnidirectional) should be separated by a certain distance and here the distance is around 20 cm. The way how the hydrophones are held also affects the pressure field and it is observed that the influence of the stems is relatively small when they are parallel with the hydrophones (Figure 17). Furthermore, the underwater parts of the metal frame are as far as possible from the hydrophones and the stems are thin to minimize their disturbance. (3) Waveform and duration of the input signal. To compare with the computed result, the received signal is modified by keeping only the direct arrival with frequency components in the range of 5-170 kHz where the frequency responses of the hydrophones are available. However, a significant difference between simulated and measured responses is still observed when the input signal takes the form of a half-sine pulse of 10 µs. A possible reason could be that the sensitivity of a hydrophone is tested using long signals of each frequency. To avoid this problem, an oscillating signal of 4 cycles of 100 kHz sine wave is used and the calculated result matches the modified received signal very well. Therefore, such an input signal is selected in this study. (4) Reflections. As shown in Figure 16, a plastic layer is used to make the wood box waterproof, and thus the surface is not really a "rigid" surface. To consider the bottom surface B to be flat (parallel to the water surface) and satisfy the Neumann boundary condition, a steel plate is placed at the bottom of the box. The impedance of water is Z water = ρc ≈ 1.48×10 6 kg m −2 s −1 , while the impedances of steel and air are, respectively, Z steel ≈ 4.63×10 7 kg m −2 s −1 and Z air ≈ 4.10×10 2 kg m −2 s −1 .
Owing to their significant differences in impedance, it is reasonable to assume that the reflections from the top and bottom surfaces are ideal.  Figure 18 shows the top and side views of the experimental setup. The x, y and z axes of the coordinate system are aligned with the length, width and height of the box, respectively. The source S(S0) and receiver R have the same y coordinate. Owing to the large size of the box, the first four reflections come from the top and bottom surfaces. The effects of these four reflections are simulated by four "virtue" sources S1-S4 (Figure 18), where the "+" or "−" denotes the sign of a virtual source. Figure 19 shows the originally received signal when there is no obstacle. The direct arrival time is denoted by t0 while the arrival times of the first four reflections are denoted by t1-t4, respectively. From Figures 18 and 19, the following system of equations can be obtained, where ts is the starting time of the input signal. With (53) and the arrival times obtained from Figure 19, the unknown distances d0-d4 can thus be calculated with a precision of 1/10 of the wavelength, i.e. With reference to the coordinate system shown in Figure 18, the coordinates of the source and receiver are (0, 0, 20.02 cm) and (21.54 cm, 0, 14.12 cm), respectively. Figure 20 plots the observed free-field pressure p F obs (t) without the obstacle and the observed pressure p obs (t) when the obstacle is present, where the recorded signals are converted to pressure using the frequency response of the receiver and only the first four reflections are considered. Since the pressure to voltage response of the receiver is available only in  With reference to (41), the evaluation of topological sensitivity requires the free field and adjoint field in the reference domain. To this end, here the direct arrival part p F 0 (t) of the observed free-field pressure p F obs (t) at the receiver and the transfer function of water are used to back-calculate the pressure wave P(t) sent from the source S, namely

Identification of the obstacle
where F is the Fourier transform operator, D denotes the distance between the source and receiver, f is the frequency and c is the wave speed in water. With P(t) in place, the free-field pressure (including the reflections) at any point of the problem domain can be computed using the transfer function of water. The simulated free-field pressure at the receiver is also shown in Figure 20. It is comparable to the observed free-field pressure but the reflections are smaller in magnitude. Given the observed pressure waves p F obs (t) and p obs (t), the adjoint field pressure p , defined as the response of the reference domain due to the time-reversal data residual p F obs (T − t) − p obs (T − t) applied at the receiver location, can be calculated directly using the transfer function of water. Note that the virtual sources due to reflections from the top and bottom surfaces should be accounted for as before.
With the treatments described above and assuming the true material properties, the topological sensitivity T (x o , m true ; E) is computed over a sampling grid {(x, y) | x = 0 : 1 : 20 cm, y = −10 : 1 : 10 cm} at z = 10 cm, with and without considering the four reflections from the top and bottom surfaces. Figure 21 displays the resulting TS distributions. They are symmetric with respect to the y axis owing to the symmetry of the experimental setup. Due to the fact that the reflections are not simulated very well ( Figure  20), the TS distribution calculated without considering the reflections (using only the first arrival) gives better results, from which, the PVC tube (indicated by white circles) is clearly located. TS distributions considering different number of reflections (1, 2 or 3) are also calculated, however, the results are worse than that without reflections. Therefore, in the following analysis, the topological sensitivity is calculated using only the first arrival. Figure 20. Simulated free-field pressure, observed free-field pressure and observed pressure when the obstacle is present. Figure 22 plots the sectional distributions of topological sensitivity at two more crosssections z = 8 cm and z = 12 cm, where the location of the PVC tube is also approximately identified. Note that at z = 8 cm, the topological sensitivity is positive at the centre of the PVC tube x c = (15 cm, 0 cm, 8 cm) and takes the most negative value at point x = (17.3 cm, 0 cm, 8 cm), denoted by a red " " on the plot.
The variations of topological sensitivity with respect to κ * /κ at the centre of the PVC tube on z = 10 cm and z = 12 cm and at x on z = 8cm are plotted in Figure 23. The topological sensitivity goes to negative as κ * /κ increases, which correctly indicates that the obstacle is "stiffer" than the background. From the plots, it can also be seen that the quality of experimental TS images is not as good as the numerical results, due to the discrepancy between simulations and experiments. However, the potential of topological sensitivity for preliminary geometric reconstruction and material characterization regarding the bulk modulus κ * is clearly demonstrated.
To further improve the geometric approximation of the tube, the iterative 3D obstacle reconstruction approach developed in Section 5 is tried. As shown in Figures 21 and 22,   the location of the tube on each x − y cross-section is indicated by relatively pronounced negative TS values, however, the magnitude of TS values on cross-sections at different z locations varies significantly, which makes 3D reconstruction of the tube with a threshold value unsuccessful even with the help of the iterative method. The main reason should be that there is only one source and one receiver. A separate parameter analysis with synthetic data (provided as supplemental material) shows that only one source and one receiver are insufficient to yield a good 3D reconstruction. Therefore, source and receiver grids covering a larger region should be used for future experimental investigations.

Conclusions
In this study, the concept of topological sensitivity is particularized to deal with the inverse acoustic scattering of penetrable solid obstacles. The asymptotic behaviour of the scattered acoustic field caused by a spherical solid obstacle is obtained based on solutions of simplified fluid-solid interaction problems resulting from preliminary limiting analysis of the boundary integral representation. With this result, the equivalence of acoustic scattering by acoustic and solid scatterers in terms of asymptotic behaviour is observed and numerically validated. Accordingly, the direct and adjoint-field topological sensitivity formulas for transient and time-harmonic excitations are expressed explicitly in terms of either the Green's functions or acoustic wave fields computed exclusively for the reference domain. The usefulness of the proposed method for providing both geometric and material information of hidden obstacles is illustrated through numerical examples and an exploratory experimental study. To further improve the geometric reconstruction of hidden obstacles, the topological sensitivity formulas are adapted to deal with acoustic reference domains with known solid obstacles and based on which an iterative 3D obstacle reconstruction approach is established. The performance of this iterative method is demonstrated with synthetic data. The details of the experimental study are described as references for further verification. Although not as good as the numerical results, the result of preliminary application to experimental data is promising and clearly demonstrates the potential of topological sensitivity for obstacle reconstruction and material characterization regarding the bulk modulus.