# Direction-of-Arrival Estimation for Coherent Sources via Sparse Bayesian Learning

Authors: Zhang-Meng Liu; Zheng Liu; Dao-Wang Feng; Zhi-Tao Huang
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

A spatial filtering-based relevance vector machine (RVM) is proposed in this paper to separate coherent sources and estimate their directions-of-arrival (DOA), with the filter parameters and DOA estimates initialized and refined via sparse Bayesian learning. The RVM is used to exploit the spatial sparsity of the incident signals and gain improved adaptability to much demanding scenarios, such as low signal-to-noise ratio (SNR), limited snapshots, and spatially adjacent sources, and the spatial filters are introduced to enhance global convergence of the original RVM in the case of coherent sources. The proposed method adapts to arbitrary array geometry, and simulation results show that it surpasses the existing methods in DOA estimation performance.


---

## Body

## 1. Introduction

The problem of direction-of-arrival (DOA) estimation with sensor arrays is common in various applications, such as radar, sonar, and wireless communications [1–4], and the requirement for finding the directions of coherent sources also widely emerges due to the factors like multipath and jamming. Most of the existing superresolution methods, for example, the subspace-based ones [1], do not adapt well to such a problem in their original form, and further research has been carried out to improve their adaptability in the coherent scenarios [2, 5–7]. The method in [2] aims only at CDMA signals with known codes, and those in [5–7] are general-purpose ones that apply to common narrowband signals. Pillai and Kwon proposed making up for the deficient rank in the array output covariance matrix via spatial smoothing [5], whose performance may deteriorate if the sources are spatially adjacent. Then Moghaddamjoo introduced the idea of spatial filtering to better decorrelate the signals [6]. After that, Delis and Papadopoulos promoted the filtering idea further by combining it with the forward and backward processing [7]. The spatial smoothing and filtering techniques require shift invariance property in the array; thus, they adapt to special array geometries such as the uniform linear arrays (ULA) only [5–7]. Moreover, after spatial preprocessing, the above methods exploit ordinary subspace-based techniques to realize DOA estimation; thus, they all lack adaptability to demanding scenarios, such as low signal-to-noise ratio (SNR), limited snapshots, and spatially adjacent sources, as good performances of the techniques rely heavily on high-precision covariance matrix and subspace estimates.The recently emerging technique of sparse reconstruction has attracted much attention from various areas, and it also provides a new perspective for DOA estimation. Existing sparsity-inducing DOA estimation methods include FOCUSS [8], L1-SVD [9], JLZA-DOA [10], CMSR [11], RVM-DOA [12–14], and so forth. Those methods combine the spatial sparsity of the incident signals when fitting the array output with an overcomplete model and show improved adaptation to the above mentioned demanding scenarios [8–14]. As they realize DOA estimation by reconstructing the array output, instead of exploiting the covariance matrix, some of them also perform well for correlated and coherent sources [8–10]. However, the way for different methods in describing the spatial sparsity of the incident signals diverges largely. The L0-norm cost function is canonical as it stands for the model dimension exactly, but it introduces too many local maxima in the objective function and makes the global maximum hardly attainable. Therefore, Lp-norm (0<p≤1) cost functions are used instead to mitigate the influence of the local maxima at the expense of possibly biased global maximum [8, 9, 12–15]. Another technique for combining model sparsity with data fitting is the relevance vector machine (RVM) [16], and conclusions have been obtained to indicate that the RVM has identical global maximum as the L0-norm based function and greatly decreased local maxima if the signals are independent [15]. It is even unimodal at high SNR when the basis amplitudes are highly scaled [15]. Based on those advantages of the RVM, our method of RVM-DOA surpasses the Lp-norm based ones in DOA estimation precision for independent sources [12]. Nonetheless, when the incident signals are correlated, the superiority of the RVM is no longer guaranteed [17], and unbiased DOA estimates may become unattainable.In this paper, we extend RVM-DOA [12], which was originally designed for independent sources, to coherent source scenarios by separating the sources with spatial filters. After such separation, the original multiple-coherent-source data is decomposed into more than one single-source components, thus improved global convergence and DOA estimation precision are expected for the RVM-based methods. The filter parameters and DOA estimates are updated via sparse Bayesian learning after reasonable initialization with the independent RVM-DOA method. The new method adapts to arbitrary array geometries, which greatly relaxes the restrictions of previous methods on the shift invariance of the array that has constrained their adaptation to special arrays like ULA only. Empirical results show that the new method obtains enhanced applicability in challenging environments and improved DOA estimation precision under moderate settings.The rest of the paper mainly consists of four parts. Section2 reviews the coherent signal model. Section 3 introduces the new method in detail. Section 4 demonstrates the performance of the proposed method via simulations. Section 5 concludes the whole paper.


## 2. Model Formulation

