# Three-Dimensional High-Resolution MIMO Radar Imaging via OFDM Modulation and Unitary ESPRIT

Authors: Jingjing Zhao; Yongxiang Liu; Kai Huo; Jiaxi Ye; Bo Xiao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Imaging and recognition of targets with complex maneuvers bring a new challenge to conventional radar applications. In this paper, the three-dimensional (3D) high-resolution image is attained in real-time by a Multiple-Input-Multiple-Output (MIMO) radar system with single Orthogonal-Frequency-Division-Multiplexing (OFDM) pulse. First, to build the orthogonal transmit waveform set for MIMO transmission, we utilize complex orthogonal designs (CODs) for OFDM subcarrier modulation. Based on the OFDM modulation, a preprocessing method is developed for transmit waveform separation without conventional matched filtering. The result array manifold is the Kronecker product of the steering vectors of subcarrier/transmit antenna/receive antenna uniform linear arrays (ULAs). Then, the high-resolution image of target is attained by the Multidimensional Unitary Estimation of Signal Parameters via Rotational Invariant Techniques (MD-UESPRIT) algorithm. The proposed imaging procedures include the multidimensional spatial smoothing, the unitary transform via backward-forward averaging, and the joint eigenvalue decomposition (JEVD) algorithm for automatically paired coordinates estimation. Simulation tests compare the reconstruction results with the conventional methods and analyze the estimation precision relative to signal-to-noise ratio (SNR), system parameters, and errors.


---

## Body

## 1. Introduction

In the past few years, the OFDM waveform has been demonstrated as a potential radar waveform [1–8]. An OFDM waveform consists of multiple subcarriers which are mathematically orthogonal in the time domain. Each subcarrier could be modulated by a phase code sequence; hence, an OFDM pulse [5] has high spectral efficiency as linear frequency modulation (LFM) pulse, along with a bit-to-bit diversity as the phase-coded pulse. Besides, the OFDM signal also offers frequency diversity to the radar system. These properties make OFDM signals particularly suitable for full orthogonal MIMO radar applications. According to the allocation schemes of multiple subcarriers, there were different applications. A step frequency technique was adopted [9, 10] that the orthogonality is fulfilled, but the range resolution is limited. Interleaved OFDM (I-OFDM) was proposed and taken in [11–16]. An equidistant and exclusive subset of subcarriers is assigned to each transmit antenna. In this way, the range resolution is preserved; however, the maximum unambiguous range interval is reduced. A novel non-equidistant subcarrier interleaving approach [17, 18] was adopted to attain a full unambiguously measurable range. However, lacking -some subcarriers may lead to an increase inthe sidelobe level. In [16, 19, 20], the complete set of subcarriers were allocated to each transmit antenna, whereas, on each subcarrier, complex orthogonal designs (CODs) [21–23] were adopted for multiple antennas to ensure orthogonality.The three-dimensional (3D) radar imaging provides more detailed and accurate descriptions of the target and is meaningful for feature extraction and identification. Along with the Inverse Synthetic Aperture Radar (ISAR) technique, Interferometric ISAR (In-ISAR) [24–26] and MIMO-ISAR [27–29] techniques have to cope with the complex motion compensation, especially for maneuvering targets. Instead of the time sampling in ISAR, MIMO radar systems utilize spatial sampling. Moreover, with an orthogonal transmit set, the equivalent virtual array could be formed in a single snapshot illumination with fewer elements than a conventional array radar. MIMO radar 3D imaging techniques [12, 14, 30–38] have been proposed in recent years. In [12, 14, 30, 31, 38], with wideband transmission for range resolution, the two-dimensional (2D) cross ranges were resolved with two mutually perpendicular uniform linear arrays (ULAs). Additional range unit alignment and compensation steps were still necessary before cross-range imaging. In [32, 33], assuming all scatterers fall into the same range bin for narrowband radar configuration, the 3D images were reconstructed with two bistatic antenna arrays. The range alignment was then avoided, but more antennas were needed compared to the wideband case. Sparse arrays [34, 36, 37] were utilized to solve the problem, and sparsity recovery-based methods were proposed for high-resolution imaging.The OFDM-MIMO radar configuration is applied in this paper to reduce antennas used for 3D imaging and remove the range alignment steps. We follow the design criterion in [20] that CODs are adopted for subcarrier coding. The uniformly spaced subcarriers form a ULA in the frequency domain for radial-range estimation. The transmit/receive antenna arrays are two ULAs which are located mutually perpendicular and form an L-shape to obtain 2D cross-range resolution. At each receiver, other than the matched filtering set for waveform separation and pulse compression, we apply Fast Fourier Transform (FFT) to the echo for subcarriers separation, and the frequency domain echoes from different transmit antennas could be separated via decoding. The demodulated results of the subcarrier/transmitter/receiver arrays are stored for 3D coordinate estimation. Therefore, no range alignment is needed.To break the Rayleigh criterion of minimum resolvable separation, we present a 3D imaging method based on the Estimation of Signal Parameters via Rotational Invariant Techniques (ESPRIT) algorithm. The ESPRIT algorithm is a commonly used high-resolution parameter estimation method [39–41]. It exploits the rotational invariance among the signal subspace and has no searching process as in the MUltiple SIgnal Classification (MUSIC) techniques. The high-resolution imaging method based on the ESPRIT algorithm was reported in [33, 36, 38, 42–46]. However, three problems have to be taken into account for the ESPRIT algorithm applied in OFDM-MIMO radar imaging. First, the preprocessed data are only one snapshot available that the frequency domain echoed from different scatterers is coherent, and the covariance matrix is no longer full rank. This leads to the failure of the ESPRIT algorithm. Second, for the 3D-ESPRIT estimation, the matrix size is usually large that the complex-valued eigenvalue decomposition (EVD) or singular value decomposition (SVD) is computationally expensive. Moreover, the standard ESPRIT algorithm does not take the utmost of the complex conjugate of the initial data. Finally, the extra pairing procedure is required for 3D reconstruction. The first aforementioned problem could be solved by applying a multidimensional spatial smoothing technique to restore the covariance matrix to full rank. For the second one, involving the backward-forward averaging method and unitary transformation, we apply a 3D Unitary-ESPRIT (3D-UESPRIT) algorithm. Although the snapshot number is doubled after backward-forward averaging, the computational burden is eased since the algorithm involves only real-valued computations afterwards. Last but not least, the joint eigenvalue decomposition (JEVD) [47–51] method is adopted for 3D estimation to diagonalize the three matrix pencils simultaneously, which leads to the 3D coordinates estimation results automatically paired.In this paper, we first establish an OFDM-MIMO radar imaging model based on OFDM's multicarrier structure; the echo is then preprocessed to construct the data vector that suits the ESPRIT form. Multidimensional spatial smoothing and backward-forward averaging techniques are introduced to build the multisnapshots data matrix. Finally, the super-resolution images are obtained by exploiting the shift-invariance property with only real-valued computations. Imaging results of the 2D EM computed data for an Airbus A320 plane and the 3D simulated data for the point scatterer model of a Boeing 777 plane are depicted and compared to prove the efficiency of our proposed imaging methods for OFDM-MIMO radar. The proposed algorithm provides better reconstruction precision compared with the conventional Fourier-based method and has lower computational complexity than the MD-MUSIC and MD-ESPRIT algorithms [42]. Reconstruction accuracy is then analyzed corresponding to the system parameters and system errors for further optimization.Notations:⋅T, ⋅H, ⋅∗, ⋅−, and ⋅† denote transpose, conjugate transpose, conjugate, inverse operation, and pseudoinverse operation, respectively. The bolded capital letters and lowercase letters represent the matrices and vectors, a denotes the modulus of the vector, and A is the norm of the matrix. ⊗, ⊙, and ° are the Kronecker product, Hadamard product, and Khatri-Rao product, respectively. Diaga is a diagonal matrix whose diagonal is the vector a, whereas diagA is the column vector of the matrix A. ℜ⋅ and ℑ⋅ denote the real and the imaginary part, respectively.


## 2. Signal Model and Preprocessing Method

In this section, we briefly introduce an OFDM-MIMO system configuration for 3D imaging, then develop the received signal model, and propose a preprocessing method for transmit waveform separation and imaging center compensation.


### 2.1. OFDM-MIMO System and Echo Model

