Near-field acoustic holography for high-frequency weak sound sources under low signal-to-noise ratio

2.1. Plane equivalent source method near-field acoustic radiation calculation model

Setting the virtual sound sources on a plane, which has better convenience and versatility [16]. As shown in Fig. 1, $S_h$ represents the measurement plane, $S_s$ represents the sound source plane, $S_v$ represents the virtual sound source plane, $n$ represents the normal direction of the sound source outside the plane, $P_h$ is the sound pressure distribution on the measurement plane, and $P_s$ is the sound source plane Sound pressure distribution, $Q_v$ is the sound intensity of the virtual sound source.

Fig. 1A schematic diagram of plane ESM

The sound pressure distribution on the measurement surface can be represented by the following integral equation:

where $r_h$ is the position vector of a point on the measurement plane, $r_v$ is the position vector of a point on the virtual sound source plane, $g\left(r_h,r_v\right)$ is the free space Green's function from the virtual point source to the measurement point, $q\left(r_v\right)$ is the virtual point sound source intensity, $\rho$ is the medium density, $c$ is the sound speed, and $k$ is the wave number.

The sound pressure on the actual sound source surface can be expressed in the following form after being discretized:

where $P_s$ is the sound source surface acoustic pressure column vector, $Q_v$ is the source strong column vector of the equivalent source sequence, and $G_{sv}$ is the sound pressure matching matrix between the equivalent source sequence and the sound source surface:

where $r_s$ is the position vector of a point on the sound source plane.

Similarly, after the sound pressure on the measurement surface is discretized, it can be expressed as:

where $P_h$ is the sound pressure column vector on the measurement surface, $Q_v$ is the source strong column vector of the equivalent source sequence, and $G_{hv}$ is the sound pressure matching matrix between the equivalent source sequence and the measurement surface:

The inversion formula for the virtual source strength is:

Further, the sound pressure of the actual sound source surface can be expressed as:

The sound pressure transfer matrix from the sound source surface to the measurement surface is expressed as:

Get the sound pressure $P_S$ on the sound source surface:

Knowing $L$ and $P_h$, the sound pressure distribution on any plane in the sound field can be reconstructed.

2.2. Array signal weight calculation

The definition of SNR is generally two kinds: one uses the energy or effective power of the original signal and the noise signal; the other is to use the effective values of the sound pressure of the original signal and the noise signal.

When the SNR of the sound array sampling signal is extremely low, not only the distribution of the sound field cannot be effectively obtained through the conventional NAH methods, but even misleading. In actual situations, the energy or sound pressure effective value of the original signal and the noise signal is often unknown, the SNR of the signal cannot be obtained, and the required weight value cannot be obtained. The EMD is used to analyze the signal obtained from each element sampling. The key to this method is the empirical mode decomposition. It can decompose the complex signal into a finite number of intrinsic mode functions (IMFs). Each IMF component contains characteristic signal of different time scales of the original signal. The basic idea can be summarized as: original waveform = $\sum IMFs$ + after wave.

Take the example of the high-frequency acoustic signal generated by the gas leakage in an open building. The noise is generated by the motor gearbox pump system, the air compressor and the lathe operation. Perform spectrum analysis on the above audio signals to obtain the following figures.

Fig. 2Frequency domain diagram of each signal

a) Gas leakage

b) Motor gearbox pump system

c) Air compressor

d) Lathe

When the above audio signals are mixed and sampled, the gas leakage sound is assumed to be the original signal, and the audio signals generated by the motor gearbox pump system, the air compressor and the lathe can be regarded as the noises. The mixed signal is analyzed to obtain the following frequency domain diagram.

Fig. 3Frequency domain diagram of mixed signal

It can be seen from each frequency domain diagram that only the frequency of the main component of the noise generated by the gas leakage is higher than 5 kHz, so the high-frequency noise in the mixed noise greater than 5 kHz can be regarded as the noise generated by the gas leakage. The high frequency part only occupies a part of the mixed signal. If the mixed signal is not preprocessed, the signal obtained by the array element measurement will have a lower SNR. Therefore, the signal can be pre-processed and multiplied by a weight $\alpha$, which takes a value between 0 and 1, which is proportional to the proportion of the high-frequency portion in this signal. This process can improve the resolving power of the array and improve the positioning accuracy of high frequency weak sound sources [17].

Using the EMD for this mixed signal, a limited number of IMF components are obtained. The result is shown in the following diagrams.

