Robust Multiple-Measurement Sparsity-Aware STAP with Bayesian Variational Autoencoder

Abstract

Due to the shortage of independent and identically distributed (i.i.d.) training samples, space−time adaptive processing (STAP) often suffers remarkable performance degradation in the heterogeneous clutter environment. Sparse recovery (SR) techniques have been introduced into STAP for the benefit of the drastically reduced training requirement, but they are incompletely robust for involving the tricky selection of hyper−parameters or the undesirable point estimation for parameters. Given this issue, we incorporate the Multiple−measurement Complex−valued Variational relevance vector machines (MCV) to model the space−time echoes and provide a Gibbs−sampling−based method to estimate posterior distributions of parameters accurately. However, the Gibbs sampler require quantities of iterations, as unattractive as traditional Bayesian type SR−STAP algorithms when the real−time processing is desired. To address this problem, we further develop the Bayesian Autoencoding MCV for STAP (BAMCV−STAP), which builds the generative model according to MCV and approximates posterior distributions of parameters with an inference network pre−trained off−line, to realize fast reconstruction of measurements. Experimental results on simulated and measured data demonstrate that BAMCV−STAP can achieve suboptimal clutter suppression in terms of the output signal to interference plus noise ratio (SINR) loss, as well as the attractive real−time processing property in terms of the convergence rate and computational loads.

1. Introduction

Airborne phased−array radar is widely used to detect moving targets because of its mobility and freedom from the radius of the earth’s curvature, but it is always plagued by the existence of the ‘nuisance’ clutter coupling in the space−time domain. To ensure the detection of weak targets, space-time adaptive processing (STAP) methods [1,2,3,4] are developed to mitigate the clutter and achieve better performance compared with traditional non−adaptive methods. For classical STAP methods, the performance is substantially determined by the accuracy of the clutter plus noise covariance matrix (CCM). Typically, the CCM of the cell under test (CUT) is estimated by training samples adjacent to the CUT, which are supposed to be independent and identically distributed (i.i.d.) [5]. According to the Reed−Mallett−Brennan (RMB) rule [6], the required training samples should be at least twice the system degree of freedom (DOF) to achieve steady performance. Nevertheless, since the clutter background is usually non−homogeneous and the DOF is usually high in practice, obtaining adequate qualified training samples is challenging, leading to dramatically degraded performance of CCM estimation and clutter suppression. A large number of methods to improve the performance of STAP with finite training samples have been proposed in the past few years, such as the reduced−dimension (RD) methods [7,8,9,10], the reduced−rank methods (RR) [11,12,13], the direct data domain (D3) methods [14] and the knowledge−aided STAP (KA−STAP) methods [15,16,17,18,19]. However, the number of required samples is still large in the RD and RR methods, the DOF of D3 is significantly reduced resulting in severe performance loss, and the KA−STAP methods depend heavily on prior knowledge. By contrast, the sparse recovery STAP (SR−STAP) methods [20,21,22,23], which are not only more tolerant with the number of training samples but also more independent from the prior knowledge, have received significant attention. In SR−STAP, every atom of the dictionary is a space−time steering vector denoting a grid of the space−time plane, which is uniformly over−sampled in both the angle and Doppler. The space−time profile, recovered by solving the formulated optimization problem, will be leveraged to reconstruct the CCM for effective clutter suppression. Since the clutter suppression performance mainly relies on the accuracy of recovery results, many efforts have been put into the problem formulation of sparse recovery. The pioneering work, such as greedy algorithms [24,25,26,27], considers the problem of SR−STAP as an $l_0$−norm minimization problem which is known as NP−hard [28,29]. To tackle the NP−hard problem, appealing approaches transform it into a convex problem [20,30,31,32], but the strong coherence among the atoms of sparse dictionary often leads to the severe performance degradation of sparse reconstruction. Besides, the recovery results obtained by convex optimization are particularly sensitive to the selection of regularization parameters.

Due to the superior ability to deal with the highly coherent dictionary situation and the tricky selection of regularization parameters, Bayesian type sparse recovery (SR) algorithms [22,33,34,35], represented by the sparse Bayesian learning (SBL) methods [36] and its extension for the multiple measurement vector (MMV) case, i.e., M−SBL [37], have received much attention in the field of STAP. However, the potential need for a large number of iterations per testing in these methods is a fatal drawback for STAP where real−time processing is desired. Some methods [33,34,35,38] are devoted to modifying the time−consuming process and accelerating the convergence, but they are at the price of either significant signal loss or expensive computational loads, which also hinder the real−time processing. In addition, all these Bayesian type SR methods are incompletely robust as a result of only providing point estimates for some parameters, such as the parameter controlling the variance of the space−time profile and the parameter denoting the noise power. In summary, the completely robust Bayesian type SR−STAP algorithms that ensure both the stable performance of sparse reconstruction and the feasibility of real−time processing are rarely studied.

To model the space−time echoes with a completely robust Bayesian framework, this paper introduces variational relevance vector machines (VRVM) [39] into STAP due to its attractive advantage of providing distributional estimates for all parameters without any user parameters (known as the parameter−free property) and sufficient ability to model the sparsity. However, the original VRVM is proposed in the real domain for the single measurement vector (SMV) case, which is not suitable for STAP where the complex−valued multiple measurement vector (MMV) case exists. As a result, we first develop the multiple−measurement complex−valued variational relevance vector machines (MCV), which can address the condition where complex−valued multiple measurements occur. To achieve accurate inference, we employ the Gibbs sampler as a general technique to estimate all parameters for MCV. Then, we present the novel STAP algorithm via MCV, named MCV−STAP, to design the adaptive filters. However, a large number of iterations, as well as expensive computational loads, are required in the Gibb sampler, which limits the efficiency of the model for out−of−sample prediction. To render the real−time processing of the MCV−STAP, inspired by the structure of the variational autoencoder (VAE) [40,41,42,43], the powerful models for unsupervised learning in the field of deep learning, we introduce a VAE into our model to provide the variational inference and expand it to a VAE−structured model named Bayesian Autoencoding MCV−STAP (BAMCV−STAP), which is scalable in the training phase and fast in the testing phase. Then, our model can be optimized by minimizing the evidence lower bound (ELBO) via a gradient ascent scheme. We note that due to the gradient ascent scheme without involving inverse operations when learning the parameters, BAMCV−STAP can approximate posterior densities of space−time profiles with low computational loads. Moreover, since the inference network can be pre−trained off−line, only a few iterations are needed before testing. Finally, the space−time profiles recovered via BAMCV−STAP are adopted to estimate the CCM for STAP. To the best of our knowledge, BAMCV−STAP is the first work that attempts a VRVM−with−VAE strategy in the context of STAP, and our method appears to outperform several other sparsity−enhanced methodologies in terms of clutter suppression performance and real−time processing, as will be demonstrated in the experimental results.

The main contributions of this paper are summarized as follows:

• Generalized from the original VRVM to the multiple measurements case existing in the complex domain, a parameter−free probabilistic model called MCV is derived to recover space−time profiles via the Gibbs sampling method for STAP.
• Since all parameters are estimated based on their posterior distributions in MCV, the robustness to the number of training samples and noise power estimation is significantly improved compared with other SR−STAP methods for the MMV case.
• Incorporating a suitable VAE into MCV, a novel method called BAMCV is developed to accelerate the convergence of iterative procedures for estimating parameters. As the inference network is pre−trained off−line, BAMCV−STAP can realize the sparse reconstruction with lower computational loads and much fewer iterations compared with conventional SR−STAP methods.
• As demonstrated on both simulated and measured data, the final proposed method BAMCV−STAP can process space−time echoes in real−time without degrading clutter suppression performance.

The rest of this paper is organized as follows. Section 2 introduces the signal model for airborne radar. Section 3 details the proposed algorithm. Numerical experiments with both simulated and measured data are shown in Section 4 before discussion in Section 5. Section 6 presents the conclusion and future work of findings.

2. System Model and Derivation of Optimal Filters

