Adaptive Support-Driven Sparse Recovery STAP Method with Subspace Penalty

Abstract

Detecting a moving target is an attractive topic in many fields, such as remote sensing. Space-time adaptive processing (STAP) plays a key role in detecting moving targets in strong clutter backgrounds for airborne early warning radar systems. However, STAP suffers serious clutter suppression performance loss when the number of training samples is insufficient due to the inhomogeneous clutter environment. In this article, an efficient sparse recovery STAP algorithm is proposed. First, inspired by the relationship between multiple sparse Bayesian learning (M-SBL) and subspace-based hybrid greedy algorithms, a new optimization objective function based on a subspace penalty is established. Second, the closed-form solution of each minimization step is obtained through the alternating minimization algorithm, which can guarantee the convergence of the algorithm. Finally, a restart strategy is used to adaptively update the support, which reduces the computational complexity. Simulation results show that the proposed algorithm has excellent performance in clutter suppression, convergence speed and running time with insufficient training samples.

1. Introduction

An important task of an airborne early warning radar is to detect slow moving targets. However, since the airborne early warning radar is in the downward-looking working mode, it often faces ground clutter with strong energy, and the target is often submerged in the strong clutter. How to effectively suppress the ground clutter is the key problem to be solved in the airborne radar system. Space-time adaptive processing (STAP) has played a key role in airborne radar systems because of its superior clutter suppression performance [1]. The traditional STAP algorithm uses the training samples adjacent to the cell under test (CUT) to obtain the clutter pulse noise covariance matrix (CNCM) estimation through the maximum likelihood method. The well-known Reed–Mallet–Brennan (RMB) criterion [2] shows that the following conditions must be met if the performance loss is not to exceed 3 dB. (1) Training samples must be independent and identically distributed (i.i.d); (2) the number of training samples should not be less than twice of the system’s degree of freedom (DOF); (3) and there should be no target in the training sample. However, it is difficult to obtain enough i.i.d training samples in the actual environment due to the ground environment, antenna configuration and so on.

In order to reduce the demand for the i.i.d training samples, many effective algorithms have been proposed in the past few decades. Reduced-dimension (RD) [3,4] and reduced-rank (RR) [5,6,7] algorithms can reduce the requirements of i.i.d training samples to twice of the reduced-dimension and twice of the clutter rank, respectively. However, due to the huge number of antenna elements and pulses, the number of i.i.d samples after RD or RR is still difficult to meet. In [8], an algorithm named direct data domain (D3) using only the CUT is proposed to solve the problem of insufficient training samples, but the D3 algorithm brings aperture loss. In recent years, knowledge-aided (KA) algorithms [9,10] have been proposed to improve the clutter suppression performance in inhomogeneous environments by using the prior information of the system and the surface. However, the performance of the KA algorithm is restricted by the accuracy of prior knowledge. Inaccurate prior knowledge will lead to performance degradation.

In recent years, with the rise of compressed sensing (CS) technology, many scholars have begun to focus on using sparse recovery (SR) technology to solve the problem of clutter suppression under the condition of insufficient samples. Many SR-based STAP (SR-STAP) algorithms have been proposed and achieved performance improvements. SR-STAP exploits the inherent sparsity of clutter in the angle-Doppler plane to achieve the accurate recovery of the clutter spectrum. In [11,12,13], some SR-STAP algorithms based on ${\ell }_1$ norm minimization are proposed. However, these algorithms with single measurement vectors (SMV) only utilize the information provided by a single sample, and do not consider the joint sparse recovery of multiple measurement vectors (MMV). In [14,15,16,17], many SR-STAP algorithms with MMV are proposed, which result in better recovery performance. Moreover, a hybrid greedy algorithm called subspace augmented MUSIC (SA-MUSIC) is proposed [18,19], which can achieve a better recovery performance. This algorithm combines greedy steps with generalized MUSIC criteria to obtain a complete subspace.

In addition to the above-mentioned ${\ell }_1$ norm-based and subspace-based greedy algorithms, there is also a class of sparse recovery algorithms with a great performance called sparse Bayesian learning (SBL). SBL does not need to set regularization parameters and is robust to noise, so it has attracted the attention of many scholars. SBL was proposed by Tipping in 2001 [20] and extended to the case of MMV in [21]. SBL was introduced into STAP (M-SBL) and it obtained a superior clutter suppression performance by Duan in [22]. However, M-SBL has large computational complexity and slow convergence speed, which limits its practical application. In [23], a fast-converging SBL (M-FCSBL) algorithm is proposed. However, the global convergence of the M-FCSBL algorithm is not guaranteed and it is sensitive to errors.

Although SBL differs greatly from the traditional sparse recovery algorithms in form, these two types of algorithms are intrinsically related. Wipf et al. pointed out that SBL can be transformed into the standard sparse restoration form by adding an additional non-separable sparsity inducing the prior, i.e., the log determinant penalty in [24]. In addition, the authors analyzed the benefits of the non-separable penalty; they found that SBL can recover true data from observations without noise. In [25], the authors reveal the connection between SBL and subspace-based hybrid greedy algorithms.

In this paper, we use the connection of SBL to subspace-based hybrid greedy algorithms to formulate a subspace penalty-based optimization problem. In order to ensure the convergence of the algorithm, we solve the above optimization problem by using the alternating minimization algorithm, and obtain the closed-form solution of the parameter update. In order to reduce the computational complexity, we adaptively update the support set in the iterative process which is similar to [26], thereby reducing the dimension of the dictionary and improving the computational efficiency. The major contributions of this paper are as follows:

