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

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


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 , which is the ratio of the fluid volume to the total sample volume; the tortuosity ϱ , 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 ⌳ and ⌳Ј, 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 to the fluid permeability f .6][7][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,100,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 9,10 by means of the wave splitting and transmission Green's functions method.13 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͒ 14 ; another one concerns the simultaneous measurement of the pressure drop accross a sample and the rate of fluid flow through a porous sample ͑flowmeter for the determination of the flow resitivity͒ 15 ; more recent methods deal with ultrasonic measurements in both time and frequency domains to retrieve the following parameters: the tortuosity, 16 the viscous 17 and thermal characteristic lengths ͑the Q␦ method͒, 11 the porosity and the tortuosity simultaneously, 18,19 and finally the porosity, the tortuosity, and the viscous characteristic length simultaneously. 20The 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 R f 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.][23][24][25][26] More recently, similar methods, together with an optimization approach, were developed in the frequency domain. 27,28he 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. 24,27The 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." 29The 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. 9,10,12][3][4] 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 K e of a so-called equivalent fluid. 7he model of rigid frame has been recently extended to macroscopically inhomogeneous porous media. 9,10,12In the frequency domain, the Helmholtz equation in terms of the fluid pressure p inside the equivalent inhomogeneous fluid is 9,10 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 , ͒ = ͱ K e ͑x , ͒ / e ͑x , ͒ and Z e = e ͑x , ͒c e ͑x , ͒.
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 1-7 extended to macroscopically inhomogeneous porous materials: 9,10

͑3͒
wherein ␥ is the specific heat ratio, f the saturating fluid density, P 0 the atmospheric pressure and P r the Prandtl number.The well-defined correction functions of this model are [5][6][7] F͑x,͒ = ͱ1+j4 f ϱ 2 ͑x͒ The five so-called acoustical parameters ͑x͒, ϱ ͑x͒, ⌳͑x͒, ⌳Ј͑x͒, and R f ͑x͒ described in the introduction become here space-dependent functions. 9,10Without 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. 6,7he 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͒: 9,10 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 onedimensional ͑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.
,13,27 In order to avoid convolution products in time and fractional derivatives, 18,20 the direct and inverse scattering problems are treated in the frequency domain.Also, the Fourier transform p͑x , ͒ of the measurable pressure field p͑x , t͒ is used.The inverse Fourier transform is chosen such that p͑x,t͒ = ͵ −ϱ ϱ p͑x,͒exp͑jt͒d.

͑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 tortuosity 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 9,10,12,27 , 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 = 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.23][24]32 The transformation is gener- 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 .In a homogeneous medium, p − are the backward propagating waves and p + are the forward propagating waves.alized for inhomogeneous media since it can be seen as an algebraic transformation from the base of functions ͓p ; ͑x͒V͔ to the base of functions ͓p + ; p − ͔.Several transformations are available.In the present case, the so-called vacuum wave splitting tranformation 9,10,27,13 is chosen,

