An Unambiguous Super-Resolution Algorithm for TDM-MIMO-SAR 3D Imaging Applications on Fast-Moving Platforms

1. Introduction

In recent years, the Multiple Input Multiple Output (MIMO) technique, which achieves high Degree of Freedom (DOF) for radar systems, has drawn an increasing amount of attention in various applications, such as Autonomous Driving [1], Ground Penetrating Radar [2], Remote Sensing [3], Through-wall Imaging [4], etc. Orthogonal waveforms are radiated by transmitters (TXs), thus a MIMO system can separate different signals transmitted by different TXs. For example, a MIMO radar which has $N_T$ transmitters and $N_R$ receivers can generate $N_T\times N_R$ virtual antennas. There are three multiplexing schemes to obtain waveform diversity, namely Time-Division Multiplexing (TDM), Frequency-Division Multiplexing (FDM) and Code-Division Multiplexing (CDM). Due to the advantages of low hardware complexity and high signal separation, TDM-MIMO system is easy to implement, thus it has a wide range of applications [5]. A large amount of research has focused on the algorithms and implementations for TDM-MIMO radar [6,7,8]. However, a TDM-MIMO system has an inherent drawback which limits its performance. As transmitters take turns to transmit their signal, phase errors occur in nonstationary scenarios [9]. In addition, more transmitters mean lower Pulse Repetition Frequency (PRF) and reduced maximum velocity that can be detected [5].

Synthetic Aperture Radar (SAR) is an important technique for remote sensing and radar imaging. With the help of MIMO technology, a new architecture called MIMO-SAR not only increases spatial resolution [10], but also provides novel imaging capability such as 3D imaging [11] and Ground-Moving Target Indication (GMTI) [12]. Due to the great potential of MIMO-SAR architecture, many researchers have contributed to various aspects of MIMO-SAR. In [10], it is shown that the mere use of orthogonal waveforms is insufficient for most SAR applications, and thus a new class of waveform is introduced. Yanik et al. addressed the issues of the design process of MIMO-SAR systems and presents three versions of MIMO-SAR testbeds with different costs and performances [13]. As in many cases, such as Unmanned Aerial Vehicle (UAV) imaging, irregular scanning geometries and computational complexity limit the capability of MIMO-SAR systems, Smith et al. proposed a new algorithm to leverage SAR and MIMO while achieving high resolution with high efficiency near field MIMO-SAR imaging [14]. The effect of virtual array approximation in near field applications is studied by Hu et al., and a new design criterion is presented [15]. An efficient image reconstruction algorithm which avoids time-consuming wavenumber accumulation by decomposing the range-wavenumber based on binomial approximation theory is proposed for MIMO-SAR 3D imaging [16]. A beat frequency division (BFD) frequency-modulated continuous wave (FMCW) MIMO-SAR system and the corresponding signal process method which are suitable for UAV environments are presented [17]. With appropriate transformation, MIMO-SAR 3D imaging can be achieved by inverse Fourier transform (IFT) [18].

Apart from the imaging algorithms and array design of MIMO-SAR, it is the realization of MIMO systems that significantly affect both the performance and costs of MIMO-SAR. There are numerous studies discussing the problem of MIMO-SAR associated with a specific MIMO realization. Jin et al. proposed a novel MIMO-SAR scheme based on segmented-phase-code waveform which can relieve the short-term shift-orthogonality condition without losing imaging performance [19]. In [20], a MIMO-SAR system using chirp diverse waveform is discussed to achieve high resolution for wide-swath remote sensing, as chirp diverse waveform can easily provide high bandwidth which is necessary for this radar imaging. It is acknowledged that a large time-bandwidth product, low cross-correlation interferences and a low peak-average ratio have significant importance for MIMO radar imaging. To reach these aims, an orthogonal frequency division multiplexing (OFDM) chirp waveform with random matrix modulation is proposed for MIMO-SAR applications [21]. As mentioned before, TDM-MIMO enjoys many advantages with respect to hardware complexity and cost. A forward-looking MIMO-SAR enhanced imaging scheme is presented, and its performance is validated with a TDM-MIMO radar [22]. In [23], the authors proposed a hierarchical high-resolution MIMO-SAR imaging algorithm which is derived in case of TDM-MIMO realization and the algorithm is tested on a TDM-MIMO system. There are also several researches focused on near-field imaging [24], Ground Based SAR(GB-SAR) applications [25]. To tackle the problem of ghost targets caused by sparse MIMO array, an algorithm called Sequential Spatial Masking (SSM) is proposed for forward-looking TDM-MIMO-SAR imaging [26]. As discussed above, the nature of TDM-MIMO brings challenges to imaging algorithms in nonstationary scenarios because different transmitters are activated separately. Since Synthetic Aperture Radar (SAR) systems rely on the relative movement between platforms and targets, TDM-MIMO brings difficulties to using TDM-MIMO for SAR imaging [27]. Besides, to achieve higher resolution with TDM-MIMO array, it is often required to have more transmitters, which leads to a lower PRF. However, if a MIMO-SAR system is installed on a fast-moving platform, such as a car or an aircraft, a low PRF means grating lobes and ghost targets appears in SAR imaging process. Notably, it is the time-division nature of TDM-MIMO that brings obstacles to its applications MIMO-SAR 3D imaging on a fast-moving platform.