Fig. 4The time-frequency diagrams of IMF components

a) IMF1

b) IMF2

c) IMF3

d) IMF4

e) IMF5

f) IMF6

g) IMF7

h) IMF8

i) IMF9

j) IMF10

k) IMF11

l) IMF12

From the time-frequency diagrams of each IMF component, it can be known that the high-frequency part of the mixed signal ($f>$ 5 kHz) is mainly distributed in IMF1 and IMF2. In the frequency domain diagram, the larger the ordinate value, the greater the contribution of the signal portion representing this frequency in the component. The amplitude corresponding to each frequency point on the frequency domain diagram corresponding to each IMF component is accumulated, and a total of 12 values are obtained, and then it is determined whether the corresponding frequency of the main component of each component is greater than 5 kHz. In this signal, the main component of IMF1 and IMF2 corresponds to a frequency greater than 5 kHz. The 12 values obtained are accumulated as denominators, and the corresponding values of IMF1 and IMF2 are accumulated as molecules. The weight value required for the mixed signal is obtained, $\alpha =$ 0.5529, and the sound pressure value of the high frequency portion of the signal is obtained.

2.3. Sparse representation model of acoustic pressure array signals

Fig. 5 shows the microphone array signal sampling model. On the measurement plane $S_h$, there is a regular array of M-element microphones. The number of sound source points is usually limited, so the sound source points are sparsely distributed on the sound source surface $S_s$. The sound source surface $S_s$ can be evenly divided into grid points so that the sound source points can be sparsely represented on the sound source surface [18].

As shown in Fig. 6, the sound source surface $S_s$ is divided into $R$ areas. A solid area indicates that there is a sound source at that location, and a hollow area indicates that there is no sound source at that location.

The sound pressure signals received by the regularly distributed M-element microphone array are sparsely represented as:

In the formula, $v$ is an $M$×$1$ dimensional matrix, which represents the signal received by the measuring surface microphone array. $A$ is an $M$×$R$ dimensional matrix that represents the sensing matrix. $u$ is an $R$×1 dimensional matrix and represents a sound source signal containing position information, of which there are only $N\left(N\ll R\right)$ non-zero values. $\theta$ is the received noise signal.

Fig. 5Microphone array signal sampling model

Fig. 6Sound source surface sparse representation

2.4. A plane equivalent source near-field acoustic holography combining compressive sensing

The actual number and distribution of sound sources are often sparse relative to the space plane to be measured, which is in accordance with the requirements for signal sparsity in compressive sensing theory. Choosing the appropriate observation matrix can achieve sound source location under the condition of few measurement values [19]. In the planar equivalent source near-field acoustic radiation model, the virtual sound source surface $S_v$ is discretized and is similar to the sparse representation model of the sound pressure matrix signal, so the two can be combined. The virtual sound source points are sparsely distributed on the virtual sound source plane. The transfer matrix in the plane ESM is directly used as a sensing matrix, and the matrix conforms to the Restricted Isometry Property (RIP) [20]. Firstly, the position of the virtual sound source on the virtual sound source surface is obtained by orthogonal matching pursuit (OMP) algorithm based on compressive sensing theory, and then the holographic surface sound pressure distribution is reconstructed in combination with the known transfer matrix. Due to the inconvenient calculation of the area of the small polygons of the non-uniform grid, a uniform grid is generally taken. In the plane rectangular coordinate system, the virtual sound source plane may be discretized into $I\times I$ uniform squares, and the microphones on the measurement surface may also be evenly spaced to facilitate calculation.

The sound pressure transfer matrix $L$ from the virtual sound source plane to the holographic plane has been given by Eq. (8). This is the sensing matrix. Substituting it into Eq. (10) yields:

The presence of a noise signal will affect the accuracy of sound source localization, where $\theta$ is the noise signal, which is directly recorded by the measurement array, so Eq. (11) can be written as:

In the flow of the OMP reconstruction algorithm [21, 22]: $r_t$ indicates residual, $t$ indicates the number of iterations, $Ø$ indicates an empty set, ${\mathrm{\Lambda }}_t$ indicates the index (column number) found by the $t$th iteration, $a_j$ indicates the $j$th column of the matrix $L$, $L_t$ represents the set of columns of the matrix $L$ selected by the index ${\mathrm{\Lambda }}_t$ (as a matrix of $M\times t$), ${\theta }_t$ is a column vector of $t\times$1, $K$ is the sparsity of the signal. The process is shown in Fig. 7.

