Clear imaging method based on far-field beamforming for substation noise source identification

Introduction

In recent years, the problem of substation noise nuisance has gradually gained attention, and the demand for substation noise control has become increasingly strong.^1,2 The premise of noise control is noise source identification, yet few studies have been reported on the substation noise identification. One way to achieve acoustic imaging of substation equipment is sound intensity method.^3,4 It obtains the sound intensity distribution on the measurement surface by scanning point by point, and thus determines the location of the main noise sources. To further determine the noise contribution from each vibration source of the structure, selective sound intensity method combines the vibration acceleration of the structure as the reference signal. Another approach is near-field acoustic holography, ⁵ collecting sound field information by arranging microphone arrays near substation equipment. Near-field acoustic holography is especially accurate at recognizing low-frequency sound sources. However, the above methods can only be applied for small substation equipment and cannot measure multiple devices at once. As a result, none of them are applicable to the identification of noise source in large outdoor substations.

By virtue of fast measurement speed and suitability for medium and long distance measurement, the beamforming noise source identification technology based on microphone arrays has been widely used in the fields of automobiles, airplanes and high-speed railways.^6,7 Microphone array form has also gradually developed from the early linear array to planar array and spherical array. Among them, the solid spherical microphone array cannot only record panoramic sound scenes but also acquire high signal-to-noise ratio, which is especially suitable for noise source identification in large outdoor substations.

Spherical harmonics beamforming (SHB) is a classical post-processing algorithm matched to spherical microphone arrays. ⁸ SHB enables rapid sound source localization, but its imaging clarity is limited. This is mainly due to the wide mainlobes at low frequencies and the serious sidelobe contamination at high frequencies. To improve the clarity of acoustic imaging with spherical array beamforming, filter and sum (FAS), ⁹ functional delay and sum (FDAS) ¹⁰ and deconvolution methods^11,12 have been developed. FAS enhances sidelobe suppression by imposing constraints on the finite impulse response filters set at each focusing point. FDAS combines functional beamforming and the delay-and-sum algorithm to construct the output of spherical array beamforming. It not only has significant sidelobe suppression, but also has high computational efficiency. The deconvolution approaches represent the output of SHB as a convolution of the acoustic source intensity and the point spread function, and gradually removes the influence of the point spread function corresponding to the peaks of mainlobes. It thus simultaneously considers the reduction of mainlobe width and the suppression of sidelobes. However, the above algorithm researches are mainly carried out in the field of near-field noise source identification such as the cabins of automobiles or airplanes,^13,14 which are not suitable to be directly applied to the far-field noise source identification in large outdoor substations.

Therefore, in order to realize the panoramic sound source identification in large substations, this paper firstly derives the output of far-field SHB based on the plane-wave assumption, and then introduces the correction factor to make its output in the direction of the sound source equal to the source intensity. Subsequently, based on the idea of deconvolution methods and combined the corrected far-field SHB output, a clear imaging method for far-field SHB, namely clear SHB (C-SHB), is proposed. C-SHB iteratively searches and locates the sound sources one by one, and finally reconstructs the output of far-field SHB to obtain clear panoramic imaging results.

Theory

Output of far-field SHB

This paper focuses on noise source identification in outdoor trial environments, such as large substations, which belongs to far-field sound source identification. Therefore, the plane wave model is used to construct the far-field SHB output. Assuming that the far-field sound source is located in the direction Ω S and the source intensity is s , the theoretical output of the classical SHB ⁸ can be expressed as

b ( k , Ω F ) = s ∑ n = 0 ∞ ∑ m = − n n Y n m * ( Ω S ) Y n m ( Ω F )

(1)