For tractability purpose, a side−looking airborne phased−array radar with a uniformly linear antenna (ULA) [34,44] consisting of M half wavelength spaced elements on the airborne radar platform is considered. As shown in Figure 1, the platform is at altitude $h_p$ and moving with constant velocity $v_p$. $\phi$ and $\theta$ are the azimuth angle and the elevation angle, respectively. The radar transmits N pulses at a constant pulse repetition frequency (PRF) in a coherent processing interval (CPI). The model for the space−time clutter plus noise snapshot $\mathit{x}\in {\mathbb{C}}^{MN\times 1}$ from the CUT is [45]

$$\begin{array}{c}\begin{array}{cc}\hfill \mathit{x}& ={\mathit{x}}_c+\mathit{n}\hfill \\ \hfill & =\sum _{i=1}^{N_c}{\alpha }_i\mathit{v}\left(f_{d,i},f_{s,i}\right)+\mathit{n},\hfill \end{array}\end{array}$$

(1)

where ${\mathit{x}}_c\in {\mathbb{C}}^{MN\times 1}$ is the space−time clutter vector and $\mathit{n}\in {\mathbb{C}}^{MN\times 1}$ is the Gaussian white thermal noise whose variance is ${\sigma }^2$; $N_c$ denotes the number of independent clutter patches in the CUT; ${\alpha }_i,f_{d,i},f_{s,i},\mathit{v}\left(f_{d,i},f_{s,i}\right)\in {\mathbb{C}}^{MN\times 1}$ are the complex amplitude, the Doppler frequency, the spatial frequency and the space−time steering vector of the ith clutter patch, respectively. Considering the criterion of linearly constrained minimum variance (LCMV), the optimal STAP weight vector is [45]

$${\mathit{S}}_w=\kappa {\mathit{R}}^{-1}{\mathit{v}}_t,$$

(2)

where $\kappa ={\left({\mathit{v}}_t^H{\mathit{R}}^{-1}{\mathit{v}}_t\right)}^{-1}$ is the normalization constant, ${\mathit{v}}_t$ is the space−time steering vector of the target, and $\mathit{R}$ denotes the CCM, which is normally unknown and estimated from i.i.d training samples. ${(·)}^H$ denotes the conjugate operation. In conventional STAP algorithms, the CCM can be calculated by [45]

$$\mathit{R}=\frac{1}{L}\sum _{l=1}^L{\mathit{x}}_l{\mathit{x}}_l^H,$$

(3)

where L is the number of training samples and ${\mathit{x}}_l\in {\mathbb{C}}^{MN\times 1}(l=1,2,\dots ,L)$ denotes the space−time clutter plus noise vector of lth training sample. In order to achieve stable performance, L is required to be at least twice of $MN$. However, it is difficult to obtain sufficient i.i.d. training samples in practice. As a result, SR−STAP algorithms, which exhibit better performance using far fewer training samples, stand out from the crowd. They start by dividing the whole space−time plane uniformly into $K=N_s{N}_t$ bins, where $N_s={\rho }_{s}M\left({\rho }_s>1\right)$ and $N_t={\rho }_{t}N\left({\rho }_t>1\right)$ are the number of spatial frequency grid points and Doppler frequency grid points, respectively; ${\rho }_s$ and ${\rho }_d$ determine the resolution of the spatial frequency and Doppler frequency axes. As each grid point can be represented as a space−time steering vector ${\mathit{v}}_k(k=1,2,3,\dots ,K)$, STAP dictionary $\mathsf{\Phi }\in {\mathbb{C}}^{MN\times K}$ can be denoted as a collection of all these space−time steering vectors, i.e.,

$$\begin{array}{c}\hfill \mathsf{\Phi }=\left[{\mathit{v}}_1,{\mathit{v}}_2,\dots ,{\mathit{v}}_K\right].\end{array}$$

(4)

For the SMV case, the signal model in Equation (1) can be re−expressed as

$$\mathit{x}=\mathsf{\Phi }\mathit{w}+\mathit{n},$$

(5)

where $\mathit{w}\in {\mathbb{C}}^{K\times 1}$ denotes the space−time profile with non−zero elements representing the presence of the clutter. The formulation of the CCM for the SMV case can be expressed as

$$\mathit{R}=\mathsf{\Phi }\mathit{w}{\mathit{w}}^H{\mathsf{\Phi }}^H+{\sigma }^2\mathbf{I}.$$

(6)

For the MMV case, the signal model of L range cells $\mathit{X}=\left[{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_L\right]\in {\mathbb{C}}^{MN\times L}$ can be expressed as

$$\mathit{X}=\mathsf{\Phi }\mathit{W}+\mathit{N},$$

(7)

where $\mathit{W}=\left[{\mathit{w}}_1,{\mathit{w}}_2,\dots ,{\mathit{w}}_L\right]\in {\mathbb{C}}^{K\times L}$ is an unknown solution matrix, $\mathit{N}=\left[{\mathit{n}}_1,{\mathit{n}}_2,\dots ,{\mathit{n}}_L\right]\in {\mathbb{C}}^{MN\times L}$ is a Gaussian white thermal noise matrix. Inconsistent with the SMV case, row sparsity is encouraged simultaneously in the MMV case, i.e., each training sample in $\mathit{X}$ has the same indices of direction of arrival (DOA). The calculation of the CCM for the MMV case can be expressed as

$$\mathit{R}=\frac{1}{L}\mathsf{\Phi }\mathit{W}{\mathit{W}}^H{\mathsf{\Phi }}^H+{\sigma }^2\mathbf{I}.$$

(8)

The optimal solution to $\mathit{W}$ is discussed in the next section.

3. Proposed MCV−STAP Algorithm and BAMCV−STAP Algorithm

In this section, we first introduce the MCV approach by generalizing the complex−valued VRVM from the SMV case to the MMV case. Then we propose an MCV−STAP algorithm for designing the adaptive filters to suppress clutter and develop a Gibbs sampler to estimate parameters. To remedy the unattractive property requiring quantities of iterations in the Gibbs sampler of the MCV−STAP algorithm, we provide an efficient autoencoding variational inference network for the generative model built on the MCV approach and develop a VAE−structured method named BAMCV−STAP.

3.1. Derivation of MCV

We now derive the MCV approach by generalizing the complex−valued VRVM from the SMV case to the MMV case. Firstly, for the SMV case, assuming the noise is white complex Gaussian distribution with unknown power ${\sigma }^2$, we can obtain the Gaussian likelihood function of the signal model in Equation (5) as

$$p\left(\mathit{x}\mid \mathit{w},{\sigma }^2\right)\triangleq \mathcal{CN}\left(\mathsf{\Phi }\mathit{w},{\sigma }^2\mathbf{I}\right).$$

(9)

Drawing on the concept of automatic relevance determination (ARD) [46], we assign the Gaussian prior on the kth row of $\mathit{w}\in {\mathbb{R}}^K$ as

$$p\left(w_k\mid {\gamma }_k\right)\triangleq \mathcal{CN}\left(0,{\gamma }_k^{-1}\right)$$

(10)

where ${\gamma }_k$ is an unknown variance parameter. Combining these row priors, a full prior of variable $\mathit{w}$ can be expressed as

$$p\left(\mathit{w}\mid \mathit{\gamma }\right)\triangleq \prod _{k=1}^{K}p\left(w_k\mid {\gamma }_k\right)$$

(11)

where $\mathit{\gamma }={\left[{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_K\right]}^T\in {\mathbb{R}}_{+}^K$ is a vector of hyper−parameters. In this way, every weight has an individual hyper−parameter to moderate the strength of the prior. To ensure the parameter−free property of our method and avoid point estimation of $\mathit{\gamma }$ and ${\sigma }^2$ [36], the Gamma prior [47] over ${\gamma }_k$ as well as over the noise variance ${\sigma }^2$ in the signal model are defined, i.e.,

$$p\left(\mathit{\gamma }\right)=\prod _{k=1}^K\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left({\gamma }_k\mid a_k,b_k\right),$$

(12)

$$p\left(\beta \right)=\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\beta \mid c,d\right),$$

(13)

where $\beta ={\sigma }^{-2}$ and

$$\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left({\gamma }_k\right|a_k,b_k)=\Gamma {\left(a_k\right)}^{-1}{b_k}^{a_k}{{\gamma }_k}^{a_k-1}exp(-b_k{\gamma }_k).$$

(14)

$\Gamma \left(a\right)$ is the gamma function expressed as

$$\Gamma \left(a\right)={\int }_0^{\infty }{t}^{a-1}exp(-t)dt.$$

(15)

By providing posterior distribution estimates for $\mathit{\gamma }$ and $\beta$, our method increases the robustness of the noise power estimates while requiring fewer training samples.

Extending the above model for the SMV case to the MMV case to obtain the proposed model MCV, the Gaussian likelihood function of the signal model can be rewritten as

$$p\left(\mathit{X}\mid \mathit{W},\beta \right)\triangleq \mathcal{CN}\left(\mathsf{\Phi }\mathit{W},{\beta }^{-1}\mathit{I}\right).$$

(16)

The prior of the kth row of $\mathit{W}$ in Equation (7) denoted by ${\mathit{w}}_{k·}$ is

$$p\left({\mathit{w}}_{k·}\mid {\gamma }_k\right)\triangleq \mathcal{CN}\left(\mathbf{0},{{\gamma }_k}^{-1}\mathit{I}\right).$$

(17)

So the prior of the lth column of $\mathit{W}$ denoted by ${\mathit{w}}_{·l}$ is

$$p\left({\mathit{w}}_{·l}\mid \mathit{\gamma }\right)\triangleq \mathcal{CN}(\mathbf{0},{\mathbf{A}}^{-1})$$

(18)

where $\mathbf{A}=diag\left(\gamma \right)$ and $diag(·)$ denotes the diagnonal operation. Hence the full prior of $\mathit{W}$ is

$$p\left(\mathit{W}\right|\mathit{\gamma })\stackrel{\Delta }{=}\prod _{l=1}^{L}p\left({\mathit{w}}_{·l}\right|\mathit{\gamma }).$$

(19)

The priors of $\mathit{\gamma }$ and $\beta$ are the same as the case of SMV.

3.2. Proposed MCV−STAP Algorithm

The Gibbs sampler to estimate all parameters for MCV are derived from the SMV case in this subsection. Firstly, for the SMV case, the conditional posterior distribution over $\mathit{w}$ can be written as

$$\begin{array}{c}\begin{array}{cc}\hfill p\left(\mathit{w}\mid \mathit{x},\mathit{\gamma },\beta \right)& =\frac{p\left(\mathit{x}|\mathit{w},\beta \right)p\left(\mathit{w}\mid \mathit{\gamma }\right)}{p\left(\mathit{x}\mid \mathit{\gamma },\beta \right)}\hfill \\ \hfill & =\frac{{\left|\Sigma \right|}^{-1}}{{\pi }^K}\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\left(\mathit{w}-\mathit{\mu }\right)}^H{\Sigma }^{-1}\left(\mathit{w}-\mathit{\mu }\right)\right]\hfill \end{array}\end{array}$$

