N Totaro

D Vigoureux

Q Leclère

J Lagneaux

J L Guyader

Sound elds separation and reconstruction of irregularly shaped sources

Keywords: source identication, underdetermined inverse problem, source separation, Tikhonov regularization, Generalized Cross Validation 1

Nowadays, the need of source identication methods is still growing and application cases are more and more complex. As a consequence, it is necessary to develop methods allowing to reconstruct sound elds on irregularly shaped sources in reverberant or conned acoustic environment. The inverse Patch Transfer Functions (iPTF) method is suitable to achieve these objectives. Indeed, as the iPTF method is based on Green's identity and double measurements of pressure and particle velocity on a surface surrounding the source, it is independent of the acoustic environment. In addition, the nite element solver used to compute the patch transfer functions permits to handle sources with 3D irregular shapes. In the present article, two experimental applications on a at plate and an oil pan have been carried out to show the performances of the method on real applications. As for all ill-posed problem, it is shown that the crucial point of this method is the choice of the parameter of the Tikhonov regularization, one of the most widely used in the literature. The classical L-curve strategy sometimes fails to choose the best solution. This issue is clearly explained and an adapted strategy combining L-curve and acoustic power conservation is proposed. The eciency of this strategy is demonstrated on both applications and compared to results obtained with Generalized Cross Validation (GCV) technique.

Introduction

In acoustics, inverse methods cover a wide eld of applications. Nevertheless, they all rely on the same concept : identifying the inputs of a linear system by observing the outputs. If the observed outputs are in most cases the pressure and/or particle velocity elds radiated by the source, the objectives of the identication processes can be very dierent. Indeed, the term identication methods encompasses goals as dierent as source localization in wide acoustic spaces, source eld reconstruction on academic or complex geometries, noise source identication in enclosed sound elds, source separation or source ranking. A unique method to treat all these applications does not exist and several specialized methods have been developed since decades [START_REF] Magalhães | Sound sources reconstruction techniques : a review of their evolution and new trends[END_REF]. The Near-eld Acoustic Holography (NAH) is one of the most widely used method. Firstly introduced by Williams et al [START_REF] Williams | Sound source reconstructions using microphone array[END_REF] and Maynard et al [START_REF] Maynard | Neareld acoustic holography : I. Theory of generalized holography and the development of NAH[END_REF], this Fourier approach reconstructs the sound sources using a spatial 2D Fourier transform of the Helmholtz equation. Initially, NAH was developped to be used with a planar array of microphones but has been extended to cylindrical [START_REF] Williams | Generalized neareld acoustical holography for cylindrical geometry : Theory and experiment[END_REF] or spherical geometries [START_REF] Jacobsen | Near eld acoustic holography with microphones on a rigid sphere[END_REF]. Many alternative methods, all based on NAH, can be found in the literature. One can cite, for example, the works done by Chardon et al [START_REF] Chardon | Near-eld acoustic holography using sparse regularization and compressive sampling principles[END_REF] who use the property of sparsity of signals in NAH. NAH can also be coupled to force identication method [START_REF] Pézerat | Identication of vibration excitations from acoustic measurements using near eld acoustic holography and the force analysis technique[END_REF] to provide a complete pipeline from radiated noise to structural excitation. A second category of identication methods is based on equivalent sources. In this type of methods, the source (eventually irregularly shaped) are replaced by a set of monopoles and dipoles [START_REF] Koopmann | A method for computing acoustic elds based on the principle of wave superposition[END_REF][START_REF] Pereira | Acoustic imaging in enclosed spaces[END_REF] and the unknowns are the strengths of equivalent sources. A third category is based on solving the Helmholtz-Kirchho equation by Boundary Element Method (BEM) [START_REF] Bai | Application of BEM (boundary element method)-based acoustic holography to radiation analysis of sources with arbitrarily shaped geometries[END_REF][START_REF] Kim | On the reconstruction of the vibro-acoustic eld over the surface enclosing an interior space using the boundary element method[END_REF][START_REF] Langrenne | Boundary element method for the acoustic characterization of a machine in bounded noisy environment[END_REF] or Finite Element Method (FEM) [START_REF] Aucejo | Identication of source velocities on 3D structures in non-anechoic environments : Theoretical background and experimental validation of the inverse patch transfer functions method[END_REF][START_REF] Kletschkowski | Identication of noise source in an aircraft fuselage using an inverse method based on a nite element model[END_REF]. The identication method used in the present article belongs to the latter. The inverse Patch Transfer Functions (iPTF) method, rstly introduced by Aucejo [START_REF] Aucejo | Identication of source velocities on 3D structures in non-anechoic environments : Theoretical background and experimental validation of the inverse patch transfer functions method[END_REF] and Vigoureux [START_REF] Vigoureux | Déconnement de sources acoustiques par utilisation d'une méthode holographique à double information[END_REF] is aimed to reconstruct the velocity eld of a radiating structure using acoustic measurements. The basic idea of the iPTF method is based on the use of the Green identity on a nite acoustic volume. It can be considered as a data completion method following the denition given by Langrenne and Garcia [START_REF] Langrenne | Data completion method for the characterization of sound sources[END_REF] : "the acoustic gradients and pressures are known on a surface surrounding the source, but are unknown on its structure".

In addition, iPTF can also be considered as a separation method [START_REF] Fernandez-Grande | Sound eld separation with sound pressure and particle velocity measurements[END_REF][START_REF] Braikia | Evaluation of a separation method for source identication in small spaces[END_REF] as it permits to identify a source in small spaces or in non-anechoic acoustic environment. The iPTF method has been applied with success on baed pistons with the presence of external correlated sources [START_REF] Aucejo | Identication of source velocities in presence of external correlated sources with the inverse Patch Transfer Functions (iPTF) method[END_REF], on L-shaped plate [START_REF] Aucejo | Identication of source velocities with Inverse Patch Transfer Functions method[END_REF] and on irregularly shaped structures such as a guitar [START_REF] Vigoureux | Mesure de vitesse vibratoire par la méthode iPTF : simulation et expérimentation[END_REF], a precatalyst of a diesel engine [START_REF] Totaro | Identify velocity of a complex source with iPTF method[END_REF] or a reduction gearbox [START_REF] Totaro | Source identication on a reduction gearbox using acoustical measurements in a non-anechoic environment[END_REF]. However, as all inverse methods [START_REF] Hansen | Rank-decient and discrete ill-posed problems[END_REF][START_REF] Williams | Regularization methods for near-eld acoustical holography[END_REF], the iPTF suers from problems associated with the inversion of an ill-posed problem. The objective of the present article is to understand the physical origin of that ill-conditioning by deeply analysing the Tikhonov regularization process, one of the most widely used in the literature. As usual, the choice of the regularization parameter is crucial. It is demonstrated that, at some frequencies, the well-known L-curve strategy fails. This is due to the fact that L-curve exhibits more than one corner. The presence of these multiple corners is explained and a modied Lcurve strategy, called weighted curvature strategy, is proposed and compared to results obtained with Generalized Cross Validation (GCV) technique. Two experiments are used to illustrate the regularization process and to evaluate capabilities of the iPTF method : a at plate and an oil pan.

Inverse Patch Transfer Function (iPTF) method

In the example of Fig. 1, a vibrating irregularly shaped surface radiates in any acoustic environment. The objective of the iPTF method is to separate the sound eld radiated by the source and the sound eld coming from the acoustic environment (disturbing sources, reections, etc.). This is achieved by dening a nite acoustic volume Ω surrounding this vibrating surface as presented in Fig. 1. This virtual volume is delimited by surfaces Σ (surface of the vibrating object), Γ (an arbitrarily chosen virtual surface surrounding the object) and Λ (a physically rigid surface). In that case, the Green's identity writes

Ω (Ψ∆Φ -Φ∆Ψ)dΩ = Σ Φ ∂Ψ ∂n -Ψ ∂Φ ∂n dΣ + Γ Φ ∂Ψ ∂n -Ψ ∂Φ ∂n dΓ + Λ Φ ∂Ψ ∂n -Ψ ∂Φ ∂n dΛ. (1)
In Eq. ( 1), the functions Φ and Ψ can be arbitrarily chosen (provided that they are dened and twice dierentiable on the domain). In the following, the function Φ will be the acoustic pressure p dened by the Helmholtz equation (considering there is no source in the volume)

∆p(M ) + k 2 p(M ) = 0, ∀M ∈ Ω, (2)
and boundary conditions on surfaces Σ, Γ and Λ are given by :

∂p(N ) ∂ n = 0, ∀N ∈ Λ, (3)
∂p(Q) ∂ n = -iωρ 0 V (Q), ∀Q ∈ Γ (4)
and

∂p(Q ) ∂ n = -iωρ 0 V (Q ), ∀Q ∈ Σ (5)
where V (Q) is the normal velocity at each point Q of the surface and n is the outer normal of the volume Ω, ρ 0 is the density of air, ω is the angular frequency and k the acoustic wavenumber. The function Ψ of the Green's identity is chosen to be the eigenmodes φ n of the virtual acoustic volume respecting Helmholtz equation :

