Three-dimensional sparse recovery space-time adaptive processing for airborne radar

: Conventional 3-dimensional (3D) space-time adaptive processing (STAP) has achieved good performance for non- stationary clutter suppression. However, the performance of 3D STAP rapidly degrades when applied to heterogeneous clutter environment due to the requirement of a large number of independent and identically distributed training samples to estimate the clutter covariance matrix. This study applies an efficient sparse recovery algorithm, i.e. multiple sparse Bayesian learning (MSBL), with tuning, to solve the limited sample problem of 3D STAP in airborne radar. Differing with the conventional 2D sparsity-based STAP, the proposed method utilises the sparsity of clutter in elevation-azimuth-Doppler domain and recovers the 3D clutter spectrum. Whereas in its’ large computational load, a fast algorithm extending the relevance vector machine is also considered.


Introduction
As an efficient technique for clutter suppression in airborne radar, space-time adaptive processing (STAP) has been put forward and developed over 40 years [1]. As it is well known, the performance of STAP filter is dependent on the estimation accuracy of clutter covariance matrix (CCM) corresponding to the range cell under test (CUT). For classical STAP methods, the CCM can be obtained using the training samples via maximum likelihood estimate. According to the Reed-Mallett-Brennan rule [2], full-dimension STAP requires at least twice of the system degrees of freedom (DOFs) of the independent and identically distributed (i.i.d.) training samples to achieve the steady performance. However, it is hardly met in practical non-homogeneous clutter environment which results the degraded performance of the adaptive processor compared to theoretical predictions. Clutter non-homogeneity mainly includes two kinds [3]: clutter heterogeneity and geometryinduced angle-Doppler variation (i.e. clutter non-stationarity).
So far, much effort has been put into clutter non-stationarity problem. Compensation methods [4][5][6] work well when range ambiguity does not exist. Once the radar works at the medium or high pulse repletion frequency (PRF), it cannot identically compensate the ambiguity clutter. In fact, the short-range clutter which is induced from the array elevation sidelobe is the dominant component contributed to clutter non-stationarity. The threedimensional (3D) STAP [7] with elevation elements or sub-arrays can cancel the short-range clutter effectively in theory. Unfortunately, it is hard to get enough i.i.d. sample data and requires more computation load because of its large system DOFs. Sub-array synthesis algorithm with non-adaptive pre-filtering in elevation [8] can effectively eliminate the short-range clutter when the element errors are not existed. However, the array antenna errors are inevitable in practice, thereby significantly degrading the performance. To overcome this problem, an adaptive pre-filter in elevation for eliminating the short-range clutter following is proposed in [9]. Its performance is not near optimal because only two system DOFs are available in elevation. Cross beam forming (CBF) STAP [10] utilises the multi-beams along the cross-line to cancel the clutter induced from azimuth and elevation sidelobe simultaneously. Although it obtains robust performance for nonstationary clutter suppression to some degree, part performance loss is resulted by the deficiency of the chosen beams. Most importantly, above methods are not low-sample support and, therefore, are unavailable to heterogeneous clutter situation.
In real clutter environment, the heterogeneity and nonstationarity often exist together, and hence robust STAP method suited for both heterogeneous and non-stationary clutter background is needed. In this paper, a 3D multiple sparse Bayesian learning (MSBL) STAP method, which can effectively eliminate the non-stationary clutter no matter what the PRF is and only needs a few training samples, is proposed and discussed.

Characteristics of non-stationary clutter
In this section, the range dependence of clutter for non-side looking airborne radar (non-SLAR) is discussed in detail. Consider a pulse Doppler airborne radar with a uniform planar array consisting of M and N elements in elevation and azimuth, respectively, as shown in Fig. 1. The interval between two elements is equal to half wavelength. The x-axis is aligned with the normal of the array antenna, and the z-axis points vertically up. The counterclockwise angle between the flight direction and y-axis is θ a ; the angles θ, φ, and ψ denote the azimuth angle, elevation angle, and cone angle, respectively. The platform is at altitude H and moving with constant velocity v. R c denotes the slant range from the clutter patch to array antenna. Assume that a CPI contains K pulses.
According to Fig. 1, the relationship between the spatial frequency (i.e. cosine of cone angles) and the Doppler frequency can be described as the following well-known formula [4] where denotes radar wavelength. When θ a ≠ 0, i.e. the array antenna is mounted as a non-SLAR, clutter Doppler varies with the slant range and thus clutter is non-stationary. A non-SLAR is considered in the paper and the system parameters are shown in Table 1.
The clutter spectrum distribution of non-SLAR is depicted in Fig. 2, where Figs. 2a-c display the corresponding spectrum in azimuth-Doppler plane, elevation-Doppler plane, and range-Doppler plane, respectively. Fig. 2a shows that the clutter in azimuth-Doppler domain consists of the short-range clutter and the range ambiguity clutter. It is evident that the short-range clutter and range ambiguity clutter cannot be distinguished both in azimuth domain and Doppler domain. However, from Fig. 2b, we can see that the short-range clutter and the range ambiguity clutter can be J. Eng distinguished along the elevation domain. As shown in Fig. 2c, the curved clutter (i.e. non-stationary clutter component) is mainly composed by range-short clutter, and the range ambiguity clutter is nearly stationary. According to the conventional planar array pattern, we know that the non-stationary clutter is induced from elevation sidelobe. Therefore, 3D STAP, i.e. elevation-azimuth-Doppler STAP, can effectively cancel the non-stationary clutter component in theory.