Suppose thatK coherent signals impinge onto an M-element array from the directions of ϑ=[ϑ1,…,ϑK] simultaneously; the array output at time t is
(1)x(t)=A(ϑ)s(t)+v(t),
where A(ϑ)=[a(ϑ1),…,a(ϑK)], a(ϑk)=[ejφk,1,…,ejφk,M]T is the array responding vector to the kth incident signal, φk,m is the phase shift of this signal when propagating from the reference point to the mth sensor, s(t)=[s1(t),…,sK(t)]T is the source waveform vector, and v(t) is the zero-mean white Gaussian noise with known unified covariance matrix Q and unknown power σ2. In order to simplify notation, we omit the time identifier (t) in (1) on the variables in the case of single snapshot when no confusion emerges and use a(·) and A(·) to represent the mappings from a direction and direction set to the corresponding array responding vector and matrix, respectively.Due to the coherence between the incident signals, the source covariance matrixRs=E{ssH} is singular, and the signal subspace of the array output covariance matrix RX=E{xxH}=A(ϑ)RsAH(ϑ)+σ2Q is rank deficient. Thus preprocessing techniques, such as spatial smoothing and filtering, are required to make up for the deficient rank by exploiting the shift invariance property of special arrays [5–7]. Although those techniques help to enhance the adaptability of the ordinary subspace-based methods to coherent signal scenarios to some extent, they perform poorly when the covariance matrix and subspace estimates are greatly inaccurate in cases like low SNR, limited snapshots, and spatially adjacent sources.Some of the sparsity-inducing DOA estimation methods apply to coherent sources directly [8–10]. However, in such scenarios, the N snapshots can be rewritten as X=A(ϑ)(αs0)+V, where X=[x(t1),…,x(tN)], V=[v(t1),…,v(tN)], s0=[s0(t1),…,s0(tN)] is the waveform of the reference signal, and α=[α1,…,αK]T contains the amplitudes and phases of the incident signals with respect to the reference one. This form of the array output implies that, when the sources are coherent, more snapshots only help increase the SNR but do not enhance global convergence of the sparsity-inducing methods by smoothing away some of the local maxima [15]. Therefore, although those methods still work in the coherent scenarios, they may not obtain unbiased and high-precision DOA estimates.


## 3. Spatial Filtering-Based Relevance Vector Machine (RVM) for Coherent DOA Estimation

In this section, we first review the overcomplete model for the sparsity-inducing DOA estimation methods and analyze how the previous independent RVM-based method performs in the case of coherent sources. Then in Section3.2, we introduce the basic idea of the spatial filters to be used to separate the coherent signal components. After that, we discuss in detail how the spatial filter parameters are updated and the DOA estimates are refined iteratively via Bayesian learning. Finally, we outline the proposed method by joining the main procedures together and complementing them with reasonable initializations in Section 3.4.


### 3.1. Overcomplete Model and Problem Formulation

In the sparse reconstruction-based DOA estimation methods, the potential space of the incident signals is sampled discretely to obtain a direction setΘ with cardinality of card(Θ)=L to contain all the possible signal directions. The corresponding array responding matrix is denoted by A-=A(Θ)∈ℂM×L. Generally, L is much larger than M to meet the requirement in superresolution; thus, A- is overcomplete. The array output given by (1) can then be rewritten as
(2)x=A-s-+v,
where s- is the extension of s from ϑ to Θ and has nonzero amplitudes only for the ϑ-associated bases. The model mismatch between (1) and (2) due to the discretization of Θ is skipped over to facilitate notation unless otherwise stated; such approximation does not influence the implementation of those methods. The signal directions can be estimated according to the locations of the nonzero magnitudes in s- if it can be estimated. However, the equation given by (2) is ill posed and it is impossible to estimate s- unless further sparsity constraints on s-, which indicates that only a few signals impinge simultaneously, are added. When multiple snapshots are collected, one needs to estimate S-=[s-(t1),…,s-(tN)] from X=[x(t1),…,x(tN)].In RVM, the columns ofS- are assumed to be independent and Gaussian distributed; that is, s-~𝒩(0,Γ), with Γ=diag⁡{γ} and γ=[γ1,…,γL]. Such an assumption leads to the following probability function of X with respect to γ and σ2:
(3)p(X∣γ,σ2)=∫p(X∣S-,σ2)p(S-∣γ)dS-=|πΣX|-Nexp⁡(-tr⁡(XHΣX-1X)),
where ΣX=A-ΓA-H+σ2Q. Its logarithmic form is used to estimate γ and σ2; that is,
(4)ℒ(γ,σ2)=ln⁡p(X∣γ,σ2)≡-Nln⁡|ΣX|-tr⁡(XHΣX-1X).
The function of ℒ(γ,σ2) can also be rewritten as ℒ(Γ,σ2) if intersource correlation exists and the diagonal structure of Γ is deviated.Then the expectation-maximization (EM) algorithm [18], which has already been used for wideband DOA estimation [19, 20], can be used to maximize ℒ(γ,σ2) and estimate the source directions [12]. We have shown in [12] that the RVM-based method, if well designed, surpasses its counterparts with significant margin in DOA estimation performance when the incident signals are independent. However, how it will behave in the coherent source settings still remains unknown. Fortunately, some conclusions available in the related neural source localization area have provided some clues about the behavior of the independent RVM when intersource correlation exists [17].Denotecpq=〈s-p•,s-q•〉/N, with s-p• standing for the pth row of S-; then, according to Theorem  1 in [17], the following propositions can be obtained to indicate the behavior of RVM when applied to independent and correlated DOA estimation.Proposition 1.
Suppose that the SNR is high and the number of the incident sources is smaller than that of the sensors, thenℒ(γ,σ2) and ℒ(Γ,σ2) own the following properties.(1)
Ifcpq=0for each two sources, thenℒ(γ,σ2)is unimodal and no local maxima exists.(2)
For arbitrarycpq, the unique global maximum ofℒ(Γ,σ2), denoted byΓ*(whileσ2is negligible), produces the true source matrixS-*=E(S-;X,Γ*); that is, it achieves perfect source reconstruction.Proposition1 is a straightforward extension of Theorem  1 in [17]. Property (1) of Proposition 1 implies that well-designed RVM is able to achieve consistent DOA estimation for independent sources, as our RVM-DOA method has done [12]. Property (2) tells that accurate source covariance matrix estimation, which is associated with accurate DOA estimates, is necessary for precise source matrix reconstruction in the case of coherent sources. However, as some of the off-diagonal elements of the coherent source covariance matrix have significant nonzero magnitudes, the original RVM should be modified to take all possible intersource correlations into consideration, so as to make global convergence achievable. Unfortunately, all the L2 elements of Γ should be estimated under such a modified model, and the great increase in the number of unknown parameters not only aggravates the computational burden of the RVM prohibitively, but also increases the risk of overfitting [21].In addition, it should be noted that, although the independent RVM-based DOA estimation method is not able to perform well enough in the case of coherent sources, it still converges to the vicinity of the global maximum [22]. Therefore, we use the independent RVM-DOA method in [12] to obtain coarse DOA estimates in this paper.


