CN113687307B - Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment - Google Patents

Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment Download PDF

Info

Publication number
CN113687307B
CN113687307B CN202110957118.3A CN202110957118A CN113687307B CN 113687307 B CN113687307 B CN 113687307B CN 202110957118 A CN202110957118 A CN 202110957118A CN 113687307 B CN113687307 B CN 113687307B
Authority
CN
China
Prior art keywords
sound
sound source
signal
matrix
calculating
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110957118.3A
Other languages
Chinese (zh)
Other versions
CN113687307A (en
Inventor
郭文勇
潘兴隆
夏菁
张文群
伍哲
曹承昊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Naval University of Engineering PLA
Original Assignee
Naval University of Engineering PLA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Naval University of Engineering PLA filed Critical Naval University of Engineering PLA
Priority to CN202110957118.3A priority Critical patent/CN113687307B/en
Publication of CN113687307A publication Critical patent/CN113687307A/en
Application granted granted Critical
Publication of CN113687307B publication Critical patent/CN113687307B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/22—Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

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 low signal-to-noise ratio and reverberation environments. The method 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 application, the propagation rule of sound waves in the reverberant field is described by constructing the impulse response function, so that the algorithm positioning accuracy is greatly improved compared with the conventional algorithms such as the conventional beamforming method, the convolution beamforming method based on orthogonal matching tracking and the like. Can be used as effective supplement of sound source localization method under low signal-to-noise ratio and reverberation environment, and has practical engineering application value.

Description

Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment
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 low signal-to-noise ratio and reverberation environments.
Background
When sound source localization is carried out in a complex sound field environment (such as an actual ship cabin), a large amount of noise signals exist in sound field sound pressure data acquired by an array, reflected sound signals exist in a reverberation environment, the 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 ambient noise and the disturbance caused by 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 solves the sound pressure distribution of the sound source surface in an inverse way by introducing a point propagation (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 the reverberant 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 of low signal-to-noise ratio. Nevertheless, in the conventional compressed sensing technology at present, the CoSaMP method has high dependence on the sparseness K value, and the SAMP method has long running time, has various drawbacks or defects, and is long in time consumption, low in anti-interference capability and the like when being applied to sound source positioning under low signal-to-noise ratio and reverberation environments, so that further application of the method for sound source positioning under more complex sound field environments is limited.
Disclosure of Invention
The application aims to provide the self-adaptive beam forming method with shorter time consumption, which can effectively remove the influence of RIR on a received signal to eliminate reverberation interference under the conditions of low signal to noise ratio and reverberation environment and has strong anti-interference capability under the conditions of low signal to noise ratio and reverberation environment.
In order to achieve the above purpose, the present application adopts the following technical scheme.
Aiming at the problems of sound source positioning and identification under the low signal-to-noise ratio and reverberation environment, based on an improved beamforming method based on adaptive compressed sensing, the adaptive beamforming method under the low signal-to-noise ratio and reverberation environment is provided, which comprises the following steps:
step A, acquiring sound field sound pressure data based on an acoustic microphone array;
specifically, in the measurement plane S h Acquiring a sparse sound source surface S containing sound source points by using an M-element microphone array which is arranged according to rules s The method comprises the steps of carrying out a first treatment on the surface of the Sound source surface S s The method comprises the steps of equally dividing R areas, wherein the R areas comprise N areas with sound sources; 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 receiving signal of the M multiplied by 1 dimension microphone array; a is an M x N dimension transfer matrix between a sound source and a measurement microphone array; s (t) represents an N1-dimensional sound source signal; n (t) represents an mx 1-dimensional noise signal;
and taking each focusing point divided by the sound source surface as a potential sound source, wherein the sound pressure obtained on the measurement surface is the sum of products of the source intensity of a single sound source and a transmission matrix, and the mathematical expression is expressed as: p=gq;
wherein: p represents the sound pressure obtained by the measuring surface of the microphone array, and M is multiplied by 1; g represents a transfer matrix between a sound source surface and a microphone array measurement surface, and M is multiplied by N; q represents the source intensity of the sound source, and N is multiplied by 1; and is also provided with
r mn Represents the distance between the nth focus point and the mth array element, r n Representing the distance between the nth focal point and the origin of coordinates;
step B, calculating indoor impulse response based on mirror image source method
Specifically, an indoor reverberation mirror image source model is built based on a mirror image source method, and an indoor microphone array receiving signal is 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;
fourier transforming the above formula to a frequency domain to obtain Y (ω) =x (ω) H (ω) +n (ω);
thereby constructing an indoor impulse response function
Wherein beta is x,1 、β y,1 And beta z,1 Representing the reflection coefficient of the wall near the origin of coordinates in each direction; beta x,2 、β y,2 And beta z,2 Representing the reflection coefficient of the wall away from the origin of coordinates in each direction; r is R p Representing the distance of the actual or mirror image sound source to each array element；R r Representing the virtual space size corresponding to the multi-order reflection; u is the sound velocity; q (Q) 1 、Q 2 、Q 3 R is 0 or 1 p 8 positions can be taken;
the conversion of the room impulse response function to the frequency domain is:
for the actual sound source point coordinates (x, y, z), the sensor element coordinates (x ', y ', z ') are:
R p ＝(x-x′+2Q 1 x′,y-y′+2Q 2 y′,z-z′+2Q 3 z′)
R r ＝2(n x L x ,n y L y ,n z L z )
wherein: n is n x ,n y ,n z Is 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: c=pp H ＝Gqq H G H The method comprises the steps of carrying out a first treatment on the surface of the H refers to conjugation treatment; g is a transfer matrix G;
when the sound source is incoherent, then the frequency of qq H Simplifying off-diagonal elements in qq H Simplified into
The sound pressure cross-spectrum matrix C is expressed as further expressed as:
wherein: g n For the corresponding column vector in the transfer matrix G;
step D, calculating a point propagation function;
the traditional beam forming output result based on the cross spectrum function is as follows: b=ω H Cω＝ω H pp H ω；
Where b represents the acoustic power output by each grid point, ω represents the steering vector, ω= [ ω ] 1 ，ω 2 …ω M ]，Obtain->
Based on the foregoing, a point spread function is obtainedSound source surface acoustic power->
Step E, a step of original signal reconstruction
1) Initializing sparsity K 0 Initializing a sparsity estimation step l=s, s being a search step; the support set f=Φ, Φ 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 >I, j=1, 2,..n }, where j is the j-th column of the sensing matrix and r is the residual;
extracting K from the correlation coefficient matrix 0 The index value corresponding to the maximum value is stored in F;
3) Constraint equidistant conditions based on one of the compressed sensing preconditions, ifAdaptive dilution k=k 0 +L, turning to step (2);
if it isThen go to step 4;
4) Solving for initial margin by least square method
5) The initial estimation signal x=0, the number of initialization stages stage=1, the number of initialization iterations k=1, the initialization index value set s=phi, and the candidate set c=phi;
6) Using u= { u j |u j ＝|<r,φ j >New correlation coefficient matrices are calculated for j=1, 2,..n } and are calculated as per u j The corresponding index value is stored into the index set S according to the standard of more than or equal to 0.5 max|u|;
7) The index value set t=f u S is merged by the formula u= { u j |u j ＝|<r,φ j >Calculating the correlation coefficient of the index value corresponding to the atom and the allowance in T, and extracting K 0 The index value corresponding to the maximum value is stored in F new In the method, based on the least square method, the method is adoptedCalculating an estimated signal x new And utilize r new ＝y-Φ F x updating the allowance;
8) Calculating the error between the two iterations and judging whether the error is smaller than the specified error, if so, the error between the two iterations is x new -x|| 2 Stopping iteration if epsilon is less than or equal to epsilon, otherwise turning to the step (9), wherein epsilon refers to a specified error, and x refers to a previous error value;
9) If r new || 2 ≥||r|| 2 Adding 1 to the number of stages, stage=stage+1, l=stage×s, and skipping step (6); if r new || 2 ＜||r|| 2 F=f new ，r＝r new K=k+1, go to step (6).
Further improvements and optimizations to the foregoing adaptive beamforming method in low snr and reverberant environments include that each wall acoustic characteristic is represented by a sound absorption coefficient γ, a reflection coefficientWherein gamma is the sound absorption coefficient of the wall surface.
The method for self-adaptive beam forming under the low signal-to-noise ratio and reverberation environment further comprises the following steps that in the process of calculating the indoor impulse response based on a mirror image source method, reflected sound waves are regarded as direct sound waves from a mirror image sound source to a sensor when the indoor impulse response function is processed, and sound paths are calculated based on the position relation between the mirror image sound source and the sensor; and meanwhile, the energy attenuation of the sound wave in the reflection process is considered in the calculation process based on the reflection coefficient of each wall surface.
The beneficial effects are that:
(1) The improvement method gets rid of the dependence of the CoSaMP method on the sparsity K value when the compressed sensing reconstruction is carried out.
(2) Under the low signal-to-noise ratio and reverberation environment, the constructed impulse response function is used for replacing the free field transfer function of the signal, so that the influence of RIR on a 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 stronger 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 microphone array signal acquisition model;
FIG. 2 is a schematic representation of a sparse representation of a sound source on a sound source surface;
FIG. 3 is a schematic diagram of a calculation based on a mirror source method;
FIG. 4 is a flow chart of a method of adaptive beamforming in low signal-to-noise and reverberant environments;
FIG. 5 is a model of reverberant field sound radiation for a ray acoustic module
FIG. 6 is a schematic diagram of an array position;
FIG. 7 is a ray trace diagram;
FIG. 8 is a measurement face acoustic pressure level distribution plot;
FIG. 9 is a measurement surface sound pressure level distribution in a noise-free environment;
fig. 10 is a measurement plane sound pressure level distribution in an snr=40 dB environment;
fig. 11 is a measurement plane sound pressure level distribution in an snr=20 dB environment;
fig. 12 is a measurement surface sound pressure level distribution in an snr=10 dB environment;
fig. 13 is a graph of imaging results obtained based on CBF at f=1000 Hz;
fig. 14 is a graph of imaging results obtained based on OMP-DAMAS at f=1000 Hz;
fig. 15 is a graph of imaging results obtained based on the present application at f=1000 Hz;
fig. 16 is a graph of imaging results obtained based on CBF at f=2000 Hz;
fig. 17 is a graph of imaging results obtained based on OMP-DAMAS at f=2000 Hz;
fig. 18 is a graph of imaging results obtained based on the present application at f=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 f=4000 Hz;
fig. 21 is a graph of imaging results obtained based on the present application at f=4000 Hz;
FIG. 22 is a waveform of standard signal source sound pressure in a reverberant low signal to noise ratio background;
fig. 23 is a graph of standard signal source sound pressure spectrum in a reverberant low signal to noise ratio background.
Detailed Description
The present application will be described in detail with reference to specific examples.
Based on compressed sensing theory, the reconstruction of signals is the key of solving, but in the common reconstruction method, some sparsity K values need to be predicted, and some running time is too long, so the application designs the self-adaptive beamforming method under the low signal-to-noise ratio and reverberation environment which does not depend on the sparsity K values and has higher processing speed by fusing compressed sampling matching tracking and sparsity self-adaptive matching and based on backtracking thinking 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, a is the point propagation function psf, and x is solved in reverse in the case of both being known. The essence of compressed sensing is that when the observation matrix phi and the measurement vector y are known, the inverse solution is carried out on x through a corresponding reconstruction method. Both are thought to be known two quantities to solve for the unknown quantity in the reverse direction, which has high similarity, and compressed sensing provides multiple reconstruction methods to solve in the reverse direction. Therefore, b and A can be regarded as compressed sensing input quantity, and then the original signal x is obtained through reduction by a compressed sensing reconstruction method.
Based on the above, the self-adaptive beam forming method under the low signal-to-noise ratio and reverberation environment 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, in the measurement plane S h Acquiring a sparse sound source surface S containing sound source points by using an M-element microphone array which is arranged according to rules s The method comprises the steps of carrying out a first treatment on the surface of the Sound source surface S s The method comprises the steps of equally dividing R areas, wherein the R areas comprise N areas with sound sources; as shown in fig. 2, the sound pressure signals received by the regularly distributed M-ary microphone array are sparsely represented as:
X(t)＝A·S(t)+N(t)；
wherein, at a certain time t, X (t) represents the receiving signal of the M multiplied by 1 dimension microphone array; a is an M x N dimension transfer matrix between a sound source and a measurement microphone array; s (t) represents an N1-dimensional sound source signal; n (t) represents an mx 1-dimensional noise signal;
and taking each focusing point divided by the sound source surface as a potential sound source, wherein the sound pressure obtained on the measurement surface is the sum of products of the source intensity of a single sound source and a transmission matrix, and the mathematical expression is expressed as: p=gq;
wherein: p represents the sound pressure obtained by the measuring surface of the microphone array, and M is multiplied by 1; g represents a transfer matrix between a sound source surface and a microphone array measurement surface, and M is multiplied by N; q represents the source intensity of the sound source, and N is multiplied by 1; and is also provided with
r mn Represents the distance between the nth focus point and the mth array element, r n Representing the distance between the nth focal point and the origin of coordinates;
step B, calculating indoor impulse response based on mirror image source method
Specifically, an indoor reverberation mirror image source model is built based on a mirror image source method, and an indoor microphone array receiving signal is 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; fourier transforming the above formula to a frequency domain to obtain Y (ω) =x (ω) H (ω) +n (ω); thereby constructing an indoor impulse response function
Wherein beta is x,1 、β y,1 And beta z,1 Representing the reflection coefficient of the wall near the origin of coordinates in each direction; beta x,2 、β y,2 And beta z,2 Representing the reflection coefficient of the wall away from the origin of coordinates in each direction; r is R p Representing the distance from the actual sound source or the mirror image sound source to each array element; r is R r Representing the virtual space size corresponding to the multi-order reflection; u is the sound velocity; q (Q) 1 、Q 2 、Q 3 R is 0 or 1 p 8 positions can be taken;
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 during the processing of the indoor impulse response function, 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, and the acoustic characteristics of each wall surface are represented by sound absorption coefficient gamma, and the reflection coefficientWherein gamma is the sound absorption coefficient of the wall surface.
The conversion of the room impulse response function to the frequency domain is:
for the actual sound source point coordinates (x, y, z), the sensor element coordinates (x ', y ', z ') are:
R p ＝(x-x′+2Q 1 x′,y-y′+2Q 2 y′,z-z′+2Q 3 z′)
R r ＝2(n x L x ,n y L y ,n z L z )
wherein: n is n x ,n y ,n z Is 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: c=pp H ＝Gqq H G H The method comprises the steps of carrying out a first treatment on the surface of the H refers to conjugation treatment; g is a transfer matrix G; when the sound source is incoherent, then the frequency of qq H Simplifying off-diagonal elements in qq H Simplified intoThe sound pressure cross-spectrum matrix C is expressed as further expressed as: />Wherein: g n For the corresponding column vector in the transfer matrix G;
step D, calculating a point propagation function;
the traditional beam forming output result based on the cross spectrum function is as follows: b=ω H Cω＝ω H pp H ω；
Where b represents the acoustic power output by each grid point, ω represents the steering vector, ω= [ ω ] 1 ，ω 2 …ω M ]，Obtain->Based on the foregoing, a point spread function is obtainedSound source surface acoustic power->
Step E, a step of original signal reconstruction
1) Initializing sparsity K 0 Initializing a sparsity estimation step l=s, s being a search step; the support set f=Φ, Φ 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 >I, j=1, 2,..n }, where j is the j-th column of the sensing matrix and r is the residual;
extracting K from the correlation coefficient matrix 0 The index value corresponding to the maximum value is stored in F;
3) Constraint equidistant conditions based on one of the compressed sensing preconditions, ifAdaptive dilution k=k 0 +L, turning to step (2);
if it isThen go to step 4;
4) Solving for initial margin by least square method
5) The initial estimation signal x=0, the number of initialization stages stage=1, the number of initialization iterations k=1, the initialization index value set s=phi, and the candidate set c=phi;
6) Using u= { u j |u j ＝|<r,φ j >New correlation coefficient matrices are calculated for j=1, 2,..n } and are calculated as per u j The corresponding index value is stored into the index set S according to the standard of more than or equal to 0.5 max|u|;
7) The index value set t=f u S is merged by the formula u= { u j |u j ＝|<r,φ j >Calculating the correlation coefficient of the index value corresponding to the atom and the allowance in T, and extracting K 0 The index value corresponding to the maximum value is stored in F new In the method, based on the least square method, the method is adoptedCalculating an estimated signal x new And utilize r new ＝y-Φ F x updating the allowance;
8) Calculating the error between the two iterations and judging whether the error is smaller than the specified error, if so, the error between the two iterations is x new -x|| 2 Stopping iteration if epsilon is less than or equal to epsilon, otherwise turning to the step (9), wherein epsilon refers to a specified error, and x refers to a previous error value;
9) If r new || 2 ≥||r|| 2 Adding 1 to the number of stages, stage=stage+1, l=stage×s, and skipping step (6); if r new || 2 ＜||r|| 2 F=f new ，r＝r new K=k+1, go to step (6).
The schematic flow diagram of the beam forming method based on the above-mentioned step improvement is shown in fig. 4:
based on the steps and the method, verification and explanation are carried out in a simulation mode;
constructing a reverberant field area with the size of 3m multiplied by 4m multiplied by 3m in a simulation environment, setting a monopole source at (1.5 m,3m and 1.5 m) as a sound source, and setting the power of the sound source to be 10 -5 W, the medium of the whole area is set as air. The fluid model was set to atmospheric attenuation, the relative humidity was set to 50%, and the initial ray count was set to 100000. The array measuring surface is arranged at a position 1m away from the sound source surface, the array element distance is 0.1m, the array element number is 11 multiplied by 11, and the array measuring surface is 1m multiplied by 1m.
The reverberant field sound radiation model and the array position schematic diagrams of the ray acoustic module are shown in fig. 5 and 6, reflection and absorption are generated when sound waves propagate to each wall surface, 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 takes a closed room as a reference, and the closed room can be approximately regarded as a cube closed space, and the wall surface materials are as follows: the front and left sides are large glass, the right side and the ceiling are walls, the ground is paved with ceramic tiles, the back is a half window half cabinet, and the values of the sound absorption coefficients of the wall surfaces under different frequencies are obtained as shown in table 1:
table 1 table of sound absorption coefficients of walls of rooms
The sound absorption coefficients of the walls at different frequencies were set in a simulated environment according to table 1. When the research is carried out, the designated time step is calculated, each mode taking 0.01s as the step length is calculated within 0-0.1 s, and the parameterized scanning is carried out on the sound source frequency. Taking the case of a sound source frequency of 2000Hz and a time of 0.1s as an example, the ray trace and the corresponding measurement surface sound pressure level distribution are shown in fig. 7 and 8.
To obtain low signal-to-noise ratio conditions in the reverberant field, gaussian white noise with relative amplitudes of 0.01 times, 0.1 times and 0.3 times is added to the obtained measurement surface sound pressure level data, and the obtained measurement surface sound pressure level data are used as the measurement surface sound pressure level data in the environments of snr=40 dB, 20dB and 10dB respectively. The measurement surface sound pressure level distribution under the conditions of low signal-to-noise ratio and reverberation environment is obtained by simulation by adding Gaussian white noise with different relative amplitudes based on the measurement surface sound pressure level data obtained when the sound source frequency is 2000Hz and the time is 0.1s as shown in figures 9, 10, 11 and 12:
the comparative analysis of the above-described fig. 9 to 12 can be performed: in the reverberation background, because of the reflection and absorption of sound waves, the sound pressure level distribution of a measuring surface is more messy than that of a free field, and the sound pressure level distribution is not regular round any more, and is more scattered in a block shape. After adding the factor of the ambient noise, the phenomenon remains the same as when the free field is present. The method comprises the following steps: under the noiseless condition, the sound pressure level distribution of the measuring surface is relatively concentrated, and the approximate position of the sound source can be judged. After adding noise, the measurement surface sound pressure level distribution becomes gradually blurred or even confused, and the sound source position cannot be judged through the measurement surface sound pressure level distribution condition.
The measured sound field data are processed in the snr=0 dB environment by using CBF, OMP-DAMAS and the method of the present application, respectively, to obtain sound pressure distribution of the sound source surface at different frequencies as shown in fig. 13 to 21:
the specific results of sound source localization performed by the three methods are shown in table 2:
table 2 comparison of sound source localization results
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)
The application is that (1.5m,1.5m) (1.5m,1.5m) (1.5m,1.45m)
The imaging results of the three methods are compared in Table 3:
table 3 comparison of imaging results
Sidelobe Predicting the number of sound sources Application range
CBF Has the following components Whether or not Is not suitable for
OMP-DAMAS Without any means for Is that Is not suitable for
The application is that Without any means for Whether or not Broadband light source
As can be seen from a comprehensive analysis of fig. 12, 13 and 14 and table 2: under the conditions of the same signal-to-noise ratio and the same sound source frequency, the sound source positioning result of the CBF method has a larger gap from the actual sound source position, and an interference sound source with certain intensity appears in the figure. Although no disturbing sound source appears in the OMP-DAMAS imaging result, the problem of poor positioning accuracy exists at the middle and high frequencies (f=2000 Hz and 4000 Hz). Only the method of the application keeps higher positioning precision and has no side lobe influence, and the performance is obviously superior to that of 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 gradually reduces along with the increase of the sound source frequency, but the positioning result is inaccurate and a certain number of interference sound sources can appear, and the sound source position cannot be judged according to the imaging result. With the increase of frequency, the OMP-DAMAS method imaging result has no obvious side lobe, but has the problem of reduced precision. The frequency change has little influence on the method of the application, the imaging result has higher precision, no side lobe influence and good stability and reliability.
As shown in Table 3, the CBF method and the OMP-DAMAS method have certain limitations in application, and the method of the application overcomes the defects of the CBF method and the OMP-DAMAS method and has good practical application value.
In order to verify the sound source localization and imaging effect of the method in the actual reverberation environment, the method is designed and an acoustic imaging experiment is carried out by taking standard signal source noise as a research object and air compressor noise as background noise. The room is 5.4mX3mX3.7m enclosed space, the model of standard signal generator is HD1910, and the model and parameters of air compressor are shown in Table 4.
Table 4 air compressor parameter table
Model number ZBM-0.1/8 type Volumetric flow rate 0.1m 3 /min
Working pressure 0.8Mpa Power of 1.5KW
Weight of whole machine 22kg External dimension 630×275×625mm
Rated rotational speed 2850r/min Air storage tank volume 30L
Before the experiment starts, the standard signal source and the working noise sound pressure level of the air compressor are measured to determine the signal to noise ratio of the whole experimental environment. Firstly, a standard HD1910 type signal generator is used for generating a section of sound wave formed by mixing sine waves of 1200Hz, 2800Hz and 4000Hz, the sound pressure level of the signal generator is adjusted to be maximum, and the sound pressure level is recorded to be 75.7dB at the moment. And then the signal generator is closed, the right rear air compressor is started, and the sound pressure level of the air compressor in normal operation is recorded to be 88.8dB. Calculations show that when the two are operated simultaneously, the signal to noise ratio of the signal measured by the acoustic array is-13.1 dB. The room was a 5.4mX3mX3.7mclosed space, and the specific sound absorption coefficients of the wall surfaces are shown in Table 1. In summary, the whole experimental field can be regarded approximately 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 horizontal distance between the signal generator and the center of the array is 1.5m, and the vertical distance between the signal generator and the center of the array is 0.6m. Simultaneously starting the signal generator and the air compressor, and then starting the sound array equipment to collect sound field data.
After the test is completed, 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.
From the spectrogram analysis, it can be seen that: at this time, the main characteristic frequencies of the sound waves in the experimental environment are 7, wherein 1200Hz, 2800Hz and 4000Hz are characteristic frequencies of the sound waves of the signal generator, and the remaining 47.4Hz, 189.8Hz, 569.7Hz and 1325Hz are characteristic frequencies of the working noise of the air compressor. The noise amplitude is larger under the corresponding characteristic frequency of the air compressor, which indicates that the working noise sound pressure level of the air compressor is higher than the standard signal source noise, and the negative signal-to-noise ratio condition obtained above is corresponding.
Corresponding sound source imaging is obtained through the similar steps, and specific results of sound source positioning performed by the three methods are shown in Table 5:
TABLE 5 Sound source localization results 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)
The application is that (0.4m,0m) (0.05m,0.05m) (0m,0m)
The comparison of imaging results of the three methods is shown in Table 6:
table 6 comparison of imaging results
Sidelobe Predicting the number of sound sources Application range
CBF Has the following components Whether or not Is not suitable for
OMP-DAMAS Without any means for Is that Is not suitable for
The application is that Without any means for Whether or not Broadband light source
The analysis can be obtained: under the conditions of the same signal-to-noise ratio and the same sound source frequency, the CBF method and the OMP-DAMAS method are poor in precision, and the CBF method is also affected 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 position the sound source position gradually close to the real sound source position along with the rising of the sound source frequency, but the precision is still poor. The method has certain error, but the error is gradually reduced and is within an acceptable range.
Under the condition that the sound source frequency is 1200Hz, the three methods all have larger errors, and the analysis is that in the whole section of sound wave shown in fig. 22, the amplitude of the sound wave of 1200Hz is the smallest in the amplitudes of the sound waves of several characteristic frequencies, and related information is easily covered by the sound waves of other frequency bands, so that larger errors appear in the imaging result.
As shown in Table 6, both the CBF method and the OMP-DAMAS method have certain limitations in application, and the method of the application overcomes the defects of the CBF method and the OMP-DAMAS method and has good practical application value.
The results obtained by the experiment are basically consistent with the simulation experiment results, and can prove that the method can still keep better effectiveness under the low signal-to-noise ratio and reverberation environment.
Finally, it should be noted that the above embodiments are only for illustrating the technical solution of the present application, and not for limiting the scope of the present application, and although the present application has been 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 to the technical solution of the present application without departing from the spirit and scope of the technical solution of the present application.

