Inverse Scheme for Acoustic Source Localization using Microphone Measurements and Finite Element Simulations

We propose an inverse scheme for acoustic source localization based on solving the corresponding partial differential equation in the frequency domain (Helmholtz equation) by applying the Finite Element (FE) method. This allows us to fully take into account the actual boundary conditions as given in the measurement setup. To recover the source locations, an inverse scheme based on a sparsity promoting Tikhonov functional to match measured (microphone signals) and simulated pressure is proposed. The properties of this inverse scheme and its applicability to source localization in the low frequency range will be demonstrated.


Introduction
Acoustic source localization is amain task in the development of newp roducts.In the last years, considerable improvements have been achievedi na coustic source localization using microphone arrays.Forthis purpose, usually acoustic beamforming is used to determine source locations and distributions, measure acoustic spectra for complete models and subcomponents, and project results from the array to farfi eld points.The fundamental processing method, FrequencyD omain Beamforming (FDBF) [ 1] is robust, fast, and renders continuous source distributions as continuous images.To this end, the beamforming map a is computed by a(g) = ḡT Cg (1) with g denoting the steering vector [2], C the cross spectral matrix (CSM)o fthe microphone array signals and T the transposed.Here, abar denotes complexconjugation.
The CSM is computed from the measured microphone signals and is modeled by Received20September 2017, accepted 25 June 2018.
with N the assumed number of sound sources.Fixing only the source strength σ k to be nonzero results in and the source computes by By normalizing ḡT k g k = 1, one obtains the correct source location with strength σ k .T he theoretical resolution is givenb yt he Rayleigh limit [3] as well as by the Sparrowl imit [4].Drawbacks are caused by the Point Spread Function (PSF)o ft he array,w hich is the convolution of the spatial array response with as ingle point source.The PSF has as trong main lobe as well as strong side lobes for lowfrequencies, so that weaker sources may be hidden.Therefore, as econd signal processing step is required to convert ar aw FDBF map into as ource map.This can be done by application of an overall scaling factor,k nown as the integration technique [1] or deconvolution by,e.g., DAMAS [5], Clean-SC [6], SC-DAMAS [7], L1-Generalized Inverse Beamforming (L1-GIB) [8].Adetailed comparison and arecent overviewondeconvolution techniques and application to 2D and 3D sound source localization can be found in [9,10].There, the different techniques are compared with respect to position detection, source levelestimation and computational time.The main findings for source localization in afree radiation environment can be summarized as follows: (1) SC-DAMAS provides the best source map at highest computational costs; (2) Clean-SC has the best trade-off between fast computation and correct source detection; (3) for acase of twoorthogonal microphone arrays only L1-GIB results in good 3D source maps at reasonable computational costs.Despite these advances in beamforming techniques, it has to be mentioned that major limitations arise from the source model.In the standard approach, the sources are modeled as monopoles or/and dipoles and the steering vector g describing the transfer function between source and microphone signal is modeled by Green'sf unction for free radiation.Adifferent choice of steering vectors can improve the results as demonstrated in [2].Thereby,ar easonable enhancement in the correct estimation of source location could be obtained.
Af urther challenge for these methods are source localization at lowf requencies and in environments with partially or fully reflecting surfaces, for which beamforming techniques do not provide physically reasonable source maps.In such cases, the steering vectors have to be adapted to takethe reverberant environment into account.In arecent publication, twoapproaches have been investigated [11]: (1) modeling reflections by aset of monopoles located at the image source positions; (2) experimentally based identification of Green'sfunction.Thereby,the best results could be obtained by using the formulation with the experimentally obtained Green'sfunctions.
The second approach for source localization can be attributed to the class of time reversal (TR) techniques.Originally,TRhas been invented as an experimental technique [12,13] and has then be applied to manya pplications involving, e.g., underwater acoustics [14], geophysics [15], biomedical imaging [16], non-destructive testing (NDT) [ 17] and also to acoustic sound source localization [18,19] as well as aeroacoustic source localization [20].TR relies on the linearity and areciprocity property of wave propagation phenomena in non-dissipative media.Therefore, recording wave signals and propagating them back in time will result in ar e-focusing of the time-reversed wavesa tt he source positions.This technique has been shown to be quite robust, sometimes even in the presence of very high noise levels [17].In [18], the time-reversal acoustic sink (TRAS)approach has been experimentally applied to the audible frequencyr ange using different transient sound sources in ar everberation room.There, af ocal spot of as eventh of aw avelength wasa chievedb ya pplying TRAS.Hence, it can be concluded that this approach in combination with numerically solving the wave equation may be apromising method in acoustic imaging and localization of acoustic and vibrational sources.In the presence of attenuation, TRAS does not perform well, though.
The third category of methods for source localization is provided by solving the inverse-source problem based on ac onstrained minimization problem, where ac ost functional is minimized under the constraint that the forward problem (waveequation with source term)isfulfilled.This approach assumes that the material and geometric properties of the domain of interest are known, and aims at find-ing the position and strength of all sources.This approach can be performed both in time and in frequencyd omain.Main application towards acoustics is the reconstruction of surface velocity distributions using the boundary element method for numerically solving the wave equation [21,22,23].Therefore, this approach is often called the inverse boundary element method.Since apart of the sound field decays quickly (evanescent waves),t he reconstruction of particular vibration patterns is ill-posed and needs appropriate regularization.In [24] different regularization methods such as Tikhonovregularization with the L-curve criterion or generalized cross-validation, have been discussed and applied to stationary and rolling tire noise reconstruction.In particular,t he L-curvec riterion provided realistic source reconstruction, and especially for measured data, where spatially correlated measurement errors such as transducer phase and position error are present, showed to be very robust.
Adetailed overviewofall different approaches including au nified formalism, ac ategorization and very interesting comparisons based on numerical and experimental data is provided in [25].According to this categorization, our approach belongs to inverse methods based on iterative Tikhonovr egularization including sparsity of the reconstruction.It relies on solving the wave equation in the frequencydomain (Helmholtz equation)b yapplying the Finite Element (FE) method.Forthis purpose, anyother numerical method, e.g. the Boundary Element (BE) method can be applied, which may be even more efficient with respect to CPU time, depending on the particular scenario.In doing so, we can consider the actual boundary conditions as giveninthe measurement setup.To recoverthe source locations, we apply an inverse scheme based on asparsity promoting Tikhonovf unctional to match measured (microphone signals)a nd simulated pressure.Aclear advantage of such inverse methods is its ability to fully takeinto account realistic geometry and boundary condition scenarios, as well its straightforward generalizability to situations with convection and/or attenuation.
In summary,itcan be concluded that acoustic source localization is astrong research topic in various application fields and alarge number of scientists from different communities have contributed by asubstantial amount of special methods.As stated in [25], there is no general source localization approach and the most appropriate choice depends on the combination of pre-knowledge and application case.
The rest of the paper is organized as follows.In Section 2t he physical/mathematical model, which is able to consider sources both inside the domain and on its boundary,ispresented.The optimization approach based on the adjoint method and its numerical scheme are discussed in Section 3. Numerical results on asimplified car model in 2D are presented in Section 4. Finally,wesummarize our approach and provide an outlook on further research in Section 5.

