CN113687306A - Multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification method - Google Patents

Multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification method Download PDF

Info

Publication number
CN113687306A
CN113687306A CN202110580351.4A CN202110580351A CN113687306A CN 113687306 A CN113687306 A CN 113687306A CN 202110580351 A CN202110580351 A CN 202110580351A CN 113687306 A CN113687306 A CN 113687306A
Authority
CN
China
Prior art keywords
sound source
sound
source
coordinates
estimation
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.)
Granted
Application number
CN202110580351.4A
Other languages
Chinese (zh)
Other versions
CN113687306B (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.)
Chongqing University
Original Assignee
Chongqing University
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 Chongqing University filed Critical Chongqing University
Priority to CN202110580351.4A priority Critical patent/CN113687306B/en
Publication of CN113687306A publication Critical patent/CN113687306A/en
Application granted granted Critical
Publication of CN113687306B publication Critical patent/CN113687306B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H17/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Circuit For Audible Band Transducer (AREA)

Abstract

The invention discloses a multi-frequency synchronization two-dimensional off-network compressed beam forming sound source identification method, which comprises the following steps: 1) establishing a multi-frequency sound pressure signal measurement model; 2) establishing a maximum likelihood estimation model with the sound source position and the source intensity as parameters; 3) performing preliminary estimation on the sound source position to obtain a preliminary estimation coordinate of the sound source position; 4) carrying out local optimization on the preliminary estimation coordinates of the sound source position by utilizing a Newton method to obtain optimized estimation coordinates of the sound source position; 5) carrying out global cyclic optimization on all the recognized sound source position coordinates by utilizing a Newton method; 6) performing orthogonal solution on all the recognized sound source intensities, and updating residual errors; the method provided by the invention can effectively overcome the problem of base mismatch, is compatible with a planar array with randomly arranged microphones, is suitable for both steady-state and unsteady-state sound sources, and has high spatial resolution and strong anti-noise interference capability.

Description

Multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification method
Technical Field
The invention relates to the field of sound source identification, in particular to a multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification method.
Background
The beam forming sound source identification technology based on the remote measurement in the microphone array is widely applied to the fields of automobiles, airplanes, high-speed trains, wind power generation and the like. The two-dimensional compressed beam forming method using the planar microphone array and the compressed sensing theory is concerned with due to the advantages of wide sound source identification space range, applicability to coherent and incoherent sound sources, clear imaging and the like, and is a current research hotspot. The traditional compressed beam forming method discretizes a sound source imaging area, establishes an underdetermined equation set between source intensity distribution and microphone measurement sound pressure, and solves to obtain position and intensity estimation of a sound source. However, this method has a fundamental mismatch problem and a problem of poor sound source identification performance at small separations or low signal-to-noise ratios. The base mismatch problem refers to: the sound source recognition performance is degraded when the sound source is deviated from the grid points. Although the encryption grid can alleviate the problem to some extent, the cost is increased, and the excessive encryption grid causes excessive coherence of the sensing matrix column (steering vector), which causes algebraic aliasing and degrades the sound source identification performance. Poor sound source recognition performance at small separations or low signal-to-noise ratios is also associated with the perceptual matrix array.
In order to solve the problem of base mismatch, many scholars make efforts, and three methods, namely a non-grid method, a dynamic grid method and an off-grid method are proposed to effectively solve the problem of base mismatch. To improve the sound source identification performance under small separation or low signal-to-noise ratio conditions, many scholars have attempted to use a multi-snapshot approach. The multi-snapshot method is based on time invariant hypothesis sampling of sound source signals to obtain a plurality of array signal data blocks (snapshots), all the snapshots share the same sound source sparsity, and source intensity distribution is obtained through joint sparse recovery processing. The joint processing effectively promotes the common sparsity in space and time, improves the anti-noise interference capability, and increases the spatial resolution, the sound source position estimation and the source intensity quantization precision due to the increase of the microphone sound pressure matrix rank. Whether conventional fixed-grid, non-grid, or off-grid compressed beamforming, corresponding multi-snapshot versions have been developed. However, the multi-shot method is limited by the time-invariant assumption, and is only applicable to a stationary sound source, and is not applicable to an unsteady sound source.
In summary, there has not been a method that can simultaneously overcome the above two problems, and can be applied to a planar array with arbitrarily arranged microphones and can also be applied to various sound sources (steady state, unsteady state).
Disclosure of Invention
The invention aims to provide a multi-frequency synchronization two-dimensional off-network compressed beam forming sound source identification method, which comprises the following steps:
1) and establishing a multi-frequency sound pressure signal measurement model.
The m-th microphone has the coordinate of
Figure BDA0003085921650000021
The coordinates of the s-th sound source are
Figure BDA0003085921650000022
Wherein the content of the first and second substances,
Figure BDA0003085921650000023
the x-axis coordinate and the y-axis coordinate of the mth microphone.
Figure BDA0003085921650000024
The x-axis coordinate and the y-axis coordinate of the s-th sound source are shown. h is the distance between the sound source plane and the array plane. M is 1,2, 3 …, M. S is 1,2, 3 …, S. Theoretical sound pressure at the location of the mth microphone
Figure BDA0003085921650000025
As follows:
Figure BDA0003085921650000026
in the form of matrix
Figure BDA0003085921650000027
Matrix array
Figure BDA0003085921650000028
The superscript T denotes transpose. f. oflThe frequency corresponding to the ith frequency point. q. q.ss,lIndicating the intensity of the s-th sound source at the l-th frequency point.
Figure BDA0003085921650000029
Is an imaginary unit. And c represents the speed of sound. dm,sIndicating the distance of the mth microphone from the mth sound source. I | · | purple wind2Representing the two-norm of the vector.
The multi-frequency sound pressure signal measurement model is as follows:
p＝A(ΩS)qS+n (2)
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500000210
and a source intensity vector formed by the source intensities of S sound sources at each frequency.
Figure BDA00030859216500000211
The source intensity vector for all frequencies of the s-th source.
Figure BDA00030859216500000212
Representing the microphone measurement noise vector. n is·,l＝[n1,l,n2,l,...,nM,l]Representing the measured noise vector of the microphone at the point corresponding to the l-th frequency point.
Figure BDA00030859216500000213
Representing the transfer matrix between S sound sources to all microphones at all L frequencies.
Figure BDA00030859216500000214
A coordinate matrix composed of coordinates of all sound source positions,
Figure BDA00030859216500000215
a complex set is represented.
Figure BDA00030859216500000216
Representing a set of real numbers. Wherein the transfer matrix between the s-th sound source to all microphones at all L frequencies
Figure BDA00030859216500000217
As follows:
Figure BDA00030859216500000218
in the formula, the frequency flIs first ofTransfer vector between s sound sources to all microphones
Figure BDA00030859216500000219
2) And establishing a maximum likelihood estimation model with the sound source position and the source intensity as parameters.
The step of establishing a maximum likelihood estimation model with the sound source position and the source intensity as parameters comprises the following steps:
2.1) minimizing the residual error of the multi-frequency Acoustic pressure Signal measurement model
Figure BDA00030859216500000220
Obtaining maximum likelihood estimates of sound source position and source intensity
Figure BDA0003085921650000031
Namely:
Figure BDA0003085921650000032
in the formula (I), the compound is shown in the specification,
Figure BDA0003085921650000033
the representation takes the real part of the complex number. The superscript H denotes the hermitian conjugate transpose operator.
2.2) establishing a likelihood function
Figure BDA0003085921650000034
Namely:
Figure BDA0003085921650000035
2.3) Source intensity vector q formed by the Source intensities of S Sound sources at respective frequenciesSFor variables, a least squares solution is established, i.e.:
qS＝(A(ΩS)HA(ΩS))-1A(ΩS)Hp (6)
2.4) establishing a maximum likelihood ratio test cost function
Figure BDA0003085921650000036
Namely:
Figure BDA0003085921650000037
2.5) establishing an estimation equation of the sound source position coordinates, namely:
Figure BDA0003085921650000038
in the formula (I), the compound is shown in the specification,
Figure BDA0003085921650000039
coordinates are estimated for the sound source position.
3) And performing preliminary estimation on the sound source position to obtain a preliminary estimation coordinate of the sound source position.
The step of preliminary estimating the sound source position includes:
3.1) spatially discretizing the sound source detection into N grid points. N > M. Wherein the coordinates of the nth grid point are
Figure BDA00030859216500000310
Set of grid point coordinates
Figure BDA00030859216500000311
Wherein, the matrix
Figure BDA00030859216500000312
3.2) establishing a perceptual matrix of all frequency correspondences between the N grid points to the M microphones
Figure BDA00030859216500000313
Wherein, the sensing matrix corresponding to all frequencies from the nth grid point to the M microphones
Figure BDA00030859216500000314
3.3) the actual measured microphone sound pressure vector p can be expressed as:
p＝A(ΩG)qG+n (9)
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500000315
and a source intensity vector formed by source intensities of all frequencies of the N grid point positions.
Figure BDA00030859216500000316
Source strength vector for all frequencies of nth grid point
Figure BDA0003085921650000041
Representing the frequency f at the nth grid pointlThe source is strong.
3.4) updating the likelihood function, the maximum likelihood ratio test cost function and the least square function into the identification formula of the single sound source to obtain the function
Figure BDA0003085921650000042
Sum function
Figure BDA0003085921650000043
Namely:
Figure BDA0003085921650000044
Figure BDA0003085921650000045
Figure BDA0003085921650000046
3.5) from the set ΩGIn selecting a maximum likelihood ratio test cost function
Figure BDA0003085921650000047
Maximum mesh node as sound source positionPreliminary estimated coordinates, denoted
Figure BDA0003085921650000048
Preliminary estimation coordinates of sound source position
Figure BDA0003085921650000049
As follows:
Figure BDA00030859216500000410
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500000411
representing the residual vector of the microphone sound pressure signal obtained after the first s-1 sound source estimates.
Figure BDA00030859216500000412
And
Figure BDA00030859216500000413
respectively the estimated coordinates and source strength of the jth sound source.
3.6) preliminary estimation of coordinates of the sound source position
Figure BDA00030859216500000414
Substituting the following formula to obtain the source intensity estimation vector of all frequencies of the s sound source, namely:
Figure BDA00030859216500000415
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500000416
vectors are estimated for the source strengths of all frequencies of the s-th source.
4) And carrying out local optimization on the preliminary estimation coordinates of the sound source position by utilizing a Newton method to obtain the optimized estimation coordinates of the sound source position.
The method for carrying out local optimization on the initial estimation coordinates of the sound source position by utilizing the Newton method comprises the following steps:
4.1) establishing a two-dimensional Newton optimization function, namely:
Figure BDA00030859216500000417
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500000418
and
Figure BDA00030859216500000419
respectively representing a jacobian matrix and a hessian matrix.
Figure BDA00030859216500000420
The estimated coordinates are optimized for the sound source position.
Jacobian matrix
Figure BDA00030859216500000421
Heisen matrix
Figure BDA00030859216500000422
Respectively as follows:
Figure BDA00030859216500000423
Figure BDA0003085921650000051
wherein the parameter item
Figure BDA0003085921650000052
Parameter item
Figure BDA0003085921650000053
Parameter item
Figure BDA0003085921650000054
Parameter item
Figure BDA0003085921650000055
Parameter item
Figure BDA0003085921650000056
Parameter item
Figure BDA0003085921650000057
Respectively as follows:
Figure BDA0003085921650000058
Figure BDA0003085921650000059
Figure BDA00030859216500000510
Figure BDA00030859216500000511
Figure BDA00030859216500000512
Figure BDA00030859216500000513
in the formula, AfRepresenting a transfer matrix at a frequency f;
4.2) optimizing and estimating coordinates of sound source position
Figure BDA00030859216500000514
Substituting the obtained value into the following formula to obtain the sound source intensity optimized estimation
Figure BDA00030859216500000515
Namely:
Figure BDA00030859216500000516
4.3) updating the microphone pressure residual vector prNamely:
Figure BDA00030859216500000517
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500000518
the transfer matrix at all frequencies between the jth source to all M microphones.
5) And carrying out global cyclic optimization on the coordinates of all the identified sound source positions by utilizing a Newton method to obtain global optimization estimated coordinates of all the identified sound source positions.
The method for carrying out global cyclic optimization on all the recognized sound source position coordinates by utilizing the Newton method comprises the following steps:
5.1) establishing a source intensity estimation vector
Figure BDA0003085921650000061
Transmission matrix
Figure BDA0003085921650000062
And harmonic source location estimation matrices
Figure BDA0003085921650000063
5.2) carrying out the t-th iterative optimization on the kth sound source, wherein the steps are as follows:
5.2.1) computing the transition residual vector that preserves the kth source component at the t-1 th iteration
Figure BDA0003085921650000064
Namely:
Figure BDA0003085921650000065
in the formula, t is initially 1.
5.2.2) residual transition
Figure BDA0003085921650000066
Substituting into the following formula to obtain the transition source intensity
Figure BDA0003085921650000067
Namely:
Figure BDA0003085921650000068
5.2.3) calculating the global optimization estimation coordinates of the sound source position
Figure BDA0003085921650000069
Namely:
Figure BDA00030859216500000610
updating microphone sound pressure residual vector
Figure BDA00030859216500000611
Namely:
Figure BDA00030859216500000612
5.3) determining the residual energy change of two continuous cycles
Figure BDA00030859216500000613
If yes, entering step 6), otherwise, making t equal to t +1, and returning to step 5.2);
6) the steps of solving all the identified sound source intensities orthogonally and updating the residual errors include:
6.1) optimizing the estimation coordinate set of the sound source position
Figure BDA00030859216500000614
Substituting into the following formula to obtain the sound source intensity estimation
Figure BDA00030859216500000615
Namely:
Figure BDA00030859216500000616
6.2) updating the microphone pressure residual vector prNamely:
Figure BDA00030859216500000617
6.3) judging whether k is more than or equal to S, if so, ending the circulation and outputting a set
Figure BDA00030859216500000618
And
Figure BDA00030859216500000619
otherwise, let k be k +1 and return to step 3.5).
The technical effect of the invention is undoubted, the method firstly establishes a maximum likelihood estimation model taking the sound source position and the intensity of each frequency source as parameters, and then utilizes the NOMP algorithm and the block sparse thought to form a BNOMP algorithm for solving to obtain the estimation of the sound source position and the intensity of each frequency source. The method provided by the invention can effectively overcome the problem of base mismatching, is compatible with a planar array with randomly arranged microphones, and is suitable for both steady-state and unsteady-state sound sources. In addition, the performance of the test piece is analyzed based on simulation and verification tests, and the result shows that: the method has high spatial resolution and strong anti-noise interference capability due to suppression of mutual interference of sound sources by global loop feedback optimization, promotion of multi-frequency joint sparse constraint on sparsity and suppression of noise interference.
Drawings
FIG. 1 is a planar microphone array measurement model;
fig. 2 shows the result of gaussian pulse signal identification (SNR 20 dB);
FIG. 2(a) is a sound source localization cloud diagram obtained by superimposing the energy of the estimated positions of the 1000-and 5000-Hz sound source according to the conventional orthogonal matching pursuit compressed beam forming method with single frequency processing;
FIG. 2(b) is a sound source localization cloud diagram obtained by superimposing the energy of the estimated positions of the 1000-5000Hz sound source according to the conventional single-frequency-processing Newton orthogonal matching tracking compressed beam forming method;
FIG. 2(c) is a sound source localization cloud diagram of an orthogonal matching pursuit compression beamforming method of multi-frequency synchronization processing;
FIG. 2(d) is a sound source localization cloud of the present invention;
fig. 2(e) shows the source strength quantization result of the conventional single-frequency orthogonal matching pursuit compressed beam forming method at each frequency;
FIG. 2(f) is a diagram illustrating the source strength quantization results of the conventional single-frequency Newton orthogonal matching pursuit compression beamforming method at each frequency;
FIG. 2(g) is a diagram illustrating the source-intensive quantization result of the conventional orthogonal matching pursuit compressed beamforming method with multi-frequency synchronization processing;
FIG. 2(h) is the source intensity quantization result at each frequency according to the present invention;
fig. 3 shows the result of gaussian pulse signal identification (SNR 5 dB);
FIG. 3(a) is a sound source localization cloud diagram obtained by superimposing the energy of the estimated positions of the 1000-and 5000-Hz sound source according to the conventional orthogonal matching pursuit compressed beam forming method with single frequency processing;
FIG. 3(b) is a sound source localization cloud chart obtained by superimposing the energy of the estimated positions of the 1000-5000Hz sound source according to the conventional single-frequency-processing Newton orthogonal matching tracking compressed beam forming method;
FIG. 3(c) is a sound source localization cloud diagram of the orthogonal matching pursuit compression beam forming method of multi-frequency synchronization processing;
FIG. 3(d) is a sound source localization cloud of the present invention;
fig. 3(e) shows the source strength quantization result of the conventional single-frequency orthogonal matching pursuit compressed beamforming method at each frequency;
FIG. 3(f) is a diagram illustrating the source strength quantization results of the conventional single-frequency Newton orthogonal matching pursuit compressed beamforming method at each frequency;
FIG. 3(g) is a diagram illustrating the source-intensive quantization result of the conventional orthogonal matching pursuit compressed beamforming method with multi-frequency synchronization processing;
FIG. 3(h) is the source intensity quantization result at each frequency according to the present invention;
fig. 4 shows the test arrangement and the amplitude spectra of the measured signals of the microphones (a) the test arrangement and (b) the amplitude spectra of the measured signals of the microphones.
Fig. 5 shows the test results.
FIG. 5(a) is a sound source localization cloud diagram obtained by superimposing the energy of the estimated positions of the 1000-and 5000-Hz sound source according to the conventional orthogonal matching pursuit compressed beam forming method with single frequency processing;
FIG. 5(b) is a sound source localization cloud chart obtained by superimposing the energy of the estimated positions of the 1000-5000Hz sound source according to the conventional single-frequency-processing Newton orthogonal matching tracking compressed beam forming method;
FIG. 5(c) is a sound source localization cloud diagram of the orthogonal matching pursuit compression beamforming method of multi-frequency synchronization processing;
FIG. 5(d) is a sound source localization cloud of the present invention;
Detailed Description
The present invention is further illustrated by the following examples, but it should not be construed that the scope of the above-described subject matter is limited to the following examples. Various substitutions and alterations can be made without departing from the technical idea of the invention and the scope of the invention is covered by the present invention according to the common technical knowledge and the conventional means in the field.
Example 1:
referring to fig. 1 to 5, the multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification method includes the following steps:
1) and establishing a multi-frequency sound pressure signal measurement model.
The m-th microphone has the coordinate of
Figure BDA0003085921650000081
The coordinates of the s-th sound source are
Figure BDA0003085921650000082
Wherein the content of the first and second substances,
Figure BDA0003085921650000083
the x-axis coordinate and the y-axis coordinate of the mth microphone.
Figure BDA0003085921650000084
The x-axis coordinate and the y-axis coordinate of the s-th sound source are shown. h is the distance between the sound source plane and the array plane. M is 1,2, 3 …, M. S is 1,2, 3 …, S. Theoretical sound pressure at the location of the mth microphone
Figure BDA0003085921650000085
As follows:
Figure BDA0003085921650000086
in the form of matrix
Figure BDA0003085921650000091
Matrix array
Figure BDA0003085921650000092
The superscript T denotes transpose. f. oflThe frequency corresponding to the ith frequency point. q. q.ss,lIndicating the intensity of the s-th sound source at the l-th frequency point.
Figure BDA0003085921650000093
Is an imaginary unit. And c represents the speed of sound. dm,sIndicating the distance of the mth microphone from the mth sound source. I | · | purple wind2Representing the two-norm of the vector.
The multi-frequency sound pressure signal measurement model is as follows:
p＝A(ΩS)qS+n (2)
in the formula (I), the compound is shown in the specification,
Figure BDA0003085921650000094
and a source intensity vector formed by the source intensities of S sound sources at each frequency.
Figure BDA0003085921650000095
The source intensity vector for all frequencies of the s-th source.
Figure BDA0003085921650000096
Representing the microphone measurement noise vector. n is·,l＝[n1,l,n2,l,...,nM,l]Representing the measured noise vector of the microphone at the point corresponding to the l-th frequency point.
Figure BDA0003085921650000097
Representing the transfer matrix between S sound sources to all microphones at all L frequencies.
Figure BDA0003085921650000098
A coordinate matrix composed of coordinates of all sound source positions,
Figure BDA0003085921650000099
a complex set is represented.
Figure BDA00030859216500000910
Representing a set of real numbers. Wherein the transfer matrix between the s-th sound source to all microphones at all L frequencies
Figure BDA00030859216500000911
As follows:
Figure BDA00030859216500000912
in the formula, the frequency flTransfer vector between the s-th sound source to all microphones
Figure BDA00030859216500000913
2) And establishing a maximum likelihood estimation model with the sound source position and the source intensity as parameters.
The step of establishing a maximum likelihood estimation model with the sound source position and the source intensity as parameters comprises the following steps:
2.1) minimizing the residual error of the multi-frequency Acoustic pressure Signal measurement model
Figure BDA00030859216500000914
Obtaining maximum likelihood estimates of sound source position and source intensity
Figure BDA00030859216500000915
Namely:
Figure BDA00030859216500000916
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500000917
the representation takes the real part of the complex number. The superscript H denotes the hermitian conjugate transpose operator.
2.2) establishing a likelihood function
Figure BDA00030859216500000918
Namely:
Figure BDA00030859216500000919
2.3) Source intensity vector q formed by the Source intensities of S Sound sources at respective frequenciesSFor variables, a least squares solution is established, i.e.:
qS＝(A(ΩS)HA(ΩS))-1A(ΩS)Hp (6)
2.4) establishing a maximum likelihood ratio test cost function
Figure BDA0003085921650000101
Namely:
Figure BDA0003085921650000102
2.5) establishing an estimation equation of the sound source position coordinates, namely:
Figure BDA0003085921650000103
in the formula (I), the compound is shown in the specification,
Figure BDA0003085921650000104
coordinates are estimated for the sound source position.
3) And performing preliminary estimation on the sound source position to obtain a preliminary estimation coordinate of the sound source position.
The step of preliminary estimating the sound source position includes:
3.1) spatially discretizing the sound source detection into N grid points. N > M. Wherein the coordinates of the nth grid point are
Figure BDA0003085921650000105
Set of grid point coordinates
Figure BDA0003085921650000106
Wherein, the matrix
Figure BDA0003085921650000107
3.2) establishing a perceptual matrix of all frequency correspondences between the N grid points to the M microphones
Figure BDA0003085921650000108
Wherein, the sensing matrix corresponding to all frequencies from the nth grid point to the M microphones
Figure BDA0003085921650000109
3.3) the actual measured microphone sound pressure vector p can be expressed as:
p＝A(ΩG)qG+n (9)
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500001010
and a source intensity vector formed by source intensities of all frequencies of the N grid point positions.
Figure BDA00030859216500001011
Source strength vector for all frequencies of nth grid point
Figure BDA00030859216500001012
Representing the frequency f at the nth grid pointlThe source is strong.
3.4) updating the likelihood function, the maximum likelihood ratio test cost function and the least square function into the identification formula of the single sound source to obtain the function
Figure BDA00030859216500001013
Sum function
Figure BDA00030859216500001014
Namely:
Figure BDA00030859216500001015
Figure BDA00030859216500001016
Figure BDA00030859216500001017
3.5) from the set ΩGIn selecting a maximum likelihood ratio test cost function
Figure BDA00030859216500001018
The largest grid node is taken as the initial estimation coordinate of the sound source position and is recorded as
Figure BDA00030859216500001019
Preliminary estimation coordinates of sound source position
Figure BDA0003085921650000111
As follows:
Figure BDA0003085921650000112
in the formula (I), the compound is shown in the specification,
Figure BDA0003085921650000113
representing the residual vector of the microphone sound pressure signal obtained after the first s-1 sound source estimates.
Figure BDA0003085921650000114
And
Figure BDA0003085921650000115
respectively the estimated coordinates and source strength of the jth sound source.
3.6) preliminary estimation of coordinates of the sound source position
Figure BDA0003085921650000116
Substituting the following formula to obtain the source intensity estimation vector of all frequencies of the s sound source, namely:
Figure BDA0003085921650000117
in the formula (I), the compound is shown in the specification,
Figure BDA0003085921650000118
vectors are estimated for the source strengths of all frequencies of the s-th source.
4) And carrying out local optimization on the preliminary estimation coordinates of the sound source position by utilizing a Newton method to obtain the optimized estimation coordinates of the sound source position.
The method for carrying out local optimization on the initial estimation coordinates of the sound source position by utilizing the Newton method comprises the following steps:
4.1) establishing a two-dimensional Newton optimization function, namely:
Figure BDA0003085921650000119
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500001123
and
Figure BDA00030859216500001110
respectively representing a jacobian matrix and a hessian matrix.
Figure BDA00030859216500001111
The estimated coordinates are optimized for the sound source position.
Jacobian matrix
Figure BDA00030859216500001112
Heisen matrix
Figure BDA00030859216500001113
Respectively as follows:
Figure BDA00030859216500001114
Figure BDA00030859216500001115
wherein the parameter item
Figure BDA00030859216500001116
Parameter item
Figure BDA00030859216500001117
Parameter item
Figure BDA00030859216500001118
Parameter item
Figure BDA00030859216500001119
Parameter item
Figure BDA00030859216500001120
Parameter item
Figure BDA00030859216500001121
Respectively as follows:
Figure BDA00030859216500001122
Figure BDA0003085921650000121
Figure BDA0003085921650000122
Figure BDA0003085921650000123
Figure BDA0003085921650000124
Figure BDA0003085921650000125
in the formula, AfRepresenting a transfer matrix at a frequency f;
4.2) optimizing and estimating coordinates of sound source position
Figure BDA0003085921650000126
Substituting the obtained value into the following formula to obtain the sound source intensity optimized estimation
Figure BDA0003085921650000127
Namely:
Figure BDA0003085921650000128
4.3) updating the microphone pressure residual vector prNamely:
Figure BDA0003085921650000129
in the formula (I), the compound is shown in the specification,
Figure BDA00030859216500001210
the transfer matrix at all frequencies between the jth source to all M microphones.
5) And carrying out global cyclic optimization on the coordinates of all the identified sound source positions by utilizing a Newton method to obtain global optimization estimated coordinates of all the identified sound source positions.
The method for carrying out global cyclic optimization on all the recognized sound source position coordinates by utilizing the Newton method comprises the following steps:
5.1) establishing a source intensity estimation vector
Figure BDA00030859216500001211
Transmission matrix
Figure BDA00030859216500001212
And harmonic source location estimation matrices
Figure BDA00030859216500001213
5.2) carrying out the t-th iterative optimization on the kth sound source, wherein the steps are as follows:
5.2.1) computing the transition residual vector that preserves the kth source component at the t-1 th iteration
Figure BDA00030859216500001214
Namely:
Figure BDA00030859216500001215
in the formula, t is initially 1.
5.2.2) residual transition
Figure BDA0003085921650000131
Substituting into the following formula to obtain the transition source intensity
Figure BDA0003085921650000132
Namely:
Figure BDA0003085921650000133
5.2.3) calculating the global optimization estimation coordinates of the sound source position
Figure BDA0003085921650000134
Namely:
Figure BDA0003085921650000135
the arrow in the formula indicates that the value of the left symbol is updated to the right expression value.
5.2.4) updating microphone pressure residual vector
Figure BDA0003085921650000136
Namely:
Figure BDA0003085921650000137
3) determining the residual energy change of two consecutive cycles
Figure BDA0003085921650000138
If yes, entering step 6), otherwise, making t equal to t +1, and returning to step 5.2);
6. the steps of solving all the identified sound source intensities orthogonally and updating the residual errors include:
6.1) optimizing the estimation coordinate set of the sound source position
Figure BDA0003085921650000139
Substituting into the following formula to obtain sound source intensityDegree estimation
Figure BDA00030859216500001310
Namely:
Figure BDA00030859216500001311
6.2) updating the microphone pressure residual vector prNamely:
Figure BDA00030859216500001312
6.3) judging whether k is more than or equal to S, if so, ending the circulation and outputting a set
Figure BDA00030859216500001313
And
Figure BDA00030859216500001314
otherwise, let k be k +1 and return to step 3.5).
Example 2:
a multi-frequency synchronization two-dimensional off-network compressed beam forming sound source identification method comprises the following steps:
1) near-field measurements were performed on S broadband sound sources using a random microphone array with M microphones, the layout being shown in fig. 1, where "·" represents a microphone and "·" represents a point sound source. Let the cartesian coordinate system origin be at the array center and the XOY plane coincide with the array plane, assuming that the sound source plane is parallel to the array plane and at a distance h. The m-th microphone has the coordinate of
Figure BDA00030859216500001315
And order
Figure BDA00030859216500001316
The coordinates of the s-th sound source are
Figure BDA00030859216500001317
And order
Figure BDA00030859216500001318
The superscript "T" denotes the transpose operator. Considering a plurality of frequencies, the ith frequency point is denoted as flThe corresponding frequency f at the mth microphone positionlTheoretical sound pressure p ofm,lCan be expressed as:
Figure BDA00030859216500001319
in the formula, qs,lThe intensity of the ith frequency of the s-th sound source (expressed by the sound pressure generated by the sound source at 1m from it),
Figure BDA0003085921650000141
is an imaginary unit, c represents the speed of sound, dm,sRepresents the distance between the mth microphone and the mth sound source, | · | | survival2Representing the two-norm of the vector. From equation (1), the frequency f can be constructedlTransfer vector between the s-th sound source to all microphones
Figure BDA0003085921650000142
(ii) a Further, the transfer matrix between the s-th sound source to all microphones at all L frequencies may be constructed
Figure BDA0003085921650000143
The transfer matrix between all S sound sources to all M microphones at all L frequencies can be expressed as:
Figure BDA0003085921650000144
wherein the content of the first and second substances,
Figure BDA0003085921650000145
a coordinate matrix composed of all S sound source position coordinates,
Figure BDA0003085921650000146
the complex set is represented as a complex set,
Figure BDA0003085921650000147
representing a set of real numbers, superscript "ML×1、ML×SL"etc. represent matrix dimensions.
In consideration of noise interference, sound pressure components of sound pressure signals received by all M microphones at L frequencies form microphone sound pressure quantities in actual measurement
Figure BDA0003085921650000148
(wherein,
Figure BDA0003085921650000149
indicating that all microphones are at frequency flMeasured sound pressure, pm,lIndicating the m-th microphone at frequency flThe actual measured sound pressure signal) can be expressed as:
p＝A(ΩS)qS+n (2)
wherein the content of the first and second substances,
Figure BDA00030859216500001410
a source intensity vector representing the source intensities of all S sound sources at the respective frequencies,
Figure BDA00030859216500001411
the source intensity vectors for all frequencies of the s-th source,
Figure BDA00030859216500001412
representing the microphone measurement noise vector. Defining a frequency of flSignal to noise ratio of time is SNRl＝20log10(||p·,l-n·,l||2/||n·,l||2)。
2) Mathematical model building
By minimizing the residual error of the multi-frequency measurement model (equation (2)))
Figure BDA00030859216500001413
Maximum likelihood estimates of sound source position and source intensity can be obtained:
Figure BDA00030859216500001414
in the formula (I), the compound is shown in the specification,
Figure BDA0003085921650000151
representing the real part of the complex number, superscript "H"denotes the Hermite conjugate transpose operator. Defining the likelihood function as:
Figure BDA0003085921650000152
synchronous solving is carried out so that
Figure BDA0003085921650000153
Maximum omegaSAnd q isSDifficult to realize, so only qSAs variables, their least squares solution can be obtained:
qS＝(A(ΩS)HA(ΩS))-1A(ΩS)Hp (5)
substituting equation (6) to obtain the maximum likelihood ratio test cost function
Figure BDA0003085921650000154
Figure BDA0003085921650000155
Estimation of the coordinates of the sound source position
Figure BDA0003085921650000156
Comprises the following steps:
Figure BDA0003085921650000157
estimating coordinates
Figure BDA0003085921650000158
Substituting equation (7) can obtain the corresponding source strength estimation
Figure BDA0003085921650000159
3) Preliminary estimation of sound source position coordinates
Direct simultaneous estimation in the entire two-dimensional continuous plane makes all sound sources with the largest cost function difficult to implement. Therefore, the solution process is simplified in two ways: (1) discretizing the continuous plane to obtain a rough sound source position through preliminary estimation; (2) all sound sources are solved one by one (steps 2, 3, 4, 5 below all describe the identification and updating for the s-th sound source).
For this purpose, the target sound source area is discretized into N (N > M) grid points, and the coordinates of the nth grid point are
Figure BDA00030859216500001510
Order to
Figure BDA00030859216500001511
Constructing a set of grid point coordinates
Figure BDA00030859216500001512
Combining the sensing matrices of all L frequencies from all N grid points to all M microphones
Figure BDA00030859216500001513
(wherein,
Figure BDA00030859216500001514
then, the actual measured microphone sound pressure vector p can be expressed as:
p＝A(ΩG)qG+n (7)
wherein the content of the first and second substances,
Figure BDA00030859216500001515
a source intensity vector formed by the source intensities of all L frequencies representing all N grid point positions (where,
Figure BDA00030859216500001516
the source intensity vectors for all L frequencies at the nth grid point,
Figure BDA00030859216500001517
representing the frequency f at the nth grid pointlStrong source of time).
In addition, to solve all sound sources one by one, the likelihood function, the cost function, and the least square function for source intensity estimation are respectively degenerated into identification formulas for a single sound source (s-th sound source as an example):
Figure BDA0003085921650000161
Figure BDA0003085921650000162
Figure BDA0003085921650000163
through two simplifications, the problem then turns into: from set omega based on multi-frequency synchronous networked compressed beamforming model (equation (4))GIn selecting a maximum likelihood ratio test cost function
Figure BDA0003085921650000164
Preliminary estimation of maximum mesh node as sound source position coordinate
Figure BDA0003085921650000165
I.e. from the block containing N a priori blocks A (r)Gn) N1.., N, a perceptual matrix a (Ω)G) In search of cost function
Figure BDA0003085921650000166
Maximum one prior block:
estimating the position coordinates of the sound source:
Figure BDA0003085921650000167
wherein the content of the first and second substances,
Figure BDA0003085921650000168
representing the residual vector of the microphone sound pressure signal obtained after the first s-1 sound source estimates.
Figure BDA0003085921650000169
And
Figure BDA00030859216500001610
respectively the estimated coordinates and source strength of the jth sound source. Preliminarily estimating the coordinates of the sound source position
Figure BDA00030859216500001611
Substituting into formula (10) can obtain the source intensity estimation vector corresponding to all frequencies of the sound source
Figure BDA00030859216500001612
4) The Newton method carries out local optimization on the position coordinate of the s-th sound source which is roughly estimated, and carries out Newton optimization on the initial estimation in a local continuous two-dimensional plane, so that the likelihood function can be further converged to the current local maximum value, and more accurate sound source position and intensity estimation can be further obtained. Specifically, two-dimensional newton optimization is expressed as follows:
Figure BDA00030859216500001613
wherein the content of the first and second substances,
Figure BDA00030859216500001614
and
Figure BDA00030859216500001615
the jacobian matrix and the hessian matrix respectively have the following specific expressions:
Figure BDA00030859216500001616
Figure BDA0003085921650000171
wherein the content of the first and second substances,
Figure BDA0003085921650000172
Figure BDA0003085921650000173
Figure BDA0003085921650000174
Figure BDA0003085921650000175
Figure BDA0003085921650000176
Figure BDA0003085921650000177
obtaining an optimized estimate of the location coordinates of a sound source
Figure BDA0003085921650000178
The sound source intensity estimates may then be updated in sequence according to equation (10)
Figure BDA0003085921650000179
By using
Figure BDA00030859216500001710
Updating microphone sound pressure residual vector pr。
5) Global circulation Newton optimization for improving positioning precision
Global loop optimization provides feedback on all identified sound sources according to the current residuals, giving them the opportunity to re-optimize.Initialization based on all identified sound sources
Figure BDA00030859216500001711
And
Figure BDA00030859216500001712
(where k 1, 2.. and s are the identified sound source indices), a sound source position estimation matrix is constructed
Figure BDA00030859216500001713
Source strength estimation vector
Figure BDA00030859216500001714
And a transfer matrix
Figure BDA00030859216500001715
Cyclically subjecting it to Newton optimization in turn and updating the set
Figure BDA00030859216500001716
Where t represents the number of cycles. The process of optimizing the kth sound source in the t-th cycle is described as follows:
to optimize the kth source, a transition residual vector is first defined that retains the last kth source component by subtracting the removed component
Figure BDA00030859216500001717
The projection on the corresponding basis of the kth sound source is obtained back by adding, i.e.
Figure BDA00030859216500001718
Residual of transition
Figure BDA0003085921650000181
Substituting equation (10) can obtain the transition source intensity
Figure BDA0003085921650000182
At this time, the optimized sound source position can be obtained by further using the Newton's optimization formula (14)
Figure BDA0003085921650000183
Source strengths can be updated using equation (10)
Figure BDA0003085921650000184
And updating the residual
Figure BDA0003085921650000185
And the above-mentioned set
Figure BDA0003085921650000186
Figure BDA0003085921650000187
Circularly and sequentially optimizing each identified sound source until residual energy changes of two continuous cycles
Figure BDA0003085921650000188
6) Source-strong orthogonal solution
Because the group formed by the atoms corresponding to each sound source is not a group of orthogonal groups, the obtained residual error always has a component which is not orthogonal to the atoms corresponding to the optimized identified sound source. To remove these components, an orthogonal solution update of all frequency source intensities of the identified sound source is required. Based on the set formed by terminating after undergoing T times of global loop optimization in the step (4)
Figure BDA0003085921650000189
The source intensity of all frequency components of all identified sound sources can be updated using equation (5) and using equation
Figure BDA00030859216500001810
Updating residual pr。
All S sound source position estimates can be obtained through S times of four-step circulation
Figure BDA00030859216500001811
And all frequency source strength estimation thereof
Figure BDA00030859216500001812
Example 3:
simulation and experimental verification of the multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification method comprise the following steps:
1) in order to analyze and verify the sound source identification performance of the method (marked as BNOMP-CB), simulation tests are carried out on the method, and the method is compared with the existing orthogonal matching pursuit compressed beam forming method (marked as OMP-CB) of single-frequency processing, the Newton orthogonal matching pursuit compressed beam forming method (marked as NOMP-CB) and the orthogonal matching pursuit compressed beam forming method (marked as BOMP-CB) of multi-frequency synchronous processing.
The measurement array is the Bruel&
Figure BDA00030859216500001813
(
Figure BDA00030859216500001814
Denmark) 36-channel sector wheel array with a diameter of 0.65m and an average microphone spacing of 0.1m, geometrically arranged as in fig. 1. The array is parallel to the sound source plane and at a distance of 1m, the sound source plane is divided into grids with a pitch of 0.05m, and the number of grid points is 21 × 21 ═ 441. Gaussian pulse sound sources are used as sound sources to be identified, and the sound source identification performance of the four methods is contrastively analyzed. Suppose that there are three Gaussian pulse sound sources S1, S2 and S3 at (-0.04, -0.2) m, (0.08, -0.2) m and (-0.2,0.32) m, respectively, with half-peak bandwidths of 0.0001S and a peak time interval of 0.001S between the sources. The sound pressure signals at each microphone of the array are subjected to Fast Fourier Transform (FFT) to obtain Fourier spectra of the signals, wherein the sampling frequency is 102400Hz and the frequency resolution is 50 Hz. FIG. 2 shows the sound source identification results obtained by processing 81 frequencies from 1000Hz to 5000Hz with the above method when the SNR is 20 dB. Fig. 3 shows the recognition result at a corresponding signal-to-noise ratio of 5 dB. In the two groups of results, subgraphs (a) - (d) are cloud images for each method, and subgraphs (e) - (h) are results of each method for the respective frequency source of all sound sources.
For sound source S3, fig. 2(a) -2(d) show: the OMP-CB method has relatively dispersed positioning at each frequency and the worst positioning precision; the NOMP-CB method has relatively centralized frequency positioning and higher positioning precision; the BOMP-CB method sound source is positioned to nearby grid points, namely the positioning precision of the grid (base) is limited due to the existence of mismatch of the constraint bases; and the BNOMP-CB obtains accurate positioning. The source intensity curves in fig. 2(e) -2(h) show: the 4 methods can accurately quantize the intensity of each frequency source of the sound source S3. For the two sound sources S1 and S2 that are closer together, the OMP-CB method only achieves correct quantization at frequencies above 4500Hz, which is related to the fact that the OMP-CB method correctly separates the sound sources S1 and S2 at frequencies above 4500 Hz; the NOMP-CB method can accurately separate and correctly quantize the two sound sources when the frequency is higher than 2500 Hz; the BOMP-CB method does not obtain accurate quantification of the source intensities of the two sound sources within the frequency band of interest of 1000-; the BNOMP-CB method realizes the accurate separation and positioning of the two sound sources and obtains the accurate source strength estimation of the two sound sources in the whole frequency band of the concerned frequency band of 1000-plus 5000 Hz. In summary, the BNOMP-CB method has a spatial resolution far superior to that of the other 3 methods, has good applicability to transient sound source identification, and can accurately locate a sound source and obtain accurate quantification of the intensity of each frequency source of the sound source.
Comparing the recognition results of the methods in the case of the signal-to-noise ratio of 20dB in FIG. 2 and the recognition results of the methods in the case of the signal-to-noise ratio of 5dB in FIG. 3, it can be seen that the positioning results of the sound source at each frequency by the OMP-CB and NOMP-CB methods become more dispersed as the signal-to-noise ratio decreases; the positioning results obtained by the BOMP-CB and BNOMP-CB methods are unchanged and slightly changed. In addition, due to the introduction of strong noise, the source intensity quantization accuracy of all sound sources is reduced by the 4 methods. As can be seen from the positioning and quantization results of the sound sources S1 and S2, the NOMP-CB method still fails to accurately separate the two sound sources when the frequency is higher than 2500Hz when the signal-to-noise ratio is 5dB, but obtains more accurate strong source quantization when the frequency is close to 4500Hz, and the quantization curve presents a rule similar to that of the OMP-CB method. This phenomenon indicates that strong noise interference reduces the spatial resolution of the NOMP-CB method. While the quantization cases of the BOMP-CB and BNOMP-CB methods are similar to each other at a signal-to-noise ratio of 20dB except for the increase of fluctuation, the introduction of the multi-frequency synchronization processing makes the sound source identification performance of the two methods relatively robust to noise interference, and therefore, the proposed BNOMP-CB method also enjoys good sound source identification performance at a low signal-to-noise ratio.
2) Verification of comparative experiment
Three small speaker sources emitting random noise (Abramtek Type M5, Shenzhen, China) were tested for identification outdoors. The array and sound source arrangement is shown in FIG. 4(a), where the array is a Bruel with a 0.65m diameter microphone average pitch of 0.1m&
Figure BDA0003085921650000201
(Type 8608,
Figure BDA0003085921650000202
Denmark)36 channel planar microphone array with three loudspeakers arranged in a plane parallel to and 1m from the array plane. During the test, a PULSE data acquisition and analysis system (Bruel) is adopted&
Figure BDA0003085921650000203
Type 3660C,
Figure BDA0003085921650000204
Denmark) synchronously acquires sound pressure signals at 36 microphones, and the sampling frequency is 16384 Hz. And performing Fast Fourier Transform (FFT) on the sound pressure signal acquired by each microphone by using BKConnect to obtain a Fourier spectrum of the signal, wherein the frequency resolution is 50 Hz. Fig. 4(b) shows the amplitude spectrum of each microphone measurement signal.
The method comprises the following steps of adopting an orthogonal matching pursuit compression beam forming method (marked as OMP-CB) of single frequency processing, a Newton orthogonal matching pursuit compression beam forming method (marked as NOMP-CB) of single frequency processing, an orthogonal matching pursuit compression beam forming method (marked as BOMP-CB) of multi-frequency synchronous processing and the method invention (marked as BNOMP-CB) to carry out positioning identification on a small loudspeaker sound source, and comparing identification results:
fig. 5 is a sound source localization cloud image at [1000,1050., 5000] Hz (L ═ 81) broadband. As can be seen from fig. 5(a) (b), the respective frequency localization results of the OMP-CB method are scattered, while the localization results obtained by the NOMP-CB method are slightly more concentrated than those of the OMP-CB method, but both methods have low localization accuracy at many frequencies and estimate more false sound sources far from the sound source position. The BOMP-CB and BNOMP-CB methods are viewed reversely, and both methods can correctly position each sound source. In particular, for the left two sources, which are relatively close together, the OMP-CB and NOMP-CB methods fail to separate them correctly at many frequencies, whereas both the BOMP-CB and BNOMP-CB methods succeed in separating the two sources and obtaining accurate position estimates. In addition, the positioning results of the BOMP-CB and the BNOMP-CB on the left lower sound source and the right upper sound source are compared, the BNOMP-CB effectively overcomes the problem of base mismatching due to the introduction of Newton optimization, and positioning superior to that of the BOMP-CB is obtained. In conclusion, test cases show that the proposed BNOMP method has good sound source identification performance due to the overcoming of Newton optimization to the base mismatching problem and the promotion of multi-frequency joint sparse constraint to sparsity.

