A Sparse Bayesian Learning-Based DOA Estimation Method With the Kalman Filter in MIMO Radar

: The direction of arrival (DOA) estimation problem as an essential problem in the radar system is important in radar applications. In this paper, considering a multiple-input and multiple-out (MIMO) radar system, the DOA estimation problem is investigated in the scenario with fast-moving targets. The system model is ﬁrst formulated, and then by exploiting both the target sparsity in the spatial domain and the temporal correlation of the moving targets, a sparse Bayesian learning (SBL)-based DOA estimation method combined with the Kalman ﬁlter (KF) is proposed. Moreover, the performances of traditional sparse-based methods are limited by the off-grid issue, and Taylor-expansion off-grid methods also have high computational complexity and limited performance. The proposed method breaks through the off-grid limit by transforming the problem in the spatial domain to that in the time domain using the movement feature. Simulation results show that the proposed method outperforms the existing methods in the DOA estimation problem for the fast-moving targets.


Introduction
Multiple-input and multiple-out (MIMO) radar systems [1][2][3] transmit orthogonal waveforms, and the corresponding matched filters are used in the receiver to distinguish the orthogonal waveforms. Then, the virtual aperture is provided in the MIMO radar system and improves the radar performance in the target detection, estimation, and tracking [4,5]. The MIMO radar systems usually have two types: the colocated MIMO radar with the distance between antennas being comparable with the wavelength [6,7], and the distributed MIMO radar having larger distance between antennas [8][9][10]. Compared with the distributed MIMO radar, the colocated MIMO radar has better characteristics in terms of system synchronization and waveform diversity [11]. Therefore, in this paper, the colocated MIMO radar is used to estimate the direction of arrival (DOA).
In the DOA estimation problem, many methods have been proposed for decades [12,13]. For only one target, interferometer methods using different phases among antennas achieve excellent performances for DOA estimation. Notably, the interferometer methods have lower computational complexity and can be realized in the field-programmable gate array (FPGA) easily [14,15]. However, the interferometer methods can only estimate one target [16]. In the scenario with more than one target, early methods are based on the discrete Fourier transform (DFT) [17,18], where the DOA estimation problem is similar to the frequency estimation from the samples in the time domain. However, the resolution of the DFT-based methods is limited by Rayleigh criteria, and the targets in the same beam cannot be distinguished. To break through the Rayleigh criteria, the super-resolution methods have been proposed. The most famous methods are the ones based on the subspace theory. For example,the multiple signal classification (MUSIC) method is based on the noise subspace [19,20] using the eigenvalue decomposition [21]; Root-MUSIC [22] avoids the peak searching to find the DOA; the estimation of signal parameters via the rotational invariant techniques (ESPRIT) method uses the rotational invariant feature to estimate the DOA from the signal subspace [23][24][25]. In the subspace-based methods, the subspaces are estimated from the received signals with multiple snapshots. For the MUSIC method, more than 3000 snapshots are used to estimate the covariance matrix, so with the pulse repetition interval (PRI) being 1 ms, more than 3 s are needed to estimate the DOA. For a fast-moving target, the target will move more than 1000 m in 3 s. Therefore, the subspace-based methods are not suitable for the DOA estimation in the scenario with fast-moving targets.
To improve the DOA estimation performance with fewer samples, the compressed sensing (CS)-based methods have been proposed [26,27], where the target sparsity in the spatial domain is exploited and transforms the problem of the DOA estimation to that of sparse reconstruction [28][29][30]. In the sparse reconstruction methods, the orthogonal matching pursuits (OMP) [31], stagewise OMP (StOMP), etc., are widely used. To further improve the performance, a sparse Bayesian learning (SBL)-based method has been proposed in [27,32]. However, the SBL-based methods usually have high computational complexity. In the sparse reconstruction method, a dictionary matrix is formulated by discretizing the spatial domain into grids, which will introduce an off-grid error and limit the improvement of the DOA estimation. Hence, the off-grid methods have been proposed by a one-order Taylor expansion, such as in [3,33]. Moreover, an atomic norm-based method is also given in [34,35], where a semidefinite programming (SDP)-based method is formulated to solve the atomic norm minimization problem.
In this paper, the DOA estimation problem is investigated in the colocated MIMO radar system, and the fast-moving targets are considered. The system model with moving targets is first formulated. Then, To exploit the target sparsity in the spatial domain, an SBL-based method is proposed to estimate the DOA. Since the performance of the SBL-based method is limited by the off-grid problem, we combine the SBL with the Kalman filter and propose an SBL Kalman filter (SBLKF) method, where the temporal correlation of DOA is exploited for the moving targets. Furthermore, simulation results are also given for the proposed method and compared with the existing methods.
The remainder of this paper is organized as follows. The MIMO radar system model for moving targets is formulated in Section 2. A DOA estimation method, combined with the SBL and KF for the fast-moving targets, is proposed in Section 3. Simulation results are given in Section 4. Section 5 concludes the paper.