∆φ n (M ) + k 2 n φ n (M ) = 0, ∀M ∈ Ω, (6)
where k n is the natural wavenumber of mode n. The boundary conditions on surfaces Σ, Γ and Λ are given by :

∂φ n (N ) ∂n = 0, ∀N ∈ Λ, (7)
∂φ n (Q) ∂n = 0, ∀Q ∈ Γ (8)
and

∂φ n (Q ) ∂n = 0, ∀Q ∈ Σ. (9)
These chosen eigenmodes satisfy the homogeneous Neumann's problem. This is an arbitrary choice and does not correspond to the boundary conditions of the real problem. The eigenmodes are only used here as a mathematical basis of functions to solve the homogeneous Helmholtz equation. This is why the acoustic volume delimited by surfaces Γ, Σ and Λ is called virtual volume.

It is not a closed volume but just a part of the acoustic environment of the source.

The velocity of the virtual surface will represent the presence of the acoustic environment outside the virtual volume even if it is reverberant, geometrically complex or even if it contains disturbing sources (sources outside the volume of interest).

Solving the Helmholtz Eq. ( 2) using decomposition on modes φ n (M ), one can express the pressure at any point of the volume Ω as a function of normal velocities of surfaces Γ and Σ

p(M ) = iωρ 0 ∞ n=1 φ n (M ) χ n (k 2 -k 2 n ) Γ φ n (Q)V (Q)dQ + +iωρ 0 ∞ n=1 φ n (M ) χ n (k 2 -k 2 n ) Σ φ n (Q )V (Q )dQ , (10)
where χ n is the norm of mode n.

Due to the boundary conditions given by Eqs. ( 3) and ( 7), the presence of the physically rigid surface Λ does not appear in the expression of the pressure. Indeed, the presence of the rigid wall is already taken into account in the modal basis of modes φ n (M ). This can be an additional advantage of using eigenmodes of the virtual volume : this volume can be inhomogeneous and reections on some rigid structures (a heavy frame for example) placed inside the virtual volume are taken into account by eigenmodes.

Dividing the surfaces Γ and Σ into elementary surfaces -called patches Γ i (or Γ j ) and Σ k (or Σ l ) -the integrals are replaced by sums in Eq. [START_REF] Pereira | Acoustic imaging in enclosed spaces[END_REF]. The averaged pressure on one patch of the surface Γ can be expressed as a function of velocities of patches Γ j and Σ k (of areas S j and S k )

p i = iωρ 0 N j j=1 V j ∞ n=1 S j φ n i φ n j χ n (k 2 -k 2 n ) + iωρ 0 N k k=1 V k ∞ n=1 S k φ n i φ n k χ n (k 2 -k 2 n ) , (11)
where • j = 1 S j S j •dS represents the space average on a patch j. The pressures on patches j of the surface Γ can be written in a matrix form

p i = Z ij V j + Z ik V k ( 12
)
where Z ij (resp. Z ik ) is an impedance matrix where each term is dened by the ratio between the mean pressure on patch i and the mean normal velocity on patch j (resp. k)

Z ij = p i V j . ( 13
)
As a consequence by inverting Eq. ( 12) one can express the source velocities V k (on patches Σ k ) as a function of pressures p i and velocities V j measured on the virtual surface Γ.

Finally, the source velocities are obtained solving

V k = Z -1 ik (p i -Z ij V j ) . (14)
Even if the source velocity is the primary unknown, it is possible with no additional inversion of matrix, to compute the boundary pressure p k on source surface. Indeed, the pressure on the surface of the source is due to the contribution of the velocity of the source V l (identied during the previous step) and to the contribution of the particle velocity V j of the virtual surface (already measured) :

p k = Z kj V j + Z kl V l . (15)
The normal intensity is simply deduced from :

I k = 1 2 (p k • V * k ) , (16)
where (•) and • * denote the real part and the complex conjugate of a complex number. • represents the Hadamard's product. Thus, the acoustic source is completely characterized by the identication process. At the end, global indicators like active and reactive acoustic power and radiation eciency of the source can also be produced with no eort (neither numerically nor experimentally).

3. Classical Tikhonov's regularization strategies

Ill-posed and underdetermined problems

The primary unknown of the problem is the source velocity on surface Γ. To obtain this source velocity, the inversion of an ill-conditioned matrix is needed. Eq. ( 14) has the classical form Ax = b where pi = p i -Z ij V j . The impedance matrix Z ik is computed using modes of the virtual volume with rigid walls extracted by standard nite element methods. This matrix represents the model used to represent the system. The vector pi is measured with pU probes [27] on the virtual surface Γ. This vector represents the output of the system. This vector is obviously aected by measurement noise. The objective here is to deduce the input of the system (the source velocities) knowing the model and measuring the output of the system. This problem is a classical ill-posed problem as dened by Hansen [START_REF] Hansen | Rank-decient and discrete ill-posed problems[END_REF]. If N k (the number of patches on the source surface) is higher than N j (the number of patches on the virtual surface), then the problem is underdetermined and many regularization techniques can be used. As this ill-posed problem is classic in acoustic holography, these techniques have been extensively studied in the literature. The most known method is the Tikhonov regularization but alternative approaches like the Truncated Singular Values Decomposition [START_REF] Hansen | Rank-decient and discrete ill-posed problems[END_REF], the Bayesian regularization [START_REF] Antoni | A Bayesian approach to sound source reconstruction : optimal basis, regularization, and focusing[END_REF][START_REF] Venkatesh | On Tikhonov regularisation[END_REF], the sparse regularization [START_REF] Chardon | Near-eld acoustic holography using sparse regularization and compressive sampling principles[END_REF] or the l 1 -norm minimization [START_REF] Simard | Acoustic source identication : experimenting the l 1 minimization approach[END_REF] also exist. The Tikhonov regularization is based on the minimization of a functional J dened as the sum of the residual norm and a penalty term proportional to the norm of the regularized solution

J = ||Z ik V k -pi || 2 + β||LV k || 2 (18)
where || • || 2 is the l 2 -norm, β is the regularization parameter and L is a non-negative denite matrix often assumed to be the identity matrix. The crucial point of the Tikhonov regularization is the choice of this regularization parameter. A small β value (less ltering) gives more importance to the minimization of the residual norm so that the solution ts well with measurements. However, the norm of this solution is generally high in that case. Conversely, a high β value (more ltering) gives more importance to the minimization of the norm of the solution but the solution deviates from the measurements and the residual norm increases. The good choice of β is a balance between these two opposite behaviors.

Using the Singular Values Decomposition (SVD) of the matrix Z ik , the regularized solution V R k that minimizes the functional J is given by

V R k = N V S m=1 σ 2 m σ 2 m + β u H m pi σ m v m , (19)
where σ m are singular values in a decreasing order, u m and v m are the left and right singular vectors of matrix Z ik and N V S is the number of singular values (equal to N j the number of patches on the virtual surface).

In fact, the Tikhonov regularization acts like a low pass lter on singular values given by lter factor

f m = σ 2 m σ 2 m + β . ( 20
)
For TSVD, the lter factor is step-like (1 or 0). The Tikhonov lter factor is smoother and continuously weights the singular values as a function of β parameter. The choice of β is, in turn, the choice of a threshold on the acceptable singular values [START_REF] Gauthier | Acoustical inverse problems regularization : Direct denition of lter factors using signal-tonoise Ratio[END_REF].

The search of the optimal regularization parameter β for Tikhonov regularization has produced a lot of interesting studies in the framework of acoustical inverse problems. For instance, the performances of the L-curve and the Generalized Cross Validation methods (two known regularization parameterdetermination methods) have been studied by Williams [START_REF] Williams | Regularization methods for near-eld acoustical holography[END_REF], Leclère [START_REF] Leclère | Acoustic imaging using under-determined inverse approaches : Frequency limitations and optimal regularization[END_REF] and Kim and Nelson [START_REF] Kim | Spatial resolution limits for the reconstruction of acoustic source strength by inverse methods[END_REF].

Curvature of the L-curve

The L-curve is a widely used technique to select the optimal regularization parameter β. The L-curve plots the evolution of the two terms of Eq. ( 18) : the norm of the regularized solution as a function of the residual norm. The compromise between both terms is located at the corner of the curve which is generally L-shaped. To easily nd the β value corresponding to the corner of the L-curve, the curvature function is generally used [START_REF] Gauthier | Acoustical inverse problems regularization : Direct denition of lter factors using signal-tonoise Ratio[END_REF][START_REF] Leclère | Acoustic imaging using under-determined inverse approaches : Frequency limitations and optimal regularization[END_REF][START_REF] Kim | Optimal regularisation for acoustic source reconstruction by inverse methods[END_REF][START_REF] Choi | Comparison of methods for parameter selection in Tikhonov regularization with application to inverse force determination[END_REF]. A corner of the L-curve corresponds to a maximum value of the curvature function.

Generalized Cross Validation