In the MIMO radar system,P transmit and Q receive antennas form two ULAs. As shown in Figure 1, the transmit array is parallel to the x-axis and uniformly spaced by dt; that is, coordinate of the p-th transmit antenna is tp=t0+pdt,0,0T=t0+pdtx, where p=0,…,P−1 indicates the transmit antenna index, t0=xt,0,yt,0,zt,0T denotes the start of transmit array which is adjacent to the coordinate system center or, and x=1,0,0T is the unit vector of x-axis. Similarly, the receive array is parallel to y-axis with the space being dr, q-th of which is denoted by rq=r0+qdry, where r0=xr,0,yr,0,zr,0T and y=0,1,0T.Figure 1
MIMO radar system configuration.Each transmit element in the radar system emits an OFDM signal modulated at the initial carrier frequencyf0. The baseband OFDM signal consists of N mutually orthogonal subcarriers with a uniform frequency spacing Δf; hence, the total bandwidth is Bw=NΔf. The n-th subcarrier is modulated with a K-bit coding sequence cn,p=an,p0,…,an,pK−1T. Complex orthogonal designs (CODs) [21–23] which fulfill(1)cn,pHcn,p′=1,p=p′,0,p≠p′,are applied [20] to build the orthogonal waveform set. Construct the matrix Cn=cn,0,…,cn,P−1∈ℂK×P and rewrite (1) as(2)CnHCn=IP.The bit widthtb meets the orthogonal condition that tbΔf=1. As a cyclic prefix (CP), tail of the original OFDM bit is added at the head of each transmitted OFDM symbol to remove the Intersymbol Interference (ISI) [7, 20]. Thus, the symbol duration is ts=tb+tc, with the CP duration being tc=αtb. The p-th baseband OFDM pulse is given by(3)upt=∑n=0N−1∑k=0K−1an,pkexpj2πnΔft−k+1tcrectt−ktsts−12,where p=0,…,P−1 and rectx=1,−1/2≤x<1/2,0,otherwise is the rectangular window. The baseband signal upt is then modulated at the initial carrier frequency f0 for transmission.Consider a far-filed target composed ofI ideal point-like scatterers, the echo model is given. The i-th scatterer coordinate is si=xi,yi,ziT, and the relative scattering amplitude is σi, which is a constant within an OFDM pulse duration and for multiple transmit/receive paths. Suppose the target velocity has been fully estimated and compensated before preprocessing [7], the Doppler effect is ignored in the echo model. At the q-th q=0,…,Q−1 receive antenna, let t=τmin+kts+tc+t0, with τmin being the start of the sampling window and t0∈0,tb. The relative round-way delay of the p,q-th transmit/receive path and i-th scatterer is(4)τp,qi=si−tp+rq−sic−τmin,where c is the velocity of light; then the echo after downconversion is expressed as(5)xqkt0=∑p=0P−1∑i=1Iσiexp−j2πf0τp,qi∑n=0N−1expj2πnΔft0−τp,qi⋅∑k′=0K−1an,pk′rectk−k′ts+tc+t0−τp,qits−12+wqkt0,where wqkt0 is the noise term, which is the zero-mean complex Additive White Gaussian Noise (AWGN). Assume τp,qi∈0,tc,∀p,q,i; then, the window function goes as(6)rectk−k′ts+tc+t0−τp,qits−12=δk−k′,where δ⋅ is the impulse function. Hence, there is no intersymbol interference (ISI) in the echo. Furthermore, set the sampling frequency as fs=Bw=NΔf; thus, t0=ltb/N,l=0,…,N−1; the discrete-time echo is written as(7)xl,qk=∑n=0N−1expj2πnlN∑p=0P−1an,pk∑i=1Iexp−j2πf0+nΔfτp,qiσi+wl,qk,where xl,qk and wl,qk denote the sampled echo and noise term, respectively.


### 2.2. Preprocessing Method

At theq-th receive element, for the k-th OFDM bit, compute the DFT of (7) towards the fast-time sampling index l; the frequency-domain echo is denoted by(8)Xn,qk=1N∑l=0N−1exp−j2πnlNxl,qk=∑p=0P−1an,pk∑i=1Iexp−j2πf0+nΔfτp,qiσi+Wn,qk,where Xn,qk is the n-th point of the DFT results and also denotes frequency-domain echo on the n-th subcarrier and Wn,qk is the frequency-domain noise. N subcarriers are separated without any Intercarrier Interference (ICI) in (8), but echoes from P transmit antennas are still mixed. For transmit separation and subcarrier decoding, for each subcarrier index n, we stack the DFT results of K bits in a column vector xn,q=Xn,q0,…,Xn,qK−1T and give the matrix form expression of xn,q as(9)xn,q=CnBn,qσ+wn,q∈ℂK×1,with Cn∈ℂK×P defined in (2) and(10)Bn,q=βn,0,q1⋯βn,0,qI⋮⋱⋮βn,P−1,q1⋯βn,P−1,qI∈ℂP×I,where the p,i-th element(11)βn,p,qi=exp−j2πf0+nΔfτp,qi,is the channel response of the i-th scatterer, on n-th subcarrier, for the p,q-th transmit/receive path, and σ=σ1,…,σIT∈ℂI×1 and wn,q=Wn,q0,…,Wn,qK−1T∈ℂK×1 are the scattering coefficient vector and the frequency-domain noise vector, respectively.Based on (1), left multiplying (9) with CnH yields(12)yn,q=CnHxn,q=Bn,qσ+wn,q′=yn,0,q,…,yn,P−1,qT∈ℂP×1,where yn,p,q denoted the demodulated echo and wn,q′ is the modified noise vector without change of power. Clearly, the P transmit waveforms are separated in yn,q.As is shown in Figure1, let o=x0,y0,z0T be the center of the imaging area, let d0=o be the radial distance, and let d0=o/d0 be the unit vector of the radar Line of Sight (LOS). For each OFDM subcarrier and each transmit/receive pair, we define the following compensation term:(13)βn,p,q0=expj2πf0+nΔfτp,q0,where τp,q0=o−tp+rq−o/c−τmin. yn,p,q is multiplied with βn,p,q0 for initial phase compensation, and(14)sn,p,q=βn,p,q0yn,p,q=∑i=1Iexpj2πf0+nΔfτp,q0−τp,qiσi+w˜n,p,q,where w˜n,p,q denotes the noise after compensation, the power of which is also unchanged. Since τp,q0−τp,qi=o−tp+rq−o−si−tp−rq−si/c, following the lemma in [31] and ignoring the quadratic and higher order terms,(15)o−tp+rq−o−si−tp−rq−si≈−2s¯iTd0+1d0s¯i−s¯iTd0d0tp+rq=−2ui+1d0pdtx˜i+qdry˜i+ξi,where s¯i=si−o, ui=s¯iTd0 denotes the radial range, which is the projection of vector s¯i on the direction of radar LOS d0, and vi=s¯i−s¯iTd0d0≐x˜i,y˜i,z˜iT is the residual vector of s¯i with respect to the radar LOS component. Then x˜i=viTx and y˜i=viTy are 2D cross-range coordinates, and ξi=viTt0+r0 is the initial range term with respect to starts of the transmit/receive ULAs. Conversely, the scatterer coordinates could be reconstructed with ui,x˜i,y˜i by(16)x^i=x˜i+uid01+x0,y^i=y˜i+uid02+y0,z^i=ui−x^i−x0d01−y^i−y0d02d03+z0.With the narrowband and far-filed assumption, that is,f0≫Bw and x˜i,y˜i,dt,dr≪d0, ignore cross-range offsets on multiple subcarriers and rewrite (14) as(17)sn,p,q≈∑i=1Iσ˜iexp−j2πnuiNΔuexpj2πpx˜iPΔx˜expj2πqy˜iQΔy˜=∑i=1Iσ˜iejnω1,iejpω2,iejqω3,i+w˜n,p,q,where σ˜i=σiexpj2π−2ui+ξi/λ denotes the phase modulated scattering amplitude, which is also constant for multiple transmit/receive paths,(18)Δu=c2NΔf=c2Bw,Δx˜=cd0f0Pdt=λ0d0Lt,Δy˜=cd0f0Qdr=λ0d0Lr,which are the resolutions of the radial/cross ranges, respectively, with λ0=c/f0 being the system wavelength and Lt=Pdt and Lr=Qdr denoting the transmit/receive array lengths. Therefore, sn,p,q is the uniformly spaced sampling in the 3D spatial spectrum with(19)ω1,i=−2πuiNΔu,ω2,i=2πx˜iPΔx˜,ω3,i=2πy˜iQΔy˜being the 3D spatial frequencies. In Algorithm 1, we summarize the proposed procedures of OFDM-MIMO radar echo preprocessing.Algorithm 1: OFDM-MIMO radar echo preprocessing.
Input: sampled echo without CP xqk, N phase code matrix Cn, center point compensation term βn,p,q0Output: spatial spectrum sampling sn,p,q(1)
Functionpreprocxqk(2)
for receive index q=0toQ−1do(3)
forOFDM symbol indexk=0toK−1(4)
Xn,qk←FFTxqk;(5)
end(6)
forsubcarrier indexn=0toN−1do(7)
xn,q=Xn,q0,…,Xn,qK−1T;(8)
yn,q=CnHxn,q=yn,0,q,…,yn,P−1,qT;(9)
fortransmit indexp=0toP−1do(10)
sn,p,q=βn,p,q0yn,p,q;(11)
end(12)
end(13)
end(14)
returnsn,p,q;(15)
end


## 2.1. OFDM-MIMO System and Echo Model

