A Novel High-Resolution Imaging Method Using Reduced-Dimension Beamspace Unitary MUSIC for OFDM-MIMO Radar

Imaging of maneuvering targets has become a challenging topic recently. This study proposes an improved MUSIC-based single-snapshot imaging method for OFDM-MIMO radar. At the transmitter part, for orthogonal waveform generation, the interleaved OFDM structure is developed that each transmitter only radiates at an equidistant set of OFDM subcarriers unique to itself. At the receiver part, a novel preprocessing scheme is presented to separate multiple transmit/receive paths and compensate the reference point. Preprocessed echo matrix can fit the model of two-Dimensional (2D) MUSIC algorithm, except for the coupling between radial/cross-range in the virtual antenna array. Correspondingly, we modify the 2D ElementSpace (ES)-MUSIC. In the meantime, a 2D Reduced-Dimension BeamSpace Unitary-MUSIC (RD-BS-UMUSIC) algorithm is proposed for computational complexity reduction and estimation performance enhancement. Numerical results of a simulated target consists of limited number of scatterers and a complex target (a Boeing 737-800 plane) verify that the proposed algorithm can improve the imaging quality and reduce computational complexity effectively compared with the traditional method.


I. INTRODUCTION
The radar imaging technique provides sufficient information for target feature extraction and recognition. Traditionally, a radar system obtains a two-dimensional (2D) image from multiple wideband returns (multiple real receive antennas or slow-time samples). However, real aperture imaging radar is immobile and expensive. The Inverse Synthetic Aperture Radar (ISAR) technique requires a comparatively long Coherent Processing Interval (CPI). And even worse, for non-cooperative targets, especially those traveling at high-speed or with maneuvering, the motion compensation procedures are considerably complex, and the specific cross-range scaling is hard to attain. Multiple Input Multiple Output (MIMO) radar technique [1]- [4] has the potential to solve the problems of traditional radar imaging methods and a new research hotspot in recent years. A MIMO radar system consists of M transmit (Tx) and N receive (Rx) antennas. By adding diversities of the waveform, frequency, polarization, or spatial, up to MN , independent observation The associate editor coordinating the review of this manuscript and approving it for publication was Maurizio Magarini . channels could be formed, and the degrees of freedom are greatly improved.
Several approaches for MIMO radar system configuration were proposed in [5]- [12] to fulfill imaging with the echo in a single snapshot duration, hence the complex motion compensation could be avoided. Similar to the traditional radar imaging techniques, the range resolution was achieved by transmitting a wideband signal in [5], [7], [8]. The azimuth resolution increased with the widened virtual antenna array, which was the spatial convolution of the monostatic Tx/Rx Uniform Linear Arrays (ULAs). However, the range alignment and compensation were still necessary before the azimuth imaging. In [6], [10], [12], narrowband bistatic MIMO systems were developed. In [10], the 2D cross-range images measured by bistatic Tx/Rx antenna arrays were used to reconstruct a 3D image. The hardware cost of each transmitter was reduced with the narrowband transmission, and the range alignment steps were removed since the radial-range resolution is low. However, the array size and element number should be large enough for an acceptable resolution. To further expand array aperture size while reducing element number, wideband MIMO radars with sparse arrays were investigated [9], [11]. In [9], the Rx antenna array is sparse and non-uniform. In [11], the fast-time sampling was also sparse, and a 2D (range and azimuth) sparse model was established for simultaneous 2D imaging.
In the literatures mentioned above, Tx antennas emit waveforms which are mutually orthogonal, whereas at each Rx antenna, a matched filter set associated with multiple transmitted waveforms was utilized for waveform separation and pulse compression. With the waveform diversity, the virtual array aperture is extended for higher resolution. Thus, the imaging performance relies on the orthogonal transmit waveform set considerably. However, the ideal orthogonality is hard to realize in practice, especially with arbitrary Delay and Doppler frequency offset. In the past few years, Orthogonal Frequency Division Multiplexing (OFDM) technique, inspired by its application in communication, has attracted attention in radar community. Applications of OFDM as radar waveforms were reported in [13]- [20], and the advantages of OFDM radar are listed as followed: extra frequency diversity, no coupling between delay and Doppler, capacity of Fast Fourier Transform (FFT) processing with low hardware cost, flexibility of waveform parameters design, potential of integrated radar and communication systems, and fine performance of anti-interference.
In this paper, we combine OFDM with the MIMO scheme for high-resolution 2D imaging. Investigations of such applications were reported in [21]- [32]. Interleaved OFDM (I-OFDM) signals were firstly introduced in [22] for Simultaneous Polarimetric Measurements (SPM), the twice I-OFDM (I-OFDM) signals were verified theoretically orthogonal to each other and could significantly increase the efficiency of SPM. Inspired by the twice I-OFDM, a new way of orthogonal waveform set generation was proposed in [23]- [26] by further interleaving conventional OFDM signal. In [30], with sufficient Cyclic Prefix (CP), OFDM waveforms emitted by multiple transmitters shared the same frequency band and were mutually orthogonal on each subcarrier in the discrete frequency domain. Complex Orthogonal Designs (CODs) were applied for subcarriers coding.
The conventional imaging method is the Delay-and-Sum (DAS) algorithm implemented via the Fast-Fourier-Transform (FFT). However, the FFT spectrum has comparatively wide mainlobes and high sidelobes, whereas the resolutions were limited by the bandwidth and array aperture. The MUltiple SIgnal Classification (MUSIC) algorithms [33]- [37], most frequently used for angle estimation, were recently applied in Through-the-Wall Imaging (TWI) [38], MIMO radar [23], [25], [26] and ISAR [39] for azimuth imaging due to its super-resolution ability. In these works, the radial-range profile was firstly retrieved computing FFT. In [25], [26], [39], within each range unit, the cross-range was estimated using MUSIC afterwards by exploiting the phase differences among the virtual elements (virtual antennas in MIMO or Doppler virtual elements in ISAR). The transmitted waveform bandwidth needs to be high enough to meet the radial-range resolution requirement, whereas the computation and storage cost is heavy for these methods. In [38], the MUSIC algorithm was applied to the Fourier spatial beams instead of the raw sensor data. The proposed algorithm was referred to as BeamSpace MUSIC (BS-MUSIC) and outperformed the conventional ElementSpace MUSIC (ES-MUSIC) in a low Signal-to-Noise Ratio (SNR) scenarios due to the processing gain offered by beamforming.
In the proposal, we establish a 2D imaging model for a MIMO radar system with monostatic Tx/Rx ULAs and I-OFDM sets. Abandoning the matched filtering set, we separate echoes from multiple Tx antennas and compensate the reference point in the preprocessing procedures. The matrix after preprocessing contains echo data from multiple Tx/Rx paths and on multiple subcarrier frequencies. We adopt the 2D MUSIC algorithm for the simultaneous high-resolution 2D spectrum estimation. Ignore the physical coupling among real antenna arrays, we consider the coupling between radial/cross-range in the virtual element array brought by the I-OFDM configuration. Modified 2D ES-MUSIC is introduced to deal with the coupling. Then in 2D RD-BS-MUSIC algorithm, we employ the conjugate centrosymmetric Discrete Fourier Transformation (DFT) matrix for Fourier beamforming and propose a target beam selection method. The MUSIC estimation works in Fourier beamspace and involves only real-valued computations with reduced matrix dimension. The numerical experiments use synthetic data of a simulated target composed of several ideal point scatterers and a real Boeing 737-800 plane. Simulation results show that the proposed method provides better imaging quality compared to DAS, OMP, 1D ES-MUSIC, and 2D ES-MUSIC, with comparatively low computational complexity. Estimation error analyses concerning SNR and algorithm parameters are also given.
The remaining sections are organized as follows. Section II introduces the 2D imaging model of I-OFDM MIMO radar and preprocessing steps. Next, Section III details the imaging method, including the modified 2D ES-MUSIC and 2D RD-BS-UMUSIC algorithm, followed by the analysis of computational complexity. The numerical experiments and simulation results are presented in Section IV. Finally, this paper is concluded in Section V.
Notation: (·) T , (·) H , (·) − , (·) † denote transpose, conjugatetranspose, inverse and pseudo-inverse operation, respectively. The bold capital and lowercase letters represent the matrices and vectors, respectively, whereas the euclidean vector from A to B is represented by − → AB, and − → AB denotes modulus of the vector. ⊗, * and are the Kronecker, Hadamard and Khatri-Rao product, respectively. diag(a) is a diagonal matrix whose diagonal is the vector a. mod (·) denotes the remainder, and · is taking integers downwardly.