The Generalized Cross Validation [START_REF] Gauthier | Acoustical inverse problems regularization : Direct denition of lter factors using signal-tonoise Ratio[END_REF][START_REF] Leclère | Acoustic imaging using under-determined inverse approaches : Frequency limitations and optimal regularization[END_REF][START_REF] Kim | Optimal regularisation for acoustic source reconstruction by inverse methods[END_REF][START_REF] Choi | Comparison of methods for parameter selection in Tikhonov regularization with application to inverse force determination[END_REF] is an alternative to the L-curve to determine the optimal parameter β. This technique is based on the minimization of the GCV function given by

GCV (β) = ||I -B(β) pi || 2 2 [T r(I -B(β))] 2 (21)
where the inuence matrix is

B(β) = Z ik (Z H ik Z ik + βI) -1 Z H ik 4.
2D underdetermined problem and regularization process

System under study

To evaluate the ability of iPTF method in the reconstruction of source eld, two experimental validations have been set up : one 2D application on an academic structure and one 3D application on an industrial structure. The rst system under study consists in an aluminium rectangular plate glued to a wood frame as can be seen in Fig. 2. The free surface of the plate is 60 cm long and 30 cm large. The plate is excited by an electro-dynamic shaker fed with a white noise between 100 and 1000 Hz and radiates in a quiet room. The measurements have thus not been done in free eld conditions. However, the iPTF method has been developed to isolate the object from its acoustic environment as previously explained. The iPTF method is considered as a eld separation technique thanks to the use of the double information pressure/particle velocity on a surface surrounding the source. A pU probe has been used to measure the acoustic pressure and the particle velocity on : a 17 × 10 grid mesh (see Fig. 3(a)) to have reference elds of source velocity, boundary pressure and normal intensity. These reference measurements will be used to evaluate the iPTF capabilities. a virtual surface (a 5 cm high rectangular box) presented in Fig. 3(b) divided into 535 patches. The acoustic pressure and the particule velocity have been measured at center of each patch. These measurements are used to compute vector pi in Eq. (17). As the measurements on the grid meshes are sequential, a phase reference (an accelerometer) has been glued on the surface of the plate. The iPTF method is applied here to reconstruct the source elds on a 45×25 identication grid mesh, see Fig. 3(c). The problem is then highly underdetermined as there are 1125 patches of identication (on the surface of the plate) and only 535 patches of measurements (on the virtual surface). To compute the matrix Z ik in Eq. ( 17) a nite element model, shown in Fig. 3(d), has been used to extract 1883 modes of the virtual volume up to 11kHz (respecting a λ/6 criterion for the FE mesh size). The condition number of the matrix Z ik is plotted in Fig. 4 as a function of frequency and some additional information about numerical computations are given in Tab. 1. It is interesting to notice that the condition number exhibits peaks in the whole frequency range. As presented in Fig. 4 exactly correspond to the eigenfrequencies of the virtual volume (plotted in dashed lines). This particularity is discussed in the following.

Identication results

In Fig. 5, the identied mean square velocity, mean square pressure and acoustic power are compared to the reference measurements. The measures on the virtual surface are also plotted to show how dierent are the radiated and the reference elds.

Both the Curv (maximum of the curvature of the L-curve) and the gcv strategies have been applied. Both regularization parameter-determination methods compare well in the whole frequency range except at two particular frequencies. Indeed, at 250Hz and 880Hz, the results obtained using Curvstrategy seem to deviate from reference measurements as shown in Fig. 5(b).

To understand why the Curv-strategy fails to identify the good solution for some frequencies, an insight into regularization process has been carried out at two frequencies : at 137.5Hz (where both Curv-strategy and gcv-strategy give satisfactory results) and at 250Hz (where Curv-strategy fails). Fig. 6 shows L-curve, curvature function, GCV function and power ratio Π(β) as a function of regularization parameter β. The power ratio is the ratio between the acoustic power Π id Σ (β) identied on the surface Σ of the vibrating object for each β value and the acoustic power Π meas Γ measured on the virtual surface Γ as dened in Eq. [START_REF] Vigoureux | Mesure de vitesse vibratoire par la méthode iPTF : simulation et expérimentation[END_REF].

Π(β) = Π id Σ (β) Π meas Γ ( 22
)
The power ratio is obviously independent of the reference acoustic power and is obtained using quantities needed (pressure and particule velocity on the virtual surface) or produced (pressure and velocity on the identication surface) by iPTF method. If the power conservation is achieved, the power ratio tends to unity.

At 137.5Hz, both Curv and gcv strategies produce good results. At this frequency, the L-curve is indeed L-shaped. The optimal β values (represented by a cross for Curv-strategy and by a square for gcv-strategy) correspond to the (unique) corner of the L-curve even if they are a little bit dierent (β Curv = 712.2 and β GCV = 2270). As a consequence, β GCV does not correspond exactly to the maximum of curvature of the L-curve and β Curv does not correspond exactly to the minimum of the GCV function. It can be noticed that both strategies imply a power ratio close to 1. Fig. 7 demonstrates and conrms that both strategies produce satisfactory identied elds compared to reference. Precise details can be observed with both strategies.

At 250Hz, the Curv-strategy fails as can be seen in Fig. 5(b). This is mainly due to the non-uniqueness of the corner of the L-curve. Indeed, in that case, Fig. 6 shows that the L-curve is rather a W-shaped curve and exhibits two corners (multiple corners of L-curve have been also observed by Pereira [START_REF] Pereira | A theoretical and experimental comparison of the equivalent source method and a Bayesian approach to noise source identication[END_REF][START_REF] Pereira | Acoustic imaging in enclosed spaces[END_REF] using Equivalent Source Method). These two corners result in two curvature maxima. As can be seen in Fig. 6(d), the Curv-strategy chooses the highest curvature value that corresponds here to the second peak where β Curv = 1.416 × 10 8 . For the rst peak of curvature function (which is not automatically chosen by Curv-strategy), the β value is 179.5. At this frequency, the β Curv value is completely dierent from the one dened by gcv-strategy (β GCV = 600) which is closer to the value of the rst peak of curvature function. It is interesting to notice that it is dicult to distinguish two clearly dierent possibilities for β GCV as with curvature function. One can at most distinguish inection of the GCV curve around β Curv = 1.416 × 10 8 . gcvstrategy seems to be less aected than Curv-strategy by the high condition number of the matrix to be inverted at this frequency. An important remark raises from the analysis of the power ratio in Fig. 6(h). For β GCV , the power ratio is still close to 1 but for β Curv , the solution identied has a power ratio less than 0.3 indicating that the power conservation is not achieved in that case. This is due to an excessive regularization (Eq. 18). On the contrary, the power conservation is achieved for the rst peak of the curvature function. Analysing maps obtained with this β value (noted wCurv in Fig. 6), one can see on Fig. 8 that this solution produces more reliable results.

The wCurv-strategy

To overcome this issue, it is possible to weight the curvature function by the power ratio function to lter the solution for which the power conservation is not achieved. The wCurv-strategy (weighted curvature strategy) only consists in nding the maximum of function wCurv(β) dened by

wCurv(β) = Curv(β)Π(β) if Π(β) > 0 0 if Π(β) < 0 (23)
where Curv(β) is the curvature function and Π(β) is the power ratio function given by Eq. ( 22). The wCurv function is plotted in Fig. 6(d) (in dashed line). The optimal β wCurv values obtained with wCurv-strategy are represented by triangles in Fig. 6. Fig. 9 presents the eect of the weighting on the curvature function as a function of frequency and β value. Fig. 9(a) shows that wrong solutions may appear at frequencies corresponding to the eigenfrequencies of the virtual volume with rigid walls (the associate problem) (vertical dashed lines). The weighting lters out second peak of curvature that produces wrong solutions.

wCurv-strategy takes advantage of information available during the iPTF identication process.

The iPTF method can be experimentally constraining because it needs pressure and particle velocity on a virtual surface surrounding the source, but it provides a complete characterization of the source from the velocity eld to the acoustic power, the active and reactive intensity and also the radiation eciency as presented in Fig. 10. The very low radiation eciency of the plate in this frequency band is validated by the Crocker and Price's model [START_REF] Crocker | Sound transmission using Statistical Energy Analysis[END_REF].

4.4. Interpretation of wrong solutions of Curv-strategy wCurv-strategy permits to overcome wrong choices of Curv-strategy by promoting solutions that respects the power conservation. However, this does not explain why such solutions appear with Curv-strategy.

Analysing the lter factors (Eq. ( 20)) obtained with β dened by the rst and the second peak of curvature observed at 250Hz, one can see that Curvstrategy (β Curv = 1.416 × 10 8 ) strongly lter singular values. Indeed, with such a regularization parameter, only two singular values have a lter factor higher than 0.1 which means that almost all the singular values are ltered as shown in Fig. 11. On the contrary, with β value obtained with wCurvstrategy (β wCurv = 179.5), a lot of singular values participate (187 singular values have a lter factor higher than 0.1).

