Fast Variational Bayesian Inference for Space-Time Adaptive Processing

Abstract

Space-time adaptive processing (STAP) approaches based on sparse Bayesian learning (SBL) have attracted much attention for the benefit of reducing the training samples requirement and accurately recovering sparse signals. However, it has the problem of a heavy computational burden and slow convergence speed. To improve the convergence speed, the variational Bayesian inference (VBI) is introduced to STAP in this paper. Moreover, to improve computing efficiency, a fast iterative algorithm is derived. By constructing a new atoms selection rule, the dimension of the matrix inverse problem can be substantially reduced. Experiments conducted on the simulated data and measured data verify that the proposed algorithm has excellent clutter suppression and target detection performance.

1. Introduction

Space-time adaptive processing (STAP) is widely studied and applied in airborne phased array radar for clutter suppression and target detection [1,2,3,4,5]. The performance of STAP is related to the clutter plus noise covariance matrix (CNCM), which is estimated by the training samples adjacent to the cell under test (CUT). More than two times the system degrees of freedom (DOFs) of independent and identically distributed (IID) training samples are required to estimate the CNCM [6,7]. However, it is difficult to obtain enough training samples in a practical environment.

Various algorithms have been proposed to improve the STAP performance with a small number of training samples. Reduced-rank (RR) approaches make use of the structural characteristics of the covariance matrix to reduce the requirement for samples, such as principal components (PC), the multistage winner filter (MWF), and cross-spectral metric (CSM) [8,9,10,11]. Reduced-dimension (RD) approaches reduce the requirement for samples by selecting space-time or angle-Doppler channels, such as joint-domain localization (JDL), the extend factor approach (EFA), and auxiliary channel processor (ACP) [12,13,14,15]. However, the performance of the RR or RD approaches is affected due to complex practical environments. Knowledge-aided (KA) [16,17,18] approaches use prior knowledge to improve clutter suppression performance. Nevertheless, it is hard to obtain precise prior knowledge. Additionally, some covariance estimation algorithms in limited samples have been proposed, such as the fast maximum likelihood (ML) estimation [19], rank-constrained ML estimation [20], knowledge-aided ML estimation [21], and distribution-free covariance estimation [22]. These approaches are sensitive to model mismatch.

Motivated by the compressed sensing (CS) theory, sparse recovery STAP (SR-STAP) methods have attracted much attention over the last few years [23,24,25,26]. Using the norm minimization method and the sparsity of the clutter spectrum, SR-STAP approaches can improve the clutter suppression performance with few samples or even a single sample. Currently, several SR methods have been proposed. In [27], a convex optimization algorithm named the focal underdetermined system solver (FOCUSS) was proposed, which uses the ${\mathsf{\ell }}_p$ norm penalty to relax the cost function. This kind of algorithm requires a suitable regularization parameter. However, the parameter chosen empirically is not robust to environmental change. A greedy pursuit approach, an orthogonal matching pursuit, is presented in [28]. It performs well when the signal sparsity is known. The difficulty in obtaining signal sparsity in practice limits the application of this type of approach. For its self-regularizing nature, sparse Bayesian learning (SBL) has received much attention and has been applied to STAP [29,30,31]. However, the problems of heavy computational burden and slow convergence speed prevent its practical application.

Variational Bayesian inference (VBI) can be seen as a Bayesian approximate method that deals with the SR problem by maximizing a lower bound of likelihood function [32,33,34]. Compared with SBL, VBI replaces the point estimation of parameters with approximate posterior distribution estimates, which gives it a fast convergence speed. The VBI has been widely used in the fields of engineering and signal processing, some of which include model degradation analysis [35], statistical properties estimation [36], and Orthogonal Time Frequency Space (OTFS) modulation [37]. In this paper, the VBI algorithm with multiple measurement vectors (MMV) is introduced in STAP (M-VBI-STAP). Because of the matrix inverse operation, the computational burden is still heavy. To reduce the computational complexity, a fast STAP algorithm based on VBI is proposed. The following are the main contributions of this paper:

(1)

The frame of VBI with multiple measurement vectors is generalized to the STAP application. This allows for an increase in convergence speed.

(2)

In order to improve computational efficiency, a fast STAP algorithm based on VBI (M-FVBI-STAP) is proposed. A new atoms selection rule is constructed, which is able to substantially reduce the dimension of the matrix inverse problem.

(3)