Signal Model in MIMO Radar
Considering the colocated MIMO radar system, M antennas are adopted to transmit the orthogonal waveforms, and N antennas are for receiving, as shown in Figure 1. In the m-th (m = 0, 1, . . . , M − 1) transmitting antenna, the transmitted waveforms are denoted as s m (t) with the PRI being T p . Since the orthogonal waveforms are transmitted, we have where T is the pulse duration, and * denotes the conjugate. Therefore, in the n-th (n = 0, 1, . . . , N − 1) receiving antenna, we use M matched filters corresponding to the transmitted waveforms to distinguish the orthogonal waveforms, where the filter banks are the same among antennas. For the m-th mathed filter, we have h m (t) = s * m (T − t). With K far-field targets having the same distance between the radar system and the targets, we consider the DOA estimation problem for these targets. With the scattering coefficient of the k-th (k = 0, 1, . . . , K − 1) target during the p-th pulse being α [p] k and the DOA being θ [p] k , the received signal in the n-th antenna can be expressed as where λ denotes the wavelength, d is the distance between the adjacent antennas, τ is the delay from the transmitter to the receiver, and w n (t) is the additive white Gaussian noise (AWGN). In the n-th receiving antenna, the signal r n (t) passes the m-th matched filter, and sampled at t = T + τ as n,m .
Upon collecting all the sampled signals into a vector, we have where (·) T denotes the transpose operation, and the signal y [p] contains the signals from all the receiving antennas (y n , p = 0, 1, . . . , P − 1), and for the n-th antenna, we define the signal y The noise vector w [p] has the same form, so the received signal can be rewritten as where w [p] denotes AWGN with the variance being σ 2 w , and we define the steering vectors of transmitter and receiver respectively as The steering vectors are obtained from the structure of the MIMO radar system [7,36]. As shown in Figure 1, the distance between the adjacent antennas is d, and the DOA for the k-th target during the p-th pulse is denoted as θ k . Therefore, the delay between the adjacent antennas can be obtained via d sin θ [p] k /c, where c is the speed of electromagnetic wave, so the phase change is 2π d λ sin θ [p] k ; and (3) can be obtained from the phase change. Then, collect all the phase change into a vector, and the corresponding steering vectors in (7) and (8) can be obtained.

Target Movement Model
For the movement targets, the DOA can be measured by a movement state. For the k-th target, the state during the p-th pulse is described by a vector as whereθ [p] k denotes the acceleration of DOA θ k . Then, with the AWGN in DOA being n k,1 and that in the acceleration of DOA being n k,2 , the target movement is shown as which can be rewritten as the following vector form where we define F 1, T p 0, 1 and n k n k,1 n k,2 .
Note that the movement model is for the DOA instead of the range. Since the movement of a target in range can be easily transformed into the movement of DOA, we directly describe the target movement as a DOA model.