In the MIMO radar system,P transmit and Q receive antennas form two ULAs. As shown in Figure 1, the transmit array is parallel to the x-axis and uniformly spaced by dt; that is, coordinate of the p-th transmit antenna is tp=t0+pdt,0,0T=t0+pdtx, where p=0,…,P−1 indicates the transmit antenna index, t0=xt,0,yt,0,zt,0T denotes the start of transmit array which is adjacent to the coordinate system center or, and x=1,0,0T is the unit vector of x-axis. Similarly, the receive array is parallel to y-axis with the space being dr, q-th of which is denoted by rq=r0+qdry, where r0=xr,0,yr,0,zr,0T and y=0,1,0T.Figure 1
MIMO radar system configuration.Each transmit element in the radar system emits an OFDM signal modulated at the initial carrier frequencyf0. The baseband OFDM signal consists of N mutually orthogonal subcarriers with a uniform frequency spacing Δf; hence, the total bandwidth is Bw=NΔf. The n-th subcarrier is modulated with a K-bit coding sequence cn,p=an,p0,…,an,pK−1T. Complex orthogonal designs (CODs) [21–23] which fulfill(1)cn,pHcn,p′=1,p=p′,0,p≠p′,are applied [20] to build the orthogonal waveform set. Construct the matrix Cn=cn,0,…,cn,P−1∈ℂK×P and rewrite (1) as(2)CnHCn=IP.The bit widthtb meets the orthogonal condition that tbΔf=1. As a cyclic prefix (CP), tail of the original OFDM bit is added at the head of each transmitted OFDM symbol to remove the Intersymbol Interference (ISI) [7, 20]. Thus, the symbol duration is ts=tb+tc, with the CP duration being tc=αtb. The p-th baseband OFDM pulse is given by(3)upt=∑n=0N−1∑k=0K−1an,pkexpj2πnΔft−k+1tcrectt−ktsts−12,where p=0,…,P−1 and rectx=1,−1/2≤x<1/2,0,otherwise is the rectangular window. The baseband signal upt is then modulated at the initial carrier frequency f0 for transmission.Consider a far-filed target composed ofI ideal point-like scatterers, the echo model is given. The i-th scatterer coordinate is si=xi,yi,ziT, and the relative scattering amplitude is σi, which is a constant within an OFDM pulse duration and for multiple transmit/receive paths. Suppose the target velocity has been fully estimated and compensated before preprocessing [7], the Doppler effect is ignored in the echo model. At the q-th q=0,…,Q−1 receive antenna, let t=τmin+kts+tc+t0, with τmin being the start of the sampling window and t0∈0,tb. The relative round-way delay of the p,q-th transmit/receive path and i-th scatterer is(4)τp,qi=si−tp+rq−sic−τmin,where c is the velocity of light; then the echo after downconversion is expressed as(5)xqkt0=∑p=0P−1∑i=1Iσiexp−j2πf0τp,qi∑n=0N−1expj2πnΔft0−τp,qi⋅∑k′=0K−1an,pk′rectk−k′ts+tc+t0−τp,qits−12+wqkt0,where wqkt0 is the noise term, which is the zero-mean complex Additive White Gaussian Noise (AWGN). Assume τp,qi∈0,tc,∀p,q,i; then, the window function goes as(6)rectk−k′ts+tc+t0−τp,qits−12=δk−k′,where δ⋅ is the impulse function. Hence, there is no intersymbol interference (ISI) in the echo. Furthermore, set the sampling frequency as fs=Bw=NΔf; thus, t0=ltb/N,l=0,…,N−1; the discrete-time echo is written as(7)xl,qk=∑n=0N−1expj2πnlN∑p=0P−1an,pk∑i=1Iexp−j2πf0+nΔfτp,qiσi+wl,qk,where xl,qk and wl,qk denote the sampled echo and noise term, respectively.


## 2.2. Preprocessing Method

At theq-th receive element, for the k-th OFDM bit, compute the DFT of (7) towards the fast-time sampling index l; the frequency-domain echo is denoted by(8)Xn,qk=1N∑l=0N−1exp−j2πnlNxl,qk=∑p=0P−1an,pk∑i=1Iexp−j2πf0+nΔfτp,qiσi+Wn,qk,where Xn,qk is the n-th point of the DFT results and also denotes frequency-domain echo on the n-th subcarrier and Wn,qk is the frequency-domain noise. N subcarriers are separated without any Intercarrier Interference (ICI) in (8), but echoes from P transmit antennas are still mixed. For transmit separation and subcarrier decoding, for each subcarrier index n, we stack the DFT results of K bits in a column vector xn,q=Xn,q0,…,Xn,qK−1T and give the matrix form expression of xn,q as(9)xn,q=CnBn,qσ+wn,q∈ℂK×1,with Cn∈ℂK×P defined in (2) and(10)Bn,q=βn,0,q1⋯βn,0,qI⋮⋱⋮βn,P−1,q1⋯βn,P−1,qI∈ℂP×I,where the p,i-th element(11)βn,p,qi=exp−j2πf0+nΔfτp,qi,is the channel response of the i-th scatterer, on n-th subcarrier, for the p,q-th transmit/receive path, and σ=σ1,…,σIT∈ℂI×1 and wn,q=Wn,q0,…,Wn,qK−1T∈ℂK×1 are the scattering coefficient vector and the frequency-domain noise vector, respectively.Based on (1), left multiplying (9) with CnH yields(12)yn,q=CnHxn,q=Bn,qσ+wn,q′=yn,0,q,…,yn,P−1,qT∈ℂP×1,where yn,p,q denoted the demodulated echo and wn,q′ is the modified noise vector without change of power. Clearly, the P transmit waveforms are separated in yn,q.As is shown in Figure1, let o=x0,y0,z0T be the center of the imaging area, let d0=o be the radial distance, and let d0=o/d0 be the unit vector of the radar Line of Sight (LOS). For each OFDM subcarrier and each transmit/receive pair, we define the following compensation term:(13)βn,p,q0=expj2πf0+nΔfτp,q0,where τp,q0=o−tp+rq−o/c−τmin. yn,p,q is multiplied with βn,p,q0 for initial phase compensation, and(14)sn,p,q=βn,p,q0yn,p,q=∑i=1Iexpj2πf0+nΔfτp,q0−τp,qiσi+w˜n,p,q,where w˜n,p,q denotes the noise after compensation, the power of which is also unchanged. Since τp,q0−τp,qi=o−tp+rq−o−si−tp−rq−si/c, following the lemma in [31] and ignoring the quadratic and higher order terms,(15)o−tp+rq−o−si−tp−rq−si≈−2s¯iTd0+1d0s¯i−s¯iTd0d0tp+rq=−2ui+1d0pdtx˜i+qdry˜i+ξi,where s¯i=si−o, ui=s¯iTd0 denotes the radial range, which is the projection of vector s¯i on the direction of radar LOS d0, and vi=s¯i−s¯iTd0d0≐x˜i,y˜i,z˜iT is the residual vector of s¯i with respect to the radar LOS component. Then x˜i=viTx and y˜i=viTy are 2D cross-range coordinates, and ξi=viTt0+r0 is the initial range term with respect to starts of the transmit/receive ULAs. Conversely, the scatterer coordinates could be reconstructed with ui,x˜i,y˜i by(16)x^i=x˜i+uid01+x0,y^i=y˜i+uid02+y0,z^i=ui−x^i−x0d01−y^i−y0d02d03+z0.With the narrowband and far-filed assumption, that is,f0≫Bw and x˜i,y˜i,dt,dr≪d0, ignore cross-range offsets on multiple subcarriers and rewrite (14) as(17)sn,p,q≈∑i=1Iσ˜iexp−j2πnuiNΔuexpj2πpx˜iPΔx˜expj2πqy˜iQΔy˜=∑i=1Iσ˜iejnω1,iejpω2,iejqω3,i+w˜n,p,q,where σ˜i=σiexpj2π−2ui+ξi/λ denotes the phase modulated scattering amplitude, which is also constant for multiple transmit/receive paths,(18)Δu=c2NΔf=c2Bw,Δx˜=cd0f0Pdt=λ0d0Lt,Δy˜=cd0f0Qdr=λ0d0Lr,which are the resolutions of the radial/cross ranges, respectively, with λ0=c/f0 being the system wavelength and Lt=Pdt and Lr=Qdr denoting the transmit/receive array lengths. Therefore, sn,p,q is the uniformly spaced sampling in the 3D spatial spectrum with(19)ω1,i=−2πuiNΔu,ω2,i=2πx˜iPΔx˜,ω3,i=2πy˜iQΔy˜being the 3D spatial frequencies. In Algorithm 1, we summarize the proposed procedures of OFDM-MIMO radar echo preprocessing.Algorithm 1: OFDM-MIMO radar echo preprocessing.
Input: sampled echo without CP xqk, N phase code matrix Cn, center point compensation term βn,p,q0Output: spatial spectrum sampling sn,p,q(1)
Functionpreprocxqk(2)
for receive index q=0toQ−1do(3)
forOFDM symbol indexk=0toK−1(4)
Xn,qk←FFTxqk;(5)
end(6)
forsubcarrier indexn=0toN−1do(7)
xn,q=Xn,q0,…,Xn,qK−1T;(8)
yn,q=CnHxn,q=yn,0,q,…,yn,P−1,qT;(9)
fortransmit indexp=0toP−1do(10)
sn,p,q=βn,p,q0yn,p,q;(11)
end(12)
end(13)
end(14)
returnsn,p,q;(15)
end


