阅读说明：本技术 低信噪比及混响环境下的自适应波束成形方法 (Self-adaptive beam forming method under low signal-to-noise ratio and reverberation environment ) 是由 郭文勇 潘兴隆 夏菁 张文群 伍哲 曹承昊 于 2021-08-19 设计创作，主要内容包括：本申请属于复杂声场环境下声源定位方法技术领域,尤其涉及一种低信噪比及混响环境下的自适应波束成形方法。其包括如下步骤：基于声传声器阵列获取声场声压数据的步骤；基于镜像源法计算室内脉冲响应；计算声压互谱矩阵的步骤；计算点传播函数的步骤；原信号重构的步骤；本申请通过构造脉冲响应函数来描述混响场中声波的传播规律,使得算法定位精度相较于常规波束成形法、基于正交匹配追踪的卷积波束成形法等常规算法有较大提高。可以作为低信噪比及混响环境下声源定位方法的有效补充,具有实际工程应用价值。(The application belongs to the technical field of sound source positioning methods in complex sound field environments, and particularly relates to a self-adaptive beam forming method in a low signal-to-noise ratio and reverberation environment. Which comprises the following steps: acquiring sound field sound pressure data based on the sound microphone array; calculating an indoor impulse response based on a mirror image source method; calculating a sound pressure cross-spectrum matrix; calculating a point spread function; reconstructing an original signal; according to the method, the propagation rule of the sound wave in the reverberation field is described by constructing the impulse response function, so that the positioning accuracy of the algorithm is greatly improved compared with the conventional algorithms such as a conventional beam forming method and a convolution beam forming method based on orthogonal matching pursuit. The method can be used as an effective supplement of a sound source positioning method in a low signal-to-noise ratio and reverberation environment, and has practical engineering application value.)

1. A method for adaptive beamforming in a low signal-to-noise ratio and reverberant environment, comprising the steps of:

step A, acquiring sound field sound pressure data based on an acoustic microphone array;

in particular, in the measuring plane S_hObtaining a sound source surface S which contains sound source points and is sparsely distributed by utilizing an M-element microphone array which is regularly arranged_s(ii) a Sound source surface S_sThe method comprises the steps of dividing R areas into equal parts, wherein N areas with sound sources exist; then the sound pressure signal received by the regularly distributed M-ary microphone array is sparsely expressed as: x (t) a · s (t) + n (t);

wherein, at a certain time t, x (t) represents the M × 1 dimension microphone array receiving signal; a is the M multiplied by N dimensional transfer matrix between the sound source and the measuring microphone array; (t) represents an N × 1-dimensional sound source signal; n (t) represents an mx 1-dimensional noise signal;

regarding each focus point of each grid line obtained in the process of equally dividing the sound source surface as a potential sound source, the sound pressure obtained on the measurement surface is the sum of the products of the source intensity of a single sound source and the transfer matrix, and the mathematical expression is recorded as: p is Gq;

in the formula: p represents sound pressure obtained by a measuring surface of the microphone array, and the dimension of M is multiplied by 1; g represents a transfer matrix between a sound source surface and a microphone array measuring surface, and the dimension is M multiplied by N; q represents the source intensity of the sound source, dimension Nx 1; and is

r_mnRepresenting the distance, r, between the nth focus point and the mth array element_nRepresenting the nth focus point and the coordinate sourceThe distance between points i is in imaginary units, r_mnThe distance between a focus point on a sound source surface and an array element on a measuring plane is obtained;

step B, calculating indoor impulse response based on mirror image source method

Specifically, an indoor reverberation mirror image source model is constructed based on a mirror image source method, and signals received by an indoor microphone array are determined

Wherein y (t) is a signal received by the microphone array, h (t) is an indoor impulse response function, x (t) is a sound source signal, n (t) is a noise signal, tau represents a delay amount generated by sound wave reflection, and x represents convolution operation;

performing fourier transform on the above formula to obtain Y (ω) ═ X (ω) H (ω) + N (ω) after transforming into a frequency domain;

thereby constructing an indoor impulse response function

In the formula, beta_x,1、β_y,1And beta_z,1Representing the reflection coefficients of walls close to the origin of coordinates in each direction; beta is a_x,2、β_y,2And beta_z,2Representing the reflection coefficient of the wall away from the origin of coordinates in each direction; r_pRepresenting the distance from the actual sound source or the mirror image sound source to each array element; r_rRepresenting the virtual space size corresponding to the multi-order reflection; u is the speed of sound; q₁、Q₂、Q₃Each 0 or 1 is R_pCan take 8 positions;

converting the indoor impulse response function to the frequency domain has:

for the actual source point coordinates (x, y, z), the sensor array element coordinates (x ', y ', z '), there are:

R_p＝(x-x′+2Q₁x′,y-y′+2Q₂y′,z-z′+2Q₃z′) R_r＝2(n_xL_x,n_yL_y,n_zL_z)

in the formula: n is_x,n_y,n_zIs an integer related to the reflection order;

step C, calculating a sound pressure cross-spectrum matrix;

the sound pressure cross-spectrum matrix C expression is as follows: c ═ pp^H＝Gqq^HG^H(ii) a H refers to conjugation treatment; g is a transfer matrix G; when the sound source is an incoherent sound source, the pair qq^HThe off-diagonal elements in (1) are simplified, qq^HSimplified toThen the expression of the sound pressure cross spectrum matrix C is further expressed as:in the formula: g_nIs the corresponding column vector in the transfer matrix G;

step D, calculating a point spread function;

the conventional beamforming output result based on the cross-spectrum function is: b ═ ω^HCω＝ω^Hpp^Hω；

Where b denotes the acoustic power output from each focusing point, ω denotes a steering vector, and ω ═ ω [ ω ]₁，ω₂…ω_M]，To obtain

Based on the foregoing, a point spread function is obtainedSound source surface acoustic power

Step E, original signal reconstruction step

1) Initializing sparsity K₀1, initializing a sparsity estimation step length L as s, wherein s is a search step length; the supporting set F is equal to phi, and phi is an empty set;

2) calculating the product of the residual error and each column of the sensing matrix to obtain a correlation coefficient matrix u ═ u_j|u_j＝|<r,φ_j>1,2,. N }, where j is the jth column of the sensing matrix and r is the residual;

extracting K from correlation coefficient matrix₀Storing the index value corresponding to the maximum value into F;

3) constrained isometry conditions based on one of the compressive sensing preconditions, ifThen the adaptive dilution K ═ K₀+ L, go to step (2);

wherein y refers to sound field data measured by the acoustic sensor array; delta_KIs a constant between 0 and 1, δ_K0.3; t is a transposed matrix; wherein Φ refers to a sensing matrix in compressed sensing;

if it isGo to step 4;

4) solving initial margin by least square method

5) The initial estimation signal x is equal to 0, the number of initialization stages is equal to 1, the number of initialization iterations k is equal to 1, the set of initialization index values S is equal to phi, and the candidate set C is equal to phi;

6) using the formula u ═ u_j|u_j＝|<r,φ_j>N, and calculating a new correlation coefficient matrix according to u_jThe criterion of > 0.5max | u |, stores the corresponding index valueIn an index set S;

7) merging the index value set T ═ FU S, using the formula u ═ u-_j|u_j＝|<r,φ_j>Calculating a correlation coefficient of an index value corresponding to an atom in the T and the residue, and extracting K₀Storing index value corresponding to maximum value into F_newIn the formula based on the least square methodCalculating an estimated signal x_newAnd use of r_new＝y-Φ_Fx updating the allowance;

8) calculating the error between two iterations and judging whether the error is less than the specified error, if the error between the two iterations is | | x_new-x||₂Stopping iteration if the error value is less than or equal to epsilon, otherwise, turning to the step (9), wherein epsilon refers to a designated error value, and x refers to a previous error value;

9) if r_new||₂≥||r||₂If yes, adding 1 to the number of stages, and jumping to the step (6); if r_new||₂＜||r||₂If F is equal to F_new，r＝r_newAnd k is k +1, and the step (6) is carried out.

Outputting the latest estimation signal x after stopping iteration_new(ii) a And after the algorithm iteration is finished and the result is output, obtaining the finally-desired beam data.

2. The method of claim 1, wherein the acoustic characteristics of each wall surface are expressed by an acoustic absorption coefficient γ and a reflection coefficient γWherein gamma is the sound absorption coefficient of the wall surface.

3. The adaptive beam forming method under the low signal-to-noise ratio and reverberation environment of claim 1, wherein in the step B, in the process of calculating the indoor impulse response based on the mirror image source method, the reflected sound wave is regarded as the direct sound wave from the mirror image sound source to the sensor when the indoor impulse response function is processed, and the sound path is calculated based on the position relationship between the mirror image sound source and the sensor; meanwhile, based on the reflection coefficient of each wall surface, the energy attenuation of the sound wave in the reflection process is considered in the calculation process.

Technical Field

The application belongs to the technical field of sound source positioning methods in complex sound field environments, and particularly relates to a self-adaptive beam forming method in a low signal-to-noise ratio and reverberation environment.

Background

When sound source positioning is performed in a complex sound field environment (such as an actual ship cabin), a large number of noise signals and reflected sound signals in a reverberation environment exist in sound field sound pressure data acquired by an array, capturing and extracting of target sound source signal information are difficult, and the performance of a conventional algorithm is greatly reduced. In this case, to determine the position of the target sound source, it is necessary to consider the interference caused by the ambient noise and the reverberant environment. Beamforming is a noise source identification technique based on microphone array measurements, which is essentially a spatial filtering technique.