Physical and mathematical model
Assume that the original geometry of the setup and Fourier-transformed acoustic pressure signals p ms i (ω)(ωbeing the angular frequency, i = 1,...,M)measured by microphones at positions x i are given.
Since the identification is performed for each fixed frequency ω,t he dependencyo nωis neglected in the notation.To mimic the Sommerfeld radiation condition on a finite domain Ω ⊆ R d , d ∈{ 2 , 3 } ,ap erfectly matched layer (PML)i su sed and the following partial differential equation (PDE)k nown as the Helmholtz equation on with the wave number k ∈ R,t he searched for interior acoustic sources σ in (x)inthe domain Ω sc ⊆ Ω acou and with appropriately chosen complexv alued functions η x , η y ,η z (for details see [26]).Note that the real part of We want to stress that the acoustic pressure p is acomplex valued function.First of all, the weak form of ( 5),( 6) is derivedbytesting with an arbitrary complexv alued function v ∈ V and integrating by parts where the function space V and its dual V * will be specified later on.More generally,i no rder to include delta pulses as well, we consider not only sound sources as regular functions of the space variable, butaselements of the dual V * ,sothat (7) becomes The inverse problem under consideration is to reconstruct σ in , σ bd (oro ne of them)f rom microphone array measurements of the pressure at certain points x 1 ,...,x M ∈ Ω \ Ω sc ,s ot he measurements takep lace at ac ertain distance from the sources.
Here σ in , σ bd may be (a) point sources, i.e. sums of delta pulses δ x j for x j ∈ Ω sc ; (b) sums of smooth peak functions g j (x), e.g., Gaussian peaks e −|x−x j | 2 /s 2 j with givenvariances s j for x j ∈ Ω sc ; (c) smoothly distributed sources; as well as combinations of these cases.All these cases can be covered (case (c) in an approximate sense)b ythe following ansatz, that will be used throughout this paper a j e iϕ j µ j (10) with ap riori givenf unctionals µ j ∈ V * as well as searched for amplitudes a 1 ,a 2 ,...,a N ∈ R and phases Here, N = N in + N bd with N in the number of possible sources in the domain Ω sc ⊂ Ω and N bd their number on the boundary.This can be seen as follows: (a) Using the delta distributions δ x j concentrated at the points x j with x j ∈ Ω sc ,j = 1 ,...,N in (11) results in provided functions in V are continuous at the points x j so that their point evaluations makesense.To this end, we use the Sobolevspace V = W 1,r (Ω)with r>d .(b) Smooth functions (c) The use of finite element nodal basis functions N j (x) on Ω satisfying N j (x k ) = δ jk ,and allows to approximate smooth interior and boundary sources a j e iϕ j N j with a j = a(x j ),ϕ j = ϕ ( x j ) ,j = 1 ,...,N in , a j = ã(x j ),ϕ j = φ ( x j ) ,j = N in + 1,...,N, and use V = W 1,2 (Ω).In all cases, elliptic regularity provides W 2,2 smoothness of p away from the source domain, so that point evaluations ( 9) makesense and we end up with the solution space loc (Ω \ Ωsc )} p .I ndeed, well-posedness of the forward problem for k 2 outside a discrete set of eigenvalues accumulating only at infinity, follows from spectral theory for compact operators.