Numerous experiment results on simulated data and measured data illustrate that the proposed algorithm can achieve better clutter suppression and target detection performance.

The remainder of this paper is outlined as follows. The signal model and Bayesian model are introduced in Section 2. In Section 3, the VBI algorithm is extended to STAP and proposes a fast STAP method based on VBI. Computational complexity analysis is detailed in Section 4. Section 5 shows the experimental results. The conclusions are presented in Section 6.

Notation: Matrices, vectors, and scalar quantities are denoted by boldface uppercase, boldface lowercase, and lightface letters, respectively. The symbols $ℂ$ and $I$ represent the complex filed and the identity matrix. For a matrix $X$, ${\left(X\right)}^{-1}$, ${\left(X\right)}^T$, and ${\left(X\right)}^H$ denote the inverse, transpose and conjugate transpose of $X$, respectively. $\Vert \bullet {\Vert }_2$, $\Vert \bullet {\Vert }_{2,0}$, and $\Vert \bullet {\Vert }_F$ denote the ${\mathsf{\ell }}_2$, ${\mathsf{\ell }}_{2,0}$ and Frobenius norms, respectively.$diag\left(x\right)$ denotes a diagonal matrix whose diagonal element is $x$. $\otimes$ stands for the kronecker product. $\langle \bullet \rangle$ denotes the expectation.

2. Signal Model and Bayesian Model

2.1. Signal Model

The system under consideration is a side-looking airborne-phased array radar with a uniform linear array (ULA), which consists of N elements spaced at half a wavelength. The radar transmits M pulses at a constant pulse repetition frequency (PRF) in a coherent processing interval (CPI). Ignoring the range ambiguity, the received space-time signal is modeled as

$$\begin{array}{cc}x& =x_c+x_t+n\\ & ={\displaystyle \sum _{i=1}^{N_c}{\alpha }_{i}s\left(f_{d,i},f_{s,i}\right)+{\alpha }_t{a}_t+n}\end{array}$$

(1)

where $x_c$ is the clutter component; $N_c$ denotes the number of clutter patches for each range cell; $x_t$ stands for the target component; $n$ denotes the complex white Gaussian noise; $f_{d,i}$ and $f_{s,i}$ are the normalized Doppler frequency and normalized spatial frequency; ${\alpha }_t$ and ${\alpha }_i$ are the complex amplitude of the target and the i-th clutter patch, respectively; $a_t$ represents the space-time steering vector of the target and $s\left(f_{d,i},f_{s,i}\right)$ is the space-time steering vector of the clutter with the form

$$s\left(f_{d,i},f_{s,i}\right)=s_t\left(f_{d,i}\right)\otimes s_s\left(f_{s,i}\right)$$

(2)

where $s_t\left(f_{d,i}\right)$ and $s_s\left(f_{s,i}\right)$ are the temporal steering vector and spatial steering vector

$$s_t\left(f_{d,i}\right)={[1,\exp \left(j2\pi f_{d,i}\right),\cdots ,\exp \left(j2\pi \left(M-1\right)f_{d,i}\right)]}^T$$

(3)

$$s_s\left(f_{s,i}\right)={[1,\exp \left(j2\pi f_{s,i}\right),\cdots ,\exp \left(j2\pi \left(N-1\right)f_{s,i}\right)]}^T$$

(4)

The optimum weight of the STAP filter is obtained from the problem,