## 3. OFDM-MIMO Radar High-Resolution Imaging via MD-UESPRIT

In this section, we use the Multidimensional Unitary-ESPRIT (MD-UESPRIT) algorithm to simultaneously estimate the 3D spatial frequenciesωd,id=1,2,3,i=1,…,I with the spatial spectrum echo sn,p,q and further reconstruct the scatterers. Since sn,p,q is only one snapshot available, we first employ a Multidimensional Spatial Smoothing (MD-SS) technique to restore the full-rank echo matrix. Including the backward-forward averaging, the unitary transformation then takes advantage of the complex conjugate of the initial data matrix and doubles the available snapshots. Finally, the 3D spatial frequencies are automatically paired, calculating the joint eigenvalue decomposition (JEVD) of the estimated matrix pencils.


### 3.1. Multidimensional Spatial Smoothing

MD-SS is introduced in this section to restore the full-rank matrix. LetD∈ℤ+ as the number of dimensions, stack the data sl1,…,lDd=1,…,D,ld=0,…,Ld−1 into an L1×⋯×LD array, and scan it by an R1×⋯×RD window. The sliding step is 1; hence, there are M=∏d=1DMd windows altogether, with Md=Ld−Rd+1. Set the scanning permutation as 1,…,D and denote m1,…,mD (md∈0,Md−1) as the initial index of the m=mD∏d=1D−1Md+mD−1∏d=1D−2Md+⋯+m1-th window; an R=∏d=1DRd-length vector sm is then formed with the elements in the window, which is organized by(20)sm=sm1,…,mD,…,sm1+R1−1,…,m3,sm1,m2+1,…,mD,…,sm1,m2+R2−1,…,mD,…,sm1+R1−1,…,mD+RD−1T=Aσm+nm,where A=AD⊕⋯⊕A1∈ℂR×I is the final steering matrix with R=∏d=1DRd and(21)Ad=1⋯1ejωd,1⋯ejωd,I⋮⋱⋮ejRd−1ωd,1⋯ejRd−1ωd,I∈ℂRd×Iis the steering matrix for the d-th dimension. σm=σ˜1,m,…,σ˜i,m,…,σ˜I,mTis the m-th signal vector with σ˜i,m=σ˜iexpj∑d=1Dmdωd,i, and nm is the noise vector within the m-th sliding window. The final matrix after spatial smoothing could be expressed as(22)S=A∑+N∈ℂR×M,with S=s0,…,sM−1, ∑=σ0,…,σM−1∈ℂI×M, and N=n0,…,nM−1∈ℂR×M. Specifically, in this 3D imaging application for OFDM-MIMO radar, set D=3, L1=N, L2=P, and L3=Q.


### 3.2. Unitary Transformation via Backward-Forward Averaging

Taking advantage of the complex conjugate of the data matrix, define the extended data matrix including the backward-forward averaging as(23)S˜=S,∏RS∗∏M∈ℂR×2M,which doubles the available number of snapshots, with ∏m=1⋰1m×m. It could be deduced that S˜ is centro-Hermitian that(24)∏RS˜∗∏2M=S˜.According to the theorem in [39], with the unitary matrix(25)Q2m+1=12Im0jIm0T20TΠm0−jΠm(the center column and row vector are dropped for even case), the following matrix(26)Z=TS=QRHS˜Q2M∈ℝR×2Mis real-valued. To simplify the computation, set S=S1TsS2TT, where submatrices S1 and S2 are of same size (the row vector sT is dropped when R is even); the real-valued matrix Z in (26) could be straightforwardly calculated as(27)Z=ℜS1+∏S2∗−ℑS1−∏S2∗2ℜsT−2ℑsTℑS1+∏S2∗ℜS1−∏S2∗,and only R×2M real additions are needed. Estimate the scatterer number I [52–54] (the specific procedures of which are ignored in this paper), and the signal subspace Es consists of the I dominant left singular vectors of Z. Correspondingly, the I dominant left singular vectors of S are QREs based on the corollary in [39].


### 3.3. Matrix Pencils Construction

D matrix pencils are constructed in this section exploiting the shift-invariance property of the uniform linear array manifold Ad in (20). For the d-th (d=1,…,D) steering matrix, we have(28)Jd,10AdΦd0=Jd,20Ad,where the selection matrices are Jd,10=IRd−10Rd−1×1, Jd,20=0Rd−1×1IRd−1=∏Rd−1Jd,10∏Rd, and Φd0=diagexpjωd,1,…,expjωd,I which is a diagonal matrix.Given the selection matrices(29)Jd,1=I∏d1=d+1DRd1⊗Jd,10⊗I∏d1=1d−1Rd1,Jd,2=I∏d1=d+1DRd1⊗Jd,20⊗I∏d1=1d−1Rd1,note that I∏d1=d+1DRd1=1 when d=D, and I∏d1=1d−1Rd1=1 when d=1. Based on the Kronecker product properties A⊗BC⊗D=AC⊗BD,A⊗B⊗C=A⊗B⊗C, the equations (30)Jd,1AΦd0=Jd,2A(31)Jd,2=∏R−RdJd,1∏Rare established. Since A and QREs span the same vector space, that is, there exists a nonsingular matrix T that A=QREsT, (30) could be rewritten as(32)Jd,1QREsTΦd0=Jd,2QREsT.Furthermore, based on (31),(33)QR−RdHJd,2QR=QR−RdH∏R−RdJd,1∏RQR=QR−RdHJd,1QR∗.Therefore, combining (32) and (33) yields(34)Kd,1EsΨd=Kd,2Es,where(35)Kd,1=ℜQR−RdHJd,2QR,Kd,2=ℑQR−RdHJd,2QR,and Ψd=TΦdT−1, with Φd=Diagtanωd,1/2,…,tanωd,I/2. Basically, the d-th pencil Ψd could be estimated via least square (LS) technique as(36)Ψ^d=Kd,1Es†Kd,2Es.


### 3.4. Location Estimation via JEVD