A. SIGNAL MODEL
The signal model of the proposed OFDM-MIMO radar is built in this section. As is depicted in Fig. 1, the radar system VOLUME 8, 2020  Rx and the m-th Tx antenna, respectively. All antennas are omnidirectional. The system configuration, along with an orthogonal transmitted waveform set, leads to an equivalent ULA of Q = MN virtual single-input-single-output (SISO) elements [1].
An original OFDM waveform involves P subcarriers transmitted simultaneously, each of which is modulated by a K -length phase code sequence. To ensure subcarriers orthogonality while maximizing the spectrum efficiency, set ft b = 1, where f is the frequency spacing between two adjacent subcarriers, and t b denotes the OFDM bit duration. In real applications, the tail of the OFDM bit is included at the beginning of each transmitted symbol [30] as a cyclic prefix (CP). Therefore, the final OFDM symbol width t s = t b + t c , with t c = αt b being the CP width.
For inter-symbol-interference (ISI) removal, set t c ≥ τ max , where τ max is the maximum delay of the imaging area. Bandwidth and pulse width of the waveform are Bw = P f and t p = Kt s = K (1 + α) / f , respectively.
The multiplex of Tx antennas via subcarrier interleaving is performed to combine OFDM technique with MIMO radar [23]. As the time-frequency structure of the transmitted waveforms shown in Fig. 2, an equidistant and unique set of subcarriers is assigned to each Tx antenna. The m-th transmitted waveform has L = P/M (w.l.o.g., P is the integral multiple of M ) subcarriers which are equally spaced by M f . The (l, m)-th carrier frequency is f l,m = f 0 + (lM + m) f , with f 0 denoting the system initial carrier frequency. Hence, the k-th symbol of the m-th I-OFDM waveform is specifically given by l,m is the payload subcarrier weighting and phase code modulated onto (l, m)th subcarrier and k-th symbol, satisfying for waveform energy normalization, where E k {·} denote the expected value toward k, trace (·) is the trace of the matrix which is the sum of the diagonal elements, l,m ∈ C L×K is the code matrix in the k-th bit. The design and analysis of the payload codes were reported in our previous work [15] and not detailed in this paper. With the scheme, for one thing, the multiple transmit waveforms have a frequency diversity and hence are mutually orthogonal; for another, they have the same bandwidth without loss of the range resolution. Now consider a far-field target modeled as a collection of I ideal point scattering centers S i . Generally speaking, the OFDM pulse width t p is a short duration, that the target could be assumed still, or the velocity could be estimated and compensated based on our previous work in [17], hence the Doppler effect is not considered in this work. At each Rx antenna, the received signal is the superimposition of echo emitted from M Tx antenna and reflected by I scatterers n (t) is the additive white Gaussian noise (AWGN). (2) is down-converted to baseband, sampled from the starting point T min , and the sampling frequency

