A Sound Source Identification Algorithm Based on Bayesian Compressive Sensing and Equivalent Source Method

Abstract

Near-field acoustic holography (NAH) based on equivalent source method (ESM) is an effective method for identifying sound sources. Conventional ESM focuses on relatively low frequencies and cannot provide a satisfactory solution at high frequencies. So its improved method called wideband acoustic holography (WBH) has been proposed, which has high reconstruction accuracy at medium-to-high frequencies. However, it is less accurate for coherent sound sources at low frequencies. To improve the reconstruction accuracy of conventional ESM and WBH, a sound source identification algorithm based on Bayesian compressive sensing (BCS) and ESM is proposed. This method uses a hierarchical Laplace sparse prior probability distribution, and adaptively adjusts the regularization parameter, so that the energy is concentrated near the correct equivalent source. Referring to the function beamforming idea, the original algorithm with order v can improve its dynamic range, and then more accurate position information is obtained. Based on the simulation of irregular microphone array, comparisons with conventional ESM and WBH show that the proposed method is more accurate, suitable for a wider range of frequencies, and has better reconstruction performance for coherent sources. By increasing the order v, the coherent sources can be located accurately. Finally, the stability and reliability of the proposed method are verified by experiments.

1. Introduction

Acoustic holography is an advanced technology for sound source identification and sound field visualization. Williams and Maynard put forward the generalized theory of acoustic holography in 1980s [1]. The technique requires holographic data to be measured in the near-field of the sound source, so it is also called near-field acoustical holography (NAH). Since this method based on two-dimensional spatial Fourier transform is mainly applied to regular measurement surfaces, it has limitations in engineering applications.

To reduce the limitation of measurement conditions, in the past few decades various NAH methods have been developed, such as the boundary element method (BEM) [2,3], the statistically optimized near-field acoustic holography (SONAH) [4,5], and the equivalent source method (ESM) [6,7]. ESM, also known as wave superposition method (WSM), is an efficacious method for identifying sound sources. For ESM, it assumes that the sound field is a superposition of some equivalent sound sources and approximates the boundary conditions by placing several simple sound sources. The calculation error and robustness analysis results of ESM show that it is an effective method with high precision [6]. Conventional ESM usually uses $l_2$ norm regularization method, represented by Tikhonov regularization method (TRESM), to find the least energy solution in the solution subspace. It provides good resolution at low frequencies, but the reconstruction accuracy is low at mid-and-high frequencies. To overcome this problem, J. Hald proposed the wideband acoustic holography (WBH) method [8,9], which uses the steepest descent method combining dynamic filtering to obtain a sparse solution of equivalent source strength, with high reconstruction accuracy at mid-and-high frequencies. It has high computational efficiency and can be applied in irregular planar arrays. However, for coherent sound sources, the resolution of WBH at low frequencies is restricted due to the average spacing of array microphones.

In recent years, as an important tool for improving resolution, many sparse regularization methods have attracted the attention of scholars and been widely used in the field of image processing [10,11]. The sparse regularization method has also been applied to the field of sound source identification. Compressive sensing (CS) [12,13,14] based on sparse representation is a technique of signal acquiring and processing whose basic assumption is that data usually has sparse representations on certain bases and can be described with a few non-zero coefficients. In acoustics, it is often meaningful to obtain a sparse representation of the sound field. Chardon et al. [15] combined the sparsity and CS theory with NAH to identify the sound source with fewer microphones. Since the assumption of CS technology is the sparsity of the signal, the CS-based NAH can reconstruct high-resolution sound sources with less measured values. N. Chu et al. [16] combined aeroacoustic imaging with Bayesian compressive sensing (BCS) and made a comparison with the conventional beamforming method, such as DAMAS, CLEAN, etc., showing that the Bayesian method has a wider dynamic display range and it is more robust. In addition, there are some other methods [17,18,19,20,21] to prove the importance of regularization and CS framework.

The purpose of this paper is to further improve the reconstruction accuracy of WBH and conventional ESM methods, and an alternative ESM-based algorithm is proposed to reconstruct sound field information more accurately. Compared with the conventional ESM model, a near-field acoustic holography model based on ESM and BCS is established firstly. Then the independent Laplace prior is utilized for each equivalent source strength in a hierarchical manner. In the process of finding optimal sources, an effective greedy construction algorithm [14] is used to continuously update the regularization parameter in the iteration to broaden the equivalent source dynamic range. Since the parameters required in the model are completely estimated automatically during the calculation, the algorithm is fully automated. Finally, the equivalent source strength is estimated by the distribution of the equivalent source strength. For the sake of brevity, this BCS-based method using Laplace priors is referred to as LPBCS in the following article. In addition, referring to the function beamforming idea [22,23], the original algorithm is extended to the high order v form (LPBCS-v for brevity), and its dynamic display range is significantly improved, thereby obtaining more accurate information of positioning.

The outline of this paper is as follows. Firstly, the application of NAH and CS in the field of sound source identification are simply introduced. Then, Section 2 presents the basic theory of ESM, TRESM, and WBH. Section 3 presents the theoretical description of LPBCS and its extension form LPBCS-v. In Section 4, numerical simulations are performed to evaluate the proposed method by comparing it with the TRESM and WBH. Section 5 further validates the effectiveness of the proposed method through experiments. Finally, in Section 6 the conclusions are summarized.

2. Theoretical Background

2.1. Equivalent Source Model