In [42], the matrix T is estimated by computing the EVD of the matrix Ψ^=∑d=1DαdΨ^d, where αd is the coefficient satisfying ∑d=1Dαd=1,αd>0. However, the diagonalizing performance relies heavily on the choice of coefficients αd, which are difficult to choose; meanwhile, the method has poor robustness with perturbations. In our work, the joint eigenvalue decomposition (JEVD) technique is adopted to jointly diagonalize the D matrix pencils, as the following way:(37)Ψ^d=T^Φ^dT^−1,where Ψ^d, Φ^d, and T^ are real. The uniqueness condition of the JEVD problem [55] is as follows:Construct the matrix(38)Ω=ϕ11⋯ϕ1D⋮⋱⋮ϕI1⋯ϕID∈ℂI×D,where ϕid is the i-th diagonal term in Φ^d. JEVD is unique if and only if all the rows of Ω are distinct.JEVD algorithms such as the “sh-rt” algorithm [47], JUST algorithm [48], JDTM algorithm [49], and JET algorithm [50] are Jacobi-like methods based on different matrix decomposition techniques. “sh-rt” has the poorest eigenvector estimation performance against noise [56]. JUST is slow to converge. JDTM requires only a few iterations before convergence but has high computation complexity per round. Conversely, JET decreases the numerical complexity per iteration but needs more iterations to converge. The JEVD based on a Taylor Expansion (JDTE) for real matrices is proposed in [51], which is also iterative but has no sweeping steps like Jacobi approaches. JDTE decreases the overall computational cost by limiting both the average number of iterations and the cost per iteration. Here, we choose JDTE for the estimation of T^ and Φ^d.In JDTE, an iteration procedure is used for JEVD. Initializing a nonsingular matrixH0 yields(39)Φd0=H0−1Ψ^dH0,whereas at each iteration, a matrix Hm,m=1,2,… is computed to update Φdm−1 as(40)Φdm=Hm−1Φdm−1Hm.The diagonalization parameter corresponding toHm is built as(41)cHm=∑d=1DΦdm−Λdm22=∑d=1DOdm22,where ⋅2 is the Euclidean norm of a real-valued matrix, Λdm=DiagdiagΦdm is a diagonal matrix consisting of the diagonal terms of Φdm, and Odm=Φdm−Λdm. Easily, if the algorithm achieves convergence after M iterations, then the final estimations are T^=∏m=0MHm and Φ^d=ΦdM.The indexm is afterwards dropped to simplify the notations. Rewrite the matrices H=I+Z, Φd=Λd+Od, and take them into (40) and (41). Suppose that we are close to the final solution; hence, Z,Od≪1, the first order Taylor approximation of (40), gives(42)I−Z−1ΦdI+Z≈I−ZΦdI+Z≈Φd−ZΦd+ΦdZ≈Λd+Od−ZΛd+ΛdZ,which is taken into (41) that the diagonalization parameter c towards Z is(43)cZ=∑d=1DOd−ZΛd+ΛdZ2=∑i=1I∑j=1j≠I−1∑d=1Doi,jd+λi,id−λj,jdzi,j2≐∑i=1I∑j=1j≠iIfzi,j,where oi,jd, λi,id, and zi,j are the corresponding element of matrices Od, Λd, and Z, respectively, and fzi,j=∑d=1Doi,jd+λi,id−λj,jdzi,j2.MinimizingcZ equals to minimizing each fzi,j. Let(44)∂fzi,j∂zi,j=2∑d=1Doi,jd+λi,id−λj,jdzi,jλi,id−λj,jd=0,and consequently, the analytical solution of zi,j is(45)zi,j=−∑d=1Doi,jdλi,id−λj,jd∑d=1Dλi,id−Λj,jd2.The algorithm is stopped if the iteration numberm reaches the maximum Mmax or the following condition is fulfilled:(46)Δc=cHm−cHm−1cHm−1≤ε.The performance relies significantly on the choice of initial inputH0. Here, we choose(47)H0=minevdΨ^dcevdΨ^d,where evd⋅ returns the matrix of eigenvectors and min⋅ denotes the minimum.With the final estimationΦ^d, extracting the diagonal elements yields the estimation of the spatial frequencies,(48)ω^d,i=2tan−1ϕ^i,id,where ϕ^i,id is the i-th diagonal element in Φ^d. So far, the spatial frequencies are estimated and automatically paired via the MD-UESPRIT algorithm.The scatterers’ locations in the radar coordinate system are then reconstructed with the estimated spatial frequencies based on (16) and (19), whereas the scattering intensities are estimated through the LS algorithm [39]. Procedures of the MD unitary-ESPRIT algorithm are summarized in Algorithm 2, where svd⋅,I returns the dominant I left singular vectors and JDTE⋅ is JEVD via JDTE algorithm, which is listed in Algorithm 3.Algorithm 2: Imaging via the MD unitary-ESPRIT algorithm.
Input: D-dimension spatial spectrum samplings sl1,…,lDd=1,…,D,ld=0,…Ld−1, window widths Rd, selection matrix Kd,1 and Kd,2, maximum iteration number Mmax, stopping threshold εOutput: spatial frequencies ω^1,i,…,ω^D,i(1)
Function MD-SS (sl1,…,lD,Rd)(2)
stacksl1,…,lD into an L1×⋯×LD array;(3)
ford=1toDdo(4)
Md=Ld−Rd+1;(5)
end(6)
m=0;(7)
formD=0toMD−1do(8)
⋯(9)
form1=0toM1−1do(10)
buildsm as in (20);(11)
m←m+1;(12)
end(13)
⋯(14)
end(15)
returnS=s0,…,sM−1;(16)
end(17)
FunctionMD-UESPRIT (S,Kd,1,Kd,2,Mmax,ε)(18)
Z=TS as in (27);(19)
estimate the numberI of scatterers [52–54];(20)
Es←svdZ,I;(21)
ford=1toD(22)
Ψ^d=Kd,1Es†Kd,2Es;(23)
end(24)
T^,Φ^d←JDTEΨ^d,Mmax,ε;(25)
ford=1toDdo(26)
fori=1toIdo(27)
ω^d,i=2tan−1ϕ^i,id;(28)
end(29)
end(30)
returnω^1,i,…,ω^D,i;(31)
endAlgorithm 3: JEVD via JDTE.
Input: matrix pencils Ψ^d, maximum iteration number Mmax, stopping threshold εOutput: eigenvectors matrix T^, diagonal matrices Φ^d(1)
FunctionJDTEΨ^d,Mmax,ε:(2)
m=0;(3)
initializeH0 with (47), Φdm with (39);(4)
computecH0 as (41), Δc=cH0;(5)
whilem≤Mmax and Δc≥εdo(6)
m=m+1;(7)
fori=1toIdo(8)
forj=1toIdo(9)
ifi≠jthen(10)
calculateZi,j with (44);(11)
end(12)
end(13)
end(14)
Hm=I+Z;(15)
ford=1toDdo(16)
updateΦdm−1 with (40)(17)
end(18)
calculateΔc using (41) and (46);(19)
end(20)
Φ^d=Φdm;(21)
T^=∏Hm;(22)
returnT^, Φ^d;(23)
end


### 3.5. Computational Complexity Analysis

We list the real multiplications required for each step in the imaging procedures in Table1. For the OFDM-MIMO echo preprocessing, the subcarriers separation via FFT needs 4×OQKNlnN flops, where N, K, and Q signify the number of the sampling points within one OFDM bit, the OFDM bits in one pulse, and the receive antennas, respectively; the decoding and center point compensation need 4×ONPQK+1 flops, with P denoting the transmit antenna number. For the imaging algorithm, the spatial smoothing and unitary transformation in (27) need no real multiplications, whereas the signal subspace and 3D spatial frequencies estimation involve only real-value computations after unitary transformation. Denote L as the iteration number of JDTE algorithm, and the total numerical complexities are shown in Table 1.Table 1
Computational complexity analysis.
StepReal multiplicationsPreprocessing4×OQKNlnN+QNPK+1Spatial smoothing—Unitary transformation—SVDOR3Matrix pencils computationO3R−∑Rd2RI+I2JEVD via JDTEL×OII−12D+1+2I3D+1


## 3.1. Multidimensional Spatial Smoothing

MD-SS is introduced in this section to restore the full-rank matrix. LetD∈ℤ+ as the number of dimensions, stack the data sl1,…,lDd=1,…,D,ld=0,…,Ld−1 into an L1×⋯×LD array, and scan it by an R1×⋯×RD window. The sliding step is 1; hence, there are M=∏d=1DMd windows altogether, with Md=Ld−Rd+1. Set the scanning permutation as 1,…,D and denote m1,…,mD (md∈0,Md−1) as the initial index of the m=mD∏d=1D−1Md+mD−1∏d=1D−2Md+⋯+m1-th window; an R=∏d=1DRd-length vector sm is then formed with the elements in the window, which is organized by(20)sm=sm1,…,mD,…,sm1+R1−1,…,m3,sm1,m2+1,…,mD,…,sm1,m2+R2−1,…,mD,…,sm1+R1−1,…,mD+RD−1T=Aσm+nm,where A=AD⊕⋯⊕A1∈ℂR×I is the final steering matrix with R=∏d=1DRd and(21)Ad=1⋯1ejωd,1⋯ejωd,I⋮⋱⋮ejRd−1ωd,1⋯ejRd−1ωd,I∈ℂRd×Iis the steering matrix for the d-th dimension. σm=σ˜1,m,…,σ˜i,m,…,σ˜I,mTis the m-th signal vector with σ˜i,m=σ˜iexpj∑d=1Dmdωd,i, and nm is the noise vector within the m-th sliding window. The final matrix after spatial smoothing could be expressed as(22)S=A∑+N∈ℂR×M,with S=s0,…,sM−1, ∑=σ0,…,σM−1∈ℂI×M, and N=n0,…,nM−1∈ℂR×M. Specifically, in this 3D imaging application for OFDM-MIMO radar, set D=3, L1=N, L2=P, and L3=Q.


## 3.2. Unitary Transformation via Backward-Forward Averaging

Taking advantage of the complex conjugate of the data matrix, define the extended data matrix including the backward-forward averaging as(23)S˜=S,∏RS∗∏M∈ℂR×2M,which doubles the available number of snapshots, with ∏m=1⋰1m×m. It could be deduced that S˜ is centro-Hermitian that(24)∏RS˜∗∏2M=S˜.According to the theorem in [39], with the unitary matrix(25)Q2m+1=12Im0jIm0T20TΠm0−jΠm(the center column and row vector are dropped for even case), the following matrix(26)Z=TS=QRHS˜Q2M∈ℝR×2Mis real-valued. To simplify the computation, set S=S1TsS2TT, where submatrices S1 and S2 are of same size (the row vector sT is dropped when R is even); the real-valued matrix Z in (26) could be straightforwardly calculated as(27)Z=ℜS1+∏S2∗−ℑS1−∏S2∗2ℜsT−2ℑsTℑS1+∏S2∗ℜS1−∏S2∗,and only R×2M real additions are needed. Estimate the scatterer number I [52–54] (the specific procedures of which are ignored in this paper), and the signal subspace Es consists of the I dominant left singular vectors of Z. Correspondingly, the I dominant left singular vectors of S are QREs based on the corollary in [39].


## 3.3. Matrix Pencils Construction

