# Complex Multisnapshot Sparse Bayesian Learning for Offgrid DOA Estimation

Authors: Qinghua Liu; Yuanxin He; Kai Ding; Quanmin Xie
Journal: International Journal of Antennas and Propagation
 (2022)

---

## Abstract

Direction of arrival (DOA) estimation has recently been developed based on sparse signal reconstruction (SSR). Sparse Bayesian learning (SBL) is a typical method of SSR. In SBL, the two-layer hierarchical model in Gaussian scale mixtures (GSMs) has been used to model sparsity-inducing priors. However, this model is mainly applied to real-valued signal models. In order to apply SBL to complex-valued signal models, a general class of sparsity-inducing priors is proposed for complex-valued signal models by complex Gaussian scale mixtures (CGSMs), and the special cases correspond to complex versions of several classical priors are provided, which is helpful to analyze the connections with different modeling methods. In addition, the expression of the SBL form of the real- and complex-valued model is unified by parameter values, which makes it possible to generalize and improve the properties of the SBL methods. Finally, the SBL complex-valued form is applied to the offgrid DOA estimation complex-valued model, and the performance between different sparsity-inducing priors is compared. Theoretical analysis and simulation results show that the proposed algorithm can effectively process complex-valued signal models and has lower algorithm complexity.


---

## Body

## 1. Introduction

Direction of arrival (DOA) estimation of spatial signals is an important content of array signal processing research. In recent years, compressive sensing (CS) [1] and sparse signal reconstruction (SSR) [2] have been introduced into the field of DOA estimation owing to their ultra-high resolution, good robustness to noise, and low dependence on the number of snapshots. One representative method is the L1-SVD algorithm proposed in [3], which uses L1-norm to construct an SSR model, reduces the complexity of the algorithm by singular value decomposition (SVD) under multisnapshots, and has a high resolution for the correlated sources as well. Sparse Bayesian learning (SBL) [4–7] or Bayesian compressed sensing (BCS) [8, 9] is another representative method for the SSR in CS. In SBL, the signal recovery problem is formulated from a Bayesian perspective, while the sparsity information is exploited by assuming a sparse prior for the signal of interest. One merit of SBL is its flexibility in modeling sparse signals, which can improve the sparsity of its solution [10]. Therefore, research work on DOA estimation based on SBL has been gaining momentum in recent years [11, 12]. While these methods have shown advantages over the conventional ones, however, there exists still difficulties in practical applications.For these SSR methods mentioned above, which are called ongrid methods, true DOAs are assumed to lie exactly on a set of fixed sampling grids, and the existing sparse representation techniques can be directly applied. When the assumption that the true DOAs are located on a sampling grid fails, the performance of on-grid methods deteriorates due to the problem of mismatch. Although the offgrid (the distance from the true DOA to the nearest grid) error can be reduced by a dense sampling grid, the computational complexity will increase significantly. Furthermore, a dense grid may lead to the correlation between two steering vectors becoming high, which means the performance cannot be improved by the theory of CS.Accordingly, an offgrid DOA estimation model was first studied in [13], where the true DOAs are no longer constrained in the sampling grid, and they proposed a sparsity cognizant total least-squares (S-TLS) method for the perturbed compressive sensing under sparsity constraints. It has been shown in [13] that the S-TLS can yield a MAP (maximum a posteriori) optimal estimate if the matrix perturbation caused by the basis mismatch is Gaussian. But the Gaussian condition cannot be satisfied in all DOA estimations. In [14], a new algorithm named offgrid sparse Bayesian inference (OGSBI) is proposed for the offgrid case. Firstly, the true DOA steering vector is approximated by a first-order Taylor expansion to handle the off-grid problem, and it is assumed to have a uniform distribution for off-grid parameters. Based on the OGSBI method, a linear interpolation between two adjacent grids is adopted in [15] to approximate the true DOA steering vector, and they proposed a perturbed sparse Bayesian learning (PSBL) algorithm to solve the offgrid DOA estimation problem.Both ongrid and offgrid methods are grid-based methods. Another method to deal with the offgrid problem is the gridless method [16–18], which directly operates in the continuous domain so that can completely avoid the grid mismatch problem. However, this kind of method needs strong theoretical guarantees. In addition, it can only be applied to uniform or sparse linear arrays [19].The SBL algorithm mainly deals with real-valued signals. In practical applications, such as DOA estimation, the experimental signals are often in complex-valued form. The results obtained in the real-valued model cannot be directly applied to the complex-valued model. In order to apply SBL to the complex-valued model, in [20–22], the real and imaginary parts of the complex-valued signal were separated to construct a new real-valued model, and then the SBL real-valued model was applied for processing. But in [20, 21], the real and imaginary parts of the complex-valued signal were treated as two independent variables, and the different hyperparameters were specified by ignoring the correlation between the real and imaginary parts of the complex-valued signal, which has obvious theoretical defects, and the reconstruction performance of the algorithm is low. It was improved in [22], where the performance of the algorithm was improved by assigning the same hyperparameters to the real and imaginary parts of the complex-valued signal by taking advantage of the fact that they have the same sparse structures. However, by decomposing the complex-valued signal into real and imaginary parts, the dimensionality of the signal is doubled, and the sensing matrix is tripled, so the complexity of the algorithm will increase significantly.In this paper, the extension of SBL from the real-valued model to the complex-valued model in conjunction with the DOA estimation complex-valued model is considered, and the complex multisnapshot SBL (CSBL) offgrid DOA estimation algorithm is proposed. The expression of the SBL form is unified under the two models by parameter values, which enables us to generalize and improve the properties of the SBL methods. Theoretical analysis and simulation results show that the proposed algorithm can effectively process complex-valued signals with lower complexity.A comment on notation: we use boldface lowercase letters for vectors and boldface uppercase letters for matrices.⋅T, ⋅−1, ⋅∗, and ⋅H denote the transpose, inverse, conjugation, and conjugation-transpose operations, respectively. ·F and ·2 are the F-norm and 2-norm of a matrix, respectively. tr⋅ denotes the trace operation of a matrix. ⊙ describes Hadamard product operator. p⋅ represents the probability density function of random variables.


