CN114325567A - A High Accuracy Estimation Method for Angle of Arrival - Google Patents

Abstract

The invention discloses a high-precision estimation method of a DOA (angle of arrival), which sequentially comprises the following steps of: a: firstly, establishing a lattice DOA bilinear system model, and then carrying out factorization on the lattice DOA bilinear system model by using a Bayes formula; b: for sparse incoming wave signalSCarrying out estimation; c: for lattice error parameter vectorβCarrying out estimation; d: according to the obtained sparse incoming wave signalSSum-off-grid error parameter vectorβAnd D, performing high-precision calculation on the arrival angle through the off-grid DOA bilinear system model obtained in the step A. The method can solve the contradiction between complexity and robustness in the off-grid DOA estimation, effectively improves the calculation precision of the DOA, and has obvious performance advantages in the scene with less snapshots such as airborne radar.

Description

A High Accuracy Estimation Method for Angle of Arrival

technical field

The invention relates to the field of array signal processing, in particular to a high-precision estimation method of the angle of arrival in array signal processing.

Background technique

Estimating the angle of arrival (DOA) of multiple far-field signals using an antenna array is a fundamental problem in array signal processing, and has been widely used in radar, sonar, and earthquake prediction. For the DOA estimation problem, domestic and foreign researchers have proposed a variety of different types of methods, such as subspace decomposition, support vector machine, matrix eigenspace decomposition, etc. Each method has its own advantages and disadvantages. With the emergence of compressed sensing technology, researchers have generalized DOA estimation as a sparse recovery problem, and proposed a variety of different estimation algorithms, such as orthogonal matching pursuit algorithm. In the usual sparse estimation method, the range of the angle of arrival is divided into a uniform grid, and it is assumed that the actual angle of arrival is exactly on a certain grid. Compared with methods such as subspace division and feature space decomposition, compressed sensing methods can usually achieve better performance, especially when the number of snapshots is small. Sparse Bayesian learning is a class of classical sparse recovery algorithms that can achieve better estimation performance with fewer observations by adding sparse bootstrap priors to variables. But when the true incoming wave direction is not precisely located on the grid, the above scheme will lead to serious estimation bias. On the one hand, finer grid division will lead to higher computational complexity, and on the other hand, it will lead to a higher correlation between the column vectors of the perception matrix, resulting in poor convergence of the estimation algorithm.

In recent years, many teams at home and abroad have proposed a variety of sparse reconstruction algorithms based on off-grid signals. Assuming that the real angle is not limited to the grid, the error can be calculated and the estimated value can be corrected. For example, the root sparse Bayesian learning method (RSBL) can achieve good estimation performance in linear array scenarios. Due to its low complexity, the classical approximate message passing algorithm (AMP) is widely used in the sparse Bayesian learning framework, but the ordinary AMP algorithm will cause problems in the iterative process when encountering a less Gaussian distribution of observation matrices diverge.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a high-precision estimation method for the angle of arrival, which can solve the contradiction between complexity and robustness in the estimation of the off-grid angle of arrival, effectively improve the calculation accuracy of the angle of arrival, and can be used in airborne radar and other fast There is a clear performance advantage in scenes with fewer shots.

The present invention adopts following technical scheme:

A high-precision estimation method for the angle of arrival, which sequentially includes the following steps:

A: First, establish an off-grid DOA bilinear system model Y=(A+∑ _n β _n E _n )S+W=Φ{β}S+W;

n＝1,…,N，s^t∈C^N×1表示长度为N的待求稀疏来波信号向量；L为快拍数量，w^t表示均值为0且方差为1/λ的高斯白噪声向量；矩阵Y＝[y¹,y²,…,y^L]、S＝[s¹,s²,…,s^L]和W＝[w¹,w²,…,w^L]分别表示接收信号向量、来波信号向量和高斯白噪声向量的集合矩阵，矩阵S具有相同的稀疏特征；A和E分别为阵列流型矩阵A(θ)和误差矩阵E(θ)的简记；θ为来波信号的角度空间；Among them, Φ{β}=A+Ediag(β)=A+∑ _n β _n E _n , E _n =[0,…,e(θ _n ),…,0]∈C ^M×N , off-grid DOA double The linear system model is a shorthand for the received signal vector Y=(A(θ)+E(θ)diag(β))S+W; Y=(A(θ)+E(θ)diag(β))S +W is the multiple observation form of y ^t =(A(θ)+E(θ)diag(β))s ^t +w ^t ; Suppose there is a uniform linear array composed of M non-directional array elements, there are K Narrowband far-field incoming signal, y ^t represents the received signal vector in the t-th snapshot, A(θ)=[a(θ ₁ ),...,a(θ _N )]∈C ^M×N , C ^M×K represents a complex-valued matrix with dimension M×K; suppose the angle space θ of the incoming wave signal is divided into N grids, expressed in the form of a vector θ=[θ ₁ ,...,θ _N ], each grid The size is Δ _θ = π/N; the matrix E(θ)=[e(θ ₁ ),…,e(θ _N )] represents the error matrix, and the vector β=[β ₁ ,…β _N ] ^T represents the off-grid error parameter vector, off-grid error parameter

n=1,...,N, s ^t ∈C ^N×1 represents the sparse incoming wave signal vector of length N; L is the number of snapshots, w ^t represents white Gaussian noise with mean 0 and variance 1/λ vector; matrices Y=[y ¹ , y ² ,...,y ^L ], S=[s ¹ ,s ² ,...,s ^L ] and W=[w ¹ ,w ² ,...,w ^L ] represent receiving, respectively The set matrix of signal vector, incoming wave signal vector and Gaussian white noise vector, the matrix S has the same sparse characteristics; A and E are the abbreviations of the array manifold matrix A(θ) and the error matrix E(θ) respectively; θ is The angle space of the incoming signal;