where b ( k , Ω F ) denotes the output of the classical SHB in the focus direction Ω F , k = 2 π f / c is the wavenumber corresponding to the frequency f and the speed of sound c , Y n m ( Ω ) denotes the spherical harmonic corresponding to the direction Ω , both n and m are the degrees of the spherical harmonic, and “ ( ⋅ ) * ” represents the conjugation. According to the orthogonality of the spherical harmonics, the output of equation (1) equals to infinite when Ω F = Ω S , otherwise it equals to zero, thus realizing the sound source localization.

Due to the sound source direction Ω S is unknown, Y n m ( Ω S ) in equation (1) cannot be solved directly. To address the problem, the spherical Fourier transform coefficients of the pressure signal sampled with spherical microphone array p n m are substituted into equation (1).

b ( k , Ω F ) = ∑ n = 0 ∞ ∑ m = − n n p n m R n ( k a ) Y n m ( Ω F )

(2)

where R n ( k a ) is n th order radial function and a is the radius of the spherical array. For the far-field SHB, R n ( k a ) is based on the plane wave assumption instead of the spherical wave assumption in. ⁸ Therefore, R n ( k a ) takes the form

R n ( k a ) = − 4 π j n ( j n ( k a ) − j n ′ ( k a ) h n ( 2 ) ′ ( k a ) h n ( 2 ) ( k a ) )

(3)

where h n ( 2 ) ( ⋅ ) and h n ( 2 ) ′ ( ⋅ ) are the spherical Hankel function of the second kind and its derivative respectively, j n ( ⋅ ) and j n ′ ( ⋅ ) are the spherical Bessel function of the first kind and its derivative respectively, and j = − 1 is the imaginary unit.

The spherical Fourier transform coefficients p n m are theoretically obtained by integrating the sound pressure over the entire surface of the spherical array. In practical, the number of microphones cannot be infinite, thus it is necessary to approximate the result of the integration with the numerical summation of a finite number of microphone signals.

p n m = S { p ( k a , Ω M q ) } = ∫ Ω q ∈ S 2 = 4 π p ( k a , Ω M q ) Y n m * ( Ω M q ) d Ω M q ≈ ∑ q = 1 Q α q p ( k a , Ω M q ) Y n m * ( Ω M q )

(4)

where p ( k a , Ω M q ) denotes the sound pressure in direction Ω M q on a solid sphere of radius a , “ S ” represents the spherical Fourier transform, “ S 2 ” represents the surface area of a unit sphere, Q is the total number of microphones, and α q is the weight of the q th microphone. In addition, the order of the spherical harmonic n cannot be taken to infinity, and the truncated order N is usually used instead of ∞ . ¹⁵ Finally, the uncorrected output of far-field SHB is obtained as follows.

b ( k , Ω F ) = ∑ n = 0 N ∑ m = − n n ∑ q = 1 Q α q p ( k a , Ω M q ) Y n m * ( Ω M q ) R n ( k a ) Y n m ( Ω F ) = y F N B F N − 1 Y M N H Γ p ⋆

(5)

where y F N is the vector consisting of spherical harmonics corresponding to the direction Ω F , Y M N is the matrix consisting of spherical harmonics corresponding to the directions of all microphones, B F N is the diagonal matrix consisting of radial functions, Γ is the diagonal matrix consisting of the weights of all microphones, p ⋆ is the vector of the microphone array measured sound pressure, and “ ( ⋅ ) − 1 ” and “ ( ⋅ ) H ” represent the inverse of a matrix and the Hermitian transpose respectively. The specific definitions of above variables are as follows.

① Vector/matrix of spherical harmonics

y F N = [ Y 0 0 Ω F ⏟ n = 0 Y 1 − 1 Ω F Y 1 0 Ω F Y 1 1 Ω F ⏟ ⋯ n = 1 Y N − N Ω F ⋯ Y N N Ω F ⏟ n = N ] ∈ ℂ 1 × N + 1 2

(6)