## 2. Offgrid DOA Estimation Model

Consider a uniform linear array (ULA) withN isotropic sensors and interelement spacing Δd as shown in Figure 1. Suppose that K far-field narrowband signals of frequency f0 illuminate on this ULA from directions θk,k=1,2,…,K. By defining the first sensor of the array as the reference, the received signals can be expressed as(1)Y=AθS+E,where Y=y1,y2,…,yM, S=s1,s2,…,sM, E=e1,e2,…,eM, and M is the number of snapshots. ym=y1m,y2m,…,yNmT∈ℂN×1, Aθ=aθ1,aθ2,…,aθK∈ℂN×K, aθk=1,expj2πf0Δdsinθk/c,…,expjN−12πf0Δdsinθk/cT∈ℂN×K, sm=s1m,s2m,…,sKmT∈ℂK×1, em=e1m,e2m,…,eNmT∈ℂN×1 is the additive complex white Gaussian noise vector with zero mean and variance σ2I. When M>1, equation (3) is the case of multiple measurement vector (MMV).Figure 1
The configuration of a ULA.Letθ˜=θ˜1,θ˜2,…,θ˜L be a fixed sampling grid in the DOA range -90°,90°, where L denotes the grid number and typically satisfies L≫N>K. Without loss of generality, let θ˜ be a uniform grid with a grid interval r=θ˜2−θ˜1∝L−1. In practical situations, the target locations are randomly distributed in space, which leads to the offgrid problem.Supposeθk∉θ˜k∈1,2,…,K and that θ˜l˜l˜∈1,2,…,L is the nearest grid to θk, which satisfies θk−θ˜l˜≤r/2. So, the steering vector aθk can be approximated by the first-order Taylor [23],(2)aθk≈aθ˜l˜+bθ˜l˜θk−θ˜l˜,where bθ˜l˜=a′θ˜l˜ is the first-order derivative of aθ˜l˜. Let Aθ˜=aθ˜1,aθ˜2,…,aθ˜L, Bθ˜=bθ˜1,bθ˜2,…,bθ˜L, and grid parameter βl=θk−θ˜l˜,l=l˜,k=1,2,…,K. By considering the approximation error into the noise, the offgrid DOA estimation sparse model can be expressed as(3)Y=ΦβX+E,where Φβ=Aθ˜+Bθ˜diagβ is abbreviated as Φ in the following paragraphs, Y∈ℂN×M, E∈ℂN×M, X=x1,x2,…,xM∈ℂL×M, Aθ˜, Bθ˜, and Φβ∈ℂN×L.In SBL, the sparse representation ofΦX is induced by designing a prior px for x. Instead of using the prior px directly, the SBL typically uses a two-layer hierarchical prior model, which involves a conditional prior px|α and a hyperprior pα with the hyperparameters α=α1,α2,…,αL, which controls the estimation accuracy and sparsity of x.


## 3. CSBL Offgrid DOA Estimation Algorithm

### 3.1. Two-Layer Hierarchical Prior Model

The extension of SBL from the real-valued model is considered to the complex-valued model, and the SBL form is unified, which can be applied to both real-valued and complex-valued models with different values of parameters. Assumed that the signals are independent at different snapshots, i.e.,xm,m=1,2,…,M are independent of each other, there exist the following conditional prior equations [24],(4)pX|α=∏m=1Mpxm|α=∏m=1M∏i=1Lpxim|αi.or(5)pX|γ=∏m=1Mpxm|γ=∏m=1M∏i=1Lpxim|γi,where pxi|αi=p/παipexp−pxi2/αi, and pxi|γi=pγi/πpexp−pγixi2. αi is the variance of xi, and γi=1/αi is the precision of xi. The parameter p=1 when xi is complex, and p=1/2 when xi is real. Next, the hyperprior pαi for variance αi and hyperprior pγi for precision γi are modeled, respectively.The generalized inverse Gaussian (GIG) distribution is chosen [25, 26] as the hyperprior pαi for variance αi with the following expression with three parameters λ, a, and b.(6)pα|λ,b,a=∏iGIGαi|λ,b,a=∏ia/bλ/22Kλabαiλ−1exp−12aαi+bαi,where Kλ⋅ denotes the modified Bessel function of the second kind of order λ, and the value range of parameters λ, a, and b are described later. The GIG distribution includes three distributions as special cases: the Gamma distribution, the inverse Gamma distribution, and the IG distribution according to the different values of the three parameters shown in Figure 2.Figure 2
The relationship between different two-layer hierarchical priors.For the hyperpriorpγi, from the inverse relationship between precision and variance γi=1/αi, the probability density function transformation formula [26] is used, and the hyperprior pγi is given by(7)pγ|λ,a,b=∏iGIGγi|−λ,a,b=∏ib/a−λ/22K−λabγi−λ−1exp−12bγi+aγi.In equations (6) and (7), the parameters satisfy the relationship as follows:(8)b≥0,a>0ifλ>0,b>0,a>0ifλ=0,b>0,a≥0ifλ<0.The marginal distribution ofx can be obtained by the two-layer hierarchical prior model, px|λ,b,a=∏ipxi|λ,b,a, where(9)pxi|λ,b,a=∫pxi|αipαi|λ,b,adαi=pπpa/bλ/2aλ/2−p/2Kλabb+2pxi2λ/2−p/2Kλ−pab+2pxi2.or(10)pxi|λ,a,b=∫pxi|γipγi|λ,a,bdαi=pπpb/a−λ/2aλ−p/2K−λabb+2pxi2λ/2−p/2K−λ−pab+2pxi2.The abovementioned setting of independent hyperparameters for each parameter to be estimated is the most significant feature of the SBL model, and it is also the fundamental reason for the sparsity of the model [27, 28]. When xi∈ℝ takes the value p=1/2, the equation (9) corresponds to the generalized hyperbolic (GH) distribution of the real-valued model, where equation (10) is equivalent to (9). In this paper, p=1 corresponds to the GH distribution of the complex-valued model.As mentioned above, the GIG distribution is selected as the hyperpriorpαi and pγi, respectively, as it includes a fairly broad class of distributions commonly used as hyperpriors, and the resulting marginal distribution, the GH distribution, again covers a large number of distributions as special cases. Due to this generalization, the connections between different modeling strategies can be analyzed, and several special cases that correspond to standard priors commonly used in sparse modeling are summarized in Figure 2.


