- Email: [email protected]
Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar
Author’s Accepted Manuscript Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar Zetao Wang, Wenchong Xie, Keqing Duan, Yongliang Wang www.elsevier.com/locate/sigpro
PII: DOI: Reference:
S0165-1684(16)30134-7 http://dx.doi.org/10.1016/j.sigpro.2016.06.023 SIGPRO6183
To appear in: Signal Processing Received date: 29 January 2016 Revised date: 21 June 2016 Accepted date: 23 June 2016 Cite this article as: Zetao Wang, Wenchong Xie, Keqing Duan and Yongliang Wang, Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar, Signal Processing, http://dx.doi.org/10.1016/j.sigpro.2016.06.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar Zetao Wanga*, Wenchong Xieb, Keqing Duanb, Yongliang Wangb a
College of Electronic Science and Engineering, NUDT, Changsha 410073, China
b
Wuhan Radar Academy, Wuhan 430019, China
[email protected] [email protected] [email protected] [email protected] *
Corresponding author. Tel.: +86 2785965760.
Abstract Adapting the space-time adaptive processing (STAP) filter with finite number of secondary data is of particular interest for airborne phased-array radar clutter suppression. Sparse representation (SR) technique has been introduced into the STAP framework for the benefit of drastically reduced training requirement. However, most SR algorithms need the fine tuning of one or more user parameters, which affect the final results significantly. Sparse Bayesian learning (SBL) and multiple sparse Bayesian learning (M-SBL) are robust and user parameter free approaches, but they converge quite slowly. To remedy this limitation, a fast converging SBL (FCSBL) approach is proposed based on Bayesian inference along with a simple approximation term, then, it is extended to the multiple measurement vector case, and the resulting approach is termed as M-FCSBL. To improve the performance of STAP in finite secondary data situation, the M-FCSBL is utilized to estimate the clutter plus noise covariance matrix (CCM) from a limited number of secondary data, and then the resulting CCM is adopted to devise the STAP filter and suppress the clutter. Numerical experiments with both simulated and Mountain-Top data are carried out. It is shown that the proposed algorithm has superior clutter suppression performance in finite secondary data situation. Keywords: space-time adaptive processing, clutter suppression, sparse representation, Bayesian 1
inference. 1. Introduction Space-time adaptive processing (STAP) holds the inherent capacity of detecting slowly moving targets in a strong clutter background for airborne phased-array radar. The optimum STAP weight vector is formulated as the product of the inverse of clutter plus noise covariance matrix (CCM) and the target space-time steering vector [1]. When testing one range cell for the existence of a target, the CCM is usually estimated from secondary/training data that are selected from the range cells adjacent to the cell under test (CUT) [2]. In such a regime, the secondary data are assumed to be target-free and contain interference with the same statistical characteristics as the CUT. However, as a result of both geometric configuration and terrain variations, the qualified secondary data may be rather limited [3]. And the performance of STAP degrades drastically due to the inaccuracy of the estimated CCM. To increase the performance of STAP in finite secondary data situation, many different types of algorithms have been developed. Earlier achievements include reduced-dimension (RD) and reduced-rank (RR) algorithms which reduce the training requirement to twice of the reduced-dimension and twice of the clutter rank, respectively [4]. Recent developments include the parametric adaptive matched filter (PAMF) [5], direct data domain (D3) algorithms [6], and knowledge-aided (KA) algorithms [7, 8]. Although the training requirements are further reduced for these algorithms, they suffer from several drawbacks. The order of the autoregressive model for the PAMF is difficult to be determined, the system degrees of freedom (DOFs) are significantly reduced for the D3 algorithms, and the performances of the KA algorithms depend heavily on the accuracy of a priori knowledge. In most recent years, sparse representation/recovery (SR) technique has received continuous attention fueled by advanced theoretical developments and practical applications [9, 10]. It provides a general signal model that represents or approximates a signal under an overcomplete dictionary with a sparse coefficient vector. The dictionary in the context of STAP is a collection of space-time steering vectors each associated with a sampling point in the normalized space-time plane. And “sparse” means that 2
most of the elements in the coefficient vector are equal to zero. The approach of applying SR technique to STAP has been well studied in [11] and it was pointed out that the approach tries to recover the clutter subspace using as few space-time steering vectors as possible. STAP algorithms based on SR technique (SR-STAP) can be summarized as two-step procedures: first, recovering the clutter space-time profile by a certain SR algorithm; second, constructing the CCM with the result obtained in the first step and calculating the STAP weight vector [12-19]. For example, the algorithm proposed in [19] first estimates the clutter profile in the CUT by considering only the clutter entries within the learned clutter support (the clutter support is estimated from secondary data by SR technique), and then devise the STAP filter based on the estimated clutter profile. The most significant benefit brought by SR-STAP is that the training requirement can be drastically reduced without incorporating any a priori knowledge of clutter. As can be inferred from the two-step scheme of SR-STAP, the final clutter suppression performance depends heavily on the recovery result which pertains to the SR algorithm adopted. Unfortunately, most SR algorithms need the fine tuning of one or more user parameters which affect the recovery results significantly. For example, the Homotopy-STAP algorithm proposed in [17] requires the fine tuning of the regularization parameter, and the final clutter suppression performance is affected by this parameter. Bayesian type algorithms [20-26] have developed the reputation of favorable performance when the dictionary is of a high coherence. Among them, the sparse Bayesian learning (SBL) [20, 21], which was proposed for the single measurement vector (SMV) case, and its counterpart for the multiple measurement vector (MMV) case, i.e., M-SBL [22] have received much attention in the field of SR for their excellent performances. In [27], a modified SBL approach is proposed in the context of structured noise. The SBL approach is also extended to recover block sparse signals in [28]. The origin of the high performance of the M-SBL has been analyzed in [29]. One attractive property of the SBL and M-SBL is that they do not require any user parameter. Empirical results also show that the SBL and M-SBL are quite robust, but they converge very slowly. The Bayesian compressive sensing (BCS) [23], multitask 3
BCS (MT-BCS) [24], and variational Bayesian group-sparse estimation [25] are also Bayesian type SR algorithms, and they have their fast implementations, however they are proposed in the real domain and not applicable to the complex radar data. The complex multitask BCS (CMT-BCS) extends the MT-BCS to complex domain by exploiting the fact that the real and imaginary components of a complex value share the same sparsity pattern [26]. However, it is quite computational expensive. To improve the convergence of SBL, a fast converging SBL (FCSBL) approach is proposed for the SMV case. Specifically, the derivation of the FCSBL follows the same framework as the SBL but starts with a modified signal model. The key step in the derivation is the incorporation of a simple approximation term, which enables the proposed approach to achieve a favorable recovery result even with only a few times of iteration. The extension of the FCSBL to the MMV case is rather straightforward, and the resulting approach is named as M-FCSBL. Then the M-FCSBL approach is utilized to devise a novel STAP algorithm (M-FCSBL-STAP). In particular, the clutter profile and noise power are recovered by the M-FCSBL from a limited number of secondary data. Then the CCM is constructed based on the recovery results and the STAP weight vector is calculated by utilizing the CCM. Numerical experiments with both simulated and Mountain-Top data show that the proposed algorithm exhibit excellent clutter suppression performance in limited secondary data situation. The remainder of this paper is organized as follows. Section 2 provides a description of the signal model in airborne phased-array radar. In section 3, the FCSBL and M-FCSBL approaches are derived first, and then their convergences are analyzed, finally the M-FCSBL-STAP algorithm is proposed. Numerical experiments with both simulated and Mountain-Top data are carried out in Section 4 to assess the performance of the proposed algorithm. Concluding remarks are given in Section 5. Notation: Scalar quantities are denoted with lightface letters. Vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively.
and
represent the real filed and complex
field, respectively. The transpose, conjugate transpose, and inverse of a matrix are denoted by the superscripts
T
, , and H
1
respectively. The expectation operator is represented by E . The 4
trace of a square matrix is represented by Tr . stands for the absolute value of a scalar or the determinant of a square matrix. The 2 norm of a vector and the Frobenius norm of a matrix are denoted by
2
and F , respectively.
2. Signal model For tractability purpose, we consider an airborne pulse Doppler radar with a uniform linear array (ULA) consisting of N e elements. Each element receives N p pulses at a constant pulse repetition frequency (PRF) in a coherent processing interval (CPI). Let M N e N p . A general model for the space-time clutter plus noise snapshot x
M 1
is [30]
Nc
x n v f d ,n , f s ,n n
(1)
n 1
where N c is the number of independent clutter patches, n , f d ,n , f s ,n , and v f d ,n , f s ,n
M 1
are
the complex amplitude, normalized Doppler frequency, normalized spatial frequency, and space-time steering vector of the n th clutter patch, respectively, n is the thermal noise, modeled as zero-mean complex Gaussian vector with covariance matrix 2I , 2 is the noise power and I is the identity matrix. The clutter returns from all range gates are considered to be circularly symmetric zero-mean complex Gaussian process. Let us discretize the normalized space-time plane uniformly into K N s N d grid points, where N s (
N e ) is the number of normalized spatial frequency bins and N d (
N p ) is the number of
normalized Doppler frequency bins. Then each grid point is associated with a space-time steering vector v k ( k 1, ,K ). The STAP dictionary D
M K
is defined as the collection of these space-time
steering vectors, i.e.,
D v1 , , v K It is straightforward to recast the signal model (1) as the following form 5
(2)
x Dβ n
where β
K1
(3)
is the space-time profile with non-zero elements indicating the presence of clutter.
Actually, (3) is just the canonical signal model considered in the SR problem, which can be interpreted as estimating a sparse coefficient vector β , i.e., as few non-zero elements as possible, from noise contaminated measurement x . The approach of applying SR technique to STAP has been well studied in [11], and the authors argue that the approach does not recover the positions and amplitudes of clutter patches but recovers the clutter subspace using as small number of space-time steering vectors as possible. According to the linearly constrained minimum variance (LCMV) criterion, the optimal STAP weight vector is [30]
wOPT R 1s
(4)
where R is the CCM of the CUT, s is the target space-time steering vector, and sH R 1s
1
is a
normalization scalar to guarantee a constant response to s . Since R is unknown in practice, it has to be estimated from homogeneous secondary data. 3. Proposed M-FCSBL-STAP algorithm In this section, we first derive the FCSBL approach based on Bayesian inference and then extend it to the MMV case. After an investigation of their convergence issues, we propose the M-FCSBL-STAP algorithm and address some implementation issues. 3.1. Derivation of FCSBL Let Di
M K 1
denote the matrix with v i removed from D . Let βi
K 11
denote the
vector with i , i.e., the i th element of β , removed from β . The signal model (3) can be recast as
x i vi Di βi n
(5)
As mentioned above, n is modeled as zero-mean complex Gaussian vector with covariance matrix 6
2I . For ease of derivation, 2 is supposed to be known temporarily. To begin, we assign to each i a zero-mean Gaussian prior of the form i 2 1 p i i exp i π i
where i 0 ( i
(6)
) is the hyper-parameter controlling the prior variance of i . And the
hyper-parameters i ( i 1, ,K ) are assumed to be statistically independent. Let Mi denote the covariance matrix of
Di βi n , it can be formulated as Mi k v k v kH 2I
(7)
k i
The likelihood of x has the following form
p x i
1 H exp x i vi Mi1 x i vi π Mi M
(8)
The hyper-parameter i can be estimated from measured data by performing maximum likelihood (ML) optimization on the marginal likelihood p x i , which is defined as
p x i p x i p i i di
(9)
M Mi i vi viH
(10)
Let us define
Substituting (6) and (8) into (9) and integrating over i yields
p x i
1 exp x H M 1x π M M
(11)
Maximizing p x i is equivalent to minimizing ln p x i , this gives the cost function
CSMV ln M x H M1x Based on the Bayesian theorem, the posterior density of coefficient i for a fixed value of i is 7
(12)
p i x, i
p x i p i i p x i
(13)
Substituting (6), (8), and (11) into (13) yields i i 2 1 p i x, i exp π i i
(14)
where
i viH Mi1vi i1
1
(15)
and
i i viH Mi1x
(16)
The form of (14) indicates that the posterior density of i is Gaussian with mean i and variance i . Thus the estimate of i can be taken as the value of i . Rewriting (10) and applying the matrix inversion lemma results in
i M 1vi viH M 1 1 i viH M 1vi
Mi1 M 1
(17)
Substituting (17) into (15), after some algebra we have
i i 1 i viH M1vi
(18)
Substituting (17) and (18) into (16), we arrive at
i i viH M1x
(19)
i i viH M1x
(20)
Thus the estimate of i is
Or in a compact form, the estimate of β is given by
β ΓDH M1x where Γ is the diagonal matrix constructed by all the hyper-parameters i ( i 1, 8
(21)
,K ). Obviously,
the task remains in estimating the hyper-parameters. To find i , we resort to the expectation maximization (EM) algorithm to maximize p x i . First, treating the unknown coefficient i as a hidden variable, and calculating the expected value of
p i i with respect to the conditional distribution p i x, i
j
estimate of i , and superscript
j
, where
i j denotes the current
stands for the j th iteration. Then, find the parameter that
maximizes this quantity. Mathematically, the EM formulation can be described by the following steps. E STEP:
j E x , p i i p i i p i x, i di
(22)
i j 1 arg max E x, p i i 0
(23)
j
i
i
M STEP: j
i
i
i
Substituting (6) and (14) into (22), by integrating out i we have
E x , p i i j π i i
j
i
j j j j where i i 1 i viH M
i
1
1
j 2 i exp i i j
vi , i j i j viH M j
1
(24)
K
x , and M j k j v k v kH 2I . Thus k 1
the M STEP is equivalent to
i j 1 arg min ln E x , p i i 0 i
j
i
i
j
arg min ln i i i 0
i j i
2
(25)
i j
Taking the gradient of the target function on the right hand side of the second equality in (25) with respect to i , equating that to zero and solving the equation results in
i j 1 i j i j 2
9
(26)
Substituting the expressions of i j and i
j
into (26) yields
i j 1 i j viH M j
1
2
2
x i j i j viH M j
By the properties of standard Capon beamformer [31], when
i
1
vi
(27)
is accurate we have the
approximation
i viH M1vi
1
(28)
Applying (28) into the last two terms in (27) we arrive at the update formula of FCSBL, i.e.,
i j 1 i j viH M j
1
2
x
(29)
In above, 2 is supposed to be known, here we address the unknown case. To estimate 2 , the same procedure as in [21] can be utilized. For the sake of simplicity, we only present the final result. Let
μ j 1 j , , K j , 2 is updated by T
2
j 1
x Dμ
j
2 2
2
j
v M K
j
k
H k
j
1
k 1
M
Since the FCSBL is an iterative approach, the hyper-parameters i ( i 1,
i1 viH x
2
vk
(30)
,K ) can be initialized as
2
viH vi . Generally, 2 can be initialized with any positive scalar, but it would be more
preferable if a rough estimate, i.e., ˆ 2 , is provided. Concerning the termination, two stopping criteria are listed below: (a) Halt after a fixed number of iterations, i.e., J . (b) Halt when the relative change in γ 1 , γ j 1 γ j
2
γ j 1
2
, K
T
is smaller than a predetermined threshold , i.e.,
.
For clarity, the pseudo-code of the FCSBL is provided in Table 1.
10
Table 1 Pseudo-code of FCSBL Initialization:
2
1
ˆ 2
i1 viH x
2
2
viH vi , ( i 1,
,K )
M k v k v kH 2 I K
1
1
1
k 1
Repeat:
μ Γ DH M j
j
2
j 1
j
j 1
x
2 j K j j j 1 M x Dμ 2 k v kH M 2 k 1
i j 1 i j viH M j M
1
K
k k 1
j 1
1
1
vk
2
x , ( i 1, , K )
v k v kH 2
j 1
I
Until a stopping criterion is satisfied
3.2. Extension to the MMV case When multiple snapshots are available, say x l ( l 1,
, L ), let X x1 , ,x L , and define B and
N similarly, the SR signal model for the MMV case is X DΒ N
(31)
In this case, the estimate of B is a matrix form extension of (21), i.e.,
B ΓDH M1X
(32)
Let R ML XX H L , (29) can be straightforwardly extended to the MMV case, yields the update formula of M-FCSBL, i.e.,
i j 1 i j
2
viH M j
Let Γ j denote the estimate of Γ at the
1
R ML M j
1
vi
(33)
j th iteration, and let Ψ Γ DH M j
j
j
1
X.
Following the same procedure as in [22], 2 is updated by
2 j 1
X DΨ j
2 F
K j j L M k v kH M k 1
11
1
vk
(34)
The hyper-parameters i ( i 1,
L
,K ) can be initialized as i1 1 L viH x l
2
2
viH vi .
l 1
For clarity, the pseudo-code of the M-FCSBL is presented in Table 2. Also for completeness, the pseudo-code of the M-SBL in provided in the appendix. Table 2 Pseudo-code of M-FCSBL Initialization:
2
1
ˆ 2 L
i1 1 L viH x l
2
2
viH vi , ( i 1,
,K )
l 1
M k v k v kH 2 I K
1
1
1
k 1
R ML XX H L Repeat:
Ψ Γ DH M j
2
j
j 1
j 1
1
X
1 L X DΨ
i j 1 i j M
j
K
2
k k 1
viH M
j 1
j
1
j
2 F
R ML M
v k v kH 2
K j H j M k v k M k 1
j 1
j
1
1
vk
vi , ( i 1, , K )
I
Until a stopping criterion is satisfied
3.3 Convergence analysis It is difficult to make rigorous mathematical convergence proofs for the proposed approaches. Thus we make the following explanations instead. Equation (29) can be recast as
j
i j 1
i
v M H i
1 j
M
v iH M j
2
1
vi
2
v
H i
j
1
1
2
x vi
2
(35)
Let us define the following scalars
ci j
1
viH M j
12
1
(36)
vi
v
d i j
v
H i
H i
M M j
1
1
j
2
x
(37)
vi
Then (35) can be recast as i j j j di ci 2
j 1
i
Note that M j is a positive definite Hermitian matrix, ci estimator power spectrum [1], and d i
j
j
ci j
j
ci j
Before proceeding, a property of i
is the well-known minimum variance
can be recognized as the output power (corresponding to the
M j
v i component) by applying the adaptive processor easily observed from (38), the term i
j
(38)
1
vi
viH M j
1
vi
to x [1]. As can be
2
plays the role of a rescaling factor.
2
should be clarified.
Property: For the particular structure of M j , the relationship ci j i j will always hold true. And
we have i
j
ci j
2
1.
Proof: Let A k j v k v kH 2I , then M j can be reformulated as M j i j vi viH A . Note that k i
A is a positive definite matrix, and it is invertible. Using the matrix inversion lemma we get
M j
1
A 1
i j A 1vi viH A 1 1 i j viH A 1vi
(39)
Substituting (39) into (36), after some algebra yields
ci i j
As A 1 is a positive definite matrix, we have Note that d i
j
j
v
H i
1 v A 1vi
(40)
H i
j j A 1vi 0 , thus i ci 1
2
1 .□
is the estimated power of the v i component that exists in x , and its accuracy 13
depends on the weight vector of the adaptive processor. Since can be written as
M U
ui U1 2 vi
y U1 2x , then
and
j
1
1 2 H
of
H i
j
M j
H
1
H
j , we have x H M
1
j
1
is a positive definite matrix, it
U 1 2 , where U 1 2 is a square root matrix of
d i j
can be reformulated as j
2
u y u ui y y , then d i y y , i.e., d i x H i
M
2
H
M j
1
M j
di j uiH y uiH ui
1
. Let
2
. As
2
x . Let be the maximum eigenvalue
x x H x , note that 2 , then di j 2 x H x . Thus, d i 2
2
j
is upper bounded. The value of i j 1 can be discussed in the following cases. Case 1. i j 0 . In this case i j 1 will also be zero.
Case 2. i
j
ci j
2
0 . In this case i j 0 , and from (40) we have
v
H i
A 1vi
1
i j , which
means A contains many v k ( k 1, ,K and k i ) components that are correlated with v i (by saying this we mean the associated k j s have large values). As d i
j
is upper bounded, the value of
i j 1 will tend to be zero. The same result also suits to i j 2 as i j 1 ci j 1
Case 3. i
j
ci j
2
1 . In this case i j 0 , and from (40) we have
v
H i
2
0 will hold.
A 1vi 0 , which means 1
A hardly contains any v k ( k 1, ,K and k i ) components that are correlated with v i (by
saying this we mean the associated k j s have very small values or zero values). Substituting (39) into (37), after some algebra yields di viH A 1x viH A 1vi . Let x k v k n , the signal model for 2
j
k i
x
can be recast as
x i vi x . Substituting it into the formulation of
di j i viH A 1x viH A 1vi . Since 2
v
H i
A 1vi 0 , d i 1
14
j
d i j
yields
tends to be i . Thus i j 1 will 2
approximate to the accurate value. The same result also suits to i j 2 as i
Case 4. i
j
ci j
j 1
ci j 1
2
1 will hold.
2
neither approximates to be zero nor to be one. In this case, it is not necessary to
know what value will i j 1 have. But in the next iteration, the value of i j 2 will fall into one of the four cases. From above discussions we can see that if M j contains many v k ( k 1, ,K ) components that are correlated with each other, the case 2 will happen to some of them. Once the value of i j 1 fall into one of the first three cases, the value of i j 2 will fall into the same case. Finally the sequence of
i j will converge to a steady state. The particular structure of (38) promotes the FCSBL to converge much faster than the SBL. The analysis for the M-FCSBL follows in a similar way, and we omit it to save space. Additionally, the approximation term (28) indicates lim i j 0 , thus lim i j is the exact value of j
j
i . By taking the limit on both side of (29), we obtain
i
1
viH M
1
2
x
(41)
Applying (28) to (41) gives
viH M
1
vi viH M
1
xx H M
1
vi
(42)
Equation (42) indicates that M E xx H (in (42) x should be treated as the SR signal model). This relationship implies that the proposed approach can be utilized to get an effective CCM estimation, which is just desired by a STAP algorithm. 3.4. M-FCSBL-STAP filter Once the space-time profiles and the noise power have been estimated by the M-FCSBL, the CCM can be calculated by 15
K 2 L R 1 L i ,l vi viH 2I i 1 l 1
(43)
where i ,l denotes the element at the i th row and l th column of B , 1 is a real constant, and
2 is the estimated noise power. In (43), diagonal loading technique is adopted for avoiding distorted beam shapes and high sidelobes, and controls the loading level (satisfactory results are obtained in our experiments by considering 5 10 dB ). Another approach for calculating the CCM is K
R i vi viH 2I
(44)
i 1
where i denotes the estimate of i . Virtually, there is hardly any difference between (43) and (44) 2 L j j since lim i lim i,l L (this relationship can be easily drawn from (26) by the fact j j l 1
lim i j 0 ). Finally, the M-FCSBL-STAP weight vector is given by j
w M-FCSBL
R 1s H 1 s R s
(45)
The STAP dictionary D must be determined before applying the M-FCSBL-STAP algorithm. If N s and N d are too small, it would result in serious mismatch between the space-time steering vectors in D and those consistent with the clutter manifold. If N s and N d are too large, the computational
complexity would increase seriously. From our empirical study, it always leads to approving performance by considering N s 4 N e and N d 4 N p . Sometimes the noise power can be measured previously, for example by placing the radar in receive-only mode. If the noise power is known, the procedure for estimating it can be omitted from the M-FCSBL, and the computational load of the M-FCSBL can be further reduced. Even if the noise power is unknown, the M-FCSBL can estimate it effectively from measured data. Although the M-FCSBL-STAP algorithm is derived for ULA and constant PRF, it also suits for other array geometries and random slow-time samples. Another problem concerning the implementation of 16
the M-FCSBL-STAP is the selection of qualified secondary data. The recently proposed secondary data selection algorithm based on the clutter spectrum similarity may be a preferable candidate [2]. 4. Numerical experiments In this section, we first verify the advantages of the proposed FCSBL and M-FCSBL approaches, and then assess the performance of the M-FCSBL-STAP with both simulated and Mountain-Top data. We use the signal to interference plus noise ratio (SINR) loss [30] as a measure of performance, which is defined as the ratio of output SINR to the SNR can be achieved by a matched filter in an interference-free environment, i.e.,
LSINR
H 2 w s
2
M w H Rw
(46)
where w denotes the STAP weight vector. The simulated data are generated by considering the following scenario: a side-looking ULA with half wavelength inter-element spacing, N e 8 elements with uniform transmit pattern, wavelength is 0.23 m, PRF is 2434.8 Hz, N p 8 , 2 10 , platform height is 8 km, platform velocity is 140 m/s, N c 180 clutter patches uniformly distributed from azimuth π 2 to π 2 , the amplitude for each
clutter patch follows a complex Gaussian distribution, clutter to noise ratio (CNR) is 40 dB. 4.1 Advantages of FCSBL and M-FCSBL Fig. 1(a) provides the Capon spectrum which is calculated with the clairvoyant CCM. Figs. 1(b)–(e) illustrate the estimated clutter profiles for different approaches. Herein, 10 snapshots are used as secondary data for the M-SBL and M-FCSBL, the noise power is initialized as ˆ 2 1 , the SBL and M-SBL are terminated after 400 iterations, the FCSBL and M-FCSBL are terminated after 20 iterations. Figs. 1(b)–(e) are plotted according to the estimated hyper-parameters. Specifically, each of the estimated hyper-parameters is associated with a sampling point in the normalized space-time plane, and their values are depicted. 17
0.3 30 0.1 20 -0.1
10
-0.3
0 -0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
40 0.3 30 0.1 20 -0.1
0
-0.5 -0.5
-0.3
-0.1
0.1
40 0.3 30 0.1
20
-0.1
10
-0.3
0 -0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
0.5 40 0.3
30
0.1
20
-0.1
10
-0.3
0
-0.5 -0.5
-0.3
0.5
-0.1
0.3
0.5
(c) M-SBL 0.5 40
0.3
30 0.1 20 -0.1
10
-0.3
0
-0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(d) FCSBL
0.1
normalized Doppler frequency
(b) SBL
0.5
-0.5 -0.5
0.3
normalized Doppler frequency
(a) Clairvoyant Capon spectrum normalized spatial frequency
10
-0.3
normalized spatial frequency
-0.5 -0.5
0.5
normalized spatial frequency
40
normalized spatial frequency
normalized spatial frequency
0.5
0.5
(e) M-FCSBL
Fig. 1. Clairvoyant Capon spectrum and estimated clutter profiles
As shown in Figs. 1(b)–(e), the actual continuous clutter distribution converges into several disconnected clutter scatters for the SR approaches. It is common in the field of SR since the SR approaches do not try to recover the positions and amplitudes of clutter patches but recover the clutter subspace using as small number of space-time steering vectors as possible. We can also find that there are a few false and weak peaks induced by noise in Figs. 1(b) and 1(d); however the clutter profiles in Figs. 1(c) and 1(e) are quite favorable (the significant values are strictly located along the clutter ridge). This is because the M-SBL and M-FCSBL are more robust to noise compared with their SMV counterparts. In addition, we observe that the number of significant values in Fig. 1(e) equals to the clutter rank, which is 15 according to the Brennan’s rule [32]. This equality is in accordance with the results presented in [11]. Although Figs. 1(b)–(e) are plotted based on one arbitrary realization, the same conclusions can be drawn for any other realizations. For the SR-STAP algorithms, it is the estimated CCM that determines the clutter suppression performance. Thus it would be helpful to study the accuracy of the estimated CCM. For this purpose, we compare the Capon spectrums calculated with different CCM estimations. Particularly, the CCMs 18
0.5
normalized spatial frequency
normalized spatial frequency
are estimated based on (44) for different algorithms. The results are provided in Fig. 2.
40 0.3 30 0.1 20 -0.1
10
-0.3
0
-0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
0.5 40 0.3 30 0.1
20
-0.1
10
-0.3
0
-0.5 -0.5
-0.3
40 0.3 30 0.1
20
-0.1
10
-0.3
0 -0.3
-0.1
0.1
0.1
0.3
0.5
(b) M-SBL normalized spatial frequency
normalized spatial frequency
(a) SBL 0.5
-0.5 -0.5
-0.1
normalized Doppler frequency
0.3
normalized Doppler frequency
0.5
0.5 40
0.3
30 0.1 20 -0.1 10 -0.3
0
-0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(c) FCSBL
0.5
(d) M-FCSBL
Fig. 2. Capon spectrums calculated with different CCM estimations
By comparing Figs. 2(a)–(d) with Fig. 1(a), it is obvious that the Capon spectrums of SBL and FCSBL deviate a lot from the clairvoyant one. On the contrary, the Capon spectrums of M-SBL and M-FCSBL are very close to the clairvoyant one. By comparing Fig. 1(e) and Fig. 2(d), we can find that even though there are only 15 significant elements in the estimated clutter profile for the M-FCSBL algorithm, the Capon spectrum calculated with the estimated CCM approximates to the clairvoyant one very well. This result verifies that the proposed M-FCSBL approach can be utilized to get an effective CCM estimation. The cost function can be utilized to investigate the convergence performance of an approach. The cost functions for the SBL and M-SBL are the same as those for the FCSBL and M-FCSBL, respectively, with the latter being CMMV ln M Tr M1R ML . Since the SBL and M-SBL cannot accurately estimate the noise power, we use the known value for them. And 10 snapshots are used for the M-SBL and M-FCSBL. The results based on 100 Monte-Carlo trials are shown in Fig. 3. Some of the preceding values are very significant for the SBL and M-SBL, and these values are not shown for 19
the sake of easy comparison.
value of cost function (dB)
28
SBL M-SBL FCSBL M-FCSBL
27 26 25 24 0
100
200
300
number of iterations
400
500
Fig. 3. Value of cost function versus the number of iterations
As shown in Fig. 3, the four approaches converge to almost the same value (actually, the SBL converges to an even lower value compared with others, but this difference can be ignored since it is rather small). Observe that with as few as 20 iterations, the FCSBL and M-FCSBL converge to their steady-state values. However, the SBL and M-SBL converge quite slowly, and more than 400 iterations are needed for them to achieve their steady-state values. Thus the proposed FCSBL and M-FCSBL approaches converge much faster compared with the SBL and M-SBL approaches, respectively. The computational loads of the SBL, M-SBL, FCSBL, and M-FCSBL for a single iteration are summarized in Table 3. Herein the computational load is measured in term of the number of complex multiplications. Table 3 Comparison of computational load for a single iteration Approach SBL M-SBL FCSBL M-FCSBL
Computational load
M M M M
2 K 2 KL 2 K M KL 3K 2 6KM 3K M 1 2 KL 4 K L M 3K L
3
2 KM K 4 K 1 M K 3K 2
3
2 KM
2
3
4 KM
2
3
5KM
2
2
2
2
Fig. 4 gives a more direct illustration of the computational loads for the four approaches. It shows the computational load as a function of the system DOFs for a single iteration, furthermore it is supposed that K 16M and L 10 .
20
computational load
10
11
10
10
10
9
10
8
10
7
10
6
10
5
SBL M-SBL FCSBL M-FCSBL
0
100
200
M
300
400
500
Fig. 4. Computational load versus system DOFs for a single iteration
Clearly, for a single iteration, the FCSBL and M-FCSBL have lower computational loads than the SBL and M-SBL, respectively. In addition, as the FCSBL and M-FCSBL converge much faster than the SBL and M-SBL, respectively, the overall computational loads of the FCSBL and M-FCSBL are even lower than those of the SBL and M-SBL, respectively. 4.2 Performance of M-FCSBL-STAP A. Simulated data In Fig. 5, we assess the performance of the M-FCSBL-STAP by investigating the SINR loss in the main beam direction. The M-FCSBL-STAP is compared with the loaded sample matrix inversion (LSMI) algorithm, the SR-STAP algorithm based on M-SBL (M-SBL-STAP), and the optimum performance which is obtained using the clairvoyant CCM. The stopping criterion (b) in Section 3.1 is adopted for the two SR-STAP algorithms. To guarantee convergence, the threshold is set as 103 for the M-FCSBL-STAP while 105 for the M-SBL-STAP (a smaller implies more iterations). The loading factor for the LSMI is 10 dB to the noise power. The number of snapshots used for training is 10 for the two SR-STAP algorithms and 30 for the LSMI. All the presented results are averaged over 100 Monte Carlo trials.
21
SINR loss (dB)
0
-3
-6
-9 -0.5
LSMI M-SBL-STAP M-FCSBL-STAP Optimum
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
Fig. 5. SINR loss versus normalized Doppler frequency in main beam direction
As illustrated in Fig. 5, the LSMI approximates the -3dB loss since the number of training samples is twice the clutter rank. Note that the plain sample matrix inversion algorithm does not work for such few samples. On the other hand, the M-FCSBL-STAP can achieve comparable performance to the M-SBL-STAP, and their performances are very close to the optimum performance even though only 10 snapshots are used for training. Additionally, we found in the experiment, 480 iterations on average are used for the M-SBL-STAP, but this number is only 17 for the M-FCSBL-STAP. Fig. 6 presents the curves of average SINR loss versus the number of snapshots. Herein, the average SINR loss is defined as the mean of the SINR loss values in the entire normalized Doppler frequency range. The thresholds for the two SR-STAP algorithms are the same as above.
average SINR loss (dB)
0
-3
-6
-9 0
LSMI M-SBL-STAP M-FCSBL-STAP Optimum
10
20
30
40
number of snapshots
50
60
Fig. 6. Average SINR loss versus the number of snapshots
As depicted in Fig. 6, the LSMI converges quite slowly. On the contrary, the M-SBL-STAP and M-FCSBL-STAP achieve both rapid convergence rates and favorable steady-state performances. The results in Figs. 5 and 6 indicate that the proposed M-FCSBL-STAP can achieve comparable 22
performance to the M-SBL-STAP. However, as the M-FCSBL requires far less iterations than the M-SBL to obtain a favorable estimated clutter profile, the M-FCSBL-STAP is more suitable for application. B. Measured data Then we apply the M-FCSBL-STAP to the public available Mountain-Top data set [33]. For this data set, the receive antenna array consists of 14 elements, echoes from 16 coherent pulses are organized into a CPI data. A 500 kHz linear frequency modulated pulse was used for transmitting, the pulse width was 100 microseconds, and the PRF was 625 Hz. The data file t38pre01v1 CPI6 is used in our experiment. In this data file, there are 403 range snapshots, the clutter locates around 245° relative to true North and the target locates at 275° relative to true North with the same normalized Doppler frequency 0.25. The
normalized spatial frequency
clutter power spectrum estimated by all 403 range cells is depicted in Fig. 7. 0.5
0
0.3
-5
0.1
-10
-0.1 -15 -0.3 -0.5 -0.5
-20 -0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
Fig. 7. Estimated clutter spectrum using all 403 range snapshots
We adopt a symmetrical sliding window to detect the target. For each CUT, the eight snapshots around the CUT are treated as guard cells, and 10 snapshots out of the 20 snapshots located in the window are selected as secondary data based on the selection strategy proposed in [2]. The parameters for the M-SBL-STAP and M-FCSBL-STAP are the same as above. Fig. 8 displays the adapted patterns for the LSMI, M-SBL-STAP, and M-FCSBL-STAP at the target range cell. As expected, deep clutter notches are formed at the clutter ridge region for the M-SBL-STAP and M-FCSBL-STAP. These pattern notches will suppress the clutter at the output to well below thermal noise. However, the clutter notch 23
0
0.3
-20
0.1 -40 -0.1 -60
-0.3 -0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
-80
0.5
0
0.3
normalized spatial frequency
0.5
normalized spatial frequency
normalized spatial frequency
formed by the LSMI is not so deep.
-20
0.1 -40 -0.1 -60
-0.3 -0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(a) LSMI
0.5
-80
0.5
0
0.3
-20
0.1 -40 -0.1 -60
-0.3 -0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(b) M-SBL-STAP
0.5
-80
(c) M-FCSBL-STAP
Fig. 8. Adapted patterns at the target range cell
Fig. 9 shows the STAP output power for different algorithms in the range of 148–160 km. We can observe that peaks are formed at the target (at 154 km) range cell for all of the algorithms. However, there are some clutter residuals for the LSMI (at the few range cells around 151.8 km, 152.8 km, and 154.9 km). The target is clearly visible for both SR-STAP algorithms even though only 10 snapshots are used as secondary data. Indeed, the two SR-STAP algorithms show nearly the same output powers. However, the M-FCSBL-STAP is more favorable for practical application as it needs far fewer
output power (dB)
iterations and lower computation cost at each iteration. LSMI M-SBL-STAP M-FCSBL-STAP
80 70 60 50 148
150
152
154
range (km)
156
158
160
Fig. 9. Range plots using 10 snapshots as secondary data
5. Conclusions We derived the FCSBL approach to improve the convergence performance of the SBL approach in this paper. We also extended the proposed FCSBL approach to the MMV case. And the proposed M-FCSBL approach is utilized to devise a SR-STAP algorithm for the benefit of drastically reduced training requirement. Numerical results confirm that the M-FCSBL-STAP can achieve comparable 24
clutter suppression performance to the M-SBL-STAP. As the M-FCSBL approach enjoys several favorable properties such as user parameter free and fast converging, the resulting M-FCSBL-STAP algorithm is more preferable for application. Acknowledgment This work was supported in part by National Natural Science Foundation of China under Grant 61501506. The authors would like to thank the anonymous reviewers for their valuable comments. Appendix For completeness, the pseudo-code of the M-SBL is provided in Table 4. Table 4 Pseudo-code of M-SBL Initialization:
2
1
ˆ 2 L
i1 1 L viH x l
2
2
viH vi , ( i 1,
,K )
l 1
M k v k v kH 2 I K
1
1
1
k 1
Repeat: j
j
Ψ Γ DH j
2
M
Γ Γ DH M j
j 1
j
2
K
k k 1
j 1
1
DΓ
j
2
j
X
1 L X DΨ
i j 1 1 L Ψi,j j 1
1
j
2
M
j
j i ,i
F
K M K k 1
, ( i 1,
v k v kH 2
j 1
j k ,k
k j
,K )
I
Until a stopping criterion is satisfied
In Table 4,
j k ,k
denotes the k th diagonal elements of
j
, Ψi ,j stands for the i th row of Ψ j .
References [1]
R. Klemm, Principles of space-time adaptive processing, 3d ed., IET, London, 2006.
[2]
Y. Wu, T. Wang, J. Wu, J. Duan, Robust training samples selection algorithm based on spectral similarity for space–time adaptive processing in heterogeneous interference environments, IET Radar, Sonar and Navigation. 25
9 (7) (2015) 778–782. [3]
S. Kang, J. Ryu, J. Lee, J. Jeong, Analysis of space–time adaptive processing performance using K-means clustering algorithm for normalisation method in non-homogeneity detector process, IET Signal Processing. 5 (2) (2011) 113–120.
[4]
W. L. Melvin, A STAP overview, IEEE Aerospace and Electronic Systems Magazine. 19 (1) (2004) 19–35.
[5]
J. R. Roman, M. Rangaswamy, D. W. Davis, Q. Zhang, B. Himed, J. H. Michels, Parametric adaptive matched filter for airborne radar applications, IEEE Transactions on Aerospace and Electronic Systems. 36 (2) (2000) 677–692.
[6]
S. Burintramart, T. K. Sarkar, Y. Zhang, M. C. Wicks, Performance comparison between statistical-based and direct data domain STAPs, Digital Signal Processing. 17 (4) (2007) 737–755.
[7]
J. R. Guerci, E. J. Baranoski, Knowledge-aided adaptive radar at DARPA: an overview, IEEE Signal Processing Magazine. 23 (1) (2006) 41–50.
[8]
X. M. Zhu, J. Li, P. Stoica, Knowledge-aided space-time adaptive processing, IEEE Transactions on Aerospace and Electronic Systems. 47 (2) (2011) 1325–1336.
[9]
J. A. Tropp, S. J. Wright, Computational methods for sparse solution of linear inverse problems, Proceedings of the IEEE. 98 (6) (2010) 948–958.
[10] S. Foucart, H. Rauhut, A mathematical introduction to compressive sensing, Springer, New York, 2013. [11] Z. Yang, X. Li, H. Wang, W. Jiang, On clutter sparsity analysis in space-time adaptive processing airborne radar, IEEE Geoscience and Remote Sensing Letters. 10 (2013) 1214–1218. [12] K. Sun, H. Meng, Y. Wang, X. Wang, Direct data domain STAP using sparse representation of clutter spectrum, Signal Processing. 91 (9) (2011) 2222–2236. [13] K. Sun, H. D. Meng, F. D. Lapierre, X. Q. Wang, Registration based compensation using sparse representation in conformal array STAP, Signal Processing. 91 (10) (2011) 2268–2276. [14] Z. Ma, Y. Liu, H. Meng, X. Wang, Jointly sparse recovery of multiple snapshots in STAP, in: Proceedings of IEEE Radar Conference, April 29 2013–May 3 2013, pp. 1–4. [15] S. Sen, OFDM radar space-time adaptive processing by exploiting spatio-temporal sparsity, IEEE Transactions on Signal Processing. 61 (1) (2013) 118–130. [16] Z. Yang, X. Li, H. Wang, W. Jiang, Adaptive clutter suppression based on iterative adaptive approach for 26
airborne radar, Signal Processing. 93 (12) (2013) 3567–3577. [17] Z. Yang, X. Li, H. Wang, L. Nie, Sparsity-based space-time adaptive processing using complex-valued Homotopy technique for airborne radar, IET Signal Processing. 8 (5) (2014) 552–564. [18] Z. C. Yang, R. Fa, Y. L. Qin, X. Li, H. Q. Wang, Direct data domain sparsity-based STAP utilizing subaperture smoothing techniques, International Journal of Antennas and Propagation. 2015 (2015) 1–10. [19] Q. Wu, Y. D. Zhang, M. G. Amin, B. Himed, Space–time adaptive processing and motion parameter estimation in multistatic passive radar using sparse Bayesian learning, IEEE Transactions on Geoscience and Remote Sensing. 54 (2) (2016) 944–957. [20] M. E. Tipping, Sparse Bayesian learning and the relevance vector machine, Journal of Machine Learning Research. 1 (2001) 211–244. [21] D. P. Wipf, B. D. Rao, Sparse Bayesian learning for basis selection, IEEE Transactions on Signal Processing. 52 (8) (2004) 2153–2164. [22] D. P. Wipf, B. D. Rao, An empirical Bayesian strategy for solving the simultaneous sparse approximation problem, IEEE Transactions on Signal Processing. 55 (7) (2007) 3704–3716. [23] S. Ji, Y. Xue, L. Carin, Bayesian compressive sensing, IEEE Transactions on Signal Processing. 56 (6) (2008) 2346–2356. [24] S. Ji, D. Dunson, L. Carin, Multitask compressive sensing, IEEE Transactions on Signal Processing. 57 (1) (2009) 92–106. [25] S. D. Babacan, S. Nakajima, N. N. Do, Bayesian group-sparse modeling and variational inference, IEEE Transactions on Signal Processing. 62 (11) (2014) 2906–2921. [26] Q. Wu, Y. D. Zhang, M. G. Amin, B. Himed, Complex multitask Bayesian compressive sensing, in: Proceedings of 2014 IEEE International Conference on Acoustic, Speech and Signal Processing, 4–9 May 2014, pp. 3375– 3379. [27] M. Hurtado, C. H. Muravchik, A. Nehorai, Enhanced sparse Bayesian learning via statistical thresholding for signals in structured noise, IEEE Transactions on Signal Processing. 61 (21) (2013) 5430–5443. [28] Z. L. Zhang, B. D. Rao, Extension of SBL algorithms for the recovery of block sparse signals with intra-block correlation, IEEE Transactions on Signal Processing. 61 (8) (2013) 2009–2015. [29] J. C. Ye, J. M. Kim, Y. Bresler, Improving M-SBL for joint sparse recovery using a subspace penalty, IEEE 27
Transactions on Signal Processing. 63 (24) (2015) 6595–6605. [30] J. Ward, Space-time adaptive processing for airborne radar, Technical report 1015, MIT Lincoln Lab, Lexington, MA, USA, 1994. [31] J. Li, P. Stoica, Z. S. Wang, On robust Capon beamforming and diagonal loading, IEEE Transactions on Signal Processing. 51 (7) (2003) 1702–1715. [32] L. E. Brennan, F. M. Staudaher, Subclutter visibility demonstration, RL-TR-92-21, Adaptive Sensors, Inc., Santa Monica, CA, 1992. [33] G. W. Titi, D. F. Marshall, The ARPA/NAVY mountaintop program: adaptive signal processing for airborne early warning radar, in: Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing, May 1996, pp. 1165–1168.
Highlights
We study the problem of clutter suppression in STAP with finite training samples. (c)
Fast converging sparse Bayesian learning approaches are derived. (d)
A novel STAP algorithm named as M-FCSBL-STAP is proposed. (e)
The M-FCSBL-STAP has superior performance in low training support situation.
28
PII: DOI: Reference:
S0165-1684(16)30134-7 http://dx.doi.org/10.1016/j.sigpro.2016.06.023 SIGPRO6183
To appear in: Signal Processing Received date: 29 January 2016 Revised date: 21 June 2016 Accepted date: 23 June 2016 Cite this article as: Zetao Wang, Wenchong Xie, Keqing Duan and Yongliang Wang, Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar, Signal Processing, http://dx.doi.org/10.1016/j.sigpro.2016.06.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar Zetao Wanga*, Wenchong Xieb, Keqing Duanb, Yongliang Wangb a
College of Electronic Science and Engineering, NUDT, Changsha 410073, China
b
Wuhan Radar Academy, Wuhan 430019, China
[email protected] [email protected] [email protected] [email protected] *
Corresponding author. Tel.: +86 2785965760.
Abstract Adapting the space-time adaptive processing (STAP) filter with finite number of secondary data is of particular interest for airborne phased-array radar clutter suppression. Sparse representation (SR) technique has been introduced into the STAP framework for the benefit of drastically reduced training requirement. However, most SR algorithms need the fine tuning of one or more user parameters, which affect the final results significantly. Sparse Bayesian learning (SBL) and multiple sparse Bayesian learning (M-SBL) are robust and user parameter free approaches, but they converge quite slowly. To remedy this limitation, a fast converging SBL (FCSBL) approach is proposed based on Bayesian inference along with a simple approximation term, then, it is extended to the multiple measurement vector case, and the resulting approach is termed as M-FCSBL. To improve the performance of STAP in finite secondary data situation, the M-FCSBL is utilized to estimate the clutter plus noise covariance matrix (CCM) from a limited number of secondary data, and then the resulting CCM is adopted to devise the STAP filter and suppress the clutter. Numerical experiments with both simulated and Mountain-Top data are carried out. It is shown that the proposed algorithm has superior clutter suppression performance in finite secondary data situation. Keywords: space-time adaptive processing, clutter suppression, sparse representation, Bayesian 1
inference. 1. Introduction Space-time adaptive processing (STAP) holds the inherent capacity of detecting slowly moving targets in a strong clutter background for airborne phased-array radar. The optimum STAP weight vector is formulated as the product of the inverse of clutter plus noise covariance matrix (CCM) and the target space-time steering vector [1]. When testing one range cell for the existence of a target, the CCM is usually estimated from secondary/training data that are selected from the range cells adjacent to the cell under test (CUT) [2]. In such a regime, the secondary data are assumed to be target-free and contain interference with the same statistical characteristics as the CUT. However, as a result of both geometric configuration and terrain variations, the qualified secondary data may be rather limited [3]. And the performance of STAP degrades drastically due to the inaccuracy of the estimated CCM. To increase the performance of STAP in finite secondary data situation, many different types of algorithms have been developed. Earlier achievements include reduced-dimension (RD) and reduced-rank (RR) algorithms which reduce the training requirement to twice of the reduced-dimension and twice of the clutter rank, respectively [4]. Recent developments include the parametric adaptive matched filter (PAMF) [5], direct data domain (D3) algorithms [6], and knowledge-aided (KA) algorithms [7, 8]. Although the training requirements are further reduced for these algorithms, they suffer from several drawbacks. The order of the autoregressive model for the PAMF is difficult to be determined, the system degrees of freedom (DOFs) are significantly reduced for the D3 algorithms, and the performances of the KA algorithms depend heavily on the accuracy of a priori knowledge. In most recent years, sparse representation/recovery (SR) technique has received continuous attention fueled by advanced theoretical developments and practical applications [9, 10]. It provides a general signal model that represents or approximates a signal under an overcomplete dictionary with a sparse coefficient vector. The dictionary in the context of STAP is a collection of space-time steering vectors each associated with a sampling point in the normalized space-time plane. And “sparse” means that 2
most of the elements in the coefficient vector are equal to zero. The approach of applying SR technique to STAP has been well studied in [11] and it was pointed out that the approach tries to recover the clutter subspace using as few space-time steering vectors as possible. STAP algorithms based on SR technique (SR-STAP) can be summarized as two-step procedures: first, recovering the clutter space-time profile by a certain SR algorithm; second, constructing the CCM with the result obtained in the first step and calculating the STAP weight vector [12-19]. For example, the algorithm proposed in [19] first estimates the clutter profile in the CUT by considering only the clutter entries within the learned clutter support (the clutter support is estimated from secondary data by SR technique), and then devise the STAP filter based on the estimated clutter profile. The most significant benefit brought by SR-STAP is that the training requirement can be drastically reduced without incorporating any a priori knowledge of clutter. As can be inferred from the two-step scheme of SR-STAP, the final clutter suppression performance depends heavily on the recovery result which pertains to the SR algorithm adopted. Unfortunately, most SR algorithms need the fine tuning of one or more user parameters which affect the recovery results significantly. For example, the Homotopy-STAP algorithm proposed in [17] requires the fine tuning of the regularization parameter, and the final clutter suppression performance is affected by this parameter. Bayesian type algorithms [20-26] have developed the reputation of favorable performance when the dictionary is of a high coherence. Among them, the sparse Bayesian learning (SBL) [20, 21], which was proposed for the single measurement vector (SMV) case, and its counterpart for the multiple measurement vector (MMV) case, i.e., M-SBL [22] have received much attention in the field of SR for their excellent performances. In [27], a modified SBL approach is proposed in the context of structured noise. The SBL approach is also extended to recover block sparse signals in [28]. The origin of the high performance of the M-SBL has been analyzed in [29]. One attractive property of the SBL and M-SBL is that they do not require any user parameter. Empirical results also show that the SBL and M-SBL are quite robust, but they converge very slowly. The Bayesian compressive sensing (BCS) [23], multitask 3
BCS (MT-BCS) [24], and variational Bayesian group-sparse estimation [25] are also Bayesian type SR algorithms, and they have their fast implementations, however they are proposed in the real domain and not applicable to the complex radar data. The complex multitask BCS (CMT-BCS) extends the MT-BCS to complex domain by exploiting the fact that the real and imaginary components of a complex value share the same sparsity pattern [26]. However, it is quite computational expensive. To improve the convergence of SBL, a fast converging SBL (FCSBL) approach is proposed for the SMV case. Specifically, the derivation of the FCSBL follows the same framework as the SBL but starts with a modified signal model. The key step in the derivation is the incorporation of a simple approximation term, which enables the proposed approach to achieve a favorable recovery result even with only a few times of iteration. The extension of the FCSBL to the MMV case is rather straightforward, and the resulting approach is named as M-FCSBL. Then the M-FCSBL approach is utilized to devise a novel STAP algorithm (M-FCSBL-STAP). In particular, the clutter profile and noise power are recovered by the M-FCSBL from a limited number of secondary data. Then the CCM is constructed based on the recovery results and the STAP weight vector is calculated by utilizing the CCM. Numerical experiments with both simulated and Mountain-Top data show that the proposed algorithm exhibit excellent clutter suppression performance in limited secondary data situation. The remainder of this paper is organized as follows. Section 2 provides a description of the signal model in airborne phased-array radar. In section 3, the FCSBL and M-FCSBL approaches are derived first, and then their convergences are analyzed, finally the M-FCSBL-STAP algorithm is proposed. Numerical experiments with both simulated and Mountain-Top data are carried out in Section 4 to assess the performance of the proposed algorithm. Concluding remarks are given in Section 5. Notation: Scalar quantities are denoted with lightface letters. Vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively.
and
represent the real filed and complex
field, respectively. The transpose, conjugate transpose, and inverse of a matrix are denoted by the superscripts
T
, , and H
1
respectively. The expectation operator is represented by E . The 4
trace of a square matrix is represented by Tr . stands for the absolute value of a scalar or the determinant of a square matrix. The 2 norm of a vector and the Frobenius norm of a matrix are denoted by
2
and F , respectively.
2. Signal model For tractability purpose, we consider an airborne pulse Doppler radar with a uniform linear array (ULA) consisting of N e elements. Each element receives N p pulses at a constant pulse repetition frequency (PRF) in a coherent processing interval (CPI). Let M N e N p . A general model for the space-time clutter plus noise snapshot x
M 1
is [30]
Nc
x n v f d ,n , f s ,n n
(1)
n 1
where N c is the number of independent clutter patches, n , f d ,n , f s ,n , and v f d ,n , f s ,n
M 1
are
the complex amplitude, normalized Doppler frequency, normalized spatial frequency, and space-time steering vector of the n th clutter patch, respectively, n is the thermal noise, modeled as zero-mean complex Gaussian vector with covariance matrix 2I , 2 is the noise power and I is the identity matrix. The clutter returns from all range gates are considered to be circularly symmetric zero-mean complex Gaussian process. Let us discretize the normalized space-time plane uniformly into K N s N d grid points, where N s (
N e ) is the number of normalized spatial frequency bins and N d (
N p ) is the number of
normalized Doppler frequency bins. Then each grid point is associated with a space-time steering vector v k ( k 1, ,K ). The STAP dictionary D
M K
is defined as the collection of these space-time
steering vectors, i.e.,
D v1 , , v K It is straightforward to recast the signal model (1) as the following form 5
(2)
x Dβ n
where β
K1
(3)
is the space-time profile with non-zero elements indicating the presence of clutter.
Actually, (3) is just the canonical signal model considered in the SR problem, which can be interpreted as estimating a sparse coefficient vector β , i.e., as few non-zero elements as possible, from noise contaminated measurement x . The approach of applying SR technique to STAP has been well studied in [11], and the authors argue that the approach does not recover the positions and amplitudes of clutter patches but recovers the clutter subspace using as small number of space-time steering vectors as possible. According to the linearly constrained minimum variance (LCMV) criterion, the optimal STAP weight vector is [30]
wOPT R 1s
(4)
where R is the CCM of the CUT, s is the target space-time steering vector, and sH R 1s
1
is a
normalization scalar to guarantee a constant response to s . Since R is unknown in practice, it has to be estimated from homogeneous secondary data. 3. Proposed M-FCSBL-STAP algorithm In this section, we first derive the FCSBL approach based on Bayesian inference and then extend it to the MMV case. After an investigation of their convergence issues, we propose the M-FCSBL-STAP algorithm and address some implementation issues. 3.1. Derivation of FCSBL Let Di
M K 1
denote the matrix with v i removed from D . Let βi
K 11
denote the
vector with i , i.e., the i th element of β , removed from β . The signal model (3) can be recast as
x i vi Di βi n
(5)
As mentioned above, n is modeled as zero-mean complex Gaussian vector with covariance matrix 6
2I . For ease of derivation, 2 is supposed to be known temporarily. To begin, we assign to each i a zero-mean Gaussian prior of the form i 2 1 p i i exp i π i
where i 0 ( i
(6)
) is the hyper-parameter controlling the prior variance of i . And the
hyper-parameters i ( i 1, ,K ) are assumed to be statistically independent. Let Mi denote the covariance matrix of
Di βi n , it can be formulated as Mi k v k v kH 2I
(7)
k i
The likelihood of x has the following form
p x i
1 H exp x i vi Mi1 x i vi π Mi M
(8)
The hyper-parameter i can be estimated from measured data by performing maximum likelihood (ML) optimization on the marginal likelihood p x i , which is defined as
p x i p x i p i i di
(9)
M Mi i vi viH
(10)
Let us define
Substituting (6) and (8) into (9) and integrating over i yields
p x i
1 exp x H M 1x π M M
(11)
Maximizing p x i is equivalent to minimizing ln p x i , this gives the cost function
CSMV ln M x H M1x Based on the Bayesian theorem, the posterior density of coefficient i for a fixed value of i is 7
(12)
p i x, i
p x i p i i p x i
(13)
Substituting (6), (8), and (11) into (13) yields i i 2 1 p i x, i exp π i i
(14)
where
i viH Mi1vi i1
1
(15)
and
i i viH Mi1x
(16)
The form of (14) indicates that the posterior density of i is Gaussian with mean i and variance i . Thus the estimate of i can be taken as the value of i . Rewriting (10) and applying the matrix inversion lemma results in
i M 1vi viH M 1 1 i viH M 1vi
Mi1 M 1
(17)
Substituting (17) into (15), after some algebra we have
i i 1 i viH M1vi
(18)
Substituting (17) and (18) into (16), we arrive at
i i viH M1x
(19)
i i viH M1x
(20)
Thus the estimate of i is
Or in a compact form, the estimate of β is given by
β ΓDH M1x where Γ is the diagonal matrix constructed by all the hyper-parameters i ( i 1, 8
(21)
,K ). Obviously,
the task remains in estimating the hyper-parameters. To find i , we resort to the expectation maximization (EM) algorithm to maximize p x i . First, treating the unknown coefficient i as a hidden variable, and calculating the expected value of
p i i with respect to the conditional distribution p i x, i
j
estimate of i , and superscript
j
, where
i j denotes the current
stands for the j th iteration. Then, find the parameter that
maximizes this quantity. Mathematically, the EM formulation can be described by the following steps. E STEP:
j E x , p i i p i i p i x, i di
(22)
i j 1 arg max E x, p i i 0
(23)
j
i
i
M STEP: j
i
i
i
Substituting (6) and (14) into (22), by integrating out i we have
E x , p i i j π i i
j
i
j j j j where i i 1 i viH M
i
1
1
j 2 i exp i i j
vi , i j i j viH M j
1
(24)
K
x , and M j k j v k v kH 2I . Thus k 1
the M STEP is equivalent to
i j 1 arg min ln E x , p i i 0 i
j
i
i
j
arg min ln i i i 0
i j i
2
(25)
i j
Taking the gradient of the target function on the right hand side of the second equality in (25) with respect to i , equating that to zero and solving the equation results in
i j 1 i j i j 2
9
(26)
Substituting the expressions of i j and i
j
into (26) yields
i j 1 i j viH M j
1
2
2
x i j i j viH M j
By the properties of standard Capon beamformer [31], when
i
1
vi
(27)
is accurate we have the
approximation
i viH M1vi
1
(28)
Applying (28) into the last two terms in (27) we arrive at the update formula of FCSBL, i.e.,
i j 1 i j viH M j
1
2
x
(29)
In above, 2 is supposed to be known, here we address the unknown case. To estimate 2 , the same procedure as in [21] can be utilized. For the sake of simplicity, we only present the final result. Let
μ j 1 j , , K j , 2 is updated by T
2
j 1
x Dμ
j
2 2
2
j
v M K
j
k
H k
j
1
k 1
M
Since the FCSBL is an iterative approach, the hyper-parameters i ( i 1,
i1 viH x
2
vk
(30)
,K ) can be initialized as
2
viH vi . Generally, 2 can be initialized with any positive scalar, but it would be more
preferable if a rough estimate, i.e., ˆ 2 , is provided. Concerning the termination, two stopping criteria are listed below: (a) Halt after a fixed number of iterations, i.e., J . (b) Halt when the relative change in γ 1 , γ j 1 γ j
2
γ j 1
2
, K
T
is smaller than a predetermined threshold , i.e.,
.
For clarity, the pseudo-code of the FCSBL is provided in Table 1.
10
Table 1 Pseudo-code of FCSBL Initialization:
2
1
ˆ 2
i1 viH x
2
2
viH vi , ( i 1,
,K )
M k v k v kH 2 I K
1
1
1
k 1
Repeat:
μ Γ DH M j
j
2
j 1
j
j 1
x
2 j K j j j 1 M x Dμ 2 k v kH M 2 k 1
i j 1 i j viH M j M
1
K
k k 1
j 1
1
1
vk
2
x , ( i 1, , K )
v k v kH 2
j 1
I
Until a stopping criterion is satisfied
3.2. Extension to the MMV case When multiple snapshots are available, say x l ( l 1,
, L ), let X x1 , ,x L , and define B and
N similarly, the SR signal model for the MMV case is X DΒ N
(31)
In this case, the estimate of B is a matrix form extension of (21), i.e.,
B ΓDH M1X
(32)
Let R ML XX H L , (29) can be straightforwardly extended to the MMV case, yields the update formula of M-FCSBL, i.e.,
i j 1 i j
2
viH M j
Let Γ j denote the estimate of Γ at the
1
R ML M j
1
vi
(33)
j th iteration, and let Ψ Γ DH M j
j
j
1
X.
Following the same procedure as in [22], 2 is updated by
2 j 1
X DΨ j
2 F
K j j L M k v kH M k 1
11
1
vk
(34)
The hyper-parameters i ( i 1,
L
,K ) can be initialized as i1 1 L viH x l
2
2
viH vi .
l 1
For clarity, the pseudo-code of the M-FCSBL is presented in Table 2. Also for completeness, the pseudo-code of the M-SBL in provided in the appendix. Table 2 Pseudo-code of M-FCSBL Initialization:
2
1
ˆ 2 L
i1 1 L viH x l
2
2
viH vi , ( i 1,
,K )
l 1
M k v k v kH 2 I K
1
1
1
k 1
R ML XX H L Repeat:
Ψ Γ DH M j
2
j
j 1
j 1
1
X
1 L X DΨ
i j 1 i j M
j
K
2
k k 1
viH M
j 1
j
1
j
2 F
R ML M
v k v kH 2
K j H j M k v k M k 1
j 1
j
1
1
vk
vi , ( i 1, , K )
I
Until a stopping criterion is satisfied
3.3 Convergence analysis It is difficult to make rigorous mathematical convergence proofs for the proposed approaches. Thus we make the following explanations instead. Equation (29) can be recast as
j
i j 1
i
v M H i
1 j
M
v iH M j
2
1
vi
2
v
H i
j
1
1
2
x vi
2
(35)
Let us define the following scalars
ci j
1
viH M j
12
1
(36)
vi
v
d i j
v
H i
H i
M M j
1
1
j
2
x
(37)
vi
Then (35) can be recast as i j j j di ci 2
j 1
i
Note that M j is a positive definite Hermitian matrix, ci estimator power spectrum [1], and d i
j
j
ci j
j
ci j
Before proceeding, a property of i
is the well-known minimum variance
can be recognized as the output power (corresponding to the
M j
v i component) by applying the adaptive processor easily observed from (38), the term i
j
(38)
1
vi
viH M j
1
vi
to x [1]. As can be
2
plays the role of a rescaling factor.
2
should be clarified.
Property: For the particular structure of M j , the relationship ci j i j will always hold true. And
we have i
j
ci j
2
1.
Proof: Let A k j v k v kH 2I , then M j can be reformulated as M j i j vi viH A . Note that k i
A is a positive definite matrix, and it is invertible. Using the matrix inversion lemma we get
M j
1
A 1
i j A 1vi viH A 1 1 i j viH A 1vi
(39)
Substituting (39) into (36), after some algebra yields
ci i j
As A 1 is a positive definite matrix, we have Note that d i
j
j
v
H i
1 v A 1vi
(40)
H i
j j A 1vi 0 , thus i ci 1
2
1 .□
is the estimated power of the v i component that exists in x , and its accuracy 13
depends on the weight vector of the adaptive processor. Since can be written as
M U
ui U1 2 vi
y U1 2x , then
and
j
1
1 2 H
of
H i
j
M j
H
1
H
j , we have x H M
1
j
1
is a positive definite matrix, it
U 1 2 , where U 1 2 is a square root matrix of
d i j
can be reformulated as j
2
u y u ui y y , then d i y y , i.e., d i x H i
M
2
H
M j
1
M j
di j uiH y uiH ui
1
. Let
2
. As
2
x . Let be the maximum eigenvalue
x x H x , note that 2 , then di j 2 x H x . Thus, d i 2
2
j
is upper bounded. The value of i j 1 can be discussed in the following cases. Case 1. i j 0 . In this case i j 1 will also be zero.
Case 2. i
j
ci j
2
0 . In this case i j 0 , and from (40) we have
v
H i
A 1vi
1
i j , which
means A contains many v k ( k 1, ,K and k i ) components that are correlated with v i (by saying this we mean the associated k j s have large values). As d i
j
is upper bounded, the value of
i j 1 will tend to be zero. The same result also suits to i j 2 as i j 1 ci j 1
Case 3. i
j
ci j
2
1 . In this case i j 0 , and from (40) we have
v
H i
2
0 will hold.
A 1vi 0 , which means 1
A hardly contains any v k ( k 1, ,K and k i ) components that are correlated with v i (by
saying this we mean the associated k j s have very small values or zero values). Substituting (39) into (37), after some algebra yields di viH A 1x viH A 1vi . Let x k v k n , the signal model for 2
j
k i
x
can be recast as
x i vi x . Substituting it into the formulation of
di j i viH A 1x viH A 1vi . Since 2
v
H i
A 1vi 0 , d i 1
14
j
d i j
yields
tends to be i . Thus i j 1 will 2
approximate to the accurate value. The same result also suits to i j 2 as i
Case 4. i
j
ci j
j 1
ci j 1
2
1 will hold.
2
neither approximates to be zero nor to be one. In this case, it is not necessary to
know what value will i j 1 have. But in the next iteration, the value of i j 2 will fall into one of the four cases. From above discussions we can see that if M j contains many v k ( k 1, ,K ) components that are correlated with each other, the case 2 will happen to some of them. Once the value of i j 1 fall into one of the first three cases, the value of i j 2 will fall into the same case. Finally the sequence of
i j will converge to a steady state. The particular structure of (38) promotes the FCSBL to converge much faster than the SBL. The analysis for the M-FCSBL follows in a similar way, and we omit it to save space. Additionally, the approximation term (28) indicates lim i j 0 , thus lim i j is the exact value of j
j
i . By taking the limit on both side of (29), we obtain
i
1
viH M
1
2
x
(41)
Applying (28) to (41) gives
viH M
1
vi viH M
1
xx H M
1
vi
(42)
Equation (42) indicates that M E xx H (in (42) x should be treated as the SR signal model). This relationship implies that the proposed approach can be utilized to get an effective CCM estimation, which is just desired by a STAP algorithm. 3.4. M-FCSBL-STAP filter Once the space-time profiles and the noise power have been estimated by the M-FCSBL, the CCM can be calculated by 15
K 2 L R 1 L i ,l vi viH 2I i 1 l 1
(43)
where i ,l denotes the element at the i th row and l th column of B , 1 is a real constant, and
2 is the estimated noise power. In (43), diagonal loading technique is adopted for avoiding distorted beam shapes and high sidelobes, and controls the loading level (satisfactory results are obtained in our experiments by considering 5 10 dB ). Another approach for calculating the CCM is K
R i vi viH 2I
(44)
i 1
where i denotes the estimate of i . Virtually, there is hardly any difference between (43) and (44) 2 L j j since lim i lim i,l L (this relationship can be easily drawn from (26) by the fact j j l 1
lim i j 0 ). Finally, the M-FCSBL-STAP weight vector is given by j
w M-FCSBL
R 1s H 1 s R s
(45)
The STAP dictionary D must be determined before applying the M-FCSBL-STAP algorithm. If N s and N d are too small, it would result in serious mismatch between the space-time steering vectors in D and those consistent with the clutter manifold. If N s and N d are too large, the computational
complexity would increase seriously. From our empirical study, it always leads to approving performance by considering N s 4 N e and N d 4 N p . Sometimes the noise power can be measured previously, for example by placing the radar in receive-only mode. If the noise power is known, the procedure for estimating it can be omitted from the M-FCSBL, and the computational load of the M-FCSBL can be further reduced. Even if the noise power is unknown, the M-FCSBL can estimate it effectively from measured data. Although the M-FCSBL-STAP algorithm is derived for ULA and constant PRF, it also suits for other array geometries and random slow-time samples. Another problem concerning the implementation of 16
the M-FCSBL-STAP is the selection of qualified secondary data. The recently proposed secondary data selection algorithm based on the clutter spectrum similarity may be a preferable candidate [2]. 4. Numerical experiments In this section, we first verify the advantages of the proposed FCSBL and M-FCSBL approaches, and then assess the performance of the M-FCSBL-STAP with both simulated and Mountain-Top data. We use the signal to interference plus noise ratio (SINR) loss [30] as a measure of performance, which is defined as the ratio of output SINR to the SNR can be achieved by a matched filter in an interference-free environment, i.e.,
LSINR
H 2 w s
2
M w H Rw
(46)
where w denotes the STAP weight vector. The simulated data are generated by considering the following scenario: a side-looking ULA with half wavelength inter-element spacing, N e 8 elements with uniform transmit pattern, wavelength is 0.23 m, PRF is 2434.8 Hz, N p 8 , 2 10 , platform height is 8 km, platform velocity is 140 m/s, N c 180 clutter patches uniformly distributed from azimuth π 2 to π 2 , the amplitude for each
clutter patch follows a complex Gaussian distribution, clutter to noise ratio (CNR) is 40 dB. 4.1 Advantages of FCSBL and M-FCSBL Fig. 1(a) provides the Capon spectrum which is calculated with the clairvoyant CCM. Figs. 1(b)–(e) illustrate the estimated clutter profiles for different approaches. Herein, 10 snapshots are used as secondary data for the M-SBL and M-FCSBL, the noise power is initialized as ˆ 2 1 , the SBL and M-SBL are terminated after 400 iterations, the FCSBL and M-FCSBL are terminated after 20 iterations. Figs. 1(b)–(e) are plotted according to the estimated hyper-parameters. Specifically, each of the estimated hyper-parameters is associated with a sampling point in the normalized space-time plane, and their values are depicted. 17
0.3 30 0.1 20 -0.1
10
-0.3
0 -0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
40 0.3 30 0.1 20 -0.1
0
-0.5 -0.5
-0.3
-0.1
0.1
40 0.3 30 0.1
20
-0.1
10
-0.3
0 -0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
0.5 40 0.3
30
0.1
20
-0.1
10
-0.3
0
-0.5 -0.5
-0.3
0.5
-0.1
0.3
0.5
(c) M-SBL 0.5 40
0.3
30 0.1 20 -0.1
10
-0.3
0
-0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(d) FCSBL
0.1
normalized Doppler frequency
(b) SBL
0.5
-0.5 -0.5
0.3
normalized Doppler frequency
(a) Clairvoyant Capon spectrum normalized spatial frequency
10
-0.3
normalized spatial frequency
-0.5 -0.5
0.5
normalized spatial frequency
40
normalized spatial frequency
normalized spatial frequency
0.5
0.5
(e) M-FCSBL
Fig. 1. Clairvoyant Capon spectrum and estimated clutter profiles
As shown in Figs. 1(b)–(e), the actual continuous clutter distribution converges into several disconnected clutter scatters for the SR approaches. It is common in the field of SR since the SR approaches do not try to recover the positions and amplitudes of clutter patches but recover the clutter subspace using as small number of space-time steering vectors as possible. We can also find that there are a few false and weak peaks induced by noise in Figs. 1(b) and 1(d); however the clutter profiles in Figs. 1(c) and 1(e) are quite favorable (the significant values are strictly located along the clutter ridge). This is because the M-SBL and M-FCSBL are more robust to noise compared with their SMV counterparts. In addition, we observe that the number of significant values in Fig. 1(e) equals to the clutter rank, which is 15 according to the Brennan’s rule [32]. This equality is in accordance with the results presented in [11]. Although Figs. 1(b)–(e) are plotted based on one arbitrary realization, the same conclusions can be drawn for any other realizations. For the SR-STAP algorithms, it is the estimated CCM that determines the clutter suppression performance. Thus it would be helpful to study the accuracy of the estimated CCM. For this purpose, we compare the Capon spectrums calculated with different CCM estimations. Particularly, the CCMs 18
0.5
normalized spatial frequency
normalized spatial frequency
are estimated based on (44) for different algorithms. The results are provided in Fig. 2.
40 0.3 30 0.1 20 -0.1
10
-0.3
0
-0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
0.5 40 0.3 30 0.1
20
-0.1
10
-0.3
0
-0.5 -0.5
-0.3
40 0.3 30 0.1
20
-0.1
10
-0.3
0 -0.3
-0.1
0.1
0.1
0.3
0.5
(b) M-SBL normalized spatial frequency
normalized spatial frequency
(a) SBL 0.5
-0.5 -0.5
-0.1
normalized Doppler frequency
0.3
normalized Doppler frequency
0.5
0.5 40
0.3
30 0.1 20 -0.1 10 -0.3
0
-0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(c) FCSBL
0.5
(d) M-FCSBL
Fig. 2. Capon spectrums calculated with different CCM estimations
By comparing Figs. 2(a)–(d) with Fig. 1(a), it is obvious that the Capon spectrums of SBL and FCSBL deviate a lot from the clairvoyant one. On the contrary, the Capon spectrums of M-SBL and M-FCSBL are very close to the clairvoyant one. By comparing Fig. 1(e) and Fig. 2(d), we can find that even though there are only 15 significant elements in the estimated clutter profile for the M-FCSBL algorithm, the Capon spectrum calculated with the estimated CCM approximates to the clairvoyant one very well. This result verifies that the proposed M-FCSBL approach can be utilized to get an effective CCM estimation. The cost function can be utilized to investigate the convergence performance of an approach. The cost functions for the SBL and M-SBL are the same as those for the FCSBL and M-FCSBL, respectively, with the latter being CMMV ln M Tr M1R ML . Since the SBL and M-SBL cannot accurately estimate the noise power, we use the known value for them. And 10 snapshots are used for the M-SBL and M-FCSBL. The results based on 100 Monte-Carlo trials are shown in Fig. 3. Some of the preceding values are very significant for the SBL and M-SBL, and these values are not shown for 19
the sake of easy comparison.
value of cost function (dB)
28
SBL M-SBL FCSBL M-FCSBL
27 26 25 24 0
100
200
300
number of iterations
400
500
Fig. 3. Value of cost function versus the number of iterations
As shown in Fig. 3, the four approaches converge to almost the same value (actually, the SBL converges to an even lower value compared with others, but this difference can be ignored since it is rather small). Observe that with as few as 20 iterations, the FCSBL and M-FCSBL converge to their steady-state values. However, the SBL and M-SBL converge quite slowly, and more than 400 iterations are needed for them to achieve their steady-state values. Thus the proposed FCSBL and M-FCSBL approaches converge much faster compared with the SBL and M-SBL approaches, respectively. The computational loads of the SBL, M-SBL, FCSBL, and M-FCSBL for a single iteration are summarized in Table 3. Herein the computational load is measured in term of the number of complex multiplications. Table 3 Comparison of computational load for a single iteration Approach SBL M-SBL FCSBL M-FCSBL
Computational load
M M M M
2 K 2 KL 2 K M KL 3K 2 6KM 3K M 1 2 KL 4 K L M 3K L
3
2 KM K 4 K 1 M K 3K 2
3
2 KM
2
3
4 KM
2
3
5KM
2
2
2
2
Fig. 4 gives a more direct illustration of the computational loads for the four approaches. It shows the computational load as a function of the system DOFs for a single iteration, furthermore it is supposed that K 16M and L 10 .
20
computational load
10
11
10
10
10
9
10
8
10
7
10
6
10
5
SBL M-SBL FCSBL M-FCSBL
0
100
200
M
300
400
500
Fig. 4. Computational load versus system DOFs for a single iteration
Clearly, for a single iteration, the FCSBL and M-FCSBL have lower computational loads than the SBL and M-SBL, respectively. In addition, as the FCSBL and M-FCSBL converge much faster than the SBL and M-SBL, respectively, the overall computational loads of the FCSBL and M-FCSBL are even lower than those of the SBL and M-SBL, respectively. 4.2 Performance of M-FCSBL-STAP A. Simulated data In Fig. 5, we assess the performance of the M-FCSBL-STAP by investigating the SINR loss in the main beam direction. The M-FCSBL-STAP is compared with the loaded sample matrix inversion (LSMI) algorithm, the SR-STAP algorithm based on M-SBL (M-SBL-STAP), and the optimum performance which is obtained using the clairvoyant CCM. The stopping criterion (b) in Section 3.1 is adopted for the two SR-STAP algorithms. To guarantee convergence, the threshold is set as 103 for the M-FCSBL-STAP while 105 for the M-SBL-STAP (a smaller implies more iterations). The loading factor for the LSMI is 10 dB to the noise power. The number of snapshots used for training is 10 for the two SR-STAP algorithms and 30 for the LSMI. All the presented results are averaged over 100 Monte Carlo trials.
21
SINR loss (dB)
0
-3
-6
-9 -0.5
LSMI M-SBL-STAP M-FCSBL-STAP Optimum
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
Fig. 5. SINR loss versus normalized Doppler frequency in main beam direction
As illustrated in Fig. 5, the LSMI approximates the -3dB loss since the number of training samples is twice the clutter rank. Note that the plain sample matrix inversion algorithm does not work for such few samples. On the other hand, the M-FCSBL-STAP can achieve comparable performance to the M-SBL-STAP, and their performances are very close to the optimum performance even though only 10 snapshots are used for training. Additionally, we found in the experiment, 480 iterations on average are used for the M-SBL-STAP, but this number is only 17 for the M-FCSBL-STAP. Fig. 6 presents the curves of average SINR loss versus the number of snapshots. Herein, the average SINR loss is defined as the mean of the SINR loss values in the entire normalized Doppler frequency range. The thresholds for the two SR-STAP algorithms are the same as above.
average SINR loss (dB)
0
-3
-6
-9 0
LSMI M-SBL-STAP M-FCSBL-STAP Optimum
10
20
30
40
number of snapshots
50
60
Fig. 6. Average SINR loss versus the number of snapshots
As depicted in Fig. 6, the LSMI converges quite slowly. On the contrary, the M-SBL-STAP and M-FCSBL-STAP achieve both rapid convergence rates and favorable steady-state performances. The results in Figs. 5 and 6 indicate that the proposed M-FCSBL-STAP can achieve comparable 22
performance to the M-SBL-STAP. However, as the M-FCSBL requires far less iterations than the M-SBL to obtain a favorable estimated clutter profile, the M-FCSBL-STAP is more suitable for application. B. Measured data Then we apply the M-FCSBL-STAP to the public available Mountain-Top data set [33]. For this data set, the receive antenna array consists of 14 elements, echoes from 16 coherent pulses are organized into a CPI data. A 500 kHz linear frequency modulated pulse was used for transmitting, the pulse width was 100 microseconds, and the PRF was 625 Hz. The data file t38pre01v1 CPI6 is used in our experiment. In this data file, there are 403 range snapshots, the clutter locates around 245° relative to true North and the target locates at 275° relative to true North with the same normalized Doppler frequency 0.25. The
normalized spatial frequency
clutter power spectrum estimated by all 403 range cells is depicted in Fig. 7. 0.5
0
0.3
-5
0.1
-10
-0.1 -15 -0.3 -0.5 -0.5
-20 -0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
Fig. 7. Estimated clutter spectrum using all 403 range snapshots
We adopt a symmetrical sliding window to detect the target. For each CUT, the eight snapshots around the CUT are treated as guard cells, and 10 snapshots out of the 20 snapshots located in the window are selected as secondary data based on the selection strategy proposed in [2]. The parameters for the M-SBL-STAP and M-FCSBL-STAP are the same as above. Fig. 8 displays the adapted patterns for the LSMI, M-SBL-STAP, and M-FCSBL-STAP at the target range cell. As expected, deep clutter notches are formed at the clutter ridge region for the M-SBL-STAP and M-FCSBL-STAP. These pattern notches will suppress the clutter at the output to well below thermal noise. However, the clutter notch 23
0
0.3
-20
0.1 -40 -0.1 -60
-0.3 -0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
0.5
-80
0.5
0
0.3
normalized spatial frequency
0.5
normalized spatial frequency
normalized spatial frequency
formed by the LSMI is not so deep.
-20
0.1 -40 -0.1 -60
-0.3 -0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(a) LSMI
0.5
-80
0.5
0
0.3
-20
0.1 -40 -0.1 -60
-0.3 -0.5 -0.5
-0.3
-0.1
0.1
0.3
normalized Doppler frequency
(b) M-SBL-STAP
0.5
-80
(c) M-FCSBL-STAP
Fig. 8. Adapted patterns at the target range cell
Fig. 9 shows the STAP output power for different algorithms in the range of 148–160 km. We can observe that peaks are formed at the target (at 154 km) range cell for all of the algorithms. However, there are some clutter residuals for the LSMI (at the few range cells around 151.8 km, 152.8 km, and 154.9 km). The target is clearly visible for both SR-STAP algorithms even though only 10 snapshots are used as secondary data. Indeed, the two SR-STAP algorithms show nearly the same output powers. However, the M-FCSBL-STAP is more favorable for practical application as it needs far fewer
output power (dB)
iterations and lower computation cost at each iteration. LSMI M-SBL-STAP M-FCSBL-STAP
80 70 60 50 148
150
152
154
range (km)
156
158
160
Fig. 9. Range plots using 10 snapshots as secondary data
5. Conclusions We derived the FCSBL approach to improve the convergence performance of the SBL approach in this paper. We also extended the proposed FCSBL approach to the MMV case. And the proposed M-FCSBL approach is utilized to devise a SR-STAP algorithm for the benefit of drastically reduced training requirement. Numerical results confirm that the M-FCSBL-STAP can achieve comparable 24
clutter suppression performance to the M-SBL-STAP. As the M-FCSBL approach enjoys several favorable properties such as user parameter free and fast converging, the resulting M-FCSBL-STAP algorithm is more preferable for application. Acknowledgment This work was supported in part by National Natural Science Foundation of China under Grant 61501506. The authors would like to thank the anonymous reviewers for their valuable comments. Appendix For completeness, the pseudo-code of the M-SBL is provided in Table 4. Table 4 Pseudo-code of M-SBL Initialization:
2
1
ˆ 2 L
i1 1 L viH x l
2
2
viH vi , ( i 1,
,K )
l 1
M k v k v kH 2 I K
1
1
1
k 1
Repeat: j
j
Ψ Γ DH j
2
M
Γ Γ DH M j
j 1
j
2
K
k k 1
j 1
1
DΓ
j
2
j
X
1 L X DΨ
i j 1 1 L Ψi,j j 1
1
j
2
M
j
j i ,i
F
K M K k 1
, ( i 1,
v k v kH 2
j 1
j k ,k
k j
,K )
I
Until a stopping criterion is satisfied
In Table 4,
j k ,k
denotes the k th diagonal elements of
j
, Ψi ,j stands for the i th row of Ψ j .
References [1]
R. Klemm, Principles of space-time adaptive processing, 3d ed., IET, London, 2006.
[2]
Y. Wu, T. Wang, J. Wu, J. Duan, Robust training samples selection algorithm based on spectral similarity for space–time adaptive processing in heterogeneous interference environments, IET Radar, Sonar and Navigation. 25
9 (7) (2015) 778–782. [3]
S. Kang, J. Ryu, J. Lee, J. Jeong, Analysis of space–time adaptive processing performance using K-means clustering algorithm for normalisation method in non-homogeneity detector process, IET Signal Processing. 5 (2) (2011) 113–120.
[4]
W. L. Melvin, A STAP overview, IEEE Aerospace and Electronic Systems Magazine. 19 (1) (2004) 19–35.
[5]
J. R. Roman, M. Rangaswamy, D. W. Davis, Q. Zhang, B. Himed, J. H. Michels, Parametric adaptive matched filter for airborne radar applications, IEEE Transactions on Aerospace and Electronic Systems. 36 (2) (2000) 677–692.
[6]
S. Burintramart, T. K. Sarkar, Y. Zhang, M. C. Wicks, Performance comparison between statistical-based and direct data domain STAPs, Digital Signal Processing. 17 (4) (2007) 737–755.
[7]
J. R. Guerci, E. J. Baranoski, Knowledge-aided adaptive radar at DARPA: an overview, IEEE Signal Processing Magazine. 23 (1) (2006) 41–50.
[8]
X. M. Zhu, J. Li, P. Stoica, Knowledge-aided space-time adaptive processing, IEEE Transactions on Aerospace and Electronic Systems. 47 (2) (2011) 1325–1336.
[9]
J. A. Tropp, S. J. Wright, Computational methods for sparse solution of linear inverse problems, Proceedings of the IEEE. 98 (6) (2010) 948–958.
[10] S. Foucart, H. Rauhut, A mathematical introduction to compressive sensing, Springer, New York, 2013. [11] Z. Yang, X. Li, H. Wang, W. Jiang, On clutter sparsity analysis in space-time adaptive processing airborne radar, IEEE Geoscience and Remote Sensing Letters. 10 (2013) 1214–1218. [12] K. Sun, H. Meng, Y. Wang, X. Wang, Direct data domain STAP using sparse representation of clutter spectrum, Signal Processing. 91 (9) (2011) 2222–2236. [13] K. Sun, H. D. Meng, F. D. Lapierre, X. Q. Wang, Registration based compensation using sparse representation in conformal array STAP, Signal Processing. 91 (10) (2011) 2268–2276. [14] Z. Ma, Y. Liu, H. Meng, X. Wang, Jointly sparse recovery of multiple snapshots in STAP, in: Proceedings of IEEE Radar Conference, April 29 2013–May 3 2013, pp. 1–4. [15] S. Sen, OFDM radar space-time adaptive processing by exploiting spatio-temporal sparsity, IEEE Transactions on Signal Processing. 61 (1) (2013) 118–130. [16] Z. Yang, X. Li, H. Wang, W. Jiang, Adaptive clutter suppression based on iterative adaptive approach for 26
airborne radar, Signal Processing. 93 (12) (2013) 3567–3577. [17] Z. Yang, X. Li, H. Wang, L. Nie, Sparsity-based space-time adaptive processing using complex-valued Homotopy technique for airborne radar, IET Signal Processing. 8 (5) (2014) 552–564. [18] Z. C. Yang, R. Fa, Y. L. Qin, X. Li, H. Q. Wang, Direct data domain sparsity-based STAP utilizing subaperture smoothing techniques, International Journal of Antennas and Propagation. 2015 (2015) 1–10. [19] Q. Wu, Y. D. Zhang, M. G. Amin, B. Himed, Space–time adaptive processing and motion parameter estimation in multistatic passive radar using sparse Bayesian learning, IEEE Transactions on Geoscience and Remote Sensing. 54 (2) (2016) 944–957. [20] M. E. Tipping, Sparse Bayesian learning and the relevance vector machine, Journal of Machine Learning Research. 1 (2001) 211–244. [21] D. P. Wipf, B. D. Rao, Sparse Bayesian learning for basis selection, IEEE Transactions on Signal Processing. 52 (8) (2004) 2153–2164. [22] D. P. Wipf, B. D. Rao, An empirical Bayesian strategy for solving the simultaneous sparse approximation problem, IEEE Transactions on Signal Processing. 55 (7) (2007) 3704–3716. [23] S. Ji, Y. Xue, L. Carin, Bayesian compressive sensing, IEEE Transactions on Signal Processing. 56 (6) (2008) 2346–2356. [24] S. Ji, D. Dunson, L. Carin, Multitask compressive sensing, IEEE Transactions on Signal Processing. 57 (1) (2009) 92–106. [25] S. D. Babacan, S. Nakajima, N. N. Do, Bayesian group-sparse modeling and variational inference, IEEE Transactions on Signal Processing. 62 (11) (2014) 2906–2921. [26] Q. Wu, Y. D. Zhang, M. G. Amin, B. Himed, Complex multitask Bayesian compressive sensing, in: Proceedings of 2014 IEEE International Conference on Acoustic, Speech and Signal Processing, 4–9 May 2014, pp. 3375– 3379. [27] M. Hurtado, C. H. Muravchik, A. Nehorai, Enhanced sparse Bayesian learning via statistical thresholding for signals in structured noise, IEEE Transactions on Signal Processing. 61 (21) (2013) 5430–5443. [28] Z. L. Zhang, B. D. Rao, Extension of SBL algorithms for the recovery of block sparse signals with intra-block correlation, IEEE Transactions on Signal Processing. 61 (8) (2013) 2009–2015. [29] J. C. Ye, J. M. Kim, Y. Bresler, Improving M-SBL for joint sparse recovery using a subspace penalty, IEEE 27
Transactions on Signal Processing. 63 (24) (2015) 6595–6605. [30] J. Ward, Space-time adaptive processing for airborne radar, Technical report 1015, MIT Lincoln Lab, Lexington, MA, USA, 1994. [31] J. Li, P. Stoica, Z. S. Wang, On robust Capon beamforming and diagonal loading, IEEE Transactions on Signal Processing. 51 (7) (2003) 1702–1715. [32] L. E. Brennan, F. M. Staudaher, Subclutter visibility demonstration, RL-TR-92-21, Adaptive Sensors, Inc., Santa Monica, CA, 1992. [33] G. W. Titi, D. F. Marshall, The ARPA/NAVY mountaintop program: adaptive signal processing for airborne early warning radar, in: Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing, May 1996, pp. 1165–1168.
Highlights
We study the problem of clutter suppression in STAP with finite training samples. (c)
Fast converging sparse Bayesian learning approaches are derived. (d)
A novel STAP algorithm named as M-FCSBL-STAP is proposed. (e)
The M-FCSBL-STAP has superior performance in low training support situation.
28