Space-time sparse signal model
The received echoes of each range cell in airborne radar can be modelled as the superposition of echo signals corresponding to all spatial frequency and Doppler frequency. Considering a multichannel-phased array radar with M received channels in elevation, N received channels in azimuth, and K coherent pulses, we discretise the whole elevation-azimuth-Doppler plane uniformly into sufficiently small grid points, e.g. N e = ρ e M elevation angle bins, N a = ρ a N azimuth angle bins and N d = ρ d K Doppler bins, where ρ e , ρ a , and ρ d determine the smoothness of the elevationazimuth-Doppler images. The corresponding sets of space-time steering vectors of all grid points are formulated as where ψ ∈ ℂ MNK × N e N a N d is also called space-time steering dictionary, '⊗' denotes Kronecker product, 1 , s e, 2 , …s e, N e ] ∈ ℂ M × N e , , s a, n , and s t, k are the spatial and temporal steering vectors, given by where () T denotes the transposition operation, f e, m is the spatial frequency of the mth elevation angle grid point, f a, n is the spatial frequency of the nth azimuth angle grid point, and f d, k is the Doppler frequency of the kth Doppler grid point. The received clutter plus noise snapshots of L range cells (1) , …, x (L) ] ∈ ℂ MNK × L can be expressed by where A = [α (1) , …, α (L) ] ∈ ℝ N e N a N d × L is an unknown clutter source matrix (or called a solution matrix) with each row representing a possible clutter patch, N = [n (1) , …, n (L) ] ∈ ℂ MNK × L is a Gaussian white thermal noise matrix, whose elements are i.i.d. and have the   variance σ 2 . A key assumption in the multiple measurement vectors (MMV) model is that the support (i.e. indexes of non-zero entries) of every column in A is identical, which means there are the same clutter degrees of sparsity (DOSs) among the training samples. When the noise is taken into consideration, the sparsest representation in MMV is the solution to the following optimisation problems where ∥ ⋅ ∥ 2, 0 is a mixed norm defined as the sum of the l 2norms of the row-vectors, ∥ ⋅ ∥ F is the Frobenius norm (Euclidean norm of a matrix), r s ∈ ℝ + is the DOSs, and λ ∈ ℝ + is the so-called regularisation parameter, which provides a trade-off between the clutter sparsity and the output residual clutter power.
Of course, problems (P1) and (P2) are equivalent. The position and power of clutter source in elevation-azimuth-Doppler plane can be obtained by solving the problems, and then the space-time CCM of the lth range cell is calculated by where () H denotes the conjugate transposition operation, ∑ = diag(β (l) ) with diagonal elements representing the power of the clutter sources, β (l) = [β 1, 1, 1 (l) , β 1, 1, 2 (l) , …, β 1, 2, 1 (l) , β 1, 2, 2 (l) , …, β N e , N a , N d ] T and β m, n, k Under the criterion of linearly constrained minimum variance, the optimal STAP weight vector of the lth range cell is where is the normalisation constant, s 0 = s t0 ⊗ s a0 ⊗ s e0 is spacetime steering vector of the target.

Algorithm description of 3D MSBL
We first construct a convex relaxation of (8) and attempt to solve the MMV optimisation problem From a Bayesian perspective, (11) is equivalent to MAP using the Laplacian prior on the l 2 norm of each row by applying an exponent transformation [11]. Analogously, a higher sparse prior, e.g. the Jeffreys prior, is also chosen for the row norms in [12], which is easily trapped into local optima. The MSBL strategy assumes a prior distribution modulated by a vectors of hyperparameters controlling the prior variance of each row of A, the values of which can be learned from the measured clutter data.
In the presence of noise, the MSBL has been demonstrated a robust, sparse enough, parameter-blind algorithm, with no local minima [13]. Therefore, we mainly focus on exploiting MSBL theory to improve the ability for clutter suppression. Assuming the noise is white complex Gaussian distribution with unknown power σ 2 , we can obtain the Gaussian likelihood function of measurement model in (6) as Suppose each column of A is a complex Gaussian prior a (l) ∼ N (0, Γ) (13) where 0 ∈ ℂ MNK × 1 is a zero vector, Γ = diag(γ), are unknown variance parameters corresponding to all grid point in angel-Doppler plane and can be viewed as the power β of the clutter sources. By combining each column prior, we obtain a full prior of A By combining likelihood and prior, the posterior density of A then becomes (15) which is modulated by the hyperparameter vectors γ and σ 2 . Actually, estimating the sparsity profile of A is conveniently shifted to estimating the hyperparameter vector γ with the correct number, value, and location of non-zero elements. The latter can be utilised to estimate the CCM and can be effectively accomplished via expectation maximisation (EM) algorithm [14]. According to EM algorithm, the hyperparameter vectors γ and σ 2 of the q + 1 step can be calculated as where The iteration procedure for updating γ and σ 2 will not stop until a predefined criterion is satisfied. For example, if the q + 1 iteration satisfies ∥ γ q + 1 − γ q ∥ 2 / ∥ γ q ∥ 2 ≤ δ or σ 2 ≤ (σ*) 2 where δ and (σ*) 2 are small enough positive threshold, then the iteration is over. While the iterative algorithm described above has been demonstrated to yield a highly accurate sparse linear-regression representation, the practical limitation is observed when applied to large-scale problems, e.g. 3D situation. When evaluating (18) and (19), one must multiply matrices of size (N e N a N d ) × (MNK), an O((N e N a N d ) 2 (MNK)) operation , thereby making this approach relatively slow. There we apply a fast algorithm [15] and reduce the computation complexity to O(L(N e N a N d )r s 2 ) where r s is the DOSs and has r s ≪ N e N a N d .