As a consequence, the wrong solution of Curv-strategy is highly dominated by the rst singular value of the matrix. Comparing this wrong solution to the projection of the second mode (eigenfrequency = 250 Hz) of the virtual volume on the surface of the plate, it is easy to conclude that the rst singular value of the matrix is highly dominated by the eigenmode of the virtual cavity, see Fig. 12. The same conclusion holds at 880 Hz. The identied velocity eld highly correlates to the pressure mode shape of the eighth modes (eigenfrequency = 873.5 Hz).

In the case of single layer pressure-velocity (as it is the case for iPTF),

Fernandez-Grande et al. [START_REF] Fernandez-Grande | Sound eld separation with sound pressure and particle velocity measurements[END_REF] propose to use a weighted least squares solution.

This is necessary beacause the solution vector is, in that case, composed by pressures and velocities

p h u h (24)
with completely dierent amplitudes. The weighting is supposed to reduce the problem, so that the energies of the pressure and velocity observations are similar. This is not the case in iPTF method. Indeed, the solution vector is equivalent to a pressure vector even if it is the sum of contributions of measured pressure and velocity as shown in Eq. ( 17). Thus, the presence of multiple corners in L-curve is not the consequence of the combination of pressure and velocity measurements. In addition, Pereira et al. [START_REF] Pereira | A theoretical and experimental comparison of the equivalent source method and a Bayesian approach to noise source identication[END_REF] have observed the same phenomenon using only pressure measurements. As the multiple corners problem appears only at eigenfrequencies of the virtual volume, we suspect that the wrong solutions found by Curv-strategy correspond to the presence of solutions of the associate problem that is to say to the eigenmodes of the virtual volume.

Application to a 3D underdetermined problem : the oil pan

One of the main advantages of iPTF method is its ability to handle complex 3D source surfaces. To illustrate that crucial point, the iPTF method has been applied to identify the velocity, pressure and intensity elds of an oil pan presented in Fig. 13 Table 2: Information on the FE mesh of the virtual volume used to compute impedance matrix Z ik in the case of the oil pan.

a virtual rectangular box shown in Fig. 14(b). The reference measurements have been acquired on the smallest box surrounding the oil pan, see Fig. 14(a). The reference surface does not correspond to the exact geometry of the source because the measurements would have been to complicated. The identication surface (see Fig. 14(c)) follows the exact geometry of the source thanks to the use of FE mesh of the virtual volume presented in Fig. 14(d).

The oil pan was excited by an electro-dynamic shaker in a frequency band from 100 to 5000 Hz. The acoustic pressure and the particle velocity have been measured at the center of each patch of (i) the reference surface and (ii) the virtual surface. An accelerometer glued on the surface of the oil pan has been used as a phase reference. The impedance matrices have been computed up to 10kHz. The number of modes used to compute the impedance matrix is enough to ensure that the rank of the matrix is driven by the number of patches on the virtual surface and not by the number of modes. The condition number of this matrix is plotted in Fig. 15. As previously remarked, the condition number exhibits peaks at eigenfrequencies of the virtual volume. Some additional information can be found in Tab. 2.

The mean square pressures are plotted in Fig. 16. As can be seen, results obtained with wCurv and GCV strategies are in good agreement with reference measurement. Again, as presented in Fig. 17, for some frequencies, the L-curve is not Lshaped. This is the case, for example, at 432Hz (Fig. 17(a)). The Curv chooses the β value that gives the maximum of curvature but by doing this chooses the wrong solution at this frequency. Again, this solution does not respect the power conservation (power ratio close to 0). Using wCurv-strategy, this solution is weighted by the power ratio and so is highly soften. The rst peak which respects the power conservation is chosen by the wCurv-strategy.

At 904Hz, this problem does not appear and the solutions found by Curv and wCurv-strategy are equal.

The velocity elds obtained by both strategies at these two frequencies are presented in Fig. 18. When the L-curve is not L-shaped, the Curv-strategy produces results that are highly dominated by mode shapes of the virtual cavity as visible in Fig. 18(b).

The GCV-strategy produces slightly under-regularized solutions in medium and high frequency as presented in Figs. [START_REF] Vigoureux | Déconnement de sources acoustiques par utilisation d'une méthode holographique à double information[END_REF], 17 and 18. The norm of the solution is a little higher than the one obtained by wCurv-strategy. However, this conclusion must be nuanced because of the dierence between identication and reference surfaces. In any case, a combinaison between GCV and wCurv strategies may produce satisfactory results in a lot of dierent applications.

Finally, to demonstrate the accuracy of the elds identied by the iPTF method associated to the wCurv-strategy of regularization, one has plotted the phase of the pressure eld in high frequency (2980Hz). Very ne details of the highly variating eld are reproduced with very good accuracy compared to reference measurements.

Conclusion

Several acoustic source identication methods exist in the literature. Each of them has its particularities, its advantages and its drawbacks. The inverse Patch Transfer Function method has been developped some years ago to identify the velocity, pressure and intensity eld of a 3D vibrating surface by measuring the radiated acoustic elds (pressure and particle velocity) on an arbitrary virtual surface surrounding the source. The solver used here is based on the computation of an impedance matrix from eigenmodes of the virtual volume (dened by the union of the source surface and the virtual surface of measurements). Neumann's boundary conditions have been chosen even if they don't reect any reality. As possible in Green's identity, they have been chosen arbitrarily to solve the problem on a basis of orthonormal functions. However, the matrix to be inverted is ill-conditioned (as in the major part of inverse methods), especially at frequencies corresponding to the eigenfrequencies of the associate problem (eigenmodes of the virtual volume). This results in the appearance of alternative solutions during the Tikhonov re-gularization process. At these frequencies, the strategy based on the choice of the regularization parameter corresponding to the maximum of curvature of the L-curve sometimes fails. A wrong solution highly dominated by the eigenmodes of the virtual volume is produced. It has been demonstrated on two real examples (a at plate and an oil pan) that these solutions don't respect the power conservation. This permits to introduce a weighted curvature function that overcomes the problems observed at frequencies close to eigenfrequencies of the virtual volume. The GCV technique seems to be less aected by the high condition number but can, in some applications, induce a slightly under-regularized solution. A combination of both GCV and wCurv strategies can then be useful.

The iPTF produces then very accurate source elds and gives a lot of information about the source : velocity, wall pressure, active and reactive intensity, acoustic power and radiation eciency.

  by measuring the acoustic radiated elds on Number of measurement patches N j 672 Number of identication patches N k 930 Maximum frequency of the modal basis 10kHz Number of modes 1516 Number of DoF 63224 Rank of the matrix Z ik 672

Figure 1 :

 1 Figure 1: Denition of the nite acoustic volume Ω (virtual volume) and the boundary surfaces Σ (surface of the vibrating object), Γ (virtual surface surrounding the object) and Λ (physically rigid surface).

Figure 2 :Figure 3 :

 23 Figure 2: Photographies of the baed plate under study. The measures of the pressures and the acoustic velocities have been done using a pU probe. The plate is excited by a shaker fed with white noise.

Figure 4 :Figure 5 :

 45 Figure 4: Condition number of the matrix to be inverted. The dashed lines represent the positions of the eigen-frequencies of the virtual cavity.

Figure 6 :

 6 Figure 6: (a)-(b) L-curve, (c)-(d) curvature of the L-curve, (e)-(f) GCV and (g)-(h) power ratio as a function of regularization parameter at frequencies 137.5Hz and 250Hz. Crosses : Curv-strategy ; Triangles : wCurv-strategy ; Squares : GCV-strategy. The dashed lines in (c) and (d) represent the weighted curvature function dened by Eq. (23)

Figure 7 :

 7 Figure 7: Comparison between reference measurements ((a)-(d)) and identication results obtained with GCV ((e)-(h)) and Curv ((i)-(l)) strategies in terms of velocity (V , m.s -1 ), pressure (P , Pa), active intensity (I ac , W.m -2 ) and reactive intensity (I re , W.m -2 ) elds at 137.5Hz.

Figure 8 :

 8 Figure 8: Comparison between reference measurements ((a)-(d)) and identication results obtained with GCV ((e)-(h)), Curv ((i)-(l)) and wCurv ((m)-(p)) strategies in terms of velocity (V , m.s -1 ), pressure (P , Pa), active intensity (I ac , W.m -2 ) and reactive intensity (I re , W.m -2 ) elds at 250Hz.

Figure 9 :

 9 Figure 9: (a) curvature of L-curve and (b) weighted curvature of L-curve as a function of frequency and regularization parameter β.

Figure 10 :Figure 11 :

 1011 Figure 10: Radiation eciency of the plate. Solid black line : identication with wCurvstrategy ; dotted red line : identication with GCV-strategy ; dashed-dotted line : Crocker and Price's model[START_REF] Crocker | Sound transmission using Statistical Energy Analysis[END_REF]

Figure 12 :

 12 Figure 12: Comparison between the identied velocity eld using Curv-strategy and projection of pressure mode shapes of the virtual cavity on the identication surface. (a) identied velocity eld at 250Hz ; (b) projection of the second mode (eigenfrequency = 250Hz) ; (c) identied velocity eld at 880Hz ; (d) projection of the eighth mode (eigenfrequency = 873.5Hz).

