Sofiane Bousabaa
email: sofiane.bousabaa@onera.fr

Jean Bulté
email: jean.bulte@onera.fr

Daniel-Ciprian Mincu
email: daniel-ciprian.mincu@onera.fr

Régis Marchiano
email: regis.marchiano@upmc.fr

François Ollivier
email: francois.ollivier@upmc.fr

Ph D Candidate

Sparse Green's Functions Estimation using Orthogonal Matching Pursuit: Application to Aeroacoustic Beamforming

The paper presents a new methodology for the numerical estimation of the Green's functions in complex external aeroacoustic configurations. Computational aeroacoustics is used to propagate multi-frequency signals from focus points to microphones. The method takes advantage of the sparsity of the Green's functions in the time-domain to minimize the simulation time. It leads to a complex sparse linear regression problem. To solve it, the Orthogonal Matching Pursuit algorithm is adapted. The method is first applied on the case of the diffraction by a rigid sphere. Results are studied both in terms of Green's function estimation and aeroacoustic beamforming. They show that the Green's functions are obtained with a good accuracy and enable to localize acoustic sources placed behind the diffracting object. The methodology is then applied on a NACA0012 2D wing in a potential flow for which the Green's function is not known analytically. The use of the reverse-flow reciprocity principle enables to reduce the complexity of the estimation problem when there are more scan points than microphones. It is shown that it is possible to take advantage of the presence of diffracting objects to improve the capability of detection of a sensor array.

= Active set at iteration q relative to source i A (q) = Global active set at iteration q * = Convolution product ⊗ = Cross-correlation product j l , y l = Spherical Bessel functions of order l h (1) l , h

I. Introduction

A source localization techniques aim to detect, localize and characterize acoustic sources on a given area of interest. Often, for practical reasons, microphones cannot be placed in this area (installation problems, intrusiveness...). Acoustic imaging is interested in monitoring this area using an array of microphones that collects spatial and temporal samples of propagating wave fields. The set of acquired microphone signals are then processed in order to extract the relevant information. Several methods can then be used to recover sound sources: holography [START_REF] Ahmed | Holography and its application to acoustic imaging[END_REF], beamforming [START_REF] Veen | Beamforming: a versatile approach to spatial filtering[END_REF][START_REF] Elias | Experimental techniques for source location[END_REF], and time reversal [START_REF] Fink | Time-reversed acoustics[END_REF][START_REF] Deneuve | A coupled time-reversal/complex differentiation method for aeroacoustic sensitivity analysis: towards a source detection procedure[END_REF][START_REF] Mimani | A simulation-based analysis of the effect of a reflecting surface on aeroacoustic time-reversal source characterization and comparison with beamforming[END_REF] among others.

Beamforming methods require knowing the Green's function (GF) of the considered medium. This GF, which represents the acoustic transfer function of the environment, must be determined between each focus point and each microphone position. Classical Beamforming methods use the free-field analytical GFs to describe the propagation of the sound from acoustic sources to the microphones. Several methods have been proposed to deal with more complex situations. An analytical GF is proposed in [START_REF] Amiet | Refraction of sound by a shear layer[END_REF] to take into account the refraction by a thin shear-layer. In [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF], the study of the propagation of sound in small rectangular rooms is investigated. The analytical solution for diffraction of acoustic waves by some simple shapes can also be found in [START_REF] Bowman | Electromagnetic and acoustic scattering by simple shapes[END_REF]. More recently, an analytical GF for annular ducts was calculated for a uniform flow [START_REF] Rienstra | An analytic Green's function for a lined circular duct containing uniform mean flow[END_REF] or for a stepwise constant flow in the radial direction [START_REF] Sijtsma | Green's Functions for In-Duct Beamforming Applications[END_REF]. However, in complex media with flows, scattering objects and reflecting surfaces, an analytic expression of the GF is difficult to obtain. In those cases, the use of beamforming methods based on approximate GFs may result in biased source positions and power estimates or even lead to spurious sources.

GFs can be obtained by doing a Computational AeroAcoustics [START_REF] Lele | Computational aeroacoustics -A review[END_REF][START_REF] Tam | Computational aeroacoustics -Issues and methods[END_REF] (CAA) simulation of the response of the medium to a pulse. However, this strategy will result in a high computational cost when the wavelength is small in comparison with the characteristic length of the problem: a minimum number of points per wavelength are required depending on the order of the solver [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF]. The simulation of high frequency signals may result in a finer mesh and in an increase of the computational cost. For helicopters for example, part of the noise is generated in the plane of the main rotor blades (cf. figure 1). However, acoustic field on the ground results from the diffraction of those sources by the fuselage. The study of shielding effects is difficult because the excited frequencies can reach several kilohertz resulting in a high ratio between the characteristic lengths of the fuselage and the wavelength. Minimizing the computational cost in such configurations is of prior interest.

Depending on the application, the number of focus points can vary from a few tens to several thousands. One may be tempted to minimize as much as possible the simulation time while estimating all the GFs from a single computation. An approach using excited CFD simulations with system identification (CFD/SI) has already been proposed by Polifke et al. [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics[END_REF] and applied to the aeroacoustic characterization of T-junctions [START_REF] Föller | Aeroacoustic Characterization of T-Junctions Based on Large Eddy Simulation and System Identification[END_REF][START_REF] Föller | Advances in Identification Techniques for Aero-Acoustic Scattering Coefficients from Large Eddy Simulation[END_REF]. A good overview of the acoustic modelling using linear time invariant systems is Sovardi and Polifke [START_REF] Sovardi | CFD-Based Modelling of Sound Generation in Ducted Discontinuities[END_REF] and especially the application of system identification [START_REF] Ljung | System Identification[END_REF] techniques for their characterization. More recently, the method was applied by Sovardi et al. [START_REF] Sovardi | Parametric LES/SI Based Aeroacoustic Characterization of Tandem Orifices in Low Mach Number Flows[END_REF][START_REF] Sovardi | Concurrent identification of aero-acoustic scattering and noise sources at a flow duct singularity in low Mach number flow[END_REF] for the concurrent identification of both scattering and noise sources at flow duct singularities. In their approach, the auto-and cross-correlations of the pressure fluctuations are estimated. Then, the unit impulse responses of the multi-port are computed by solving the Wiener filter equation. In this paper, it is also proposed to use CAA simulations and system identification techniques for the estimation of the GFs for external aeroacoustics. However, the cost in terms of CAA simulation is highly decreased by taking into account the sparsity in time of the GFs in such configurations and by estimating all the GFs from a single simulation. The paper is organized as follows. The GFs estimation problem is presented in Sec. II. It consists in the propagation of multi-frequency excitation signals from the focus points to the microphone positions. Focus points are excited as point monopole sources. The direct problem is computed using a single simulation. The GFs, which represent the impulse responses of the system, are then derived from inputs (excitation signals) and outputs (microphone signals). We show that the GFs appear to be solution of an ill-posed deconvolution problem. To regularize the problem, a constraint on the l 0 -norm is added. This enables taking into account the sparsity of the GFs in external aeroacoustical configurations. A methodology for solving the regularized estimation problem is then presented. In this approach, the Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF] algorithm is adapted. The overfitting phase is avoided by using a cross-validation [START_REF] Stone | Cross-Validatory Choice and Assessment of Statistical Predictions[END_REF][START_REF] Hastie | The elements of statistical learning[END_REF] stopping criterion. Results are then presented in Sec. III. The GFs obtained by this method are compared with both the analytical and experimental solutions in the case of a diffracting sphere in a potential flow (Sec. III.A). In terms of aeroacoustic imaging, we show that the GF estimates enable the characterization of both the amplitude and location of an acoustic source hidden below the diffracting sphere. Then, the methodology is applied to extract GFs in the presence of a NACA0012 2D wing in Sec. III.B in a potential flow. The very high number of focus points in this case increases the difficulty to solve the estimation problem. However, the use of the reverse-flow reciprocity principle will allow to apply the method. The GFs are estimated using the proposed methodology and the ability to perform acoustic imaging in presence of flow is studied numerically. We show that the proposed method provides GF estimates with good accuracy and enables a decrease of the computational cost.

II. Green's function estimation problem

In this section, the methodology for the numerical estimation of the GFs is presented. First, a CAA simulation is performed in order to propagate multi-frequency excitation signals from focus points to the microphone positions. Pressure perturbations extracted from this simulation are used to extract the GFs representing the acoustic impulse response of the medium. GFs appear to be solution of an ill-posed system identification problem that is regularized by taking into account the sparsity of the GFs for external aeroacoustical configurations. A methodology for solving the regularized estimation problem is then proposed. It consists of an adaptation of the OMP algorithm with a cross-validation stopping criterion.

A. Simulation of the direct problem

Acoustic model

