The factorization method and Capon's method for random source identification in experimental aeroacoustics

Experimental aeroacoustics is concerned with the estimation of acoustic source power distributions, which are for instance caused by fluid structure interactions on scaled aircraft models inside a wind tunnel, from microphone array measurements of associated sound pressure fluctuations. In the frequency domain aeroacoustic sound propagation can be modelled as a random source problem for a convected Helmholtz equation. This article is concerned with the inverse random source problem to recover the support of an uncorrelated aeroacoustic source from correlations of observed pressure signals. We show that a variant of the factorization method from inverse scattering theory can be used for this purpose. We also discuss a surprising relation between the factorization method and a commonly used beamforming algorithm from experimental aeroacoustics, which is known as Capon's method or as the minimum variance method. Numerical examples illustrate our theoretical findings.


Introduction
In experimental aeroacoustic testing a solid object (e.g., a model of an aircraft component) is placed inside a wind tunnel, and the fluid structure interaction between the flow field and the object generates sound pressure fluctuations, i.e., aeroacoustic noise. The raw acoustic time signal is recorded by an array of microphones and further post-processed to obtain correlation data in the frequency domain. Based on these correlation data one then seeks to localize and quantify the power distribution of the aeroacoustic sources (see, e.g., [4,27,34]).
In this work we restrict the discussion to subsonic homogeneous unidirectional flow fields in free space, and we use the convected Helmholtz equation to model the propagation of timeharmonic aeroacoustic pressure waves. These waves and the associated sources are usually considered as random functions. Following [17] we model the aeroacoustic acoustic pressure signal as a Hilbert space process with zero mean and a covariance operator that acts as a multiplication operator. The inverse source problem then amounts to reconstructing the source power function from the corresponding covariance operator of the aeroacoustic pressure signal on the microphone array. In practise the latter can be estimated from microphone array measurements by Welch's method [36]. In our analysis we assume that the covariance operator corresponding to an idealized continuum model for the microphone array is available. It has been shown in [17] that this inverse random source problem has a unique solution. For further results on inverse random source problems for time-harmonic acoustic waves, which are not directly related to aeroacoustic imaging, we refer, e.g., to [2,3,12,24].
Various reconstruction procedures have been discussed for correlation based random source identification in aeroacoustics. Covariance fitting (see, e.g., [5,37]) estimates source powers directly from correlation data of the observed acoustic random pressure signal by minimizing a suitably regularized output least squares functional. A faster and therefore more popular reconstruction technique in experimental aeroacoustics is beamforming (see, e.g., [6,11,30,31,35]). Instead of solving the inverse source problem for all source positions at once, beamforming estimates the source powers at individual source positions separately. In particular DAMAS [7] and CLEAN-SC [33], which combine beamforming methods with suitable postprocessing schemes to improve the spatial resolution of the reconstruction, have become standard tools. Both, covariance fitting and beamforming, have recently been reviewed from a continuous perspective in [17].
In this work we focus on the localization of extended aeroacoustic source power functions. We show that a variant of the factorization method from inverse scattering theory can be used to recover the support of a random source from correlations of aeroacoustic pressure fluctuations. The factorization method has been introduced in the framework of inverse obstacle scattering [19] and inverse medium scattering [20] by Kirsch. It has subsequently been attracting a considerable amount of attention over the past twenty-five years. We will show that the mathematical structure of the covariance operator of the aeroacoustic pressure signal is closely related to the structure of the Born approximation of the far field operator for the inverse medium scattering problem (see, e.g., [21]). This will be used to establish the theoretical foundation of the factorization method for the aeroacoustic inverse source problem. On the other hand, we will see that the inf-criterion of the factorization method and also the traditional imaging functional that is obtained from Picard's criterion (see, e.g, [22]) is closely related to another well-established beamforming algorithm that is known as Capon's method or as the minimum variance method (see, e.g., [9,25,26]). In particular, our results give a mathematically rigorous theoretical interpretation of the reconstructions obtained by Capons's method. We show for the first time that (for our idealized measurement setup and in the absence of measurement errors) Capon's method recovers the correct support of locally strictly positive source power functions.
This article is organized as follows. In Section 2 we briefly recall some basic facts on solutions to the convected Helmholtz equation, and we introduce the stochastic model for the aeroacoustic source problem with uncorrelated extended sources. In Section 3 we establish the main result of this work, which is a theoretical justification of the factorization method for reconstructing the support of the source power function from the covariance operator corresponding to the radiated sound pressure fluctuations. In Section 4 we discuss the relation between Capon's method and the factorization method. Some numerical results on experimental data are presented in Section 5.