### 3.2. Basic Ideas for Designing the Spatial Filters

According to Proposition1 and the analysis that follows, the global maximum of ℒ(Γ,σ2) can hardly be achieved by RVM in the case of coherent sources, and the corresponding DOA estimation method gets trapped in local maxima easily. In order to improve the DOA estimation precision, we design a group of spatial filters to separate the incident signals from the array output. After the separation, the original multiple-coherent-source array output decomposes into a group of single-source components, and the convergence of RVM based on each component obeys the unimodality described in Property (1) of Proposition 1.Define the required filters byT1,…,TK, with Tk used to separate the kth source from the array output, which satisfies
(5)TkA(ϑ)s=a′(ϑk)sk,k=1,…,K,
where a′(ϑk) is a certain function of ϑk. As s and sk are random and unknown, (5) can be transformed into the following equations for certain Tk:
(6)Tka(ϑi)=a′(ϑk)δ(k-i),i=1,…,K,
where δ(k-i) is the indicative function that takes nonzero value of 1 only when k-i=0.Equation (6) indicates that Tk is orthogonal to the column subspace of [a(ϑ1),…,a(ϑk-1),a(ϑk+1),…,a(ϑK)]; thus, the dimension of Tk is (M-K+1)×M and its explicit value is determined by Tka(ϑk)=a′(ϑk). An intuitive selection of a′(ϑk) is
(7)a′(ϑk)=[I(M-K+1)×(M-K+1),0(M-K+1)×(K-1)]a(ϑk),
which means reserving the array output from the first M-K+1 sensors. In the following, we present a way for designing such a Tk that satisfies (6) with a′(ϑk) given by (7).Denote the responding matrix of the subarray consisting of themth to (m+K-1)th elements of the original array to sources impinging from directions of ϑ by A(m)(ϑ), and also define B(m)(ϑ)=[A(m)(ϑ)]-1, with its kth row denoted by bk(m); then,
(8)B(m)(ϑ)A(m)(ϑ)s=Φms,m=1,…,M-K+1,bk(m)A(m)(ϑ)s=ejφk,msk,k=1,…,K,
where Φm=diag⁡([ejφ1,m,…,ejφK,m]).Based on the above property of the subfilters, we extend them in dimension to act on the whole array output to separate the incident signals. DefineTk∈ℂ(M-K+1)×M, whose mth row is (Tk)m•=[01×(m-1),bk(m),01×(M-K+1-m)]; then,
(9)yk≜Tkx=[ejφk,1,…,ejφk,M-K+1]Tsk+Tkv≜a′(ϑk)sk+Tkv,k=1,…,K,
where a′(ϑk)=[ejφk,1,…,ejφk,M-K+1]T. It can be concluded from (9) that the kth signal can be extracted from the array output x provided that the spatial filter Tk is well-designed. However, the parameters of the filters, that is, the signal directions, are not available in practice. Therefore, we introduce an iterative procedure to refine the signal directions and update the filter parameters simultaneously.


### 3.3. Bayesian Learning for Filter Update and DOA Refinement

