MIMO Radar Imaging With Multiple Probing Pulses for 2D Off-Grid Targets via Variational Sparse Bayesian Learning

Spatial sparsity of the target space has been successfully exploited to provide accurate range-angle images by the methods based on sparse signal reconstruction (SSR) in Multiple-Input Multiple-Output (MIMO) radar imaging applications. The SSR based method discretizes the continuous target space into finite grid points and generates an observation model utilized in image reconstruction. However, inaccuracies in the observation model may cause various degradations and spurious peaks in the reconstructed images. In the process of the image formation, the off-grid problem frequently occurs that the true locations of targets that do not coincide with the computation grid. In this article, we consider the case that the true location of a target has both range and angle-varying two-dimensional (2D) off-grid errors with a noninformative prior. From a variational Bayesian perspective, an iterative algorithm is developed for joint MIMO radar imaging with orthogonal frequency division multiplexing (OFDM) linear frequency modulated (LFM) waveforms and 2D off-grid error estimation of off-grid targets. The targets during multiple probing pulses are modeled as Swerling II case and a unified generalized inverse Gaussian (GIG) prior is adopted for the target reflection coefficient variance at all snapshots. Furthermore, an approach to reducing the computational workload of the signal recovery process is proposed by using singular value decomposition. Experimental results show that the proposed algorithm is insensitive to noise and has improved accuracy in terms of mean squared estimation error under different computation grid interval.


I. INTRODUCTION
Multiple-Input Multiple-Output (MIMO) radar is an emerging field of radar research which has attracted intensive researches. MIMO radars with waveform diversity offer unique advantages over conventional phased-array counterparts, such as extra degrees of freedom, higher resolution [1], [2] and better parameter identifiability [3]. Besides, MIMO radar which employs transmit Linear frequency modulated (LFM) waveforms can achieve a better range resolution [4]. Different conventional phased-array, MIMO radar adopts orthogonal waveforms transmitted from different antenna elements, which avoids transmitting waveform interference The associate editor coordinating the review of this manuscript and approving it for publication was Zaharias D. Zaharis . from each other and enables each received signal separated by using a bank of matched filters [5]. Recently, an orthogonal frequency division multiplexing (OFDM) LFM waveform has been applied to MIMO radar for high-resolution remote sensing [6]. The OFDM scheme is introduced in [7], which is one of the ways to accomplish simultaneous use of several carriers. With each transmitter of the MIMO radar operating at one of the OFDM carrier frequencies, OFDM-LFM MIMO radar employing OFDM LFM waveforms can generate rangeangle images of targets. This article aims at colocated OFDM-LFM MIMO radar imaging.

A. REVIEW OF EXISTING LITERATURE AND MOTIVATION
The goal of the MIMO radar imaging is to provide an estimate of the radar cross sections (RCS) of targets at precise VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ angular locations and distances relative to the radar. One of the most important problems in MIMO radar imaging is to obtain the two-dimensional (2D) range-angle images. In most radar imaging applications, the number of targets in the region of interest (ROI) is sufficiently less than the number of potential source locations [8]. The sparsity of the target scene has been successfully exploited to provide accurate range-angle images in MIMO radar imaging applications, owing to the development of methods based on sparse signal reconstruction (SSR) or compressive sensing (CS) [9]. The problem of synthetic aperture radar (SAR) image degradation at low signal-to-noise ratio (SNR) has been addressed in [10] where spatial sparsity has been utilized. Sparsity in transform domain has also been exploited for image reconstruction in [11] where smart sensing technique is proposed for online and automatic evaluation of carbon fiber reinforced polymer structure integrity. In these methods, e.g., [8], [10]- [12], the angle or the range domain which is continuous space in practice has to be discretized into a finite set, or equivalently, the angles or the ranges of interest are constrained on a fixed grid. The SSR problem can be formulated, where target locations of interest constitute the support of the sparse signal to be reconstructed. By assuming that all true target locations are exactly on the selected grid, an observation model is established for image reconstruction by the SSR based method. Under the assumption that the reflection coefficients of the targets are constant varying from pulse to pulse, which is also referred to as the single-measurement-vector (SMV) case, 1 optimization is a favourable approach in MIMO radar imaging [13] due to its guaranteed recovery accuracy [9]. Further, q optimization based on maximum a posteriori (MAP) approach is proposed in [8] with the small values of q promoting sparser solutions. However, these methods are based on the conventional CS model which adopts the fixed dictionary. Recently, the models combined with parametric dictionary are proposed in [10], [14]- [17]. They have achieved improvements in SAR imaging [10], SAR autofocusing [14], sparse array imaging [15], moving target imaging [16] and time-frequency analysis [17], by optimizing both the off-grid parameters and the sparse solution. Motived by this idea, in this work, we focus on the off-grid problem of MIMO radar imaging during multiple probing pulses, with the reflection coefficients of the targets varying independently from pulse to pulse where Swerling II model is assumed [18]. It should be noted that the off-grid problem is in multiple-measurement-vector (MMV) case. The MMV case is often encountered in practical applications, such as source localization and direction-of-arrival estimation (DOA) [19]. In the MMV case, the sparse signals at all snapshots share the same support [20]. Such joint sparsity has been exploited to improve the probability of successful recovery [21]. Furthermore, sparse Bayesian learning (SBL) algorithm [22], [23] is another popular method for the SSR problems, which induces less structural error and convergence error [24] and flexibly models the possible structures of signals in order to obtain improved performance [25], [26].
The MIMO radar imaging can be regarded as the estimation of the target amplitude response on a two-dimensional (2D) grid, with one dimension being the angular bin and the other dimension being the range bin. It is worth noting, however, that the SSR algorithms aforementioned involve division of the 2D range-angle space of the potential tagets into a mass of grids, and they are based on the on-grid assumption that all targets should fall right on the predefined computation grid so that reliable estimations can be guaranteed [12]. Unfortunately, such on-grid scenario can rarely be met in reality, where there is a gap between the true location and its nearest computation grid which the estimated range-angle points are constrained on. The employment of a coarse grid set often leads to unsatisfactory estimation accuracy because of the off-grid error. Although accurate solution can be achieved by denser grid search following the coarse grid search, this would bring about much more computation burden. Furthermore, too dense grid set may result in unacceptable degradation of computational efficiency as well as highly coherent matrix that violates the condition for the SSR. Another solution is to establish an off-grid model and the improved estimation performance can be expected. With the basis mismatch of the measurement matrix in the observation model caused by the off-grid error taken into account, an off-grid model with Gaussian prior on off-grid error has been successfully applied to DOA estimation [27]. In fact, the noninformative off-grid error is more reasonable. The uniformly distributed off-grid angluar error has been well addressed in [28]. Furthermore, we focus on the off-grid target imaging with the 2D noninformative range-angle error in MMV case from variational Bayesian perspective.