Claims (7)

1. The method for identifying the sound source formed by the multi-frequency synchronous two-dimensional off-network compressed wave beams is characterized by comprising the following steps of:
1) and establishing the multi-frequency sound pressure signal measurement model.
2) Establishing a maximum likelihood estimation model with the sound source position and the source intensity as parameters;
3) performing preliminary estimation on the sound source position to obtain a preliminary estimation coordinate of the sound source position;
4) carrying out local optimization on the preliminary estimation coordinates of the sound source position by utilizing a Newton method to obtain optimized estimation coordinates of the sound source position;
5) carrying out global cyclic optimization on all the recognized sound source position coordinates by utilizing a Newton method to obtain global optimization estimation coordinates of all the recognized sound source positions;
6) and carrying out orthogonal solving on all the recognized sound source intensities, and updating the residual error according to the obtained sound source intensities.
2. The method of claim 1, wherein the multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification model is as follows:
p＝A(ΩS)qS+n (1)
in the formula (I), the compound is shown in the specification,
Figure FDA0003085921640000011
source intensity vectors representing source intensities of the S sound sources at each frequency;
Figure FDA0003085921640000012
source intensity vectors of all frequencies of the s-th sound source;
Figure FDA0003085921640000013
representing a microphone measurement noise vector; n is·,l＝[n1,l,n2,l,...,nM,l]Representing a measured noise vector corresponding to the microphone at the l-th frequency point;
Figure FDA0003085921640000014
a transfer matrix representing all L frequencies between S sound sources to all microphones;
Figure FDA0003085921640000015
a coordinate matrix composed of coordinates of all sound source positions,
Figure FDA0003085921640000016
representing a complex set;
Figure FDA0003085921640000017
representing a set of real numbers; p is the actually measured microphone sound pressure vector;
wherein the mth microphone has the coordinate of
Figure FDA0003085921640000018
The coordinates of the s-th sound source are
Figure FDA0003085921640000019
Wherein the content of the first and second substances,
Figure FDA00030859216400000110
the x-axis coordinate and the y-axis coordinate of the mth microphone;
Figure FDA00030859216400000111
the x-axis coordinate and the y-axis coordinate of the s-th sound source are obtained; h is the distance between the sound source plane and the array plane; m is 1,2, 3 …, M; s is 1,2, 3 …, S; m, S are the number of microphones and sound sources, respectively;
distance d between mth microphone and s sound sourcem,sAs follows:
Figure FDA00030859216400000112
in the formula, vector
Figure FDA00030859216400000113
Vector
Figure FDA00030859216400000114
Superscript T denotes transpose; i | · | purple wind2A two-norm representation of a vector;
Figure FDA0003085921640000021
a transfer matrix representing all L frequencies between the s-th source to all microphones is shown below:
Figure FDA0003085921640000022
in the formula (I), the compound is shown in the specification,
Figure FDA0003085921640000023
is a frequency flTransmission of the s-th sound source to all microphonesVector increment quantity
Figure FDA0003085921640000024
3. The method of claim 2, wherein the step of establishing a maximum likelihood estimation model with the sound source position and the source intensity as parameters comprises:
1) minimizing residual errors of a multi-frequency acoustic pressure signal measurement model
Figure FDA0003085921640000025
Obtaining maximum likelihood estimates of sound source position and source intensity
Figure FDA0003085921640000026
Namely:
Figure FDA0003085921640000027
in the formula (I), the compound is shown in the specification,
Figure FDA0003085921640000028
representing the real part of the complex number; superscript H represents the hermitian conjugate transpose operator;
2) establishing a likelihood function
Figure FDA0003085921640000029
Namely:
Figure FDA00030859216400000210
3) source intensity vector q formed by source intensities of S sound sources at various frequenciesSFor variables, a least squares solution is established, i.e.:
qS＝(A(ΩS)HA(ΩS))-1A(ΩS)Hp (6)
4) establishing a maximum likelihood ratio test cost function
Figure FDA00030859216400000211
Namely:
Figure FDA00030859216400000212
5) establishing a sound source position coordinate estimation equation, namely:
Figure FDA00030859216400000213
in the formula (I), the compound is shown in the specification,
Figure FDA00030859216400000214
coordinates are estimated for the sound source position.
4. The method of claim 3, wherein the step of performing the preliminary estimation of the sound source location comprises:
1) dispersing a sound source detection space into N grid points; n > M; wherein the coordinates of the nth grid point are
Figure FDA00030859216400000215
Set of grid point coordinates
Figure FDA00030859216400000216
Wherein, the vector
Figure FDA00030859216400000217
2) Establishing a sensing matrix of all frequency correspondences between N grid points and M microphones
Figure FDA0003085921640000031
Wherein, the sensing matrix corresponding to all frequencies from the nth grid point to the M microphones
Figure FDA0003085921640000032
3) The actual measured microphone sound pressure vector p is represented as follows:
p＝A(ΩG)qG+n (9)
in the formula (I), the compound is shown in the specification,
Figure FDA0003085921640000033
a source intensity vector formed by source intensities of all frequencies of N grid point positions is represented;
Figure FDA0003085921640000034
source strength vector for all frequencies of nth grid point
Figure FDA0003085921640000035
Representing the frequency f at the nth grid pointlThe source of time is strong;
4) updating likelihood function, maximum likelihood ratio test cost function and least square function as identification formula for single sound source to obtain function
Figure FDA0003085921640000036
Sum function
Figure FDA0003085921640000037
Namely:
Figure FDA0003085921640000038
Figure FDA0003085921640000039
Figure FDA00030859216400000310
5) from the set omegaGIn selecting a maximum likelihood ratio test cost function
Figure FDA00030859216400000311
The largest grid node is taken as the initial estimation coordinate of the sound source position and is recorded as
Figure FDA00030859216400000312
Preliminary estimation coordinates of sound source position
Figure FDA00030859216400000313
As follows:
Figure FDA00030859216400000314
in the formula (I), the compound is shown in the specification,
Figure FDA00030859216400000315
a residual vector representing the microphone sound pressure signal obtained after the first s-1 sound source estimation;
Figure FDA00030859216400000316
and
Figure FDA00030859216400000317
respectively the estimated coordinates and source strength of the jth sound source.
6) Preliminarily estimating coordinates of sound source position
Figure FDA00030859216400000318
Substituting the following formula to obtain the source intensity estimation vector of all frequencies of the s sound source, namely:
Figure FDA00030859216400000319
in the formula (I), the compound is shown in the specification,
Figure FDA00030859216400000320
vectors are estimated for the source strengths of all frequencies of the s-th source.
5. The method of claim 1, wherein the step of locally optimizing the preliminary estimated coordinates of the sound source location by newton's method comprises:
1) establishing a two-dimensional Newton optimization function, namely:
Figure FDA0003085921640000041
in the formula (I), the compound is shown in the specification,
Figure FDA0003085921640000042
and
Figure FDA0003085921640000043
respectively representing a Jacobian matrix and a Hessian matrix;
Figure FDA0003085921640000044
optimizing and estimating coordinates for the sound source position;
jacobian matrix
Figure FDA0003085921640000045
Heisen matrix
Figure FDA0003085921640000046
Respectively as follows:
Figure FDA0003085921640000047
Figure FDA0003085921640000048
wherein the parameter item
Figure FDA0003085921640000049
Parameter item
Figure FDA00030859216400000410
Parameter item
Figure FDA00030859216400000411
Parameter item
Figure FDA00030859216400000412
Parameter item
Figure FDA00030859216400000413
Parameter item
Figure FDA00030859216400000414
Respectively as follows:
Figure FDA00030859216400000415
Figure FDA00030859216400000416
Figure FDA00030859216400000417
Figure FDA00030859216400000418
Figure FDA00030859216400000419
Figure FDA00030859216400000420
in the formula, AfRepresenting a transfer matrix at a frequency f;
2) optimizing and estimating coordinates of sound source position
Figure FDA0003085921640000051
Substituting the obtained value into the following formula to obtain the sound source intensity optimized estimation
Figure FDA0003085921640000052
Namely:
Figure FDA0003085921640000053
3) updating microphone sound pressure residual vector prNamely:
Figure FDA0003085921640000054
in the formula (I), the compound is shown in the specification,
Figure FDA0003085921640000055
the transfer matrix at all frequencies between the jth source to all M microphones.
6. The method of claim 1, wherein the step of performing a global loop optimization on all the identified sound source position coordinates by using newton's method comprises:
1) establishing a source strength estimation vector
Figure FDA0003085921640000056
Transmission matrix
Figure FDA0003085921640000057
And harmonic source location estimation matrices
Figure FDA0003085921640000058
2) And (3) carrying out the t iterative optimization on the kth sound source, wherein the steps are as follows:
2.1) computing the transition residual vector of the kth sound source component at the time of retaining the t-1 th iteration
Figure FDA0003085921640000059
Namely:
Figure FDA00030859216400000510
in the formula, the initial value of t is 1;
2.2) residual transition
Figure FDA00030859216400000511
Substituting into the following formula to obtain the transition source intensity
Figure FDA00030859216400000512
Namely:
Figure FDA00030859216400000513
2.3) calculating the global optimization estimated coordinates of the sound source position
Figure FDA00030859216400000514
Namely:
Figure FDA00030859216400000515
2.4) updating microphone pressure residual vector
Figure FDA00030859216400000516
Namely:
Figure FDA00030859216400000517
3) determining the residual energy change of two consecutive cycles
Figure FDA00030859216400000518
And if so, performing orthogonal solution on all the recognized sound source intensities and updating the residual error, otherwise, making t equal to t +1, and returning to the step 2).
7. The method of claim 1, wherein the method for recognizing a sound source by multi-frequency synchronous two-dimensional off-network compressed beam forming comprises: the steps of solving all the identified sound source intensities orthogonally and updating the residual errors include:
1) optimizing and estimating coordinate set of sound source position
Figure FDA0003085921640000061
Substituting into the following formula to obtain the sound source intensity estimation
Figure FDA0003085921640000062
Namely:
Figure FDA0003085921640000063
2) updating microphone sound pressure residual vector prNamely:
Figure FDA0003085921640000064
3) judgmentIf the K is greater than or equal to S, if so, ending the circulation and outputting the set
Figure FDA0003085921640000065
And
Figure FDA0003085921640000066
otherwise, let k be k +1 and from the set ΩGIn selecting a maximum likelihood ratio test cost function
Figure FDA0003085921640000067
And the largest grid node is used as a preliminary estimation coordinate of the sound source position.
CN202110580351.4A 2021-05-26 2021-05-26 Multi-frequency synchronous two-dimensional off-grid compressed beam forming sound source identification method Active CN113687306B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110580351.4A CN113687306B (en) 2021-05-26 2021-05-26 Multi-frequency synchronous two-dimensional off-grid compressed beam forming sound source identification method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110580351.4A CN113687306B (en) 2021-05-26 2021-05-26 Multi-frequency synchronous two-dimensional off-grid compressed beam forming sound source identification method