Then the out-of-grid DOA bilinear system model is decomposed by Bayesian formula;

和

均表示稀疏来波信号向量，数学表达为

δ(·)为delta函数；

表示在第t个快拍中，接收信号向量y^t的似然函数，似然函数表示为复高斯形式

其中

表示均值为

协方差矩阵为λ^-1I的复高斯分布，I表示对角阵；概率分布

代表向量

的先验分布，γ表示超先验参数向量；所有来波信号向量

t＝1:L都具有共同的稀疏特性，

共享同一个超先验参数向量γ；离格误差参数向量β服从均匀分布p(β)＝U[-Δ_θ/2,Δ_θ/2]，噪声方差的先验分布假设为p(λ)＝1/λ；Among them, the sparse incoming wave signal vector s ^t is defined as two identical equivalent vectors

and

Both represent the sparse incoming wave signal vector, and the mathematical expression is

δ( ) is a delta function;

Represents the likelihood function of the received signal vector y ^t in the t-th snapshot, and the likelihood function is expressed as a complex Gaussian

in

means that the mean is

The covariance matrix is a complex Gaussian distribution of λ ^-1 I, where I represents a diagonal matrix; probability distribution

representative vector

The prior distribution of , γ represents the super-prior parameter vector; all incoming wave signal vectors

t=1:L all have the same sparse property,

Share the same super-prior parameter vector γ; the off-grid error parameter vector β obeys the uniform distribution p(β)=U[ _-Δθ /2, _Δθ /2], the prior distribution of noise variance is assumed to be p(λ) =1/λ;

B: Estimate the sparse incoming wave signal S;

C: Estimate the off-grid error parameter vector β;

D: According to the obtained sparse incoming wave signal S and the off-grid error parameter vector β, the off-grid DOA bilinear system model obtained in step A is used to calculate the high-precision angle of arrival.

In the step B, the approximate message passing algorithm is applied to estimate the sparse incoming wave signal S.

In the described step C, the application adopts the expectation maximization method to estimate the off-grid error parameter vector β

In the described step A, the establishment steps of the off-grid DOA bilinear system model are as follows,

A1: Suppose there is a uniform linear array composed of M non-directional array elements, and there are K narrowband far-field incoming wave signals, then the receiving vector of the t-th snapshot is expressed as

为阵列流型矩阵，符号C^M×K表示维度为M×K的复值矩阵，

表示第t个快拍下的真实来波信号向量，C^K×1表示维度为K×1的复值向量，w^t表示均值为0且方差为1/λ的高斯白噪声向量，λ表示噪声精度；Among them, y ^t ∈ C ^M×1 is the received signal vector of all M array elements in the t-th snapshot, and the symbol C ^M×1 represents a complex-valued vector with dimension M×1,

is an array manifold matrix, the symbol C ^M×K represents a complex-valued matrix of dimension M×K,

Represents the real incoming wave signal vector in the t-th snapshot, C ^K×1 represents a complex-valued vector with dimension K×1, w ^t represents a Gaussian white noise vector with mean 0 and variance 1/λ, λ represents noise precision;

分解为向量形式

A2: The array manifold matrix

Decompose into vector form

表示第k个来波信号的导向矢量，分解表示为

j表示虚数单位，

表示第k个来波信号的真实波达角度，k＝1,…,K；来波信号的角度空间θ被划分为N个网格，表示为向量形式θ＝[θ₁,…,θ_N]，每个网格大小为Δ_θ＝π/N；每个角度θ_n均对应于一个潜在的来波信号源

n＝1,…,N；来波信号方向在角度域是稀疏分布的，即N＞＞K，构成稀疏来波信号向量

where the vector

Represents the steering vector of the kth incoming signal, decomposed and expressed as

j represents the imaginary unit,

Represents the true angle of arrival of the kth incoming wave signal, k=1,...,K; the angle space θ of incoming wave signal is divided into N grids, which are expressed in vector form θ=[θ ₁ ,...,θ _N ], the size of each grid is Δ _θ = π/N; each angle θ _n corresponds to a potential incoming signal source

n=1,...,N; the direction of the incoming wave signal is sparsely distributed in the angle domain, that is, N >> K, which constitutes a sparse incoming wave signal vector

重构为A(θ)＝[a(θ₁),...,a(θ_N)]∈C^M×N；A3: The array manifold matrix

Reconstructed as A(θ)=[a(θ ₁ ),...,a(θ _N )]∈C ^M×N ;

与任何一个网格都不完全重合，即假定存在离格误差参数

根据离格模型，将导向矢量

一阶泰勒展开为：A4: Introduce the off-grid model, assuming the true angle of arrival of the signal source

Does not exactly coincide with any grid, i.e. assumes that there is an off-grid error parameter

According to the out-of-cell model, the steering vector

The first-order Taylor expansion is:

向量

表示误差向量，其中

表示对导向矢量

依波达角度

求导数；下标n_k表示距离真实波达角度

最近的网格，

表示网格n_k对应的波达角度；where the out-of-frame error parameter β _n is defined as

vector

represents the error vector, where

Represents the pair steering vector

Ipodar angle

Derivative; subscript n _k represents the distance from the true arrival angle

the nearest grid,

represents the arrival angle corresponding to grid n _k ;

A5: According to the first-order Taylor expansion, rewrite the received signal vector of the t-th snapshot as:

y ^t =(A(θ)+E(θ)diag(β))s ^t +w ^t ;

Among them, the matrix E(θ)=[e(θ ₁ ),…,e(θ _N )] represents the error matrix, the vector β=[β ₁ ,…β _N ] ^T represents the off-grid error parameter vector, s ^t ∈ C ^N×1 represents the sparse incoming wave signal vector of length N; when there are L snapshots in the DOA estimation model,

