Abstract
The well-known Space-Alternating Generalized Expectation Maximisation (SAGE) algorithm has been recently considered for multipath mitigation in Global Navigation Satellite System (GNSS) receivers. However, the implementation of SAGE in a GNSS receiver is a challenging issue due to the numerous number or parameters to be estimated and the important size of the data to be processed. A new implementation of the SAGE algorithm is proposed in this paper in order to reach the same efficiency with a reduced complexity. This paper focuses on the trade-off between complexity and performance thanks to the Cramer Rao bound derivation. Moreover, this paper shows how the proposed algorithm can be integrated with a classical GNSS tracking loop. This solution is thus a very promising approach for multipath mitigation.
1. Introduction
In Global Navigation Satellite System (GNSS) applications, multipath (MP) errors are still one of the major error sources for conventional receivers. The additional signal replicas due to reflections on the local environment introduce a bias in the delay lock loops (DLLs), which finally leads to a positioning error [1, 2]. Several techniques have been developed for multipath mitigation. One of the most popular approaches is the Narrow Correlator Spacing [3], which reduces the chip spacing between the early and late correlators in order to mitigate the impact of multipath. However, this technique suffers from high sensitivity to noise and cannot perform with short delay multipath (<0.1 chip). Based on the Maximum Likelihood (ML) estimation, the Multipath Estimating Delay-Lock-Loop (MEDLL) [4] algorithm has also been proposed to estimate the delay and the power of all the paths by studying the shape of the cross-correlation function. This approach shows better performances than the Narrow Correlator Spacing technique, but short delay multipath mitigation is still an issue [4]. More recently, Bayesian approaches have been proposed [5–7]. Indeed, most of the time, prior information could be used in order to improve the delays estimations. However in practice, it is difficult to get correct prior information. Measurement campaigns can be used to build a first-order Markov process for a sequential estimation, but the performance will consequently be strongly dependent on the measured environment (design of the city…).
Last, the use of array antenna algorithms has been proposed for multipath mitigation [8, 9]. Array antennae enable a spatial sampling that makes it possible to distinguish different sources in the spatial domain. Therefore, mitigation techniques based on such an array are independent of the relative delays of the MP. Consequently, the rejection of any kind of MP seems possible (and especially short-time delay MP). In conventional mobile receivers, the room available to incorporate an array antenna is reduced. Only a small number of elements can be integrated. This study will focus on a square elements antenna. The choice of the shape of the array antenna is out of the scope of this paper, and we will only focus on the array processing.
Using appropriate weights on each channel, beamforming technics perform a spatial filtering in order to concentrate the energy beam towards the signal of interest, while trying to minimise the gain towards the interferences. However, in the MP mitigation context, these classical approaches present limited performances for the following reasons [10].(i)Near angle signals are difficult to separate because of the low antenna resolution [10, 11].(ii)Line-of-sight signal (LOSS) and multipaths are strongly correlated, which implies a severe degradation of adaptive beamforming algorithms.(iii)Adaptive algorithms should operate after the correlation step in order to work with positive SNR signals. That implies that only a small number of samples are available to estimate the array covariance matrix. Bad covariance matrix estimation can lead to inversion or eigenvalues decomposition instabilities.
To overcome these problems, several solutions have been proposed. In [12, 13], for example, the authors propose to use additional Choke ring techniques, GPS microstrip array antenna, and/or angle constraint for negative elevation in order to reject ground multipaths. However, for a generic application, we cannot assume any specific MP DOA.
In order to improve the MP mitigation, another approach proposes to include the different DOA paths parameters in the estimation procedure [10, 14–18]. In other words, we estimate a set of parameters (amplitudes, times/delays, Doppler shifts, elevations, and azimuths) for all the incoming sources. The main difference with the beamforming approach is that, instead of filtering the sources in the spatial domain only, the different incoming paths are jointly identified in the space, time, and frequency domains. In order to estimate the parameters of all impinging signals, the Space Alternating Generalized Expectation Maximisation (SAGE) algorithm [16], which is a low-complexity generalization of the Maximum Likelihood (ML) algorithm, has been considered. The SAGE algorithm is usually used in communication systems [16], but the potential of SAGE in a navigation context has also been proven [17, 18]. Nevertheless, the computational cost increases due to the number of unknown parameters. Moreover, the memory size is also a challenging issue. Thus, the SAGE algorithm can hardly be directly implemented in real time for GNSS receivers. Last, the SAGE algorithm provides an estimation of the time delay and phase of the LOSS, as do the DLL and PLL (Phase Lock Loop). Thus, the use of SAGE implies to switch off the DLL/PLL, and, consequently, we lose the “smoothing effect” of the tracking loops.
In order to reach the efficiency of the SAGE algorithm with a reduced complexity, and to keep the compatibility with conventional GNSS tracking loops (DLL/PLL), we proposed a new implementation of SAGE in [10]. The main idea is to apply the SAGE algorithm after the local correlation step. Indeed, the cross-correlation between the received signal and the local code can be seen as a compression step of the baseband signal. Thus, by reducing the size of the input signal, the complexity of the algorithm will reduce in the same proportion. By applying the ML estimation on the postcorrelated signal, we estimate the relative delay and Doppler of the paths. These estimations are then used to drive the GNSS tracking loops. We named this new algorithm the SAGE/STAP multicorrelator algorithm.
Comparisons based on Monte Carlo simulations between the SAGE/STAP algorithm, the SAGE algorithm, and beamformers have been done in [10]. This paper focuses on the Cramer Rao Bound (CRB) derivation in order to present the theoretical trade-off between the complexity reduction and the estimation accuracy of the SAGE/STAP algorithm. Moreover, this paper presents an implementation of the proposed algorithm hybridized with the GNSS tracking loops
This paper is organized as follows. The signal model and the main assumptions are outlined in Section 2. A review of the SAGE algorithm, the SAGE/STAP multicorrelator algorithm, and its implementation in a GNSS receiver is given in Section 3. The simulation results in Section 4 show the trade-off between the complexity and the accuracy of the SAGE/STAP multicorrelator algorithm and present the performances of the proposed algorithm in realistic scenarios. Finally, in Section 5 we present our conclusions.
2. Signal Model
Let us assume that we receive narrowband planar wave fronts of wavelength on an array of isotropic sensors. The M-sized vectorized received signal, , can be written as where is given by and is an additional complex white Gaussian noise. Let us note that the number of incoming paths (LOSS included) and are, respectively, the complex amplitude, elevation, azimuth, time delay, and Doppler shift of the th path. Note that the index corresponds to the LOSS. Here, denotes the pseudorandom-noise (PRN) sequence that consists of a Gold code with a code period ms, 1023 chips per period, and a rectangular chip shape. represents the steering vector of a square array antenna. The channel parameters are assumed constant during the observation time.
Collecting the samples of the complex baseband signal leads to a , matrix, where denotes the number of SnapShots. Typically in GNSS, is larger than 2000, and, thus, a direct processing of the baseband signal is hardly implemented in real time. In order to compress the signal, we propose to use the scheme presented in the Figure 1. The aim of the array antenna is to sample the wave fronts and to get access to the spatial domain. The banc of correlators is used to compress the signal and to get access to the relative delays of the sources, and finally the postcorrelation time taps enable to get access to the frequency domain. We can see that the first and last parts of this architecture are equivalent to a STAP (space time adaptive process) array applied to the postcorrelated signal. Thus, we call this architecture the STAP multicorrelator array. This array is then defined by 4 parameters: the number of antennas, the number of postcorrelation taps, the number of correlators, and the time spacing between each correlator. Note that the sampling time of the postcorrelated signal is equal to . Therefore, we have the following relation between the number of SnapShots and the number of postcorrelated Taps: , where denotes the sampling frequency of the baseband signal. The th output of one correlator delayed by a time spacing of is given by with , if is even and if is odd . We introduce the relative delay of the DLL estimation and the relative Doppler of the FLL (Frequency Lock Loop) estimation with respect to the th path parameters: . We can then approximate the integral as with the autocorrelation function of the PRN code and . The amplitude and phase components are independent of and ; they can be inserted in the modified complex amplitude of the incoming path: Then, with
The outputs of the correlators are rearranged in a column vector. The cross-correlation functions obtained are then concatenated in order to get the temporal evolution of the postcorrelation signal:
Finally, the output signal of the STAP multicorrelator array can be collected in a column vector: where is the output noise and is given by
where denotes the Kronecker product, and the parameters are now , where denote the relative delay and Doppler of the path .
3. SAGE Algorithm
3.1. Concept of SAGE
The problem associated with the signal model (2) consists in estimating the parameters for all the paths. The estimation of is not discussed in this paper. Usually, is fixed to a value large enough to capture all the dominant impinging waves. Classical information theory methods for model size selection like Akaike’s and Rissanen’s [19] criteria can be used. The ML estimation is given by , where is the sampled complex baseband signal and is the likelihood function. The direct maximization of the likelihood function is a computationally prohibitive task since there is no analytical solution. Moreover, is generally not a concave function of , and is usually high. To perform this optimization, we use the iterative process of the SAGE algorithm [16]. The basic concept of the SAGE algorithm is to us a hidden data space. Instead of estimating the parameters of all impinging waves in one search, the SAGE algorithm sequentially estimates the parameters of each signal. The SAGE algorithm breaks down the multidimensional optimization problem into several smaller problems. In spite of this complexity reduction, SAGE is still hardly implemented in real time due to the size of the baseband signal.
3.2. The SAGE/STAP Multicorrelator Algorithm
The main idea to take advantages of SAGE with more reasonable hardware requirements is to process the data at the output of the STAP multicorrelator array. Using the signal model (10), the likelihood function is then where contains the superimposition of the post-correlation signals. denotes the covariance matrix of the postcorrelation noise. The noise is no longer white after the correlation step, and it is shown in [20] that the postcorrelation noise covariance matrix is given by where denotes the noise power after the correlation step, and is define thanks to the auto-correlation function of the PRN code :
The first step of the SAGE algorithm, so-called expectation step (E-STEP), consists in estimating the hidden data space with
The second step, so-called maximization step (M-STEP), carries out the maximization of the log-likelihood function which is associated with the estimated hidden data space. In the case of the STAP multicorrelator signal, the log-likelihood function is
The maximization of Λ with respect to can be concentrated, as the dependence to is linear: where is given by
Finally, the reduced likelihood function to maximize is
In order to better understand the expression (18), let us assume no noise is present, and we can develop (10) as where
Then, inserting (19), (20), and (21) in (18), we can express the reduced likelihood function as
As we can see, the first two terms are similar to space and frequency Fourier transforms. Thus, the two parameters and will mainly influence the accuracy of the DOA and the Doppler estimations for both algorithms (SAGE and SAGE/STAP). However, the couple influences only the performances of the SAGE/STAP algorithm with respect to the paths delay estimation. This point will be detailed in more detail in Section 4.
3.3. Implementation in a GNSS Receiver
To implement the SAGE/STAP algorithm in a GNSS receiver, we propose the architecture illustrated in Figure 2.
The first front end bloc contains the RF filters, the downconversion step, the ADC and a calibration stage [21]. A calibration stage aims to compensate the technological defects such as the RF dispersion of the filters and antenna defects.
The baseband signal is the output of this blocx. We must therefore use the STAP multicorrelator array for each tracked PRN code in order to compress the data. Afterward, we use the SAGE algorithm in the processing bloc. Finally, the SAGE estimation of the relative delay and Doppler is used to drive the DLL and PLL/FLL. Thus, the key idea of our solution is to substitute the conventional DLL and PLL/FLL discriminators by the SAGE/STAP estimation.
4. Simulation Results
The aim of this section is to study the impact of the compression step (realized by the STAP multicorrelator array) on the SAGE algorithm performances. Thus, this section focuses on the comparison between the classical SAGE and the SAGE/STAP multicorrelator algorithms. The estimation performances of both algorithms were numerically evaluated with Monte Carlo simulations and theoretically through the CRB derivation. The derivation of the CRB for the SAGE/STAP algorithm is given in the annexe. For the conventional SAGE algorithm, the CRB can be found in [17]. The simulations parameters are as follows: dBW, dBW, , , Hz, and we use blocks of 1 ms of integration. The sampling frequency of the baseband signal is MHz. The square array antenna contains 4 isotropic identical sensors spaced by .
In a first time, we analyze the influence of (the number of correlators) and (the correlators’ delay spacing) on the performances of the SAGE/STAP algorithm. In order to cover all the cross-correlation function, we should have . This condition is necessary if we want to take into account all the multipaths with a relative delay . Therefore, we have the following relationship , where the operator denotes the nearest integer greater than or equal to . We now focus our attention on the parameters.
To study the influence of , we calculate the root-mean-square error (RMSE) of the time delay based on 100 Monte Carlo simulations for both algorithms (SAGE and SAGE/STAP multicorrelator). In this first simulation, we assume no multipath, and the RMSE is plotted on Figure 3 as a function of , for different signal bandwidths: MHz and MHz, where denotes the two-sided RF bandwidth. When comparing the RMSE of SAGE and SAGE/STAP algorithms, we can see that both algorithms present the same performances in the case . This behavior is confirmed theoretically with the curve of the CRB as a function of . We can see that the standard deviation increases when becomes higher than the inverse of the signal bandwidth. We can deduce from this observation the following empirical rule which results from a tradeoff between accuracy and complexity:
This conclusion can also be observed on the shape of the likelihood function. In Figures 4 and 5, we plot a section of the likelihood function through the delay for the SAGE and SAGE/STAP multicorrelator algorithms. In other words, we compare and for different signal bandwidths ( MHz on Figure 4 and MHz on Figure 5) and different correlator spacings . We can see that if the condition (23) is fulfilled, the likelihood functions are identical for the SAGE and SAGE/STAP multicorrelator algorithms. If not, the likelihood function of the SAGE/STAP algorithm becomes wider than the classical SAGE likelihood function, which will consequently increase the variance of the estimation.
In a second time, we evaluate the performances of both algorithms in the case of one multipath. The reflected multipath and the LOSS are considered to be in-phase which corresponds to one of the worst possible cases. In Figure 6, we plot the RMSE of the estimated time delay for both algorithms as a function of the relative time delay and relative azimuth between the LOSS and the MP (denoted and ). The multipath parameters are dBW, , Hz.
We assume that the receiver is perfectly synchronized with the LOSS at the beginning of the simulation. Therefore the initial relative parameters are Hz, Hz, Chip. In order to choose the couple , we use the results of the first simulation together with the condition (23): MHz, Chip, and . Here again, we can see that both algorithms present the same performances whatever the multipath position (in space and delay). We also plot in Figure 7 the CRB for the SAGE and SAGE/STAP algorithms. The behaviour of both CRB shows that the algorithms are theoretically equivalent. Thus, the SAGE/STAP multicorrelator can reach the efficiency of the SAGE algorithm but with a strongly reduced complexity. Indeed, the size of the baseband signal is compared to for the postcorrelated signal in the optimal configuration. The compression factor is .
We compare now the RMSE and the CRB of the SAGE/STAP algorithm in Figure 8 in the same configuration. We can see that the RMSE reaches the CRB if the sources are spaced by more than 50°. This result is particularly interesting because in this situation, close and far time delay multipaths are completely mitigated by the algorithm, which is mainly due to the multiantenna contribution. In the case where the relative azimuth is weaker, the algorithms mainly use the time and frequency domains to mitigate the multipath. In this condition, the array antenna is not useful anymore, and the rejection performances are similar to single antenna mitigation techniques (e.g., MEDLL). In Figure 9, we plot the RMSE and the CRB for the LOSS elevation estimation. Here again, the RMSE reaches the CRB, excepted in the case of short-time delay ( Chip) and close space sources (). This shows the limitation of the SAGE approach for cases where the MP is close to the LOSS on the 4 dimensions considered in the estimation process.
The efficiency of the SAGE/STAP algorithm was illustrated in the previous simulations. Now, we want to address the impact of this algorithm on the performances of conventional GNSS tracking loops.
First of all, we evaluate the error envelope for a conventional DLL enhanced by a conventional beamformer (CBF) and for a DLL driven by SAGE (denoted DLL/SAGE). The loops parameters are loop order = 2, damping factor = 0.707, Hz, Hz. In the case of classical DLL, we use the early-minus-late discriminator (EMLD) with a chips spacing equal to 0.1 Chip and the arctangent discriminator for the PLL. In Figure 10, we plot the error envelope as a function of the relative azimuth and delay of the MP. Although the DLL is enhanced by a CBF, the SAGE/STAP algorithm provides the best multipath mitigation whatever the MP position.
In the last simulation, we propose to test the SAGE/STAP algorithm with the DLR channel model available online (http://www.kn-s.dlr.de/satnav/) [22]. The DLR channel model generates the number of paths (LOSS + MP), their complex amplitudes, relative delays, and Doppler shifts. In the free version available online, the DOA are not generated. Thus, we computed the DOA of the MP thanks to the position of the scatterers, which are randomly positioned based on statistics established during measurement campaigns [22].
The LOSS parameters are , , Hz, and we simulated a CBOC signal [23] on the L1 band with MHz. Thus, according to (23), we have and ns. The speed of the vehicle has been chosen to 0 m/s in order to tackle the problem of stationary MP. The other simulation parameters (, DLL and PLL) are the same as in the previous simulations. In Figure 11, we plot the pseudorange error and the Doppler estimation of the conventional DLL/PLL and the DLL/PLL driven by the SAGE/STAP multicorrelators algorithm. As we can see, the SAGE approach provides a real improvement in the time delay and Doppler estimation with respect to the conventional multiantenna receiver architecture in a tracking scenario. More tracking simulations with this solution can be found in [24].
(a)
(b)
5. Conclusion
In this work we have addressed the problem of estimating the propagation time delay of the LOSS in a GNSS receiver under severe multipath conditions. To reduce the influence of the multipaths, we investigated the use of array antenna algorithms. Previous studies suggest using the SAGE algorithm with an array antenna in order to reduce the estimation error of the delay of the LOSS. However, SAGE is hardly implemented in real time due to the memory requirements and computation cost. Moreover, SAGE is hardly compatible with classical GNSS tracking loops.
In order to take the advantages of SAGE with more reasonable hardware requirements, we proposed a new implementation based on a STAP multicorrelator array. The CRB derivation and the Monte Carlo simulations were used in this paper to study the trade-off between accuracy and complexity. This trade-off appears by the influence of the parameters (number of correlators) and (correlator space) on the SAGE/STAP multicorrelator performances. The size of the data to be processed is then drastically reduced, and the CRB and Monte Carlo simulations attest that there is no performance loss between the classical SAGE algorithm and the SAGE/STAP multicorrelator algorithm. Moreover, the SAGE/STAP multicorrelator algorithm can constitute the discriminator of GNSS tracking loops. Thus, with this new implementation, the high-resolution performances of SAGE are now available for GNSS receiver, which is a very promising approach for the multipath mitigation problem in GNSS.
Appendix
Cramer Rao Bound
The likelihood function for the signal at the output of the STAP multicorrelator array is given in (11). Therefore, the log-likelihood function is where contains the superimposition of the postcorrelated signals, denotes the covariance matrix of the postcorrelated noise, and is the parameters vector with the modified complex amplitude vector of the impinging wave fronts, the elevation vector, the azimuth vector, the relative Doppler vector, and the vector of relative delay.
The CRB is found as the element of the inverse of the so-called Fisher Information Matrix (FIM) :
The definition of the FIM is
As we assume that the noise covariance is independent of the parameters, the FIM for a complex multivariate Gaussian process is [6, 11]
Note that the FIM can be a bloc partitioned symmetric matrix. Thus, we can write where the bloc matrices can be obtained by with . Last, we need to calculate the differential of the STAP multicorrelator signal model given in (10) with respect to . To that end, we use the following signal model:
The differential of the Kronecker product of 2 matrixes and with respect to a parameter is given by [6] with the permutation matrix. Note that in our model, we are working with column vectors, and each subvector is depending on different independent parameters. Therefore, the differential of the signal is simply given by
Acknowledgments
This work was supported by the French Space Agency (CNES), French Aerospace Lab (ONERA), and Thales Alenia Space—France.