Simulations
In this section, we assess the proposed 3D MSBL method using simulated radar data. The parameters of the simulated radar platform are shown in Table 1. In the following examples, we assume the grid points ρ e = 2, ρ a = 4 and ρ d = 4, the stopping criterion of the iteration procedure is (σ*) 2 = 10 −20 .

3D sparse spectrum
In the first experiment, we examined the sparse clutter spectrum of 3D MSBL STAP via the comparison with the conventional capon clutter spectrum, in which the CCM is assumed to be exactly known. We chose the 200th range cell as the CUT and only four training samples are utilised for sparse recovery. Figs. 3a and b show the capon clutter spectrum in elevation-azimuth-Doppler domain and two slices along the elevation angles corresponding to short range and mainbeam, respectively. Similarly, the sparse clutter spectrum is shown in Fig. 4. From Fig. 4, we see that the location and power of clutter can be well recovered via 3D MSBL compared with the capon clutter spectrum in Fig. 3. Meanwhile, we also see that favourable sparse properties in elevation-azimuth-Doppler domain are displayed, and the short-range clutter and mainbeam clutter are clearly distinguished along the elevation direction.

Comparison with other STAP algorithms
The 200th range cell, in which the non-stationary clutter is severe, is also chosen to be the CUT in this experiment. Four correlative STAP methods are considered for comparison. The conventional 2D loaded sample matrix inversion (LSMI) [16] method with 11 spatial channels in azimuth is first considered. The second is the 3D LSMI [7] with synthesised four spatial channels in elevation and eight spatial channels in azimuth. Then, the STAP method based on adaptive spatial pre-filtering in elevation [9] is investigated. Finally, the recent CBF STAP [10] with three beams and five beams in elevation and azimuth, respectively, is also taken into account. We measure the improvement factor (IF) curve, where the IF is defined as follows.
s 0 is the space-time adaptive weight, R and R are the estimated and exact CCM, respectively, and tr(·) is the trace of a matrix.
In the first part, as shown in Fig. 5, we present the IF performance against the normalised Doppler frequency for our proposed method and above referred methods using four training samples. The performance of optimum 3D STAP is also shown here to give the upper bound. The plots show that the proposed method provides a suboptimum IF performance, but it outperforms the other STAP methods and no obvious performance loss in the short-range clutter Doppler bins. The results shown in Fig. 5 also display that, compared with existing STAP methods, the 3D MSBL STAP can achieve good performance for non-stationary clutter cancellation with only a few training samples, which means it also available to the heterogeneous clutter background.
In the second part, we pay our attention to the convergence speed of our proposed method. As shown in Fig. 6, the plots depict that the convergence speed of our proposed method distinctly outperform the other methods. Although the IF performance of 3D LSMI STAP is better than our proposed method when the number of the available training samples is >30, the IF steady-state performance of our proposed method is obtained especially when the training samples is insufficient.

Conclusion
In this paper, we have applied the MSBL algorithm to 3D STAP for airborne radar system to solve the clutter suppression problem in non-stationary and heterogeneous clutter background. The improved 3D STAP method exploits the sparsity in elevationazimuth-Doppler domain and imposes the separableness of shortrange clutter and ambiguity clutter in elevation. To overcome the large computation problem, we utilised the fast algorithm to estimate the hyperparameters. The simulation results have shown that the proposed method outperforms the conventional STAP   methods for non-stationary clutter elimination with a few training samples, which means that the method proposed in this paper can work well in non-stationary and heterogeneous clutter environment. However, the problem of off-grid clutter in elevationazimuth-Doppler domain is not mentioned in this paper. In fact, the pixels corresponding to clutter are not always perfectly aligned on the clutter ridge, especially when the number of antennas does not equal the number of pulses. Therefore, further work should be devoted to the off-grid problem of 3D MSBL STAP.

Acknowledgments
This work was supported by the China National Science Foundation (Grant No. 61501506).