A6: When there are L snapshots in the DOA estimation model, write the received signal vector from the t-th snapshot in step A5 as multiple observations:

Y=(A(θ)+E(θ)diag(β))S+W;

Among them, the matrices Y=[y ¹ , y ² ,...,y ^L ], S=[s ¹ ,s ² ,...,s ^L ] and W=[w ¹ ,w ² ,...,w ^L ] represent the receiving The set matrix of the signal vector, the incoming wave signal vector and the Gaussian white noise vector, the matrix S has the same sparse feature, that is, the position of the non-zero elements in each column is the same, and the angle space θ of the wave signal is a fixed value;

A7: The array manifold matrix A(θ) and the error matrix E(θ) are abbreviated as A and E respectively, and the multiple observation form in step A6 is abbreviated as

Y=(A+Ediag(β))S+W; (1)

A8: Since the expression (A+Ediag(β)) in formula (1) has only the outlier error parameter vector β as unknown, it is redefined as

Φ{β}=A+Ediag(β)=A+∑ _n β _n E _n ; (2)

Among them, β _n represents the off-grid error parameter, and E _n is defined as E _n =[0,…,e(θ _n ),…,0]∈C ^M×N , the DOA estimation model in equation (1) is transformed into bilinear form

Y=(A+∑ _n β _n E _n )S+W=Φ{β}S+W (3)

时，β_n取非零值，向量β和矩阵S具有完全相同的稀疏特性；Among them, S=[s ¹ ,...,s ^L ] represents L sparse incoming wave signal vectors to be obtained; only when

When β _n takes a non-zero value, the vector β and the matrix S have exactly the same sparse characteristics;

A9: According to the off-grid DOA bilinear system model shown in formula (3), the joint distribution function of all known and unknown variables in the system is factorized using the Bayesian formula by using the constraint relationship between the data ,

Described step B includes the following specific steps:

为长度为N的全零向量，中间变量

初始化为1，噪声精度初始化为λ＝1，超先验参数向量γ初始化为长度N的全1向量；然后进入步骤B2；B1: Initialize the off-grid error parameter vector β and intermediate variables

is an all-zero vector of length N, the intermediate variable

It is initialized to 1, the noise accuracy is initialized to λ=1, and the super-priority parameter vector γ is initialized to an all-1 vector of length N; then enter step B2;

的后验概率为高斯形式

其中均值

和方差

分别计算为B2: Calculated to sparse wave signal vector

The posterior probability of is in Gaussian form

where the mean

and variance

are calculated as

Then enter step B3;

In the formula, the operator symbol <·> represents the average of the vectors, I represents the diagonal matrix, and y represents the received signal vector;

和方差

以及步骤B1中初始化的中间变量

和

计算中间变量

和

分别为B3: Use the mean value obtained in step B2

and variance

and the intermediate variables initialized in step B1

and

Calculate intermediate variables

and

respectively

Then enter step B4;

和

计算中间变量

的后验分布为高斯形式

其中均值

和方差

分别计算为B4: Utilize the super-prior parameters γ initialized in step B1 and the intermediate variables obtained in step B3

and

Calculate intermediate variables

The posterior distribution of is in Gaussian form

where the mean

and variance

are calculated as

Then enter step B5;

的均值

和方差

计算超先验参数γ为B5: Use the intermediate variable obtained in step B4

mean of

and variance

Calculate the hyper-prior parameter γ as

Then enter step B6; wherein L represents the number of snapshots in DOA estimation;

的均值

和方差

以及步骤B3得到的中间变量

和

计算中间变量

和

B6: Use the intermediate variable obtained in step B4

mean of

and variance

and the intermediate variable obtained in step B3

and

Calculate intermediate variables

and

Then enter step B7;

B7: Determine whether the set number of iterations is reached, if so, end the iteration; if not, return to step B2; finally determine the sparse incoming wave signal vector

The steps B2 to B6 are iterative processes, and the number of iterations is set to 15 times.

Described step C includes following concrete steps:

C1: Initialize the out-of-grid error parameter β as an all-zero vector of length N, and substitute β into formula (2) to construct a matrix Φ{β}; then enter step C2;

t＝1:L计算中间变量矩阵D为C2: Use the variance obtained in step B2

t=1:L calculates the intermediate variable matrix D as

Then enter step C3;

和公式(2)中定义的E_n,n＝1,2,...,N分别构造矩阵H∈C^N×N和向量α∈C^N×1为C3: Use the sparse incoming wave signal vector obtained in step B7

and E _n , n=1,2,...,N defined in formula (2), respectively construct matrix H∈C ^N×N and vector α∈C ^N×1 as

α＝(α₁…α_i…α_N)^T；

α=(α ₁ ...α _i ...α _N ) ^T ;

Among them, the elements H _ij ,i,j=1,2,...,N in the matrix and vector, α _i ,i=1,2,...,N,

i, j represent the element labels of the matrix H and the vector α, the matrices E _i , E _j have the same definition as the matrix _En in the formula (2); then enter step C4;

C4: find the estimated value of the off-grid error parameter vector β, β=H ^-1 α; then enter step C5;

C5: Determine whether the set number of iterations has been reached, and if so, end the iteration; if not, return to step C2; and finally obtain the off-grid error parameter vector β.

As mentioned above, steps C2 to C4 are iterative processes, and the number of iterations is set to 5 times.

The invention converts the DOA estimation model into a bilinear problem, and uses the bilinear VAMP algorithm to solve the problem. In the present invention, the off-grid DOA estimation model is first converted into a bilinear problem, and the bilinear VAMP algorithm is used to deduce and operate the bilinear model; then a DOA estimation simulation model is established to compare the proposed algorithm with the existing methods in the literature. , the results show that the algorithm proposed in the present invention has obvious performance advantages in scenarios with a small number of snapshots such as airborne radar.

