A Fast Space-Time Adaptive Processing Algorithm Based on Sparse Bayesian Learning for Airborne Radar

Abstract

Space-time adaptive processing (STAP) plays an essential role in clutter suppression and moving target detection in airborne radar systems. The main difficulty is that independent and identically distributed (i.i.d) training samples may not be sufficient to guarantee the performance in the heterogeneous clutter environment. Currently, most sparse recovery/representation (SR) techniques to reduce the requirement of training samples still suffer from high computational complexities. To remedy this problem, a fast group sparse Bayesian learning approach is proposed. Instead of employing all the dictionary atoms, the proposed algorithm identifies the support space of the data and then employs the support space in the sparse Bayesian learning (SBL) algorithm. Moreover, to extend the modified hierarchical model, which can only apply to real-valued signals, the real and imaginary components of the complex-valued signals are treated as two independent real-valued variables. The efficiency of the proposed algorithm is demonstrated both with the simulated and measured data.

1. Introduction

Space-time adaptive processing (STAP) has the capability to detect slow-moving targets that might otherwise be swallowed up by the strong sidelobe clutter. The performance of STAP filter is dependent on the accuracy of the clutter plus noise covariance matrix (CNCM) of the cell under test (CUT) [1]. According to the well-known Reed-Mallet-Brennan (RMB) rule [2], to achieve the steady performance, the number of the independent and identically distributed (i.i.d) secondary range cells used to estimate CNCM is no less than twice the system degrees of freedom (DOF). However, it is hard to obtain enough i.i.d samples in practice because of array geometry structures, nonhomogeneous clutter environment and so on [1]. Knowing how to improve the performance of STAP with limited samples has been a hot topic until now.

Reduced-dimension (RD) [3] and reduced-rank (RR) [4] algorithms are proposed to improve the performance of STAP with limited secondary samples. Even though they are easy to implement, their performance gets worse when the number of secondary samples is smaller than twice the rank of clutter [5].

Apart from RD and RR algorithms, some other algorithms [6,7,8,9,10] are proposed and succeed in suppressing clutter in theory. However, they face some disadvantages in practice. Sufficiently, (1) direct data domain (DDD) algorithms in [6] that achieve enough samples from the CUT suffer from aperture loss; and (2) knowledge-aided (KA) algorithms in [7,8,9,10] state that the accurate prior knowledge must be required in advance. Since the cost to achieve the accurate prior knowledge is expensive and the prior knowledge changes with time, KA algorithms are not widely used in practical applications.

Over the past twenty years, sparse recovery/representation (SR) algorithms have received continuous attention in STAP [11,12,13,14,15] because they have enormous potential to estimate the clutter spectrum with limited samples. The sparse Bayesian learning (SBL)-type algorithms are robust and have drawn continuous attention in the last five years.

The SBL algorithm was proposed by Tipping in [16]. Wipf enhanced it in the single measurement vector (SMV) case in [17] and then the multiple measurement vectors (MMV) case in [18]. Due to the excellent performance of SBL algorithms, the fast marginal likelihood maximization (FMLM) algorithm [19], the Bayesian compressive sensing (BCS) algorithm [20], the multitask BCS (M-BCS) algorithm [21], and other Bayesian algorithms [22,23,24] were proposed by researchers in the following years. SBL was firstly introduced into STAP with MMV, defined as the M-SBL-STAP algorithm in [25], by Duan to estimate CNCM and Wang improved the fast convergence of the algorithm in [26]. SBL has also been used to solve some common problems in STAP in the last five years, such as off-grid in [27] and discrete interference in [28]. However, the Bayesian algorithms aforementioned have one or more of the following disadvantages: (a) high computational cost, (b) inaccurate estimation of the noise variance and (c) inapplicability to complex-valued signals.

In this paper, to improve the computational efficiency of the M-SBL-STAP algorithm, we extend the FMLM algorithm to the conventional M-SBL-STAP algorithm. The real and imaginary components of the complex-valued signals are treated as two independent real-valued variables in order to extend the modified hierarchical model. Simulation results with both simulated and Mountain-Top data demonstrate that the proposed algorithm can achieve high computational efficiency and good performance.

The main contributions of this paper can be listed as follows:

• We extend the FMLM algorithm into M-SBL-STAP for the purpose of identifying the support space of the data, i.e., the atoms whose corresponding hyper-parameters are non-zero. After support space identification, the dimensions of the effective problems are drastically reduced due to sparsity, which can reduce computational complexities and alleviate memory requirements.
• According to [18], we have no access to obtain the accurate value of the noise variance. Instead of estimating the noise variance, we extend the modified hierarchical model, introduced in Section 4, to the SBL framework.
• Although the hierarchical models apply to the real-valued signals, they cannot be extended directly to the complex-valued signal according to [29,30]. The data needed to be dealt with in STAP are all complex-valued. To solve the problem, we transform sparse complex-valued signals into group sparse real-valued signals.

Notation: In this article, scalar quantities are denoted with italic typeface. Boldface small letters are reserved for vectors, and boldface capital letters are reserved matrices. The $i-\mathrm{th}$ entry of a vector $x$ is denoted by $x_i$, while $A_i$ and ${\mathrm{A}}_{ij}$ denote the $i-\mathrm{th}$ row and $\left(i,j\right)-\mathrm{th}$ element of a matrix $A$, respectively. The symbols ${\left(•\right)}^{\mathrm{T}}$ and ${\left(•\right)}^{\mathrm{H}}$ denote the matrix transpose and conjugate transpose, respectively. The symbols ${‖•‖}_1$, ${‖•‖}_2$ and ${‖•‖}_{\mathrm{F}}$ are reserved for ${\ell }_1$, ${\ell }_2$ and Frobenius (${\ell }_{\mathrm{F}}$) norms, respectively. ${‖•‖}_0$ is reserved for ${\ell }_0$ pseudo-norm. ${‖•‖}_{2,0}$ stands for a mixed norm defined as the number of non-zero elements of the vector formed by the ${\ell }_2$ norm of each row vector. The symbol $\left|•\right|$ is reserved for the determinant. The notations $I$, $0$ and $1$ represent the identity matrix, the all zero matrix and the all one matrix, respectively. The expectation of a random variable is denoted as $\mathsf{{\rm E}}\left(•\right)$.

2. Background and Problem Formulation

2.1. STAP Signal Model for Airborne Radar

Consider an airborne pulsed-Doppler radar system with a side-looking uniform linear array (ULA) consisted of $N$ elements. A coherent burst of $M$ pulses is transmitted at a constant pulse repetition frequency (PRF) in a coherent processing interval (CPI). The complex sample at the CUT from the $n-\mathrm{th}$ element and the $m-\mathrm{th}$ pulse is denoted as $x_{mn}$, and the data for the CUT can be written as a $MN\times 1$ vector $x$, termed a space-time snapshot.

$$x={\left[x_{11},x_{12},\cdots ,x_{1N},x_{21},\cdots ,x_{MN}\right]}^{\mathrm{T}}$$

(1)

Radar systems need to ascertain whether the targets are present in the CUT; thus, target detection is formulated into a binary hypothesis problem: the hypothesis ${\mathrm{H}}_0$ represents target absence and the other hypothesis ${\mathrm{H}}_1$ represents target presence.