Y M N = [ Y 0 0 ( Ω M 1 ) Y 1 − 1 ( Ω M 1 ) Y 1 0 ( Ω M 1 ) Y 1 1 ( Ω M 1 ) ⋯ Y N − N ( Ω M 1 ) ⋯ Y N N ( Ω M 1 ) Y 0 0 ( Ω M 2 ) Y 1 − 1 ( Ω M 2 ) Y 1 0 ( Ω M 2 ) Y 1 1 ( Ω M 2 ) ⋯ Y N − N ( Ω M 2 ) ⋯ Y N N ( Ω M 2 ) ⋮ ⋮ ⋮ ⋮ ⋯ ⋱ ⋯ ⋮ Y 0 0 ( Ω M Q ) ุ n = 0 Y 1 − 1 ( Ω M Q ) Y 1 0 ( Ω M Q ) Y 1 1 ( Ω M Q ) ุ n = 0 ⋯ Y N − N ( Ω M Q ) ⋯ Y N N ( Ω M Q ) ุ n = N ] ∈ C Q × ( N + 1 ) 2

(7)

where ℂ represents the set of complex numbers.

② Diagonal matrix of radial functions

B F N Diag ( [ R 0 k a ⏟ n = 0 R 1 k a R 1 k a R 1 k a ⏟ n = 1 ⋯ R N k a ⋯ R N k a ⏟ n = N ) ] ∈ ℂ N + 1 2 × N + 1 2

(8)

where “ Diag ( ⋅ ) ” represents the construction of a diagonal matrix with the vector in the bracket as diagonals.

③ Diagonal matrix of microphone weights

Γ = Diag α 1 α 2 ⋯ α Q ∈ ℝ Q × Q

(9)

where ℝ represents the set of real numbers.

④ Vector of sound pressure measured by the array

p ⋆ = p k a , Ω M 1 p k a , Ω M 2 ⋯ p k a , Ω M Q H ∈ ℂ Q × 1

(10)

The SHB output calculated according to equation (5) only enables the estimation of the source direction and cannot accurately quantify the source intensity. For this reason, the correction factor w , which is derived as follows, is therefore introduced. In the case of single sound source, if the focus direction is the same as the source direction ( Ω F = Ω S ), the theoretical output of the classical SHB in equation (1) with the truncated order N can be expressed as

b ( k , Ω S ) = s ( k , Ω S ) ∑ n = 0 N ∑ m = − n n Y n m * ( Ω S ) Y n m ( Ω S ) = ( N + 1 ) 2 4 π s ( k , Ω S )

(11)

In order to make the corrected output of far-field SHB b ( k , Ω S ) = s ( k , Ω S ) , we take the value of the correction factor as w = 4 π / ( N + 1 ) 2 . The corrected output is as follows.

b c ( k , Ω F ) = w ⋅ y F N B F N − 1 Y M N H Γ p ⋆ = 4 π ( N + 1 ) 2 y F N B F N − 1 Y M N H Γ p ⋆

(12)

Both equations (5) and (11) adopt the sound pressure vector measured by the array as input. To facilitate subsequent acoustic imaging clarification, the corrected output of far-field SHB based on the measured cross-spectral matrix is constructed.

b c k , Ω F = E b ^ k , Ω F b ^ * k , Ω F = 4 π N + 1 2 y F N B F N − 1 Y M N H Γ C Γ Y M N B F N − 1 H y F N H

(13)

where “ E ⋅ ” denotes the expectation and C = E p ⋆ p ⋆ H ∈ ℂ Q × Q is the cross-spectral matrix of the array-measured sound pressure.

Clear SHB

To obtain clear imaging of far-field SHB, multiple iterations are required to determine the source locations and source intensities, which in turn reconstruct the output of C-SHB. The iteration process from γ to γ + 1 is as follows.

(1) Search for the maximum output of equation (13) after γ iterations, and consider its location Ω max ( r + 1 ) and the maximum value b c , max ( γ + 1 ) = max ( b c ( γ ) ) as the sound source direction and intensity respectively, where “ max ( ⋅ ) ” denotes taking the maximum value.