Figure 13 :Figure 14 :

 1314 Figure 13: Photographies of the baed oil pan under study. The measures of the pressures and the acoustic velocities have been done using a pU probe. The oil pan is excited by a shaker fed with white noise.

Figure 15 :

 15 Figure 15: Condition number of the matrix to be inverted. The dashed lines represent the positions of the ten rst eigen-frequencies of the virtual cavity.

Figure 16 :

 16 Figure 16: Mean square pressure. Red dashed line : reference measurements, black solid line : identication with Curv-strategy, green dashed-dotted line : identication with GCV-strategy.

  Power ratio 432Hz (h) Power ratio 904Hz

Figure 17 :

 17 Figure 17: (a)-(b) L-curve, (c)-(d) curvature of the L-curve, (e)-(f) GCV and (g)-(h) power ratio as a function of regularization parameter at frequencies 432Hz and 904Hz.Crosses : Curv-strategy ; Triangles : wCurv-strategy ; Squares : GCV-strategy.

Figure 18 :

 18 Figure 18: Amplitude of the source velocity eld (m.s -1 ). Comparison between reference measurements and identication results obtained with Curv, wCurv and GCV strategies at 432Hz (a)-(d) and at 904Hz (e)-(h). Scale of subgure (b) is not the same, it is 0.0109 × 10 -6 to 0.7860 × 10 -6 .

Figure 19 :

 19 Figure 19: Phase (rad) of the pressure eld at 2980Hz. (a) reference measurements ; (b) wCurv strategy.

Table 1 :

 1 Information on the virtual volume used to compute impedance matrix Z ik in the case of the rectangular plate.

	, these peaks

Z ik V k = pi ,(17)

Acknowledgements

This work was performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007).


Dorian Vigoureux

Nicolas Totaro
email: nicolas.totaro@insa-lyon.fr

Jonathan Lagneaux

Jean-Louis Guyader

Inverse Patch Transfer Functions method as a tool for source field identification

Many methods to detect, quantify or reconstruct acoustic sources exist in the literature and are widely used in industry (Near-field acoustic holography, inverse Boundary Element Method, etc.). However, the source identification in a reverberant or non-anechoic environment on an irregularly shaped structure is still an open issue.

In this context, the inverse Patch Transfer Functions Method firstly introduced by Aucejo et al.

[1] can be a suitable method. Indeed, the iPTF method has been developed to identify source velocity on complex geometries and in a non-anechoic environment. However, to obtain good results, the application of the method must follow rigorous criteria that were not fully investigated yet. In addition, as it was firstly defined, the iPTF method only provides source velocity while wall pressure or intensity should also give useful information to engineers.

In the present article, a procedure to identify wall pressure and intensity of the source without any additional measurement is proposed. This procedure only needs simple numerical post-processing. Using this new intensity identification, the influence of background noise, evanescent waves and mesh discretization are illustrated on numerical examples. Finally an experiment on a vibrating plate is shown to illustrate the iPTF procedure.

Introduction

Many methods, based on acoustic measurements, have been developed to detect and quantify noise sources. Those methods are particularly used when a direct measure on an object is not possible (because of the complexity of its geometry for example) or if one wishes to predict the noise produced by the object. Among all these methods allowing to solve the inverse problem in acoustics, the Near-field Acoustical Holography has to be cited as an example, as well as the inverse Boundary Element Method. The Near-field Acoustical Holography (NAH) was developed by J. D. Maynard et al. [START_REF] Maynard | Nearfield acoustic holography: I. theory of the generalized holography and the development of NAH[END_REF]. It permits to identify the velocity field of the source, using the pressure field measured on a hologram located near the source.

Initially, this method could only be applied for the reconstruction of a planar structure and was brought to various academic shapes by E. G. Williams et al. [START_REF] Williams | Generalized nearfield acoustical holography for cylindrical geometry: Theory and experiment[END_REF][START_REF] Williams | Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography[END_REF]. One of its limitations is that it has to be applied in the near-field to prevent effects due to vanishing waves or finite 2D Fourier transform.

The inverse Boundary Elements Method (iBEM) is based on the numerical computation and inversion of transfer matrices. This method is suitable for irregularly shaped sources but suffers from a high computational cost, especially when the complexity of the structure and of the acoustical domain increases.

A hybrid method [START_REF] Wu | Hybrid nearfield acoustic holography[END_REF] has been set up using elements from both the NAH and the iBEM in order to tackle their limitations. This hNAH method is based on the least square solution of the Helmholtz equation to express the radiated sound pressure as a combination of inbound and outbound spherical waves.

The NAH can also be coupled to force distribution identification technique as presented by Pézerat et al. [START_REF] Pézerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF]. The objective is here to use the velocity field identified by NAH method to deduce the force distribution that produces vibration and sound radiation. The source of vibration is then quantified and localized by acoustic measurements. Another method has been recently introduced by Langrenne et al. [START_REF] Langrenne | Boundary element method for the acoustic characterization of a machine in bounded noisy environment[END_REF][START_REF] Langrenne | Measurement of confined acoustic sources using near-field acoustic holography[END_REF]. This method, known as Deconfined Acoustic Holography (DAH) is based on the decomposition of the acoustic field measured on a surface surrounding the source. It states that the measured source is the composition of the free field pressure radiated by the primary source, the incoming field radiated by all secondary sources and reflected by the enclosure and the field scattered by the surface of the tested source. The separation of the outgoing field from the incoming field is performed using a Helmholtz standard integral formulation and the removing of the scattering of the incoming field is possible as the primary source is considered to be a rigid body.

We present in this paper the inverse Patch Transfer Functions method (iPTF), introduced to identify the vibration velocity on irregularly shaped sources [1] in a non-controlled acoustic environment (reverberant room, stationary disturbing sources, etc.) even in the presence of rigid obstacles (non-homogeneous acoustic medium). This inverse formulation is based on the direct formulation firstly developed by Ouisse et al. [9].

In its direct formulation (from the source to the radiated sound pressure), the PTF consists in calculating acoustic impedances of uncoupled sub-domains (cavities, structures or semi-infinite media) to predict the sound pressure or the structural velocity when all the sub-domains are coupled. The whole domain can contain sources, absorbing materials, rigid surfaces, openings, etc. To apply the PTF method, one thus proceeds to the decomposition of the global domain into sub-domains to reduce the computation time. These acoustical sub-domains are then coupled by their common surfaces divided into elementary surfaces called patches. The acoustic impedances are evaluated numerically on these patches using the suitable method (analytic solution, finite element method, Rayleigh's equation, measurements, etc.). Using the coupling conditions on the coupling surfaces, one obtains a system of linear equations where unknowns are the coupling velocities.

Solving this system of equations leads to the possibility of computing the response at any point of the whole domain. This method, developped in the framework of the VisPeR and SILENCE european projects, has been successfully applied on industrial applications [START_REF] Ouisse | Application of a substructuring approach for linear acoustics in automotive context[END_REF]11] and can be a useful tool to detect efficient positions of absorbing materials to reduce noise in the domain [12].

In the direct formulation, the coupling velocity is thus deduced from a system of equations when the source is known (monopoles, velocity field for example). In the inverse formulation presented here, the source is unknown and the coupling velocity is measured on an arbitrarily defined surface surrounding the source. Using this method makes possible to identify the velocity field of the source even if the considered surface is unreachable thanks to the finite element solver used here [START_REF] Vigoureux | Identification des champs de vitesse, pression et intentié à l'aide de la méthode iPTF[END_REF].

In this article, the direct formulation will be associated to the inverse formulation to give access to information about boundary pressure and intensity on the source surface. Then, the limitations of the method will be investigated. A particular attention will be paid to the distance to the object as the radiated field partly consists of evanescent waves that cannot be measured beyond the near-field. A notion of penetration depth to estimate the bandwidth limitation of the iPTF method will be introduced.

In addition, the influence of the discretization of the identification and the measurement surfaces will be investigated and rigorous criteria will be defined. The sensitivity of iPTF to noise measurement will also be evaluated.

Finally, an experiment on the vibration of a flat plate excited by a shaker will be shown to illustrate the iPTF procedure.

Theoretical background

The method used hereby combines the direct and inverse formulation of the PTF method. We will then remind the equations used in both cases, starting with the direct formulation.

Solving the direct problem

As a general matter, the problem here is to compute the pressure and velocity fields at any point in a cavity due to one or more sources. Such a cavity is illustrated in Fig. 1. The source is represented by the surface Σ vibrating with a known velocity field. In the example of Fig. 1, Γ is a rigid surface and Σ is an opening towards an arbitrary acoustic environment. In this cavity, the acoustical problem can be solved using the Helmholtz equation

∆p(N) + k 2 p(N) = 0 ∀N ∈ A, (1)
where p(N) is the acoustic pressure at point N and k = ω/c is the wavenumber with ω the angular frequency and c the sound speed in the acoustic medium. The boundary conditions are defined in Eqs. ( 2) and ( 3)

∂p(Q) ∂n = 0 ∀Q ∈ Γ (2)
and