$$\{\begin{array}{c}\underset{w}{\min }\begin{array}{cc}{w}^H{R}_{c+n}w& \end{array}\\ s.t.\begin{array}{cc}{}^{}{w}^H{a}_t=1& \end{array}\end{array}.$$

(5)

The optimum weight of the STAP filter can be acquired by

$$w_{opt}=\frac{R_{c+n}^{-1}{a}_t}{a_t^H{R}_{c+n}^{-1}{a}_t}.$$

(6)

This assumes that the clutter and noise components are mutually uncorrelated. Then, the clutter plus noise covariance matrix (CNCM) can be written as

$$\begin{array}{cc}{R}_{c+n}& =R_c+R_n\\ & =\langle x{x}^H\rangle +{\sigma }_{}^{2}I\end{array}$$

(7)

where $R_c\in {ℂ}^{MN\times MN}$ and $R_n\in {ℂ}^{MN\times MN}$ represent the CCM and noise covariance matrix and ${\sigma }_{}^2$ is the power of the noise.

For the SR-STAP algorithms, the angle-Doppler plane is uniformly discretized into $K=N_s{N}_d$ grids. Here, $N_s={\rho }_{s}N$ and $N_d={\rho }_{d}M$ are the number of angle bins, and Doppler bins, respectively. The received signal $X\in {ℂ}^{NM\times L}$ from L IID training samples can be written as

$$X=\Phi A+N$$

(8)

where $\Phi =[s_1,s_2\cdots ,s_K]\in {ℂ}^{MN\times K}$ is the dictionary matrix, $A=[{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_L]\in {ℂ}^{K\times L}$ denotes the sparse coefficient matrix, and $N=[n_1,n_2,\cdots ,n_L]\in {ℂ}^{MN\times L}$ represents the noise matrix.

The angle-Doppler profile can be estimated by solving the following problem:

$$A=\arg \underset{A}{\min }{\Vert A\Vert }_{2,0}\begin{array}{cc}& s.t{.}_{}^{}\end{array}{\Vert X-\Phi A\Vert }_F^2\le \epsilon$$

(9)

where $\epsilon$ is the noise error tolerance. Thus, the CNCM can be estimated by

$${\tilde{R}}_{c+n}=\frac{1}{L}{\displaystyle \sum _{l=1}^L\Phi diag({\alpha }_l^2)}{\Phi }^H+{\sigma }_{}^2{I}_{MN}.$$

(10)

2.2. Bayesian Model

From a Bayesian perspective, noise is assumed to be a zero-mean complex of Gaussian distribution. Then, the complex Gaussian likelihood for the $X$ can be modeled as

$$p\left(X|A,{\sigma }^2\right)={\displaystyle \prod _{l=1}^L\mathcal{C}\mathcal{N}\left(x_l|\Phi {\alpha }_l,{\delta }^{-1}{I}_{MN}\right)}$$

(11)

where $\delta ={\sigma }^{-2}$ is the noise precision. According to the conjugate prior principle [38], the $\delta$ follows a Gamma distribution,

$$p\left(\delta \right)=Gamma\left(\delta |c,d\right)$$

(12)

where $Gamma\left(\delta |c,d\right)=\frac{d^c}{\gamma \left(c\right)}{\delta }^{c-1}\exp \left(-d\delta \right)$, $\gamma \left(c\right)$ is the Gamma function and $c$ and $d$ represent the shape parameter and scale parameter.

To model the sparse coefficient, a complex Gaussian prior is employed first for each column of $A$

$$p\left({\alpha }_l|\lambda \right)={\displaystyle \prod _{k=1}^K\mathcal{C}\mathcal{N}\left(a_{k,l}|0,{\lambda }_k^{-1}\right)}$$

(13)

where $\lambda ={\left[{\lambda }_1,\cdots ,{\lambda }_k,\cdots ,{\lambda }_K\right]}^T$. A Gamma prior is assigned on the hyper-parameter ${\lambda }_k$

$$p\left({\lambda }_k\right)=Gamma\left({\lambda }_k|a,b\right)$$

(14)

In a word, the Bayesian model of the received signal has been established. Figure 1 shows the corresponding graphical model. Each node represents a conditional probability density.

3. Proposed Algorithm

In this section, a VBI-based STAP algorithm is derived for the MMV case. Then, a fast iterative structure is proposed to improve the computational efficiency of M-VBI-STAP.

3.1. VBI Algorithm

The set of hiding variables for the Bayesian model denotes $\mathsf{\Omega }=\left[A,\lambda ,\delta \right]$. The marginal likelihood function $p\left(X\right)$ needs to be calculated to obtain the closed form of the posterior distribution $p\left(\mathsf{\Omega }|X\right)$. Since $p\left(X\right)$ is intractable, it can be decomposed as

$$\ln p\left(X\right)=\mathcal{L}\left(q\right)+\mathcal{K}\mathcal{L}\left(q|p\right)$$

(15)

with

$$\mathcal{L}\left(q\right)={\displaystyle \int q\left(\mathsf{\Omega }\right)}\ln \left(\frac{p\left(X,\mathsf{\Omega }\right)}{q\left(\mathsf{\Omega }\right)}\right)d\mathsf{\Omega }$$

(16)

and

$$\mathcal{K}\mathcal{L}\left(q|p\right)=-{\displaystyle \int q\left(\mathsf{\Omega }\right)}\ln \left(\frac{p\left(\mathsf{\Omega }|X\right)}{q\left(\mathsf{\Omega }\right)}\right)d\mathsf{\Omega }$$

(17)

where $q\left(\mathsf{\Omega }\right)$ denotes a variational distribution, $\mathcal{L}\left(q\right)$ is the lower bound of $\ln p\left(X\right)$, and $\mathcal{K}\mathcal{L}\left(q|p\right)$ represents the Kullback–Leibler divergence between the $p\left(\mathsf{\Omega }|X\right)$ and $q\left(\mathsf{\Omega }\right)$ regularization parameters. The VBI algorithm approximates the posterior distribution using the factorized variable distribution. Assuming the variables are independent of each other, the factorized variable distribution can be expressed as

$$q\left(\mathsf{\Omega }\right)=q\left(X\right)q\left(\lambda \right)q\left(\delta \right).$$

(18)

Maximizing $\mathcal{L}\left(q\right)$, the optimal factorized variable distribution is

$$\ln q^{\ast }\left(\mathsf{\Omega }\right)={\langle \ln p\left(X,\mathsf{\Omega }\right)\rangle }_{q\left(\mathsf{\Omega }\backslash {\mathsf{\Omega }}_i\right)}+\mathrm{const}$$

(19)

where ${\mathsf{\Omega }\backslash \mathsf{\Omega }}_i$ means removing ${\mathsf{\Omega }}_i$ from $\mathsf{\Omega }$, and

$$p\left(X,\mathsf{\Omega }\right)=p\left(X|A,\lambda ,\delta \right)p\left(A|\lambda \right)p\left(\lambda \right)p\left(\delta \right)$$

(20)

Combining Equations (19) and (20), the following gives the iterative update procedure of the variational parameters:

(1)

update $A$

$$\begin{array}{cc}\ln q^{\ast }\left(A\right)& ={\langle \ln p\left(X|A,\lambda ,\delta \right)+\ln p\left(A|\lambda \right)\rangle }_{q\left(\mathsf{\Omega }\backslash X\right)}+\mathrm{const}\\ & =-\delta {\displaystyle \sum _{l=1}^L{\Vert x_l-\Phi {\alpha }_l\Vert }_2^2}+{\displaystyle \sum _{l=1}^L{\alpha }_l^{H}diag\left(\lambda \right)}{\alpha }_l+\mathrm{const}\end{array}$$

(21)

It can be found that $q^{\ast }\left(A\right)$ obeys the Gaussian distribution

$$q^{\ast }\left(A\right)={\displaystyle \prod _{l=1}^L\mathcal{C}\mathcal{N}\left({\alpha }_l|{\mu }_l,\Sigma \right)}$$

(22)

where

$${\mu }_l=\delta \Sigma {\Phi }^H{x}_l$$

(23)

$$\Sigma =diag\left(\lambda \right)-diag\left(\lambda \right){\Phi }^H{\left(\delta I_{MN}+\Phi diag\left(\lambda \right){\Phi }^H\right)}^{-1}\Phi diag\left(\lambda \right)$$

(24)

(2)

Update $\lambda$

$$\begin{array}{cc}\ln q^{\ast }\left(\lambda \right)& ={\langle \ln p\left(A|\lambda \right)+\ln p\left(\lambda \right)\rangle }_{q\left(\mathsf{\Omega }\backslash \lambda \right)}+\mathrm{const}\\ & ={\displaystyle \sum _{k=1}^K\left[\left(a+L-1\right)\ln {\lambda }_k-\left(b+{\displaystyle \sum _{l=1}^L\langle a_{k,l}{a}_{k,l}^H\rangle }\right){\lambda }_k\right]}\end{array}$$

(25)

From Equation (25), $q^{\ast }\left(\lambda \right)$ follows Gamma distribution and ${\lambda }_k$ can be obtained by

$${\lambda }_k=\frac{a+L}{b+{\displaystyle \sum _{l=1}^L\langle a_{k,l}{a}_{k,l}^H\rangle }}$$

(26)

where $\langle a_{k,l}{a}_{k,l}^H\rangle ={\mu }_{k,l}^2+{\Sigma }_{k,k}$, ${\mu }_{k,l}$ is the k-th element of ${\mu }_l$ and ${\Sigma }_{k,k}$ denotes the k-th main diagonal element of $\Sigma$.

(3)

Update $\delta$:

$$\begin{array}{cc}\ln q^{\ast }\left(\delta \right)& ={\langle \ln \ln p\left(X|A,\lambda ,\delta \right)+\ln p\left(\delta \right)\rangle }_{q\left(\mathsf{\Omega }\backslash \delta \right)}+\mathrm{const}\\ & =\left(MNL+c-1\right)\ln \delta -\left(\langle {\displaystyle \sum _{l=1}^L{\Vert x_l-\Phi {\alpha }_l\Vert }_2^2}\rangle +d\right)\delta \end{array}$$

(27)

Then, the result of $\delta$ is

$$\delta =\frac{MNL+c}{{\Vert X-\Phi \mu \Vert }_F^2+trace\left(\Sigma {\Phi }^H\Phi \right)+d}$$

(28)

When the maximum number of iterations or the predefined stop condition is reached, the iterative procedure is terminated. Finally, the CNCM can be estimated by

$${\tilde{R}}_{c+n}=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{k=1}^K{\Vert {\mu }_{k,l}\Vert }^2{s}_k}}{s}_k^H+{\sigma }_{}^2{I}_{MN}.$$