### 3.2. Complex Sparse Bayesian Learning

As mentioned above, the model of the signal matrixX has been defined by different two-layer hierarchical priors, and the SBL model characterization will be completed in the following by modeling the observation matrix Y in equation (3). Firstly, according to the different values in the model, the noise E is modeled as Gaussian (real-valued model) or complex Gaussian (complex-valued model) with independent and identical distribution as follows:(11)pE|τ=∏m=1M∏n=1Npenm|τ=pπτ−1pNMexp−pτEF2,where τ=1/σ2 is the noise precision with variance σ2. The parameter p=1 when enm is complex, and p=1/2 when enm is real. Therefore, the following expression is obtained:(12)pY|X,β,τ=pπτ−1pNMexp−pτY−ΦXF2.The noise precisionτ is unknown, and a Gamma hyperprior is assumed for τ with pτ|c,d=Γτ|c,d, where τ obeys the Gamma distribution of parameters c and d. The offgrid parameter β obeys a uniform distribution with β∼U−r/2,r/2L. Taking the two-layer hierarchical prior CNxi|0,αi, Γγi|λ,a/2 as an example, the relationship between the parameters of the SBL model is shown in Figure 3.Figure 3
SBL off-grid DOA estimation parameter model.For equation (3) of the off-grid DOA estimation model, when the signal variance α is used to model the two-layer hierarchical prior, the joint probability distribution of the complete SBL model is represented as(13)pY,X,α,β,τ=pY|X,β,τpX|αpαpβpτ.According to Bayesian theory, the posterior distribution of the parameters is estimated as(14)pX,α,β,τ|Y=pY|X,β,τpX|αpαpβpτpY.However, the exact posteriorpX,α,β,τ|Y cannot be given in closed form due to the fact that pY is not calculated analytically [25]. Therefore, SBL needs an effective approximation pX,α,β,τ|Y. The sparsity-inducing property of the resulting estimator X does not only depend on the two-layer hierarchical prior, but also the approximation method used.In SBL, two widespread approximate approaches referred to as type I and type II estimation have been used. Type II is considered in this paper. In type II estimation [4–9], the impact of hyperparameter αi or γi is concerned on the model, and pX,α,β,τ|Y can be decomposed as(15)pX,α,β,τ|Y=pX|Y,α,β,τpα,β,τ|Y,where pα,β,τ|Y=pY|α,β,τpαpβpτ/pY. Therefore, type II estimation needs to effectively approximate the posterior distribution pα,β,τ|Y, by using MAP to estimate hyperparameter α as follows:(16)α^MAP=argmaxαpY|α,β,τlikelihoodpαprior=argmaxαpα∫pY|X,β,τpX|αdX=argmaxαlnpα+ln∫pY|X,β,τpX|αdX,and type II estimation is equivalent to maximizing the cost function of α, then(17)ℒα=−p×trYHΣY−1Y−pMlnΣY+lnpα,where ΣY is the covariance matrix of the likelihood function pY|α,β,τ. In the current widely studied real-valued model, it is usually assumed that a uniform hyperprior pα∝1 for α, which leads to maximizing α^MAP equivalently to maximizing the likelihood function pY|α,β,τ. Therefore, type II estimation is also called as type II maximum likelihood (ML) estimation.In [26], the results show that type II estimation has superior reconstruction performance than type I estimation. In subsequent paragraphs, type II estimation combined with the above two-layer hierarchical prior extension from the real-valued model is changed to the complex-valued model for the SBL derivation for offgrid DOA estimation.In type II estimation, equation (15) is used to decompose the exact posterior pX,α,β,τ|Y, and the following expression can be acquired:(18)pX|Y,α,β,τ=pY|X,β,τpX|αpY|α,β,τ=pY|X,β,τpX|α∫pY|X,β,τpX|αdX=pπpLMΣX−pMexp−p×trX−μXHΣX−1X−μX,where ΣX and μX are the covariance matrix and the mean of pX|Y,α,β,τ, respectively,(19)ΣX=τΦHΦ+Λα−1−1μX=τΣXΦHY.The signal varianceα and noise precision τ are estimated by using MAP as follows:(20)pα,β,τ|Y∝pY|α,β,τpαpβpτ=∫pY|X,β,τpX|αpαpβpτdX=pπpNMΣY−pMexp−p×trYHΣY−1Ypαpβpτ,where ΣY is the covariance matrix of pα,β,τ|Y, and ΣY=τ−1IN+ΦΛαΦH. If α and τ are known, type II estimated value of X is obtained as(21)X^=μX=ΛαΦHΣY−1Y,where Λα=diagα is a diagonal matrix composed of signal variance, and its diagonal elements control the row-sparsity of X^. As αi tends to a particularly small value (theoretically zero), the corresponding i-th row of X^ becomes 0T. DOA parameters θ˜k^ (the closest grid to the true DOA θ˜k, k=1,2,…,K) can be estimated from the position of nonzero rows in X^.When signal precisionγ is used to model two-layer hierarchical prior, similar expression results are acquired as(22)ΣX=τΦHΦ+Λγ−1,ΣY=τ−1IN+ΦΛγ−1ΦH,X^=μX=τΣxΦHYΛγ−1ΦHΣY−1Y.where Λγ=diagγ is a diagonal matrix composed of signal precision, as γi tends to a particularly large value (theoretically infinity), the corresponding i-th row of X^ becomes 0T.In the above, the complete SBL model is analyzed with signal modeling varianceα or precision γ separately. Next, the iterative process of α will be mainly discussed here, the iterative process of γ is similar and will not be repeated.