In the prior art, the deconvolution beam forming method reversely solves sound source surface sound pressure distribution by introducing a point spread (psf) function, and the spatial resolution is obviously improved compared with the conventional beam forming; the impulse response function is a function for describing the propagation rule of sound waves in a reverberation field, and a common construction method is a mirror image source method; and the compressed sensing can still keep good signal reconstruction accuracy under the environment with low signal-to-noise ratio. Nevertheless, in the conventional compressed sensing technology at present, the CoSaMP method has a high dependency on the sparsity K value, and the SAMP method has a long running time, and has various drawbacks or defects, and when the CoSaMP method is applied to sound source positioning in a low signal-to-noise ratio and reverberation environment, the time consumption is long, the anti-interference capability is low, and the like, so that the further application of the method to sound source positioning in a more complex sound field environment is limited.

Disclosure of Invention

The application aims to provide a self-adaptive beam forming method which is short in time consumption, can effectively remove the influence of RIR on received signals and eliminate reverberation interference in a low signal-to-noise ratio and reverberation environment, and is high in anti-interference capability.

In order to achieve the purpose, the following technical scheme is adopted in the application.

Aiming at the problems of sound source positioning and identification in the low signal-to-noise ratio and reverberation environment, the method for self-adaptive beam forming in the low signal-to-noise ratio and reverberation environment is provided on the basis of an improved beam forming method based on self-adaptive compressed sensing, and comprises the following steps:

step A, acquiring sound field sound pressure data based on an acoustic microphone array;

in particular, in the measuring plane S_hObtaining a sound source surface S which contains sound source points and is sparsely distributed by utilizing an M-element microphone array which is regularly arranged_s(ii) a Sound source surface S_sThe method comprises the steps of dividing R areas into equal parts, wherein N areas with sound sources exist; then the sound pressure signal received by the regularly distributed M-ary microphone array is sparsely expressed as:

X(t)＝A·S(t)+N(t)；

wherein, at a certain time t, x (t) represents the M × 1 dimension microphone array receiving signal; a is the M multiplied by N dimensional transfer matrix between the sound source and the measuring microphone array; (t) represents an N × 1-dimensional sound source signal; n (t) represents an mx 1-dimensional noise signal;

regarding each focus point divided by the sound source surface as a potential sound source, the sound pressure obtained on the measurement surface is the sum of the products of the source intensity of a single sound source and the transfer matrix, and the mathematical expression is recorded as: p is Gq;

in the formula: p represents sound pressure obtained by a measuring surface of the microphone array, and the dimension of M is multiplied by 1; g represents a transfer matrix between a sound source surface and a microphone array measuring surface, and the dimension is M multiplied by N; q represents the source intensity of the sound source, dimension Nx 1; and is

r_mnRepresenting the distance, r, between the nth focus point and the mth array element_nRepresenting the distance between the nth focus point and the coordinate origin;

step B, calculating indoor impulse response based on mirror image source method

Specifically, an indoor reverberation mirror image source model is constructed based on a mirror image source method, and signals received by an indoor microphone array are determined

Wherein y (t) is a signal received by the microphone array, h (t) is an indoor impulse response function, x (t) is a sound source signal, n (t) is a noise signal, t represents a delay amount generated by sound wave reflection, and x represents convolution operation;

performing fourier transform on the above formula to obtain Y (ω) ═ X (ω) H (ω) + N (ω) after transforming into a frequency domain;

thereby constructing an indoor impulse response function

In the formula, beta_x,1、β_y,1And beta_z,1Representing the reflection coefficients of walls close to the origin of coordinates in each direction; beta is a_x,2、β_y,2And beta_z,2Representing the reflection coefficient of the wall away from the origin of coordinates in each direction; r_pRepresenting the distance from the actual sound source or the mirror image sound source to each array element; r_rRepresenting the virtual space size corresponding to the multi-order reflection; u is the speed of sound; q₁、Q₂、Q₃Each 0 or 1 is R_pCan take 8 positions;

converting the indoor impulse response function to the frequency domain has:

for the actual source point coordinates (x, y, z), the sensor array element coordinates (x ', y ', z '), there are:

R_p＝(x-x′+2Q₁x′,y-y′+2Q₂y′,z-z′+2Q₃z′)

R_r＝2(n_xL_x,n_yL_y,n_zL_z)

in the formula: n is_x,n_y,n_zIs an integer related to the reflection order;

step C, calculating a sound pressure cross-spectrum matrix;

the sound pressure cross-spectrum matrix C expression is as follows: c ═ pp^H＝Gqq^HG^H(ii) a H refers to conjugation treatment; g is a transfer matrix G;

when the sound source is an incoherent sound source, the pair qq^HThe off-diagonal elements in (1) are simplified, qq^HSimplified to

Then the expression of the sound pressure cross spectrum matrix C is further expressed as:

in the formula: g_nIs the corresponding column vector in the transfer matrix G;

step D, calculating a point spread function;

the conventional beamforming output result based on the cross-spectrum function is: b ═ ω^HCω＝ω^Hpp^Hω；

Where b denotes the acoustic power output at each grid point, ω denotes a steering vector, and ω ═ ω₁，ω₂…ω_M]，To obtain

