Robust Space-Time Adaptive Processing Method for GNSS Receivers in Coherent Signal Environments

: In the coherent signal environments caused by multipath propagation, the interference suppression performance of the global navigation satellite systems (GNSS) receivers decreases sharply. In this paper, a robust space-time adaptive processing (STAP) method for GNSS receivers is proposed to suppress interferences in coherent signal environments, by using the modiﬁed space-time two-dimensional iterative adaptive approach (ST2D-IAA) spectrum estimation. This method applies the IAA algorithm to the ST2D signal model of GNSS receivers, and further modiﬁes the ST2D-IAA algorithm to accurately estimate the power spectrum and noise power simultaneously. The space-time interference-plus-noise covariance matrix (STINCM) is reconstructed by using the estimated power spectrum and noise power in the interference angle region. Based on the reconstructed STINCM, we construct the STAP beamforming optimization problem for the space-time steering vector (STSV) error vector, and further correct the STSV of GNSS signal. Finally, the weight vector of STAP beamforming is calculated by using the reconstructed STINCM and the corrected STSV of GNSS signal. Simulation results show that the proposed method can suppress interferences in coherent signal environments and outperforms the current methods.


Introduction
The global navigation satellite systems (GNSS) use the advantages of satellites to provide real-time, high-precision navigation, positioning and timing services, which play a vital role in the economic and military fields [1][2][3]. The current GNSS mainly includes the United States global positioning systems (GPS), China's Beidou Systems, Europe's GALILEO systems and Russia's GLONASS systems. However, the accuracy and reliability of GNSS receivers can be easily affected by all kinds of intentional and unintentional interferences [4]. With the significant increase in the types and quantities of interference, GNSS applications have posed a great threat, especially the rapid increase in cheap jamming devices [5]. Therefore, it is urgent to improve the interference suppression ability of GNSS receivers. An effective method to suppress interferences is to use antenna array processing technology, which can form a nulling in the direction of interference [6].
A large number of interference suppression algorithms for GNSS receivers based on antenna arrays have been proposed. Among these algorithms, space only processing (SOP) and space-time adaptive processing (STAP) are two typical representatives [7]. The SOP algorithms [8] use only one spatial filter to process the observed signals, mainly divided into power inversion algorithms and adaptive beamforming algorithms. The power inversion algorithms [9] do not need to know the array structure and satellite signal information in advance, and utilizes the characteristics that the intensity of the satellite signal is lower than the noise and interferences when it reaches the receivers. However, this kind of algorithms cannot guarantee that the mainlobe of the beampattern always points to the satellite signal.
Adaptive beamforming algorithms are proposed to enhance desired satellite signals while suppressing interference signals which are based on the criteria of minimum variance distortionless response (MVDR), minimum power distortionless response (MPDR) and maximum signal-to-interference-plus-noise ratio (SINR). In practice, there may be direction errors, local scattering, wavefront distortion and limited number of sampling snapshots, which lead to the covariance matrix mismatches and steering vector errors. In order to improve the beamforming performance, a large number of robust adaptive beamforming methods have been proposed to eliminate these errors, such as diagonal loading [10], eigenspace-based [11], uncertainty set-based [12], interference-plus-noise covariance matrix (INCM) reconstruction [13][14][15][16] methods. Although these SOP algorithms can effectively suppress multiple interferences, the degrees-of-freedom of interference suppression is limited by the number of antenna array elements.
In order to overcome the shortcomings of the SOP algorithms, the STAP algorithms are proposed to increase the number of interference suppression degrees-of-freedom in GNSS receivers. The STAP algorithms work by placing a finite impulse response (FIR) filter behind each sensor without the need to increase the number of physical array antennas [17][18][19]. The distortionless STAP method can guarantee the distortionless response of GNSS signal and suppress the interference signals [20]. The essence of this method is to ensure the linearity of the phase by constraining the coefficients of the FIR filter. However, the imprecise prior information caused by the direction errors will result in the mismatch of the steering vector, which further deteriorates the performance of STAP method [21]. Therefore, it is of great significance to study a robust STAP method that can deal with steering vector errors. In the vicinity of GNSS receivers, due to the obstruction of buildings or obstacles, GNSS signal will be refracted and reflected, resulting in the coherent interference, which seriously reduces the positioning performance of GNSS receivers [22].
The multipath signal can be used for remote sensing purposes in GNSS interferometric reflectometry, such as the estimation of lake ice thickness in [23]. In general, the GNSS multipath signal is coherent with the direct line-of-sight signal. As a result, the coherent interference causes the spatial spectrum diffusion in the angle region, which reduces the accuracy of the angle estimation algorithms. In order to deal with the coherent interference problems, the spatial smoothing technique [24] resorves the correlation by restoring the rank of the signal covariance matrix. Forward/backward spatial smoothing techniques [25] construct an average forward/backward covariance matrix by dividing the uniform linear array into multiple coincident subarrays, thereby overcoming the defect of rank deficiency. However, this kind of techniques deal with the coherence problems at the cost of the number of effective degrees-of-freedom, and still cannot deal with the coherent interference completely in the low signal-to-noise ratios (SNR) environments. Since the parametric maximum likelihood beamforming technique [26] can handle the coherent interference, it can be considered as an alternative to the methods described above. This method parameterizes all the incident signals, but the derived beamforming optimization problem is a complex high-dimensional expression, which dramatically increases the computational complexity of the system. In [27], a polarization smoothing method is proposed to restore the rank of the matrix by averaging the array data across the polarization domain. Another polarization smoothing method is developed in [28] by performing polarization differential smoothing on the array data, which can work with colored noise. However, both of these methods cannot estimate the polarization parameter. In [29], the 2D direction-of-arrival (DOA) estimation problem of coherent sources is revisited via a polarized uniform rectangular array, where the rank-deficiency of the source matrix is solved by the rearrangement of the array data.
In recent years, various deformation forms of beamformers based on sparse methods [30] have been extensively studied, also known as iterative adaptive approaches (IAA) algorithms, showing good robustness to coherent interference and array perturbations. The IAA algorithms use the weighted minimum variance algorithm to update variables iteratively, and has been widely used in many fields such as multiple-input multiple-output radar imaging, source localization, blood flow velocity estimation, synthetic aperture radar imaging and channel prediction. A generalized IAA [31] is proposed for robust spectral analysis, which utilizes a weighted l p norm to achieve high accuracy and outlier suppression. The IAA-based adaptive beamforming algorithms [32] have robust performance for model errors. The idea of these algorithms is to iteratively estimate the spatial spectrum, which is then used to reconstruct INCM and correct the direction of the desired signal. A virtual uniform linear IAA [33] is proposed for the design of a coprime array robust adaptive beamformer that can directly process snapshots data in the longest virtual uniform linear array. A covariance matrix reconstruction method based on IAA is proposed in [34], which can overcome the direction errors and deal with coherent interference when linear arrays are mutually coupled. So far, there is no literature on the application of IAA algorithm in space-time two-dimensional (ST2D) signal model of GNSS receivers.
In order to deal with the interference suppression of GNSS receivers in coherent signal environments, a robust STAP method based on modified ST2D-IAA spectrum estimation is proposed. First, the IAA algorithm is applied to the ST2D signal model of GNSS receivers, and the ST2D-IAA algorithm is further modified to accurately estimate the power spectrum and the noise power simultaneously. Secondly, the estimated power spectrum and the noise power of the interference angle region are used to reconstruct the space-time interferenceplus-noise covariance matrix (STINCM) of GNSS receivers. Then, the STAP beamforming optimization problem is constructed, and the space-time steering vector (STSV) error vector is obtained after solving this optimization problem, which is used to correct the STSV of GNSS signal. Finally, the weight vector of STAP beamforming is calculated. Simulation results show that the proposed method can suppress interferences in coherent signal environments and is robust to model errors.
Notations:The lower-case (upper-case) bold characters are used to represent the vectors (matrices). C denotes the set of the complex numbers. The superscripts (·) T , (·) H , (·) * and (·) −1 stand for the transpose, Hermitian transpose, conjugate and matrix inversion, respectively. j denotes the imaginary unit with j 2 = −1. ⊗ stands for the Kronecker product. · and | · | denote the Euclidean norm and absolute value, respectively. diag{·} denotes the diagonalization operator. E{·} represents the mathematical expectation. I N denotes the N × N identity matrix. Figure 1 shows a general STAP filter structure for GNSS receiver. Consider a uniform linear array composed of M antenna elements, and the spacing between adjacent elements is d. A temporal filter with N taps and the sampling period T s is connected behind each antenna element. The received signal of each antenna element is down-converted to the baseband by the radio frequency (RF) front end for the analog-to-digital conversion. After passing through the RF front end, the digitized signal is processed by a FIR filter. Finally, the output signal of the STAP filter is sent to each channel of the GNSS receivers for signal acquisition, tracking, and navigation solution. Consider that the direction of GNSS signal reaching the antenna array is θ d , and the directions of G interferences are θ i g (g = 1, 2, · · · , G). The received signal vector of the antenna array can be expressed as

