Forward and inverse acoustic scattering problems involving the mass density

TL;DR

This paper investigates the inverse problem of retrieving the mass density or wavespeed of a layered medium in 2D acoustic scattering, using explicit solutions and perturbation techniques to improve retrieval accuracy.

Contribution

It introduces a perturbation-based series solution for the inverse problem, enabling more accurate retrieval of layer properties in acoustic scattering.

Findings

Perturbation series effectively estimates mass density contrast.

Explicit solutions facilitate inverse problem analysis.

Numerical results confirm the method's accuracy.

Abstract

This investigation is concerned with the 2D acoustic scattering problem of a plane wave propagating in a non-lossy fluid host and soliciting a linear, isotropic, macroscopically-homogeneous, lossy, flat-plane layer in which the mass density and wavespeed are different from those of the host. The focus is on the inverse problem of the retrieval of either the layer mass density or the real part of the layer wavespeed. The data is the transmitted pressure field, obtained by simulation (resolution of the forward problem) in exact, explicit form via separation of variables. Another form of this solution, which is exact and more explicit in terms of the mass-density contrast (between the host and layer), is obtained by a domain-integral method. A perturbation technique enables this solution to be cast as a series of powers of the mass density contrast, the first three terms of which are employed as the trial models in the treatment of the inverse problem. The aptitude of these models to retrieve the mass density contrast and real part of the layer wavespeed is demonstrated both theoretically and numerically. Keywords: 2D acoustics, forward scattering, inverse scattering, domain integral equation, constant mass density assumption, small density contrast, retrieval accuracy

Equations

$${\nabla }_{\mathbf{x}}\cdot {\nabla }_{\mathbf{x}}p(\mathbf{x},t)-\frac{1}{c(\mathbf{x}{)}^2}\frac{{\partial }^{2}p(\mathbf{x},t)}{\partial t^2}-\frac{{\nabla }_{\mathbf{x}}\rho (\mathbf{x})}{\rho (\mathbf{x})}\cdot {\nabla }_{\mathbf{x}}p(\mathbf{x},t)=\rho (\mathbf{x}){\nabla }_{\mathbf{x}}\cdot \mathbf{f}(\mathbf{x},t)=0;\forall \mathbf{x}\in {\mathbb{R}}^2,$$

$$\frac{1}{c(\mathbf{x}{)}^2}:=\frac{\rho (\mathbf{x})}{K(\mathbf{x})}=\rho (\mathbf{x})\kappa (\mathbf{x}),$$

$$p(\mathbf{x},t)={\int }_{-\infty }^{\infty }p(\mathbf{x},\omega )\exp (-i\omega t)d\omega ,$$

$$\mathbf{f}(\mathbf{x},t)={\int }_{-\infty }^{\infty }\mathbf{f}(\mathbf{x},\omega )\exp (-i\omega t)d\omega ,$$

$${\nabla }_{\mathbf{x}}\cdot {\nabla }_{\mathbf{x}}p(\mathbf{x},\omega )+k^2(\mathbf{x},\omega )p(\mathbf{x},\omega )-\frac{{\nabla }_{\mathbf{x}}\rho (\mathbf{x},\omega )}{\rho (\mathbf{x},\omega )}\cdot {\nabla }_{\mathbf{x}}p(\mathbf{x},\omega )=\rho (\mathbf{x},\omega ){\nabla }_{\mathbf{x}}\cdot \mathbf{f}(\mathbf{x},\omega ;\forall \mathbf{x}\in {\mathbb{R}}^2,$$

$$k(\mathbf{x},\omega ):=\frac{\omega }{c(\mathbf{x},\omega )}$$

$$p(\mathbf{x},\omega)=p^{r}(\mathbf{x},\omega)+\\ \int_{\mathbb{R}^{2}}G^{[0]}(\mathbf{x};\mathbf{x^{\prime}};\omega)\left(\left[\big{(}k(\mathbf{x^{\prime}},\omega)\big{)}^{2}-\big{(}k^{[0]}\big{)}^{2}\right]p(\mathbf{x^{\prime}},\omega)-\frac{\nabla_{\mathbf{x^{\prime}}}\rho(\mathbf{x^{\prime}},\omega)}{\rho(\mathbf{x^{\prime}},\omega)}\cdot\nabla_{\mathbf{x^{\prime}}}p(\mathbf{x^{\prime}},\omega)\right)d\Omega(\mathbf{x^{\prime}})~{}~{};~{}~{}\forall~{}\mathbf{x}\in\mathbb{R}^{2}~{},$$

$${\nabla }_{\mathbf{x}}\cdot {\nabla }_{\mathbf{x}}{G}^{[0]}(\mathbf{x};{\mathbf{x}}^{\prime };\omega )+{\left(k^{[0]}\right)}^2{G}^{[0]}(\mathbf{x};{\mathbf{x}}^{\prime };\omega )=-\delta (\mathbf{x}-{\mathbf{x}}^{\prime });\forall \mathbf{x},{\mathbf{x}}^{\prime }\in {\mathbb{R}}^2,$$

$${\nabla }_{\mathbf{x}}\rho (\mathbf{x},\omega )=0,{\nabla }_{\mathbf{x}}c(\mathbf{x},\omega )=0;\forall \mathbf{x}\in {\Omega }_1,$$

$$\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+{\left(k^{[j]}(\omega )\right)}^2\right)p^{[j]}(\mathbf{x},\omega )=0;\forall \mathbf{x}\in {\Omega }_j,j=0,1,2$$

$$p^i(\mathbf{x},\omega )=a^i(\omega )\exp (i{\mathbf{k}}^{\mathbf{i}}\cdot \mathbf{x})=a^i(\omega )\exp [i(k_x^{i[0]}x-k_y^{i[0]}y)],$$

$$p^d(\mathbf{x},\omega )\sim \text{outgoing waves};\Vert \mathbf{x}\Vert \to \infty .$$

$$p^{[0]}(x,0,\omega )-p^{[1]}(x,0,\omega )=0;\forall \mathbf{x}\in \mathbb{R},$$

$$\frac{1}{{\rho }^{[0]}}{p}_{,y}^{[0]}(x,0,\omega )-\frac{1}{{\rho }^{[1]}}{p}_{,y}^{[1]}(x,0,\omega )=0;\forall \mathbf{x}\in \mathbb{R},$$

$$p^{[1]}(x,-h,\omega )-p^{[2]}(x,-h,\omega )=0;\forall \mathbf{x}\in \mathbb{R},$$

$$\frac{1}{{\rho }^{[1]}}{p}_{,y}^{[1]}(x,0,\omega )-\frac{1}{{\rho }^{[0]}}{p}_{,y}^{[2]}(x,0,\omega )=0;\forall \mathbf{x}\in \mathbb{R}.$$

$$p^{d[j]}(\mathbf{x},\omega )=a^{[j]}(\omega )\exp [i(k_x^{i[0]}x-k_y^{i[j]}y)]+b^{[j]}(\omega )\exp [i(k_x^{i[0]}x+k_y^{i[j]}y)],j=0,1,2;$$