(20)

where the posterior covariance and mean are respectively

$$\Sigma ={\left(\beta \sum _{i=1}^{MN}{{\mathit{\varphi }}_i}^H{\mathit{\varphi }}_i+\mathbf{A}\right)}^{-1},$$

(21)

$$\mathit{\mu }=\beta \Sigma {\mathsf{\Phi }}^H\mathit{x},$$

(22)

with ${\mathit{\varphi }}_i$ denoting the i-th row vector of $\mathsf{\Phi }$.

Similarly, the conditional posterior distribution over $\mathit{\gamma }$ can be expressed as

$$\begin{array}{c}\begin{array}{cc}\hfill p\left(\mathit{\gamma }|\mathit{w},{\mathit{a}}_0,{\mathit{b}}_0\right)& =\frac{p\left(\mathit{w}\mid \mathit{\gamma }\right)p\left(\mathit{\gamma }\mid {\mathit{a}}_0,{\mathit{b}}_0\right)}{p\left(\mathit{w}\mid {\mathit{a}}_0,{\mathit{b}}_0\right)}\hfill \\ \hfill & =\prod _{k=1}^K\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left({\gamma }_k\mid a_k,b_k\right)\hfill \end{array}\end{array}$$

(23)

where the shape parameter and the scale parameter can be expressed as

$$\begin{array}{cc}\hfill a_k& =a_{k,0}+1,\hfill \\ \hfill b_k& =b_{k,0}+{w_k}^2.\hfill \end{array}$$

(24)

${\mathit{a}}_0={[a_{1,0},a_{2,0},...,a_{K,0}]}^T$ denotes the vector of initial values for shape parameters of $\mathit{\gamma }$, and ${\mathit{b}}_0={[b_{1,0},b_{2,0},...,b_{K,0}]}^T$ denotes the vector of initial values for scale parameters of $\mathit{\gamma }$. The conditional posterior distribution over $\beta$ is

$$\begin{array}{c}\begin{array}{cc}\hfill p\left(\beta \mid \mathit{x},\mathit{w},c_0,d_0\right)& =\frac{p\left(\mathit{x}\mid \mathit{w},\beta \right)p\left(\beta \mid c_0,d_0\right)}{p\left(\mathit{x}\mid \mathit{w},c_0,d_0\right)}\hfill \\ \hfill & =\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\beta \mid c,d\right),\hfill \end{array}\end{array}$$

(25)

where the shape parameter and the scale parameter respectively become

$$\begin{array}{c}\hfill \begin{array}{c}c=c_0+MN,\\ d=d_0+{\mathit{x}}^H\mathit{x}-2{\mathit{w}}^H{\mathsf{\Phi }}^H\mathit{x}+{\displaystyle \sum _{i=1}^{MN}}{\mathit{\varphi }}_i\left(\mathit{w}{\mathit{w}}^H\right){\mathit{\varphi }}_i^H.\end{array}\end{array}$$

(26)

$c_0$ denotes the initial value for shape parameters of $\beta$ and $d_0$ denotes the initial values for scale parameters of $\beta$.

Extending the above process to the MMV case, we develop the Gibbs sampler for MCV here. The conditional posterior distribution over the latent variable $\mathit{W}$ is thus given by

$$\begin{array}{c}\begin{array}{cc}\hfill p\left(\mathit{W}\mid \mathit{X},\mathit{\gamma },\beta \right)& =\frac{p\left(\mathit{X}|\mathit{W},\beta \right)p\left(\mathit{W}\mid \mathit{\gamma }\right)}{p\left(\mathit{X}\mid \mathit{\gamma },\beta \right)}\hfill \\ \hfill & =\frac{{\left|{\Sigma }_{\mathbf{M}}\right|}^{-1}}{{\pi }^K}\mathrm{e}\mathrm{x}\mathrm{p}\left\{\mathrm{t}\mathrm{r}\left[-{\left(\mathit{W}-{\mathit{\mu }}_{\mathbf{M}}\right)}^H{{\Sigma }_{\mathbf{M}}}^{-1}\left(\mathit{W}-{\mathit{\mu }}_{\mathbf{M}}\right)\right]\right\}\hfill \end{array}\end{array}$$

(27)

where the covariance and mean are respectively

$${\Sigma }_{\mathbf{M}}={\left(\mathbf{fi}\beta \sum _{\mathbf{i}=\mathbf{1}}^{\mathbf{MN}}{{\mathit{\varphi }}_{\mathbf{i}}}^{\mathbf{H}}{\mathit{\varphi }}_{\mathbf{i}}+\mathbf{A}\right)}^{-\mathbf{1}},$$

(28)

$${\mathit{\mu }}_{\mathbf{M}}=\beta {\Sigma }_{\mathbf{M}}{\mathsf{\Phi }}^H\mathit{X},$$

(29)

and $tr(·)$ denotes the trace of a square matrix. The conditional posterior distribution over $\mathit{\gamma }$ can be expressed as

