DOA estimation for monostatic MIMO radar using enhanced sparse Bayesian learning

This study discusses the problem of direction-of-arrival estimation (DOA) estimation for a monostatic multiple-input multipleoutput (MIMO) radar system, and a novel sparse Bayesian learning (SBL) framework is presented. To lower the computational load, the matched array data is firstly compressed via reduced-dimension transformation. Then the problem of DOA estimation is linked to a sparse inverse problem. Finally, a forgotten factor-based root SBL algorithm is derived from hyperparameters learning, which can solve the offgrid problem by finding the roots of a polynomial. The proposed algorithm does not require the prior of the source number, and it can apply to the scenario with a small snapshot as well as coarse grid, thus it has a blind and robust characteristic. Numerical simulations verify the effectiveness of the proposed algorithm.


Introduction
The topic of multiple-input multiple-output (MIMO) radar system has aroused extensive attention in the past decade [1][2][3][4][5].MIMO radar is characterised by the features that emit mutual orthogonal waveforms with multiple antennas and receive the echoes with multiple antennas.The spatial diversity enables a MIMO radar to achieve an enhanced performance in parameter estimation, interference and jamming suppression and so on.Generally, the MIMO radar can be divided into statistical MIMO and colocated MIMO radar in term of its antenna configuration [6].The former utilises distribute antennas to combat radar cross section (RCS) coefficient fluctuation.The later explore colocated antenna array to obtain higher-resolution direction estimation.In this paper, we focus on the uniform linear array (ULA)-based monostatic MIMO radar, which belong to the later.
Direction estimation plays an important role in a colocated MIMO radar.Many excellent algorithms have been developed, such as multiple signal classification (MUSIC) [7], estimation method of signal parameters via rotational invariance techniques (ESPRIT) [8], propagator method (PM), tensor-based methods [9-14] and so on.Unlike the angle estimation problem in a bistatic MIMO radar, only one-dimensional angle needs to be determined in a monostatic scenario.Although the estimation algorithms in a bistatic configuration are also suitable for a monostatic case, they are usually too complicating to implement.To lower the computation burden, several reduced-complexity (RC) strategies have been presented, e.g.RC-Capon [15], RC-MUSIC [16], RC-ESPRIT [17] and RC-root-MUSIC [18].It has been shown that typical subspacebased methods, such as RC-MUSIC and RC-ESPRIT, provide accurate estimation performance.However, all the methods mentioned above require the prior information of source number.Moreover, they may fail to work with the small snapshot scenario, especially when the snapshot number is smaller than source number.The optimisation-based methods can overcome such drawbacks [19], among which the sparse Bayesian learning (SBL) methods are more attractive from the perspective of estimation accuracy [20].SBL was first proposed by Tipping at the beginning of this century [21].It is shown that the global minimum of SBL is always the sparest solution to the inverse problem [22].Besides, SBL has much fewer local minimum than many existing algorithms.Usually, a fixed grid is set in SBL, and the DOA estimation problem is resolved via the maximum-likelihood criterion or the maximum-posterior probability principle.Unfortunately, it is usually very hard to choose a proper grid, as a coarse grid may lead to imprecise DOA estimation while a refined grid will bring extensive computational complexity.
Recently, the off-grid DOA estimators are developed to balance the above tradeoff, in which a coarse grid is initiated as the base matrix.The mismatch between the adopted grid and the real direction is considered as a perturbation in the sparse model.The grid can be adaptively updated along with algorithm iteration.The off-grid DOA estimation problem was first stressed in [23], where the perturbation is approximated by the first-order Taylor expansion of the real direction at its nearest grid.The off-grid gap is assumed to satisfy a Gaussian prior, and a sparse total least squares (STLS) solver is investigated.Since the off-grid error is within the grid interval, the Gaussian prior is replaced by the uniform hypothesis, and an off-grid sparse Bayesian inference (OGSBI) algorithm is presented in [24].Moreover, a polynomial rooting-based off-grid sparse Bayesian learning (ROGSBL) scheme is derived in [25], where the grid is updated by finding the roots of a polynomial.It has been shown that the ROGSBL method is more suitable for a coarse grid, and it has more accurate DOA estimation performance than the OGSBI method.Since the rooting processing in ROGSBL is computationally inefficient, not all the grid points are updated in the iteration.Instead, a hard threshold is set to determine the active grid point.Unfortunately, the updating of the grid is still inefficient and can be improved further.
In this paper, an improved DOA estimation algorithm is proposed for a ULA-based monostatic MIMO radar.A reduced-dimension transformation is firstly applied to obtain a Vandermonde structure of the direction matrix as well as lower the computation burden.Then a sparse inverse model is established that links the DOA estimation problem to an off-grid optimisation problem.An improved root-SBL algorithm is derived, in which a forgotten factor model is utilised to update the grid dynamically.Numerical experiments show that the proposed method provides a more accurate DOA estimation than the existing off-grid SBL methods.
The remaining of the paper is organised as follows.The data model for the monostatic MIMO radar is presented in Section 2. The details of the proposed scheme are given in Section 3. Simulation results are illustrated in Section 4. The paper is ended by a brief conclusion in Section 5.
Notation, bold capital letters, e.g.X , and bold lowercase letters, e.g.x denote matrices and vectors, respectively.The M × M identity matrix is denoted by I M .The superscript • ( ) T , • ( ) H and • ( ) −1 stand for the operations of transpose, Hermitian transpose and inverse, respectively; ⊗ and ⊙ represent, respectively, the Kronecker product and the Khatri-Rao product (column-wise Kronecker product).The subscript • F and • 2 represent the Frobenius norm and 2-norm of a matrix; diag • ( ) and vec • ( ) denote the diagonalisation and the vectorisation operations, respectively.tr • ( ) and E • ( ) return the trace and the expectation, respectively.X •n and X m• denote the nth column and the mth row of X , respectively.