Echo described in
i.e., all the relative delays lie within the CP duration. A suitable system parametrization can achieve this. Hence, after removing CP in front of each OFDM symbol, there is no interference from adjacent symbols, and what remains is the cyclic shift of the original OFDM bit. We store the rest sampling points in a column vector r (k) n , the w-th element is given by n is a P-length measurement column vector, p-th element of which is n is the discrete noise vector and ϒ (k) n is its DFT vector element-wise multiplied with α where

B. PREPROCESSING STEPS
According to (4), multiple subcarriers could be separated via DFT towards the term w and be demodulated as following Based on the interleaved structure, an L × P partial identity matrix m could be used to obtain (m, n)-th spatial channel, i.e. the Tx-Rx pair, from R n : l-th row vector of m is the lM + m-th row vector in the P-order unity matrix I P . R m,n is the measurement matrix with noise of the spatial channel between m-th Tx and n-th Rx antenna, involving frequency-domain echoes at L uniformly spaced subcarriers.
The measurement matrix R m,n is then left-multiplied bȳ A m,n for reference point compensation. Set O s as the reference point (as depicted in Fig. 1), which is usually the centre of the imaging area, define an L × L diagonal compensation matrix, the l-th diagonal element is where d  [8], it could be deduced that The secondary and higher order terms in (8) are ignored for far-field assumption. Further, u i =s i u T 0 denotes the radial-range, hence u i u 0 is the projection vector ofs i on u 0 . Correspondingly, define v i =s i − u i u 0 as the discrepancy vector. Since the m-th transmit antenna − − → OT m = mNdx, and the n-th receive antenna The compensation result is Z m,n =Ā m,n R m,n ∈ C L×K , the (l, k)-th term of which is specifically given by where λ = c/f 0 and u = c/ (2Bw) denote the system wavelength and the radial-range resolution, respectively, l,m,n are the modified scattering coefficient and noise term without change in intensity. Far-filed and narrowband approximation is adopted that the cross-range frequency offsets among multiple subcarriers are ignored in (9). Denote (µ i , ν i ) as radial/cross-range frequency, expressed as VOLUME 8, 2020 terms in the brackets are dropped to simplify the expression. Then we concatenate the MN matrices row-wisely as definitions of (µ i , ν i ) are shown in (10). For one thing, the virtual element array consists of Q = MN SISO elements [5], whose steering matrix is A v = A r ⊗ A t . Hence, the Degrees of Freedom (DOFs) in cross-range improve by M compared with the traditional phased-array radars. For another, there are both radial/cross-range terms (u, v) coupled in A t due to the frequency diversity among Tx array for I-OFDM transmission. Consequently, the coupling of (u, v) also exists in A v . Modified estimation methods should be introduced to process the coupling.

