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

Sparse Bayesian learning (SBL) has been successfully applied in solving the problem of direction-of-arrival (DOA) estimation. However, SBL needs multiple snapshots to ensure accuracy and costs huge computational workload. To reduce the requirement of snapshot and computational burden, a DOA estimation method based on the randomize-then-optimize (RTO) algorithm is first time introduced. RTO algorithm uses the optimization and Metropolis-Hastings approach to avoid the “learning” process of SBL in updating hyperparameters. And in order to apply RTO algorithm in the circumstance of signal with Laplace prior, a prior transformation technique is first induced. Compared with conventional CS based DOA methods, the proposed method has a better accuracy with single snapshot and shorter processing time. Some simulations are conducted to demonstrate the effectiveness of the proposed method.


Introduction
The Direction-of-Arrival (DOA) estimation problem is an important research area in radar and antenna which mainly concerns on recovering the signals and getting the impinging angles.In recent years, compressive sensing (CS) [1] theory has been successfully applied in DOA estimation.Its key idea is to reconstruct sparse signal.
Based on CS theory, sparse signal reconstruction method is extended to the Bayesian compressive sensing (BCS) [5] which is formulated from a Bayesian perspective based on the sparse prior assumption of signal and noise.
Under the frame of BCS, signal reconstruction is mainly achieved by sparse Bayesian learning (SBL) [6].The SBL has obtained great development in the recent researches, such as root SBL [7], variational SBL [8], jointly SBL [9] and grid evolution method [15].SBL has a multilayered assumption frame which is designed to iterate to "learn" new information and update the hyperparameters.
In the practical applications of DOA estimation problems, the snapshot is commonly short or even single because of the need of real-time processing.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 [17] and the spatial smoothing technique [18] are widely used.The influence of snapshot on estimation accuracy still exists in sparse signal reconstruction methods.Usually for SBL based DOA methods, the accuracy of estimated DOAs strongly depends on the number of snapshots.It causes poor performance if the number of snapshots is not large enough for retrieving the DOAs.This phenomenon also occurs in the l1-SVD method as it requires sufficient data for the singular value decomposition and dimensionality reduction.The studies on single snapshots for sparse signal reconstruction are limited.
Thus, aim to relax the requirement of snapshots and improve the performance of sparse signal reconstruction for single snapshot, this letter proposes a novel DOA estimation method based on the Randomize-Then-Optimize (RTO) [10] 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 [13].This MH method is based on Markov chain Monte Carlo (MCMC) [14], which has been widely used to evaluate the posterior distribution in the Bayesian inverse problem.Different from traditional SBL, the "learning" 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" the hyperparameters in the assumption frame.The RTO has no more requirement on the hyperparameters so that the processing time is significantly reduced.Furthermore, the RTO approach proposed in [10] is based on Gaussian prior, while the Laplace prior has better parametric sparsity for DOA estimation problem [11].Therefore, a prior transformation method [12] is adopted.Simulation results show that good accuracy can be achieved even with one single snapshot.
( ) ( ) is the array manifold matrix.( ) The model given in Eq. ( 1) can be regarded as a Bayesian inverse problem.The manifold matrix () θ A is a nonlinear parameter-to-measurement mapping.Commonly, solving Bayesian inverse problem mainly focuses on characterising the posterior density.For traditional SBL method, the process of obtaining the posterior density requires the hyperparameters γ.According to the literature [16], solving the Bayesian inverse problem leads to minimize the cost function of γ.The likelihood that SBL will converge to the global minimum of cost function is increased along with the increasing of snapshot number.This is important because, the maximally sparse representation is guaranteed due to the globally minimizing SBL hyperparameters and increasing snapshot number improves the probability that these hyperparameters are found.Too few snapshots number produces poor performance of SBL.
From another perspective of Bayesian compressive sensing, RTO method is freed from the dependency of snapshots 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 will be used in Metropolis-Hastings (MH) approach as a Metropolis independence proposal.