$$\left\{\begin{array}{cc}x=x_{\mathrm{c}}+n& {\mathrm{H}}_0\\ x={\alpha }_{\mathrm{t}}{S}_{\mathrm{t}}+x_{\mathrm{c}}+n& {\mathrm{H}}_1\end{array}\right.$$

(2)

where ${\alpha }_{\mathrm{t}}$ is the target amplitude and $S_{\mathrm{t}}\in {ℂ}^{MN\times 1}$ is the space-time vector of the target. The vector $n\in {ℂ}^{MN\times 1}$ is the thermal noise, which is uncorrelated both spatially and temporally. A general model for the space-time clutter $x_{\mathrm{c}}\in {ℂ}^{MN\times 1}$ is

$$x_{\mathrm{c}}={\displaystyle \sum _{k=1}^{N_{\mathrm{c}}}{\alpha }_k\upsilon \left(f_{\mathrm{d},k},f_{\mathrm{s},k}\right)}$$

(3)

$${\upsilon }_{\mathrm{st}}\left(f_{\mathrm{d},k},f_{\mathrm{s},k}\right)={\upsilon }_{\mathrm{t}}\left(f_{\mathrm{d},k}\right)\otimes {\upsilon }_{\mathrm{s}}\left(f_{\mathrm{s},k}\right)$$

(4)

$${\upsilon }_{\mathrm{t}}\left(f_{\mathrm{d},k}\right)={\left[1,\exp \left(j\pi f_{\mathrm{d},k}\right),\dots ,\exp \left(j\pi \left(M-1\right)f_{\mathrm{d},k}\right)\right]}^{\mathrm{T}}$$

(5)

$${\upsilon }_{\mathrm{s}}\left(f_{\mathrm{s},k}\right)={\left[1,\exp \left(j\pi f_{\mathrm{s},k}\right),\dots ,\exp \left(j\pi \left(N-1\right)f_{\mathrm{s},k}\right)\right]}^{\mathrm{T}}$$

(6)

where ${\alpha }_k$ denotes the random amplitude of the echo from the corresponding clutter patch; $N_{\mathrm{c}}$ denotes the number of clutter patches in a clutter ring; ${\upsilon }_{\mathrm{st}}$, ${\upsilon }_{\mathrm{t}}$ and ${\upsilon }_{\mathrm{s}}$ denote space-time steering vector, temporal steering vector and spatial steering vector, respectively; and $f_{\mathrm{d},k}$ and $f_{\mathrm{s},k}$ denote the corresponding normalize temporal frequency and spatial frequency, respectively.

According to the linearly constrained minimum variance (LCMV) criterion, the optimal STAP weight vector is

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

(7)

where the CNCM $R_{\mathrm{c}+\mathrm{n}}$ is expressed as

$$R_{\mathrm{c}+\mathrm{n}}=R_{\mathrm{c}}+R_{\mathrm{n}}=\mathrm{E}\left\{x_{\mathrm{c}}{x}_{\mathrm{c}}^{\mathrm{H}}\right\}+{\sigma }^{2}I$$

(8)

2.2. SR-STAP Model and Principle

Discretize the space-time plane uniformly into $K=N_S{N}_d$ grids, where $N_S={\phi }_{S}N$ $\left({\phi }_S>1\right)$ denotes the number of normalized spatial frequency bins and $N_d={\phi }_{d}M$ $\left({\phi }_d>1\right)$ denotes the number of normalized Doppler frequency bins. Each grid corresponds to a space-time steering vector ${\upsilon }_k , \left(k=1,2,\cdots , K\right)$. The dictionary $D\in {ℂ}^{MN\times K}$ used in STAP is the collection of space-time steering vectors of all grids.

$$D=\left[{\upsilon }_1,{\upsilon }_2,\cdots ,{\upsilon }_K\right]$$

(9)

The signal model in STAP can be cast in the following form

$$X=DA+n$$

(10)

where $X=\left[x_1,x_2,\cdots ,x_L\right]$; $A\in {ℂ}^{K\times L}$ denotes sparse coefficient matrix; non-zero elements indicate the presence of clutter on the space-time profile; $L$ denotes the number of the range gates in the data; and $n\in {ℂ}^{MN\times L}$ denotes zero-mean Gaussian noise.

In sparse signal recovery algorithms with MMV, our goal is to represent the measurement $X$, which is contaminated by noise as a linear combination of as few dictionary atoms as possible. Therefore, the objective function is expressed as

$$A^{\ast }=\underset{A}{\arg \min }{‖A‖}_{2,0}, \mathrm{s}.\mathrm{t}. {‖X-DA‖}_{\mathrm{F}}^2\le \epsilon$$

(11)

where $\epsilon$ is the noise error allowance.

3. M-SBL-STAP Algorithm

3.1. Sparse Bayesian Learning Formulation

In the SBL framework, ${\left\{x_l\right\}}_{l=1}^L$ are assumed as i.i.d training snapshots. The noise is submitted to white complex Gaussian distribution with unknown power ${\sigma }^2$. The likelihood function of the measurement vectors can be denoted as

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

(12)

Since the training snapshots are i.i.d, each column in $A$ is independent and shares the same covariance matrix. Assign to each column in $A$ with a zero-mean Gaussian prior

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

(13)

where $0\in {ℂ}^{K\times 1}$ is a zero vector and $\mathsf{\Gamma }=\mathrm{diag}\left(\gamma \right)$. $\gamma =\left\{{\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_K\right\}$ are unknown variance parameters corresponding to sparse coefficient matrix. The prior of $A$ can be expressed as

$$p\left(A|\mathsf{\Gamma }\right)={\pi }^{-KL}{\mathsf{\Gamma }}^{-L}\exp ⟮-{\displaystyle \sum _{l=1}^L{a}_l^{\mathrm{H}}{\mathsf{\Gamma }}^{-1}{a}_l}⟯$$

(14)

Combining likelihood with prior, the posterior density of $A$ can be easily expressed as

$$p\left(A|X,\mathsf{\Gamma },{\sigma }^2\right)=\frac{p\left(X|A,{\sigma }^2\right)p\left(A|\mathsf{\Gamma }\right)}{{\displaystyle \int p\left(X|A,{\sigma }^2\right)p\left(A|\mathsf{\Gamma }\right)dA}}$$

(15)

modulated by the hyper-parameters $\gamma$ and ${\sigma }^2$. To find the hyper-parameters $\gamma$ and ${\sigma }^2$, which are enough accurate to estimate the CNCM, the most common method is the expectation maximization (EM) algorithm. The EM algorithm is divided into two parts: E-step and M-step. $t$ represents the sequence number of the current iteration. At the E-step: treat $A$ as hidden data, and the posterior density can be described with hyper-parameters $\gamma$ and ${\sigma }^2$.

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

(16)

with covariance and mean given by

$${\Sigma }^{\left(t+1\right)}={\left({\left\{{\sigma }^{-2}\right\}}^{\left(t\right)}{D}^{\mathrm{H}}D+{\Lambda }^{\left(t\right)}\right)}^{-1}$$

(17)

$${\mathsf{\mu }}^{\left(t+1\right)}={\left\{{\sigma }^{-2}\right\}}^{\left(t\right)}{\Sigma }^{\left(t+1\right)}{D}^{\mathrm{H}}X$$

(18)

where

$${\Lambda }^{\left(t\right)}={\left\{{\mathsf{\Gamma }}^{\left(t\right)}\right\}}^{-1}=\mathrm{diag}⟮\left[\frac{1}{{\gamma }_1^{\left(t\right)}},\frac{1}{{\gamma }_2^{\left(t\right)}},\cdots ,\frac{1}{{\gamma }_K^{\left(t\right)}}\right]⟯$$

(19)

$${\mathsf{\mu }}^{\left(t+1\right)}=\left[{\mathsf{\mu }}_1^{\left(t+1\right)},{\mathsf{\mu }}_2^{\left(t+1\right)},\cdots ,{\mathsf{\mu }}_L^{\left(t+1\right)}\right]$$

(20)

At the M-step, we update the hyper-parameters by maximizing the expectation of $p\left(X,A^{\left(t+1\right)}|{\mathsf{\Gamma }}^{\left(t\right)},{\left\{{\sigma }^2\right\}}^{\left(t\right)}\right)$.

$$\left[{\gamma }^{\left(t+1\right)},{\left\{{\sigma }^2\right\}}^{\left(t+1\right)}\right]=\underset{{\mathsf{\Gamma }}^{\left(t\right)},{\left\{{\sigma }^2\right\}}^{\left(t\right)}}{\mathrm{argmax}}\mathrm{E}\left[\ln p\left(X,A^{\left(t+1\right)}|{\mathsf{\Gamma }}^{\left(t\right)},{\left\{{\sigma }^2\right\}}^{\left(t\right)}\right)\right]$$

(21)

where

$${\gamma }_i^{\left(t+1\right)}=\frac{1}{L}{\displaystyle \sum _{l=1}^L{‖{\mathsf{\mu }}_l^{\left(t+1\right)}\left(i\right)‖}_2^2}+{\Sigma }_{ii}^{\left(t+1\right)}$$

(22)

$${\left\{{\sigma }^2\right\}}^{\left(t+1\right)}=\frac{\frac{1}{L}{\displaystyle \sum _{l=1}^L{‖x_l-D{\mathsf{\mu }}_l^{\left(t+1\right)}‖}_2^2}}{NM-K+{\displaystyle \sum _{k=1}^K\frac{{\Sigma }_{kk}^{\left(t+1\right)}}{{\gamma }_k^{\left(t\right)}}}}$$

(23)

The M-SBL-STAP algorithm is shown in Algorithm 1.

Algorithm 1: Pseudocode for the M-SBL-STAP algorithm.

Step 1: Input: the clutter data X, the dictionary $D$

Step 2: Initialization: initial the values of $\gamma$ and ${\sigma }^2$.

Step 3: E-step: update the posterior moments ${\Sigma }^{\left(t+1\right)}$ and ${\mathsf{\mu }}^{\left(t+1\right)}$ using (17) and (18).

Step 4: M-step: update ${\gamma }^{\left(t+1\right)}$ and ${\left\{{\sigma }^2\right\}}^{\left(t+1\right)}$ using (22) and (23).

Step 5: Repeat step 3 and step 4 until a stopping criterion is satisfied.

Step 6: Estimate the CNCM $R$ by the formula      $R=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{k=1}^K\left({‖{\mathsf{\mu }}_{l,k}^{\ast }‖}^2\right){\upsilon }_k{\upsilon }_k^{\mathrm{H}}}}+\beta {\rm I}$  where $\beta$ is a load factor and the symbol * represents the stopping criterion.

Step 7: Compute the space-time adaptive weight $w$ using (7).

Step 8: The output of the M-SBL-STAP algorithm is $Z=w^{\mathrm{H}}X$.

3.2. Problem Statement of the M-SBL-STAP Algorithm

At the E-step in each iteration, the inversion of a $K\times K$ matrix, which brings high computational complexities in the order of $o\left(K^3\right)$, is required to be calculated when update covariance $\Sigma$ (17). $K$ is the number of the atoms in the dictionary, and the value of $K$ is usually large. If we avoid calculating the matrix inversion of a $K\times K$ matrix, the computational complexities can be reduced a lot.