Publications (2)

Publication Number Publication Date
CN113687306A true CN113687306A (en) 2021-11-23
CN113687306B CN113687306B (en) 2023-07-21

Family

ID=78576500

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110580351.4A Active CN113687306B (en) 2021-05-26 2021-05-26 Multi-frequency synchronous two-dimensional off-grid compressed beam forming sound source identification method

Country Status (1)

Country Link
CN (1) CN113687306B (en)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109061727A (en) * 2018-08-14 2018-12-21 中国石油大学(华东) A kind of viscous acoustic medium full waveform inversion method of frequency domain
CN109870669A (en) * 2019-02-19 2019-06-11 重庆工业职业技术学院 How soon a kind of two dimension claps mesh free compression Wave beam forming identification of sound source method
CN111830465A (en) * 2020-07-27 2020-10-27 重庆大学 Two-dimensional Newton orthogonal matching tracking compressed beam forming method
CN111965599A (en) * 2020-07-03 2020-11-20 重庆大学 Sound source identification method for two-dimensional dynamic grid compressed beam forming
CN112198517A (en) * 2020-09-03 2021-01-08 浙江工业大学 High-precision ultrasonic distance measurement method adaptive to high-speed movement of sound source and high environmental wind

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109061727A (en) * 2018-08-14 2018-12-21 中国石油大学(华东) A kind of viscous acoustic medium full waveform inversion method of frequency domain
CN109870669A (en) * 2019-02-19 2019-06-11 重庆工业职业技术学院 How soon a kind of two dimension claps mesh free compression Wave beam forming identification of sound source method
CN111965599A (en) * 2020-07-03 2020-11-20 重庆大学 Sound source identification method for two-dimensional dynamic grid compressed beam forming
CN111830465A (en) * 2020-07-27 2020-10-27 重庆大学 Two-dimensional Newton orthogonal matching tracking compressed beam forming method
CN112198517A (en) * 2020-09-03 2021-01-08 浙江工业大学 High-precision ultrasonic distance measurement method adaptive to high-speed movement of sound source and high environmental wind

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
XIN ZHANG: "An Alternative Hybrid Time-Frequency Domain Approach Based on Fast Iterative ShrinkageThresholding Algorithm for Rotating Acoustic Source Identification", 《 IEEE ACCESS》, vol. 7, pages 59797 - 59805, XP011725195, DOI: 10.1109/ACCESS.2019.2915001 *