Description of drawings

Fig. 1 is the schematic flow chart of the present invention;

Fig. 2 is the comparison diagram of estimated value and real value in a Monte Carlo simulation of the present invention;

Fig. 3 is the change curve diagram of the present invention and normalized mean square error with the number of snapshots;

4 is a graph showing the variation of estimation performance with signal-to-noise ratio when the number of snapshots L is 2 in the present invention and an existing algorithm;

5 is a graph showing the variation of the estimated performance with the signal-to-noise ratio of the present invention and the existing algorithm when the number of snapshots L is 4;

6 is a graph showing the variation of the estimated performance with the signal-to-noise ratio of the present invention and the existing algorithm when the number of snapshots L is 6;

FIG. 7 is a graph showing the variation of the estimation performance with the signal-to-noise ratio when the number of snapshots L is 8 in the present invention and the existing algorithm.

Detailed ways

Below in conjunction with accompanying drawing and embodiment, the present invention is described in detail:

As shown in FIG. 1, the high-precision estimation method of the angle of arrival according to the present invention includes the following steps:

A: Establish an off-grid DOA bilinear system model Y=(A+Σ _n β _n E _n )S+W=Φ{β}S+W, and use Bayesian formula for factorization;

A1: Suppose there is a uniform linear array composed of M non-directional array elements, and there are K narrowband far-field incoming wave signals, then the receiving vector of the t-th snapshot is expressed as

为阵列流型矩阵，符号C^M×K表示维度为M×K的复值矩阵，

表示第t个快拍下的真实来波信号向量，C^K×1表示维度为K×1的复值向量，w^t表示均值为0且方差为1/λ的高斯白噪声向量，其中λ表示噪声精度。A2：将阵列流型矩阵

分解为向量形式：Among them, y ^t ∈ C ^M×1 is the received signal vector of all M array elements in the t-th snapshot, and the symbol C ^M×1 represents a complex-valued vector with dimension M×1,

is an array manifold matrix, the symbol C ^M×K represents a complex-valued matrix of dimension M×K,

Represents the real incoming wave signal vector of the t-th snapshot, C ^K×1 represents the complex-valued vector with dimension K×1, w ^t represents the Gaussian white noise vector with mean 0 and variance 1/λ, where λ represents Noise Accuracy. A2: The array manifold matrix

Decompose into vector form:

表示第k个来波信号的导向矢量，可以分解表示为

j表示虚数单位，

表示第k个来波信号的真实波达角度，k＝1,…,K；由于DOA估计的目标是利用接收向量y^t估计所有K个来波信号的波达角度

和幅度，因此在本发明中，将来波信号的角度空间θ划分为N个网格，表示为向量形式θ＝[θ₁,…,θ_N]，每个网格大小为Δ_θ＝π/N。假设每个角度θ_n均对应于一个潜在的来波信号源

n＝1,…,N；由于来波信号方向在角度域是稀疏分布的，即N＞＞K，从而构成稀疏来波信号向量

where the vector

Represents the steering vector of the kth incoming signal, which can be decomposed and expressed as

j represents the imaginary unit,

Represents the true angle of arrival of the k-th incoming signal, k=1,...,K; since the goal of DOA estimation is to use the receiving vector y ^t to estimate the angle of arrival of all K incoming signals

and amplitude, so in the present invention, the angular space θ of the future wave signal is divided into N grids, which are expressed in the vector form θ=[θ ₁ ,...,θ _N ], and the size of each grid is Δ _θ =π/ N. Assume that each angle θ _n corresponds to a potential incoming wave signal source

n=1, .

重构为A(θ)＝[a(θ₁),...,a(θ_N)]∈C^M×N；A3: The array manifold matrix

Reconstructed as A(θ)=[a(θ ₁ ),...,a(θ _N )]∈C ^M×N ;

与任何一个网格都不完全重合，即假定存在离格误差参数

根据离格模型，将导向矢量

一阶泰勒展开为：A4: To improve the accuracy of DOA, the most direct method is to divide the angular space θ of the future wave signal into smaller grids, that is, to increase the number of grids N. However, the increase of the number of grids N will lead to an increase in complexity, and the correlation of the matrix A(θ) will increase, and the matrix A(θ) with high correlation will lead to poor estimation performance, causing iterative classes such as message passing. The divergence of the algorithm causes the estimation to fail. In order to avoid the above problems, an off-grid model is introduced in the present invention, that is, the real angle of arrival of the signal source is assumed

Does not exactly coincide with any grid, i.e. assumes that there is an off-grid error parameter

According to the out-of-cell model, the steering vector

The first-order Taylor expansion is:

向量

表示误差向量，其中

表示对导向矢量

依波达角度

求导数；下标n_k表示距离真实波达角度

最近的网格，

表示网格n_k对应的波达角度。where the out-of-frame error parameter β _n is defined as

vector

represents the error vector, where

Represents the pair steering vector

Ipodar angle

Derivative; subscript n _k represents the distance from the true arrival angle

the nearest grid,

represents the angle of arrival corresponding to the grid n _k .

A5: According to the first-order Taylor expansion, rewrite the received signal vector of the t-th snapshot as:

y ^t =(A(θ)+E(θ)diag(β))s ^t +w ^t ;

Among them, the matrix E(θ)=[e(θ ₁ ),…,e(θ _N )] represents the error matrix, the vector β=[β ₁ ,…β _N ] ^T represents the off-grid error parameter vector, s ^t ∈ C ^N×1 represents the sparse incoming wave signal vector of length N to be obtained.