Randomize-Then-Optimize-Metropolis-Hastings (RTO-MH) Algorithm
Considering the Bayesian inverse problem with Gaussian measurement errors and the Gaussian prior, linear transformations are used to "whiten" the prior and the error model so that the inverse problem has the following model, where θ0 is the prior mean, In and Im are identity matrices with size n and m, respectively.And y is the measurement vector, f is the forwarding function (also called as the "parameter-to-observable mapping") and ε is the measurement error.
RTO obtains candidate samples by repeatedly optimizing a randomly perturbed cost function.In order to obtain the unknown parameter θ, RTO requires target distribution to have a specific form.This distribution allows the RTO samples being used in the MH process.The MH approach corrects these samples which contain the information of spatial spectrum.Especially, it needs the target density (usually is the posterior density of θ) to be written as where ||•|| denotes the norm.The right side in the symbol of norm is defined as a vector-valued function of the parameters θ which has the form of ( ) ( ) We illustrate the steps of sampling from the posterior as Eq. ( 3) by RTO.Firstly, find a linearization point, denoted as θ , and fix it.Usually, this point is set to be posterior mode.The following equation is used to obtain the posterior mode, ( ) Secondly, calculate the Jacobian of F(θ) at the linearization point, denoted as ( ) The orthonormal basis, denoted as Q for the column space of ( ) F θ J , is evaluated through a thin QR factorization of ( ) Thirdly, compute independent samples ξ according to an n-dimensional standard Gaussian distribution.The proposal points prop θ are found by solving the following optimization problem, ( ) According to the previous researches 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 as an independence proposal in the Metropolis-Hastings approach.For the update from a point ( 1)   i θ − to the proposed point ( )   i prop θ , the Metropolis-Hastings acceptance ratio is where w(θ) are The process mentioned above is called as RTO-Metropolis-Hastings (RTO-MH) algorithm which combines RTO algorithm and MCMC method [10].The samples obtained by using RTO can be used within this framework of the Metropolis-Hastings (MH) algorithm.This process produces the samples from a posterior based on an arbitrary measurement model.As described in Eq. ( 2), the RTO-MH is proceeded with a Gaussian prior.
However, in the BCS based DOA estimation methods, a Laplace prior is preferred.From this consideration, a prior transformation technique is induced in RTO-MH in this letter.

RTO-MH With A prior Transformation
Firstly, consider the single parameter with Laplace prior of the model in Eq. ( 2) Then the posterior distribution is where obs σ is the error standard deviation.
To satisfy the posterior form in Eq. ( 3), an invertible mapping function T1D which connects a Gaussian reference random variable v with 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 and L is the cdf of the Laplace distribution.
Based on the single parameter prior transformation, extend this technique to the multiple parameters case.
Assume the prior on θ is where (Dθ)i is the ith element of vector Dθ and D is an invertible matrix.Thus, the posterior can be derived ) Reference random Guassian distributed variables can be converted to each Laplace-distributed element of Dθ by means of the one-dimensional transformation T1D defined in Eq. (11).So Dθ = T(v), where : ,, where vi is the corresponding each reference variables, resulting in the prior transformation, ( ) The Jacobian of the transformation is D -1 JT, where JT is the Jacobian of T noted as is the derivative of Eq. ( 11), ( ) where ( ) is the pdf (probability density function) of the standard Gaussian distribution.
Based on the transformation step as in Eq. ( 15), the posterior density of v can be derived as

Implementation
As the observation model stated in Eq. ( 1), the measurement vector is complex.However, the RTO-MH with prior Equation ( 19) still satisfies the Bayesian inverse problem model as Eq. ( 2) for the corresponding parameters.
Therefore, the DOAs can be retrieved by 3 steps.Firstly, make the measurement matrix transformation from complex to real.Secondly, input the real measurement matrix to the RTO-MH with prior transformation.Finally, calculate the normalized power of the results and select the main lobes as DOAs.
The RTO-MH with prior transformation algorithm for DOA estimation is summarized in Table 1.
Table 1.DOA Estimation Algorithm Based On RTO

DOA estimation algorithm based on RTO-MH with prior transformation
1. Do an observation model data transformation from complex value to real value using Eq.(19).
2. Use the prior transformation function as Eq. ( 11) so that ( ) vT has a standard Gaussian distribution.
3. Find the mode v of transformed posterior density ( ) p vy in the form of Eq. ( 18) by Eq. (4).
4. Evaluate Q through a thin QR factorization of ( ) F θ J as defined after Eq. (4).

