Ground Clutter Mitigation for Slow-Time MIMO Radar Using Independent Component Analysis

Abstract

The detection of low, slow and small (LSS) targets, such as small drones, is a developing area of research in radar, wherein the presence of ground clutter can be quite challenging. LSS targets, because of their unusual flying mode, can be easily shadowed by ground clutter, leading to poor radar detection performance. In this study, we investigated the feasibility and performance of a ground clutter mitigation method combining slow-time multiple-input multiple-output (st-MIMO) waveforms and independent component analysis (ICA) in a ground-based MIMO radar focusing on LSS target detection. The modeling of ground clutter under the framework of st-MIMO was first defined. Combining the spatial and temporal steering vector of st-MIMO, a universal signal model including the target, ground clutter, and noise was established. The compliance of the signal model for conducting ICA to separate the target was analyzed. Based on this, a st-MIMO-ICA processing scheme was proposed to mitigate ground clutter. The effectiveness of the proposed method was verified with simulation and experimental data collected from an S-band st-MIMO radar system with a desirable target output signal-to-clutter-plus-noise ratio (SCNR). This work can shed light on the use of ground clutter mitigation techniques for MIMO radar to tackle LSS targets.

1. Introduction

Small drone detection using radar has attracted enormous attention in recent years [1,2,3]. With the rapid growth of the consumer drone market, unmanned aerial vehicles (UAVs) have become a significant threat to civil aviation, anti-terrorism, and private security. As a powerful sensor that can operate regardless the time and weather, radar plays an important role in tackling these low, slow and small (LSS) targets. Many systems [4,5,6,7,8,9,10,11,12,13] and techniques [14,15,16,17,18,19] focusing on LSS target detection in the field of radar have been researched and developed.

The early LSS target surveillance radar systems were modified from the maritime radar systems such as the Merlin^TM Radar System from DeTect Inc. [20] and the initial product of Robin Radar Systems Inc. [21]. They consist of two maritime radar antennas that rotate along the azimuth and elevation to achieve quasi-3D detection. With the development of phased array radar, more LSS target detection systems have turned to utilize antenna array and digital beamforming technology to obtain better target detection performance [4,5]. In recent years, some relatively new radar concepts, such as multiple-input multiple-output (MIMO) radar [6,7,8], multistatic radar [9,10], and ubiquitous radar (which is also named holographic or staring radar) [11,12,13], have also been introduced to the field of LSS target detection.

There are two main challenges in detecting LSS targets. The first challenge is the poor target signal-to-clutter-plus-noise ratio (SCNR). Due to the low flying altitude of the small drones, the radar beam must have a rather small grazing angle. This significantly raises the ground clutter energy in the received signal. Considering the small radar cross-section (RCS) of LSS targets, the target SCNR is significantly reduced. Secondly, the slow flying, or hovering, velocity of the target renders the detection of LSS in slow-moving ground clutter via conventional techniques such as moving target indication (MTI) or moving target detector (MTD) ineffective. This is because ground clutter mainly consists of buildings, trees, cars, etc. The velocity spectrum spread of this clutter ranges from zero to tens of meters per second, which overshadows the speed of LSS targets. Thus, it is worth investigating clutter mitigation methods for LSS target detection.

Independent component analysis (ICA) [22,23,24] has aroused worldwide research interest in the field of signal processing since the 1990s. It has extensive applications in many fields such as communication, radar, image processing, acoustic processing, biomedical signal processing, and even financial data analysis [25]. Based on a MIMO system, the purpose of ICA is to simultaneously separate independent non-Gaussian components from observed multi-channel signals. This process can be used to help solve the ground clutter mitigation problem.

MIMO radar, as a typical multi-channel system, utilizes omnidirectional antennas and an orthogonal waveform in the transmitting stage, and then it forms a synthesized MIMO beam in the receiving stage. MIMO radar leverages waveform diversity to further increase the scale of the virtual array and spatial diversity [25,26,27]. A slow-time MIMO (st-MIMO) waveform [28,29,30] expands a conventional radar waveform by phase-coding the pulses of different channels to achieve orthogonal transmission and MIMO demodulation after pulse-Doppler (PD) processing. Although the unambiguous Doppler speed is divided by the number of Doppler sub-bands [31], st-MIMO is still feasible for slow-moving target detection with acceptable orthogonality and bandwidth efficiency.

Over the last few years, various techniques have been developed to tackle LSS targets from the radar detection stage to the parameter estimation stage and the target classification stage. Regarding methods used to enhance LSS target detection performance at the range-Doppler (r-D) level, one study [14] used the stationary point concentration technique to reduce the noise floor caused by transmitter leakage and increase the signal-to-noise ratio (SNR) of a UAV. Another study [15] adopted an iterative adaptive approach to enhance the Doppler resolution in an r-D map, which led to improved target detectability. For the parameter estimation stage, the authors of [16] proposed a long-time coherent integration method for maneuvering LSS targets to improve target estimation accuracy. Regarding the target classification stage, many works [2,17,18] focused on the micro-Doppler signature (m-DS) of UAVs. The numbers of m-DSs have been thoroughly analyzed among different types of UAVs and utilized to conduct target classification via multiple classifiers or neural networks.

There have been many studies in the field of clutter mitigation involving the aforementioned st-MIMO or ICA methods separately. Nevertheless, the combination of these two techniques while focusing on ground clutter mitigation in detecting LSS targets has been limited. For clutter mitigation involving st-MIMO, the authors of [28] developed spatial-time adaptive processing (STAP) in st-MIMO radar to mitigate multipath clutter. The spatial-time structure of st-MIMO was established, and a data covariance model, including target, direct path clutter, multipath clutter, jamming, and white noise, was defined. In [32], the beamspace st-MIMO, which can form virtual transmit nulls in directions that would result in multipath clutter returns in the main lobe of a radar system, was investigated. In [33], ground clutter was removed via principal components analysis (PCA) to enhance micro-Doppler feature extraction. For clutter mitigation involving the ICA or BSS methods, many works have focused on the multi-mode clutter suppression of over-the-horizon (OTH) radar. In [34,35], the spread-Doppler clutter caused by multi-mode propagation in OTH radar was suppressed via the second-order blind identification (SOBI) method. In addition, some works applied the BSS method in the field of main lobe jamming suppression [36].