B. CONTRIBUTIONS OF PRESENT WORK
A variational SBL algorithm and an acceleration scheme are proposed for joint OFDM-LFM MIMO radar imaging and 2D off-grid error estimation of targets. The off-grid errors are considered as model errors that are estimated and compensated during image formation, which enables our technique to alleviate degradations and spurious peaks in the reconstructed image. At Bayesian formulation stage, a hierarchical model is developed with a unified generalized inverse Gaussian (GIG) prior adopted for all snapshots. The GIG is referred to as a family of distributions including Gaussian, Laplacian and many other heavy-tailed distributions, which provides more flexibility to model the prior information for target. This is distinguished from simply assigning Laplacian prior for coefficient variance [29], [30]. At Bayesian learning stage, a variational Bayesian technique is incorporated in hyperparameters inference, where the sparsity regularizing parameter is learned rather than fixed. Since the reconstruction process with measurement matrix based on the 2D imaging model is usually a large computational burden. An acceleration scheme is also incorporated in our algorithm by using sequential singular value decomposition (SVD) to reduce the rows and the columns of measurement matrix rather than only the columns [29]. The signal and the noise subspace induced by SVD can be exploited to improve the image quality at low SNR [10].
The proposed algorithm is named as 2D off-grid sparse Bayesian learning with two SVDs (2DOGSBL-2SVDs). We will compare our method with the widely used on-grid MMV model based approaches, such as 1 optimization based algorithms [12], [31] and q optimization based algorithms [32], [33], and the similar work of off-grid model based method OGSBI [29] which has been extended to 2D parameter estimation [30]. Note that it is difficult for these typical methods to solve the problem to achieve good performance, so they are combined with the proposed acceleration scheme in order to more fairly compete with the proposed method. It is shown that the proposed method is more accurate in terms of mean squared error (MSE) under different computation grid interval as well as more insensitive to Gaussian noise.
The remaining sections are organized as follows. Section 2 presents the OFDM-LFM MIMO radar signal model. Next, off-grid OFDM-LFM MIMO sparse imaging model is studied in Section 3. Section 4 details the proposed 2DOGSBL-2SVDs algorithm with 2SVDs strategy. Our simulation results are shown in Section 5. Finally, the conclusions on this article are provided in Section 6.
Notation: Bold-face letters are reserved for vectors and matrices. (·) * , (·) T and (·) H denote conjugate, transpose and conjugate transpose operations, respectively. · 2 and · F denote the 2 norm and Frobenious norm, respectively. I N represents a N × N identity matrix. A ∈ R M ×N and R ∈ C M ×N denote a real and complex-valued M × N matrix A, respectively. diag(A) denotes a column vector composed of the diagonal elements of a matrix A,and diag (x) is a diagonal matrix with x being its diagonal elements. x i denotes the ith entry of x. (A) i , (A) j and A i,j are the ith column, jth row and (i, j)th entry of a matrix A. (·) and (·) take the real and imaginary parts of a complex variable respectively.