At the M-step in each iteration, the noise variance ${\sigma }^2$ needs to be estimated. However, it has been demonstrated in [18] that ${\sigma }^2$ estimated by (23) can be extremely inaccurate when the dictionary is structured and $K\ge NM$. We can see in [21] that ${\sigma }^2$ is a nuisance parameter in the iterative procedure and an inappropriate value may contaminate the convergence performance of the algorithm.

We extend a modified hierarchical model to the SBL framework, which aims at integrating ${\sigma }^2$ out instead of estimating ${\sigma }^2$. However, the modified model, which applies to the real-valued signals, cannot be directly extended to complex-valued signals.

In the next sections, we introduce the proposed algorithm to solve the above problems. The proposed algorithm is defined as M-FMLM-STAP algorithm.

4. The Proposed M-FMLM-STAP Algorithm

4.1. Modified Hierarchical Model

In [16], the scholars define the following hierarchical hyperpriors over $\mathsf{\Gamma }$, as well as the noise variance ${\sigma }^2$.

$$p\left(\mathsf{\Gamma }|c,d\right)={\displaystyle \prod _{k=1}^K\mathrm{Gamma}⟮\frac{1}{{\gamma }_k}|c,d⟯}$$

(24)

$$p\left({\sigma }^2|a,b\right)=\mathrm{Gamma}\left({\sigma }^{-2}|a,b\right)$$

(25)

where $\mathrm{Gamma}\left(\alpha |a,b\right)=\Gamma {\left(a\right)}^{-1}{b}^a{\alpha }^{a-1}{e}^{-b\alpha }$ and the ‘gamma function’ $\Gamma \left(a\right)={\displaystyle {\int }_0^{\infty }{t}^{a-1}{e}^{-t}dt}$. It has been demonstrated in [16] that an appropriate choice of the shape parameter $a$ and the scale parameter $b$ encourages the sparsity of the coefficient matrix.

In the STAP framework, CNCM can be expressed as

$$R_{\mathrm{c}+\mathrm{n}}=D\mathsf{\Gamma }D+{\sigma }^{2}I$$

(26)

We can translate STAP weight vector in the form

$$\begin{array}{ll}{w}_{\mathrm{opt}}\hfill & =\mu R_{\mathrm{c}+\mathrm{n}}^{-1}{S}_{\mathrm{t}}\hfill \\ & =\lambda {⟮D\frac{\mathsf{\Gamma }}{{\sigma }^2}D+I⟯}^{-1}{S}_{\mathrm{t}}\hfill \\ & =\lambda {\tilde{R}}^{-1}{S}_{\mathrm{t}}\hfill \end{array}$$

(27)

where $\lambda$ and $\mu$ are constants.

The above analysis shows that $\tilde{R}$ is equivalent to $R_{\mathrm{c}+\mathrm{n}}$ in the performance of clutter suppression. Thus, including ${\sigma }^2$ in the prior of $\mathsf{\Gamma }$, each component of $A$ is defined as a zero-mean Gaussian distribution and the modified hierarchical model follows (28) and (29).

$$p\left(A|{\mathsf{\Gamma }}^{\prime },{\sigma }^2\right)={\displaystyle \prod _{l=1}^L{\displaystyle \prod _{k=1}^K\mathcal{C}\mathcal{N}\left(a_{l,k}|0,{\gamma }_k\right)}}={\displaystyle \prod _{l=1}^L{\displaystyle \prod _{k=1}^K\mathcal{C}\mathcal{N}\left(a_{l,k}|0,{\sigma }^2{\gamma }_k^{\prime }\right)}}$$

(28)

$$p\left({\sigma }^2|a,b\right)=\mathrm{Gamma}\left({\sigma }^{-2}|a,b\right)$$

(29)

where ${\mathsf{\Gamma }}^{\prime }=\mathrm{diag}\left({\gamma }^{\prime }\right)=\mathrm{diag}\left({\gamma }_1^{\prime },{\gamma }_2^{\prime },\cdots ,{\gamma }_{2K}^{\prime }\right)=\frac{\mathsf{\Gamma }}{{\sigma }^2}$.

4.2. Application of the Modified Hierarchical Model to Complex-Valued Signals

However, although both of the original and the modified models apply to real-valued signals, they cannot be directly extended to complex-valued signals. To solve the above problem, the real and imaginary components of the complex-valued signals are treated as two independent real-valued variables.

$$y_l=\left[\begin{array}{c}\mathrm{Re}\left(x_l\right)\\ \mathrm{Im}\left(x_l\right)\end{array}\right]$$

(30)

The new signal model is expressed as

$$Y=\Psi B+N$$

(31)

where $Y=\left[y_1,y_2,\cdots ,y_L\right]\in {ℝ}^{2MN\times L}$ is the new data model, $B=\left[b_1,b_2,\cdots ,b_L\right]\in {ℝ}^{2K\times L}$ is the new sparse coefficient matrix and $N\in {ℝ}^{2MN\times L}$ is the new noise matrix in new model. The new dictionary $\Psi \in {ℝ}^{2MN\times 2K}$ is expressed as

$$\Psi =\left[\begin{array}{ccccc}\mathrm{Re}\left({\upsilon }_1\right)& -\mathrm{Im}\left({\upsilon }_1\right)& \cdots & \mathrm{Re}\left({\upsilon }_K\right)& -\mathrm{Im}\left({\upsilon }_K\right)\\ \mathrm{Im}\left({\upsilon }_1\right)& \mathrm{Re}\left({\upsilon }_1\right)& \cdots & \mathrm{Im}\left({\upsilon }_K\right)& \mathrm{Re}\left({\upsilon }_K\right)\end{array}\right]$$

(32)

while the new covariance, mean and hyper-parameter in the modified model are expressed as

$$\Sigma ={\left({\sigma }^{-2}{\Psi }^{\mathrm{T}}\Psi +{\mathsf{\Gamma }}^{-1}\right)}^{-1}$$

(33)

$$\mathsf{\mu }={\sigma }^{-2}\Sigma {\Psi }^{\mathrm{T}}y$$

(34)

$$\mathsf{\Gamma }=\mathrm{diag}\left({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_K,{\gamma }_{K+1},{\gamma }_{K+2},\cdots ,{\gamma }_{2K}\right)$$

(35)

where ${\gamma }_{2i-1}={\gamma }_{2i}, i=1,2,\cdots ,K$. The signal $y_l$ is group sparse because the entries in each group, saying ${\gamma }_{\left\{i\right\}}=\left\{{\gamma }_{2i-1},{\gamma }_{2i}\right\}, i=1,2,\cdots ,K$, are either zero or nonzero simultaneously.

In order to avoid calculating ${\sigma }^2$ in (33) and (34), we extend the modified hierarchical model to the SBL framework. The new real-valued likelihood function can be expressed as

$$\begin{array}{ll}p\left(Y|B,{\sigma }^2\right)\hfill & ={\left(2\pi {\sigma }^2\right)}^{-\frac{2MNL}{2}}\exp \left[-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{l=1}^L{‖y_l-\Psi b_l‖}^2}\right]\hfill \\ & ={\left(2\pi {\sigma }^2\right)}^{-MNL}\exp \left[-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{l=1}^L{‖y_l-\Psi b_l‖}^2}\right]\hfill \end{array}$$

(36)

and the new prior of $B$ can be expressed as

$$\begin{array}{ll}p\left(B|{\mathsf{\Gamma }}^{\prime },{\sigma }^2\right)\hfill & ={\left(2\pi \right)}^{-\frac{2KL}{2}}{\mathsf{\Gamma }}^{-\frac{L}{2}}\exp ⟮-\frac{1}{2}{\displaystyle \sum _{l=1}^L{b}_l^{\mathrm{T}}{\mathsf{\Gamma }}^{-1}{b}_l}⟯\hfill \\ & ={\left(2\pi \right)}^{-KL}{\left({\sigma }^2{\mathsf{\Gamma }}^{\prime }\right)}^{-\frac{L}{2}}\exp ⟮-\frac{1}{2}{\displaystyle \sum _{l=1}^L{b}_l^{\mathrm{T}}{\left({\sigma }^2{\mathsf{\Gamma }}^{\prime }\right)}^{-1}{b}_l}⟯\hfill \end{array}$$

(37)

Combining the new likelihood and prior, the new posterior density of $B$ can be expressed as

$$p\left(B|Y,{\mathsf{\Gamma }}^{\prime },{\sigma }^2\right)=\frac{p\left(Y|B,{\sigma }^2\right)p\left(B|{\mathsf{\Gamma }}^{\prime },{\sigma }^2\right)}{{\displaystyle \int p\left(Y|B,{\sigma }^2\right)p\left(B|{\mathsf{\Gamma }}^{\prime },{\sigma }^2\right)dB}}$$

(38)

Integrate ${\sigma }^2$ out, and the posterior density function of $B$ is