Also Published As

Publication number Publication date
CN113687306B (en) 2023-07-21

Similar Documents

Publication Publication Date Title
US10429488B1 (en) System and method for geo-locating and detecting source of electromagnetic emissions
CN112526451B (en) Compressed beam forming and system based on microphone array imaging
CN107247251B (en) Three-dimensional sound source positioning method based on compressed sensing
CN106021637B (en) DOA estimation method based on the sparse reconstruct of iteration in relatively prime array
Chu et al. Deconvolution for three-dimensional acoustic source identification based on spherical harmonics beamforming
CN111123192A (en) Two-dimensional DOA positioning method based on circular array and virtual extension
CN113064147B (en) Novel matching field passive positioning method under low signal-to-noise ratio
CN115236594B (en) Under-ice sound source positioning method suitable for polar impulse noise environment
US20120109563A1 (en) Method and apparatus for quantifying a best match between series of time uncertain measurements
CN108398659B (en) Direction-of-arrival estimation method combining matrix beam and root finding MUSIC
CN111812581B (en) Spherical array sound source direction-of-arrival estimation method based on atomic norms
CN116068493A (en) Passive sound source positioning method for deep sea large-depth vertical distributed hydrophone
CN115032591A (en) Broadband multi-sound-source positioning asynchronous measurement method and device and related medium
CN113866718A (en) Matching field passive positioning method based on co-prime matrix
CN111830465B (en) Two-dimensional Newton orthogonal matching pursuit compressed beam forming method
Ping et al. Compressive spherical beamforming for acoustic source identification
CN113687306A (en) Multi-frequency synchronous two-dimensional off-network compressed beam forming sound source identification method
Yin et al. Super-resolution compressive spherical beamforming based on off-grid sparse Bayesian inference
CN113381793B (en) Coherent information source estimation-oriented non-grid direction-of-arrival estimation method
Oertwig et al. Extension of the source localization method SODIX for the determination of partially coherent sound sources
CN105652264B (en) Multipath propagation acoustic signal based on Higher Order Cumulants
Pan Spherical harmonic atomic norm and its application to DOA estimation
Rani et al. High-resolution spectral analysis methods for self-potential data inversion
CN115267673A (en) Sparse sound source imaging method and system considering reconstruction grid offset
Mao et al. Scanning Measurement Based on FISTA Phase Calibration

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