II. OFDM-LFM MIMO RADAR SIGNAL MODEL A. BASIC ASSUMPTION
Consider the scenario of the targets in the far field and a narrowband MIMO radar system composed by N transmitters and M receivers. Due to the attractive ambiguity characteristics, such as range resolution, adjacent band interference and matched filtering sidelobe performance [6], OFDM LFM waveforms are utilized to the MIMO radar. In fact, the OFDM LFM waveform employed here can be regarded as a kind of frequency stepped LFM waveform.
Assuming each transmitter transmit LFM signal with the equal chirp rate γ c and operates at one of the OFDM carrier frequency f n . Frequency intervals between the arbitrary two adjacent transmitters f are assumed to be same. The transmitted LFM waveform for the n-th transmitter can be expressed as where f n = f c + (n − 1) f , T P denotes the pulse width and f c is the carrier frequency of the signal on the first transmit channel. It can be easily shown that when f = κ 0 /T s with κ 0 ∈ N † holds, the requirement of orthogonal multitransmission is satisfied with For simplicity's sake, a uniform linear array is applied to both at the transmitter and receiver. The steering vectors can be written as and where d t and d r denote the inter-element spacing of the transmit and receive arrays, respectively, λ is the wavelength. In this article, the targets in our scenario are modeled according to Swerling II case [18], and are assumed to have no translation relative to the radar so that a target can be modeled as a peak in final range-angle image.. However, the model allows for target spinning relative to the radar, with the RCS fluctuations of targets fixing during a pulse and varying independently pulse to pulse. The received signal at the mth receive element can be expressed as where c is the light velocity, b m (θ ) denotes the nth element of b(θ ), a n (θ ) denotes the nth element of a(θ), R(τ ) denotes the instantaneous distance from the scattering center to the radar, σ (τ ) represents the complex reflection coefficients of the targets which fix during a pulse and vary independently pulse to pulse, t represents the fast time, t s = (τ − 1)T s represents the slow time, T s denotes the pulse repetition interval, τ = 1, . . . , T denotes pulse index and T denotes the number of pulses.

B. MATCH FILTERING
For the τ th pulse, D(m, τ, t) can be filtered out as the mth signal by the matched filter set [6] with multiplying by a phase factor e j2πf n t s . The matched filtering of each receive channel with respect to the (n, m)th channel can be represented as where B = γ c T P , t s = (τ − 1)T s represents the slow time and the received signal D MF (n, m, t, τ ) corresponding to each target has the same time delay 2R/c with the instantaneous distance fixing from pulse to pulse, for the assumption that target has no translation relative to the radar during T pulses. We consider P range bins and K angular bins relative to the radar, so that the grids R p P p=1 and {θ k } K k=1 divide the ROI. Let σ p,k (τ ) represent the complex reflection coefficients for the targets in the ROI. Assuming the pth range cell contains K scattering centers with K angle locations, the received signal in the pth range cell corresponding to t = 2R p /c during the τ th pulse can be derived from (6) as follows: by inserting B = 0. Thus, (7) can be rewritten as With respect to (n, m)th channel, the received signal from the ROI can be formed by a superposition of the received signal from the (p, k)th position.
The received signal can be represented in vector form which leads to the on-grid observation model where e(τ ) denotes the measurement noise vector, σ (τ ) = [σ 1,1 (τ ), · · · , σ K ,1 (τ ), σ 1,2 (τ ), · · · , σ K ,P (τ )] T is a vector representing the sampled and column-stacked version of the range-angle image during the τ th pulse, τ = 1, · · · , T and G ∈ C NMT ×PK is the true observation matrix composed of the element in (9), which can be expressed as which is the ( It is evident that the observation element can map any target location in terms of the range-angle plane to the measurement, and the target is continously distributed at the 2D space. As a result, the goal of the radar imaging is to find σ p,k (τ ) given y(τ ), where the key point is the accurate modeling of the observation matrix. For many radar imaging applications like air surveillance, the target of ROI is spatially sparse.
The on-grid assumption based SSR algorithms like 1 -SVD and M-FOCUSS has successfully produce ange-angle image by exploiting the sparsity in σ (τ ). However, these algorithms obtain the accurate results only if the targets exactly fall on the fixed computation grid. Any grid mismatch leads to ambiguities in the estimation due to the leakage of energy over all the grid cells which depends on the kernel used for the construction of the observation matrix. To overcome this problem, we develop the off-grid model for MIMO radar imaging and show its relationship with the on-grid one.

III. OFF-GRID OFDM-LFM MIMO SPARSE IMAGING MODEL
The on-grid model adopts the discretization of the continous range-angle plane of potential targets into PK computation grids. P and K denote the number of points in the range and the angle dimension, respectively. Let (R,θ ) = {(R 1 ,θ 1 ), . . . , (R P ,θ K )} be the set of the fixed computation grid points in range-angle plane. Without loss of generality, letR andθ be the uniform grid with the grid interval of range cell and angle cell being R and θ , respectively, as is shown in Fig.1. By contrast, off-grid assumption is adopted in this section that the true targets are distributed in a continuous manner and off the fixed computation grid points. Suppose the lth location of the true 1), . . . , (P, K )}. Note that the grid point number should be far greater than the true target number L, which is a kind of sparse description of ROI. However, the off-grid error of the true target leads to the fact that y(τ ) is not truly sparse with respect to the observation matrix. To alleviate this problem, the observation element g R p , θ k n,m,τ can be approximated via first-order partial derivative: denote the two fixed kernels, and the true off-grid location is modeled with the fixed grid (the first term) and the 2D off-grid distance (the last two terms). It is evident that, the true off-grid location in terms of the observation during the τ th pulse is regarded as a two-variable function.
Substitute (13) into (11), the off-grid observation model can be obtained as where where G = [G 0 B C], B and C ∈ C NMT ×PK are the matrices with the fixed kernels, and G (α, β) ∈ C NMT ×PK is parameterized by the off-grid range dependent vector (16) and the off-grid anlge dependent vector which ensures the joint sparsity between the σ (τ ) and the two off-grid distance dependent vectors.
It should be noted that the off-grid observation model is closely related to the on-grid one. By setting (α, β) = (0, 0) and G (0, 0) = G 0 , the on-grid observation model can be obtained. The on-grid observation model only exploits the observation matrix with the fixed kernel, which suffers from the quantization error caused by the 2D discretization of potential space. By contrast, the off-grid model incorporates the 2D parameterized observation matrix with first-order approximation of the true kernel. It is observed that off-grid model has a smaller modeling error and higher accuracy than the on-grid one under the same computation grid. Note that higher order approximation of the kernel can be adopted in the off-grid model to achieve less modeling error, but this would lead to heavier computational burden and is out of the scope of this research.
In order to form the image of MIMO radar, we need to obtain not only the the support of coefficient vector, but also the 2D off-grid error (α, β). Aimed at the joint estimation of them, an iterative algorithm is developed from a Bayesian perspective in the following section.