III. SUPER-RESOLUTION IMAGING METHODS
The traditional imaging algorithm via 2D-FFT suffers from limited resolution (defined by the Rayleigh criterion [40]) and comparatively high sidelobe level. In [23], the radial-range image was obtained by FFT, whereas in each radial-range unit, the MUSIC algorithm was applied for high azimuth resolution that breaks the Rayleigh criterion. This method may lose some inherent coupling information between two directions. Therefore, we modify the 2D ES-MUSIC algorithm to fit the model in (11)- (12) and estimate the 2D spatial spectrum simultaneously in this section. Additionally, to lessen the high computational burden in the modified 2D ES-MUSIC and improve the estimation accuracy, we propose a 2D BS-UMUSIC (U-MUSIC) algorithm with the reduced dimension (RD) technique.

A. MODIFICATION OF 2D ES-MUSIC ALGORITHM
For the data matrix in (11), the adoption of the MUSIC algorithm has to meet the following basic assumptions: 1) Noise covariance matrix 2) Each column of the steering matrix A is linear independent; 3) Signal covariance matrix We then construct two mutually orthogonal subspaces by computing the Singular Value Decomposition (SVD) of Z. One consists of the I dominant singular vectors and denotes the signal subspace E s , whereas the other of dimension (R − I ) is called the noise subspace E n . Based on the orthogonality between target steering vector and noise subspace, the 2D ES-MUSIC spectrum is given by where (u, v) denote radial/cross-range grids, and a (u, v) ∈ C R×1 is the steering vector.