(29)

For clarification, the steps of the M-VBI-STAP algorithm are described in Algorithm 1.

Algorithm 1: Pseudocode of the M-VBI.

step 1: Input: data $X$ and dictionary $\Phi$; step 2: Initialize: $a=b=c=d={10}^{-6}$, ${\lambda }^{\left(0\right)}=$ 1; step 3: While if it does not converge:    1. Calculate ${\mu }_l$, and $\Sigma$ by Equations (23) and (24);    2. Calculate $\lambda$ by Equation (26);    3. Calculate $\delta$ by Equation (28);    End step 4: Obtain CNCM ${\tilde{R}}_{c+n}$ by Equation (29); step 5: Output: the STAP weight.

3.2. Variational Fast Solution

In a previous subsection, the updated criteria for the variational parameters were derived. However, the real-time performance of the algorithm is also crucial due to the matrix inverse operation in Equation (24). In addition, the overcomplete dictionaries used in the traditional SR-STAP algorithms also affect the computational efficiency. Here, a fast algorithm to reduce the computational complexity is introduced by analyzing the behavior of the distribution $q^{\ast }\left(\lambda \right)$. The high dimensional matrix inverse operation and a full dictionary are not needed in the fast algorithm.

The ${\mu }_{k,l}^2+{\Sigma }_{k,k}$ in Equation (26) can be rewritten as $e_k^H\left({\mu }_l{\mu }_l^H+\Sigma \right)e_k$, where $e_k$ is the k-th column of the unit matrix with the dimension K. By noting ${\mu }_l{\mu }_l^H={\delta }^2\Sigma {\Phi }^H{x}_l{x}_l^H\Phi {\Sigma }^H$ and $diag\left(\lambda \right)={\displaystyle \sum _{k=1}^K{\lambda }_k}{e}_k{e}_k^H$, Equation (24) is equivalent to