Claims (3)

1. An adaptive beamforming method 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;
specifically, in the measurement plane S h Obtaining the content by using M element microphone arrays arranged according to rulesSparse sound source surface S with sound source points s The method comprises the steps of carrying out a first treatment on the surface of the Sound source surface S s The method comprises the steps of equally dividing R areas, wherein the R areas comprise N areas with sound sources; 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 receiving signal of the M multiplied by 1 dimension microphone array; a is an M x N dimension transfer matrix between a sound source and a measurement microphone array; s (t) represents an N1-dimensional sound source signal; n (t) represents an mx 1-dimensional noise signal;
taking each focusing point of each grid line obtained in the process of equally dividing the sound source surface as a potential sound source, and then obtaining sound pressure on the measuring surface, namely the sum of products of single sound source intensity and a transfer matrix, wherein mathematical expression is recorded as: p=gq;
wherein: p represents the sound pressure obtained by the measuring surface of the microphone array, and M is multiplied by 1; g represents a transfer matrix between a sound source surface and a microphone array measurement surface, and M is multiplied by N; q represents the source intensity of the sound source, and N is multiplied by 1; and is also provided with
r mn Represents the distance between the nth focus point and the mth array element, r n Indicating that the distance i between the nth focus point and the origin of coordinates is an imaginary unit, r mn Distance from the focusing point on the sound source surface to the array element on the measuring plane;
step B, calculating indoor impulse response based on mirror image source method
Specifically, an indoor reverberation mirror image source model is built based on a mirror image source method, and an indoor microphone array receiving signal is 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, τ represents a delay amount generated by sound wave reflection, and x represents convolution operation;
fourier transforming the above formula to a frequency domain to obtain Y (ω) =x (ω) H (ω) +n (ω);
thereby constructing an indoor impulse response function
Wherein beta is x,1 、β y,1 And beta z,1 Representing the reflection coefficient of the wall near the origin of coordinates in each direction; beta x,2 、β y,2 And beta z,2 Representing the reflection coefficient of the wall away from the origin of coordinates in each direction; r is R p Representing the distance from the actual sound source or the mirror image sound source to each array element; r is R r Representing the virtual space size corresponding to the multi-order reflection; u is the sound velocity; q (Q) 1 、Q 2 、Q 3 R is 0 or 1 p 8 positions can be taken;
the conversion of the room impulse response function to the frequency domain is:
for the actual sound source point coordinates (x, y, z), the sensor element coordinates (x ', y ', z ') are:
R p ＝(x-x′+2Q 1 x′,y-y′+2Q 2 y′,z-z′+2Q 3 z′)R r ＝2(n x L x ,n y L y ,n z L z )
wherein: n is n x ,n y ,n z Is 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: c=pp H ＝Gqq H G H The method comprises the steps of carrying out a first treatment on the surface of the H refers to conjugation treatment; g is a transfer matrix G; when the sound source is incoherent, the q is the pair H Simplifying off-diagonal elements in qq H Simplified intoThe sound pressure cross-spectrum matrix C is expressed as further expressed as: />Wherein: g n For the corresponding column vector in the transfer matrix G;
step D, calculating a point propagation function;
the traditional beam forming output result based on the cross spectrum function is as follows: b=ω H Cω＝ω H pp H ω；
Where b represents the acoustic power output by each focus point, ω represents the steering vector, ω= [ ω ] 1 ，ω 2 …ω M ]，Obtain->
Based on the foregoing, a point spread function is obtainedSound source surface acoustic power->
Step E, a step of original signal reconstruction
1) Initializing sparsity K 0 Initializing a sparsity estimation step l=s, s being a search step; the support set f=Φ, Φ 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 >I, j=1, 2, … N, where j is the j-th column of the sensing matrix and r is the residual;
extracting K from the correlation coefficient matrix 0 Corresponding to maximum valueThe index value is stored in F;
3) Constraint equidistant conditions based on one of the compressed sensing preconditions, ifAdaptive dilution k=k 0 +L, turning to step (2);
wherein y refers to sound field data measured by the acoustic sensor array; delta K Is a constant between 0 and 1, delta K =0.3; t refers to the transpose matrix; wherein Φ refers to the sensing matrix in compressed sensing;
if it isThen go to step 4;
4) Solving for initial margin by least square method
5) The initial estimation signal x=0, the number of initialization stages stage=1, the number of initialization iterations k=1, the initialization index value set s=phi, and the candidate set c=phi;
6) Using u= { u j |u j ＝|<r,φ j >New correlation coefficient matrix is calculated for j=1, 2, … N and is calculated as per u j The corresponding index value is stored into the index set S according to the standard of more than or equal to 0.5 max|u|;
7) The index value set t=f u S is merged by the formula u= { u j |u j ＝|＜r,φ j Calculating the correlation coefficient of the index value corresponding to the atom and the allowance in T according to the ratio of > |, j=1, 2 and … N, and extracting K 0 The index value corresponding to the maximum value is stored in F new In the method, based on the least square method, the method is adoptedCalculating an estimated signal x new And utilize r new ＝y-Φ F x updating the allowance;
8) Calculating the error between two iterations and determining if it is less than a specified errorIf the difference is the error between two iterations ||x new -x|| 2 Stopping iteration if epsilon is less than or equal to epsilon, otherwise turning to the step (9), wherein epsilon refers to a specified error, and x refers to a previous error value;
9) If r new || 2 ≥||r|| 2 Adding 1 to the number of stages, stage=stage+1, l=stage×s, and skipping step (6); if r new || 2 ＜||r|| 2 F=f new ，r＝r new K=k+1, go to step (6);
after stopping iteration, outputting the latest obtained estimated signal x new The method comprises the steps of carrying out a first treatment on the surface of the And after the algorithm iterates, outputting a result to obtain the beam data which is finally wanted to be obtained.
2. The method of claim 1, wherein the acoustic properties of each wall are represented by an acoustic absorption coefficient y, a reflection coefficientWherein gamma is the sound absorption coefficient of the wall surface.
3. The method for adaptive beamforming in a low signal-to-noise ratio and reverberant environment according to claim 1, wherein in the step B, in calculating an indoor impulse response based on a mirror image source method, the reflected sound wave is regarded as a direct sound wave from a mirror image sound source to a sensor when the indoor impulse response function is processed, and the sound path is calculated based on a positional relationship between the mirror image sound source and the sensor; and meanwhile, the energy attenuation of the sound wave in the reflection process is considered in the calculation process based on the reflection coefficient of each wall surface.
CN202110957118.3A 2021-08-19 2021-08-19 Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment Active CN113687307B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110957118.3A CN113687307B (en) 2021-08-19 2021-08-19 Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110957118.3A CN113687307B (en) 2021-08-19 2021-08-19 Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment

Publications (2)

Publication Number Publication Date
CN113687307A CN113687307A (en) 2021-11-23
CN113687307B true CN113687307B (en) 2023-08-18

Family

ID=78580789

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110957118.3A Active CN113687307B (en) 2021-08-19 2021-08-19 Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment

Country Status (1)

Country Link
CN (1) CN113687307B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114333876B (en) * 2021-11-25 2024-02-09 腾讯科技（深圳）有限公司 Signal processing method and device
CN114779203B (en) * 2022-06-20 2022-09-09 中国电子科技集团公司第五十四研究所 Target positioning method for wide-area random sparse array airspace scanning energy information matching
CN116148770B (en) * 2023-04-21 2023-07-07 湖南工商大学 Sound source positioning method, device and system based on array signal processing

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102183435A (en) * 2011-01-25 2011-09-14 中国船舶重工集团公司第七一五研究所 Method for measuring submarine density and sound velocity based on multi-path reflection theory
US9869752B1 (en) * 2016-04-25 2018-01-16 Ocean Acoustical Services And Instrumentation Systems, Inc. System and method for autonomous joint detection-classification and tracking of acoustic signals of interest
CN110082725A (en) * 2019-03-12 2019-08-02 西安电子科技大学 Auditory localization delay time estimation method, sonic location system based on microphone array
CN112180329A (en) * 2020-09-07 2021-01-05 黑龙江工程学院 Automobile noise source acoustic imaging method based on array element random uniform distribution spherical array deconvolution beam forming

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6594201B2 (en) * 2001-10-23 2003-07-15 Lockheed Martin Corporation System and method for localizing targets using multiple arrays

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102183435A (en) * 2011-01-25 2011-09-14 中国船舶重工集团公司第七一五研究所 Method for measuring submarine density and sound velocity based on multi-path reflection theory
US9869752B1 (en) * 2016-04-25 2018-01-16 Ocean Acoustical Services And Instrumentation Systems, Inc. System and method for autonomous joint detection-classification and tracking of acoustic signals of interest
CN110082725A (en) * 2019-03-12 2019-08-02 西安电子科技大学 Auditory localization delay time estimation method, sonic location system based on microphone array
CN112180329A (en) * 2020-09-07 2021-01-05 黑龙江工程学院 Automobile noise source acoustic imaging method based on array element random uniform distribution spherical array deconvolution beam forming

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
张宏宇 ; 郭文勇 ; 韩江桂 ; 陈汉涛 ; 夏菁 ; .基于镜像源方法抑制混响干扰的声全息算法改进.科学技术与工程.2019,(第33期),全文. *