B. IMAGING VIA 2D RD-BS-UMUSIC
The 2D ES-MUSIC has high computational complexity in complex-valued high-order matrix SVD and 2D searching. In this section, matrix Z is firstly transformed into the real-valued Fourier beamspace. Then, only the target beams are selected out for the MUSIC spectrum calculation to reduce computational complexity and improve estimation accuracy. The proposed method is denoted as the 2D RD-BS-UMUSIC algorithm.
We modify the conjugate centro-symmetrized DFT matrix in [41] for Fourier beamforming. Denote W H s and W H v as the beamforming matrices of subcarrier and virtual element arrays. Based on (12), the l-th row vector of W H s is given by where µ (f) l = 2π [l − (L − 1)/2] /L denotes the l-th radial-range spatial frequency. Similarly, the q-th row vector from W H v is expressed as where Q = MN denotes the virtual element number, ν (f) q = 2π [q − (Q − 1) /2] /Q is the q-th cross-range spatial frequency.
Correspondingly, w H s,l and w H v,q represent the Fourier beams steered at the spatial frequencies µ Left multiplied with the beamforming matrix W H s , the subcarrier array manifold is transformed into the Fourier beamspace as D s , the (l, i )-th element of which is as following Based on the Kronecker product property (A ⊗ B) (C ⊗ D) = (AC) ⊗ (BD) with the right matrix size, the (q, i)-th element of the beamspace virtual antenna array manifold D v is Let the final Fourier beamforming matrix be and the steering matrix in 2D Fourier beamspace can be expressed as The complex steering matrix A is transformed to a real-valued one D.
Build the matrix Z = W H Z , W H Z , the noise subspace E n is then estimated through computing the SVD of Z, The 2D BS-UMUSIC spatial spectrum is formed as d (u, v) ∈ R R×1 is the beamspace steering vector towards the 2D coordinates (u, v), and the (qL + l)-th element is given by where l = 0, . . . , L − 1, q = 0, . . . , Q − 1. To note that, in (17), (18) and (22), (µ, ν) are all calculated following (10). In order to take advantage of the beamforming gain and reduce the matrix dimension, we examine Z to choose a small set of target beams and drop the noise beams, before the noise subspace and MUSIC spectrum estimation.
Basically, calculate the incoherent mean energy of the r-th beam, i.e., the r-th row vector in Z as where · 1 is the 1 norm of the vector. If E r ≥ γ , with γ denoting a specific threshold, then the r-th beam will be picked out as the target beam. Accordingly, the data matrix after beams selection is given bỹ J ∈ RR ×R is a sparse selection matrix. For r = 0, . . . , R − 1, if E r ≥ γ , the r-th row vector in the R-order identity matrix I R will be included in J.
The noise subspace is estimated with the RD matrixZ asẼ n , and the corresponding RD-BS-UMUSIC spatial spectrum is defined similar as (21) with d (u, v) the same as in (22). The degree of matrix dimension reduction in RD-BS-UMUSIC is determined by the threshold γ for target beams selection. Here, we employ an adaptive threshold defined as where η = R−1 r=0 E r /R is the total mean energy, ε is an adjustment parameter determined by system Signal-to-Noise Ratio (SNR).
Peaks off BS (u, v) are pick out as û i ,v i |i = 1, . . . , I . The estimated 2D coordinates are automatically paired in this procedure. The MUSIC spectrum amplitude does not represent the corresponding scattering intensity of each grid [33]. Therefore, the scattering coefficients need to be further estimated by least-square (LS) algorithm after peaks searching, the explicit process of which is ignored in this paper. Fig. 3 shows the proposed imaging procedures for I-OFDM signals with a MIMO radar configuration, and the detailed steps are listed as follows: 1) For each Rx antenna, sample the echo from T min with f s = Bw, divide the sampling points into K segments, remove the CP in front of each segment, and store the remaining part in a matrix R n as in (4); 2) for each Rx antenna, compute the DFT of R n for subcarrier separation, demodulate phase coding matrix, and select the row vectors of each Tx antenna to obtain the measurement matrix Z m,n for each spatial channel as described in (5), (6); 3) compensate the initial phase term of the reference point O s specific to each spatial channel, and concatenate the MN compensated matrices row-wise, attain the data matrix Z via (7)-(11); 4) employ the conjugate centro-symmetrized DFT matrices (15)- (20) for Fourier beamforming, yield the real-valued data matrix in beamspace Z; 5) calculate the threshold defined in (25), examine Z to construct the selection matrix J, and get the subsetZ following (23); 6) compute the SVD ofZ and build the noise space E n with the estimated scatterer number I [42]- [44]; 7) form the 2D RD-BS-UMUSIC spatial spectrum f BS (u, v) based on (24); 8) search the peaks off BS (u, v) and estimate the scattering coefficients using LS algorithm [33].