ST2D Signal Model of GNSS Receivers
where d 1 (t), i 1 (t), and n 1 (t) are respectively the GNSS signal, interference and noise, a s (θ d ) and a s (θ i g ) are respectively the steering vectors of the GNSS signal and the g-th interference, s d (t) and s i g (t) are respectively the complex waveforms of GNSS signal and the g-th interference. The steering vector a s (θ) ∈ C M×1 is defined as where λ is the signal wavelength. Stack MN time domain snapshots signals into a vector as where T are the complex envelope sets of the GNSS signal and the g-th interference, respectively.
x mn (t) is the time domain sample at the n-th tap of the temporal filter after the m-th element. Define the weight vector as w = [w 11 , · · · , w M1 , · · · , w 1N , · · · , w MN ] T ∈ C MN×1 , where w mn is the STAP weight at the n-th tap of the temporal filter after the m-th element. The STAP beamformer output can be written as To generate the weight vector, the well-known Power Inversion (PI) criterion is proposed to minimize the output power of the STAP beamformer. In the PI algorithm, an antenna element is selected as the reference element and the remaining antennas are then viewed as auxiliary elements. The weight vector of the PI algorithm is calculated as where

Proposed Robust STAP Method for GNSS Receivers in Coherent Signal Environments
In this section, a robust STAP MVDR beamforming method for coherent signals in GNSS receivers is proposed. We firstly develop a ST2D-IAA algorithm and modify it to obtain the power spectrum of 2D grid points and the noise power. Then we reconstruct the STINCM for GNSS receivers by using the estimated power spectrum of the interference angle region and the noise power. Furthermore, the STSV of GNSS signal is corrected. Finally, the weight vector of the STAP beamforming is calculated based on the reconstructed STINCM and the corrected STSV of GNSS signal.

