Direction-of-ArrivaI Method Based on Randomize-Then-Optimize Approach

Cai-Yi Tang | Sheng Peng | Zhi-Qin Zhao | Bo Jiang

Abstract—The direction-of-arrival (DOA) estimation problem can be solved by the methods based on sparse Bayesian learning (SBL).To assure the accuracy,SBL needs massive amounts of snapshots which may lead to a huge computational workload.In order to reduce the snapshot number and computational complexity,a randomizethen-optimize (RTO) algorithm based DOA estimation method is proposed.The “l(fā)earning” process for updating hyperparameters in SBL can be avoided by using the optimization and Metropolis-Hastings process in the RTO algorithm.To apply the RTO algorithm for a Laplace prior,a prior transformation technique is induced.To demonstrate the effectiveness of the proposed method,several simulations are proceeded,which verifies that the proposed method has better accuracy with 1 snapshot and shorter processing time than conventional compressive sensing (CS) based DOA methods.

1.Introduction

The direction-of-arrival (DOA) estimation problem is an important research area in radar and antenna,which mainly concerns on recovering signals and getting impinging angles.In recent years,the compressive sensing (CS)[1]theory has been successfully applied in DOA estimation,aiming to reconstruct sparse signals.Several sparse signal construction algorithms have been developed,such as orthogonal matching pursuit(OMP) and basis pursuit (BP).OMP[2]is originated from matching pursuit (MP)[3].To reduce the computational burden of BP methods,[4] proposed a dimensionality reduction method,named asl1-SVD.Based on the CS theory,the sparse signal reconstruction method was extended to Bayesian compressive sensing (BCS)[5]which was formulated from a Bayesian perspective based on the sparse prior assumption of signal and noise.Under the frame of BCS,signal reconstruction was mainly achieved by sparse Bayesian learning (SBL)[6].SBL has obtained great developments in recent studies,such as root SBL[7],variational SBL[8],jointly SBL[9],and the grid evolution method[10].SBL has a multilayered assumption frame which is designed to iterate to“l(fā)earn” new information and update the hyperparameters.

In the practical applications of DOA estimation problems,the snapshot is usually short or even single.However,the limited data results in poor accuracy of DOA estimation problems.To lower the demand of snapshots,most of the existed methods are developed to improve the traditional spatial spectrum estimation methods.The pseudo covariance matrix method[11]and the spatial smoothing technique[12]are widely used.The influence of snapshots on estimation accuracy still exists in sparse signal reconstruction methods.Usually for SBL based DOA methods,the accuracy of estimated DOA strongly depends on the number of snapshots.It causes poor performance if the number of snapshots is not large enough for retrieving DOA.This phenomenon also occurs in thel1-SVD method as it requires sufficient data for the singular value decomposition and dimensionality reduction.To sum up,the studies on single snapshot for sparse signal reconstruction are limited.

Thus,aim to relax the requirement of snapshots and improve the performance of sparse signal reconstruction using 1 snapshot,this paper proposes a DOA estimation method based on the randomizethen-optimize (RTO)[13]approach.The RTO based method solves the Bayesian nonlinear inverse problem by using a process of optimization to generate the proposal samples and correcting these samples by the Metropolis-Hastings (MH) approach[14].This MH method is based on the Markov chain Monte-Carlo(MCMC)[15],which has been widely used to evaluate the posterior distribution in the Bayesian inverse problem.Different from traditional SBL,the “l(fā)earning” process in updating hyperparameters is no longer needed in this RTO based method,reducing the required number of snapshots.Also,an intrinsic shortcoming of SBL is that it consumes relative long time to “train” hyperparameters in the assumption frame.RTO has no more requirement on hyperparameters,so that the processing time is significantly reduced.Furthermore,the RTO approach[13]is based on a Gaussian prior,while a Laplace prior has better parametric sparsity for the DOA estimation problem[16].Therefore,a prior transformation method[17]is adopted.Simulation results show that good accuracy can be achieved even with 1 snapshot.

2.Proposed Method

2.1.ProbIem FormuIation

Assume that there areKnarrowband far field signalssk(t),wherek=1,2,···,K,impinging on a linear array ofMsensors fromKdirectionsθk,wherek=1,2,···,K.In the DOA estimation problem,the time delays at different sensors can be represented by simple phase shifts,leading to an observation model of

where y(t)=[y1(t),y2(t),…,yM(t)]T,θ=[θ1,θ2,…,θK]T,s(t)=[s1(t),s2(t),…,sK(t)]T,and e(t)=[e1(t),e2(t),…,eM(t)]T.ym(t) andem(t) (m=1,2,···,M) are the measurement and noise of themth sensor at timet,respectively.And A(θ)=[a(θ1),a(θ2),…,a(θK)] is the array manifold matrix,a(θk) is the steering vector of thekth source whose entryam(θK) contains the delay information of thekth source to themth sensor.To locate the direction of the sparse signal,the spatial domain is sampled by a grid,whereNdenotes the grid number and usuallyN＞＞M＞K.The sampling is usually uniform because there is no prior information of sources.The measurement vector y(t) and array manifold matrix A () are known,but s (t) needs to be estimated.