C. COMPUTATIONAL COMPLEXITY
In this section, we analyze the number of real multiplications involved in echo preprocessing and imaging via 2D RD-BS-UMUSIC.
The computational complexity of the preprocessing procedure is 4 × O (KNP log (P) + 2KNP), with P = LM denoting the total number of subcarriers. In Table 1, we compare the computational complexity of the method in [23], the modified 2D ES-MUSIC, and the proposed 2D RD-BS-UMUSIC. The three methods are executed with the same searching grid numbers in both dimension. For the method in [23], for each virtual antenna,P =LM -points zero-padded FFT is executed for radial-range pulse compression, thenL radial-range bins are selected because of the I-OFDM scheme. In each radial-range bin, One-Dimension (1D) ES-MUSIC estimation is applied to a Q × K matrix, with the cross-range grid number beingQ. The 2D MUSIC algorithms have no pulse compression step, but need a 2D searching process ofL ×Q grids. Adding the Fourier beamforming and the target beams selection, the 2D RD-BS-UMUSIC involves only real-valued computations afterwards, and the matrix dimension decreases toR, hence the computing burden is released significantly.

IV. SIMULATION RESULTS
In this section, simulation results of I-OFDM MIMO radar imaging are presented to illustrate the algorithm performance of the proposed methods. Target location errors are then analyzed in terms of SNR and various algorithm parameters. We have to mention that, the simulation parameters in this section are set for theoretical analysis rather than the real-word measuring.

