An Improved Algorithm for Noise Source Localization Based on Equivalent Source Method

Near-field acoustical holography (NAH) based on the equivalent source method (ESM) is an efficient method applied for noise source identification. As l2-norm-based regularization cannot produce a satisfactory solution of the ill-conditioned problem in high frequency, the conventional ESM is restricted to relatively low frequency, and the resolution of conventional ESM inmiddle to high frequency remains a limitation open to investigation. This article presents an algorithm known as improved functional equivalent source method (IFESM), designed to enhance the resolution of conventional ESM. This method is developed in the framework of wideband acoustical holography (WBH) combining with functional beamforming (FB). Through numerical simulations, it is proved that the proposed method can localize noise with higher resolution compared with WBH and conventional ESM, and the ghosts on noise source map can be suppressed effectively.The validity and the feasibility of the proposed method are manifested by experiments including single-source and coherent-source localization.


Introduction
Since near-field acoustic holography (NAH) was first proposed by Williams et al. [1][2][3] in 1980s, it has been widely studied by the worldwide scholars.The equivalent source method (ESM), which is a kind of NAH reconstruction algorithm for sound sources with arbitrary shape, has some particular advantages such as robust calculation, fast speed, and easy implementation [4,5].Its core idea is to substitute the real source with a series of virtual equivalent sources.Then the equivalent source amplitude could be reconstructed with the transfer matrix and acoustic pressure vector, in which the pressure on reconstruction surface can be calculated for noise source identification [6,7].However, in the application of the ESM-based NAH, the number of microphones is always less than that of equivalent sources.This leads to a preliminary key problem: ill-posedness, which is often encountered in the inverse problems [8,9].Typically, classical  2 -norm-based regularization methods such as Tikhonov regularization [10][11][12] and truncated singular value decomposition (TSVD) [7,13,14] are employed to find the optimal solution of the equivalent source amplitude.As the solution obtained is a least square solution, the ESM with  2 -norm regularization focused on relatively low frequency, and the resolution of ESM-based NAH in middle to high frequency remains a limitation open to investigation.
Recently, sparse regularization based on minimizing a  1norm cost function has emerged in the field of noise source identification [15,16].The  1 -norm-based regularization has also received much attention, particularly motivated by compressive sensing theory under sparsity condition, where the minimum  1 -norm solution is the sparsest solution [17][18][19].A sparse sequence means the number of nonzero components is quite small, compared with its dimension.Under some conditions, it is possible to reconstruct the sparse signals with significantly less measurements (i.e., microphones) than classically required.In ESM-based NAH, the equivalent source amplitude vector can be thought as sparse because most components equal zero.In a sparsity framework, the sparse set of equivalent source amplitude is calculated by solving an objective function containing a penalty on the  1norm of equivalent source vector.For example, Chardon et al. [20] used sparsity and compressive sensing principles for near-field acoustic holography to identify vibrating source velocity with fewer microphones than Tikhonov regularization required.Suzuki [16]  This paper focuses on improving the resolution of ESM in middle to high frequency.Firstly, a different filter function and a modification of the step length are applied into WBH to improve its reconstruction accuracy, and the improved WBH is known as iWBH in the following sections.Next, cross spectrum of the equivalent source amplitude vector is decomposed into a specific form.Finally, the high order matrix function [22,23] is introduced to suppress the ghosts and reduce the beam-width.Through numerical simulation, it is proved that iWBH can reconstruct noise source with better accuracy than WBH, and the proposed method can localize noise sources more accurately in middle to high frequency.

