Random Nyström Approach to Clutter Subspace Estimation for Polarimetric Space-Time Adaptive Processing

: A random Nyström (R-Nyström) scheme for clutter subspace estimation is proposed in the context of polarimetric space-time adaptive processing (pSTAP). Unlike the standard Nyström scheme making use of only partial columns of the clutter plus noise covariance matrix (CNCM), R-Nyström exploits full CNCM information with a properly designed selection procedure under the newly developed random ridge cross leverage score (RRCLS) criterion. With R-Nyström, sup-ported by the complete CNCM columns, upgraded clutter subspace estimation can be achieved at the expense of an insigniﬁcant increase in computational complexity, in contrast to the standard Nyström. The R-Nyström-based pSTAP, termed pR-Nyström, is shown to be superior over the current eigendecomposition-free subspace pSTAP in the signal to clutter plus noise loss and computational complexity. The e ﬃ cacy of R-Nyström / pR-Nyström is validated by the simulation results.


Introduction
Polarimetric space-time adaptive processing (pSTAP) plays an important role in improving the detection performance of airborne array radar operating in clutter environments [1][2][3][4]. The purpose of pSTAP is to suppress clutter by maximizing the signal to clutter plus noise ratio (SCNR) [5] in the cell under test (CUT), with the direct use of the clutter plus noise covariance matrix (CNCM). In pSTAP, the estimation accuracy of CNCM determines the performance of the processor. In general, CNCM can be estimated by using the independent and identically distributed (i.i.d.) training samples obtained from the adjacent range bins of CUT. However, in a real nonhomogeneous environment, the clutter characteristic in different range cells would change quickly and dynamically, which caused difficulty in collecting enough training samples [6].
Subspace-based methods have also been studied in the literature by using the principal subspace of CNCM instead. For example, the eigen-canceller method based on the direct eigenvalue decomposition (EVD) of CNCM was developed in [7]. This method, however, is not attractive for real-time processing due to its high computational complexity. For this reason, some fast subspace methods without EVD of CNCM have been proposed further. For instance, the projection approximation subspace tracking (PAST) scheme in [8] and the fast approximated power iteration (FAPI) scheme in [9]. Both schemes are implemented iteratively for the principal subspace estimation with some training samples available. A variant of PAST, called L1-PAST, was developed in [10], and incorporates an additional sparse constraint on the filter. These mentioned EVD-free subspace methods are computationally more efficient than the EVD method but suffer a severe performance loss in the estimation of the principal subspace.
The Nyström scheme is a popular technique for the fast estimation of the matrix representation of the principal subspace [11][12][13][14] while it was originally designed for solving the integral equations [15].
The main idea of Nyström is to select a subset called landmark points for the calculation of the kernel eigenfunction approximately. The most favored sampling method for subset selection is random sampling, which leads to the fast versions of kernel machines [16] and spectral clustering [17].
Although the random sampling method can work in practice, it will fail for many natural kernel matrices where the relative "importance" of points is not uniform across the dataset. A ridge leverage score (RLS) method is proposed for measuring the importance of points [18][19][20][21] such that the kernel matrices can be selected more reasonably. When the number of selected columns is fixed, the corresponding number of rows of the selection matrix is also fixed. The performance of the existing Nyström-based subspace estimation method is proportional to the number of selected columns. Therefore, it is worth studying improving the estimation accuracy of the principal subspace without increasing the selection number [22].
Most current Nyström variants employ only partial column information of CNCM. There is still room to improve the application of Nyström for estimating the clutter subspace. In this paper, we propose a new method called random Nyström (R-Nyström) which makes full use of CNCM under the newly developed "random" ridge cross leverage score (RRCLS) rule. In the RRCLS rule, a new selection matrix is constructed for obtaining the low-dimensional matrix in Nyström. The R-Nyström scheme can be directly applied for implementing high computational efficient pSTAP.
The major contributions of the paper are the following: (1) Enabling fast principal subspace identification within the Nyström framework while making full use of the CNCM information. (2) The computational burden in pSTAP can be reduced with high clutter suppression performance compared with the existing methods.

