Method for direction of arrival estimation based on robust sparse bayesian learning with variational inference

Abstract

This invention provides a novel robust sparse Bayesian learning-based direction-of-arrival (DOA) estimation method. Traditional sparse Bayesian learning methods based on the Gaussian noise assumption degrade in impulsive noise environments. This invention transforms the ODA problem in array signal processing into a sparse reconstruction problem. It inputs a perception matrix composed of the array observation vector and the array steering vector, initializes parameters, iteratively updates the parameters and the posterior distribution of the sparse signal, and finally uses the approximate posterior expectation of the sparse signal as the reconstruction result. The ODA result is obtained based on the support set of the reconstructed sparse signal. A robust sparse Bayesian learning algorithm based on variational inference is designed to obtain higher accuracy reconstruction and ODA results in non-Gaussian noise environments by inferring the approximate posterior distribution. This invention exhibits more robust estimation performance in impulsive noise environments.

Description

Technical Field

[0001] This invention relates to array signal processing technology, particularly to direction-of-arrival estimation technology for communication, radar and other systems in impact noise environments. Technical Background

[0002] Direction-of-arrival (DOA) estimation, as an important research topic in array signal processing, has received widespread attention and in-depth study in radar and communication fields. To avoid the major shortcomings of traditional subspace methods in terms of sampling quantity and coherent sources, research on DOA estimation methods based on sparse reconstruction has been gradually carried out. Existing sparse reconstruction algorithms can handle Gaussian noise well, but impulse noise, due to its significant energy characteristics over a short period, can have a global impact on observations within the traditional compressed sensing framework, leading to a deterioration in the reconstruction performance of traditional sparse reconstruction algorithms. Robust algorithms based on the traditional sparse Bayesian learning framework assume that all observations are affected by impulse noise; however, in real-world environments, impulse noise often has a short duration, affecting only a portion of the observations. The robust sparse Bayesian learning algorithm based on variational inference sets a flag variable to indicate whether the observation is affected by impact noise, and then establishes different conditional distribution models for different types of measurement data. Among them, BP-RBCS (Beta-Bernoulli Prior model-based Robust Bayesian Compressed Sensing) removes the identified outliers from the observation vector. The model is relatively simple but loses some observation data. Q.Wan,HPDuan,J.Fang,HBLi,ZLXing.Robust Bayesian compressed sensing with outliers[J].Signal Processing,2017,5(17):104-109.; Mix-RSBL algorithm (Robust Sparse Bayesian Learning using Mixture The model makes a prior probability distribution assumption on the impact component in the noise, identifies and processes the observation data affected by the impact noise, and applies it to reconstruction and estimation. It makes fuller use of the observation data, improves the accuracy of estimation, and avoids waste of resources. However, the number of parameters that need to be updated increases significantly, so the computational load is greater. R.Zheng,X.Xu,ZFYe,etal.Robust sparse Bayesian learning for DOA estimation in impulsive noise environments[J].Signal Processing,2020,171(5):1-6. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to provide a method for obtaining more accurate sparse reconstruction results under the influence of impact noise without significantly increasing the amount of computation, thereby completing the direction of arrival estimation.

[0004] The technical solution adopted by this invention to solve the above-mentioned technical problems is a direction-of-arrival estimation method based on variational inference robust sparse Bayesian learning, which includes the following steps:

[0005] 1. A direction-of-arrival estimation method based on robust sparse Bayesian learning with variational inference, characterized by the following steps:

[0006] Input steps: Input the sensing matrix A composed of the array observation vector y and the array steering vector, determine the total number of directions of arrival N, the total number of arrays M, and the set of directions of arrival.

[0007] Initialization steps: iteration number L, iteration threshold ε, hyperparameters α, β1, β2, z, a, b, c1, d1, c2, d2, e, f; where α is the inverse variance of the support set for the direction-of-arrival estimation, and each element of α follows a gamma distribution; z is a flag vector indicating whether the observations are affected by impulse noise, and each element of z follows a Bernoulli distribution; β1 is the inverse variance of the conditional likelihood function of the array observation vectors unaffected by impulse noise; β2 is the inverse variance of the likelihood function of the array observation data affected by noise; a, b, c1, d1, c2, d2, e, f are all parameters used when iteratively updating the posterior distribution of the sparse signal;

[0008] Iterative steps: The specific steps for iteratively updating the hyperparameters and the posterior distribution of the sparse signal are as follows:

[0009] (1) Calculate the covariance matrix Φ of the sparse signal using the sensing matrix A, the currently updated flag vector z, and the inverse variance α of the support set of the direction-of-arrival estimation. x =(A T DA+D α ) -1 ; T This represents the matrix transpose, where matrix D = <β1>D z +<β2>D 1-z D z =diag( <z>), D 1-z =ID z , diag means extracting the diagonal matrix, < > means the estimated value, D α =diag(<α>);

[0010] (2) Utilizing the covariance matrix Φ of the sparse signal x The reconstructed vector update value μ is calculated using the perception matrix A and the array observation vector y. x =Φ x A T Dy, then the latest direction-of-arrival estimation reconstruction vector x new Updated to μ x The reconstruction vector of the previously updated direction-of-arrival estimate is x. old ;

[0011] (3) Update the estimated value corresponding to each element in α. Direction of arrival (DOA) indices n = 1, ..., N; where the intermediate values ​​are... x n Let x be the nth element of the support set x for direction-of-arrival estimation, and x = x new , Φ x (n, n) is the covariance matrix Φ of the sparse signal. x The element in the nth row and nth column;

[0012] (4) Update the estimate of β1 Among them, the median value <(y-Ax) T D z (y-Ax)>=(y-Aμ x ) T D z (y-Aμ x )+tr(A T D z AΦ x ), tr represents the trace of the matrix, m is the array index, and z m This is the m-th element of the flag vector z;

[0013] (5) Update the estimate of β2 Among them, the median value

[0014] (6) Update the estimated value corresponding to each element in z. Where p is the probability density function, p(z) m =1)=Cexp[-0.5<β1><(y m -a m x) 2 >+(lnπ m >],π m Let a be the m-th element of probability π that follows different noise distributions. m p(z) represents the m-th row of the array guiding vector matrix A; m =0)=Cexp[-0.5<β2><(y m -a m x) 2 >+(ln(1-π m )>], where C is a constant,

[0015] (7) Calculate ||x new -x old || 2 If ||x new -x old || 2 If ||x|| ≤ ε or k = L, then proceed to the following output step; if ||x|| ≤ ε or k = L, then proceed to the next output step; new -x old || 2 If x > ε and k < L, then x old The value is updated to x new Return to step 1;

[0016] Output steps: for x new A binarized vector is obtained through thresholding, thus yielding the estimated support set I for the direction of arrival (DOA). The estimated DOA vector is constructed using the angles corresponding to each element in I within the DOA set p. And output it.

[0017] This invention models noise as a Bernoulli-Gaussian distribution, forming a Bayesian probability model controlled by multiple parameters. Likelihood functions are established for array observation data affected and unaffected by impact noise using two variance parameters. A novel variational inference sparse Bayesian learning algorithm is designed. This algorithm iteratively updates parameters and the approximate posterior distribution of the sparse signal, ultimately using the approximate posterior expectation of the sparse signal as its estimate. This achieves higher accuracy in sparse reconstruction results without a significant increase in computational cost, thereby obtaining a support set and ultimately realizing direction-of-arrival estimation. Attached Figure Description

[0018] Figure 1 This is a comparison chart of the normalized mean square error of the BG-RSBL algorithm of this invention with the existing BP-RBCS and Mix-RSBL algorithms as a function of signal-to-noise ratio.

[0019] Figure 2 This is a comparison chart of the estimation success rate of the present invention BG-RSBL with that of existing BP-RBCS and Mix-RSBL as a function of signal-to-noise ratio. Detailed Implementation

[0020] The present invention proposes a direction-of-arrival estimation method based on variational inference robust sparse Bayesian learning, hereinafter referred to as BG-RSBL. BG-RSBL: Robust Sparse Bayesian Learning based on the Bernoulli-Gaussian noise model.

[0021] The flowchart of the direction-of-arrival estimation method based on BG-RSBL is shown in the table below:

[0022]

[0023]

[0024]

[0025] principle:

[0026] 1. Impulse noise and non-impulse noise are expressed using two Gaussian distribution models with different parameters, using the indicator variable z. i For the observed value y i To distinguish whether it is affected by impact noise, a conditional likelihood function for the array observation vector is established.

[0027] Assume the array observation vector is y = [y1, y2, ..., y3]. M ] T for:

[0028]

[0029] in,

[0030] (1)a i. Let the i-th row be the array guide vector matrix A;

[0031] (2) x is the support set for the direction-of-arrival estimation. Assume that each element follows a Gaussian distribution with zero mean and unknown variance, and that they are independent of each other. Its probability density function is: The reciprocal of the variance of x, α = [α1, α2, ..., α...]. N ] T Assuming that all elements follow a gamma distribution, i.e.

[0032] (3) For Gaussian components with small variance, the reciprocal of the variance For Gaussian components with large variance, the reciprocal of the variance

[0033] (4) z = [z1, z2, ..., z M ] T Let the vector be a label, and its elements follow a Bernoulli distribution: When z i When = 1, it means that the noise component in the i-th row of observation data is n. 1i When z i When = 0, it means that the noise component in the i-th row of observation data is n. 2i The elements of π follow a beta distribution, and their joint probability density is: π ~ [Beta(e, f)] M .

[0034] Based on the above assumptions, the conditional likelihood function of the observation vector can be expressed as follows:

[0035]

[0036] 2. Under the variational inference framework, the approximate posterior probability distributions of hyperparameters and sparse signals are obtained by minimizing the KL divergence, and the relevant parameters of each distribution are iteratively updated.

[0037] The approximate posterior probability density can be obtained through lnq. j =<lnp(y,θ)> i≠j +const is obtained, where θ={θ1, θ2,…,θ n } represents all unknown parameters, q j The <> represents the approximate posterior probability density of each component, and <> represents the estimated value, i.e., the expectation. A This represents the expectation of the probability with respect to index A. Decomposing the joint probability density p(y, θ) and ignoring the parts unrelated to the parameters, we get:

[0038] (1)q x (x)

[0039]

[0040] Substituting the conditional likelihood function of y and the prior probability density of x, we get:

[0041]

[0042] Among them, D z =diag( <z>), D 1-z =ID z D α =diag(<α>), so the approximate posterior distribution of x is q. x (x)~N(μ x , Φ x ),in:

[0043] Φ x =(A T DA+D α ) -1

[0044] μ x =Φ x A T Dy

[0045] D=<β1>D z +<β2>D 1-z

[0046] (2)q α (α)

[0047]

[0048] Substituting the prior probability density of x and the prior probability density of α, we get:

[0049]

[0050] Therefore, the approximate posterior distribution of α is: in,

[0051]

[0052] (3)

[0053]

[0054] therefore, It is a gamma distribution, that is The parameters are as follows:

[0055]

[0056] (4)

[0057]

[0058] therefore, It is a gamma distribution, that is The parameters are as follows:

[0059]

[0060] (5)q z (z)

[0061]

[0062] Substituting the conditional likelihood function of y and the prior probability density of z, we get:

[0063]

[0064] Therefore, q π (π) follows a Bernoulli distribution, with probabilities of taking the values ​​1 and 0 respectively:

[0065] p(z m =1)=Cexp[-0.5<β1><(y m -a m x) 2 >+(lnπ m >]

[0066] p(z m =0)=Cexp[-0.5<β2><(y m -a m x) 2 >+ <ln(1-π m )>]

[0067] (6)q π (π)

[0068]

[0069] Substituting the prior probability density of z, we get:

[0070]

[0071] Therefore, q π (π) represents the beta distribution, i.e. in:

[0072]

[0073] Based on the approximate posterior probability derived above, the hyperparameter α... n , β1, β2, z m Iterative updates are performed based on its approximate posterior expectation:

[0074]

[0075] 3. Put μ x The results of two consecutive iterations are represented as x. old With x new If ||x new -x old || 2 If the iteration count reaches the upper limit (<ε), then x will be... new Use this as an estimate of the direction-of-arrival support set x; otherwise, continue iterating.

[0076] 4. Regarding x new After thresholding, the support set I is obtained, and the actual direction of arrival is estimated as follows: p is the direction vector corresponding to the guide vector. i The direction corresponding to the i-th element in support set I.

[0077] Performance simulation verification:

[0078] The simulation experiment assumes a uniform linear array with M = 90 elements, k = 5 sources, and a spatial search range of 0–179° with a search interval of 1° (N = 180°). The source variance α is initialized as a vector of all 1s, following a gamma distribution with constant coefficients a = 10. -10 b = 10 -10 The initial noise variance β1 = β2 = 1 follows a gamma distribution with constant coefficients c1 = 100, d1 = 1, and c2 = d2 = 10. -10 The algorithm terminates when the number of iterations reaches 250, or when the l2 norm of the difference between two iterations is less than 10. -6 L = 300 Ment Carlo experiments were conducted, using Bernoulli-Gaussian noise to simulate the impact noise environment. The probability of noise component generation was p1 = 0.8, p2 = 0.2, and the variance ratio was... The performance of the algorithm is analyzed using the normalized mean square error (NMSE) and the estimation success rate. The formula for calculating NMSE is: The method for estimating the success rate is as follows: calculate the success rate after each experiment. If the value is less than 1°, the estimation is considered successful. The percentage of successful experiments out of the total number of experiments is the success rate of direction estimation.

[0079] The effects of this invention are demonstrated through simulation. Figure 1 and 2 Please provide an explanation. Figure 1 This is a comparison chart of the normalized mean square error of the present invention with the BP-RBCS algorithm and the Mix-RSBL algorithm as a function of signal-to-noise ratio. Figure 2 This is a comparison graph showing the estimation success rate of the present invention, the BP-RBCS algorithm, and the Mix-RSBL algorithm as a function of signal-to-noise ratio. The simulation results show that the present invention has higher estimation accuracy and more robust estimation performance in impulsive noise environments compared to existing algorithms.< / z> < / z>

Claims

1. A direction-of-arrival estimation method based on robust sparse Bayesian learning with variational inference, characterized in that, Includes the following steps: Input steps: Input the sensing matrix A composed of the array observation vector y and the array steering vector, determine the total number of directions of arrival N, the total number of arrays M, and the set of directions of arrival. Initialization steps: iteration number L, iteration threshold ε, hyperparameters α, β1, β2, z, a, b, c1, d1, c2, d2, e, f; where α is the inverse variance of the support set for the direction-of-arrival estimation, and each element of α follows a gamma distribution; z is a flag vector indicating whether the observations are affected by impulse noise, and each element of z follows a Bernoulli distribution; β1 is the accuracy of the conditional likelihood function of the array observation vectors unaffected by impulse noise; β2 is the accuracy of the likelihood function of the array observation data affected by noise; a, b, c1, d1, c2, d2, e, f are all parameters used when iteratively updating the posterior distribution of the sparse signal. Iterative steps: The specific steps for iteratively updating the hyperparameters and the posterior distribution of the sparse signal are as follows: (1) Calculate the sparse signal Φ using the sensing matrix A, the currently updated flag vector z, and the inverse variance α of the support set of the direction-of-arrival estimation. x =(A T DA+D α ) -1 T denotes matrix transpose, and matrix D = <β1>D z +<β2>D 1-z D z =diag( <z>), D 1-z =ID z ,diag indicates extracting the diagonal matrix, 〈〉 indicates the estimated value, D α =diag(〈α〉);< / z> (2) Utilizing sparse signal Φ x The reconstructed vector update value μ is calculated using the perception matrix A and the array observation vector y. x =Φ x A T D y Then, the reconstructed vector x of the latest direction-of-arrival estimation is... new Updated to μ x The reconstruction vector of the previously updated direction-of-arrival estimate is x. old ; (3) Update the estimated value corresponding to each element in α. The direction of arrival (DOA) numbers are n = 1, ..., N; where the intermediate values ​​are... x n Let x be the nth element of the support set x for direction-of-arrival estimation, and x = x new ,Φ( n ,n) is a sparse signal Φ x The element in the nth row and nth column; (4) Update the estimate of β1 Among them, the median value <(y-Ax) T D z (y-Ax)>=(y-Aμ) x ) T D z (y-Aμ x )+tr(A T D z AΦ x ), tr represents the trace of the matrix, m is the array index, and z m This is the m-th element of the flag vector z; (5) Update the estimate of β2 Among them, the median value <(y-Ax) T D 1-z (y-Az)>=(y-Aμ) x ) T D z (y-Aμ x )+tr(A T D z AΦ x ), (6) Update the estimated value corresponding to each element in z. Where p is the probability density function, p(z) m =1)=Cexp[-0.5<β1><(y m -a m x) 2 >+ <lnπ m >],π m Let a be the m-th element of probability π that follows different noise distributions. m p(z) represents the m-th row of the array guiding vector matrix A; m =0)=Cexp[-0.5<β2><(y m -a m x) 2 >+ <ln(1-π m )>], where C is a constant, y m Let m be the m-th element of the array observation vector y; (7) Calculate ||x new -x old || 2 If ||x new -x old || 2 If ||x|| ≤ ε or the current iteration number k = L, then proceed to the following output step; if ||x|| ≤ ε, then proceed to the next output step; new -x old || 2 If x > ε and k < L, then x old The value is updated to x new Return to step 1; Output steps: for x new A binarized vector is obtained through thresholding, thus yielding the estimated support set I for the direction of arrival (DOA). The estimated DOA vector is constructed using the angles corresponding to each element in I within the DOA set p. And output it.

Images