Advances in Signal Tracking for GNSS Receivers: Theory and Implementation
View this Special IssueResearch Article  Open Access
Nonlinear Bayesian Tracking Loops for Multipath Mitigation
Abstract
This paper studies Bayesian filtering techniques applied to the design of advanced delay tracking loops in GNSS receivers with multipath mitigation capabilities. The analysis includes tradeoff among realistic propagation channel models and the use of a realistic simulation framework. After establishing the mathematical framework for the design and analysis of tracking loops in the context of GNSS receivers, we propose a filtering technique that implements RaoBlackwellization of linear states and a particle filter for the nonlinear partition and compare it to traditional delay lock loop/phase lock loopbased schemes.
1. Introduction
Global Navigation Satellite Systems (GNSS) are the general concept used to identify those systems that allow user positioning based on a constellation of satellites. Specific GNSS are the wellknown American GPS, the Russian GLONASS, or the forthcoming European Galileo. All those systems rely on the same principle: the user computes its position by means of measured distances between the receiver and the set of inview satellites. These distances are calculated estimating the propagation time that synchronously transmitted signals take from each satellite to the receiver. Therefore, GNSS receivers are only interested in estimating the delays of signals which are received directly from the satellites, referred to as lineofsight signal (LOSS), since they are the ones that carry information of direct propagation time. Hence, reflections distort the received signal in a way that may cause a bias in delay and carrierphase estimations. Multipath is probably the dominant source of error in highprecision applications, especially in urban scenarios, since it can introduce a bias up to a hundred of meters when employing a 1chip wide (standard) delay lock loop (DLL) to track the delay of the LOSS, which is a common synchronization method used in spreadspectrum receivers. This error might be unacceptable in many applications.
Sophisticated synchronization techniques estimate not only LOSS parameters but those of multipath echoes. This results in enhanced, virtually biasfree pseudorange measurements. In this paper, we investigate multipath estimating tracking loops in realistic scenarios, where this effect is known to be severe. The analysis is driven in two directions. Firstly, a review of statistical characterization of the channel model in such situations is performed and a commercial signal simulator. Secondly, a novel multipath estimating tracking loop is discussed, providing details on the implementation, as well as comparisons to stateoftheart techniques when different channel characteristics are considered. This tracking loop resorts to the Bayesian nonlinear filtering framework, sequentially estimating the unknown states of the system (i.e., parameters of the LOSS and echoes) and providing robust pseudorange estimates, subsequently used in the positioning solution. The socalled multipath estimating particle filter (MEPF) considers RaoBlackwellization of signal amplitudes and the use of a suitable nonlinear filter for the rest of nonlinear states, for example, timedelays and their rate. More precisely, RaoBlackwellization involves marginalization of linear states and the use of a standard Kalman filter to track signal amplitudes with the goal of reducing the estimation variance, since (i) the dimensionality of the problem that nonlinear filters solve is reduced and (ii) linear states are optimally tackled. For the nonlinear part of the state space we consider sequential MonteCarlo methods (specifically, the standard particle filtering) as one of the most promising alternatives in advanced GNSS receiver designs. Realistic computer simulation results are presented using the GRANADA FCM signal simulator and the performance of the MEPF is evaluated.
The remainder of the paper is organized as follows. Section 2.1 provides a brief overview of the fundamentals of GNSS, their signal structure, available channel models, and receivers’ architecture and describes a realistic simulation platform. Section 3 sketches the basics of particle filters, and Section 4 is devoted to their application to GNSS signal synchronization in the presence of multipath. Section 5 presents computer simulations, and finally Section 6 concludes the paper. For the sake of completeness, the paper shows in the Appendix the equivalence between precorrelation and postcorrelation processing of GNSS signals. Notice that in this paper, the MEPF method operates after correlation is performed in order to operate at a lower data rate.
2. Fundamentals of Global Navigation Satellite Systems
GNSS space vehicles broadcast a lowrate navigation message that modulates continuous repetitions of pseudorandom spreading codes, that in turn are modulating a carrier signal allocated in the L band. The navigation message, after proper demodulation, contains among other information the socalled ephemeris, a set of parameters that allow the computation of the satellite position at any time. These positions, along with the corresponding distance estimations, allow the receiver to compute its own position and time, as we will see hereafter. Basically, a GNSS receiver performs trilateration, a method for determining the intersections of three or more sphere surfaces given the centers and radii of the spheres. In this case, the centers of the spheres are the satellites, whose position can be computed from the navigation message, and the radii of the spheres are the distances between the satellites and the receiver, estimated from the time of flight.
The distance between the receiver and a given satellite can be computed by where m/s is the speed of light, is the receiving time in the receiver’s clock, and the time of transmission for a given satellite . Receiver clocks are inexpensive and not perfectly in sync with the satellite clock, and thus this time deviation is another variable to be estimated. The clocks on all of the satellites belonging to the same system , where , are in sync with each other, so the receiver’s clock will be out of sync with all satellites belonging to the same constellation by the same amount . In GNSS, the term pseudorange is used to identify a range affected by a bias, directly related to the bias between the receiver and satellite clocks. There are other factors of error: since propagation at speed is only possible in the vacuum, atmospheric status affects the propagation speed of electromagnetic waves modifying the propagation time and thus the distance estimation. For instance, the ionosphere, that is the part of the atmosphere above km until km of the Earth surface, is a plasmatic medium that causes a slowdown in the group velocity and a speed up of the phase velocity, having an impact in code and phase delays and, thus, impeding precise navigation when its effects are not mitigated. Actually, errors can be on the order of tens of meters in geomagnetic storm episodes [1].
For each inview satellite of system , we can write where is the satellite’s position (known from the navigation message), the receiver’s position, and gathers other sources of error. Since the receiver needs to estimate its own 3D position (three spatial unknowns) and its clock deviation with respect to the satellites’ time basis, at least satellites must be seen by the receiver at the same time, where is the number of different navigation systems available (inview) at a given time. Each received satellite signal, once synchronized and demodulated at the receiver, defines one equation such as the one defined in (2), forming a set of nonlinear equations that can be solved algebraically by means of the Bancroft algorithm [2] or numerically, resorting to multidimensional NewtonRaphson and weighted least square methods [3]. When a priori information is added we resort to Bayesian estimation, a problem that can be solved recursively by a Kalman filter or any of its variants. The problem can be further expanded by adding other unknowns (for instance, parameters of ionospheric and tropospheric models), sources of information from other systems, mapping information, and even motion models of the receiver. In the design of multiconstellation GNSS receivers, the vector of unknowns can also include the receiver clock offset with respect to each system in order to take advantage of a higher number of inview satellites and using them jointly in the navigation solution, therefore increasing accuracy.
2.1. Signal Model
A general signal model for most navigation systems consists of a directsequence spreadspectrum (DSSS) signal [4], synchronously transmitted by all the satellites in the constellation. This type of signals enables code division multiple access (CDMA) transmissions, that is, satellite signals are distinguished by orthogonal (or quasiorthogonal) codes. At a glance, these signals consists of two main components: a ranging code (the PRN spreading sequence) and a low rate data link (broadcasting necessary information for positioning such as satellites orbital parameters and corrections). The complex baseband model of the signal transmitted by a GNSS space vehicle reads as where being the transmitting power, a parameter controlling the power balance, the data symbols, the bit period, the number of repetitions of a full codeword that spans a bit period, the codeword period, a chip of a spreading codeword of length chips, the transmitting chip pulse shape, which is considered energy normalized for notation clarity, and is the chip period. Figure 1 aims at clarifying the relation between those bits/chips parameters. Subindex refers to the inphase component, and all parameters are equivalently defined for the quadrature component, referred to with the subindex . This signal model describes all GNSS’s signalsinspace, for instance GPS L1, GPS L5, Galileo E1, and Galileo E5. Refer to [5] for the details.
2.2. Propagation Channel Model
A key aspect in the definition of the propagation channel model between satellites’ antenna and the user’s receiver antenna is whether it can be considered narrowband or wideband, which depends on the bandwidth of the propagation channel in which a given signal is transmitted, being assessed with respect to the channel coherence bandwidth. The coherence bandwidth is defined as the frequency band within which all frequency components are equally affected by fading due to multipath. In narrowband systems, all the components of the signal are equally influenced by multipath, while in wideband systems the various frequency components of the signal are differently affected by fading. Narrowband systems, therefore, are affected by nonselective fading, whereas wideband systems are affected by selective fading. The coherence bandwidth depends on the environment and is given by where is the delay spread, which is the time span between the arrival of the first and the last multipath signals that can be sensed by the receiver. In a fading environment, a propagated signal arrives at the receiver through multiple paths. For a typical GNSS multipath propagation channel in which s (the limit can be greater in nonurban areas, but in general it is not lower), we obtain that the system is wideband if transmitted signals are wider than kHz, which is the case for GNSS waveforms (in the order of MHz). Hence, we conclude that we need to define propagation channel models considering wideband systems. Another important definition within this context concerns coherence time. The coherence time, , is defined as the time interval during which the characteristics of the propagation channel remain approximately constant, and it is given as where is the maximum Doppler shift. The Doppler shift is given as , where is the radial speed of the mobile terminal with respect to the satellite and is the signal wavelength. A channel is considered WSSUS (widesense stationary with uncorrelated scatterers) during the coherence time.
In the following, we describe four of the most relevant satellite channel models found in the literature.
2.2.1. Jahn’s Channel Characterization
Jahn et al. provided a wideband channel model for land mobile satellite services [6]. The model was derived from a channel measurements campaign performed in the L band at 1820 MHz. An aircraft transmitted a spread spectrum signal of 30 MHz, being received by a mobile receiver (handheld or car terminal). From those measurements, authors characterized the channel assuming WSSUS and modeling it as a filter structure with delay taps. Then, they provided statistical models for LOS (Rician probability density function for the amplitude of the direct path), shadowing (ray amplitude following a Raileigh distribution with a lognormal distributed mean power), near echoes (the number of the near echoes follows a Poisson distribution, with delays being exponentially distributed and amplitudes following a Rayleigh distribution), and far echoes (same distributions than near echoes but with other parameters). Table 1 summarizes the main features of Jahn’s statistical channel model.