IV. 2DOGSBL AND 2DOGSBL-2SVDs
Since the complex valued reflection coefficients varies from pulse to pluse, our variational SBL algorithm is derived in the MMV case. For multiple snapshots, the off-grid OFDM-LFM MIMO imaging model is given as follows: (1), . . . , e(T )] ∈ C NMT ×T stands for noise matrix.

A. SPARSE BAYESIAN FORMULATION
Some sparse prior has been frequently utilized for SSR based parameter estimation. The variational SBL method is introduced and developed to obtain the uncertainty information during estimation, where the signal is hierarchically modeled to impose a prior that promotes sparsity as well as learning. Instead of using too dense grid, we adopt a set of coarse grids to reduce the cost of computation, which is more practical.

1) GAMMA PRIOR FOR NOISE MODELING
From pulse to pulse, multiple reflection coefficient vectors σ (τ ), τ = 1, . . . , T and noise vectors e (τ ), τ = 1, . . . , T , are assumed to be independently and identically complex Gaussian distributed, i.e., where γ −1 0 denotes the noise variance with γ 0 being the noise precision, I is an identity matrix with a proper dimension and CN (µ, ) denotes the complex Gaussian probability density function (PDF) of a random variable with mean µ and covariance [34].
Since γ 0 is unknown, and we embed this unknown parameter in the hierarchical framework. For Gamma prior is a conjugate prior of the Gaussian distribution, the perturbation precision is modeled as a Gamma distribution

2) UNIFIED GIG FOR REFLECTION COEFFICIENT MODELING
Each column of is jointly sparse. In this case, the sparse reflection coefficient vectors σ (τ ) share the same support. The spatial distribution of the nonzero coefficients is substantially sparse, so the hierarchical prior is adopted which favors VOLUME 8, 2020 most elements of being zeroes. First, we append a complex Gaussian prior to : where = diag (γ ) and γ = [γ 1 , . . . , γ PK ] T ∈ R PK with γ η being the variance of σ η (τ ), η = 1, 2, . . . , PK . Depending on the choice of the prior distribution for the variance parameters in γ at the second level of hierarchy, various sparsity promoting distributions may arise for σ (τ ). To unify these distributions in a single model and take the nonnegative value of γ into consideration, the variance parameters in γ are assumed to follow a GIG distribution [35] p γ η |ρ, κ, λ = (ρ/κ) λ/2 where K λ ( ) is the modified Bessel function of the second kind, the hyperparameters λ, κ and ρ can be selected so as to formulate the widely used heavy-tailed priors. According to (23), the marginal distribution of σ (τ ) can be obtained by integrating out the latent variables as the generalized hyperbolic distribution: And the resulting marginal distribution, the generalized hyperbolic distribution, covers a large number of distributions as special cases. Heavy-tailed prior distributions have been exploited for sparsity inducing, due to the fact that their heavy tails favor a few large coefficients and severely attenuate the substantial noise. Instead of fixing the sparsity regularizing parameter ρ, we infer it from the data and assume it to follow a Gamma distribution with parameters a and b at the third level, i.e.
which provides more adaptiveness regarding ρ during the learning process. Thus, the proposed modeling constitutes a the three-level hierarchical prior that is equivalent to assigning an almost Jeffrey's prior distribution over γ . The first two levels of this hierarchical prior (22) and (23), and the last level (25) is embedded to calculate ρ. Note that the hyperparameters a, b, c and d at the highest level are set to be trivial values to induce non-informative prior on ρ and γ , respectively, which means their posterior depends only on the data and not the prior knowledge. In this article, a = b = c = d = 10 −6 . Although various other heavy-tailed distributions can be generated by the unified GIG function, the related deep study is not covered in this work.