To further improve imaging capabilities of MIMO-SAR technique, the super-resolution algorithms are also applied and modified for MIMO-SAR imaging. Traditional imaging methods, which lack the super-resolution ability, are well-researched in MIMO imaging scenarios. For instance, back projection (BP) [28] and rang migration (RM) [29] have been implemented in high-resolution MIMO imaging. However, the performances of these methods highly depend on the aperture of virtual array. To further improve imaging ability, the compressed sensing theory, which has been strictly proved [30], are widely used in radar imaging, including MIMO imaging and ISAR imaging [31]. Therefore, numerous research has been focused on modifying and improving CS-based methods for MIMO and MIMO-SAR 3D imaging. In [32], the author proposed a CS-based method combined with a 2-step motion compensation process for MIMO-SAR imaging. In [33], a L1-norm minimization method is used to do spatial filtering for deep neural network (DNN) for MIMO-SAR 3D imaging in low signal-to-noise ratio (SNR) scenarios. Z. Zhou et al. [34] made a comparison of BP and sparse Bayesian recovery via iterative minimum (SBRIM) and concluded that generally SBRIM has better imaging ability but suffers from high computational complexity. A tensor compressed sensing algorithm is proposed for 3D imaging with Cross-MIMO array [35]. An efficient MIMO-SAR imaging algorithm based on alternating direction method of multipliers (ADMM) is conducted and validated by experimental results [36]. Most existing researches focus on computational efficiency, motion compensation and multipath effect. However, few researchers pay attention to super-resolution MIMO-SAR 3D imaging while considering the unique nature of TDM-MIMO.

This article is mainly contributed to utilizing the property of TDM-MIMO to achieve unambiguous MIMO-SAR 3D imaging on a fast-moving platform. As discussed in paragraphs above, if a TDM-MIMO radar employs more transmitters to achieve high resolution, it will also have a lower PRF which makes MIMO-SAR imaging ambiguous on a fast-moving platform, such as a car or an aircraft. In this article, a novel CS-based super-resolution algorithm is proposed to achieve unambiguous MIMO-SAR 3D imaging by making good use of the phase shift caused by motion of platform. Through both simulation and real data experiments, we show that our algorithm perform well on a fast moving platform, while other methods fail.

The rest of this article is organized as follows. Details of the theory and the algorithm are shown in Section 2. Simulation and experimental results are listed and discussed in Section 3. Finally, a conclusion is presented in Section 4.

2. Signal Model

In this section, the signal model of TDM-MIMO-SAR 3D imaging is deduced and simplified. Moreover, our deductions are verified by simulation data.

The geometry of TDM-MIMO-SAR 3D imaging is shown in Figure 1. Figure 1(a) demonstrates the imaging scenario of a typical airborne sensing application with a linear MIMO array of 4 transmitters and 8 receivers. For simplicity, an equivalent geometry is demonstrated in Figure 1(b) where only one single point target is shown with its coordinates. In the following discussion, we will follow the notation and geometry depicted in Figure 1(b), where the x axis, y axis and z axis stand for the cross-track direction, the track direction and the height direction. The linear TDM-MIMO array is distributed along the z axis, and thus the height resolution is provided by the MIMO array. The platform is moving along the y axis with a velocity $vm/s$. The MIMO array composed of $N_t$ transmitters and $N_r$ receivers can be equivalent to $N_t\times N_r$ virtual elements. The inter-element space of transmitters and receivers are $d_t=M_r{d}_r$ and $d_r=\lambda /2$. $\lambda$ stands for the wavelength of the transmitted signal. For a typical time-division working mode, it can be assumed that all receivers are activated during the entire TDM period, while transmitters take turns to transmit radar waves. A single TDM period contains $N_t$ pulses and each pulse lasts for $T$ seconds. So, a full TDM period lasts for $N_{t}T$ seconds. To complete a MIMO-SAR imaging working circle, it takes $N_{sar}$ TDM periods to transmit $N_{sar}{N}_t$ pulses.

It is well known that virtual array approximation is useful and sufficient for most MIMO applications. In a TDM-MIMO-SAR working mode shown in Figure 1(b), the equivalent virtual array is distributed in a two-dimensional aperture shown in Figure 2. As discussed before, each TDM period costs for $N_{t}T$ seconds and the entire working circle which contains $L$ TDM period lasts for $L{N}_{t}T$ seconds. With a platform velocity of $v$ m/s along y axis, the aperture length along y axis is $(L-1)N_{t}Tv$ m, while the aperture length along z axis is $\left(N_t{N}_r-1\right)d_r$ m. However, it is shown in Figure 2 that due to the property of time-division working mode, virtual elements of the same TDM period are not exactly distributed on the same line. For instance, if we assume the coordinates of the first sub-aperture of $N_r$ virtual elements in the first TDM period are $\{\left(\mathrm{0,0},0\right),\left(\mathrm{0,0},d_r\right),\dots \dots ,(\mathrm{0,0},(N_r-1)d_r\left)\right\}$, the coordinates of the second sub-aperture in the first TDM period are $\{\left(0,Tv,N_r{d}_r\right),\left(0,Tv,\left(N_r+1\right)d_r\right),\dots \dots ,(0,Tv,(2N_r-1)d_r\left)\right\}$. There are several motion compensation methods developed to achieve 3D imaging for TDM-MIMO radar. However, these methods cannot exploit the full potential of TDM-MIMO-SAR. In the following sections, we will demonstrate a novel algorithm which can exploit this property of TDM-MIMO to achieve unambiguous 3D imaging.

For simplicity, we assume that the coordinates of the platform at time $t=0$ are $(\mathrm{0,0},0)$, thus the $N_t$ transmitters are at $\{\left(\mathrm{0,0},0\right),\left(\mathrm{0,0},d_t\right),\dots ,(\mathrm{0,0},\left(N_t-1\right)d_t\left)\right\}$ and the $N_r$ receivers are at $\{\left(\mathrm{0,0},0\right),\left(\mathrm{0,0},d_r\right),\dots ,(\mathrm{0,0},\left(N_r-1\right)d_r\left)\right\}$. With a velocity $v$ towards y axis, the platform moves to $(0,l{N}_{t}Tv,0)$ at the $lth$ TDM period. Suppose that there are $P$ targets with corresponding coordinates $\{\left(x_1,y_1,z_1\right),\dots ,\left(x_p,y_p,z_p\right),\dots ,(x_P,y_P,z_P\left)\right\}$. For a virtual element synthesized by the $n_{t}th$ transmitter and the $n_{r}th$ receiver, the baseband receiver signal $y\left(t,n_t,n_r,p,l\right)$ reflected by the $pth$ target is written as

$\begin{array}{c}y\left(t,n_t,n_r,p,l\right)=A_{p}s\left(t-{\displaystyle \frac{R_{n_{t}pl}+R_{n_{r}pl}}{c}}\right)\left(1\right)\end{array}$

$\begin{array}{c}{R}_{n_{t}pl}=\sqrt{x_p^2+{\left(y_p-\left(\left(l-1\right)N_t+n_t-1\right)Tv\right)}^2+{\left(z_p-\left(n_t-1\right)d_t\right)}^2}\left(2\right)\\ R_{n_{r}pl}=\sqrt{x_p^2+{\left(y_p-\left(\left(l-1\right)N_t+n_t-1\right)Tv\right)}^2+{\left(z_p-\left(n_r-1\right)d_r\right)}^2}\left(3\right)\end{array}$

Where $t$ denotes the fast time of one pulse, $A_p$ is the reflection coefficient of the$pth$ targets and $s\left(t\right)$ is the baseband signal transmitted by antennas. For simplicity, we assume that a chirp signal which is widely used in radar imaging is employed in the MIMO-SAR application.

$\begin{array}{c}s\left(t\right)=exp(j2\pi f_{c}t+j\pi K{t}^2)\left(4\right)\end{array}$

$\begin{array}{c}y\left(t,n_t,n_r,p,l\right)=A_p\exp \left(j2\pi f_c\left(t-{\displaystyle \frac{R_{n_{t}pl}+R_{n_{r}pl}}{c}}\right)+j\pi K{\left(t-{\displaystyle \frac{R_{n_{t}pl}+R_{n_{r}pl}}{c}}\right)}^2\right)\left(5\right)\end{array}$

Where $K$ represents the slope of the frequency change and $f_c$ is the carrier frequency. For a de-chirping process, the received signal is multiplied by the complex conjugate of the transmitted signal.

$\begin{array}{c}y\left(t,n_t,n_r,p,l\right)s^{*}\left(t\right)=A_p\exp \left(-j2\pi f_c{\tau }_{n_t{n}_{r}pl}-j2\pi K{\tau }_{n_t{n}_{r}pl}t+j\pi K{\tau }_{n_t{n}_{r}pl}^2\right)\left(6\right)\end{array}$

$\begin{array}{c}{\tau }_{n_t{n}_{r}pl}={\displaystyle \frac{R_{n_{t}pl}+R_{n_{r}pl}}{c}}\left(7\right)\end{array}$

A Fast Fourier Transform (FFT) can be conducted to complete the de-chirping process.

$\begin{array}{c}x\left(k,n_t,n_r,p,l\right)=FT\left(y\left(t,n_t,n_r,p,l\right)s^{*}\left(t\right)\right)\left(8\right)\\ x\left(k,n_t,n_r,p,l\right)=A_p\exp \left(-j2\pi f_c{\tau }_{n_t{n}_{r}pl}+j\pi K{\tau }_{n_t{n}_{r}pl}^2\right)FT\left(\exp \left(-j2\pi K{\tau }_{n_t{n}_{r}pl}t\right)\right)\left(9\right)\end{array}$

Where $k$ represents the $kth$ range bin. The equation above can be further simplified and rewritten as

$\begin{array}{c}x\left(k,n_t,n_r,p,l\right)\approx {\widehat{A}}_p\delta \left(k-r_p\right)\exp \left(-j2\pi f_c{\tau }_{n_t{n}_{r}pl}+j\pi K{\tau }_{n_t{n}_{r}pl}^2\right)\left(10\right)\end{array}$