• We formulate a new optimization problem based on a subspace penalty. This optimization problem considers the connection between SBL and subspace-based algorithms;
• To solve the above optimization problem, we use an alternating minimization algorithm. For each minimization step, a closed-form solution is obtained, which guarantees the convergence of the algorithm;
• The restart strategy is used to adaptively estimate the support set in each iteration, thereby reducing the computational complexity;
• The proposed algorithm (ASDSPL) is compared with other common SR-STAP algorithms in terms of the output signal to interference plus noise ratio (SINR) loss, improvement factor (IF), Capon spectrum, computational complexity and convergence speed using the simulation data.

The rest of this paper is organized as follows. In Section 2, the SR-STAP model is established. In addition, we also review some related knowledge of SBL. The proposed algorithm is presented in Section 3. In Section 4, the computational complexity analysis is evaluated. In Section 5, the performance of the proposed method is evaluated with simulation data. The conclusion is drawn in Section 6.

Notation: Throughout this paper, scalar quantities are denoted with italic typeface and lowercase bold and uppercase bold represent vectors and matrices, respectively. The $l-\mathrm{th}$ column of the matrix $A$ is represented as $a_l$. $a_{l,k}$ represents the $k-\mathrm{th}$ element of $a_l$. The matrix transpose and conjugate transpose are represented by ${\left(·\right)}^{\mathrm{T}}$ and ${\left(·\right)}^{\mathrm{H}}$, respectively. The symbol ${‖·‖}_{\mathrm{F}}$ and ${‖·‖}_0$ are reserved for Frobenius (${\ell }_{\mathrm{F}}$) norms and ${\ell }_0$ pseudo-norms. ${‖·‖}_{2,0}$ represents a mixed norm, which is defined as the non-zero element of the column vector consisting of the ${\ell }_2$ norm of each row vector. The determinant of a matrix is denoted by $\left|·\right|$. The expectation of a random variable is denoted by the symbol ${\rm E}\left(·\right)$. The symbol $\otimes$ is used to represent the Kronecker product. The symbol $CN(x|c,d)$ represents that the variable $x$ obeys the complex Gaussian distribution with the mean $c$ and variance $b$. The symbol $\mathrm{tr}(A)$ denotes the trace of the matrix $A$, the sum of the diagonal elements of $A$.

2. SR-STAP Model and SBL Review

2.1. SR-STAP Model

In this paper, an airborne radar system with a uniform line array (ULA) is considered. The antenna array consists of $N$ array elements with an array element spacing $d$ of half of the wavelength $\lambda$. The platform flies at a constant velocity $v_p$. $M$ pulses are transmitted at a constant pulse repetition frequency (PRF) $f_r$ within a coherent processing interval (CPI).

According to the Ward clutter model [1], the ground clutter is divided into many range rings, and each range ring is evenly divided into $N_c$ clutter patches, as shown in Figure 1. The clutter on each range ring can be viewed as a superposition of $N_c$ clutter patches. The spatial frequency and normalized Doppler frequency of the $i-\mathrm{th}$ clutter patch can be expressed as:

$$f_{s,i}=\frac{d}{\lambda }\cos ({\theta }_i)\sin ({\phi }_i)$$

(1)

$$f_{d,i}=\frac{2v_p}{\lambda f_r}\cos ({\theta }_i)\sin ({\phi }_i)$$

(2)

where ${\theta }_i$ and ${\phi }_i$ are the elevation angle and the azimuth angle corresponding to the $i-\mathrm{th}$ clutter patch, respectively.

Then the $MN\times 1$ space-time snapshot of the $l-\mathrm{th}$ range ring can be expressed as:

$$x_l=x_c+x_n={\displaystyle \sum _{a=1}^{N_a}{\displaystyle \sum _{i=1}^{N_c}{\alpha }_{a,i}s(f_{s,a,i},f_{d,a,i})}}+x_n$$

(3)

where $N_a$ indicates the number of ambiguous ranges, ${\alpha }_{a,i}$ denotes the random complex amplitude corresponding to the $i-\mathrm{th}$ clutter patch of the $a-\mathrm{th}$ ambiguous range and the $x_n$ is assumed to be zero-mean Gaussian white noise. The space-time steering vector of the $i-\mathrm{th}$ clutter patch $s(f_{s,i},f_{d,i})$ can be written as:

$$s(f_{s,i},f_{d,i})=b(f_{d,i})\otimes a(f_{s,i})$$

(4)

$$b(f_{d,i})={[1,\exp (j2\pi f_{d,i}),...,\exp (j(M-1)2\pi f_{d,i})]}^{\mathrm{T}}$$

(5)

$$a(f_{s,i})={[1,\exp (j2\pi f_{s,i}),...,\exp (j(N-1)2\pi f_{s,i})]}^{\mathrm{T}}$$

(6)

here $b(f_{d,i})$ and $a(f_{s,i})$ are the corresponding temporal and spatial steering vectors.

The ideal CNCM $R_{c+n}$ can be expressed as:

$$R_{c+n}=\mathrm{E}[x_l{x}_l{}^H]=R_c+{\sigma }^2{I}_{MN}$$

(7)

where $R_c$ is the covariance matrix of clutter and ${\sigma }^2$ denotes the noise power. The optimal STAP weight vector can be obtained by the linearly constrained minimum variance (LCMV) criterion:

$$w_{\mathrm{opt}}=\frac{R_{c+n}^{-1}{s}_{\mathrm{t}}}{s_{\mathrm{t}}^{\mathrm{H}}{R}_{c+n}^{-1}{s}_{\mathrm{t}}}$$

(8)

where $s_{\mathrm{t}}$ is the space-time steering vector pointing in the direction of the main beam.

In traditional SR-STAP algorithms, the angular-Doppler plane is uniformly discretized into $K=N_s{N}_d$ grid points, where $N_s={\rho }_{s}N$ (${\rho }_s>1$) represents the number of angle bins and $N_d={\rho }_{d}M$ (${\rho }_d>1$) represents the number of Doppler bins. Then the SR signal model with MMV can be written as:

$$X=DA+X_n$$

(9)

where $X=[x_1,x_2,...,x_L]$; $D$ is the space-time overcomplete dictionary consisting of the space-time steering vectors corresponding to all grid points; $A=[a_1,a_2,...,a_K]$ denotes the sparse coefficient matrix; $L$ denotes the number of snapshots; and $X_n$ denotes zero-mean Gaussian white noise.

Then the SR-STAP problem can be formulated in the form of basis pursuit denoising (BPDN):

$$\begin{array}{l}\tilde{A}=\arg \underset{A}{\min }\left|\right|A|{|}_{2,0}\\ \mathrm{s}.\mathrm{t}.\left|\right|X-DA|{|}_{\mathrm{F}}^2\le \epsilon \end{array}$$

(10)

where the ${\ell }_{2,0}$ norm is used to measure the sparsity and the ${\ell }_F$ norm is used to indicate the fitness of the data. $\epsilon$ denotes the fitting error tolerance.

2.2. SBL Review

In SBL framework, assume that the noise follows a complex Gaussian distribution with a mean of zero and the variance $\sigma$. Then the likelihood function of the observation data $X$ can be expressed as:

$$p\left(X|A,\sigma \right)={\left(\pi \sigma \right)}^{-MNL}\exp \left[-{\sigma }^{-1}{\displaystyle \sum _{l=1}^L{‖x_l-D{a}_l‖}^2}\right]$$

(11)

In the MMV case, it is assumed that multiple measurement vectors satisfy the i.i.d assumption and have the same sparse property. A Gaussian prior with zero-mean is assigned to each column of $A$

$$a_l~\mathcal{C}\mathcal{N}\left(0,\mathsf{\Gamma }\right),l=1,2,\cdots ,L$$

(12)

where $\mathsf{\Gamma }=\mathrm{diag}\left(\mathsf{\gamma }\right)$, $\mathsf{\gamma }=\left\{{\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_K\right\}$ represents the hyper-parameter vector.

According to the likelihood function (11) and the prior distribution (12), the posterior distribution can be expressed as:

$$\begin{array}{l}p\left(A|X,\mathsf{\Gamma },\sigma \right)=\frac{p\left(X,A|\mathsf{\Gamma },\sigma \right)}{p\left(X|\mathsf{\Gamma },\sigma \right)}\\ ={\pi }^{-KL}{\left|\mathsf{\Sigma }\right|}^{-L}\exp \left[-{\displaystyle \sum _{l=1}^L{\left(a_l-{\mu }_l\right)}^{\mathrm{H}}{\left(\mathsf{\Sigma }\right)}^{-1}\left(a_l-{\mathsf{\mu }}_l\right)}\right]\end{array}$$

(13)

with the mean and covariance given by:

$$\mathsf{\mu }=\mathsf{\Gamma }{D}^{\mathrm{H}}{\mathsf{\Sigma }}_x^{-1}X$$

(14)

$$\mathsf{\Sigma }=\mathsf{\Gamma }-\mathsf{\Gamma }{D}^{\mathrm{H}}{\mathsf{\Sigma }}_x^{-1}D\mathsf{\Gamma }$$

(15)

where ${\mathsf{\Sigma }}_x=\sigma I+D\mathsf{\Gamma }{D}^{\mathrm{H}}$.

The type-II maximum likelihood approach obtains the updated formula of the hyper-parameters by maximizing the marginal likelihood function in the $\mathsf{\gamma }$ space by integrating the unknown sparse coefficients out [24]. The cost function is expressed as:

$$G(\mathsf{\gamma })=\mathrm{tr}({\mathsf{\Sigma }}_x^{-1}X{X}^{\mathrm{H}})+L\log \left|{\mathsf{\Sigma }}_x\right|$$

(16)

where $\log (•)$ represents the natural logarithm.

Then the SBL can also be expressed as a standard sparse recovery problem form [24]:

$$\underset{A}{\min }\left|\right|X-DA|{|}_F^2+\sigma g(A)$$

(17)

where $g(A)$ represents a penalty term as follows:

$$g(A)=\underset{\mathsf{\gamma }\ge 0}{\min }G\left(A,\mathsf{\gamma }\right)$$

(18)

where

$$G(A,\mathsf{\gamma })=\mathrm{tr}(A^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}A)+L\log |\sigma I+D\mathsf{\Gamma }{D}^{\mathrm{H}}|$$

(19)

The penalty term in Equation (17) corresponds to a non-separable sparsity promoting penalty. Additionally, the optimization problem in Equation (17) can be solved by the alternating minimization method.

Step 1: Minimization with respect to $A$. In $t-\mathrm{th}$ iteration, we are given the estimate ${\mathsf{\gamma }}^{(t)}$, then get an updated formula for $A^{(t+1)}$:

$$\begin{array}{l}{A}^{(t+1)}={\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}}{(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}X\\ {\mathsf{\Gamma }}^{(t)}=\mathrm{diag}({\mathsf{\gamma }}^{(t)})\end{array}$$

(20)

Step 2: Minimization with respect to $\mathsf{\gamma }$. After being given the estimate $A^{(t)}$, we obtain the updated formula of ${\mathsf{\gamma }}^{(t+1)}$ by solving the optimization problem as follows:

$${\mathsf{\gamma }}^{(t+1)}=\underset{\mathsf{\gamma }\ge 0}{\arg \min }\mathrm{tr}({\left\{A^{(t+1)}\right\}}^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}{\left\{A^{(t+1)}\right\}}^{\mathrm{H}})+L\log |\sigma I+D\mathsf{\Gamma }{D}^{\mathrm{H}}|$$