A6: When there are L snapshots in the DOA estimation model, write the above formula in the form of multiple observations

Y=(A(θ)+E(θ)diag(β))S+W;

Among them, the matrices Y=[y ¹ , y ² ,...,y ^L ], S=[s ¹ ,s ² ,...,s ^L ] and W=[w ¹ ,w ² ,...,w ^L ] represent the receiving The set matrix of the signal vector, the incoming wave signal vector and the Gaussian white noise vector, and the matrix S has the same sparse feature, that is, the position of the non-zero elements in each column is the same. The angle space θ of the incoming wave signal in the above formula is a fixed value.

A7: The array manifold matrix A(θ) and the error matrix E(θ) are abbreviated as A and E respectively, then the DOA estimation model is abbreviated as

Y=(A+Ediag(β))S+W; (1)

It can be seen from the observation that the unknown parameters in the above DOA estimation model are S and β. In the field of signal processing, the above DOA estimation model belongs to a bilinear model that uses a set of observations to simultaneously estimate two sets of independent parameters.

A8: Since the expression (A+Ediag(β)) in formula (1) only has the out-of-frame error parameter vector β as unknown, it can be redefined as

Φ{β}=A+Ediag(β)=A+∑ _n β _n E _n ; (2)

Among them, β _n represents the off-grid error parameter, and E _n is defined as E _n =[0,…,e(θ _n ),…,0]∈C ^M×N , then the DOA estimation model in equation (1) is transformed into in bilinear form

Y=(A+Σ _n β _n E _n )S+W=Φ{β}S+W (3)

时β_n取非零值，从而可确定向量β和矩阵S具有完全相同的稀疏特性，即非零元素位置一样。Wherein, S=[s ¹ ,...,s ^L ] represents L sparse incoming wave signal vectors to be obtained. According to the physical meaning of the vector β in equation (1), it can be seen that only when

When β _n takes a non-zero value, it can be determined that the vector β and the matrix S have exactly the same sparse characteristics, that is, the positions of the non-zero elements are the same.

A9: According to the off-grid DOA bilinear system model shown in formula (3), the joint distribution function of all known and unknown variables in the system is factorized using the Bayesian formula by using the constraint relationship between the data ,

和

三者完全相等，即

均表示稀疏来波信号向量，数学表达为

此处δ(·)为delta函数；

表示在第t个快拍中，接收信号向量y^t的似然函数，由于噪声服从方差为1/λ的复高斯分布，则似然函数可以表示为复高斯形式

其中

表示均值为

协方差矩阵为λ^-1I的复高斯分布，I表示对角阵；概率分布

代表向量

的先验分布，γ表示超先验参数向量。由于本发明假定对于所有来波信号向量

t＝1:L都具有共同的稀疏特性，则

共享同一个超先验参数向量γ。离格误差参数向量β服从均匀分布p(β)＝U[-Δ_θ/2,Δ_θ/2]，噪声方差的先验分布假设为p(λ)＝1/λ。根据(4)式所示双线性模型的因式分解，可以利用双线性近似消息传递和期望最大化算法对双线性模型进行算法推导。具体步骤分为稀疏来波信号S和离格误差参数β估计两部分。In the above formula, the sparse incoming wave signal vector s ^t is defined as two identical equivalent vectors

and

The three are completely equal, namely

Both represent the sparse incoming wave signal vector, and the mathematical expression is

Here δ( ) is the delta function;

Represents the likelihood function of the received signal vector y ^t in the t-th snapshot. Since the noise follows a complex Gaussian distribution with a variance of 1/λ, the likelihood function can be expressed as a complex Gaussian form

in

means that the mean is

The covariance matrix is a complex Gaussian distribution of λ ^-1 I, where I represents a diagonal matrix; probability distribution

representative vector

The prior distribution of , γ denotes the hyper-prior parameter vector. Since the present invention assumes that for all incoming signal vectors

t=1:L all have common sparse characteristics, then

share the same hyper-prior parameter vector γ. The outlier error parameter vector β obeys the uniform distribution p(β)=U[ _-Δθ /2, _Δθ /2], and the prior distribution of noise variance is assumed to be p(λ)=1/λ. According to the factorization of the bilinear model shown in equation (4), the bilinear model can be algorithmically derived using the bilinear approximation message passing and expectation maximization algorithm. The specific steps are divided into two parts: the sparse incoming wave signal S and the estimation of the off-grid error parameter β.

B: Estimate the sparse incoming wave signal S;

进行估计，具体估计过程如下：B1: If the off-grid error parameter vector β is known, the matrix Φ{β} is also known, so that the bilinear model shown in formula (3) degenerates into a single linear model. When solving a single linear problem, the approximate message passing algorithm is applied to the sparse incoming signal vector

The specific estimation process is as follows:

为长度为N的全零向量，中间变量

初始化为1，噪声精度初始化为λ＝1，超先验参数向量γ初始化为长度N的全1向量；然后进入步骤B2；Initialize the out-of-frame error parameter vector β and intermediate variables

is an all-zero vector of length N, the intermediate variable

It is initialized to 1, the noise accuracy is initialized to λ=1, and the super-priority parameter vector γ is initialized to an all-1 vector of length N; then enter step B2;

的后验概率为高斯形式

其中均值

和方差

分别计算为B2: Calculated to sparse wave signal vector

The posterior probability of is in Gaussian form

where the mean

and variance

are calculated as

Then enter step B3;

In the formula, the operator symbol <·> represents the average of the vectors, I represents the diagonal matrix, and y represents the received signal vector. The average of variance vectors is used in the present invention, which can improve the robustness of the algorithm.

和方差

以及步骤B1中初始化的中间变量

和

计算中间变量

和

分别为B3: Use the mean value obtained in step B2