### 3.3. Two-Layer Hierarchical Priorpxi|αi, IGαi|−λ,b/2 Type II Estimation

Maximizing the cost function ofα is equivalent to maximizing the following function:(23)ℒ=−pMlnΣy−p×trYHΣy−1Y+∑i=1Lλ−1lnαi−b2αi+c−1lnτ−dτ.Then the partial derivative ofℒ with respect to αi is(24)∂ℒ∂αi=pMΣxi,iαi−2−pMαi+pμXi22αi−2+λ−1αi+b2αi−2.Let the partial derivative be zero, it will result in the hyperparameter with the following expression:(25)αi=pMΣxi,i+pμXi22+b/2pM−λ+1.In this paper, the Gamma hyperprior for the noise precisionτ is considered, the update of τ can be similarly obtained as follows:(26)τ=pM×N+c−1pM×trΦHΦΣx+pY−ΦμXF2+d.


### 3.4. Estimate of the Offgrid Parameterβ

As mentioned above, the updates of signal variance and noise precision by SBL mode are derived. In the following, the update of the important parameterβ is explained in the offgrid DOA estimation model. Considering that the only information β∼U−r/2,r/2L is bounded, the update of β cannot be obtained by MAP, so the maximum expectation algorithm is used to estimate the offgrid parameter β.The quantityX is treated as a hidden variable and maximize EpX|Y,α,β,τlnpY,X,α,β,τ, where EpX|Y,α,β,τ⋅ denotes an expectation with respect to the posterior of X as given in equation (18) using the current estimates of the hyperparameters, and is abbreviated as E⋅. For the offgrid parameter β, the irrelevant items are ignored, and the value of β is computed by(27)βnew=argmaxβElnpY|X,β,τpβ.According to equations (12) and (14), maximizing ElnpY|X,β,τpβ is equivalent to minimizing EY−ΦXF2=βTPβ+2vTβ+C, where C is a constant independent of β, P=ReBHB∗⊙M×ΣX+μXμXH is a positive semidefinite matrix, v=ReM×diagBHAΣX−Re∑m=1MdiagμX∗mBHym−AμXm. The optimization problem of β is given by(28)βnew=argminβ∈−r/2,r/2LβTPβ+2vTβ.which is a linear least-squares problem with boundary constraints. Let its partial derivative with respect to β be zero as follows:(29)∂∂ββTPβ+2vTβ=2Pβ+v=0.IfP is invertible, equation (29) has a unique solution, that is(30)β˜=−P−1v.Considering the boundedness, whenβ˜∈−r/2,r/2L, then βnew=β˜. While one of the above two conditions is not satisfied, β is updated element by element. When updating βl, other elements are fixed and the l-th equation with the partial derivative is zero, and the expression is given by(31)β^l=−vp+Pl−lTβ−lPl,l=1,2,…,L,where ⋅−l is the vector to remove the l-th element and Pl is the l-th column in Pl. The calculation formula of β is described as(32)βlnew=β˜lifP−1existsandβ˜∈−r/2,r/2L1β^lif1doesn'tholdandβ^l∈−r/2,r/22−r/2if1and2don'tholdandβ^l<−r/23r/2if1and2don'tholdandβ^l>r/24.Finally, theK DOAs will be estimated in the following equation:(33)θ^k=θ˜k^+βk^,k=1,2,…,K,where θ˜k^ is the closest grid to the true DOA θ˜k, and βk^ is the offgrid parameter corresponding to θ˜k.


### 3.5. Proposed Algorithm Flow

Through the above derivation, the algorithm flow is concluded in Table1 for complex sparse Bayesian learning for offgrid DOA estimation, where the two-layer hierarchical prior pxi|αi=CNxi|0,αi, IGαi|−λ,b/2 as an example.Table 1
Complex sparse Bayesian learning algorithm for offgrid DOA estimation.
Input quantities: the received signals Y, and the sampling grid in the DOA range θ˜Output quantities: offgrid θInitialization: giving the initial values of parameters λ,b,c,d,α,τ,β←0, and setting the initial values of the convergence threshold ε and the maximum number of iterations max_iter, computing the array steering matrices Aθ˜ and Bθ˜.(1) Matrix construction: Computing Φβ=A+Bdiagβ from Aθ˜, Bθ˜, and β.(2) Parameters update:1)Updating matricesΣX and μX based on equation (19) from the current values of α,τ,Φ.2)Acquiring the estimated valuesθ˜k^ (the closest grid to the true DOA θ˜k, k=1,2,…,K) from the position of nonzero rows in X^=μX.3)Updating quantitiesα,τ from the current values of ΣX and μX with equations (25) and (26).4)Updating the offgrid parameterβ from the current values Aθ˜, Bθ˜, ΣX, and μX according to equation (32).(3) Iteration termination: Calculating the residual value Y−ΦβμXF2. If the residual value is less than ε or the number of iterations is more than max_iter, terminating the iteration process, otherwise, jumping to (1) to start the iteration, and updating the matrices Aθ˜, Bθ˜, and Φβ based on equation (33).(4) DOA detection:K DOAs will be estimated in equation (33).In this paper, for the two-layer hierarchical priorpxi|γi=CNxi|0,γi−1, Γγi|−λ,b/2 and pxi|αi=CNxi|0,αi, IGαi|−λ,b/2, there are λ=−1,b⟶0+. In the models of the two-layer hierarchical prior pxi|αi=CNxi|0,αi, Γαi|λ,a/2 and pxi|γi=CNxi|0,γi−1, IGγi|λ,a/2, with the complex-valued DOA estimation model, in order to obtain the L1-norm of xi, there are λ=3/2, a=1. The pτ=Γτ|c,d for the noise precision, there are c=1, d⟶0+.