(21)

Then we can get the updated formula of ${\mathsf{\gamma }}^{(t+1)}$ as:

$${\gamma }_i^{(t+1)}={\left(\frac{\frac{1}{L}{\displaystyle \sum _{j=1}^L|a_{i,j}{|}^2}}{d_i^{\mathrm{H}}{(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}{d}_i}\right)}^{1/2}$$

(22)

For the non-separable penalty in Equation (19), a new interpretation can be obtained through some simple transformations. Because of the non-negative constraint on $\mathsf{\gamma }$ in Equation (18), the solution of Equation (18) needs to satisfy the Karush–Kuhn–Tucker (KKT) condition [27]:

$${\gamma }_i\frac{\partial P(A,\mathsf{\gamma })}{\partial {\gamma }_i}={\gamma }_i\left(-\frac{{\displaystyle \sum _{l=1}^L|a_{i,l}{|}^2}}{{\gamma }_i^2}+L{d}_i^{\mathrm{H}}{(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}{d}_i\right)=0$$

(23)

Then, we can get the following equation:

$$\frac{\left|\right|A_{i,:}|{|}^2}{{\gamma }_i}=\left\{\begin{array}{cc}L{\gamma }_i{d}_i^{\mathrm{H}}{(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}{d}_i,& {\gamma }_i\ne 0\\ 0,& {\gamma }_i=0\end{array}\right.$$

(24)

where $A_{i,:}$ represents the $i-\mathrm{th}$ row of matrix $A$.

Therefore, the first item in Equation (19) can be rewritten as:

$$\begin{array}{l}\mathrm{tr}(A^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}A)={\displaystyle \sum _{i\in S}L{\gamma }_i{d}_i^{\mathrm{H}}{(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}{d}_i}\\ =L\mathrm{tr}({(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}{D}_S{\mathsf{\Gamma }}_S{D}_S^{\mathrm{H}})\\ =L\mathrm{tr}({(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}D\mathsf{\Gamma }D)\end{array}$$

(25)

where $S$ denotes the support set where ${\gamma }_i\ne 0$.

Let $q_j,(j=1,...,MN)$ denote the singular values of $D\mathsf{\Gamma }{D}^{\mathrm{H}}$, then Equation (25) can be rewritten as:

$$\mathrm{tr}(A^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}A)=L{\displaystyle \sum _{j=1}^{MN}\frac{q_j}{\sigma +q_j}}$$

(26)

In the same way, the second item in Equation (19) can be expressed as:

$$L\log |\sigma I+D\mathsf{\Gamma }{D}^{\mathrm{H}}|=L{\displaystyle \sum _{j=1}^{MN}\log (\sigma +q_j)}$$

(27)

In many rank norm minimization optimization problems, the matrix singular value is often used as a new sparsity metric to approximate the rank norm [28,29,30]. Therefore, we can treat the functions in Equations (26) and (27) as non-convex rank surrogates. Substituting Equations (26) and (27) into Equation (18) yields:

$$g(A)=L{f}_1(D{\mathsf{\Gamma }}^{1/2})+L{f}_2(D{\mathsf{\Gamma }}^{1/2})$$

(28)

where $f_1(•)$ and $f_2(•)$ represent two different non-convex rank surrogates, respectively.

We consider the situation where $D_S{\mathsf{\Gamma }}_S^{1/2}$ is a full column rank. Then we have:

$$\mathrm{rank}(D_S{\mathsf{\Gamma }}_S{D}_S^{\mathrm{H}})=\left|S\right|=\left|\right|\mathsf{\gamma }|{|}_0$$

(29)

Therefore, the role of the above rank surrogates in SBL can be understood as a non-separable approximation of the $\left|\right|A|{|}_{2,0}$.

3. Proposed Method

3.1. Subspace Penalty

Although the above non-separable penalty brings a lot of advantages, as shown in Equation (28), two different non-convex rank surrogates make the analysis difficult and do not use any subspace information. In sparse recovery, the Schatten $p-quasi$ norm is a special norm, which is defined on the singular values of a matrix and can be used to solve various low-rank problems. Compared with the nuclear norm, the Schatten $p-quasi$ norm has one more parameter $p$, which acts on the power of the singular value of the matrix. Therefore, the Schatten $p-quasi$ norm can be viewed as a weighted nuclear norm. In addition, the Schatten $p-quasi$ norm has a more flexible form than the nuclear norm and can reduce the weight of large singular values in the objective function. Therefore, the Schatten $p-quasi$ norm is more suitable for describing the rank minimization problem [31,32].

Inspired by the ${\ell }_0$ norm and rank norm similarity, a simpler and more intuitive Schatten $p-quasi$ norm rank surrogate is used to replace the penalty in SBL in this paper. This new penalty, named the subspace penalty (SPL), can be expressed as:

$$g_{SPL}(A)=2L\mathrm{tr}|U^{\mathrm{H}}D{\mathsf{\Gamma }}^{1/2}{|}^p,(0<p\le 1)$$

(30)

where $U$ denotes the incomplete subspace of the clutter. The Schatten $p-quasi$ norm rank surrogate is defined as:

$$\mathrm{tr}|Q{|}^p=\mathrm{tr}({(Q{Q}^{\mathrm{H}})}^{p/2})={\displaystyle \sum _{j=1}^{MN}{q}_j^p(Q)}$$

(31)

Similar to the derivation of Equation (18) in Section 2, we can rewrite Equation (29) as follows:

$$g_{SPL}(A)=\underset{\mathsf{\gamma }\ge 0}{\min }{G}_{SPL}(A,\mathsf{\gamma })$$

(32)

where

$$G_{SPL}(A,\mathsf{\gamma })=\mathrm{tr}(A^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}A)+L\mathrm{tr}({(U^{\mathrm{H}}A\mathsf{\Gamma }{A}^{\mathrm{H}}U)}^{p/2})$$

(33)

Then a new optimization problem is established based on the new penalty in Equation (32),

$$\underset{A}{\min }\left|\right|X-DA|{|}_F^2+\sigma g_{SPL}(A)$$

(34)

It is worth noting that $\mathrm{tr}({(U^{\mathrm{H}}A\mathsf{\Gamma }{A}^{\mathrm{H}}U)}^{p/2})$ is a non-convex function with respect to its singular values, and thus replace the SPL penalty in Equation (32) with its convex conjugate [33].

$$G_{SPL}(A,\mathsf{\gamma })=\underset{\mathsf{\Psi }\in \mathrm{T}}{\min }{F}_{SPL}(A,\mathsf{\gamma },\mathsf{\Psi })$$

(35)

where $F_{SPL}(A,\mathsf{\gamma },\mathsf{\Psi })$ can be expressed as:

$$F_{SPL}(A,\mathsf{\gamma },\mathsf{\Psi })=\mathrm{tr}(A^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}A)+Lp(\mathrm{tr}((U^{\mathrm{H}}A\mathsf{\Gamma }{A}^{\mathrm{H}}U)\mathsf{\Psi })-\frac{2}{s}\mathrm{tr}({\mathsf{\Psi }}^{s/2}))$$

(36)

where $\mathsf{{\rm T}}$ denotes the set of symmetric positive semi-definite matrices and $p,s$ satisfies

$$\frac{2}{p}+\frac{2}{s}=1$$

(37)

Although Equation (35) is not jointly convex, its greatest advantage over Equation (32) is its convexity for a single variable, while keeping the rest of the variables constant. Therefore, we can use an alternating minimization algorithm to obtain a closed-form solution at each minimization step. Similar to the derivation in Section 2.2, we can prove the equivalence of Equations (29) and (34).

Now, the final form of the optimization problem based on the SPL penalty can be written as:

$$\ell (A,\mathsf{\gamma },\mathsf{\Psi })=\left|\right|X-DA|{|}_F^2+\sigma \left\{\mathrm{tr}(A^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}A)+Lp(\mathrm{tr}((U^{\mathrm{H}}A\mathsf{\Gamma }{A}^{\mathrm{H}}U)\mathsf{\Psi })-\frac{2}{s}\mathrm{tr}({\mathsf{\Psi }}^{s/2}))\right\}$$

(38)

Since Equation (37) is convex for $A$, $\mathsf{\gamma }$ and $\mathsf{\Psi }$, respectively, an alternating minimization algorithm is used to solve Equation (37). The solution of Equation (37) also should satisfy the KKT conditions, as we have:

$$\frac{\partial \ell (A,\mathsf{\gamma },\mathsf{\Psi })}{\partial A}=-2D^{\mathrm{H}}(X-DA)+\sigma {\mathsf{\Gamma }}^{-1}A=0$$

(39)

$$\frac{\partial \ell (A,\mathsf{\gamma },\mathsf{\Psi })}{\partial \mathsf{\Psi }}=Lp(U^{\mathrm{H}}A\mathsf{\Gamma }{A}^{\mathrm{H}}U)-Lp{\mathsf{\Psi }}^{s/2-1}=0$$

(40)

$$\frac{\partial \ell (A,\mathsf{\gamma },\mathsf{\Psi })}{\partial {\gamma }_i}-{\mu }_i=-\frac{{\displaystyle \sum _{l=1}^L|a_{i,l}{|}^2}}{{\gamma }_i^2}+L{d}_i^{\mathrm{H}}U\mathsf{\Psi }{U}^{\mathrm{H}}{d}_i-{\mu }_i=0$$

(41)

$$\begin{array}{l}-{\mu }_i=0,\forall i\\ {\mu }_i{\gamma }_i=0,{\mu }_i\ge 0,{\gamma }_i\ge 0,\forall i\end{array}$$

(42)

According to Equations (38)–(41), we can get the following fixed-point iterations:

Step 1: Minimization with respect to $A$. In the $t-\mathrm{th}$ iteration, we are given the estimate ${\mathsf{\gamma }}^{(t)}$, then get an updated formula for $A^{(t+1)}$:

$$\begin{array}{l}{A}^{(t+1)}={\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}}{(\sigma I+D{\mathsf{\Gamma }}^{(t)}{D}^{\mathrm{H}})}^{-1}X\\ {\mathsf{\Gamma }}^{(t)}=\mathrm{diag}({\mathsf{\gamma }}^{(t)})\end{array}$$

(43)

Step 2: Update of $\mathsf{\Psi }$. The closed-form solution for ${\mathsf{\Psi }}^{(t+1)}$ can be expressed as:

$${\mathsf{\Psi }}^{(t+1)}={(U^{\mathrm{H}}A{\mathsf{\Gamma }}^{(t+1)}{A}^{\mathrm{H}}U)}^{\frac{1}{s/2-1}}$$

(44)

Step 3: Update of $\mathsf{\gamma }$. For the given ${\mathsf{\Psi }}^{(t)}$ and $A^{(t)}$, according to Equations (40) and (41), we have

$${\gamma }_i\frac{\partial \ell (A,\mathsf{\gamma },\mathsf{\Psi })}{\partial {\gamma }_i}={\gamma }_i\left(-\frac{\left|\right|A_{i,:}^{(t)}|{|}^2}{{\gamma }_i^2}+L{d}_i^{\mathrm{H}}U{\mathsf{\Psi }}^{(t)}{U}^{\mathrm{H}}{d}_i\right)=0$$

(45)

then the update formula for $\mathsf{\gamma }$ can be written as:

$${\gamma }_i=\left\{\begin{array}{cc}\infty & \mathrm{if}{d}_i^{\mathrm{H}}U{\mathsf{\Psi }}^{(t)}{U}^{\mathrm{H}}{d}_i=0\mathrm{and}\left|\right|A_{i,:}^{(t)}|{|}^2\ne 0\\ {\left(\frac{\frac{1}{L}{\displaystyle \sum _{j=1}^L|a_{i,j}{|}^2}}{d_i^{\mathrm{H}}U{\mathsf{\Psi }}^{(t)}{U}^{\mathrm{H}}{d}_i}\right)}^{1/2},& \mathrm{if}{d}_i^{\mathrm{H}}U{\mathsf{\Psi }}^{(t)}{U}^{\mathrm{H}}{d}_i\ne 0\\ const,& \mathrm{if}{d}_i^{\mathrm{H}}U{\mathsf{\Psi }}^{(t)}{U}^{\mathrm{H}}{d}_i=0\mathrm{and}\left|\right|A_{i,:}^{(t)}|{|}^2=0\end{array}\right.$$

(46)

After obtaining the sparse coefficient $A$, we can get the estimated CNCM $\stackrel{⌢}{R}$ by:

$$\stackrel{⌢}{R}=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{k=1}^{\mathrm{K}}\left|\right|a_{l,k}|{|}^{{}^2}}}{d}_k{d}_k^H+\alpha I_{MN}$$