Signal model
Consider a monostatic MIMO radar system equipped with M transmitters and N receivers, both of which are colocated ULAs with a half-wavelength spacing, as depicted in Fig. 1.Suppose that there are K far-field sources appearing in the same range bin of the radar system, and the kth (k = 1, 2, . . ., K) DOA is denoted by u k .An additional assumption is that the transmitters emit mutual orthogonal waveforms.The transmitted signals are reflected by the sources, and the echoes are collected by the receivers.The output of the matched filters can be expressed as [16,17] where A T , A R , S and E denote, respectively, the transmit loading matrix, the receive loading matrix, the echo loading matrix and array noise matrix, and where a t (u k ), a r (u k ) and s f k are the kth (k = 1, 2, . . ., K) column vectors of A T , A R and S, respectively; a k , f k and f s stand for the RCS amplitude, the Doppler frequency and the pulse repeat frequency, respectively; L is the number of snapshots, E [ C MN ×L is a white Gaussian noise matrix.From (1), we can observe that the received array data R is redundant, since the dimension of virtual steering vector a t u k ⊗ a r u k is MN, while only R = M + N − 1 unique elements in it.The redundancy in (1) can be expressed as where ã(u k ) = 1, e −jp sin (u k ) , . . ., e −jp(R−1) sin (u k ) T , and However, G cannot be applied on R directly, otherwise, the coloured noise will appear and results in degraded estimation performance.Instead, a new weight matrix which is a full rank matrix.Therefore, multiplying W with R yielding where N = WE is the compressed noise matrix.Obviously, the compressed array data Y is non-redundant, and the associated loading matrix is Ã.
3 Proposed scheme

Sparse representation-based DOA estimation
Generally, sources can be considered to be sparse in the spatial domain.Usually, a fixed grid w 1 , w 2 , . . ., is the virtual loading matrix, and the associate qth (q = 1, 2, . . ., Q) steering vector is b w q = 1, e −jp sin (w q ) , . . ., e −jp(R−1) sin (w q ) T [ C R×1 .Thus the array data Y in (4) can be rewritten as Clearly, X is a row sparse matrix.If we can recovery the support of X , the DOAs are tantamount obtained.However, this is an NP-hard problem.Alternatively, approximate solutions are achieved via solving the following relaxed optimisation problem: where 1 is a tradeoff parameter balancing noise with row sparsity, d X ( ) is a regularisation term encouraging row sparsity.Different strategies pursue various d X ( ), e.g. 5) and ( 6) are established on the assumption that sources are on the grids.Unfortunately, there is always a mismatch between the true DOAs and the fixed grids, which is well known as the off-grid problem.To overcome this drawback, a perturbed sparse representation model is utilised, in which [26] where E is the perturbation to B, and it can be approximated by the first-order Taylor expansion of the off-grid DOAs as d q represents the nearest grid interval with d q = w q − u k if w q ≤ u k ≤ w q+1 and zeros elsewhere (q [ 1, 2, . . ., Q , k [ {1, 2, . . ., K}).F = B + E is the perturbed measurement matrix.Accordingly, the optimisation problem in ( 6) is modified into arg min To obtain the real DOA, one need to find out the support of X as well as the intervals d q (q = 1, 2, . . ., Q).In what follows, an enhanced SBL algorithm is developed for the DOA estimation.

Bayesian model
As a statistical method, SBL is based on the assumption of white Gaussian noise and linear system.In this paper, we assume the compressed noise measurements fulfil a circular symmetric complexity Gaussian distribution with zero mean and unknown variance s 2 .Moreover, we suppose the noise is uncorrelated with the source matrix X .According to the Bayesian principle, we have the following Gaussian likelihood: where w = w 1 , w 2 , . . ., w Q is the grid that determine F. Generally, the rows of X are assigned with an L-dimensional Gaussian prior, i.e.
where g q is a non-negative hyperparameter that controls the row sparsity of X .To obtain a two-stage hierarchical prior to X q• , g q (q = 1, 2, . . ., Q) are modelled with independent Gamma hyperpriors in the empirical Bayesian strategy, which is given by where g = [g 1 , g 2 , . . ., g Q ], r is a small constant.In addition, the unknown hyperparameter s 2 is also modelled with a Gamma hyperprior with Usually, a, b are initialised with a, b 0 to achieve a board hyperprior.

Sparse Bayesian learning
It can be easily concluded from the Bayes rule that the posterior density p X |Y ( ) is Gaussian, i.e.
where the mean m = m 1 , m 2 , . . ., m L and the covariance S are given by where L = diag g .In this paper, the well-known expectationmaximisation (EM) method is explored to find the hyperparameters V = g, s 2 , w , which can be accomplished via maximising p Y ; V ( ).Equivalently, one need to minimise − ln p Y ; V ( ).By treating X as a hidden variable, the EM method try to maximise where V (old) is the estimated hyperparameters in the previous iteration.
To update g, Q V ( ) is firstly simplified by treating s 2 and w as hidden variables.The learning rule for g is therefore obtained by setting the derivative of Q V ( ) with respect to g q to zeros, leading to where J l = m l m H l + S, and • [ ] q,q is to get the qth diagonal element of a matrix.By treating g and w as hidden variables and forcing the gradient of Q V ( ) over s 2 to zero yields to

Grid update
Now we need to determine the precise grid.In [24], the grid interval d q (q = 1, 2, . . ., Q) is assumed to fulfil a uniform distribution between two adjacent grids.By setting the gradient of Q V ( ) over d q to zero, one can get the update rule for d q .Since the coarse dictionary matrix B is known, then the measurement matrix can be obtained via F = B + E. Alternatively, the accurate grid w can be directly achieved via a polynomial root method.By setting Q V ( ) over v w q = exp −jp sin w q to zero yields to [25] a H w q a w q o q + p q = 0 (18 where o q and p q are given by Also, ( 18) can be formulated as According to (20), the grid w q can be updated by finding the roots of a known polynomial.Therefore, the root method can work with a very coarse grid.It has been suggested in [25] that the updated w (new) where the superscripts 'old' and 'new' denote, respectively, the obtained grid in the last iteration and in this iteration.
Since the polynomial rooting is computationally inefficient, thus it is not recommended to update every w n in each iteration.Instead, a hard threshold is set to select h (h ≥ R) 'proper' active grid points, which is implemented by finding the first h maxima of mq• 2 (q = 1, 2, . . ., Q), where m denotes the estimated mean matrix m.However, this scheme still suffers from the slowness of convergence.In fact, the number of active points would decrease with the increasing iteration number, and it would barely change after several iterations, which coincides a forgotten curve.An illustration is given in Fig. 2, where three independent trails were carried out (SNR = 5 dB, other experiment conditions are the same to that in the simulation section), the forgotten ratio is figured by calculating the ratio of dominate grids (recognised if the associated elements in m is bigger than a threshold) versus iteration number.To accelerate the polynomial rooting, a soft threshold is set to update the grid dynamically.The following forgetting factor f is calculated to determine the active grid points: where h records the grid index, the qth element of which would add one if the 2-norm of mq• 2 is larger than h t , iter represents the iteration number.If f (q) = 1, the corresponding grid point will be updated, otherwise it will invariant in the iteration.To guarantee the ratio of the active grid points, h t should be increased with the increasing iter.In this paper, h t is initialised with 0 and increase with 0.0001 per iteration.
Up to now, we have achieved the proposal of our DOA estimation algorithm for monostatic MIMO radar, which can be summarised as following steps: Step 1: Arrange the matched measurements into R according to (1).Construct the dimensionality reduction matrix according to (3), and compute the non-redundant measurement matrix Y .
Step 3: Update the means m l and the covariance matrix S via (14), and update g q and s 2 via ( 16) and (17), respectively.
Step 4: Compute the forgotten factor f according to (21), and refine the active grid points via (20).
Step 5: Repeat steps 3 to 4 until algorithm convergence, i.e. iteration number reaches max_item or the relative residual is smaller than tol.

Advantages of our algorithm
Compared with the traditional subspace-based methods, e.g.RC-MUSIC method, RC-ESPRIT method and the off-grid SBL methods, the advantages of the proposed method are given as follows: (a) The proposed algorithm does not require the prior of sources number, which means the proposed method has a blind characteristic, while the subspace methods are not blind; (b) The proposed method can work with a small snapshot as well as coherent sources while the traditional subspace methods may invalid in such backgrounds.(c) The proposed method can achieve an accurate DOA estimation while the traditional peak-searching methods can only find the on-grid solutions.(d) The proposed require fewer iteration steps than the existing off-grid SBL methods, and it has better estimation accuracy than the existing off-grid SBL methods.] with an interval 6°.Fig. 3 shows the result of 200 independent trials.Clearly, the proposed method is effective for one snapshot scenario.

In
In the second simulation, we test the DOA estimation performance versus grid interval, and we compare the proposed method ] , respectively, where the resolution is set to 0.01°.Other conditions are the same to that in the first simulation.Fig. 4 depicts the resultant RMSEs at the different grid interval.It can be observed that the RMSE performance of the proposed method gradually improves with the growing grid interval, and it will keep stable after the interval is larger than 4°.Besides, the proposed method provides more accurate DOA estimation performance than OGSBI and ROGSBL.Fig. 5 gives the average iteration number versus grid interval, from which we can observe the proposed method requires nearly the same iteration number than ROGSBL.Note that in each iteration, the rooting time of the proposed method is no more than that in the ROGSBL, thus the proposed is more efficient than ROGSBL.
In the third simulation, we test the DOA estimation versus SNR, where L = 40 and other conditions in the first simulation keep unchanged.Fig. 6 illustrates the comparison, from which we can see that RMSE of the proposed method and ROGSBL gradually improved with the increasing SNR.Moreover, the proposed has lower RMSR performance than ROGSBL and OGSBI, since the dynamic determination of the updated grid number results in a more accurate dictionary matrix.

Conclusion
In this paper, we have developed a reduced-dimensional SBL-based DOA estimation method for monostatic MIMO radar.It can deal with the off-grid problem in DOA estimation, and it has an estimation better accuracy than the existing off-grid SBL methods and subspace-based methods.The proposed method does not require the prior information of source number, and it supports small snapshot scenario, which will lead to a brighter prospect in applications.Finally, numerical simulations are given to verify the improvement of the proposed method.

Acknowledgment
This work was supported by China NSF grants (61701046, 61471191 and 61501233).

Fig. 1
Fig. 1 Bistatic MIMO radar configuration this section, 200 Monte Carlo trials were performed to verify the improvement of the proposed method (Matlab code is available on https://pan.baidu.com/s/1dF6MyTb).The ULA-based monostatic MIMO radar is configured with M = 6 transmit elements and N = 4 receive elements.Assume that there exist K uncorrelated sources located at the angles u k , (k [ 1, 2, . . ., K { } ).Assume the RCS coefficients fulfil the Swirling II model and L snapshots are collected.The signal-to-noise ratio (SNR) is defined as SNR = 10 log 10 R − E 2 2 / E 2 2 dB [ ], where R and E are the matrices in (1).The following root mean square error (RMSE) is used for performance k − u k 2 where ûi,k represents the estimates of u k for the ith Monte Carlo trial.In the first simulation, we assume K = 2 sources located at directions u 1 = −25.4°,u 2 = 7.54°, L = 1 and SNR = 10 dB are considered.The coarse grid in the proposed method is set to −90°, 90°[

Fig. 2
Fig. 2 Forgotten character of active grid points