ST2D-IAA Spectrum Estimation
Considering the STAP beamformer, the ST2D parameters of impinging signals is denoted as where θ represents the pitching angle, f represents the frequency, and f h represents the cutoff frequency of the low-pass filters in GNSS receivers. The whole space domain Θ = {θ|θ ∈ [0, π/2]} can be uniformly discretized asΘ = {θ 1 , θ 2 , · · · , θ K } by K grid points, and the time domain for k = 1, 2, · · · , K and l = 1, 2, · · · , L. The temporal steering vector a t ( f l ) is defined as where f s represents the sampling frequency. The space-time steering matrix of all 2D grid points inΦ is denoted as The complex envelope vector of all 2D grid points inΦ is written as The power vector of all 2D grid points inΦ is given by In order to apply the IAA algorithm to the STAP beamforming of GNSS receivers, we propose a ST2D-IAA algorithm to estimate the power spectrum. Similar to the IAA algorithm, we build the weighted least squares cost function for ST2D-IAA algorithm as x H Q −1 (k, l)x and the matrix Q(k, l) is denoted as with Minimizing (13) with respect to s k,l (t) yields Using (14) and the matrix inversion lemma (see, e.g., [35]), (17) can be rewritten as The steps of the ST2D-IAA algorithm are shown in the Appendix A.1. It should be noted that this algorithm has to be performed in an iterative way because it requires R which contains the unknown signal powers. The initialization of this algorithm is implemented using a standard delay-and-sum beamformer. Let p (i) denote the power vector at the i-th iteration in the ST2D-IAA algorithm. The iteration will be terminated when the relative change p (i) − p (i−1) is less than a specified tolerance. Extensive simulations have shown that the number of iterations is less than 12 for the ST2D-IAA algorithm to converge. One can obtain an accurate power spectrum of all 2D grid points from the results of ST2D-IAA algorithm. In addition, the presence of coherent sources has little influence on the performance of ST2D-IAA algorithm.
Although the ST2D-IAA algorithm can obtain accurate power spectrum, it cannot estimate the noise power iteratively. In order to estimate the power spectrum and the noise power simultaneously in the iteration, we further modify the ST2D-IAA algorithm. Based on the weighted least squares method, we modify the optimization problem for estimating s k,l (t) as below min s k,l (t) where the matrixR is denoted asR = Adiag{p}A H +σ 2 n I MN with the power vector p = [p 1,1 , · · · ,p 1,L , · · · ,p K,1 , · · · ,p K,L ] T and the noise powerσ 2 n . By iteratively solving (19), s k,l (t) in the (i + 1)-th iteration is estimated as whereR where the scaling coefficientφ (i) is obtained from the following optimization problem min ϕ (i) The solution of (26) is given bȳ The criterion for the termination of the iteration process (20)- (27) is given by p (i) −p (i−1) < ς. Through extensive simulations, the modified ST2D-IAA algorithm will converge after at most 12 iterations. The steps of the modified ST2D-IAA algorithm are shown in Appendix A.2.