Based on the foregoing, a point spread function is obtainedSound source surface acoustic power

Step E, original signal reconstruction step

1) Initializing sparsity K₀1, initializing a sparsity estimation step length L as s, wherein s is a search step length; the supporting set F is equal to phi, and phi is an empty set;

2) calculating the product of the residual error and each column of the sensing matrix to obtain a correlation coefficient matrix u ═ u_j|u_j＝|<r,φ_j>1,2,. N }, where j is the jth column of the sensing matrix and r is the residual;

extracting K from correlation coefficient matrix₀Storing the index value corresponding to the maximum value into F;

3) constrained isometry conditions based on one of the compressive sensing preconditions, ifThen the adaptive dilution K ═ K₀+ L, go to step (2);

if it isGo to step 4;

4) solving initial margin by least square method

5) The initial estimation signal x is equal to 0, the number of initialization stages is equal to 1, the number of initialization iterations k is equal to 1, the set of initialization index values S is equal to phi, and the candidate set C is equal to phi;

6) using the formula u ═ u_j|u_j＝|<r,φ_j>N, and calculating a new correlation coefficient matrix according to u_jStoring the corresponding index value into an index set S according to the standard of more than or equal to 0.5max | u |;

7) merging the index value set T ═ FU S, using the formula u ═ u-_j|u_j＝|<r,φ_j>Calculating a correlation coefficient of an index value corresponding to an atom in the T and the residue, and extracting K₀Storing index value corresponding to maximum value into F_newIn the formula based on the least square methodCalculating an estimated signal x_newAnd use of r_new＝y-Φ_Fx updating the allowance;

8) calculating the error between two iterations and judging whether the error is less than the specified error, if the error between the two iterations is | | x_new-x||₂Stopping iteration if the error value is less than or equal to epsilon, otherwise, turning to the step (9), wherein epsilon refers to a designated error value, and x refers to a previous error value;

9) hollow if |)r_new||₂≥||r||₂If yes, adding 1 to the number of stages, and jumping to the step (6); if r_new||₂＜||r||₂If F is equal to F_new，r＝r_newAnd k is k +1, and the step (6) is carried out.

The further improvement and optimization of the self-adaptive beam forming method under the low signal-to-noise ratio and reverberation environment also comprises that the acoustic characteristics of each wall surface are expressed by a sound absorption coefficient gamma, and a reflection coefficientWherein gamma is the sound absorption coefficient of the wall surface.

The further improvement and optimization of the self-adaptive beam forming method under the low signal-to-noise ratio and reverberation environment further comprises the following steps that in the step B, in the process of calculating the indoor impulse response based on the mirror image source method, the reflected sound wave is regarded as the direct sound wave from the mirror image sound source to the sensor when the indoor impulse response function is processed, and the sound path is calculated based on the position relation between the mirror image sound source and the sensor; meanwhile, based on the reflection coefficient of each wall surface, the energy attenuation of the sound wave in the reflection process is considered in the calculation process.

The beneficial effects are that:

(1) the improved method gets rid of the dependence of the CoSaMP method on the sparsity K value when carrying out compressed sensing reconstruction.

(2) Under the environment with low signal-to-noise ratio and reverberation, the constructed impulse response function is used for replacing the free field transfer function, so that the influence of RIR on the received signal can be effectively removed, and the interference of reverberation is eliminated.

The self-adaptive beam forming method under the low signal-to-noise ratio and reverberation environment has strong anti-interference capability, can be used as effective supplement of a sound source positioning method under the low signal-to-noise ratio and reverberation environment, and has practical engineering application value.

Drawings

FIG. 1 is a model of microphone array signal acquisition;

FIG. 2 is a schematic view of sparse representation of sound sources on a sound source plane;

FIG. 3 is a schematic diagram of a mirror source based calculation;

FIG. 4 is a flow diagram of a method of adaptive beamforming in a low signal-to-noise and reverberant environment;

FIG. 5 is a model of reverberant field acoustic radiation for a ray-acoustic module

FIG. 6 is a schematic of an array position;

FIG. 7 is a ray tracing diagram;

FIG. 8 is a measured surface sound pressure level distribution plot;

FIG. 9 is a measured surface sound pressure level distribution in a noise-free environment;

FIG. 10 is a measured surface sound pressure level distribution in an SNR-40 dB environment;

FIG. 11 is a measured surface sound pressure level distribution in an SNR-20 dB environment;

FIG. 12 is a measured surface sound pressure level distribution in a SNR-10 dB environment;

fig. 13 is a graph of imaging results obtained based on CBF at 1000 Hz;

fig. 14 is a graph of imaging results obtained based on OMP-DAMAS at 1000 Hz;