Sparse Bayesian Learning-Based DOA Estimation
To estimate the DOA during the p-th pulse, the target sparsity in the spatial domain can be exploited to improve the estimation performance. First, the sparse system model will be formulated. A dictionary containing all possible steering vectors for target DOAs can be expressed as where Q denotes the number of columns in the dictionary matrix D. The spatial is descretized into Q angles, and formulates a vector Hence, the q-th column in the dictionary matrix D is where φ q is the q-th discrete DOA. Then, the system model in (6) can be rewritten as the following sparse model where the vector x [p] is a sparse vector denoting the scattering coefficients, and the positions of non-zero entries indicate the corresponding DOAs. For example, if the q-th discrete angle φ q is equal to the DOA k . To reconstruct the sparse vector x [p] in (16), a SBL-based method is proposed in this paper, and we first show the probability density function (PDF) of AWGN (w [p] ) as where the noise variance σ 2 w is usually unknown, and we define the precision of noise variance as β = σ −2 w . Then, β follows a Gamma distribution as where β 1 and β 2 are the hyperparameters, and the Gamma function is defined as To simplify the analysis, as the Gamma distribution is a conjugate prior of the Gaussian distribution, we assume that the target scattering coefficients x [p] follow the following Gaussian distribution: where the covariance matrix Σ is a diagonal matrix Σ diag{σ 2 x,0 , σ 2 x,1 , . . . , σ 2 x,Q−1 }. Similarity, the variances of the entries in the sparse vector x can also described by a Gamma distribution where we define γ q = σ −2 x,q . g and h are the hyperparameters of γ γ 0 , γ 1 , . . . , γ Q−1 T . Therefore, the DOA estimation problem with the unknown noise variance can be described by a problem to maximize the posterior probability, which is expressed as {x [p] ,β,γ} = arg max f (x [p] , β, γ|y [p] ).
Then, for the sparse vector x [p] , the estimated mean and variance can be obtained by the following theorem Theorem 1. With the received signal in (16) and the Gaussian distribution assumption in (20), the mean and variance of x [p] are, respectively, Proof. To estimate the mean and variance of the sparse vector, we first formulate the joint distribution Then, with the joint distribution f (y [p] , x [p] , β, γ), the posterior for x [p] is obtained as , γ) can be obtained as and Substitute (27) and (28) into (26), and we can find that the sparse vector x [p] follows the Gaussian distribution. The mean and variance can be obtained as From Theorem 1, we find how to estimate the mean and variance of the sparse vector x, but the unknown parameters including β and γ are still unknown. The expectation-maximum (EM)-based method is used to estimate these unknown parameters, and can be obtained by the following theorem Theorem 2. The unknown parameters β and γ in Theorem 1 can be obtained by the EM-based method aŝ Proof. The likelihood function for the unknown parameters β and γ can be expressed as Therefore, for the noise precision β, we can simplify the likelihood function by ignoring terms independent thereof in (33) as Then, the noise precision can be estimated by maximizing the likelihood function L(β). Set ∂L(β) β = 0, and the estimated noise precision iŝ Similarity, for the precision γ, we have the following likelihood function By setting ∂L(γ) γ = 0, we can estimate the q-th entry in γ aŝ The details for the sparse reconstruction using the SBL-based method are shown in Algorithm 1. Using the proposed method, the DOA can be estimated. Moreover, the movements of targets can be exploited to improve the estimation performance further.

The KF-Combined SBL Estimation Method
For the far-field targets, the movement of targets are stationary, and cannot change significantly. Thus, the KF can be used to improve further the DOA estimation performance when combined with the SBL-based estimation method, the movement model in (12), and the estimated DOAθ k obtained by Algorithm 1.
During the p-th pulse, the estimated target state is denoted aŝ where the estimated stateŝ k−1 for the k-th target has the DOA and acceleration information. Then, the transition matrix for all the targets can be expressed as

