1. Introduction
Moving target detection is a fundamental function of radar systems for military surveillance and reconnaissance. Space-time adaptive processing (STAP) has become an effective and mature clutter suppression technique for airborne early warning phased array radar systems [
1,
2,
3]. The optimal STAP filter is constructed using the ideal clutter plus noise covariance matrix (CNCM). In practice, the ideal CNCM is unknown a priori and is typically estimated using independent and identically distributed (i.i.d.) samples around the cell under test (CUT). According to the well-known Reed–Mallet–Brennan criterion [
4], the number of i.i.d. samples with at least twice the system degree of freedom (DOF) is required to ensure a performance loss under 3 dB. Unfortunately, sufficient i.i.d. samples are often challenging to obtain due to various terrains, artificial structures and array configurations.
Over the past few decades, many STAP methods have been developed to improve clutter suppression performance in heterogeneous environments, including classical data-dependent reduced-rank (RR) methods [
5,
6,
7,
8] and data-independent reduced-dimension (RD) methods [
9,
10,
11,
12]. Although the required number of training samples is reduced to twice the reduced dimension, the requirement is still difficult to meet under severe non-stationary clutter environments for modern large-scale radar systems.
During the recent fifteen years, many sparse recovery-based STAP (SR-STAP) algorithms have been proposed by exploiting the underlying connection between the compressed sensing technique and the intrinsic sparsity prior of the clutter spectrum [
13,
14,
15,
16,
17,
18,
19,
20,
21]. With the help of sparse recovery theory, these SR-STAP algorithms have shown excellent clutter suppression performance with limited samples. Depending on the sparse recovery algorithms, SR-STAP can be generally classified into several categories, such as convex relaxation [
13,
14,
15], greedy algorithms [
16], Bayesian inference [
17,
18,
19], and other methods [
20,
21].
Convex relaxation substitutes the
norm with the
norm as the sparse penalty. The
norm optimization problem has shown that sparse solutions can be stably obtained under certain conditions [
22]. It has been extensively applied to the least absolute shrinkage and selection operator (LASSO) and basis pursuit (BP) based STAP methods [
23]. Most convex relaxation methods require careful tuning of the regularization parameter, and inappropriate parameter selection will jeopardize the performance of clutter suppression as well as slow-moving target detection [
14,
15]. However, choosing an appropriate regularization parameter is quite challenging in practice.
In recent years, sparse Bayesian learning (SBL) has drawn much effort due to its preferable advantages, such as automatic self-regulation [
18] and flexibility in exploiting the potential signal structure [
24]. SBL was first proposed by Tipping in 2001 [
25] and introduced to the field of STAP by Duan in 2017, termed M-SBL-STAP [
17]. Numerous empirical results indicate that the SBL based SR-STAP can provide satisfactory performance and is quite robust to noise and high coherence dictionary [
19]. However, M-SBL-STAP faces an overwhelming computational burden and large memory requirements, hindering its implementation in large-scale radar systems. Many efficient SBL based STAP algorithms have been developed to address this issue. In [
18], Wang proposed a fast-converging SBL algorithm by combining an approximation term, but the global convergence property of the algorithm is not guaranteed. An iterative reweighted based M-SBL (M-SBL-IR
) STAP algorithm was proposed by Liu [
19], which exhibits better reconstruction accuracy and has a favourable convergence speed.
Our experience with numerous simulation experiments demonstrates that the majority of the space-time coefficients to be recovered are zero or close to zero, and only a few have nontrivial values, owing to the sparse nature of the clutter. Therefore, in this paper, we are inspired to propose an improved Bayesian probabilistic model for SR-STAP, exploiting the aforementioned sparsity feature of the clutter. Instead of assigning a separate hyperparameter to each space-time coefficient in the conventional Bayesian model, only the significant coefficients are assigned independent hyperparameters, while the remaining irrelevant coefficients share a common hyperparameter. As a result, the parameter space to be updated will be dramatically reduced, and the computational efficiency will naturally be promoted. In [
26], a heuristic expansion-compression variance-component based method (ExCoV) has been proposed to guide us to classify those hyperparameters.
Specifically, the main contributions of this paper are summarized as follows:
By exploiting the inherent sparsity nature of the clutter, the space-time coefficients are divided into two disjoint groups, i.e., the significant and irrelevant groups, thus reducing redundancy, scaling down the parameter space and yielding an improved Bayesian probabilistic model for SR-STAP with reduced computational complexity and reduced memory requirements.
The space-time coefficients are partitioned into the significant and irrelevant groups according to the generalized maximum likelihood (GML) criterion, which preserves both accuracy and efficiency, unlike the conventional SBL cost function.
We extend and modify the real-value ExCoV method to complex-value STAP applications to approximately maximize the GML objective function. Using the ExCoV scheme, it is unnecessary to specify convergence thresholds and maximum iteration times.
Extensive experiments as well as detailed comparative analyses are presented, such as clutter suppression performance and target detection performance, etc.
Notations used in this paper are as follows: vectors, matrices and scalars are denoted by bold lowercase, bold uppercase and italic letters, respectively. , and stand for conjugate, transpose and conjugate transpose. and are the Kronecker and Hadamard (elementwise) product. is the trace operator. represents the set of complex values. and are respectively defined as the Frobenius norm and mixed norm, which is the number of non-zero elements of the vector formed by the norm of each row. denotes the expectation operator.
The remainder of the paper is organized as follows. In
Section 2, the airborne radar signal model is established. The formulation of SR-STAP model is presented in
Section 3. The proposed algorithm is introduced in
Section 4. The computational complexity analysis is presented in
Section 5. Experiments and analyses are carried out on
Section 6. Finally,
Section 7 discuss the conclusions.
4. Proposed Algorithm
Following the conventional sparse Bayesian learning framework, all the unknowns are treated as stochastic variables with assigned independent probability distributions. First of all, the noise matrix is modelled as white complex Gaussian noise with unknown power
; thus the observed Gaussian likelihood function of the measurements
for the MMV case can be expressed as.
Since measurements are mutually independent and identically distributed, we suppose each column of the space-time coefficient matrix
is assigned with the same zero-mean complex Gaussian prior
, governing the prior variance of each unknown space-time coefficient.
where
is a zero vector,
is the unknown variance hyperparameters corresponding to the rows in
, and covariance matrix
is a diagonal matrix, where
is its diagonal elements, i.e.,
.
Further, following the Bayesian theorem, the posterior probability density
can be easily calculated by
The posterior probability
obeys a multivariate complex Gaussian distribution
with mean and covariance respectively given by
Thus, with a fixed
, the estimated sparse recovery solution
of M-SBL is
Accordingly, for , the corresponding row of the space-time coefficient matrix will be zeros as well. In other words, if is sparse, the corresponding space-time coefficient estimation will also be sparse.
The hyperparameters vector can be estimated by performing the type-II maximum likelihood procedure or evidence maximization in
space [
29]. Mathematically, the cost function can be expressed as minimizing the marginal likelihood function
with respect to
:
After the unknown coefficients
A have been integrated out,
stands for the covariance of the measurements
with the hyperparameters
and
.
This minimization can be performed by an iterative reweighted SBL based algorithm [
19,
30], which will be modified and described in the following part of this paper.
As we can observe from (11) and (12), the computational bottleneck is mainly manifested in large-scale matrix multiplication and matrix inversion operations. The storage requirement is also heavy. Moreover, during each iteration, a hyperparameter space of dimension needs to be updated. These factors make SBL based STAP algorithms considerably slower than other types of sparse recovery algorithms, even though excellent clutter suppression performance can be achieved with limited training samples.
In fact, due to the intrinsic sparsity of the clutter spectrum on the spatial-temporal plane, only a few rows of
have nontrivial magnitudes, and the remaining elements are strictly zero (or close to zero). It implies that most of the hyperparameters in
will converge to zeros, and they are redundant. Therefore, we can naturally consider partitioning the space-time coefficients into two parts: the significant elements part and the complementary irrelevant elements part. Instead of assigning independent hyperparameters to all coefficients in the conventional SBL framework, we can only allocate independent hyperparameters to significant coefficients part, while the remaining irrelevant coefficients share an identical hyperparameter. Consequently, the dimension of the hyperparameter space can be reduced to be proportional to the sparsity level of the clutter, which is much smaller than
[
23]. As a result, the computation complexity will be greatly mitigated, and the memory consumption will be reduced, particularly in large-scale modern radar systems.
Next, define as the complete index set. The index set of significant coefficients is denoted by with size , and the complementary index set of irrelevant coefficients is denoted by with cardinality . Accordingly, we will also divide dictionary matrix and space-time coefficients matrix into two submatrices, i.e., , , and , respectively. For more details:
● is the submatrix of corresponding to the significant coefficients index set , e.g., if , then , where is the ith column of .
● is the submatrix of corresponding to the index set , e.g., if , then , where is the ith row of .
Then (9) can be reformulated as
where
and
denote the covariance matrix of the significant part
and irrelevant part
, respectively.
If all the coefficients belong to the
index set and
, then the above probabilistic model will degenerate to the original SBL full model in (9). Similarly, according to the index sets
and
, (11) and (12) can be respectively partitioned into submatrices
where the covariance matrix
in (15) now can be rewritten as
For the sake of expression simplicity, we define the set of all unknowns:
Now, let us consider estimating the most superior and efficient coefficient index sets
and
, i.e., the optimal representation. According to [
26], the generalized maximum-likelihood (GML) criterion can be used to assign the hyperparameters, and its equivalence with the canonical optimization problem in (7) is demonstrated. The GML criterion maximizes:
where the first term is the same marginal likelihood function in (14) with
fixed, enforcing the estimate fit the measurements; the second term
is the Fisher information matrix (FIM) for the hyperparameters. Based on the well-known FIM for the Gaussian measurement model [
31], we can deduce the block-partitioned FIM result for the hyperparameters vector (the detailed derivation is shown in
Appendix A).
with each block computed as
Since the hyperparameter space is nested, namely any
is a subset of
, the more parameters in the set
, the larger the value
will be, which reduces the model mismatch. Therefore, there will always be a tendency to choose the complete set
[
31]. However, the second term increases as the set
grows, penalising the growth of the set. The GML criterion is, therefore, able to maintain a balance between underfitting and overfitting of the model, ensuring both accuracy and efficiency.
Nevertheless, directly applying the GML criterion (23) to obtain the optimal parameter index sets still requires an exhaustive search and is hard to implement in practical applications. Subsequently, we will employ an ExCoV based method to maximize the objective function approximately [
32]. The fundamental idea is interleaving the expansion/compression step and the parameter update step. In each expansion or compression procedure, we modify the current estimate of
by one element per step to obtain a larger
. Then followed by the parameter update step, the hyperparameters
are updated for a fixed index set
. Subsequently, we introduce the three main iteration steps involved in ExCoV: expansion step, compression step and parameter update step.
A. Expansion step. In this step, we determine the index
corresponding to the row of
with the largest
norm
where the superscript
represents the iteration number. Then the index
q is moved from
to
, i.e.,
and
, yielding
,
. The new expanded vector
can be expressed as
.
B. Compression step. Here, we determine the index
corresponding to the smallest element of
Then the index t is moved from to , i.e., and , yielding , . The new compressed vector can be constructed as .
C. Parameters update step. In [
26], the original ExCoV algorithm employs the expectation-maximization (EM) method for Bayesian inference to estimate the hyperparameters, which usually requires extensive iterations to converge [
18]. To overcome this drawback and thus further improve code running speed, an iterative reweighted M-SBL strategy is modified to update hyperparameters.
In this step, we firstly assume that the index sets
and
are known and fixed. We next present the derivation of the parameters update step based on the iterative reweighed
M-SBL algorithm from [
19] modified by the previous proposed Bayesian probabilistic model (details can be seen in
Appendix B), which can remarkably accelerate the convergence speed. Then, the hyperparameters
can be updated by
where
.
We term the proposed algorithm as an improved iterative reweighted sparse Bayesian learning algorithm based on expansion-compression variance components (ExCoV-IIR-MSBL).
Figure 2 illustrates the processing flowchart of the proposed algorithm and its procedures are summarized as follows.
Step1 Initialization: Calculate the minimum norm solution of the space-time coefficients
Then we can utilize the prior knowledge provided by the inertial navigation system and radar system to initialize
where
represents the slope of the clutter ridge [
33] and
denotes rounding-up operation; or simply set
, since the ExCoV method is insensitive to the initial value of
. Hence, the index set
is constructed using the rows of
corresponding to the first
largest
norm. Subsequently,
and
.
Step2 Cycle Initialization: Set the initial values for the hyperparameters
.
Generally, these parameters can be initialized to an arbitrary positive vector, but it would be more beneficial if a rough estimate could be given. Then the initial values of the optimal record are set to and .
Step3 Expansion: Apply the aforementioned expansion step, yielding the updated , and .
Step4 Parameter Update: With the aforementioned modified iterative reweighted M-SBL update step, the updated parameters set and space-time coefficients estimate are obtained.
Step5 Optimal Estimates Update: During the iterations, we shall record the optimal unknown parameters set
that maximizes the GML, and the corresponding space-time coefficients estimates
and
in the latest iteration cycle. After the unknown parameters
are updated, the new
can be calculated using (23). Next, we will verify whether the following condition holds
If the inequality holds, then we set , and ; instead, keep , and unchanged.
Step6 Termination Condition Check: In this step, we would verify whether the expansion is still required to continue [
26].
where
T is the length of a sliding average window. This condition will help us to determine if there is still a need for the expansion operation and prevents premature termination. If the inequality is satisfied, then the expansion operation is terminated, and Step7 begins; otherwise, the expansion continues and returns to Step3.
Step7 Compression: Apply the aforementioned compression step, yielding the updated , and .
Step8 Parameter Update: With the aforementioned iterative reweighted M-SBL update step, the updated parameters set and space-time coefficients estimate are obtained.
Step9 Optimal Estimates Update: Verify whether the condition (38) holds. If the inequality holds, then we set , and ; instead, keep , and unchanged.
Step10 Termination Condition Check: In this step, we would verify whether the compression is still required to continue and check the same condition (39) in Step6. If the inequality is satisfied, then the compression operation is terminated, and Step11 begins; otherwise, the compression continues and returns to Step7.
Step11 Globally Optimal Estimates Update: Moreover, we keep a record of the globally optimal parameter set
and the corresponding coefficients
in the entire iteration cycles, verifying the condition
If the inequality holds, then the globally optimal estimates and can be updated with and ; instead, keep them unchanged.
Step12 Algorithm Iteration Termination Check: If the globally optimal index set is consistent between two consecutive cycles. The proposed algorithm is next terminated and the final globally optimal space-time coefficients and noise power are output.
Otherwise, the globally optimal index set
is updated, then the algorithm continues. We reset the iteration index
,
, and calculate the initial values for the space-time coefficients
Step13 STAP Weight Calculation: Once the estimated sparse recovery solution
and
is obtained, the CNCM can be reconstructed from
Then the STAP weight can be obtained by using (5).
Step14 Output: Give the filtered output of the cell under test .
From the above steps, it is worth noting that the proposed algorithm is fully automatic and does not require setting any convergence threshold as well as the maximum iteration number.
6. Numerical Experiment
In this section, numerous experiments are performed with simulated data to evaluate the performance of the proposed algorithm. The simulated data are generated by the well-known Ward clutter model introduced in the previous section. The main simulation parameters of a side-looking ULA are listed in
Table 2. The dictionary resolution scales are set to be
. Furthermore, all the simulation results are averaged over 100 independent Monte-Carlo trails.
First of all, we examine the average code running time of the proposed algorithm and compare it with other state-of-the-art SR-STAP algorithms, as listed in
Table 3. M-CVX, M-IAA, M-OMP, M-FOCUSS, M-SBL, M-FCSBL and M-SBL-IR
are employed as benchmarks. Note that all the simulations are operated on the same workstation with Intel Xeon 4114 CPU @2.2GHz and 128 GB RAM.
As seen from
Table 3, since the parameter space is dramatically shrunk, the average running time of the proposed ExCoV-IIR-SBL is remarkably faster than the other compared SBL based STAP algorithms. It is even better than M-FOCUSS and M-IAA since the proposed algorithm can reach the steady state with fewer iterations with the help of a modified iterative reweighed parameter update step. Thus, the efficiency of the proposed algorithm is demonstrated.
In the following experiments, we assess the clutter suppression performance via the metric of signal to interference plus noise ratio (SINR) loss [
27], defined as the ratio of output SINR to the optimum output SNR in a noise-only environment, i.e.,
where
is the clairvoyant CNCM of the CUT, and
is the STAP weight.
Next, the clutter suppression performance in the ideal case is considered. The SINR loss performance of the proposed algorithm is evaluated and compared with other state-of-the-art SR-STAP algorithms, including M-CVX, M-IAA, M-OMP, M-FOCUSS, M-SBL, M-FCSBL and M-SBL-IR
. The number of i.i.d. training samples are set to be
L = 9.
Figure 4 depicts the curves of the SINR loss against the normalized Doppler frequency. As clearly illustrated in
Figure 4, the proposed algorithm achieves the same near-optimal performance as M-IAA and the other three SBL based algorithms, revealing that the novel Bayesian probabilistic model proposed in this paper for STAP can accurately reconstruct the CNCM with a smaller parameter space. The performance of M-FOCUSS is found to be slightly inferior to theirs. Meanwhile, the proposed ExCoV-IIR-MSBL consumes the least running time among these three SBL based algorithms on this basis.
To more explicitly demonstrate the clutter suppression performance,
Figure 5 plots the spatial-temporal adapted responses, given the STAP weight formed by the abovementioned algorithms. The presumed main lobe is located at a normalized spatial frequency of 0 with a normalized Doppler frequency of 0.2. From
Figure 5, it is evident that the proposed method is able to maintain high gain at the presumed target location while precisely forming notches at the clutter ridge, suppressing both main lobe and sidelobe of the clutter component. Poorly performing algorithms have distorted two-dimensional responses and cannot accurately form notches at the clutter ridge.
Note that because of the massive computational complexity of the M-CVX, it will not be discussed in the subsequent experiments. Then,
Figure 6 illustrates the curves of the average SINR loss versus different a number of training samples. It should be noted that the average SINR loss is calculated by the mean of SINR loss for
in this paper. As demonstrated in
Figure 6, all algorithms have a certain degree of performance loss with one snapshot and increase with the number of training samples. ExCoV-IIR-MSBL apparently exhibits promising performance similar to the M-IAA and other SBL algorithms under small sample conditions. The proposed algorithm achieves a near-optimum performance with merely three training samples.
Subsequently, we verify the target detection performance via the probability of detection (
Pd) versus the target SNR curves, which are acquired by employing the cell-average constant false alarm rate (CA-CFAR) detector [
37]. One hundred targets within the boresight are randomly added to the entire Range-Doppler plane, and the false alarm probability (
Pfa) is fixed to 10
−3.
Figure 7 shows that the proposed algorithm and other SBL based algorithms have noticeable improvements compared to M-OMP and M-FOCUSS. This result indicates that the proposed algorithm has superior target detection performance.
Besides, we also consider the SINR loss performance under the spatial mismatch case in the presence of angle-dependent gain-phase error. Specifically, the gain errors among the antennas are generated by a Gaussian distribution with a standard deviation of 0.03, and the phase errors are uniformly distributed within [0°, 2°].
Figure 8 shows the SINR loss versus normalized Doppler frequency. From
Figure 8, the curves show that all the SR-STAP algorithms suffer severe performance degradation due to the model mismatch between the dictionary and the actual clutter component. The M-IAA has the most serious performance loss, despite the fact that it can approach near-optimal as the SBL based algorithms in the ideal case.
Figure 9 shows the SINR loss versus the number of samples of above referred algorithms in the non-ideal case. As depicted in
Figure 9, the proposed algorithm can yield a steady-state with only a few training samples as in the ideal case.
From the above experiments, it has been demonstrated that the proposed algorithm ensures a satisfactory clutter suppression and target detection performance while keeping a remarkable computational efficiency, thereby making it more suitable for implementation in modern large-scale radar systems.