Results and Discussion
Firstly, the performance of the SBL method and the proposed method are compared.In the simulations, single 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 The spatial spectrum of the BCS and the proposed method.
Figure 1 shows the simulations results of the proposed method with single snapshot and the results of the SBL method with single and 200 snapshots, respectively.The vertical dotted lines indicate the true locations of DOAs.For single snapshot, the SBL method can only retrieve the DOA at 30 o .And its sidelobes are so high, implying the poor performance of the SBL for single snapshot.When the snapshots go to 200, the performance of the SBL improves.The proposed method and SBL method both retrieve the signal at 30 o exactly and have an 1 o deviation at 60 o .Though the main lobe at 60 o of the proposed method is wider that of SBL method, the side lobes of the proposed method are much lower.
The SBL method is simulated for 1 snapshot and 200 snapshots respectively.The l1-SVD algorithm is simulated for 1 snapshot and 50 snapshots, respectively.Fig. 2 The RMSE of OMP, l1-SVD and the proposed method versus SNR.
Figure 2 shows that the proposed method has an obvious lower RMSE level than OMP, SBL with 1 snapshot and l1-SVD with 1 snapshot algorithms.For the SBL algorithm with 200 snapshots, it has an almost same accuracy with the proposed method.This indicates that the proposed method has advantage over snapshot number comparing with SBL algorithm.For the l1-SVD with 50 snapshots, its accuracy exceeds all the mentioned methods, including the proposed method.TABLE 2 gives a comparison on the processing times of the above methods at SNR=10 dB.It can be found that the proposed method costs the shortest time compared with SBL and l1-SVD.OMP method has a similar time cost.However, as shown in Fig. 2, the proposed method has better accuracy than O M P.

Conclusion
In this letter, the RTO-MH with prior transformation algorithm is proposed to improve the performance of BCS method in the DOA estimation.Compared with conventional BCS methods, such as SBL, whose accuracy highly depends on the snapshot numbers, the proposed method doesn't require many samples to update the hyperparameters.Simulation results demonstrate that the proposed method has better estimation accuracy than SBL, OMP and l1-SVD when the number of snapshots is single and the processing time is reduced effectively.Its computational burden is close to OMP algorithm and lower than l1-SVD and SBL methods.However, a shortcoming of the proposed method is the limited resolution.The resolution limit is 5 o based on several experiments.In the future work, this proposed method will be extended to off-grid signals under the circumstance of wideband sources and its resolution limit will be reduced to 1 o .
The spatial spectrum of the BCS and the proposed method.
The RMSE of OMP, l1-SVD and the proposed method versus SNR.
vector of the kth source whose entry ( ) mK θ a contains the delay information of the kth source to the mth sensor.To locate the direction of the sparse signal, the spatial domain is sampled by a grid N denotes the grid number and usually N>>M>K.The sampling is usually uniform because there is no prior information of sources.The measurement vector () t y and array manifold matrix () θ A are known.() t s needs to be estimated.

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

3 Experimental
Some simulations are conducted to verify the effectiveness of the proposed method.Two uncorrelated sources are from 30 o and 60 o , respectively.The signal sources are narrowband and impinge on a uniform linear array with 16 sensors.The spatial range of interested signal is discretized by means of a uniform grid sample in range [-90 o ,90 o ]with the grid interval of 1 o .Also, to illustrate how the accuracy performs, the proposed method is compared with O M P, l1-SVD and SBL in RMSE.The signal-to-noise ratio (SNR) in the comparison is [-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10].200 independent Monte-Carlo trials for each SNR are implemented.Finally, to state the superiority of the proposed method over processing time, the processing times of the methods mentioned above are measured at SNR=10dB.All experiments are carried out in MATLAB on a desktop with a Windows 10 system and eight 2.2GHz CPUs.
transformation can only operate for real values.In order to apply the algorithm in DOA estimation, the observational model is transformed into real form, is the real noise matrix.Φ is the real manifold matrix which has the form

Table 2 .
Processing Time Comparison