The aeroacoustic inverse source problem
, be relatively open domains such that Ω ∩ Σ 0 = ∅. In the following M represents an idealized (d − 1)dimensional continuous measurement array, and Ω is supposed to be a region in space that contains all possible aeroacoustic sources.

The convected Helmholtz equation
The basic sound propagation model that is used in experimental aeroacoustics to describe time-harmonic sound waves inside a subsonic homogeneous flow field u ∈ R d is the convected Helmholtz equation. Given a source term Q ∈ L 2 (Ω), the associated sound pressure field p satisfies ∆p where k := ω/c is the wave number, ω the frequency, and c the speed of sound. Here, subsonic means that the Mach vector m := u/c satisfies |m| < 1. In the following we also use the notation β := 1 − |m| 2 . Throughout, | · | denotes the Euclidean norm on R d . We will assume that the convective field u is aligned with the i.e., w 0 is a weak solution to a standard Helmholtz equation with wavenumber k/β.
Proof. This may be verified by direct calculation.
Using Proposition 2.1, the Sommerfeld radiation condition, which determines outgoing solutions to the standard Helmholtz equation on unbounded domains (see, e.g., [10, p. 18]), can be transferred to the convective Helmholtz equation. Let U ⊂ R d be a bounded domain, and let p ∈ C 2 (R d \ U ) be a solution to Then we call p radiating if it satisfies the radiation condition uniformly with respect to all directions x/ |x| ∈ S d−1 . Similarly, using Proposition 2.1 the fundamental solution of the convected Helmholtz equation can be obtained from the fundamental solution for the standard Helmholtz equation (see, e.g., [10, p. 19 and p. 89]). To simplify the notation, we define the Mach norm on R d by Therewith, the fundamental solution of the convected Helmholtz equation is given by 0 denotes the Hankel function of the first kind of order zero. For later reference, we note that Proof. This follows from the one-to-one correspondence between radiating solutions to the stan- The next proposition gives an integration by parts formula that is a consequence of Green's second theorem (see, e.g., [10, p. 19]). A complete proof can be found in Appendix A of [29].
Finally, we transfer the Helmholtz representation formula for radiating solutions of the standard Helmholtz equation (see, e.g., [10, Thm. 2.5]) to radiating solutions of the convected Helmholtz equation. Again a complete proof, which employs Proposition 2.1, can be found in Appendix A of [29].
Then, for any x ∈ R d \ U , we have (2.5)