STINCM Reconstruction
Since the powers of interferences are much higher than those of GNSS signal and noise, the ST2D parameters associated with the interferences can be easily estimated through Pauta criterion. The ST2D parameter of the q-th interference can be denoted as where θ and f denote the spatial and temporal widths occupied by a certain basic 2D subregion, respectively. Φ g corresponds to the area marked in yellow as shown in Figure 2, which contains the following discretized 2D grid points (θ i g ,2 , f i g ,1 ) (θ i g ,2 , f i g ,2 ) · · · (θ i g ,2 , f i g ,V ) . . . . . . . . . . . .
(θ i g ,U , f i g ,1 ) (θ i g ,U , f i g ,2 ) · · · (θ i g ,U , f i g ,V ) Figure 2. Division of the joint space-time domain for the STAP architecture.

STSV Estimation
In order to suppress the interferences with high powers while remaining a distortionless response for the GNSS signal, we formulate the optimization problem of the STAP beamforming based on the MVDR criterion as where R in = E{(i(t) + n(t))(i(t) + n(t)) H } ∈ C MN×MN represents the theoretical STINCM. The constraints in (35) are imposed to guarantee the linearity of the STAP beamforming response so that no biases are introduced into the code and carrier phase measurements [12]. The optimal solution of (35) is given by Since R in is unavailable in practice, we replace it with the reconstructed STINCMR in in (34). Besides, the STSV a is not accurate because it may have STSV error caused by many factors in practical circumstances. By taking the STSV error e into account, the actual STSV can be written as a + e. Based on the reconstructed STINCMR in and the actual STSV a + e, the STAP beamformer output power can be denoted aŝ P(a + e) = 1 (a + e) HR −1 in (a + e) (37) The STSV error e can be estimated by maximizingP(a + e) or, equivalently, by minimizing the denominator ofP(a + e). The STSV error e can be further decomposed into the orthogonal component e ⊥ and the parallel component e . Since the parallel component e does not have influence on the STAP beamforming performance, we dismiss e and only take e ⊥ into computation. The optimization problem for estimating e ⊥ can be formulated as SinceR −1 in 0 with the symbol denoting the positive-definite operation, the convex optimization problem (38) can be easily solved to obtain the solutionê ⊥ . With the estimated e ⊥ , the STSV a can be corrected asâ ). The optimal weight vector is obtained aŝ It should be mentioned that the GNSS signal is generated with the direct-sequence spread-spectrum (DSSS) technique, which shows that it lies in the whole frequency band after passing through low-pass filters in GNSS receivers. In other words, the STSV of the GNSS signal does not have the form a(θ d , f d ) = a s (θ d ) ⊗ a t ( f d ) with a certain frequency f d . Therefore, the constraints of the optimization problem (35) have to be revised for the GNSS signal. The revised optimization problem can be obtained from min w w HR in w subject to N is odd By solving (41), the final weight vector can be given bȳ Based on the weight vector in (42), the array output carrier-to-noise ratio C/N 0 of the GNSS signal can be denoted as where B denotes the signal bandwidth and σ 2 d denotes the power of the GNSS signal. The work [36] assesses the performance of receivers by measuring a zero baseline and a short baseline in the field. In the short baseline test, the multipath error is compensated by a set of harmonic functions, which can guarantee the property of unbiasedness of the least-squares estimators. However, The least-squares method provided in [36] may be less agile than the proposed method of this paper in terms of iterations.

