# Reconstruction of material properties profiles in one-dimensional macroscopically inhomogeneous rigid frame porous media in the frequency domain

Laurent de Ryck, Walter Lauriks, Philippe Leclaire, Jean-Philippe Groby, Armand Wirgin, Claude Depollier

► To cite this version:

Laurent de Ryck, Walter Lauriks, Philippe Leclaire, Jean-Philippe Groby, Armand Wirgin, et al. Reconstruction of material properties profiles in one-dimensional macroscopically inhomogeneous rigid frame porous media in the frequency domain. Journal of the Acoustical Society of America, 2008, 10.1121/1.2959734. [hal-01330348](https://hal.science/hal-01330348)

HAL Id: hal-01330348

https://hal.science/hal-01330348

Submitted on 10 Jun 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Reconstruction of material properties profiles in one-dimensional macroscopically inhomogeneous rigid frame porous media in the frequency domain

L. De Ryck and W. Lauriks

Laboratorium voor Akoestiek en Thermische Fysica, Katholieke Universiteit Leuven, Celestijnenlaan 200D, B-3001 Heverlee, Belgium

P. Leclaire

Laboratoire de Recherche en Mécanique et Acoustique, ISAT/Université de Bourgogne, 49 rue Mademoiselle Bourgeois, B.P.31-58027 Nevers cedex, France

J. P. Groby

Département de Recherche en Electromagnétisme, Laboratoire des Signaux et Systèmes, UMR8506 CNRS/Supelec/UPS11, 3 rue Joliot-Curie, 91192 Gif-sur-Yvette, France

A. Wirgin

Laboratoire de Mécanique et d'Acoustique, UPR 7051 CNRS, 31 chemin Joseph Aiguier, 13009 Marseille, France

C. Depollier

Laboratoire d'Acoustique de l'Université du Maine, UMR 6613 CNRS, Avenue Olivier Messiaen, 72085 cedex 09 Le Mans, France

The present paper deals with the inverse scattering problem involving macroscopically inhomogeneous rigid frame porous media. It consists of the recovery, from acoustic measurements, of the profiles of spatially varying material parameters by means of an optimization approach. The resolution is based on the modeling of acoustic wave propagation in macroscopically inhomogeneous rigid frame porous materials, which was recently derived from the generalized Biot's theory. In practice, the inverse problem is solved by minimizing an objective function defined in the least-square sense by the comparison of the calculated reflection (and transmission) coefficient(s) with the measured or synthetic one(s), affected or not by additive Gaussian noise. From an initial guess, the profiles of the x-dependent material parameters are reconstructed iteratively with the help of a standard conjugate gradient method. The convergence rate of the latter and the accuracy of the reconstructions are improved by the availability of an analytical gradient.

# I. INTRODUCTION

Porous materials are encountered in a large domain of application in building acoustics (concrete walls, rockwools, plastic foams), in aeronautic, in geophysics, in petroleum prospection, in medicine (osteoporosis diagnose)... The resolution of the inverse scattering problem for the characterization of the parameters of rigid frame porous materials is of large importance, particularly when they are spatially varying. The parameters of interest in rigid frame porous materials are the porosity $\phi$, which is the ratio of the fluid volume to the total sample volume; the tortuosity $\tau_\infty$, which is a parameter describing the change in magnitude and direction of the fluid microvelocity due to the curliness of the pores; the characteristic viscous and thermal lengths $\Lambda$ and $\Lambda'$, which are illustrating the geometry of the pores through the

viscous and thermal losses; and the flow resistivity $R_f$, which is the ratio of the fluid viscosity $\eta$ to the fluid permeability $\kappa_f$.

Wave propagation in homogeneous porous materials was mainly studied after Biot's theory$^{1-4}$ and latter contributions.$^{5-8}$ However, the equations of motion for macroscopically inhomogeneous porous materials were only recently derived from the alternative formulation of Biot's theory.$^{3,4,9,10}$ In particular, when the assumption of rigid frame is valid,$^{5-7,11}$ i.e., when the saturating fluid is light, such as air, and the frame is not moving, the equations of motion reduce to those of an equivalent inhomogeneous fluid with spatial and frequency dependent effective density and bulk modulus.$^{9,10,12}$ The frequency band suitable to this approximation is bounded at high frequency by the diffusion limit (when the wavelength is of the order of, or smaller than the pore size), and at low frequency by the Biot characteristic frequency, below the which the skeleton may vibrate. Under these conditions, the reflected and transmitted pressure fields
as well as the internal pressure field can be numerically determined⁹,¹⁰ by means of the wave splitting and transmission Green's functions method.¹³ This method is employed in the present paper to generate synthetic data to feed the optimizer solving the inverse scattering problem. For multilayered materials, the analytic calculation of the scattered fields is preferred.

Several methods have been developed to characterize homogeneous porous samples, or layered homogeneous porous samples. One is based on the measurement of the change in the fluid volume of a chamber without or with the sample (porosimeter for the determination of the porosity)¹⁴; another one concerns the simultaneous measurement of the pressure drop across a sample and the rate of fluid flow through a porous sample (flowmeter for the determination of the flow resistivity)¹⁵; more recent methods deal with ultrasonic measurements in both time and frequency domains to retrieve the following parameters: the tortuosity,¹⁶ the viscous¹⁷ and thermal characteristic lengths (the Qδ method),¹¹ the porosity and the tortuosity simultaneously,¹⁸,¹⁹ and finally the porosity, the tortuosity, and the viscous characteristic length simultaneously.²⁰ The complexity of the problem is hardened by the different sensitivity of each parameter on the physical quantities that are measured and to the frequency of the solicitation. For example, both density and bulk modulus do not depend on Rf in the asymptotic high frequency model.

The aim of the following material is to reconstruct the depth profiles of all or some material parameters from the measurement of the acoustic field scattered (more specifically: reflected and/or transmitted) by a sample of known thickness. To solve such inverse problems, the Wave Splitting method was extensively studied in the time domain.²¹⁻²⁶ More recently, similar methods, together with an optimization approach, were developed in the frequency domain.²⁷,²⁸ The latter are more suitable to model wave propagation in rigid frame macroscopically inhomogeneous porous materials in all its complexity, i.e., over a large frequency range for which the five parameters have to be accounted for. Moreover, this allows us to work with band-limited data, which is more appropriate when one has to deal with measured data. The inverse scattering problem considered here is unidimensional (depth profiling) and is solved iteratively via a minimization approach based on a classical conjugate gradient method (CG). This method minimizes an objective function which compares measured (or synthetic) reflection coefficients with numerically estimated ones obtained by the frequency domain wave splitting and invariant imbedding technique.²⁴,²⁷ The availability of the analytical gradient of the objective function allows to accelerate the rate of convergence of the optimization method. The choice of the invariant imbedding technique instead of the transmission Green's function method is motivated by the use of the reflection coefficient and optionally the transmission coefficient, but also by the necessity to avoid computational traps and errors due to the so-called "back marching effect," also known as the "inverse crime."²⁹ The latter arises when the synthetic data are generated with a direct solver identical or very similar to the one employed in the minimization method. When

the invariant imbedding (inverse problem) and transmission Green's functions technique (direct problem) are used together with the wave splitting method, Gaussian noise is added to the data to avoid the inverse crime.

The last part of this paper deals with profile reconstructions of one or several material parameters. The accuracy and computation time of the reconstructions are discussed.

## II. MACROSCOPICALLY INHOMOGENEOUS RIGID FRAME POROUS MATERIALS

In this section, the acoustic wave propagation model in macroscopically inhomogeneous rigid frame porous materials is presented.⁹,¹⁰,¹² The theory of wave propagation in homogeneous porous materials was initially given by Biot.¹⁻⁴ In most of the plastic foams saturated by a light fluid like air, the rigid frame assumption is valid so that an acoustic excitation impinging on a porous sample induces wave propagation only in the fluid phase. Therefore the viscothermal effects taking place in the pore channels are included in an effective density ρe and an effective bulk modulus Ke of a so-called equivalent fluid.⁷

The model of rigid frame has been recently extended to macroscopically inhomogeneous porous media.⁹,¹⁰,¹² In the frequency domain, the Helmholtz equation in terms of the fluid pressure p inside the equivalent inhomogeneous fluid is⁹,¹⁰

$$ \nabla \cdot \left( \frac{1}{\rho_e(x, \omega)} \nabla p \right) + \frac{\omega^2}{K_e(x, \omega)} p = 0, \quad (1) $$

where ω is the angular frequency and ∇ the nabla operator. Attenuation, viscothermal losses and dispersion are accounted for in the effective density and bulk modulus. Effective sound speed and characteristic impedance are $c_e(x, \omega) = \sqrt{K_e(x, \omega)/\rho_e(x, \omega)}$ and $Z_e = \rho_e(x, \omega)c_e(x, \omega)$.

In all of this paper, inhomogeneities occur along the x direction and so a unidimensional approach is considered. According to the Biot–Johnson–Allard model¹⁻⁷ extended to macroscopically inhomogeneous porous materials:⁹,¹⁰

$$ \rho_e(x, \omega) = \rho_f \frac{\tau_{\infty}(x)}{\phi(x)} \left[ 1 - j \frac{R_f(x) \phi(x)}{\rho_f \tau_{\infty}(x) \omega} F(x, \omega) \right], \quad (2) $$

$$ K_e(x, \omega) = \frac{\gamma P_0 / \phi(x)}{\gamma - (\gamma - 1) \left[ 1 - j \frac{R_f(x) \phi(x)}{\rho_f \tau_{\infty}(x) P_r \omega} G(x, P_r \omega) \right]^{-1}}, \quad (3) $$

wherein γ is the specific heat ratio, ρf the saturating fluid density, P0 the atmospheric pressure and Pr the Prandtl number. The well-defined correction functions of this model are⁵⁻⁷

$$ F(x, \omega) = \sqrt{1 + j4 \frac{\eta \rho_f \tau_{\infty}^2(x)}{R_f(x)^2 \phi(x)^2 \Lambda^2(x)} \omega}, \quad (4) $$

$$ G(x, P_r \omega) = \sqrt{1 + j4 \frac{\eta \rho_f \tau_{\infty}^2(x)}{R_f(x)^2 \phi(x)^2 \Lambda^{\prime 2}(x)} P_r \omega}. \quad (5) $$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Slab of porous material of thickness L. $p^i$ is the incident signal, $p^r$ is the reflected signal, and $p^t$ is the transmitted signal. The angle of incidence is denoted by $\theta$. In a homogeneous medium, $p^-$ are the backward propagating waves and $p^+$ are the forward propagating waves.</figcaption></figure>

The five so-called acoustical parameters $\phi(x)$, $\tau_\infty(x)$, $\Lambda(x)$, $\Lambda'(x)$, and $R_f(x)$ described in the introduction become here space-dependent functions.⁹,¹⁰ Without addition of shape factors, the thermal and viscous permeabilities are assumed to be equal so that the flow resistivity $R_f$ appears in both correction functions.⁶,⁷

The wave splitting and invariant imbedding method chosen to solve both the direct and inverse scattering problem is solely based on the equations of motion (6):⁹,¹⁰

$$ \frac{d}{dx} \begin{pmatrix} p \\ \phi(x)V \end{pmatrix} = \begin{pmatrix} 0 & -j\omega\rho_e(x, \omega) \\ -\frac{j\omega}{K_e(x, \omega)} & 0 \end{pmatrix} \begin{pmatrix} p \\ \phi(x)V \end{pmatrix}, \quad (6) $$

where $V$ is the fluid particle velocity. This set of equations is equivalent to the Helmholtz equation (1).

### III. STATEMENTS OF THE PROBLEM AND DIRECT SCATTERING SOLUTION

In all the following paper, a macroscopically one-dimensional (1D) inhomogeneous porous slab of thickness $L$, solicited by an incident plane wave propagating initially in the air (the saturating fluid here), is considered. This is depicted in Fig. 1.

The wave splitting method was mainly used to solve both direct and inverse scattering problems in several fields of physics, especially in electromagnetism, initially in the time domain<sup>21–23,30</sup> and then in the frequency domain.<sup>9,10,13,27</sup> In order to avoid convolution products in time and fractional derivatives,<sup>18,20</sup> the direct and inverse scattering problems are treated in the frequency domain. Also, the Fourier transform $p(x, \omega)$ of the measurable pressure field $p(x, t)$ is used. The inverse Fourier transform is chosen such that

$$ p(x,t) = \int_{-\infty}^{\infty} p(x, \omega) \exp(j\omega t) d\omega. \quad (7) $$

The inverse scattering problem consists in the retrieval of the profile of the five acoustical parameters from the measurement of the pressure waves reflected by an inhomogeneous porous sample. The solution of the direct scattering

problem is used to improve each profile at each step of the iterative optimization approach. In order to be able to simultaneously reconstruct more than one depth-dependent material property, enough information must be provided. Several frequencies are used in order (I) to regularize the inverse scattering problem, (II) to use a large enough frequency band for most of the acoustic parameters to have an influence on the scattered field, and (III) to be able to recover from the “low frequency” content of the scattered field the average values of the acoustic parameters and from its “high frequency” content their spatial variations, which can be sharp in case of layered medium. On the other hand, the problem is much more complicated than the one usually treated in homogeneous case because, for instance, the porosity and totuosity are space dependent functions and cannot be “simply” recovered from the instantaneous reflection at the first interface of the sample. The additional information available from the wave front during its propagation inside the whole sample has to be caught. To increase the amount of information, several operations are applicable: (I) For the simultaneous reconstruction of two profiles, a two-sided normal reflection measurement can be used, i.e., the recording of both the reflected data at $x=0$ after an excitation from a source at $x<0$ and the reflected data at $x=L$ after an excitation in the opposite direction from $x>L$; (II) For the simultaneous reconstruction of one to four depth profiles, several two-sided reflection measurements at normal incidence can be done by imposing or not a rigid wall condition on the opposite side of the solicitation; (III) Another classical way to solve the problem of reconstructing many depth profiles is to use acoustic waves impinging on the porous sample at various oblique incidence. The last operation is followed and as many profiles as angles of incidence should normally be reconstructed.

The direct scattering problem consists of determining the scattered, i.e., reflected and transmitted, pressure field due to the presence of a macroscopically inhomogeneous porous slab. Therefore, the aim is to solve the system of equations (6) for an incident acoustic excitation from a position $x<0$. The solution is obtained through the wave splitting and invariant imbedding approach<sup>9,10,12,27</sup>, which requires boundary conditions. They are found via the characteristic impedances of both surrounding media for $x<0$ and $x>L$. Here, for all $x<0$, the medium is a homogeneous free fluid (air). Its characteristic impedance is denoted by $Z_0=\rho_f c_f$ with $c_f$ the sound speed in the free fluid. For all $x>L$, the medium can be either air or a rigid wall. Its impedance is denoted by $Z_L$, which is respectively $Z_0$ in the first case and the infinity in the second case.

#### A. Wave splitting transformation

Wave splitting is a method based on the mathematical properties of the wave equation. It had been demonstrated since d'Alembert<sup>31</sup> that the solution of the wave equation in terms of the pressure field $p(x, t)$, which is the inverse Fourier transform of the Helmholtz Eq. (1), can be exactly decomposed into two components $p^+(x, t)$ and $p^-(x, t)$ corresponding to forward and backward propagating waves in a homogeneous medium.<sup>13,21–24,32</sup> The transformation is gener-
alized for inhomogeneous media since it can be seen as an algebraic transformation from the base of functions [p; $\phi(x)V$] to the base of functions [p$^+$; $p^-$]. Several transformations are available. In the present case, the so-called vacuum wave splitting transformation$^{9,10,27,13}$ is chosen,

$$ \begin{pmatrix} p^+ \\ p^- \end{pmatrix} = \frac{1}{2} \begin{pmatrix} 1 & Z_0 \\ 1 & -Z_0 \end{pmatrix} \begin{pmatrix} p \\ \phi(x)V \end{pmatrix}. \quad (8) $$

This choice is possible because the wave splitting is coupled with an invariant imbedding technique. In a general wave splitting transformation, the characteristic material impedance $Z_e(x, \omega)$ should be used.$^{21-23}$ However, the invariant imbedding technique, which will be detailed in the next section, allows a transformation that uses the known characteristic impedance $Z_0=\rho_f c_f$ of the saturating fluid instead.$^{27,28}$ As a schematic explanation, one considers that at the initial state, and for all $x<L$, there is only the ambient fluid of characteristic impedance $Z_0$. Then, layer after layer, the inhomogeneous porous material is constructed from $x=L$ to $x=0$. Therefore, in the present formulation, each nth layer of infinitesimal thickness $dx$ added at $x=L-ndx$ is taken into account as a perturbation of the impedance $Z_0$ at the point $x=L-ndx$. Continuity conditions from layer to layer are implicitly accounted for.

Applying the transformation (8) to the equations of motion (6) in the frequency domain straightforwardly yields$^{9,10}$

$$ \frac{d}{dx} \begin{pmatrix} p^+ \\ p^- \end{pmatrix} (x, \omega) = \begin{pmatrix} -A^+ & -A^- \\ A^- & A^+ \end{pmatrix} \begin{pmatrix} p^+ \\ p^- \end{pmatrix} (x, \omega), \quad (9) $$

with

$$ A^{\pm}(x, \omega) = \frac{j\omega}{2} \left[ \frac{Z_0}{K_e(x, \omega)} \pm \frac{\rho_e(x, \omega)}{Z_0} \right]. \quad (10) $$

For all $x \le 0$, the pressures waves $p^+(x, \omega)$ and $p^-(x, \omega)$ are respectively the incident and reflected pressure waves. Similarly, for all $x \ge L$, $p^+(x, \omega)$ corresponds to the transmitted pressure wave. Therefore, the splitted pressure waves $p^+$ are easily connected to the reflection and transmission coefficients $R(\omega)$ and $T(\omega)$. These coefficients contain information on the inhomogeneities inside the sample, but an extension of these data for each depth $x$ is necessary to take accurately into account the material properties variations between $x=0$ and $x=L$. A continuous approach is done in the following by taking the limit $dx \to 0$.

## B. Invariant imbedding method

In the invariant imbedding method, an imbedding geometry is considered, i.e., a subslab between $x$ and $L$ of the total material slab of thickness $L$. For each subslab, it is assumed that the medium filling the domain between 0 and $x$ is the same as the medium chosen for all $x<0$, e.g., air in the treated case. This allows us to define the $x$-dependent reflection and transmission coefficients $R(x, \omega)$ and $T(x, \omega)$ as follows:

$$ p^{-}(x, \omega) = R(x, \omega)p^{+}(x, \omega), \quad (11) $$

$$ p^{+}(L, \omega) = T(x, \omega)p^{+}(x, \omega). \quad (12) $$

From the above relations, one deduces that $R(0, \omega) = R(\omega)$ and $T(0, \omega) = T(\omega)$ are the usual reflection and transmission coefficients for the total porous slab of length $L$. These coefficients link the measurable reflected and transmitted pressure fields, respectively denoted by $p^r(\omega)=p^-(0, \omega)$ and $p^i(\omega)=p^+(L, \omega)$, to the incident pressure field $p^i(\omega) = p^+(0, \omega)$ illuminating the porous slab at $x=0$, as in Fig. 1.

Combining Eqs. (11) and (12) with (9) and (10) leads to the imbedding equations for $R(x, \omega)$ and $T(x, \omega)$,

$$ \frac{d}{dx}R = 2A^+R + A^-(1 + R^2), \quad (13) $$

$$ \frac{d}{dx}T = (A^+ + A^-R)T. \quad (14) $$

The first equation is a nonlinear Riccati differential equation for $R$, which can be solved with a classical fourth order Runge Kutta algorithm. The second one is a linear differential equation for $T$, which can be solved with the help of the method described in Appendix A, once the Riccati equation is solved for $R$. The resolution of these equations requires boundary conditions to be satisfied. According to the definition of the imbedding geometry, there is no porous subslab between $x$ and $L$ when $x=L$. This means that there is just a simple interface between the two homogeneous media of characteristic impedances $Z_0$ and $Z_L$. These impedances are constants so that the boundary conditions at $x=L$ for the imbedding geometry are

$$ R(L, \omega) = \frac{Z_L - Z_0}{Z_L + Z_0}, \quad (15) $$

$$ T(L, \omega) = \frac{2Z_L}{Z_L + Z_0}. \quad (16) $$

If the homogeneous medium for all $x>L$ is the same free fluid as the one for all $x<0$, then $Z_L=Z_0$ and the boundary conditions reduce to

$$ R(L, \omega) = 0, \quad (17) $$

$$ T(L, \omega) = 1. \quad (18) $$

If the homogeneous medium for all $x>L$ is a rigid wall, then $Z_L=\infty$. There is no transmitted waves, and only the imbedding Eq. (13) for $R(x, \omega)$ is used. The boundary condition reduces to

$$ R(L, \omega) = 1. \quad (19) $$

Equations (11) and (12) together with the boundary conditions (15) and (16) constitute the direct solver.

## C. Notes about the oblique incidence

The oblique incidence implies some slight changes in Eqs. (6)–(14). In the latter, for an angle of incidence $\theta$, the effective bulk modulus $K_e$, the characteristic impedances $Z_0$ and $Z_e$ must be replaced by$^{10}$

$$ \bar{K}_e = K_e \left( 1 - \frac{c_e^2}{c_f^2} \sin^2 \theta \right)^{-1}, \quad (20) $$
$$ \bar{Z}_0 = \frac{Z_0}{\cos \theta}, \qquad (21) $$

$$ \bar{Z}_e = Z_e \left( 1 - \frac{c_e^2}{c_f^2} \sin^2 \theta \right)^{-1/2}. \qquad (22) $$

These expressions are extracted from the projection of the wave vector onto the x axis.

In the next section, the optimization method that is employed to solved the inverse problem of the recovery of the profiles is presented. The transmission data are not often used²⁸ and may even trap the optimizer.²⁷ Therefore, in the following, only the minimization of the reflection data is detailed.

## IV. OPTIMIZATION APPROACH OF THE INVERSE SCATTERING PROBLEM

Let us introduce a five-element vector

$$ \mathbf{p} = (\phi(x), \tau_{\infty}(x), \Lambda(x), \Lambda'(x), R_f(x))^T, \qquad (23) $$

where the superscript T denotes the transpose operation. An objective function $J(\mathbf{p})$ is defined such that²⁷

$$ J(\mathbf{p}) = \sum_{\theta_{\min}}^{\theta_{\max}} \sum_{\omega_{\min}}^{\omega_{\max}} W_R(\omega) |R(0, \omega; \theta) - R_m(\omega; \theta)|^2, \qquad (24) $$

where $R_m$ is the measured reflection coefficient, $R(0, \omega; \theta)$ the reflection coefficient calculated at $x=0$ with the help of our direct solver presented in the previous section, and $W_R$ a non-negative frequency-dependent weighting function.

The objective function is defined such that a summation is performed over the angles of incidence between $\theta_{\min}$ and $\theta_{\max}$, and over a frequency band [$\omega_{\min}, \omega_{\max}$]. The larger the number of angles of incidence and the number of discrete frequencies, the more regularized the minimization problem is.

The optimal values of the elements of **p** are obtained when the objective function $J(\mathbf{p})$ is minimal. These optimal values correspond to the profiles that are to be reconstructed. Therefore the optimization approach consists of minimizing the function $J(\mathbf{p})$ from an initial guess of the elements of the vector **p**. The resolution of our problem, involving the minimization of the function $J(\mathbf{p})$, is hopeless if no other information is available. If the gradient of the objective function is known for all x—analytically or numerically—the convergence to the global minimum is much improved and accelerated.²⁵,²⁷ This is due to the fact that the gradient provides the best directions to follow in order to converge to a minimum. Moreover, if the gradient is available for all x, the profiles of each parameter, along the sample thickness, are optimizable. Then, well-known iterative minimization routines can be applied.

### A. Analytical gradient of the objective function

A small perturbation $\delta\mathbf{p}$ is applied to the material property vector **p**. By definition, the perturbation of the reflection coefficient is $\delta R = \tilde{R} - R$ where $\tilde{R} = R(x, \omega; \mathbf{p} + \delta\mathbf{p})$ is the solution of the perturbed Eq. (13), which becomes:²⁷

$$ \frac{d}{dx} \delta R - 2(A^+ + A^- R) \delta R = 2R \delta A^+ + (1 + R^2) \delta A^- \qquad (25) $$

with the perturbed condition (15) reducing to

$$ \delta R(L, \omega) = 0 \qquad (26) $$

because the impedances $Z_0=\rho_f c_f$ and $Z_L$ are independent of the material properties, which are solely included in the vector **p**.

The perturbations of the coefficients $A^\pm$ are deduced from Eq. (10) and are expressed, for any oblique angle of incidence, as

$$ \delta A^{\pm} = \frac{j\omega}{2} [Z_0 \delta \bar{K}_e^{-1} \pm Z_0^{-1} \delta \rho_e], \qquad (27) $$

wherein the perturbation of $\rho_e$ [Eq. (2)] and of $\bar{K}_e^{-1}$ [Eq. (20)] yield

$$ \delta \rho_e = \frac{\partial \rho_e}{\partial \phi} \delta \phi + \frac{\partial \rho_e}{\partial \tau_{\infty}} \delta \tau_{\infty} + \frac{\partial \rho_e}{\partial \Lambda} \delta \Lambda + \frac{\partial \rho_e}{\partial R_f} \delta R_f, \qquad (28) $$

$$ \delta \bar{K}_e^{-1} = \delta K_e^{-1} + \frac{\sin^2 \theta}{\rho_e^2 c_f^2} \delta \rho_e \qquad (29) $$

with

$$ \delta K_e^{-1} = \frac{\partial K_e^{-1}}{\partial \phi} \delta \phi + \frac{\partial K_e^{-1}}{\partial \tau_{\infty}} \delta \tau_{\infty} + \frac{\partial K_e^{-1}}{\partial \Lambda'} \delta \Lambda' + \frac{\partial K_e^{-1}}{\partial R_f} \delta R_f. \qquad (30) $$

From Eqs. (2) and (3), the partial derivatives in Eqs. (28) and (30) are reduced to functions in terms of $\phi(x)$, $\tau_{\infty}(x)$, $\Lambda(x)$, $\Lambda'(x)$ and $R_f(x)$.

The increment of the objective function (24) induced by the perturbation $\delta\mathbf{p}$ takes the form²⁵,²⁷

$$ \delta J(\mathbf{p}) = J(\mathbf{p} + \delta\mathbf{p}) - J(\mathbf{p}) = 2 \operatorname{Re} \sum_{\theta, \omega} u(0, \omega; \theta) \delta R(0, \omega; \theta), \qquad (31) $$

where $\operatorname{Re}$ represents for the real part and wherein the auxiliary function $u(x, \omega; \theta)$ is defined such that

$$ u(0, \omega; \theta) = W_R(\omega)(R(0, \omega; \theta) - R_m(\omega; \theta))^*, \qquad (32) $$

in which the superscript * denotes the complex conjugate.

Using the boundary conditions (26) at $x=L$, the following integration rule can be applied:

$$ -u(0, \omega; \theta) \delta R(0, \omega; \theta) = \int_{0}^{L} \left[ \frac{d}{dx} u(x, \omega; \theta) \delta R(x, \omega; \theta) \right]_{x=z} dz, \qquad (33) $$

the integrand of which, from Eq. (25), is developed as follows:

$$ \frac{d}{dx}(u \delta R) = \frac{du}{dx} \delta R + u \frac{d \delta R}{dx} = \left( \frac{du}{dx} + 2(A^+ + A^- R)u \right) \delta R + 2Ru \delta A^+ + (1 + R^2)u \delta A^- \qquad (34) $$

$u(x, \omega; \theta)$ is an arbitrary function, except at $x=0$, so that we can choose it in order to eliminate the dependence on the
unknown perturbation $\delta R$ in Eq. (34). Thus, the following linear first order differential equation must be satisfied:

$$ \frac{du}{dx} + 2(A^+ + A^- R)u = 0, \quad (35) $$

with the boundary conditions at $x=0$ given in Eq. (32). This equation can be solved iteratively, as explained in Appendix A.

Combining Eqs. (33) and (31) with the condition imposed in Eq. (35), the perturbation of the objective function reduces on one hand to

$$ \begin{aligned} \delta J(\mathbf{p}) &= -2\mathrm{Re} \sum_{\theta, \omega} \int_0^L \{2Ru \delta A^+ + (1+R^2)u \delta A^-\} dx \\ &= -\mathrm{Re} \int_0^L \sum_{\theta, \omega} j\omega \{Z_0(1+R)^2 u \delta K_e^{-1} \\ &\quad - Z_0^{-1} (1-R)^2 u \delta \rho_e\} dx, \end{aligned} \quad (36) $$

and on the other hand, the same perturbation can be expressed as the inner product

$$ \delta J(\mathbf{p}) = \int_0^L \left[ \frac{\partial J}{\partial \phi} \delta \phi + \frac{\partial J}{\partial \tau_\infty} \delta \tau_\infty + \frac{\partial J}{\partial \Lambda} \delta \Lambda + \frac{\partial J}{\partial \Lambda'} \delta \Lambda' \right. \\ \left. + \frac{\partial J}{\partial R_f} \delta R_f \right] dx. \quad (37) $$

Therefore, with the help of Eqs. (28) and (30), the identification of (36) with (37) leads to the exact expression of the gradient of the objective function. The components of the gradient in terms of each of the parameters to reconstruct are

$$ \frac{\partial J}{\partial \phi}(x) = - \mathrm{Re} \sum_{\theta, \omega} j\omega \left\{ Z_0(1+R)^2 u \frac{\partial K_e^{-1}}{\partial \phi} - Z_0^{-1}(1-R)^2 u \frac{\partial \rho_e}{\partial \phi} \right\}, \quad (38) $$

$$ \frac{\partial J}{\partial \tau_{\infty}}(x) = - \mathrm{Re} \sum_{\theta, \omega} j\omega \left\{ Z_0(1+R)^2 u \frac{\partial K_e^{-1}}{\partial \tau_{\infty}} - Z_0^{-1}(1-R)^2 u \frac{\partial \rho_e}{\partial \tau_{\infty}} \right\}, \quad (39) $$

$$ \frac{\partial J}{\partial \Lambda}(x) = - \mathrm{Re} \sum_{\theta, \omega} -j\omega \left\{ Z_0^{-1}(1-R)^2 u \frac{\partial \rho_e}{\partial \Lambda} \right\}, \quad (40) $$

$$ \frac{\partial J}{\partial \Lambda'}(x) = - \mathrm{Re} \sum_{\theta, \omega} j\omega \left\{ Z_0(1+R)^2 u \frac{\partial K_e^{-1}}{\partial \Lambda'} \right\}, \quad (41) $$

$$ \frac{\partial J}{\partial R_f}(x) = - \mathrm{Re} \sum_{\theta, \omega} j\omega \left\{ Z_0(1+R)^2 u \frac{\partial K_e^{-1}}{\partial R_f} - Z_0^{-1}(1-R)^2 u \frac{\partial \rho_e}{\partial R_f} \right\}. \quad (42) $$

In the above relations, all partial derivatives are functions of the components of the vector

$\mathbf{p}=[\phi(x), \tau_\infty(x), \Lambda(x), \Lambda'(x), R_f(x)]^T$ via the definitions of $\rho_e$ and $K_e$ in Eqs. (2) and (3). The available analytical expression of the gradient will help the iterative minimization scheme to converge faster to the optimal depth profiles of the five parameters of interest. Classically, one of the most suitable iterative algorithms for our problem is the conjugate gradient method.²⁵,²⁷,³³⁻³⁵

A discussion of the numerical calculations, and a presentation of some reconstructions, are given in the following section.

## V. PROFILE RECONSTRUCTION OF x-DEPENDENT POROUS MATERIAL PROPERTIES: SIMULATIONS

In this section, the computation of the direct scattering solver from Sec. III and the iterative method chosen to solve the optimization problem from Sec. IV are described. The differences with similar problems in electromagnetism²⁵ or transmission line theory²⁷ mainly concern the numerical aspects, particularly the choice of the frequency range and the number of points for discretizing our samples. Physical aspects, such as the specific dispersion and absorption processes and properties of porous materials linked to the relations (2)–(5), in which the five material parameters we deal with are introduced, are also treated.

### A. Conjugate gradient method

Due to the fact that the objective function $J(\mathbf{p}(x))$ (24) is not a quadratic form and continuous in $\mathbf{p}(x)$, the optimization problem we deal with is nonlinear. The minimization should give accurate enough results when the analytical gradient of $J$ is available. The classical and widely used conjugate gradient method (CG) is one of the most suitable iterative algorithms²⁵,³³⁻³⁵ for minimization. For clear and rigorous details about the method, the reader can refer to a tutorial (with provided algorithms).³³ For CG algorithms, as well as other optimization techniques, the reader can refer to Polak³⁵ and Press et al.³⁴ Herein, we describe the CG algorithm for the reconstruction of a single parameter $p(x)$, e.g., the porosity $\phi(x)$. The reconstruction of the other parameters can be simply derived in closed form. The gradient $\partial J / \partial p$ is denoted by $G(p)$. The implementation is as follows:²⁵,³³,³⁴

* Step 0: Initial approximation: For all $x$ in $[0, L]$, $p(x)=1 = p^{(0)}$.

* Step 1: Set $i=0$; Solve the direct problem for all $x$ in $[0, L]$ using Eq. (13)–(16); Calculate the gradient of $J$ with respect to $\phi$ with the help of Eq. (38); For all $x$ in $[0, L]$, $G^{(0)}=G(p^{(0)})=\partial J / \partial \phi(p^{(0)})$; Set $D^{(0)}=G^{(0)}$. For any iteration order $i$, $D^{(i)}$ is the search direction.

* Step 2: The line search consists in computing a scalar $\lambda_i$ such that

$$ \forall x \in [0, L], J(p^{(i)} - \lambda_i D^{(i)}) = \min_{\lambda \in \mathbb{R}} \{J(p^{(i)} - \lambda D^{(i)}); \lambda \ge 0\}. \quad (43) $$

* Step 3: Improve the reconstruction by setting $p^{(i+1)}=p^{(i)} - \lambda_i D^{(i)}$.

* Step 4: Compute a new gradient $G^{(i+1)}=G(p^{(i+1)})$ with the help of Eq. (38).
<table><caption>TABLE I. Values of the porous material at $x=0$.</caption><thead><tr><th>Φ</th><th>τ<sub>∞</sub></th><th>Λ<br>(μm)</th><th>Λ'<br>(μm)</th><th>R<sub>f</sub><br>(Ns m<sup>-4</sup>)</th><th>L<br>(cm)</th></tr></thead><tbody><tr><td>0.96</td><td>1.07</td><td>273</td><td>672</td><td>2843</td><td>2</td></tr></tbody></table>

* Step 5: if $G^{(i+1)}<ε$, with $ε \ll 1$ the convergence criterion, STOP. Else, if $i>MaxIter$, with MaxIter the maximal number of iterations, STOP. Otherwise, update the direction search with a Polak-Ribière condition:<sup>33,35</sup>

$$D^{(i+1)} = G^{(i+1)} + \gamma_i D^{(i)} \quad (44)$$

with

$$i = \max\left( \frac{G^{(i+1)T} G^{(i+1)} - G^{(i+1)T} G^{(i)}}{G^{(i)T} G^{(i)}}, 0 \right). \quad (45)$$

Set $i=i+1$ and GO TO Step 2.

The minimization over the parameter $λ$ is easily achieved by means of the Nelder and Mead Simplex method.<sup>34,36</sup> This is a fast way to solve an unconstrained minimization problem with only one or a few variables. In case of the simultaneous reconstruction of more than one parameter, $λ$ becomes a vector with as many components as analytical gradients (or material properties to reconstruct).

## B. Numerical settings

The general procedure for carrying out the reconstruction of the material profiles is now described.

The simulations are performed for an acoustic plane wave excitation coming from $x<0$, for several angles of incidence ranging from $θ=0$ rad to $θ=π/3$ rad, and either for $Z_L=Z_0=ρ_f c_f$ (same fluids on each side of the slab) or $Z_L=∞$ (rigid wall at $x=L$). The advantage of the rigid wall is to detect a reflected signal with large amplitude despite the strongly absorptive material through which the waves travel two times.

## 1. Generating profiles and scaling

For each parameter to reconstruct, a depth profile is created. The value of each parameter at the boundary $x=0$ is given in Table I. These values were measured on a porous polyurethane foam.<sup>11</sup> When performing the numerical reconstruction of some of the five parameters along the depth $x$, the others remain constant.

The parameters in Table I can be classified as: (i) those without dimensions ($Φ$ and $τ_∞$), which vary around one, their value in the ambient fluid, (ii) those with dimensions of length ($Λ$ and $Λ'$), which are of the order of several hundreds of micrometers, and (iii) the last one ($R_f$), whose value is of the order of one to ten thousand.

In order to optimize the detection of variations in the profiles that are of similar range for each parameter, a scaling can be operated. While the parameters keep their physical values when the cost function is evaluated by solving the direct problem at each iteration, they are scaled while the

minimization is running by means of the CG algorithm. When a scaling is applied, some new parameters are defined,

$$\forall x \in [0, L], \quad \begin{cases} \Lambda_S = 10\,000\Lambda, \\ \Lambda'_S = 10\,000\Lambda', \\ R_{fS} = R_f/1000. \end{cases} \quad (46)$$

## 2. Generating synthetic data

Once the profiles are selected for numerical tests, the effective density $ρ_e$ and $K_e$ are calculated using Eqs. (2) and (3). From the effective characteristic impedance $Z_e=\sqrt{\rho_e K_e}$, the synthetic reflection coefficient $R_m$ is generated.

In a multilayer case, the reflection coefficients $R_m(ω; θ)$ are calculated by a recursive method such as the determination of local surface impedances<sup>7</sup> $Z(x)$ in an imbedding geometry,

$$Z(x_i) = \bar{Z}_e(x_i) \frac{-jZ(x_{i+1})\cotan(kdx) + \bar{Z}_e(x_i)}{Z(x_{i+1}) - j\bar{Z}_e\cotan(kdx)}, \quad (47)$$

$$R_m = R(x_1) = \frac{Z(x_1) - \bar{Z}_0}{Z(x_1) + \bar{Z}_0} \quad (48)$$

wherein $x_i$, $i=1...N_x-1$ is the discretization of $x$ from $x_1=0$ to $x_N=L$, $dx=x_{i+1}-x_i$ the thickness of an infinitesimal layer, $k=ω/\bar{c}_e(x, ω)$ the component of the wave vector along the $x$ axis inside the porous slab in which $\bar{c}_e=ρ_e/\bar{Z}_e$ with $\bar{Z}_e$ from Eq. (22); and $Z(x_N)=\bar{Z}_L$ is the known impedance at $x=L$. Another advantage of generating data by these methods is that the inverse problem can be solved without adding artificial noise. Indeed, the direct solver of the optimization process uses a completely different technique from the solver for generating synthetic data, and thus the “inverse crime”<sup>29</sup> is avoided.

For other inhomogeneous profiles, continuous or not, the synthetic data $R_m$ can be generated by a wave splitting and transmission Green's functions method<sup>9,10</sup> (WS-TGF). The latter is a direct solver that accurately predicts the reflected and transmitted pressure fields, with explicit expressions of the reflection and transmission coefficients, for all the considered frequencies. The WS-TGF method is therefore suitable to generate data in the present work. This method is not similar to the direct solver, except that it uses forward- and backward-propagating waves, e.g., the wave splitting transformation (WS). As explained in Lesselier and Tabbara,<sup>26</sup> the invariant imbedding<sup>21-23,25,27</sup> technique is similar to layer stripping<sup>37</sup> and both are based on iterative resolutions of a Riccati equation as in Eq. (13), whereas the transmission Green's functions (TGFs) technique are linked to the Cholesky factorization of the matrix equation as in (9), which is also similar to the so-called downward continuation.<sup>26,38</sup> Despite these differences between the method of generating data and calculating the field during the optimization, Gaussian noise is added to the synthetic data to simulate experimental error.
<table><caption>TABLE II. Initial values for the reconstructions.</caption><thead><tr><th>φ</th><th>τ<sub>∞</sub></th><th>Λ<br>(μm)</th><th>Λ'<br>(μm)</th><th>R<sub>f</sub><br>(Ns m<sup>-4</sup>)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>100</td><td>300</td><td>1000</td></tr></tbody></table>

### 3. Initial guess of the parameters

The optimization approach consists of retrieving the material parameter profiles by minimizing the objective function $J(p(x))$ from Eq. (24). Unfortunately, in general, we do not know a priori much about the sample. Then, initial values must be chosen with respect to some general knowledge about rigid frame porous media. In such materials, the values of the porosity $φ(x)$ and the tortuosity $τ_∞$ are close to their values in the surrounding fluid, which is one. As to the characteristic lengths, in plastic foams entering into the framework of the rigid frame assumption, we can assume initially that $Λ'(x) ≈ 3Λ(x)$, with $Λ(x)$ varying potentially around 100 µm. The flow resistivity $R_f(x)$ has a wide range of variation, but to be able to use the equivalent fluid model, this parameter must not be too large. It is assumed here that $R_f(x)$ can vary from 1000 N s m⁻⁴ to 20 000 N s m⁻⁴. On account of these remarks, but in somewhat arbitrary fashion, we chose the initial values in Table II for the profile reconstructions.

In practice, the instantaneous reflection at the first interface can provide the material parameter values at $x=0^+$. This corresponds to the first pressure wave recorded in reflection in the time domain, due to the sudden addition of the solid frame at the interface at $x=0$. This method is valid as concerns the porosity and the tortuosity in the high frequency asymptotic model.⁷,¹⁸,²⁰ Then, the values so obtained can be employed as more suitable initial guesses.

To bound the parameter values into the minimizer is similar to add constraints. However, the unconstrained minimization scheme can be kept when an appropriate nonlinear transformation is used. The latter is given by the relation $p_i = f(ψ_i, p_{min_i}, p_{max_i})$, $i=1...5$, such that $p_i$ is the ith component of the vector p, i.e., one of the five acoustical parameters of the porous medium, with $p_{min_i}$ and $p_{max_i}$ the related lower and upper bounds. The parameter $ψ_i$ is the transformed unconstrained variable on which is applied the conjugate gradient algorithm. An infinite number of nonlinear transformations exists³⁹ and one is selected and detailed in Appendix B. The use of bounds can help to retrieve the profiles of the characteristic lengths $Λ(x)$ and $Λ'(x)$ simultaneously with one or several other parameters since only their squared values appear in the equivalent fluid model, in Eqs. (1)–(5). The bounds force the minimizer to seek positive values within a physical meaning.

### 4. Choice of the weighting function $W_R(ω)$ and filtering

The weighting function $W_R(ω)$ is chosen to emphasize the high frequency components in the minimization procedure because these components can reproduce fast variations of the material properties (i.e., because the wavelengths are small). Unfortunately, reconstructions based only on high

frequency components are not accurate enough because of unwanted variations in the retrieved profiles.²⁷ Moreover, the precision of the Runge-Kutta method to solve the differential Eq. (13) decreases when frequency increases. Low frequency components seem to stabilize and smooth the reconstructed profiles. A suitable weighting function is $W_R(ω) = 100^{-(ω/ω_{max})^2}$ or $W_R(ω) = 1000^{-(ω/ω_{max})^2}$ if one wants to reduce the high-frequency contribution.

For most of the reconstructions, the band-limited data induce oscillations in the gradient of the objective function, and therefore in the reconstructed profiles. To limit these artifacts, a sliding average window filter is applied to the gradient along the depth $x$. The length of the window is chosen to be floor($N_x/N_f$) with $N_x$ the number of points to discretize the sample thickness, $N_f=128$ the desired number of windows to apply along the thickness and floor the floor function that converts any real number $a$ to the highest integer less than or equal to $a$. The gradient is averaged in this window and then the window is slid from $x=0$ to $x=L$ with a step $dx$. The result of this procedure is a smoothed gradient and smoothed material parameter profiles, without a noticeable change in the physical information.

### 5. Convergence criterion

In the given conjugate gradient algorithm, the convergence criterion is based on the value of the gradient of the cost function $J$. When the minimum of $J$ is found, the gradient should have converged to zero. However, in practice, a small value $ε$ is chosen such that the algorithm stops when the value of the gradient for all $x$ is below the threshold $ε$. It is difficult to select an optimal criterion $ε$ because the influence of each parameter on the reflection data is different. This means that to screen the variations in the material parameter profiles that do not substantially influence the reflection coefficient a strict criterion is necessary, while a weak criterion is sufficient for the other parameters. To overcome this problem, another convergence criterion is added to the previous one. For all the reconstructions performed in the following, the evolution of the objective function $J$ with the number of iterations $N_{iter}$ is obtained as in Fig. 2. After a certain number of iterations, the cost function $J$ is not decreasing anymore, and its value remains the same for the next iterations. The applied convergence criterion consists of evaluating the variations of $J(p(x))$ for several consecutive iterations. If the variations of $J$ between five consecutive iterations are below $ε$, the global minimum might have been found and the algorithm stops. Five iterations were chosen because it can occur that two or three consecutive evaluations of $J$ are equal. For all the reconstructions, the algorithm stops if the convergence criterion on the gradient, normalized by the number of frequency points $N_ω$, is below $ε=10^{-7}$ or $5 \cdot 10^{-7}$, or if the value of $J$ does not vary more than $ε=10^{-12}$ after five successive iterations, or if the number of iterations exceeds a given number $N_{MaxIter}=100$.

With these ingredients, some reconstructions are performed.
<figure><img src="image_3.png" /><figcaption>FIG. 2. (Left panel) porosity profile reconstruction after 17 iterations for a two-layered medium surrounded and saturated by air. The computation time was 5 min and 11 s. The data were noiseless. The initial guess is $\phi(x)=1$ for all x in $[0, L]$. (Right panel) decrease of the objective function $J(\phi(x))$ as a function of the number of iterations $N_{iter}$. The other parameters are constants from Table I.</figcaption></figure>

### C. Single profile reconstruction

Each of the five acoustical parameters of a virtual mac- roscopically inhomogeneous rigid frame porous material are reconstructed independently and successively.

The frequency band [500 Hz, 500 kHz] is discretized into $N_\omega=200$ frequency points with a logarithmic spacing. In practice, according to the value of $\Lambda'$ from the first layer in Table I, the bandwidth of the incident spectrum signal should not exceed the frequency $f=500$ kHz ($\omega=2\pi f$). This verifi- cation is essential to be sure to avoid the diffusion effect that appears when the wavelength of the acoustic excitation is of the order of the average pore radius, which can be related in a rough approximation to $\Lambda'$. The lower bound is chosen small enough to have a relatively large frequency range. This can help to improve the reconstruction procedure since the objective function $J$ in Eq. (24) and its gradient of compo- nents Eqs. (38)–(42) are summed over a range of frequen- cies.

The differential Eq. (13) is solved using the fourth order Runge–Kutta method. However, the thickness of the sample must be discretized with a small-enough spatial step $dx$ so that, from $x=0$ to $x=L$, at least $N_x=400$ points are needed in order for the direct solver to converge with enough precision. This choice is a bit empirical but, to allow us to change the frequency band and the material thickness, $N_x=800$ points are used.

As a first test, a four-layered material in which only the porosity is varying is simulated with the help of Eqs. (47) and (48). The optimal profile after 17 iterations and 5 min, 11 s of computation time is presented in Fig. 2 together with the evolution of the objective function $J$ with the number of iterations. The optimization was performed on a pentium 4 PC, with a CPU frequency around 2.99 GHz and 1 Gb of RAM memory.

To test the stability of the algorithm, a larger amount of Gaussian noise was added to the synthetic reflection coeffi- cient. Both the porosity profile and the reflection coefficients are plotted in Fig. 3 for a level of Gaussian noise with a standard deviation set to 0.05. For the three-layer material in panel in Fig. 3(a), the maximum number of iterations was set to 50. One of the reconstruction process stopped after 37

iterations because the last five evaluations of the cost func- tion $J(\phi(x))$ gave the same values. This shows that from one computation to another, there is good stability and that the stopping criteria involved are not always the same from one reconstruction to another. The reconstruction of a continuous profile of the porosity with addition of two air layers is given in panel 3(b) with a filtered contaminated reflection coeffi- cient. The convergence is faster, after 21 iterations and 6 min and 34 s. From these minimization results, one remarks that the reconstructed reflection coefficients are smoother than the “measured” ones. The algorithm therefore appears to be stable with respect to additive Gaussian noise.

The other four parameters, tortuosity, viscous and ther- mal characteristic lengths and flow resistivity, can also be reconstructed from data pertaining to single-sided normal re- flection. However, for a given number of iterations, the ac- curacy of the reconstructions is not the same from one pa- rameter to another. For instance, the tortuosity is retrieved with an accuracy similar to the one of the porosity, but using the stronger weighting function $W_R=1000^{-(\omega/\omega_{max})^2}$ and a sliding average window filter to avoid oscillations due to the band-limited data. The use of double-sided reflection data helps to improve the reconstruction, but makes the minimi- zation process slower. Using the same weighting function and filter as the ones used to reconstruct the tortuosity en- ables the characteristic viscous length $\Lambda$ is to be well recon- structed with the use of a rigid wall backing. A slightly better reconstruction is obtained by scaling the parameter $\Lambda_s =10 000\Lambda$. The reconstruction of the characteristic thermal length requires more computational time, but the conver- gence to the true profile in the studied case is not optimal when single-sided reflection data are employed. With a rigid wall backing, the minimization converged only for plane wave excitation at oblique incidence. The difficulty of recon- structing the characteristic thermal length might come from the fact that thermal losses are much less important than viscous losses in most of the porous materials in which the rigid frame assumption is valid. These losses were even ne- glected for a long time in the theory of wave propagation in porous materials.¹⁻⁵ The fifth parameter, the flow resistivity, is also difficult to reconstruct despite the fact that the objec-
<figure><img src="image_4.png" /><figcaption>FIG. 3. Convergence and stability of the conjugate gradient method as a function of the noise level. Each row shows the reconstructed profile of the porosity and the corresponding real part of the reflection coefficient. The maximum number of iterations was set to 50. The added noise is Gaussian, with a standard deviation 0.05. The noise added to the reflection coefficient is not filtered in panel (a) and filtered in panel (b). The other parameters are constants from Table I.</figcaption></figure>

tive function is substantially diminished after a few itera-
tions. This result is mainly due to the poor influence of the
flow resistivity on the reflection data. This parameter mostly
influences low frequency wave propagation attenuation. A
more favorable reconstruction can be carried out when the
lower bound of the frequency band is decreased.

Examples of reconstructed profiles for the tortuosity, the
viscous and thermal characteristic lengths, and the flow re-
sistivity are given in Fig. 4. Results and details of the com-
putational choices are gathered in Table III. According to
these results, the combined use of scaling and no backwall
(Z<sub>L</sub>=Z<sub>0</sub>) at normal incidence seems appropriate to recon-
struct several parameters, except R<sub>f</sub>.

The reconstructed profiles of properties that are discon-
tinuous exhibit large jumps, which amount to singularities
that are the signatures of numerical difficulties. Continuous
parameter profiles, even with discontinuous gradients, induce
a less ill-posed optimization problem allowing successful
and fast reconstructions. To illustrate the last remark, some
examples of single continuous profile (inside the porous
slab) reconstructions are presented in Fig. 5. These profiles,
which are constructed in a manner similar to that of the po-
rosity in Fig. 3(b), yield more satisfactory reconstructions
that are easier to carry out than the ones in Table III and Fig.
4, as we expected. We also studied the case of two air layers,
corresponding to the existence of two jumps in the material

parameter profiles. These jumps were found to be easier to
retrieve than the ones in the profiles of Fig. 4 because air is
a homogeneous fluid whose physical properties are known.
They are introduced in the computation as known param-
eters. The reflection coefficients were generated with the
help of the wave splitting and transmission Green's functions
method (WS-TGF).<sup>9,10</sup> Since this method is considered too
similar to the one (WS-II) of the direct problem solver used
in the minimization process, a small random Gaussian noise
with standard deviation of 0.02 was added to the synthetic
data. This corresponds to a large signal to noise ratio (SNR)
of 100.

A final remark about the reconstruction of a single pa-
rameter: The use of several angles of incidence θ, and of
both backings, increases the available data and therefore the
rate of convergence and the accuracy of the reconstructions.
Several attempts were made to recover Λ'(x) (which is dif-
ficult to reconstruct) without scaling, for two or three angles
of incidence θ=0, θ=π/6 and θ=π/3. The results are shown
in Fig. 6.

D. Simultaneous reconstruction of two profiles

The reconstruction of several material property profiles
simultaneously is more difficult to carry out. Since our aim is
<figure><img src="image_5.png" /><figcaption>FIG. 4. Reconstruction of (a) the tortuosity $\tau_\infty$, (b) the viscous characteristic length $\Lambda$, (c) the thermal characteristic length $\Lambda'$, and (d) the flow resistivity $R_f$. A small amount of Gaussian noise with a standard deviation of 0.01 is added to the data employed for each reconstruction. The performances of the reconstruction procedures are gathered in Table III. For each reconstruction, the other parameters are constants from Table I.</figcaption></figure>

ideally to retrieve five profiles, several records of reflected signals are necessary. Two very different angles of incidence are chosen, $\theta_1=0$ rad and $\theta_2=\pi/3$ rad.

The porosity and the tortuosity are two parameters without dimensions, that have an important influence on the reflection coefficient. For continuous profiles, the simultaneous

reconstruction of these two parameters is quite satisfactory, as depicted in Fig. 7. The reconstruction is obtained after 41 iterations and a bit more than 1/2 h of computation time. The backing is such that $Z_L=Z_0$.

Attempts of simultaneous reconstruction of other couples of material parameters were also made. It was

<table><caption>TABLE III. Convergence of the reconstruction procedure as a function of the number of iterations $N_{iter}$ and computation time. Simulations are performed with or without a rigid wall backing, with or without scaling, and at normal or oblique incidence. N.A. stands for Non Accurate, when reconstructions oscillate too much, or when there are giant values at jump positions. Fails designates reconstructions that do not vary much from the initial guess. "..." means that no computation was done.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Z<sub>L</sub>=Z<sub>0</sub></th><th colspan="2">Z<sub>L</sub>=∞</th><th rowspan="2">Scaling</th><th rowspan="2">θ</th><th rowspan="2">ε</th></tr><tr><th>N<sub>iter</sub></th><th>CPU time</th><th>N<sub>iter</sub></th><th>CPU time</th></tr></thead><tbody><tr><td>φ(x)</td><td>20</td><td>6 min 31 s</td><td>13</td><td>4 min 37 s</td><td>No</td><td>0</td><td>5e<sup>-7</sup></td></tr><tr><td>τ<sub>∞</sub>(x)</td><td>24</td><td>5 min 31 s</td><td>22</td><td>8 min 11 s</td><td>No</td><td>0</td><td>5e<sup>-7</sup></td></tr><tr><td>Λ(x)</td><td>58</td><td>91 min 26 s</td><td>20</td><td>43 min 12 s</td><td>Yes</td><td>0</td><td>5e<sup>-7</sup></td></tr><tr><td></td><td>N.A.</td><td>N.A.</td><td>19</td><td>19 min 52 s</td><td>No</td><td>0</td><td>5e<sup>-7</sup></td></tr><tr><td>Λ'(x)</td><td>21</td><td>67 min 38 s</td><td>N.A.</td><td>N.A.</td><td>Yes</td><td>0</td><td>1e<sup>-7</sup></td></tr><tr><td></td><td>6</td><td>2 min 38 s</td><td>N.A.</td><td>N.A.</td><td>No</td><td>0</td><td>5e<sup>-7</sup></td></tr><tr><td></td><td>...</td><td>...</td><td>12</td><td>23 min 13 s</td><td>Yes</td><td>π/3</td><td>5e<sup>-7</sup></td></tr><tr><td></td><td>...</td><td>...</td><td>11</td><td>9 min 53 s</td><td>No</td><td>π/3</td><td>5e<sup>-7</sup></td></tr><tr><td>R<sub>f</sub>(x)</td><td>Fails</td><td>...</td><td>Fails</td><td>...</td><td>Yes</td><td>0</td><td>5e<sup>-7</sup></td></tr><tr><td></td><td>8</td><td>13 min 26 s</td><td>32</td><td>91 min 40 s</td><td>No</td><td>0</td><td>5e<sup>-7</sup></td></tr><tr><td></td><td>49</td><td>96 min 41 s</td><td>...</td><td>...</td><td>No</td><td>0</td><td>1e<sup>-7</sup></td></tr><tr><td></td><td>...</td><td>...</td><td>31</td><td>99 min 44 s</td><td>No</td><td>π/3</td><td>5e<sup>-7</sup></td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Reconstruction of (a) the tortuosity $\tau_\infty$ after 30 iterations and 10 min, 44 s of computation time, (b) the viscous characteristic length $\Lambda$ after 30 iterations and 26 min, 41 s, (c) the thermal characteristic length $\Lambda'$ after 80 iterations and 48 min, 55 s, and (d) the flow resistivity $R_f$ after 30 iterations and 55 min, 11 s. For each reconstruction, a small amount of Gaussian noise with a standard deviation of 0.02 is added to the data. For each reconstruction, the other parameters are constants from Table I.</figcaption></figure>

pointed out previously that the characteristic lengths $\Lambda$ and $\Lambda'$ require more computation time to sufficiently diminish the cost function $J$. Moreover, as noticed for the reconstruction of the porosity-tortuosity couple, the computation time needed to retrieve several profiles simultaneously increases drastically. The way to cope with this problem is to decrease

<figure><img src="image_7.png" /><figcaption>FIG. 6. Reconstruction of $\Lambda'(x)$ with the use of two angles of incidence ($\theta=0$ and $\pi/3$) or three angles of incidence ($\theta=0, \pi/6$ and $\pi/3$) with a rigid backing. The number of iterations is respectively 10 and 11, which required 15 min, 2 s and 17 min, 3 s. The data included a small amount of Gaussian noise with a standard deviation of 0.01. The other parameters are constants from Table I.</figcaption></figure>

the number of points $N_x$ from 800 to 400 for discretizing the material thickness, and decreasing the higher bound of the frequency range. Thus, only 150 frequency points with a logarithmic spacing are chosen to describe the frequency band [400 Hz, 400 kHz]. In Fig. 8 is presented the simultaneous reconstruction of the porosity $\phi$ and the viscous characteristic length $\Lambda$ when the medium is backed by a rigid wall ($Z_L = \infty$).

<figure><img src="image_8.png" /><figcaption>FIG. 7. Simultaneous reconstruction of the porosity and the tortuosity. The number of frequency points is 200, the number of depth points is 800. The reconstruction is given after 40 iterations that took 30 min, 11 s. Addition of a weak Gaussian noise with a standard deviation of 0.01. The other parameters are constants from Table I.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 8. Simultaneous reconstruction of the porosity $\phi$ and the characteristic viscous length $\Lambda$ after 100. The computation time was 55 min and 11 s. The number of points was reduced and the upper bound of the frequency band set to 400 kHz. A small amount of Gaussian noise with a standard deviation of 0.01 was added to the synthetic data. The other parameters are constants from Table I.</figcaption></figure>

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

<figure><img src="image_12.png" /><figcaption>FIG. 9. Simultaneous reconstruction of the porosity $\phi$, tortuosity $\tau_\infty$ and characteristic length $\Lambda$ after 100 iterations that took 118 min and 19 s. 150 frequency points and 400 equally spaced points were employed in a frequency band [500 Hz, 500 kHz]. A small amount of Gaussian noise with a standard deviation of 0.01 was added to the synthetic data. The other parameters are constants from Table I.</figcaption></figure>

Satisfactory reconstructions can be obtained in a reasonable amount of time with reduced sets of points to discretize the frequency range as well as the material thickness. However, the optimization procedure is sensitive to these discretizations. The initial guesses of the minimizer are also very important to accelerate the convergence. The reconstructions are faster when they are closer to their true values.

## E. Simultaneous reconstruction of more than two material parameter profiles

Employing the reflection data from acoustical excitations at only three different angles of incidence, $\theta_1=0$ rad, $\theta_2=\pi/6$ rad and $\theta_3=\pi/3$ rad, a rigid wall at $x=L$ and a wide-enough frequency band [400 Hz, 400 kHz], only 150 frequency points and 400 equally spaced points along the material thickness are needed to reconstruct three profiles simultaneously, as is shown in Fig. 9. The chosen parameters are the porosity, tortuosity, and viscous characteristic length. They are the most sensitive parameters with respect to the reflection coefficient. Moreover, they are also the parameters of interest in the high frequency asymptotic model for rigid frame porous plastic foams.⁷,¹⁸,²⁰ The flow resistivity $R_f(x)$ has such a small influence that it vanishes, and the thermal characteristic length $\Lambda'(x)$ is chosen such that $\Lambda'(x) = 3\Lambda(x)$ or $\Lambda'(x) = 2\Lambda(x)$ according to the foam.

The computation time exceeds 1 h, but is reasonable for such a type of optimization. The results are presented after 100 iterations. Nevertheless, the main variations of the material parameters are already visible after 50 iterations. With experimental data, the expected precision of the reconstructions can be expected to be much less satisfactory.

The result of an attempt to simultaneously reconstruct the same three parameters for a material with strong discontinuities, i.e., a three-layered material, is depicted in Fig. 10. The reconstruction is carried out by minimizing data pertaining to four reflection coefficients. The “measured” ones are obtained by the recursive method from Eqs. (47) and (48). Three of them are simulated without backwall ($Z_L=Z_0$) at three different angles of incidence ($\theta=\pi/6$ rad, $\pi/4$ rad and $\pi/3$ rad) and the fourth one is simulated with a backwall ($Z_L=\infty$) at the angle $\theta=\pi/4$ rad. The frequency band is [1 kHz, 400 kHz], discretized into $N_\omega=300$ equispaced points. The depth is discretized into $N_x=400$ points. The reconstruction took 5 h, 49 min, 47 s and 105 iterations. The result is a bit less accurate and is more difficult to carry out than with a continuous profile. There is a good enough reproductibility but not as satisfactory as with continuous profiles (due to the directions taken by the CG iterations).

The reconstruction of more than three profiles simultaneously can be achieved in a similar way. To improve con-
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 10. Simultaneous reconstruction of the porosity $\phi$, tortuosity $\tau_\infty$ and characteristic length $\Lambda$ after 105 iterations that took 5 h, 49 min. and 47 s; 300 frequency points and 400 equally spaced points were used in a frequency band [1 kHz, 400 kHz]. The data were noiseless. The other parameters are constants from Table I.</figcaption></figure>

vergence, one can use more sets of data than the number of parameters we aim to reconstruct. For instance, to satisfactorily reconstruct four profiles simultaneously, it can happen that four reflections coefficient spectra at four different angles of incidence are not sufficient. One can then employ data relative to more than four angles of incidence, or one can add measurements with air for all $x > L$ instead of a rigid wall, or one can also impose a rigid wall at $x=0$ while acoustically exciting the sample at $x=L$.

## VI. CONCLUSION

An optimization approach to the frequency domain inverse problem is described in this paper in order to reconstruct simultaneously the profiles of several material parameters such as the porosity, the tortuosity, the viscous and thermal characteristic lengths, and the flow resistivity in any unidimensional rigid frame macroscopically inhomogeneous porous materials.

The method is based on the minimization of a suitable objective function, chosen such that its gradient can be analytically obtained. From this amount of information the direct problem is solved iteratively by replacing at each step the values of the parameters to reconstruct by their optimized values. The direct solver is based on the wave splitting and invariant imbedding methods. The optimization approach is based on the conjugate gradient method. Reconstructions of

one, two, and even three profiles are obtained with accuracy after only few iterations. The computation time is reasonable in most of the cases. The simultaneous reconstruction of the five profiles is also discussed. The use of several impedances at one end point ($x=L$) and the use of several angles of incidence helps to carry out simultaneous reconstructions based on the reflection data only.

Better optimized results might be obtained when a preconditioner $M$ is used.33 The latter is a diagonal matrix whose components are those of the main diagonal of the Hessian matrix $H=(\partial^2 J(p(x))/\partial p_i \partial p_j)$, where $p_i, i=1...5$, refers to the five material parameters to reconstruct. When a preconditioner is used, the direction search $D^{(i)}$ in the conjugate gradient algorithm is replaced by $MD^{(i)}$. When the signal to noise ratio is weak and the variations in the properties to reconstruct are important (i.e., jumps), the optimization problem can become highly ill posed and regularization methods such as Tikhonov regularization40 might be applied.

In a further work, some questions about the reconstruction of material profiles from measured data need to be solved. For instance, the frequency bandwidths of most of the transducers are not as large as the frequency band used to reconstruct profiles in this article. Therefore, in order to apply the technique detailed in the present paper, numerical methods to extract the necessary large band reflection coefficient are under study.

## ACKNOWLEDGMENT

The authors are grateful to Professor M. Norgren for his advice about the present optimization approach.

## APPENDIX A: ITERATIVE SOLUTIONS OF LINEAR FIRST ORDER DIFFERENTIAL EQUATIONS

While the Riccati Eq. (13) is solved with a Runge-Kutta method, the other first order differential Eqs. (14), (35) are linear and then easier and faster to solve. An analytical solution of Eq. (14) is easy to carry out. Nevertheless, it may be not the faster way to compute the solution, because for each $x$, an integral is evaluated. A fast and easy to compute iterative method extracted from the analytical solution is preferred.28 The technique is general and just a detailed example is given for Eq. (14). The equation is separated as follows in order to have only $dx$ on the right hand side:

$$ \frac{dT}{(A^+ + A^- R)T} = dx. \qquad (A1) $$

Integrating from $x$ to $x+dx$, and denoting $x_i=x$, $x+dx = x_{i+1}$, $f(x)=f_i$ and $f(x+dx)=f_{i+1}$, with $f$ any function of $x$, one gets

$$ \int_{T_i}^{T_{i+1}} \frac{dT}{(A_i^+ + A_i^- R_i)T} = \int_{x_i}^{x_{i+1}} dx, \qquad (A2) $$

$$ [\ln T]_{T_i}^{T_{i+1}} = (A_i^+ + A_i^- R_i) \Delta x, \qquad (A3) $$

$$ T_i = T_{i+1} e^{-(A_i^+ + A_i^- R_i) \Delta x}, \qquad (A4) $$

where $\Delta x = x_{i+1} - x_i$.
The recurrence relation to solve Eq. (35) is

$$u_{i+1} = u_i e^{-2(A_i^+ + A_i^- R_i) \Delta x}. \quad (A5)$$

## APPENDIX B: NONLINEAR TRANSFORMATION TO ADD BOUNDS TO THE UNCONSTRAINED MINIMIZER

Instead of performing an unconstrained minimization, the addition of bounds on the model parameters may help the optimizer to converge faster to the desired profiles. The bounds in the present case must at least force the algorithm to seek positive real values for each component of the vector $p(x)$. This gives a lower bound vector $p_{\min}$. Higher bounds, in $p_{\max}$, can be guessed according to the physics of porous materials in the rigid frame assumption. However, they should be taken large enough to ensure the accounting of unexpected values to inhomogeneities.

Calling $p_i$, $i=1, \dots, 5$, the ith component of the vector $p(x)=[\phi(x), \tau_{\infty}(x), \Lambda(x), \Lambda'(x), R_f(x)]^T$, an unconstrained vector of components $\psi_i$ in $]-\infty, \infty[$ is defined such that

$$p_i = f(\psi_i, p_{\min i}, p_{\max i}) \\ = \frac{p_{\max i} + p_{\min i}}{2} + \frac{p_{\max i} - p_{\min i}}{2} \sin(\psi_i), \quad (B1)$$

with

$$p_i \rightarrow p_{\min i}, \quad \text{as} \quad \sin(\psi_i) \rightarrow -1, \\ p_i \rightarrow p_{\max i}, \quad \text{as} \quad \sin(\psi_i) \rightarrow +1.$$

The conjugate gradient algorithm is applied to the transformed unconstrained vector $\psi$. At the kth iteration, the real values of the parameters must be used in the model, and

$$p_i^{(k+1)} = \frac{p_{\max i} + p_{\min i}}{2} + \frac{p_{\max i} - p_{\min i}}{2} \sin(\psi_i^{(k)} - \lambda_k D_i^{(k)}), \quad (B2)$$

with $\lambda_k$ a positive real value minimizing the objective function $J(p)$, the gradient of which, in terms of $\psi_i$, is linked to the search direction $D_i^{(k)}$. Straightforwardly, the components of the gradient are

$$\frac{\partial J}{\partial \psi_i} = \sqrt{(p_{\max i} - p_i)(p_i - p_{\min i})} \frac{\partial J}{\partial p_i}. \quad (B3)$$

**Remark:** In Eq. (B3), it is obvious that the initial guess cannot be set equal to one of the bound. This avoids a non-desired null gradient that implies errors in the Polak-Ribière direction update (45) in the CG algorithm.

<footnote n="1">¹M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated porous solid - I. low frequency range and II. high frequency range.” J. Acoust. Soc. Am. **28**, 168–178 (1956); and **28**, 179–191 (1956).</footnote>

<footnote n="2">²M. A. Biot and D. G. Willis, “The elastic coefficients of the theory of consolidation.” J. Appl. Mech. **24**, 594–601 (1957).</footnote>

<footnote n="3">³M. A. Biot, “Mechanics of deformation and acoustic propagation in porous media.” J. Appl. Phys. **33**, 1482–1498 (1962).</footnote>

<footnote n="4">⁴M. A. Biot, “Generalized theory of acoustic propagation in porous dissipative media.” J. Acoust. Soc. Am. **34**, 1254–1264 (1962).</footnote>

<footnote n="5">⁵D. L. Johnson, J. Koplik, and R. Dashen, “Theory of dynamic permeability and tortuosity in fluid-saturated porous media.” J. Fluid Mech. **176**,</footnote>

<footnote>379–402 (1987).</footnote>

<footnote n="6">⁶Y. Champoux and J. F. Allard, “Dynamic tortuosity and bulk modulus in air-saturated porous media.” J. Appl. Phys. **70**, 1975–1979 (1991).</footnote>

<footnote n="7">⁷J. F. Allard, *Propagation of Sound in Porous Media: Modeling Sound Absorbing Materials* (Chapman and Hall, London, 1993).</footnote>

<footnote n="8">⁸K. Attenborough, “Acoustical characteristics of porous materials.” Phys. Rep. **82**, 179–227 (1982).</footnote>

<footnote n="9">⁹L. De Ryck, J. P. Groby, P. Leclaire, W. Lauriks, A. Wirgin, Z. E. A. Fellah, and C. Depollier, “Acoustic wave propagation in a macroscopically inhomogeneous porous medium saturated by a fluid.” Appl. Phys. Lett. **90**, 181901 (2007).</footnote>

<footnote n="10">¹⁰L. De Ryck, W. Lauriks, Z. E. A. Fellah, A. Wirgin, J. P. Groby, P. Leclaire, and C. Depollier, “Acoustic wave propagation and internal fields in rigid frame macroscopically inhomogeneous porous media.” J. Appl. Phys. **102**, 024910 (2007).</footnote>

<footnote n="11">¹¹P. Leclaire, L. Kelders, W. Lauriks, N. R. Brown, M. Melon, and B. Castagnède, “Determination of the viscous and thermal characteristic lengths of plastic foams by ultrasonic measurements in helium and air.” J. Appl. Phys. **80**, 2009–2012 (1996).</footnote>

<footnote n="12">¹²J. P. Groby, L. De Ryck, P. Leclaire, A. Wirgin, W. Lauriks, R. P. Gilbert, and Y. S. Xu, “Use of specific Green’s functions for solving direct problems involving a heterogeneous rigid frame porous medium slab solicited by acoustic waves.” Math. Methods Appl. Sci. **30**, 91–122 (2007).</footnote>

<footnote n="13">¹³O. Forslund and S. He, “Electromagnetic scattering from an inhomogeneous grating using a wave splitting approach,” Prog. Electromagn. Waves **19**, 147–171 (1998).</footnote>

<footnote n="14">¹⁴L. L. Beranek, “Acoustic impedance of porous materials,” J. Acoust. Soc. Am. **13**, 248–260 (1942).</footnote>

<footnote n="15">¹⁵R. L. Brown and R. H. Bolt, “The measurement of flow resistance of porous acoustic materials,” J. Acoust. Soc. Am. **13**, 337–344 (1942).</footnote>

<footnote n="16">¹⁶J.-F. Allard, B. Castagnede, M. Henry, and W. Lauriks, “Evaluation of tortuosity in acoustic porous materials saturated by air,” Rev. Sci. Instrum. **65**, 754–755 (1994).</footnote>

<footnote n="17">¹⁷N. R. Brown, M. Melon, B. Castagnede, V. Montembault, W. Lauriks, and P. Leclaire, “Evaluation of the viscous characteristic length of air-saturated porous materials from the ultrasonic dispersion curve,” C. R. Acad. Sci. URSS **322**, 121–127 (1996).</footnote>

<footnote n="18">¹⁸Z. E. A. Fellah, S. Berger, W. Lauriks, C. Depollier, C. Aristegui, and J-Y. Chapelon, “Measuring the porosity and the tortuosity of porous materials via reflected waves at oblique incidence,” J. Acoust. Soc. Am. **114**, 2424–2432 (2003).</footnote>

<footnote n="19">¹⁹O. Umnova, K. Attenborough, H. C. Shin, and A. Cummings, “Deduction of tortuosity and porosity from reflection and transmission measurements on thick samples of porous materials,” Appl. Acoust. **66**, 607–624 (2005).</footnote>

<footnote n="20">²⁰Z. E. A. Fellah, C. Depollier, S. Berger, W. Lauriks, C. Aristegui, P. Trompette, and J-Y. Chapelon, “Determination of transport parameters in air-saturated porous materials via reflected ultrasonic waves,” J. Acoust. Soc. Am. **114**, 2561–2569 (2003).</footnote>

<footnote n="21">²¹G. Kristensson and R. J. Krueger, “Direct and inverse scattering in the time domain for a dissipative wave equation. i. scattering operators and ii. simultaneous reconstruction of dissipation and phase velocity profiles,” J. Math. Phys. **27**, 1667–1682 and **27**, 1683–1693 (1986).</footnote>

<footnote n="22">²²G. Kristensson and R. J. Krueger, “Direct and inverse scattering in the time domain for a dissipative wave equation. iii. scattering operators in the presence of a phase velocity mismatch,” J. Math. Phys. **28**, 360–370 (1987).</footnote>

<footnote n="23">²³G. Kristensson and R. J. Krueger, “Direct and inverse scattering in the time domain for a dissipative wave equation. iv. use of phase velocity mismatches to simplify inversions,” Inverse Probl. **5**, 375–388 (1989).</footnote>

<footnote n="24">²⁴J. P. Corones, M. E. Davison, and R. J. Krueger, “Direct and inverse scattering in the time domain via invariant imbedding equations,” J. Acoust. Soc. Am. **74**, 1535–1541 (1983).</footnote>

<footnote n="25">²⁵M. Gustafsson and S. He, “An optimization approach to multi-dimensional time domain acoustic inverse problems,” J. Acoust. Soc. Am. **108**, 1548–1556 (2000).</footnote>

<footnote n="26">²⁶D. Lesselier and W. Tabbara, “Problèmes unidimensionnels de diffraction inverse en acoustique et électromagnétisme. étude bibliographique (Unidimensional inverse scattering problems in acoustic and electromagnetism. Bibliographical study),” J. Acoust. **1**, 363–384 (1988).</footnote>

<footnote n="27">²⁷M. Norgren and S. He, “An optimization approach to the frequency-domain inverse problem for a nonuniform lcg transmission line,” IEEE Trans. Microwave Theory Tech. **44**, 1503–1507 (1996).</footnote>

<footnote n="28">²⁸C. MacKay, D. Hayward, S. McKee, A. J. Mulholland, and R. A. Pethrick, “An inverse problem of reconstructing the electrical and geometrical pa-</footnote>
rameters characterizing airframe structures and connector interfaces,” In-
verse Probl. Eng. **15**, 177–190 (2007).

<footnote n="29"><sup>29</sup>A. Wirgin, *The Inverse Crime*. Available at http://hal.archives-ouvertes.fr/ hal-00001084/fr/, last viewed online on 04/03/2008.</footnote>

<footnote n="30"><sup>30</sup>R. J. Krueger and L. Ochs, “A Green’s function approach to the determination of internal fields,” *Wave Motion* **11**, 525–543 (1989).</footnote>

<footnote n="31"><sup>31</sup>J. le Rond d'Alembert, *Opuscules Mathématiques ou Mémoires sur Différents Sujets de Géométrie, de Méchanique, d'Optique, d'Astronomie (Mathematical Opus or Manuscripts on Several Subjects About Geometry, Mechanics, Optics and Astronomy)* (Chez David, Paris, 1761), pp. 1761–1780. Available at http://gallica.bnf.fr/ark:/12148/bpt6k62394p/ last viewed online on 7/31/08.</footnote>

<footnote n="32"><sup>32</sup>S. He, S. Ström, and V. H. Weston, *Time Domain Wave-Splittings and Inverse Problems* (Oxford University Press, Oxford, 1998).</footnote>

<footnote n="33"><sup>33</sup>J. R. Shewchuk, *An Introduction to the Conjugate Gradient Method Without the Agonizing Pain*. Available at www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf, last viewed online on 04/04/2008.</footnote>

<footnote n="34"><sup>34</sup>W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in Fortran* 77, 2nd ed. (Cambridge University Press, Cambridge, 1992).</footnote>

<footnote n="35"><sup>35</sup>E. Polak, *Computational Methods in Optimization* (Academic, New York, London, 1971).</footnote>

<footnote n="36"><sup>36</sup>J. A. Nelder and R. Mead, “A simplex method for function minimization,” *Comput. J.* **7**, 308–313 (1965).</footnote>

<footnote n="37"><sup>37</sup>A. E. Yagle, “A layer stripping fast algorithm for reconstructing a lossy layered medium,” *J. Acoust. Soc. Am.* **82**, 927–936 (1987).</footnote>

<footnote n="38"><sup>38</sup>K. P. Bube and R. Burridge, “The one-dimensional inverse problem of reflection seismology,” *SIAM Rev.* **25**, 497–559 (1983).</footnote>

<footnote n="39"><sup>39</sup>T. M. Habashy and A. Abubakar, “A general framework for constraint minimization for the inversion of electromagnetic measurements,” *Electromagn. Waves* **46**, 265–312 (2004).</footnote>

<footnote n="40"><sup>40</sup>A. Franchois and C. Pichot, “Microwave imaging-complex permittivity reconstruction with a Levenberg-Marquardt method,” *IEEE Trans. Antennas Propag.* **45**, 203–215 (1997).</footnote>