(47)

where $\alpha$ is the load factor.

As shown in Equation (45), even in the iterations, we mistakenly set the $\left|\right|A_{i,:}^{(t)}|{|}^2$ to zero, but the corresponding hyperparameter still has a probability of non-zero when $d_i^{\mathrm{H}}U{\mathsf{\Psi }}^{(t)}{U}^{\mathrm{H}}{d}_i=0$, which makes the proposed algorithm have the ability to correct this error. This is different from M-SBL, where once we mistakenly set the $\left|\right|A_{i,:}^{(t)}|{|}^2$ to zero, the mistake will continue in subsequent iterations. Therefore, the proposed algorithm has a stronger applicability and maintains a good performance even in the case of SMV.

3.2. Adaptive Update Support Set and Calculate $U$

Although the algorithm based on the SPL penalty proposed in the previous section had a good recovery performance, it needed to calculate the inverse of the matrix in the iterative process, which leads to a huge computational burden. In particular, the degree of freedom of modern airborne early warning radar systems is extremely high, which severely limits the application of algorithms in practice.

Due to the inherent sparsity of clutter signals, there are already a large number of redundant atoms in the dictionary $D$ that do not contribute to the model. In addition, when the radar is working, the energy of the area irradiated by the main beam is much higher than the area irradiated by the side lobes. Inspired by this, we considered setting a threshold adaptively in each iteration to update the support set $S$. In the next iteration, we ignored the atoms outside the supported set, so that the dimension of the dictionary $D$ was gradually reduced, then the computational complexity could be effectively reduced.