The specific flow is as follows:

(1) Initialization: $r_0=v$, ${\mathrm{\Lambda }}_0=Ø$, $L_0=Ø$, $t=$ 1.

(2) Find the index ${\lambda }_t$ so that: ${\lambda }_t=\mathrm{a}\mathrm{r}\mathrm{g}\underset{j=\mathrm{1,2},\cdots ,N}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|⟨r_{t-1},{\alpha }_j⟩\right|$.

(3) Order: ${\mathrm{\Lambda }}_{\mathrm{t}}={\mathrm{\Lambda }}_{\mathrm{t}-1}\cup \left\{{\lambda }_t\right\}$, $L_t=L_{t-1}\cup {\alpha }_{{\lambda }_t}$.

(4) Find the least squares solution of $v=L_t{\theta }_t$: ${\widehat{\theta }}_t=\mathrm{a}\mathrm{r}\mathrm{g}\underset{{\theta }_t}{\mathrm{m}\mathrm{i}\mathrm{n}}‖v-L_t{\theta }_t‖$.

(5) Update residual: $r_t=v-L_t{\widehat{\theta }}_t$.

(6) If $t\le K$, returns to step (2), $t=t+1$, otherwise stop iterating into step (7).

(7) Obtain the estimated value of the original signal $u_t$.

After finding the sound source plane information $u\text{'}$ with sparse information and knowing the transfer matrix $L$, the holographic plane sound pressure distribution map of any near-field plane with the distance $d$ can be obtained from Eq. (13):

Fig. 7The process of OMP reconstruction algorithm

3.1. Discs vibration model

Using Comsol simulation software, two metal discs excited by high-frequency points were set up to simulate the sound field distribution caused by high-frequency vibration. Both discs have a radius of 0.05 m and a thickness of 0.005 m. Both materials are alloy structural steels. The central axis of the two discs is parallel to the $y$-axis, and the two discs are equally distributed along the $y$-axis. The coordinates of the center of the discs distributed on the $x$-$z$ plane are (0.2, 0.2) and (–0.2, –0.2), as shown in Fig. 8(a). A sinusoidal excitation force of $F=\text{5}\mathrm{s}\mathrm{i}\mathrm{n}\left(\text{6000}\pi t\right)+\text{5}$ is applied at the same time in the center of the two metal discs and the unit is Newton. A fixed constraint is applied around the two metal discs to simulate the sound field distribution produced by the high frequency sound source.

A cuboid free sound field space of 0.6 m×1 m×1 m is established, and a perfect matching layer with a thickness of 0.2 m is added around. The presence of a perfect matching layer allows incident sound waves to pass through the interface without reflection, so that no reflected wave interference occurs. The airspace is divided by using a more refined free-splitting tetrahedral mesh, and the two metal discs are divided by a regular free-splitting tetrahedral mesh. Matlab is used to write related programs, and the data obtained by simulation is used to perform operations.

The frequency of the sound source for this simulation is set at 3 kHz. Fig. 8(b) shows the distribution of sound pressure isosurfaces in the space at time $t=$ 0.2 s.

Fig. 8Sound field

a) Sound field three-dimensional model

b) Sound pressure isosurface distribution ($t=$ 0.2 s)

3.2. Sound field holographic simulation imaging analysis

Set two planes, 0.1 m and 0.2 m away from the metal discs. The sound array samples a signal of 0.005 s duration. Calculate the obtained data to obtain the ideal sound pressure distribution of the two planes as shown in Fig. 9(a) and Fig. 9(b). Using an acoustic array formed by 11×11 uniformly arrayed elements, the sound pressure distribution shown in Fig. 9(b) is sampled, as shown in Fig. 10(a) and Fig. 10(b).

Fig. 9Sound pressure distribution maps

a) Distance 0.1 m

b) Distance 0.2 m

Fig. 10Microphones arrangement diagram

a)

b)

Acoustic array sampling is performed in the time domain. Each microphone samples a signal for a certain period of time. Selecting a microphone in the sound array, the sound pressure changes of the signal sampled in the time domain is shown in Fig. 11. It can be seen from the figure that the change in sound pressure at this point tends to be stable after 0.4 s. Therefore, the following are selected after the 0.4 s data calculation.

Fig. 11Time-domain map of a single element signal