## 3.1. Two-Layer Hierarchical Prior Model

The extension of SBL from the real-valued model is considered to the complex-valued model, and the SBL form is unified, which can be applied to both real-valued and complex-valued models with different values of parameters. Assumed that the signals are independent at different snapshots, i.e.,xm,m=1,2,…,M are independent of each other, there exist the following conditional prior equations [24],(4)pX|α=∏m=1Mpxm|α=∏m=1M∏i=1Lpxim|αi.or(5)pX|γ=∏m=1Mpxm|γ=∏m=1M∏i=1Lpxim|γi,where pxi|αi=p/παipexp−pxi2/αi, and pxi|γi=pγi/πpexp−pγixi2. αi is the variance of xi, and γi=1/αi is the precision of xi. The parameter p=1 when xi is complex, and p=1/2 when xi is real. Next, the hyperprior pαi for variance αi and hyperprior pγi for precision γi are modeled, respectively.The generalized inverse Gaussian (GIG) distribution is chosen [25, 26] as the hyperprior pαi for variance αi with the following expression with three parameters λ, a, and b.(6)pα|λ,b,a=∏iGIGαi|λ,b,a=∏ia/bλ/22Kλabαiλ−1exp−12aαi+bαi,where Kλ⋅ denotes the modified Bessel function of the second kind of order λ, and the value range of parameters λ, a, and b are described later. The GIG distribution includes three distributions as special cases: the Gamma distribution, the inverse Gamma distribution, and the IG distribution according to the different values of the three parameters shown in Figure 2.Figure 2
The relationship between different two-layer hierarchical priors.For the hyperpriorpγi, from the inverse relationship between precision and variance γi=1/αi, the probability density function transformation formula [26] is used, and the hyperprior pγi is given by(7)pγ|λ,a,b=∏iGIGγi|−λ,a,b=∏ib/a−λ/22K−λabγi−λ−1exp−12bγi+aγi.In equations (6) and (7), the parameters satisfy the relationship as follows:(8)b≥0,a>0ifλ>0,b>0,a>0ifλ=0,b>0,a≥0ifλ<0.The marginal distribution ofx can be obtained by the two-layer hierarchical prior model, px|λ,b,a=∏ipxi|λ,b,a, where(9)pxi|λ,b,a=∫pxi|αipαi|λ,b,adαi=pπpa/bλ/2aλ/2−p/2Kλabb+2pxi2λ/2−p/2Kλ−pab+2pxi2.or(10)pxi|λ,a,b=∫pxi|γipγi|λ,a,bdαi=pπpb/a−λ/2aλ−p/2K−λabb+2pxi2λ/2−p/2K−λ−pab+2pxi2.The abovementioned setting of independent hyperparameters for each parameter to be estimated is the most significant feature of the SBL model, and it is also the fundamental reason for the sparsity of the model [27, 28]. When xi∈ℝ takes the value p=1/2, the equation (9) corresponds to the generalized hyperbolic (GH) distribution of the real-valued model, where equation (10) is equivalent to (9). In this paper, p=1 corresponds to the GH distribution of the complex-valued model.As mentioned above, the GIG distribution is selected as the hyperpriorpαi and pγi, respectively, as it includes a fairly broad class of distributions commonly used as hyperpriors, and the resulting marginal distribution, the GH distribution, again covers a large number of distributions as special cases. Due to this generalization, the connections between different modeling strategies can be analyzed, and several special cases that correspond to standard priors commonly used in sparse modeling are summarized in Figure 2.


## 3.2. Complex Sparse Bayesian Learning