and variance

and the intermediate variables initialized in step B1

and

Calculate intermediate variables

and

respectively

Then enter step B4;

和

计算中间变量

的后验分布为高斯形式

其中均值

和方差

分别计算为B4: Utilize the super-prior parameters γ initialized in step B1 and the intermediate variables obtained in step B3

and

Calculate intermediate variables

The posterior distribution of is in Gaussian form

where the mean

and variance

are calculated as

Then enter step B5;

的均值

和方差

计算超先验参数γ为B5: Use the intermediate variable obtained in step B4

mean of

and variance

Calculate the hyper-prior parameter γ as

Then enter step B6; wherein L represents the number of snapshots in the DOA estimation.

的均值

和方差

以及步骤B3得到的中间变量

和

计算中间变量

和

B6: Use the intermediate variable obtained in step B4

mean of

and variance

and the intermediate variable obtained in step B3

and

Calculate intermediate variables

and

Then enter step B7;

B7: Determine whether the set number of iterations is reached, if so, end the iteration; if not, return to step B2; finally determine the sparse incoming wave signal vector

In this embodiment, steps B2 to B6 are iterative processes, and the number of iterations is set to 15 times.

C: Estimate the off-grid error parameter vector β;

The present invention adopts the expectation maximization method for the estimation of the off-grid error parameter vector β, and the calculation steps are as follows:

C1: Initialize the out-of-grid error parameter β as an all-zero vector of length N, and substitute β into formula (2) to construct a matrix Φ{β}; then enter step C2;

t＝1:L计算中间变量矩阵D为C2: Use the variance obtained in step B2

t=1:L calculates the intermediate variable matrix D as

Then enter step C3;

和公式(2)中定义的E_n,n＝1,2,...,N分别构造矩阵H∈C^N×N和向量α∈C^N×1为C3: Use the sparse incoming wave signal vector obtained in step B7

and E _n , n=1,2,...,N defined in formula (2), respectively construct matrix H∈C ^N×N and vector α∈C ^N×1 as

α＝(α₁…α_i…α_N)^T；

α=(α ₁ ...α _i ...α _N ) ^T ;

Among them, the elements H _ij ,i,j=1,2,...,N in the matrix and vector, α _i ,i=1,2,...,N,

Wherein i, j represent the element labels of matrix H and vector α, where matrix E _i , E _j have the same definition as matrix E _n in formula (2); then enter step C4;

C4: Find the estimated value of the off-grid error parameter vector β, β=H ^-1 α; then go to step C5;

C5: Determine whether the set number of iterations has been reached, and if so, end the iteration; if not, return to step C2; and finally obtain the off-grid error parameter vector β.

In this embodiment, steps C2 to C4 are iterative processes, and the number of iterations is set to 5 times.

D: Using the sparse incoming wave signal vector S obtained in step B, and the off-grid error parameter β obtained in step C, combined with the DOA estimation model established in step A Y=(A+∑ _n β _n E _n )S +W=Φ{β}S+W, and finally a high-precision estimate of the angle of arrival is obtained.

B步的其余运算均为标量形式，复杂度为

其中符号

表示复杂度的阶数。C步骤需要计算矩阵求逆，复杂度为

所以本发明的整体复杂度可写为

In terms of complexity, the present invention can be divided into two steps B and C, wherein step B2 requires matrix multiplication Φ ^T {β}Φ{β}, which has the complexity

The rest of the operations in step B are in scalar form, and the complexity is

where the symbol

Indicates the order of complexity. Step C needs to calculate the matrix inversion, and the complexity is

So the overall complexity of the present invention can be written as

In order to illustrate the performance of the present invention, a simulation environment can be established to carry out numerical simulation, and a comparison is made with the existing domestic and foreign algorithms. The simulation environment selects a uniform linear array, the number of array elements is M=30, then the number of virtual grids is N=30, that is, the grid width is 6°, the number of snapshots is set to L=2-12, and the threshold _γmax is set = 10 ³ . In order to measure the effectiveness of the present invention, the existing gridded sparse Bayesian learning (OGSBL), root-seeking SBL (RTSBL) and Gauss-Selder root (GSROOT) are selected for comparison in the simulation, and Kramer is also selected. Luojie (CRB) as a reference.

Figure 2 is a comparison diagram of the estimated value Estimated and the true value of True Angle in a Monte Carlo simulation. The true angle selected in the simulation is [-28.6,-18.6,3.5,15.6,31.7]. In Fig. ₂ , the abscissa is the angle, and the ordinate is the absolute value of the estimated vector s2. In this simulation, the signal-to-noise ratio SNR=30dB is set. It can be seen from FIG. 2 that the present invention can accurately capture the off-grid angle and accurate power of five incoming wave signals, and the angle estimation result approaches 0 when there is no incoming wave signal. ,

FIG. 3 is a graph showing the variation of the normalized mean square error (NMSE) of the present invention and the existing algorithm with the number of snapshots. The numerical simulation results are the average of 500 Monte Carlo simulations. In each simulation, the direction of arrival is based on [-30,-18,6,18,30], and a uniform random angular offset of (-3,3) is superimposed. In Fig. 3, the signal-to-noise ratio SNR=20dB is set, and the number of incoming signals is K=5. It can be seen from the simulation results that the performance of the present invention, namely the Bad-VAMP shown in the figure, is close to that of the existing RTSBL and GSROOT when the number of snapshots is large. However, with the reduction of the number of snapshots, the performance of the RTSBL and GSROOT methods is obviously deteriorated, so that the present invention shows a large performance gain. It can also be seen from Figure 3 that the performance of the OGSBL algorithm is poor because the grid bias is not considered, which proves the necessity of the out-of-grid estimation algorithm.