The propagation of acoustic waves through a moving medium can be modeled by the linearized isentropic Euler equations [START_REF] Bailly | Numerical solution of acoustic propagation problems using linearized Euler's equations[END_REF]:

         ∂p ∂t + u 0 • ∇p + u • ∇p 0 + γp 0 ∇ • u + γp ∇ • u 0 = c 2 0 N S i=1 Θ i ∂u ∂t + u 0 • ∇u + u • ∇u 0 + 1 ρ 0 ∇p -p ρ 2 0 c 2 0 ∇p 0 = 0 p = c 2 0 ρ and c 2 0 = γp 0 ρ 0 (1)
The density ρ , the velocity u and the pressure p designate small perturbations superimposed on a mean steady flow of density ρ 0 , velocity u 0 and pressure p 0 . The coefficient γ designates the ratio of specific heats, and is taken as γ = 1.4 for air. The domain is excited by N S point monopole sources situated at the positions of the focus points:

Θ i (x, t) = δ (x -x i ) × θ i (t) (2)
where δ represent the Dirac distribution, x i is the position of the focus point i and θ i its mass injection source term.

Computational aeroacoustic simulation

ONERA's code sabrina_v0 [START_REF] Redonnet | Numerical simulation of propagation of small perturbations interacting with flows and solid bodies[END_REF] is used to solve the linearized Euler equations in the non-conservative form with disturbances on body fitted structured mesh. This code uses high order finite difference schemes (up to the sixth order), high order filtering schemes (up to the 14th order) and an explicit (third order) time scheme. In the frequency domain, the GFs need to be determined only on a given frequency band [ f min , f max ]. Low-dispersive and low-dissipative acoustic propagation is ensured using at least 10 grid points per wavelength at the highest frequency of interest, f max . To ensure outflow boundary conditions, a combination of a smooth stretching grid and Tam conditions [START_REF] Tam | Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow[END_REF] is used at the boundaries of the domain. To evacuate acoustic waves accurately, sources are placed at least at a distance of one wavelength from the outflow boundaries (relatively to the smallest frequency of interest, f min ). The propagation of waves outside the frequency range [ f min , f max ] may results in higher numerical dissipation and/or spurious reflections at the domain boundaries.

For a given source number i, the forcing term is approached by a normalized Gaussian-distribution:

c 2 0 ∂Θ i ∂t = 1 3 exp -π x -x i 2 2 × s i (t) (3)
where represents the spatial extent of the acoustic source and the {s i } i are the input excitation signals. The parameter must be taken sufficiently large with respect to the cell size in order to ensure computational stability, while being sufficiently small with respect to the smallest acoustic wavelength to reproduce an acoustically compact noise source. For the simulation performed in the following, = 2∆x where ∆x is the mesh step.

Excitation signals

In their survey of excitation signals for FFT based signal analyzers, Schoukens et al. [START_REF] Schoukens | Survey of excitation signals for FFT based signal analyzers[END_REF] explain that, to contribute constructively to the measurement, the energy must lie in the frequency range of interest. Moreover, the use of periodic signals measured over an integral number of periods enables to avoid leakage effects. If a non-integral number of periods are measured or if aperiodic signals are used, leakage effects cannot be avoided leading to a an increase of the needed measurement time to get a specified accuracy. In the following, multisine signals exciting this frequency range with harmonic components every ∆ f are used:

s i (t n ) = N F u=1 a u sin 2π f exc u t n + φ u,i , (4)
where N F is the number of frequencies and

f exc 1 = f min ≤ • • • ≤ f exc N F = f max
are the excited frequencies. In order not to favour any particular frequency, the amplitude of all the harmonics is taken unitary a u = 1 ∀u. The signal is thus T-periodic with T = 1/∆ f . The sampling period ∆t is constant and the discrete time steps are given by t n = n∆t. The number of time steps for one period of signal is called N T and T = N T ∆t.