It has been noted at the end of Section3.1 that the independent RVM-DOA can be used to obtain coarse DOA estimates for coherent sources, which are located in the vicinity of the true source directions with somewhat significant biases. We then refine those estimates iteratively via Bayesian learning.To keep in line with the iterative format of the learning techniques, we use(q) to identify the updating variables in the qth iteration. Denote the output of Tk(q) in the case of multiple snapshots by Yk(q)=[yk(q)(t1),…,yk(q)(tN)] and the extension of sk(q) from ϑk to Θ by s-k(q), which is zero-mean and Gaussian distributed as s-k(q)~𝒩(0,Γk(q-1)) with Γk(q-1)=diag⁡{γk(q-1)}, and the noise contained in Yk(q), that is, Tk(q)V, is also zero-mean and Gaussian distributed with variance of σ2Tk(q)Q(Tk(q))H; then, the posterior probability of S-k(q)=[s-k(q)(t1),…,s-k(q)(tN)] with Yk(q) given is
(10)p(S-k∣Yk(q),(σk2)(q-1))=p(Yk(q)∣S-k,(σk2)(q-1))p(S-k∣γk(q-1))∫p(Yk(q)∣S-k,(σk2)(q-1))p(S-k∣γk(q-1))d(S-k)=|πΣS-k(q)|-N×exp⁡{-tr⁡[(S-k-ℳk(q))H(ΣS-k(q))-1(S-k-ℳk(q))]},
where the source-intrinsic noise variance of σk2 is introduced because they will be updated independently and may differ from each other,
(11)ℳk(q)=Γk(q-1)A-H×(A-Γk(q-1)A-H+(σk2)(q-1)Tk(q)Q(Tk(q))H)-1Yk(q),ΣS-k(q)=Γk(q-1)-Γk(q-1)A-H×(A-Γk(q-1)A-H+(σk2)(q-1)Tk(q)Q(Tk(q))H)-1A-Γk(q-1).
The likelihood function of Yk(q) with respect to γk and σk2 is
(12)p(Yk(q)∣γk,σk2)=∫p(Yk(q)∣S-k,σk2)p(S-k∣γk)d(S-k)=|πΣYk(q)|-Nexp⁡{-tr⁡((Yk(q))H(ΣYk(q))-1Yk(q))},
where ΣYk(q)=A-Γk(q-1)A-H+(σk2)(q-1)TkQTkH.The logarithmic likelihoodℒk(γk,σk2)=ln⁡p(Yk(q)∣γk,σk2) can then be maximized by maximizing ℒk(γk;σk2) and ℒk(σk2;γk) alternatively with respect to γk and σk2. Firstly, ℒk(γk;σk2) can be maximized with the EM algorithm [18] as follows:
(13)(γk(q))l=argmax⁡(γk)l≥0ln⁡p(Yk(q)∣γk;(σk2)(q-1))=∥(ℳk(q))l•∥22N+(ΣS-k(q))l,l,
where ℳk(q) and ΣS-k(q) are computed according to (11), (·)l, (·)l•, and (·)l,l that represent the lth vector element, the lth matrix row, and the (l,l)th matrix element, respectively, and ∥·∥2 stands for the 2-norm of a vector. Based on γk(q) and Yk(q), the kth signal direction ϑ^k can be obtained and the noise variance is updated by projecting the filter output onto the noise subspace [12]; that is,
(14)(σk2)(q)=tr⁡(Pa′(ϑ^k)⊥R^Yk(q))tr⁡(Pa′(ϑ^k)⊥Tk(q)QTk(q)H),
where R^Yk(q)=Yk(q)(Yk(q))H/N, and Pa′(ϑ^k)⊥=I(M-K+1)×(M-K+1)-(a′H(ϑ^k)a′(ϑ^k))-1a′(ϑ^k)a′H(ϑ^k) is the orthogonal subspace of a′(ϑ^k). The way for estimating ϑ^k is presented in the following.Theγk(q) obtained in the above iteration reflects the power distribution of the kth filter output on the predefined discrete direction set, and further efforts should be made to obtain refined DOA estimate of the kth source, denoted by ϑ^k, from it. Assume that the peak cluster (the peak corresponding to each signal generally consists of two significant magnitudes on the left and right hand side of the true signal direction; we call them apeak cluster) of γk(q) is located at (ϑ_k,ϑ-k), and denote Θ-k=Θ∖{ϑ_k,ϑ-k}, that is, the direction set by removing ϑ_k and ϑ-k from Θ, Α¯-k′=A′(Θ-k), ΣYk-=Α¯-k′Γk-Α¯-k′H+σk2TkQTkH with Γk- equaling Γk=diag⁡(γk) by removing the two rows and columns that associate with (ϑ_k,ϑ-k); then, ΣYk- can be deemed to be the covariance matrix of Yk by removing the kth signal component. Complementing ΣYk- with the kth signal according to its power ηk and direction ϑk, instead of a peak cluster, yields ΣYk=ΣYk-+ηka′(ϑk)a′H(ϑk), and a refined direction estimate can be obtained by maximizing ℒk based on this ΣYk form according to ϑ^k=argmax⁡ηk,ϑkℒk as follows [12]:
(15)ϑ^k=argmax⁡θ∈[ϑ_k,ϑ-k)×|Re[a′H(θ)(ΣYk-)-1mmmm×(a′(θ)a′H(θ)(ΣYk-)-1R^Yk(q)mmmmmm-R^Yk(q)(ΣYk-)-1a′(θ)a′H(θ))mmmm×(ΣYk-)-1d′(θ)]|-1,
where d′(θ)=∂a′(θ)/∂θ. The signal direction estimate is derived from (15) according to the peak location by scanning the scope of [ϑ_k,ϑ-k) with required precision.After updating the source directions with (15), the spatial filters can be redesigned with them according to B(m)(ϑ^)=[A(m)(ϑ^)]-1 and (Tk)m•=[01×(m-1),bk(m),01×(M-K+1-m)]; the filter outputs can then be recomputed to prepare for the next iteration. Such iterations are repeated until a predefined convergence criterion is satisfied.We note here that the spatial filters used to separate the incident signals do not require any special property in the array geometry such as shift invariance [5–7]; hence, the new method can be applied more comprehensively than the existing counterparts.


### 3.4. Outline of the New Method