$$\begin{array}{c}\begin{array}{cc}\hfill p\left(\mathit{\gamma }|\mathit{W},{\mathit{a}}_0,{\mathit{b}}_0\right)& =\frac{p\left(\mathit{W}\mid \mathit{\gamma }\right)p\left(\mathit{\gamma }\mid {\mathit{a}}_0,{\mathit{b}}_0\right)}{p\left(\mathit{W}\mid {\mathit{a}}_0,{\mathit{b}}_0\right)}\hfill \\ \hfill & =\prod _{k=1}^K\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left({\gamma }_k\mid a_k,b_k\right)\hfill \end{array}\end{array}$$

(30)

where the shape parameter and the scale parameter can be expressed as

$$\begin{array}{cc}\hfill a_k& =a_{k,0}+L,\hfill \\ \hfill b_k& =b_{k,0}+{∥{\mathit{w}}_k.∥}_2^2.\hfill \end{array}$$

(31)

The conditional posterior distribution over $\beta$ is

$$\begin{array}{c}\begin{array}{cc}\hfill p\left(\beta \mid \mathit{X},\mathit{W},c_0,d_0\right)& =\frac{p\left(\mathit{X}\mid \mathit{W},\beta \right)p\left(\beta \mid c_0,d_0\right)}{p\left(\mathit{X}\mid \mathit{W},c_0,d_0\right)}\hfill \\ \hfill & =\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\beta \mid c,d\right),\hfill \end{array}\end{array}$$

(32)

where the shape parameter and the scale parameter respectively become

$$\begin{array}{c}\hfill \begin{array}{c}c=c_0+MNL,\\ d=d_0+{\displaystyle \sum _{l=1}^L}{{\mathit{x}}_l}^H{\mathit{x}}_l-2{\displaystyle \sum _{l=1}^L}{\mathit{W}}^H{\mathsf{\Phi }}^H{\mathit{x}}_l+{\displaystyle \sum _{i=1}^{MN}}{\mathit{\varphi }}_i\left(\mathit{W}{\mathit{W}}^H\right){\mathit{\varphi }}_i^H.\end{array}\end{array}$$

(33)

As conditional distributions over all parameters are specified clearly in the above, the space−time profiles estimated accurately by the Gibbs sampler of MCV will be adopted to design filters to suppress the clutter. For clarity, the pseudo−codes of the abovementioned MCV−STAP algorithm are provided in Algorithm 1.

Algorithm 1 MCV−STAP algorithm.

Step1: Give initial values $a_{k,0}=b_{k,0}=c_0=d_0={10}^{-6},k=1,2,3,...,K$.

Step2: Sample ${\mathit{\gamma }}^{\left(0\right)}$ from Equation (12) and sample ${\beta }^{\left(0\right)}$ from Equation (13).

Step3: For $t=1,2,...$ do

Calculate ${\mathit{\mu }}_{\mathbf{M}}^{\left(\mathbf{t}\right)}$ and ${\Sigma }_{\mathbf{M}}^{\left(\mathbf{t}\right)}$ via Equations (28) and (29);

Sample ${\mathit{W}}^{\left(t\right)}$ from Equation (20);

Calculate ${\mathit{a}}^{\left(t\right)}$ and ${\mathit{b}}^{\left(t\right)}$ via Equation (31);

Sample ${\mathit{\gamma }}^{\left(t\right)}$ from Equation (30);

Calculate $c^{\left(t\right)}$ and $d^{\left(t\right)}$ via Equation (33);

Sample ${\beta }^{\left(t\right)}$ from Equation (32);

Check for convergence.

End For

Denote the last iter number as T.

Step4: Assume $\mathit{W}={\mathit{\mu }}_M^{\left(T\right)}$, calculate the CCM via Equation (8) and the space−time adaptive optimal

weight vector via Equation (1).

Step5: Denote the CUT as ${\mathit{x}}_{vt}$, the output of MCV−STAP algorithm solved by Gibbs sampling

is $\mathit{y}={{\mathit{S}}_{vt}}^H{\mathit{x}}_{vt}$.

While having closed−form updating equations and providing accurate estimation for parameters, the Gibbs sampler is still not attractive since it takes a large number of iterations to infer the sparse representation during the testing phase, which hinders real−time processing of the incoming CUT and motivates us to construct a VAE with fast testing, as described below.

3.3. BAMCV−STAP Algorithm

In this section, to ensure the real−time processing of MCV−STAP, we extend it into a VAE−structured method named BAMCV−STAP to map the observation $\mathit{x}$ directly to the latent representation for out−of−sample prediction. Specifically, the generative model (decoder) of BAMCV−STAP is built based on MCV as shown in Figure 2a, and the inference network (encoder) is accomplished by the neural network with the structure shown in Figure 2b. Instead of needing expensive iterative inference schemes with Gibbs samplers introduced above, the variational inference network of BAMCV allows us to efficiently realize sparse reconstructions by fitting an approximate inference model to the unknown posterior using standard stochastic gradient methods. According to the theory of mean−field variational Bayes [48,49,50], the ELBO of BAMCV can be expressed as

$$\begin{array}{c}\hfill L=-{\mathbb{D}}_{KL}\left[q\left(\mathit{\gamma }\right)\left|\right|p\left(\mathit{\gamma }\right)\right]-{\mathbb{D}}_{KL}\left[q\left(\mathit{w}\right|\mathit{x}\left)\right|\left|p\right(\mathit{w}\left|\mathit{\gamma }\right)\right]-{\mathbb{D}}_{KL}\left[q\left(\beta \right)\left|\right|p\left(\beta \right)\right]\\ \hfill +{\mathbb{E}}_{q(\mathit{w},\mathit{\gamma },\beta |\mathit{x})}\left[p\left(\mathit{x}\right|\mathit{w},\mathit{\gamma },\beta )\right].\end{array}$$

(34)

where the approximate conditional posterior distributions $q\left(\mathit{\gamma }\right)$, $q\left(\mathit{w}|\mathit{x}\right)$ and $q\left(\beta \right)$ can be written as

$$\begin{array}{c}\hfill q\left(\mathit{\gamma }\right)=\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left({\mathit{a}}_h,{\mathit{b}}_h\right),\\ \hfill q\left(\mathit{w}|\mathit{x}\right)=\mathcal{CN}\left({\mathit{\mu }}_h,{\Sigma }_h\right),\\ \hfill q\left(\beta \right)=\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(c_h,d_h\right).\end{array}$$

(35)

Noticeably, ${\mathbb{D}}_{KL}\left[q(·)\right|\left|p(·)\right]$ denotes the Kullback–Leibler (KL) divergence between $q(·)$ and $p(·)$.

Different from traditional Gaussian distribution based VAEs, the gamma distribution, modeling the sparsity and satisfying non−negative constraint, is defined over $\mathit{\gamma }$ and $\beta$ here. However, it is hard to compute the gradient of the ELBO to $\mathit{\gamma }$ and $\beta$ due to the difficulty of reparameterizing gamma distributed random variables [51,52]. To reparameterize $\mathit{\gamma }$ and $\beta$ easily without deviating from gamma distributions, the Weibull distribution is considered an ideal choice to replace the gamma distribution, as demonstrated in [53]. Consequently, gamma distributed conditional posteriors, $q\left(\mathit{\gamma }\right)$ and $q\left(\beta \right)$, are approximated by Weibull distributions expressed as

$$\begin{array}{cc}\hfill q\left(\mathit{\gamma }\right)& =\prod _{i=1}^{K}q\left({\gamma }_i\right)=\prod _{i=1}^K\mathrm{W}\mathrm{e}\mathrm{i}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{l}\left(k_{{\gamma }_i},{\lambda }_{{\gamma }_i}\right),\hfill \\ \hfill q\left(\beta \right)& =\mathrm{W}\mathrm{e}\mathrm{i}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{l}\left(k_{\beta },{\lambda }_{\beta }\right).\hfill \end{array}$$

(36)