͑8͒
This choice is possible because the wave splitting is coupled with an invariant imbedding technique.2][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 = f c f of the saturating fluid instead. 27,28s 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 implicitely accounted for.
Applying the transformation ͑8͒ to the equations of motion ͑6͒ in the frequency domain straightforwardly yields 9,10

͑10͒
For all x Յ 0, the pressures waves p + ͑x , ͒ and p − ͑x , ͒ are respectively the incident and reflected pressure waves.Similarly, for all x Ն L, p + ͑x , ͒ corresponds to the transmitted pressure wave.Therefore, the splitted pressure waves p Ϯ are easily connected to the reflection and transmission coefficients R͑͒ and T͑͒.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 → 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 , ͒ and T͑x , ͒ as follows: From the above relations, one deduces that R͑0,͒ = R͑͒ and T͑0,͒ = T͑͒ 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 ͑͒ = p − ͑0,͒ and p t ͑͒ = p + ͑L , ͒, to the incident pressure field p i ͑͒ = p + ͑0,͒ 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 , ͒ and T͑x , ͒, 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 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 If the homogeneous medium for all x Ͼ L is a rigid wall, then Z L = ϱ.There is no transmitted waves, and only the imbedding Eq. ͑13͒ for R͑x , ͒ is used.The boundary condition reduces to 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 latters, for an angle of incidence , the effective bulk modulus K e , the characteristic impedances Z 0 and Z e must be replaced by 10 . ͑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 28 and may even trap the optimizer. 27Therefore, 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 where the superscript T denotes the transpose operation.An objective function J͑p͒ is defined such that 27 where R m is the measured reflection coefficient, R͑0, ; ͒ 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 min and max , and over a frequency band ͓ min , 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͑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͑p͒ from an initial guess of the elements of the vector p.The resolution of our problem, involving the minimization of the function J͑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. 25,27This 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 ␦p is applied to the material property vector p.By definition, the perturbation of the reflection coefficient is ␦R = R ˜− R where R ˜= R͑x , ; p + ␦p͒ is the solution of the perturbed Eq. ͑13͒, which becomes: 27 d dx with the perturbed condition ͑15͒ reducing to ␦R͑L,͒ = 0 ͑26͒ because the impedances Z 0 = 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 Ϯ are deduced from Eq. ͑10͒ and are expressed, for any oblique angle of incidence, as wherein the perturbation of e ͓Eq.͑2͔͒ and of K ¯e −1 ͓Eq.͑20͔͒ yield

͑31͒
where Re represents for the real part and wherein the auxiliary function u͑x , ; ͒ is defined such that in which the superscript Ã denotes the complex conjugate.Using the boundary conditions ͑26͒ at x = L, the following integration rule can be applied: the integrand of which, from Eq. ͑25͒, is developed as follows: u͑x , ; ͒ is an arbitrary function, except at x = 0, so that we can choose it in order to eliminate the dependence on the unknown perturbation ␦R in Eq. ͑34͒.Thus, the following linear first order differential equation must be satisfied: 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 and on the other hand, the same perturbation can be expressed as the inner product

͑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

͑42͒
In the above relations, all partial derivatives are functions of the components of the vector p = ͓͑x͒ , ϱ ͑x͒ , ⌳͑x͒ , ⌳Ј͑x͒ , R f ͑x͔͒ T via the definitions of 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.][35] 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 25 or transmission line theory 27 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͑p͑x͒͒ ͑24͒ is not a quadratic form and continuous in 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 25,[33][34][35] for minimization.For clear and rigorous details about the method, the reader can refer to a tutorial ͑with provided algorithms͒. 33For CG algorithms, as well as other optimization techniques, the reader can refer to Polak 35 and Press et al. 34 Herein, we describe the CG algorithm for the reconstruction of a single parameter p͑x͒, e.g., the porosity ͑x͒.The reconstruction of the other parameters can be simply derived in closed form.The gradient ‫ץ‬J / ‫ץ‬p is denoted by G͑p͒.The implementation is as follows: 25,33,34 • 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 with the help of Eq. ͑38͒; For all x in ͓0,L͔, For any iteration order i, D ͑i͒ is the search direction. with The minimization over the parameter is easily achieved by means of the Nelder and Mead Simplex method. 34,36This 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.

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. 11When 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, ∀x ͓0,L͔, Ά ⌳ S = 10 000⌳,

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 = ͱ 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 7 Z͑x͒ in an imbedding geometry, x i the thickness of an infinitesimal layer, k = / c ¯e͑x , ͒ the component of the wave vector along the x axis inside the porous slab in which c ¯e = e / Z ¯e with Z ¯e from Eq. ͑22͒; and Z͑x N x ͒ = 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" 29 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 9,10 ͑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, 26 the invariant imbedding [21][22][23]25,27 technique is similar to layer stripping 37 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. 26,38Despite 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.

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 −4 to 20 000 N s m −4 .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. 7,18,20Then, 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 39 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.

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. 27Moreover, 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.

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.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.