As mentioned above, the model of the signal matrixX has been defined by different two-layer hierarchical priors, and the SBL model characterization will be completed in the following by modeling the observation matrix Y in equation (3). Firstly, according to the different values in the model, the noise E is modeled as Gaussian (real-valued model) or complex Gaussian (complex-valued model) with independent and identical distribution as follows:(11)pE|τ=∏m=1M∏n=1Npenm|τ=pπτ−1pNMexp−pτEF2,where τ=1/σ2 is the noise precision with variance σ2. The parameter p=1 when enm is complex, and p=1/2 when enm is real. Therefore, the following expression is obtained:(12)pY|X,β,τ=pπτ−1pNMexp−pτY−ΦXF2.The noise precisionτ is unknown, and a Gamma hyperprior is assumed for τ with pτ|c,d=Γτ|c,d, where τ obeys the Gamma distribution of parameters c and d. The offgrid parameter β obeys a uniform distribution with β∼U−r/2,r/2L. Taking the two-layer hierarchical prior CNxi|0,αi, Γγi|λ,a/2 as an example, the relationship between the parameters of the SBL model is shown in Figure 3.Figure 3
SBL off-grid DOA estimation parameter model.For equation (3) of the off-grid DOA estimation model, when the signal variance α is used to model the two-layer hierarchical prior, the joint probability distribution of the complete SBL model is represented as(13)pY,X,α,β,τ=pY|X,β,τpX|αpαpβpτ.According to Bayesian theory, the posterior distribution of the parameters is estimated as(14)pX,α,β,τ|Y=pY|X,β,τpX|αpαpβpτpY.However, the exact posteriorpX,α,β,τ|Y cannot be given in closed form due to the fact that pY is not calculated analytically [25]. Therefore, SBL needs an effective approximation pX,α,β,τ|Y. The sparsity-inducing property of the resulting estimator X does not only depend on the two-layer hierarchical prior, but also the approximation method used.In SBL, two widespread approximate approaches referred to as type I and type II estimation have been used. Type II is considered in this paper. In type II estimation [4–9], the impact of hyperparameter αi or γi is concerned on the model, and pX,α,β,τ|Y can be decomposed as(15)pX,α,β,τ|Y=pX|Y,α,β,τpα,β,τ|Y,where pα,β,τ|Y=pY|α,β,τpαpβpτ/pY. Therefore, type II estimation needs to effectively approximate the posterior distribution pα,β,τ|Y, by using MAP to estimate hyperparameter α as follows:(16)α^MAP=argmaxαpY|α,β,τlikelihoodpαprior=argmaxαpα∫pY|X,β,τpX|αdX=argmaxαlnpα+ln∫pY|X,β,τpX|αdX,and type II estimation is equivalent to maximizing the cost function of α, then(17)ℒα=−p×trYHΣY−1Y−pMlnΣY+lnpα,where ΣY is the covariance matrix of the likelihood function pY|α,β,τ. In the current widely studied real-valued model, it is usually assumed that a uniform hyperprior pα∝1 for α, which leads to maximizing α^MAP equivalently to maximizing the likelihood function pY|α,β,τ. Therefore, type II estimation is also called as type II maximum likelihood (ML) estimation.In [26], the results show that type II estimation has superior reconstruction performance than type I estimation. In subsequent paragraphs, type II estimation combined with the above two-layer hierarchical prior extension from the real-valued model is changed to the complex-valued model for the SBL derivation for offgrid DOA estimation.In type II estimation, equation (15) is used to decompose the exact posterior pX,α,β,τ|Y, and the following expression can be acquired:(18)pX|Y,α,β,τ=pY|X,β,τpX|αpY|α,β,τ=pY|X,β,τpX|α∫pY|X,β,τpX|αdX=pπpLMΣX−pMexp−p×trX−μXHΣX−1X−μX,where ΣX and μX are the covariance matrix and the mean of pX|Y,α,β,τ, respectively,(19)ΣX=τΦHΦ+Λα−1−1μX=τΣXΦHY.The signal varianceα and noise precision τ are estimated by using MAP as follows:(20)pα,β,τ|Y∝pY|α,β,τpαpβpτ=∫pY|X,β,τpX|αpαpβpτdX=pπpNMΣY−pMexp−p×trYHΣY−1Ypαpβpτ,where ΣY is the covariance matrix of pα,β,τ|Y, and ΣY=τ−1IN+ΦΛαΦH. If α and τ are known, type II estimated value of X is obtained as(21)X^=μX=ΛαΦHΣY−1Y,where Λα=diagα is a diagonal matrix composed of signal variance, and its diagonal elements control the row-sparsity of X^. As αi tends to a particularly small value (theoretically zero), the corresponding i-th row of X^ becomes 0T. DOA parameters θ˜k^ (the closest grid to the true DOA θ˜k, k=1,2,…,K) can be estimated from the position of nonzero rows in X^.When signal precisionγ is used to model two-layer hierarchical prior, similar expression results are acquired as(22)ΣX=τΦHΦ+Λγ−1,ΣY=τ−1IN+ΦΛγ−1ΦH,X^=μX=τΣxΦHYΛγ−1ΦHΣY−1Y.where Λγ=diagγ is a diagonal matrix composed of signal precision, as γi tends to a particularly large value (theoretically infinity), the corresponding i-th row of X^ becomes 0T.In the above, the complete SBL model is analyzed with signal modeling varianceα or precision γ separately. Next, the iterative process of α will be mainly discussed here, the iterative process of γ is similar and will not be repeated.


## 3.3. Two-Layer Hierarchical Priorpxi|αi, IGαi|−λ,b/2 Type II Estimation

Maximizing the cost function ofα is equivalent to maximizing the following function:(23)ℒ=−pMlnΣy−p×trYHΣy−1Y+∑i=1Lλ−1lnαi−b2αi+c−1lnτ−dτ.Then the partial derivative ofℒ with respect to αi is(24)∂ℒ∂αi=pMΣxi,iαi−2−pMαi+pμXi22αi−2+λ−1αi+b2αi−2.Let the partial derivative be zero, it will result in the hyperparameter with the following expression:(25)αi=pMΣxi,i+pμXi22+b/2pM−λ+1.In this paper, the Gamma hyperprior for the noise precisionτ is considered, the update of τ can be similarly obtained as follows:(26)τ=pM×N+c−1pM×trΦHΦΣx+pY−ΦμXF2+d.


## 3.4. Estimate of the Offgrid Parameterβ