The acoustic array is used to sample the sound pressure in the 0.2 m plane from the metal discs in the time domain. The signal durations of all microphones are 0.005 s. The time domain signals obtained by sampling are calculated to obtain the sound pressure values of 121 points.

Set the virtual sound source plane on the actual sound source plane and evenly divide it into 11×11 grids. Uniformly divide the holographic surface into 21×21 grids for sound pressure reconstruction. The acoustic array on the measuring plane 0.2 m from the discs. The sound pressure distribution obtained by sampling is shown in Fig. 12(a).

In the actual ship cabin, most of the noises come from low-frequency noises generated by regular vibrations of other types of equipment or systems, and the noise sources that emit noise are relatively far away from the measurement surface and can be regarded as far-field sound sources. When the noise reaches the measurement plane, it can be approximated as a plane wave, so the noise on the acoustic array plane can be regarded as evenly distributed on each sampling point. Analog noise interference requires the addition of 0.005 s duration noise to the data at each sample point. The sampled values on the measurement surface are calculated, and the average sound pressure amplitude of each sampling point is 0.038 Pa.

Fig. 12Sound pressure sampling data with different noise in the measurement plane

a) Noiseless

b) The noise amplitude is 0.019 Pa

c) The noise amplitude is 0.038 Pa

d) The noise amplitude is 0.19 Pa

The experiment is divided into three groups and each group is added with different noise signals. The specific conditions are:

(1) The adding duration is 0.005 s, and the expression is the noise signal of:

(2) The adding duration is 0.005 s, and the expression is the noise signal of:

(3) The adding duration is 0.005 s, and the expression is the noise signal of:

On the measuring surface of 0.2 m from the discs, three sampled sound pressure distribution diagrams as shown in Fig. 12(b), Fig. 12(c) and Fig. 12(d) are obtained. From the comparison of the above diagrams, it can be seen that there is no significant difference in the visual comparison of the graphs. The main difference is that the range of sound pressure values is different.

When the noise is added to the case (1), the 0.1 m plane sound pressure reconstruction is performed using the NAH based on Fourier transform, the traditional BEM, and the method proposed in this paper. In the following text, the method proposed in this paper is simply referred to as weight method. The results obtained are as follows.

Fig. 13Sound pressure reconstruction maps

a) NAH based on Fourier transform

b) Traditional BEM

c) Weight method

When the noise is added to the case (2), the 0.1 m plane sound pressure reconstruction is performed using the NAH based on Fourier transform, the traditional BEM, and the weight method. The results obtained are as follows.

Fig. 14Sound pressure reconstruction maps

a) NAH based on Fourier transform

b) Traditional BEM

c) Weight method

When the noise is added to the case (3), the 0.1 m plane sound pressure reconstruction is performed using the NAH based on Fourier transform, the traditional BEM, and the weight method. The results obtained are as follows.

Fig. 15Sound pressure reconstruction maps.

a) NAH based on Fourier transform

b) Traditional BEM

c) Weight method

3.3. Results and discussion

In this simulation experiment, three groups of experiments are performed, and each group uses three different sound pressure reconstruction methods. In the three groups of experiments, the average SNR of the sampled signals is 6.02 dB, 0 dB, and –13.98 dB, respectively. When the SNR is zero or negative, it indicates that the signal power is less than the noise power, and the noise is so large that the original signal is submerged.

Under the condition that the noise signal expression is added as:

Compare and analyze Fig. 13. The reconstruction sound pressure distribution obtained by using only the weight method is in good agreement with Fig. 9(a), and the number and location of sound sources can be clearly distinguished. The effect of the reconstruction sound pressure distribution map obtained by the NAH based on Fourier transform and the traditional BEM is not ideal, and only the area where the sound source point is located can be roughly determined, and the number of sound sources cannot be clearly distinguished.

Under the condition that the noise signal expression is added as:

Compare and analyze Fig. 14. The reconstruction sound pressure distribution of each method is basically the same as when the noise signal expression is added as:

Under the condition that the noise signal expression is added as:

Compare and analyze Fig. 15. The reconstruction sound pressure distribution obtained by using only the weight method is in good agreement with Fig. 9(a), and the number and location of sound sources can be clearly distinguished. The reconstructed sound pressure distribution map obtained by adopting the NAH based on Fourier transform and the traditional BEM is greatly affected by the noise at this time, and no longer has a reference value for locating the sound source.

