阅读说明：本技术 一种低信噪比条件下基于自适应压缩感知的波束成形方法 (Beam forming method based on self-adaptive compressed sensing under low signal-to-noise ratio condition ) 是由 陈汉涛 韩江桂 余丽 施冠羽 黄家宁 余良武 于 2021-08-19 设计创作，主要内容包括：本申请属于波束成形方法技术领域,尤其涉及一种低信噪比条件下基于自适应压缩感知的波束成形方法。包括如下步骤：基于声阵列获取声场声压数据的步骤；计算声压互谱矩阵的步骤；计算点传播函数的步骤；原信号重构的步骤；本申请的低信噪比条件下基于自适应压缩感知的波束成形方法,与常规波束形成法(CBF)、基于正交匹配追踪方法的反卷积波束成形法(OMP-DAMAS)相比,在低信噪比条件下对声源定位效果较好,具备较强的抗干扰能力,可作为低信噪比条件下进行声源定位的一种有效补充方法。(The application belongs to the technical field of beam forming methods, and particularly relates to a beam forming method based on adaptive compressed sensing under the condition of low signal-to-noise ratio. The method comprises the following steps: acquiring sound field sound pressure data based on the sound array; calculating a sound pressure cross-spectrum matrix; calculating a point spread function; reconstructing an original signal; compared with a conventional beam forming method (CBF) and an inverse convolution beam forming method (OMP-DAMAS) based on an orthogonal matching tracking method, the beam forming method based on the adaptive compressed sensing under the condition of low signal-to-noise ratio has a good sound source positioning effect under the condition of low signal-to-noise ratio, has strong anti-jamming capability, and can be used as an effective supplement method for sound source positioning under the condition of low signal-to-noise ratio.)

1. A beamforming method based on adaptive compressed sensing under the condition of low signal-to-noise ratio is characterized by comprising the following steps:

step A, acquiring sound field sound pressure data based on a sound 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 The sound source surface S_sAre divided into R regions including N regions where 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 dimensional array receiving signal; a is an M x N dimensional transfer matrix between the sound source and the 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 the sound pressure obtained by the array measuring surface, and the dimension is M multiplied by 1; g represents a transfer matrix between the sound source surface and the array measurement surface, and the dimension is M multiplied by N; q represents the source intensity of the sound source, dimension Nx 1; and is

g_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 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 When the sound source is an incoherent sound source, the pair qq^HOfSimplification of the diagonal elements, qq^HSimplified toThe 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 C, calculating a point spread function;

traditional beam forming output result b-omega based on cross-spectrum function^HCω＝ω^Hpp^Hω；

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

Step D, 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 phi_jFor the jth column of the sensing matrix, r is the residual, and K is extracted from the correlation coefficient matrix₀Storing the index values corresponding to the maximum values into an index set F;

(3) based on a constrained isometry condition that is one of the compressive sensing preconditions,

if it is notThen 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_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 determining whether the error is less than the fingerDetermining the error if the error between 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; 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.

(9) 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 beamforming method based on adaptive compressed sensing under low snr condition according to claim 1, wherein the beamforming method is performed in a measurement plane S_hThe upper future wave is considered to be a plane wave.

Technical Field

The application belongs to the technical field of beam forming methods, and particularly relates to a beam forming method based on adaptive compressed sensing under the condition of low signal-to-noise ratio.

Background

The beam forming is one of noise source identification technologies based on microphone array measurement, and is essentially a spatial filtering technology, the beam forming idea is to firstly discretize a sound source surface to form a focusing grid, then to use a beam forming method to perform focusing reconstruction on the sound source surface on sound field data acquired by a sound array, and finally to realize sound source positioning in a mode of enhancing output energy of a focus point where the sound source is located. However, in practical application or engineering practice, the target sound source signal is often interfered by the environmental noise signal and even annihilated therein, which makes it difficult to measure and extract the target sound source signal information, thereby affecting the positioning and imaging of the target sound source, resulting in poor effect of the beam forming method and affecting the subsequent processing.

Disclosure of Invention