In this study, we investigated the feasibility and performance of a ground clutter mitigation method combining the st-MIMO and ICA techniques in a ground-based MIMO radar focusing on LSS target detection. The main contributions of this work are as follows. Firstly, we propose a way to modulate ground clutter under the framework of st-MIMO. The ground clutter covariance was derived based on a Gaussian-shaped power spectrum. Combining the spatial and temporal steering vector of st-MIMO, a universal signal modelling including the target, ground clutter, and noise is provided. Secondly, we propose a st-MIMO-ICA processing scheme to separate the target signal from the received data including ground clutter plus noise. Mathematical proof that the signal model fits the framework of the ICA problem is provided. The multi-targets can be simultaneously separated via the proposed method. Finally, we validated our proposed st-MIMO-ICA method and evaluated its performance with both simulations and experiments using an S-band st-MIMO radar system developed in our previous work [6]. Comparisons with the PCA and adaptive techniques were performed, and the st-MIMO-ICA method showed the highest target output SCNR of the tested approaches.

The remainder of this paper is organized as follows. Section 2 describes the signal modeling of the st-MIMO radar and the ground clutter modulation in st-MIMO. Section 3 describes the signal modeling for conducting ICA and proposes the st-MIMO-ICA method. Section 4 presents the simulation and experimental results for the performance of the st-MIMO-ICA method. Finally, Section 5 concludes the paper and outlines possible future work.

2. Signal Modeling of St-MIMO Radar

In this section, the signal modeling of st-MIMO is described. Furthermore, the modulation of ground clutter under the framework of st-MIMO is derived.

2.1. St-MIMO Waveform Modeling and Processing

Consider a co-located 1D linear antenna array with $M$ transmitting elements and $N$ receiving elements that are omnidirectional. The element distance to the reference antenna of the $m^{th},m=0,\cdots ,M-1$ element in the transmitting array is $d_m$. Likewise, the element distance to the reference antenna of the $n^{th},n=0,\cdots ,N-1$ element in the receiving array is $d_n$. Note that this definition of the array is general regardless of whether the array is uniform or not. Figure 1 shows the geometry of the signal modeling setup. The operating frequency is $f_0$ and the operating wavelength is ${\lambda }_0$. There are $K$ pulses in one coherent processing interval (CPI) and the pulse repetition frequency (PRF) $f_r=1/T_r$, where $T_r$ refers to pulse repetition interval (PRI).

The slow-time MIMO approach split the whole Doppler PRF into $M$ orthogonal Doppler sub-bands with a bandwidth of $\Delta f_{sub}=f_r/M$ via slow-time phase coding. Here, we set the number of the Doppler sub-bands as the same as the number of transmitting antennas to simplify the modeling. Furthermore, the Doppler sub-bands can be redundant, leading to some empty Doppler sub-bands that can be utilized to enlarge the velocity measurement range. The baseband pulse waveform $u_p(t)$ of each transmitting element is an identical linear frequency modulated (LFM) signal but with varying starting phases $\phi (m,k)$, which is a function of the transmitting element index $m$ and the pulse index (slow-time) $k$. Let ${\rho }_t$ represent the constant transmit amplitude in each antenna without beamforming transmission. Then, the transmitting waveform of the $m^{th}$ element is:

$${\mathrm{s}}_m(t)={\rho }_t{\displaystyle \sum _{k=0}^{K-1}{u}_p\left(t-k{T}_r\right)e^{j2\pi \left(f_{0}t+\phi (m,k)\right)}}$$

(1)

Herein, let $\phi (m,k)$ have a linear form:

$$\phi (m,k)={\alpha }_{m}k{T}_r$$

(2)

In this way, the Doppler domain is divided into $M$ identical sub-bands. The instantaneous Doppler frequency $f_d^m$ of the $m^{th}$ transmitting antenna is the derivative of the slow time variable $k{T}_r$:

$$f_d^m=\frac{\partial \phi (m,k)}{\partial k{T}_r}={\alpha }_m$$

(3)

Let the Doppler frequency of each transmitting antenna be evenly distributed in the Doppler domain and consider that the sign of the Doppler, ${\alpha }_m$, has the form of:

$${\alpha }_m=\frac{m}{M}\cdot f_r$$

(4)

One of the advantages of the slow-time MIMO approach is the good hardware compatibility to traditional phased array radar systems. Orthogonal transmission is realized in the Doppler domain via time-varying starting phases, while the carrier frequency remains $f_0$ for each channel. This makes it easy to implement to existing array radar systems. Note that the time-varying starting phases between antennas also presumably steer the beam as a function of the slow-time pulse $k$, which means that the main lobe direction of the transmitted beampattern sweeps the angle domain with a period of $M$ pulses [30]. The ratio $f_r/\Delta f_{sub}$, i.e., the number of the Doppler sub-bands, can be modified to adjust the sweep rate of the main lobe.

Next, consider a far-field moving target with target speed $v_t$ and corresponding Doppler shift $f_t{=2\mathrm{v}}_t/{\lambda }_0$. The target is located at an angle ${\varphi }_t$ related to the array boresight. The amplitude of the echo is ${\rho }_r$. Thus, the backscattered signal of the $n^{th}$ receiving element via the $m^{th}$ transmitting element is shown in (5).

$${\mathrm{s}}_{mn}(t)={\rho }_r{\displaystyle \sum _{k=0}^{K-1}{u}_p\left(t-{\tau }_{mn}-k{T}_r\right)e^{j2\pi \left(f_0+f_t\right)\left(t-{\tau }_{mn}\right)}}{e}^{j2\pi {\alpha }_{m}k{T}_r}$$

(5)