Algorithm 1
The SBL-based DOA estimation method.
1: Input: the received signal y [p] during the p-th pulse, the dictionary matrix D, the hyperparameters β 1 , β 2 , g, and h. 2: Initialization: set the values of hyperparameters as 10 −4 , and the number of iterations is N iter = 100. 3: for iter = 1 to N iter do 4: Estimate the mean and variance as Estimate the precision of noise variance aŝ Estimate the q-th precision of sparse vector aŝ 7: end for 8: Estimate the DOA by peak-searching the spatial spectrum using the mean µ [p] and variance Σ x . 9: Output: the estimated DOAθ k .
The predicted target state using the KF method can be obtained as and the estimated state during the (p + 1)-th pulse iŝ Then, with the KF theory, the details for the sparse DOA estimation method are shown in Algorithm 2.
In Figure 2, the flow chart of the proposed method is shown, and two main steps, including the SBL-based sparse reconstruction and the KF-based prediction, are combined in the proposed method. Hence, better DOA estimation performance can be achieved, especially for the fast-moving targets.

Algorithm 2
The KF-based sparse DOA estimation method. Additionally, the proposed method can also deal with the off-grid problem caused by the discretized grids in the sparse reconstruction method. In the SBL-based reconstruction step, the on-grid DOAs are estimated. To further improve the estimation performance, the KF-based method is combined to predict the DOAs, so the off-grid errors can be reduced. For example, when the ground-truth DOA is 30.2 • , the traditional SBL method cannot breakthrough the resolution of grid size. The estimated DOA using the SBL method can be 30 • or 31 • with the grid size being 1 • . However, when the proposed method is used, the resolution can be less than 1 • when the KF is combined with the SBL method. Therefore, the proposed method can solve the off-grid problem.

Computational Complexity
The computational complexity can be estimated from the two main steps of the proposed method. In the SBL-based reconstruction step, we can find that the computational complexities Since the number of antennas is much lower than that of discretized grids in the practically MIMO radar system, i.e., MN Q, and the number of targets is also much lower than that of Q, i.e., K Q, we can find that the computational complexity of the proposed method is O(Q 3 ).