Data Model
Consider a side-looking uniform linear polarimetric radar array system with N elements, each of which contains two independent polarimetric receiving channels. The inner-element spacing of the radar array d is half of the radar wavelength. The geometry of the system is shown in Figure 1, where (θ, φ) is the beam direction, and v is the platform velocity. The polarimetric radar transmits uniformly M pulses. The pulse repetition frequency (PRF) in a coherent processing interval (CPI) is ι. methods are computationally more efficient than the EVD method but suffer a severe performance loss in the estimation of the principal subspace.
The Nyström scheme is a popular technique for the fast estimation of the matrix representation of the principal subspace [11][12][13][14] while it was originally designed for solving the integral equations [15]. The main idea of Nyström is to select a subset called landmark points for the calculation of the kernel eigenfunction approximately. The most favored sampling method for subset selection is random sampling, which leads to the fast versions of kernel machines [16] and spectral clustering [17].
Although the random sampling method can work in practice, it will fail for many natural kernel matrices where the relative "importance" of points is not uniform across the dataset. A ridge leverage score (RLS) method is proposed for measuring the importance of points [18][19][20][21] such that the kernel matrices can be selected more reasonably. When the number of selected columns is fixed, the corresponding number of rows of the selection matrix is also fixed. The performance of the existing Nyström-based subspace estimation method is proportional to the number of selected columns. Therefore, it is worth studying improving the estimation accuracy of the principal subspace without increasing the selection number [22].
Most current Nyström variants employ only partial column information of CNCM. There is still room to improve the application of Nyström for estimating the clutter subspace. In this paper, we propose a new method called random Nyström (R-Nyström) which makes full use of CNCM under the newly developed "random" ridge cross leverage score (RRCLS) rule. In the RRCLS rule, a new selection matrix is constructed for obtaining the low-dimensional matrix in Nyström. The R-Nyström scheme can be directly applied for implementing high computational efficient pSTAP.
The major contributions of the paper are the following: (1) Enabling fast principal subspace identification within the Nyström framework while making full use of the CNCM information.
(2) The computational burden in pSTAP can be reduced with high clutter suppression performance compared with the existing methods.

Data Model
Consider a side-looking uniform linear polarimetric radar array system with N elements, each of which contains two independent polarimetric receiving channels. The inner-element spacing of the radar array d is half of the radar wavelength. The geometry of the system is shown in Figure 1   The radar data vector under the target free hypothesis H 0 and the target presence hypothesis H 1 can be written by H 0 : x = u H 1 : where u = c + n is the clutter plus noise vector which is, as usual, a random vector drawn from a complex Gaussian process with zero-mean and unknown positive definite covariance matrix R; t = Ap Electronics 2020, 9,1488 3 of 13 is the target polarization-space-time steering vector (or the target vector), p is a complex random vector characterizing the unknown amplitude, phase, and polarization parameters of target return, and A = I 2 ⊗ d f D ⊗ s f S , f D = 2v cos θ cos ϕ/(λι) and f S = d cos θ cos ϕ/λ denote the target normalized Doppler frequency and normalized spatial frequency, respectively; where ⊗ denotes the Kronecker product, I m denotes the m × m identity matrix, λ is the radar wavelength, and where superscript "T" denotes transpose. Following Ward's clutter model [5] and polarization characterize [1], the 2MN × 1 vector c can be written as follows where N C is the number of the independent clutter patches in one range cell, the MN × 1 zero-mean complex Gaussian random vectors c H and c V , given by represent the horizontal and vertical clutter components whose covariance matrices are in which E{·} denotes expectation and superscript "H" denotes conjugate transpose; f D,n and f S,n are, respectively, the Doppler frequency and the spatial frequency of the nth clutter patch, The 2 × 1 vector p n can be written as p n = P 1/2 ρ n , where E ρ n ρ H n = σ 2 H,n I 2 , and in which ζ and γ are the correlation coefficient and the power ratio between the two polarimetric channels.