It is worth noting that the selection of the threshold affects the performance of the proposed algorithm. When setting the threshold, it is necessary to balance the computational efficiency and the recovery accuracy. In this paper, we set the threshold as:

$${\eta }^{(t+1)}=\kappa \max (|\left|A^{(t)}\right|{|}_F^2)$$

(48)

where $\kappa$ denotes an artificially set parameter used to balance the computational efficiency and recovery accuracy.

Then we could get the update of support set $S$ as follows:

$$S^{(t+1)}=\left\{i\right|{\gamma }_i>{\eta }^{(t)}\},i=1,...,K$$

(49)

Moreover, we needed to obtain the estimation methods for the incomplete clutter subspace $U$. Usually, we use the eigen-decomposition of the covariance matrix to obtain the clutter subspace. However, such estimates are inaccurate when the number of samples is insufficient. In this paper, we only needed to estimate an incomplete partial clutter subspace $U$.

First, the sample covariance matrix $\tilde{R}$ can be obtained by maximum likelihood estimation.

$$\tilde{R}=\frac{X{X}^{\mathrm{H}}}{L}$$

(50)

Then, the biased matrix $\widehat{R}$ can be computed by [18]:

$$\widehat{R}=\tilde{R}-{\lambda }_{MN}(\tilde{R})I$$

(51)

where ${\lambda }_{MN}(\tilde{R})$ denotes the $MN-\mathrm{th}$ eigenvalue of matrix $\tilde{R}$. The matrix $\tilde{R}$ and $\widehat{R}$ have the same eigenvectors.

Next, the incomplete partial clutter subspace $U$ can be expressed as:

$$U=[v_1,v_2,...,v_r]$$

(52)

where $v_j,j=1,2...,r$ denotes the $j-\mathrm{th}$ eigenvector of $\widehat{R}$. The $r$ can be determined by the following criteria:

$$\frac{{\lambda }_r(\tilde{R})-{\lambda }_{r+1}(\tilde{R})}{{\lambda }_1(\tilde{R})}\ge \tau \ge \frac{{\lambda }_k(\tilde{R})-{\lambda }_{k+1}(\tilde{R})}{{\lambda }_1(\tilde{R})},\forall k>r$$

(53)

where $\tau$ is a non-negative threshold.

The pseudocode of the proposed algorithm is shown in Algorithm 1.

Algorithm 1. Pseudocode for ASDSPL-STAP algorithm.

Input: training samples $X$, dictionary matrix $D$.

Initialize: $t=0$, $\kappa =0.01$, $\epsilon =\tau ={10}^{-6}$, $p=0$, ${\mathrm{iter}}_{\max }=150$.