∂p(Q) ∂n = -iωρ 0 V (Q) ∀Q ∈ Σ, Σ , (3)
where ρ 0 is the density of the acoustic medium and V (Q) is the normal velocity on surfaces Σ and Σ . To solve Eq. (1) one can use Green's theorem:

A [φ n (N)∆p(N) -p(N)∆φ n (N)] dV = Σ∪Γ φ n (N) ∂p(N) ∂n -p(N) ∂φ n (N) ∂n dS, (4)
where φ n (N) could be an arbitrary function (defined and twice differentiable on the domain). To simplify Eq. ( 4), it has been chosen to use the eigen-modes of the cavity with all boundary surfaces (Γ, Σ and Σ ) considered to be rigid. This decomposition is only a mathematical tool used to solve the problem even if this is not in accordance with the physical boundary conditions. Such a decomposition is valid in PTF approach as it is a substructuring approach where sub-domains have to be solved independently from each other. In that case, for example, the opening is not represented by the modal basis of the cavity but by coupling conditions between sub-domains. An illustration of that method on two coupled cavities is shown in [9].

As φ n (N) are eigen-modes of the cavity with rigid walls, they respect Helmholtz equation

∆φ n (N) + k 2 n φ n (N) = 0 ∀N ∈ A, (5)
∂φ n (Q) ∂n = 0 ∀Q ∈ Σ ∪ Σ ∪ Γ. (6)
Writing p(N) as a sum of contributions of eigen-modes and considering the boundary conditions (Eq. ( 2) and ( 3)), one can express the pressure at any point in the cavity:

p(N) = iωρ 0 ∞ ∑ p=1 φ p (N) Λ p k * 2 -k 2 p Σ φ p (Q)V (Q)dQ + Σ φ p (Q )V (Q )dQ (7)
where Λ p is the norm of the p-th mode, k * = k(1 + iη), k is the acoustic wavenumber and η is the damping loss factor of the fluid. As can be seen in Eq. ( 7), the surface Γ, which is physically rigid, does not appear in the calculation of the pressure at point N. This is due to the rigid wall boundary conditions we used to express the pressure that actually corresponds to real boundary conditions.

To compute Eq. ( 7), it is necessary to express it in a discrete form

p(N) = iωρ 0 ∞ ∑ p=1 φ p (N) Λ p k 2 -k 2 p N Σe ∑ Σ e =1 φ p Σ e V Σ e A Σ e + N Σ e ∑ Σ e =1 φ p Σ e V Σ e A Σ e     (8)
where χ S is the mean of the variable χ on surface S, and A S is the area of surface S. Eq. ( 8) can be reorganized as a function of velocity information. Indeed, if one isolates the terms depending on the velocity from the rest of the equation, one can write

p(N) = N Σe ∑ Σ e =1 V Σ e • Z Σ e →N + N Σ e ∑ Σ e =1 V Σ e • Z Σ e →N (9)
where

Z Σ e →N = iωρ 0 ∞ ∑ p=1 A Σ e φ p (N) φ p Σ e Λ p k * 2 -k 2 p . (10)
The term Z Σ e →N is called acoustic impedance. It permits to calculate the pressure at a point N using the velocities of a series of surfaces Σ e called patches. The notation Σ e → N indicates that we consider the transmission path from the patch Σ e to the point N. Eq. (9) shows that the radiated pressure can be computed using the velocity of the vibrating surface Σ and the velocity of the opening surface Σ . To solve the direct problem, the velocity on Σ is supposed to be known and the velocity on Σ is deduced from coupling conditions between sub-domains as presented in [9]. The basic idea of the inverse Patch Transfer Function is included in Eq. ( 9). Indeed, if the coupling velocity V Σ e is not unknown but measured all the informations of the outer sub-domain are taken into account and the pressure p(N) can be computed at any point of the sub-domain. If the effective velocity on coupling surface V Σ e is known, the modelization of acoustic environment outside the acoustic volume is not necessary but intrinsically taken into account by measurements.

Solving the inverse problem

Eq. ( 9) permits to compute the pressure at any point of the sub-domain and particularly on the coupling surface Σ .

Averaging this pressure on an elementary surface Σ e , one can express the following system of equation

P Σ e = Z Σ e →Σ e V Σ e + Z Σ e →Σ e V Σ e (11)
where the patch to patch impedance is defined as

Z Σ e →Σ e = iωρ 0 ∞ ∑ p=1 A Σ e φ p Σ e φ p Σ e Λ p k * 2 -k 2 p ( 12
)
The general purpose of the inverse problem is to compute the velocity field of the source based on measurements of the radiated field. Starting from Eq. ( 11), the velocities of elementary surfaces Σ e of the source can be expressed as a function of acoustic patch impedance matrices and patch pressure and patch acoustic velocity on elementary surfaces Σ e of the opening:

V Σ e = Z Σ e →Σ e -1 P Σ e -Z Σ e →Σ e V Σ e (13)
The acoustic patch impedance matrices of the volume defined by surfaces Σ and Σ are computed using Eq. ( 12) and eigen-modes extracted with standard finite element solver. The pressure and velocity vectors of Eq. ( 13) are measured on the virtual surface Σ using pU probe for example.

The geometry of the acoustic volume is defined by the shape of the object under study (Σ), by the presence of rigid walls (Γ) and by the measurement surface (Σ ). If no rigid wall physically exists, the measurement surface Σ has to completely surround the source. In that case, Σ can be chosen arbitrarily and the acoustic volume is virtual (not defined by a physical surface). In the following, it is assumed that the exact geometry of the vibrating object is known and that the surface Σ corresponds to it. However, there is no theoretical need for this and the surface Σ can be defined arbitrarily as, for example, a sphere or an ellipsoid around the source. Surface Σ is a surface of identification of particle velocity in the acoustic medium.

If the surface of identification corresponds to the surface of the source, one considers, thanks to continuity conditions, that the particle velocity is equal to the source velocity.