Algorithm Description
The ESM-based NAH uses a source model in terms of a set of elementary sources and solves an inverse problem to reconstruct the complex amplitude of all elementary sources.According to the theory of ESM, the arrangements of equivalent source surface, measurement surface, and reconstruction surface are shown in Figure 1.
Assume that there are  monopole point sources on the equivalent source surface and  microphones set up on the measurement surface.In the framework of ESMbased NAH, the sound field excited by the real sources is considered to be equivalent to the sound field excited by these equivalent sources.Then the relationship between the measured pressure and the amplitude of these equivalent sources can be written in matrix vector notation as where p (size  × 1) denotes acoustic pressure measured by the microphone array.q (size  × 1) represents the complex amplitude of these equivalent sources.A (size  × ) denotes the transfer matrix between acoustic pressure at the microphone positions and equivalent source positions.Each element of the transfer matrix is computed using the monopole radiation formula: with  being the distance between the considered sourcemicrophone couple (in m) and  the acoustic wave number (in rad/m).When using ESM to reconstruct the sound field, Tikhonov regularization is one of the efficient methods applied to find the optimal solution of the equivalent source amplitude.The Tikhonov regularization calculates the vector Shock and Vibration 3 q by making a tradeoff between residual norm ‖p − Aq‖ 2 and the solution norm ‖q‖ 2 .It can be described as minimize where ‖ ⋅ ‖ 2 2 denotes the  2 -norm and  ∈ [0, +∞) is the regularization parameter.Note that the transfer matrix is A underdetermined, because the number of microphones is limited in most acoustic studies.With this additional term ‖q‖ 2  2 , the ill-posed inverse problem becomes well-posed.Generally, the parameter  is chosen to be a fraction of the greatest eigenvalue of AA  or determined by L-curve method.For a fixed parameter , Tikhonov regularization always has an explicit analytic solution to formula (3).
where I is a unit diagonal matrix and  represents Hermitian transpose.The main limitation of the ESM with Tikhonov regularization is that the reconstructed sound pressure is much lower than the theoretical value, especially in middle to high frequency.
As mentioned before, the distribution of equivalent source amplitude is known to be sparse in the sense that most components of q are close to zero.The limitation encountered in Tikhonov regularization can be alleviated by replacing the  2 -norm in the penalty term by a  1 -norm.minimize q ‖q‖ 1 ; where ‖q‖ 1 denotes the  1 -norm and  denotes noise bound.Formula (5) can be written as a minimizer of a convex unconstrained objective function.minimize From formulas ( 5) and ( 6), one can find that the sparsity is enforced by use of  1 -norm penalty on the solution vector.The sparse solution of the objective function is typically obtained in an iterative manner as there is no analytic solution.Besides, the regularization parameter  in formula ( 6) is determined empirically.
Different from  1 -norm-based regularization, WBH algorithm uses a dedicated iterative solver including a dynamic filter to find the sparse solution of equivalent source amplitude.Firstly, the residual vector r is defined related to formula (1).
The residual quadratic function F to be minimized is defined as Denoting by q  the approximate solution after the th iteration step, the optimal solution of the residual function F(q) is computed with the steepest descent method.Besides, the least squares solution q 0 is used as the initial value of the iteration.The step length Δ(q  ) that minimizes F(q) in the steepest direction is computed first: Here,  is a constant vector related to the initial value q 0 , instead of using a relaxation factor in WBH.As the solution obtained with formula ( 4) is stable, the parameter  in formula (10) can modify the step length to avoid big fluctuation.g  and   are search direction and optimal search step length of the th iteration in the steepest decent method, respectively.
To suppress the ghosts, all components in q +1 below a certain threshold   are set to zero.Different from the filter process in WBH, a new filter function including natural logarithm is used here, in which the filter function is smoother than that in WBH.Although there is no proof of how formula (11) performs better, the new filter function can improve reconstruction accuracy of WBH in the simulations, particularly for noise sources at middle to high frequency.The threshold value in each step can be expressed as Here,   is updated during the iteration as shown in formula (13), and |q +1,max | is the amplitude of the largest element in q +1 .Accordingly, the approximate solution q +1 after ( + 1)th iteration step can be written as where the range of   is determined by initial value  0 , maximum value  max and the constant Δ.In formula (13),   → ∞ for  → ∞, it could lead the filter process to be invalid.In simulations, the parameter in formula ( 13) is found to have a significant influence on the iteration result.When comparing iWBH with WBH, the value of  0 ,  max , and Δ in WBH is set to be the same as in [21].For iWBH, the following values have been found to work very well: max = 40 dB; Suppose that final solution of F(q) is calculated after the th iteration step.
As shown in Figure 1, R is the transfer matrix between equivalent source surface and reconstruction surface.The sound pressure p  on the reconstruction surface can be calculated by Due to the fact that iterative algorithm can converge to the optimal solution, the resolution of WBH would reach a certain level when increasing iterative steps.To further improve the resolution, the output of iWBH p  is modified with a coefficient matrix obtained with FB.
Firstly, the cross spectrum of q  is defined as follows: where C is nonnegative and Hermitian matrix.The high order matrix function with the order 1/V is introduced, and it is applied to C by computing the spectral decomposition of C and applying the function to the eigenvalues: where U is a unitary matrix of C of dimension  × .Σ = diag( 1 , . . .,   ) is the diagonal eigenvalue matrix.The output of FB at th grid can be expressed as follows: with ‖ ⋅ ‖ 2 being the  2 -norm and R  the th row of the transfer matrix R. Note that the components in matrix Z are not the sound intensity on the reconstruction surface but the coefficients related to the cross spectral matrix of reconstruction sound intensity.
In order to explain how the FB algorithm focuses more on main lobe and suppress side lobe, a simple model is constructed as shown in Figure 2. In the model, there is only a single sound source and only a corresponding nonzero equivalent source.Under this condition, we define  as the singlesource amplitude and w 0 as propagation vector between the sound source and equivalent source.w 1 represents the propagation vector from the only equivalent source to side lobe, and w 2 denotes the propagation vector to the main lobe.
The relative level of the side lobe can be expressed by ratio between the amplitude of Z 1 and Z 2 .
where w 0 , w 1 , and w 2 are assumed to be unit vectors. 1 is the included angle of w 1 and w 0 , and  2 is the included angle between w 2 and w 0 .As the distribution of equivalent source amplitude and the sound intensity on reconstruction surface are almost consistent in terms of position, w 0 is approximately parallel to w 2 .Then,  in formula (20) can be simplified as a function that depends on  1 and  2 .For a fixed V, the reconstruction sound intensity in side lobe direction would decay faster than that in the main lobe direction as 0 < cos( 1 ) < cos( 2 ) = 1.Finally, the output of FB would focus more on the main lobe with the increase of V.However, the output of FB is just the coefficient matrix corresponding to the amplitude of sound pressure matrix produced by iWBH.The final output of IFESM is generated by modifying the iWBH output with Z.
In formula (21), every component in p  is modified with Z to suppress the side lobe, and the amplitude remains accurate due to the good reconstruction accuracy of iWBH.In other words, IFESM combines the advantages of iWBH and FB into one algorithm to produce better resolution than iWBH and conventional ESM.