The random source process
In experimental aeroacoustics sources are usually considered as random functions. Following [17], we use a Hilbert space process, i.e., a bounded linear operator where (X, A, P) is the underlying probability space, to model the source problem. Then, the associated random pressure signal is given by where g is the fundamental solution from (2.2). Using (2.3) we see that g(x, · ) is square integrable for any x ∈ M, and thus (2 and the covariance operator Cov[Q] : L 2 (Ω) → L 2 (Ω) is the unique self-adjoint and positivesemidefinite operator that satisfies It is commonly assumed in experimental aeroacoustics that the random source has zero mean and is spatially uncorrelated.
We note that in the special case when q ≡ 1, a process Q that satisfies Assumption 1 is called a white noise process. Since Cov[Q] is symmetric and positive-semidefinite, the source power function q is real-valued and nonnegative a.e. in Ω. For any x ∈ M, the pressure signal p(x) is a scalar, complex random variable with E[p(x)] = 0, and the correlation between two observation positions x, y ∈ M satisfies Cov(p(x), p(y)) = Cov (Q(g(x, ·)), Q(g(y, ·))) = Cov[Q]g(x, ·), g(y, ·) L 2 (Ω) Accordingly, the covariance operator of the aeroacoustic pressure signal C(q) : L 2 (M) → L 2 (M) is given by Using (2.3) it follows that, for any q ∈ L ∞ (Ω), the covariance operator C(q) is a Hilbert-Schmidt operator (see [17,Pro. 2.2]).
In experimental aeroacoustics finite dimensional approximations of the covariance operator C(q) are obtained from microphone array measurements by estimating the covariance matrix of the microphone signals. This estimation is usually carried out by Welch's method [36]. We are interested in the inverse source problem to reconstruct the support of the source power function q ∈ L ∞ (Ω) from observations of C(q) ∈ HS(L 2 (M)). We note that in [17] it has been established that in fact even q is uniquely determined by C(q). In this work we will show that a variant of the factorization method from inverse scattering theory can be utilized to recover the support of q from C(q).

The factorization method in aeroacoustic source imaging
From now on we let q ∈ L ∞ (Ω) with q ≥ 0 a.e. in Ω be a fixed source power function, and we denote by C(q) the associated covariance operator. Following [16,23] we distinguish the support, the inner support, and the outer support of q. These notions will be used in the characterization of the support of q in terms of C(q) below.

Remark 3.2.
In the following we denote by D ⊂ R d the interior of supp(q). Following [15] we say that the source power function q is locally strictly positive on D, if for each x ∈ D there exist ε x , r x > 0 such that B rx (x) ⊂ D and where B rx (x) denotes the ball of radius r x centered at x. If this is the case, and if R d \ supp(q) is connected, then inn supp(q) = supp(q) = out supp(q) (see [16,Cor. 2.5]). In general, the outer support is basically the support plus the holes that cannot be connected to infinity.
Our goal is to reconstruct D from the covariance operator C(q) under minimal assumptions on q. The techniques that we use have been developed for time-harmonic inverse scattering problems in [19,20,21,22], and we further apply ideas that have been proposed in [8,14,15]. We define the operator H D : Then the adjoint H * D : L 2 (M) → L 2 (D) of H D is given by Therewith, the covariance operator C(q) can be decomposed as The following range identities are the first ingredient of our reconstruction method.
The second ingredient of our reconstruction method is the following characterization of the support of the source power function in terms of the point sources g( · , z)| M , z ∈ Ω, and the range of the operator H D .
the volume integral on the right hand side of (3.2) vanishes. Moreover, the function w is a radiating solution of the homogeneous convected Helmholtz equation on R d \ B ε (z). Hence, (2.5) can be applied to conclude that This yields the assertion.
(b) Suppose that z ∈ Ω \ out supp(q), and that g( · , z)| M ∈ ran(H D ). Then there is ψ ∈ L 2 (D) such that is an entire radiating solution to the convected Helmholtz equation Thus, v must vanish identically on R d (see [10, p. 28] for the corresponding result for the standard Helmholtz equation, and use one-to-one correspondence between radiating solutions to the standard Helmholtz equation and radiating solutions to the convected Helmholtz equation by means of the Lorentz transformation). Therefore, v vanishes on R d − , and we find by analytic continuation that v is zero on R d \ (out supp(q) ∪ {z}). Here we used that R d \ out supp(q) is connected. This means that However, the left hand side of (3.3) is unbounded on B 1 (z) ∩ (R d \ (out supp(q) ∪ {z})), while for the right hand side we obtain that D
Combining Theorems 3.3-3.4 gives the following result.
(a) If z ∈ inn supp(q), then Since C(q) : L 2 (M) → L 2 (M) is compact, self-adjoint, and positive-semidefinite, it has a complete orthonormal eigensystem. We assume that the (possibly finite) sequence of positive eigenvalues (λ j ) j∈N is in decreasing order such that each eigenvalue is repeated according to its multiplicity, and we denote by (ψ j ) j∈N the corresponding sequence of orthonormal eigenfunctions. Accordingly, the nonzero eigenvalues and the corresponding eigenvectors of C(q) 1/2 are given by ( λ j ) j∈N and (ψ j ) j∈N , respectively.  Proof. We first show that g( · , y) ∈ ran(C(q) 1/2 ) for any y ∈ Ω . (3.6) To see this, let φ ∈ ker(H * D ), i.e., By assumption there is an open subset B ⊂ D such that ess inf(q| B ) > 0. Thus, M g(x, y)φ(x) dx = 0 for any y ∈ B , and by analytic continuation this holds even for any y ∈ Ω. Accordingly, g( · , y) ∈ ker(H * D ) ⊥ = ran(H D ) for any y ∈ Ω .
Usually in practice, only a finite number of microphones at positions x 1 , . . . , x M ∈ M is available to measure the random pressure fluctuations. A self-adjoint, positive-semidefinite correlation matrix C ∈ C M ×M , which approximates the covariance operator C(q), can be obtained from these observations using Welch's method [36]. We denote by (λ j , ψ j ) 1≤j≤M an orthonormal eigensystem of C such that the eigenvalues are in decreasing order and counted with multiplicity. Let 0 < M 0 ≤ M be the number of positive eigenvalues. Then we define the imaging functional I fac : Ω → R of the factorization method by Denoting by C † and (C 1/2 ) † the pseudoinverses of C and C 1/2 , respectively, (3.7) can be rewritten as According to Theorem 3.7, the values of I fac (z) should be much smaller for z ∈ Ω \ out supp(q) than for z ∈ inn supp(q). The imaging functional in (3.8) is closely related to Capon's method [9] from seismic imaging. In the context of correlation based aeroacoustic source mapping this method is also known as the minimum variance method (see, e.g., [25,26]). This relationship is quite surprising as Capon's method was originally derived from a totally different viewpoint. In the next subsection we discuss this observation in some more detail.

Capon's method
In aeroacoustic source identification imaging functionals I : Ω → R are usually defined on a source region Ω ⊂ R d as introduced at the beginning of Section 2. Imaging procedures that map focus points z ∈ Ω in the source region directly to an image value I(z) independently of all other focus points z ∈ Ω, z = z, are called beamforming methods. As they do not require evaluations of the source problem, a main advantage of beamforming methods is that they are usually very fast. On the other hand, these methods typically rely on heuristic arguments and can only capture the main features of the source power function rather than providing an exact reconstruction.
Following the usual presentation in the field (see, e.g., [32]) a beamforming imaging functional is defined by with a steering vector w(z) ∈ C M that depends on the focus point z and is assumed to satisfy the constraint w(z) * g(z) = 1 . The latter is often called unit gain. A particular beamforming method is therefore fully determined by its steering vector. The steering vector of Capon's method is given by This yields the imaging functional I Cap : Ω → R, which coincides with the discrete imaging functional I fac of the factorization method in (3.8).
In the traditional derivation of Capon's method (see, e.g., [18, p. 358]) it is assumed that the correlation matrix C is positive-definite, and the steering vector w Cap (z) is obtained, for any z ∈ Ω, as the solution of the constrained optimization problem, If g(z) ∈ ran(C), which is always the case when C is positive-definite, then w Cap (z) from (4.2) is a solution to (4.4) (see, e.g., [28, pp. 443-447]). The minimization problem (4.4) is usually motivated as follows. According to our model in Section 2 the pressure signals at the microphone positions x 1 , . . . , x M ∈ M are zero-mean, complex random variables. We collect them in a vector-valued random variable p := [p(x 1 ), . . . , p(x M )] ∈ C M with zero mean. Then Cov(p) = E(pp * ) = C, and considering the inner product with the steering vector w(z) * p, one seeks to reduce noise as well as signals coming from other focus points z ∈ Ω, z = z, whereas the signal originating at the focus point z should not be dampened. The latter requirement is ensured by the unit gain constraint (4.1). The first requirement is enforced by minimizing the variance of w(z) * p. Therefore, Capon's method is also known as minimum variance method. Minimizing the variance yields which explains the cost functional in (4.4). Finally, we note that Capon's beamformer can equivalently be written as whenever g(z) ∈ ran(C). This is the discrete analogue of the infimum in the inf-criterion of the factorization method in Corollary 3.6.

Numerical examples
We conclude our investigations with some numerical results for the factorization method on experimental data, and we compare these reconstructions to results that are obtained using two commonly used conventional beamforming schemes. The dataset was measured at the cryogenic wind tunnel in Cologne (DNW-KKK) on a 1 : 9.24 scaled Dornier 728 half model [1]. Figure 5.1 shows the setup of this experiment. The measurement array (on the right hand side of the picture) consists of 134 microphones, which are flush-mounted at the wall of the wind tunnel. The Mach number of the flow field is |m| = 0.125 (i.e., the wind speed is |u| = 43 m/s), the angle of attack (i.e., the inclination angle of the wing's cross section plane) is 9 • , and the temperature is 11 • C. The raw output data of the experiment consists of time series of acoustic pressure fluctuations for each microphone with a total measurement time interval of 30 s and a sampling frequency of 120 kHz. These time series are then post-processed to obtain an estimated correlation matrix C ∈ C 134×134 using Welch's method [36] with a Hann weighting window, a block size of 1024 time samples and an overlap factor of 0.5. We evaluate the imaging functional I fac of the factorization method from (3.7), which coincides with the imaging functional of Capon's method from (4.3), on a two-dimensional plane Ω ⊂ Ω that is aligned to the cross-section of the aircraft wing. The map size is 1.05 m × 1.45 m and the grid spacing is 1 cm. We compare these results to source maps obtained using two conventional beamformers with and without diagonal removal (cf., e.g., [32]), which are defined by Diagonal removal is often used in experimental aeroacoustics to lower the effect of wind noise due to turbulent boundary layers directly at the microphone array (see, e.g., [32]). The imaging functionals in (5.1)-(5.2) can also be written as where · F denotes the Frobenius norm.
To further reduce noise effects, the imaging outputs I(f, z) of each imaging functional for single frequencies f are averaged over a frequency band B, i.e., we evaluate the sum Here we consider third octave bands with center frequency f1 /3Oct , which are defined by With a frequency resolution ∆f = 120 kHz 1024 ≈ 117 Hz, the number of discrete frequencies that are contained in a frequency band B = [f 1 , f 2 ] is given by We note that in the notation of Section 2 this corresponds to ω = 2πf , i.e., to a wave num-  (5.3) those frequency bands contain 16, 23, and 32 frequencies. At 8 kHz the factorization method provides a significant improvement in spatial resolution when compared to conventional beamforming with and without diagonal removal. On the other hand, the reconstructions of the factorization method contain more low frequent artifacts in regions apart from the wing, where no sources are to be expected (e.g. in the top left corner of the source maps). At 12 kHz and 16 kHz only one or two dominating sources are recovered by the conventional beamformers, while the factorization method reconstructs regularly spaced sources on the leading edge of the wing and a localized source at the end of the wing flaps. All main source mechanisms are visible in the source maps of the factorization method. The processing time for the factorization method is comparable to that of the conventional beamformers.
Recalling the equivalence of (3.8) and (4.3) we conclude from (4.4) that the imaging functional of the factorization method (or equivalently of Capon's method) gives the smallest values among all beamformers maintaining the unit gain constraint (4.1). This cannot be seen in Comparing this with (3.7) shows that source components corresponding to large eigenvalues of the correlation matrix dominate the reconstruction that is obtained by the conventional beamformers, while the factorization method emphasizes on source components related to smaller eigenvalues of the correlation matrix. This, and our theoretical results from Section 3, might be used to explain the larger number of reconstructed source components that is obtained by the factorization method in Figure 5.2 at 12 and 16 kHz. On the other hand, small eigenvalues of the correlation matrix do not affect the stability of the conventional beamformers, while they lead to instability of the factorization method (without further regularization), which yields artifacts in reconstructions from noisy data.

Conclusions
In this article, we have demonstrated that a variant of the factorization method from inverse scattering theory can be used to reconstruct the support of aeroacoustic random sources from correlations of observed pressure fluctuations. We established a rigorous characterization of the support of the source power function in terms of the correlation data. Moreover we have shown that the factorization method is closely related to Capon's method, which is a well-established beamforming method in experimental aeroacoustics. This unexpected relationship gives a new theoretically rigorous interpretation of the reconstructions that are obtained by Capon's method. Our results basically say that Capon's method recovers the correct support of the source power function, at least when the latter is locally strictly positive.