A. RESULTS WITH A SIMPLE TARGET COMPOSED OF SEVERAL SCATTERERS
The simulated system parameters are listed as followed: the initial carrier frequency f 0 = 10 GHz and the wavelength λ = 0.03 m, N = 8 Rx and M = 4 Tx ULAs are adopted, with d = 1.5 m and Nd = 12 m being the corresponding element spacing. Therefore, the equivalent array consists of Q = MN = 32 virtual antennas.
Each transmitted waveform is composed of L = 32 subcarriers, i.e., P = 128 subcarriers altogether. The frequency spacing between two adjacent subcarriers is f = 0.5 MHz, with the total bandwidth Bw = P f = 64 MHz. The bit duration t b = 1/ f = 2 µs, the CP duration is set t c = 0.25t b = 0.5 µs, hence the total transmitted OFDM symbol width t s = t b + t c = 2.5 µs. K = 64 symbols are transmitted consequently as a pulse, the width of which is T p = Kt s = 0.16 ms.
A simulated target model with I = 16 ideal independent scatterers depicted in Fig. 4, size of each point denotes the scattering intensity. Set the reference center O s = [5000 m, 5000 m, 6000 m], the radial/cross-range coordinates (u i , v i ) (computed following (8)) and the corresponding scattering intensities |σ i | are listed in Table 2. {|σ i |} are chosen randomly.  Denote S n and ϒ n as the noiseless echo and AWGN matrices at the n-th Rx antenna, i.e., R n = S n + ϒ n is the received signal. The Signal-to-Noise Ratio (SNR) is then defined as with · being the Frobenius norm of the matrix. In our simulation, set ∀n, SNR n = SNR = −5 dB, the estimated spatial spectrum, reconstructed scatterers and running time are shown in Fig. 5. The proposed imaging methods is compared with 2D-FFT, 2D-OMP, and the method in [23]. To implement 2D-FFT, we reshape each column vector z k of Z as an L ×Q matrix. The zero-padded lengths are set 512 and 128 for u and v, respectively. Without regard to windowing, rectangular window is chosen in both dimensions. Amplitudes of K 2D-FFT outputs are summed up for incoherent integration. The radial/cross-range resolutions for 2D-FFT are where Bw = P f is the transmitted signal bandwidth, L v = Qd/2 is the virtual array aperture, d 0 is the radial range from coordinate center to imaging area center. Thus u = 2.3438 m and v = λd 0 / (Qd) = 5.7960 m in the simulation, whereas the FFT unit width is u 0 = 0.5859 m and v 0 = 1.4490 m, respectively. The 2D-FFT spectrum in Fig. 5(a) has wide mainlobes and high sidelobes, but owns a good reconstruction accuracy, due to the robustness of DAS method to the noise.
According to (10) and (12), the maximum unambiguous intervals of radial/cross-range are In radial-range, u max is determined by the frequency interval f , as for the terms in the brackets, the first one is for ISI avoidance [30], whereas the second one is due to I-OFDM transmission, which is a topic need to be further solved. In this paper, we set M relatively small to meet the application requirements. In cross-range, v max is decided by λ, element spacing d and d 0 . In our settings, u max = 75 m, v max = 185.47 m. Therefore, in Fig. 5 In Fig. 5(b), we apply 2D-OMP for each OFDM symbol and integrate K estimation results incoherently. Comparing Fig. 5(b) with Fig. 5(a), the mainlobe width and sidelobe level decrease drastically. However, on one hand, the estimation process is time-consuming for 2D-OMP of K symbols. On the other hand, the closely spaced scatterers cannot be resolved in the spectrum and the reconstruction error is high, as is depicted in the red square.
In Fig. 5(c), the radial-range image is obtained via 1024-point zero-padded FFT, hencethe radial-range bing width is also u 1 = 0.2930 m; whereas the cross-range image is achieved calculating the 1D ES-MUSIC spectrum with v 1 . The mainlobe width and sidelobe level in cross-range decreases compared with Fig. 5(c). However, as shown in the red square, the high radial-range sidelobe becomes a fake peak whereas the correct scatterer's coordinates cannot be reconstructed.  in Fig. 5(d) takes more time than in Fig. 5(c) because of the 2D searching, but has higher radial-range resolution with smaller reconstruction deviation (shown in the red square). 2D BS-UMUSIC algorithm in Fig. 5(e) consumes less time compared to Fig. 5(c) and Fig. 5(d), with better estimation accuracy. However, the sidelobe level is still high because of the noise. In Fig. 5(f), set the threshold adjustment parameter ε = 0.2, andR = 51 beams are selected for estimation. Taking advantage of the beamforming gain, the proposed method is robust to low SNR with narrow mainlobes and low sidelobes. Meanwhile the computational cost decreases significantly compared to all the other methods.
Scatterers reconstruction results and comparison under SNR = −5dB Then we analyze the reconstruction accuracy of proposed method with respect to SNR and the algorithm parameters.
In this test, definition of SNR is the same as (26). We conduct 200 independent Monte Carlo trials to compute the Root-Mean-Square Errors (RMSEs), in order to measure the reconstruction precision, defined as where û i,h ,v i,h are the auto-paired coordinates of the peaks in 2D RD-BS-UMUSIC spectrum. The RMSE curves corresponding to SNR are depicted in Fig. 6, with SNR varying from −15 dB to 15 dB. Set the searching steps ( u, v) = (u 1 , v 1 ) as in Fig. 5(f), Fig. 6(a) and Fig. 6(b) compare Clearly as depicted in the four subfigures, in low SNR condition (SNR ≤ −3 dB), noise has a great influence on the estimation accuracy, i.e., RMSE decreases fast with the increase of SNR. When SNR varies from −3 dB to 15 dB, the reconstruction precision is more related to searching step and off-grid error, showing little relevance with noise.
In Fig. 6(a) and Fig. 6(b), RMSE u and RMSE v both show a slight reduction with the decrease of ε in low SNR circumstances (SNR ≤ −3 dB), but are not very sensitive to ε with higher SNRs. Since those "weak" target beams may be covered by noise beams in the high-noisy environment. Consequently, the threshold should be set lower with smaller ε to avoid the omission of some target beams. Whereas with high SNR, the noise coverage effect in Fourier beamspace has little influence on the target beam selection. However, the selected beam numberR may increase with lower threshold (small ε), then the following MUSIC spectrum estimation will consume more time (consulting Table 1). Therefore, higher threshold with larger ε is preferred considering the algorithm performance and complexity in high SNR circumstances.
As is seen in Fig. 6(c) and Fig. 6(d), the final stable value of RMSE u decreases with the searching step u, similar as RMSE v . On the other hand, algorithm complexity increases linearly with the decrease of u and v as shown in Table 1. Hence, trade-off needs to be made between the estimation precision and computational complexity for the selection of searching steps.