C. Single profile reconstruction
Each of the five acoustical parameters of a virtual macroscopically inhomogeneous rigid frame porous material are reconstructed independently and successively.
The frequency band ͓500 Hz, 500 kHz͔ is discretized into N = 200 frequency points with a logarithmic spacing.In practice, according to the value of ⌳Ј from the first layer in Table I, the bandwidth of the incident spectrum signal should not exceed the frequency f = 500 kHz ͑ =2f͒.This verification 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 ⌳Ј.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 components Eqs.͑38͒-͑42͒ are summed over a range of frequencies.
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 coefficient.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 function J͑͑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 coefficient.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 thermal characteristic lengths and flow resistivity, can also be reconstructed from data pertaining to single-sided normal reflection.However, for a given number of iterations, the accuracy of the reconstructions is not the same from one parameter 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 −͑/ 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 minimization process slower.Using the same weighting function and filter as the ones used to reconstruct the tortuosity enables the characteristic viscous length ⌳ is to be well reconstructed with the use of a rigid wall backing.A slightly better reconstruction is obtained by scaling the parameter ⌳ S = 10 000⌳.The reconstruction of the characteristic thermal length requires more computational time, but the convergence 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 reconstructing 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.2][3][4][5] The fifth parameter, the flow resistivity, is also difficult to reconstruct despite the fact that the objec- tive function is substantially diminished after a few iterations.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 resistivity are given in Fig. 4. Results and details of the computational choices are gathered in Table III.According to these results, the combined use of scaling and no backwall ͑Z L = Z 0 ͒ at normal incidence seems appropriate to reconstruct several parameters, except R f .
The reconstructed profiles of properties that are discontinuous 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 porosity 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 parameters.The reflection coefficients were generated with the help of the wave splitting and transmission Green's functions method ͑WS-TGF͒. 9,10Since 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 parameter: 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 difficult 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 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.
ideally to retrieve five profiles, several records of reflected signals are necessary.Two very different angles of incidence are chosen, 1 = 0 rad and 2 = / 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 FIG. 4. Reconstruction of ͑a͒ the tortuosity ϱ , ͑b͒ the viscous characteristic length ⌳, ͑c͒ the thermal characteristic length ⌳Ј, 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. pointed out previously that the characteristic lengths ⌳ and ⌳Ј 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 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 and the viscous characteristic length ⌳ when the medium is backed by a rigid wall ͑Z L = ϱ ͒.I 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, 1 = 0 rad, 2 = / 6 rad and 3 = / 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. 7,18,20The flow resistivity R f ͑x͒ has such a small influence that it vanishes, and the thermal characteristic length ⌳Ј͑x͒ is chosen such that ⌳Ј͑x͒ =3⌳͑x͒ or ⌳Ј͑x͒ =2⌳͑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 ͑ = / 6 rad, / 4 rad and / 3 rad͒ and the fourth one is simulated with a backwall ͑Z L = ϱ ͒ at the angle = / 4 rad.The frequency band is ͓1 kHz, 400 kHz͔, discretized into N = 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- FIG. 9. Simultaneous reconstruction of the porosity , tortuosity ϱ and characteristic length ⌳ 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. 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. 33The latter is a diagonal matrix whose components are those of the main diagonal of the Hessian matrix H = ‫ץ͑‬ 2 J͑p͑x͒͒ / ‫ץ‬p i ‫ץ‬ 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 regularization 40 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. 28The 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: 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 dx, ͑A2͒ ͓ln T͔ T i T i+1 = ͑A i + + A i − R i ͒⌬x, ͑A3͒ T i = T i+1 e −͑A i where ⌬x = x i+1 − x i .I.
The recurrence relation to solve Eq. ͑35͒ is ͑A5͒

FIG. 2 .
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 ͑x͒ = 1 for all x in ͓0,L͔.͑Right panel͒ decrease of the objective function J͑͑x͒͒ as a function of the number of iterations N iter .The other parameters are constants from TableI.

FIG. 5 .FIG. 6 .
FIG.5.Reconstruction of ͑a͒ the tortuosity ϱ after 30 iterations and 10 min, 44 s of computation time, ͑b͒ the viscous characteristic length ⌳ after 30 iterations and 26 min, 41 s, ͑c͒ the thermal characteristic length ⌳Ј 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 TableI.

FIG. 8 .
FIG.8.Simultaneous reconstruction of the porosity and the characteristic viscous length ⌳ 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 TableI.

FIG. 10 .
FIG.10.Simultaneous reconstruction of the porosity , tortuosity ϱ and characteristic length ⌳ 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 TableI.

TABLE I .
Values of the porous material at x =0.

TABLE II .
Initial values for the reconstructions.

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. .
φ 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 TableI.