3) NONINFORMATIVE PRIOR FOR 2D OFF-GRID ERROR MODELING
Since there is usually no additional a priori knowledge of (α, β) with the two independent variables α and β. We use its boundedness as the prior information. A 2D uniform distribution is defined on the 2D off-grid error By combining the stages of the hierarchical model, the posterior distribution is given with α and β being independent of each other. The problem in (27) is often intractable to calculate explicitly, since a multidimensional integral is required. In most cases, it is either impossible or very difficult to compute the integral in closed form. Instead of solving it directly, a deterministic approximation method is utilized through variational Bayeisan approximation.

B. BAYESIAN LEARNING
According to the variational Bayeisan technique [36], the posterior in (28) is calculated with , γ , ρ and γ 0 treated as latent variables, and α and β being the parameters. The approximated posterior is factorizable under the meanfield assumption. Furthermore, by minimizing the KL divergence between the true posterior and the approximated one, the approximated posterior w.r.t. the individual variable can be updated via alternate operation. The piteratively updating rocedures of the individual variables and the parameter as follows.

2) UPDATING FOR COEFFICIENT VARIANCE
The approximated posterior w.r.t. γ is given by Substituting (22) and (23) into (32), the approximated posterior for each γ η obeys a GIG distribution [19] and the hth moment of the GIG distrbution is obtained as where η 2

3) UPDATING FOR ρ
The approximated posterior w.r.t. ρ is given by Substituting (23) and (25) into (34), the approximated posterior for ρ obeys a Gamma distribution and the mean of ρ is given by

5) Updating FOR 2D OFF-GRID PARAMETERS
Since no informative prior for the off-grid difference (α, β) is available, the inference of (α, β) can be obtained by minimizing the negative expected log-likelihood function as where ϕ = α T , β T T , C is a constant term independent of ϕ, where G = [B, C]. Next, the joint sparsity of ϕ with σ (τ ) is exploited so that the dimension of ϕ can be reduced according to the number of true targets L for efficiently computing ϕ.
The sparse prior information is exploited, only L entries of α and β that correspond to locations of the L maximum entries of α and β are calculated. It should be noted that we estimate the 2D off-grid error denoted as α and β simultaneously rather than separately in our frame. Because ϕ is obtained by stacking α and β, 2L entries of ϕ that correspond to the maximum 2L entries of γ should be preserved and others be set to zeroes. As a result, ϕ, v and P can be truncated into dimension of 2L × 1 or 2L × 2L. Their truncated versions are still denoted by ϕ, v and P hereafter just for simplicity. Specifically, if P is invertible andφ = P −1 v ∈ R 2L×1 . Otherwise,φ is updated element wise, i.e., at each step we update one ϕ i by fixing up the other entries of ϕ. For i = 1, . . . , 2L, where α l ,β l ∈ ([− R/2, R/2], [− θ /2, θ /2]) should be satisfied, and the subscript −i of the vector stands for the vector with the ith entry being zero.
Thus, the proposed 2DOGSBL algorithm is implemented as follows. First, the hyperparameters γ , γ 0 , α and β are initialized, and and µ (τ ), τ = 1, . . . , T , are calculated by (30) and (31) respectively. Then γ , ρ, γ 0 , α and β are updated by (33), (35), (37) and (38), respectively. The process is repeated until the convergence criterion that < tol or the maximum number of iterations is reached is satisfied, where tol is a user-defined tolerance and the superscript iter refers to the iteration. According to the property of expectation-maximization algorithm [37], the convergence of the algorithm is guaranteed since the function p (γ , γ 0 , α, β|Y ) increases with iteration.