Step1: Calculate $U$ using (49)-(52); Step2: While $t\le {\mathrm{iter}}_{\max }$ $\mathrm{ t}\leftarrow t+1$;                 Update $A^{(t)}$ using (42);                 Update ${\mathsf{\Psi }}^{(t)}$ using (43);                 Update ${\mathsf{\gamma }}^{(t)}$ using (45);                 Update support set $S^{(t)}$ using (47) and (48);                 Update dictionary matrix $D^{(t)}=D_S^{(t-1)}$;                 If $\left|\right|{\mathsf{\gamma }}^{(t)}-{\mathsf{\gamma }}^{(t-1)}|{|}_2/\left|\right|{\mathsf{\gamma }}^{(t-1)}|{|}_2<\epsilon$                    break               end end

Step2: Estimate the CNCM $\stackrel{⌢}{R}$ by (46)

Step3: Calculate the STAP weight $w$ by (8);

Output: $Y=w^{\mathrm{H}}X$.

4. Computational Complexity Analysis

In this section, we discuss the computational complexity of the proposed algorithm in detail and compare it with some common SR-STAP algorithms, such as M-CVX, M-OMP, M-FOCUSS, M-IAA, M-SBL and M-FCSBL. The number of multiplications is used to evaluate the computational complexity of the algorithm. For simplicity, the low-order terms of multiplications are omitted.

The calculation amount of the proposed algorithm and M-SBL is mainly concentrated in the matrix inversion of matrix $U^{\mathrm{H}}A\mathsf{\Gamma }{A}^{\mathrm{H}}U$ and matrix $D\mathsf{\Gamma }{D}^{\mathrm{H}}$. Compared with M-SBL, the proposed algorithm reduces the computational complexity from two aspects: (1) with the SPL penalty, the dimension of matrix $U^{\mathrm{H}}A\mathsf{\Gamma }{A}^{\mathrm{H}}U$ is $r\times r$, and the $r$ satisfies $r<MN$. (2) By updating the atoms corresponding to the support set in each iteration, the dimension of the dictionary $K_{ASD}$ gradually decreases as the number of iterations increases. The relationship between the cardinality of the support set and iteration sequence number t is shown in Figure 2. It can be seen from Figure 2 that the cardinality of the support set drops sharply in the first few iterations, and then tends to converge.

From the above analysis, it can be seen that the ASDSPL algorithm has lower computational complexity than M-SBL. The computational complexities of different algorithms are shown in

Table 1. Computational complexity comparison.

Algorithm	Computational Load for a Single Iteration
M-OMP	$O\left(NMKL+r_s{}^3+NM{r}_s{}^2+2NM{r}_s\right)$
M-FOCUSS	$O\left(NMKL+{\left(NM\right)}^3+2{\left(NM\right)}^2{N}_s{N}_d\right)$
M-IAA	$O\left(NMKL+{\left(NM\right)}^3+2{\left(NM\right)}^2{N}_s{N}_d\right)$
M-SBL	$O\left({(NM)}^3+2{(NM)}^2{N}_s{N}_d+3K^{2}NM+NM{N}_s{N}_{d}L\right)$
M-FCSBL	$O\left({(NM)}^3+5K{(NM)}^2+3K+(2KL+4K+L)NM+L\right)$
ASDSPL	$O\left({(r)}^3+2{(NM)}^2{K}_{ASD}+3K_{ASD}^{2}NM+NM{K}_{ASD}L\right)$

, where we assumed that the sparsity $r_s$ is equal to the clutter rank in M-OMP.

5. Numerical Simulation

In this section, we evaluate the performance of the proposed M-F-LSBL algorithm using simulated data and compare it with the existing STAP algorithms. The IF and SINR loss are used as the evaluation metrics of the clutter suppression performance.

$$\mathrm{IF}=\frac{|w^H{s}_t{|}^2\mathrm{tr}(R)}{w^{H}Rw{s}_t^H{s}_t}$$

(54)

$$L_{SINR}=\frac{{\sigma }^2}{MN}\frac{|w^H{s}_t{|}^2}{w^{H}Rw}$$

(55)

where $R$ is the known CNCM of the CUT, and $w$ denotes the STAP weight vector.

A side-looking phased array radar was considered in our experiment. The main simulation parameters are shown in

Table 2. List of simulation parameters.

Symbols	Parameters	Value
$\lambda$	Wavelength	0.3 m
$d$	Element spacing	0.15 m
$\vartheta$	Platform velocity	150 m/s
$H$	height of platform	9000 m
$M$	Number of pulses	8
$N$	Number of array elements	8
$f_r$	PRF	2000 Hz
$f_s$	Sampling frequency	2.5 MHz
$\theta$	Azimuth angle	90°
$\mathrm{CNR}$	Clutter to noise ratio	50 dB

. The CUT is $\mathit{600}-\mathrm{th}$ range gate. We set the ${\rho }_s=5$, ${\rho }_d=5$ and $L=6$. Moreover, all presented results are averaged over 100 independent Monte Carlo runs.

5.1. Performance Analysis in the Ideal Case

In this subsection, an ideal error free simulation scenario is considered. First, we give an intuitive view of the reconstructed Capon spectrum of ASDSPL, and compare the Capon spectrum obtained by the ideal CNCM, M-OMP algorithm, M-FOCUSS algorithm and other MMV SR-STAP algorithms. From Figure 3, we can observe that the Capon spectrums in Figure 3d–g obtained from the M-IAA, M-SBL, M-FCSBL and M-F-LSBL algorithms are similar to each other and are close to the ideal Capon spectrum. The value of the color bar corresponds to the energy value in the figure. We can obtain a more realistic clutter Capon spectrum via the MIAA, M-SBL, M-FCSBL and ASDSPL algorithms. The M-OMP selects the support set and then uses the least square to calculate the power, which brings additional errors.