Optimization based source identification
The number N of sources is unknown and actually the mentioned identification problem can be considered as ad iscretization of an infinite dimensional inverse problem with af orward operator mapping from the space of regular Radon measures M(Ω) = C 0 (Ω) * (particularly in case (a) above)t oR M via the parameter-to-state map S : M(Ω) → U , σ → p,a nd the measurement operator O : U → R M , p → (p(x 1 ),...,p(x M )),see also [27] for the wave equation in time domain.By compactness of the observation operator,the inverse problem under consideration is ill-posed.
Fitting of the parameters by means of Tikhonovregularization amounts to solving the constrained optimization problem min p∈U,a∈R N a j e iϕ j µ j ,v V * ,V ∀v ∈ V, a = (a 1 ,...,a N ), ϕ = (ϕ 1 ,...,ϕ N ), and In here, the exponent q ∈ (1, 2] enhances sparsity of the reconstruction, when chosen close to one.Moreover, the regularization parameters α, β have to be chosen appropriately.According to the sequential discrepancyprinciple [28], we choose β = α = α 0 2 −m with m the smallest exponent such that following inequality is fulfilled, where ε denotes the measurement error.I ndeed, according to [29], we can expect this to lead to a convergent regularization method.
The box constraints on the phases ϕ i can be realized, e.g., by abarrier term in the cost functional, i.e., by replacing the original Tikhonovfunctional by with some penalty parameter ρ>0.This also helps to avoid phase wrapping artifacts.
To derive the first order optimality conditions and the gradient of the reduced cost function j(a, ϕ) = J p(a, ϕ),a,ϕ , where p = p(a, ϕ)issupposed to satisfy the constraint the following Lagrange functional is considered L(a, ϕ, p, z) = J (p, a, ϕ) + A(p, z) Due to regularity of the constraint, aminimizer has to satisfy the following optimality conditions, with some adjoint state z: = αq a j q−1 sign(a j ) + Re e iϕ j µ j , z V * ,V Here ( 19) is just the state equation ( 15),w hereas ( 18) is the adjoint equation for z = z(a, ϕ), whose strong form (interms of z)is Therewith, in order to carry out, e.g., some gradient or quasi Newton method for solving this optimality system, the gradient of the reduced cost functional is computed via (and analogously for the derivative w.r.t ϕ i ), where the first equality holds by definition, the second equality by the fact that p(a, ϕ)satisfies the constraint, and the last equality by the state and the adjoint equation.This allows us to compute all 2N components of the gradient by just solving twoPDEs, the state and the adjoint equation.
Evaluation of these gradients according to ( 16),( 17) requires the quantities µ j , z V * ,V for z solving the adjoint equation (18).After discretization of (18),zwill be afinite element function with N a the finite element basis functions.Thus, with precomputed values m ja = µ j ,N a V * ,V evaluation of the above mentioned quantities can be simply carried out by amatrix vector multiplication Note that in case (c),the matrix M = (m ja ) j,a=1,...,n eq is the FE mass matrix, whereas in case (a),itisjust the identity matrix.
In our special parametrization of the sources (10), a j are actually signed amplitudes, hence the true physical quantities compute from a and ϕ as follows.
Amplitude : Forour practical realization, aminimization by agradient method with Armijo line search is applied.Since the differential operators for the state equation (see (15)) and the adjoint equation (see (20))a re the same, the FE system matrix for both PDEs is the same, which results in ac omputational highly efficient solution process.Furthermore, the computational time does not depend on the number of microphones nor on the assumed number of possible sources.Finally,t he inverse scheme results in a source map both for amplitude and phase and in addition the reconstructed acoustic field is provided.Nevertheless, it should be mentioned that the applicability of the scheme towards computational time is mainly restricted to the low frequencyrange.The discretization effort and therefore the number of degree of freedoms in 2D is of the order O(h −2 ) (in3DO(h −3 )),where h is the mesh size being determined by (assuming linear finite element basis functions) In (21) c is the speed of sound and f max the highest frequencyo fa coustic sources.Ad etailed description of the implemented scheme including apseudo-code is presented in the Appendix.

