Nonlinear Bayesian Tracking Loops for MultipathMitigation

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 Rao-Blackwellization of linear states and a particle filter for the nonlinear partition and compare it to traditional delay lock loop/phase lock loop-based schemes.


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 well-known 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 in-view 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 line-of-sight 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 carrier-phase estimations.Multipath is probably the dominant source of error in high-precision applications, especially in urban scenarios, since it can introduce a bias up to a hundred of meters when employing a 1-chip wide (standard) delay lock loop (DLL) to track the delay of the LOSS, which is a common synchronization method used in spread-spectrum 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 bias-free 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 state-ofthe-art 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 Rao-Blackwellization of signal amplitudes and the use of a suitable nonlinear filter for the rest of nonlinear states, for example, time-delays and their rate.More precisely, Rao-Blackwellization 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 Monte-Carlo methods (specifically, the standard particle filtering) as one of the most promising alternatives in advanced International Journal of Navigation and Observation 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.

Fundamentals of Global Navigation Satellite Systems
GNSS space vehicles broadcast a low-rate 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 c = 299792458 m/s is the speed of light, t Rx i is the receiving time in the receiver's clock, and t Tx i the time of transmission for a given satellite i. 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 s, where s = {GPS, Galileo, GLONASS, . ..}, 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 Δt (s) .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 c 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 60 km until 2000 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 in-view satellite i of system s, we can write where (x Tx i , y Tx i , z Tx i ) is the satellite's position (known from the navigation message), (x, y, z) the receiver's position, and σ e 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 3 + N s satellites must be seen by the receiver at the same time, where N s is the number of different navigation systems available (in-view) 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 Newton-Raphson 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 multi-constellation 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 in-view satellites and using them jointly in the navigation solution, therefore increasing accuracy.