where $k_{{\gamma }_i}$ and ${\lambda }_{{\gamma }_i}$ are the parameters of ${\gamma }_i$, and ${\lambda }_{\beta }$ and $k_{\beta }$ are the parameters of ${\gamma }_i$. Latent variables ${\gamma }_i$ and $\beta$ can be easily reparameterized as [53,54]

$$\begin{array}{c}\hfill {\gamma }_i={\lambda }_{{\gamma }_i}{(-ln(1-{\epsilon }_{{\gamma }_i}))}^{\frac{1}{k_{{\gamma }_i}}},{\epsilon }_{{\gamma }_i}\sim \mathrm{U}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}(0,1),\\ \hfill \beta ={\lambda }_{\beta }{(-ln(1-{\epsilon }_{\beta }))}^{\frac{1}{k_{\beta }}},{\epsilon }_{\beta }\sim \mathrm{U}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}(0,1).\end{array}$$

(37)

In addition, the Gaussian distribution based variable $\mathit{w}$ can be simply reparameterized as [40]

$$\begin{array}{c}\hfill {\mathit{w}}^{\mathit{r}}={\mathit{\mu }}_h^r+{\Sigma }_h\ast {\mathit{z}}^{\mathit{r}},{\mathit{z}}^{\mathit{r}}\sim \mathcal{N}(\mathbf{0},\mathit{I}),\\ \hfill {\mathit{w}}^i={\mathit{\mu }}_h^i+{\Sigma }_h\ast {\mathit{z}}^{\mathit{i}},{\mathit{z}}^{\mathit{i}}\sim \mathcal{N}(\mathbf{0},\mathit{I})\end{array}$$

(38)

where the superscript r of all parameters represents their real part and the superscript i of all parameters represents their imaginary part. $(\ast )$ denotes matrix multiplication. According to the illustration in [53], the KL divergence from the gamma distribution to the Weibull distribution, which is the first term of the ELBO in Equation (34), can be expressed analytically as

$$\begin{array}{cc}\hfill {\mathbb{D}}_{KL}\left[q\left(\mathit{\gamma }\right)\right|\left|p\left(\mathit{\gamma }\right)\right]& ={\mathbb{D}}_{KL}[\mathrm{W}\mathrm{e}\mathrm{i}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{l}({\mathit{k}}_{\mathit{\gamma }},{\mathit{\lambda }}_{\mathit{\gamma }})\left|\right|\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}({\mathit{a}}_0,{\mathit{b}}_0)]\hfill \\ \hfill & =\sum _{i=1}^K\frac{{\gamma }_e{a}_{i,0}}{k_{{\gamma }_i}}-a_{i,0}ln{\lambda }_{{\gamma }_i}+ln{k}_{{\gamma }_i}+b_{i,0}{\lambda }_{{\gamma }_i}\Gamma (1+\frac{1}{k_{{\gamma }_i}})-{\gamma }_e-1-a_{i,0}ln\left(b_{i,0}\right)+ln\Gamma \left(a_{i,0}\right).\hfill \end{array}$$

(39)

where ${\gamma }_e$ is the Euler–Mascheroni constant. The second term of the ELBO can be expressed analytically as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbb{D}}_{KL}\left[q\left(\mathit{w}\right|\mathit{x})\right|\left|p\left(\mathit{w}\right|\mathit{\gamma })\right]& ={\mathbb{D}}_{KL}[\mathcal{CN}({\mathit{\mu }}_h,{\Sigma }_h)\left|\right|\mathcal{CN}(\mathbf{0},{\mathbf{A}}^{-\mathbf{1}})]\hfill \\ \hfill & =\sum _{i=1}^K(1+log{\epsilon }_{kk}-|{\mu }_{h,i}{|}^2-{\epsilon }_{kk}).\hfill \end{array}\end{array}$$

(40)

${\epsilon }_{kk}$ means the k-th diagonal element of ${\Sigma }_h$, and ${\mu }_{h,i}$ denotes the i-th element of ${\mathit{\mu }}_h$. $|·|$ means the modulus operation. The third term and the last term of the ELBO can be expressed analytically as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbb{D}}_{KL}\left[q\left(\beta \right)\right|\left|p\left(\beta \right)\right]& ={\mathbb{D}}_{KL}[\mathrm{W}\mathrm{e}\mathrm{i}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{l}(k_{\beta },{\lambda }_{\beta })\left|\right|\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}(c_0,d_0)]\hfill \\ & =\frac{{\gamma }_e{c}_0}{k_{\beta }}-c_{0}ln{\lambda }_{\beta }+ln{k}_{\beta }+d_0{\lambda }_{\beta }\Gamma (1+\frac{1}{k_{\beta }})-{\gamma }_e-1-c_{0}ln\left(d_0\right)+ln\Gamma \left(c_0\right).\hfill \end{array}\end{array}$$

(41)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbb{E}}_{q(\mathit{w},\mathit{\gamma },\beta |\mathit{x})}\left[p\left(\mathit{x}\right|\mathit{w},\mathit{\gamma },\beta )\right]={∥\mathit{x}-\mathsf{\Phi }\mathit{w}∥}_2^2\end{array}\end{array}$$

(42)

As illustrated in Figure 2a, all the parameters relative to conditional distributions are transformed from the observation ${\mathit{x}}_{input}$ with the neural network. Specifically,

$$\begin{array}{c}{\mathit{\mu }}_h^r=\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left({\mathbf{C}}_{\mathit{\mu }}^{\mathit{r}}{\mathit{h}}_1+g_{\mathit{\mu }}^{\mathit{r}}\right),\hfill \\ {\mathit{\mu }}_h^i=\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left({\mathbf{C}}_{\mathit{\mu }}^{\mathit{i}}{\mathit{h}}_1+g_{\mathit{\mu }}^{\mathit{i}}\right),\hfill \\ {\mathit{\mu }}_h={\mathit{\mu }}_h^r+j{\mathit{\mu }}_h^i,\hfill \end{array}$$

(43)

where ${\mathit{\mu }}_h$ is the mean vector of $\mathit{w}$,

$$\begin{array}{c}{\mathit{e}}_h=\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left({\mathbf{C}}_{\mathit{e}}{\mathit{h}}_1+{\mathit{g}}_{\mathit{e}},\right)\hfill \end{array}$$

(44)

where ${\mathit{e}}_h=diag\left({\Sigma }_h\right)$, and

$$\begin{array}{c}{\mathit{h}}_1=\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left({\mathbf{C}}_{1h}{\mathit{x}}_{input}+{\mathit{g}}_{1h}\right)\hfill \end{array}$$

(45)

where $\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}(·)$ denotes the non−linear activation function. Note that ${\mathit{x}}_{input}$ denotes the concatenation of the real part and the imaginary part of training data, as $\mathit{x}$ and latent representation $\mathit{w}$ are complex−valued. Similarly,

$$\begin{array}{c}{\mathit{k}}_{\gamma }=\mathrm{S}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({\mathit{C}}_a{\mathit{h}}_2+{\mathit{g}}_a\right),{\mathit{\lambda }}_{\mathit{\gamma }}=\mathrm{S}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({\mathit{C}}_b{\mathit{h}}_2+{\mathit{g}}_b\right)\hfill \\ k_{\beta }=\mathrm{S}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({\mathit{C}}_c{\mathit{h}}_{\beta }+g_c\right),{\lambda }_{\beta }=\mathrm{S}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({\mathit{C}}_d{\mathit{h}}_{\beta }+g_d\right)\hfill \\ \hfill \end{array}$$

(46)

$$\begin{array}{c}{\mathit{h}}_2=\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left({\mathbf{C}}_{2h}{\mathit{h}}_1+{\mathit{g}}_{2h}\right)\hfill \\ {\mathit{h}}_{\beta }=\mathrm{R}\mathrm{e}\mathrm{L}\mathrm{U}\left({\mathbf{C}}_{\beta }{\mathit{x}}_{input}+{\mathit{g}}_{\beta }\right)\hfill \end{array}$$

(47)