The basic idea of the NAH based on ESM is that the sound field radiated by actual sound sources can be modeled by the superposition of a series of simple monopoles. The spatial sound field can be reconstructed by solving the strength of each equivalent source, as shown in Figure 1. The equivalent source surface is constructed near the actual sound sources. It assumes that the equivalent source radiates the sound pressure outward by the free-field radiation law. The sound pressure measured by the microphone array is used as the boundary condition to establish the equation. The equivalent source strength can be obtained by solving the inverse problem equation so that it is further possible to reconstruct the sound field information on the reconstruction surface. To simplify the calculation, the equivalent source and reconstruction surfaces are set to be rectangular planes. Assuming that there are H microphones on the measurement surface and L virtual equivalent sources on the equivalent source surface, the theoretical calculation model can be obtained

$$\left[\begin{array}{c}{p}_1\\ p_2\\ ⋮\\ p_H\end{array}\right]=\left[\begin{array}{ccccccc}{G}_{11}& & G_{12}& & \cdots & & G_{1L}\\ G_{21}& & G_{22}& & \cdots & & G_{2L}\\ ⋮& & ⋮& & \ddots & & ⋮\\ G_{H1}& & G_{H2}& & \cdots & & G_{HL}\end{array}\right]\left[\begin{array}{c}{q}_1\\ q_2\\ ⋮\\ q_L\end{array}\right].$$

(1)

Complex numbers $q_l$ and $p_h$ are the strength of the equivalent source and measured pressure, respectively. $G_{hl}$ can describe the connection between $q_l$ and $p_h$. It can be expressed as

$$G_{hl}=\frac{\exp (-\mathrm{j}k{∥{\mathit{r}}_h-{\mathit{r}}_l∥}_2)}{4\pi {∥{\mathit{r}}_h-{\mathit{r}}_l∥}_2},$$

(2)

where k is the wave number, ${\mathit{r}}_h$ and ${\mathit{r}}_l$ are the Cartesian coordinate of the h-th microphone and the l-th equivalent source, respectively. ${∥{\mathit{r}}_h-{\mathit{r}}_l∥}_2$ is the Euclidean distance of the vector difference. Equation (1) can be simplified as matrix-vector form

$$\mathit{p}=\mathit{G}\mathit{q}.$$

(3)

Here, $\mathit{G}$ is called the Green matrix. When the calculation of equivalent source strength is completed, the reconstruction surface sound pressure can be expressed as

$$\tilde{\mathit{p}}=\tilde{\mathit{G}}\mathit{q}.$$

(4)

Similar to the Green matrix $\mathit{G}$, $\tilde{\mathit{G}}$ represents the transfer matrix between the equivalent source surface grids and the reconstruction surface grids. $\tilde{\mathit{p}}$ indicates the reconstruction surface sound pressure. It can be seen that the accuracy of the calculated strength of equivalent sources will greatly affect the accuracy of the reconstructed sound pressure of the reconstruction surface, so the method used to calculate the equivalent source strength is crucial.

However, in general engineering applications, the number of measurement points is often much less than the number of virtual equivalent sources due to the limitation of measurement conditions. The problem is severely ill-conditioned and usually underdetermined and has a non-unique solution subspace. It is necessary to constrain the equivalent source strength vector so that it has a meaningful solution.

2.2. Tikhonov Regularization and Wideband Acoustical Holography

In general, limiting the solution vector $\mathit{q}$ requires regularization methods. The Tikhonov regularization method is one of the effective methods for calculating the optimal solution of the equivalent source strength. TRESM usually adds an $l_2$ norm penalty to the solution vector $\mathit{q}$ when minimizing the residual norm ${∥\mathit{p}-\mathit{G}\mathit{q}∥}_2^2$. Since the result given by the Tikhonov regularization method is only a least squares solution, some ghost sources appear near the real sound source. Therefore, it causes the sound source image to become blurred at the mid-and-high frequencies. To suppress the number of ghost sources, WBH uses the steepest descent iterative algorithm to solve the problem. In each iteration, the threshold filtering process is used to set elements whose values are smaller than the threshold value to zero, while retaining elements above the threshold value. At the same time, the threshold value decreases as the number of iterations increases, thereby reducing its own effect. The threshold filtering process can effectively reduce the number of equivalent sources, so that the energy is more concentrated near the true equivalent source. However, WBH can accurately identify the coherent sound sources only when the microphone spacing in the array is appropriate. It is recommended to identify the coherent sound sources above 0.7 times the frequency of a half wavelength average array microphones spacing.

Therefore, in this paper, we focus on developing a method to solve the ESM-based sound sources identification problem in a Bayesian framework, which can better identify the sound source. This method enables the coherent sound sources to be accurately identified at relatively low frequencies while maintaining high frequencies accuracy. In addition, taking advantage of the function beamforming properties, for the problem with low resolution of coherent sound sources, the cross-spectral matrix is used to greatly increase the dynamic display range by increasing the order, which makes the sound source location more accurate.

3. Reconstruction Algorithm

3.1. Compressive Sensing Model for Sound Field Reconstruction