D matrix pencils are constructed in this section exploiting the shift-invariance property of the uniform linear array manifold Ad in (20). For the d-th (d=1,…,D) steering matrix, we have(28)Jd,10AdΦd0=Jd,20Ad,where the selection matrices are Jd,10=IRd−10Rd−1×1, Jd,20=0Rd−1×1IRd−1=∏Rd−1Jd,10∏Rd, and Φd0=diagexpjωd,1,…,expjωd,I which is a diagonal matrix.Given the selection matrices(29)Jd,1=I∏d1=d+1DRd1⊗Jd,10⊗I∏d1=1d−1Rd1,Jd,2=I∏d1=d+1DRd1⊗Jd,20⊗I∏d1=1d−1Rd1,note that I∏d1=d+1DRd1=1 when d=D, and I∏d1=1d−1Rd1=1 when d=1. Based on the Kronecker product properties A⊗BC⊗D=AC⊗BD,A⊗B⊗C=A⊗B⊗C, the equations (30)Jd,1AΦd0=Jd,2A(31)Jd,2=∏R−RdJd,1∏Rare established. Since A and QREs span the same vector space, that is, there exists a nonsingular matrix T that A=QREsT, (30) could be rewritten as(32)Jd,1QREsTΦd0=Jd,2QREsT.Furthermore, based on (31),(33)QR−RdHJd,2QR=QR−RdH∏R−RdJd,1∏RQR=QR−RdHJd,1QR∗.Therefore, combining (32) and (33) yields(34)Kd,1EsΨd=Kd,2Es,where(35)Kd,1=ℜQR−RdHJd,2QR,Kd,2=ℑQR−RdHJd,2QR,and Ψd=TΦdT−1, with Φd=Diagtanωd,1/2,…,tanωd,I/2. Basically, the d-th pencil Ψd could be estimated via least square (LS) technique as(36)Ψ^d=Kd,1Es†Kd,2Es.


## 3.4. Location Estimation via JEVD

In [42], the matrix T is estimated by computing the EVD of the matrix Ψ^=∑d=1DαdΨ^d, where αd is the coefficient satisfying ∑d=1Dαd=1,αd>0. However, the diagonalizing performance relies heavily on the choice of coefficients αd, which are difficult to choose; meanwhile, the method has poor robustness with perturbations. In our work, the joint eigenvalue decomposition (JEVD) technique is adopted to jointly diagonalize the D matrix pencils, as the following way:(37)Ψ^d=T^Φ^dT^−1,where Ψ^d, Φ^d, and T^ are real. The uniqueness condition of the JEVD problem [55] is as follows:Construct the matrix(38)Ω=ϕ11⋯ϕ1D⋮⋱⋮ϕI1⋯ϕID∈ℂI×D,where ϕid is the i-th diagonal term in Φ^d. JEVD is unique if and only if all the rows of Ω are distinct.JEVD algorithms such as the “sh-rt” algorithm [47], JUST algorithm [48], JDTM algorithm [49], and JET algorithm [50] are Jacobi-like methods based on different matrix decomposition techniques. “sh-rt” has the poorest eigenvector estimation performance against noise [56]. JUST is slow to converge. JDTM requires only a few iterations before convergence but has high computation complexity per round. Conversely, JET decreases the numerical complexity per iteration but needs more iterations to converge. The JEVD based on a Taylor Expansion (JDTE) for real matrices is proposed in [51], which is also iterative but has no sweeping steps like Jacobi approaches. JDTE decreases the overall computational cost by limiting both the average number of iterations and the cost per iteration. Here, we choose JDTE for the estimation of T^ and Φ^d.In JDTE, an iteration procedure is used for JEVD. Initializing a nonsingular matrixH0 yields(39)Φd0=H0−1Ψ^dH0,whereas at each iteration, a matrix Hm,m=1,2,… is computed to update Φdm−1 as(40)Φdm=Hm−1Φdm−1Hm.The diagonalization parameter corresponding toHm is built as(41)cHm=∑d=1DΦdm−Λdm22=∑d=1DOdm22,where ⋅2 is the Euclidean norm of a real-valued matrix, Λdm=DiagdiagΦdm is a diagonal matrix consisting of the diagonal terms of Φdm, and Odm=Φdm−Λdm. Easily, if the algorithm achieves convergence after M iterations, then the final estimations are T^=∏m=0MHm and Φ^d=ΦdM.The indexm is afterwards dropped to simplify the notations. Rewrite the matrices H=I+Z, Φd=Λd+Od, and take them into (40) and (41). Suppose that we are close to the final solution; hence, Z,Od≪1, the first order Taylor approximation of (40), gives(42)I−Z−1ΦdI+Z≈I−ZΦdI+Z≈Φd−ZΦd+ΦdZ≈Λd+Od−ZΛd+ΛdZ,which is taken into (41) that the diagonalization parameter c towards Z is(43)cZ=∑d=1DOd−ZΛd+ΛdZ2=∑i=1I∑j=1j≠I−1∑d=1Doi,jd+λi,id−λj,jdzi,j2≐∑i=1I∑j=1j≠iIfzi,j,where oi,jd, λi,id, and zi,j are the corresponding element of matrices Od, Λd, and Z, respectively, and fzi,j=∑d=1Doi,jd+λi,id−λj,jdzi,j2.MinimizingcZ equals to minimizing each fzi,j. Let(44)∂fzi,j∂zi,j=2∑d=1Doi,jd+λi,id−λj,jdzi,jλi,id−λj,jd=0,and consequently, the analytical solution of zi,j is(45)zi,j=−∑d=1Doi,jdλi,id−λj,jd∑d=1Dλi,id−Λj,jd2.The algorithm is stopped if the iteration numberm reaches the maximum Mmax or the following condition is fulfilled:(46)Δc=cHm−cHm−1cHm−1≤ε.The performance relies significantly on the choice of initial inputH0. Here, we choose(47)H0=minevdΨ^dcevdΨ^d,where evd⋅ returns the matrix of eigenvectors and min⋅ denotes the minimum.With the final estimationΦ^d, extracting the diagonal elements yields the estimation of the spatial frequencies,(48)ω^d,i=2tan−1ϕ^i,id,where ϕ^i,id is the i-th diagonal element in Φ^d. So far, the spatial frequencies are estimated and automatically paired via the MD-UESPRIT algorithm.The scatterers’ locations in the radar coordinate system are then reconstructed with the estimated spatial frequencies based on (16) and (19), whereas the scattering intensities are estimated through the LS algorithm [39]. Procedures of the MD unitary-ESPRIT algorithm are summarized in Algorithm 2, where svd⋅,I returns the dominant I left singular vectors and JDTE⋅ is JEVD via JDTE algorithm, which is listed in Algorithm 3.Algorithm 2: Imaging via the MD unitary-ESPRIT algorithm.
Input: D-dimension spatial spectrum samplings sl1,…,lDd=1,…,D,ld=0,…Ld−1, window widths Rd, selection matrix Kd,1 and Kd,2, maximum iteration number Mmax, stopping threshold εOutput: spatial frequencies ω^1,i,…,ω^D,i(1)
Function MD-SS (sl1,…,lD,Rd)(2)
stacksl1,…,lD into an L1×⋯×LD array;(3)
ford=1toDdo(4)
Md=Ld−Rd+1;(5)
end(6)
m=0;(7)
formD=0toMD−1do(8)
⋯(9)
form1=0toM1−1do(10)
buildsm as in (20);(11)
m←m+1;(12)
end(13)
⋯(14)
end(15)
returnS=s0,…,sM−1;(16)
end(17)
FunctionMD-UESPRIT (S,Kd,1,Kd,2,Mmax,ε)(18)
Z=TS as in (27);(19)
estimate the numberI of scatterers [52–54];(20)
Es←svdZ,I;(21)
ford=1toD(22)
Ψ^d=Kd,1Es†Kd,2Es;(23)
end(24)
T^,Φ^d←JDTEΨ^d,Mmax,ε;(25)
ford=1toDdo(26)
fori=1toIdo(27)
ω^d,i=2tan−1ϕ^i,id;(28)
end(29)
end(30)
returnω^1,i,…,ω^D,i;(31)
endAlgorithm 3: JEVD via JDTE.
Input: matrix pencils Ψ^d, maximum iteration number Mmax, stopping threshold εOutput: eigenvectors matrix T^, diagonal matrices Φ^d(1)
FunctionJDTEΨ^d,Mmax,ε:(2)
m=0;(3)
initializeH0 with (47), Φdm with (39);(4)
computecH0 as (41), Δc=cH0;(5)
whilem≤Mmax and Δc≥εdo(6)
m=m+1;(7)
fori=1toIdo(8)
forj=1toIdo(9)
ifi≠jthen(10)
calculateZi,j with (44);(11)
end(12)
end(13)
end(14)
Hm=I+Z;(15)
ford=1toDdo(16)
updateΦdm−1 with (40)(17)
end(18)
calculateΔc using (41) and (46);(19)
end(20)
Φ^d=Φdm;(21)
T^=∏Hm;(22)
returnT^, Φ^d;(23)
end


## 3.5. Computational Complexity Analysis

We list the real multiplications required for each step in the imaging procedures in Table1. For the OFDM-MIMO echo preprocessing, the subcarriers separation via FFT needs 4×OQKNlnN flops, where N, K, and Q signify the number of the sampling points within one OFDM bit, the OFDM bits in one pulse, and the receive antennas, respectively; the decoding and center point compensation need 4×ONPQK+1 flops, with P denoting the transmit antenna number. For the imaging algorithm, the spatial smoothing and unitary transformation in (27) need no real multiplications, whereas the signal subspace and 3D spatial frequencies estimation involve only real-value computations after unitary transformation. Denote L as the iteration number of JDTE algorithm, and the total numerical complexities are shown in Table 1.Table 1
Computational complexity analysis.
StepReal multiplicationsPreprocessing4×OQKNlnN+QNPK+1Spatial smoothing—Unitary transformation—SVDOR3Matrix pencils computationO3R−∑Rd2RI+I2JEVD via JDTEL×OII−12D+1+2I3D+1