where ${\mathit{k}}_{\gamma }=[k_{{\gamma }_1},k_{{\gamma }_2},...,k_{{\gamma }_K}]\in {\mathbb{R}}^{K\times 1}$, ${\mathit{\lambda }}_{\gamma }=[{\lambda }_{{\gamma }_1},{\lambda }_{{\gamma }_2},...,{\lambda }_{{\gamma }_K}]\in {\mathbb{R}}^{K\times 1}$ and $\mathrm{S}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}(·)$ applies $\mathrm{l}\mathrm{o}\mathrm{g}[1+\mathrm{e}\mathrm{x}\mathrm{p}(·\left)\right]$ non-linearity to each element to ensure positive Weibull shape and scale parameters.

Obviously, ${\mathit{C}}_{\mathit{\mu }}^{\mathit{r}}$, ${\mathit{C}}_{\mathit{\mu }}^{\mathit{i}}$, ${\mathit{C}}_{\mathit{e}}$, ${\mathit{C}}_{\mathbf{1}\mathbf{h}}$, ${\mathit{g}}_{\mathit{\mu }}^{\mathit{r}}$, ${\mathit{g}}_{\mathit{\mu }}^{\mathit{i}}$, ${\mathit{g}}_{\mathit{e}}$, ${\mathit{g}}_{\mathbf{1}\mathbf{h}}$, ${\mathit{C}}_a$, ${\mathit{C}}_b$, ${\mathit{C}}_c$, ${\mathit{C}}_d$, ${\mathit{C}}_{2h}$, ${\mathit{C}}_{\beta }$, ${\mathit{g}}_a$, ${\mathit{g}}_b$, $g_c$, $g_d$, ${\mathit{g}}_{2h}$ and ${\mathit{g}}_{\beta }$ are parameters in the neural network. The superscript r of all parameters represents the real part, and the superscript i of all parameters represents the imaginary part. For clarity, the dimensions of parameters in the inference network are summarized in

Table 1. Parameters in the inference network.

Parameters	Dimensions	Parameters	Dimensions	Parameters	Dimensions
${\mathit{C}}_{\mathit{\mu }}^{\mathit{r}}$	$K\times K$	${\mathit{g}}_{1h}$	$2MN\times 1$	${\mathit{g}}_a$	$K\times 1$
${\mathit{C}}_{\mathit{\mu }}^{\mathit{i}}$	$K\times K$	${\mathit{C}}_a$	$K\times K$	${\mathit{g}}_b$	$K\times 1$
${\mathit{C}}_{\mathit{e}}$	$K\times K$	${\mathit{C}}_b$	$K\times K$	$g_c$	$1\times 1$
${\mathit{C}}_{1h}$	$2MN\times K$	${\mathit{C}}_c$	$1\times 2MN$	$g_d$	$1\times 1$
${\mathit{g}}_{\mathit{\mu }}^{\mathit{r}}$	$K\times 1$	${\mathit{C}}_d$	$1\times 2MN$	${\mathit{g}}_{2h}$	$K\times 1$
${\mathit{g}}_{\mathit{\mu }}^{\mathit{i}}$	$K\times 1$	${\mathit{C}}_{2h}$	$K\times K$	${\mathit{g}}_{\beta }$	$2MN\times 1$
${\mathit{g}}_{\mathit{e}}$	$K\times 1$	${\mathit{C}}_{\beta }$	$2MN\times 2MN$

.

Above all, due to the simple reparameterization tricks, our proposed VAE−structured MCV−STAP methods can obtain the latent representations directly with the aid of an autoencoding inference network.

Moreover, to make the proposed method further attractive in real−time processing, we pre−train the inference network off−line with simulated data preparing for sequential fine−tuning facing with realistic clutter data to be processed [55]. The main steps of the proposed BAMCV for STAP (BAMCV−STAP) algorithm are summarized in Algorithm 2. For clarity, the flow chart of BAMCV is plotted in Figure 3.

Algorithm 2 BAMCV−STAP algorithm.

Step1: Simulate data for pre−training using radar system parameters of realistic clutter data under

test. Denote the simulated dataset as $D_{pre}$ and the realistic dataset as $D_{test}$.

Step2: Pre−train the inference network with dataset $D_{pre}$ off−line until the ELBO converges.

Step3: Select CUT in $D_{test}$ and choose L training samples around CUT.

Step4: Fine−tune the inference network with the selected L training samples until the ELBO

converges again. Assume $\mathit{W}={\mathit{\mu }}_h$.

Step6: Estimate the CCM via Equation (8) and the space−time adaptive optimal weight vector via Equation (2).

Denote the weight vector as ${\mathit{S}}_{vt}.$

Step7: Denote data in CUT as ${\mathit{x}}_{cut}$, the output of BAMCV−STAP solved by the inference network is

$y={\mathit{S}}_{vt}^H{\mathit{x}}_{cut}$.

Compared with MCV, BAMCV has two explicit advantages which can be summarized as follows:

• Decreasing the computational loads. Inspired by the gradient ascent scheme for parameters optimization, the parameters of BAMCV are updated by the backpropagation of the gradient, which only involves some linear operations, exponential operations and logarithmic operations, instead of involving complex inverse operations and multiple samplings appeared in MCV.
• Improving the convergence rate for testing. As the inference network is pre−trained off−line, only a few iterations are taken in BAMCV to obtain the recovery results of observations rather than a large number of time−consuming iterations per testing in MCV.

Experimental results shown in the following section will verify these advantages.

4. Numerical Results

In this section, we evaluate the proposed algorithms, MCV−STAP and BAMCV−STAP, using both simulated data and measured data in terms of clutter suppression performance, time efficiency and computational complexity. Significantly, the metric to evaluate clutter suppression performance for the simulated data is the signal to interference plus noise ratio (SINR) loss [45], which is defined as follows

$$L_{\mathrm{SINR}}=\frac{{\sigma }_n^2}{MN}\frac{{\left|{\mathit{S}}_w^H{\mathit{v}}_t\right|}^2}{{\mathit{S}}_w^H\mathit{R}{\mathit{S}}_w}.$$

(48)

where ${\mathit{S}}_w$ denotes the STAP weight vector, $\mathit{R}$ denotes the true CCM, ${\mathit{v}}_t$ denotes the target steering vector and ${\sigma }_n^2$ denotes the actual noise power. In addition, as the true CCM is unable to be obtained for the measured data, we evaluate the target detection performance for targets with normalized Doppler frequency of 0.1 by the probability of detection (PD) versus signal to noise ratio (SNR) curves, which are achieved by utilizing the adaptive matched filter (AMF) detector [56] and averaged over ${10}^5$ Monte Carlo trials. The probability of the false alarm rate (PFA) is set as ${10}^{-6}$, the PD is estimated based on ${10}^6$ randomly sampled CUTs and the PD to SNR curves are averaged over ${10}^5$ Monte Carlo trials. SNR is defined as the ratio of signal power (for a single pulse per antenna element) to the rough estimate of noise power (the average energy of echoes within the support region) here. In addition, to evaluate the time efficiency, we evaluate the time spent on the convergence rate per testing. The computational complexity is also evaluated by the comparison of computational loads per testing. For clarity, the performance metrics used in this paper are summarized in

Table 2. Summary of Performance Metrics.

	Performance Metric	Simulated Data	Measured Data
Clutter Suppression	SINR loss	√
	PD		√
Real−time Processing	Convergence rate	√	√
	Computational loads	√	√

, in which the check mark means the metric is used in the data and the blank means not.