In the compressive sensing problem, it assumes that for a signal $\mathit{d}\in {\mathit{R}}^N$, it can be represented by a set of orthogonal basis $\mathit{d}={\sum }_i{\phi }_i{x}_i=\Phi \mathit{x}$, where $\Phi$ is the basis matrix and $\mathit{x}\in {\mathit{R}}^N$ is a sparse vector, i.e., most of its coefficients are zero. Then the M-dimensional observation vector $\mathit{y}\in {\mathit{R}}^M$ is obtained by the projection matrix $\Psi$, i.e.$\mathit{y}=\Psi \mathit{d}=\Psi \Phi \mathit{x}=\mathit{A}\mathit{x}$. When the measurements are noisy, the form is as follows

$$\mathit{y}=\mathit{A}\mathit{x}+\mathit{n}.$$

(5)

Since $\mathit{x}$ is compressible, its reconstruction can be defined as a sparse regularization problem as follows

$$\widehat{\mathit{x}}=arg\underset{\mathit{x}}{min}\{{∥\mathit{y}-\mathit{A}\mathit{x}∥}_2^2+\rho {∥\mathit{x}∥}_1\},$$

(6)

where the term $\rho$ called regularization parameter is significant and it controls the relative importance applied to the Euclidian error and the sparsity of $\mathit{x}$.

In Bayesian framework, the sparse vector $\mathit{x}$ can be estimated from the CS-based measured value $\mathit{y}$. The posterior probability distribution function of the parameters is first calculated, and then the statistics in the distribution are used to estimate the specific required values.

In order to apply the Bayesian compressive sensing method to sound field reconstruction, the strategy is to convert the complex vector and matrix into real values, that is, to decompose the complex values into real and imaginary parts. Both parts have the same sparsity. Then the sound field reconstruction problem is the same as Equation (5) with

$$\mathit{y}=\left[\begin{array}{c}\mathrm{Re}\left(\mathit{p}\right)\\ \mathrm{Im}\left(\mathit{p}\right)\end{array}\right],$$

(7)

$$\mathit{x}=\left[\begin{array}{c}\mathrm{Re}\left(\mathit{q}\right)\\ \mathrm{Im}\left(\mathit{q}\right)\end{array}\right],$$

(8)

$$\mathit{A}=\left[\begin{array}{cc}\mathrm{Re}\left(\mathit{G}\right)& -\mathrm{Im}\left(\mathit{G}\right)\\ \mathrm{Im}\left(\mathit{G}\right)& \mathrm{Re}\left(\mathit{G}\right)\end{array}\right].$$

(9)

3.2. Reconstruction via Bayesian Method

In the reconstruction process, the observation model and the signal model are first introduced, and then sparse Bayesian inference is performed which includes a fast greedy algorithm.

In this the observation model, it assumes that the observed noise is independent and a Gaussian distribution with zero mean and the variance $s_0$. Therefore, the probability distribution of the measured vector $\mathit{y}$ is

$$\mathcal{P}\left(\mathit{y}\right|\mathit{x},s_0^{-1})=\mathcal{N}\left(\mathit{y}\right|\mathit{A}\mathit{x},s_0),$$

(10)

where $s_0^{-1}$ is the inverse-variance or precision of observed noise. Because the conjugate probability distribution of Gaussian distribution precision is Gamma distribution, there is

$$\mathcal{P}\left(s_0^{-1}\right|\alpha ,\beta )=\Gamma \left(s_0^{-1}\right|\alpha ,\beta ).$$

(11)

The conjugate prior makes the functional form of the posterior probability the same as the prior probability, which greatly simplifies the Bayesian analysis.

Then the signal model is introduced. Due to the sparsity of the equivalent sources, the unknown parameter $\mathit{x}$ is a sparse vector. According to the Bayesian method, a hierarchical prior is introduced due to the sparse prior probability of $\mathit{x}$. As the first layer of the hierarchical model, the Gaussian distribution is employed on $\mathit{x}$

$$\mathcal{P}\left(\mathit{x}\right|\mathit{s})=\prod _{i=1}^N\mathcal{N}\left(x_i\right|0,s_i).$$

(12)

As a second layer of the hierarchical model, the following Gamma distribution is used on the hyper-parameter $s_i$,

$$\mathcal{P}\left(s_i\right|\lambda )=\Gamma \left(s_i\right|1,\frac{\lambda }{2})=\frac{\lambda }{2}\exp \left(\frac{-\lambda s_i}{2}\right).$$

(13)

The overall prior on $\mathit{x}$ is evaluated as

$$\mathcal{P}\left(\mathit{x}\right|\lambda )=\int \mathcal{P}\left(\mathit{x}\right|\mathit{s})\mathcal{P}\left(\mathit{s}\right|\lambda )d\mathit{s}=\frac{{\lambda }^{N/2}}{2^N}\exp (-\sqrt{\lambda }\sum _{i=1}^N|x_i|).$$

(14)

In the last layer, as a hyper-parameter, $\lambda$ has the following hyper-prior distribution

$$\mathcal{P}\left(\lambda \right|\theta )=\Gamma \left(\lambda \right|\frac{\theta }{2},\frac{\theta }{2}).$$

(15)

The parameter selection of the prior distribution of $\lambda$ is very flexible. When $\theta$ approaches to zero, the distribution provides a very ambiguous message. The hierarchical prior described above has a multi-layer combination. The first two layers constitute the Laplace prior distribution, which has many advanced properties, such as heavy tails, the sharp peak at zero and the shape regulated by parameters. And the last layer can be used to calculate $\lambda$.

Here, Bayesian inference is performed by using type-II maximum likelihood approach. In the process of reconstruction, it is the posterior distributions of parameters that need to be solved. According to Bayesian inference, there is the posterior distribution