FIG. 4 to FIG. 7 are graphs showing the variation of estimation performance with signal-to-noise ratio when the number of snapshots L is 2, 4, 6 and 8, respectively. From Figures 4 to 7, the same conclusion as Figure 3 can be drawn, that is, the present invention has obvious performance advantages when the number of snapshots is small. To achieve the same estimation performance, the present invention requires fewer snapshots, and is more suitable for fast time-varying scenarios such as airborne radar.

Claims (8)

1. A high-precision estimation method of the angle of arrival is characterized in that: the method sequentially comprises the following steps:

a: firstly, establishing a lattice-separated DOA bilinear system model Y (A + ∑ Y)_nβ_nE_n)S+W＝Φ{β}S+W；

Where Φ { β } ═ a + Ediag (β) ═ a + ∑_nβ_nE_n，E_n＝[0,…,e(θ_n),…,0]∈C^M×NThe de-lattice DOA bilinear system model is a shorthand for a received signal vector Y ═ a (θ) + E (θ) diag (β)) S + W; y ═ a (θ) + E (θ) diag (β)) S + W is Y^t＝(A(θ)+E(θ)diag(β))s^t+w^tMultiple observation modalities of (1); the device is provided with a uniform linear array formed by M non-directional array elements, and K narrow-band far-field incoming wave signals, y^tDenotes a received signal vector at the t-th snapshot, where a (θ) ═ a (θ)₁),...,a(θ_N)]∈C^M×N，C^M×KRepresenting a complex matrix with dimension M × K; let the angle space theta of the incoming wave signal be divided into N grids and expressed as a vector form theta ═ theta₁,…,θ_N]Each grid size is Δ_θpi/N; matrix E (θ) ═ E (θ)₁),…,e(θ_N)]Representing an error matrix, vector β ═ β₁,…β_N]^TRepresenting vectors of, parameters of, off-grid errors

s^t∈C^N×1Representing a sparse incoming wave signal vector to be solved with the length of N; l is the number of snapshots, w^tRepresenting a Gaussian white noise vector with a mean of 0 and a variance of 1/λ; matrix Y ═ Y¹,y²,…,y^L]、S＝[s¹,s²,…,s^L]And W ═ W¹,w²,…,w^L]Respectively representing a set matrix of a received signal vector, an incoming wave signal vector and a Gaussian white noise vector, wherein the matrix S has the same sparse characteristic; a and E are shorthand for array flow matrix A (theta) and error matrix E (theta) respectively; theta is the angle space of the incoming wave signal;

then carrying out factorization on the off-grid DOA bilinear system model by using a Bayesian formula;

wherein, a sparse incoming wave signal vector s is obtained^tDefined as two identical equivalent vectors

And

all represent sparse incoming wave signal vectors, expressed mathematically as

δ (·) is a delta function;

indicating that in the t-th snapshot, the received signal vector y^tOf the likelihood function expressed in complex gaussian form

Wherein

Represents a mean value of

Covariance matrix of λ^-1A complex gaussian distribution of I, I representing a diagonal matrix; probability distribution

Representative vector

γ represents a vector of the prior parameters; all incoming wave signal vectors

All have a common sparsity characteristic that is,

sharing the same super-prior parameter vector gamma; the off-grid error parameter vector beta is subject to a uniform distribution of p (beta) ═ U [ -delta ]_θ/2,Δ_θ/2]The prior distribution of noise variance is assumed to be p (λ) ═ 1/λ;

b: estimating a sparse incoming wave signal S;

c: estimating an off-grid error parameter vector beta;

d: and D, according to the obtained sparse incoming wave signal S and the off-grid error parameter vector beta, performing high-precision calculation on the arrival angle through the off-grid DOA bilinear system model obtained in the step A.

2. The high-precision estimation method of the angle of arrival according to claim 1, characterized in that: in the step B, an approximate message transfer algorithm is applied to estimate the sparse incoming wave signal S.

3. The high-precision estimation method of the angle of arrival according to claim 1, characterized in that: in the step C, the off-grid error parameter vector beta is estimated by using an expectation maximization method.

4. The high-precision estimation method of the angle of arrival according to claim 1, characterized in that: in the step A, the off-grid DOA bilinear system model is established as follows,

a1: the receiving vector under the t-th snapshot is expressed as

Wherein, y^t∈C^M×1For the received signal vectors of all M array elements in the t-th snapshot, symbol C^M×1Representing a complex-valued vector of dimension mx 1,

for an array flow pattern matrix, symbol C^M×KRepresenting a complex-valued matrix of dimension M x K,

representing the true incoming wave signal vector, C, at the t-th snapshot^K×1Representing a complex-valued vector of dimension K x 1, w^tRepresenting a Gaussian white noise vector with a mean value of 0 and a variance of 1/lambda, lambda representing the noise precision;

a2: array flow pattern matrix

Decomposition into vector form

Wherein the vector

A steering vector representing the k-th incoming wave signal, decomposed as

j represents the unit of an imaginary number,

represents the true arrival angle of the kth incoming wave signal, K is 1, …, K; the angle space theta of the incoming wave signal is divided into N grids and expressed as a vector form theta ═ theta₁,…,θ_N]Each grid size is Δ_θpi/N; each angle theta_nAll correspond to a potential incoming wave signal source

The directions of incoming wave signals are sparsely distributed in an angle domain, namely N > K, and a sparse incoming wave signal vector is formed

A3: array flow pattern matrix

Reconstructed as a (θ) ═ a (θ)₁),...,a(θ_N)]∈C^M×N；

A4: introducing a lattice-separating model, and assuming the true angle of arrival of the signal source

Not completely coincident with any of the meshes, i.e. assuming the presence of a misstep parameter

According to the off-grid model, the guide vector is divided into