(2) Reconstruct the residual sound pressure cross-spectral matrix of the ( γ + 1 ) th iteration C ( r + 1 ) as follows.

{ C γ + 1 = C γ − G γ + 1 G γ + 1 = b max γ + 1 2 t k , Ω max r + 1 H t k , Ω max r + 1 ∈ ℂ Q × Q

(14)

where G ( γ + 1 ) is the cross-spectral matrix corresponding to the source identified in the ( γ + 1 ) th iteration, t ( k , Ω max ( r + 1 ) ) is the focus vector from the source in the direction Ω max ( r + 1 ) to the microphone array, which can be expressed as follows.

{ t k , Ω max r + 1 = t k a , Ω M 1 | Ω max r + 1 t k a , Ω M 2 | Ω max r + 1 ⋯ t k a , Ω M q | Ω max r + 1 ∈ ℂ 1 × Q t k a , Ω M q | Ω max r + 1 = ∑ n = 0 ∞ ∑ m = − n n R n k a Y n m * Ω max r + 1 Y n m Ω M q

(15)

where t ( ( k a , Ω M q ) | Ω max ( r + 1 ) ) is the acoustic transfer function from the source in the direction Ω max ( r + 1 ) to the q th microphone.

(3) Calculate the residual output of far-field SHB after γ + 1 iterations.

{ b c γ + 1 = b c γ + 1 k , Ω F , Ω F ∈ B b c γ + 1 k , Ω F = 4 π N + 1 2 y F N B F N − 1 Y M N H Γ C γ + 1 Γ Y M N B F N − 1 H y F N H

(16)

where b c ( γ + 1 ) ( k , Ω F ) is the residual SHB output in the focus direction Ω F , b c ( γ + 1 ) is the residual SHB output vector consisting of all focus directions after γ + 1 iterations, and B is the set of focus directions.

(4) Reconstruct the C-SHB output of the sound source identified in the ( γ + 1 ) th iteration as follows.

{ b s γ + 1 = b s γ + 1 k , Ω F ⋅ Ψ , Ω F ∈ B b s γ + 1 k , Ω F = 4 π N + 1 2 y F N B F N − 1 Y M N H Γ G γ + 1 Γ Y M N B F N − 1 H y F N H

(17)

where b s ( γ + 1 ) ( k , Ω F ) is the reconstructed C-SHB output in the focus direction Ω F , b s ( γ + 1 ) is the reconstructed C-SHB output vector consisting of all focus directions during the ( γ + 1 )th iteration. Meanwhile, set the beamwidth function Ψ to clarify the acoustic imaging results. When ∠ ( Ω F , Ω max ( r + 1 ) ) ≤ Φ , 0 < Ψ ≤ 1 and Ψ ( 0 ) ≡ 1 ; otherwise, Ψ ≡ 0 . It should be noted that Φ denotes the preset beamwidth and ∠ ( Ω F , Ω max ( r + 1 ) ) represents the angle between the focus direction Ω F and the sound source direction Ω max ( r + 1 ) .

(5) After the termination of iterations, the total beamforming output of C-SHB is calculated as follows.

b s = ∑ γ = 1 γ max b s ( γ )

(18)

Clear SHB can use max ( b s ( γ + 1 ) ) ≥ max ( b s ( γ ) ) as the criterion for iteration termination, while setting the maximum number of iterations γ max .

Simulation and experiment

The simulation is conducted with a 36-channel solid spherical microphone array with the radius of 0.0975 m, as shown in Figure 1, which is also used in the subsequent experiment. The origin of the coordinate system is located at the center of the spherical microphone array. An arbitrary direction in 3D space can be indicated as Ω = ( θ , φ ) , where θ and φ are the elevation and azimuth angles respectively. The symbol “★” indicates the far-field sound source, and the symbols “ • ” indicate microphones arranged on the surface of the solid spherical array.OPEN IN VIEWER

Three sound sources with the directions ( θ , φ ) of ( 80 ∘ , 140 ∘ ) , ( 80 ∘ , 230 ∘ ) and ( 120 ∘ , 170 ∘ ) are simulated. The theoretical sound source intensities of the three sources are 2 Pa, 2 Pa and 2 / 2 Pa respectively, which can be reformulated as 100 dB, 97 dB and 91 dB in decibels. The maximum number of iterations is set to γ max = 50 and the beamwidth is preset to Φ = 20 ∘ .

Figure 2 shows the imaging colormaps of output-uncorrected SHB, output-corrected SHB and C-SHB, with the columns from left to right corresponding to 500 Hz, 1500 Hz, and 4500 Hz. The symbols “+” represent the directions of true sound sources. As can be seen from Figure 2, the mainlobe width of SHB is large at low frequencies, resulting in the weak source being covered, while at high frequencies, the sidelobe amplitudes of SHB are higher than the weak source, resulting in the weak source being masked. Compared with the results of SHB, the imaging clarity of C-SHB is significantly improved by reducing the mainlobe width and suppressing the sidelobe contamination, thus effectively enhancing the spatial resolution and weak source identification capability. In addition, according to the quantization results labeled near the true sound sources, except the output-uncorrected SHB, both the output-corrected SHB and C-SHB can accurately quantify the source intensities, with the C-SHB being more advantageous for quantifying the weak source at low frequency.OPEN IN VIEWER

Figure 3 shows the imaging colormaps of C-SHB at different signal-to-noise ratios (SNRs). The columns from left to right correspond to SNRs of 10 dB, 20 dB, and 30 dB, and the rows from top to bottom correspond to incoherent and coherent sound sources. The frequency is set to 500 Hz, which is the lowest frequency of the simulation in Figure 2. As is shown in Figure 3, C-SHB has good noise resistance and is not affected by source coherence. This is particularly important for outdoor substation noise source identification, because outdoor testing is susceptible to external interference and coherence may exist between sound sources.OPEN IN VIEWER

In order to verify the effectiveness of C-SHB in practical applications, noise source identification based on spherical microphone array was performed in a large substation. The experiment setup is shown in Figure 4. The sound pressure signals are collected by a 36-channel solid spherical array from Brüel & Kjær with the radius of 0.0975 m, whose microphone coordinates are the same as those used in the simulation. The arithmetic average of the A-weighted auto-spectrum of all channels is shown in Figure 5, where frequency band within 1000 Hz is the dominant part. The most significant noise is seen at 500 Hz, 600 Hz and 700 Hz. The acoustic imaging colormaps corresponding to the above three frequencies are further shown in Figure 6. It is obvious that the sound source localization of SHB and C-SHB is basically the same, but the imaging clarity of C-SHB is significantly better than that of SHB, which is consistent with the simulation conclusion.OPEN IN VIEWER

OPEN IN VIEWER

Figure 5.

Auto-spectrum of experimental measurement.

OPEN IN VIEWER

Figure 6.

Experimental imaging result. (a) 500 Hz, SHB, (b) 600 Hz, SHB, (c) 700 Hz, SHB, (d) 500 Hz, C-SHB, (e) 600 Hz, C-SHB, (f) 700 Hz, C-SHB.

Conclusions

In order to realize the identification of noise sources in large substations, we first derive the corrected output of far-field SHB based on the plane-wave assumption, and further propose a clear acoustic imaging method, C-SHB. The simulation results show that the corrected far-field SHB is able to roughly locate the sound source and accurately quantify the source intensity. Subsequently, the newly proposed C-SHB approach provides clearer acoustic imaging with the following advantages: (1) Smaller mainlobes of sound sources and fewer sidelobes, resulting in better spatial resolution of sound source identification and clearer acoustic imaging; (2) Stronger weak source identification capability, including localization and quantization. The experimental result further verified the effectiveness of C-SHB for noise source identification in the large outdoor substation.