C. 2DOGSBL-2SVDs
In the practical implement of the 2DOGSBL algorithm, the computation burden is often huge, which is mainly caused by the large number of pulses T . In [29], SVD is used to reduce the columns of the measurement matrix Y . In this section, another SVD is used to reduce its rows. We refer VOLUME 8, 2020 to this scheme of use of SVD twice as 2SVDs. The dimensionality reduction of the measurement matrix Y ∈ C NMT ×T preserves the most significant components while it reduces the computation of the signal reconstruction process and the sensitivity to the measurement noise.
Consider the imaging model in (14) where G (α, β) = G I PK α β T . There is a correlation among the row elements used for constructing G with respect to τ and that Rank (G) ≤ NM < NMT can be satisfied when T > 1.
To some extent, the information in G does not significantly increase with a larger T . It calls for the dimensionality reduction of G. Performing the SVD of the observation matrix G = G 0 B C ∈ C NMT ×3PK without unknown α and β, we have The optimal rank-k r approximation of G, denoted as G r , is computed by According to the theory of principal component analysis (PCA) [38], G r can be obtained by preserving most energy of G in [39]. In our work, we consider that the leakage of the energy of G is so small that it can be negligible. Given λ i > 0, k r can be obtained by approximat- λ i to 1, where k r represents the number of the largest k r singular values of G. Thus, the left singular values can be discarded as noise components. Then, with the projection matrix U H G , the dimensionality reduction version of G can be obtained as Most of the energy of Y can be preserved as a good representation with fewer rows, and we have where Y r ∈ C k r ×T and E L = U H G E ∈ C k r ×T . Another dimensionality reduction is implemented on the number of columns of Y r . A similar subspace-based analysis in [29] is adopted. In the noise-free case of (41), we have where V L and V TL are matrices that consist of the first L and the other T −L columns of V c , respectively. Then we have that Y rc = Y r V L ∈ C k r ×L preserves all signal information. Let c = V L and E rc = E L V L . From (41), we have totally new version of (19) which can be expressed as where Y rc , G S r , c and E rc are all new. The joint sparsity still holds in c and the proposed 2DOGSBL algorithm can be applied. As the resulting algorithm, 2DOGSBL-2SVDs is used to solve (42) to estimate c with α and β.
The complexity for the computation of the proposed method is analyzed. It can be shown that 2DOGSBL-2SVDs has a computational complexity in order of O(k r P 2 K 2 ) per iteration while that for 2DOGSBL is O(max(NMT P 2 K 2 , N 2 M 2 T 2 PK )) per iteration. The additional computational workloads of O(max(k 2 r T , k r T 2 )) and O(max(9NMT P 2 K 2 , 3N 2 M 2 T 2 PK )) are for the SVDs of Y r and G in 2DOGSBL-2SVDs, respectively. Given PK > NM ≥ k r , it can be shown that the whole computational workload of 2DOGSBL-2SVDs is less than that of 2DOGSBL.

D. RANGE-ANGLE IMAGE FORMATION
Since the reflection coefficients of targets are assumed to be complex-valued, the power estimate of in range-angle plane is used to form the range-angle image. The locations of the highest peaks in the range-angle image with off-grid error (α,β) are used to estimate the target locations. Let . . , PK . Let ψ be the vector representation of the gridded map (R,θ) of the range-angle plane. ψ η denotes the power of the ηth location which corresponds to the (p,k) = (ceil(η/K ), mod(η, K ))th element of (R,θ ), where (p,k) = (ceil(η/K ), mod(η, K )) , mod(η, K ) = 0 (ceil(η/K ), K ) , others where ceil(·) and mod (·), which are matlab functions, denote the operations of rounding toward positive infinity and returning modulus after division, respectively. The estimate of the power ψ η can be obtained by computing the expectation Thus, the locations of the nearest grid points are estimated using the locations of the highest peaks of the spectrum. Suppose that the grid point indices of the highest L peaks ofψ = [ψ 1 , . . . , ψ PK ] T are (p l ,k l ), l = 1, . . . , L, and the estimated locations are (R p ,θ k ) = (Rp l ,θk l ) + (αp l ,βk l ).