The application aims to provide an improved beam forming method combined with compressed sensing, based on self-adaption and backtracking thought, the existing scheme is optimized, so that the operation speed of the compressed sensing reconstruction method is higher and does not depend on the sparsity K value, and then the method is combined with a deconvolution beam forming method, so that a better positioning effect and stronger anti-jamming capability are realized.

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

A beamforming method based on adaptive compressed sensing under the condition of low signal-to-noise ratio comprises the following steps:

step A, acquiring sound field sound pressure data based on a sound 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 The sound source surface S_sAre divided into R regions including N regions where 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 dimensional array receiving signal; a is an M x N dimensional transfer matrix between the sound source and the 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 the sound pressure obtained by the array measuring surface, and the dimension is M multiplied by 1; g represents a transfer matrix between the sound source surface and the array measurement surface, and the dimension is M multiplied by N; q represents the source intensity of the sound source, dimension Nx 1; and is

In the formula, g_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 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 When the sound source is an incoherent sound source, the pair qq^HThe off-diagonal elements in (1) are simplified, qq^HSimplified toThe 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 C, calculating a point spread function;

traditional beam forming output result b-omega based on cross-spectrum function^HCω＝ω^Hpp^Hω；

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

Step D, 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 phi_jFor the jth column of the sensing matrix, r is the residual, and K is extracted from the correlation coefficient matrix₀Storing the index values corresponding to the maximum values into an index set F;

(3) based on a constrained isometry condition that is one of the compressive sensing preconditions,

if it is notThen 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_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; 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.

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

For the beam forming method based on the self-adaptive compressed sensing under the condition of the low signal-to-noise ratioFurther development or optimization method, in the measuring plane S_hThe upper future wave is considered to be a plane wave.

The beneficial effects are that:

compared with a conventional beam forming method (CBF) and an inverse convolution beam forming method (OMP-DAMAS) based on an orthogonal matching tracking method, the beam forming method based on the adaptive compressed sensing under the condition of low signal-to-noise ratio has a good sound source positioning effect under the condition of low signal-to-noise ratio, has strong anti-jamming capability, and can be used as an effective supplement method for sound source positioning under the condition of low signal-to-noise ratio.

Drawings

FIG. 1 is a flow chart of a beamforming method based on adaptive compressed sensing under low SNR condition;

FIG. 2 is a schematic diagram of microphone array signal acquisition;

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

FIG. 4 is a time diagram of reconstruction in a noise-free condition according to the prior art method and the method of the present application;

FIG. 5 is a graph illustrating reconstruction error rates for a prior art method and a method of the present application under noiseless conditions;

FIG. 6 is a schematic diagram of reconstruction time of the prior art method and the present application method under the condition of SNR (signal to noise ratio) of 30 dB;

FIG. 7 is a graph illustrating the error rate of prior art methods and the error rate of the present application under the condition of SNR (signal to noise ratio) of 30 dB;

FIG. 8 is a graph of reconstruction times for the prior art method and the present application method under the condition of 0dB SNR;

FIG. 9 is a graph illustrating the reconstruction error rate for the prior art method and the present application method for a signal-to-noise ratio SNR of 0 dB;

Detailed Description

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

The application discloses a beamforming method based on adaptive compressed sensing under a low signal-to-noise ratio condition, which comprises the following steps:

step A, acquiring sound field sound pressure data based on a sound 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 The sound source surface S_sIs divided into R areas including N areas with sound source, and the hollow area indicates that no sound source exists in the area, as shown in FIG. 1 and FIG. 2, and in the implementation process, in the measuring plane S_hThe upper future wave is considered to be a plane wave.

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 dimensional array receiving signal; a is an M x N dimensional transfer matrix between the sound source and the microphone array; (t) represents an N × 1-dimensional sound source signal; n (t) represents an mx 1-dimensional noise signal;

the deconvolution beam forming method takes DAMAS and CLEAN methods as main modes, reversely solves sound source surface sound pressure distribution through introducing a point spread (psf) function, and obviously improves the spatial resolution compared with the conventional beam forming.

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 the sound pressure obtained by the array measuring surface, and the dimension is M multiplied by 1; g represents a transfer matrix between the sound source surface and the array measurement surface, and the dimension is M multiplied by N; q represents the source intensity of the sound source, dimension Nx 1; and is