Moreover, the proposed algorithm is compared with some other classical methods especially for the MMV case, such as loading sample matrix inversion (LSMI) for STAP [57,58], multiple orthogonal matching pursuit (M−OMP) for STAP [59], M−SBL for STAP and multiple fast converging SBL (M−FCSBL) for STAP [34]. The LSMI for STAP [57,58], which is the representative of classical STAP algorithm, is an adaptive algorithm imposing diagonal loading (DL) to improve the robustness of conventional Capon beamformers but fails to achieve satisfactory performance unless the number of training samples is more than twice the clutter rank. The M−OMP for STAP, orthogonal matching pursuit (OMP) [25] for the MMV case, obtains sparse representations through a greedy algorithm but is particularly unfriendly to the highly coherent dictionary. As the classical technique to solve sparse recovery issues utilizing the greedy strategy, M−OMP is chosen as one of comparison methods in this paper. The M−SBL for STAP, suppressing the clutter imposing M−SBL algorithm based on Bayesian theory as well as ours, receives much attention in the field of STAP for the excellent performance with a highly coherent dictionary but potentially requires a large number of iterations per testing leading to unsuccessful real−time processing. The M−FCSBL for STAP, a fast converging M−SBL algorithm, improves the convergence by incorporating a simple approximation term and achieves a favorable recovery result with only a few iterations. However, the experimental results reveal that this approximation term leads to sub−optimal recovery results, as shown in its convergence curve.

Furthermore, the elements of all weight matrices in the neural network employed in this paper are initialized with Gaussian distributions with standard deviations set to 0.1, and all bias terms are set to zero. During the pre−training stage, we set the mini−batch size to 32. For optimization, the Adam optimizer [60] with a learning rate of ${10}^{-4}$ is used. On a Pentium PC with a 3.7 GHz CPU and 64 GB RAM, the inference network is pre−trained and fine−tuning in non−optimized Python software before being moved to MATLAB for testing. The amount of samples in the pre−training data set is ten times the number of parameters in the neural network to ensure the generalization of the pre−trained neural network. The settings of each layer in the neural network have been stated in Table 1.

4.1. Simulated Data

In this paper, a side−looking airborne phased array radar with a uniformly linear antenna (ULA), whose inter−element spacing is half wavelength, is considered. The parameters of radar system are shown in

Table 3. Radar System Parameters of Simulated data.

Parameters	Value	Parameters	Value
Carrier frequency (Hz)	1.25 G	Platform velocity (m/s)	125
Bandwidth (Hz)	2.5 M	Platform height (m)	6000
Mainbeam azimuth (∘)	0	Pulse number in one CPI	8
Mainbeam elevation (∘)	0	Antenna elements number	8
Pulse repetition frequency (Hz)	2000	Range cell number	400

. In the simulations, we set $N_s=4M$, $N_t=4N$, ${\sigma }_n^2=1$ and 180 clutter patches uniformly distributed from $-\pi /2$ to $\pi /2$. According to Algorithm 2, the data set $D_{pre}$ and $D_{test}$ for BAMCV−STAP are simulated with the same radar system parameters, respectively. Obviously, MCV−STAP is executed on $D_{test}$ directly without need of $D_{pre}$. To make the inference network easily generalized to the training samples of CUT, We simulated as many training samples as ten times the number of parameters with random clutter to signal ratio (CNR) sampled from 30 dB to 70 dB in $D_{pre}$. As for $D_{test}$, to investigate the performance of all methods under the heterogeneous clutter scenario conveniently, we assume the 150−th range cell as CUT which is located at the interested area with the average CNR = 45 dB and around by only 10 training samples satisfying the i.i.d constraint, as shown in Figure 4. Compared with the clairvoyant spectrum calculated by known CCM in Figure 5a, the clutter spectrum estimated by 128 training samples, as many as twice the DOF, exhibits too large value along the clutter ridge due to the heterogeneity of samples, as depicted in Figure 5b. Noticeably, the clutter suppression performance with the different number of i.i.d training samples in the interested area will also be investigated later.

4.1.1. Analysis of Clutter Suppression Performance

The clutter suppression performance of the proposed method is evaluated here. The clutter spectrum of clutter data under test utilizing various approaches are depicted in Figure 6. The number of training samples is fixed to 10 since the number of i.i.d training samples is set as 10 according to Figure 4. As is customary in LSMI, the load factor is set at 1. The clutter spectrum computed by MCV in Figure 6f and BAMCV in Figure 6h is the closest to the ideal one in Figure 6a in terms of both position and value. As the posterior distributions are strictly derived via Bayesian theory, the CCM computed using Gibbs sampling is as close to the clairvoyant one as BAMCV−STAP. Because of limited i.i.d. training samples, the approximated spectrum by M−OMP in Figure 6c and LSMI in Figure 6b expose their powerlessness in clutter subspace estimation. Instead, the clutter spectrum obtained by M−SBL, as shown in Figure 6d, is extremely close to the ideal one, with a small spreading and a few noise peaks due to imprecise noise power assessment. In Figure 6e, M−FCSBL [34] achieves an accurate clutter subspace and performs no noise peaks. However, upon further investigation of these SR−STAP approaches, as shown in Figure 7, we discover that values along the clutter ridge in the clutter spectrum calculated by M−FCSBL are not as well as expected. The clutter ridge formed by our method, on the other hand, is closer to the ideal result.

Moreover, to demonstrate the importance of the pre−training to the inference network, we utilize the inference network after pre−training to recover the space−time profile of clutter data under test, as shown in Figure 6g. Although the significant values along the clutter ridge are not accurate enough, the sparsity of the clutter spectrum and the shape of the clutter ridge are initialized successfully as the preparation for sequential fine−tuning. Apart from the sparse subspace, the noise power estimation is also an important issue for the CCM, according to Equation (8). Consequently, we enumerate the final results estimated by all comparison methods in

Table 4. Noise power estimation.

True Noise Power	M−SBL	M−FCSBL	MCV	BAMCV
0.01	$2.49\times {10}^{-10}$	0.0111	0.00973	0.0101
0.1	$4.15\times {10}^{-9}$	0.1197	0.1061	0.1103
1	$7.92\times {10}^{-8}$	1.0117	0.9937	0.9925
5	$2.98\times {10}^{-6}$	4.7692	4.9793	4.9001
10	$7.57\times {10}^{-6}$	11.1382	9.9992	9.9856

. Clearly, our methods outperform M−SBL and M−FCSBL, proving the robustness to the estimation for noise power.

When the number of training samples is 10, the clutter suppression performance measured by SINR loss is evaluated in Figure 8. The optimal SINR loss curve is completed using the clairvoyant CCM. In comparison to previous Bayesian SR−STAP approaches, LSMI and M−OMP display severe performance loss in the sidelobe region, as depicted in Figure 8a. Specifically, the output SINR loss of MCV−STAP is the closest to the optimal one among various Bayesian SR−STAP methods, as shown in Figure 8b (which is an enlargement of Figure 8a). BAMCV−STAP achieves the lower SINR loss than MCV−STAP as shown in Figure 8 but still surpasses other comparison methods.

Besides, we compare the effect of reducing the i.i.d training samples on the SINR loss of each method, and the results are shown in Figure 9. The average SINR loss is defined in this section as the mean of the SINR loss values across the entire normalized Doppler frequency range. According to Figure 9, on average sense, when L is a large value, all the methods have almost equivalent performance. But with an decrease in L, the Bayesian type SR−STAP methods, including M−SBL−STAP, MCV−STAP and BAMCV−STAP, have the larger SINR loss. Notably, our methods, both MCV−STAP and BAMCV−STAP, maintain stable performance as the number of training samples decreases and achieve a 3 dB higher SINR loss than M−SBL−STAP with only two training samples, demonstrating the robustness of the proposed methods to the number of training samples.

4.1.2. Analysis of Computational Loads and Convergence Rate

To study the real−time processing of several SR−STAP algorithms with similarly superior performance i.e., M−SBL, M−FCSBL and our methods, we detail the analysis of the convergence rate and computational loads in Figure 10.

For the on−line phase, convergence curves, utilizing the negative log−likelihood as the unified loss function [34] for fair comparison, are shown in Figure 10a, i.e.,

$$\begin{array}{c}{L}_{\mathrm{cost}}=ln\left|\mathit{M}\right|+\mathrm{t}\mathrm{r}\left({\mathit{M}}^{-1}{\mathit{R}}_{\mathrm{cost}}\right)\hfill \end{array}$$

(49)

where