$$\begin{array}{ll}p\left(B|Y,{\mathsf{\Gamma }}^{\prime }\right)\hfill & ={\displaystyle \int p\left(B|Y,{\mathsf{\Gamma }}^{\prime },{\sigma }^2\right)p\left({\sigma }^2|a,b\right)d{\sigma }^2}\hfill \\ & =\frac{\Gamma \left(a+LK\right){\left[1+\frac{1}{2b}{\displaystyle \sum _{l=1}^L{\left(a_l-{\mathsf{\mu }}_l\right)}^{\mathrm{H}}{\Sigma }^{-1}\left(a_l-{\mathsf{\mu }}_l\right)}\right]}^{-\left(a+LK\right)}}{\Gamma \left(a\right){\left(2\pi b\right)}^{LK}{\left|\Sigma \right|}^{\frac{L}{2}}}\hfill \end{array}$$

(39)

where

$$\Sigma ={\left({\Psi }^{\mathrm{T}}\Psi +{\mathsf{\Gamma }}^{\prime -1}\right)}^{-1}$$

(40)

$${\mathsf{\mu }}_l=\Sigma {\Psi }^{\mathrm{T}}{y}_l$$

(41)

From (39), we can draw the conclusion that the modified formulation induces a heavy-tailed Student-t distribution on the residual noise, which improves the robustness of the algorithm.

The logarithm of the marginal likelihood function is

$$\begin{array}{ll}\mathcal{L}\left({\mathsf{\Gamma }}^{\prime }\right)\hfill & =\log p\left(Y|{\mathsf{\Gamma }}^{\prime }\right)\hfill \\ & =\log {\displaystyle \int p\left(Y|B,{\sigma }^2\right)p\left(B|{\mathsf{\Gamma }}^{\prime },{\sigma }^2\right)p\left({\sigma }^2|a,b\right)dBd{\sigma }^2}\hfill \\ & =-\frac{1}{2}{\displaystyle \sum _{l=1}^L\left[\left(2MN+2a\right)\log \left(y_l^{\mathrm{T}}{C}^{-1}{y}_l+2b\right)+\log \left|C\right|\right]}+\mathrm{constant}\hfill \end{array}$$

(42)

with

$$C=I_{2MN}+{\Psi \mathsf{\Gamma }}^{\prime }{\Psi }^{\mathrm{T}}$$

(43)

Since the signals $y_l, l=1\cdots ,L$ is group-sparse, we define

$${\alpha }_i={\gamma }_{2i-1}^{\prime }={\gamma }_{2i}^{\prime },i=1,2,\cdots ,K$$

(44)

and then

$${\mathsf{\Gamma }}^{\prime }=\mathrm{diag}\left(\alpha \right)\otimes I_2$$

(45)