fig. 15 is a graph of imaging results obtained based on the present application at 1000 Hz;

fig. 16 is a graph of imaging results obtained based on CBF at 2000 Hz;

FIG. 17 is a graph of imaging results obtained based on OMP-DAMAS at 2000 Hz;

fig. 18 is a graph of imaging results obtained based on the present application at 2000 Hz;

fig. 19 is a graph of imaging results obtained based on CBF at f 41000 Hz;

FIG. 20 is a graph of imaging results obtained based on OMP-DAMAS at 4000 Hz;

fig. 21 is a graph of imaging results obtained based on the present application at 4000 Hz;

FIG. 22 is a waveform of the sound pressure of a standard signal source in a background with low signal-to-noise ratio of reverberation;

fig. 23 is a graph of sound pressure spectrum of a standard signal source against a background of low signal-to-noise ratio in reverberation.

Detailed Description

The present application will be described in detail with reference to specific examples.

Based on a compressed sensing theory, signal reconstruction is the key of solution, but in a common reconstruction method, some sparsity K values need to be predicted, and some reconstruction methods have overlong running time, so that the method integrates compressed sampling matching pursuit and sparsity self-adaptive matching, and designs a low signal-to-noise ratio and self-adaptive beamforming method under a reverberation environment, which does not depend on the sparsity K values and has a faster processing speed, based on a backtracking thought and self-adaptive processing. The essence of deconvolution beamforming is the process of solving the system of equations b ═ Ax in reverse, where b is the conventional beamforming output and a is the point spread function psf, in both known cases solving x in reverse. The essence of compressed sensing is that when the observation matrix phi and the measurement vector y are known, x is solved reversely by a corresponding reconstruction method. Both are conceptually known to solve back for the unknown quantities, which has a high degree of similarity, and compressed sensing provides a variety of reconstruction methods that can be solved back. Therefore, b and A can be regarded as input quantities of compressed sensing, and then the original signal x is obtained through reduction by a compressed sensing reconstruction method.

Based on the above basis, the adaptive beamforming method under the low signal-to-noise ratio and reverberation environment of the present application mainly comprises the following steps:

step A, acquiring sound field sound pressure data based on an acoustic microphone array;

specifically, as shown in FIG. 1, the measurement plane S is_hObtaining a sound source surface S which contains sound source points and is sparsely distributed by utilizing an M-element microphone array which is regularly arranged_s(ii) a Sound source surface S_sThe method comprises the steps of dividing R areas into equal parts, wherein N areas with sound sources exist; as shown in fig. 2, the sound pressure signal received by the regularly distributed M-ary microphone array is sparsely represented as:

X(t)＝A·S(t)+N(t)；

wherein, at a certain time t, x (t) represents the M × 1 dimension microphone array receiving signal; a is the M multiplied by N dimensional transfer matrix between the sound source and the measuring microphone array; (t) represents an N × 1-dimensional sound source signal; n (t) represents an mx 1-dimensional noise signal;

regarding each focus point divided by the sound source surface as a potential sound source, the sound pressure obtained on the measurement surface is the sum of the products of the source intensity of a single sound source and the transfer matrix, and the mathematical expression is recorded as: p is Gq;

in the formula: p represents sound pressure obtained by a measuring surface of the microphone array, and the dimension of M is multiplied by 1; g represents a transfer matrix between a sound source surface and a microphone array measuring surface, and the dimension is M multiplied by N; q represents the source intensity of the sound source, dimension Nx 1; and is

r_mnRepresenting the distance, r, between the nth focus point and the mth array element_nRepresenting the distance between the nth focus point and the coordinate origin;

step B, calculating indoor impulse response based on mirror image source method

Specifically, an indoor reverberation mirror image source model is constructed based on a mirror image source method, and signals received by an indoor microphone array are determinedWherein y (t) is a signal received by the microphone array, h (t) is an indoor impulse response function, x (t) is a sound source signal, n (t) is a noise signal, t represents a delay amount generated by sound wave reflection, and x represents convolution operation; performing fourier transform on the above formula to obtain Y (ω) ═ X (ω) H (ω) + N (ω) after transforming into a frequency domain; thereby constructing an indoor impulse response function

In the formula, beta_x,1、β_y,1And beta_z,1Representing the reflection coefficients of walls close to the origin of coordinates in each direction; beta is a_x,2、β_y,2And beta_z,2Representing the reflection coefficient of the wall away from the origin of coordinates in each direction; r_pRepresenting the distance from the actual sound source or the mirror image sound source to each array element; r_rRepresenting the virtual space size corresponding to the multi-order reflection; u is the speed of sound; q₁、Q₂、Q₃Each 0 or 1 is R_pIs convenient to use8 kinds of positions;