## 4. Simulation Results

### 4.1. Imaging Results Comparison

In this section, we use the simulated data for the 127-point scatterer model of a Boeing 777 airplane to verify the performance of the proposed 3D imaging method. The target scatterer model is depicted in Figure2, including the projection in each plane and the 3D image. For ease of presentation and comparison, the scattering amplitudes of the scatterers are set uniformly distributed, with the dynamical range being 10 dB. The center o of the imaging area is o=1767.8,3061.9,3535.5T m, with the initial radial range d0 = 6  km.Figure 2
Scatterer model of a Boeing 777 plane.The MIMO radar system consists ofP=8-transmit ULA with element spacing dt = 4 m and Q=16-receive ULA with element spacing dr = 2 m. Each transmit antenna radiates the OFDM waveform isotropically in space. The transmitted OFDM waveform consists of N=128 subcarriers, the frequency interval Δf = 1 MHz, the OFDM bit/CP/symbol durations are tb=1µs, tc=0.5µs, and ts=1.5µs, respectively, and K=16; hence, the bandwidth is Bw=NΔf=128MHz and the pulse width is tp=Kts=24µs . After preprocessing, the size of the spatial spectrum data array is 128×8×16.Figure3 compares the imaging results obtained, respectively, from the 3D-FFT, MD-MUSIC, and MD-ESPRIT in [42] and the proposed MD-UESPRIT algorithm. Each image is normalized by the values of the strongest points. All the tests are conducted with SNR=15dB, where SNR is defined as(49)SNR=20lg∑q=0Q−1Es,q∑q=0Q−1En,q,with Es,q and En,q being the signal/noise power, respectively, at the q-th receiver.Figure 3
Imaging results of a Boeing 777 plane.
(a)(b)(c)(d)The 3D-FFT result is computed with the zero-padded data array of size256×64×128. In Figure 3(a), the first three subfigures show the superimposed 2D images in u−x˜, u−y˜, and x˜−y˜ planes, respectively. The grid widths of the radial range and 2D cross ranges are u0=0.5859m, x˜0=0.7031m, and y˜0=0.7031m. The resolution of FFT spectrum is so low that we could hardly identify the target from the results. It can be seen in the last subfigure that the target could hardly be identified from the reconstructed image.Set the smoothing window sizes as64,4,8; the snapshots after MD-SS are M=∏d=1DMd=2925. The MD-MUSIC algorithm is adopted in Figure 3(b) with the same searching steps as in Figure 3(a). Mainlobe width and sidelobe level decrease significantly in the MUSIC spectrum compared with the FFT spectrum. In the last subfigure, peaks in the MUSIC spectrum are selected for target reconstruction, whereas the scattering amplitudes are not directly equal to the spectrum value but are also calculated via the LS algorithm. However, the reconstruction error is still high; meanwhile, it consumes a long time for 3D searching.As can be seen in Figures3(c) and 3(d), the imaging results of the ESPRIT-class methods are the specific locations and intensities of distinct scatterers without any sidelobes. The proposed method in Figure 3(d) has best imaging quality with high reconstruction accuracy and takes less time than MD-MUSIC and MD-ESPRIT.The performance of the proposed algorithm is further demonstrated with 2D EM simulated data of an Airbus A320 airplane. The transmitted OFDM waveform consists ofN=256 subcarriers with frequency interval Δf=1MHz; the bandwidth is hence Bw=256MHz with a carrier frequency of f0 = 10 GHz. 64-element ULA is adopted for receiving; the incident angle for each receiver is increased by 0.009°. The spatial spectrum echo is 256×64 after preprocessing. Set the initial incident/azimuth angle as (90°, 0°) and (90°, 45°); the normalized imaging results are depicted in Figures 4(a) and 4(b), respectively. In each subfigure, the results of 2D-FFT, 2D-MUSIC, and 2D-ESPRIT are also shown as a comparison. Figure 4 indicates that the high-resolution spatial frequency estimation methods such as MUSIC and ESPRIT provide better resolution images compared with the Fourier-based method. Among them, MUSIC is time-consuming with multidimensional searching. The proposed method has better reconstruction performance with fewer fake peaks and a dynamical range closer to that of the FFT method.Figure 4
Imaging results of an Airbus A320 plane.
(a)(b)


### 4.2. Reconstruction Accuracy Analysis

The reconstruction accuracy of 3D-UESPRIT estimation is analyzed in this section, corresponding toSNR, system parameters, and errors. As shown in Figure 5, we choose a target with fewer scatterers for ease of analysis with initial radial distance being d0=4km. The size of scatterers indicates the scattering amplitude. For each SNR and parameters combination, we carry out M=200 Monte Carlo simulations. The Root-Mean-Square Error (RMSE) is computed to measure the reconstruction precision and is expressed as(50)RMSE1=1I∑i=1I1M∑m=1Mx^i−xi2,RMSE2=1I∑i=1I1M∑m=1My^i−yi2,RMSE3=1I∑i=1I1M∑m=1Mz^i−zi2,RMSE=1D∑d=1DRMSEd.Figure 5
Simulated 3D target model with the red rectangular representing the imaging center.In Figure6, SNR is varied from −18 dB to 18 dB with 3 dB as the step. The four subfigures in Figure 6 show that RMSE decreases with the increase of SNR, especially when SNR≤0dB. Keeping the system parameters, the same with those in Figure 3(d), we utilize different window sizes for MD-SS in Figure 6(a). For each dimension, the window size decides the order of the final steering matrix and the incoherence property of multiple scatterers in matrix S. Figure 6(a) shows that, within a certain range, the estimation precision improves with the increase of the total matrix order R=∏d=1DRd, while after that, the estimation accuracy decreases on the contrary. This is because when R exceeds some value, the coherence among scatterers increases that does not meet the criterion of the ESPRIT-class algorithm. Meanwhile, the identifiability of the algorithm [39] is defined as(51)Imax=minR1−1∏d=2DRd,R1R2−1∏d=2DRd,…,∏d=1D−1RD−1,2M.Therefore, considering both the estimation accuracy and the maximum estimable scatterer number, we choose Rd=Ld/2d=1,…,D for the MD-UESPRIT algorithms.Figure 6
Reconstruction accuracy corresponding toSNR and parameters.
(a)(b)(c)(d)Suppose other parameters be unchanged as in Figure3(d); the subcarrier/transmit/receive array parameters are analyzed, respectively, in Figures 6(b) and 6(d). For each array, with the same element number, RMSE decreases with the increase of interelement spacings (Δf, dt, dr), which is equivalent to the increase of bandwidth Bw and array lengths, while the maximum unambiguous intervals are(52)umax=αc2Δf,x˜max=λd0dt,y˜max=λd0dr,which are inversely proportional to the interelement spacings. Hence, the interelement spacings cannot increase without limitation, and the array lengths are subsequently restricted. With the bandwidth/array lengths unchanged, the estimation accuracy ameliorates with the increase of element numbers (N,P,Q). However, both the system cost and computational complexity improve with more elements. Thus, we have to reach a compromise between the element number and interelement spacing when designing the system parameters.Last but not least, reconstruction accuracy is analyzed concerning the system errors. Position errors in the transmit/receive antenna array are considered. The coordinate ofp-th transmit antenna is set as tp=p+Δpdt,0,0T, where Δp is a random value uniformly distributed in −σd/2,σd/2. Similarly, for q-th receive antenna, rq=0,q+Δqdr,0T, and Δq∼U−σd/2,σd/2. For simplification, we also model the Doppler frequency compensation residual error on each subcarrier, a uniformly distributed random, i.e., the n-th subcarrier frequency after Doppler compensation is fn=f0+n+ΔnΔf and Δn∼U−σf/2,σf/2. RMSE curves with respect to σd and σf are depicted in Figure 7, with SNR=5dB. For one thing, the estimation accuracy is sensitive to antenna position error; hence, effective antenna array calibration method is required before imaging. For another, with the system setting in the tests, the method has a tolerance of Doppler compensation error. However, we ignore the Intercarrier Inference (ICI) for subcarrier separation and decoding in our Doppler perturbation model. Therefore, more detailed considerations towards Doppler compensation error are required in the future work.Figure 7
Reconstruction accuracy corresponding to antenna arrays and frequency errors.


## 4.1. Imaging Results Comparison