We name the new method spatial filtering-based relevance vector machine for DOA estimation (SF RVM-DOA for short), and outline its main procedures in Table1.Table 1
Outline of SF RVM-DOA.
Initialization: Use the independent RVM-DOA [12] to estimate source directions coarsely until primary convergence criterion is satisfied (and jointly detect their number if necessary).
DOA estimate refinement:
(1) Design spatial filters based onϑ^1,…,ϑ^K and compute Y1,…,YK.
(2) Updateγ^k and σ^k2 for k=1,…,K according to (13) and (14).
(3) Estimateϑ1,…,ϑK according to (15).
(4) Stop if final convergence criterion is satisfied; otherwise, go back to (1).Some supplements listed as follows should be added to complete the proposed method.Remark 2.
When turning from theInitialization procedure to theRefinement procedure, the γ^k’s are derived from γ^ by setting the other K-1 peak clusters to zeroes.Remark 3.
The EM procedure is generally much slow in its original form and is usually speeded up with the fixed-point iteration of(γk(q))l=(1/N)∥(ℳk(q))l•∥22/[1-(ΣS-k(q))l,l/(γk(q-1))l]+ς, where ς is a small positive value, such as 10-10 [12–14, 16].Remark 4.
The primary and final convergence criteria are application dependent. We set them as∥γ(q+1)-γ(q)∥2/∥γ(q+1)∥2<ε1 and ∥[γ1(q+1),…,γK(q+1)]-[γ1(q),…,γK(q)]∥F/∥[γ1(q+1),…,γK(q+1)]∥F<ε2 in this paper, where ε1≥ε2 and ∥·∥F stands for the Frobenius-norm of a matrix. In the simulation section, we set ε1=10-3 and ε2=10-4.


## 3.1. Overcomplete Model and Problem Formulation

In the sparse reconstruction-based DOA estimation methods, the potential space of the incident signals is sampled discretely to obtain a direction setΘ with cardinality of card(Θ)=L to contain all the possible signal directions. The corresponding array responding matrix is denoted by A-=A(Θ)∈ℂM×L. Generally, L is much larger than M to meet the requirement in superresolution; thus, A- is overcomplete. The array output given by (1) can then be rewritten as
(2)x=A-s-+v,
where s- is the extension of s from ϑ to Θ and has nonzero amplitudes only for the ϑ-associated bases. The model mismatch between (1) and (2) due to the discretization of Θ is skipped over to facilitate notation unless otherwise stated; such approximation does not influence the implementation of those methods. The signal directions can be estimated according to the locations of the nonzero magnitudes in s- if it can be estimated. However, the equation given by (2) is ill posed and it is impossible to estimate s- unless further sparsity constraints on s-, which indicates that only a few signals impinge simultaneously, are added. When multiple snapshots are collected, one needs to estimate S-=[s-(t1),…,s-(tN)] from X=[x(t1),…,x(tN)].In RVM, the columns ofS- are assumed to be independent and Gaussian distributed; that is, s-~𝒩(0,Γ), with Γ=diag⁡{γ} and γ=[γ1,…,γL]. Such an assumption leads to the following probability function of X with respect to γ and σ2:
(3)p(X∣γ,σ2)=∫p(X∣S-,σ2)p(S-∣γ)dS-=|πΣX|-Nexp⁡(-tr⁡(XHΣX-1X)),
where ΣX=A-ΓA-H+σ2Q. Its logarithmic form is used to estimate γ and σ2; that is,
(4)ℒ(γ,σ2)=ln⁡p(X∣γ,σ2)≡-Nln⁡|ΣX|-tr⁡(XHΣX-1X).
The function of ℒ(γ,σ2) can also be rewritten as ℒ(Γ,σ2) if intersource correlation exists and the diagonal structure of Γ is deviated.Then the expectation-maximization (EM) algorithm [18], which has already been used for wideband DOA estimation [19, 20], can be used to maximize ℒ(γ,σ2) and estimate the source directions [12]. We have shown in [12] that the RVM-based method, if well designed, surpasses its counterparts with significant margin in DOA estimation performance when the incident signals are independent. However, how it will behave in the coherent source settings still remains unknown. Fortunately, some conclusions available in the related neural source localization area have provided some clues about the behavior of the independent RVM when intersource correlation exists [17].Denotecpq=〈s-p•,s-q•〉/N, with s-p• standing for the pth row of S-; then, according to Theorem  1 in [17], the following propositions can be obtained to indicate the behavior of RVM when applied to independent and correlated DOA estimation.Proposition 1.
Suppose that the SNR is high and the number of the incident sources is smaller than that of the sensors, thenℒ(γ,σ2) and ℒ(Γ,σ2) own the following properties.(1)
Ifcpq=0for each two sources, thenℒ(γ,σ2)is unimodal and no local maxima exists.(2)
For arbitrarycpq, the unique global maximum ofℒ(Γ,σ2), denoted byΓ*(whileσ2is negligible), produces the true source matrixS-*=E(S-;X,Γ*); that is, it achieves perfect source reconstruction.Proposition1 is a straightforward extension of Theorem  1 in [17]. Property (1) of Proposition 1 implies that well-designed RVM is able to achieve consistent DOA estimation for independent sources, as our RVM-DOA method has done [12]. Property (2) tells that accurate source covariance matrix estimation, which is associated with accurate DOA estimates, is necessary for precise source matrix reconstruction in the case of coherent sources. However, as some of the off-diagonal elements of the coherent source covariance matrix have significant nonzero magnitudes, the original RVM should be modified to take all possible intersource correlations into consideration, so as to make global convergence achievable. Unfortunately, all the L2 elements of Γ should be estimated under such a modified model, and the great increase in the number of unknown parameters not only aggravates the computational burden of the RVM prohibitively, but also increases the risk of overfitting [21].In addition, it should be noted that, although the independent RVM-based DOA estimation method is not able to perform well enough in the case of coherent sources, it still converges to the vicinity of the global maximum [22]. Therefore, we use the independent RVM-DOA method in [12] to obtain coarse DOA estimates in this paper.


## 3.2. Basic Ideas for Designing the Spatial Filters