$$\begin{array}{cc}\Sigma & ={\left(\delta {\Phi }^H\Phi +{\displaystyle \sum _{j\ne k}{\lambda }_k}{e}_k{e}_k^H+{\lambda }_j{e}_j{e}_j^H\right)}^{-1}\\ & ={\overline{\Sigma }}_j-\frac{{\overline{\Sigma }}_j{e}_j{e}_j^H{\overline{\Sigma }}_j}{{\lambda }_j^{-1}+e_j^H{\overline{\Sigma }}_j{e}_j}\end{array}$$

(30)

where the Woodbury Matrix Identity is used in the second expression, and ${\overline{\Sigma }}_j={\left(\delta {\Phi }^H\Phi +{\displaystyle \sum _{j\ne k}{\lambda }_k}{e}_k{e}_k^H\right)}^{-1}$.

Defining ${\varsigma }_j=e_j^H{\overline{\Sigma }}_j{e}_j$ and ${\rho }_{j,l}=\delta e_j^H{\overline{\Sigma }}_j{\Phi }^H{x}_l$, the modified version of Equation (26) can be obtained

$${\lambda }_j^{\left(t+1\right)}=\frac{\left(a+L\right)\left(1+{\lambda }_j^{\left(t\right)}{\varsigma }_j\right)}{b+{\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2+{\varsigma }_j+{\lambda }_j^{\left(t\right)}{\varsigma }_j^2\right)}}$$

(31)

where t is the number of iterations. For brevity, when setting $a=b=0$, Equation (31) can be seen as nonlinear maps ${\lambda }_j^{\left(t+1\right)}=\Gamma \left({\lambda }_j^{\left(t\right)}\right)$. Calculating the stationary points of this map [39], the optimal ${\lambda }_j^{\ast }$ is