where $\alpha =\left[{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K\right]$.

4.3. Maximization of $\mathcal{L}\left({\mathsf{\Gamma }}^{\prime }\right)$ to Estimate $\alpha$

A most probable way to point estimate hyper-parameters $\alpha$ may be found via a type-II maximum likelihood procedure. Mathematically, maximize the marginal likelihood function or its logarithm (42) with respect to $\alpha$. In the following substance, the FMLM algorithm is applied to maximize (42) to estimate $\alpha$. Unlike the EM algorithm, the FMLM algorithm reduces the computational complexities by identifying the support space of data, i.e., the atoms in the dictionary whose corresponding values in $\alpha$ is non-zero.

For notational convenience, we ignore the symbol ${ }^{\prime }$. $\gamma$ and $\mathsf{\Gamma }$ in the following substance all represent ${\gamma }^{\prime }$ and ${\mathsf{\Gamma }}^{\prime }$, respectively.

Define $\Omega$ as the set of the non-zero values in $\mathsf{\Gamma }$ and $\mathsf{\psi }$ as the support space of data.

$$\Omega =\mathrm{diag}\left({\alpha }_{{\omega }_1},{\alpha }_{{\omega }_2},\cdots ,{\alpha }_{{\omega }_J}\right)\otimes I_2$$

(46)

$$\mathsf{\psi }=\left[{\mathsf{\Phi }}_{{\omega }_1},{\mathsf{\Phi }}_{{\omega }_2},\cdots ,{\mathsf{\Phi }}_{{\omega }_J}\right]$$

(47)

where

$$\left\{{\omega }_1,{\omega }_2,\cdots ,{\omega }_J\right\}=\left\{i|1\le i\le K,{\alpha }_i\ne 0\right\}$$

(48)

$${\mathsf{\Phi }}_i=\left[\begin{array}{cc}\mathrm{Re}\left({\upsilon }_i\right)& -\mathrm{Im}\left({\upsilon }_i\right)\\ \mathrm{Im}\left({\upsilon }_i\right)& \mathrm{Re}\left({\upsilon }_i\right)\end{array}\right],i=1,2,\cdots ,K$$

(49)

and $J$ is the number of non-zeros in $\alpha$.

At the beginning of the FMLM algorithm, we initialize ${\alpha }^{\left(0\right)}=0$, namely, ${\Omega }^{\left(0\right)}=\varnothing$ and ${\mathsf{\psi }}^{\left(0\right)}=\varnothing$. Then, we need to identify the support space of data in each iteration until that $\alpha$ converges to a steady point.

The matrix $C$ can be decomposed into two parts.

$$\begin{array}{ll}C\hfill & =I_{2MN}+{\displaystyle \sum _{k\ne i,1\le k\le K}{\alpha }_k{\mathsf{\Phi }}_k{\mathsf{\Phi }}_k^{\mathrm{T}}}+{\alpha }_i{\mathsf{\Phi }}_i{\mathsf{\Phi }}_i^{\mathrm{T}}\hfill \\ & =C_{-i}+{\alpha }_i{\mathsf{\Phi }}_i{\mathsf{\Phi }}_i^{\mathrm{T}}\hfill \end{array}$$

(50)

The first term $C_{-i}=I_{2MN}+{\displaystyle \sum _{k\ne i,1\le k\le K}{\alpha }_k{\mathsf{\Phi }}_k{\mathsf{\Phi }}_k^{\mathrm{T}}}$ contains all terms that are independent of ${\alpha }_i$, and the second term includes all the terms related to ${\alpha }_i$.

Using the Woodbury Matrix Identity, the matrix inversion and matrix determinate lemmas are expressed as

$$C^{-1}=C_{-i}^{-1}-C_{-i}^{-1}{\mathsf{\Phi }}_i{\left({\alpha }_i^{-1}{I}_2+{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}{\mathsf{\Phi }}_i\right)}^{-1}{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}$$

(51)

$$\left|C\right|=\left|C_{-i}\right|\left|{\alpha }_i{I}_2\right|\left|{\alpha }_i^{-1}{I}_2+{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}{\mathsf{\Phi }}_i\right|$$

(52)

Then, (42) can be expressed as

$$\begin{array}{l}\mathcal{L}\left(\mathsf{\Gamma }\right)=-\frac{1}{2}\sum _{l=1}^L\left[\left(2MN+2a\right)log\left({\mathbf{y}}_l^{\mathrm{T}}{\mathbf{C}}_{-i}^{-1}{\mathbf{y}}_l+2b\right)+log\left|{\mathbf{C}}_{-i}\right|\right]\\ -\frac{1}{2}{\displaystyle \sum _{l=1}^L}\left[\begin{array}{cc}\hfill & \left(2MN+2a\right)log⟮1,-,\frac{{\mathbf{y}}_l^{\mathrm{T}}{\mathbf{C}}_{-i}^{-1}{\mathsf{\Phi }}_i{\left({\alpha }_i^{-1}{\mathbf{I}}_2+{\mathsf{\Phi }}_i^{\mathrm{T}}{\mathbf{C}}_{-i}^{-1}{\mathsf{\Phi }}_i\right)}^{-1}{\mathsf{\Phi }}_i^{\mathrm{T}}{\mathbf{C}}_{-i}^{-1}{\mathbf{y}}_l}{{\mathbf{y}}_l^{\mathrm{T}}{\mathbf{C}}_{-i}^{-1}{\mathbf{y}}_l+2b}⟯\\ \hfill & +log\left|{\alpha }_i{\mathbf{I}}_2\right|+log\left|{\alpha }_i^{-1}{\mathbf{I}}_2+{\mathsf{\Phi }}_i^{\mathrm{T}}{\mathbf{C}}_{-i}^{-1}{\mathsf{\Phi }}_i\right|\end{array}\right]+\mathrm{constant}\hfill \\ ={\mathcal{L}}_{-i}+\mathcal{L}\left({\alpha }_i\right)\hfill \end{array}$$

(53)

where ${\mathcal{L}}_{-i}$ contains all terms that are independent of ${\alpha }_i$ and $\mathcal{L}\left({\alpha }_i\right)$ includes all the terms related to ${\alpha }_i$.

Define

$${\widehat{S}}_i\triangleq {\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}{\mathsf{\Phi }}_i,{\widehat{q}}_{l,i}\triangleq {\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}{y}_l,{\widehat{g}}_{l,i}\triangleq y_l^{\mathrm{T}}{C}_{-i}^{-1}{y}_l+2b$$

(54)

Then, $\mathcal{L}\left({\alpha }_i\right)$ can be expressed as

$$\mathcal{L}\left({\alpha }_i\right)=-\frac{1}{2}{\displaystyle \sum _{l=1}^L\left[\begin{array}{l}\left(2MN+2a\right)\log ⟮1-\frac{{\widehat{q}}_{l,i}^{\mathrm{T}}{\left({\alpha }_i^{-1}{I}_2+{\widehat{S}}_i\right)}^{-1}{\widehat{q}}_{l,i}}{{\widehat{g}}_{l,i}}⟯\\ +2\log \left|{\alpha }_i\right|+\log \left|{\alpha }_i^{-1}{I}_2+{\widehat{S}}_i\right|\end{array}\right]}$$

(55)

The eigenvalue decomposition (EVD) of ${\widehat{S}}_i$

$${\widehat{S}}_i=V_i{\Lambda }_i{V}_i^{\mathrm{T}}={\displaystyle \sum _{j=1}^2{\widehat{s}}_{i,j}{\upsilon }_{i,j}{\upsilon }_{i,j}^{\mathrm{T}}}$$

(56)

where ${\widehat{s}}_{i,j}$ and ${\upsilon }_{i,j}$ denote the $j-\mathrm{th}$ eigenvalue and eigenvector of ${\widehat{S}}_i$ respectively. $V_i=\left[{\upsilon }_{i,1},{\upsilon }_{i,2}\right]$ is the eigen-matrix of ${\widehat{S}}_i$.

Define $c_{l,i}=V_i^{\mathrm{T}}{\widehat{q}}_{l,i}=\left[\begin{array}{c}{c}_{l,i,1}\\ c_{l,i,2}\end{array}\right]$, and the formula (55) can then be expressed as

$$\mathcal{L}\left({\alpha }_i\right)=-\frac{1}{2}{\displaystyle \sum _{l=1}^L\left[\begin{array}{l}\left(2MN+2a\right)\log \left(1-\frac{\frac{{\alpha }_i{c}_{l,i,1}^2}{1+{\alpha }_i{\widehat{s}}_{i,1}}+\frac{{\alpha }_i{c}_{l,i,2}^2}{1+{\alpha }_i{\widehat{s}}_{i,2}}}{{\widehat{g}}_{l,i}}\right)\\ +2\log \left({\alpha }_i\right)+{\displaystyle \sum _{j=1}^2\log \left(\frac{1+{\alpha }_i{\widehat{s}}_{i,j}}{{\alpha }_i}\right)}\end{array}\right]}$$

(57)

The next step involves maximizing (57) to estimate the hyper-parameters ${\alpha }_i ,\forall i$, and we choose only one candidate from ${\alpha }_i ,\forall i$ that can maximize (53).

Differentiate $\mathcal{L}\left({\alpha }_i\right)$ with respect to ${\alpha }_i$, and set the result to zero

$$\frac{\partial \mathcal{L}\left({\alpha }_i\right)}{\partial {\alpha }_i}=-\frac{1}{2}{\displaystyle \sum _{l=1}^L\left\{\left(\frac{-\left(2MN+2a\right)⟮\frac{c_{l,i,1}^2}{{\left(1+{\alpha }_i{\widehat{s}}_{i,1}\right)}^2}+\frac{c_{l,i,2}^2}{{\left(1+{\alpha }_i{\widehat{s}}_{i,2}\right)}^2}⟯}{{\widehat{g}}_{l,i}-\frac{{\alpha }_i{c}_{l,i,1}^2}{1+{\alpha }_i{\widehat{s}}_{i,1}}-\frac{{\alpha }_i{c}_{l,i,2}^2}{1+{\alpha }_i{\widehat{s}}_{i,2}}}\right)+\frac{{\widehat{s}}_{i,1}}{1+{\alpha }_i{\widehat{s}}_{i,1}}+\frac{{\widehat{s}}_{i,2}}{1+{\alpha }_i{\widehat{s}}_{i,2}}\right\}}$$

(58)

which is equivalent to solve the following polynomial function

$$f\left({\alpha }_i\right)=d_3{\alpha }_i^3+d_2{\alpha }_i^2+d_1{\alpha }_i^1+d_0=0$$

(59)

where

$$\begin{array}{cc}\hfill & d_0=\sum _{l=1}^L\left\{{\widehat{g}}_{l,i}\left({\widehat{s}}_{i,1}+{\widehat{s}}_{i,2}\right)-\left(2MN+2a\right)\left(c_{l,i,1}^2+c_{l,i,2}^2\right)\right\}\hfill \\ \hfill & d_1=\sum _{l=1}^L\left\{\begin{array}{c}4{\widehat{g}}_{l,i}{\widehat{s}}_{i,1}{\widehat{s}}_{i,2}+{\widehat{g}}_{l,i}\left({\widehat{s}}_{i,1}^2+{\widehat{s}}_{i,2}^2\right)-\left(c_{l,i,1}^2+c_{l,i,2}^2\right)\left({\widehat{s}}_{i,1}+{\widehat{s}}_{i,2}\right)\hfill \\ -2\left(2MN+2a\right)\left(c_{l,i,1}^2{\widehat{s}}_{i,2}+c_{l,i,2}^2{\widehat{s}}_{i,1}\right)\hfill \end{array}\right\}\hfill \\ & d_2=\sum _{l=1}^L\left\{\begin{array}{c}3{\widehat{g}}_{l,i}{\widehat{s}}_{i,1}^2{\widehat{s}}_{i,2}+3{\widehat{g}}_{l,i}{\widehat{s}}_{i,1}{\widehat{s}}_{i,2}^2-3c_{l,i,1}^2{\widehat{s}}_{i,1}{\widehat{s}}_{i,2}-3c_{l,i,2}^2{\widehat{s}}_{i,1}{\widehat{s}}_{i,2}-c_{l,i,1}^2{\widehat{s}}_{i,2}^2\hfill \\ -c_{l,i,2}^2{\widehat{s}}_{i,1}^2-\left(2MN+2a+1\right)c_{l,i,1}^2{\widehat{s}}_{i,2}^2-\left(2MN+2a+1\right)c_{l,i,2}^2{\widehat{s}}_{i,1}^2\hfill \end{array}\right\}\hfill \\ & d_3=\sum _{l=1}^{L}2{\widehat{g}}_{l,i}{\widehat{s}}_{i,1}^2{\widehat{s}}_{i,2}^2-2c_{l,i,1}^2{\widehat{s}}_{i,1}{\widehat{s}}_{i,2}^2-2c_{l,i,2}^2{\widehat{s}}_{i,1}^2{\widehat{s}}_{i,2}\hfill \end{array}$$

(60)

It has at most 3 roots found by standard root-finding algorithms. Considering that ${\alpha }_i\ge 0$ and the corresponding ${\mathsf{\Phi }}_i$ is not in the support space when ${\alpha }_i=0$, the set of the positive roots is

$$\Theta =\left\{{\alpha }_i|{\alpha }_i\in {ℝ}^{+},f\left({\alpha }_i\right)=0\right\}$$

(61)

The optimal ${\alpha }_i^{\ast }$ is given by

$${\alpha }_i^{\ast }=\left\{\begin{array}{ccc}0& & if\underset{{\alpha }_i\in \mathsf{\Theta }}{ \max }\mathcal{L}\left({\alpha }_i\right)<0\\ \underset{{\alpha }_i\in \mathsf{\Theta }}{\arg \max }\mathcal{L}\left({\alpha }_i\right)& & \mathrm{else}\end{array}\right.$$

(62)

If ${\alpha }_i^{\ast }=0$, it means the set $\left\{{\alpha }_i,{\mathsf{\Phi }}_i\right\}$ has no contribution to $\mathcal{L}\left(\mathsf{\Gamma }\right)$, that is to say, $\mathcal{L}\left({\alpha }_i^{\ast }\right)=0$. If $\underset{{\alpha }_i\in \mathsf{\Theta }}{\max }\mathcal{L}\left({\alpha }_i\right)<0$, we need to set ${\alpha }_i^{\ast }=0$ because $\mathcal{L}\left(0\right)=0>\underset{{\alpha }_i\in \mho }{ \max }\mathcal{L}\left({\alpha }_i\right)$.

In the $\left(t+1\right)-\mathrm{th}$ iteration, we choose only one ${\alpha }_j,1\le j\le K$ that can maximize (53), while fix the other $\left\{{\alpha }_i|1\le i\le K,i\ne j\right\}$. The sequence number $j$ is expressed as

$$j=\underset{i,1\le i\le K}{\arg \max }\mathcal{L}\left(\mathsf{\Gamma }\right)\underset{i,1\le i\le K}{=\arg \max }{\mathcal{L}}_{-i}+\mathcal{L}\left({\alpha }_i^{\ast }\right)$$

(63)

In order to avoid calculating ${\mathcal{L}}_{-i},\forall i$, we define

$$\Delta \mathcal{L}\left(i\right)\triangleq {\mathcal{L}}_{-i}+\mathcal{L}\left({\alpha }_i^{\ast }\right)-\mathcal{L}\left({\mathsf{\Gamma }}^{\left(t\right)}\right)={\mathcal{L}}_{-i}+\mathcal{L}\left({\alpha }_i^{\ast }\right)-{\mathcal{L}}_{-i}-\mathcal{L}\left({\alpha }_i^{\left(t\right)}\right)=\mathcal{L}\left({\alpha }_i^{\ast }\right)-\mathcal{L}\left({\alpha }_i^{\left(t\right)}\right)$$

(64)

Since $\mathcal{L}\left({\mathsf{\Gamma }}^{\left(t\right)}\right)$ in the $\left(t+1\right)-\mathrm{th}$ iteration is a fixed constant, the sequence number $j$ can also be expressed as

$$j=\underset{i,1\le i\le K}{\mathrm{argmax}}\Delta \mathcal{L}\left(i\right)=\underset{i,1\le i\le K}{\mathrm{argmax}}\mathcal{L}\left({\alpha }_i^{\ast }\right)-\mathcal{L}\left({\alpha }_i^{\left(t\right)}\right)$$

(65)

The next step is to change the value of ${\alpha }_j$ while fixing the other $\left\{{\alpha }_i|1\le i\le K,i\ne j\right\}$.

$${\alpha }_j^{\left(t+1\right)}={\alpha }_j^{\ast }$$

(66)

$${\alpha }_i^{\left(t+1\right)}={\alpha }_i^{\left(t\right)},1\le i\le K,i\ne j$$

(67)

If ${\alpha }_j^{\left(t+1\right)}>0$ and ${\alpha }_j^{\left(t\right)}=0$, add ${\mathsf{\Phi }}_j$ to ${\mathsf{\psi }}^{\left(t\right)}$ (i.e., ${\mathsf{\Phi }}_j$ is not in ${\mathsf{\psi }}^{\left(t\right)}$, ${\mathsf{\psi }}^{\left(t+1\right)}=[{\mathsf{\psi }}^{\left(t\right)},{\mathsf{\Phi }}_j]$, ${\Omega }^{\left(t+1\right)}=\left[\begin{array}{cc}{\Omega }^{\left(t\right)}& \\ & {\alpha }_j^{\left(t+1\right)}\otimes I_2\end{array}\right]$); if ${\alpha }_j^{\left(t+1\right)}>0$ and ${\alpha }_j^{\left(t\right)}>0$ (i.e., ${\mathsf{\Phi }}_j$ is already in ${\mathsf{\psi }}^{\left(t\right)}$, ${\mathsf{\psi }}^{\left(t+1\right)}={\mathsf{\psi }}^{\left(t\right)}$), replace ${\alpha }_j^{\left(t\right)}$ with ${\alpha }_j^{\left(t+1\right)}$ in ${\Omega }^{\left(t\right)}$; if ${\alpha }_j^{\left(t+1\right)}=0$ and ${\alpha }_j^{\left(t\right)}>0$, delete ${\mathsf{\Phi }}_j$ from ${\mathsf{\psi }}^{\left(t\right)}$ and delete ${\alpha }_j^{\left(t\right)}$ from ${\Omega }^{\left(t\right)}$; and if ${\alpha }_j^{\left(t+1\right)}=0$ and ${\alpha }_j^{\left(t\right)}=0$ (i.e., ${\mathsf{\Gamma }}^{\left(t+1\right)}={\mathsf{\Gamma }}^{\left(t\right)}$), stop iteration because the procedure converges to steady state. Finally, ${\mathsf{\psi }}^{*}$ is the exact support space of the data and ${\Omega }^{\ast }$ is the set of the non-zero values in the exact $\mathsf{\Gamma }$ where the symbol * represents the stopping criterion.

4.4. Fast Computation of $\left\{{\widehat{S}}_i,{\widehat{q}}_{l,i},{\widehat{g}}_{l,i}\right\}$

In each iteration, we need to calculate the matrix inversion of all $C_{-i}, \forall i$ when updating $\left\{{\widehat{S}}_i,{\widehat{q}}_{l,i},{\widehat{g}}_{l,i}\right\}$. In order to reduce computation complexities, we need a fast means to update $\left\{{\widehat{S}}_i,{\widehat{q}}_{l,i},{\widehat{g}}_{l,i}\right\}$.

Define

$$S_i\triangleq {\mathsf{\Phi }}_i^{\mathrm{T}}{C}^{-1}{\mathsf{\Phi }}_i,q_{l,i}\triangleq {\mathsf{\Phi }}_i^{\mathrm{T}}{C}^{-1}{y}_l,g_l\triangleq y_l^{\mathrm{T}}{C}^{-1}{y}_l+2b$$

(68)

If we can calculate $\left\{{\widehat{S}}_i,{\widehat{q}}_{l,i},{\widehat{g}}_{l,i}\right\}$ with $\left\{S_i,q_{l,i},g_l\right\}$, the computational complexities can be reduced a lot because we only calculate the matrix inversion of $C$.

Substituting (51) into (68), we can arrive at the following formula:

$$\begin{array}{ll}{S}_i\hfill & ={\mathsf{\Phi }}_i^{\mathrm{T}}{C}^{-1}{\mathsf{\Phi }}_i\hfill \\ & ={\mathsf{\Phi }}_i^{\mathrm{T}}\left(C_{-i}^{-1}-C_{-i}^{-1}{\mathsf{\Phi }}_i{\left({\alpha }_i^{-1}{I}_2+{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}{\mathsf{\Phi }}_i\right)}^{-1}{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}\right){\mathsf{\Phi }}_i\hfill \\ & ={\widehat{S}}_i-{\widehat{S}}_i{\left({\alpha }_i^{-1}{I}_2+{\widehat{S}}_i\right)}^{-1}{\widehat{S}}_i\hfill \\ & =V_i\left[\begin{array}{cc}\frac{{\widehat{S}}_{i,1}}{1+{\alpha }_i{\widehat{S}}_{i,1}}& \\ & \frac{{\widehat{S}}_{i,2}}{1+{\alpha }_i{\widehat{S}}_{i,2}}\end{array}\right]V_i^{\mathrm{T}}\hfill \end{array}$$

(69)

We can draw the conclusion that the eigen-matrix of $S_i$ is the same as that of ${\widehat{S}}_i$ from (56) and (69). The EVD of $S_i$ can be also expressed as

$$S_i=V_i\left[\begin{array}{cc}{S}_{i,1}& \\ & S_{i,2}\end{array}\right]V_i^{\mathrm{T}}$$

(70)

where $S_{i,j}$ denotes the $j-\mathrm{th}$ eigenvalue of $S_i$.

$$S_{i,j}=\frac{{\widehat{S}}_{i,j}}{1+{\alpha }_i{\widehat{S}}_{i,j}},j=1,2$$

(71)

Thus, ${\widehat{S}}_i$ can be obtained via the EVD of $S_i$

$${\widehat{S}}_i=V_i\left[\begin{array}{cc}\frac{S_{i,1}}{1-{\alpha }_i{S}_{i,1}}& \\ & \frac{S_{i,2}}{1-{\alpha }_i{S}_{i,2}}\end{array}\right]V_i^{\mathrm{T}}$$

(72)

${\widehat{q}}_{l,i}$ and ${\widehat{g}}_{l,i}$ can be also computed with $q_{l,i}$ and $g_l$, respectively, via the EVD of $S_i$.

$$\begin{array}{ll}{q}_{l,i}\hfill & ={\mathsf{\Phi }}_i^{\mathrm{T}}{C}^{-1}{y}_l\hfill \\ & ={\mathsf{\Phi }}_i^{\mathrm{T}}\left(C_{-i}^{-1}-C_{-i}^{-1}{\mathsf{\Phi }}_i{\left({\alpha }_i^{-1}{I}_2+{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}{\mathsf{\Phi }}_i\right)}^{-1}{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}\right)y_l\hfill \\ & ={\widehat{q}}_{l,i}-{\widehat{S}}_i{\left({\alpha }_i^{-1}{I}_2+{\widehat{S}}_i\right)}^{-1}{\widehat{q}}_{l,i}\hfill \\ & =\left(I_2-{\widehat{S}}_i{\left({\alpha }_i^{-1}{I}_2+{\widehat{S}}_i\right)}^{-1}\right){\widehat{q}}_{l,i}\hfill \end{array}$$

(73)

Thus,

$${\widehat{q}}_{l,i}={\left(I_2-{\widehat{S}}_i{\left({\alpha }_i^{-1}{I}_2+{\widehat{S}}_i\right)}^{-1}\right)}^{-1}{q}_{l,i}=V_i\left[\begin{array}{cc}\frac{1}{1-{\alpha }_i{s}_{i,1}}& \\ & \frac{1}{1-{\alpha }_i{s}_{i,2}}\end{array}\right]V_i^{\mathrm{T}}{q}_{l,i}$$

(74)

Similarly,

$$\begin{array}{ll}{g}_l\hfill & =y_i^{\mathrm{T}}\left(C_{-i}^{-1}-C_{-i}^{-1}{\mathsf{\Phi }}_i{\left({\alpha }_i^{-1}{I}_2+{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}{\mathsf{\Phi }}_i\right)}^{-1}{\mathsf{\Phi }}_i^{\mathrm{T}}{C}_{-i}^{-1}\right)y_l+2b\hfill \\ & ={\widehat{g}}_{l,i}-{\widehat{q}}_{l,i}^{\mathrm{T}}{\left({\alpha }_i^{-1}{I}_2+{\widehat{S}}_i\right)}^{-1}{\widehat{q}}_{l,i}\hfill \\ & ={\widehat{g}}_{l,i}-{\widehat{q}}_{l,i}^{\mathrm{T}}{V}_i\left[\begin{array}{cc}{\alpha }_i\left(1-{\alpha }_i{S}_{i,1}\right)& \\ & {\alpha }_i\left(1-{\alpha }_i{S}_{i,2}\right)\end{array}\right]V_i^{\mathrm{T}}{\widehat{q}}_{l,i}\hfill \end{array}$$

(75)

Using (74), we can obtain

$${\widehat{g}}_{l,i}=g_l+q_{l,i}^{\mathrm{T}}{V}_i\left[\begin{array}{cc}\frac{{\alpha }_i}{1-{\alpha }_i{s}_{i,1}}& \\ & \frac{{\alpha }_i}{1-{\alpha }_i{s}_{i,2}}\end{array}\right]V_i^{\mathrm{T}}{q}_{l,i}$$

(76)

Therefore, $\left\{{\widehat{S}}_i,{\widehat{q}}_{l,i},{\widehat{g}}_{l,i}\right\}$ can be obtained from $\left\{S_i,q_{l,i},g_l\right\}$ using (72), (74) and (76).

Compared with (68), there is another approach that requires fewer computational complexities to calculate $\left\{S_i,q_{l,i},g_l\right\}$. The formula (43) is equivalent to

$$C=I_{2MN}+\mathsf{\psi }\Omega {\mathsf{\psi }}^{\mathrm{T}}$$

(77)

With matrix inversion lemmas and (77), it is fast and more convenient to update $\left\{S_i,q_{l,i},g_l\right\}$ with the following formulae compared with (68).

$$\begin{array}{l}{S}_i={\mathsf{\Phi }}_i^{\mathrm{T}}{\mathsf{\Phi }}_i-{\mathsf{\Phi }}_i^{\mathrm{T}}\mathsf{\psi }\Sigma {\mathsf{\psi }}^{\mathrm{T}}{\mathsf{\Phi }}_i\\ q_{l,i}={\mathsf{\Phi }}_i^{\mathrm{T}}{y}_l-{\mathsf{\Phi }}_i^{\mathrm{T}}\mathsf{\psi }\Sigma {\mathsf{\psi }}^{\mathrm{T}}{y}_l\\ g_l=y_l^{\mathrm{T}}{y}_l-y_l^{\mathrm{T}}\mathsf{\psi }\Sigma {\mathsf{\psi }}^{\mathrm{T}}{y}_l+2b\end{array}$$

(78)

where $\Sigma$ herein represents the covariance of the non-zeros in $\mathsf{\Gamma }$.

$$\Sigma ={\left({\mathsf{\psi }}^{\mathrm{T}}\mathsf{\psi }+{\left\{\Omega \right\}}^{-1}\right)}^{-1}$$

(79)

Additionally, the mean of the non-zeros in $\mathsf{\Gamma }$ herein is expressed as

$${\mathsf{\mu }}_l=\Sigma {\mathsf{\psi }}^{\mathrm{T}}{y}_l,l=1,\cdots ,L$$

(80)

The computational complexities are measured in terms of the number of multiplications. Assuming that the dimension of $\mathsf{\psi }$ is $2MN\times r$, $r$ is not a fixed value in each iteration but satisfies the condition $r\ll MN<K$ because the measurements are sparse. When we calculate $S_i$ with (68), the computation complexities are in order of $o\left({\left(2MN\right)}^3+8{\left(MN\right)}^2+8MN\right)$. With (78), the computation complexities of calculating $S_i$ are in order of $o\left(r^3+16MN+8MNr+2r^2\right)$. It is apparent that the latter operation is faster and more convenient.

Since the measurements are sparse, the dimension of the matrix ${\mathsf{\psi }}^{\mathrm{T}}\mathsf{\psi }+{\left\{\Omega \right\}}^{-1}$ in (79) is far smaller than the dimension of ${\Psi }^{\mathrm{T}}\Psi +{\mathsf{\Gamma }}^{\prime -1}$ in (41). The proposed M-FMLM-STAP algorithm identifies the support space of data to reduce the effective problem dimensions drastically due to sparsity. Therefore, the proposed algorithm has tremendous potential for real- time operation.

The proposed M-FMLM-STAP algorithm is shown in Algorithm 2. The detailed update formulae are shown in the Appendix A.

Algorithm 2: Pseudocode for M-FMLM-STAP algorithm.

Step 1: Input: the original data X, the original dictionary $D$ and $a=b=0$.

Step 2: $y=\left[\mathrm{Re}\left(X\right);\mathrm{Im}\left(X\right)\right]$ and $\Psi$.

Step 3: Initialize:$C=C_{-i}=I$ $\alpha =0$ and $\mathsf{\psi }=\varnothing$.

Step 4: while not converged do     Choose only one candidate $j\in \left\{1,2,\cdots ,K\right\}$ and find optimal ${\alpha }_j$ using (65).     If ${\alpha }_j^{\left(t+1\right)}>0$ and ${\alpha }_j^{\left(t\right)}=0$, then       ${\mathsf{\psi }}^{\left(t+1\right)}=[{\mathsf{\psi }}^{\left(t\right)},{\mathsf{\Phi }}_j]$, and ${\Omega }^{\left(t+1\right)}=\left[\begin{array}{cc}{\Omega }^{\left(t\right)}& \\ & {\alpha }_j^{\left(t+1\right)}\otimes I_2\end{array}\right]$.     Otherwise, if ${\alpha }_j^{\left(t+1\right)}>0$ and ${\alpha }_j^{\left(t\right)}>0$, then       ${\mathsf{\psi }}^{\left(t+1\right)}={\mathsf{\psi }}^{\left(t\right)}$, and replace ${\alpha }_j^{\left(t\right)}$ with ${\alpha }_j^{\left(t+1\right)}$.     Otherwise, if ${\alpha }_j^{\left(t+1\right)}>0$ and ${\alpha }_j^{\left(t\right)}>0$, then       delete ${\mathsf{\Phi }}_j$ from ${\mathsf{\psi }}^{\left(t\right)}$, and delete ${\alpha }_j^{\left(t\right)}$ from ${\Omega }^{\left(t\right)}$.     end     Update $\Sigma$,$\mathsf{\mu }$,$S_i$,$q_{l,i}$ and $g_l$ referring to Appendix A.    end while

Step 5: Estimate the CNCM $R$ by     $R=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{i=1}^J\left({‖{\mathsf{\mu }}_{l,2i-1}^{\ast }‖}^2+{‖{\mathsf{\mu }}_{l,2i}^{\ast }‖}^2\right){\upsilon }_{{\omega }_i}{\upsilon }_{{\omega }_i}^{\mathrm{H}}}}+\beta {\rm I}$   where the vector ${\upsilon }_{{\omega }_i}$ is the ${\omega }_i-\mathrm{th}$ column of $D$, $J$ is the number of non-zeros in ${\alpha }^{\ast }$ and $\beta$ is a load factor. The symbol * represents the stopping criterion.

