KUNGL TEKNISKA HÖGSKOLANInstitutionen förSignaler, Sensorer & SystemSignalbehandling100 44 STOCKHOLMROYAL INSTITUTEOF TECHNOLOGY
Department of
Signals, Sensors & Systems
Signal Processing
S-100 44 STOCKHOLM
ALG
WSF
1020304050607080
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
DOA (deg)
b
i
a
s
(
d
e
g
)
ALG
WSF
1012141618202224262830
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
DOA (deg)
B
i
a
s
(
d
e
g
)
第47卷第6期世界地震译丛V S.47，N o.62016 年11 月（543 — 548) TRANSLATED WORLD SEISMOLOGY n〇〇"2016N. M. S h a p iro»M. C am pillo. 2004. E m ergence of b ro ad ban d R ayleigh w aves fro m correlation s of th Geophys ,Res, Lett.31(7) , L07614 , d o i； 10. 1029/2004G L019491赵爰华译.2016.从环境地震噪声的互相关中产生的宽带瑞利波.世界地震译丛.47(6)543 — 548. d i10. 16738,.cnki.issn. 1003 — 3238. 201606006从环境地震噪声的互相关中产生的宽带瑞利波N.M.Shapiro M.Campillo摘要证明了从环境地震噪声中可以提取地球结构的有关信息。
对间隔从约100千米至2000多千米不同的台站对，计算了垂直分量上数天噪声的互相关，清楚地得到了相干的宽带频散波列，其群速度与根据全球瑞利波层析成像图（以弹道式面波约束）所预测的相似。
计算的结果表明可以从环境噪声中提取相干的瑞利波，并能在很宽的周期范围内测量频散特性。
这为那些不能对弹道式波进行采样的大量路径提供了 一种新的面波测量之源，从而可显著地提高地震图像的分辨率。
0引言地震成像传统上基于爆炸或地震等震源 产生的直达波。
弹道式波则用于测量体波的 走时、面波的频散曲线和地球简正模式的频 率。
测量结果的反演提供了有关地球内部的 基础性信息，从而得以认识地球的结构、球 对称的分层以及地震性质在三维及各向异性 上的变化。
地震图像分辨率的进一步提高主要受限 于震源及检波器的非均匀分布，它们对弹道 式地震波测量的分辨率有些根本性的限制。
在本文中，我们集中探讨地震面波，它是地 壳与上地幔结构的主要信息来源之一。
Sparsity-Inducing Direction Finding for Narrowband and Wideband Signals Based on Array Covariance Vectors Zhang-Meng Liu,Member,IEEE,Zhi-Tao Huang,and Yi-Yu ZhouAbstract—Among the existing sparsity-inducing direction-of-arrival(DOA)estimation methods,the sparse Bayesian learning (SBL)based ones have been demonstrated to achieve enhanced precision.However,the learning process of those methods con-verges much slowly when the signal-to-noise ratio(SNR)is relatively low.In this paper,wefirst show that the covari-ance vectors(columns of the covariance matrix)of the array output of independent signals share identical sparsity profiles corresponding to the spatial signal distribution,and their SNR exceeds that of the raw array output when moderately many snapshots are collected.Thus the SBL technique can be used to estimate the directions of independent narrowband/wideband signals by reconstructing those vectors with high computational efficiency.The method is then extended to narrowband correlated signals after proper modifications.In-depth analyses are also provided to show the lower bound of the new method in DOA estimation precision and the maximal signal number it can separate in the case of independent signals.Simulation results finally demonstrate the performance of the proposed method in both DOA estimation precision and computational efficiency.Index Terms—Direction-of-arrival(DOA)estimation,sparse reconstruction,relevance vector machine(RVM),covariance vector.I.I NTRODUCTIONT HE array output covariance matrix contains the direc-tional information of the incident signals and well con-centrates the signal energy distributed in all the snapshots, thus it is widely exploited in the direction-of-arrival(DOA) estimation methods for both narrowband[1]and wideband [2],[3]signals.In particular applications,the calculation of the covariance matrix may even greatly facilitate the measurement formulation of wideband signals[4].However,most of the existing covariance matrix-based DOA estimation methods require the prior information of the incident signal number[1]-[3],and the spectral decomposition and focusing procedures of the wideband array outputs may introduce extra imperfections into the measurements[2],[3].The sparse reconstruction techniques have attracted much interest recently in various areas including wireless communi-cations,with special attention in thisfield paid to channel esti-Manuscript received August29,2012;revised December18,2012and February27,2013;accepted April15,2013.The associate editor coordinating the review of this paper and approving it for publication was A.Zajic.This work was supported in part by the National Natural Science Founda-tion(NO.61072120).The authors are with the School of Electronic Science and Engineering, National University of Defense Technology,Changsha,410073,China(e-mail:zm liur@;taldcn@).Digital Object Identifier10.1109/TWC.2013.071113.121305mation[5]-[7],data gathering[8],interference estimation[9], multiuser detection[10],etc..Another important application of those techniques lays in DOA estimation with antenna arrays [4],[11]-[14].The sparsity-inducing DOA estimation methods make use of the spatial sparsity of the incident signals,and they succeed to estimate the signal directions by reconstructing the array outputs on a directional overcomplete dictionary. Previous simulation results have demonstrated the superiority of those methods in superresolution and robustness,especially in much demanding scenarios with low signal-to-noise ratio (SNR),limited snapshots and spatially adjacent signals[4], [11]-[14].The existing sparsity-inducing DOA estimation methods mainly divide into two categories,the p-norm(0≤p≤1) based ones[4],[11]-[13]and the sparse Bayesian learn-ing(SBL)based ones[14].It has been demonstrated both theoretically and empirically that,the SBL technique[15] induces less structural error(biased global minimum)and convergence error(failure in achieving the global minimum) than the p-norm(0≤p≤1)based ones[16],[17],and it also performs well in capturing local signal properties to facilitate the refined DOA estimation process[14].Several other literatures have also introduced the Bayesian idea to thefield of array signal processing[18]-[20].Nonetheless, the beamformers generally do not perform satisfyingly enough in superresolution[18],[19],and further research is required to apply the periodic cost function[20]in scenarios of multiple signals with unknown waveforms.A more relevant work with the method to be proposed is the relevance vector machine(RVM)-based DOA estimator[14].The estimator has been shown to surpass its subspace-based and p-norm-based counterparts in adaptation to much demanding scenarios and DOA estimation precision.However,the simulation results in [14]also indicate a significant drawback of the RVM-DOA method,i.e.,its computational efficiency deteriorates rapidly as the signal-to-noise ratio(SNR)decreases due to the slowed down convergence rate of the reconstruction process.In this paper,we calculate the covariance matrix of the array outputfirst as is usually done in the subspace-based DOA estimators,and realize DOA estimation of independent narrowband/wideband signals and correlated narrowband sig-nals by reconstructing the covariance vectors(columns of the covariance matrix)sparsely in the spatial domain,rather than exploiting the orthogonality between the signal-and noise-subspaces.The method to be proposed is also RVM-based1536-1276/13$31.00c 2013IEEEas the RVM-DOA method [14]and we name it CV-RVM with CV denoting Covariance Vectors.The major motivation for us to resort to this strategy from that in [14]is that,the SNR of those vectors is higher than that of the raw array outputs when moderately many snapshots are collected,while the vectors still share identical sparsity profiles as the outputs.The enhanced SNR is expected to make up for the drawback of the RVM-DOA method in computational efficiency.Moreover,due to the structure differences of the covariance vectors and the raw array outputs,it is not an easy task to derive the CV-based DOA estimator straightforwardly from the work in [14],and special attentions should be paid to the implementation and properties of CV-RVM.The idea of estimating narrowband and wideband signal directions by reconstructing the covariance vectors can also be found in [12]and [4].However,those methods can hardly obtain DOA estimates with satisfying precision due to the shortcomings of the p -norm reconstruction techniques they use,as we will show in the simulations of this paper.The rest of the paper mainly consists of six parts.Section II reviews the formulation of the array covariance vectors.In Section III,we analyze the first-and second-order statistics of the covariance vector estimation error in the case of limited snapshots,and show that the SNR of the vectors exceeds that of the raw array outputs if moderately many snapshots are collected.In Section IV ,we introduce the SBL technique to estimate the directions of independent narrowband and wideband signals by reconstructing those vectors,and then extend the proposed method to correlated narrowband signals.The theoretical lower bound of the proposed method in DOA estimation precision and the maximal source number that it is able to separate in the case of independent signals are analyzed in Section V .Numerical examples are carried out in Section VI to demonstrate the performance of the proposed method.Section VII concludes the whole paper.II.M ODEL F ORMULATIONSuppose that K stochastic Gaussian signals impinge onto an M -element array from directions of ϑ=[ϑ1,···,ϑK ]simultaneously,the output of the m th sensor at time t isx m (t )=K k =1s k (t +τk,m )+v m (t ),(1)where s k (t )is the k th signal waveform,τk,m is the prop-agating time-delay of this signal between the m th sensorand the reference,v m (t )is the independent white Gaussian noise with variance σ2.Suppose that the array sensors are located in the same 2-D plane,then τk,m =d Tm g k νwithg k =[sin ϑk ,cos ϑk ]T,d m being the location of the m th sensor and νthe propagation velocity of the waveforms.If sufficient snapshots are collected,one can obtain the perturbation-free covariance matrix R asR =limN →+∞1NNn =1x (t n )x H (t n ),(2)where x (t )=[x 1(t ),···,x M (t )]T,(•)Hand (•)Tare the conjugate transpose and transpose operators,respectively.Thesampling rate of the array receiver is assumed to keep constantduring the observation time,i.e.,t n =(n −1)T s with T s being the sampling interval.The (m 1,m 2)th element of R derived from (1)is given as follows,R m 1,m 2=K k =1K k =1E [s k (t +τk,m 1)s ∗k (t +τk ,m 2)]+σ2δ(m 1−m 2),m 1,m 2=1,···,M,(3)where E (•)and (•)∗stand for the expectation and conjugate operators,respectively,and δ(•)is the indicator function.The structure of R can be simplified further by neglecting the correlation items in (3)in the case of independent narrowband/widebandsignals,where E [s k (t +τk,m 1)s ∗k (t +τk ,m 2)]=ηk r k (τk,m 1−τk,m 2)δ(k −k)with ηk =E |s k (t )|2 being the power of the k th signal,r k (τ)being the unified self-correlation of this signal at time-delay τthat satisfies r k (0)=1.Then the m th column of R for independent signals can be formulated as follows,y m =[R 1,m ,···,R M,m ]T=K k =1ηk b k,m +σ2e m=B m η+σ2e m ,(4)where b k,m =[r k (τk,1−τk,m ),···,r k (τk,M −τk,m )]T,B m =[b 1,m ,···,b K,m ],η=[η1,···,ηK ]T,and e m is a M ×1vector with the m th element being 1and the others being 0.The connection of b k,m to ϑk can be made more explicitly by combining τk,m =d T m g k νto conclude in b k,m =r k(d 1−d m )T g k ν ,···,r k (d M −d m )Tg k ν T ,with d m ,g k and νdefined in the same way as those in (1).Eq.(4)indicates that each column of R ,after removing the item of σ2e m ,is a weighted summation of the K column vectors in B m ,and the vector formed by aligning the M columns of R one-after-another has the following expression,y =vec (R )=Bη+σ2˜e,(5)where vec (•)is the vectorization operator that forms a vector satisfying [y ](m 2−1)×M +m 1=R m 1,m 2,B = B T 1,···,B T M T and ˜e=vec (I M ).Eq.(5)gives a re-shaped equation of the M equations given in (4)for m =1,···,M .Since the vectors of b 1,m ,···,b K,m and matrices of B 1,···,B M and B rely on the time-delays of the signals through the array,they are signal direction-dependent,we can thus add the directional label to them as b k,m (ϑk ),B m (ϑ)and B (ϑ)when necessary.Suppose that the incident signals are narrowband ones,or wideband ones with identical and known modulations as is assumed in [4]and [21],the columns of B rely only on the signal directions.Therefore,those directions can be estimated by recovering the K directional components from y .In this paper,we name both y 1,···,y M and y as covariance vectors,and seek to estimate the signal directions by decomposing them in the spatial domain.In the case of correlated signals,the columns of R also contain inter-signal correlation items and they cannot be simplified as that in (4).If the incident signals are wide-band ones,those cross-correlation items rely on both the signal directions and multipath delays,which complicates the covariance vectors significantly.We skip over the corre-lated wideband scenarios due to space limitation and leave it for future research.However,when the incident signals are correlated narrowband ones,the cross-correlation items can be simplified as E [s k (t +τk,m 1)s ∗k (t +τk ,m 2)]=αk,k √ηk ηk ξk,m 1ξ∗k ,m 2,where αk,k is the correlation co-efficient between the k th and k th signals,f is the frequency shared by all the signals and ξk,m =exp (j 2πfτk,m ).Then the m th column of R can be rewritten asy m =K k =1 K k =1αk,k √ηk ηk ξ∗k ,ma (ϑk )+σ2e m =A (ϑ)u m +σ2e m ,(6)where u m =[g m,1,···,g m,K ]T,g m,i =K k =1αi,k √ηi ηk ξ∗k,m ,A (ϑ)=[a (ϑ1),···,a (ϑK )]and a (ϑk )=[ξk,1,···,ξk,M ]T.The signal directions are used as identifiers in A (ϑ)and a (ϑk )because the time delays are direction-dependent.It can be concluded from (6)that,after removing the item of σ2e m ,the covariance vectors of narrowband correlated signals also consist of K directional components,while the weight vector of u m varies with the column index.Therefore,the directions of correlated narrowband signals can also be estimated by reconstructing the covariance vectors.Based on the above formulation of the covariance vectors,we propose a DOA estimator named covariance vector-based relevance vector machine,CV-RVM for short,to estimate the directions of independent narrowband/wideband signals and correlated narrowband signals.In the case of independent signals,the estimation error-contaminated counterpart of y is taken as the single measurement to accomplish the SMV (single measurement vector)implementation of CV-RVM.When the incident signals are correlated narrowband ones,the vectors of y 1,···,y M are reconstructed jointly for DOA estimation,which forms the multiple measurement vector (MMV)implementation of CV-RVM.In order to distinguish those two implementations of the new method,we name them SMV CV-RVM and MMV CV-RVM,respectively.III.P ROPERTIES OF THE A RRAY O UTPUT C OVARIANCEV ECTORS In this part,we first analyze the first-and second-order statistics of the estimation error of the covariance vectors caused by finite sampling.Then we compare the SNR of the covariance vectors with that of the raw array output,so as to partially verify the motivation for us to resort from the DOA estimator in [14]to the one in this paper.A.first-and second-order statistics of the covariance vector estimation errorIn practical applications,the covariance matrix can only be estimated using the N snapshots collected at time instants oft =t 1,···,t N as follows,ˆR =1N Nn =1x (t n )x H (t n ).(7)The covariance matrix estimate is estimation error-contaminated due to finite sampling.Denote E =ˆR−R ,E =[ε1,···,εM ]and ε=vec (E ),then the covariancevectors can be formulated as ˆy m =y m +εm and ˆy =y +ε,where ˆym =ˆRe m ,ˆy =vec ˆRand the expressions of y m and y are given in (4)-(6).In order to better distinguishv (t )=[v 1(t ),···,v M (t )]Tcontained in the raw array output and the εm ’s in the covariance vector estimates,we call v (t )”noise”and ε”perturbation”or ”estimation error”in the rest of the paper.By combining (1)and (7),one can obtain the explicitexpression of ˆRm 1,m 2,i.e.,the (m 1,m 2)th element of ˆR ,by taking the effect of finite sampling into account as follows,ˆRm 1,m 2=K k =1K k =11N N n =1s k (t n +τk,m 1)s ∗k (t n +τk ,m 2)+Kk =11N Nn =1s k (t n +τk,m 1)v ∗m 2(t n )+K k =11NN n =1s ∗k (t n +τk,m 2)v m 1(t n )+1NNn =1v m 1(t n )v ∗m 2(t n ).(8)Thus the expression of E m 1,m 2=[E ]m 1,m 2can be obtainedby combining (8)and (3)as follows,E m 1,m 2=K k =1K k =11NN n =1s k (t n +τk,m 1)s ∗k (t n +τk ,m 2)−E [s k (t +τk,m 1)s ∗k (t +τk ,m 2)]]+K k =11N N n =1s k (t n +τk,m 1)v ∗m 2(t n )+K k =11N N n =1s ∗k (t n +τk,m 2)v m 1(t n )+1NNn =1v m 1(t n )v ∗m 2(t n )−σ2δ(m 1−m 2) .(9)Denote the four items on the right hand side of (9)by υ1,···,υ4,they are stochastic due to the randomicity of the signal and noise amplitudes.Then it can be concluded from the zero-mean property and mutual independence of the signal and noise sequences that,E (E m 1,m 2)=4 i =1E (υi )=0.(10)It can also be concluded that E m 1,m 2is Gaussian distributeddue to the Gaussian distribution of the signal and noise amplitudes according to the law of large numbers when N is moderately large,and the second order statistic of E m 1,m 2can be obtained via straightforward calculation that (with detailedderivation provided in the Appendix),E E m 1,m 2E ∗m 1,m2=1N Δt =nT sR m 1,m 1(Δt )R ∗m 2,m 2(Δt )+O 1N 2 ,(11)where R (Δt )=E x (t +Δt )x H(t ) ,R m 1,m 2(Δt )=[R (Δt )]m 1,m 2,and O 1N 2 stands for a constant having ascaled magnitude of 1N 2with the scale being 0when the incident signals are narrowband.ThusE εm εH m =1N Δt =nT s[R (Δt )]m,m R (Δt )+O 1N 2Δ=Q m ,(12)andE εεH =1NΔt =nT s R (Δt )T ⊗R (Δt )+O1N 2Δ=Q ,(13)where ⊗represents the Kronecker product.It should be noted that,R (Δt )=δ(Δt )R and (13)can be simplified to E εεH =(1/N )R T⊗R for narrowband signals,which tallies with the results in [22,Chp.4],while R (Δt )has nonzero values for small Δt ’s in wideband scenarios and the simplification no longer holds,e.g.,it is nonzero for all |Δt |<1/B if the signals are PN ones with code rate B [4],[21].Nonetheless,Q m and Q have nonzero off-diagonal ele-ments for both narrowband and wideband signals according to (12)and (13),thus it can be concluded that the estimation errors of different covariance vector elements are correlated.Strict analysis of such correlations is necessary to facilitate the implementation of any maximum likelihood or Bayesian method.Another point that should be noted is that,the perturbation-free entities of Q m ,Q ,R and σ2are actually unknown beforehand,but we use them directly during the theoretical analyses without any remark for notational conve-nience.The way for estimating those entities will be provided in the simulation section.B.SNR of the covariance vectorsIn the following,we take independent narrowband signals for example to analyze the SNR of the covariance vectors theoretically,and compare it with that of the raw array output.In such scenarios,the array SNR (ASNR)of the k th signal in the raw array outputs can be calculated based on (1)asASNR k =ηkσ2.(14)The formulation of the m th covariance vector can be derived from (4)as follows,ˆym =B m η+σ2e m +εm =A (ϑ)w m +σ2e m +εm ,(15)where B m =A (ϑ)Φm for independent narrowband signals,Φm =diag ξ∗1,m ,···,ξ∗K,mT ,w m =Φm η,A (ϑ),ξk,m and e m are defined in the same way as those in (6),and εm is the estimation error.As the estimation errors of different covariance elements are correlated,a decorrelation process should be introduced tofacilitate the calculation of the SNR of the covariance vectors.The decorrelated vector isˆym =Q −1/2m ˆy m =Q −1/2m A (ϑ)w m +σ2Q −1/2m e m +Q −1/2m εm ,(16)where Q m is given by (12),and the estimation error of ˆy m satisfies E Q −1/2m εm Q −1/2m εm H=I M .The component corresponding to the k th signal in the covariancevector is Q −1/2m a (ϑk )ηk ξ∗k,m ,and the stochastic perturbationis Q −1/2m εm ,thus the ASNR of the signal component with respect to the estimation perturbation isASNR k =η2k M a H (ϑk )Q −1m a (ϑk )=Nη2k MR m,m a H (ϑk )R−1a (ϑk ).(17)The ratio of the ASNR of the covariance vector to that ofthe raw array output can then be calculated based on (14)and (17)as follows,ρ=ASNR kASNR k =Nηk σ2R m,m a H (ϑk )R −1a (ϑk )M.(18)When multiple signals impinge simultaneously,it is difficult to give the explicit value of the ratio directly from (18),but the proportionality of it to the snapshot number holds without doubt.That is to say,when moderately many snapshots are collected,the ASNR of the covariance vector will surpass that of the raw array output.In order to make the degree of the ASNR improvement clearer,we further simplify the scenario to the single-signal case.Based on such simplification,one can easily conclude that R =η1a (ϑ1)a H (ϑ1)+σ2I M and R m,m =σ2+η1,thus the ratio given in (18)can be rewritten asρ=N η1σ2(σ2+η1)(σ2+Mη1).(19)As (σ2+η1)(σ2+Mη1)η1σ2≥ √M +1 2and the equivalence holds when σ2η1=√M ,one can conclude from Eq.(19)that,improved ASNR is obtained in the covariance vector when N >√M +1 2.The minimal value of N for such improvement may vary with σ2η1,but it exists for certain due to the proportionality between ρand N .During the above analysis,we have taken the independent narrowband scenario for example for convenience.But as the magnitudes of the estimation errors of both narrowband and wideband covariance vectors are inversely proportional to √N ,which is indicated by (12),similar conclusions on the ASNR improvement for different kinds of signals can be obtained as that for independent narrowband ones.IV.C OVARIANCE V ECTOR -B ASED DOA E STIMATION In this part,we propose the SMV CV-RVM method to estimate the directions of independent narrowband/widebandsignals by reconstructing ˆy,and propose the MMV CV-RVM method to estimate the directions of correlated narrowbandsignals by reconstructing ˆy1,···,ˆy M jointly.A.independent narrowband/wideband DOA estimation In the independent signal scenarios,we remove the noisecomponent from ˆyto make the relationship between the vector and the signal directions clearer.Denote the noise-removedcounterpart of ˆyby ˆz ,i.e.,ˆz=ˆy −σ2˜e =K k =1ηk b (ϑk )+ε.(20)where b (ϑk )=b T k,1,···,b Tk,M T with b k,m defined inthe same way as that in (4).Eq.(20)indicates that ˆzis a weighted summation of the K signal-components of b (ϑk )for k =1,···,K besides the estimation error,thus the signal directions can be estimated if those vectors are recoveredfrom ˆz.In the following,we introduce the SBL technique [15]to recover those vectors and estimate the directions of independent narrowband/wideband signals.In order to recover the signal components from ˆz,we first sample the potential space of the incident signals discretely to yield a direction set Θ=[θ1,···,θI ]and forms the corre-sponding manifold dictionary B (Θ)=[b (θ1),···,b (θI )],with b (θi )constructed according to the formulation of b (ϑk )in (20)by replacing ϑk with θi .For example,sampling the [-90o 90o ]space with interval Δθ=1o forms a set Θ=[−90◦,−89◦···,90◦]and a dictionary accordingly,the di-rection set covers the possible space of the incident signals.Then ˆzcan be rewritten in the following form,ˆz=I i =1¯ηi b (θi )+ε=B (Θ)¯η+ε,(21)where ¯η=[¯η1,···,¯ηI ]Tis a zero-padded extension ofη=[η1,···,ηK ]Tfrom ϑto Θ,i.e.,it has nonzero values only for θ∈ϑ.Actually,the discrete set Θis never dense enough to include all the possibilities of the true source directions,but the grid mismatch will not deteriorate the validity of the overcomplete model significantly [4],[11]-[14].In most practical array processing problems,it is satisfied that I M >K ,thus the dictionary B (Θ)is overcomplete and (21)is a sparse model.The SBL technique is introduced to extract the basis set of [b (ϑ1),···,b (ϑK )]from B (Θ)toapproximate ˆzunder a model parsimony constraint.We refer the interested readers to [14]for more detailed explanations for the predominance of the SBL technique in DOA estimation by exploiting the spatial sparsity of the incident signals.Based on the above overcomplete formulation of the co-variance vector,we assume that ¯ηis Gaussian distributed as ¯η∼N (0,Γ),with Γ=diag (γ)and γ=[γ1,···,γI ],then the probability of ˆzwith respect to γis given as p (ˆz ;γ)=p (ˆz |¯η)p (¯η;γ)d ¯η=|πΣˆz |−1exp −ˆz H Σ−1ˆz ˆz × |πΣ¯η|−1exp − (¯η−μ)H Σ−1¯η(¯η−μ) d ¯η=|πΣˆz |−1exp −ˆz H Σ−1ˆzˆz ,(22)whereμ=ΓB H (Θ) Q +B (Θ)ΓB H (Θ)−1ˆz ,(23)Σ¯η=Γ−ΓB H (Θ) Q +B (Θ)ΓB H(Θ)−1B (Θ)Γ,(24)Σˆz =Q +B (Θ)ΓB H(Θ).(25)The coefficient vector of ¯ηin (21)should be non-negative as it represents the spatial power distribution of the incident signals.However,since non-Gaussian assumptions generally greatly block Bayesian parameter estimation processes,we abandon this prior information and follows the guideline ofSBL to append a Gaussian distribution to ¯η.Eq.(22)reveals the relationship between ˆzand γ,and γcan be optimized by maximizing the likelihood function.After that,¯ηis calculated according to (23)and the indexes of its nonzero elements indicate the signal directions.Taking the logarithm of (22)and neglecting the constants yields the following objective function for optimizing γ,L (γ)=ln |Σˆz |+ˆz H Σ−1ˆzˆz .(26)The EM algorithm [23]can then be used to estimate γbyminimizing this objective function.During each EM iteration,the first-and second-order posterior moments of ¯ηare calcu-lated with (23)and (24)in the E-step,and γis updated by minimizing L (γ)according to ∂L (γ)/∂γ=0in the M-step,which results in the following update strategy of γ,γ(q )i = μ(q )i 2+ Σ(q )¯η i,i,(27)where the superscript (•)(q )represents the q th iteration,μ(q )and Σ(q )¯ηare calculated with (23)and (24)in the q th iteration,μi is the i th element of μ,and (Σ¯η)i,i is the (i,i )th element of Σ¯η.The update strategy given in (27)can be substituted with a fixed-point iteration as γ(q )i = μ(q )i 2 1−Σ(q )¯η i,iγ(q −1)i +ςto speed up the con-vergence of the EM algorithm,with ςbeing a small positivevalue [14],[15].The initialization and termination criteria of the EM algorithm are also set similarly as those in [14].It should be noted that,the noise component of σ2˜ecan also be recovered from the covariance vector by taking ˆyas the measurement and combining ˜einto B (Θ)to form a dictionary of [B (Θ),˜e],the reconstruction process follows the same guideline as the one presented above.In the re-construction result,the coefficient of ˜estands for the noise variance estimate.However,as such a model extension com-plicates the sparsity profile,we have concluded from sufficient empirical evidences that the extension deteriorates both the computational efficiency and the reconstruction performance of the method in most of the cases.Therefore,we choose to estimate the noise variance directly from the array output (with the estimation method given in Section VI),and remove thenoise component from ˆybefore the reconstruction process.Since the predefined direction set Θis formed via discrete spatial sampling,notable quantization errors may be intro-duced into the DOA estimates if the peak locations of ˆγare taken as the source directions directly.Therefore,we introduce a refined scanning process similar as that in [14]to improve the DOA estimation precision based on the reconstruction result.Denote the estimates of γand Σˆz when the iterative reconstruction process is terminated by γ#and Σ#ˆz,respec-tively,the direction sets in Θcorresponding to the spec-tral lines associated with each signal by θ1,···,θK ,the hyper-parameter set associated with θk by γk ,and define Θ−k =Θ\θk ,i.e.,removing θk from Θyields Θ−k ,γ−k =γ#\γk ,Γ−k =diag (γ−k ),Γk =diag (γk ),and Σ−k =Q +B (Θ−k )Γ−k B H (Θ−k ),then Σ−k can be deemed as the covariance matrix of the estimation error andthe other K −1signal components in ˆzexcept the k th one.In order to obtain refined DOA estimates for the incident signals,we use a single spectral line of βk b (θ)b H (θ)to substitute the spectral peak of B (θk )Γk B H (θk )to denote the k th signal component in the covariance matrix.By introducing the new formulation into (26),one can obtain the following objective function for estimating ϑk ,L (βk ,θ)=ln Σ−k +βk b (θ)b H (θ)+ˆz H Σ−k +βk b (θ)b H (θ) −1ˆz .(28)The estimate of βk can be derived according to∂L (βk ,θ)/∂βk =0as follows,ˆβk =b H (θ)Σ−1−k ˆz ˆz H −Σ−k Σ−1−k b (θ) b H (θ)Σ−1−k b (θ) 2.(29)Then substituting (29)into ∂L (βk ,θ)/∂θ=0yields thefollowing equality,g (θ)Δ=Re b (θ)H Σ−1−k b (θ)b H (θ)Σ−1−k ˆzˆz H −ˆz ˆz H Σ−1−k b (θ)b H(θ) Σ−1−k d [b (θ)]dθ =0.(30)where b (θ)is constructed according to the formulation of b (ϑk )in (20)by replacing ϑk with θ.In practice,this equality may not hold due to measurement perturbation or inaccurate signal reconstruction,and the refined DOA estimate of the k th signal should be obtained via 1-D scanning by checking the distance of g (θ)from 0,i.e.,ˆϑk =arg max θ∈Ωk|g (θ)|−1,(31)where Ωk represents the peak scope of the k th signal.B.correlated narrowband DOA estimationIn the case of correlated narrowband signals,the signal-components have different weights in y 1,···,y M accordingto (6),thus they cannot be reconstructed uniformly in the same way as that for independent signals.However,as they share identical bases of a (ϑ1),···,a (ϑK ),a joint reconstruction procedure can be introduced to recover the signal-components in those vectors and estimate the signal directions.When one reconstructs the covariance vectors ofˆy1,···,ˆy M jointly for DOA estimation,their estimation errors should be dealt with carefully.That is because the second-order statistic given in (13)indicates that theestimation errors of ˆy1,···,ˆy M are correlated with each other,thus they should be decorrelated as a whole before the reconstruction procedure.We use Q to decorrelateˆy= ˆy T1,···,ˆy T M T after removing the noise component to obtain the following covariance vector,ˆz =Q −1/2ˆy −σ2˜e =z +Q −1/2ε,(32)where z =Q −1/2vec (R )is the perturbation-free counterpartof ˆz,and the decorrelated perturbation satisfies E Q −1/2ε Q −1/2ε H =I M 2,(33)which means that the M M ×1subvectors of Q −1/2εare independent to each other.Moreover,as Q −1/2=√N R −1/2 T⊗R −1/2,the m th M ×1subvector of z can be formulated aszm=√N M i =1c m,i R −1/2 y i −σ2e i =√N R −1/2A (ϑ)M i =1c m,i u m ,(34)where c m,i is the (m,i )th element ofR −1/2 T .Eq.(34)indicates that the zm’s share identical sparsity profiles when decomposed on the overcomplete dictionary of R −1/2A (Θ),with A (Θ)=[a (θ1),···,a (θI )]defined similarly as B (Θ)in (21),thus the directions of narrow-band correlated signals can be estimated by reconstructingˆz 1,···,ˆz M jointly.Denote u m =M i =1c mi u m ,A (Θ)=√N R −1/2A (Θ),¯u m is the zero-padded extension of u mfrom ϑto Θ,and assume ¯u m∼N (0,Γ),then the update strategy for estimating γcan be derived similarly as that given in (27)as follows,γ(q )i = M (q )i • 22M + Σ(q )Mi,i,(35)where M (q )and Σ(q )M are the first-and second-order posterior moments of ¯U =[¯u 1,···,¯u M ],with M (q )=Γ(q −1)A (Θ)H Ψ−1ˆZ and Σ(q )M=Γ(q −1)−Γ(q −1)A (Θ)H Ψ−1A (Θ)Γ(q −1),in which Ψ=I M +A (Θ)Γ(q −1)A (Θ)H ,ˆZ =[ˆz 1,···,ˆz M ]and ˆzm is the m th M ×1subvector of ˆz =Q −1/2 ˆy −σ2˜e .When the predefined terminating criterion of the EM algo-rithm is satisfied,the joint reconstruction result can be used to realize refined DOA estimation according to (31)with g (θ)given asg (θ)=Re a(θ)H Σ−1−ka (θ)a (θ)H Σ−1−k ˆZ ˆZ H −ˆZ ˆZ H Σ−1−k a (θ)a (θ)H Σ−1−k d [a(θ)]dθ ,(36)where a(θ)=√N R −1/2a (θ)and Σ−kis defined similarlyas its counterpart in (31).It should be noted that,when all the incident narrowband signals are coherent to each other,the matrices of R and Q are nearly singular when the SNR is high,thus the decorrelation process in (32)becomes unstable.In such scenarios,the matrix inverse lemma should be applied when calculating the inverse of R and Q .Take R for example,its inverse can be calculated asR −1=σ2I M −σ2tr (R )R −σ2I M ,(37)。