2.2.2. Loo’s Channel Characterization
The Loo’s land mobile satellite channel model [7] is a statistical model that assumes that the LOS component under foliage attenuation (shadowing) is lognormally distributed and that the multipath effect is Rayleigh distributed. This model provides complete statistical descriptions for different shadowing and multipath conditions based on an extensive measurement campaign for different frequency bands. For the L band, the “Inmarsat’s Marecs A” satellite was used as transmitter, while a mobile laboratory was considered for signal reception, resulting in a fixed elevation. Many more investigations on Lband measurements are also referred to in [8], obtaining results for other elevation angles. Table 1 summarizes the main features of Loo’s statistical channel model.
2.2.3. PérezFontán’s Channel Characterization
The model presented by Fontán et al. in [9] addressed the statistical modeling of shadowing and multipath effects in land mobile satellite applications for a wide range of environments with different clutter densities (from open to dense urban areas) and elevation angles (from 5° to 90°) at L, S, or Ka Bands, using a comprehensive experimental database to extract the model parameters for the different bands, environments, and elevations. One of its main contributions consists of producing time series of any channel parameter whose study is required, instead of just cumulative distribution functions. These ones may be computed later from the generated series. The model uses a firstorder Markov chain to describe the slow variations of the direct signal, basically due to shadowing/blockage effects. The overall signal variations due to shadowing and multipath effects within each individual Markov state are assumed to follow a Loo distribution with different parameters for each shadowing condition (Markov state). Up to this point the model is of the narrowband type since it does not account for time dispersion effects. These effects are introduced by using an exponential distribution to represent the excess delays of the different echoes. Table 1 summarizes the main features of PérezFontán’s channel model.
2.2.4. Steingass/Lehner’s Channel Characterization
The Steingass/Lehner land mobile channel model presented in [10] was developed using data recorded in a highresolution measurement campaign carried out in Munich in 2002. Different types of environments (urban, suburban, and rural) were measured for car and pedestrian applications. It has been approved as standard by the ITU [11]. For the measurements, a 100 MHz signal near the GPS L1 band was used. This signal provided a time resolution of about 10 ns. The received signal was processed using a superresolution algorithm to extract the single reflections. With this information, the probability density distribution of the parameters of the reflected rays, such as Doppler shift, power of echoes, duration of a reflector, and number of echoes, were extracted. In urban environments, three major obstacles influence the propagation of the LOS signal: house fronts, trees, and lamp posts. The model is comprised of a deterministic part with a generated scenery, which computes geometrically the LOS signal shadowing and knifeedge diffraction for house fronts, lamp posts, and trees. The other observables like the number of coexisting echoes, life span of reflectors, and the mean power of the echoes are generated stochastically, using the probability density distribution extracted from the measurements. The output of the model is a complex timevariant channel impulse response recalculated each time step. Table 1 summarizes the main features of Steingass/Lehner’s channel model.
2.3. A Realistic Signal/Channel Simulator
When transmitted, satellite’s signals travel through a propagation channel which modifies its amplitude, phase, and delay. Indeed, many replicas of the same transmitted signal can reach the receiver’s antenna due to multipath propagation. In general, these replicas are caused by reflections of the direct signal in surrounding obstacles (e.g., buildings, trees, and ground etc.). As shown above, such propagation channel is generically modeled by a linear timevarying impulse response with propagation paths: where , and are the amplitude, phase, and delay of the th propagation path for the th satellite, is the multipath delay axis and the index stands for the lineofsight signal. These channel parameters can be seen as realizations of random processes with underlying probability density functions , , and , respectively, whose shape and parameters are approximated by the models outlined above.
Therefore, considering visible satellites, the signal received at the receiver’s antenna is the superposition of the transmitted signals, as propagated through the corresponding channel, and corrupted by additive noise, . This reads as where is the transmitted signal corresponding to the th satellite.
As shown in [12], the term can be approximated by its firstorder Taylor expansion as . Hence, the general baseband equivalent model that will be used along this paper is
The first element in the receiver RF chain is a right hand circularly polarized (RHCP) antenna, usually with nearly hemispherical gain coverage, with the mission to receive the radionavigation signals of all the satellites in view. The RF signals collected by the antenna are immediately amplified by a low noise amplifier (LNA), a key element which is the most contributing block to the noise figure of the receiver. The LNA also acts as a filter, minimizing outofband RF interferences and setting the sharpness of the received code. After the LNA, the amplified and filtered RF signals are then downconverted to an intermediate frequency (IF) using signal mixing frequencies from local oscillators (LOs). These LOs are derived from a receiver reference oscillator, often an ovenstabilized clock with typical accuracies of . There is a need for one LO per downconversion stage. Two or three downconversion stages are commonly devoted to reject mirror frequencies or large out of band jamming signals, in particular the 900 MHz used by the GSM mobile communication system. However, depending on the subsequent analogtodigital converter (ADC) characteristics, a onestage downconversion or even a direct Lband sampling is also possible [13]. The lower sideband generated by the mixer process is selected, while the upper sideband is filtered by a postmixer bandpass filter. It is important to point out that signal Doppler’s and PRN codes are preserved after the mixing stage, only the carrier frequency is lowered.
In the sequel, we focus on the contribution of a single satellite and thus omit the dependence with of the signal model. Considering a generic data sequence , chip code , chipshaping pulse , chip period , full codes in a whole bit, and data period , the baseband equivalent received signal for a channel model as in (7) but particularized to (i.e., only one line of sight signal) can be put in the form where is the pulse received at the antenna and then filtered by a precorrelation filter (usually the LNA), is the filtered version of , and the term stands for the filtered thermal noise and other unmodeled terms. The objective of a synchronization method is to estimate the time delay , Doppler shift and the carrier phase information embedded into the phase of the complex amplitude .
The analogtodigital conversion and the automatic gain control (AGC) processes take place at IF or baseband, where all the signals from GNSS satellites in view are buried in thermal noise. Once the received signal is digitized, it is ready to feed each of the digital receiver channels. Every receiver channel is intended to acquire and track the signal of a single GNSS satellite; typical receivers are equipped with channels. The multiplication of the IF digitized signal by a local replica of its carrier frequency allows to produce the inphase (I) and quadraturephase (Q) components of the digitized signal.
Assuming as additive white Gaussian noise (AWGN), at least in the band of interest, it is well known that the optimum receiver is the code matched filter, expressed as where are local estimates of the time delay, Doppler shift, and carrier phase of the received signal, and stands for the complex conjugate operator. Theoretically , but actual implementations make use of approximated versions: while is a rectangular pulse filtered at the satellite, is digitally generated at the receiver and therefore not filtered. In addition, is usually filtered again by a precorrelation filter before the matched filter, as expressed in (10) with . The code matched filter output can be written in the form
Notice that, in the matched filter, we have substituted the estimates , , and for trial values obtained from previous (in time) estimates of these parameters which we have defined as , , and , respectively. This is the usual procedure in GNSS receivers, since the estimates are not really available, but to be estimated after correlation.
In DSSS terminology, the matched filter is often referred to as correlator, while the processing it performs is called despreading. Since the correlators perform accumulation of the sampled signal during a period and then release an output, we can write the discrete version of the signal as where is the sampling period, is the integration time (usually, ) and stands for the nearest integer towards zero.
Equation (13) can be expressed more conveniently by solving the convolution in (12), which yields [14] where we defined , and (i.e., the estimation errors), stands for the nearest integer toward zero, and means the integer part of , being the navigation bit period, and is the correlation function. An equivalent derivation for the arm leads to
Terms , , and should be regarded as the average local phase error over the integration interval, that is, , assuming a frequency rate error (i.e., a phase acceleration error) equal to zero. In case of inclusion of such effect in the model, the average phase error can be expanded as In this expression, the terms , , and are referred to the error values at the beginning of the integration interval.
In the following, we will consider as the integer number of samples collected in an accumulation. This number will not be integer in receiver configurations having a sample rate incommensurable with the chip rate, and thus some integration blocks will have samples instead of . This effect can be considered negligible for the analysis presented in this paper.
In the case of (i.e., in the presence of multipath), (12) becomes a sum of all the replicas convoluted with a filter matched to the line of sight signal, whose estimated parameters are possibly biased by the presence of multipath. Since the convolution is a linear operator, the correlator output will be a linear combination of the contributions made by each signal path.
Note that an arbitrary number of correlators (very early, early, prompt, late, very late, etc.) can be used in the filter update, just adding or subtracting the correlator offset to the argument of (i.e., , , etc.). The correlators’ output can be stacked in a vector , which will be the measurements used in next section.
In the context of this work, we used the GRANADA (Galileo Receiver ANAlysis and Design Application) simulation platform to simulate realistic channel and receiver scenarios. The GRANADA Factored Correlator Model (FCM) blockset (see Figure 2) is a MATLAB/Simulink (MATLAB and Simulink are registered trademarks of The MathWorks, Inc.) library that provides a swift, flexible, and realistic way of simulating different signal processing architectures, either of standalone GNSS receivers or multisystem solutions. The FCM was included in a Simulink blockset, which, since 2007, has been commercially available as part of the GRANADA product family, whose remaining products were developed by DEIMOS Space in the frame of the Galileo Receiver Development activities (GARDA), funded by the Galileo Joint Undertaking (now European GNSS Agency, GSA) under the 6th Framework Program of the European Union.
The FCM separates the effects of carrier and code Doppler and misalignment on a GNSS receiver’s correlator outputs into several multiplicative factors and allows the inclusion (or not) of each factor independently. Since it is an analytical model, the computation rate can be as low as the tracking loop rate, dramatically increasing simulation speed: the FCM provides directly the correlators’ output, precluding the need of simulating the lowerlevel signal processing stages, significantly reducing the computational load and hence decreasing processing and memory requirements, while still accounting for various effects (as filtering, carrier phase and frequency errors, code delay error, code Doppler, noise, and multipath), thus keeping a high level of realism [15]. Since, statistically speaking, it is equivalent to work with samples before or after the correlation process (proof in the Appendix), we take advantage of working at the correlator output since it considerably reduces the computational load.
Once configured (type of signal, propagation channel, user dynamics, sampling frequency before correlation, number of correlators and their spacing, integration period, environment, etc., see Figure 3), FCM provides the measurements used in the simulations presented in Section 5.
3. Particle Filtering
Bayesian filtering involves the recursive estimation of states given measurements at time based on all available measurements, . To that aim, we are interested in the filtering distribution , which can be recursively expressed as with and referred to as the likelihood and the prior distributions, respectively. Unfortunately, (18) can only be obtained in closedform in some special cases. For instance, when the model is linear and Gaussian, the Kalman Filter (KF) [16] provides the optimal solution. In more general setups—nonlinear and/or nonGaussian—we should resort to more sophisticated methods [17]. In this paper we consider particle filters (PFs) [18, 19].
PFs approximate the filtering distribution by a set of weighted random samples, forming the socalled set of particles . These random samples are drawn from the importance density distribution, , and weighted according to the general formulation
Algorithm 1 outlines the operation of the Standard PF (SPF) when a new measurement becomes available. After particle generation, weighting, and normalization, a minimum mean square error (MMSE) estimate can be obtained by a weighted sum of particles. A typical problem of PFs is the degeneracy of particles, where all but one weight tend to zero. This situation causes the particle to collapse to a single state point. To avoid the degeneracy problem, we apply resampling, consisting in eliminating particles with low importance weights and replicating those in highprobability regions [20, 21]. In this work, we consider a multinomial sampling scheme for the resampling step.