Numerical results
The numerical example consists of asimplified car as displayed in Figure 1a nd Figure 2, where we position two acoustic sources with equal intensity,one near aside mirror and one near awheel housing.The chosen frequencyi s5 00 Hz, so that the wavelength is about 0.68 ma ssuming as peed of sound of 340 m/s in air.Asdisplayed in Figure 1, the microphones of setup 1are positioned at the top and at both sides with at otal number of 23 measurement points.The three line arrays were positioned according to [30].The microphone spacing wasset to 0.3 m, satisfying the spatial aliasing theorem [31].In setup 2( see Figure 2) just 11 microphones have been used.The microphones were placed closer to the car and the microphone spacing has been increased to 0.6 m, knowing that the design rules for beamforming arrays do not hold.This microphone arrangement is just used for the inverse scheme to demonstrate that the method is sensitive to the positioning of the microphones.
The floor as well as the boundary around the car is sound hard (fully reflective)a nd on the remaining three  sides free radiation is modeled using the mentioned PML technique.To realistically simulate the measured values at the microphone positions, we perform computations on afi ne grid with 250.514 degrees of freedom for setup 1 (respectively 80.990 for setup 2) with aPML region which is twice as thick as the one used in inverse computations later on.Moreoverweadd random noise to this simulated data.
In afi rst step, we demonstrate that the identification with beamforming algorithms does not work well for such ac ase.To this end, we use the simulated acoustic pressure values at the microphone positions (without noise)as input.Figure 3d isplays the source maps of setup 1o btained with the classical beamforming algorithm and with DAMAS using Green'sfunction of free radiation for computing the steering vectors.Thereby the amplitude of the sources were normalized to the maximum of the original source.
Classical beamforming findsalot of fictious sources, which cannot be distinguished from the true sources.Due to the coherence between the twos ources, DAMAS will have problems to findthe correct source position (see Fig- ure 3b).However, activating for each computation just one source, DAMAS performs quite well as expected.
Thereby,the position of the source near the side mirror is well identified (see Figure 3c), whereas the position of the source near the wheel housing does not match the correct one (see Figure 3d).Here, modified steering vectors, e.g. based on measured Green'sfunctions as described in [11], may improve the result.Summarizing, it can be stated that deconvolution approaches such as DAMAS are robust and perform well.However, instead of the real source distribution (see Figure 4a)s ingle point sources are identified.Therefore, the obtained source amplitudes are higher.Forthe identification with our inverse scheme, we construct adifferent grid with just 62.854 degrees of freedom in the case of setup 1( respectively 20.431 for setup 2).In addition, the thickness of the PMLo nt his grid is just half the size as in the data simulation.We set the starting values of our regularization and penalty parameters to α = 1, β = 1a nd ρ = 1.Furthermore, for the exponent q,w ec hoose the value 1.1.All other parameters of the inverse scheme are chosen as described in the Appendix.The searched for sources are modeled by delta peaks for each of the 2.881 nodes within the source region.The total elapsed CPUtime (stand-alone PC withanIntel Xeon E5-2640, 2.50 GHz processor)f or getting the source map wasabout 45 minutes for setup 1(respectively 15 minutes for setup 2).Thereby,15outer loops have been performed (see Algorithm 1i nt he Appendix)r educing the regularization parameters α, β, ρ to 6.10352 • 10 −5 .
To mimic the real situation, we add noise to obtain different signal-to-noise ratios (SNR)for the simulated pressure values at the microphone positions computed on the finegrid.These values are then used in our inverse scheme, which performs the identification of the sources on the coarser grid.The results of the reconstructed sources are displayed in Figure 4a  of the sources were normalized to the maximum of the original source.One can clearly see that we get areasonable reconstruction of the sources for both setups and for all cases of SNR.The amplitude for each source is also giveninT able I.It should be noted that some artificial sources not present in the forward simulation are reconstructed, an effect which gets more pronounced as the SNR decreases.The main advantage of this approach is that ad etailed source distribution both in amplitude and phase is identified, and finally with these information we can perform a numerical simulation to obtain the acoustic field.
Applying nowt he identified source distributions and performing sound field computations result in the acoustic fields as displayed in Figure 6and      The acoustic field computed by the reconstructed sources compares well to the original acoustic field (see Figure 6a and Figure 7a), and just for the lowest SNR deviations are clearly observed (see Figure 6f and Figure 7f).Here, the measurement positions of setup 2provides better results compared to setup 1.
Table II summarizes the achievederror value between the acoustic pressure values at the microphone positions computed according to (14) and normalized to the L 2norm of the measured pressure values.
Here, we can observethat the relative L 2 error is quite low, even in the case of SNR being 2dB.This fact clearly indicates that we should increase the number of microphones or use other microphone positions with ah igher sensitivity to improve our identification of the acoustic sources.
These results demonstrate the applicability of our inverse scheme.Af urther improvement can most probably be obtained by using microphone positions with ah igher sensitivity as can be seen by using the different microphone positions in setup 1a nd 2. Therefore, the question of optimal positioning of the microphones arises, which will be subject of our further research.