Step 6: Compute the space-time adaptive weight $w$ using (7).

Step 7: The output of M-FMLM-STAP is $Z=w^{\mathrm{H}}X$.

The main symbols aforementioned are listed in

Table 1. The main symbols.

Symbols	Parameters	Symbols	Parameters
$X$	The original data	$Y$	The new data
$D$	The original dictionary	$\Psi$	The new dictionary
$A$	The original coefficient matrix	$B$	The new coefficient matrix
$\mathsf{\Gamma }$	The original hyper-parameter	${\mathsf{\Gamma }}^{\prime }$	The new hyper-parameter
$\Sigma$	The covariance of coefficient	$\mathsf{\mu }$	The mean of coefficient
$\Omega$	The set of the non-zero values in ${\mathsf{\Gamma }}^{\prime }$	$\mathsf{\psi }$	The support space of data
$\left\{{\widehat{S}}_i,{\widehat{q}}_{l,i},{\widehat{g}}_{l,i}\right\}$	See (54)	$\left\{S_i,q_{l,i},g_l\right\}$	See (68)
${\mathsf{\Phi }}_i$	See (49)

.

5. Complexity Analysis and Convergence Analysis

5.1. Complexity Analysis

Here, we compare the computational complexities of the proposed M-FMLM-STAP algorithm and the M-SBL-STAP algorithm for a single iteration. The computational complexities are measured in terms of the number of multiplications.