$\delta \left(n\right)=\left\{\begin{array}{c}1;n=0\\ 0;n\ne 0\end{array}\right.$

$n_p$ means that the $pth$ target is on the $r_{p}th$ range bin. Sum up reflected signals of all targets and gaussian noise. After the de-chirping process, the four-dimensional sampled data can be expressed as

$\begin{array}{c}x\left(k,n_t,n_r,l\right)=\sum _{p=1}^{P}x\left(k,n_t,n_r,p,l\right)+n\left(k,n_t,n_r,l\right)\left(11\right)\end{array}$

Where $n\left(k,n_t,n_r,l\right)$ represents an additive white gaussian noise. According to the virtual array approximation, equation (7) can be approximated to a more concise form. The coordinates of virtual element of the transmitter and receiver discussed above, which is the$\left(n_t-1\right)N_r+n_{r}th$ virtual element, should be $\left(0,\left(\left(l-1\right)N_t+n_t-1\right)Tv,\left(n_t-1\right)d_t/2+\left(n_r-1\right)d_r/2\right)$. For simplicity, this virtual element is called the $mth$ virtual element and the distance between the $mth$ virtual element and the $pth$ target is denoted as $R_{mpl}$.

$\begin{array}{c}{R}_{n_{t}pl}\approx \sqrt{\begin{array}{c}{x}_p^2+{\left(y_p-\left(\left(l-1\right)N_t+n_t-1\right)Tv\right)}^2+{\left(z_p-\left(\left(n_t-1\right)d_t/2+\left(n_r-1\right)d_r/2\right)\right)}^2\\ +{\left(z_p-\left(n_t-1\right)d_t\right)}^2-{\left(z_p-\left(\left(n_t-1\right)d_t/2+\left(n_r-1\right)d_r/2\right)\right)}^2\end{array}}\\ =\sqrt{\begin{array}{c}{R}_{mpl}^2+{\left(z_p-\left(n_t-1\right)d_t\right)}^2-{\left(z_p-\left(\left(n_t-1\right)d_t/2+\left(n_r-1\right)d_r/2\right)\right)}^2\end{array}}\\ =R_{mpl}{\left(1+{\displaystyle \frac{{\left(z_p-\left(n_t-1\right)d_t\right)}^2-{\left(z_p-\left(\left(n_t-1\right)d_t/2+\left(n_r-1\right)d_r/2\right)\right)}^2}{R_{mpl}^2}}\right)}^{{\displaystyle \frac{1}{2}}}\\ \approx R_{mpl}+{\displaystyle \frac{{\left(z_p-\left(n_t-1\right)d_t\right)}^2-{\left(z_p-\left(\left(n_t-1\right)d_t/2+\left(n_r-1\right)d_r/2\right)\right)}^2}{{2R}_{mpl}}}\end{array}\left(12\right)$

Similarly, $R_{n_{r}pl}$ can also be approximated according to the virtual element approximation.

$\begin{array}{c}{R}_{n_{r}pl}\approx R_{mpl}+{\displaystyle \frac{{\left(z_p-\left(n_r-1\right)d_r\right)}^2-{\left(z_p-\left(\left(n_t-1\right)d_t/2+\left(n_r-1\right)d_r/2\right)\right)}^2}{{2R}_{mpl}}}\left(13\right)\end{array}$

By combining the equation (12) and equation (13), we get an approximation of equation (7). The second approximation in equation (14) is done by the far field assumption.

$\begin{array}{c}{\tau }_{n_t{n}_{r}pl}\approx {\displaystyle \frac{2R_{mpl}}{c}}+{\displaystyle \frac{{\left(\left(n_t-1\right)d_t-\left(n_r-1\right)d_r\right)}^2}{2c{R}_{mpl}}}\\ \approx {\displaystyle \frac{2R_{mpl}}{c}}\end{array}\left(14\right)$

With the virtual array approximation, the four-dimensional data can be reshaped to a three-dimensional cube shown in Figure 3. For simplicity, we assume that $L$ is small enough that range cell migration does not happen in one MIMO-SAR working circle and both the synthetic aperture and MIMO aperture are significantly smaller than the rang distance of the target and the radar, such that the far field assumption can be established.

By dividing the $N_s\times N_t{N}_r\times L$ data into $N_s$ separate 2d $N_t{N}_r\times L$ matrix, we can further simplify the 3D imaging problem to a 2D imaging problem in every range bins. Suppose that the $pth$ target is in the $kth$ range bin. The corresponding 2d data matrix can be denoted as

$\begin{array}{c}{X}_{2d}\left(m,l,p\right)={\widehat{A}}_p\exp \left(-j2\pi f_c{\tau }_{mpl}+j\pi K{\tau }_{mpl}^2\right)\left(15\right)\end{array}$

Where $m$ represents the $mth$ virtual element and $l$ represents the $lth$ TDM-period. The coordinates of this virtual element are $(0,y\left(m,l\right),z(m\left)\right)$. By introducing the far-field assumption, the equation (15) can be simplified.

$\begin{array}{c}{\tau }_{mpl}={\displaystyle \frac{2R_{mpl}}{c}}\approx {\displaystyle \frac{2\sqrt{x_p^2+y_p^2+z_p^2}\left(1-\frac{{(y}_{p}y\left(m,l\right)+z_{p}z\left(m\right)}{\sqrt{x_p^2+y_p^2+z_p^2}}\right)}{c}}\left(16\right)\end{array}$

$\begin{array}{c}{X}_{2d}\left(m,l,p\right)\approx {\widehat{A}}_p\exp \left(j4\pi {\displaystyle \frac{\left(\left(l-1\right)N_t+n_t-1\right)Tv{y}_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\exp \left(j2\pi {\displaystyle \frac{\left(\left(n_t-1\right)d_t+\left(n_r-1\right)d_r\right)z_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\left(17\right)\end{array}$

Equation (17) is a result of ignoring constant term and quadratic term while substituting (16) for (15). The echo signal can be expressed as:

$\begin{array}{c}{X}_{2d}\left(m,l\right)=\sum _{p\in \mathsf{\Omega }\left(k\right)}{\widehat{A}}_p\exp \left(j4\pi {\displaystyle \frac{\left(\left(l-1\right)N_t+n_t-1\right)Tv{y}_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}+j2\pi {\displaystyle \frac{\left(\left(n_t-1\right)d_t+\left(n_r-1\right)d_r\right)z_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\left(18\right)\end{array}$

It is shown clearly in equation (18) that the TDM-MIMO-SAR 3D imaging is different from a naïve 2d imaging problem, as the property of time-division working mode introduces additional phase shift. Moreover, the phase shift is proportional to the exact $y_p$ of the target location. By simulation, our deduction can be validated in Figure 4. Figure 4 shows a simulation of a TDM-MIMO-SAR echo signal data after pulse compression. Details about the simulation settings are listed in Table 1.

3. Algorithm

In this section, we will firstly demonstrate a matched filter-based method with phase compensation for TDM-MIMO-SAR 3D imaging. Although this matched filter-based method performs well when the platform velocity is low, simulations show its limits when the platform moves fast due to gating lobes. Then we will show our proposed CS-based super-resolution method and its ability to suppress gating lobes for a TDM-MIMO-SAR on a fast-moving platform. To reduce the computational complexity, a 2D CS method is employed to reduce the computational complexity after regrouping the echo signal data.

3.1. The Matched Filter-based Imaging Scheme

To make a comparison, a conventional matched-filter based imaging scheme is demonstrated in this section. Following the deduction shown in Section 2, a straightforward matched filter-based imaging scheme with phase compensation is shown in Figure 5.

According to equation (18), the phase compensation term $h\left(m\right)$ is expressed as

$\begin{array}{c}h\left(m\right)=\exp \left(j4\pi {\displaystyle \frac{floor\left(\frac{m-1}{N_r}\right)Tv{y}_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right);m=\mathrm{1,2},\dots ,N_t{N}_r\left(19\right)\end{array}$

The MF-based algorithm provides imaging results without super-resolution capability. Moreover, it does not utilize the phase shift caused by the TDM working mode to deal with gating lobes which occur when the platform moves fast. A simulation result shown in Figure 6 will demonstrate how gating lobes affect the imaging process of MIMO-SAR as the velocity of the platform increases. The detailed parameters of this simulation are the same as parameters listed in Table 1. The number of TDM-periods $L$ is 16. The geometry of imaging process is shown in Figure 2.

For real world MIMO-SAR applications, a fast-moving platform, such as a car or an aircraft, is a common situation. Although it is possible to avoid the gating lobes problem by reducing PRI. However, the PRI of TDM-MIMO mainly depends on the number of transmitters and reducing the number of transmitters significantly affects the number of virtual elements, which is vital for MIMO radar ability. For TDM-MIMO enjoys the advantage of low hardware complexity and costs, it is important to develop a new algorithm for TDM-MIMO-SAR applications which can deal with the gating lobes problem so that combining the powerful capability of MIMO-SAR imaging scheme and the low complexity advantage of TDM-MIMO can be feasible on fast-moving platforms.

3.2. Proposed 2D-CS based Imaging Scheme

In this section, we demonstrate our proposed imaging algorithm based on regrouping echo data and 2D-CS super-resolution method. A brief mathematical proof is provided, and numerical simulation is conducted to compare the proposed method and other methods.

3.1.1. Compressed Sensing Theory and Methods

Let $X\in {\mathbb{C}}^N$ denotes a signal of finite length $N$. There exists a basis matrix $\mathsf{\Psi }=\{{\psi }_1,{\psi }_2,\dots ,{\psi }_N\}$ such that $X=\mathsf{\Psi }\mathsf{\Theta }$, where $\mathsf{\Theta }=\left\{{\mathsf{\theta }}_{\mathrm{i}}\right\}$ is a K-sparse vector with length $N$. For a K-sparse vector, $\mathsf{\Theta }$ can be approximated by its $K$ largest coefficients. Let $\mathsf{\Phi }\in {\mathbb{C}}^{M\times N}$ denotes the measurement matrix. The measured signal $Y\in {\mathbb{C}}^M$ can be expressed as

$\begin{array}{c}Y=\Phi X=\Phi \Psi \Theta \left(19\right)\end{array}$

Recovering $X$ from $Y$ is an ill-posed problem for $M\ll N$. However, according to the CS theory, the recovery is possible as long as $\mathsf{\Phi }\mathsf{\Psi }$ has the Restricted Isometry Property (RIP). The detailed expression of RIP property is given as follows:

$\begin{array}{c}1-\epsilon K\le {\displaystyle \frac{{‖\mathsf{\Phi }\mathsf{\Psi }\mathsf{\Theta }‖}_2}{{‖\mathsf{\Theta }‖}_2}}\le 1+\epsilon K\left(20\right)\end{array}$

Although it is difficult to validate whether a not a measurement matrix satisfies the RIP directly, RIP is usually satisfied if the rows of $\mathsf{\Phi }$ and columns of $\mathsf{\Psi }$ are incoherent. With the help of RIP, the recovery problem is equal to solving the optimization problem which can be expressed as

$\begin{array}{c}\min {‖\mathsf{\Theta }‖}_{0}s.t.Y=\Phi \Psi \Theta \left(21\right)\end{array}$

However, the discontinuousness of$L0$ norm leads to the difficulty of solving optimization problem (21). Therefore, there are many CS methods solving the optimization problem by applying relaxation of $L0$ norm. For instance, the smoothed $l0$ (SL0) method and $l1$ method. In this article, we mainly focus on the $l1$ method which can be denoted as

$\begin{array}{c}\min {‖\mathsf{\Theta }‖}_{1}s.t.Y=\Phi \Psi \Theta \left(22\right)\end{array}$

Moreover, the constrained optimization problem (22) can be converted to the unconstrained convex optimization problem shown in (23).

$\begin{array}{c}\widehat{\mathsf{\Theta }}={\underset{\mathsf{\Theta }}{\arg \min }‖Y-A\mathsf{\Theta }‖}_2^2+\lambda {‖\mathsf{\Theta }‖}_1\left(23\right)\end{array}$

Where $A\in {\mathbb{C}}^{M\times N},A=\mathsf{\Phi }\mathsf{\Psi }$. An efficient approach to solve the problem (23) is the fast iterative shrinkage-thresholding algorithm (FISTA). By applying proximal gradient and Nesterov’s acceleration, the FISTA enjoys the advantage of both efficiency and accuracy. The detailed procedure of the FISTA is listed in Table 3.

By applying the two-dimensional model of $L1$ norm optimization problem, a 2D-CS model can be expressed as

$\begin{array}{c}\widehat{\mathsf{\Theta }}={\underset{\mathsf{\Theta }}{\arg \min }‖Y-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}‖}_F^2+\lambda {‖\mathsf{\Theta }‖}_1\left(24\right)\end{array}$

Where $Y\in {\mathbb{C}}^{M\times M}$, $A\in {\mathbb{C}}^{M\times N}$, $B\in {\mathbb{C}}^{M\times N}$, $\mathsf{\Theta }\in {\mathbb{C}}^{N\times N}$. Therefore, after applying the proximal gradient with the optimization problem (24), we can get

$\begin{array}{c}L\left(\mathsf{\Theta }\right)={‖Y-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}‖}_F^2+\lambda {‖\mathsf{\Theta }‖}_1=F\left(\mathsf{\Theta }\right)+G\left(\mathsf{\Theta }\right)\left(25\right)\end{array}$

Where $F\left(\mathsf{\Theta }\right)={‖Y-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}‖}_F^2$ and $G\left(\mathsf{\Theta }\right)=\lambda {‖\mathsf{\Theta }‖}_1$. It is obvious that $G\left(\mathsf{\Theta }\right)$ is convex and $F\left(\mathsf{\Theta }\right)$ follows the Lipschitz continuity.

$\begin{array}{c}F\left({\mathsf{\Theta }}_2\right)\le F\left({\mathsf{\Theta }}_1\right)+<\nabla F\left({\mathsf{\Theta }}_1\right),\left({\mathsf{\Theta }}_2-{\mathsf{\Theta }}_1\right)>+{\displaystyle \frac{L}{2}}{‖{\mathsf{\Theta }}_2-{\mathsf{\Theta }}_1‖}_F^2\left(26\right)\end{array}$

Then the approximation of $L\left(\mathsf{\Theta }\right)$ can be expressed as

$\begin{array}{c}Q\left({\mathsf{\Theta }}_2,{\mathsf{\Theta }}_1\right)=F\left({\mathsf{\Theta }}_1\right)+<\nabla F\left({\mathsf{\Theta }}_1\right),\left({\mathsf{\Theta }}_2-{\mathsf{\Theta }}_1\right)>+{\displaystyle \frac{L}{2}}{‖{\mathsf{\Theta }}_2-{\mathsf{\Theta }}_1‖}_F^2+G\left({\mathsf{\Theta }}_2\right)\left(27\right)\end{array}$

The proximal gradient of $L\left(\mathsf{\Theta }\right)$ is the gradient of $Q\left({\mathsf{\Theta }}_2,{\mathsf{\Theta }}_1\right)$.

$\begin{array}{c}{\nabla }_{{\mathsf{\Theta }}_2}Q\left({\mathsf{\Theta }}_2,{\mathsf{\Theta }}_1\right)=\nabla F\left({\mathsf{\Theta }}_1\right)+L\left({\mathsf{\Theta }}_2-{\mathsf{\Theta }}_1\right)+\nabla G\left({\mathsf{\Theta }}_2\right)\left(28\right)\end{array}$

The optimal solution of the proximal gradient can be achieved when the equation (28) is equal to $0$. Therefore, the optimal solution $P_L\left({\mathsf{\Theta }}_1\right)$ is

$\begin{array}{c}{P}_L\left({\mathsf{\Theta }}_1\right)=soft\left({\mathsf{\Theta }}_1-{\displaystyle \frac{1}{L}}\nabla F\left({\mathsf{\Theta }}_1\right),{\displaystyle \frac{\lambda }{L}}\right)\left(29\right)\end{array}$

$\begin{array}{c}soft\left(x,t\right)=\left\{\begin{array}{c}0,\left|x\right|\le t\\ sign\left(x\right)\left(\left|x\right|-t\right),\left|x\right|>t\end{array}\right.\left(30\right)\end{array}$

Now, we substitute (25) into (29).

$\begin{array}{c}\nabla F\left(\mathsf{\Theta }\right)=\nabla {‖Y-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}‖}_F^2\\ =\nabla tr\left({\left(\mathrm{Y}-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}\right)}^{\mathrm{H}}\left(\mathrm{Y}-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}\right)\right)=2{\mathrm{A}}^{\mathrm{H}}\left(\mathrm{Y}-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}\right)B^H\end{array}\left(31\right)$

With the help of Nesterov’s acceleration, the detailed procedure of 2D-FISTA is demonstrated in Table 4.

3.1.2. The Proposed Algorithm

Following the expression in the equation (18), a novel algorithm which can utilize the property of TDM-MIMO to achieve unambiguous MIMO-SAR 3D imaging. In this section, first we analyze the motion-induced phase shift along the virtual array, then the full procedure of the proposed algorithm is demonstrated. Finally, the performance of the proposed algorithm is verified by simulation data.

After slow time compressing of the 2d data shown the equation (18) by MF methods or super-resolution methods, the 2d data $Z_{2d}(m,y)$ is denoted as

$\begin{array}{c}Z\left(m,y\right)=FT\left(X_{2d}\left(m,l\right)\right)\left(31\right)\end{array}$For a point target $(x_p,y_p,z_p)$, the corresponding signal after ignoring constant term is

$\begin{array}{c}{Z}_{y_p}\left(n_t,n_r\right)=\exp \left(j4\pi {\displaystyle \frac{\left(n_t-1\right)Tv{y}_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\exp \left(j2\pi {\displaystyle \frac{\left(\left(n_t-1\right)d_t+\left(n_r-1\right)d_r\right)z_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\left(32\right)\end{array}$

With proper phase compensation shown in (19), the compensated 2d data is

$\begin{array}{c}{Z^{\text{'}}}_{y_p}\left(n_t,n_r\right)=h\left(n_t,n_r,y_p\right)Z_{y_p}\left(n_t,n_r\right)=\exp \left(j2\pi {\displaystyle \frac{\left(\left(n_t-1\right)d_t+\left(n_r-1\right)d_r\right)z_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\left(33\right)\end{array}$

However, as discussed previously, for a TDM-MIMO radar, it needs more transmitters to obtain more degree of freedom for MIMO imaging, while the increasing number of transmitters directly leads to the decreasing of PRF. And for MIMO-SAR applications, low PRF usually leads to ambiguous imaging results. Assume that for the point target $(x_p,y_p,z_p)$, it has a gating lobe located at $y_{pg}$. After the phase compensation in (33), the corresponding signal is

$\begin{array}{c}{Z^{\text{'}}}_{y_{pg}}\left(n_t,n_r\right)=h\left(n_t,n_r,y_{pg}\right)Z_{y_p}\left(n_t,n_r\right)\\ =\exp \left(j4\pi {\displaystyle \frac{\left(n_t-1\right)Tv\left(y_p-y_{pg}\right)}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}+j2\pi {\displaystyle \frac{\left(\left(n_t-1\right)d_t+\left(n_r-1\right)d_r\right)z_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\\ =\exp \left(j2\pi {\displaystyle \frac{\left(n_t-1\right)d_t{z}_{pg}+\left(n_r-1\right)d_r{z}_p}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right)\\ \end{array}\left(34\right)$

Where $z_{pg}={\displaystyle \frac{2Tv\left(y_p-y_{pg}\right)}{d_t}}+z_p$ .The main idea of our proposed method is to distinguish the real target and its gating lobe, so that an unambiguous MIMO-SAR 3D image can be achieved. Now, we apply the 2D-FISTA demonstrated in Table 4 to the phase-compensated 2D signal ${Z^{\text{'}}}_{y_p}\left(n_t,n_r\right)$.

$\begin{array}{c}\widehat{\mathsf{\Theta }}={\underset{\mathsf{\Theta }}{\arg \min }‖Z^{\text{'}}-A\mathsf{\Theta }{\mathrm{B}}^{\mathrm{H}}‖}_F^2+\lambda {‖\mathsf{\Theta }‖}_1\left(35\right)\end{array}$

Where $Z^{\text{'}}\in {\mathbb{C}}^{N_t\times N_r}$, $Z^{\text{'}}={\left\{Z^{\text{'}\left(n_t,n_r\right)}\right\}}_{ij}$ , $\mathsf{\Theta }\in {\mathbb{C}}^{N_g\times N_g}$, $A\in {\mathbb{C}}^{N_t\times N_g}$, $A=\left\{a_i\right\}$, $B\in {\mathbb{C}}^{N_t\times N_g}$, $B=\left\{b_i\right\}$. Measurement matrix $A$ and $B$ consists of steering vectors $a_i$ and $b_i$, $a_i\in {\mathbb{C}}^{N_t}$, $b_i\in {\mathbb{C}}^{N_r}$. $N_g$ is the number of atoms in the dictionary, which contains $\{z_1,z_2,\dots ,z_{N_g}\}$.

$\begin{array}{c}{a}_i\left(n\right)=\exp \left(j2\pi {\displaystyle \frac{\left(n-1\right)d_t{z}_i}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right);n=\mathrm{0,1},2,\dots ,n_t\left(36\right)\end{array}$

$\begin{array}{c}{b}_i\left(n\right)=\exp \left(-j2\pi {\displaystyle \frac{\left(n-1\right)d_r{z}_i}{\lambda \sqrt{x_p^2+y_p^2+z_p^2}}}\right);n=\mathrm{0,1},2,\dots ,n_r\left(37\right)\end{array}$

After 2D-FISTA process, a 2D result $\widehat{\mathsf{\Theta }}\in {\mathbb{C}}^{N_g\times N_g}$ is obtained.

For a real target and its corresponding signal ${Z^{\text{'}}}_{y_p}\left(n_t,n_r\right)$, its imaging result appear on the main diagonal of $\widehat{\mathsf{\Theta }}$. On the other hand, for a gating lobe, its imaging result is located off the main diagonal, as it is shown in the equation (38). Based on the equation (33) and (34), we can obtain

$\begin{array}{c}\left\{\begin{array}{c}{Z^{\text{'}}}_{y_p}=a_{z_p}{b}_{z_p}^H\\ {Z^{\text{'}}}_{y_{pg}}=a_{z_{pg}}{b}_{z_p}^H\end{array}\right.\left(38\right)\end{array}$

By extracting the main diagonal of $\widehat{\mathsf{\Theta }}$, we manage to achieve unambiguous MIMO-SAR 3D imaging. A simulation result shown in Figure 7 demonstrates how the real target differs from gating lobes in the 2D result $\widehat{\mathsf{\Theta }}$ obtained by equation (35). The simulation settings are the same as in Figure 6.

As is clearly shown in Figure 7, the proposed scheme manages to eliminate gating lobes which distribute along y axis. After the phase compensation and 2D-CS process, the gating lobes should be filtered by extracting the main diagonal. Therefore, although the imaging result along y axis is ambiguous, the gating lobes would not appear in the final 3D imaging result as gating lobes are filtered and eliminated.

The full imaging scheme is demonstrated in Figure 8.

3.3. Simulation Results

In this section, the imaging effects of the proposed method in the TDM-MIMO-SAR 3D imaging are analyzed by simulation results. A comparison of the ordinary MF-based method and the proposed method is made to demonstrate the difference between ambiguous imaging results and unambiguous ones. The major parameters of this simulation test are listed in Table 6. The geometry of the MIMO-SAR imaging is shown in Figure 2.

As can be seen from simulation results shown in Figure 9, the MF-based method suffers from gating lobes when the platform velocity is high. The increasing of the platform velocity severely degrades the imaging result of the MF-based method, meanwhile the proposed method manages to achieve unambiguous imaging.

As is demonstrated in Table 5, the process of the proposed method relies on information about the velocity of the platform. In real world applications, the measured velocity may not be the exact velocity. Therefore, it is vital to analyze the effect of velocity error. In this section, the effect of velocity error is analyzed by simulation results under the same setup as shown in Table 6. The simulation of the proposed method under velocity error is demonstrated in Figure 10.

It is indicated in simulation results that although the velocity error indeed degrades the imaging performance of the proposed method, the proposed method still maintains unambiguous imaging capability with a velocity error less than 10% of the exact velocity.

4. Experimental and Analysis

In this section, a real world TDM-MIMO-SAR 3D imaging experiment is conducted to verify the proposed method. An AWR2243 cascaded MIMO radar is used to collect measured data. Developed by Texas Instruments (TI), AWR2243 is a multichannel array radar with short-range MIMO imaging capability. The key parameters of AWR2243 in MIMO mode are shown in Table 7.

Among all 12 transmitters of AWR2243, 9 of these transmitters are placed in a uniform linear array. Detailed structure of AWR2243 is shown in Figure 11.

As shown in Figure 12(a)(b), the MIMO radar system is installed on a moving vehicle which can reach 3 m/s. The platform moves along a straight road as indicated in Figure 12(c). Three corner reflectors are set by the road as shown in Figure 12(d). It is clearly shown in Figure 12(e) that gating lobes occur in the experiment circumstance. Figure 12(f) demonstrates the proposed method’s ability for unambiguous MIMO-SAR imaging.

To further demonstrate the performance of the proposed method under more realistic circumstances, we installed the AWR2243 radar on an actual car which runs at a speed of 20km/h as shown in Figure 13. The main target is a stationary car by the road which the platform is moving along. Figure 13(c, d) demonstrates the 2D projection of the MF-based method and the proposed method and it is clearly shown that the proposed method effectively eliminates the gating lobes and achieves unambiguous MIMO-SAR imaging.

5. Conclusions

In this paper, an unambiguous MIMO-SAR 3D imaging algorithm is proposed for applications on a fast-moving platform with TDM-MIMO. For TDM-MIMO radar, an increase of the number of transmitters significantly decreases the PRF of each transmitter, which leads to ambiguous MIMO-SAR imaging for a fast-moving platform, making it necessary to develop an unambiguous TDM-MIMO-SAR imaging algorithm. By analyzing the detailed virtual array geometry and the corresponding signal model of TDM-MIMO-SAR, we manage to utilize the phase shift induced by the time division property to solve the ambiguity of MIMO-SAR imaging. A 2D-CS method based on FISTA is used in the proposed methods for super-resolution imaging. Compared with the MF-based method, the proposed method successfully eliminates gating lobes and thus achieve unambiguous MIMO-SAR imaging. The performance of the proposed method is verified by both simulation data and real-world experiments. In order to fulfill the requirements of real world applications, the effect of velocity error is analyzed by simulation results. The main contribution of this paper is to propose a super-resolution method which can utilize the time-division property of TDM-MIMO radar to solve ambiguity of MIMO-SAR imaging.