Clutter Subspace Based pSTAP
The clutter covariance matrix is According to [23], the effective rank of R HH is N + β(M − 1) , where . denotes rounding to the nearest integer and β = f D / f S = 2v/(dι) denotes the slope of the clutter ridge. From (12), the effective rank of R c is rank(R c ) = rank(P) × rank(R HH ) Because the rank of P is 2, the effective rank of the clutter covariance matrix is As the clutter and noise are independent, we have Let U c = [v 1 , · · · , v Q ] whose column span is the so-called clutter subspace, Λ c = diag(γ 1 , · · · , γ Q ), which consists of the Q principal eigenvalues of R, in which diag(·) denotes a diagonal matrix; and σ 2 is the variance of the Gaussian white noise.
The optimal pSTAP filter vector is given by [3], Applying the matrix inversion lemma to (15) yields Since the diagonal elements in Σ are much larger than σ 2 , (17) can be approximated as Therefore, the pSTAP filter vector can be written by In the following, we give a new method for estimating U c and p.

Revisit of Standard Nyström
The CNCM R can be approximated as follows Electronics 2020, 9, 1488 5 of 13 where U is an 2MN × Q orthogonal matrix whose columns are the Q principal eigenvectors of R, Λ is a diagonal matrix whose diagonal elements are the Q principal eigenvalues of R; U 11 ∈ C K×Q , and where S is a 2MN × K selection matrix, and S 11 ∈ C K×K , and S 21 ∈ C (2MN−K)×K . In the standard Nyström, S is composed of K different but arbitrary columns of I 2MN . Define thenÛ 11 , in whichÛ 11 is an K × Q orthogonal matrix whose columns are the eigenvectors ofĈ 11 , andΛ 11 is the Q × Q diagonal matrix whose diagonal elements are the eigenvalues ofĈ 11 .
It can be verified that rank(Û) = Q, and It follows thatÛ This implies that span(Û) = span(U), which is exactly the Q-dimensional clutter subspace to be identified.
The computational complexity of the standard Nyström is about O(2MNK 2 ). If K << 2MN, it is computationally much simpler than the classical eigendecomposition method whose computational complexity is about O((2MN) 3 ). However, making use of a relatively small K in the method would result in poor identification span(U) since only K columns of R are employed. In the sequel, a R-Nyström method is developed which fulfills a good compromise between the computational complexity and identification performance. In the new method, the dimension of the selection matrix S is still of 2MN × K. However, each column of S is designed by a new criterion based on R. Note that all Nyström-based methods require EVD operation, except that R-Nyström is based on a low-dimensional matrix for EVD.

The Random Nyström Method
To show the estimated performance of the clutter subspace from R, a cost function (projection error) which denotes the difference between the estimated and the real projection operator of the clutter subspace is given. The smaller the function, the more accurate the estimation. The cost function is given by where U e and U r denote the estimated and real matrix representation of the clutter subspace, respectively; ||·|| F denotes the Frobenius norm. Because the estimated matrix U e is not a column orthogonal matrix, the projection operator is used by U e (U H e U e ) −1 U H e in (25). Observed from the standard Nyström method in part 4.1, U e is determined by the selection matrix S. Therefore, f (U e ) is also determined by S. In the following, a new criterion is proposed to determine the selection matrix S.
Define a 2MN × K matrix S, and assume all the elements in S are determined with B(i, j) as the probability distribution function, namely the value at (i, j) depends on B(i, j). Let S(i, j) be the function at (i, j) of matrix S whose definition is represented as follows [24] where E(i, j) is a uniformly distributed random number between [0,1]. The derivation of B(i, j) is given in the following. First, perform singular value decomposition (SVD) for R, The jth column of R can be represented as follows where v q,j is the jth element of the qth right singular vector, and " * " denotes a conjugate. Observed from (28), R j is a linear combination of the left singular vectors and right singular values, and the coefficients are the corresponding jth row of V . Hence, R j can be approximated with a linear combination of the larger K left singular vectors and its corresponding singular values All the coefficients are extracted and used to construct the normalized statistical leverage score, which is represented as follows where ε j > 0, 2MN j=1 ε j = 1 and 1 ≤ j ≤ 2MN. The whole leverage scores form a probability distribution over the 2MN columns [18]. A new RRCLS based on related works on RLS [19][20][21][22] can be defined by in which D is a K × 2MN matrix whose elements follow a standard normal distribution, κ is a constant and 1 ≤ i ≤ 2MN, 1 ≤ j ≤ K. B(i, j) can be obtained by where δ is an adjustable factor. After obtaining the probability distribution function B(i, j), the selection matrix S can be determined by (26). Compared to the standard Nyström, the number of non-zero elements in each column of S is more than one. The selected low dimension matrix can be represented as Similar to (22), define where U 21 = R 21 U 11 Λ −1 11 , in which U 11 is an K × Q orthonormal matrix whose columns are the eigenvectors of R 11 , and Λ 11 is the Q × Q diagonal matrix whose diagonal elements are the eigenvalues of R 11 . Note that all the columns of U are not orthogonal to each other. The R-Nyström algorithm for the matrix representation of the clutter subspace is summarized as follows: (a) determine the selection matrix S by (26), (31) and (32) (b) obtain the low-dimensional matrix by (33) (c) estimate the matrix representation of the clutter subspace by (34)