First, we analyze the computational complexities of the M-SBL-STAP algorithm. It is apparent that the main computational complexities of the M-SBL-STAP algorithm are related to the formulae (17) and (18). Noting that $\Lambda$ in (17) is a diagonal matrix, the computational complexities of (17) are in the order of $o\left(K^3+K^{2}MN\right)$. The computational complexities of (18) are in the order of $o\left(K^{2}MN+KMNL\right)$. Thus, the computational complexities of the M-SBL-STAP algorithm are in the order of $o\left(K^3+2K^{2}MN+KMNL\right)$.

Second, we analyze the computational complexities of the M-FMLM-STAP algorithm. Noting that the dimension of $S_i$ is $2\times 2$, the computational complexities of the EVD of $S_i$ are small enough to be ignored, that is to say, the computational complexities to calculate $\left\{{\widehat{S}}_i,{\widehat{q}}_{l,i},{\widehat{g}}_{l,i}\right\}$ with $\left\{S_i,q_{l,i},g_l\right\}$ can be ignored. Meanwhile, the computational complexities to solve the polynomial function (59), $\forall i,1\le i\le K$ and to find the sequence number (65) are in direct proportion to $K$, which are also small enough to be ignored. Thus, the main computational complexities of the M-FMLM-STAP algorithm are used to update $\left\{S_i,q_{l,i},g_l,\Sigma ,\mathsf{\mu }\right\}$. Assuming that the dimension of $\mathsf{\psi }$ is $2MN\times r$, $r$ is not a fixed value in all iterations. However, the condition $r\ll MN<K$ is satisfied because the measurements are sparse. The computational complexities of the term ${\mathsf{\Phi }}_i^{\mathrm{T}}-{\mathsf{\Phi }}_i^{\mathrm{T}}\mathsf{\psi }\Sigma {\mathsf{\psi }}^{\mathrm{T}}$ are in the order of $o\left(8MNr+2r^2\right)$. Thus, $\forall i$ and $\forall l$, the computational complexities of (78), are in the order of $o\left(\left(8MNr+2r^2\right)K+8MNK+4MNLK+\left(4MNr+r^2+4MN\right)L\right)$. The computational complexities of (79) and (80) are in the order of $o\left(2MN{r}^2+r+r^3\right)$ and $o\left(2MN{r}^2+2MNrL\right)$, respectively. To sum up, ignore the low-order computational complexities, and the computational complexities of the M-FMLM-STAP algorithm for a single iteration are in the order of $o\left(8MNrK+4MNLK+6MNrL\right)$.