In order to better illustrate the effectiveness of the proposed algorithm, the eigen-spectra of CNCMs obtained by different STAP methods are shown in Figure 4. Under the simulation conditions set in this subsection, the clutter rank is equal to 15 according to Brennan’s rule [1]. From Figure 4, we can see that the eigen-spectra obtained by the ideal CNCM decreases significantly when the eigenvalue index is equal to 15. Moreover, the eigen-spectra obtained by the M-IAA, M-SBL, M-FCSBL and ASDSPL algorithms is close to the optimal eigen-spectra. Therefore, the M-IAA, M-SBL, M-FCSBL and ASDSPL algorithms can achieve better clutter suppression performance compared with the M-OMP algorithm and M-FOCUSS algorithm.

In addition, the clutter suppression performance of different SR-STAP algorithms is also evaluated by the IF against the normalized Doppler frequency, as shown in Figure 5. We found that the IF curve obtained by the M-IAA, M-SBL, M-FCSBL and ASDSPL algorithms is close to the optimal IF curve regardless of the depth or shape of the notch.

In this paper, the relationship between the cost function $C=\ln \left|\stackrel{⌢}{R}\right|+\mathrm{Tr}\left({\stackrel{⌢}{R}}^{-1}R\right)$ in [23] and the number of iterations is used to evaluate the convergence performance of different algorithms. $R$ is the exact CNCM. In addition, we assume that the noise power can be obtained by the silence of the receiver. As shown in Figure 6, we only considered five algorithms with a better performance: the M-IAA, M-SBL, M-FCSBL, M-FOCUSS and the ASDSPL algorithms. From Figure 6, we can observe that the M-IAA, M-FCSBL and the ASDSPL algorithms converge faster and only need 20 iterations to converge to the steady state, while the M-SBL and M-FOCUSS algorithms converge slower and need 150 iterations to reach the steady state.

In the next experiment, the average running time was used to assess the computational efficiency of different algorithms. The experiment results are presented in

Table 3. Average running time comparison.

Algorithm	Running Time
M-OMP	0.025 s
M-FOCUSS	11.294 s
M-IAA	4.368 s
M-SBL	16.511 s
M-FCSBL	4.213 s
ASDSPL	0.119 s

. Note that all simulations results are implemented by MATLAB 2021b on a computer equipped with Intel Xeon E5-2620 CPU @2.10GHz×2. Table 3 shows that the average running time of the ASDSPL algorithm is significantly less than other algorithms. The great computational efficiency of the ASDSPL algorithm comes from its low computational complexity and fast convergence speed.

Finally, the average SINR, which is defined as the mean of the $L_{SINR}$ for $f_d\in (-0.5,-0.1)\cup (0.1,0.5)$ against the number of snapshots, is shown in Figure 7. It can be observed that the M-IAA, M-SBL, M-FCSBL and the ASDSPL algorithms achieved the best performance with only three snapshots. The experimental results clearly prove the great performance of the proposed algorithm when the number of i.i.d training samples is small.

5.2. Performance Analysis in the Non-Ideal Case

In this subsection, we consider the non-ideal case with gain-phase (GP) errors. The GP error is caused by the inconsistent amplitude and phase characteristics of the RF amplifier part of the array channel. The GP error manifests itself as an inconsistency in the amplitude and phase of the individual channels. It can be described by introducing an error matrix $T$ into the steering vector modeling. Then, the signal model with the GP error can be expressed as:

$$X_{GP}=TDA+X_n$$

(56)

where error matrix $T$ is:

$$T=I_M\otimes \mathrm{diag}({[{\rho }_1\exp (j{\xi }_1),...,{\rho }_N\exp (j{\xi }_N)]}^{\mathrm{T}})$$

(57)

where ${\rho }_i,{\xi }_i,i=1,...,N$ represents the amplitude error and phase error of $i-\mathrm{th}$ element.

The standard deviation of the gain error and phase error is 0.03 and 2°. Figure 8 shows the Capon spectrum recovered by the different SR-STAP algorithms in the case of a GP error. We can find that the existence of the GP error leads to the broadening of the Capon spectrum. Meanwhile, the M-SBL algorithm and ASDSPL algorithm can achieve a better recovery performance of the Capon spectrum, both in terms of location and power. Although, there is some performance loss compared to the error-free case, which is caused by the mismatch between the data and the dictionary.

Next, we illustrate the influence of the GP error by the eigen-spectra of CNCMs obtained by different STAP methods. As shown in Figure 9, the eigen-spectra obtained by the M-SBL and ASDSPL algorithms is close to the optimal eigen-spectra. The IF curve is shown in Figure 10; we find that the performance of the M-IAA and M-FCSBL algorithms is degraded due to the use of the standard Capon beamformer in each iteration. The standard Capon beamformer exhibits a poor performance in the presence of GP errors due to the mismatch of steering vectors. Moreover, compared with the ideal IF curve, the M-SBL and ASDSPL algorithms had a slight performance loss. The main lobe notch of the IF curve obtained by the M-SBL and ASDSPL algorithms has slightly widened, and the side lobe has decreased by 3 dB.

6. Conclusions

In this work, an adaptive support-driven sparse recovery STAP method with a subspace penalty was proposed. First, we developed a new optimization objective function-based subspace penalty. The objective function combines the information of SBL and the subspace, and achieves a great recovery performance. Second, the alternating minimization algorithm was used to get the closed-form solution of the optimization function. Due to the convexity of the objective function for each individual variable, the convergence of the proposed algorithm is guaranteed. Finally, an adaptive updating support set method was used to improve computational efficiency. The simulation results show that the proposed algorithm has an excellent clutter suppression performance and great computational efficiency.