$${\lambda }_j^{\ast }=\{\begin{array}{cc}L{\left[{\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}\right]}^{-1},& {\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}>0\\ \infty ,& {\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}\le 0\end{array}$$

(32)

This means that:

(1)

If $s_j$ is excluded from the model and ${\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}>0$, add $s_j$ to the model;

(2)

If $s_j$ is in the model and ${\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}>0$, re-estimate ${\lambda }_j$;

(3)

If $s_j$ is in the model and ${\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}\le 0$, remove $s_j$ from the model.

To efficiently compute ${\varsigma }_j$ and ${\rho }_{j,l}$, the dictionary and sparse coefficient are set as ${\Psi }^{\left(0\right)}\in \varnothing$, ${\alpha }^{\left(0\right)}=0$. Thus, the ${\overline{\Sigma }}_j$ can be written as

$${\overline{\Sigma }}_j={\left(\begin{array}{cc}\delta {\Psi }^H\Psi +diag(\alpha )& \delta {\Psi }^H{s}_j\\ \delta s_j^H\Psi & \delta s_j^H{s}_j\end{array}\right)}^{-1}$$

(33)

Using the block matrix inversion, Equation (35) is equivalent to

$${\overline{\Sigma }}_j=\left(\begin{array}{cc}\Sigma +{\gamma }_j{\delta }^2\Sigma {\Psi }^H{s}_j{s}_j^H\Psi \Sigma & -{\gamma }_j\delta \Sigma {\Psi }^H{s}_j\\ -{\gamma }_j\delta s_j^H\Psi \Sigma & {\gamma }_j\end{array}\right)$$

(34)

where ${\gamma }_j={\left(\delta s_j^H{s}_j-{\delta }^2{s}_j^H\Psi \Sigma {\Psi }^H{s}_j\right)}^{-1}$, $\Sigma ={\left(\delta {\Psi }^H\Psi +diag\left(\alpha \right)\right)}^{-1}$ is the variance and the mean is ${\mu }_l=\delta \Sigma {\Psi }^H{x}_l$.

Then, ${\varsigma }_j$ and ${\rho }_{j,l}$ can be updated as

$${\varsigma }_j={\left(\delta s_j^H{s}_j-{\delta }^2{s}_j^H\Psi \Sigma {\Psi }^H{s}_j\right)}^{-1}$$

(35)

$${\rho }_{j,l}=\delta {\varsigma }_j{s}_j^H{x}_l-{\delta }^2{\varsigma }_j{s}_j^H\Psi \Sigma {\Psi }^H{x}_l$$

(36)

The computation complexity of sequential selecting all dictionary atoms is still large. Inspired by the marginal likelihood maximization algorithm [40], a function $f\left({\lambda }_j\right)$ can be constructed which has a unique maximum value for ${\lambda }_j$. Combining Equations (25) and (31), we have

$$f\left({\lambda }_j\right)=L\ln {\lambda }_j-{\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}{\lambda }_k$$

(37)

In the $\left(t+1\right)$ iteration, we maximize the incremental $f\left({\lambda }_j\right)$ to select a candidate atom. The serial number $z$ is obtained by

$$z=\underset{j,1\le j\le N_s{N}_d}{\arg \max }\left[f\left({\left({\lambda }_j^{\ast }\right)}^{\left(t+1\right)}\right)-f\left({\left({\lambda }_j^{\ast }\right)}^{\left(t\right)}\right)\right]$$

(38)

Utilizing Equation (38), the computational complexity decreases because the number of the selected atoms is smaller than the size of dictionary. In summary, an iterative process to update parameter ${\lambda }_j$ can be obtained. The update formulas are detailed in Appendix A. The pseudocode of the proposed M-F-VBI algorithm is described in Algorithm 2.

Algorithm 2: Pseudocode of the M-F-VBI.

step 1: Input: data $X$ and dictionary $\Phi$; step 2: Initialize: ${\lambda }^{\left(0\right)}=$ 0, ${\Psi }^{\left(0\right)}=\varnothing$, and $E^0=I$; step 3: While, if it does not converge:    1. Calculate all ${\lambda }_j^{\ast }$, $\forall j$ by (31) and choose the atomic index by (38);    2. If $s_j$ is in the model and ${\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}\le 0$, use (44)–(47) to update parameters; If $s_j$ is in the model and ${\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}>0$, use (48)–(51) to update parameters; If $s_j$ is excluded from the model and ${\displaystyle \sum _{l=1}^L\left({\rho }_{j,l}^2-{\varsigma }_j\right)}>0$, use (53)–(56) to update parameters.    3. Update $\delta$ by (28).    End step 4: Obtain CNCM ${\tilde{R}}_{c+n}$ by (29); step 5: Output: the STAP weight.

4. Computational Complexity Analysis

The computational complexity of the proposed algorithm is investigated via the number of complex multiplications, and the lower-order terms of the multiplication are omitted. In the proposed M-FVBI-STAP algorithm, the computational complexity is related to the dimension of $\Psi \in {ℂ}^{MN\times P}$. $P$ is the number of selected atoms that is assumed to be twice the rank of the clutter. The most computational part of M-VBI-STAP is Equation (24).

Table 1. Computational complexity.

Algorithm	Computational Load
M-SBL	$o\left(\left({\left(MN\right)}^3+2K{\left(MN\right)}^2+K^{2}MN+KMNL\right)T_{SBL}\right)$
M-FOCUSS	$o\left(\left({\left(MN\right)}^3+2{\left(MN\right)}^{2}K+MNLK\right)T_{FOCUSS}\right)$
M-VBI	$o\left(\left({\left(MN\right)}^3+2K{\left(MN\right)}^2+K^{2}MN+KMNL\right)T_{VBI}\right)$
M-FVBI	$o\left(\left(MNK+MNL+4MNKP+MNLP\right)T_{FVBI}\right)$

shows the computational complexity of different approaches. $T_{SBL}$, $T_{FOCUSS}$, $T_{VBI}$ and $T_{FVBI}$ represent the number of iterations of M-FOCUSS-STAP, M-SBL-STAP [24], M-VBI-STAP and M-FVBI-STAP. The results are plotted in Figure 2. As the number of system DOFs increases, it is observed that the computational complexity of the M-VBI-STAP is lower than M-SBL-STAP and M-FOCUSS-STAP. Compared to M-VBI-STAP, M-FVBI-STAP can improve computational efficiency.

5. Numerical Experiments

5.1. Simulated Data

In this section, numerical experiments are conducted to illustrate the STAP and target detection performance of the proposed algorithm on simulated data. The main parameters of the simulations are listed in

Table 2. Parameters of radar system.

Parameter	Value	Unit
Number of elements	8	/
Number of pulses	8	/
Wavelength	0.3	m
Bandwidth	2.5	MHz
Height of platform	150	m/s
Velocity of platform	9000	m
Pulse repetition frequency (PRF)	2000	Hz
Clutter to noise ratio (CNR)	40	dB

and ${\rho }_s={\rho }_d=5$. For comparison, M-SBL-STAP and M-FOCUSS-STAP are presented. In addition, the optimal STAP with known CNCM and diagonal loading SMI-STAP (LSMI-STAP) are also included. The results are averaged over 100 Monte Carlo trials, and the number of training samples used in all experiments is ten. To measure the performance of clutter suppression, the improvement factor (IF) and signal-to-interference plus noise ratio (SINR) loss are used, which are denoted as

$$IF=\frac{\left|w^H{a}_t\right|}{w^H{R}_{c+n}w}\frac{tr\left(R_{c+n}\right)}{a_t^H{a}_t}$$

(39)

$$L_{SINR}=\frac{{\sigma }^2}{MN}\frac{{\left|w^H{a}_t\right|}^2}{w^H{R}_{c+n}w}$$

(40)

The STAP performance in terms of Capon spectra is first investigated. Figure 3a plots the optimal spectrum with known CNCM. From Figure 3b–f, it can be noted that the estimation results of the proposed M-FVBI, M-VBI algorithm, and M-SBL-STAP are close to the optimal spectrum. In addition, the spectra of M-FOCUSS-STAP have an extension in the main-lobe region. It can also be seen that the spectrum is inaccurate when estimated by LSMI-STAP due to the lack of training samples. Figure 4 plots the Capon spectra with gain–phase (GP) error, which are set as 0.05 and 2°. Compared to Figure 3a, it can be seen that the results of the proposed M-FVBI, M-VBI algorithm, and M-SBL-STAP are close to the optimal spectrum.

Then, the IF curves of different approaches against the normalized Doppler frequency in the main beam direction are presented in Figure 5. It is shown in Figure 5 that M-FVBI-STAP, M-VBI-STAP, and M-SBL-STAP have a similar performance, which is near the optimal performance. In the main-lobe area, the curves of the M-FOCUSS-STAP are broader than others.

For the next experiment, the average SINR loss curves versus the number of training samples are plotted in Figure 6, which were calculated as the mean SINR loss curves for $f_d=\left(-0.5,-0.1\right)\cup \left(0.1,0.5\right)$. It can be seen from Figure 6 that all the approaches achieved a relatively stable performance when the number of training samples was more than three, and the stable performance of M-FOCUSS-STAP was worse than others. The M-VBI-STAP and M-FVBI-STAP had the best performance.

For a clearer illustration, the target detection performance of the proposed algorithm was accessed. According to [41], the probability of detection (PD) versus the probability of false alarm (PFA) is defined as

$$PD=PF{A}^{\frac{1}{1+\xi }}$$

(41)

where $\xi$ is the output signal-to-clutter-plus-noise ratio (SCNR). We set the target located in the main beam direction.

As shown in Figure 7, the receiver operating characteristics (ROC) (i.e., PD versus PFA) are plotted with the target signal-to-noise ratio (SNR), 0 dB and −4 dB, respectively. The normalized Doppler frequency is $f_d=0.1$. It can be seen that the detection performances of the M-SBL-STAP, M-VBI-STAP, and M-FVBI-STAP are close to optimal. Another notable result is that the detection performances increased with the growth of the SNR value. The ROC results with $f_d=0.3$ are shown in Figure 8. The detection performance increased compared to $f_d=0.1$. In both cases, M-VBI-STAP and M-FVBI-STAP are considered optimal.

Finally, the average running time for different cases is summarized in

Table 3. The average running time (second).

Algorithms	Running Time
M-SBL	34.88 s
M-FOCUSS	11.16 s
M-VBI	3.70 s
M-FVBI	0.05 s

. The experiments were conducted by MATLAB 2017b on the computer equipped with Intel(R) Xeon(R) E5-2620 CPU @ 2.10Hz 2.10Hz and 64G RAM. It was found that the average running time of M-VBI-STAP was less than M-SBL-STAP. M-FVBI-STAP significantly accelerated the average running time.

5.2. Measured Data

In this subsection, the performance of the proposed algorithm is accessed with publicly available data from Multi-Channel Airborne Radar Measurement (MCARM) data. [42]. The airborne radar of the MCARM can be viewed as side looking. In CPI data, the number of coherent pluses, array elements, and range cells are 128, 22, and 630, respectively. For 12 pluses, data from the first space channels are selected for experimental analysis. There is a target added in the 320th range cell. Since the first 150 range cells cannot be used, the clutter power spectrum shown in Figure 9 was estimated using ten training samples selected from ye 151st range cell to the 630th range cell. It can be seen that all the approaches can recover the clutter ridge. The spectra of the M-FOCUSS-STAP have an extension in the main-lobe region. In addition, M-FVBI-STAP estimated CNCM by selecting the atoms. Thus, the clutter ridge of the M-FVBI-STAP is concentrated on a diagonal line.

The STAP output power from the 280th range cell to the 360th range cell is depicted in Figure 10. The generalized inner product (GIP) method [43] was used to select ten range cells of the twenty range cells around the CUT as the training samples. Additionally, four range cells on each side of the CUT were chosen as the guard cells. It can be observed that all algorithms clearly detected the target. To further illustrate the detection performance, the output signal to clutter plus noise ratio (SCNR) was calculated at the target range cell. The result corresponding to the algorithms is listed in

Table 4. The output SCNR (dB).

Algorithms	Results
M-SBL	16.43 dB
M-FOCUSS	16.89 dB
M-VBI	18.86 dB
M-FVBI	19.58 dB

. The proposed algorithm clearly had the best detection performance of all approaches.

6. Conclusions

In this paper, the framework of VBI into the STAP algorithm has been developed, which can improve the convergence speed compared to SBL. Due to the heavy computational burden caused by the matrix inverse operation, a fast iterative algorithm was derived. The numerical results of the simulated data and measured data demonstrate that the proposed algorithm has superior STAP and target detection performance and can effectively reduce computational complexity. VBI still belongs to SR algorithms, which require dictionary construction. When the grid mismatch exists, the performance of VBI decreases. The direction for future research is mitigating the grid mismatch effect on SR-STAP.