3.1. RaoBlackwellized Particle Filter
In this paper, we analyze a way to alleviate the dimensionality problem based on the marginalization of linear states. The basic idea is that a KF can optimally deal with these states, while reducing the dimension of the state space that the nonlinear filter has to explore. The procedure was proposed in [23, 24] for the case of dealing with the nonlinear states with a PF. The algorithm was termed Marginalized particle filter (MPF), although the same concept is also referred to as RaoBlackwellized PF (RBPF) in other works [25, 26]. The latter nomenclature is because marginalization resorts to a general result due to [27, 28] referred to as the RaoBlackwell theorem, which shows that the performance of an estimator can be improved by using information about conditional probabilities. The RaoBlackwell theorem states let represent any unbiased estimator for and be a sufficient statistic for under . Then the conditional expectation is independent of , and it is the uniformly minimum variance unbiased estimator (cf. [29, 30] for the details) The result of a corollary points out that the use of a RaoBlackwellized estimator effectively reduces the variance of the estimation error. Therefore, when possible, it is desirable to apply marginalization procedures. Corollary: let be an unbiased estimator and let be the RaoBlackwell estimator, then
Final remarks on RaoBlackwellization are worth mentioning. (i)RaoBlackwellization is a procedure suitable when linear substructures are present in the dynamical model. (ii)It is a variance reduction technique, in the sense that the estimation variance of a filter considering this marginalization procedure is less than a filter estimating the complete state space. (iii)Filtering linear states with a Kalman filter has twofold benefits: (1) linear states are optimally filtered and (2) the system coped by the nonlinear filter has reduced dimensionality (with large benefits in terms of computational resources).
4. Joint Filtering of LOSS and Multipath Parameters
The technique herein investigated attempts to estimate the synchronization parameters of both the LOSS and multipath components. We refer to the algorithm as the multipath estimating particle filtering, or MEPF for short. Here the term Bayesian means that the algorithm is using some sort of a priori information regarding these parameters (such as interdependencies and time evolution models). This approach was first introduced in [31] and further refined in [32], although other papers might be found following the same scheme [33] with more complex timeevolving models. The application of Bayesian filtering techniques becomes straightforward when one describes the problem at hand in terms of a measurement equation and a process equation (i.e., how unknowns evolve randomly over time).
4.1. Observations
A receiver implementing such Bayesian tracking loops typically processes each satellite independently, and most of the work in the literature discusses architectures using IF signal. Here we are interested in operating at the output of the bank of correlators.
Observations for the th satellite are gathered into a random vector , where we omitted the subindex for the sake of clarity. The th element in corresponds to the sample of the th correlator, and it is expressed as accounting that corresponds to the point where the th early/late sample is evaluated. As usual, denotes LOSS. Here we consider a noncoherent tracking architecture that operates with the squared outputs. This scheme avoids the estimation of carrier phases, and thus it reduces the statespace dimension. In our implementation, a conventional PLL/FLL network is used in parallel to the MEPF. Therefore, the observations are the parallel outputs of the correlation bank, which we denote as where is the total number of correlators used at the receiver. We made apparent the dependence of measures on unknown states: real amplitude and time delay of each replica of the signal.
4.2. Process Dynamics
The state space is composed of the unknown parameters of the model, namely, delay, delay rate, and real amplitude of the LOSS and its multipath replica: where is the delay rate of the th component, related to the Doppler shift. We have introduced this delay rate to better capture the dynamics of the timeevolving delay of the signals.
One could adopt many alternatives to specify the timeevolving processes for each state, ranging from the simplistic (although effective in some situations) autoregressive model to more sophisticated models. Here, we adopt a channel state model based on that presented in [34], adapted to the noncoherent scheme. This model was motivated by channel modeling work for multipath prone environments such as the urban satellite navigation channel [35].
The dynamics of time delay and delay rate for the LOSS (i.e., ) are described by where is the integration period and the process noise is an uncorrelated zeromean Gaussian random variable with diagonal entries and .
The evolution of and for the echoes is modeled with a truncated Gaussian distribution as in [31], which allows us to introduce the fact that due to physical reasons in outdoor propagation channels [6, 11, 36]. Taking (26) into account, we force this situation using the evolution with and being zeromean Gaussian random variables with variances and , respectively. For the evolution of each we consider independent autoregressive models with variance . The overall covariance matrix of the process is denoted as and is constructed with , , , , , and in its diagonal.
4.3. Algorithm Implementation
From the previous modeling, we realize that the state space can be partitioned into linear and nonlinear subspaces. Clearly, these can be identified as
By the chain rule of probability, linear states can be analytically marginalized out from : and, taking into consideration that generates a linear Gaussian statespace, can be updated analytically via a KF conditional on and only the nonlinear part of needs to be estimated with a nonlinear filter. In the proposed scheme, an SPF is run to characterize and a KF is executed to obtain .
Notice that both linear and nonlinear states are interdependent, thus the algorithm has to be aware of this coupling. The details might be consulted in [23] for the general algorithm and in [12] for the specific GNSS setup considered here. At a glance, each particle in the PF has an associated KF that tracks amplitudes. Then, before particle generation, KF prediction is run and the results are used in the particle filter. Similarly, once particles are weighted this information is used in the update step of the KF.
5. Results in Realistic Scenarios
We used the GRANADA FCM blockset of Simulink to simulate the GPS L1 C/A signal, the propagation channel, and the inaccuracies of the receiver front end. An initial set of controlled scenarios is simulated to analyze the method. Then, from the set of reviewed channel models, we have selected Jahn’s to show simulation results in a realistic environment. The GPS signal is spread spectrum with a code length of chips and a chip rate of Mchips/s (notice that a chip of the signal corresponds to approximately 300 meters in length and the duration of an entire codeword is one millisecond). The carrier frequency of the transmitted signal was MHz and the receivers precorrelation bandwidth was MHz. Estimates of time delay were performed at a rate of Hz, which corresponds to an integration time of milliseconds, assuming bit synchronization. The carriertonoise density ratio () of the simulated satellite was dBHz. The dynamics of the scenario were due to the relative motion of the satellitereceiver, which is completely simulated by the GRANADA FCM blockset, and the receiver performed a pedestrianlike trajectory at m/s. Simulation time was seconds.
We compared the performance of the MEPF with the results of a narrow chip spacing DLL (stateoftheart in GNSS receivers) with an equivalent noise bandwidth of Hz. This architecture uses correlators. Also, the benchmark receiver implements a coherent phase lock loop (PLL) carrier phase discriminator using a secondorder filter and an error accumulator with equivalent noise bandwidth Hz. The initial timedelay ambiguity at which the filter was initialized was drawn from , with the chip period.
It has been reported in [37] that the number of correlators () used in the PF plays an important role. For instance, in AWGN on the order of correlators are required to obtain stable results. Also, the algorithm improves its performance with the number of particles although this improvement saturates at particles.
Figures 4–7 show the behavior of the classical DLL/PLL scheme and the proposed MEPF, respectively, in a multipath scenario. In this experiments, we used correlators for the MEPF in order to span correlators along regions of interest in terms of multipath estimation and mitigation. The results are organized as follows. Top figure represents the obtained pseudorange error. Central figure is the relative delay between the LOSS and the multipath replica, in the first representative interval () it has been set to chips and in the second interval () to chips. Bottom figure plots the signaltomultipath ratio (SMR) in linear scale of the simulated scenario. During the first interval, the SMR was abruptly kept constant to and during the second interval it grew linearly from to . Since the MEPF is very sensitive to the tuning of process covariance matrix—as many Bayesian filtering solutions,—we have investigated three different setups with particles. Namely, (i) in Figure 5 we used standard deviations , , , , and ; (ii) in Figure 6 we used , , , ,, and ; and finally (iii) in Figure 7 we used , , , ,, and . At the light of the results, the latter configuration provided a good performance as it allowed for sufficient delay excursions to explore the state space and fast variations in multipath amplitude were coped. A summary of results in terms of bias, variance, and RMSE over the entire simulation can be consulted in Table 2. We can observe that, compared to DLL schemes, a remarkable performance improvement can be obtained after properly adjusting the covariances.