Figure 1 illustrates the complexity requirements of two algorithms for a single iteration. It shows the computational complexities to the number of pluses $M$ for a single iteration. Suppose that ${\phi }_S={\phi }_d=4$, $N=10$ and $L=6$. Although the value $r$ is unknown and unfixed in each iteration, we can suppose that it is twice as much to the rank of clutter in the current iteration. We can draw the conclusion that the computational complexities of the M-SBL-STAP algorithm are far more than that of the M-FMLM-STAP algorithm for a single iteration.

5.2. Convergence Analysis

It is obvious that the logarithm of the marginal likelihood function $\mathcal{L}\left(\mathsf{\Gamma }\right)$ has upper bound [31]. From (62), in the $t-\mathrm{th}$ iteration, $\mathcal{L}\left({\alpha }_i^{\ast }\right)$ is the maximal value of $\mathcal{L}\left({\alpha }_i\right)$, that is to say, $\mathcal{L}\left({\alpha }_i^{\ast }\right)\ge \mathcal{L}\left({\alpha }_i^{\left(t\right)}\right)$. Therefore, the condition $\mathcal{L}\left({\mathsf{\Gamma }}^{\left(t+1\right)}\right)\ge \mathcal{L}\left({\mathsf{\Gamma }}^{\left(t\right)}\right),\forall t$ is promised in each iteration referring to (65). Thus, the convergence performance of our proposed algorithm is guaranteed.

6. Performance Assessment

In this section, we firstly verify the performance of the loading sample matrix inversion (LSMI) algorithm, the multiple orthogonal matching pursuit (M-OMP-STAP) algorithm [32], the M-SBL-STAP algorithm and the proposed M-FMLM-STAP algorithm with the simulated data of side-looking ULA in the ideal case and the array errors case. The first two algorithms are common approaches in STAP. We then assess the performance of the M-SBL-STAP algorithm and the proposed M-FMLM-STAP algorithm with the Mountain-Top data. We utilize the improvement factor (IF) and signal to interference plus noise ratio (SINR) loss as two measurements of performance, which are expressed as

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

(81)

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

(82)

where $R$ is the exact CNCM of the CUT.

6.1. Simulated Data

The parameters of a side-looking phased array radar are listed in

Table 2. Radar system parameters.

Symbols	Parameters	Value
$\lambda$	Wavelength	0.3 m
$d$	Distance between elements	0.15 m
$\vartheta$	Platform velocity	150 m/s
$H$	Platform height	9000 m
$M$	Number of pulses	8
$N$	Number of channels	8
$f_{\mathrm{PRF}}$	Pulse repetition frequency	2000 Hz
$f_S$	Range sampling frequency	2.5 MHz
$\theta$	Perspective angle	90°
$\mathrm{CNR}$	Clutter to noise ratio	30 dB

. The 600th range gate is chosen to be CUT. According to the parameters of the radar system, the slope of the clutter ridge is 1. We set the resolution scales ${\phi }_S=5$ and ${\phi }_d=5$. Total of six training samples are utilized for two algorithms and the parameters $a=b=0$.

In Figure 2, there are five clutter power spectrum figures in the framework of the side-looking ULA in the ideal case. Figure 2a–e show the clutter power spectrums estimated by the exact CNCM, the LSMI algorithm, the M-OMP-STAP algorithm, the M-SBL-STAP algorithm and the M-FMLM-STAP algorithm, respectively. We note that the clutter spectrums obtained by the M-SBL-STAP algorithm and the M-FMLM-STAP algorithms are closer to the exact spectrum both on the location and power than the other algorithms in the ideal case. As shown in Figure 3, the IF factor curve achieved by the M-SBL-STAP algorithm is a litter better than that achieved by the M-FMLM-STAP algorithm because the latter algorithm is sensitive to noise fluctuation. However, the latter algorithm needs much fewer computational complexities, and the loss on performance can be offset by the decrease on the computational complexities. Using a computer equipped with two Inter(R) Xeon(R) Gold 6140 CPU @2.30 GHz 2.29 GHz processors, the former spends more than 40 s and the latter spends less than 4 s at average. The larger the value of $K$, the bigger gap of computational complexities between the M-SBL-STAP algorithm and the M-FMLM-STAP algorithm.

In Figure 4, we apply four approaches in the non-ideal case with amplitude Gaussian error (standard deviation 0.03) and phase random error (standard deviation $2^{\circ }$). The error is the same at all direction. Figure 4a–e show the clutter power spectrums estimated by the exact CNCM, the LSMI algorithm, the M-OMP-STAP algorithm, the M-SBL-STAP algorithm and the M-FMLM-STAP algorithm, respectively. We note that the clutter spectrums obtained by the M-SBL-STAP algorithm and the M-FMLM-STAP algorithm are much closer to the exact spectrum, both in terms of the location and power in the non-ideal case. As shown in Figure 5, the notch of the IF factor curve achieved by the M-SBL-STAP algorithm is close to that achieved by the M-FMLM-STAP algorithm, which means the performance of two algorithms is about the same in the non-ideal case. However, the running time of the latter is much less than that of the former.

In Figure 6, the average SINR loss, defined as the mean of the SINR loss values in the entire normalized Doppler frequency range to the number of training samples of the M-SBL-STAP algorithm and the M-FMLM-STAP algorithm, is presented. The curve of the average SINR loss achieved by the M-SBL-STAP algorithm is a litter better than that achieved by the M-FMLM-STAP algorithm within 0.5 dB, which can be ignored in practice. However, the gap on computational complexities of two algorithms will be larger as the DOF of the radar system increases. All presented results are averaged over 100 independent Monte Carlo runs.

6.2. Measured Data

We apply the M-SBL-STAP algorithm and the M-FMLM-STAP algorithm to the public available Mountain-Top set, i.e., t38pre01v1 CPI6 [33]. In this data file, the number of array elements and coherent pulses are 14 and 16, respectively. There are 403 range cells, and the clutter locates around 245° relative to the true North. The target is located in 147th range cell, and the azimuth angle is 275° to the true North. The normalized target Doppler frequency is 0.25. The estimated clutter Capon spectrum utilizing all 403 range cells is provided in Figure 7.

Figure 8 depicts the STAP output power of the M-SBL-STAP algorithm and the M-FMLM-STAP algorithm in the range cells from 130 to 165, and 10 range cells out of 20 range cells located next to the CUT are selected as training samples. The target locates at 147th range cell and can be detected by two algorithms. Obviously, the detection performance of the proposed M-FMLM-STAP algorithm is close to that of M-SBL-STAP algorithm. However, the operation time of the former algorithm is much less than that of the latter algorithm. Therefore, the proposed algorithm is applicable in practice.

7. Conclusions

In this paper, we have extended the real-valued multitask compressive sensing technique to suppress complex-valued heterogeneous clutter for airborne radar. Unlike the conventional M-SBL-STAP algorithm, the proposed algorithm avoids the inversion of a $K\times K$ matrix at each iteration to guarantee real-time operations. We integrate the noise ${\sigma }^2$ out instead of estimating ${\sigma }^2$ to overcome the problem that we have no access to obtain the accurate value of ${\sigma }^2$. The complex-valued multitask clutter data are translated into group real-valued sparse signals in the article to suppress clutter because the modified hierarchical model is not suitable to complex-valued signals. At the end, simulation results demonstrate the high computational efficiency and the great performance of the proposed algorithm.