As mentioned above, the updates of signal variance and noise precision by SBL mode are derived. In the following, the update of the important parameterβ is explained in the offgrid DOA estimation model. Considering that the only information β∼U−r/2,r/2L is bounded, the update of β cannot be obtained by MAP, so the maximum expectation algorithm is used to estimate the offgrid parameter β.The quantityX is treated as a hidden variable and maximize EpX|Y,α,β,τlnpY,X,α,β,τ, where EpX|Y,α,β,τ⋅ denotes an expectation with respect to the posterior of X as given in equation (18) using the current estimates of the hyperparameters, and is abbreviated as E⋅. For the offgrid parameter β, the irrelevant items are ignored, and the value of β is computed by(27)βnew=argmaxβElnpY|X,β,τpβ.According to equations (12) and (14), maximizing ElnpY|X,β,τpβ is equivalent to minimizing EY−ΦXF2=βTPβ+2vTβ+C, where C is a constant independent of β, P=ReBHB∗⊙M×ΣX+μXμXH is a positive semidefinite matrix, v=ReM×diagBHAΣX−Re∑m=1MdiagμX∗mBHym−AμXm. The optimization problem of β is given by(28)βnew=argminβ∈−r/2,r/2LβTPβ+2vTβ.which is a linear least-squares problem with boundary constraints. Let its partial derivative with respect to β be zero as follows:(29)∂∂ββTPβ+2vTβ=2Pβ+v=0.IfP is invertible, equation (29) has a unique solution, that is(30)β˜=−P−1v.Considering the boundedness, whenβ˜∈−r/2,r/2L, then βnew=β˜. While one of the above two conditions is not satisfied, β is updated element by element. When updating βl, other elements are fixed and the l-th equation with the partial derivative is zero, and the expression is given by(31)β^l=−vp+Pl−lTβ−lPl,l=1,2,…,L,where ⋅−l is the vector to remove the l-th element and Pl is the l-th column in Pl. The calculation formula of β is described as(32)βlnew=β˜lifP−1existsandβ˜∈−r/2,r/2L1β^lif1doesn'tholdandβ^l∈−r/2,r/22−r/2if1and2don'tholdandβ^l<−r/23r/2if1and2don'tholdandβ^l>r/24.Finally, theK DOAs will be estimated in the following equation:(33)θ^k=θ˜k^+βk^,k=1,2,…,K,where θ˜k^ is the closest grid to the true DOA θ˜k, and βk^ is the offgrid parameter corresponding to θ˜k.


## 3.5. Proposed Algorithm Flow

Through the above derivation, the algorithm flow is concluded in Table1 for complex sparse Bayesian learning for offgrid DOA estimation, where the two-layer hierarchical prior pxi|αi=CNxi|0,αi, IGαi|−λ,b/2 as an example.Table 1
Complex sparse Bayesian learning algorithm for offgrid DOA estimation.
Input quantities: the received signals Y, and the sampling grid in the DOA range θ˜Output quantities: offgrid θInitialization: giving the initial values of parameters λ,b,c,d,α,τ,β←0, and setting the initial values of the convergence threshold ε and the maximum number of iterations max_iter, computing the array steering matrices Aθ˜ and Bθ˜.(1) Matrix construction: Computing Φβ=A+Bdiagβ from Aθ˜, Bθ˜, and β.(2) Parameters update:1)Updating matricesΣX and μX based on equation (19) from the current values of α,τ,Φ.2)Acquiring the estimated valuesθ˜k^ (the closest grid to the true DOA θ˜k, k=1,2,…,K) from the position of nonzero rows in X^=μX.3)Updating quantitiesα,τ from the current values of ΣX and μX with equations (25) and (26).4)Updating the offgrid parameterβ from the current values Aθ˜, Bθ˜, ΣX, and μX according to equation (32).(3) Iteration termination: Calculating the residual value Y−ΦβμXF2. If the residual value is less than ε or the number of iterations is more than max_iter, terminating the iteration process, otherwise, jumping to (1) to start the iteration, and updating the matrices Aθ˜, Bθ˜, and Φβ based on equation (33).(4) DOA detection:K DOAs will be estimated in equation (33).In this paper, for the two-layer hierarchical priorpxi|γi=CNxi|0,γi−1, Γγi|−λ,b/2 and pxi|αi=CNxi|0,αi, IGαi|−λ,b/2, there are λ=−1,b⟶0+. In the models of the two-layer hierarchical prior pxi|αi=CNxi|0,αi, Γαi|λ,a/2 and pxi|γi=CNxi|0,γi−1, IGγi|λ,a/2, with the complex-valued DOA estimation model, in order to obtain the L1-norm of xi, there are λ=3/2, a=1. The pτ=Γτ|c,d for the noise precision, there are c=1, d⟶0+.


## 4. Simulation and Performance Analysis