as shown in fig. 3, in the process of calculating the indoor impulse response based on the mirror image source method, the reflected sound wave is regarded as the direct sound wave from the mirror image sound source to the sensor when the indoor impulse response function is processed, and the sound path is calculated based on the position relationship between the mirror image sound source and the sensor; meanwhile, based on the reflection coefficient of each wall surface, the energy attenuation of sound waves in the reflection process is considered in the calculation process, the acoustic characteristics of each wall surface are expressed by sound absorption coefficients gamma, and the reflection coefficientsWherein gamma is the sound absorption coefficient of the wall surface.

Converting the indoor impulse response function to the frequency domain has:

for the actual source point coordinates (x, y, z), the sensor array element coordinates (x ', y ', z '), there are:

R_p＝(x-x′+2Q₁x′,y-y′+2Q₂y′,z-z′+2Q₃z′)

R_r＝2(n_xL_x,n_yL_y,n_zL_z)

in the formula: n is_x,n_y,n_zIs an integer related to the reflection order;

step C, calculating a sound pressure cross-spectrum matrix;

the sound pressure cross-spectrum matrix C expression is as follows: c ═ pp^H＝Gqq^HG^H(ii) a H refers to conjugation treatment; g is a transfer matrix G; when the sound source is an incoherent sound source, the pair qq^HThe off-diagonal elements in (1) are simplified, qq^HSimplified toThen the expression of the sound pressure cross spectrum matrix C is further expressed as:in the formula: g_nIs the corresponding column vector in the transfer matrix G;

step D, calculating a point spread function;

the conventional beamforming output result based on the cross-spectrum function is: b ═ ω^HCω＝ω^Hpp^Hω；

Where b denotes the acoustic power output at each grid point, ω denotes a steering vector, and ω ═ ω₁，ω₂…ω_M]，To obtainBased on the foregoing, a point spread function is obtainedSound source surface acoustic power

Step E, original signal reconstruction step

1) Initializing sparsity K₀1, initializing a sparsity estimation step length L as s, wherein s is a search step length; the supporting set F is equal to phi, and phi is an empty set;

2) calculating the product of the residual error and each column of the sensing matrix to obtain a correlation coefficient matrix u ═ u_j|u_j＝|<r,φ_j>1,2,. N }, where j is the jth column of the sensing matrix and r is the residual;

extracting K from correlation coefficient matrix₀Storing the index value corresponding to the maximum value into F;

3) constrained isometry conditions based on one of the compressive sensing preconditions, ifThen the adaptive dilution K ═ K₀+ L, go to step (2);

if it isGo to step 4;

4) solving initial margin by least square method

5) The initial estimation signal x is equal to 0, the number of initialization stages is equal to 1, the number of initialization iterations k is equal to 1, the set of initialization index values S is equal to phi, and the candidate set C is equal to phi;

6) using the formula u ═ u_j|u_j＝|<r,φ_j>N, and calculating a new correlation coefficient matrix according to u_jStoring the corresponding index value into an index set S according to the standard of more than or equal to 0.5max | u |;

7) merging the index value set T ═ FU S, using the formula u ═ u-_j|u_j＝|<r,φ_j>Calculating a correlation coefficient of an index value corresponding to an atom in the T and the residue, and extracting K₀Storing index value corresponding to maximum value into F_newIn the formula based on the least square methodCalculating an estimated signal x_newAnd use of r_new＝y-Φ_Fx updating the allowance;

8) calculating the error between two iterations and judging whether the error is less than the specified error, if the error between the two iterations is | | x_new-x||₂Stopping iteration if the error value is less than or equal to epsilon, otherwise, turning to the step (9), wherein epsilon refers to a designated error value, and x refers to a previous error value;

9) if r_new||₂≥||r||₂If yes, adding 1 to the number of stages, and jumping to the step (6); if r_new||₂＜||r||₂If F is equal to F_new，r＝r_newAnd k is k +1, and the step (6) is carried out.

The schematic flow chart of the improved beam forming method based on the foregoing steps is shown in fig. 4:

based on the steps and the method, verification is performed in a simulation mode;

constructing a reverberation field area with the size of 3m multiplied by 4m multiplied by 3m in a simulation environment, arranging a single-pole source as a sound source at (1.5m,3m,1.5m), and setting the power of the sound source to be 10^-5W, setting the medium in the whole area as air. The fluid model was set to atmospheric attenuation, the relative humidity was set to 50%, and the initial number of rays was set to 100000. The array measuring surface is arranged at a position 1m away from the sound source surface, the array element interval is 0.1m, the array element number is 11 multiplied by 11, and the size of the array measuring surface is 1m multiplied by 1 m.

Fig. 5 and 6 show a reverberation field acoustic radiation model and an array position schematic diagram of the ray acoustic module, wherein when sound waves propagate to each wall surface, reflection and absorption are generated, and reflection and absorption amounts are different due to different sound absorption coefficients of each wall surface. The wall surface condition of the room model in the simulation is referred to a closed room, the closed room can be approximately regarded as a cubic closed space, and the wall surface materials are as follows: the front and left side are bulk glass, the right side and the ceiling are walls, the floor is paved with ceramic tiles, the back is a half-window and half-cabinet, and the values of the sound absorption coefficient of each wall under different frequencies are shown in table 1:

TABLE 1 Sound absorption coefficient of each wall of room

The sound absorption coefficients of the walls at different frequencies were set in the simulated environment according to table 1. In research, the specified time step is calculated, each mode with 0.01s as the step length in 0-0.1 s is calculated, and the parametric scanning is carried out on the sound source frequency. Taking the sound source frequency of 2000Hz and time of 0.1s as an example, the ray trace and the corresponding measured surface sound pressure level distribution are shown in fig. 7 and 8.

In order to obtain the condition of low signal-to-noise ratio in the reverberation field, Gaussian white noise with relative amplitude of 0.01 times, 0.1 times and 0.3 times is added to the obtained measurement surface sound pressure level data respectively, and the Gaussian white noise is used as the measurement surface sound pressure level data obtained under the environment of SNR (signal to noise ratio) of 40dB, 20dB and 10dB respectively. The measurement surface sound pressure level distribution under the low signal-to-noise ratio and reverberation environment is obtained by simulation based on the measurement surface sound pressure level data obtained when the sound source frequency is 2000Hz and the time is 0.1s and adding gaussian white noise with different relative amplitudes as shown in fig. 9, 10, 11 and 12:

the comparative analysis of the above-mentioned FIGS. 9 to 12 can be found: in the reverberation background, because of the reflection and absorption of sound waves, the distribution of the sound pressure level of a measuring surface is more messy than that of a free field, and does not present a regular circle, and more is block-shaped scattered distribution. After the addition of the factors of ambient noise, the phenomenon is consistent with that of free field. Namely: under the condition of no noise, the distribution of the sound pressure level of the measuring surface is relatively concentrated, and the approximate position of the sound source can be judged. After the noise is added, the measurement surface sound pressure level distribution gradually becomes fuzzy or even disordered, and the sound source position can not be judged by measuring the surface sound pressure level distribution condition.

The CBF, OMP-DAMAS, and the method of the present application are now used to process the measured sound field data in an environment with SNR of 0dB, respectively, and the sound pressure distribution of the sound source plane at different frequencies is obtained as shown in fig. 13 to 21:

the specific results of sound source localization by the three methods are shown in table 2:

TABLE 2 sound source localization result comparison table

	1000Hz	2000Hz	4000Hz
				CBF	(1.55m,1.5m)	(1.7m,1.95m)	(1.6m,1.75m)
OMP-DAMAS	(1.5m,1.45m)	(1.8m,1.75m)	(1.15m,1.6m)
				This application	(1.5m,1.5m)	(1.5m,1.5m)	(1.5m,1.45m)

The comparison of the imaging results of the three methods is shown in table 3:

TABLE 3 comparison of imaging results

	Side lobe	Predicting the number of sound sources	Application scope
				CBF	Is provided with	Whether or not	Not applicable to
OMP-DAMAS	Is free of	Is that	Not applicable to
				This application	Is free of	Whether or not	Broad frequency band

As can be seen from the comprehensive analysis of fig. 12, 13, 14 and table 2: under the conditions of the same signal-to-noise ratio and the same sound source frequency, the CBF sound source positioning result has a larger difference with the actual sound source position, and an interference sound source with certain intensity appears in the map. The OMP-DAMAS method has a problem that the positioning accuracy is not good at medium and high frequencies (f 2000Hz, 4000Hz), although no interfering sound source is present in the imaging results. Only the method of the application keeps higher positioning precision and has no side lobe influence, and the performance has obvious advantages compared with the CBF method and the OMP-DAMAS method.

Under the conditions of the same signal-to-noise ratio and different sound source frequencies, the imaging range of the CBF method is gradually reduced along with the increase of the sound source frequency, but the positioning result is not accurate, a certain number of interference sound sources appear, and the judgment of the sound source position cannot be carried out according to the imaging result. As the frequency is increased, although the OMP-DAMAS imaging result has no obvious side lobe, the accuracy is also reduced. The influence of the change of the frequency on the method is small, the imaging result has high precision, no side lobe influence and good stability and reliability.

As can be seen from Table 3, both the CBF method and the OMP-DAMAS method have certain limitations in application, but the method of the present application overcomes the disadvantages of the CBF method and the OMP-DAMAS method and has good practical application value.

In order to verify the sound source positioning and imaging effects of the method in the practical reverberation environment, the noise of a standard signal source is used as a research object, the noise of an air compressor is used as background noise to design and develop an acoustic imaging experiment. The room is a closed three-dimensional space of 5.4m × 3m × 3.7m, the model of the standard signal generator is HD1910, and the models and parameters of the air compressors are shown in Table 4.

TABLE 4 air compressor parameter table