Finally, we tested the algorithm in a more realistic scenario. We selected the Jahn’s channel model with the same receiver parameters as before. Particularly, the considered channel was that of a satellite at an elevation angle of in an urban scenario with an average of dBHz. The results can be consulted in Figure 8, where it can be observed that MEPF requires an initial convergence time (depending on the covariance matrix set) larger than DLL schemes. Conversely, it appears more robust to channel impairments. Numerically, the RMSE in the overall simulation is of m and m for DLL and MEPF, respectively. For the MEPF we used paths, particles, and , , , , , and .
6. Conclusions
In this paper we have analyzed an advanced tracking loop for timedelay and carrierphase estimation in a GNSS receiver based on sequential MonteCarlo methods. The algorithm builds upon previous work by the authors on RaoBlackwellized particle filtering while introducing more realistic process dynamics and the usage of postcorrelation observations, that reduce the computational burden at the receiver. The paper presents the general signal model, GNSS concept, and tradeoffs the most common propagation channel models. A realistic scenario simulator based on the FCM blockset of Simulink was used Section 5. Results point out the need for properly setting not only the number of particles but the number of correlation outputs used as observations. Also, degradation of conventional DLL/PLL schemes in multipathrich scenarios became clear. Nevertheless, the correct selection of a process covariance matrix was seen to affect significantly the performance of the MEPF and future work should be devoted in selfadjustment of such matrix.
Appendix
Equivalence of Pre/Postcorrelation Receiver Architectures
In this appendix we establish a basic result showing the equivalence between processing pre and postcorrelation signals. That is to say, from a statistical point of view, an estimator of a given parameter (e.g., time delay) computed using a bunch of snapshots taken at the IF signal level () is the same as that which is derived using the output of the correlators (). It is a wellknown result in statistical signal processing that both signals are sufficient statistics, and thus one is able to derive an estimator of using either. However, we will see that this equivalence becomes evident when one examines the likelihood distribution (the density where the information from measurements is gathered) for each approach.
If we analyze first the case of using the IF signal we should be aware of the following.(i)This approach does not force an implementation based on early, prompt, and late samples; as observations are directly the baseband signal at the sampling frequency. (ii)It is necessary to use a sufficiently large set of IF data to be able to infer any parameter from it. That is, one has to integrate over a certain integration time, , since the signaltonoise ratio of GNSS signals is typically well below the noise level.
The term stands for the vector of snapshots of the IF signal, as gathered for the th integration interval, defined as using the same notation conventions used along the document. Then, the likelihood can be decomposed as the independent contribution of each snapshot and assuming Gaussianity for the noise term, we could identify that where stands for the precorrelation signal model, which was defined earlier as
Further manipulation of the loglikelihood yields to with the latter step being clear if one accounts for the definition of as the output of a correlator. Recall that is the true unknown parameter of the signal. An ML estimator of could be obtained after maximizing the latter equation.
Drawbacks of this approach are twofold. (i)It might be computationally expensive as large data sets need to be processed to increase the signaltonoise ratio, and thus might be large depending on . (ii)There is a requirement for performing signal processing operations at a high rate, since it operates at the sampling frequency.
If we turn our attention to the conventional approach in which one uses samples at the output of a bank of correlators, we should see the following. (i)This approach forces an implementation based on early, prompt, and late samples; this means that samples are taken assuming a previous estimation (prompt) of the parameters, denoted as . (ii)Few samples are sufficient to infer estimates of . After correlation an integration over a certain interval is already done, , and therefore the signaltonoise ratio is relatively high.
In this case, measurements can be expressed as at the output of the th integration interval. In this measurement we explicitly expressed that samples are taken with respect to the error between true and prompt parameters, . Notice that we considered that only the prompt is used for the sake of clarity. It is easy to obtain a similar result, as the one shown here, when one accounts for several early and late samples.
Then, the loglikelihood under the Gaussian assumption is with being the postcorrelation signal model and the unknown parameter we want to estimate at .
If we set , we can identify that
From the latter mathematical derivations, we can conclude an important result: for a given integration interval considering snapshots. As said, similar results apply for larger integration and more early/late samples.
As a consequence, we can state the following: the ML estimator of computed from the data sets and is equivalent.
To sum up, from a statistical point of view, both approaches are equivalent and the choice should be made considering implementation aspects. For instance, it is clear that using precorrelation measurements involves larger computational burden than using postcorrelation samples. Another important conclusions is that since in the precorrelation approach we also need to integrate in order to increase the signaltonoise ratio, effects happening faster than will not be captured by the estimation algorithm. The same happens in the postcorrelation case. Therefore, the limitation of which phenomena could be tracked is inherent to the GNSS signal, instead of the way it is processed (i.e., preor postcorrelated samples).
Acknowledgment
P. Closas and C. FernándezPrades were supported by the European Commission under COST Action IC0803 (RFCSET).
References
 M. HernándezPajares, J. M. J. Zornoza, J. S. Subirana, R. Farnworth, and S. Soley, “EGNOS test bed ionospheric corrections under the October and November 2003 storms,” IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 10, pp. 2283–2293, 2005. View at: Publisher Site  Google Scholar
 S. Bancroft, “An algebraic solution of the GPS equations,” IEEE Transactions on Aerospace and Electronic Systems, vol. 21, no. 1, pp. 56–59, 1985. View at: Google Scholar
 G. Strang and K. Borre, Linear Algebra, Geodesy, and GPS, Wellesley Cambridge Press, 1997.
 B. HofmannWellenhof, H. Lichtenegger, and E. Wasle, GNSS Global Navigation Satellite Systems: GPS, GLONASS, Galileo & More, SpringerVerlag, Wien, Austria, 2008.
 C. FernándezPrades, L. L. Presti, and E. Falletti, “Satellite radiolocalization from GPS to GNSS and beyond: novel technologies and applications for civil mass market,” Proceedings of the IEEE, vol. 99, no. 11, pp. 1882–1904, 2011. View at: Publisher Site  Google Scholar
 A. Jahn, H. Bischl, and G. Heiss, “Channel characterization for spread spectrum satellite communications,” in Proceedings of the 4th International Symposium on Spread Spectrum Techniques & Applications (ISSSTA '96), pp. 1221–1226, September 1996. View at: Google Scholar
 C. Loo and J. S. Butterworth, “Land mobile satellite channel measurements and modeling,” Proceedings of the IEEE, vol. 86, no. 7, pp. 1442–1462, 1998. View at: Google Scholar
 M. A. V. Castro, F. P. Fontan, A. A. Villamarín, S. Buonomo, P. Baptista, and B. Arbesser, “Lband Land Mobile Satellite (LMS) amplitude and multipath phase modeling in urban areas,” IEEE Communications Letters, vol. 3, no. 1, pp. 12–14, 1999. View at: Google Scholar
 F. P. Fontán, M. VázquezCastro, C. E. Cabado, J. P. García, and E. Kubista, “Statistical modeling of the LMS channel,” IEEE Transactions on Vehicular Technology, vol. 50, no. 6, pp. 1549–1567, 2001. View at: Publisher Site  Google Scholar
 A. Steingass and A. Lehner, “A channel model for land mobile satellite navigation,” in Proceedings of the the European Navigation Conference, pp. 2132–2138, German Institute of Navigation (DGON), July 2005. View at: Google Scholar
 Recommendation ITUR P.6817, “Propagation data required for the design of Earthspace land mobile telecommunication systems,” 2009, http://www.itu.int/rec/RRECP.6817200910I/en. View at: Google Scholar
 P. Closas, Bayesian signal processing techniques for GNSS receivers: from multipath mitigation to positioning [Ph.D. dissertation], Universitat Politècnica de Catalunya (UPC), Department of Signal Theory and Communications, Barcelona, Spain, 2009.
 D. M. Akos, M. Stockmaster, J. B. Y. Tsui, and J. Caschera, “Direct bandpass sampling of multiple distinct RF signals,” IEEE Transactions on Communications, vol. 47, no. 7, pp. 983–988, 1999. View at: Publisher Site  Google Scholar
 B. Parkinson and J. Spilker, Eds., Global Positioning System: Theory and Applications, vol. 1 of Progress in Astronautics and Aeronautics, American Institute of Aeronautics, Washington, DC, USA, 1996.
 J. S. Silva, P. F. Silva, A. Fernández, J. Diez, and J. F. M. Lorga, “Factored correlator model: a solution for fast, flexible, and realistic GNSS receiver simulations,” in Proceedings of the 20th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS '07), pp. 2676–2686, Fort Worth, TX, USA, September 2007. View at: Google Scholar
 R. E. Kalman, “A new approach to linear filtering and prediction problems,” Transactions of the ASMEJournal of Basic Engineering, vol. 82, pp. 35–45, 1960. View at: Google Scholar
 Z. Chen, “Bayesian filtering: from Kalman filters to particle filters, and beyond,” Tech. Rep., Adaptive Systems Laboratory, McMaster University, Ontario, Canada, 2003. View at: Google Scholar
 M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/nonGaussian Bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174–188, 2002. View at: Publisher Site  Google Scholar
 P. M. Djurić, J. H. Kotecha, J. Zhang et al., “Particle filtering,” IEEE Signal Processing Magazine, vol. 20, no. 5, pp. 19–38, 2003. View at: Publisher Site  Google Scholar
 M. Bolić, P. M. Djuric, and S. Hong, “Resampling algorithms for particle filters: a computational complexity perspective,” Eurasip Journal on Applied Signal Processing, vol. 2004, no. 15, pp. 2267–2277, 2004. View at: Publisher Site  Google Scholar
 R. Douc, O. Cappé, and E. Moulines, “Comparison of resampling schemes for particle filtering,” in Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis (ISPA '05), pp. 64–69, Zagreb, Croatia, September 2005. View at: Google Scholar
 GRANADA Galileo Receiver ANalysis And Design Application. The Reference Galileo Simulation Toolkit for GNSS Receiver Research And Development. Factored Correlator Model Blockset v2.0 User Manual, Deimos Engenharia, S.A., 2009.
 T. Schön, F. Gustafsson, and P. J. Nordlund, “Marginalized particle filters for mixed linear/nonlinear statespace models,” IEEE Transactions on Signal Processing, vol. 53, no. 7, pp. 2279–2289, 2005. View at: Publisher Site  Google Scholar
 R. Karlsson, Particle filtering for positioning and tracking applications [Ph.D. dissertation], Linköping University, Linköping, Sweden, 2005.
 R. Chen and J. S. Liu, “Mixture Kalman filters,” Journal of the Royal Statistical Society B, vol. 62, no. 3, pp. 493–508, 2000. View at: Google Scholar
 A. Doucet, N. de Freitas, and N. Gordon, Eds., Sequential Monte Carlo Methods in Practice, Springer, 2001.
 C. Rao, “Information and the accuracy attainable in the estimation of statistical parameters,” Bulletin of Calcutta Mathematical Society, vol. 37, pp. 81–91, 1945. View at: Google Scholar
 D. Blackwell, “Conditional expectation and unbiased sequential estimation,” The Annals of Mathematical Statistics, vol. 18, no. 1, pp. 105–110, 1947. View at: Google Scholar
 E. Lehmann, Theory of Point Estimation. Probability and Mathematical Statistics, John Wiley & Sons, 1983.
 A. Papoulis and S. U. Pillai, Probability, Random Variables and Stochastic Processes, McGrawHill, New Delhi, India, 4th edition, 2001.
 P. Closas, C. FernándezPrades, and J. A. FernándezRubio, “Bayesian DLL for multipath mitigation in navigation systems using particle filters,” in Proceedings of the IEEE (ICASSP '06), Toulouse, France, May 2006. View at: Google Scholar
 P. Closas, C. FernándezPrades, and J. A. FernándezRubio, “A Bayesian approach to multipath mitigation in GNSS receivers,” IEEE Journal on Selected Topics in Signal Processing, vol. 3, no. 4, pp. 695–706, 2009. View at: Publisher Site  Google Scholar
 M. Lentmaier, B. Krach, and P. Robertson, “Bayesian time delay estimation of GNSS signals in dynamic multipath environments,” International Journal of Navigation and Observation, vol. 2008, Article ID 372651, 11 pages, 2008. View at: Publisher Site  Google Scholar
 B. Krach, P. Robertson, and R. Weigel, “An efficient twofold marginalized Bayesian filter for multipath estimation in satellite navigation receivers,” Eurasip Journal on Advances in Signal Processing, vol. 2010, Article ID 287215, 2010. View at: Publisher Site  Google Scholar
 A. Steingass and A. Lehner, “Measuring the navigation multipath channel—a statistical analysis,” in Proceedings of the 17th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS '04), pp. 1157–1164, Long Beach, Calif, USA, September 2004. View at: Google Scholar
 M. Irsigler, J. A. ÁvilaRodríguez, and G. W. Hein, “Criteria for GNSS multipath performance assessment,” in Proceedings of the International Technical Meeting of the Institute of Navigation(ION GPS/GNSS '05), Long Beach, Calif, USA, September 2005. View at: Google Scholar
 P. Closas, C. FernándezPrades, J. Diez, and D. de Castro, “Multipath estimating tracking loops in advanced GNSS receivers with particle filtering,” in Proceedings of the IEEE Aerospace Conference, Big Sky, Mont, USA, March 2012. View at: Google Scholar
Copyright
Copyright © 2012 Pau Closas et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.