$$\mathcal{P}(\mathit{x},\mathit{s},\lambda ,s_0^{-1}|\mathit{y})=\frac{\mathcal{P}(\mathit{x},\mathit{s},\lambda ,s_0^{-1},\mathit{y})}{\mathcal{P}\left(\mathit{y}\right)}.$$

(16)

Here, the dependency on $\theta$ is dropped for conciseness. Due to the hierarchical prior structure, the joint distribution can be defined as

$$\mathcal{P}(\mathit{x},\mathit{s},\lambda ,s_0^{-1},\mathit{y})=\mathcal{P}\left(\mathit{y}\right|\mathit{x},s_0^{-1})\mathcal{P}\left(\mathit{x}\right|\mathit{s})\mathcal{P}\left(\mathit{s}\right|\lambda )\mathcal{P}\left(\lambda \right)\mathcal{P}\left(s_0^{-1}\right).$$

(17)

The denominator of the posterior formulation is not easy to solve because it requires complicated integration and has no analytical solution. So Equation (16) needs to be properly decomposed. The decomposition used in this inference procedure is as follows

$$\mathcal{P}\left(\mathit{x},\mathit{s},\lambda ,s_0^{-1}\right|\mathit{y})=\mathcal{P}\left(\mathit{x}\right|\mathit{y},\mathit{s},\lambda ,s_0^{-1})\mathcal{P}(\mathit{s},\lambda ,s_0^{-1}|\mathit{y}),$$

(18)

$\mathcal{P}\left(\mathit{x}\right|\mathit{y},\mathit{s},\lambda ,s_0^{-1})$ is the posterior probability distribution on $\mathit{x}$ which can be expressed as a Gaussian distribution $\mathcal{N}\left(\mathit{x}\right|\mu ,\Sigma )$ with mean and variance

$$\begin{array}{c}\mu =\Sigma s_0^{-1}{\mathit{A}}^{\mathrm{T}}\mathit{y},\end{array}$$

(19)

$$\begin{array}{c}\Sigma ={({\mathit{S}}^{-1}+s_0^{-1}{\mathit{A}}^{\mathrm{T}}\mathit{A})}^{-1},\end{array}$$

(20)

where $\mathit{S}=diag(s_1,s_2,\cdots ,s_N)$. Since $\mathcal{P}(\mathit{s},\lambda ,s_0^{-1}|\mathit{y})\propto \mathcal{P}(\mathit{y},\mathit{s},\lambda ,s_0^{-1})$, the hyper-parameters $\mathit{s}$, $\lambda$, $s_0^{-1}$ can be estimated by maximizing the joint distribution $\mathcal{P}(\mathit{y},\mathit{s},\lambda ,s_0^{-1})$, which can be obtained from Equation (17) by integrating out $\mathit{x}$

$$\mathcal{P}(\mathit{y},\mathit{s},\lambda ,s_0^{-1})={\left(2\pi \right)}^{-\frac{N}{2}}{\left|\mathit{C}\right|}^{-\frac{1}{2}}\exp \left(\frac{-{\mathit{y}}^{\mathrm{T}}{\mathit{C}}^{-1}\mathit{y}}{2}\right)\mathcal{P}\left(\mathit{s}\right|\lambda )\mathcal{P}\left(\lambda \right)\mathcal{P}\left(s_0^{-1}\right),$$

(21)

where $\mathit{C}=s_0\mathit{I}+{\mathit{A}\mathit{S}\mathit{A}}^{\mathrm{T}}$. To facilitate calculations, maximize the logarithm of Equation (21) equivalently

$$\begin{array}{cc}\hfill \mathcal{L}=& -\frac{N}{2}\log \left(2\pi \right)-\frac{\log \left|\mathit{C}\right|}{2}-\frac{{\mathit{y}}^{\mathrm{T}}{\mathit{C}}^{-1}\mathit{y}}{2}+N\log \frac{\lambda }{2}-\sum _{i=1}^N\frac{\lambda s_i}{2}+\frac{\theta }{2}\log \frac{\theta }{2}\hfill \\ \hfill & -\log \Gamma \left(\frac{\theta }{2}\right)+(\frac{\theta }{2}-1)\log \lambda -\frac{\theta \lambda }{2}-(\alpha -1)\log s_0-\frac{\beta }{s_0}.\hfill \end{array}$$

(22)

To reduce the computational complexity, a greedy construction manner can be used to decompose the $\mathit{C}$ into

$$\mathit{C}=s_0\mathit{I}+\mathit{A}\mathit{S}{\mathit{A}}^{\mathrm{T}}=s_0\mathit{I}+\sum _{j\ne i}{s}_j{\mathit{a}}_j{\mathit{a}}_j^{\mathrm{T}}+s_i{\mathit{a}}_i{\mathit{a}}_i^{\mathrm{T}}={\mathit{C}}_{-i}+s_i{\mathit{a}}_i{\mathit{a}}_i^{\mathrm{T}}.$$

(23)

where ${\mathit{a}}_i$ is the i-th column basis in $\mathit{A}$, and ${\mathit{C}}_{-i}$ is $\mathit{C}$ with ${\mathit{a}}_i$ removed. By using wood-bury identity and determinant identity, ${\mathit{C}}^{-1}$ and $\left|\mathit{C}\right|$ can be obtained

$$\begin{array}{c}{\mathit{C}}^{-1}={\mathit{C}}_{-i}^{-1}-\frac{{\mathit{C}}_{-i}^{-1}{\mathit{a}}_i{\mathit{a}}_i^{\mathrm{T}}{\mathit{C}}_{-i}^{-1}}{s_i^{-1}+{\mathit{a}}_i^{\mathrm{T}}{\mathit{C}}_{-i}^{-1}{\mathit{a}}_i},\end{array}$$

(24)

$$\begin{array}{c}\left|\mathit{C}\right|=|{\mathit{C}}_{-i}\left|\right|1+s_i{\mathit{a}}_i^{\mathrm{T}}{\mathit{C}}_{-i}^{-1}{\mathit{a}}_i|.\end{array}$$

(25)

Then $\mathcal{L}$ can be treated as a function of $\mathit{s}$ only

$$\begin{array}{cc}\hfill \mathcal{L}\left(\mathit{s}\right)& =-\frac{1}{2}(N\log \left(2\pi \right)+\log \left|{\mathit{C}}_{-i}\right|+{\mathit{y}}^{\mathrm{T}}{\mathit{C}}_{-i}^{-1}\mathit{y}+\sum _{j\ne i}\frac{\lambda s_j}{2})+\frac{1}{2}(\log \frac{1}{1+s_i{f}_i}+\frac{g_i^2{s}_i}{1+s_i{f}_i}-\lambda s_i)\hfill \\ \hfill & =\mathcal{L}\left({\mathit{s}}_{-i}\right)+l\left(s_i\right).\hfill \end{array}$$

(26)

where $f_i$ and $g_i$ are defined as

$$\begin{array}{c}{f}_i={\mathit{a}}_i^{\mathrm{T}}{\mathit{C}}_{-i}^{-1}{\mathit{a}}_i,\end{array}$$

(27)

$$\begin{array}{c}{g}_i={\mathit{a}}_i^{\mathrm{T}}{\mathit{C}}_{-i}^{-1}\mathit{y}.\end{array}$$

(28)

The term $l\left(s_i\right)$ containing parameter $s_i$ in Equation (26) is separated from others so that the derivation with respect to $s_i$ is simplified. By setting the derivation to zero, its solution can be expressed analytically as

$$s_i=\left\{\begin{array}{cc}\frac{-f_i(f_i+2\lambda )+f_i\sqrt{{(f_i+2\lambda )}^2-4\lambda (f_i-g_i^2+\lambda ))}}{2\lambda f_i^2},\hfill & \mathrm{if}{g}_i^2-f_i>\lambda \hfill \\ 0,\hfill & \mathrm{if}{g}_i^2-f_i\le \lambda .\hfill \end{array}\right.$$

(29)

The updates of other hyper-parameters are similar. By taking the derivative of Equation (26) with respect to $\lambda$ and making it equal to zero, the update formula of the parameter $\lambda$ can be obtained

$$\lambda =\frac{2N-2+\theta }{{\sum }_i{s}_i+\theta }.$$

(30)

Similarly, an estimate of the noise variance can be written as

$$s_0=\frac{{∥\mathit{y}-\mathit{A}\mu ∥}_2^2+2\beta }{N+2\alpha }.$$

(31)

The regularization parameter $\rho$ in Equation (6) is driven by the parameters $\lambda$ and $s_0$ [13], so the iteration process of $\lambda$ and $s_0$ is equivalent to iteratively updating the regularization parameter $\rho$ which can make the reconstruction of sound field more accurate. This is an important reason to introduce a third layer of the hierarchical prior to overcome the drawback of using a fixed regularization parameter in conventional ESM.

So far, through the above Bayesian inference, the updated formulas of the necessary hyper-parameters have been obtained. The parameter $\mu$ is the mean of the posterior distribution of $\mathit{x}$. It is used to estimate the vector $\mathit{x}$. Then the equivalent source strength $\mathit{q}$ can be obtained from Equation (8). The LPBCS algorithm is summarized as follows:

• Input $\mathit{y}$ and $\mathit{A}$ by Equation (7) and Equation (9).
• Initialize $s_0=0.1\mathrm{var}\left(\mathit{y}\right)$, $\mathit{s}=\mathbf{0}$, $\lambda =0$.
• Select a basis vector ${\mathit{a}}_i$ out of $\mathit{A}$ and update the basis vector,
  if $g_i^2-f_i>\lambda$ and $s_i>0$, ${\mathit{a}}_i$ is in the model, re-estimate $s_i$,
  if $g_i^2-f_i>\lambda$ and $s_i=0$, add ${\mathit{a}}_i$ to the model with updated $s_i$,
  if $g_i^2-f_i\le \lambda$, prune ${\mathit{a}}_i$ from the model and set $s_i=0$.
• Update parameters $\mu$, $\Sigma$, $g_i$, $f_i$, $\lambda$.
• The iteration terminates when the value of marginal likelihood changes less than the set threshold or when the maximum number of iterations is exceeded. Otherwise go to step 3.
• Set $\mathit{x}=\mu$, then estimate the vector $\mathit{q}$ by Equation (8).

3.3. High Dynamic Range for $\mathrm{LPBCS}$-v

To locate the sound source more accurately, the dynamic range is improved. Inspired by the function beamforming, the solution of Equation (4) is modified, and the output at the field point t on the reconstruction surface can be represented by the cross-spectral matrix of $\mathit{q}$,