Model number	ZBM-0.1/8 type	Volume flow rate	0.1m3/min
				Operating pressure	0.8Mpa	Power of	1.5KW
Weight of the complete machine	22kg	Overall dimension	630×275×625mm
				Rated speed of rotation	2850r/min	Volume of gas storage tank	30L

Before the experiment begins, the sound pressure levels of the working noise of the standard signal source and the air compressor are measured to determine the signal-to-noise ratio of the whole experimental environment. Firstly, an HD1910 type standard signal generator is utilized to generate a sound wave formed by mixing sine waves of 1200Hz, 2800Hz and 4000Hz, the sound pressure level of the signal generator is adjusted to the maximum, and the sound pressure level at the moment is recorded to be 75.7 dB. And then, closing the signal generator, starting the rear right air compressor, and recording the sound pressure level of 88.8dB when the rear right air compressor normally works. The signal-to-noise ratio of the signal measured by the acoustic array is-13.1 dB when the two work simultaneously. The room is a closed three-dimensional space of 5.4m × 3m × 3.7m, and the specific sound absorption coefficient of each wall surface is shown in table 1. In summary, the entire experimental site can be approximately considered as a low signal-to-noise ratio and reverberant environment.

During experimental measurement, the signal generator and the center of the array are arranged on the same horizontal line, the distance between the signal generator and the center of the array is 1.5m horizontally, and the distance between the signal generator and the array is 0.6m vertically. And simultaneously starting the signal generator and the air compressor, and then starting the acoustic array equipment to start to acquire the sound field data.

And after the test is finished, carrying out Fourier transform on the real-time sound field sound pressure data measured by the sound array to obtain a corresponding spectrogram. The sound pressure waveform and the corresponding spectrogram are shown in fig. 15 and 16.

According to the spectrogram analysis, the method comprises the following steps: at this time, the sound wave in the experimental environment contains 7 main characteristic frequencies, wherein 1200Hz, 2800Hz and 4000Hz are the characteristic frequencies of the sound wave of the signal generator, and the remaining 47.4Hz, 189.8Hz, 569.7Hz and 1325Hz are the characteristic frequencies of the working noise of the air compressor. The noise amplitude of the air compressor is larger under the corresponding characteristic frequency, which shows that the sound pressure level of the working noise of the air compressor is higher than that of the standard signal source noise, and corresponds to the above obtained negative signal-to-noise ratio condition.

The specific results of sound source localization obtained by the three methods, obtained by obtaining corresponding sound source images similarly through the above steps, are shown in table 5:

TABLE 5 sound source localization result comparison table

	1200Hz	2800Hz	4000Hz
				CBF	(0.1m,0.5m)	(0.05m,0.25m)	(0m,0.2m)
OMP-DAMAS	(0.5m,0.5m)	(0.1m,0.2m)	(0.1m,0.2m)
				This application	(0.4m,0m)	(0.05m,0.05m)	(0m,0m)

The comparison of the imaging results of the three methods is shown in table 6:

table 6 comparative table of imaging results

	Side lobe	Predicting the number of sound sources	Application scope
				CBF	Is provided with	Whether or not	Not applicable to
OMP-DAMAS	Is free of	Is that	Not applicable to
				This application	Is free of	Whether or not	Broad frequency band

Analysis can obtain: under the conditions of the same signal-to-noise ratio and the same sound source frequency, the accuracy of the CBF method and the OMP-DAMAS method is poor, and the CBF method is also influenced by wider side lobes. The method has certain error, but has obvious advantages in performance compared with the CBF method and the OMP-DAMAS method within an acceptable range. Under the conditions of the same signal-to-noise ratio and different sound source frequencies, the CBF method and the OMP-DAMAS method locate the sound source position gradually close to the real sound source position along with the increase of the sound source frequency, but the precision is still poor. The method has certain errors, but the errors are gradually reduced and are all within an acceptable range.

Under the condition that the sound source frequency is 1200Hz, the three methods have large errors, and the analysis reason is that the amplitude of 1200Hz sound waves is the minimum of the amplitudes of several characteristic frequency sound waves in the whole section of sound waves shown in FIG. 22, and related information is easily covered by the sound waves of other frequency bands, so that the large errors occur in the imaging result.

As can be seen from Table 6, both the CBF method and the OMP-DAMAS method have certain limitations in application, but the method overcomes the disadvantages of the CBF method and the OMP-DAMAS method and has better practical application value.

The conclusion obtained by the experiment is basically consistent with the conclusion obtained by the simulation experiment, and the method can still keep better effectiveness under the environment with low signal-to-noise ratio and reverberation.

Finally, it should be noted that the above embodiments are only used for illustrating the technical solutions of the present application, and not for limiting the protection scope of the present application, and although the present application is described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that modifications or equivalent substitutions can be made on the technical solutions of the present application without departing from the spirit and scope of the technical solutions of the present application.