$$\mathit{M}={\mathsf{\Phi }}^H{\mathbf{A}}^{-\mathbf{1}}\mathsf{\Phi }+{\beta }^{-1}\mathit{I}$$

$${\mathit{R}}_{\mathrm{cost}}=\frac{1}{L}\mathit{X}{\mathit{X}}^H.$$

Different from MCV converging as slowly as M−SBL, BAMCV, which is pre−trained with lots of simulated data off−line, exhibits the dominant superiority at the beginning of the iteration. Then, after fine−tuning with several snapshots around CUT over a few iterations, BAMCV converges to a stable and optimal state. The time spent on each iteration of each method is also listed in

Table 5. Comparison of running time per iteration using simulated data.

Approach	Running Time Per Iteration (s)
M−SBL	$0.125$
M−FCSBL	$0.055$
MCV	$0.103$
BAMCV	$0.011$

. Indeed, since it uses a gradient ascent scheme without any inverse operation or sampling, BAMCV takes much less time per iteration than the other three approaches. However, M−SBL and MCV require not only a large number of iterations but also a long time for each iteration. Furthermore, even though M−FCSBL quickly converges to its steady−state value, the final solution is not optimal, which is exactly coincident with the estimated clutter spectrum in Figure 6e and Figure 7.

In addition, the computational loads of M−SBL, M−FCSBL, MCV and BAMCV are detailed in

Table 6. Comparison of computational cost.

Approach	Computational Loads
M−SBL	$({\left(K\right)}^3+{\left(MN\right)}^3+3K^{2}MN+[2{\left(MN\right)}^2+2MNL+L+1]K+MNL+1)N_s$
M−FCSBL	$({\left(MN\right)}^3+K^{2}MN+6K{\left(MN\right)}^2+(2KL+2K)MN+1)N_f$
MCV	$({\left(MN\right)}^3+3K^{2}MN+4K{\left(MN\right)}^2+(L^2+L+1)MNK+MN{L}^2+MNL+1)N_g$
BAMCV	$2KMNL+5K^{2}L+4{\left(MN\right)}^{2}L+4MNL$

. Specifically, $N_s$ is the iteration number of the M−SBL algorithm, $N_f$ is the iteration number of the M−FCSBL algorithm and $N_g$ is the iteration number of MCV for our method. Due to the intractable inverse operations and the large number of iterations, the computational loads of M−SBL and M−FCSBL, as well as MCV are about three times of the DOF. By contrast, since the time−consuming training process of neural networks is finished before testing, BAMCV achieves fast super−resolution by directly mapping out sparse representation vectors of training samples with only twice of the DOF. To provide an explicit illustration, we also display the curve of computational loads versus DOF in Figure 10b with setting $K=16MN$, $L=10$ and $N_s=N_f=N_g=1$. Obviously, the computational complexity of BAMCV is much lower than that of other methods, and the gap between BAMCV and other methods diverges as DOF increasing.

In summary, the modified method BAMCV−STAP not only requires less online running time and lower computational loads than other comparison methods while achieving comparable clutter suppression performance, but also remains robust to the number of training samples and noise power. Consequently, BAMCV−STAP is suitable for airborne early warning (AEW) radar applications.

4.2. Measured Data

In this subsection, MCV−STAP and BAMCV−STAP are applied to the public available data set, the Multi−Channel Airborne Radar Measurements (MCARM) data set [61], to verify their attractiveness. Some indispensable parameters of the radar system are listed in

Table 7. Radar System Parameters of MCARM data.

Parameters	Value	Parameters	Value
Pulse repetition frequency (Hz)	1984	Antenna array spacing of azimuth (m)	0.1029
Wavelength (m)	0.24	Antenna array spacing of elevation (m)	0.5629
Pulse number in one CPI	128	Platform height (m)	10,188
Antenna elements number of azimuth	11	Range cell number	400
Antenna elements number of elevation	2

. The clutter spectrum estimated using 10 training samples around the selected target range cell by different methods is exhibited in Figure 11. Consistent with simulated results, limited by the number of training samples, the clutter spectrum estimated by the LSMI method in Figure 11a is as rough as the one by the M−OMP method in Figure 11b. Due to the approximation term, the clutter spectrum estimated by M−FCSBL in Figure 11d method is too sparse to reconstruct the information in the sidelobe region. M−SBL exhibits a better capability of recovering the clutter spectrum using only 10 training samples but with a slight broadening of the clutter ridge as shown in Figure 11c.

As it is difficult to judge whether our methods or M−SBL achieve the best performance according to the clutter spectrum estimated shown in Figure 11c,e,f, we explore the detection performance averaged over all range cells for the targets with the normalized Doppler frequency of 0.1 (corresponding velocity is 10 m/s) by the PD versus SNR curves as depicted in Figure 12a. The probability of the false alarm rate (PFA) is set as ${10}^{-6}$, the PD to SNR for each CUT is averaged over ${10}^5$ Monte Carlo trials. To further explore the detection performance of targets with different velocity values, we exhibit the PD versus velocity curves (SNR = 20 dB), which are averaged over ${10}^5$ Monte Carlo trials, in Figure 12b.

Moreover, the detection performance to explore the robustness to the number of training samples is provided in Figure 12c. As shown in Figure 12, similar to the LSMI method, M−FCSBL suffers a severe performance degradation as the information in the sidelobe region is ultimately given up due to the approximation term, leading to the too sparse reconstructions to obtain the accurate CCM. M−SBL method achieves sub−optimal performance among comparison methods but with very slow convergence, as shown in Figure 13 and

Table 8. Comparison of running time per iteration using MCARM data.

Approach	Running Time Per Iteration (s)
M−SBL	$1.102$
M−FCSBL	$0.333$
MCV	$0.951$
BAMCV	$0.025$

. On the contrary, BAMCV−STAP achieves the sub−optimal target detection performance. Besides that, according to the computational loads analyzed in Figure 10b and Table 6, the computational loads of BAMCV is about $\mathcal{O}\left(MN\right)$ but the other methods are about $\mathcal{O}({\left(MN\right)}^3)$. Consequently, BAMCV−STAP can achieve favorable detection of targets and keep robust to the number of training samples, with the most efficient convergence and the lowest computational loads among comparison methods.

To sum up, among these comparison approaches, our method exhibits its superiority in recovering the clutter spectrum in terms of both sparsity and accuracy on measured data, which is consistent with the results on simulated data.

5. Discussion

Along with the impressive performance achieved by our methods, the above−mentioned experimental results indicate that it is substantially more challenging for traditional SR−STAP algorithms, such as the proposed MCV−STAP, to ensure both the stable performance of sparse reconstruction and real−time processing. However, deep learning techniques, such as the proposed BAMCV−STAP, can efficiently realize parameter estimates for real−time on−line processing at the expense of nonsignificant off−line pre−training.

6. Conclusions

Traditional SR−STAP algorithms are unattractive as they only provide point estimates for several parameters and rely on a large number of iterations per testing at the testing stage. To avoid the iterative procedure, we first propose MCV−STAP, a parameter−free SR−STAP method that incorporates a modified VRVM into the field of STAP, to model the space−time echoes with a completely robust Bayesian framework. Then, to realize the real−time processing of the incoming space−time echoes, we introduce BAMCV−STAP, a VAE−structured method that builds its encoder from MCV but with an efficient inference network instead of the time−consuming Gibbs sampler in MCV−STAP. Compared with the conventional Bayesian type SR−STAP algorithms, the inference network of BAMCV−STAP can be pre−trained off−line using a gradient ascent scheme, enabling more efficient sparse reconstruction. The experimental results show that the BAMCV−STAP algorithm can achieve comparable clutter suppression with lower computational cost.

This paper focuses on the suppression of non−homogeneous clutter environments with side−looking airborne phased−array radar by means of deep learning methods. Inspired by their attractive characteristics, approaches that can solve the other tricky problems of radar signal processing, such as the suppression of the non−stationary clutter and the detection of slow−moving targets, are the promising directions that we plan to explore in the future.