Conclusion
Acoustic source localization, especially in the lowf requencyr ange, is still at opic of ongoing research with strong practical applications.Amain restriction is the assumption of Green'sf unction for free radiation.Some of currently available signal processing techniques overcome this limitation by using experimentally determined or numerical computed Green'sf unction to adapt the steering vectors.Our developed inverse scheme fully solves the Helmholtz equation with the correct boundary conditions.The first numerical results demonstrate the potential of our approach and clearly showits applicability in the lowfrequencyrange.The main advantage of this approach is that adetailed source distribution both in amplitude and phase is identified, and finally with this information we can perform an umerical simulation to obtain the acoustic field.To further improve our method, we currently investigate in the optimal positioning of the microphones.
D and b is close to unity also in the PML region.Homogeneous Neumann conditions on the whole boundary of Ω are imposed, part of them just as asimple waytoclose the outer boundary of the PML domain, part of them to model the sound-hard boundary part of the acoustic domain.Furthermore, sound sources on the surface of the obstacle are modeled by n • D∇p = σ bd on ∂Ω.

Figure 1 .
Figure 1.Computational setup 1for identification of the sources around acar.

Figure 2 .
Figure 2. Computational setup 2for identification of the sources around acar.

Figure 3 .
Figure 3. (Colour online)Source maps obtained with beamforming algorithm applied on setup 1( just the source region is displayed).The green dots showt he correct source positions, and the red spot the identified ones.(a) Classical beamforming, (b) DAMAS (both sources), (c) DAMAS (left source), (d) DAMAS (right source).

Table I .
Amplitude of the identified source at the true source position normalized to the maximum of the original source (unit dB).

Table II .
Relative L 2 error between measured and computed pressure values at the microphone positions.