According to Proposition1 and the analysis that follows, the global maximum of ℒ(Γ,σ2) can hardly be achieved by RVM in the case of coherent sources, and the corresponding DOA estimation method gets trapped in local maxima easily. In order to improve the DOA estimation precision, we design a group of spatial filters to separate the incident signals from the array output. After the separation, the original multiple-coherent-source array output decomposes into a group of single-source components, and the convergence of RVM based on each component obeys the unimodality described in Property (1) of Proposition 1.Define the required filters byT1,…,TK, with Tk used to separate the kth source from the array output, which satisfies
(5)TkA(ϑ)s=a′(ϑk)sk,k=1,…,K,
where a′(ϑk) is a certain function of ϑk. As s and sk are random and unknown, (5) can be transformed into the following equations for certain Tk:
(6)Tka(ϑi)=a′(ϑk)δ(k-i),i=1,…,K,
where δ(k-i) is the indicative function that takes nonzero value of 1 only when k-i=0.Equation (6) indicates that Tk is orthogonal to the column subspace of [a(ϑ1),…,a(ϑk-1),a(ϑk+1),…,a(ϑK)]; thus, the dimension of Tk is (M-K+1)×M and its explicit value is determined by Tka(ϑk)=a′(ϑk). An intuitive selection of a′(ϑk) is
(7)a′(ϑk)=[I(M-K+1)×(M-K+1),0(M-K+1)×(K-1)]a(ϑk),
which means reserving the array output from the first M-K+1 sensors. In the following, we present a way for designing such a Tk that satisfies (6) with a′(ϑk) given by (7).Denote the responding matrix of the subarray consisting of themth to (m+K-1)th elements of the original array to sources impinging from directions of ϑ by A(m)(ϑ), and also define B(m)(ϑ)=[A(m)(ϑ)]-1, with its kth row denoted by bk(m); then,
(8)B(m)(ϑ)A(m)(ϑ)s=Φms,m=1,…,M-K+1,bk(m)A(m)(ϑ)s=ejφk,msk,k=1,…,K,
where Φm=diag⁡([ejφ1,m,…,ejφK,m]).Based on the above property of the subfilters, we extend them in dimension to act on the whole array output to separate the incident signals. DefineTk∈ℂ(M-K+1)×M, whose mth row is (Tk)m•=[01×(m-1),bk(m),01×(M-K+1-m)]; then,
(9)yk≜Tkx=[ejφk,1,…,ejφk,M-K+1]Tsk+Tkv≜a′(ϑk)sk+Tkv,k=1,…,K,
where a′(ϑk)=[ejφk,1,…,ejφk,M-K+1]T. It can be concluded from (9) that the kth signal can be extracted from the array output x provided that the spatial filter Tk is well-designed. However, the parameters of the filters, that is, the signal directions, are not available in practice. Therefore, we introduce an iterative procedure to refine the signal directions and update the filter parameters simultaneously.


## 3.3. Bayesian Learning for Filter Update and DOA Refinement