The first order Taylor expansion is:

wherein the off-grid error parameter beta_nIs defined as

(Vector)

Represents an error vector in which

Representing a pair of steering vectors

Angle of arrival

A derivative is obtained; subscript n_kRepresenting true angle of arrival of distance

The nearest grid of the grid is,

representing a grid n_kA corresponding angle of arrival;

a5: according to the first-order Taylor expansion, the received signal vector at the t-th snapshot is rewritten as:

y^t＝(A(θ)+E(θ)diag(β))s^t+w^t；

wherein the matrix E (θ) ═ E (θ)₁),…,e(θ_N)]Representing an error matrix, vector β ═ β₁,…β_N]^TRepresenting the off-grid error parameter vector, s^t∈C^N×1Representing a sparse incoming wave signal vector to be solved with the length of N; when there are L snapshots in the DOA estimation model,

a6: when there are L snapshots in the DOA estimation model, the received signal vector at the t-th snapshot in step a5 is written as a multiple observation form:

Y＝(A(θ)+E(θ)diag(β))S+W；

wherein the matrix Y ═ Y¹,y²,…,y^L]、S＝[s¹,s²,…,s^L]And W ═ W¹,w²,…,w^L]Respectively representing a set matrix of a received signal vector, an incoming wave signal vector and a Gaussian white noise vector, wherein the matrix S has the same sparse characteristic, namely the positions of non-zero elements in each column are the same, and the angle space theta of a wave signal is a fixed value;

a7: the array flow pattern matrix A (theta) and the error matrix E (theta) are respectively abbreviated as A and E, and the multiple observation form in the step A6 is abbreviated as

Y＝(A+Ediag(β))S+W； (1)

A8: since the expression (a + Ediag (β)) in equation (1) is unknown only for the miselike error parameter vector β, it is redefined as

Φ{β}＝A+Ediag(β)＝A+Σ_nβ_nE_n； (2)

Wherein, beta_nRepresenting the off-grid error parameter, E_nIs defined as E_n＝[0,…,e(θ_n),…,0]∈C^M×NConverting the DOA estimation model in the formula (1) into a bilinear form

Y＝(A+∑_nβ_nE_n)S+W＝Φ{β}S+W (3)

Wherein S ═ S¹,...,s^L]Representing L sparse incoming wave signal vectors to be solved; only when

When is beta_nTaking a nonzero value, wherein the vector beta and the matrix S have the same sparse characteristic;

a9: according to the off-grid DOA bilinear system model shown in the formula (3), the joint distribution function of all known and unknown variables in the system is factorized by a Bayes formula by utilizing the constraint relation among data,

5. the high-precision estimation method of the angle of arrival according to claim 1, characterized in that: the step B comprises the following specific steps:

b1: initializing a de-grid error parameter vector beta and an intermediate variable

Is an all-zero vector of length N, an intermediate variable

Initializing to 1, initializing the noise precision to lambda being 1, and initializing a super-prior parameter vector gamma to a full 1 vector with the length N; then step B2 is entered;

b2 calculating sparse wave signal vector

The posterior probability of (1) is in Gaussian form

Wherein the mean value

Sum variance

Are respectively calculated as

Then step B3 is entered;

in the formula, an operation symbol < > represents averaging of vectors, I represents a diagonal matrix, and y represents a received signal vector;

b3 average value obtained in step B2

Sum variance

And intermediate variables initialized in step B1

And

calculating intermediate variables

And

are respectively as

Then step B4 is entered;

b4 using the initial gamma parameter in step B1 and the intermediate variables obtained in step B3

And

calculating intermediate variables

The posterior distribution of (A) is in the form of Gaussian

Wherein the mean value

Sum variance

Are respectively calculated as

Then step B5 is entered;

b5 intermediate variables obtained by step B4

Mean value of

Sum variance

Calculating a super-prior parameter gamma of

Then step B6 is entered; wherein L represents the number of snapshots in the DOA estimation;

b6 intermediate variables obtained by step B4

Mean value of

Sum variance

And intermediate variables obtained in step B3

And

calculating intermediate variables

And

then step B7 is entered;

b7, judging whether the set iteration times is reached, if so, ending the iteration; if not, return toGo back to step B2; finally determining sparse incoming wave signal vector

6. The high-precision estimation method of the angle of arrival according to claim 5, characterized in that: the steps B2 to B6 are iterative processes, and the number of iterations is set to 15.

7. The high-precision estimation method of the angle of arrival according to claim 5, characterized in that: the step C comprises the following specific steps:

c1, initializing a lattice error parameter beta as a full 0 vector with the length of N, and substituting the beta into the formula (2) to construct a matrix phi { beta }; then proceed to step C2;

c2 using the variance obtained in step B2

Calculating an intermediate variable matrix D of

Then proceed to step C3;

c3 using the sparse incoming wave signal vector obtained in step B7

And E defined in formula (2)_nN is 1,2, N forms a matrix H e C, respectively^N×NThe sum vector α ∈ C^N×1Is composed of

α＝(α₁…α_i…α_N)^T；

Wherein the element H in the matrix and vector_ij,i,j＝1,2,...,N，α_i,i＝1,2,...,N，

i, j denote the element indices of matrix H and vector α, matrix E_i,E_jAnd the matrix E in the formula (2)_nHave the same definition; then proceed to step C4;

c4, obtaining the estimated value of the error parameter vector beta, where beta is H^-1α; then proceed to step C5;

c5, judging whether the set iteration times are reached, if so, ending the iteration; if not, returning to the step C2; finally, the off-grid error parameter vector beta is obtained.

8. The high-precision estimation method of the angle of arrival according to claim 7, characterized in that: the steps C2 to C4 are iterative processes, and the number of iterations is set to 5.

Images