In the formula, g_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 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 When the sound source is an incoherent sound source, the pair qq^HThe off-diagonal elements in (1) are simplified, qq^HSimplified toThe 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 C, calculating a point spread function;

traditional beam forming output result b-omega based on cross-spectrum function^HCω＝ω^Hpp^Hω；

Where b denotes the acoustic power output at each grid point, ω denotes a steering vector, and ω ═ ω₁，ω₂…ω_M]，To obtainObtaining a point spread functionSound source surface acoustic powerThe above equation gives the relationship between the output of cross-spectral imaging, the point spread function and the source intensity of the sound source, which is essentially a convolution process.

Step D, original signal reconstruction step

The basic compressive sensing reconstruction methods such as Orthogonal Matching Pursuit (OMP), Regularized Orthogonal Matching Pursuit (ROMP), compressive sampling matching pursuit (CoSaMP) and the like are widely applied due to the characteristics of simple structure, low complexity and the like of the methods, but have the defect that the sparsity K value needs to be known in advance, namely, the operation can be carried out by taking the number of sound sources as prior information. In practical applications, the number of sound sources is generally unknown, which limits the application of compressed sensing to some extent. The Sparsity Adaptive Matching Pursuit (SAMP) method proposed for the problem can avoid the need of knowing the sparsity K value, but has the disadvantage of long running time. The method is characterized in that a CoSaMP method and a SAMP method are fused, the backtracking thought of the CoSaMP method and the adaptability of the SAMP method are kept, and an adaptive compression sampling reconstruction method (SAMP-CoSaMP method for short) is obtained, and the method comprises the following specific steps:

(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 phi_jFor the jth column of the sensing matrix, r is the residual, and K is extracted from the correlation coefficient matrix₀Storing the index values corresponding to the maximum values into an index set 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_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; 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.

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

In order to visually reflect the performance of the method, design simulation is carried out, and the running time and the error rate of each method are compared under the conditions of different signal-to-noise ratios. In a sparse signal, the signal observation value M is 128, and the signal length N is 512. Under the conditions of no noise and signal-to-noise ratio of 30dB and 0dB respectively, CoSaMP, SAMP and the method (SAMP-CoSaMP method for short) are adopted for reconstruction, and each group runs 100 times by adopting a Monte Carlo method. Each group compares the running time and the reconstruction error rate of each method respectively from different sparsity K values,

in order to visually reflect the difference between the method and the prior art, the running time and the reconstruction error are firstly compared under the conditions of different signal to noise ratios through simulation experiments.

Taking a certain sparse signal as a test object, the signal observation value M is 128, and the signal length N is 512. Under the conditions of no noise and signal-to-noise ratio of 30dB and 0dB respectively, the existing CoSaMP, SAMP and the method are adopted for reconstruction respectively, and each group of experiments are operated 100 times by adopting a Monte Carlo method. Respectively comparing the running time and the reconstruction error rate of each method under different sparsity K values,

the simulation results are shown in FIGS. 4-9:

as can be seen from fig. 8 and 9 in combination with fig. 7, when the sparsity is the same, the average single run time SAMP algorithm > SAMP-CoSaMP algorithm > CoSaMP algorithm.

K is 45, which is a threshold value of the error rate reconstructed by each algorithm; when K is less than 45, the error rates of the three algorithms are lower and less than 0.01; when K is greater than 45, the error rate of reconstruction of each algorithm rises sharply, and the trend is consistent.

Comparing fig. 7 and fig. 8, it can be seen that, under the noise-free condition, when K >45, the reconstruction errors of the SAMP-CoSaMP algorithm and the SAMP algorithm are consistent with the sparsity variation trend, and are smaller than those of the CoSaMP algorithm, and with the introduction of noise, the reconstruction errors of the algorithms (K >45) under the same sparsity K value are obviously increased.

Therefore, the SAMP-CoSaMP keeps the rapidity and the adaptability of the CoSaMP algorithm, the reconstruction precision of the SAMP-CoSaMP algorithm is between that of the SAMP algorithm and that of the CoSaMP algorithm, and the method can play an important role in the sound source positioning problem of actually unknown sound source number.

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.