Application of R-Nyström to pSTAP
The CNCM R is given by the sample estimate where u l is the clutter plus noise vector of the lth range cell and L is the number of training samples. Substitute (35) into (31), the new RRCLS can be written as By (36), the corresponding B(i, j) can be written as Therefore, the new selection matrix S is calculated by Because there is an error in the sample estimate, x in H 0 is processed directly for avoiding estimation error into the following calculation. The selected output is written by and the corresponding low-dimensional matrices are given by Then, performing EVD of R yy yields, where λ i (λ 1 ≥ · · · ≥ λ Q > λ Q+1 = · · · = λ K ) is the eigenvalue of R yy and u i is the corresponding eigenvector. Let U yy = [u 1 , u 2 , · · · , u Q ] and Λ yy = diag(λ 1 , λ 2 , · · · , λ Q ), and the matrix representation of the clutter subspace is given by The columns of U are also not mutually orthogonal. Based on (19), the weight vector of the pR-Nyström method is given by It notes that the polarimetric signal vector to be detected is unknown in advance, hence, p needs to be estimated by the proposed modified least square method (mLSM). In hypothesis H 1 , the radar data are represented as follows Introduce function f (p), which is written by where x = x − c denotes the data that contain no clutter components. For eliminating the clutter effect in x, x can be obtained by projecting x into the orthogonal subspace of clutter subspace, where the matrix representation of the clutter subspaceÛ can be obtained by the standard Nyström method, Substitute (47) into (46), Therefore, whenp = p, f (p) can get the minimum value. Let Taking the sub-gradient of (49) with respect to p,p can be obtained as followŝ The final filter vector is rewritten aŝ Extend the existing methods to the pSTAP case, which are called pEVD, pPAST, pFAPI, pL1-PAST, and pBiPerEFA, respectively. Note that all the methods use the estimatedp in the process of calculating the corresponding filter.

Computational Complexity Analysis
The major computational burden of the pR-Nyström method is engaged with the calculation for the selection matrix and the polarimetric signal vector to be detected. The main polarimetric radar system parameters are as shown in Table 1, and are the same as that in [2]. The approximate computational complexity and the averaged computation time of 100 independent runs are listed in Table 2, where the number of training samples is L = 130 and K = Q under the following configurations: (1) CPU: Intel(R) Core (TM) i7-6700 3.40 GHz; (2) memory: 8.00 GB; (3) system: 64-bit Windows 7; (4) MATLAB R2016a.
It can be seen from Table 2 that the average computation time of the mLSM method is denoted as t 1 , and the average computation time of the pR-Nyström method is less than that of the other existing methods except for the pNyström method and pBi-PerEFA method. The EFA referred to in this paper is with three Doppler bins and the iteration number of pBiPerEFA method [25] is nine.