$$k_y^{i[j]}=\sqrt{(k^{[j]}{)}^2-(k_x^{i[0]}{)}^2};\Re (k_y^{i[j]}\ge 0,\Im (k_y^{i[j]}\ge 0;\omega \ge 0,$$

$$b^{[2]}=0,$$

$$\begin{array}[]{cc}\left(a^{[0]}+b^{[0]}\right)-\left(a^{[l]}-b^{[l]}\right)&=0\\ \frac{ik_{y}^{i[0]}}{\rho^{[0]}}\left(-a^{[0]}+b^{[0]}\right)-\frac{ik_{y}^{i[1]}}{\rho^{[1]}}\left(-a^{[1]}+b^{[1]}\right)&=0\\ \left(a^{[1]}\exp(ik_{y}^{i[1]}h)+b^{[1]}\right)\exp(-ik_{y}^{i[1]}h)-a^{[2]}\exp(ik_{y}^{i[2]}h)&=0\\ \frac{ik_{y}^{i[1]}}{\rho^{[1]}}\left(-a^{[1]}\exp(ik_{y}^{i[1]}h)+b^{[1]}\exp(-ik_{y}^{i[1]}h)\right)-\frac{ik_{y}^{i[2]}}{\rho^{[2]}}\left(-a^{[2]}\exp(ik_{y}^{i[2]}h)\right)&=0\end{array}~{},$$

$$b^{[0]}=a^{[0]}\left[\frac{\left(({\gamma }^{[0]}{)}^2-({\gamma }^{[1]}{)}^2\right)i\sin (k_y^{[1]}h)}{2{\gamma }^{[0]}{\gamma }^{[1]}\cos (k_y^{[1]}h)-\left(({\gamma }^{[0]}{)}^2+({\gamma }^{[1]}{)}^2\right)i\sin (k_y^{[1]}h)}\right],$$

$$a^{[1]}=a^{[0]}\left[\frac{{\gamma }^{[0]}\left({\gamma }^{[1]}+{\gamma }^{[0]}\right)\exp (-i{k}_y^{[1]}h)}{2{\gamma }^{[0]}{\gamma }^{[1]}\cos (k_y^{[1]}h)-\left(({\gamma }^{[0]}{)}^2+({\gamma }^{[1]}{)}^2\right)i\sin (k_y^{[1]}h)}\right],$$

$$b^{[1]}=a^{[0]}\left[\frac{{\gamma }^{[0]}\left({\gamma }^{[1]}-{\gamma }^{[0]}\right)\exp (i{k}_y^{[1]}h)}{2{\gamma }^{[0]}{\gamma }^{[1]}\cos (k_y^{[1]}h)-\left(({\gamma }^{[0]}{)}^2+({\gamma }^{[1]}{)}^2\right)i\sin (k_y^{[1]}h)}\right],$$

$$a^{[2]}=a^{[0]}\left[\frac{2{\gamma }^{[0]}{\gamma }^{[1]}\exp (-i{k}_y^{[0]}h)}{2{\gamma }^{[0]}{\gamma }^{[1]}\cos (k_y^{[1]}h)-\left(({\gamma }^{[0]}{)}^2+({\gamma }^{[1]}{)}^2\right)i\sin (k_y^{[1]}h)}\right],$$

$$G^{[0]}(\mathbf{x};{\mathbf{x}}^{\prime })=\frac{i}{4}{H}_0^{(1)}(k^{[0]}(\mathbf{x};{\mathbf{x}}^{\prime })),$$

$$G^{[0]}(\mathbf{x};{\mathbf{x}}^{\prime })=\frac{i}{4\pi }{\int }_{-\infty }^{\infty }\exp [i(k_x(x-x^{\prime })+k_y^{[0]}|y-y^{\prime }|)]\frac{d{k}_x}{k_y^{[0]}},$$

$$k_y^{[j]}=\sqrt{(k^{[j}{)}^2-(k_x{)}^2};\Re k_y^{[j]}\ge 0,\Im k_y^{[j]}\ge ;\omega \ge 0.$$

$$p(\mathbf{x})=p^{i}(\mathbf{x})+\int_{\mathbb{R}^{2}}G^{[0]}(\mathbf{x};\mathbf{x^{\prime}})\left[\big{(}k(\mathbf{x^{\prime}})\big{)}^{2}-\big{(}k^{0}\big{)}^{2}\right]p(\mathbf{x^{\prime}})d\Omega(\mathbf{x^{\prime}})~{}~{};~{}~{}\forall~{}\mathbf{x}\in\mathbb{R}^{2}~{},$$

$$k(\mathbf{x})=k(y),\rho (\mathbf{x})=\rho (y).$$

$$p(\mathbf{x})=p(y)\exp (i{k}_x^{i[0]}x)$$

Full text

Keywords: 2D acoustics, forward scattering, inverse scattering, domain integral equation, constant mass density assumption, small density contrast, retrieval accuracy

Abbreviated title: Role of mass density in forward and inverse acoustic scattering

Corresponding author: Armand Wirgin,

e-mail: wirgin@lma.cnrs-mrs.fr

1. 1 Introduction

2. 2 The forward-scattering problem

3. 2.1 Preliminaries

4. 2.2 Description of the problem of scattering of a plane wave by a layer

5. 2.3 The domain decomposition, separation-of-variables (DD-SOV) solution via the BPDE

6. 2.4 The domain integral (DI) solution

7. 2.4.1 The assumption of invariance with respect to $x$

8. 2.4.2 The assumption of piecewise-constant density with respect to y

9. 2.4.3 The assumption of constant wavespeed within the layer

10. 2.4.4 The expansion of the transmission coefficient in a series of powers of $\epsilon$

11. 2.5 Numerical results for the comparison of the exact $a^{[2]}$ with the zeroth, first and second-order (in terms of $\epsilon$) approximations thereof

12. 2.5.1 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the density contrast $\epsilon$

13. 2.5.2 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the incident angle $\theta^{i}$

14. 2.5.3 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the frequency $f$

15. 2.5.4 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the real part of the wavespeed $\Re\left(c^{[1]}\right)$ in the layer

16. 2.5.5 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the imaginary part of the wavespeed $\Im\left(c^{[1]}\right)$ in the layer

17. 3 The inverse-scattering problem

18. 3.1 Preliminaries

19. 3.2 Retrieval of the mass density contrast

20. 3.3 Numerical results for the retrieval of the mass density contrast

21. 3.3.1 First and second-order retrievals as a function of the frequency $f$

22. 3.3.2 First and second-order retrievals as a function of the actual density contrast $\epsilon$

23. 3.3.3 First and second-order retrievals as a function of the layer thickness $h$

24. 3.4 Retrieval of the real part of the layer wavespeed

25. 3.5 Numerical results for the retrieval of the real part of the layer wavespeed

26. 3.5.1 Retrieved real part of the layer wavespeed $C^{[1]^{\prime}}$ as a function of frequency $f$

27. 3.5.2 Retrieved real part of the layer wavespeed $C^{[1]^{\prime}}$ as a function of $c^{[1]}$

28. 3.5.3 Retrieved real part of the layer wavespeed $C^{[1]^{\prime}}$ as a function of $\epsilon$

29. 4 Conclusions

1 Introduction

Scattering is produced when a material obstacle is solicited by a wave and manifests itself by the transformation of the incident wave into several waves whose amplitudes, directions, and spectral properties are different from those of the solicitation. For instance, when the solicitation is a plane wave and the obstacle is a fluid-like, macroscopically homogeneous, isotropic flat-faced layer situated within another fluid of the same nature, but of different density and/or compressibility, the scattering manifests itself by the splitting of the incident plane wave into a reflected plane wave and a transmitted plane wave in the two half-spaces above and below the layer, accompanied by one upgoing and one downgoing plane wave within the layer itself.

The task, in forward-scattering problems, is to predict the characteristics of the scattered wavefield for a given solicitation, host medium, and obstacle. This task has occupied researchers for the last two centuries [63, 56, 65, 60, 81, 3, 21, 39, 40] (and continues to occupy them [51, 37, 48, 80, 94]) because the solicitation, and/or the host medium, and/or the obstacle can (and does, in real-world situations) have complicated characteristics. For instance, in the biomedical field, predicting the acoustic wavefield within a human body [71, 49], due to an ultrasonic transducer being applied on the surface of the body, is complicated because the body is a complicated medium/structure composed of many closely-spaced heterogeneous organs of different nature and size located within a host medium containing many small-scale heterogeneities, and the ultrasonic pulse-like solicitation wavefield emitted by the transducer is not simple either. A first, crude, but very-useful, approximation is to assume that the target organ is similar to a homogeneous, isotropic fluid-like layer, the wave soliciting this organ within the body is similar to a plane wave, and that it strikes the layer-like organ at normal incidence. It ensues that the scattering manifests itself by a backward-reflected wave redirected towards the skin and then back into the transducer. A relatively-simple analysis shows how the back-scattered pulse is spectrally (i.e., affecting the shape and height of the pulse) and temporally (affecting the time-delay with respect to the incident pulse) modified, this being the information that is then exploited (a problem that is different from the forward problem we have just described) for diagnostic purposes.

The reasons why this crude paradigm is so effective are that: 1) the characteristic frequency of the incident pulse is such that the associated wavelength of the acoustic wavefield within the body is small relative to the size of the organ-like target, 2) most of the incident energy is concentrated (focused) into a beam directed on purpose towards the target and not elsewhere, and 3) the density and compressibility of the target organ are sufficiently-different from those of the host medium (within the body) for the reflected wave to have sufficient amplitude so as to be distinguishable from the background returns [71, 49] due to scattering by smaller-scale heterogeneities in the path of the incident and reflected waves.

A less-crude paradigm for the aforementioned biomedical forward problem [49] might consist in assuming that: a) the organ is not layer-like but rather potato-like [71], but still composed of an isotropic, medium, b) this medium is piecewise-homogeneous [49] and fluid-like, c) the host medium is a homogeneous, isotropic fluid, and d) the solicitation is still a plane wave, propagating in the host medium. The forward problem of predicting the acoustic wavefield outside and within this obstacle turns out to be much more complex than the one involving the layer. It is so complex that in recent years the tendency has been to not even try to solve it mathematically (such as based on low-frequency, high-frequency, canonical shape approximations) but rather numerically (e.g., by finite difference [3, 82], finite element [20, 76, 79], spectral element[54], discrete-source [21] methods, etc., which is not to say that these methods do not involve a lot of mathematics too), justified by the fact that powerful computing facilities and software are increasingly available to researchers. The results of such computations are aesthetically and scientifically-speaking quite impressive in that they reveal the striking complexity and variability of the scattered field with respect to factors such as frequency, composition and shape [81, 20, 79, 51, 61, 37, 48, 11, 50, 80, 14]. Making use of this information for diagnostic purposes such as in biomedicine is another, even more challenging, problem.

The diagnostic problem is more familiar to the layman than the forward problem because it occupies him permanently. For instance, a sound received within his auditory system is transduced, encoded and conveyed to the brain wherein it is analyzed [16] so as to provoke a reaction such as withdrawal (reaction to danger as manifested by the loud and/or nearby sound source), pleasure (reaction to a melodious sound) or other types of reaction. The diagnostic problem is therefore one of exploiting the information conveyed by a wave (e.g., sound) to analyze the source (e.g, its distance from the ear) of this wave and/or the medium/structure (e.g., a wall) which the wave encounters before being detected, the purpose of the analysis being to furnish the means for making a decision (e.g., in the biomedical sphere, this might be a decision as to whether an organ submitted to ultrasound is healthy or diseased [95, 36]).

The diagnostic problem, as it relates to the exploitation of the information contained in scattered waves, is known in the scientific community as the inverse scattering problem (because it is the other way around of the forward scattering problem, i.e., determine the source and/or host medium and/or the obstacle from the scattered wavefield). How a human being copes with such a problem in the aforementioned auditory context is nothing less than astounding and certainly enabled by the incredible power and efficiency of the brain. But, a more detailed look at how this functions would probably reveal that the judgements and reactions of a human being to a sound are often either sub-optimal (or even disastrous) and/or simply biased [9], mostly due to the ill-posed nature of inverse problems [77, 26]. Such errors and bias, which are a consequence of either the non-existence, and/or the non-uniqueness and/or the instability of the solution (i.e., the fundamental properties of ill-posed problems) are generally not admissible in the biomedical diagnostic context, but the question is: can they be avoided, or, in other terms: if one assumes so-and-so (e.g., certain aspects of the model concerning how the acoustic wave interacts with the environment), what may the consequences be for his diagnosis [33]?