In this section, we use the simulated data for the 127-point scatterer model of a Boeing 777 airplane to verify the performance of the proposed 3D imaging method. The target scatterer model is depicted in Figure2, including the projection in each plane and the 3D image. For ease of presentation and comparison, the scattering amplitudes of the scatterers are set uniformly distributed, with the dynamical range being 10 dB. The center o of the imaging area is o=1767.8,3061.9,3535.5T m, with the initial radial range d0 = 6  km.Figure 2
Scatterer model of a Boeing 777 plane.The MIMO radar system consists ofP=8-transmit ULA with element spacing dt = 4 m and Q=16-receive ULA with element spacing dr = 2 m. Each transmit antenna radiates the OFDM waveform isotropically in space. The transmitted OFDM waveform consists of N=128 subcarriers, the frequency interval Δf = 1 MHz, the OFDM bit/CP/symbol durations are tb=1µs, tc=0.5µs, and ts=1.5µs, respectively, and K=16; hence, the bandwidth is Bw=NΔf=128MHz and the pulse width is tp=Kts=24µs . After preprocessing, the size of the spatial spectrum data array is 128×8×16.Figure3 compares the imaging results obtained, respectively, from the 3D-FFT, MD-MUSIC, and MD-ESPRIT in [42] and the proposed MD-UESPRIT algorithm. Each image is normalized by the values of the strongest points. All the tests are conducted with SNR=15dB, where SNR is defined as(49)SNR=20lg∑q=0Q−1Es,q∑q=0Q−1En,q,with Es,q and En,q being the signal/noise power, respectively, at the q-th receiver.Figure 3
Imaging results of a Boeing 777 plane.
(a)(b)(c)(d)The 3D-FFT result is computed with the zero-padded data array of size256×64×128. In Figure 3(a), the first three subfigures show the superimposed 2D images in u−x˜, u−y˜, and x˜−y˜ planes, respectively. The grid widths of the radial range and 2D cross ranges are u0=0.5859m, x˜0=0.7031m, and y˜0=0.7031m. The resolution of FFT spectrum is so low that we could hardly identify the target from the results. It can be seen in the last subfigure that the target could hardly be identified from the reconstructed image.Set the smoothing window sizes as64,4,8; the snapshots after MD-SS are M=∏d=1DMd=2925. The MD-MUSIC algorithm is adopted in Figure 3(b) with the same searching steps as in Figure 3(a). Mainlobe width and sidelobe level decrease significantly in the MUSIC spectrum compared with the FFT spectrum. In the last subfigure, peaks in the MUSIC spectrum are selected for target reconstruction, whereas the scattering amplitudes are not directly equal to the spectrum value but are also calculated via the LS algorithm. However, the reconstruction error is still high; meanwhile, it consumes a long time for 3D searching.As can be seen in Figures3(c) and 3(d), the imaging results of the ESPRIT-class methods are the specific locations and intensities of distinct scatterers without any sidelobes. The proposed method in Figure 3(d) has best imaging quality with high reconstruction accuracy and takes less time than MD-MUSIC and MD-ESPRIT.The performance of the proposed algorithm is further demonstrated with 2D EM simulated data of an Airbus A320 airplane. The transmitted OFDM waveform consists ofN=256 subcarriers with frequency interval Δf=1MHz; the bandwidth is hence Bw=256MHz with a carrier frequency of f0 = 10 GHz. 64-element ULA is adopted for receiving; the incident angle for each receiver is increased by 0.009°. The spatial spectrum echo is 256×64 after preprocessing. Set the initial incident/azimuth angle as (90°, 0°) and (90°, 45°); the normalized imaging results are depicted in Figures 4(a) and 4(b), respectively. In each subfigure, the results of 2D-FFT, 2D-MUSIC, and 2D-ESPRIT are also shown as a comparison. Figure 4 indicates that the high-resolution spatial frequency estimation methods such as MUSIC and ESPRIT provide better resolution images compared with the Fourier-based method. Among them, MUSIC is time-consuming with multidimensional searching. The proposed method has better reconstruction performance with fewer fake peaks and a dynamical range closer to that of the FFT method.Figure 4
Imaging results of an Airbus A320 plane.
(a)(b)


## 4.2. Reconstruction Accuracy Analysis

The reconstruction accuracy of 3D-UESPRIT estimation is analyzed in this section, corresponding toSNR, system parameters, and errors. As shown in Figure 5, we choose a target with fewer scatterers for ease of analysis with initial radial distance being d0=4km. The size of scatterers indicates the scattering amplitude. For each SNR and parameters combination, we carry out M=200 Monte Carlo simulations. The Root-Mean-Square Error (RMSE) is computed to measure the reconstruction precision and is expressed as(50)RMSE1=1I∑i=1I1M∑m=1Mx^i−xi2,RMSE2=1I∑i=1I1M∑m=1My^i−yi2,RMSE3=1I∑i=1I1M∑m=1Mz^i−zi2,RMSE=1D∑d=1DRMSEd.Figure 5
Simulated 3D target model with the red rectangular representing the imaging center.In Figure6, SNR is varied from −18 dB to 18 dB with 3 dB as the step. The four subfigures in Figure 6 show that RMSE decreases with the increase of SNR, especially when SNR≤0dB. Keeping the system parameters, the same with those in Figure 3(d), we utilize different window sizes for MD-SS in Figure 6(a). For each dimension, the window size decides the order of the final steering matrix and the incoherence property of multiple scatterers in matrix S. Figure 6(a) shows that, within a certain range, the estimation precision improves with the increase of the total matrix order R=∏d=1DRd, while after that, the estimation accuracy decreases on the contrary. This is because when R exceeds some value, the coherence among scatterers increases that does not meet the criterion of the ESPRIT-class algorithm. Meanwhile, the identifiability of the algorithm [39] is defined as(51)Imax=minR1−1∏d=2DRd,R1R2−1∏d=2DRd,…,∏d=1D−1RD−1,2M.Therefore, considering both the estimation accuracy and the maximum estimable scatterer number, we choose Rd=Ld/2d=1,…,D for the MD-UESPRIT algorithms.Figure 6
Reconstruction accuracy corresponding toSNR and parameters.
(a)(b)(c)(d)Suppose other parameters be unchanged as in Figure3(d); the subcarrier/transmit/receive array parameters are analyzed, respectively, in Figures 6(b) and 6(d). For each array, with the same element number, RMSE decreases with the increase of interelement spacings (Δf, dt, dr), which is equivalent to the increase of bandwidth Bw and array lengths, while the maximum unambiguous intervals are(52)umax=αc2Δf,x˜max=λd0dt,y˜max=λd0dr,which are inversely proportional to the interelement spacings. Hence, the interelement spacings cannot increase without limitation, and the array lengths are subsequently restricted. With the bandwidth/array lengths unchanged, the estimation accuracy ameliorates with the increase of element numbers (N,P,Q). However, both the system cost and computational complexity improve with more elements. Thus, we have to reach a compromise between the element number and interelement spacing when designing the system parameters.Last but not least, reconstruction accuracy is analyzed concerning the system errors. Position errors in the transmit/receive antenna array are considered. The coordinate ofp-th transmit antenna is set as tp=p+Δpdt,0,0T, where Δp is a random value uniformly distributed in −σd/2,σd/2. Similarly, for q-th receive antenna, rq=0,q+Δqdr,0T, and Δq∼U−σd/2,σd/2. For simplification, we also model the Doppler frequency compensation residual error on each subcarrier, a uniformly distributed random, i.e., the n-th subcarrier frequency after Doppler compensation is fn=f0+n+ΔnΔf and Δn∼U−σf/2,σf/2. RMSE curves with respect to σd and σf are depicted in Figure 7, with SNR=5dB. For one thing, the estimation accuracy is sensitive to antenna position error; hence, effective antenna array calibration method is required before imaging. For another, with the system setting in the tests, the method has a tolerance of Doppler compensation error. However, we ignore the Intercarrier Inference (ICI) for subcarrier separation and decoding in our Doppler perturbation model. Therefore, more detailed considerations towards Doppler compensation error are required in the future work.Figure 7
Reconstruction accuracy corresponding to antenna arrays and frequency errors.


## 5. Threats to Validity

There are two threats to validity as follows.The first threat is that we ignore the Doppler effect in our echo model by assuming the target velocity is zero or fully compensated. However, this condition is hard to achieve in reality. For one thing, we model the Doppler frequency shift in 4.2 as a random error in the subcarrier frequency array to measure the reconstruction accuracy relative to the Doppler effect. For another, the Doppler offset causes the interferences among subcarriers and transmitters [8], which is not considered in this paper.The second threat is that the estimated signal subspace suffers from significant deviation with perturbations (including noise and array manifold error).The proposed method is not robust enough with SNR and system errors.


## 6. Future Work

As future work, we are planning to investigate some modifications of the proposed algorithm to make reconstruction more robust to SNR and system errors, which is significant for real applications. Secondly, we will build a complete signal model, including the Doppler effect, and propose an optimized transmitting scheme to mitigate the Doppler offset. Last but not least, we will build a physical experiment system to verify our method.


## 7. Conclusion

We present in this paper a high-resolution 3D imaging method for a MIMO-OFDM radar system with a single pulse. Orthogonal transmit and complete separation are achieved with OFDM structure and the preprocessing steps. We adopt the MD-UESPRIT algorithm to reconstruct the coordinates of scatterers simultaneously. The simulation results indicate that our proposed imaging algorithm has better reconstruction results compared with FFT, MD-MUSIC, and MD-ESPRIT. Parameters selection formulas about estimation accuracy, maximum unambiguous range, identifiability, and computational complexity are also provided in this paper.



---

*Source: 2308389-2020-06-27.xml*