V. EXPERIMENT RESULTS
The MIMO radar under consideration contains N = 5 transmit antennas and M = 5 receive antennas. The antennas are arranged to form a uniform linear array (ULA), with d t = d r = 0.5λ. Carrier frequency f c is 10GHz. Transmitted pulse width T s is 2µs. Bandwidth B is 15MHz. Frequency interval f between two adjacent channels is 15MHz with κ 0 = 30 which meets the requirement of orthogonal multitransmission. The angle grid interval is ϑ · θ , where θ = 0.1 • denotes a constant and ϑ is the factor. A uniform sampling is adopted. The uniform computation grids in angle and range dimension can be written as {−1.6 • · ϑ, (−1.6 • + θ ) · ϑ, . . . , 1.6 • · ϑ} and {998.6m, 998.6m + R, . . . , 1001.5m} with R = 0.1m. Besides, we set the reference range and frequency to 1000m and f c , respectively, in order to reduce the approximation error in (13). The number of snapshots is set to T = 100 in the MMV case. The targets are modeled as Swerling II case. We assume that the coefficients of the targets are standard complex normally distributed. For the proposed 2DOGSBL-2SVDs, we initialize ρ = 10 −4 , λ = κ = 0, α = β = 0, γ 0 = 1 and γ is initialized as a vector with all entries being 1. In our method, the large energy loss of G is not expected and δ > 0.9999 with λ i > 0.1, i = 1, . . . , k r is approximately required. For the sake of simplicity, k r = MN is set in our experiment.
We take L1 and Lq as the representatives of 1 and the more general q (0 < q ≤ 1) optimization methods based on ongrid MMV model, respectively. The value of 0 < q ≤ 1 gives sparser solution and it can be chosen according to [32]. For small values of q have been found to have a higher likelihood of getting trapped in local minima. It should be noted that the large dimensionality of G 0 in our on-grid model leads to a huge burden during signal reconstruction process. Besides, the employment of subspace-idea based SVD [12] can alleviate the sensitivity to the measurement noise in the MMV case. For fair comparison, G 0 used in on-grid model based methods in our experiment is replaced by G S r0 = U H G G 0 obtained by 2SVDs. With the off-grid error (α, β) being zero, the on-grid version of the model in (42) solved by both of L1 and Lq can be written as Then, L1 and Lq are used to solve (45) estimate , and they are referred to as L1-2SVDs and Lq-2SVDs in our experiment.
In addition, we also compare our method with the offgrid MMV model based method, i.e. OGSBI [29]. It should be noted that the large dimensionality of G (α, β) = G I PK α β T in our off-grid model of (19) leads to a great complexity for OGSBI which only models the angle error and α = diag (0) is assumed in problem (19). Although the measurement matrix can be thinner via the employment of subspace-idea based SVD in [29] which to some extent reduces the computation, it is still difficult to solve the problem in (19). For fair comparison, the dimensionality reduction has been preprocessed via the employment of 2SVDs in Section IV. Then, a similar model with (42) solved by OGSBI can be expressed as where α = diag (0). As the result, OGSBI is used to solve (46) to estimate c with only angle error β and hence it is referred to as OGSBI-2SVDs in our experiment. In subsection 5.1, we present the imaging results on three targets in both on-grid and off-grid cases with different SNRs and demonstrate the improvements in visual image quality, as compared with the on-grid model based methods. In subsection 5.2, we provide a quantitative comparison of our approach with on-grid model based methods.
Figs. 2(a) and (b) show the comparative results for the on-grid and the off-grid targets respectively with the SNR of 20dB. We observe that all the methods perform successful in on-grid case at the relatively high SNR of 20dB in Fig. 2(a). However, in off-grid case, spurious peaks appear in the results of L1-2SVDs and Lq-2SVDs in Fig. 2(b). Lq-2SVDs shows its good property of sparsity promoting but sensitivity to the off-grid error, and the off-grid error even causes its misestimation of one target location parameter. By contrast, 2DOGSBL-2SVDs method provides high resolution as well as spurious peak suppression in off-grid case. It is shown that the superior results of 2DOGSBL-2SVDs method benefits from the accurate estimation and correction of the 2D off-grid error.
Figs. 3(a) and (b) show the comparative results for the ongrid and the off-grid targets respectively with the SNR of 0dB. It is shown from Figs. 2(a) and Figs. 3(a) that almost all the methods benefit from their denoising effects, which are closely related to their sparsity promoting and the common 2SVDs. As is shown in Figs. 3(a) and (b), under the relatively low SNR of 0dB, Lq-2SVDs shows its good property of sparsity promoting in low SNR case but sensitivity to the offgrid error. As is shown in the middle column of Figs. 3(b), Lq-2SVDs misestimates the location and amplitude of one target in off-grid case. On the other hand, both in off-grid case and at relatively low SNR, 2DOGSBL-2SVDs method provides high resolution as well as relatively well spurious peak suppression. In fact, both of the noise and the off-grid error may lead to degraded quality of imaging and spurious peaks, which also degrades the accuracy of parameter FIGURE 2. Left-Images reconstructed by L1-2SVDs. Middle-Images reconstructed by Lq-2SVDs. Right-Images reconstructed by 2DOGSBL-2SVDs. The circle represents the on-grid location of a true target, and the filled square represents a target estimate. The darkness represents the normalized amplitude of the target coefficients, and the observed amplitude range is 25dB. The darker the filled square is, the larger the amplitude is.
estimation. The results demonstrate that the proposed method is more robust to Gaussian measurement noise and 2D offgrid errors than L1-2SVDs and Lq-2SVDs. All the three sparse reconstruction based methods may utilize the SVD to improve the performance at low SNR, but it is difficult for only SVD to handle well with the observation model inaccuracies caused by off-grid errors which may lead to degradations and spurious peaks in the reconstructed images. The main reason for the superior performance of the proposed method is that the 2D off-grid errors are carefully modeled with a noninformative prior and well estimated within our framework. It can be seen from Fig. 3 that 2DOGSBL-2SVDs method can handle with both of the 2D off-grid error and the low SNR.