Also Published As

Publication number Publication date
CN113687307A (en) 2021-11-23

Similar Documents

Publication Publication Date Title
CN113687307B (en) Self-adaptive wave beam forming method under low signal-to-noise ratio and reverberation environment
Varslot et al. Computer simulation of forward wave propagation in soft tissue
CN111947045B (en) GVMD parameter optimization and singular value decomposition-based fluid pipeline leakage positioning method
US20150057083A1 (en) Methods, systems, and computer readable media for simulating sound propagation in large scenes using equivalent sources
CN110109058A (en) A kind of planar array deconvolution identification of sound source method
CN113238189A (en) Sound source identification method and system based on array measurement and sparse prior information
Yan et al. Acoustic tomography system for online monitoring of temperature fields
CN111812581B (en) Spherical array sound source direction-of-arrival estimation method based on atomic norms
CN114460588B (en) High-precision imaging method based on active acoustic imager
Lingevitch et al. Parabolic equation simulations of reverberation statistics from non-Gaussian-distributed bottom roughness
Reed et al. Neural Volumetric Reconstruction for Coherent Synthetic Aperture Sonar
Kraut et al. A generalized Karhunen-Loeve basis for efficient estimation of tropospheric refractivity using radar clutter
Orofino et al. Efficient angular spectrum decomposition of acoustic sources. II. Results
Hu et al. Achieving high-resolution 3D acoustic imaging in a large-aspect-ratio cabin by the non-synchronous measurements
Waag et al. Estimates of wave front distortion from measurements of scattering by model random media and calf liver
CN110850421A (en) Underwater target detection method based on space-time adaptive processing of reverberation symmetric spectrum
Merino-Martınez et al. Three–dimensional acoustic imaging using asynchronous microphone array measurements
Dong et al. A study of multi-static ultrasonic tomography using propagation and back-propagation method
CN113126029A (en) Multi-sensor pulse sound source positioning method suitable for deep sea reliable acoustic path environment
CN116879952B (en) Calculation method, storage medium and equipment for point source elastic wave seabed reflection coefficient
Bazulin et al. Deconvolution of complex echo signals by the maximum entropy method in ultrasonic nondestructive inspection
CN117250604B (en) Separation method of target reflection signal and shallow sea reverberation
Dorausch et al. Adoption and evaluation of a multistatic Fourier-based synthetic aperture radar method for ultrasound imaging
CN114609493B (en) Partial discharge signal identification method with enhanced signal data
Ma et al. DAMAS with compression computational grid based on functional beamforming for acoustic source localization

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant