Array optimization of sparse regularization equivalent source acoustic holography algorithm

2.1. Mathematical model of acoustic imaging algorithm

The logic of ESM is that a number of virtual point sound sources are set inside the sound source area and the sound field excited by those virtual point sound sources is superimposed to equivalently fit the original sound field. The schematic diagram of ESM is shown in Fig. 1.

Fig. 1ESM schematic diagram

Assuming that $N$ equivalent sources are set on the acoustic source surface, $q\left(r_v\right)$ represent the equivalent source intensity. The acoustic pressure signal $p\left(r_h\right)$ is collected by $M$ sensors on the array. According to the propagation law of sound wave in the free field, the acoustic pressure on the measured surface can be approximated as follows:

The $g\left(r_h,r_v\right)$ represents the free space Green's function from the equivalent source to the acoustic array. Generally, it can be rewritten into the form of matrix:

where ${\mathbf{Q}}_v$ is the equivalent source intensity vector, and the elements ${\mathbf{G}}_{hv}(m,n)$ in the matrix of ${\mathbf{G}}_{hv}$ represent the transfer function between the $n$-th equivalent source and the $m$-th element:

The $k$ is the wave number. For the solution of Eq. (2), obtaining the equivalent point source intensity is the key to the equivalent source method. However, in practical applications, in order to display the detailed information in the acoustic field more precisely, more equivalent point sources will be set up to improve the resolution, and the number of array sensors is expected to be as small as possible, that is $N>M$. This leads to the underdetermined problem of Eq. (2), whose solution subspace is not unique and needs to find its optimal solution.

The optimal solution of this problem is usually obtained by regularization. The traditional solution is to use Tikhonov regularization method to deal with, through the $l_2$-norm of the solution vector to constraint, obtain the smooth minimum energy solution, and transform the solution of the equation into the minimization problem of the following function:

The acoustic pressure values obtained by this regularization method are consistent at the measurement points, but generally smaller at other locations. In the actual environment, the location of the acoustic source in space is sparse, and most of the values of equivalent source intensity vector are zero. Therefore, it is considered to obtain the optimal solution with the goal of sparseness, that is, to replace $l_2$-norm with $l_0$-norm to optimize the solution of Eq. (2):

For Eq. (5), this is a non-convex optimization problem, which is generally difficult to solve. According to the relevant theories of compressed sensing, this problem can be extended to the solution of Eq. (6) under certain conditions, that is, it can be converted into a $l_1$-norm optimization problem. For the $l_1$-norm optimization problem, convex optimization algorithm can be used to solve the problem. The specific algorithms include base tracking algorithm or gradient shadow algorithm:

By calculating the intensity of equivalent source, the acoustic quantity at any position in space can be calculated by acoustic transfer function, and the acoustic field reconstruction can be realized.

2.2. Measurement of incoherence of the sensing matrix

If the solution result of Eq. (6) is consistent with the solution of Eq. (5), Restricted Isometry Property (RIP) and Mutual Incoherence Property (MIP) [15] should be satisfied, which are two different expressions to measure the uncorrelation of the sensing matrix.

Under the theoretical framework of compressed sensing, only the sensor matrix with unrelated columns can ensure the recovery of sparse signals, and the matrix ${\mathbf{G}}_{hv}$ in Eq. (2) is the sensor matrix. The MIP condition is that any two different column vectors in the requirement are not correlated, which is generally expressed by the coherence coefficient and calculated according to Eq. (7):

where ${\mathbf{a}}_i$, ${\mathbf{a}}_j$ represents any column vector in the matrix ${\mathbf{G}}_{hv}$, $\left({\mathbf{a}}_i,{\mathbf{a}}_j\right)$ represents the inner product. The smaller $\mu \left({\mathbf{G}}_{hv}\right)$ means the stronger incoherence of ${\mathbf{G}}_{hv}$, which is more conducive to accurately calculating the intensity of the equivalent sound source.

Unlike the MIP, which only considers the correlation of several pairs of column vectors, the RIP contains S-element column vectors, which is more suitable for the evaluation of matrix ${\mathbf{G}}_{hv}$ irrelevance. The $S$ represents the sparsity of the solution vector. It can be calculated from the following inequality:

The sensor matrix ${\mathbf{G}}_{hv}$ is substituted into Eq. (8) to obtain the minimum value of ${\delta }_S$. The $\mathbf{Q}$ is any S-sparse vector. The minimum value of ${\delta }_S$ is also known as the Restricted Isometry Constant (RIC), which can be used to measure the uncorrelation of the matrix ${\mathbf{G}}_{hv}$.

For Eq. (8), the ideal is that when ${\mathbf{G}}_{hv}$ is a unitary matrix, at this time ${‖\mathbf{G}\mathbf{Q}‖}_2^2={‖\mathbf{Q}‖}_2^2$, ${\delta }_S$ is zero. It also shows that the matrix ${\mathbf{G}}_{hv}$ is completely column independent. However, due to the fact $N>M$ that ${\delta }_S$ calculated from the actual constructed sensing matrix ${\mathbf{G}}_{hv}$ is greater than zero.

In practical, because of the large dimension of sparse vector, it is very time-consuming to accurately solve the value of ${\delta }_S$, which is difficult to achieve in reality. P. Simard [16] used Monte Carlo method and probability knowledge to find the statistical solution of ${\delta }_S$, so as to quickly compare the uncorrelation degree of array sensing matrix.

Monte Carlo method is a kind of stochastic simulation calculation method, which can obtain the statistical parameters of the problem by random sampling without complicated formula calculation, so as to infer the approximate value of the parameters. First, the Eq. (8) is transformed into the form of Eq. (9), and then the sparse vector is randomly generated, which is substituted into Eq. (9) and solved. The statistical value of $\delta$ is obtained to the frequency histogram, as shown in Fig. 2, which is approximately Gaussian distribution:

Fig. 2Statistical histogram of δ

According to the statistical knowledge, the statistical values is in the interval [$\stackrel{-}{\delta }-3\sigma$, $\stackrel{-}{\delta }+3\sigma$] with high probability. $\stackrel{-}{\delta }$ is the average and $\sigma$ is the standard deviation. Therefore, the ${\delta }_S$ can be expressed through this interval, as shown in Eq. (10):

Therefore, the statistical estimation ${\tilde{\delta }}_S$ can represent the property of measurement matrix to a certain extent, and provide a certain basis for the solution of acoustic source.

4.1. Array parameters

In order to study the influence of array arrangement on acoustic imaging results, the array size and array element number should be set in a consistent manner. The array size is set as 1 m×1 m, and the number of array elements is controlled at about 40. The cross array and network array are arranged according to the regular, and the optimized pseudo-random array is obtained by iterative screening of genetic algorithm. The layout diagram of the cross array, network array and optimized pseudo-random array are shown in Fig. 5. The programmed algorithm is used to process the sound pressure values of the sound field collected by these arrays respectively, and the acoustic imaging effects of these arrays are compared.

Fig. 5Array layout data

a) Cross array of 41 arrays

b) Grid array of 36 arrays

c) Optimize pseudo-random array of 40 arrays

4.2. Establishment of simulation model

The modeling and simulation work is completed by the acoustics software COMSOL Multiphysics, and the pressure acoustics module is invoked, the double-disk interference sound field model is shown in Fig. 6. The first step in the modeling is to construct a cubic air domain with 2 m sides and add a perfectly matched layer (PML) with a thickness of 0.2 m to its exterior. The second step is to place two disks with a radius of 0.05 m and a thickness of 0.01 m as sound sources, and the center distance between them is 0.4 m. The materials are set as structural steel, and the edges of the disks are fixed and constrained. And a force of 50 N is applied to the center of these disks to excite the sound field. The third step is to place the data acquisition array, which is 0.2 m away from these disks, to collect the sound pressure data in the sound field. After modeling, free tetrahedral grids are used to divide the entire study area, and the division is conducted according to the rule of no less than 5 grids in one wavelength. The research object is the frequency of sound sources. In the frequency range of 50-2000 Hz, the sound field generated by the force exerted on disks is simulated and calculated with the step size of 50 Hz.

According to the Monte Carlo method mentioned above, the estimated values of ${\tilde{\delta }}_S$ of different arrays are calculated, and the preliminary comparison of different arrays is made. The results are shown in Table 1.

From the numerical of ${\tilde{\delta }}_S$ point of view, the optimal pseudo-random array is the most suitable array for the acoustic holography algorithm, which satisfies the limit of 0.4625 proposed by Foucault, and meets the requirements of the RIP conditions.

Fig. 6Excitation sound field model of double disk

In order to compare the acoustic imaging effect of different arrays, the data collected from different arrays are used for acoustic field reconstruction. After processing by sparse regularization equivalent source algorithm, the acoustic images of different arrays are compared. As shown in Fig. 7-9, the acoustic pressure distribution at a distance of 0.1 m from the sound source is reconstructed at 50 Hz, 500 Hz, and 1000 Hz respectively. It can be clearly seen from the figure that the data collected by regular traditional array cannot accurately reconstruct the acoustic field acoustic pressure through the algorithm under different frequencies, which will misjudge the location and quantity of acoustic sources and cannot guide the location of acoustic sources; The data collected by the optimal pseudo-random array are processed by the same algorithm to get the acoustic pressure distribution closer to the actual acoustic field, which is in line with the actual acoustic pressure distribution on the whole, and shows the acoustic source location accurately without missing or misjudgment.

Fig. 7Acoustic pressure reconstruction of different arrays (f= 50 Hz)

a) Ideal measured value

b) Grid array of 36 arrays

c) Cross array of 41 arrays

d) Optimize pseudo-random array of 40 arrays

Fig. 8Acoustic pressure reconstruction of different arrays (f= 500 Hz)

a) Ideal measured value

b) Grid array of 36 arrays

c) Cross array of 41 arrays

d) Optimize pseudo-random array of 40 arrays

Fig. 9Acoustic pressure reconstruction of different arrays (f= 1000 Hz)

a) Ideal measured value

b) Grid array of 36 arrays

c) Cross array of 41 arrays

d) Optimize pseudo-random array of 40 arrays

In order to accurately describe the difference between the reconstructed acoustic pressure and the actual situation, the reconstruction error is generally calculated by Eq. (11):

where ${\mathbf{p}}_t$ is the reconstructed acoustic pressure value matrix at 0.1 m away from the acoustic source, and ${\mathbf{p}}_r$ represents the actual acoustic pressure value matrix at 0.1 m away from the acoustic source in the simulation environment. The reconstruction errors corresponding to different arrays are calculated, and Fig. 10 is drawn. The figure shows the variation of reconstruction error from low frequency to high frequency for different arrays.

It can be seen from the Fig. 10 that the reconstruction errors of cross array and grid array are larger, and the errors under different frequencies are quite different. The reconstruction error of the optimized pseudo-random array is smaller than that of the other two arrays, and it has a good reconstruction effect in a wide frequency range, which can be used to indicate the position of acoustic source.

Fig. 10Reconstruction error graph