B. QUANTITATIVE COMPARISON
In this subsection, we compare 2DOGSBL-2SVDs with the on-grid model based methods L1-2SVDs and Lq-2SVDs and the off-grid model based method, i.e., OGSBI-2SVDs, in terms of MSE with respect to the grid interval and SNR. We assume there are L = 2 targets. The true range coordinates R 1 and where the superscript i refers to the ith trial. It should be noted that there exists a lower bound for the MSE of on-grid MMV model based methods like L1-2SVDs and Lq-2SVDs regardless of the SNR since the best locations estimate that they can obtain is the grid point nearest to the true location.
Since the angle and the range of true targets are uniformly distributed, the lower bounds can be calculated as LB_angle = (ϑ · θ ) 2 /12 and LB_range = ( R) 2 /12 which are the expectation in the case of limited trials. However, the lower bounds here are computed independently without the coupling effect between the 2D off-grid errors taken into account. This means that the lower bounds may be lower than the actual ones. Fig. 4 presents our experimental results of comparison with both on-grid and off-grid model based methods. In almost FIGURE 3. Left-Images reconstructed by L1-2SVDs. Middle-Images reconstructed by Lq-2SVDs. Right-Images reconstructed by 2DOGSBL-2SVDs. The circle represents the on-grid location of a true target, and the filled square represents a target estimate. The darkness represents the normalized amplitude of the target coefficients, and the observed amplitude range is 25dB. The darker the filled square is, the larger the amplitude is. all scenarios under consideration, 2DOGSBL-2SVDs has more accurate angle and range estimations than L1-2SVDs, Lq-2SVDs and OGSBI-2SVDs. This demonstrates that the proposed 2D off-grid model based method is superior to the on-grid model based methods; what's more, the off-grid model with 2D off-grid errors taken into account is more accurate than the one with only one-dimensional error modeled like OGSBI-2SVDs. It is seen that the estimation errors of the off-grid model based methods become smaller with decreasing grid interval, while the estimation performance of on-grid model based methods like L1-2SVDs and Lq-2SVDs degrades with the grid interval being smaller than 0.2 • . This demonstrates that although a dense computation grid is necessary for on-grid model based methods like L1-2SVDs and Lq-2SVDs to obtain an accurate estimation of the target location since the gap between the true location and its nearest grid point can be narrowed, this would lead to the sensitivity to noise of the observation matrix used for VOLUME 8, 2020 the SSR. However, this can be largely overcome through the off-grid model which has a good performance on describing the true observation model.
It is shown in Fig. 4 that 2DOGSBL-2SVDs has better precision than OGSBI-2SVDs in different grid intervals and SNRs. The main reason is that 2DOGSBL-2SVDs has modeled the 2D off-grid errors as the range-angle dependent variables of the parameterized observation matrix, rather than OGSBI-2SVDs with range dependent variables being zeros. The off-grid modeling error may result in the degraded performance of the method without range dependent offgrid error incorporated. Besides, the proposed method adopts variational Bayesian learning framework with a unified GIG prior for coefficient variance, which has more degrees of freedom to mimic a kind of heavy-tailed distribution so that better adaptive modeling of the sparsity. The lower bound can be viewed as the expectation of the off-grid error between the true location and the computation grid under the assumption that the off-grid errors are not modeled in the observation process and they are uniformly distributed. The off-grid error can be compensated with high probability using our method, because the proposed off-grid model can be viewed as the first order approximation to the true observation model with the two variables. However, there is still coupling between the two dimensional off-grid errors that is difficult to model and compensate, so the MSEs of the proposed algorithm is still difficult to exceed the lower bound.
All experiments are carried out in Matlab on a PC with 3.6 GHz CPUs. For the three off-grid target case with 500 iterations, the average computation time of L1-2SVDs, Lq-2SVDs, OGSBI-2SVDs and 2DOGSBL-2SVDs are 31.19s, 0.04s, 33.13s, 53.64s, excluding the 2SVDs preprocessing that takes about 13.51s in our case. The computational load of 2DOGSBL-2SVDs is relatively more than the on-grid model and one-dimensional off-grid based methods, but this can be justified through the benefits provided by our imaging framework, as demonstrated in our experiments.

VI. CONCLUSION
In this article, an off-grid imaging framework of OFDM-LFM MIMO radar with multiple probing pulses is presented. The orthogonal transmission property of the OFDM LFM waveform is employed in order to obtain the image in range dimension. The targets are modeled as Swerling II case that the reflection coefficients of the targets is fixed during a pulse while varying independently from pulse to pulse. Different from the mostly proposed on-grid model for SSR, we studied the off-grid sparse imaging model for OFDM-LFM MIMO radar. The true off-grid location is modeled with the firstorder partial derivative of a two-variable function, which is exploited to reduce the modeling error due to discretization of a continuous range-angle plane. From a variational Bayesian perspective, an iterative algorithm is developed for joint MIMO radar imaging and 2D position error estimation and autocorrection of off-grid targets. Besides, a subspacelike method referred to 2SVDs is used to reduce the com-putational complexity of the signal recovery process and the sensitivity to noise. Simulation results demonstrate that the proposed approach outperforms existing on-grid and off-grid model based methods.
In future, further experimental validation, including using dedicated samples can be investigated on not only imaging the static point targets with only frequency domain features, but also more complex and dynamic targets using timefrequency representations and their micro doppler characteristics [16]. His research interests include electromagnetic field and microwave technology, and network communication. VOLUME 8, 2020