The round-trip delay ${\tau }_{mn}$ in (5) takes the form of (6), which consists of the time delay from the $m^{th}$ transmitting element TX_m to the target and the time delay from the target to the $n^{th}$ receiving element TX_n [37]:

$$\begin{array}{ll}{\tau }_{mn}& ={\tau }_{tm}+{\tau }_{rn}\\ & =\left(\frac{R_t}{c}-\frac{d_m\sin {\varphi }_t}{c}\right)+\left(\frac{R_t}{c}-\frac{d_n\sin {\varphi }_t}{c}\right)\\ & =\frac{2R_t}{c}-\frac{d_m\sin {\varphi }_t}{c}-\frac{d_n\sin {\varphi }_t}{c}\end{array}$$

(6)

where $R_t$ is the target range based on the first element of the array. Under the conditions of a narrow band and slow target velocity, the following assumptions can be made:

$${\tau }_{mn}\approx {\tau }_c=\frac{2R_t}{c}$$

(7)

and $e^{j2\pi f_t{\tau }_{mn}}\approx 1$.

Then, ${\mathrm{s}}_{mn}(t)$ can be approximated as:

$${\mathrm{s}}_{mn}(t)\approx {\rho }_r{\displaystyle \sum _{k=0}^{K-1}{u}_p\left(t-{\tau }_c-k{T}_r\right)e^{j2\pi \left(f_0+f_t\right)t}{e}^{-j2\pi f_0{\tau }_{mn}}}{e}^{j2\pi {\alpha }_{m}k{T}_r}$$

(8)

The output of the $n^{th}$ receiving element is the sum of the $M$ transmit waveforms and can be expressed as:

$$\begin{array}{ll}{\mathrm{s}}_n(t)& ={\displaystyle \sum _{m=0}^{M-1}{\mathrm{s}}_{mn}(t)}\\ & ={\rho }_r{\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{k=0}^{K-1}{u}_p\left(t-{\tau }_c-k{T}_r\right)}}{e}^{j2\pi f_{0}t}{e}^{j2\pi f_{t}t}{e}^{-j2\pi f_0{\tau }_{mn}}{e}^{j2\pi {\alpha }_{m}k{T}_r}\end{array}$$

(9)

After down-converting by multiplying the echo with $e^{-j2\pi f_{0}t}$, the baseband receiving signal can be expressed as:

$$\begin{array}{ll}\mathrm{s}{\prime }_n^{}(t)& ={\mathrm{s}}_n^{}(t)\cdot e^{-j2\pi f_{0}t}\\ & ={\rho }_r{\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{k=0}^{K-1}{u}_p\left(t-{\tau }_c-k{T}_r\right)}}\hspace{0.17em}{e}^{j2\pi f_{t}t}{e}^{-j2\pi f_0\frac{2R_t}{c}}{e}^{j2\pi \frac{d_m}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi {\alpha }_{m}k{T}_r}\end{array}$$

(10)

Next, conduct matched filtering on $\mathrm{s}{\prime }_n^{}(t)$ using a baseband matched filter $h(t)=u_p^{\ast }(-t)$, where the superscript * represents the conjugate operation. The output signal after matched filtering can be derived as:

$$\begin{array}{ll}{X}_n^{}(t)& ={\displaystyle {\int }_{\infty }^{\infty }\mathrm{s}{\prime }_n^{}(x)h(t-x)dx}\\ & ={\rho }_r{e}^{-j2\pi f_0\frac{2R_t}{c}}{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{k=0}^{K-1}{e}^{j2\pi \frac{d_m}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi {\alpha }_{m}k{T}_r}}}\cdots \\ & {\displaystyle \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\int }_{\infty }^{\infty }{u}_p\left(x-{\tau }_c-k{T}_r\right)e^{j2\pi f_{t}x}{u}_p^{\ast }(-(t-x))dx}\end{array}$$

(11)

After conducting variable substitution and assuming $e^{j2\pi f_t{\tau }_c}\approx 1$, (11) can be derived as:

$$\begin{array}{ll}{X}_n^{}(t)& \approx {\rho }_r{e}^{-j2\pi f_0\frac{2R_t}{c}}{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{k=0}^{K-1}{e}^{j2\pi \frac{d_m}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi {\alpha }_{m}k{T}_r}}}\cdots \\ & {\displaystyle \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\int }_{\infty }^{\infty }{u}_p\left(\beta \right)e^{j2\pi f_t\beta }{e}^{j2\pi f_{t}k{T}_r}{u}_p^{\ast }(\beta +{\tau }_c+k{T}_r-t)d\beta }\\ & ={\rho }_r{e}^{-j2\pi f_0\frac{2R_t}{c}}{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{k=0}^{K-1}{e}^{j2\pi \frac{d_m}{{\lambda }_0}\sin {\varphi }_t}}}\cdots \\ & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{e}^{j2\pi \left({\alpha }_m+f_t\right)k{T}_r}\chi \left(t-{\tau }_c-k{T}_r,f_t\right)\end{array}$$

(12)

where

$$\chi \left(\tau ,f_t\right)={\displaystyle {\int }_{\infty }^{\infty }{u}_p\left(\beta \right)u_p^{\ast }(\beta -\tau )e^{j2\pi f_t\beta }d\beta }$$

(13)

is the ambiguity function of $u_p\left(t\right)$ with the time lag being $\tau$ and the Doppler frequency being $f_t$.

Next, we focus on the fast time of $t_k={\tau }_c+k{T}_r$ in each pulse, which corresponds to the time lag of the target. Then, the ambiguity function in (12) is $\chi \left(0,f_t\right)$. Assume $\chi \left(0,f_t\right)\approx 1$; because of the high Doppler tolerance of the LFM signal, the response of each $t_k$ in (12) can be expressed as:

$$\begin{array}{ll}{X}_{nk}^{}& =X_n^{}(t_k)\\ & ={\xi }_t{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{\displaystyle \sum _{m\prime =0}^{M-1}\hspace{0.17em}{e}^{j2\pi \frac{d_{m\prime }}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi \left({\alpha }_{m\prime }+f_t\right)k{T}_r}}\end{array}$$

(14)

where ${\xi }_t={\rho }_r{e}^{-j2\pi f_0\frac{2R_t}{c}}$.

In order to separate the response of the $m^{th}$ transmitting element and achieve MIMO demodulation, we first shift the central Doppler frequency of the $m^{th}$ Doppler sub-band to zero-Doppler by multiplying $e^{-j2\pi {\alpha }_{m}k{T}_r}$ to (14) which yields:

$$\begin{array}{ll}{X}_{nk,m}^{}& =X_{nk}^{}\cdot e^{-j2\pi {\alpha }_{m}k{T}_r}\\ & ={\xi }_t{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi f_{t}k{T}_r}{e}^{j2\pi \frac{d_m}{{\lambda }_0}\sin {\varphi }_t}\\ & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\xi }_t{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi f_{t}k{T}_r}{\displaystyle \sum _{m=0,m\prime \ne m}^{M-1}{e}^{j2\pi \frac{d_{m\prime }}{{\lambda }_0}\sin {\varphi }_t}}{e}^{j2\pi \left({\alpha }_{m\prime }-{\alpha }_m\right)k{T}_r}\end{array}$$

(15)

For certain $n$ and $m$, $X_{nk,m}^{}$ is only related to the slow-time index $k$. Next, conduct a discrete Fourier transform (DFT) of $X_{nk,m}^{}$ to obtain the Doppler spectrum, which yields:

$$X_{n,k\prime ,m}^{}={\displaystyle \sum _{k=0}^{K-1}{X}_{nk,m}^{}{e}^{-j\frac{2\pi k}{K}k\prime }}$$

(16)

where $k\prime$ represents the index in the Doppler domain. Then, apply a low-pass Doppler filter to remove all the other $M-1$ Doppler sub-bands in $X_{nk,m}^{}$, i.e., the second term in (15) [38]. Thus, the low-pass Doppler filter $H_{LP}$ has a pass-band from $-f_r/2M$ to $f_r/2M$. Conduct an inverse discrete Fourier transform (IDFT) to obtain the temporal output after Doppler filtering, which yields:

$$\begin{array}{ll}X{\prime }_{n,k,m}^{}& ={\displaystyle \sum _{k\prime =0}^{K-1}{X}_{n,k\prime ,m}^{}{H}_{LP}(k\prime )e^{j\frac{2\pi k\prime }{K}k}}\\ & \approx {\xi }_t{e}^{j2\pi \frac{d_n}{{\lambda }_0}\sin {\varphi }_t}{e}^{j2\pi f_{t}k{T}_r}{e}^{j2\pi \frac{d_m}{{\lambda }_0}\sin {\varphi }_t}\end{array}$$

(17)

Note that $X{\prime }_{n,k,m}^{}$ can be expressed as the Kronecker product of separable vectors:

$$X_t\left({\varphi }_t,f_t\right)=a_r\left({\varphi }_t\right)\otimes b\left(f_t\right)\otimes a_t\left({\varphi }_t\right)$$

(18)

where $a_r\left({\varphi }_t\right)$ and $a_t\left({\varphi }_t\right)$ are the receiving and transmitting spatial steering vectors, respectively, of the target as the function of ${\varphi }_t$ and $b\left(f_t\right)$ is the temporal steering vector of the target as the function of $f_t$:

$$a_r\left({\varphi }_t\right)={\left[1,e^{j\frac{2\pi }{{\lambda }_0}{d}_1\sin {\varphi }_t},\cdots ,e^{j\frac{2\pi }{{\lambda }_0}{d}_n\sin {\varphi }_t},\cdots ,e^{j\frac{2\pi }{{\lambda }_0}{d}_{N-1}\sin {\varphi }_t}\right]}^T$$

(19)

$$b\left(f_t\right)={\left[1,e^{j2\pi f_t{T}_r},\cdots ,e^{j2\pi f_{t}k{T}_r},\cdots ,e^{j2\pi f_t\left(K-1\right)T_r}\right]}^T$$

(20)

$$a_t\left({\varphi }_t\right)={\left[1,e^{j\frac{2\pi }{{\lambda }_0}{d}_1\sin {\varphi }_t},\cdots ,e^{j\frac{2\pi }{{\lambda }_0}{d}_m\sin {\varphi }_t},\cdots ,e^{j\frac{2\pi }{{\lambda }_0}{d}_{M-1}\sin {\varphi }_t}\right]}^T$$

(21)

Because the Doppler bandwidth of $b\left(f_t\right)$ is narrowed to $f_r/M$ after the low-pass filtering rather than the entire $f_r$, M-times decimation can be conducted on $b\left(f_t\right)$ to reduce the dimension of the vector, which yields:

$$b_{deci}\left(f_t\right)={\left[1,e^{j2\pi \cdot M\cdot f_t{T}_r},\cdots ,e^{j2\pi \cdot \left(K/M-1\right)M\cdot f_t{T}_r}\right]}^T$$

(22)

At this point, all the responses of the $M$ transmitting elements are extracted from the $N$ receiving elements. The final data cube of the target signal for one snapshot has the dimensions of $N\times \left(K/M\right)\times M$ and can be expressed as:

$${\zeta }_s=a_r\left({\varphi }_t\right)\otimes b_{deci}\left(f_t\right)\otimes a_t\left({\varphi }_t\right)$$

(23)

Furthermore, the radar return from a target-plus-clutter-plus-noise environment for one snapshot can then be expressed as:

$${\zeta }_x={\zeta }_s+{\zeta }_c+{\zeta }_n$$

(24)

where ${\zeta }_c$ and ${\zeta }_n$ are the clutter vector and noise vector with dimensions of $N\times \left(K/M\right)\times M$, respectively.

2.2. Ground Clutter Modeling under St-MIMO Framework

In this section, the modeling of the clutter matrix ${\zeta }_c$ based on the Doppler distributed clutter (DDC) model [39] is introduced. The fundamental principle of DDC modelling is to first decide the power spectrum $S_{ddc}(f)$ of the ground clutter. Then, the autocorrelation $r_{ddc}(\tau )$ and the covariance matrix $R_{ddc}$ of the clutter can be derived using an inverse Fourier transform (IFT).

The power spectrum $S_{ddc}(f)$ of the ground clutter caused by the internal motion of the clutter itself is commonly assumed to be Gaussian-shaped [19] with the form of:

$$S_{ddc}(f)=\frac{P_c}{\sqrt{2\pi {\sigma }_c^2}}\exp \left[-\frac{{\left(f-f_d\right)}^2}{2{\sigma }_c^2}\right]$$

(25)

where $P_c$ is the power of the clutter, ${\sigma }_c^{}$ is the standard deviation, and $f_d$ is the central frequency of the clutter spectrum. We chose $f_d=0$ for the ground clutter, and the 3 dB velocity spectrum width ${\sigma }_v^{}$ has the following relationship with the standard deviation ${\sigma }_c^{}$:

$${\sigma }_c^{}=\frac{2{\sigma }_v^{}}{{\lambda }_0}$$

(26)

By conducting an IFT of Equation (25), the autocorrelation function $r_{ddc}(\tau )$ of the clutter can be derived as:

$$r_{ddc}(\tau )=P_c\exp \left(j2\pi f_d\tau \right)\exp \left(-2{\pi }^2{\sigma }_c^2{\tau }^2\right)$$

(27)

Accordingly, the covariance matrix of the clutter can be represented as:

$$\begin{array}{r}{R}_{ddc}\left(k,l\right)=P_c^2{e}^{j2\pi f_d{T}_r\left(k-l\right)}{e}^{-2{\pi }^2{\sigma }_c^2{T}_r{}^2{\left(k-l\right)}^2}\\ k=1,\cdots ,K/M,\hspace{0.17em}l=1,\cdots ,K/M\end{array}$$

(28)

The temporal sequence of the clutter $x_c$ can be generated using the following equation:

$$x_c=R_{ddc}^{\frac{1}{2}}{x}_n$$

(29)

where $x_n$ is the $K/M$ white Gaussian noise vector with a variance of 1.

For the ground clutter at the target’s range bin, assume that there are $N_C$ clutter directions of arrival (DOA) in the received data. Then, the clutter snapshot ${\zeta }_c$ can be expressed as:

$${\zeta }_c={\displaystyle \sum _{i=1}^{N_C}{a}_r\left({\varphi }_{ci}\right)\otimes x_{ci}\otimes a_t\left({\varphi }_{ci}\right)}$$

(30)

where $i=1,\cdots ,N_C$ represents the $i^{th}$ clutter patch with the DOA of ${\varphi }_{ci}$ and $x_{ci}$ represents the temporal vector of the $i^{th}$ clutter patch. At this point, ground clutter modeling in st-MIMO is achieved.

3. St-MIMO-ICA Processing

In this section, the proposed st-MIMO-ICA processing method is illustrated. Firstly, the feasibility of using the ICA method to mitigate ground clutter with a st-MIMO system is demonstrated. Then, the proposed st-MIMO-ICA processing method is illustrated based on a complex fixed-point algorithm.

3.1. Signal Modeling under ICA Compliance

The previous sections illustrated the signal and clutter model under the st-MIMO framework. Now, rewrite the sequence of the terms in (23) as:

$${\zeta }_s=b_{deci}\left(f_t\right)\otimes a\left({\varphi }_t\right)$$

(31)

where $a\left({\varphi }_t\right)=a_r\left({\varphi }_t\right)\otimes a_t\left({\varphi }_t\right)$ is the $N\times M$ two-way spatial vector of the st-MIMO waveform. From this point, the virtual array of the MIMO radar is established to achieve narrower beamwidth and realize more degrees of freedom along the spatial dimension compared with a conventional phased array radar. Generally, in multiple target situations, Equation (31) can be expressed as:

$${\zeta }_s={\displaystyle \sum _{j=1}^{N_T}{b}_{deci}\left(f_{tj}\right)\otimes a\left({\varphi }_{tj}\right)}$$

(32)

where $j=1,\cdots ,N_T$ represents the $j^{th}$ target with an associated Doppler frequency of $f_{tj}$ and target DOAs of ${\varphi }_{tj}$.

Regarding ground clutter, the same reasoning can be used and the clutter snapshot ${\zeta }_c$ can be expressed as:

$${\zeta }_c={\displaystyle \sum _{i=1}^{N_C}{x}_{ci}\otimes a\left({\varphi }_{ci}\right)}$$

(33)

Further exploiting the structures of Equations (32) and (33) leads to a universal signal model, including the target and clutter in a matrix form:

$$X=A{S}_{Dop}$$

(34)

In Equation (34), the columns of matrix $A$ represent the spatial vector of the target together with the clutter:

$$A=\left[a_{t1},a_{t2},\cdots a_{t{N}_T},a_{c1},a_{c2},\cdots a_{c{N}_C}\right]$$

(35)

and the rows of matrix $S_{Dop}$ represent the frequency vectors in the Doppler domain of the target and the clutter:

$$S_{Dop}={\left[d_{t1}^H,\cdots ,d_{tj}^H,\cdots d_{t{N}_T}^H,d_{c1}^H,\cdots ,d_{ci}^H,\cdots d_{c{N}_C}^H\right]}^T$$

(36)

where $d_{tj}^{}$ and $d_{ci}^{}$ are the DFT of the temporal vectors $b_{deci}\left(f_{tj}\right)$ and $x_{ci}$, respectively.

Equation (34) shows that the signal model of the st-MIMO waveform in a ground clutter environment enables the use of the ICA method to extract targets from clutter, that is: (1) multi-channel observations of the mixed signal are obtained, (2) sources are linearly mixed and statistically independent from each other, and (3) sources have non-Gaussian distributions [40]. By exploiting the structure of Equation (34), the received signal $X$ can be regarded as multi-channel inputs for ICA. Meanwhile, the number of the input channels is enlarged from $N$ to $M\times N$ via st-MIMO demodulation. The spatial matrix $A$ can be regarded as the mixing matrix that linearly combines all the sources in the frequency matrix $S_{Dop}$, of which the vectors are non-Gaussian in the Doppler domain. Thus, using the ICA method for ground clutter mitigation with a st-MIMO system is feasible.

3.2. St-MIMO-ICA Processing

ICA can be illustrated in a general framework that consists two major parts, the separation criteria and the optimization method [22]. In the application of array radar signal processing, a number of ICA methods have been proposed based on variable criteria such as maximum likelihood (ML), information-maximization (Infomax), and the maximization of non-Gaussianity (MN). In the branch of MN, kurtosis and negentropy are two commonly used criteria used in some popular ICA methods such as JADE [24] and FastICA [40], respectively.

The famous FastICA method utilizes negentropy as the cost function and the fixed-point algorithm for optimization. Furthermore, the FastICA method has been derived from complex domains. However, the complex FastICA (c-FastICA) method does not perform well with noncircular sources [41]. In this study, we adopted the noncircular FastICA (nc-FastICA) method to separate the target signal from clutter and noise. By adding second-order information in a fixed-point update, the nc-FastICA method can provide an improved separation performance with noncircular sources.

Based on (34), the linear signal mixture model of the MIMO array with $Q=M\times N$ sensors (channels) and $N_s=N_T+N_C$ sources can be expressed as:

$$\begin{array}{ll}{\zeta }_x& ={\zeta }_s+{\zeta }_c+{\zeta }_n\\ & =A{S}_{Dop}+{\zeta }_n\end{array}$$

(37)

where ${\zeta }_x$ can be considered to be the observation matrix, which contains the aforementioned signal matrix $A{S}_{Dop}$ and the Gaussian white noise with zero-mean ${\zeta }_n$. Then, the nc-FastICA method can be conducted via the following steps.

• Whitening:
  The covariance matrix of the observed signal can be expressed as:

  $$R_Y=E\left\{{\zeta }_x{\zeta }_x{}^H\right\}$$

  (38)
  By applying eigenvalue decomposition to the covariance matrix, one can obtain:

  $$R_Y=U\Lambda U^H$$

  (39)
  Then, the whitening matrix can be denoted as:

  $$V={\Lambda }^{-\frac{1}{2}}{U}^H$$

  (40)

  and the whitened matrix $Z=V{\zeta }_x$ can be achieved.
• Formulating the optimization problem
  A cost function based on maximizing the negentropy can be express as:

  $$J(w)=E\left\{G\left({\left|w^{H}Z\right|}^2\right)\right\}$$

  (41)

  where $G:ℝ\to ℝ$ is a smooth even function and $w\in {ℂ}^Q$ is the $Q\times Q$ weight matrix used to de-mix the observed signal and obtain the estimated source matrix $e=w^{H}Z$. Then, the optimization problem is formulated as:

  $$w_{opt}=\arg \underset{{\Vert w\Vert }^2=1}{\max }E\left\{G\left({\left|w^{H}Z\right|}^2\right)\right\}$$

  (42)
  Here, we chose $G(u)=u^2/2$ as the function motivated by kurtosis.
• Fixed-point updating process
  The fixed-point algorithm is utilized to update the weight matrix $w$ in each iteration and can be expressed as:

  $$\begin{array}{ll}{w}_{n+1}& =-E\left\{g\left({\left|e\right|}^2\right)e^{*}Z\right\}+E\left\{g^{\prime }\left({\left|e\right|}^2\right){\left|e\right|}^2+g\left({\left|e\right|}^2\right)\right\}{w}_n\\ & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+E\left\{e{e}^T\right\}E\left\{g^{\prime }\left({\left|e\right|}^2\right)e^{*2}\right\}{w}_n^{*}\end{array}$$

  (43)

  where $g(u)=dG(u)/du$ and $g^{\prime }(u)=dg(u)/du$. The third term of Equation (43) includes the second-order information in terms of the pseudo-covariance matrix $E\left\{e{e}^T\right\}$, which is non-zero if the sources are noncircular. This modification ensures that the nc-FastICA method has an improved separation performance with noncircular sources.
• Obtain the weight matrix and estimated source matrix
  The estimation of the observation matrix ${\zeta }_x$ can be expressed as:

  $$e=w^{H}Z=w^{H}V{\zeta }_x$$

  (44)

  where $e=\left[e_1,e_2,\dots ,e_Q\right]$ is the estimated source matrix with $Q$ source vectors and $w={\left[w_1,w_2,\dots ,w_Q\right]}^T$ is the weight matrix with $Q$ corresponding weight vectors.
  It is worth noting that the number of sources separated by the ICA method is no more than the channel number $Q$. When the number of sources in the mixed signal is less than the number of the channels, those redundant channels will contain noise signals [22]. In this case, the PCA method can be utilized to decrease the dimensions of the data for the following ICA processing. However, this part is beyond the scope of this paper.

4. Experimental Results

For this section, simulation and field experiment were conducted based on the MIMO radar system developed by the Beijing Institute of Technology (BIT) [6,7], which is shown in Figure 2. The ground-based radar operates at the S-band and has a co-located MIMO architecture. The antenna array is organized along the elevation axis with six elements transmitting st-MIMO waveforms. Each antenna has a wide beam of 90° in elevation, where the MIMO is formed, and a narrow beam of 3° in azimuth. For the azimuth, mechanical scanning is utilized to achieve full airspace coverage. The system is mainly focused on LSS target detection and tracking, such as that of small drones or birds. The st-MIMO-ICA method for ground clutter mitigation is first demonstrated with simulation results using parameters corresponding to a real radar system. Then, the field experiment results are presented.

4.1. Simulation Results

The simulations were conducted in MATLAB. They were employed to verify the performance of the proposed method. To reflect the performance of a real radar system used in the field, consider a co-located MIMO antenna array with $M=N=6$ linear displaced elements along the elevation axis. Each element transmits a st-MIMO waveform. The simulated radar parameters are shown in

Table 1. Radar parameters.

Parameter	Value
Operating frequency	3 GHz
Bandwidth	40 MHz
Pulse repetition interval	56 us
Number of pulses in one CPI	900
Number of Doppler sub-band	6
Maximum detect range	8.4 km
Velocity measurement range	±72 m/s

.

Furthermore, the antenna array was chosen to be a non-uniform sparse array to make it cost-effective with a larger aperture and fewer elements. The element location was optimized via a genetic algorithm (GA) to suppress the grating lobes of the antenna beampattern caused by the sparse arrangement. A detailed illustration and analysis of this non-uniform and sparse arrangement can be found in [6]. The array parameters are shown in

Table 2. Array parameters.

Parameter	Value
Number of antenna	6
Number of channels	36
Array aperture	0.5 m
Virtual array aperture	1 m
Element position	[0 0.1821 0.2681 0.3555 0.4206 0.5] m
Peak sidelobe level	−20 dB
3 dB beamwidth	6.2°

, and the beam pattern of the array after GA optimization is shown in Figure 3.

The simulation experiment setup is shown in Figure 4. The radar was ground-based, and the MIMO was achieved in elevation. The ground clutter came from the horizontal direction and Target 1 to Target $N_T$ had the target DOA 1 to DOA $N_T$. A typical LSS target could be a small drone such as a DJI Phantom series drone. This type of target has a relatively slow speed and low flying altitude, which means the radar detection performance is affected by strong ground clutter. To simplify the simulation without losing generality, assume that there are three targets in the same range bin with different SNRs, speed, and DOAs. Regarding the ground clutter, we used the clutter model in Section 1 with four closely distributed DOAs to simulate a real ground environment containing trees and buildings of different heights. We set proper target SNRs to keep the target SCNRs of all the targets at −20 dB. Detailed parameters of the target and clutter setup are given in

Table 3. Target and clutter parameters.

Type	Parameter	Value
Target 1	Range	500 m
	Speed	2 m/s
	Elevation	4°
	SNR	40 dB
	SCNR	−20 dB
Target 2	Range	500 m
	Speed	−4 m/s
	Elevation	20°
	SNR	30 dB
	SCNR	−20 dB
Target 3	Range	500 m
	Speed	6 m/s
	Elevation	10°
	SNR	20 dB
	SCNR	−20 dB
Ground clutter	Spectrum center	0 m/s
	3 dB Spectrum width	±2 m/s
	DOAs	[0° 0.1° 0.2° 0.3°]
	CNR	60 dB

.

The results of signal modeling are shown in Figure 5. The range-velocity map of the received signal in antenna 1 after pulse-compression and PD processing is displayed in Figure 5a, from where the waveform diversity of st-MIMO has formed and the echoes from the six transmitting antenna, TX1 to TX6, have been shifted into different Doppler sub-bands. Moreover, the velocity spectrum at targets’ range bin is displayed in Figure 5b. where the black boxes indicate six Doppler sub-bands, DS1 to DS6. Together with the received signals in the other five antennas, all the $Q=M\times N=36$ transmitting-receiving paths, named as Channel 1 to Channel 36, can be established. The velocity spectrum of Channel 1 after MIMO demodulation is displayed in blue line in Figure 5c, and the red line represent the targets plus noise in order to indicate the targets information clearly. The CNR is 60 dB after pulse-compression and PD processing, and different target SNRs are set to keep the target SCNR as −20 dB.

The st-MIMO-ICA output of the mixed signal is shown in Figure 6. Some representative results among all the 36 output channels are displayed. Excluding the uncertainty of the sequence of the output signals [22], the three targets are separated from the mixed targets plus clutter and noise signal and located in Channel 32, 33, and 35 respectively. Channel 34 contains the combination of the ground clutter. Because the clutter DOAs are closely distributed, the method fails to separate those clutter components. Lastly, the remaining channels are noise signals. Figure 7 shows the detailed information of the three separated targets after st-MIMO-ICA processing. All the three targets are precisely separated with correct velocity and desirable output SCNR of 29.8 dB, 23.4 dB, and 33.2 dB, respectively. We compare the results of st-MIMO-ICA with two other approaches which are also frequently utilized to mitigate the ground clutter, that is PCA and adaptive digital beamforming (ADBF). Among the ADBF approaches, the sample matrix inversion (SMI) technique is chosen. The results of these two approaches are also displayed in Figure 7, and it can be seen that the ICA method enjoys the highest target output SCNR. The exact target SCNR values of the three approaches are listed in

Table 4. Comparison of simulation results.

Results (dB)		Target 1	Target 2	Target 3
PCA	SCNR	2.78	5.68	15.17
	Null level	−28.10	−34.82	−53.20
SMI	SCNR	11.79	13.24	23.82
	Null level	−34.27	−44.87	−59.90
ICA	SCNR	29.77	23.42	33.20
	Null level	−52.97	−49.97	−75.46

.

Blind beamforming could be achieved using the weight matrix generated by the st-MIMO-ICA method. Figure 8 displays the beampatterns corresponding to the three targets, with $w_{32}$, $w_{33}$, and $w_{35}$ as the weight vectors. All the beampatterns could form a deep null in the clutter DOA, which effectively mitigated the strong ground clutter. Meanwhile, when the target DOA was close to the clutter DOA, the blind beamforming also suffered from distortion in the main beam and the loss of the target SCNR, which is shown in Figure 8a. Figure 8 also displays the equivalent beampatterns of PCA that were generated from the eigenvectors corresponding to the three targets and the adaptive beampatterns of SMI with the pre-known steering vectors pointing at the three targets. The st-MIMO-ICA method showed the deepest null level at the clutter’s DOA compared with the other two approaches, which illustrates the outplayed performance of the proposed method in ground clutter mitigation. The detailed null levels of the three approaches are also listed in Table 4.

To investigate the performance boundary of the proposed st-MIMO-ICA method, a simulation of the output target SCNR under each target DOA was conducted, which is shown in Figure 9. We still focused on the three aforementioned targets: Target 1, Target 2, and Target 3. We let the DOA of the targets change from −10° to 10° while other parameters remained the same. To simplify the simulation, we also let the clutter DOA be 0°. The number of Monte Carlo simulations was chosen to be 100 for each target DOA. The output target SCNR well-corresponded to the input target SNR plus the coherent processing gain of the MIMO array, which was 15.5 dB ($Q=36$). However, when the target’s DOA was close to the clutter’s DOA, the output target SCNR was decreased due to the increase in the cross-correlation of the two components in the spatial domain. The black dashed line in the zoomed-out figure of Figure 9 indicates the sufficient detection threshold of 13 dB, and the corresponding target DOA was 0.4°, which was 6.5% of the MIMO beamwidth.

4.2. Field Experimental Results

Field experiments were conducted in a complex urban environment to verify the performance of the proposed method. The radar system is deployed in an industrial city in southeast China where factory buildings, trees, moving cars constitute the sources of the ground clutter. The satellite image of the radar location and surroundings courtesy of Google Maps is shown in Figure 10a, and the radar field of view is shown in Figure 10b. The parameters of the radar system were given in Table 1 and Table 2. The target information is listed in

Table 5. Target parameters in field test.

Type	Parameter	Value
Target 4	Range	3041 m
	Speed	2 m/s
	Height	200 m
	Elevation	3.8°
Target 5	Range	5585 m
	Speed	4 m/s
	Height	200 m
	Elevation	2°

. A small drone of DJI Phantom IV is used as a real flying target in this section. The diagonal size of the small drone is 350mm and can be treated as a point target compared to the range resolution of the radar system. There are two range bins chosen as the field experiment location, which are 3041 m and 5585 m. The small drone has a flying altitude of 200 m at both locations, resulting in different target elevations of 3.8° and 2°, respectively.

In the first experiment of 3041 m, the small drone has the speed of 2 m/s. The range-velocity map in Channel 1 after pulse-compression and PD processing is shown in Figure 11, where the data label indicates the target location, speed, and amplitude. The target is submerged in the spectrum of strong ground clutter and cannot be detected. The result of the st-MIMO-ICA processing of Target 4 is shown in Figure 12. The velocity spectrum at Target 4’s range bin before st-MIMO-ICA processing is shown in Figure 12a. The black dashed line indicates the target location in the mixed signal. The target is shadowed by strong ground clutter and cannot be detected via conventional method such as CFAR. The velocity spectrum after st-MIMO-ICA processing is shown in Figure 12b, where the target can be effectively extracted from the strong ground clutter with the highest output target SCNR of 33.9 dB. The other two approaches, PCA and SMI, are also conducted as a comparison with output target SCNR of 14.5 dB and 17.8dB, respectively. The output of the PCA approach also has a false alarm located at −1 m/s.

In the second experiment of 5585 m, the small drone has the speed of 4 m/s and can be seen in the range-velocity map together with some competing ground clutters which belong to the buildings and moving cars. The range-velocity map in Channel 1 after pulse-compression and PD processing is shown in Figure 13, where the data label indicates the target location, speed, and amplitude. The result of the st-MIMO-ICA processing of Target 5 is shown in Figure 14. The velocity spectrum at Target 5’s range bin before st-MIMO-ICA processing is shown in Figure 14a. The target had a SCNR of -10dB compared to the ground clutter located at 0 m/s and -6.5 m/s. The velocity spectrum after st-MIMO-ICA processing is shown in Figure 14b, where the target can be effectively extracted from the competing clutters with output target SCNR of 37.8 dB, and the other two remarkable ground clutters have been successfully suppressed. The results of PCA and SMI are also displayed with output target SCNR of 26.6 dB and 15.1dB, respectively. The output of the PCA approach also has a false alarm in the clutter location of −6.5 m/s.

5. Discussion

A ground clutter mitigation method for slow-time MIMO radar using independent component analysis was researched in this study. Firstly, a ground clutter model under a st-MIMO scheme was provided. The clutter covariance was derived based on a Gaussian-shaped power spectrum, and a universal signal modelling including the target, ground clutter, and noise was provided. Secondly, the compliance for conducting ICA was discussed, and the st-MIMO-ICA processing scheme was proposed. Lastly, the performance of the proposed method was verified with simulation and field experiments via an S-band MIMO radar system.

The simulation results indicated the feasibility of the proposed st-MIMO-ICA method, as the st-MIMO-ICA output showed the highest target SCNR compared with two other conventional clutter mitigation approaches of PCA and SMI. The three test targets with an input SCNR of −20 dB were precisely separated with the correct velocity and desirable output SCNRs of 29.8 dB, 23.4 dB, and 33.2 dB. The st-MIMO-ICA method also showed the deepest null level at the clutter’s DOA compared with the other two approaches. The performance boundary of the method as validated with a separable DOA difference of 6.5% of the MIMO beamwidth. The field experimental results further proved the effectiveness of the proposed method in LSS target detection with strong ground clutter compared with conventional methods.

It is worth noting that the PCA approach only utilizes the eigenvectors of the signal covariance matrix after conducting eigenvalue decomposition. As is noted in Section 1, PCA is often utilized as a preliminary step of ICA to decrease the dimensions of the data. Regarding the SMI technique, we chose to use target-plus-clutter-plus-noise signals in this study as the training data to calculate the adaptive weight vector, which limited the performance of the SMI technique. Generally, it is not feasible to obtain independent and identical distributed clutter samples in real scenes. The ICA method, however, utilizes negentropy as the cost function and the fixed-point algorithm for optimization to achieve better source separation results, which corresponds to desirable ground clutter mitigation performance.

There is still work to be conducted. It is worth further investigating the character of ground clutter, especially regarding the number of the sources of clutter. In this paper, we assumed that the DOAs of ground clutter for a certain snapshot were limited even though actual clutter characteristics are more complex. Another work will focus on the adaption of the detection threshold to further improve the detection performance for LSS targets [42,43].