Numerical Simulation
In this section, noise source maps obtained with the methods previously presented are compared.The aim is to highlight the abilities in terms of the reconstruction accuracy of iWBH and resolution of IFESM.The noise sources, measurement array, and equivalent source surface are arranged as shown in Figure 1.The simulation objects are pulsating ball sources with radius of 0.01 m, vibrating velocity is 2.5 × 10 −2 m/s, and the sound velocity is 340 m/s.In all simulations, the distance between measurement surface and origin of coordinates is 0.3 m, and the distance between reconstruction surface and origin of coordinates is 0.08 m.The equivalent source surface is located at 0.001 m from sound sources.Noise sources, measurement surface, reconstruction surface, and equivalent surface are in the same coordinates.The aperture of measurement surface is 0.1 m × 0.1 m including 8 × 8 measurement points.The mesh on reconstruction surface has 41 × 41 grids with 0.02 m spacing, and the mesh on the equivalent source surface is the same as reconstruction surface.In order to quantify the reconstruction accuracy of IWBH, the error is defined as where  true () and  cal () are, respectively, the theoretical pressure and reconstruction sound pressure at the th measurement point.

Simulation One.
The simulation one aims to analyze the reconstruction performance of iWBH for different type of sources.Figure 3 shows the reconstruction error of conventional ESM, WBH, and iWBH versus frequency.A single source is located at (0, 0, 0) m, and random noise is added to the complex measurement data at a level of 15 dB below the average sound pressure across the microphones.The reconstruction error is the average value of 20 simulations.Generally, the reconstruction accuracy of WBH and iWBH is better than conventional ESM.As conventional ESM focuses on relatively low frequency and the regularization parameter of Tikhonov method does not fit well in high frequency, the reconstruction error of conventional ESM increases rapidly versus frequency.Actually, differences between WBH and iWBH are the different filter function and modification on step length.The result confirms that iWBH performs better, especially in the frequencies higher than 600 Hz. Figure 4 presents sound pressure on reconstruction surface for two coherent noise sources at 800 Hz.The two sources with equal strength are located at (−0.2, 0, 0) m and (0.2, 0, 0) m, and random noise is also set to be 15 dB.As shown in Figure 4, WBH and iWBH have identified the peak of the sources which represent the source position and strength.It can be seen from Figure 4(d), in which a comparison along the middle row (the 21th row) is shown, that iWBH can identify the amplitude of reconstruction pressure more accurately than WBH.According to formula (22), the reconstruction error based on WBH and iWBH is −18.4 dB and −23.5 dB, respectively.Figure 5 shows the comparison of reconstruction sound pressures for two coherent noise sources at 1800 Hz, from which the same conclusion can be made that iWBH can reconstruct noise sources with better reconstruction accuracy.

Simulation Two.
To verify the better resolution of IFESM, localization results obtained with iWBH, WBH, and IFESM are compared in this section.As an important parameter, res represents the resolution of these three algorithms.It describes the ability of the algorithm to distinguish one source from another close to it and is given as res =  −15 dB  es × 100%.
Here,  −15 dB is the "diameter" of −15 dB down output of the algorithm and  es is the side length of the equivalent source surface.
Figure 6 shows the contour diagram of sound pressure on reconstruction surface.According to formula (23), res of different methods is shown in Table 1.The results show that IFESM can produce better resolution, and IFESM can further reduce the beam-width as the order V increases from 1 to 5.
Noise source maps of coherent noise sources located at (−0.2, 0, 0) m and (0.2, 0, 0) m are presented in Figure 7.The WBH noise source map is shown in Figure 7(a).Two main spots are apparent at the given positions, but strong side lobes prevent the actual detection of the source position.Figure 7(b) shows the iWBH noise source map.The noise source map obtained with iWBH exhibits two main spots at the known positions, with side lobes partially removed.The IFESM noise source map further reduces the level of side lobes and exhibits two main spots at the given positions with smaller beam-width.Figure 7(d) is IFESM noise source map with the order V set to 6, and the resolution is further improved compared with Figure 7(c) in which the order V is set to 2. Figure 8 shows the noise source maps of IFESM versus order V. Figures 8(a)-8(f) confirm that the resolution can be further improved as the order V increases appropriately.

Table 1 :
Comparison of  −15 dB and res.