The reason for obtaining the above results is that the NAH based on Fourier transform and the traditional BEM do not perform the high-frequency signal extraction processing on the sampled signals, so that the data used for the subsequent operations contain strong noise interference. The weight method uses the EMD to analyze and process the signals of each sample point. In order to extract the high frequency part of the sampled signal, the weight multiplied by each point is calculated. Therefore, under the same conditions, the weight method is more effective than the NAH based on Fourier transform and the traditional BEM in the localization and sound pressure distribution reconstruction of high-frequency weak sound sources.

Compared with the non-weighted value, the sampled signal multiplied by the weight is closer to the ideal situation without noise. When the noise is added, the weighted value of the sampling signal and the non-weighted sampling signal are respectively compared with the ideal value without noise, and the calculation formula is as follows:

$P_h\left(i\right)$ is the actual sound pressure value measured at each measurement point on the measurement plane 0.2 m away from the plane where the sound source is located. $P_f\left(i\right)$ is the theoretical sound pressure value when no interference noise is added on the corresponding plane. $P_h\left(i\right)$ and $P_f\left(i\right)$ each have 121 values. Eq. (14) is now used to quantify the difference in weighting power. Three sets of simulation experiments are performed. The results of the calculations are shown in the following table:

It can be seen from Table 1. After multiplying the weights, the data at each measurement point is closer to the ideal sound pressure value when there is no interference noise. Therefore, the sound pressure reconstruction map obtained by the weight method under low SNR conditions is closer to the actual situation.

From Fig. 13(c), Fig. 14(c), and Fig. 15(c), it can also be seen that the distribution of the sound source by the weight method is close to the actual situation, but there is a certain gap between the sound pressure distribution value and the actual situation. The reason for this situation is that the modal aliasing phenomenon exists in the EMD itself. This phenomenon can lead to the inability to separate out different modality components according to the feature scale, and the existing IMF contains different time scale components. Therefore, the weights calculated by the EMD are not accurate enough, so that the sound pressure data at each sampling point after multiplying the weights still has a certain gap with the ideal noise-free noise, and the larger the interference noise, the more obvious the gap. Ultimately, the sound pressure values at the points on the resulting sound pressure distribution reconstruction map are not accurate enough.

4.1. Composition of the experimental system

The specific composition of the experimental system is shown in Fig. 16. The power amplifier and the speakers produce a sound field environment with a low signal-to-noise ratio, the sound array performs data acquisition, and the computer calculates and analyzes the data to finally obtain an imaging result.

Fig. 16Experimental system composition

4.2. Experiment procedure

The experimental site is located in an empty, quiet room. The two speakers are placed 3 m in front of the sound array. The spectrum of the noise played by the speakers is shown in Fig. 18(a). The main frequency is around 320 Hz. $r$ is the measurement distance, $d$ is the sound array aperture, $f$ is the sound source frequency, and c is the sound speed. When $f>cr/2d^2$, the sound source can be studied as a near-field sound source. For the two speakers, it is known that $r=$ 3 m, $d=$ 0.635 m, $c=$ 340 m/s, then the sound source frequency needs to satisfy $f$ > 1265 Hz, which can be used as a near-field sound source. Obviously 320 Hz is not in this range. Therefore, the two speakers playing audio can be regarded as far-field noise with respect to the microphone elements on the sound array. The microphone at the center of the array measures the noise level of the noise produced by the two speakers at a certain moment of 62 dB.

The mini speaker is used as the target sound source, and the position is 0.2 m from the sound array. The specific placement is shown in Fig. 17. It plays 4500 Hz noise audio. The spectrum of the audio is shown in Fig. 18(b). At this time, the sound pressure level of the audio played by the mini speaker is 53.2 dB at the center of the array. The signal-to-noise ratio of the signal sampled by the microphone at the center of the array is calculated as a negative value, which can be regarded as a low SNR condition.

Fig. 17Equipments placement

Fig. 18Spectrograms

a) Speaker

b) Mini speaker

c) Mixed noise

Considering that the traditional BEM is not effective in the simulation, it is not verified here. The array acquisition data was calculated by Fourier transform method and weight method respectively, and the sound pressure inversion is performed on a plane with a frequency of 4500 Hz at a distance of 0.1 m from the array. Considering that the mini speaker aperture is small, it is better to select the smaller value of the sparsity $K$ in the weight method, so $K=$ 4 is selected.

Fig. 19Inversion sound pressure distribution map (d= 0.1 m, f= 4500 Hz)

a) Fourier transform method

b) Weight method