Simulation Results
Consider a uniform linear array consisting of 10 omnidirectional elements with half wavelength spacing and each element is followed by 8-order FIR filter. The desired GNSS signal is assumed as BeiDou-2 (BD2) signal at B3 band whose carrier frequency is 1268.52 MHz and mainlobe bandwidth is 20.46 MHz. The desired BD2 signal impinges on the array from the direction θ d = 5 • . Two interferences whose bandwidths are the same as the desired BD2 signal are respectively incident from the directions θ i 1 = 30 • and θ i 2 = −30 • and the common INR is 30 dB. Only the interference from the direction 30 • is coherent with the desired BD2 signal. In coherent signal environment, the number of samples utilized to calculate the weight vector is chosen as 2400.
Firstly, the spatial-temporal beampatterns corresponding to the proposed method are obtained to observe whether nulls are formed in the directions of interferences. Then, the output C/N 0 performance of the proposed method is mainly compared with the sample matrix inverse method in [37], the diagonal loading in [10], the reconstruct-estimate method in [38]. Although the sampling matrix inverse method is widely used to calculate the weighted vector of STAP beamformers, it can neither suppress the interferences in the coherent signal environments nor be robust to the STSV error caused by the model errors. As for the diagonal loading method, it is a widely used method to deal with space-time covariance matrix errors. The reconstruct-estimate method is also an INCM reconstruction method with STAP structure, which is robust to model errors. In order to verify the robustness of the proposed method in the coherent signal environments and model error conditions, the above methods are selected for comparison. In the reconstructestimate method, the angle region of the desired BD2 signal is [θ d − 5 • , θ d + 5 • ]. In the proposed method of this paper, the two interference angle regions are [θ i 1 − 10 • , θ i 1 + 10 • ] and [θ i 2 − 10 • , θ i 2 + 10 • ], respectively. The threshold ς is set to 10 −3 .
Example 1: Spatial-temporal beampatterns. The 3D and 2D spatial-temporal beampatterns of the proposed method are shown in Figure 3a and Figure 3b respectively. It can be seen that while the entire frequency band at the desired BD2 signal direction forms a mainlobe, nullings are formed in the entire frequency bands of the interference directions 30 • and −30 • , respectively. This shows that the proposed STAP beamforming method can effectively suppress the interferences in the coherent signal environments, because the proposed method uses the modified ST2D-IAA method to estimate the precise power spectrum.  Example 2: Output C/N 0 performance. Figure 4a shows the relationship curves of the sampling matrix inverse method, diagonal loading method, reconstruct-estimate method and the proposed method with the change of the output C/N 0 versus the input SNR. Note that although the GNSS signal input SNR is usually −20 dB, the spot beam has been or will be equipped with modern GNSS in the future, so the corresponding SNR will be increased to 27 dB, so the SNR range of BD2 signal is [−30, 10] dB. It can be seen from Figure 4a that the output C/N 0 performance of the proposed method is better than those of the other methods, because the proposed method can suppress the interferences in the coherent signal environments, while the other methods cannot deal with the interferences in the coherent signal environments. Figure 4b shows the relationship curves of the sampling matrix inverse method, diagonal loading method, reconstruct-estimate method and the proposed method with the change of the output C/N 0 versus the number of snapshots. It can be seen that the output C/N 0 performance of the proposed method is still better than the other methods, and the proposed method has good advantages in the output C/N 0 .  Example 3: The scenary of direction error. When the direction error of the desired BD2 signal is 3 • , that is, the actual direction of the BD2 signal is 8 • . In this case, Figure  5a and Figure 5b respectively show the relationship curves of the change of the output C/N 0 versus the input SNR and the number of snapshots, while Figure 5c shows the relationship curves of the output C/N 0 versus the direction error. It can be seen from these figures that both the proposed method and the reconstruct-estimate method are robust to the signal direction error, because both of them consider the signal STSV error and correct the signal STSV. However, the performance of the diagonal loading method and the sampling matrix inverse method is poor in the case of signal direction error, because these two methods do not consider the signal STSV error, and do not correct the signal STSV. In addition, the proposed method has the best performance of the output C/N 0 , which is because the proposed method not only reconstructs the STINCM, but also corrects the signal STSV. Example 4: The case of arbitrary STSV error. In the case of arbitrary STSV error of the signal, the actual STSV of the desired BD2 signal can be expressed as a = a(θ d + δθ) + δa, where the direction error δθ is fixed at 3 • , the STSV error δa is a random vector with the mean 0 and the variance σ 2 I MN×1 . The parameter σ 2 is set to 0.64. Figure 6a and Figure 6b respectively show the relationship curves of output C/N 0 versus the input SNR and versus the number of snapshots under the condition of arbitrary STSV error of the signal. It can be seen that the output C/N 0 of the proposed method is higher than that of the other methods, because the proposed method combines the accurate STINCM and the corrected STSV, which indicates that the proposed method is robust to arbitrary STSV error of signal. Example 5: The situation of small number of snapshots. We investigate the performance of the tested methods when the number of snapshots is small. The simulated curves of output C/N 0 versus the number of snapshots are displayed in Figure 7, where the number of snapshots ranges from 20 to 200. Obviously, the proposed method has the minimal performance loss, and performs better than the other methods in the situation of small number of snapshots. The reason lies in that the ST2D-IAA spectrum estimation in the proposed method can work with the limited sampling snapshots. As described in [33], the IAA spectrum estimation can even handle the single snapshot in the coprime virtual uniform linear array. Therefore, the performance of the proposed method is robust against the small number of snapshots.

Conclusions
In this paper, a robust STAP method based on the modified ST2D-IAA spectrum estimation is proposed to suppress interferences for GNSS receivers in coherent signal environments. By applying the IAA algorithm to the ST2D signal model in GNSS receivers and further modifying the ST2D-IAA algorithm, the power spectrum and noise power can be estimated accurately at the same time. The STINCM of GNSS receivers is reconstructed by selecting the estimated power spectrum and noise power corresponding to the interference angle regions. In order to obtain the error vector of the STSV, the STAP beamforming optimization problem based on the MVDR criterion is constructed to correct the STSV of GNSS signal. The weight vector of STAP beamforming is calculated by combining the reconstructed STINCM and the corrected GNSS signal STSV. Simulation results show that the proposed method can suppress the interferences in the coherent signal environments, and has better robustness than the existing methods under the conditions of model errors. The proposed method has a potential prospect in maintaining the integrity and functionality of GNSS systems against coherent interference in many applications.  The ST2D-IAA algorithm is summarized in Algorithm A1.