Our initial intention [94] was to seek answers to this question in relation to the inverse scattering problem of the retrival of one or two constitutive properties of the medium of which an obstacle of arbitrary shape (i.e., potato-like) is composed after having been solicited by a wave. If the related forward problem is difficult to solve (and, as written earlier: now usually carried out numerically), the inverse problem is even harder to deal with for several reasons: 1) it is mathematically non-linear even when the associated forward problem is linear, this being one of the causes of ill-posedness, 2) usually, it cannot be solved in a mathematically-sound (e.g., algebraic manner, as for searching for the roots of a polynomial equation) [88, 83], which means that it is treated algorithmically by what resembles a trial and error optimization technique requiring many resolutions of the associated forward-scattering problem [62, 86, 46], 3) it is not clear what aspects, and what amount, of the scattered-wave data are necessary to treat the inverse problem in the best manner, 4) in real-world situations (such as in geophysical applications [82]), one may dispose of either a very small amount of data (which may be somewhat inappropriate) or a large amount of disparate, unsynchronized data gathered by measurement methods that are not-easily controlled, 5) many of the parameters of the solicitation, the host medium and even of the obstacle (aside, from those that are searched-for) are either not at all, or poorly, known [46, 92], and 6) these parameters could also be searched-for by the retrieval scheme, but the latter becomes more difficult as the number of to-be-retrieved parameters increases [67].

The aforementioned second point is crucial in that it forces one to think twice before undertaking a purely-numerical resolution of the inverse problem. This means that the traditional, more, mathematically-oriented methods might have their say, notably in the way each of the potentially-multiple forward problems are formulated. The latter are so-called trial problems that differ from each other by the values assigned to the to-be-retrieved (in our case, constitutive) parameters. The way a trial problem is treated is called the trial model. What one usually strives for is a trial model that is as mathematically-simple as possible, while being hopefully-able to account for the most important features of the scattered wavefield. Such simple models are often the outcome of certain assumptions and approximations in what is called the domain integral formulation (treated in depth in [94]) of the forward-scattering problem. These assumptions usually have to do with:

1. treating an elastic wave problem (in a solid or porous medium) as an acoustic wave problem [13, 27, 41, 34, 42, 8] (in a so-called equivalent fluid),

2. treating a microscopically-inhomogeneous (e.g., porous) medium as a macroscopically-homogeneous (effective) medium [7, 18, 17, 19, 5, 58, 93],

3. treating the bioacoustic, marine acoustic, electromagnetic and geophysical problem as one in which the mass density, or another constitutive property, is constant everywhere (i.e., is the same and spatially-constant within the obstacle as well as in the host)[65, 30, 64, 73, 22, 23, 45, 74, 66, 27, 84, 55, 52, 15, 34, 96, 7, 29, 57, 5, 78, 36, 37, 86, 96],

4. treating a 3D problem as a 2.5D, 2D or even 1D problem [72, 24, 94] (and many other references)

5. treating the potato as a sphere, coated sphere, cylinder, coated cylinder, layer, multilayer ([38, 88, 37, 90, 4], or other body of canonical shape,

(6) treating the host and/or obstacle as being linear and isotropic [94] (and many other references),

7. invoking ad hoc expressions of the scattered field such as the so-called Born, Kirchoff, physical optics, tangent plane approximations [47, 35, 85, 13, 66, 84, 55, 75, 86, 15, 34, 96, 28, 82], etc.

Sometimes, questions have been raised and addressed as to whether the acoustic ’approximation’, the 2D ’approximation’ etc. are really justified (i.e., lead to accurate parameter retrievals), but when this question relates to the retrieval of the density it is understandable that it has not been considered to be interesting since it is not possible to retrieve the density from the scattered field with a trial model based on the constant-density assumption. Less comprehensible is the fact that the question of the influence of the constant-density assumption on the quality of the retrieval of another parameter such as the wavespeed in, or compressibility of, the medium in the obstacle has not received much attention [53, 44, 94] either, except insofar as porous media are concerned [31]. Nevertheless, it should be underlined that in the area of electromagnetic scattering, the (acoustic) constant-density assumption becomes the oft-employed assumption of constant permeability ([72, 60], p. 2), which is now routinely dropped, notably to account for the exotic properties of metamaterials [10, 69, 70, 25, 91]. Also, more and more studies of acoustic (and elastic wave) inverse scattering no longer appeal to the constant-density assumption [64, 12, 33, 52, 68, 17, 43, 53, 90, 92, 44, 59, 94]

A less-radical (than the constant-density) assumption is that: a) the density of the host is spatially-constant, b) the density of the obstacle, different from that of the host, is also spatially-constant, and c) the mass-density contrast $\epsilon$ (involving the difference of the two constant densities) is small. If the inverse problem is to retrieve the obstacle mass-density (which appears, for instance, to be of interest in certain biomedical applications [96, 43, 95, 53, 90]), an important question would appear to be: how small must $\epsilon$ (which is unknown, like the obstacle mass density, when it is the sought-for-parameter) be for the retrieval of $\epsilon$ (or of another constitutive parameter such as the wavespeed in the obstacle) to be reliable if at the outset the small-$\epsilon$ assumption is incorporated in the trial model? This is the principal question that will shall address in our investigation.

2 The forward-scattering problem

2.1 Preliminaries

The material in this section concerns the determination of the scattered pressure field due to a radiated (by a source) pressure field, initially propagating in a lossless fluid, which encounters a fluid-like obstacle placed somewhere outside of the support of the source. The characteristics of the source (or of the field it radiates), as well as the location, geometrical and physical parameters of the obstacle, are assumed to be known, which means that the only unknowns are the pressure fields outside, and within, the obstacle.

We shall assume that: 1) the source density is invariant with respect to $z$, 2) the obstacle is of infinite extension in the $z$ direction, 3) the host material surrounding the obstace is a lossless, homogeneous, isotropic fluid, 4) the fluid-like material within the obstacle is isotropic and possibly lossy, and 5) the eventual heterogeneity of the obstacle is such that the mass density and wavespeed therein may depend only on $x$ and $y$ but not on $z$. The principal consequence of these assumptions is that the incident (i.e., radiated by the source) and scattered (i.e., due to the presence of the obstacle) pressure fields do not depend on $z$, i.e., the problem is 2D and can henceforth be examined in the $x-y$ (i.e., sagittal) plane.

The time domain governing partial differential equation (PDE) for the total pressure field $p$ (in the entire sagittal plane, i.e., above, within and below the layer obstacle) is the 2D Bergman equation

[TABLE]

wherein $\mathbf{x}=(x,y)$, $t$ is the time variable, $\rho$ the mass density, $\mathbf{f}$ the applied force associated with the source, and

[TABLE]

with $c(\mathbf{x})$, $K(\mathbf{x})$ and $\kappa(\mathbf{x})$ the phase velocity (also termed wavespeed), isentropic bulk modulus and isentropic compressibility within/of the fluid medium.

By expanding the time domain pressure and force in the Fourier integrals

[TABLE]

[TABLE]

wherein $\omega=2\pi f$ is the angular frequency and $f$ the frequency, one obtains the frequency domain expression of the Bergman wave equation

[TABLE]

wherein

[TABLE]

is the wavenumber, which is a generally-complex quantity due to the fact that $c$ is assumed to be generally-complex (to account for losses in the obstacle portion of the scattering configuration). Eq. (5) is the frequency domain version of the 2D Bergman PDE (BPDE).

The BPDE implies certain properties of the pressure field:

(i) $p(\mathbf{x},\omega)$ is continuous and bounded everywhere in $\mathbb{R}^{2}$ except within the support of the source,

(ii) $\frac{\nabla_{\mathbf{x}}\rho(\mathbf{x},\omega)}{\rho(\mathbf{x},\omega)}$ is continuous and bounded everywhere in $\mathbb{R}^{2}$ except within the support of the source.

Since the BPDE applies to all of $\mathbb{R}^{2}$, the outer reaches of this infinite domain constitute a sort of boundary on which it turns out to be necessary to specify the behavior of the pressure field. Seen by an observer situated on this ’boundary’, the field is considered (on grounds of physical plausibility) to behave as a wave emerging from (rather than converging to) the sources (both applied due to $\mathbf{f}$, and induced within the obstacle due to its solicitation by the radiated wave), this being what is called the radiation (at infinity) or outgoing wave condition.

The BPDE is the differential form of the frequency domain governing equation of the scattering problem. As is well-known, it is also possible to cast this problem in integral form. The result is the so-called Bergman domain integral equation (BDIE)

[TABLE]

wherein $d\Omega$ is the differential area in the sagittal plane and $G^{[0]}(\mathbf{x};\mathbf{x^{\prime}};\omega)$ the free-space Green’s function satisfying

[TABLE]

and subject (together with $p$) to the radiation condition.

Henceforth, we shall be concerned with the resolution of the BPDE and the BDIE for a specific (i.e., layer-like) obstacle and solicitation (i.e., plane wave), keeping in mind properties (i), (ii) and the radiation condition.

2.2 Description of the problem of scattering of a plane wave by a layer

The obstacle is a fluid-like, flat-faced layer of infinite extent in the $x$ and $z$ directions. The upper and lower faces of the layer are the planes $y=0$ and $y=-h$ respectively, with $h$ the thickness of the layer. The sagittal ($x-y$) plane view of the scattering configuration is given in fig. 1, wherein $\mathbf{k^{i}}$ is the wavevector of the incident (this word is employed to distinguish it from ’radiated’, but the incident and radiated waves play the same role which is to solicit the obstacle) plane wave $p^{i}(\mathbf{x,\omega})=a^{i}(\omega)\exp(i\mathbf{k}^{i}(\omega)\cdot\mathbf{x})$ whose amplitude is $a^{i}(\omega)$.