The model given in (1) can be regarded as a Bayesian inverse problem.The manifold matrix A(θ) is a nonlinear parameter-to-measurement mapping.Commonly,solving the Bayesian inverse problem mainly focuses on characterising the posterior density.For the traditional SBL method,the process of obtaining the posterior density requires the hyperparameters γ.According to [18],solving the Bayesian inverse problem leads to minimize the cost function of γ.The likelihood that SBL will converge to the global minimum of the cost function is increased along with increasing the snapshot number.This is important,because the maximally sparse representation is guaranteed due to globally minimizing SBL hyperparameters,and increasing the snapshot number improves the probability that these hyperparameters are found.Too few snapshots induce poor performance of SBL.

From another perspective of BCS,the RTO method is freed from the dependency of the snapshot number.It solves the Bayesian inverse problem with limited snapshots.RTO produces samples by using repeated solutions of a randomly perturbed optimization problem from a proposal density,which can be used in the MH approach as a Metropolis independence proposal.

2.2.RTO-MH AIgorithm

Considering the Bayesian inverse problem with Gaussian measurement errors and the Gaussian prior,linear transformations are used to “whiten” the error model and the prior,so that the following model can be used to describe the inverse problem:

where ε～N(0,Ii),θ～N(θ0,Ij)θ0is the prior mean,and Iiand Ijare identity matrices with sizeiand sizej,respectively.The model in (2) is only used to describe the Bayesian inverse problem.And y is the measurement vector,fis the forwarding function (also called as the “parameter-to-observable mapping”),and ε is the measurement error.

By repeatedly optimizing a randomly perturbed cost function,RTO can obtain candidate samples.In order to obtain the unknown parameter θ,RTO requires the target distribution to have a specific form which is defined by (3).This distribution allows the RTO samples being used in the MH process.It is suitable for any form of the target distribution scenarios,because the model in (2) can be generally used in any inverse problem.The MH approach corrects these samples which contain the information of the spatial spectrum.Especially,it needs the target density (usually the posterior density of θ) to be in the form of

where ||·|| denotes the norm.The right side in the symbol of the norm is defined as a vector-valued function of the parameter θ which has the form of

We illustrate the steps of sampling from the posterior as (3) by RTO.

1) Find a linearization point,denoted as,and fix it.Usually,this point is set to be the posterior mode.The following equation is used to obtain the posterior mode:

2) Calculate the Jacobian ofF(θ) at the linearization point,denoted as JF().The orthonormal basis,denoted asfor the column space of JF(),is evaluated through a thin quadrature rectangle (QR)factorization approach of JF().

3) Compute independent samples ξ according to ann-dimensional standard Gaussian distribution.The proposal points θpropare evaluated by tackling the following optimization problem:

According to the previous study in [10],the proposal points are distributed in terms of the proposal density:

where |·| denotes the absolute value of the matrix determinant.This proposal density is used in the MH approach as an independence proposal.For updating a pointθ(i-1)to the proposed point,the MH acceptance ratio is

The process mentioned above is called as the RTO-MH algorithm which combines the RTO algorithm and the MCMC method[10].The framework of the MH algorithm can use the samples which are obtained by applying the RTO algorithm.This process produces the samples from a posterior based on an arbitrary measurement model.As described in (2),RTO-MH is proceeded with a Gaussian prior.However,in BCS based DOA estimation methods,a Laplace prior is preferred.From this consideration,a prior transformation technique is induced in RTO-MH in this paper.

2.3.RTO-MH with Prior Transformation

First,consider the single parameter with a Laplace prior of the model in (2):

whereλis a hyperparameter used to describe the Laplace distribution and θ is a Laplace-distributed physical parameter.

Then the posterior distribution is

whereσobsis the error standard deviation.

To satisfy the posterior form in (3),an invertible mapping functionT1Dwhich connects a Gaussian reference random variablevwith the Laplace-distributed physical parameter θ,such that θ=T1D(v) is constructed.The transformation equation can be written as

whereφis the cumulative distribution function (CDF) of the standard Gaussian distribution,andLis CDF of the Laplace distribution.

Based on the single parameter prior transformation,this technique is extended to the multiple parameters case.Assume the prior on θ is

where (Dθ)iis theith element of the vector Dθ andD is an invertible matrix.Thus,the posterior can be derived:

Random Gaussian distributed variables can be converted to each Laplace-distributed element of Dθ by means of the one-dimensional transformationT1Ddefined in (11).So Dθ=T(v),where

whereviis the corresponding reference variable,resulting in the prior transformation:

The Jacobian of the transformation is D-1JT,where JTis the Jacobian of T noted as

Based on the transformation step as in (15),the posterior density of v can be derived as

2.4.ImpIementation

As the observation model stated in (1),the measurement vector is complex.However,RTO-MH with the prior transformation can only operate for real values.In order to apply the algorithm in DOA estimation,the observational model is transformed into the real form

where tm=[Re{ym},Im{ym}]Tis the real measurement matrix,wm=[Re{si},Im{si}]Tis the real signal matrix,and nm=[Re{em},Im{em}]Tis the real noise matrix.Φ is the real manifold matrix which has the form:

Equation (19) still satisfies the Bayesian inverse problem model as (2) for the corresponding parameters.Therefore,DOA can be retrieved by three steps: 1) Make the measurement matrix transformation from complex to real;2) input the real measurement matrix to RTO-MH with prior transformation;3) calculate the normalized power of the results and select the main lobes as DOA.

Algorithm 1.

RTO-MH with the prior transformation algorithm for DOA estimation is:

1.Do an observation model data transformation from the complex value to the real value using (19).

2.Use the prior transformation function as (11),so that v=T-1(θ) has a standard Gaussian distribution.

5.Fori=1,2,…,nsamps

Draw a standard Gaussian sampleξ(i)～N(0,In).

Compute RTO samples as

6.Fori=1,2,…,nsamps

Sample v from a uniform distribution on [0,1].

7.Fori=1,2,…,nsamps

Defineθ(i)=T(v(i)),which are the corrected samples fromp(θ ∣y).

End for

8.The resultedθ(i)forms the spectrum of the spatial domain where the interested signals exist.Calculate the normalized power of the spectrum and select the main lobes as DOA.

3.SimuIations

Several simulations are conducted to verify the effectiveness of the proposed method.Two uncorrelated sources are from 30oand 60o,respectively.The signal sources are narrowband and impinge on a uniform linear array with 16 sensors.The spatial range of the interested signal is discretized by the means of a uniform grid sample in the range of [-90o,90o] with a grid interval of 1o.First,the performance of the SBL method and that of the proposed method are compared.In the simulations,1 snapshot is used for the proposed RTO based method and the SBL method.Also,the SBL method with 200 snapshots is performed for comparison.

Fig.1.Spatial spectra of BCS and the proposed method.

Fig.1 shows the simulation results of the proposed method with 1 snapshot and the SBL method with 1 and 200 snapshots,respectively.The vertical dotted lines indicate the true locations of DOA.For 1 snapshot,the SBL method can only retrieve DOA at 30o.And its side lobes are so high,implying the poor performance of SBL for 1 snapshot.When the number of snapshots increases to 200,the performance of SBL improves.Both the proposed method and the SBL method retrieve the signal at 30oexactly and have an 1odeviation at 60o.Though the main lobe at 60oof the proposed method is wider than that of the SBL method,the side lobes of the proposed method are much lower.

To illustrate how the accuracy performs,the proposed method is compared with OMP,l1-SVD,and SBL in the root mean square error (RMSE).The signal-to-noise ratio (SNR) in the comparison ranges from -10 dB to 10 dB with an interval of 2 dB.Here 200 independent Monte-Carlo trials for each SNR are implemented.The SBL method is simulated for 1 snapshot and 200 snapshots,respectively.Thel1-SVD algorithm is simulated for 1 snapshot and 50 snapshots,respectively.RMSE is used to describe the estimation accuracy,defined by

wherenis the number of Monte-Carlo simulations,is the estimated value for each simulation,andis the true value.

Fig.2 shows that the proposed method has an obvious lower RMSE level than OMP,SBL with 1 snapshot,andl1-SVD with 1 snapshot algorithms.The OMP algorithm always keeps a high RMSE level which means its estimation accuracy is much worse than SBL andl1-SVD with many snapshots.For the SBL algorithm with 200 snapshots,it has almost the same accuracy with the proposed method.This indicates that the proposed method has an advantage over the snapshot number comparing with the SBL algorithm.Forl1-SVD with 50 snapshots,its accuracy exceeds all the mentioned methods,including the proposed method.However,RMSE can be regarded as an enlarged value for the true error.Thus,the proposed method almost has the same estimation accuracy asl1-SVD with 50 snapshots.

Table 1 gives a comparison on the processing time of the above methods with SNR=10 dB.All experiments are carried out in MATLAB on a desktop with a 2.2-GHz central processing unit (CPU).It can be found that the proposed method costs the shortest time compared with SBL andl1-SVD.Although the OMP method has a similar time cost,the proposed method has better accuracy than OMP as shown in Fig.2.

Fig.2.RMSE of OMP,l1-SVD,and the proposed method(RTO-MH) versus SNR.

Table 1: Processing time comparison

4.ConcIusions

In this paper,RTO-MH with the prior transformation algorithm is proposed to improve the performance of the BCS method in DOA estimation.Compared with conventional BCS methods,such as SBL,whose accuracy highly depends on the snapshot number,the proposed method does not require many samples to update the hyperparameters.Simulation results demonstrate that the proposed method has better estimation accuracy than SBL,OMP,andl1-SVD,when the number of snapshots is 1,and the processing time is reduced effectively.Its computational burden is close to that of the OMP algorithm and lower than those ofl1-SVD and SBL methods.However,there exists a disadvantage that the proposed method shows limited resolution of 5obased on several experiments.In the future work,the proposed method will be extended to off-grid signals under the circumstance of wideband sources,and its resolution limit can be expected to be reduced to 1o.

DiscIosures

The authors declare no conflicts of interest.