Simulation Results
In this section, the simulation results are given to show the performance of the proposed method in the DOA estimation. The proposed method was carried out in a PC with Intel i7 processor and 16 GB RAM. The proposed method is also available online (https://drive.google.com/drive/folders/ 1eHSYMnwEGy2kf_2NXZSmdbuK4lyIrSEt?usp=sharing).The simulation parameters are shown in Table 1. The targets have movements with constant speed, and the DOA is described as where we can choose the speed ∆θ as 10 −3 . During the 100-th pulse, the ground-truth DOAs are 0 • and 50 • , as shown in Figure 3, and we use 3 methods to estimate the DOA: • SBL method: The method is proposed in [3,27,32,33], where the SBL theory is exploit to reconstruct the sparse vector. The better estimation performance can be achieved, but with higher computational complexity. • OMP method: The method is widely used in the sparse reconstruction theory, and has lower computational complexity. However, this method cannot have better estimation performance in the scenario with a coherent dictionary matrix.

•
Proposed method: This is the method proposed in this paper, where both the SBL and KF are combined. Hence, the off-grid problem can be solved, and better DOA estimation performance can be achieved. As shown in Figure 3, all three methods can find the two targets. The root mean square error (RMSE) in degrees is used to measure the DOA estimation performance, and is defined as We found that the RMSEs of the SBL method, the OMP method, and the proposed method were 0.82558, 0.85187, and 0.29818. Hence, the best estimation performance was achieved by the proposed method. Additionally, we found that compared with the grid size 1 • , the RMSE in degrees (0.29818) of the proposed method is much less than the gird size, so the off-grid problem is solved.
For the DOA estimation with a moving target, Figures 4 and 5 show the estimation results for the first and second targets, respectively. We use three methods to track the target movement, and the proposed method can more precisely track the moving DOA, compared with SBL and OMP methods. The corresponding RMSEs of DOA estimation with a moving target are shown in Figure 6, where the proposed method has relatively lower RMSE with the increasing index of movement. However, for the SBL and OMP methods, the performance cannot be improved using the historical information of the estimated DOA. Therefore, with more estimation information, the DOA estimation performance can be further improved using the proposed method.
In the proposed method, Algorithm 1 uses the iterations to estimate the DOA during the p-th pulse, so we show the DOA estimations with different numbers of iterations in Figure 7. We can find that when the number of iteration is greater than 20, the same DOA estimation performance is achieved. Hence, in the following simulations, the number of iterations is chosen as N iter = 20. Moreover, when SNR ≥ 15 dB, we have acceptable performance (RMSE ≤ 0.2 • ) of the DOA estimation, so we choose SNR = 15 in the following simulations.  For the sparse reconstruction method, the dictionary matrix D is formulated, so the grid size is essential to controlling the performance of sparse reconstruction. In Figure 8, the DOA estimation performance with different grid sizes is shown; we choose the grid sizes as 0.2 • , 0.5 • , 1 • , 1.5 • , 2 • , 2.5 • , and 3 • . The SBL, OMP, and proposed methods are compared in the scenario with different grid sizes. We can find that when the grid size is less than 0.5 • , the coherence between the adjacent columns in the dictionary matrix is higher, so the sparse reconstruction performance is worse. When the grid size is larger than 2 • , the DOA estimation performance is also worse with the additional off-grid problem. As a trade-off between the dictionary size and the reconstruction performance, we choose the grid size as 1 • in the following simulations. Under this condition, the RMSE of the DOA estimation using the proposed method is about 0.2 and much less than the grid size 1 • .
Then, the DOA estimation performances using the SBL, OMP, and proposed methods are shown in Figure 9 with different SNRs. As shown in Figure 9, the proposed method achieves the best estimation performance from 0 dB to 30 dB. The SBL method has better performance than the OMP method but also has higher computational complexity. Additionally, in the scenario with worse SNR (≤ 10 dB), the performance improvement is more significant.  The DOA estimation performance with different numbers of the transmitting and receiving antennas is also shown in Figure 10 and Figure 11, respectively. We can find that the estimation performance is improved with more antennas. However, with the limit of the grid size, the estimation performance cannot be improved significantly when the number of transmitting antennas is higher than 8 and that of receiving antennas greater than 7. Compared with the SBL and OMP methods, the proposed method can achieve better DOA estimation performance. When the grid size is large (i.e., greater than 2 • ), we can find that the RMSE of the DOA estimation using the proposed method is 0.25 • , which is much lower than the SBL and OMP methods, as shown in Figure 8. Additionally, the computational time with the 2 • grid size is also half of the one with the grid size being 1 • .
Finally, the computational complexity is investigated. As shown in Table 2, we show the computation time of the SBL, OMP, and proposed methods with different grid sizes. The proposed method has higher computational complexity than SBL and OMP methods but is comparable with the SBL method. For the computational time, the grid size controls the total computational time. Additionally, the relative time is also added, where we normalized the relative time of the proposed method as 100%. Hence, the computational complexity of the proposed method is acceptable.

Conclusions
The DOA estimation problem for the colocated MIMO radar has been investigated in this paper. By formulating the system model, the SBL-based DOA estimation method with KF has been proposed for the fast-moving targets. By exploiting both the target sparsity in the spatial domain and the correlation in the time domain, better performance has been achieved by the proposed method compared with the existing methods. Furthermore, the proposed method can also deal with the off-grid problem caused by the discretized dictionary matrix. Simulation results have shown both the performance improvement and the effectiveness of the proposed method. Further work will focus on the DOA estimation for fast-moving targets with imperfect MIMO radar.

Conflicts of Interest:
The authors declare no conflict of interest.