Eq. 13 has the classical form of an ill-posed problem. The solution is found using Tikhonov regularization and the maximum of curvature of the L-curve [14,[START_REF] Leclère | Acoustic imaging using under-determined inverse approaches: Frequency limitations and optimal regularization[END_REF].

Using this identified velocity field, the identification of pressure and intensity fields on the source surface is also possible.

This identification can be performed by combining both direct and inverse formulations of the PTF method.

Starting from the pressure and velocity fields measured on the opening surface Σ , the velocity field of the source can be computed using Eq. ( 13). At this step, the velocity field of the source is identified as already presented in [1].

To go further we combine the identified velocity field to Eq. ( 9). Indeed, as Eq. ( 9) is true for any point N, one can now compute the pressure on the source surface, using the velocity V Σ e measured and the velocity V Σ e computed during step 1. The wall pressure of the source can then be written

P Σ e = Z Σ e →Σ e V Σ e + Z Σ e →Σ e V Σ e . (14)
Finally, the intensity field is easily obtained using the usual equation

I Σ e = 1 2 Re P Σ e V Σ e * (15)
where P Σ e and V Σ e are the pressure and velocity fields computed during steps 1 and 2 and * represents the transpose complex conjugate.

With this new step, the source is completely characterized (velocity, wall pressure, intensity, power, radiation efficiency).

Obviously, the acoustic particle velocity is the primary unknown and its bad estimation will lead to wrong estimations of the other quantities.

Numerical validation of the intensity identification

For this numerical validation, we present the identification performed on a baffled steel plate excited by a harmonic point force as presented in Fig. are used to computed the vibro-acoustic response of the plate (427 modes).

Measurement step (numerical experiment)

The measurement step consists in the measurement of both pressure and velocity fields on the opening surface Σ . For this first application, the opening surface is defined as shown on Fig. 3. We are here in a particular case where the rigid wall Γ is reduced to null and the opening is entirely surrounding the source. It is important to notice here that the radiated field obtained at the measurement step has been provided by a numerical experiment performed using a boundary element method in infinite domain avoiding noise associated to real experiments and permitting to clearly state validity domain of the method and define rigorous criteria needed.

Identification step

We focus here on the normal intensity identified on the source surface. Fig. 4 shows two different maps of the intensity field. The first one (Fig. 4(a)) represents the reference field computed using Boundary Element Method while the second one (Fig. 4(b)) represents the intensity identified on the plate using the iPTF algorithm based on the radiated pressure and velocity fields computed during the numerical experiment. These maps show a good agreement between the reference and the identified intensity field. Besides, amplitudes and positions of maximum are well estimated. These two maps are given at a particular frequency. It is then important to look out for the intensity identification in the whole frequency range. The identified averaged intensity Īidenti f ied (ω) is compared to the reference averaged intensity Īre f (ω) as a function of frequency in Fig. 5(a). This figure shows that the intensity is properly identified although there are some discrepancies existing between 100 Hz and 1 800 Hz. This can be seen more precisely on Fig. 5(b) that represents the identification error as defined in Eq. One can see on this figure that the error is less than 1 dB for frequencies higher than 1800 Hz. Errors in low frequencies are usually due to evanescent waves in the radiated field. We will then study their influence in the identification of the intensity field using the iPTF method.

Effect of the evanescent waves

It is widely known that every structure presents an evanescent radiation below the critical frequency. In this article, we want to study how those evanescent waves are influencing the application of the iPTF method. For this purpose, we will remind here some fundamental equation of plate vibration in order to define a parameter called penetration depth.

The pressure radiated by an infinite plate in a semi-infinite medium is given by [START_REF] Crocker | Handbook of noise and vibration control[END_REF] p(x, y, z)

= -i ω 2 ρ 0 w(x, y) k z • e ik z z (17)
The penetration depth d corresponds here to the distance z at which the radiated pressure amplitude is divided by e = 2.718 compared to the boundary pressure. It is the distance at which the exponential is equal to e -1 and is expressed as the plate thanks to their large wavelengths. Then, the penetration depth quickly decreases as a function of wavelength. Thus, the ability to measure evanescent waves decreases while the frequency increases, in such a way that for ω/ω c > 0.18, being farther than 10 cm makes it difficult to measure the evanescent waves especially if the signal to noise ratio is high. When reaching the critical frequency, the penetration depth slightly increases. Above critical frequency, waves are propagative and the penetration depth is undefined.

d = 1 ω c ω c ω -1 (18)
Conversely, at a given distance z s from the plate, the reduction factor e -z s /d indicates the part of the initial amplitude (at z = 0) that is measured at the distance z s from the plate. This factor is visualized in Fig. 6(b). In this figure, one can identify two different behaviors, below and above critical frequency. Above f c , the reduction factor is close to one that indicates that the amplitude measured at z s is approximately equal to the one at z = 0. Below f c and because of the decrease of evanescent waves, their amplitudes are highly reduced. Obviously, this will have a consequence on the quality of iPTF results mainly because of the measurement noise and the regularization process. In particular, the amplitude of the real velocity of the source might be underestimated in this frequency band.

The agreement of the penetration depth criterion with results presented in Fig. 5(b) is very good. One can see that the more vanishing is the field, the greater is the identification error. As any method based on acoustical measurement, it is proven here that the iPTF method is limited by the knowledge of the evanescent part of the radiated field. Above the critical frequency, the radiated field is essentially propagative, allowing thus an identification with a low error (about 0.5 dB).

These remarks, based on an infinite plate analogy, are guidelines. In the case of a finite plate, each mode has its own radiation efficiency.

Fig. 7(a) and Fig. 7(b) show differences between the reference field and the identified field. As can be seen, both maps are slightly different, the identified one seems to be blurred and less details are visible. As explained previously, this is probably due to the evanescent part of the radiated field. To prove this conjecture, it is interesting to compute the field that could be radiated by these two different fields (reference and identification). 7(a) and 7(b) are highly different. However, these two very different velocity fields produce nevertheless the same radiated velocity at 20 cm, as presented in Fig. 7(c) and 7(d). Contrary, in the near-field (plane at 4 cm), the noticeable discrepancies are due to the influence of evanescent waves. Using the iPTF method, the identification is performed using the radiated field measured on a surface at a certain distance from the source. Due to the existence of the evanescent waves, the identified field can only represent a part of the exact field. That is the part which leads to radiations in the far-field. Even if the identified field shown on Fig. 7(b) is not the exact field that could be measured directly on the plate, it is the part responsible for the radiations that effectively contribute to the sound power in the far-field called supersonic acoustic intensity by Williams [START_REF] Williams | Supersonic acoustic intensity on planar sources[END_REF] The iPTF method is based on a discrete decomposition of the identification and measurement surfaces. Therefore, some errors may occur regarding the discretization used. The results shown here illustrate the influence of the meshes on the identification performed with the iPTF method.

Discretization of the identification surface with bending waves criterion

Fig. 8 shows the different meshes used for the identification. They are respectively made using a λ b /2, λ b /4 and λ b /6 criterion, where λ b is the wavelength of the bending waves in the plate. The measurement mesh is made using the λ a /4 criterion, where λ a is the acoustical wavelength in air. This measurement mesh is considered here as fine enough to ensure a good convergence of results (see §5.2) and evaluate the influence of discretization of identification mesh. identification is performed using the iPTF algorithm based on the simulated pressure and velocity radiated fields. Fig. 9(b) shows clearly that the λ b /2 mesh is quite insufficient to accurately represent the physical phenomena that are taking place.

In this first application, neither the localization is precise nor the amplitude is correctly evaluated, mostly due to the large areas on which the mean amplitude is calculated. The correlation between the λ b /4 identification and the reference is much better. One can easily see, on Fig. 9(c), the nodal lines of the plate. Besides, the values of the velocity are identified with a higher precision. Same as previously, the extrema are lowered by the effect of the mean but as the area of the patches is smaller, this error has a lower effect on the identification of the velocity. Regarding the increase of precision while meshing with a λ b /4 criterion instead of a λ b /2 criterion, one could think that a λ b /6 would lead to better results. This is not the case. The mesh fineness is limited to a certain point beyond which the results become incoherent. Indeed, the iPTF method consists in the inversion of a system whose size is given by the parameters of the problem. The number of unknowns is equal to the number of points where one wishes to compute the velocity, the number of equations is equal to the number of points at which the radiated pressure and velocity have been measured. Thus, the problem can be over or under determined and then treated in least square sense. To verify that the number of independent equations is effectively equal to the number of measurement points, the rank of the matrix to be inverted is computed. Tab. 1 shows the rank of the system for the three meshes. For the first two cases, one can see that the equations are linearly independent whereas the equations of the third system are dependent. The rank of this system is lower than the number of its equations and is not limited by the number of available equations. As said in §2.1, the problem is decomposed on the cavity eigen-modes basis computed up to 8000

Hz. In this frequency band, the number of modes of the virtual cavity is 427. It appears clearly that the rank of the system is limited by the number of modes taken into account in impedances computation. As a consequence, the number of modes must be at least equal to the number of identification points. To verify this hypothesis, we increased the modal basis to modes up to 11 kHz. Doing this, 963 modes are taken into account. The application of the iPTF method on the former data with this extended modal basis leads to the result shown in Fig. (9(d)). One can see on this figure that the vibration field is now identified with a good accuracy.

Discretization of the measurement surface with acoustic wave criterion

Fig. 10 shows the different meshes used for the measurement. They are respectively made using a λ a /2, λ a /4 and λ a /6 criterion, where λ a is the acoustical wavelength in air. The identification mesh shown on Fig. 10(d) is made using the λ b /4 criterion, where λ b is the wavelength of the bending waves in the plate (seen to produce good results in §5.1). Fig. 11 presents the results obtained with these three measurement meshes. see that the three maps show equivalent results. This indicates that the measurement mesh and the number of measurement points are of little importance in the accuracy of the description of the sources and an under-determined problem can be solved here.

Creating a mesh

Considering the remarks presented here about the importance of the measurement and the identification step, one can define some procedure to obtain a correct identification using the iPTF method. Let call f max the maximal frequency of identification. In order to ensure a correct visualization of the physical phenomena, it is important to choose a patch size for the identification mesh that will be small enough to render them with a sufficient precision. Thus, the identification patches should not be larger than the fourth of the bending wavelength in the structure.

Once the size of the identification patches is chosen, and knowing the identification area, it is possible to calculate the number of identification patches N Σ e . The number of modes N m that have to be taken into account in the modal basis of the cavity is directly linked to the number of identification patches as it must necessarily be greater than N Σ e to ensure the linear independence of the equations. Finally, the number of measurement patches must be equal or higher than N Σ e .

Summarizing: (i) The number of identification patches N Σ e could be arbitrarily chosen but a λ b /4 or λ b /6 criterion is a good choice; (ii) The number of modes of the virtual cavity N m has to be higher than N Σ e ; (iii) The number of measurement patches

The most common way to simulate noise during the measurement of a signal χ is to blur it using two independent gaussian variables, corresponding to the multiplicative and additive errors (respectively α and η) [START_REF] Pezerat | Force analysis technique: reconstruction of force distribution on plates[END_REF].

χ noisy = α • χ + η (19)
Most of the time, the additive part of the error is uncorrelated from the signal itself and can be suppressed using different methods such as the Principal Component Analysis [START_REF] Pearson | On line and places of closest fit to systems of points in space[END_REF] for example. For this reason, we will focus here on the modelization of the multiplicative error. We here develop α so as to let appear the percentage error p

χ noisy = χ • (p + 1) • e i∆ϕ = |χ| • (p + 1) • e i(ϕ+∆ϕ) (20)
where p and ∆ϕ are respectively the amplitude percent error and the phase error. If one consider the amplitude error ∆χ = p • |χ| instead of the percent error, it is possible to choose those two errors at random using the normal laws.

Using this modelization is to consider that the measured value is located in a circle centered on the theoretical value and whose radius characterizes the preciseness of the measurement instrument. Both parameters ∆χ and ∆ϕ can be used to adjust the severity of the error.

Fig. 12 shows the reference and the identified velocity power spectra when the error is applied on the measured velocity (Fig. 12(a)) or on the measured pressure (Fig. 12(b)) independently. For this simulation, an error of 3 % has been used. One can see on these figures that the identification error is not particularly important even if the measurement error is important (∼ 2 dB). The most important thing to notice is that beyond 3000 Hz, the identified velocity is still correct when the error is made on the measured velocity whereas it is not when the error is made on the measured pressure. This let see that it is more important to ensure a good pressure measurement. This phenomenon can ben explained from Eq. ( 13). As shown in this equation, there are two source terms:

P A Σ e and Z A Σ e ←Σ e V A Σ e
. The pressure term comes directly from the measurements and is blurred by the measurement error while the velocity term is balanced by the impedance matrix computed numerically. This balance smoothes the velocity term, leading to better results in the identification when the error is committed on the velocity.

Experimental validation

An experimental validation has been performed to illustrate the previous numerical simulation. As can be seen in Fig. 13(a), a rectangular plate made of aluminum has been glued to a wood frame on 4 cm on the edges of the plate. This 2 mm thick plate is 0.68 m long, 0.38 m wide and has thus a free surface of dimensions 0.6 × 0.3 m. The plate is excited by an electro-dynamic shaker fed with a withe noise and radiates into a quiet room. It is important to underline here that this room is not an anechoic chamber and that no particular treatement has been used to prevent reflections on facilities (on the cartesian robot used to move sensors for example). These measurements are then not in free field conditions. The vibration and boundary pressure fields on the plate surface have been measured directly with a laser vibrometer and a microphone (at 1 mm from the plate) on a regular mesh of 10 by 17 points on the whole surface of the plate. These direct measurements will be considered as the reference fields to evaluate the iPTF method. As the measurements are not synchronous, a phase reference is necessary. For sake of simplicity, the force sensor on the electro-dynamic shaker has been used as a reference. Obviously, any kind of sensor can be used in the same way (accelerometers, microphones). The reference measurements of the boundary pressure and particle velocity are represented in Fig. 14 by the solid black lines. To apply the iPTF method, an acoustic measurement surface surrounding the plate has been defined and discretized into 535 measurements points as presented in Fig. 13(b). The top surface of the acoustic measurement surface is at 5 cm from the plate. As can be seen in Fig. 13(a), the sound pressure and the particle velocity have been measured using a PU probe at the center of each patch constituting the acoustic measurement surface. It is considered here that the value of acoustic quantities at the center of the patch is representative of the mean value on the patch. The measured velocity is always normal to the considered surface and the same phase reference (the force sensor) has been used for acoustic measurements. The mean quadratic pressure and velocity measured on the acoustic measurement surface are also plot in Fig. 14 (dash-dotted green lines) to show that, at this distance of measurements, the differences between boundary and radiated fields are important.

Using measured pressure and particle velocity in the frequency domain and acoustic impedance computed with the FE model of the virtual volume, the pressure, velocity and intensity fields on the surface of the plate have been identified with Eq. (12).

For the sake of comparison, it has been chosen that the identification mesh will be the same as the reference mesh (170 patches). In addition, considering the number of measurements patches (535), the modal basis of the virtual cavity must have more than 535 modes. In the following, the modal basis of the virtual cavity contains all modes up to 11 kHz (2102 modes)

to avoid any problem of convergence.

The identification results are plotted in Fig. 14 (mean quadratic quantities) and Fig. 15 (maps at one particular frequency). Analysing these figures, it can be concluded that the method correctly predict mean quadratic velocity, mean quadratic pressure and radiated power in the whole frequency range under investigation. The particle velocity, which is the primary unknown, is very well estimated and the deduced unknowns (pressure and radiated power) show a good agreement with reference measurements. In Fig. 15 the comparison of the pressure, normal velocity and intensity fields directly measured on the plate with the laser vibrometer and the microphone at a particular frequency are presented. The identified maps shown in Fig. 15 give also a good representation of the source field in level as well as in space distribution. Indeed, the two intensity spots of the source are well estimated and localized. To quantify objectively the quality of the identification, the Frequency Domain Assurance Criterion (FDAC, Eq. ( 21)) can be computed.

FDAC(ω) = {V re f (ω)} T {V id (ω)} {V re f (ω)} T {V re f (ω)}{V id (ω)} T {V id (ω)} (21)
where {V re f (ω)} and {V id (ω)} are vectors of reference and identified velocities at angular frequency ω. FDAC values are included in [-1 1]. A value of 0 indicates no correlation between reference and identified maps, a value of 1 indicates a perfect correlation and a value of -1 a perfect correlation with a phase shift of π. In the case of Fig. 15, the FDAC value is 0.9336 that indicates a really good correlation between velocity maps.

On the whole frequency band between 20 and 1000Hz, FDAC is most of the time higher than 0.8 as can be seen in Fig. 16. Usually correlation between maps are considered to be good for FDAC values higher than 0.75. However, below 100Hz some phase inversions can be noticed (negative values). These phase inversions might be due to a wrong calibration of the PU probe in low frequency where sensitivities of sensors vary quickly (sensitivities are frequency dependent). Above 100Hz, the identification produces really good results and validates the whole procedure of iPTF method.

Conclusion

In this article, iPTF has been used to identify both pressure and velocity fields, yielding the parietal intensity. Until now, the iPTF method has mainly been used to identify the velocity field on the surface of a structure. The mathematical formulation combining the direct and the inverse formulation in order to identify the parietal pressure and intensity fields has been presented. The first results have shown a frequency band for which the iPTF method was not giving exact results mainly due to the presenec of evanescent waves. A parameter called penetration depth has been defined. It explains why identification results underestimate reality in this particular frequency band.

The influence of the discretization of the measurement and identification meshes has been studied. It has been shown that the fineness of the identification mesh is limited by the number of modes taken into account for the decomposition of the problem on the cavity eigen-modes basis.

The stability of the iPTF method regarding measurement errors has been studied. As this method is based on the measurement of both acoustical pressure and particle velocity it is interesting to know which of these two quantities was the most sensitive to noise errors. It has been demonstrated that it is more important to ensure a good pressure measurement than a good velocity measurement.

At last, the validation experiment is an example of a successful application of the iPTF method following the rigorous criteria previously defined.

Fig. 1 .

 1 Fig. 1. Example of system under study. A vibrating surface Σ radiates in an acoustic volume with a rigid surface Λ and an opening Σ .

2 .Fig. 2 .

 22 Fig.2. System under study: a plate excited by a point force radiating into a semi-infinite acoustic medium.

Fig. 3 .

 3 Fig. 3. Definition of the virtual surface Σ surrounding the rectangular plate. Top view. Heigth = 0.045m.

Fig. 4 .

 4 Fig. 4. Intensity fields at 1900 Hz. (a) Reference field; (b) Identified field.

(Fig. 5 .

 5 Fig. 5. Space averaged intensity on the surface of the plate as a function of frequency. (a) Comparison of reference and identified intensity; (b) Identification error in dB.

Fig. 6 Fig. 6 .

 66 Fig.6shows the penetration depth versus frequency. One can see that very low frequency waves can propagate far from

Fig. 7 Fig. 7 (Fig. 7 .

 777 Fig.7(d). These fields are displayed here at 710 Hz where, as can be seen on Fig.5(b), the error is the highest. Indeed, Fig.

Fig. 9 Fig. 8 .Fig. 9 .

 989 Fig. 8. Identification and measurement meshes used in the study on influence of the identification mesh. (a) Identification meshλ b /2; (b) Identification meshλ b /4; (c) Identification meshλ b /6; (d) Measurement meshλ a /4.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Identification and measurement meshes used in the study on influence of the measurement mesh. (a) Measurement meshλ a /2; (b) Measurement meshλ a /4; (c) Measurement meshλ a /6; (d) Identification meshλ b /4.

Fig. 12 .

 12 Fig. 12. Identification of the velocity when gaussian noise is added on (a) measured velocity or (b) measured pressure.

Fig. 13 .

 13 Fig. 13. System under study: an aluminum plate glued to a wood frame and excited by a shaker. (a) Photography of the system; (b) Measurement mesh surrounding the plate.

Fig. 14 .

 14 Fig. 14. Experimental comparison between reference, identified and radiated fields as function of frequency. (a) Sound pressure; (b) mean square velocity

Fig. 15 .

 15 Fig. 15. Experimental comparison between reference and identified fields at 205 Hz on the plate surface. (a) real part of the reference pressure field; (b) real part of the reference velocity field; (c) reference active intensity field; (d) real part of the idenfied pressure field; (e) real part of the identified velocity field; (f) identified active intensity field.

FDACFig. 16 .

 16 Fig. 16. Frequency Domain Assurance Criterion between reference and idenfitied velocity fields as a function of frequency. The dashed line indicates the threshold of good correlation (0.75).

Table 1 .

 1 Rank of the systems

	Mesh Number of unknowns Rank
	λ b /2	74	74
	λ b /4	316	316
	λ b /6	664	427

Acknowledgements

This work was performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007).

[1] Aucejo, M., Totaro, N., and Guyader, J.-L., 2010. "Identification of source velocities on 3d structures in non-anechoic environments: Theoretical background and experimental validation of the inverse patch transfer functions method". J.