This section validates the validity of the different two-level hierarchical priors for the application of DOA estimation. Suppose that ULA withN=21 array elements for receiving far-field narrowband sources, the sources carrier frequency is f0=10GHz with SNR=0dB and interelement spacing Δd=3×108/2f0=1.5cm in equation (1).Parameter initialization is as follows: there areαl=1,l=1,2,…L (or γl=1) and σ2=10−SNR/10YoF2/NM, where Yo is the noisefree received signal matrix. The offgrid parameter is β=0. A uniform sampling grid is θ˜=θ˜1,θ˜2,…,θ˜L with a grid interval r=1° in the DOA range -90°,90°. The number of snapshots is M=200. The proposed algorithm terminates if ε≤10−4 or the maximum number of iterations reaches max_iter = 2000.The DOA estimation results and details of the estimation results of a source located at DOA0.5° are shown in Figures 4(a) and 4(b), respectively.Figure 4
Single source of DOA estimation results: (a) single source DOA estimate and (b) single source DOA estimate details.
(a)(b)The DOA estimation results and details of the estimation results of the two sources located at DOAs−2.4° and 3.6° are shown in Figures 5(a) and 5(b), respectively.Figure 5
Two sources of DOA estimation results: (a) two sources DOA estimates and (b) two sources DOA estimates details.
(a)(b)Comparing Figures4 and 5, it can be seen that under the consideration of the offgrid DOA model, different two-layer hierarchical priors used for sparse source xi modeling can effectively estimate the DOA of the offgrid model, while the traditional MUSIC and L1-SVD algorithms have the DOA estimation results located on the spatial sampling grid nearest to the target point because the sampling grid mismatch is not considered.In order to analyze the different two-layer hierarchical priors used for sparse sourcexi modeling, the performance of the different two-layer hierarchical priors will be compared with the traditional MUSIC algorithm and L1-SVD algorithm, and the performances between the different two-level hierarchical priors are illustrated by simulation results. The root mean square error (RMSE) of the DOA estimation result is given in the form of dB as follows:(34)RMSE=10log101KT∑k=1K∑t=1Tθ^kt−θk2,where K is the total number of sources (K=3), T is the number of experiments (T=200), and θ^kt is the estimated DOA of the k-th source in the t-th experiment. Then consider 3 targets with DOAs θ1=−2.4°, θ2=0.5°, and θ3=3.6° respectively.Figure6(a) shows the RMSE of DOAs estimation result, with the variation of SNR under the condition that the number of array elements is N=21 and the number of snapshots M=200. Figure 6(b) shows the RMSE of DOAs estimation result, with the variation of array elements under the condition of SNR=0dB and the number of snapshots M=200. Figure 6(c) shows the RMSE of DOAs estimation result, with the variation of snapshots under the condition of SNR=0dB and the number of array elements N=21. The Cramer–Rao bound (CRB) for the offgrid bias parameters of the offgrid model in (2) and (3) is acquired by CRB=1/2MtrBHBSNR [29–33]. Due to the fact that the value of trBHB tends to 104, the values of CRB in Figure 6 is about 10−3∼10−4, such as CRB = [0.00107, 0.00085, 0.00067, 0.0005, 0.0004, 0.0003] in Figure 6(a), which causes that the CRB waveforms appear relatively flat by compared with RMSE value waveforms.Figure 6
RMSE versus parameters: (a) RMSE with different SNR values; (b) RMSE with different array elements; and (c) RMSE with different snapshots.
(a)(b)(c)According to the comparison results of RMSE under different conditions in Figure6, it can be seen that the RMSE of traditional MUSIC and L1-SVD both show a straight line when the SNR, the number of array elements, and the number of snapshots reach a certain level. Because they cannot estimate the offgrid parameter, and the maximum estimation precision of MUSIC and L1-SVD algorithm is −2°, 0° or 1°, 4°, respectively. In addition, the different two-layer hierarchical priors of the signal xi are considered, and their performance comparisons are shown in Figure 6. The precision of DOA estimation has been significantly improved due to the fact that the offgrid DOA model takes into account the grid offset.At the same time, the results of Figure6 are further analyzed as following, for the Gamma hyperprior of pαi and pγi as an example. The marginal prior from two-layer hierarchical prior pxi|γi and Γγi|−λ,b/2 is student’s t prior, and the marginal prior obtained by setting the values of λ for the two-layer hierarchical prior pxi|αi, Γαi|λ,a/2 is the Laplace prior. In [12], it is analyzed in detail that the Laplace prior has better performance than the student’s t prior under the condition of SMV of the real-valued model. In this paper, this conclusion is extended to the case of MMV of the complex-valued model, the effects of different hyperpriors are considered, and it is found that the two-layer hierarchical prior pxi|αi, IGαi|−λ,b/2 has better performance than pxi|γi, IGγi|λ,a/2, although the marginal prior obtained by the two-layer hierarchical prior pxi|γi, Γγi|−λ,b/2 and pxi|αi, IGαi|−λ,b/2 are the same as student’s t prior, and the marginal prior obtained by the two-layer hierarchical prior pxi|αi, Γαi|λ,a/2, and pxi|γi, IGγi|λ,a/2 are the same as Laplace prior.Through the above analysis, the performance of the two-layer hierarchical priorpxi|αi and pαi is better than pxi|αi and pγi, and the IG prior model is better than Gamma prior model.The proposed algorithm is further applied to the angle-distance positioning in the FDA (frequency diverse array) [34]. The complex multisnapshot SBL combined with FDA radar characteristics is used for target angle-distance two-dimensional localization. The log frequency offset is used at the transmitter to remove the distance periodicity of the beam pattern of the full-band frequency receiving mechanism of FDA radar for distance unambiguous localization. Through the transmitting of zero frequency offset and nonzero frequency signals, the angle-distance two-dimensional estimation is converted into two one-dimensional estimates of the azimuth and distance, respectively, thus reducing the amount of calculation. The simulation results are shown in Figure 7 by using the two-layer hierarchical prior pxi|αi, Γαi|λ,a/2. The positions of the three targets are −2.4°,9576m, 0.5°,9492m, and 3.6°,9533m, respectively. The simulation parameters are consistent with Figure 5. The log frequency offset Δfn=Δflogn,n=1,2,…,N with Δf=2MHz. The grid interval r=10m in the distance interval 9000m,10000m. The simulation results are −2.402°,9576.336m, 0.505°,9492.047m, and 3.583°,9533.193m.Figure 7
Three sources of angle-distance estimation results: (a) three sources azimuth estimates and (b) three sources distance estimates.
(a)(b)


## 5. Conclusions

A general complex signal prior construction suitable for complex-valued signal modeling is proposed. By using the complex Gaussian scale mixtures hierarchy, it has been shown that this signal model includes complex versions of a number of signal models commonly used for sparse signal modeling. In addition, the expression of the SBL form of the real-valued and complex-valued models is unified by parameter values. Our motivation is to generalize the current research work of SBL by this unified form and contribute to its generalization to complex-valued signal models. Finally, the SBL complex-valued form is applied to the off-grid DOA estimation complex-valued model and compared with the performance between different sparsity-inducing priors. Theoretical analysis and simulation results show that the proposed algorithm can effectively use SBL to process complex-valued signal models and has lower algorithm complexity.



---

*Source: 4500243-2022-02-28.xml*