$\mathbb{R}^{2}$ is initially (i.e., before the introduction of the layer) filled with a homogeneous isotropic fluid medium $M^{[0]}(\rho^{[0]},c^{[0]})$, in which $\rho^{[0]}$ and $c^{[0]}$ are position- and frequency-independent.

$M^{[0]}(\rho^{[0]},c^{[0]})$ is termed the host (i.e., to the obstacle) medium. After the introduction of the layer-like obstacle, $M^{[0]}(\rho^{[0]},c^{[0]})$ continues to occupy the spaces above ($\Omega_{0}$) and below ($\Omega_{2}$) the layer, whereas the space ($\Omega_{1}$) within the layer is occupied by the fluid-like medium $M^{[1]}(\rho^{[1]},c^{[1]})$, the latter being isotropic, possibly-heterogeneous (but only with respect to $y$) and generally-lossy, this meaning that $c^{[1]}$ (but also possibly $\rho^{[1]}$ such as in fluid-saturated porous materials) is(are) complex. Thus, $\rho(\mathbf{x},\omega)$ and $c(\mathbf{x},\omega)$ in the BPDE take on the values of $M^{[0]}$, $M^{[1]}$, and $M^{[0]}$ as one progresses from $y=-\infty$ to $y=\infty$ within $\mathbb{R}^{2}$. It may (although this is not mandatory) be convenient to assume that the total pressure field $p(\mathbf{x},\omega)$ takes on the values $p^{[2]}(\mathbf{x},\omega)$ (within $\Omega_{2}$), $p^{[1]}(\mathbf{x},\omega)$ (within $\Omega_{1}$), and $p^{[0]}\mathbf{x},\omega)$ (within $\Omega_{0}$) as one progresses from $y=-\infty$ to $y=\infty$ within $\mathbb{R}^{2}$. Finally, it is convenient, but again not mandatory, to assume that the scattered pressure field is generally equal to the total pressure field minus the incident pressure field throughout $\mathbb{R}^{2}$. Actually, we will change this rule in the so-called DD-SOV method of solution outlined hereafter.

2.3 The domain decomposition, separation-of-variables (DD-SOV) solution via the BPDE

At the heart of the DD-SOV method is the very strong assumption

[TABLE]

which means that $\rho^{[1]}$ and $c^{[1]}$ (and therefore $k^{[1]}=\omega/c^{[1]}$) are constants with respect to $\mathbf{x}$. The DD-SOV method accounts explicitly for the fact that $\mathbb{R}^{2}$ is decomposed into three subdomains, $\Omega_{0},~{}\Omega_{1},~{}\Omega_{2}$, each in which both constitutive functions $\rho(\mathbf{x},\omega)$ and $c(\mathbf{x},\omega)$ (which constitute the coefficients of the BPDE) are constants (with respect to position) so that we are faced with the problem of solving three Helmholtz equations

[TABLE]

wherein $k^{[2]}=k^{[0]}$ and $k^{[j]}=\omega/c^{[j]}~{};~{}j=0,1$ are constants with respect to $\mathbf{x}$.

For plane-wave incidence, $\mathbf{f}(\mathbf{x},t)=0$, and $p(\mathbf{x},\omega)=p^{i}(\mathbf{x},\omega)+p^{d}(\mathbf{x},\omega)$, wherein $p^{i}(\mathbf{x},\omega)$ is the incident plane-wave pressure field

[TABLE]

wherein $\mathbf{k}^{i}=(k_{x}^{i[0]},-k_{y}^{i[0]})=(k^{[0]}\sin\theta^{i},-k^{[0]}\cos\theta^{i})$, $\mathbf{x}=(x,y)$ and $\theta^{i}$ is the angle of incidence. Moreover, $p^{d}(\mathbf{x},\omega)$ is the diffracted (or scattered) pressure field which obeys the radiation condition

[TABLE]

$p^{[j]}(\mathbf{x},\omega)~{};~{}j=0,1,2$ are the total pressure fields in $\Omega_{j}~{};~{}j=0,1,2$ which (due to the assumption that $\rho^{[2]}=\rho^{[0]}$) are connected via the transmission (boundary) conditions arising from properties (i) and (ii) of sect. 2.1:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Applying SOV to the Helmholtz equations gives rise to the field representations

[TABLE]

in which $a^{[0]}(\omega)=a^{i}(\omega)$,

[TABLE]

whereas the invocation of the radiation (i.e., outgoing-wave) condition (12) entails

[TABLE]

so that the problem reduces to the determination of the four unknown amplitudes $b^{[0]},a^{[1]},b^{[1]},a^{[2]}$ by means of the four transmission conditions. This gives rise to the four relations:

[TABLE]

from which it is straightforward to extract the solutions:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

wherein $\gamma^{[j]}=k_{y}^{i[j]}/\rho^{[j]}$.

It is not evident how the mass density contrast $\epsilon=\left(\rho^{[1]}-\rho^{[0]}\right)/\rho^{[1]}$ intervenes in these expressions. For this reason, we shall search for the solution of the scattering problem by another method.

2.4 The domain integral (DI) solution

The point of departure is now the Bergman integral equation (7) wherein $G^{[0]}$ is the 2D Helmholtz operator free-space Green’s function for medium $M^{[0]}$ given by [56]

[TABLE]