B. RESULTS WITH A COMPLEX TARGET
In this section, we employ the electromagnetic synthetic echo data of a complex target ( a Boeing 737-800 plane) to verify the algorithm performance. Consider both the resolution and unambiguous interval, one Tx antenna is employed transmitting a full-subcarrier OFDM waveform, with initial carrier frequency f 0 = 10 GHz, bandwidth is Bw = 256 MHz of P = 128 subcarriers, and the receive antenna number is Q = 32. Since the text results in this section are shown for theoretical proving rather than a realistic application, we just set the observation angle step as 0.009 • instead of an exact cross-range calibration. SNR = 5 dB is computed following (26). As mentioned in Section III, amplitude for each grid in the MUSIC spectrum is incongruence with the corresponding scattering intensity. Hence all the spectrum estimation results shown in this section are normalized. Optical image of Boeing 737-800 plane Fig. 7 shows the comparison between 2D-FFT and the proposed algorithm with different searching steps. As can be seen in Fig. 7(a), 2D-FFT spectrum of P × Q-grids has the poorest resolution of both dimensions. In Fig. 7(b)-7(d), the searching grid numbers are P 1 = 128, P 2 = 256, P 3 = 512 for radial-range, Q 1 = 32, Q 2 = 128, Q 3 = 256 for cross-range, respectively, and set the threshold adjustment parameter ε = 0.2 as a constant. The proposed method  can attain images with better resolution by increasing the searching grid number. With the same grid sizes, the proposed spectrum is apparently better resolved than 2D-FFT comparing Fig. 7(a) and Fig. 7(b), due to the narrower mainlobes and lower sidelobes. Different components of the plane could be better focused and distinguished in Fig. 7(c) and Fig. 7(d) with finer resolutions. Fig. 8 exhibits the images generated by the proposed algorithm with different bandwidth Bw and subcarrier number P. The searching grid numbers are P 3 and Q 3 in the subfigures. Contrasting Fig. 8(a)-8(b) and Fig. 8(c)-8(d), it could be found that with the increase of the total bandwidth, the radial-range resolution improves. However, via analyzing Fig. 8(b) and Fig. 8(c), images generated with the same bandwidth but different subcarrier numbers, we discover that the resolution and image quality is not influenced by the subcarrier number. This conclusion is similar with that in (27).
Similarly, Fig. 9 compares the images generated with different array aperture and element number. As can be clearly seen, the cross-range resolution improves with the increase of array aperture, but has little relation to the element number, which is also analogous with (27).

V. CONCLUSION
In this paper, we have formulated the imaging model of OFDM-MIMO radar and derived super-resolution imaging method based on 2D RD-BS-UMUSIC. Different from traditional pulse compression and range unit selection, a preprocessing technique has been proposed based on FFT, and the result data matrix suits the 2D-MUSIC form except for the coupling of radial/cross-range term in the transmit antenna array due to the frequency diversity in I-OFDM structure. 2D ES-MUSIC is modified for simultaneous 2D spectrum estimation, but the algorithm has high computational complexity and the spectrum has wide mainlobes with high sidelobes. We have defined a centro-symmetrized DFT matrix for real-valued Fourier beamforming and proposed a target beams selection technique for matrix dimension reduction. The simulation results have demonstrated that the proposed method can lead to narrower mainlobes with lower sidelobes compared with conventional DAS, OMP, ES-MUSIC methods, and separate the closely-spaced scatterers effectively. The operation time also decreases significantly compared to OMP and ES-MUSIC. In the future work, we will consider how to deal with the reduction of maximum unambiguous interval in radial-range. He is also the author of more than 11 papers and holds two patents. His research interests include array signal processing, auto target detection, and waveform optimization.
BO XIAO received the B.S. degree in automation from Information Engineering University, in 2014, and the M.S. degree in modern communication technology from the National University of Defense Technology, in 2017, where he is currently pursuing the Ph.D. degree in information and communication engineering. His research concerns are radar communication integration waveform design and signal processing. VOLUME 8, 2020