It has been noted at the end of Section3.1 that the independent RVM-DOA can be used to obtain coarse DOA estimates for coherent sources, which are located in the vicinity of the true source directions with somewhat significant biases. We then refine those estimates iteratively via Bayesian learning.To keep in line with the iterative format of the learning techniques, we use(q) to identify the updating variables in the qth iteration. Denote the output of Tk(q) in the case of multiple snapshots by Yk(q)=[yk(q)(t1),…,yk(q)(tN)] and the extension of sk(q) from ϑk to Θ by s-k(q), which is zero-mean and Gaussian distributed as s-k(q)~𝒩(0,Γk(q-1)) with Γk(q-1)=diag⁡{γk(q-1)}, and the noise contained in Yk(q), that is, Tk(q)V, is also zero-mean and Gaussian distributed with variance of σ2Tk(q)Q(Tk(q))H; then, the posterior probability of S-k(q)=[s-k(q)(t1),…,s-k(q)(tN)] with Yk(q) given is
(10)p(S-k∣Yk(q),(σk2)(q-1))=p(Yk(q)∣S-k,(σk2)(q-1))p(S-k∣γk(q-1))∫p(Yk(q)∣S-k,(σk2)(q-1))p(S-k∣γk(q-1))d(S-k)=|πΣS-k(q)|-N×exp⁡{-tr⁡[(S-k-ℳk(q))H(ΣS-k(q))-1(S-k-ℳk(q))]},
where the source-intrinsic noise variance of σk2 is introduced because they will be updated independently and may differ from each other,
(11)ℳk(q)=Γk(q-1)A-H×(A-Γk(q-1)A-H+(σk2)(q-1)Tk(q)Q(Tk(q))H)-1Yk(q),ΣS-k(q)=Γk(q-1)-Γk(q-1)A-H×(A-Γk(q-1)A-H+(σk2)(q-1)Tk(q)Q(Tk(q))H)-1A-Γk(q-1).
The likelihood function of Yk(q) with respect to γk and σk2 is
(12)p(Yk(q)∣γk,σk2)=∫p(Yk(q)∣S-k,σk2)p(S-k∣γk)d(S-k)=|πΣYk(q)|-Nexp⁡{-tr⁡((Yk(q))H(ΣYk(q))-1Yk(q))},
where ΣYk(q)=A-Γk(q-1)A-H+(σk2)(q-1)TkQTkH.The logarithmic likelihoodℒk(γk,σk2)=ln⁡p(Yk(q)∣γk,σk2) can then be maximized by maximizing ℒk(γk;σk2) and ℒk(σk2;γk) alternatively with respect to γk and σk2. Firstly, ℒk(γk;σk2) can be maximized with the EM algorithm [18] as follows:
(13)(γk(q))l=argmax⁡(γk)l≥0ln⁡p(Yk(q)∣γk;(σk2)(q-1))=∥(ℳk(q))l•∥22N+(ΣS-k(q))l,l,
where ℳk(q) and ΣS-k(q) are computed according to (11), (·)l, (·)l•, and (·)l,l that represent the lth vector element, the lth matrix row, and the (l,l)th matrix element, respectively, and ∥·∥2 stands for the 2-norm of a vector. Based on γk(q) and Yk(q), the kth signal direction ϑ^k can be obtained and the noise variance is updated by projecting the filter output onto the noise subspace [12]; that is,
(14)(σk2)(q)=tr⁡(Pa′(ϑ^k)⊥R^Yk(q))tr⁡(Pa′(ϑ^k)⊥Tk(q)QTk(q)H),
where R^Yk(q)=Yk(q)(Yk(q))H/N, and Pa′(ϑ^k)⊥=I(M-K+1)×(M-K+1)-(a′H(ϑ^k)a′(ϑ^k))-1a′(ϑ^k)a′H(ϑ^k) is the orthogonal subspace of a′(ϑ^k). The way for estimating ϑ^k is presented in the following.Theγk(q) obtained in the above iteration reflects the power distribution of the kth filter output on the predefined discrete direction set, and further efforts should be made to obtain refined DOA estimate of the kth source, denoted by ϑ^k, from it. Assume that the peak cluster (the peak corresponding to each signal generally consists of two significant magnitudes on the left and right hand side of the true signal direction; we call them apeak cluster) of γk(q) is located at (ϑ_k,ϑ-k), and denote Θ-k=Θ∖{ϑ_k,ϑ-k}, that is, the direction set by removing ϑ_k and ϑ-k from Θ, Α¯-k′=A′(Θ-k), ΣYk-=Α¯-k′Γk-Α¯-k′H+σk2TkQTkH with Γk- equaling Γk=diag⁡(γk) by removing the two rows and columns that associate with (ϑ_k,ϑ-k); then, ΣYk- can be deemed to be the covariance matrix of Yk by removing the kth signal component. Complementing ΣYk- with the kth signal according to its power ηk and direction ϑk, instead of a peak cluster, yields ΣYk=ΣYk-+ηka′(ϑk)a′H(ϑk), and a refined direction estimate can be obtained by maximizing ℒk based on this ΣYk form according to ϑ^k=argmax⁡ηk,ϑkℒk as follows [12]:
(15)ϑ^k=argmax⁡θ∈[ϑ_k,ϑ-k)×|Re[a′H(θ)(ΣYk-)-1mmmm×(a′(θ)a′H(θ)(ΣYk-)-1R^Yk(q)mmmmmm-R^Yk(q)(ΣYk-)-1a′(θ)a′H(θ))mmmm×(ΣYk-)-1d′(θ)]|-1,
where d′(θ)=∂a′(θ)/∂θ. The signal direction estimate is derived from (15) according to the peak location by scanning the scope of [ϑ_k,ϑ-k) with required precision.After updating the source directions with (15), the spatial filters can be redesigned with them according to B(m)(ϑ^)=[A(m)(ϑ^)]-1 and (Tk)m•=[01×(m-1),bk(m),01×(M-K+1-m)]; the filter outputs can then be recomputed to prepare for the next iteration. Such iterations are repeated until a predefined convergence criterion is satisfied.We note here that the spatial filters used to separate the incident signals do not require any special property in the array geometry such as shift invariance [5–7]; hence, the new method can be applied more comprehensively than the existing counterparts.


## 3.4. Outline of the New Method

We name the new method spatial filtering-based relevance vector machine for DOA estimation (SF RVM-DOA for short), and outline its main procedures in Table1.Table 1
Outline of SF RVM-DOA.
Initialization: Use the independent RVM-DOA [12] to estimate source directions coarsely until primary convergence criterion is satisfied (and jointly detect their number if necessary).
DOA estimate refinement:
(1) Design spatial filters based onϑ^1,…,ϑ^K and compute Y1,…,YK.
(2) Updateγ^k and σ^k2 for k=1,…,K according to (13) and (14).
(3) Estimateϑ1,…,ϑK according to (15).
(4) Stop if final convergence criterion is satisfied; otherwise, go back to (1).Some supplements listed as follows should be added to complete the proposed method.Remark 2.
When turning from theInitialization procedure to theRefinement procedure, the γ^k’s are derived from γ^ by setting the other K-1 peak clusters to zeroes.Remark 3.
The EM procedure is generally much slow in its original form and is usually speeded up with the fixed-point iteration of(γk(q))l=(1/N)∥(ℳk(q))l•∥22/[1-(ΣS-k(q))l,l/(γk(q-1))l]+ς, where ς is a small positive value, such as 10-10 [12–14, 16].Remark 4.
The primary and final convergence criteria are application dependent. We set them as∥γ(q+1)-γ(q)∥2/∥γ(q+1)∥2<ε1 and ∥[γ1(q+1),…,γK(q+1)]-[γ1(q),…,γK(q)]∥F/∥[γ1(q+1),…,γK(q+1)]∥F<ε2 in this paper, where ε1≥ε2 and ∥·∥F stands for the Frobenius-norm of a matrix. In the simulation section, we set ε1=10-3 and ε2=10-4.


## 4. Simulation Results