wherein $H_{0}^{(1)}$ is the zeroth-order Hankel function of the first kind [1]. Note that we have dropped the $\omega$-dependence in this relation and shall do so from now on for all functions that depend on this variable. This Green’s function admits the following representation in cartesian coordinates (([56], p. 823), with $\mathbf{x}=(x,y)$ and $\mathbf{x^{\prime}}=(x^{\prime},y^{\prime})$:

[TABLE]

in which

[TABLE]

Note that when the mass density is everywhere the same, the Bergman integral equation reduces to what is usually termed the Lippmann-Schwinger (LS) equation

[TABLE]

which, of course, is easier to deal with than the Bergman integral equation. We shall return to the LS integral equation further on, and continue, for the moment, to concentrate our attention on the Bergman integral equation.

2.4.1 The assumption of invariance with respect to $x$

Recall that $M^{[0]}$ was assumed to be homogeneous and isotropic. Although the BPDE (7) and (25)-(27) hold for an arbitrary isotropic, inhomogeneous, although invariant with respect to $x_{3}=z$, medium occupying the ’obstacle’ (i.e., the layer), we shall henceforth assume $M^{[1]}$ to (also) be invariant with respect to $x$, i.e.,

[TABLE]

The consequence of this in (7) is, on account of the plane wave nature of the incident pressure field:

[TABLE]

which, after projection of (7) and employment of the identity

[TABLE]

(with $\delta(~{})$ the Dirac delta distribution), leads to the 1D integral equation

[TABLE]

2.4.2 The assumption of piecewise-constant density with respect to y

We now make the even more radical (but less radical than in the DD-SOV method) assumption that the mass density (but not the wavespeed) is piecewise-constant with respect to $y$. This translates to the relation

[TABLE]

in which $\rho^{[0]}$ and $\rho^{[1]}$ are constants with respect to the spatial coordinates and $H(~{})$ is the Heaviside distribution. It follows that

[TABLE]

so that (32) becomes

[TABLE]

wherein

[TABLE]

[TABLE]

Using the fact that since we assumed at the outset that $c(\mathbf{x})=c^{[0]}$ in $\Omega_{0}$ and $\Omega_{2}$ it follows that

[TABLE]

Now, let us return to $K(y)$. Using the sifting property of the delta distributions we obtain

[TABLE]

We cannot deal with $I$ and $K$ in analytic manner beyond this point unless further assumptions are made.

2.4.3 The assumption of constant wavespeed within the layer

We now make the further assumption (in agreement with what is assumed in the DD-SOV approach) that

[TABLE]

the principal consequence of which is that it seems legitimate to suppose (this is actually a consequence of SOV applied to the Helmholtz PDE in $\Omega_{1}$) that the pressure field within the layer can be represented by the sum of two plane waves:

[TABLE]

($a^{[1]}$ and $b^{[1]}$ are constants) from which it follows (in addition to the continuity properties of $\frac{1}{\rho}p_{,y}$) that

[TABLE]

This suggests that we focus our attention on the integral equation

[TABLE]

It is easy to find via (41):

[TABLE]

or, on account of (41) and the fact that $\left(k_{y}^{i[1]}\right)^{2}-\left(k_{y}^{i[0]}\right)^{2}=\left(k^{[1]}\right)^{2}-\left(k^{[0]}\right)^{2}$:

[TABLE]

It is also readily found via (42) that:

[TABLE]

It follows that (43) takes the form:

[TABLE]

from which, after cancelations of like terms, and the fact that the term multiplying \exp\big{[}-ik_{y}^{i[0]}\big{]} must be equal to $a^{[0]}$ and the terms multiplying \exp\big{[}ik_{y}^{i[0]}\big{]} must vanish if the relation (47) is to hold for every $y$ in $]-h,0[$, that we are left with the two results:

[TABLE]

[TABLE]

These two equations can be put in the form:

[TABLE]

[TABLE]

wherein

[TABLE]

is the mass density contrast. The solution of the two linear equations (50)-(51) in the two unknowns $a^{[1]}$, $b^{[1]}$ is then readily found to be:

[TABLE]

[TABLE]

Since $1-\epsilon=\rho^{[0]}/\rho^{[1]}$ it is easy to see that (53)-(54) are identical to their DD-SOV counterparts (22)-(23. Moreover, (53)-(54 provide the constant-density (i.e., $\epsilon=0$) solution of the Lippmann-Schwinger integral equation

[TABLE]

We now turn to the expressions of the pressure fields in $\Omega_{0}$ and $\Omega_{2}$ which are necessarily an outcome of the solution obtained for the field within the layer, as is expressed by the Bergman (and also LS) integral representations of these fields:

[TABLE]

wherein $I(y)$ and $K(y)$ (=0 for $\epsilon=0$, as for the LS integral equation) are functions of the previously-found $p(y)~{};~{}y\in]-h,0[$.

Let us first consider the field in $\Omega_{0}$. Proceeding as previously, we find

[TABLE]

[TABLE]

so that

[TABLE]

which, upon the introduction of (53)-(54), yields

[TABLE]

wherein

[TABLE]

Since $1-\epsilon=\rho^{[0]}/\rho^{[1]}$ it is easy to see that (61) is identical to its DD-SOV counterpart (21).

Let us next consider the field in $\Omega_{2}$. Proceeding as previously, we find

[TABLE]

[TABLE]

so that

[TABLE]

which, upon the introduction of (53)-(54), yields

[TABLE]

wherein

[TABLE]

Since $1-\epsilon=\rho^{[0]}/\rho^{[1]}$ it is easy to see that (66) is identical to its DD-SOV counterpart (24). Moreover, (60)-(61 on the one hand, and (65)-(66) on the other hand, provide the constant-density (i.e., $\epsilon=0$) solution of the Lippmann-Schwinger integral equations

[TABLE]

respectively.

The important point to note here is that the DI expressions for $b^{[0]},a^{[1]},b^{[1]},a^{[2]}$ enable, contrary to their DD-SOV counterparts, to account, in very explicit manner, for the mass density contrast (expressed by $\epsilon$) of the layer obstacle configuration. We shall see further on that this has important consequences, especially in the inverse scattering context.

2.4.4 The expansion of the transmission coefficient in a series of powers of $\epsilon$

What will be offered here for the transmission coefficient $a^{[2]}$ is applicable, as well to the other coefficients $b^{[0]},a^{[1]},b^{[1]}$.

Eq. (66) is of the form

[TABLE]

in which:

[TABLE]

so that

[TABLE]

This relation suggests that the transmission coefficient can be expressed in iterative manner as follows:

[TABLE]

It is readily-verified that the zeroth-order approximation $a^{[2](0)}$ of $a^{[2]}$ is the constant-density (i.e., the $\epsilon=0$) solution for this coefficient. Similarly, the first-order approximation $a^{[2](1)}$ of $a^{[2]}$ constitutes the first correction, presumably-valid for very small $\epsilon$, of the constant-density solution. This approximation is further corrected via the second-order approximation $a^{[2](2)}$, presumably-valid for small $\epsilon$, and so on.

The questions, to which numerical answers will be given in the next section, are (a) how reliable is the zeroth-order (constant-density) approximation (which does not depend on $\epsilon$)?, and (b) to what order of approximation must we go to get a decent approximation of $a^{[2]}$ for a given $\epsilon$ which is not necessarily small?

2.5 Numerical results for the comparison of the exact $a^{[2]}$ with the zeroth, first and second-order (in terms of $\epsilon$) approximations thereof

The aim of the computations was to compare numerically the exact acoustic transmitted amplitude $a^{[2]}$ to: (a) the zeroth-order $\epsilon$ prediction $a^{[2](0)}$ (equal to the constant-density) prediction of $a^{[2]}$, (b) the first-order $\epsilon$ prediction $a^{[2](1)}$ prediction of $a^{[2]}$, and (c) the second-order $\epsilon$ prediction $a^{[2](2)}$ prediction of $a^{[2]}$, as a function of the various parameters $\epsilon$, $f$, $\theta^{i}$, $\Re\left(c^{[1]}\right)$, $\Im\left(c^{[1]}\right)$, the other parameters being fixed at the following values: $\rho^{[0]}=1000~{}Kgm^{-3}$, $c^{[0]}=1500~{}ms^{-1}$, $h=0.2~{}m$, $a^{[0]}=1$.

2.5.1 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the density contrast $\epsilon$

What changes between the four figures figs. 2-5 (for the low frequency $f=2~{}KHz$) and the four figures 6-9 (for the higher frequency $f=20~{}KHz$) is the imaginary part of the layer wavespeed and it is observed in these figures that these changes have very little influence on the relative qualities of the three approximations of the transmitted amplitude. In fact, the constant mass density approximation of $a^{[2]}$ is manifestly both functionally and numerically incorrect, the linear approximation is numerically correct only within a very small interval around $\epsilon=0$, but the quadratic approximation is both functionally and numerically correct for $\epsilon$ in the interval $\epsilon\approx[-0.2,0.2]$ to which corresponds the rather large range $\rho^{[1]}\approx[1250~{}ms^{-1},1875~{}ms^{-1}]$ (when, as is here the case, $\rho^{[0]}=1500~{}ms^{-1})$ that includes the obstacle wavespeeds encountered in most bioacoustic applications, except those regarding teeth and bones (which, anyway, are more solid-like than fluid-like).

What changes between the three figures figs. 10-12 (for the frequency $f=20~{}KHz$) is the real part of the layer wavespeed and it is seen in these figures that the said changes have very little influence on the relative qualities of the three approximations of the transmitted amplitude, so that the same comment as previously applies to these three figures, as to fig. 13 and 14 relative to a frequency that changes from $f=20~{}KHz$ to $f=2~{}KHz$.

Finally, what changes between the four figures figs. 14-12 is the incident angle, which too is seen to have little influence on the relative qualities of the three approximations of the transmitted amplitude, so that the same comment as previously regarding the acceptable range of $\epsilon$ applies to these four figures as well.

2.5.2 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the incident angle $\theta^{i}$

What changes between the two figures figs. 18-19 (both of which concern the effect of the variation of the incident angle) is the layer mass density. It is observed in these figures that, reardless of the layer mass density, the constant mass density approximation of $a^{[2]}$ is functionally correct but numerically rather incorrect, whereas the linear approximation, and more so the quadratic approximation, are both functionally and numerically correct for most incident angles, this of course, being true for the chosen priors.

2.5.3 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the frequency $f$

What changes between the two figures figs. 20-21 is the imaginary part of the layer wavespeed which is seen to obliterate the small-scale dispersive effects for the larger $|\Im c^{[1]}|$. The zeroth-order approximation of the transmission amplitude is seen to be functionally-correct for the large-scale dispersive effects only, the first-order approximation to be both functionally and numerically rather correct for both the large and small-scale dispersive effects, and the second-order approximation to be be even more correct for both of the dispersive effects.

What changes between the two figures figs. 21-22 (both for $\rho^{[1}=1300~{}Kgm^{-3}$) is the real part of the layer wavespeed which is seen to not affect the aforementioned relative qualities of the three approximations of the transmitted amplitude.

What changes between the two figures figs. 23-24 (both for $\rho^{[1}=700~{}Kgm^{-3}$) is again the real part of the layer wavespeed which is seen to affect the aforementioned relative qualities of the three approximations of the transmitted amplitude to the extent that now even the second-order approximation has trouble of accurately-accounting for the small-scale dispersive effects, especially in the first of these two figures.

Figs. 25-26 (wherein the changes concern the real part of the layer wavespeed) apply to a lossless layer so that the small-scale dispersive effects are amplified. Now, the real part of the layer wavespeed is smaller than in the previous two figures, which fact improves the quality of the three approximations, especially the second-order one, even to the extent of accounting quite well for the small-scale dispersive effects. Fig. 27 is a zoom of fig. 26 whose purpose is to show that the small-scale dispersive effects of the transmitted amplitude are of pseudo-resonant (not resonant because the response is not infinite for a lossless configuration) nature.

2.5.4 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the real part of the wavespeed $\Re\left(c^{[1]}\right)$ in the layer

What changes between figs. 28 and 29 (for $\rho^{[1]}=1100~{}Kgm^{-3}$) is the frequency; likewise between figs. 30 and 31 (for $\rho^{[1]}=700~{}Kgm^{-3}$). In all four of these figures it is apparent that the zeroth-order appoximation of the transmitted amplitude is functionally quite correct, the first-order even more so and numerically quite correct, and the second-order approximation nearly-coincident with the exact transmitted amplitude over the considered range of $\Re c^{[1]}$.

2.5.5 Zeroth, first, and second-order $\epsilon$ approximations of $a^{[2]}$ compared to the exact expression thereof as a function of the imaginary part of the wavespeed $\Im\left(c^{[1]}\right)$ in the layer

What changes between figs. 32 and 33 (for $\Re c^{[1]}=700~{}ms^{-1}$) is the frequency; likewise between figs. 34 and 35 (for $\Re c^{[1]}=1300~{}ms^{-1}$). In all four of these figures it is apparent that the zeroth-order appoximation of the transmitted amplitude is functionally quite correct, the first-order even more so and numerically quite correct, and the second-order approximation coincident with the exact transmitted amplitude over the considered range of $\Im c^{[1]}$.

3 The inverse-scattering problem

3.1 Preliminaries

The following material concerns the retrieval of either $\epsilon$ or $\Re\left(c^{[1]}\right)$ from data relating to the measured field in $\Omega_{2}$.

Actually, we shall not measure this field, but rather obtain it by simulation, the latter appealing to the exact solution of the forward problem.

Furthermore, we shall assume that the configurations of the forward and inverse problems are exactly the same, this meaning that (a) the obstacle is the same layer (i.e., same: thickness, macroscopic homogeneity invariant with respect to the three space coordinates, location in space, but differing a priori only by the mass density contrast and/or the wavespeed), (b) the obstacle is submitted to the same radiation (i.e., same plane wave), and (c) the measurement location and device are the same (actually a device capable of detecting with perfect accuracy the amplitude and phase of the pressure on the lower face of the layer).

Inversion of this sort normally is the outcome of the comparison of two pressure fields (at the above-mentioned location): 1) the one predicted by a trial model of the layer involving a trial value of the to-be-retrieved parameter (which may initially be far-removed from the actual value of this parameter as it intervenes in the data) and 2) the simulated pressure field. If the difference (the meaning of this word will be detailed further on) of the trial field from the measured (simulated) data field is larger than some prescribed value, then the trial parameter in the trial model is changed in some hopefully-rational manner and the new trial model is again compared to the data. This comparison procedure is repeated as many times as necessary to obtain a difference equal to, or smaller than, a prescribed value. We shall employ this procedure to retrieve $\Re\left(c^{[1]}\right)$. However, inversion is simpler when the difference between the predictions of (the pressure fields) trial model and the simulated data are imposed at the outset to vanish and the equation translating this fact can be solved in explicit, mathematical manner. This will be the procedure [88] for retrieving the density contrast $\epsilon$.

In all the following material, upper-case letters will designate trial parameters (and the corresponding trial fields) and lower-case letters will designate the actual parameters (and data fields) involved in the simulation model employed to generate the data. When a parameter in the trial model is not varied (i.e., not the object of the retrieval), it is termed ’prior’, and due to the aforementioned assumptions, all the priors are considered to be equal to the corresponding parameters in the simulation model. This means, in particular, that when trying to retrieve the layer wavespeed, we assume the mass density contrast prior to be equal to its value in the simulation model, and when trying to retrieve the mass density contrast, we assume the layer wavespeed to be equal to its value in the simulation model.

The simulation model employed to generate the data appeals to the exact DD-SOV solutions. We shall employ three trial models to retrieve the density contrast: these are based on the zeroth, first and second-order $\epsilon$ approximations of the pressure field, and we shall show that these three inversions can be carried out in explicit, mathematical manner. We shall employ the same three trial models to retrieve the real part of the layer wavespeed, but these three inversions cannot be carried out in explicit mathematical manner. In fact, they appeal to the numerical, algorithmic comparison procedure described previously.

This will enable us to show that the employment of approximate trial models, such as the one based on constant-density assumption, lead to impossible or inaccurate retrievals, much as do the employment of exact trial models with one or several trial priors being different from the corresponding parameters employed in the simulation model.

3.2 Retrieval of the mass density contrast

The simulated (via the DD-SOV model) field in the transmission half-space is

[TABLE]

We assume the data is collected at $x=0,~{}y=-h$, so that the simulated pressure data is

[TABLE]

The trial fields, at the same location, for the three trial models, are of the general form

[TABLE]

wherein $E$ is the trial mass density contrast, and due to the assumption of the identity of the priors, $K_{y}^{i[0]}=k_{y}^{i[0]}$, $H=h$, so that

[TABLE]

wherein, it must be recalled:

[TABLE]

with the coefficients $A_{j}^{[2]}$ not depending on $E$.

The so-called cost function is a function of the difference of $P^{(l)}(0,-h)$ from $p(0,-h)$

[TABLE]

the object being to find the trial density contrast $E^{(l)}$ that minimizes $\kappa^{(l)}$ for each $l$. For reasons that will become obvious, we shall choose $\mathcal{F}\left((P^{(l)}(0,-h)-p(0,-h)\right)=P^{(l)}(0,-h)-p(0,-h))$ and the minimum to be zero, so that the inverse problems reduce to finding $E^{(l)}$ from

[TABLE]

or, more explicitly,

[TABLE]

which implies

[TABLE]

First consider the case $l=0$. On account of (76) the previous equation reduces to

[TABLE]

which possesses no solution for the mass density contrast since $A_{0}^{[2]}$ does not depend on this parameter. Thus, the somewhat-trivial conclusion is that it is not possible to retrieve the density contrast when the constant-density assumption is made.

Next consider the case $l=1$. On account of (76) the comparison equation (80) reduces to the linear equation (in terms of the mass density contrast)

[TABLE]

whose single solution is simply

[TABLE]

This shows that the employment of the linearized approximation $A^{[2](1)}$ of $A^{[2]}$ as the trial model enables: (1) not only to solve the inverse problem in explicit, mathematical form, (2) but also to obtain a unique solution of this problem. These are the reasons why linearization is often employed in dealing with inverse problems. A last property of this solution: it may be complex even when $\epsilon$ is real, this being due to the fact that the different terms in (83) are generally-complex.

Finally, consider the $l=2$ case. On account of (76) the comparison equation (80) reduces to the quadratic equation (in terms of the mass density contrast)

[TABLE]

whose two solutions are simply

[TABLE]

This shows that the employment of the second-order approximation $A^{[2](2)}$ of $A^{[2]}$ as the trial model enables: (1) not only to solve the inverse problem in explicit, mathematical form, (2) but also to obtain a non-unique solution of this problem. In fact, we find two solutions to this problem, and there does not appear to exist any easily-discernible relation between these two solutions on the one hand and the sole solution of the the $l=1$ problem on the other hand. A last property of the second-order solution: it may be complex even when $\epsilon$ is real, this being due to the fact that the different terms in (85) are generally-complex.

It has been suggested that to alleviate the non-uniqueness problem, and/or improve the accuracy of the retrievals, one should take advantage of the fact that a retrieved constitutive parameter should not depend on the characteristics of the solicitation because of our assumption of isotropy. We translate this remark here to: $E^{[2]}$ should not depend on the angle of incidence $\theta^{i}$, so that to possibly alleviate the non-uniqueness problem and/or increase the accuracy of the retrieval, we generalize the comparison equation (81) to $N$ realizations which differ from each other only by the choice of incident angle.

The chosen incident angles form the set

[TABLE]

and for the $n$-th realization, the simulated transmission coefficient is $a^{[2]}(\theta^{i}_{n})$ while its trial counterpart is $A^{[2]}(\theta^{i}_{n})$ so that the comparison equation now takes the form

[TABLE]

which is the same as (81) when $N=1$.

Proceeding as previously, we make the expansion (for each incident angle realization)

[TABLE]

whence (87) becomes

[TABLE]

For $l=1$ this equation yields the unique solution

[TABLE]

whereas for $l=2$ it gives rise to the two solutions

[TABLE]

wherein:

[TABLE]

As previously, the $l=1$ and $l=2$ retrievals for the mass density contrast may turn out to be complex even though $\epsilon$ is real, this being due to the fact that the different terms in (90) and (91) are generally-complex.

3.3 Numerical results for the retrieval of the mass density contrast

The aim of the computations was to compare numerically the actual density contrast $\epsilon$ to the retrievals thereof: (a) $E^{(1)})$ obtained via the first-order $\epsilon$ trial model, and (b) $E^{(2)}$ obtained via the second-order $\epsilon$ trial model, as a function of the various parameters $\Re\left(c^{[1]}\right)$, $\Im\left(c^{[1]}\right)$, $\epsilon$, $f$, $h$ and $\theta^{i}$, the other parameters being fixed at the following values: $\rho^{[0]}=1000~{}Kgm^{-3}$, $c^{[0]}=1500~{}ms^{-1}$, $a^{[0]}=1$. As concerns the set of incident angle(s), we shall designate it in condensed form by $\boldsymbol{\Theta}^{i}=(\theta^{i}_{b},\theta^{i}_{e},N)$, which means that $N$ equally-spaced values of $\theta^{i}$ are chosen for the retrieval ranging from $\theta^{i}_{b}$ to $\theta^{i}_{e}$.

3.3.1 First and second-order retrievals as a function of the frequency $f$

What changes between the six figures 36-41 (relative to the narrow low-frequency range $f=[5~{}KHz,10~{}KHz]$) and the four figures 42-45 (relative to the wider, including high-frequencies, range $f=[0~{}KHz,50~{}KHz]$) is $\boldsymbol{\Theta}^{i}$. The increase of the number of incident-angle realizations is observed to have the unexpected effect of decreasing the quality of the first-order mass-density retrievals. It is also observed that the second-order retrievals are divided into two widely-separated branches, one of which is rather close to the exact mass density $\epsilon$, with the $+$ and $-$ solutions suddenly shifting between these two branches as the frequency increases. These sudden shifts have been observed in other inverse problems and can be qualified as ’instability’ which is known to be one of the characteristics (although usually associated with data, rather than trial model, error [77, 26]) of the ’solutions’ of ill-posed mathematical problems. The increase of the number $N$ of realizations seems not to affect the height of the branch closest to $\epsilon$, which suggests that this might be a useful indicator for distinguishing the ’good branch’ from the ’bad’ branch. Also the increase of $N$ seems, for small $N$, to shift to higher frequencies the moment of shift from one branch to the other, but for larger $N$, ceases to produce this possibly-useful effect. Besides this, the ’good branches’ of the second-order retrievals are clearly-closer than the corresponding first-order retrievals to the actual mass-density $\epsilon$, as one would expect, and clearly, the first-order retrievals are way off mark. Other than these remarks , the noticeable feature of figs. 42-45, also discernible to some extent in figs. 36-41, is the dispersive nature of the retrieved mass density contrast even though the actual mass density contrast was assumed to be non-dispersive (i.e., to not depend on $f$). This dispersive nature is seen to particularly affect: 1) the first-order retrieval and 2) the particular second-order retrieval (of the two possible second-order retrievals) which is the farthest from the actual mass density contrast. This suggests that the dispersion of the retrieval is induced by trial model error and is all the larger the larger is the trial model error.

What changes between the two figures 46-47 is $\rho^{[1]}$. Other than the previous remarks which apply as well here as for a different layer mass density (and therefore different $\epsilon$), one gets the impression that the effects of changing $\boldsymbol{\Theta}^{i}$ are not obviously-beneficial. These remarks also hold for figs. 47-48 which differ from each other by the choice (but not the number) of incident angles.

3.3.2 First and second-order retrievals as a function of the actual density contrast $\epsilon$

What changes between the two figures 49 and 50 is again $\boldsymbol{\Theta}^{i}$ which is seen to have a small effect on the first-order, but no visible effect on the second-order retrievals. Otherwise, the patterns are the same as previously: first-order retrieval rather far from the actual $\epsilon$ except for very small $\epsilon$, and double-branch second-order retrievals, the ’bad’ one of which is really way off from $\epsilon$ except for very small $\epsilon$, and the ’good’ one of which is quite close to $\epsilon$, especially in the interval $\epsilon\approx[-0.2,0.2]$, with the $+$ and $-$ solutions suddenly shifting between these two branches as $\epsilon$ increases. Due to this latter behavior there does not appear to exist a reliable way to spot the ’good’ retrieval for a given $\epsilon$. This is, in fact, only possible when one disposes of rather narrow-range a priori information (such as: $E$ should lie in the interval $[E_{min},E_{max}]$, with the difference of the two values being as small as possible) as to where the retrieved parameter should lie. However, if the trial model is rather crude, as in the case of the first-order model, the single retrieved solution might lie outside the range $[E_{min},E_{max}]$ in which case one would be led to the conclusion that the retrieval (with this crude model) is impossible, this being an incorrect conclusion for a wider range $[E_{min},E_{max}]$, so that one should be cautious about the use of a priori information for the retrieval.

3.3.3 First and second-order retrievals as a function of the layer thickness $h$

What changes between fig. 51 and fig. 52 is the range of $h$, the latter figure constituting a sort of zoom to highlight the dispersive (with respect to $h$) behavior of the retrievals.

What changes between fig. 52 and fig. 56 is $\rho^{[1]}$, but this does not change the general pattern of behavior. Again, the dispersive (with respect to $h$) effects are quite apparent, more so for the first-order and second-order ’bad’ branch retrievals than for the second-order ’good’ branch retrievals, thus providing a possible manner of spotting the retrieval branch that is closest to the actual mass density contrast $\epsilon$.

In the retrievals of figs. 51-56 it should be noticed that we have renounced at considering more than one incident angle realization because the figures in the previous subsections showed that taking more than one incident angle into account does not consistently improve the accuracy nor the uniqueness of the retrievals. For this reason, we shall not, in the next section, take more than one incident angle realization into account in our attempts to retrieve the real part of the layer wavespeed.

3.4 Retrieval of the real part of the layer wavespeed

Again, the simulated (via the DD-SOV model) field in the transmission half-space is

[TABLE]

and once more, we assume the data is collected at $x=0,~{}y=-h$, so that the simulated pressure data is

[TABLE]

wherein $c^{[1]^{\prime}}=\Re\left(c^{[1]}\right)$. The trial fields, at the same location, for the three trial models, are, due to the assumption of the identity of the priors, $K_{y}^{i[0]}=k_{y}^{i[0]}$, $H=h$:

[TABLE]

wherein $C^{[1]^{\prime}}=\Re\left(C^{[1]^{\prime}}\right)$ is the trial parameter and

[TABLE]

with the coefficients $A_{j}^{[2]}$ not depending on $E$. Now, the density contrast is a prior, and due to our assumption of the identity of the priors of the data simulation and trial models, $E^{(l)}=\epsilon~{};~{}l=0,1,2$, so that the previous equation becomes

[TABLE]

wherein the coefficients $A_{j}^{[2]}$ depend on the to-be-retrieved trial wavespeed $C^{[1](l)^{\prime}}=\Re\left(C^{[1](l)}\right)$ in a complicated manner. For this reason, it is not obvious how to solve for $C^{[1](l)^{\prime}}$ by the algebraic method previously employed for retrieving the density contrast.

However, the cost functions are still a function of the difference of $P^{(l)}(0,-h)$ from $p(0,-h)$

[TABLE]

the object being to find the trial layer wavespeed $C^{[1](l)^{\prime}}$ that minimizes $\mathcal{K}^{(l)}$ for each $l$. More often than not, in the inverse problem context, $\mathcal{F}$ is chosen to be $\mathcal{F}=\|(P^{(l)}(0,-h)-p(0,-h)\|^{2}$ and the optimal trial solution is the one for which the cost function

[TABLE]

attains its absolute minimum.

Actually, searching for an absolute minimum is not as simple as it seems since the cost function usually has more than one minima and the number of these minima usually increases with the size of the search interval. Shrinking the search interval might be an appealing way to solve this problem, but, in fact, it is dangerous in that it increases the probability of excluding a priori a series of minima, one of which might be the absolute minimum one is looking for. In fact, there is no reasonable solution to this non-uniqueness problem unless one is sure at the outset that the to-be-retrieved parameter must lie between a given lower and upper bound.

A final remark: the proposed method for solving the wavespeed-retrieval problem is numerical and algorithmic in nature, in contrast to the mathematical, explicit, manner for solving the mass density contrast-retrieval problem; furthermore we have limited it to a single incident angle realization due to the inconclusive previously-found results as to the possible amelioration of the retrievals afforded by multiple incident angle realizations.

3.5 Numerical results for the retrieval of the real part of the layer wavespeed

The aim of the computations was to compare numerically the actual real part of the wavespeed $\Re\left(c^{[1]}\right)$ in the layer to the retrievals thereof: (a) $\Re\left(C^{[1](0)}\right)$ obtained via the zeroth-order $\epsilon$ (i.e., constant-density) trial model, (b) $\Re\left(C^{[1](1)}\right)$ obtained via the first-order $\epsilon$ trial model, and (c) $\Re\left(C^{[1](2)}\right)$ obtained via the second-order $\epsilon$ trial model, as a function of the various parameters $\Re\left(c^{[1]}\right)$, $\epsilon$, $f$ and $\theta^{i}$, the other parameters (recall that we are assuming identity of the priors of the trial models to those of the simulation model) being fixed at the following values: $\rho^{[0]}=1000~{}Kgm^{-3}$, $c^{[0]}=1500~{}ms^{-1}$, $h=0.2~{}m$, $a^{[0]}=1$, $\theta^{i}=0^{\circ}$, $\Im\left(c^{[1](2)}\right)=-21~{}ms^{-1}$.

In the following figures of this section the lowest and highest ordinate represent the lower and upper bounds of the search interval within which the real part of the layer wavespeed is sought.

3.5.1 Retrieved real part of the layer wavespeed $C^{[1]^{\prime}}$ as a function of frequency $f$

What changes between fig. 57 and fig. 58 is $\rho^{[1]}$. Throughout the (low) frequency range of the graphs, the zeroth order trial model is observed give inaccurate retrievals whereas the first-order trial model is more accurate and the second-order model coincident with the actual $c^{[1]^{\prime}}$. The closer is $\rho^{[1]}$ to $\rho^{[0]}$ (i.e., the closer is $\epsilon$ to zero), the more-accurate are all three retrievals.

What changes between fig. 58 and fig. 59 is $c^{[1]^{\prime}}$. Owing to the small value of $\epsilon$, this change has little effect on the observed highly-accurate retrievals, particularly those obtained via the first and second order trial models.

What changes between fig. 59 and fig. 60 is $\rho^{[1]}$. The same comments as for figs. 57-58 apply here except that the negative mass density contrast (case of fig. 60) introduces a more-pronounced dispersive behavior in the first and second-order retrievals.

What changes between fig. 60 and fig. 61 is the frequency range which in the second figure now extends from very small to high frequencies. The dispersive behavior of the retrievals, largest for the zeroth and smallest for the second-order trial models means that this dispersion is due to trial model error. Notice should be taken of the fact that all three retrievals tend asymptotically (i.e., for large frequencies) towards the same value $c^{[1]^{\prime}}$ which means that retrieval inaccuracy, especially that resulting from the use of the zeroth-order trial model, is essentially a low-frequency problem.

What changes between fig. 61 and fig. 62 is $\rho^{[1]}$. The latter figure applies to positive $\epsilon$ and the former to negative $\epsilon$ (both corresponding to the same $|\epsilon|$), and the comparison of the two indicates that the negative sign of $\epsilon$ is an aggravating factor for the inaccuracy of the (especially the zeroth-order) retrievals.

Thus, to resume: the dispersive nature of the retrievals of the real part of the layer wavespeed is due to (i.e., is induced by) trial model error rather than by what is often thought to be of physical origin [92], and these dispersive effects appear to be exacerbated by certain choices of the priors such as largely-negative $\epsilon$.

3.5.2 Retrieved real part of the layer wavespeed $C^{[1]^{\prime}}$ as a function of $c^{[1]}$

What changes between fig. 63 and fig. 64 the search interval of $C^{[1]^{\prime}}$ which is being narrowed. In fig. 63 we observe branch-hopping between two retrieval branches of the type formerly observed for the retrieval, via the second-order trial model, of the mass density contrast. But now the hopping occurs for all three trial models employed to obtain the retrievals. One of these branches appears to be ’bad’ and the second is ’good’ insofar as it is linear (as is the actual red line). Narrowing the search interval in fig. 64 results in the disappearance of the higher-lying branch, thus justifying its qualification as being ’bad’. In fig. 64 we also observe near-coincidence of the three retrievals with the actual $c^{[1]^{\prime}}$.

For a value of $\rho^{[1]}$ that is different from its value in the former two figures, what changes between fig. 65 and fig. 67 is again the search interval which is progressively narrowed. At first, it does not enable the elimination of the ’bad’ branch, but finally succeeds in doing this in fig. 67 in which we again observe near-coincidence of the three retrievals with the actual $c^{[1]^{\prime}}$.

For a value of $\rho^{[1]}$ that is different from its value in the preceding three figures, what changes in figs. 68-69 is the search interval which is being narrowed. This narrowing produces the sought-for effect of eliminating the ’bad’ retrieval branch, so that in fig. 69 we again observe near-coincidence of the three retrievals with the actual $c^{[1]^{\prime}}$.

3.5.3 Retrieved real part of the layer wavespeed $C^{[1]^{\prime}}$ as a function of $\epsilon$

What changes between fig. 70 and fig. 71 is the search interval which is being narrowed. In the first of these two figures we observe branch-hopping only for the zeroth-order retrieval, and in the second figure it is seen that the narrowing of the search interval enables the elimination of the ’bad’ branch of the zeroth-order retrieval. In fig. 69 we also observe near-concidence of the retrievals with the actual $c^{[1]^{\prime}}$ over a wide range of $\epsilon$. In this figure we also observe that for very large negative $\epsilon$, the second-order retrieval can even be less-accurate than the zeroth and first-order retrievals.

What changes between figs. 71-73 is the frequency $f$. The reduction of frequency is seen to reduce the interval of $\epsilon$ in which the three retrievals are acceptable. These figures also show that the second-order retrieval is generally more accurate than the other two, albeit only for positive $\epsilon$.

What changes between fig. 73 and fig. 74 is $c^{[1]^{\prime}}$ for the same low frequency $f=500~{}Hz$. Now, in conformity with expectations, the second-order retrieval is more accurate than the first-order retrieval which is more accurate than the zeroth-order retrieval, this being true over the whole range of $\epsilon$. In any case, the zeroth-order retrieval is significantly-inaccurate for $|\epsilon|>0.2$.

4 Conclusions

This investigation dealt with the question of how small must the mass density contrast $\epsilon$ (which is unknown when it is the sought-for-parameter) be for the retrieval of $\epsilon$, or of another constitutive parameter such as the wavespeed in the obstacle, to be reliable if at the outset the small(or zero)-$\epsilon$ assumption is incorporated in the trial model?

As written in the Introduction, the trial model is nothing but a (method of) solution of a forward problem. Since our inquiry had to do with a small-$\epsilon$ assumption, we first searched for a small-$\epsilon$, preferably mathematical, rather than numerical, solution of the forward-scattering problem. This could not be easily-done for an obstacle of arbitrary shape, so that the first idea was to treat the obstacle (as in the crude paradigm described in the Introduction for bioacoustic diagnosis) as being a flat-faced layer submitted to a plane wave. The mathematically-explicit solution to this forward-scattering problem is easily-obtained, but the way $\epsilon$ intervenes in this solution is not easy to discern. For this reason, we chose to search for the solution via the domain integral formulation, whereby we found the same solution as previously, but in a form in which $\epsilon$ clearly emerges. Moreover, this new form of the solution lends itself to a perturbation analysis by which the solution, and thus the sought-for trial model (of similar form to what is found in [88, 83, 2]) is expressed as a series of powers of $\epsilon$, the coefficients of which were found in explicit, algebraic form.

What we called the zeroth-order trial model results from neglecting all except the first terms in the series, and this model (which does not involve $\epsilon$ and therefore is as if $\epsilon$ were equal to zero) is nothing but the one resulting from the constant-density assumption. We showed that the inverse problem of the retrieval of $\epsilon$ is trivially-impossible using this zeroth-order trial model.

Next, we called the first-order trial model the one that results from neglecting all except the first two terms in the series. The employment of this linearized (in terms of $\epsilon$) model enabled, via comparison with simulated data relative to the transmitted field, the retrieval problem to be cast as a linear algebraic equation which was easily-solved for the mass density contrast. The solution was thus found to exist and be unique, these being features of inverse problems that are often strived-for.

Finally, we called the second-order trial model the one that results from neglecting all except the first three terms in the series. The employment of this model enabled, via comparison with simulated data relative to the transmitted field, the retrieval problem to be cast as a quadratic equation which was again easily and explicitly-solved for the mass density contrast. The solution was found to exist, but not to be unique since the quadratic equation possesses two solutions, both of which are equally-plausible.

Since there is no obvious relation between the solutions for the mass density contrast obtained via these three trial models, one can say that either we are not obtaining a solution, obtaining a single solution, or obtaining two solutions that all seem to differ from each other solely on account of the fact of employing different trial models to solve the inverse problem. This result is similar to what happens when one or more of the priors (i.e., the parameters that are not retrieved, but assumed to be known) in the trial model are, in fact, different from the corresponding parameters in the model employed to simulate the scattered field data, or different from what is guessed to be these parameters in actually-measured data [46, 67].

We showed that the use of approximate trial models (i.e., those corresponding to the various orders of $\epsilon$) lead (this was shown numerically) to several remarkable features concerning the retrieved mass density contrast:

1. complex retrievals of what is a real quantity,

2. frequency dispersion of this parameter although the latter is assumed at the outset to be non-dispersive,

3. branch-hopping for the second-order trial model solutions which means, for example, that the solution skips from one branch to another in apparently-haphazard fashion as the frequency (or another parameter) increases,

4. first-order (i.e., based on a very small density contrast assumption) retrievals that are way off mark,

5. one of the second-order solutions is more affected by frequency dispersion, and is farther from the actual mass density contrast, than the other solution, which fact furnishes a possible means of spotting the ’right’ retrieval, provided, of course, the to-be-retrieved parameter is known a priori to actually be non- or feebly-dispersive in the frequency range of interest.

Not unsurprisingly, we found that the ’right’ second-order retrievals were generally-closer than the first-order retrievals to the actual $\epsilon$ for a diversity of priors and frequencies. In fact, the ’right’ second-order retrieval was found to be practically-coincident with the actual $\epsilon$ for $\epsilon\in[-0.2,0.2]$. Finally, we found that increasing the amount of angle of incidence realizations did not have a systematically-beneficial effect on the quality of the retrievals.

The question of the effect of employing different-order $\epsilon$ models to retrieve the real part of the layer wavespeed $c^{[1]^{\prime}}$ could not be treated in mathematical terms and thus was undertaken numerically. Consequently, we considered the retrieval of the single real parameter $c^{[1]^{\prime}}$ to be a non-linear optimization problem which we solved via a standard least-squares algorithm (for the minimization of the corresponding cost function). The retrievals obtained by this technique were found to have the following characteristics:

(i) contrary to the case of the retrieval of $\epsilon$, the zeroth-order model enables the retrieval of the sought-for parameter,

(ii) the zeroth-order, and more so the two other trial models, enable quite-accurate retrievals of $c^{[1]^{\prime}}$ except for very low frequencies and/or large $|\epsilon|$

(iii) non-uniqueness (manifested by branch-hopping) was observed for all three types of retrievals (recall that it was only possible for the second-order retrievals of $\epsilon$) and was eliminated by narrowing the search interval which is justified only if strong a priori information is available as to the actual value of $c^{[1]^{\prime}}$,

(iv) the retrievals of the real part of the layer wavespeed were found to be dispersive (with respect to frequency) even though the actual value of this parameter did not depend on frequency, and this phenomenon was found to affect all three retrievals, although more so for the zeroth-order one, and to be particularly noticeable at low frequencies.

It follows from the preceding remarks that the mass density contrast ($\epsilon$, which equals zero when the constant density assumption is made) is a constitutive parameter that plays an important role in the forward problem solution for the scattered field transmitted by the fluid-like layer solicited by a plane acoustic wave, and consequently this role is as, if not more, important in the inverse problem context in which one strives to retrieve either the mass density contrast $\epsilon$ or the real part of the layer wavespeed $c^{[1]^{\prime}}$ from transmitted field data (the same is probably true for the other constitutive parameters and perhaps for the characteristics of the solicitation). This means that the constant density assumption is unjustified, particularly in the inverse problem context, when $\epsilon$ is poorly-known a priori and may therefore be large enough to cause large inaccuracies in the retrievals of the parameter(s) that is(are) sought-for via the constant density trial model.

Another, more general, conclusion of this study (underlined in our previous investigation [88] for a much simpler scattering configuration), has to do with another aspect of the non-uniqueness of inverse problem ’solutions’: there exist as many of such ’solutions’ as the number of mathematical translations (here related to the number of terms in the $\epsilon$ series) of the trial models one might want to employ.

Choosing the ’most accurate’ trial model does not eliminate this paradox when real, measured, data is compared to the trial model, because it is not certain a priori that this data is as ’exact’ (i.e., devoid of systematic and random errors) as the the trial model and this discrepancy can lead to a difficulty whose origin is described in the next paragraph. Moreover, the adopted trial model necessarily involves a number (which usually is larger the more exact is the trial model) of priors (i.e., the parameters which are not retrieved, but whose values are assumed to be known) and the uncertain knowledge one may have concerning these priors can induce the same type of non-uniqueness of the retrieval [67] as the variety of mathematical expressions for the trial-model. As far as we know, no method, including FWI (full-wave inversion [83]), are immune to these pathological aspects of inverse problem ’solving’.

Finally, a nagging issue in the inverse problem context concerns how to obtain what one can qualify as the ’best’ retrieval. If, for example, the sought-for parameter is the mass density contrast $\epsilon$, then the ’best’ retrieval is obviously attained when the trial parameter $E$ is equal to $\epsilon$, in which case, whatever the definition of the cost function, the latter attains the value zero. In our study, we employed exact simulated data and compared it to approximate trial models with the objective of attaining a zero cost, the result of which generally led to a solution $E\neq\epsilon$ which cannot be qualified as ’best’ even though it is associated with zero cost. The question is then: how to obtain $E=\epsilon$ and zero cost at the same time? The answer is that the simulated pressure field data must be identical to the trial model prediction of the pressure field, this meaning, that we should have simulated the data with the same $\epsilon$ series as the one employed for the trial model (with, of course, the mass-density contrast being denoted by $\epsilon$ in the data series of $l+1$ terms, and by $E$ in the trial model series of $L+1$ terms), which would have led, after imposing the zero- cost condition, to $E=\epsilon$, whatever the number ($l+1=L+1>1$) of terms in both these series. From this one must conclude that the ’best’ retrieval is obtained by tailoring the trial model to the data (this is called committing the ’inverse crime’ [89]), which means that if the data is ’bad’, the trial model should be ’bad’ and if the data is ’good’, the trial model should be ’good’. Of course, tailoring of this sort is only possible when the data is simulated, which means, that when the data is real (as opposed to simulated), it is: 1) impossible to commit the ’inverse crime’ and 2) impossible to obtain the ’best’ retrieval of the sought-for parameter(s).