Signal Model.
A general signal model for most navigation systems consists of a direct-sequence spread-spectrum (DS-SS) 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 International Journal of Navigation and Observation 3 where being P T the transmitting power, γ a parameter controlling the power balance, d I (m) ∈ {−1, 1} the data symbols, T bI the bit period, N cI the number of repetitions of a full codeword that spans a bit period, T PRNI = T bI /N cI the codeword period, c I (k) ∈ {−1, 1} a chip of a spreading codeword of length L cI chips, g T,I (t) the transmitting chip pulse shape, which is considered energy normalized for notation clarity, and T cI = T bI /N cI L cI is the chip period.Figure 1 aims at clarifying the relation between those bits/chips parameters.Subindex I refers to the in-phase component, and all parameters are equivalently defined for the quadrature component, referred to with the subindex Q.This signal model describes all GNSS's signals-in-space, for instance GPS L1, GPS L5, Galileo E1, and Galileo E5.Refer to [5] for the details.

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 T 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 T < 0.5 μ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 320 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, T coh , is defined as the time interval during which the characteristics of the propagation channel remain approximately constant, and it is given as where f m is the maximum Doppler shift.The Doppler shift is given as v/λ, where v is the radial speed of the mobile terminal with respect to the satellite and λ is the signal wavelength.A channel is considered WSSUS (wide-sense 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.

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.

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 19 • elevation.Many more investigations on L-band 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.

Pérez-Fontán's Channel Characterization
where α i,m , φ i,m (t) and τ i,m (t) are the amplitude, phase, and delay of the mth propagation path for the ith satellite, ξ is the multipath delay axis and the index m = 0 stands for the line-of-sight signal.These channel parameters can be seen as realizations of random processes with underlying probability density functions f αp (α), f φp (φ), and f τp (τ), respectively, whose shape and parameters are approximated by the models outlined above.
Therefore, considering M s visible satellites, the signal r(t) received at the receiver's antenna is the superposition of the transmitted signals, as propagated through the corresponding channel, and corrupted by additive noise, w(t).This reads as where s T,i (t) is the transmitted signal s T (t) corresponding to the i-th satellite.
As shown in [12], the term φ i,m (t) can be approximated by its first-order Taylor expansion as φ i,m (t) ≈ 2π f di,m (t)t + φ i,m,0 .Hence, the general baseband equivalent model that will be used along this paper is • e jφi,m,0 s T,i t − τ i,m (t) + w(t).(9) 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 out-of-band 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 oven-stabilized clock with typical accuracies of 10 −8 .There is a need for one LO per down-conversion stage.Two or three down-conversion 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 analog-to-digital converter (ADC) characteristics, a one-stage downconversion or even a direct L-band 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 i of the signal model.Considering a generic data sequence d, chip code c, chip-shaping pulse g T (t), chip period T c , N c full codes in a whole bit, and data period T b , the baseband equivalent received signal for a channel model as in (7) but particularized to M i = 1 (i.e., only one line of sight signal) can be put in the form where g R (t) is the pulse received at the antenna and then filtered by a precorrelation filter (usually the LNA), p(t) is the filtered version of p(t) = p I (t) + j p Q (t), and the term w(t) stands for the filtered thermal noise and other unmodeled terms.The objective of a synchronization method is to estimate the time delay τ 0 (t), Doppler shift f d (t) and the carrier phase information φ 0 embedded into the phase of the complex amplitude a 0 (t) = |a 0 (t)|e j(2π fd(t)t+φ0) .The analog-to-digital 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 N digital receiver channels.Every receiver channel is intended to acquire and track the signal of a single GNSS satellite; typical receivers are equipped with N = 12 channels.The multiplication of the IF digitized signal by a local replica of its carrier frequency allows to produce the in-phase (I) and quadrature-phase (Q) components of the digitized signal.
Assuming w(t) 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 • e − j φ0 e − j2π fd 0 (t)t = q * R (−t + τ 0 + L c T c )e − j φ0 e − j2π fd 0 (t)t , (11) where τ 0 , f d0 , φ 0 are local estimates of the time delay, Doppler shift, and carrier phase of the received signal, and (dot) * stands for the complex conjugate operator.Theoretically g R (t) = g T (t), but actual implementations make use of approximated versions: while g T (t) is a rectangular pulse filtered at the satellite, g R (t) is digitally generated at the receiver and therefore not filtered.In addition, g T (t) is usually filtered again by a precorrelation filter before the matched filter, as expressed in (10) with g R (t).The code matched filter output can be written in the form y t; τ0 , fd0 , φ0 = r 0 t; τ 0 , f d0 , φ 0 * h MF t; τ0 , fd0 , φ0 .(12) Notice that, in the matched filter, we have substituted the estimates τ 0 , f d0 , and φ 0 for trial values obtained from previous (in time) estimates of these parameters which we have defined as τ0 , fd0 , and φ0 , respectively.This is the usual procedure in GNSS receivers, since the estimates are not really available, but to be estimated after correlation.
In DS-SS 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 T int and then release an output, we can write the discrete version of the signal as where T s is the sampling period, T int is the integration time (usually, T int = T c L c ) 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 Δ f = f d0 − fd0 , Δφ = φ 0 − φ0 and Δτ 0 = τ 0 − τ0 (i.e., the estimation errors), • stands for the nearest integer toward zero, and [n] Tb/Tint means the integer part of nT int /T b , being T b the navigation bit period, and is the correlation function.An equivalent derivation for the Q arm leads to Terms Δ f , Δφ, and Δτ 0 should be regarded as the average local phase error over the integration interval, that is, Δφ = Δφ + 2πΔ f (T int /2), assuming a frequency rate error Δ ḟ International Journal of Navigation and Observation (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 Δφ, Δ f , and Δ ḟ are referred to the error values at the beginning of the integration interval.
In the following, we will consider K = (T int /T s ) 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 K + 1 samples instead of K.This effect can be considered negligible for the analysis presented in this paper.
In the case of M i > 1 (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 R pq (i.e., R pq (Δτ 0,n + δ), R pq (Δτ 0,n − δ), etc.).The correlators' output can be stacked in a vector y n , 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 lower-level 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 y n used in the simulations presented in Section 5.

Particle Filtering
Bayesian filtering involves the recursive estimation of states x n ∈ R nx given measurements y n ∈ R ny at time n based on all available measurements, y 1:n = {y 1 , . . ., y n }.To that aim, we are interested in the filtering distribution p(x n |y 1:n ), which can be recursively expressed as with p(y n |x n ) and p(x n |x n−1 ) referred to as the likelihood and the prior distributions, respectively.Unfortunately, (18) can only be obtained in closed-form 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 non-Gaussian-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 N weighted random samples, forming the so-called set of . 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 y n 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 high-probability regions [20,21].In this work, we consider a multinomial sampling scheme for the resampling step.

Rao-Blackwellized 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 High-frequency part of signal processing (complex)  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 Rao-Blackwellized 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 Rao-Blackwell theorem, which shows that the performance of an estimator can be improved by using information about conditional probabilities.The Rao-Blackwell theorem states let θ = g(x) represent any unbiased estimator for θ and T(x) be a sufficient statistic for θ under p(x, θ).Then the conditional expectation θ RB = E{g(x)|T(x)} 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 Rao-Blackwellized 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 θ RB be the Rao-Blackwell estimator, then Final remarks on Rao-Blackwellization are worth mentioning.
(i) Rao-Blackwellization is a procedure suitable when linear substructures are present in the dynamical model.
International Journal of Navigation and Observation (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).

Joint Filtering of LOSS and Multipath Parameters
The technique herein investigated attempts to estimate the synchronization parameters of both the LOSS and M − 1 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 time-evolving 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).

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 i-th satellite are gathered into a random vector y n , where we omitted the subindex i for the sake of clarity.The th element in y n corresponds to the sample of the -th correlator, and it is expressed as accounting that Δτ ,m,n = τ m,n − τ 0,n−1 + δ corresponds to the point where the -th early/late sample is evaluated.As usual, m = 0 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 state-space 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 L is the total number of correlators used at the receiver.We made apparent the dependence of measures on unknown states: real amplitude (α n ) and time delay (τ n ) of each replica m of the signal.

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 τm,n is the delay rate of the m-th component, related to the Doppler shift.We have introduced this delay rate to better capture the dynamics of the time-evolving 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., m = 0) are described by where T int is the integration period and the process noise is an uncorrelated zero-mean Gaussian random variable with diagonal entries σ 2 0,τ and σ 2 0, τ .The evolution of τ m,n and τm,n 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 τ m,n > τ 0,n ∀m ∈ {1, . . ., M − 1}, (26) in outdoor propagation channels [6,11,36].Taking (26) into account, we force this situation using the evolution with u τ m,n and u τ m,n being zero-mean Gaussian random variables with variances σ 2 m,τ and σ 2 m, τ , respectively.For the evolution of each α m,n we consider independent autoregressive models with variance σ 2 m,α .The overall covariance matrix of the process is denoted as Σ x and is constructed with σ 2 0,τ , σ 2 m,τ , σ 2 0, τ , σ 2 m, τ , σ 2 0,α , and σ 2 m,α in its diagonal.

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 p(x n | y 1:n ): and, taking into consideration that x l n generates a linear Gaussian state-space, p(x l n | x nl 0:n , y 1:n ) can be updated analytically via a KF conditional on x nl 0:n and only the nonlinear part of x n needs to be estimated with a nonlinear filter.In the proposed scheme, an SPF is run to characterize p(x nl 0:n | y 1:n ) and a KF is executed to obtain p(x l n | x nl 0:n , y 1:n ).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.

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 1023 chips and a chip rate of 1.023 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 1575.42MHz and the receivers precorrelation bandwidth was 2 MHz.Estimates of time delay were performed at a rate of 50 Hz, which corresponds to an integration time of 20 milliseconds, assuming bit synchronization.The carrier-to-noise density ratio (C/N 0 ) of the simulated satellite was 38 dB-Hz.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 pedestrian-like trajectory at 1 m/s.Simulation time was 50 seconds.
We compared the performance of the MEPF with the results of a narrow 0.125-chip spacing DLL (state-of-theart in GNSS receivers) with an equivalent noise bandwidth of 1 Hz.This architecture uses 3 correlators.Also, the benchmark receiver implements a coherent phase lock loop (PLL) carrier phase discriminator using a second-order filter and an error accumulator with equivalent noise bandwidth 10 Hz.The initial time-delay ambiguity at which the filter was initialized was drawn from N (τ 0,0 , T chip /2), with T chip the chip period.
It has been reported in [37] that the number of correlators (L) used in the PF plays an important role.For instance, in AWGN on the order of L = 11 correlators are required to obtain stable results.Also, the algorithm improves its performance with the number of particles although this improvement saturates at 300 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 L = 21 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 (t ∈ [10,20]) it has been set to 0.25 chips and in the second interval (t ∈ [25,40]) to 0.5 chips.Bottom figure plots the signal-to-multipath ratio (SMR) in linear scale of the simulated scenario.During the first interval, the SMR was abruptly kept constant to 0.75 and during the second interval it grew linearly from 0 to 0.75.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 N = 1000 particles.Namely, (i) in Figure 5 we used standard deviations σ 0,τ = .03/c,σ 0, τ = 0.03/c, σ m,τ = 100/c, σ m, τ = 0.03/c, σ 0,α = 0.0001, and σ m,α = 0.01; (ii) in Figure 6 we used σ 0,τ = 30/c, σ 0, τ = 0.3/c, σ m,τ = 30/c, σ m, τ = 0.3/c,σ 0,α = 0.0001, and σ m,α = 0.0001; and finally (iii) in Figure 7 we used σ 0,τ = 3/c, σ 0, τ = 0.3/c, σ m,τ = 30/c, σ m, τ = 0.3/c,σ 0,α = 0.0001, and σ m,α = 0.01.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 55 • in an urban scenario with an average C/N 0 of 38 dB-Hz.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 8.48 m and 4.82 m for DLL and MEPF, respectively.

Conclusions
In this paper we have analyzed an advanced tracking loop for time-delay and carrier-phase estimation in a GNSS receiver based on sequential Monte-Carlo methods.The algorithm builds upon previous work by the authors on Rao-Blackwellized 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 trade-offs 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 multipath-rich 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 self-adjustment of such matrix.

Figure 1 :
Figure 1: Relation among the parameters defining bits and spreading sequences in a generic navigation signal (in-phase component).
. The model presented by Fontán et al. in[9] addressed the statistical

Table 1 :
Comparison between channel models.

FCMFigure 2 :
Figure 2: Schematic of the tracking stage of a GNSS receiver and substitution of the high frequency stage of the receiver (correlation and carrier wipe-off) with the FCM blockset [22].

Figure 3 :
Figure 3: Configuration screen of the FCM blockset.

Table 1 :
Continued.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 first-order 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 narrow-band 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érez-Fontán's channel model.As shown above, such propagation channel is generically modeled by a linear time-varying impulse response with M i propagation paths: [11].4.Steingass/Lehner's Channel Characterization.The Steingass/Lehner land mobile channel model presented in[10]was developed using data recorded in a high-resolution 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 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.).