The phases φ u,i are randomly generated in [0, 2π[ in order to ensure decorrelation between all the source signals. Better decorrelation can be ensured by increasing the number of frequencies N F . However, it will also imply a longer period and thus an increase in the computational cost.

B. Sparse Green's functions estimation problem

Sparsity based regularization

For several sources emitting at the same time, the pressure fields at microphones can be expressed as:

p m (t) = N S i=1 s i * g i,m (t) ∀m, t (5)
where p m (t) is the pressure field measured at the position x m of microphone m, and * is the convolution product

(u * v)(t n ) = N T -1 k=0 ū(t n -t k )v(t k ) ( 6
)
where ū is the periodic extension of signal u ( ū(t nt k ) = u(T + t nt k ) for negative times t n < t k and ū = u elsewhere). The function g i,m (t) will be referred to as the GF between source i and microphone m. It represents the acoustic impulse response of the medium between those two points. For each microphone m independently, GFs between the sources and this microphone appear to be solution of a multi-input, single-output (MISO) system identification problem. As it will be shown in Sec. II.B.2, this problem is ill-posed. In the case of external aeroacoustical configurations, GFs can be considered as sparse. Adding a sparsity constraint to the GFs enables to regularize the problem. The GFs are estimated as a sum of Dirac delta functions corresponding to the direct path and the reflections. The resulting inverse problem can then be solved in the time domain if the input and output signals hold enough information. The problem, with l 0 -regularization, reads as:

g m = argmin g m p m - N S i=1 s i * g i,m 2 2 s. t.
i g i,m 0 ≤ λ 0 [START_REF] Amiet | Refraction of sound by a shear layer[END_REF] where g m stands for the concatenation of all the GFs g 1,m , • • • , g N S ,m between the microphone m and all the focus points and g m is an estimate for g m . The counting norm . 0 represents the number of non-zero components of a signal. This corresponds to the best subset selection problem. The OMP algorithm (also called forward stepwise regression) can be used in the following to obtain an approximate solution. In the following, we omit the m index as the same procedure can be applied independently to each microphone. Details on the adaptation of the OMP algorithm [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF] for GFs estimation are given in Sec. II.C.1 and in appendix B.

Condition number

Ignoring a multiplying factor, the gradient of the criterion with respect to g i (t n ) writes (see appendix A for the derivation):

c i,n = s i ⊗ p - N S j=1 s j * g j (t n ) ( 8
)
where ⊗ is the cross-correlation product

u ⊗ v(t n ) = N T -1 k=0 ū(t k -t n )v(t k ) ( 9
)
where ū is the periodic extension of signal u. Details on the derivation of this expression are presented in appendix A. A fitting descent algorithm tries to minimize this gradient as it will be seen in the next section. The problem can thus be formulated as finding the functions that cancel the gradient and g is solution of the Wiener filter equation:

s i ⊗ N S j=1 s j * g j (t n ) = (s i ⊗ p) (t n ) i ∈ 1, N S n ∈ 0, N T -1 (10)
The number of equations here is N S N T , equal to the number of unknowns (every GF at every time step {g i (t n )} i,n ). However, those equations are not linearly independent and the rank r of the system [START_REF] Rienstra | An analytic Green's function for a lined circular duct containing uniform mean flow[END_REF] is smaller than N S N T . First, all these equations are linear combinations of the equations obtained by simply cancelling the residue:

N S j=1 s j * g j (t n ) = p(t n ) n ∈ 0, N T -1 (11)
meaning that r ≤ N T . Besides, for periodic signals the rank is even lower. The direct discrete Fourier transform F and its inverse F -1 of a signal are respectively given by

U( f k ) = F {u} ( f k ) = N T -1 n=0 u(t n )e -i2πt n f k k ∈ 0, N T -1 , (12)
with

f k = k∆ f = k/T and u(t n ) = F -1 {U} (t n ) = 1 N T N T -1 k=0 U( f k )e i2πt n f k n ∈ 0, N T -1 . ( 13
)
With this expression, and because we are measuring one period of converged signals, the convolution theorem apply. It states that the Fourier transform of the convolution is the pointwise product of the Fourier transforms F {u * v} = F {u} • F {v}. The discrete Fourier transform being a linear operation, applying it to system [START_REF] Sijtsma | Green's Functions for In-Duct Beamforming Applications[END_REF] gives, at frequency f k :

N S j=1 S j ( f k )G j ( f k ) = P( f k ) k ∈ 0, N T -1 , (14)
where S j , G j and P are the discrete Fourier transforms of s j , g j and p respectively. The discrete Fourier transform is invertible and, thus, conserve the rank so that system [START_REF] Sijtsma | Green's Functions for In-Duct Beamforming Applications[END_REF] has the same rank as system [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF]. As only a small number N F < N T of frequencies are excited during the simulation, most of these equations are always true (like "0 = 0"). The non-trivial equations are those corresponding to the excited frequencies

N S j=1 S j ( f exc u )G j ( f exc u ) = P( f exc u ) u ∈ 1, N F , (15)
Thus we have shown that if multisine signals exciting frequencies f exc u , 1 ≤ u ≤ N F are used for the excitation, and assuming that one period of converged signal is recorded, the rank r of system [START_REF] Rienstra | An analytic Green's function for a lined circular duct containing uniform mean flow[END_REF] is smaller than the rank of system (15) that is itself smaller than N F because it contains only N F equations. The problem of GF estimation is thus severely ill-posed (r ≤ N F N S N T ) justifying the need of the regularization term introduced in [START_REF] Amiet | Refraction of sound by a shear layer[END_REF]). An approach for the resolution of (10) was proposed by Polifke et al. [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics[END_REF] for N S = 2 non-periodic input signals for the reconstruction of acoustic transfer matrices in thermoacoustic systems. The methodology was also applied N S = 3 non-periodic signals by Föller et al. [START_REF] Föller | Aeroacoustic Characterization of T-Junctions Based on Large Eddy Simulation and System Identification[END_REF] for the aeroacoustic characterization of T-junctions. It can be generalized for a larger number of periodic input signals. The Wiener filter Eq. ( 10) can be re-written:

Γg = b ( 16
)
where Γ is the correlation matrix between input signals s 1 , • • • , s N S and b is the cross-correlation between input signals and pressure at the microphone:

Γ in, jn = s i ⊗ s j (t n -t n ) i, j ∈ 1, N S n, n ∈ 0, N T -1 b i (t n ) = (s i ⊗ p) (t n ) i ∈ 1, N S n ∈ 0, N T -1 (17)
Auto-correlations (s i ⊗ s i ) (t nt n ) appear in the diagonal blocks of the matrix. They represent a physical limitation to the estimation of the GFs. Hence, the only excited frequencies can be estimated in the process [START_REF] Schoukens | Survey of excitation signals for FFT based signal analyzers[END_REF]. Considering off-diagonal blocks, it shows that the pairwise signals cross-correlations need to be low. The difficulty to solve the deconvolution problem also increases with the number of sources N S . The contributions of all the sources are added in Eq. ( 5) making it difficult to extract GFs between the considered microphone and all the focus points. We observe in practice that the number of frequencies N F required to estimate GFs increases with the number of sources N S .

C. Deconvolution

Algorithm outline

The purpose of the OMP algorithm, proposed by Pati et al. [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], is to find, from a family of vectors {u 1 , • • • , u N }, a sparse linear combination that approximates a given vector v. With α = (α 1 , • • • , α N ) the coefficients of this linear combination, the OMP algorithm solves the problem

α = argmin α v - N i=1 α i u i 2 2 s. t. i α 0 ≤ λ 0 ( 18
)
The vector α is the vector of regression coefficients and the family of vector {u 1 , • • • , u N } is referred to as the dictionary D. The OMP algorithm starts with the zero solution. It iterates in order to minimize the gradient. At each iteration, the component that maximizes the gradient is added to the active set (set of all the non-zero components of α) and a linear system is solved in order to determine the amplitudes of these components.

The algorithm was originally proposed for the wavelet decomposition in order to compute sparse representations of signals in possibly overcomplete dictionaries. In the case of external aeroacoustics, the algorithm can be used to exploit the sparse characteristics of the GFs. Most of the energy of the GFs is concentrated in a few echoes. It is possible to give a good approximation of the time-domain GFs as a sum of pulses, each pulse corresponding to one acoustic ray between the source and microphone positions. It is required to find a sparse representation of the pressure at a microphone p m in the dictionary D consisting in the N S N T signals representing the source signals delayed in time

D = {s (n) i } i,n with ∀(i, n) ∈ 1, N S × 0, N T -1 : s (n) i (t n ) = si (t n -t n ) n ∈ 0, N T -1 , (19)
where si is the periodic extension of signal s i . With these notations, the GF estimation problem of Eq. ( 7) can be reformulated as

g m = argmin g m p m - N S i=1 N T -1 n=0 g i,m (t n )s (n) i 2 2 s. t. i g i,m 0 ≤ λ 0 . (20)
The OMP algorithm can thus be used to solve the GF estimation problem. It requires a training dataset S = (s i ) i , p m m for extracting the GFs. Details on the adaptation of the OMP algorithm for the problem of GFs estimation are given in appendix B. To do so, instead of the naive approach that will consist of computing independently the gradient of the least squares criterion

∀(i, n) ∈ 1, N S × 0, N T -1 : ∂Q ∂g i,m (t n ) ( g m ) s (n) i T p m - N S j=1 N T -1 n =0 g j,m (t n )s (n ) j = s (n) i T µ, (21)
the expression of the gradient given in Eq. ( 8) is used. The computation of the gradient using this formula has a O N S N 2 T complexity if vector µ is computed initially and then used for the calculation of the gradient. The use of expression [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF] enables to reduce the algorithmic complexity to O N S N T log 2 (N T ) by using high speed convolution and correlation [START_REF] Stockham | High-speed Convolution and Correlation[END_REF].

As we will see, the algorithm succeeds well in giving a sparse solution. However, the l 0 -constraint can appear to be too strict in some cases resulting in a bad estimation of the GFs. This is the case when the input and output signals are too short and does not bring enough information for the algorithm to train on and to estimate relevant GFs. In the worst case, some GFs can even be skipped (estimated as zero functions while it is not the case). It appears in practice that the GFs between focus points located in hidden areas are more likely to be skipped. Considering longer multisine signals with a higher number of harmonics N F will improve the capacity of the algorithm to extract those.

Cross-validation stopping criterion

As it is shown in appendix B., the OMP algorithm requires to set a stopping criterion. Hence, the problem is ill-posed and two phases appear during the minimization process: (i) the fitting and (ii) overfitting phases. During the fitting phase, GF estimates are updated in a way that the distance between them and the exact GFs is reduced. The algorithm then reaches an overfitting phase in which the algorithm is trying to extract more information than what the training dataset holds. In their review on modelling issues, Maier and Dandy [START_REF] Maier | Neural networks for the prediction and forecasting of water resources variables: a review of modelling issues and applications[END_REF] explain that the overfitting appears due to the high ratio between the number of equations and the number of unknowns. Stone [START_REF] Stone | Cross-Validatory Choice and Assessment of Statistical Predictions[END_REF] proposed the cross-validation to differentiate between the fitting and the overfitting phases. It requires a second set of input signals and pressure at the microphones, the validation set S = ( s i ) i , p m m . The accuracy of the estimated GFs is calculated as their ability to generalize on the validation dataset that has not been utilized in the training process.

F =141 Valid. N F =281 Valid. N F =561 Train. N F =141 Train. N F =281 Train. N F =561 (i) (ii) (i) (ii) (i) (ii)
This works as follows. As it is explained in Sec. II.B.1, the minimization process is applied independently for each microphone index m. We define Q and Q the criteria on the training and validation set:

Q g m = p m -N S i=1 s i * g i,m 2 2 p m 2 2
and

Q g m = p m -N S i=1 s i * g i,m 2 2 p m 2 2 (22)
The evolution of this two criteria in function of the iteration number are shown in figure 2. In this test case, N S = 11 sources are located in the segment 0 ≤ x s ≤ 2 m at y=0 m every 0.2 m and z=0 m and the microphone is located at position (0, -2 m, 0). Multisine signals (4) in the frequency range [4 kHz,11 kHz] are used for the excitation of the sources. The procedure was applied for different number of frequencies N F = 141, 281 and 561 and the delimitation of the fitting (i) and overfitting (ii) phases relatively to the cross-validation stopping criterion has been highlighted. At each iteration of the algorithm, the GF estimates g m = g i,m i are updated in order to decrease the least squares criterion Q on the training set. Thus, the criterion Q decreases at every step making it difficult to differentiate fitting and overfitting phases. However, the criterion Q decreases, reaches minimum value, then increases when the algorithm starts overfitting. The best GF estimates are chosen as those that minimize the criterion Q on the validation set. The value of the minimum gives information on the quality of the GF estimates. In figure 2, the minimum value reached by the criterion Q is smaller when the number of frequencies is increased meaning that the estimation is better. This is because when the number of frequencies is increased, the algorithm is given more data to train on. The time domain GFs in the various phases of the minimization process is presented in figure 3. For this example, the number of frequencies is N F =141. The time aspects of the GFs are shown after 11, 100 and 200 iteration of the algorithm. The iteration 11, shown in the figure 3-(a), corresponds to the fitting phase. For this iteration number, exactly 11 non-zeros were added and the algorithm has decided to add exactly one component per GF at the time delay corresponding to the propagation from the source to the microphone. The best GF estimates, relatively to the cross-validation criterion, are shown in the figure 3-(b). Some more non-zero components have been added near the main spikes. Those components enable the compensation of the effects of the discretization in time. They appear when the time of propagation of the sound from the source to the microphone is not a multiple of the time sample ∆t. Two spurious spikes are visible at times t = 10 ms and t = 12.5 ms. The iteration 200, shown in the figure 3-(c), corresponds to the overfitting phase. Between iteration 100 and 200, the solution have not been modified too much in the time interval of interest [5.9 ms, 8.3 ms]. However the solution differs more outside this interval where numerous spurious spikes have been added to the GFs. The cross-validation appears to give a good estimation of the limit between the fitting and the overfitting phase in this case. The increase of the sampling frequency, until ∆t = 10 -6 s improves the performances of the OMP algorithm because spikes can be placed more precisely. However, as observed by Berger et al. [START_REF] Berger | Application of compressive sensing to sparse channel estimation[END_REF], it appears that further increase of the sampling frequency strongly diminish the quality of the estimation, probably because of the dropping of the condition number. For all the cases envisaged in the following, the sample time is taken as ∆t = 10 -6 s.

III. Green's functions estimation results

Two test cases representative of external shapes of the fuselage and wings are considered in this section. The first case we address is the diffraction by a rigid sphere. The acoustic source is located at the opposite side of a microphone line array. The pressure field resulting from the diffraction of a monopole source by a rigid sphere can be expressed analytically. In order to prove the advantages of CAA in terms of GF estimation, we will compare CAA estimated GFs with both experimental and analytical ones. The second test case is the diffraction by a NACA0012 wing profile. For this test case, the solution is not known analytically and the GFs obtained using the proposed method will be compared with GFs obtained experimentally. Beamforming results will also be presented and we will see that the method can be applied to real industrial configurations with several thousands of focus points. A. In presence of a diffracting sphere

Experimental setup

The data were obtained in the DLR Aeroacoustic Windtunnel Braunschweig (AWB) anechoic wind tunnel (1.2 m × 0.8 m) and were part of the HeliCopter Action Group HC/AG-24 of the Group for Aeronautical Research and Technology in EURope (GARTEUR) research collaboration organization. The experimental setup is presented in figure 4. The center of the 0.06 m radius sphere is taken as the origin of the spatial system. The sphere is made in aluminium and is held by a bar downstream for reducing the perturbation of the flow. ONERA's SPARC sound source, based on high power electric discharge, was used for the excitation. Details on the characteristics of SPARC source are presented in figure 5. The source generates sound from the cracking of an electric arc. Most of the energy of the acoustic pulse is released in 0.2 ms and its maximum energy lays around 10 kHz. In the following, the considered frequency range is [4 kHz,11 kHz] in order to maximize signal-to-noise ratio. For the present experiment a spark repetition rate of 10 Hz was chosen and 10 s of signal was recorded at a sampling frequency of 250 kHz. This means that approximately 100 samples are recorded per data point. Time signals are then averaged improving the signal-to-noise ratio and enabling a better reproducibility. This procedure enables to reduce the 99% confidence interval at less than ± 1 dB in the spectrum of the acoustic source on the frequency range [4 kHz,11 kHz]. A truncation in time was also performed in order to remove reflections on mounting objects.

The sound source is located at z s = +0.32 m. Inflow measurements were performed with a 1/4 inch microphone with a nose cap at z m = 0.3 m. The microphone was moved in the streamwise and spanwise directions. However, only measurements along the streamwise microphone traverse are used here as the effect of the flow is more important for those positions. The traverse is constituted of 29 positions going from x m = -0.2 m to x m = +0.2 m. Measurements were done at two flow speed: 0 and 45 m/s. As it will be seen, this is sufficient to observe the influence on the resulting GFs. The presence of the holding bar, visible in figure 4, is not taken into account in the simulation nor in the analytic formulation inducing a difference for GFs at microphone positions x m ≥ 0. In the following, comparison is done for position x m ≤ 0 in order to minimize the effect of the bar.

Numerical setup

The simulations were performed with ONERA's sAbrinA-v0 code on a body-fitted mesh. In the useful part of the mesh, the spatial step is 2 mm. It is small enough to ensure more than 10 points per wavelength for the maximum frequency f max = 11 kHz (leading to the minimum wavelength λ min 3.1 cm). A 3% stretching is added at the limits of the domain and combined with [START_REF] Tam | Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow[END_REF] outflow boundary conditions (BC) to evacuate outgoing waves properly. The simulation is led in the presence of a potential flow of v∞ = v∞ u x with v ∞ = 45 m/s the flow far from the sphere and u x the unit vector in the x-axis direction (Mach number M = 0.13). The analytial expression of the mean velocity field is v(x) = ∇ x φ with:

φ(r, θ) = v∞ cos(θ) a 3 2r 2 + r with r = ||x|| cos(θ) = u x • x/||x|| (23)
where r ∈ R + is the radial distance, θ ∈ [0, π] is the polar angle, a is the radius of the sphere and φ is the mean velocity potential. The mesh is body-fitted near the surface of the sphere in order to apply hard wall BC. The minimum mesh cell size is 0.6 mm and the time step is ∆t = 10 -6 s. This results in a maximum CFL number equal to 0.6. It is obtained near the surface of the sphere and is smaller than 1 to ensure the stability. Pressure perturbation signals were stored every ∆t ext = 40∆t = 4 10 -5 s. That gives a frequency of storage f ext = 25 kHz, enough to ensure Shannon criterion. This enables to reduce the calls to files during the simulation. To apply the multisource strategy, the signals are then oversampled back to ∆t = 10 -6 s to enable a good estimation of the GFs.

GF estimation results

In order to validate the CAA simulation approach, we start exciting only one monopole source. This enables to evaluate possible numerical errors and to quantify the quality of the GFs that can be extracted from the simulation. The monopole source is placed at the position z s = +0.32 m and is excited using a multisine signal in the frequency range [4 kHz, 11 kHz] every ∆ f = 125 Hz. The period of the signal extracted from the simulation is T = 8 ms after a transitory time of around t trans 4 ms. The simulation is done without flow or in the presence of a potential flow of 45 m/s as in the experiments presented in the previous section. Phases are generated randomly in [0, 2π[ and the same set of phases is used for both cases. Pressures at the microphone positions are recorded and the GFs are extracted dividing the microphone spectrum by the source spectrum at the excited frequencies.

In the following, the focus is put on 3 microphone positions: • In the line of sight x m = -0.2 m: the source is visible by the microphone. Two main rays are to be considered: one direct ray and a ray resulting from the reflection on the surface of the sphere. The delay between those two rays will imply positive or negative interferences on the considered frequency range. • In the masking region x m = -0.07 m: the source is hidden by the sphere. The main part of the energy is held by two main rays creeping on the surface clockwisely and counter-clockwisely. In the frequency range [4 kHz,11 kHz], interference between those two rays is mainly destructive and the masking effect is maximum. • Opposite to the source x m = 0 m: the source, the microphone and the center of the sphere are aligned. Creeping rays come from all the perimeter of the sphere and interfere in a constructive way creating an important increase of the amplitude. The analytical expression of the pressure field resulting from the diffraction of a monopole by a rigid sphere can be found in Bowman et al. [START_REF] Bowman | Electromagnetic and acoustic scattering by simple shapes[END_REF] book on electromagnetic and acoustic scattering by simple shapes. The holding bar, visible in figure 4, is not taken into account in the simulation nor in the analytical formulation which may imply a difference on the estimated GFs when compared with the experimental ones. The GFs obtained using the experimental, analytical and numerical approaches are compared in figure 7-(a), (c) and (e) respectively for positions x m = -0.2, -0.07 and 0 m. The comparison is performed in terms of modulus and phase of the GFs in the frequency range [4 kHz,11 kHz]. For the three microphone positions, a good agreement is obtained between analytical, numerical and experimental GFs especially on the phase. In the absence of flow, the GFs extracted from the simulation should equal the analytical solutions. However, because of a small mesh dissipation effect at high frequencies, we notice a small offset in amplitude that reaches 0.5 dB at 11 kHz. This error is neglected here but could be reduced more if a finer mesh is used. However, the consequence will be an increase of the computational cost. The error on the modulus for the experimental GFs is smaller than 2 dB on almost all the chosen frequency range. The same procedure is followed in the presence of a 45 m/s mean flow (Mach number M = 0.13). Again, several approximations are to be taken into account. In the simulation and for the analytical solution, a potential mean flow is considered. This may induce a difference with the experimental flow especially because of the presence of the holding bar. Moreover, for the analytical GFs, Taylor's transformation [START_REF] Taylor | A transformation of the acoustic equation with implications for low-speed flight and wind-tunnel tests[END_REF] valid under the hypothesis of a low Mach number M 2 1 is used to take into account the potential mean flow. In the present case, M 2 = 0.017 is in agreement with the hypothesis. The GFs obtained in presence of the flow using the experimental, analytical and numerical approaches are compared in figure 7-(b), (d) and (f) respectively for positions x m = -0.2, -0.07 and 0 m. On all the frequency range [4 kHz, 11 kHz], agreement is very good for the phase between the numerical and experimental GFs. This means that the assumption of a potential flow enables to reproduce correctly the influence of the mean flow on the acoustic diffraction. The error on the analytical solution is already visible at M 2 = 0.017. A phase drift appears for the analytical solution and the flow effect is underestimated by this formulation. This effect is more important for upstream microphone positions (Figure 7-(b)). In terms of magnitude, it appears that for position x m = -0.07 m, the flow induces an increase in the shielding effect of the sphere at 8500 Hz (Figure 7-(d)). This shielding effect is underestimated when using the analytical solutions. Finally, for position x m = 0 m, the flow induces a small increase of the amplitude visible on both the experimental and numerical GFs while the flow correction states that the flow has no influence onto the GF for this position.

!"#$%&'# !" #$%" #$&" #$'" #$(" #" )#$(" )#$'" )#$&" #$%" )!" (#$)'# *#$)'# #$*" #$(" #$!" #" )#$!" )#$(" )#$*" )#$(" )#$!" #" #$!" #$(" !"#$%&'# !" #$%" #$&" #$'" #$(" #" )#$(" )#$'" )#$&" #$%" )!" (#$)'# *#$)'# #$*" #$(" #$!" #" )#$!" )#$(" )#$*" )#$(" )#$!" #" #$!" #$("
Thus, the analytical formulation gives an accurate estimation of the GF in the absence of flow. For more complex geometries, methods such as boundary elements methods [START_REF] Agarwal | Low-Frequency Acoustic Shielding by the Silent Aircraft Airframe[END_REF] or ray tracing methods [START_REF] Agarwal | A Ray Tracing Approach to Calculate Acoustic Sheilding by the Silent Aircraft Airframe[END_REF] can be used to determine precisely the GFs of the propagation in the medium. However, as it is shown in this section, the main difficulty is on reproducing correctly the impact of a complex mean flow. The CAA approach enables to take into account complex flows. However, because it is a volumic method, this generally results in a high computational cost that limits its intensive use in a study. In the next section, we will show how the estimation of GFs from a single CAA simulation enables to reduce computational cost.

Beamforming using estimated GFs

In this subsection, we propose to use the approach presented in Sec. II for the estimation of the GFs in order to solve a problem of imaging with the classical beamforming method on a numerical test case. The setup is the same as in the previous section but only sources and microphone numbers and positions are changed. The 41 microphones representing the sensor array are uniformly placed every 10 mm at z m = -0.3 m and -0.2 m < x m < 0.2 m. On the other side of the sphere, the area of interest is divided into 101 focus points uniformly placed every 4 mm at z m = +0.32 m and -0.2 m < x s < 0.2 m.

These positions are highlighted in figure 8 where the pressure field resulting of the simultaneous excitation of all the source points is presented. In order to obtain the training and validation sets S and S introduced in Sec. II.C. The obtained GFs have been compared with those obtained using the monosource strategy in figure 7. We see a good agreement on the amplitude but the error is more important at the limits of the excited frequency range [4 kHz, 11 kHz] due to side effects. Hence, out of this interval, no information is brought to the model and the solution has no physical sense. On the other hand, there is a good agreement on the phases of the estimated GFs with those obtained with the monosource strategy.

The estimated GFs were then used to perform beamforming. The conventional beamforming operator [START_REF] Elias | Experimental techniques for source location[END_REF] for the focus point i can be expressed as:

P C BF i = m,m G * i,m ( f )Γ m,m ( f )G i,m ( f ) m,m G i,m ( f ) 2 G i,m ( f ) 2 ( 24
)
where

Γ( f ) = Γ m,m ( f ) m,m
is the array cross-spectral matrix at the frequency f , G i,m is the GF between the focus point i and the microphone m and superscript * denotes the complex conjugate operator. Beamforming maps obtained using the estimated GFs are plotted in figure 9 for three positions of the acoustic source x s = -0.1, 0 and +0.1 m and at several frequencies f = 5, 7.5 and 10 kHz. It appears that the use of the free-field GFs conduct to a bias in the estimation of the positions of the source and an underestimation of its power due to the shielding of the sphere. The shielding effect is more important at high frequencies. At 10 kHz, the maximum of beamforming based on the free-field GFs is not at the position of the source and two sources are seen coming from the two edges of the sphere. However, we see that the beamforming patterns obtained using the analytical and estimated GFs are of very good quality. Those GFs enable to identify a source masked by the sphere both in terms of position and amplitude. We can also notice that the beamforming patterns are affected by the presence of the sphere. Its presence increases the secondary lobe amplitudes but reduces the width of the main lobe. At 10 kHz, the maximum side-lobe amplitude are increased to -6 dB when the source is masked by the sphere whereas it only reaches -13 dB when the source is not masked. On the other hand, the width of the main lobe below the sphere is 0.04 m, narrower than in the absence of the sphere (0.064 m).

!"#$%&'# !"# $"# %"# &"# '"# "# ('"# (&"# (%"# ($"# (!"# (#$)'# *#$)'# ")%# ")&# ")'# "# (")'# (")&# (")%# (")&# (")'# "# ")'# ")&#
Finally, the application of the proposed method enabled to reduce the computational cost. Hence, supposing that the GFs are required every 125 Hz, the estimation of the GF for one focus point will require t trans + T = 12 ms. Because one simulation has to be done for each of the N S = 101 focus points, the total simulation time required to estimate the GFs using the monosource approach is t mono = N S (t trans + T) = 1.12 s. The proposed method needs only t multi = 0.216 s of CAA simulation i.e. 18% of the initial computational cost. If the GFs samples are required every 50 Hz, only 9% of the 2.42 s original computational time is required. Concerning the resolution of the estimation problem ( 7), the cost is negligible in comparison with the simulation. The CAA simulation required around 600 hours on 480 CPU whereas the deconvolution process required 24 hours of Matlab computation without parallel implementation.

B. Diffraction by a NACA0012 profile

We now address the case of the 2D acoustic imaging of a NACA0012 wing profile. For this test case, the solution is not known analytically and the GFs obtained using the proposed method will be compared with GFs obtained experimentally. The spatial sampling of the surfaces results in a high number of focus points. However, the use of the reverse-flow reciprocity will enable applying the method.

Experimental setup

Again, the data were recorded in the DLR AWB anechoic wind tunnel (1.2 m × 0.8 m) and the experimental setup is presented in figure 10. The center of the 0.2 m chord NACA0012 profile wing is taken as the origin of the spatial system. The wing is held from its extremities and is supposed infinite in the Y direction. Again, ONERA's SPARC source, is used for the excitation and the same procedure as the one described in Sec. III.A.1 is used to improve the signal to noise ratio. The measurements are performed at a flow speed of 55 m/s.

Numerical setup

The simulations were performed with ONERA's sAbrinA-v0 code on a body-fitted mesh. The spatial step is 2 mm in the useful part of the meshing order to ensure 10 points per wavelength for the maximum frequency f max = 11 kHz. A 3% stretching is added at the limits of the domain and combined with Tam [START_REF] Tam | Radiation and outflow boundary conditions for direct computation of acoustic and flow disturbances in a nonuniform mean flow[END_REF] outflow BC to evacuate outgoing waves properly. The global mesh contains 78.4 10 6 points. The simulation is led in the presence of a potential flow of 55 m/s in order to be compared with the experiments described in the previous section. This time, the potential mean flow is first obtained by the use of the FreeFem++ software [START_REF] Hecht | New development in freefem++[END_REF]. The mesh is body-fitted near the NACA profile in order to apply hard wall BC on the surface of the wing. The minimum mesh cell size is 1.4 mm and the time step is ∆t = 2.5 10 -6 s. The resulting maximum CFL number is 0.65. It is obtained near the surface at the maximum thickness of the wing and is smaller than 1 to ensure the stability. Pressure perturbation signals were stored every ∆t ext = 16∆t = 4 10 -5 s. That gives a frequency of storage f ext = 25 kHz. To apply the multisource strategy, the signals are then oversampled back to ∆t = 10 -6 s to enable a good estimation of the GFs.

GF estimation results

In order to have a reference solution for the CAA simulation approach, we start exciting only one monopole source at various focus points. This enables to evaluate possible numerical errors and to quantify the quality of GFs that can be extracted from the multisource simulation. The monopole source is excited using a multisine signal exciting the frequency range [4 kHz, 11 kHz] every ∆ f = 125 Hz. The period of the signal extracted from the simulation is T = 8 ms after a transitory time of around 3 ms. Phases were generated randomly in [0, 2π[. The pressures at the microphone positions are recorded and the GFs are extracted dividing the microphone spectrum by the source spectrum at the excited frequencies.

Comparisons between the experimental GFs and the ones obtained using CAA simulation are shown in figure 11. Three source and microphone positions are taken under consideration. The sources are placed in the visible upstream region (Figure 11-(b): x s = (-0.15 m, 0, -0.1 m) and x m = (-0.2 m, 0, +0.4 m)), exactly below the wing (Figure 11-(c): x s = (0, 0, -0.1 m) and x m = (-0.1 m, 0, +0.4 m)) and near the trailing edge of the wing (Figure 11-(d): x s = (+0.1 m, 0, -0.04 m) and x m = (0, 0, +0.4 m)). The GFs are obtained with a good accuracy even when the source is hidden below the profile. A small shift in phase appears in every case. This may be due to a more important intrusivity of the SPARC source in the flow which is placed only at 0.1 m from the diffracting object.

Aeroacoustic beamforming

In this subsection, we consider doing aeroacoustic imaging in the presence of a NACA0012 profile wing where all the GFs are estimated from one main simulation. Again, it is a completely numerical test case because no experimental data are used. The setup is presented in figure 12-(a). The first surface, plotted in red, is in the plane z = 0 containing the wing and follows the extrados with -0.2 m ≤ x ≤ 0.2 m and -0.2 m ≤ y ≤ 0.2 m. The second surface, plotted in blue, is in the plane z = -0.1 m below the wing with -0.2 m ≤ x ≤ 0.2 m and -0.2 m ≤ y ≤ 0.2 m. The last surface, plotted in green, is orthogonal to the wing profile at y = 0 with -0.2 m ≤ x ≤ 0.2 m and -0.2 m ≤ z ≤ 0.1 m. The three surfaces of interest are discretized respectively into 1681, 1681 and 1291 points giving a total number of 4653 focus points. A cross-shaped sensor array consisting of 81 microphones is considered. It is included in the plane z = 0.4 m and is at a distance of 0.4 m from the wing. Microphones are uniformly distributed every 0.01 m and the two braces of the array are parallel to x and y-axis.

Reverse-flow reciprocity:

This time, the number of focus points is very large and solving the estimation problem [START_REF] Amiet | Refraction of sound by a shear layer[END_REF] for N S = 4653 is not the best strategy. Hence, in this case the estimation problem is more difficult to solve because it requires a higher number of frequencies due to the fact that more parameters must be estimated for one minimization. The reciprocity [START_REF] Heaslet | Reciprocity relations in aerodynamics[END_REF] principle states that the acoustic response of a medium does not change when the source and microphone are interchanged. This principle doesn't apply directly to the case of a moving medium encountered in aeroacoustics. However, it was pointed out in [START_REF] Lamb | On Reciprocal Theorems in Dynamics[END_REF] that inverting the direction of the steady flow may restore reciprocity.

A numeric validation of this principle is shown in figure 13 in the presence of a NACA0012 diffracting wing. In this example, the estimation of the GF between a source point located near the trailing edge and a microphone is presented. The direct problem consists of exciting the source point using a multisine signal and to measure one period of the resulting pressure signal at the position of the microphone (Figure 13-(a)). The frequency domain GF between the two points is obtained by computing the ratio of the Fourier transforms of those signals. In the reverse approach, shown in figure 13-(b), the flow is inverted and the roles of the microphone and the source point are exchanged. The GF in the frequency domain obtained with the direct and reversed approach are shown in figure 13-(c) and compared with experimental ones. The direct and reversed approaches give the same GF and are in good agreement with the experiments.

X Z !!"# !!!"# !!!!"# !$%&'()*#$%&')"# ((#)+,# )-./#012# !%*#$%&')"# !%&')*#$%&%3)"# !$%&4)*#%&3)"# !$%&')*#%&3)"# !%*#%&3)"# !!"# !!!"# !!!!"# (a
In the present work, this principle enables to inject monopole sources from the microphones to the focus points. The problem of Eq. ( 7) is now solved for each focus point independently and reads as:

g i = argmin g i p i - N M m=1 g i,m * s m 2 2 s. t. i g i,m 0 ≤ λ 0 ( 25
)
where N M is the number of microphones and g i now stands for the concatenation of all the GFs g i,1 , • • • , g i, N M between the focus point i and all the microphones. This is of particular interest in the case of real aeroacoustic imaging applications for two reasons. On one hand, the number of focus points can reach several thousands whereas the number of microphones in the array is around one hundred. The use of the reciprocity conducts to the problem of the estimation of N M GFs simultaneously instead of N S for the direct problem. Less parameters are to be estimated which results in a reduction of the number of frequencies N F required, which also drives the quantity of information given to the model. This also implies a reduction of the period T of the signal and thus a reduction of the CAA simulation time. On the other hand, reciprocity also enables to extract GFs between microphones and focus points that may be located on the surface of diffracting objects because only pressure perturbations need to be recorded at the focus points. This enables to localize the acoustic sources resulting from the interaction of turbulent flows with structures.

The reverse-flow reciprocity was applied here and a 3D view of the performed CAA simulation is shown in figure 12-(b). In order to apply the reciprocity principle, the flow is reversed. Monopole sources are placed at the microphones positions and excited with multisinus (4) with N F = 1401 harmonics in the frequency range [4 kHz, 11 kHz] every 5 Hz. Pressure fluctuations are recorded at every focus point and an estimation problem of the form ( 25) is solved for each of them in order to extract the GFs. The estimated GFs are given in figure 11 and again there is a good agreement with those obtained with the monosource strategy. The more adverse case is the case II (figure 11-(c)) when the source is centered below the wing. Hence, there is an important difference of the GFs amplitude between this focus point and the various microphones making problem (25) more difficult to solve.

On the surface of the wing: The surface in the plane at z = 0 and following the extrados is discretized into 41 points in each direction giving a total of 1681 focus points. The estimated GFs with the proposed multisource strategy are used to perform aeroacoustic imaging on the surface of a NACA0012 wing profile. Results of beamforming are presented in figure 14 both in one and two dimensions. For the 1D-beamforming, only the streamwise brace of the sensor array is used in order to improve the results. Results after deconvolution have also been represented. The deconvolution was performed using the Least Angle Regression [START_REF] Efron | Least angle regression[END_REF] (LARS) model selection algorithm. The process is stopped when 99% of the energy of the acoustic map has been removed. It appears that the free-field GFs with flow is a good approximation in this case. This is because a uniform flow is a good approximation in this case. In practice, the flow cannot always be considered uniform or potential and a mixing layer often appears between focus points and microphones. The proposed method however does not make the assumption of the flow being potential but only require a stationary flow. For a source near the leading edge of the wing, the correction of the flow effect enables to localize the position of the acoustic source. It appears on the beamforming maps that the use of the free-field GFs with flow leads to a slight overestimation of the amplitude of the source. The amplification factor is due to the presence of a second acoustic ray diffracted by the leading-edge of the wing. The overestimation is more important at low frequencies for which the two rays reach the microphones with the same phase. On the surface of the wing (-0.1 m ≤ x ≤ 0.1 m), beamforming obtained with CAA-estimated GFs have levels 6 dB lower than those obtained with flow-corrected free-field GFs. The estimated GFs take into account the fact that, for focus points on the surface of the wing, two times the energy is scattered in the direction of the sensor array which is not true with the free-field GFs. For focus points outside the surface of the wing (x < -0.1 m or x > 0.1 m), the beamforming levels are close. Some differences still remains in the position of the lobes for x > 0.1 m because of the consideration of the wing. The mixing between surfacic and volumic focus points tends to favour the volumic ones. The GFs obtained for surfacic points have higher amplitudes resulting in lower beamforming power estimation. This amplification due to the reflection on the wing is also the reason why the maximum of beamforming using estimated GFs is not at the position of the source. After deconvolution, the estimated GFs enable to get the exact levels of a source located near the leading edge whereas the free-field GFs with correction of the flow lead to a slight over-estimation of 1.2 dB on its amplitude at 5 kHz. However, this error is rather small and, after deconvolution, they both enable a precise determination of the location of the acoustic source. Thus, the free-field GF with flow yield a good prediction for doing imaging in this configuration. Below the wing: We now focus on the plane 0.1 m below the wing which is discretized into 41×41 = 1681 focus points. Beamforming results for a source located below the profil are shown in figure 15. As it can be expected, the use of the free-field GFs does not enable to localize the source below the wing. When the correction of the effect of the flow is applied, the beamforming based on free-field GFs leads to two sources coming from the leading and trailing edges of the wing and the amplitude is underestimated (-20 dB). Beamforming results based on estimated GFs enables to characterize the acoustic source placed below the profile. For focus points located below the wing, only two interfering rays coming from leading and trailing edges contribute to the GF. The beamforming results mainly rely on the phase of those two rays. As a result, the beamforming lobes have smaller full-width at half maximum but the amplitude of secondary lobes is almost equal to the amplitude of the main lobe. Deconvolution enables to extract the position of the source but spurious sources are found at the positions of secondary lobes.

We also notice that the beamforming lobes downstream the position of the sources have higher amplitudes than those situated upstream (Figures 15-(b) and (e)). This is because the trailing edge has a higher masking effect than the leading edge. In the expression of the conventional beamforming operator [START_REF] Hastie | The elements of statistical learning[END_REF], the GF appears with a power of 2 in the numerator and a power of 4 in the denominator favoring masked regions. This effect is even more important when the source is placed in the visible region as shown in figure 16. For that case, secondary lobes have an amplitude greater than the level of the source (Figures 16-(b) and(e)). The correct level of 0 dB is found at the location of the source but the beamforming does not pass by a maximum. Discrepancies on the phases are not large enough to make the beamforming decrease because of the too important ratios between the GFs modulus. After deconvolution, the position and level of the source are found but artifact sources appear below the profile where beamforming levels are high. Free-field GFs with correction of the flow leads to an underestimation of the amplitude of the source and to the detection of a second source coming from the leading edge.

Orthogonal to the wing: 2D beamforming results on the plane orthogonal to the wing are shown in figure 17. Only the streamwise part of the sensor array is used. With the flow-corrected free-field GF, an image source is detected due to the reflection of acoustic waves on the surface of the wing (Figure 17-(d)). Beamforming based on CAA-estimated GFs leads to very high amplitudes in regions where GFs have small amplitude. This happens in masked regions and near the surface of the wing. For a monopole placed near the surface, interferences between direct and image source lead to minimum and maximum of beamforming (Figure 17-(a)). The effect is more important near the surface as the two sources have almost the same distance to the antenna and thus have the same amplitudes leading to increasing interferences.

Broadband beamforming enables to decrease secondary lobes. When the frequency is changed, lobes do not appear at the same position whereas the position of the main lobe remains the same. Beamforming in third octave frequency range [4.47 kHz, 5.62 kHz] is also shown for a unit power source with a flat power spectrum in this frequency range (Figures 17-(b) and (e)). The interference effect is still visible and was also observed in [START_REF] Dougherty | Virtual Rotating Microphone Imaging of Broadband Fan Noise[END_REF] where broadband fan noise beamforming was performed using steering vectors derived from annular hardwall duct modes. The same procedure was applied with broadband beamforming in the whole frequency range [4 kHz,11 kHz] in figures 17-(c) and (f). This time the amplitudes of the secondary lobes have been significantly decreased. Thus, it is better to avoid doing beamforming in volumic regions located not only in masked but also near flat reflecting surfaces.

The reflections on the surface of the wing do not have only negative effects here. Hence, it appears that the presence of the wing enables to increase the power of detection of the array in the Z-axis direction. If the area of interest is far enough from the surface to avoid this interference effect, the presence of reflecting surfaces can be an advantage. A possible application would be, in the framework of wind tunnel experiments, to place reflecting surfaces in order to improve the power of detection of a sensor array. However, those reflecting surfaces have to be positioned in a way that they do not disturb the flow in the wind tunnel. For instance, it was numerically shown in [START_REF] Mimani | A simulation-based analysis of the effect of a reflecting surface on aeroacoustic time-reversal source characterization and comparison with beamforming[END_REF] that mounting a rigid surface on the floor of an anechoic wind tunnel can improve the capability of a line array of microphones to characterize experimental flow-induced noise sources.

The method also allowed to reduce the computational cost. If the GFs are required every 125 Hz, one period T = 8 ms of signal has to be recorded after the transient regime of t trans 3 ms. Using the monosource approach, one simulation of time t trans + T will have to be performed for each of the N S = 4653 focus points leading to a simulation time of t mono = N S (t trans + T) = 51.2 s. The use of the reverse-flow reciprocity only requires to do one simulation for each of the N M microphone. Thus, only t reverse = N M (t trans + T) = 0.891 s of CAA simulation time is required i.e. 1.7% of the initial time. The proposed method needs only t multi = 0.214 s of CAA simulation decreasing the computational cost by again a factor 4. Concerning the deconvolution problem, the cost is negligible in comparison with the simulation. The CAA simulation required around 220 hours on 480 CPU whereas the deconvolution process required 2 weeks of Matlab computation with a parallelization on 4 processors. The parallelization of the deconvolution process ( 25) is trivial because it has to be solved independently on each focus point.

IV. Conclusion and future work

This paper presents an innovative approach for the estimation of Green's functions for external aeroacoustics in the presence of a steady flow. It requires the propagation of excitation signals from focus points to the microphones in one main CAA simulation. It takes advantage of the sparsity of the GFs to decrease the computational cost but leads to a complex sparse estimation problem. This problem can be solved using an adaptation of the OMP algorithm with a cross-validation stopping criterion. We have seen that the method enables the determination of the GFs with good accuracy in the case of the diffraction by a rigid sphere and in the presence of a NACA0012 2D-wing and to predict the effect of flow on those GFs. When the number of focus points is greater than the number of microphones, the use of the reverse-flow reciprocity principle enables a decrease of the complexity of the estimation problem. The estimated GFs allow to localize and characterize an acoustic source located below the diffracting sphere. In the case of the wing, the difficulty to localize acoustic sources placed below was highlighted and we saw that reflecting surfaces can be used to enhance the beamforming results if the relevant GF is used. The method is being applied to slat noise imaging in an open section wind tunnel with the presence of a mixing layer between the wind tunnel model and the sensor array. We expect that the use of the GFs obtained with high fidelity will increase the beamforming accuracy and enable better acoustic map reconstruction to locate and describe unpredicted aeroacoustic sources.

  Nomenclature t = Continuous time, s t n = Discrete time t n = n∆t, s N T = Number of time steps N S = Number of focus points N M = Number of microphones x = Vector of 3D-coordinates x = (x, y, z) s i = Signal of source i, kg/s 2 p m , p = Pressure at microphone m, Pa g i,m , g i = Green's function between source i and microphone m g = Concatenation of all the Green's functions g 1 , • • • , g N S for a given microphone m S

l=

  Spherical Hankel functions of order l δ = Dirac delta function H = Heaviside step function Γ = Cross-correlation matrix Γ A (q)= Submatrix of Γ relative to global active set A (q) c = Set E c is the complement of set E * = Complex conjugate operator

Fig. 1

 1 Fig. 1 Left: Installation effects on a helicopter. Broadband noise is generated in the plane of the main rotor blades. Resulting acoustic field is then diffracted by the fuselage. Right: Fan noise imaging on an installed turbofan. The propagation of acoustic waves is affected by the fan cowl and the flow making it difficult to localize acoustic sources.

Fig. 2

 2 Fig. 2 Cross-validation stopping criterion. The training (---) and validation (---) criteria are plotted as a function of the model complexity • 0 . Three numbers of frequencies N F =141, 281 and 561 are considered.

Fig. 3

 3 Fig. 3 Time GFs obtained at (a) iteration 11 during the fitting phase, (b) at iteration 100 at the minimum of the validation criterion, and (c) at iteration 200 during the overfitting phase. Vertical dashed lines represents the time of arrival of the first and the last echoes obtained analytically.

  Fig. 4 (a): Diffraction by a rigid sphere in the AWB wind tunnel. (b): Schematic view.

Fig. 5 ONERA

 5 Fig. 5 ONERA's SPARC acoustic source.

Fig. 6

 6 Fig. 6 Diffraction of a point monopole source by a rigid sphere without flow (left) or in presence of a potential flow of 45 m/s (right). The source is represented with the violet diamond whereas microphone positions are represented with gray squares. The monopole source is excited using a multisine signal in the frequency range [4 kHz, 11 kHz] every ∆f=125 Hz (f exc 1 =4 kHz, f exc 2 =4.125 kHz, . . . , f exc 57 =11 kHz).

  x m =0 m, 45 m/s flow

Fig. 7

 7 Fig. 7 Comparison between GFs obtained experimentally (--), analytically (---) and using CAA simulation with one point monopole excited (+). Subfigures (a), (c) and (e) correspond respectively to microphone positions x m =-0.2, -0.07 and 0 m in absence of flow. Subfigures (b), (d) and (f) show the results obtained in presence of a mean flow of 45 m/s. For those last three subfigures, the GFs obtained by applying the proposed multisource strategy have also been represented. For the experimental (--) and analytical (---) GFs, a linear interpolation was done on the phase to remove phase jumps.

  2, two simulations were performed. For the main simulation, corresponding to the training set S, excitation signals used are multisines as in Eq. ((4)) with N F = 1401 harmonics in the frequency range [4 kHz, 11 kHz] every ∆ f = 5 Hz. After a transient regime of t trans 4 ms, one period T = 200 ms of signal is recorded. For the validation set S, sources are excited in the frequency range [4 kHz, 11 kHz] every 125 Hz giving a period T = 8 ms. The total simulation time required for the estimation of the GFs using the proposed multisource approach is thus t multi = 2t trans + T + T = 216 ms.

Fig. 8

 8 Fig. 8 Presure field resulting by the simultaneous excitation of N S monopoles at the focus points (z s =+0.32 m). The perturbations of pressure are recorded at the microphones positions (z m =-0.3 m). The simulation is done in presence of a potential flow of 45 m/s. The sources are represented with the violet diamonds whereas microphone positions are represented with gray squares. The sources are excited using a multisine signals in the frequency range [4 kHz, 11 kHz] every ∆f = 5 Hz (f exc 1 =4 kHz, f exc 2 =4.005 kHz, . . . , f exc 1401 =11 kHz).

Fig. 9

 9 Fig. 9 Beamforming below a rigid sphere for three source positions (x s =-0.1, 0 and +0.1 m) and at three frequencies f=5, 7.5 and 10 kHz. The solid lines (--) correspond to beamforming done with GFs estimated by the proposed method. Dashed (---) and dotted (• • • • • • ) lines correspond respectively to beamforming based on analytical GFs with and without taking into account the presence of the sphere. A correction was applied on analytical GFs in order to reproduce the effect of the 45 m/s potential flow.

Fig. 10 (

 10 Fig. 10 (a): Diffraction by a NACA0012 profile in the AWB wind tunnel. (b): Schematic view (right).

Fig. 11 (

 11 Fig. 11 (a): positions of the source and microphone for the three cases envisaged for the comparison of the GFs. Positions are given in the XZ plane because all the points belong to the plane y=0. (b), (c) and (d): show the comparison between GFs obtained experimentally (--) and using CAA simulation with one point monopole excited (+) or by applying the proposed multisource strategy (•) for the three cases in the presence of a mean flow of 55 m/s. For the experimental (--) and analytical (---) GFs, a linear interpolation was done on the phase to remove phase jumps.

Fig. 12 (

 12 Fig. 12 (a) 3D representation of the three areas of interest and their discretization into 4653 focus points. The microphones of the cross-shaped microphone array used for acoustic imaging have also been plotted. (b) Pressure perturbation signal resulting from the excitation of 81 acoustic monopoles placed at the positions of the microphones of the array. The monopoles are excited using multisine signals in the frequency range [4 kHz, 11 kHz] every ∆f = 5 Hz (f exc 1 =4 kHz, f exc 2 =4.005 kHz, . . . , f exc 1401 =11 kHz). The mean flow of 55 m/s is reversed for the reciprocity principle to apply.

Fig. 13 Fig. 14 (

 1314 Fig. 13 Reverse-flow reciprocity principle. In the simulations, the monopole source is excited using a multisine signal in the frequency range [4 kHz, 11 kHz] every ∆f = 125 Hz (f exc 1 =4 kHz, f exc 2 =4.125 kHz, . . . , f exc 57 =11 kHz). (a) Direct approach. (b) Reverse approach. (c) Comparison of the GFs with experiments.

Fig. 15 (

 15 Fig. 15 (a), (d): Acoustic beamforming results on the surface of the NACA0012 wing profile for a source located at (x s , y s , z s )=(0, 0, -0.1 m) with unit power. The 2D acoustic maps are obtained using classical beamforming using CAA-estimated GFs. Positions of the source and microphones have also been represented. (b), (e): 1D beamforming results performed on the curve corresponding to y=0. Only the streamwise part of the array is used for 1D beamforming to improve the sensitivity. (c), (f): Results after deconvolution using the LARS algorithm. The considered frequency is 5 kHz ((a)-(c)) and 10 kHz ((d)-(f)).

Fig. 16 (

 16 Fig. 16 (a), (d): Acoustic beamforming results on the surface of the NACA0012 wing profile for a source located at (x s , y s , z s )=(-0.15 m, 0, -0.1 m) with unit power. The 2D acoustic maps are obtained using classical beamforming using CAA-estimated GFs. Positions of the source and microphones have also been represented. (b), (e): 1D beamforming results performed on the curve corresponding to y=0. Only the streamwise part of the array is used for 1D beamforming to improve the sensitivity. (c), (f): Results after deconvolution using the LARS algorithm. The considered frequency is 5 kHz ((a)-(c)) and 10 kHz ((d)-(f)).

Fig. 17

 17 Fig. 17 Beamforming maps on the plane orthogonal to the NACA0012 wing profile for a source located at (x s , y s , z s )=(+0.1 m, 0, +0.1 m). Beamforming is done using estimated GFs ((a)-(c)) or free-field GFs with correction of the effect of the flow ((d)-(f)). Only the streamwise part of the sensor array is used to improve the sensitivity. Positions of the source and microphones have also been represented. The beamforming is done at the frequency 5 kHz ((a), (d)), in the octave band [4.47 kHz, 5.62 kHz] ((b), (e)) or in the whole frequency range [4 kHz, 11 kHz] ((c), (f)). In every case, the total power of the source is 0 dB.

Acknowledgments

Experimental data for the sphere test case were part of the HeliCopter Action Group HC/AG-24 of the Group for Aeronautical Research and Technology in EURope (GARTEUR) research collaboration organization. Experimental data for the NACA 0012 profile were part of the Applied Vehicle Technology AVT-233 working group of the Science and Technology Organisation (STO) of NATO. Authors want to thank P. Lebigot and F. Desmerger for their help in the conception and realization of the experiments. The authors also thank D. Blacodon for the valuable technical discussions.

Appendix A. Gradient of the estimation criterion

The convolution product u * v between two discrete time signals u and v is defined as

where ū is the periodic extension of signal u ( ū(t nt n ) = u(T + t nt n ) for negative times t n < t n and ū = u elsewhere). Using this notation, the Least Squares (LS) criterion reads as:

Here, we introduced matrix A i that reproduce the convolution with source signal s i :

with si is the periodic extension of signal s i . With these notations, the gradient relatively to GF g i is

A i g i .

The matrix A i T is the transpose of matrix A i . It reproduces the cross-correlation with source signal s i :

where ⊗ is the cross-correlation product. More details on the convolution and cross-correlation products can be found in Gray and Goodman's [START_REF] Gray | Convolution and Correlation[END_REF] book on Fourier transforms. Dropping the multiplicative constant, the expression of the gradient c i,n is:

Appendix B. OMP algorithm

The OMP algorithm [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], also called forward stepwise selection, is used to solve the l 0 -constrained minimization problem (( 7)). The l 0 -norm, or counting norm, is the number of non zero components. Thus, for a given λ 0 we are looking for the solution that has the lowest number of non-zero components. Instead of searching through all possible subsets (combinatorially too expensive), OMP algorithm intends to find a good path through data. At each iteration q of the algorithm, one non-zero component is added to one of the GF estimates. The added component is the one which maximizes the gradient expressed in Eq. (( 8)). A projection on the active set is then performed and all the coefficients are updated in order to decrease the criterion of Eq. ((7)) as much as possible. The minimization process is applied independently for each microphone this is why no microphone index m is present (See Sec. II.B.1). For λ 0 fixed, the different steps of the procedure are described as follows:

Step 0 Initialization of the iteration index q = 0, the algorithm starts with the intercept:

where g i is the estimate for GF g i . Set A i is the active set relative to source i. It contains all the time steps n where g i (t n ) is not zero. The global active set A contains all the tuples (i, n) for which g i (t n ) is not zero. Superscript (q) is here to differentiate between the value at the previous iteration and the new one.

Step 1 Increment the iteration index q ← q + 1. Compute the gradient using (( 8)) based on GFs estimated at the previous step (q -1). Then determine the argument of its maximum in absolute value outside the active set and add it to the belonging active sets:

step [START_REF] Veen | Beamforming: a versatile approach to spatial filtering[END_REF] We can now look for GF estimates as a finite sum of Dirac delta functions:

The vector of coefficients β (q) = β (q)

i,n (i,n)∈A (q) ∈ R q , representing the amplitudes of the spikes, minimizes the LS criterion:

The minimum is obtained for β (q) solution of the system:

Where: Γ

Because we are doing a projection, the solution can completely change in one iteration of the algorithm. It appears that the matrix Γ A (q) is a submatrix of the cross-correlation matrix Γ presented in Eqs. (( 17)). The inversion of the linear system at each step can be avoided. Indeed, Γ A (q-1) is a Grammian matrix of a set of linearly independent time signals {s i (tt n )} (i,n)∈A (q) and Γ A (q) is a Grammian matrix corresponding to the previous set of vectors increased by one. The inverse matrix Γ -1 A (q-1) at iteration (q -1) can be used for constructing by blocks the inverse matrix at the next step Γ -1 A (q) . This enables to avoid matrix inversion at each step [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF].

Step 3 Finally GF estimates are updated:

Remaining step Repeat step 1-3 while i g (q) i 0 = q < λ 0 .