In this section, we carry out simulations to demonstrate the performance of the proposed method (named SF RVM-DOA) in coherent DOA estimation. Four other existing methods are also listed for comparison, including spatial smoothing MUSIC (named SS MUSIC) [5], spatial filtering MUSIC (SF MUSIC) [7], L1-SVD [9], and independent RVM-DOA [12]. Among those methods, L1-SVD, RVM-DOA, and SF RVM-DOA adapt to arbitrary array geometries, while SS MUSIC and SF MUSIC require shift invariance property in the array geometry. Therefore, in order to facilitate comparison, we choose a 10-element ULA for the simulations, whose interelement spacing equals half wavelength. Moreover, as only RVM-DOA and SF RVM-DOA realize joint detection and estimation, while the other three methods need to know the source number beforehand, we assume that the source number is known to make the DOA estimation performance comparison fairer. The three sparsity-inducing methods, that is, L1-SVD, RVM-DOA, and SF RVM-DOA, sample the [−90°90°] space with 1° interval to obtain the direction set Θ, and L1-SVD refines the grid to 0.01° afterwards to improve the DOA estimation precision. The SF MUSIC method is stopped after 5 repetitions, and the maximum iteration numbers of RVM-DOA and SF RVM-DOA are capped at 1000. Two equal-power sources are assumed to impinge onto the array simultaneously, and 300 independent trials are carried out in each scenario. The additive noise is assumed to be spatially white; that is, Q=I.In the first experiment, we fix the snapshot number at 30 and the angle separation of the two sources at 10°, and their directions are-5°+υ and 5°+υ, respectively, with υ chosen uniformly and randomly in [−0.5°0.5°] in each trial. The differential phase shift of the two signals is set randomly within [02π] in each trial. The SNR of each signal is increased from 0 dB to 20 dB, and the averaged root-mean-square errors (RMSEs) of the DOA estimates of the two sources are given in Figure 1, with the CRLB depicted by the dash dotted line. The results in Figure 1 indicate that the three sparsity-inducing methods succeed in separating the two sources at lower SNR, but among them, L1-SVD and RVM-DOA have larger DOA estimation RMSE than the two subspace-based methods when the SNR is higher than 6 dB. SF RVM-DOA achieves high-precision DOA estimation at a SNR of about 2 dB, which is 5 dB lower than that of the two subspace-based methods, and its DOA estimation precision surpasses the other four methods above this threshold and gradually approaches the CRLB as the SNR increases.Figure 1
Averaged DOA estimation RMSE of two equal-power coherent sources separated by 10° with 30 snapshots when the SNR increases.Then fix the angle separation of the two sources at 10° and their SNR at 7 dB; the directions and phase shift of the two sources are set similarly as the first experiment. When the snapshot number increases from 1 to 100, the averaged DOA estimation RMSE of the five methods is shown in Figure2. The results indicate that the three sparsity-inducing methods are able to separate the two sources with only 4 snapshots, while the two subspace-based methods require more than 30 snapshots. As the snapshot number increases, the DOA estimation precision of L1-SVD and RVM-DOA improves only slightly and RVM-DOA surpasses L1-SVD with a small margin. However, the RMSE of the other three methods decreases almost synchronously with the CRLB, and SF RVM-DOA obtains higher precision than the other two methods and is closer to the CRLB.Figure 2
Averaged DOA estimation RMSE of two 7 dB coherent sources separated by 10° when the snapshot number increases.Finally, fix the SNR of the two sources at 7 dB and the snapshot number at 30, and increase the angle separation of them from 6° to 16°. The directions of the two sources are set to be-(ϑ2-ϑ1)/2+υ and (ϑ2-ϑ1)/2+υ, with υ and the intersource phase-shift randomly chosen as that in the above two experiments. The averaged DOA estimation RMSE of the five methods in such a scenario is given in Figure 3. It can be concluded that, although L1-SVD owns enhanced superresolution than the other methods, its performance improves only slightly as the two sources are farther apart. When the two sources are separated by more than 10°, L1-SVD performs the worst among all the five methods. RVM-DOA and SF RVM-DOA surpass L1-SVD in DOA estimation precision when the angle separation is larger than 8° and improves at a higher speed than it, and SF RVM-DOA gains the highest precision. The two subspace-based methods can only separate two sources spaced by more than 10°, with their DOA precision lower than SF RVM-DOA all along, and is also exceeded by RVM-DOA when the angle separation of the two sources is larger than 14°.Figure 3
Averaged DOA estimation RMSE of two 7 dB coherent sources with 30 snapshots when their angle separation increases.


## 5. Conclusions and Discussions

In this paper, we introduced the relevance vector machine (RVM) to estimate the directions of coherent sources and designed a group of spatial filters to separate the incident signals from the array output and guide the RVM to better convergence. The proposed method adapts to arbitrary array geometry, and simulation results show that it shares similar superiority with other sparsity-inducing methods in adaptation to low SNR, limited snapshots, and spatially adjacent sources. Moreover, it performs much better than all the subspace-based and sparsity-inducing counterparts in DOA estimation precision in moderate settings.The spatial filters designed in this paper can also be used to replace the ordinary spatial filtering and smoothing techniques [5–7] and combined with the ordinary subspace-based DOA estimation methods, so as to extend them to arbitrary array geometry. However, the performance of such variants may still be constrained by the dependence of the subspace-based methods on high-precision covariance matrix and subspace estimates; thus, they may deteriorate significantly in the demanding environments of low SNR, limited snapshots, or spatially adjacent sources.



---

*Source: 959386-2014-04-27.xml*