Numerical Results
In this section, to assess SCNR loss, simulated data are employed as the metric for comparison among different methods, which is defined by where SCNR o and SCNR i denote the output and input signal to clutter plus noise ratio, respectively; R real represents the true CNCM which is assumed to be known in simulation experiments andt = Ap. The polarimetric radar system parameters are shown in Table 1, and κ = 1. The curves shown are the averaged results of 100 Monte Carlo trials. Shown in Figure 2 are the projection error curves versus training samples with different values of δ, where the target normalized Doppler frequency is 0.2. It is observed that the pR-Nyström method is robust for different values of δ between [10,40]. In the following simulation, we set δ = 20. Shown in Figure 3 are the projection error curves versus training samples, where the target normalized Doppler frequency is 0.2. It is observed that the pR-Nyström method has a smaller projection error and is closest to the pEVD method. This indicates that the pR-Nyström method has a better approximate performance than the other methods. Furthermore, in the case of the non-side-looking model, assuming that the angle between the array pointing and the direction of movement is 30 • , the SCNR loss is a little weaker than the side-looking model. Shown in Figure 4 are the SCNR loss curves versus training samples where the target normalized Doppler frequency is 0.2. It is observed that the pR-Nyström method has a better SCNR loss performance than the other methods with the increase in the training sample and is closest to the pEVD method.
Electronics 2020, 9, x FOR PEER REVIEW 11 of 14 its depth depends on the CNR. The clutter in the side lobe clutter region can be well suppressed to the noise level and is less affected by the CNR.   Electronics 2020, 9, x FOR PEER REVIEW 11 of 14 its depth depends on the CNR. The clutter in the side lobe clutter region can be well suppressed to the noise level and is less affected by the CNR.   its depth depends on the CNR. The clutter in the side lobe clutter region can be well suppressed to the noise level and is less affected by the CNR.    Shown in Figure 5 are the SCNR loss curves versus the target normalized Doppler frequency, where the number of training samples is 130. It is observed that the pR-Nyström method can achieve a better SCNR loss performance than the other methods with different target normalized Doppler frequencies and is closest to the pEVD method.  Shown in Figure 6 are the SCNR loss curves versus the normalized Doppler frequency, with different clutter degree of polarization (DOP). The number of training samples is 130, and CNR = 40 dB. It is observed that when the target Doppler frequency is large, the SNCR loss performance of the pR-Nyström method stays stable with different DOPs; when the normalized Doppler frequency value is located in [−0.1, 0.1], the SCNR loss performance of pR-Nyström method changes. When the clutter DOP is large, the clutter is a completely polarized wave and this will make a larger notch [2]. When the clutter DOP is gradually reduced, the clutter presents a partially polarized wave and this leads to deeper and deeper notch.   Shown in Figure 7 are the SCNR loss curves versus the normalized Doppler frequency, with different CNRs. The number of training samples is 130, and |ζ|= 0.99 . It is observed that the pR-Nyström method has a stable SCNR loss performance when the normalized Doppler frequency value is large.
When the normalized Doppler frequency value is located in [−0.1, 0.1], the SCNR loss performance of the pR-Nyström method changes with different CNRs. The reason is that when the clutter notch is formed in the main lobe clutter region of Doppler frequency by the clutter suppressor, its depth depends on the CNR. The clutter in the side lobe clutter region can be well suppressed to the noise level and is less affected by the CNR.

Conclusions
We have presented a random Nyström framework based on the matrix representation of the clutter subspace estimation method for pSTAP. The new method involves a new selection matrix for the complete use of CNCM, which is determined by the proposed random ridge cross leverage score SCNR loss(dB) Figure 7. SCNR loss performance versus target Doppler frequency: with different clutter to noise ratio (CNR).

Conclusions
We have presented a random Nyström framework based on the matrix representation of the clutter subspace estimation method for pSTAP. The new method involves a new selection matrix for the complete use of CNCM, which is determined by the proposed random ridge cross leverage score rule. Simulation results show that the proposed method has an appealing clutter suppression performance and computational efficiency as compared with the existing methods.
It can be seen from the investigation in this article that high subspace estimation accuracy and computational efficiency can be achieved based on the new column selection criterion within the Nyström framework. In order to reduce the influence of adjustable factors on subspace estimation, the number of adjustable factors in the column selection criterion should be as small as possible. Hence, it is important to investigate further the adjustable factors in the subsequent research. When the newly proposed method is applied in airborne radar, the clutter suppression performance of pSTAP can be improved, and the target detection probability can be increased.