$${\tilde{P}}_t={\tilde{p}}_t{\tilde{p}}_t^{\mathrm{H}}={\tilde{\mathit{G}}}_t\mathit{q}{\mathit{q}}^{\mathrm{H}}{\tilde{\mathit{G}}}_t^{\mathrm{H}}={\tilde{\mathit{G}}}_t\mathit{C}{\tilde{\mathit{G}}}_t^{\mathrm{H}},$$

(32)

where ${\tilde{\mathit{G}}}_t$ is the t-th row vector of $\tilde{\mathit{G}}$, $\mathit{C}=\mathit{q}{\mathit{q}}^{\mathrm{H}}$. The eigendecomposition of $\mathit{C}$ can be written as

$$\mathit{C}=\mathit{U}\mathsf{\Lambda }{\mathit{U}}^{\mathrm{H}},$$

(33)

where $\mathit{U}$ is a unitary matrix which represents the set of eigenvectors of $\mathit{C}$, $\mathsf{\Lambda }$ is a diagonal matrix $diag({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_N)$ whose diagonal elements show the eigenvalues. An exponential matrix function of $\mathit{C}$ is defined by

$${\mathit{C}}^{1/v}=\mathit{U}diag({\gamma }_1^{1/v},{\gamma }_2^{1/v},\cdots ,{\gamma }_N^{1/v}){\mathit{U}}^{\mathrm{H}},$$

(34)

where $1/v$ is the power of the eigenvalue. Then replace $\mathit{C}$ in Equation (32) with ${\mathit{C}}^{1/v}$ and take the v power of the output

$${\tilde{P}}_t^{\prime }={\left({\tilde{\mathit{G}}}_t{\mathit{C}}^{1/v}{\tilde{\mathit{G}}}_t^{\mathrm{H}}\right)}^v.$$

(35)

Finally, the output can be written in vector form ${\tilde{\mathit{P}}}_v={[{{\tilde{P}}^{\prime }}_1,\cdots ,{{\tilde{P}}^{\prime }}_t,\cdots ,{{\tilde{P}}^{\prime }}_N]}^{\mathrm{T}}$. By applying the order v, the output strength of the sound source may be smaller than the true value, so the amplitude is corrected to ensure the true estimation of the maximum pressure amplitude

$${\tilde{\mathit{P}}}_{out}={\tilde{\mathit{P}}}_v\frac{\max \left(\right|\tilde{\mathit{p}}\left|\right)}{\max \left(\right|{\tilde{\mathit{P}}}_v\left|\right)}.$$

(36)

This is the final output of LPBCS-v. The order $v=0$ is defined as the original output. If we want to ensure a near-perfect reconstruction of the sound field, we can apply the result when the order is zero. However, if we want to identify the sound source position with high precision and large dynamic range, the order is adjusted to 1, 2, 3, etc. Therefore, the selection of the order v is very flexible and can provide high-precision sound source localization.

4. Numerical Simulations

To illustrate the performance of the proposed method for sound source identification, numerical simulations are implemented. The specific settings of the simulation are as follows. In the same coordinate system, the distance between the measurement surface and the origin of the coordinate is 0.2 m, and the distances from the equivalent source surface and the reconstruction surface to the origin are 0.001 m and 0.02 m, respectively. These surfaces are parallel to each other and consistent with Figure 1. Assume that each sound source is a pulsating sphere with a radius of 0.01 m and vibrating velocity of $2.5\times {10}^{-2}$ m/s. The speed of sound is 340 m/s. Reconstruction surface and the equivalent source surface are both meshed into $41\times 41$ points with the grid spacing of 0.02 m. To ensure consistency between simulations and experimentations, a Combo microphone array is used for simulation measurements. The array has 18 microphones distributed in a circle, and the microphone coordinates of the measurement array are shown in Figure 2.

The influences of different factors on the reconstruction accuracy are studied under two simulation conditions. In condition one, a pulsating sphere is placed at (0, 0, 0) m, which is the middle of the sound source surface as shown in Figure 3a, and the simulation results are presented in Figure 4, Figure 5, Figure 6 and Figure 7. In condition two, to study the identification of coherent sound sources, two coherent pulsating spheres are placed at coordinates (0.2, 0, 0) m and (−0.2, 0, 0) m, respectively, and the distance D = 0.4 m. The specific settings are as shown in Figure 3b and the simulation results are shown in Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12. To quantify the reconstruction accuracy of the sound pressure on the reconstruction surface, the percentage of reconstruction error is defined as

$$\mathrm{RE}=\frac{{∥\left(\right|{\tilde{\mathit{p}}}_{cal}|-|{\tilde{\mathit{p}}}_{true}\left|\right)∥}_2}{{∥\left(\right|{\tilde{\mathit{p}}}_{true}\left|\right)∥}_2}\times 100\%,$$

(37)

where ${\tilde{\mathit{p}}}_{cal}$ is the calculated sound pressure of the reconstruction surface and ${\tilde{\mathit{p}}}_{true}$ is the theoretical pressure. The lower the reconstruction error is, the higher the reconstruction accuracy is.

4.1. Simulation for Single Sound Source

The purpose of this part is to compare the reconstruction performance of LPBCS with the conventional equivalent source method TRESM and the sparse method WBH. The reconstructed sound pressure along the middle row of reconstruction surface and two-dimensional pressure distribution on the reconstruction surface are shown in Figure 4 and Figure 5. Figure 4 mainly compares the accuracy of sound pressure reconstruction, while Figure 5 is mainly used to observe the accuracy of positioning. The frequencies of comparison are 400, 800, 1800, and 3000 Hz with a signal-to-noise ratio (SNR) of 20 dB.

Figure 4 shows that the sound pressure calculated by TRESM at each frequency is much lower than the theoretical sound pressure, which means that the reconstruction accuracy of sound pressure is poor. The sound pressure of WBH has a certain gap with the theoretical sound pressure at low frequencies, and gradually approaches the theoretical sound pressure as the frequency increases. The pressure calculated by LPBCS is similar to the theoretical sound at each frequency, so its reconstruction accuracy of sound pressure is better than other two methods. The purpose of LPBCS-v is not to completely reconstruct the sound field but to locate the sound source. LPBCS-v can reduce the main lobe width so it can provide a more accurate position than LPBCS.

As Figure 5 shows, the source map of TRESM has a large hot-spot at low frequency, such as 400 and 800 Hz, so the resolution is poor. However, the positioning does not fail. As the frequency increases, the hot-spot becomes smaller but more ghost sources appear, which affect its positioning performance to some extent. The source map of WBH has a slightly larger hot-spot than the result of LPBCS at 400 Hz, so the resolution is slightly worse than LPBCS. Its positioning effect is close to the LPBCS at 800, 1800, and 3000 Hz. These show that WBH and LPBCS have similar effects, and both of them are better than TRESM on reconstruction accuracy and positioning performance. When LPBCS-v is used, the main lobe width is reduced while the sound pressure peak is guaranteed. The resolution is greatly improved and the positioning is more accurate.

Figure 6 shows the reconstruction error curves in the frequency range of 400–4000 Hz under different measurement noise levels (i.e., 10, 20, and 30 dB). The reconstruction error for each frequency is calculated 20 times to obtain an average reconstruction error. The results show that the TRESM reconstruction error is large and increases rapidly with the increase of frequency. This is because the number of observations in the high frequency is seriously insufficient and therefore sparse promotion method is needed to enhance the reconstruction accuracy. Obviously, under different measurement noise levels, the reconstruction accuracy of WBH and LPBCS in the frequency band is higher than TRESM. While LPBCS is slightly better than the WBH method at low frequencies, the two methods are similar at high frequencies. Figure 7 shows the reconstruction error curves with different measurement distances (i.e., 0.1, 0.2, and 0.3 m). All methods almost maintain their trends of error curves at different measurement distances. Moreover, the results show that LPBCS can maintain a better reconstruction accuracy and maintain its low-frequency accuracy advantage under different measurement distances.

4.2. Simulation for Coherent Sound Sources

It is necessary to distinguish the pressure generated by the coherent sound sources. Figure 8 shows the reconstructed pressure of coherent sound sources along the middle row of reconstruction surface. Figure 9 and Figure 10 show the two-dimensional pressure distribution of the coherent sound sources. In order to strengthen the research on the reconstruction accuracy and localization accuracy of coherent sound sources, more frequencies are used for comparison. The frequencies of comparison are 400, 600, 800, 1000, 1200, 1800, 3000, and 4000 Hz.

It can be seen from Figure 8 that the sound pressure calculated by TRESM is significantly different from the theoretical sound pressure at each frequency, and the sound pressure reconstruction results are very poor. For WBH, the reconstructed and theoretical sound pressure are also different from 400 to 1000 Hz. When the frequency is increased to 1200 Hz, the reconstructed sound pressure of WBH is close to the theoretical value, but the amplitude is biased. At 3000 and 4000 Hz, the sound pressure reconstructed by WBH is very close to theoretical sound pressure. Compared with WBH and TRESM, LPBCS has better reconstruction accuracy at the corresponding frequencies. Only LPBCS can reconstruct two main lobes at a low frequency of 400 Hz. At 600 and 800 Hz, the reconstructed sound pressure of LPBCS is close to the theoretical value, although the sound pressure of the sound source is underestimated. It is consistent with the theoretical sound pressure from 1000 to 4000 Hz. For LPBCS-v, it cannot reconstruct the true sound pressure but can significantly reduce the width of two main lobes. Hence, it can provide more accurate positions with a high dynamic range.

Figure 9 and Figure 10 present the positioning performance of coherent sources more clearly. TRESM is unable to locate the sound sources at 400, 600, and 800 Hz, and the two sound sources are completely replaced by several larger hot spots. The positioning performance of TRESM is improved when the frequency of the sound source is set at 1000, 1200, and 1800 Hz, but the spacing of the two sound sources is underestimated. The sound sources positions can be located by TRESM at 3000 Hz, although there is a certain gap between the theoretical and reconstructed sound positions. When the frequency of sound sources is set at 4000 Hz, a large number of ghosts seriously affect the positioning accuracy of TRESM. For WBH, the positioning fails at 400, 600, and 800 Hz. To avoid these situations, WBH is recommended to be used above the 0.7 times the frequency of a half wavelength average inter-element spacing. At 1000 Hz, although there are some offsets in the positions of the two sources, the positions of the sound sources can already be identified by WBH. The sources can be accurately positioned by WBH from 1200 to 4000 Hz. For LPBCS, although there is a positioning deviation at the low frequency of 400 Hz, the positioning does not fail. From 800 to 4000 Hz, this method can accurately locate sound source positions, which are consistent with the theoretical sound source positions. By using LPBCS-v, the resolution is enhanced, and the positioning performance is significantly improved. It can be seen that the proposed method can significantly improve the positioning accuracy at low frequencies and is more applicable in a wider frequency band.

The reconstruction error curves of the coherent sound sources under different measurement noise levels are shown in Figure 11. TRESM has a large error in the frequency band shown. The sound pressure reconstruction accuracy is poor, and the reconstruction error increases versus frequency. For each different noise level, WBH has a large reconstruction error at low frequencies, but a better reconstruction accuracy at middle and high frequencies. The LPBCS reconstruction error is less than TRESM and WBH throughout the displayed frequency band under different measurement noise levels. Notably, the limitation of lower frequencies accuracy of the LPBCS is smaller than that of conventional WBH. The results show LPBCS reconstructs the sound pressure of the coherent sound sources more accurately than the conventional method, which makes the positioning of the sound sources more accurate. Figure 12 shows the reconstruction error curves with different measurement distances. It can be seen that as the distance changes, the reconstruction error curves of LPBCS can maintain a relatively stable state. The accuracy is better than the traditional ESM and WBH. The advantage of its low-frequency reconstruction accuracy is more obvious. However, although the low-frequency accuracy of LPBCS is better than the WBH method under a large noise level, its error curves fluctuate. There is a certain gap between its low-frequency accuracy and high-frequency accuracy. It can be seen that the improvement of low-frequency accuracy with the large background noise is still an open and challenging issue.

5. Experimental Validation

To assess the applicability of the LPBCS, experimental verification is carried out in this section. Three methods are compared by a couple of experiments, including the single sound source and coherent sound sources. The experimental equipment includes an 18-channel Combo microphone array and a B&K LAN-XI data collector module. Consistent with the simulation, the distance between the speakers and the microphone array is 0.2 m. The locations of the loudspeakers are as shown in Figure 13. When performing a single source test, the loudspeaker is placed in the center of the sound source surface. For coherent sound sources, the two loudspeakers are symmetrically arranged on the sound source plane with an interval of 0.4 m. The in-phase sine waves are excited by loudspeakers. The sampling frequency of the collector is set to 16,384 Hz and the sampling time is 3 s. The meshing and distance settings for the equivalent source and reconstruction surfaces are the same as the simulation. The pressure distribution maps of reconstruction surface are presented in Figure 14, Figure 15 and Figure 16. Their frequencies correspond to the simulation.

As can be seen from Figure 14, for the single source, all methods can identify the location of the loudspeaker. However, at 400 and 800 Hz, the hot-spots identified by the TRESM are larger. At 1800 and 3000 Hz, for TRESM, the range of hot-spots is reduced but the accompanying ghost sources increase. While the maps of WBH have slightly smaller hot-spots at the frequencies shown and provide higher frequency resolution as the frequency increases. Overall, LPBCS maintains a good frequency resolution at the indicated frequencies compared with the other two methods, and its positioning performance is better. The LPBCS-v can identify the location more accurately, which is in agreement with the simulation.

Figure 15 and Figure 16 correspond to the recognition result of the coherent sound sources. It can be seen that multiple hot-spots identified by TRESM cannot locate the sound source at low frequencies. A large number of ghosts appear at 3000 and 4000 Hz, which seriously affect the positioning accuracy of TRESM. WBH fails to identify the locations at 400, 600, and 800 Hz. As frequency increases, its positioning performance becomes better. The sources can be accurately positioned by WBH from 1200 to 4000 Hz. It is noteworthy that only LPBCS can locate the coherent sound sources at 400 Hz. Although the peak positions deviate from the actual speaker positions, the deviation is small. LPBCS is able to pinpoint the sound sources at other higher frequencies. Meanwhile, LPBCS-v can further improve the frequency resolution and accurately locate the sound sources.

Therefore, similar to the simulation results, LPBCS can effectively identify the coherent sound sources with high accuracy over a wider frequency range. Compared with conventional methods, LPBCS has obvious advantages. However, when the order $v>0$, if there is a relative energy difference between the two sound sources, the smaller strength source may be covered by the larger strength source. In actual measurements, in order to obtain the maximum dynamic range, the order v can be gradually increased to achieve the purpose of accurately identifying the sound source location. In other words, as v gradually increases, the range of hot-spots will gradually decrease, which can make the positioning range smaller. So when our goal is to reconstruct the sound field, set v to 0, and when we are more interested in the specific location of the sound source, we should increase the order appropriately.

6. Conclusions

A sound field reconstruction algorithm LPBCS is suggested to solve the inverse problem based on ESM over a wide frequency range from a Bayesian perspective. The proposed method uses a hierarchical form of Laplace prior on each equivalent sound source and the regularization parameter can change with iterations. The results of numerical simulations show that LPBCS is superior to TRESM in terms of reconstruction accuracy and positioning performance. Compared with WBH, LPBCS has higher precision in coherent sound sources identification. The robustness of the proposed method is verified by changing the frequency, source measurement distance, and noise level. By combining the advantage of function beamforming with LPBCS, LPBCS-v can accurately locate the equivalent sources with a high dynamic range. The experimental results are in agreement with the simulation results, which prove the stability and reliability of the proposed method. Valuable and meaningful areas of future research may include the improvement of low-frequency accuracy under large background noise and the development of CS technical for better accuracy.