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 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.

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 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-of-the-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 so-called 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 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 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 so-called 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 𝜌_𝑖=𝑐𝑡^Rx_𝑖−𝑡^Tx_𝑖, (1) where 𝑐=299792458m/s is the speed of light, 𝑡^Rx_𝑖 is the receiving time in the receiver’s clock, and 𝑡^Tx_𝑖 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 𝑠={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 Δ𝑡^(𝑠). 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 60km until 2000km 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 𝑖 of system 𝑠, we can write 𝜌_𝑖=𝑥^Tx_𝑖−𝑥²+𝑦^Tx_𝑖−𝑦²+𝑧^Tx_𝑖−𝑧²+𝑐Δ𝑡^(𝑠)+𝜎_𝑒, (2) where (𝑥^Tx_𝑖,𝑦^Tx_𝑖,𝑧^Tx_𝑖) 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 3+𝑁_𝑠 satellites must be seen by the receiver at the same time, where 𝑁_𝑠 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.

2.1. 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 quasi-orthogonal) 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 𝑠_𝑇(𝑡)=𝑃_𝑇𝛾^∞_𝑢=−∞𝑑_𝐼(𝑢)𝑝_𝐼𝑡−𝑢𝑇_𝑏𝐼+𝑗1−𝛾^2∞_𝑙=−∞𝑑_𝑄(𝑙)𝑝_𝑄𝑡−𝑙𝑇_𝑏𝑄, (3) where 𝑝_𝐼(𝑡)=^{𝑁𝑐𝐼−1}_𝑛=0𝑞_𝐼𝑡−𝑛𝑇_PRN𝐼,𝑞_𝐼(𝑡)=^{𝐿𝑐𝐼−1}_𝑘=0𝑐_𝐼(𝑘)𝑔_𝑇,𝐼𝑡−𝑘𝑇_𝑐𝐼, (4) being 𝑃_𝑇 the transmitting power, 𝛾 a parameter controlling the power balance, 𝑑_𝐼(𝑚)∈{−1,1} the data symbols, 𝑇_𝑏𝐼 the bit period, 𝑁_𝑐𝐼 the number of repetitions of a full codeword that spans a bit period, 𝑇_PRN𝐼=𝑇_𝑏𝐼/𝑁_𝑐𝐼 the codeword period, 𝑐_𝐼(𝑘)∈{−1,1} 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 in-phase component, and all parameters are equivalently defined for the quadrature component, referred to with the subindex 𝑄. 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.

Figure 1 

Relation among the parameters defining bits and spreading sequences in a generic navigation signal (in-phase component).

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 𝐵_𝑐=1,2𝜋𝑇 (5) 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 𝑇<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 320kHz, 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, 𝑇_coh, is defined as the time interval during which the characteristics of the propagation channel remain approximately constant, and it is given as 𝑇_coh=12𝑓_𝑚, (6) 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 (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.

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.

Table 1 

Comparison between channel models.

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 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.

2.2.3. Pérez-Fontá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 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.

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 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 super-resolution 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 knife-edge 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 time-variant 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 time-varying impulse response with 𝑀_𝑖 propagation paths: ℎ_𝑖(𝑡;𝜉)=^𝑀𝑖−1_𝑚=0𝛼_𝑖,𝑚(𝑡)𝑒^{𝑗𝜙𝑖,𝑚(𝑡)}𝛿𝜉−𝜏_𝑖,𝑚(𝑡), (7) 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 𝑚=0 stands for the line-of-sight 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 𝑟(𝑡)=^𝑀𝑠−1_𝑖=0𝑠_𝑇,𝑖(𝑡)∗ℎ_𝑖()(t)=𝑡;𝜉+𝑤_{𝑀𝑠−1𝑖=0+∞−∞}𝑠_𝑇,𝑖(𝑡−𝜉)ℎ_𝑖()(𝑡)=𝑡;𝜉𝑑𝜉+𝑤_{𝑀𝑠−1𝑀𝑖=0𝑖−1𝑚=0}𝛼_𝑖,𝑚(𝑡)𝑒^{𝑗𝜙𝑖,𝑚(𝑡)}𝑠_𝑇,𝑖𝑡−𝜏_𝑖,𝑚(𝑡)+𝑤(𝑡), (8) where 𝑠_𝑇,𝑖(𝑡) is the transmitted signal 𝑠_𝑇(𝑡) corresponding to the 𝑖-th satellite.

As shown in [12], the term 𝜙_𝑖,𝑚(𝑡) can be approximated by its first-order Taylor expansion as 𝜙_𝑖,𝑚(𝑡)≈2𝜋𝑓_{𝑑𝑖,𝑚}(𝑡)𝑡+𝜙_𝑖,𝑚,0. Hence, the general baseband equivalent model that will be used along this paper is 𝑟(𝑡)=^{𝑀𝑠−1𝑀}_𝑖=0^𝑖−1_𝑚=0𝛼_𝑖,𝑚(𝑡)𝑒^{𝑗2𝜋𝑓𝑑𝑖,𝑚(𝑡)𝑡}⋅𝑒^{𝑗𝜙𝑖,𝑚,0}𝑠_𝑇,𝑖𝑡−𝜏_𝑖,𝑚(𝑡)+𝑤(𝑡). (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⁻⁸. 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 𝑖 of the signal model. Considering a generic data sequence 𝑑, chip code 𝑐, chip-shaping 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 𝑀_𝑖=1 (i.e., only one line of sight signal) can be put in the form 𝑟₀̃(𝑡)=𝑎₀(𝑡)^∞_𝑢=−∞𝑑(𝑢)^{𝑁𝑐−1𝐿}_𝑛=0^𝑐−1_𝑘=0𝑐(𝑘)⋅𝑔_𝑅𝑡−𝜏₀(𝑡)−𝑘𝑇_𝑐−𝑛𝑇_PRN−𝑢𝑇_𝑏+𝑤(𝑡)=𝑎₀(𝑡)^∞_𝑢=−∞𝑑(𝑢)𝑝𝑡−𝜏₀(𝑡)−𝑢𝑇_𝑏+𝑤(𝑡)=𝑎₀(𝑡)𝑒^{𝑗(2𝜋𝑓𝑑(𝑡)𝑡+𝜙0)∞}_𝑢=−∞𝑑(𝑢)𝑝𝑡−𝜏₀(𝑡)−𝑢𝑇_𝑏+𝑤(𝑡), (10) 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 𝑎₀(𝑡)=|𝑎₀(𝑡)|𝑒^{𝑗(2𝜋𝑓𝑑(𝑡)𝑡+𝜙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 𝑁 digital receiver channels. Every receiver channel is intended to acquire and track the signal of a single GNSS satellite; typical receivers are equipped with 𝑁=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 𝑤(𝑡) 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 ℎ_MF𝑡;𝜏₀,𝑓_𝑑0,𝜙₀=^𝐿𝑐−1_𝑘=0𝑐^∗(𝑘)𝑔^∗_𝑅−𝑡−𝑘𝑇_𝑐+𝜏₀+𝐿_𝑐𝑇_𝑐⋅𝑒^−𝑗𝜙₀𝑒^−𝑗2𝜋𝑓_{𝑑0(𝑡)𝑡}=𝑞^∗_𝑅−𝑡+𝜏₀+𝐿_𝑐𝑇_𝑐𝑒^−𝑗𝜙₀𝑒^−𝑗2𝜋𝑓_{𝑑0(𝑡)𝑡}, (11) where 𝜏₀,𝑓_𝑑0,𝜙₀ 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 𝑦𝑡;𝜏̆₀,𝑓̆_𝑑0,𝜙̆₀=𝑟₀̃𝑡;𝜏₀,𝑓_𝑑0,𝜙₀∗ℎ_MF𝑡;𝜏̆₀,𝑓̆_𝑑0,𝜙̆₀. (12)

Notice that, in the matched filter, we have substituted the estimates 𝜏₀, 𝑓_𝑑0, and 𝜙₀ for trial values obtained from previous (in time) estimates of these parameters which we have defined as 𝜏̆₀, 𝑓̆_𝑑0, and 𝜙̆₀, 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 𝑇_int and then release an output, we can write the discrete version of the signal as 𝑦_𝑛=^{⌊𝑇int/𝑇𝑠−1⌋}_𝑠=0𝑦𝑛𝑇_int−𝑠𝑇_𝑠;𝜏̆₀,𝑓̆_𝑑0,𝜙̆₀, (13) where 𝑇_𝑠 is the sampling period, 𝑇_int is the integration time (usually, 𝑇_int=𝑇_𝑐𝐿_𝑐) and ⌊⋅⌋ stands for the nearest integer towards zero.

Equation (13) can be expressed more conveniently by solving the convolution in (12), which yields [14] 𝑦_𝑛,𝐼=𝑎₀2𝐾sin𝜋Δ𝑓𝑇_int𝜋Δ𝑓𝑇_int𝑑[𝑛]_{𝑇𝑏/𝑇int}𝑅_𝑝𝑞Δ𝜏₀⋅cos𝜋Δ𝑓𝑇_int+Δ𝜙+𝑣_𝐼(𝑛), (14) where we defined Δ𝑓=𝑓_𝑑0−𝑓̆_𝑑0, Δ𝜙=𝜙₀−𝜙̆₀ and Δ𝜏₀=𝜏₀−𝜏̆₀ (i.e., the estimation errors), ⌊⋅⌋ stands for the nearest integer toward zero, and [𝑛]_{𝑇𝑏/𝑇int} means the integer part of 𝑛𝑇_int/𝑇_𝑏, being 𝑇_𝑏 the navigation bit period, and 𝑅_𝑝𝑞(𝜉)=1𝑇_PRN^𝑇PRN_𝜉^+𝜉𝑝(𝑡)𝑞^∗(𝑡−𝜉)𝑑𝑡 (15) is the correlation function. An equivalent derivation for the 𝑄 arm leads to 𝑦_𝑛,𝑄=𝑎₀2𝐾sin𝜋Δ𝑓𝑇_int𝜋Δ𝑓𝑇_int𝑑[𝑛]_{𝑇𝑏/𝑇int}𝑅_𝑝𝑞Δ𝜏₀⋅sin𝜋Δ𝑓𝑇_int+Δ𝜙+𝑣_𝑄(𝑛). (16)

Terms Δ𝑓, Δ𝜙, and Δ𝜏₀ should be regarded as the average local phase error over the integration interval, that is, Δ𝜙=Δ𝜙+2𝜋Δ𝑓(𝑇_int/2), 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 𝑇Δ𝜙=Δ𝜙+2𝜋Δ𝑓^int2𝑇+2𝜋²6Δ𝑓.̇ (17) 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 𝐾=(𝑇_int/𝑇_𝑠) 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 𝐾+1 samples instead of 𝐾. This effect can be considered negligible for the analysis presented in this paper.

In the case of 𝑀_𝑖>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 𝑅_𝑝𝑞 (i.e., 𝑅_𝑝𝑞(Δ𝜏_0,𝑛+𝛿), 𝑅_𝑝𝑞(Δ𝜏_0,𝑛−𝛿), 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.

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].

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 𝐲_𝑛 used in the simulations presented in Section 5.

Figure 3 

Configuration screen of the FCM blockset.

3. Particle Filtering

Bayesian filtering involves the recursive estimation of states 𝐱_𝑛∈ℝ^𝑛𝑥 given measurements 𝐲_𝑛∈ℝ^𝑛𝑦 at time 𝑛 based on all available measurements, 𝐲_1∶𝑛={𝐲₁,…,𝐲_𝑛}. To that aim, we are interested in the filtering distribution 𝑝(𝐱_𝑛|𝐲_1∶𝑛), which can be recursively expressed as 𝑝𝐱_𝑛∣𝐲_1∶𝑛=𝑝𝐲_𝑛∣𝐱_𝑛𝑝𝐱_𝑛∣𝐱_𝑛−1𝑝𝐲_𝑛∣𝐲_1∶𝑛−1𝑝𝐱_𝑛−1∣𝐲_1∶𝑛−1, (18) with 𝑝(𝐲_𝑛|𝐱_𝑛) and 𝑝(𝐱_𝑛|𝐱_𝑛−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 𝑁 weighted random samples, forming the so-called set of particles {𝐱_𝑛^(𝑖),𝑤_𝑛^(𝑖)}^𝑁_𝑖=1. These random samples are drawn from the importance density distribution, 𝜋(⋅), 𝐱⁽ⁱ⁾_𝑛∼𝜋𝐱_𝑛∣𝐱^(𝑖)_0∶𝑛−1,𝐲_1∶𝑛, (19) and weighted according to the general formulation 𝑤_𝑛^(𝑖)∝𝑤^(𝑖)_𝑛−1𝑝𝐲_𝑛∣𝐱^(𝑖)_0∶𝑛,𝐲_1∶𝑛−1𝑝𝐱_𝑛^(𝑖)∣𝐱^(𝑖)_𝑛−1𝜋𝐱_𝑛^(𝑖)∣𝐱^(𝑖)_0∶𝑛−1,𝐲_1∶𝑛. (20)

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 high-probability regions [20, 21]. In this work, we consider a multinomial sampling scheme for the resampling step.

Require:  {𝐱(𝑖)𝑛−1,𝑤(𝑖)𝑛−1}𝑁𝑖=1 and 𝐲𝑛 Ensure:  {𝐱𝑛(𝑖),𝑤𝑛(𝑖)}𝑁𝑖=1 and 𝐱𝑛   1:   for  𝑖=1 to 𝑁  do   2:    Generate 𝐱𝑛(𝑖)∼𝜋(𝐱𝑛|𝐱(𝑖)0∶𝑛−1,𝐲1∶𝑛)   3:    Calculate 𝑤𝑛(𝑖)=𝑤(𝑖)𝑛−1𝑝(𝐲𝑛|𝐱(𝑖)0∶𝑛,𝐲1∶𝑛−1)𝑝(𝐱𝑛(𝑖)|𝐱(𝑖)𝑛−1)𝜋(𝐱𝑛(𝑖)|𝐱(𝑖)0∶𝑛−1,𝐲1∶𝑛)   4:    end  for   5:    for  𝑖=1 to 𝑁  do   6:      Normalize weights: 𝑤𝑛(𝑖)=𝑤𝑛(𝑖)∑𝑁𝑗=1𝑤𝑛(𝑗)   7:    end  for   8:    MMSE state estimation:  𝐱𝑛=∑𝑁𝑖=1𝑤𝑛(𝑖)𝐱𝑛(𝑖)   9:   {𝐱𝑛(𝑖),1/𝑁}𝑁𝑖=1 = Resample({𝐱𝑛(𝑖),𝑤𝑛(𝑖)}𝑁𝑖=1)

Algorithm 1 

Standard particle filtering (SPF).

3.1. 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 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 𝜃=𝑔(𝐱) represent any unbiased estimator for 𝜃 and 𝑇(𝐱) be a sufficient statistic for 𝜃 under 𝑝(𝐱,𝜃). Then the conditional expectation 𝜃^RB=𝔼{𝑔(𝐱)|𝑇(𝐱)} 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 𝔼𝜃−𝜃²≥𝔼𝜃−𝜃^RB2. (21)

Final remarks on Rao-Blackwellization are worth mentioning. (i)Rao-Blackwellization 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 𝑀−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).

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 𝑦_𝑛,ℓ𝜶_𝑛,𝜙_𝑛,𝝉_𝑛=^𝑀−1_𝑚=0𝛼_𝑚,𝑛𝑒^{𝑗Δ𝜙𝑚,𝑛}𝑅_𝑛,ℓΔ𝜏_{ℓ,𝑚,𝑛}+𝑣_𝑛,ℓ, (22) accounting that Δ𝜏_{ℓ,𝑚,𝑛}=𝜏_𝑚,𝑛−𝜏_0,𝑛−1+𝛿_ℓ corresponds to the point where the ℓ-th early/late sample is evaluated. As usual, 𝑚=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 𝐲_𝑛=𝑦_𝑛,1(𝜶_𝑛,𝝉_𝑛)²,…,𝑦_𝑛,𝐿(𝜶_𝑛,𝝉_𝑛)^2𝑇, (23) 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: 𝐱_𝑛=𝜏_0,𝑛,…,𝜏_{𝑀−1,𝑛𝝉𝑛},𝜏̇_0,𝑛,…,_{𝑀−1,𝑛̇𝝉𝑛}𝜏̇,𝛼_0,𝑛,…,𝛼_{𝑀−1,𝑛𝜶𝑛}^𝑇, (24) 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 time-evolving delay of the signals.

One could adopt many alternatives to specify the time-evolving 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., 𝑚=0) are described by 𝜏_0,𝑛𝜏̇_0,𝑛=1𝑇_int𝜏01𝜏̇^0,𝑛−1_0,𝑛−1+𝐮^𝜏_0,𝑛, (25) where 𝑇_int is the integration period and the process noise is an uncorrelated zero-mean Gaussian random variable with diagonal entries 𝜎²_0,𝜏 and 𝜎²_̇𝜏0,.

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 𝜏_𝑚,𝑛>𝜏_0,𝑛∀𝑚∈{1,…,𝑀−1}, (26) in outdoor propagation channels [6, 11, 36]. Taking (26) into account, we force this situation using the evolution 𝜏_𝑚,𝑛=𝜏_0,𝑛+𝜏_{𝑚,𝑛−1}+𝑢^𝜏_𝑚,𝑛,𝜏̇_𝑚,𝑛=𝜏̇_0,𝑛+𝑢^̇𝜏_𝑚,𝑛, (27) with 𝑢^𝜏_𝑚,𝑛 and 𝑢^̇𝜏_𝑚,𝑛 being zero-mean 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 𝜎²_0,𝜏, 𝜎²_𝑚,𝜏, 𝜎²_̇𝜏0,, 𝜎²_̇𝜏𝑚,, 𝜎²_0,𝛼, 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 𝐱^l_𝑛=𝜶_𝑛,𝐱^nl_𝑛=𝝉^𝑇_𝑛,𝝉̇^𝑇_𝑛^𝑇. (28)

By the chain rule of probability, linear states can be analytically marginalized out from 𝑝(𝐱_𝑛∣𝐲_1∶𝑛): 𝑝𝐱^l_𝑛,𝐱^nl_0∶𝑛∣𝐲_1∶𝑛=𝑝𝐱^l_𝑛∣𝐱^nl_0∶𝑛,𝐲_1∶𝑛𝑝𝐱^nl_0∶𝑛∣𝐲_1∶𝑛 (29) and, taking into consideration that 𝐱^l_𝑛 generates a linear Gaussian state-space, 𝑝(𝐱^l_𝑛∣𝐱^nl_0∶𝑛,𝐲_1∶𝑛) can be updated analytically via a KF conditional on 𝐱^nl_0∶𝑛 and only the non-linear part of 𝐱_𝑛 needs to be estimated with a nonlinear filter. In the proposed scheme, an SPF is run to characterize 𝑝(𝐱^nl_0∶𝑛∣𝐲_1∶𝑛) and a KF is executed to obtain 𝑝(𝐱^l_𝑛∣𝐱^nl_0∶𝑛,𝐲_1∶𝑛).

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 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 2MHz. Estimates of time delay were performed at a rate of 50Hz, which corresponds to an integration time of 20 milliseconds, assuming bit synchronization. The carrier-to-noise density ratio (C/N₀) of the simulated satellite was 38dB-Hz. The dynamics of the scenario were due to the relative motion of the satellite-receiver, which is completely simulated by the GRANADA FCM blockset, and the receiver performed a pedestrian-like trajectory at 1m/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-the-art in GNSS receivers) with an equivalent noise bandwidth of 1Hz. 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 10Hz. The initial time-delay ambiguity at which the filter was initialized was drawn from 𝒩(𝜏_0,0,𝑇_chip/2), with 𝑇_chip 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 𝐿=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 𝐿=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 (𝑡∈[10,20]) it has been set to 0.25 chips and in the second interval (𝑡∈[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 𝑁=1000 particles. Namely, (i) in Figure 5 we used standard deviations 𝜎_0,𝜏=.03/𝑐,  𝜎_̇𝜏0,=0.03/𝑐, 𝜎_𝑚,𝜏=100/𝑐,  𝜎_̇𝜏𝑚,=0.03/𝑐,𝜎_0,𝛼=0.0001, and 𝜎_𝑚,𝛼=0.01; (ii) in Figure 6 we used 𝜎_0,𝜏=30/𝑐, 𝜎_̇𝜏0,=0.3/𝑐, 𝜎_𝑚,𝜏=30/𝑐, 𝜎_̇𝜏𝑚,=0.3/𝑐,𝜎_0,𝛼=0.0001, and 𝜎_𝑚,𝛼=0.0001; and finally (iii) in Figure 7 we used 𝜎_0,𝜏=3/𝑐, 𝜎_̇𝜏0,=0.3/𝑐, 𝜎_𝑚,𝜏=30/𝑐, 𝜎_̇𝜏𝑚,=0.3/𝑐,𝜎_0,𝛼=0.0001, and 𝜎_𝑚,𝛼=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.

Table 2 

Comparison of pseudorange error metrics in meters.

Figure 4 

Evolution of pseudorange error in DLL/PLL scheme under severe multipath propagation.

Figure 5 

Evolution of pseudorange error in MEPF scheme under severe multipath propagation, setup number 1.

Figure 6 

Evolution of pseudorange error in MEPF scheme under severe multipath propagation, setup number 2.

Figure 7 

Evolution of pseudorange error in MEPF scheme under severe multipath propagation, setup number 3.

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₀ of 38dB-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.48m and 4.82m for DLL and MEPF, respectively. For the MEPF we used 𝑀=2 paths, 𝑁=1000 particles, and 𝜎_0,𝜏=.03/𝑐, 𝜎_̇𝜏0,=0.03/𝑐, 𝜎_𝑚,𝜏=3/𝑐, 𝜎_̇𝜏𝑚,=0.03/𝑐,  𝜎_0,𝛼=0.0001, and 𝜎_𝑚,𝛼=0.001.

Figure 8 

Pseudorange error in DLL and MEPF schemes in an urban scenario at an elevation of 55∘, as modeled by Jahn’s model.

6. 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.

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 well-known result in statistical signal processing that both signals are sufficient statistics, and thus one is able to derive an estimator of 𝝊=(𝜏₀,𝑓_𝑑0,𝜙₀)^𝑇 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, 𝑇_int, since the signal-to-noise 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 𝐫_𝑛=𝑟₀𝑛𝑇_int⋮𝑟₀𝑛𝑇_int+(𝐾−1)𝑇_𝑠, (A.1) using the same notation conventions used along the document. Then, the likelihood can be decomposed as the independent contribution of each snapshot log𝑝𝐫_𝑛∣𝝊_𝑛=^𝐾_𝑘=1𝑝𝑟₀𝑛𝑇_int+(𝑘−1_{𝑠𝑟0,𝑛(𝑘)})𝑇∣𝝊_𝑛, (A.2) and assuming Gaussianity for the noise term, we could identify that log𝑝𝐫_𝑛∣𝝊_𝑛∝−^𝐾_𝑘=1𝑟_0,𝑛(𝑘)−ℎ_𝑛⁽¹⁾𝑘;𝝊_𝑛², (A.3) where ℎ_𝑛⁽¹⁾ stands for the precorrelation signal model, which was defined earlier as ℎ_𝑛⁽¹⁾𝑡;𝝊_𝑛=𝑎₀(𝑡)𝑒^{𝑗(2𝜋𝑓𝑑0𝑡+𝜙0)}𝑑_𝑛𝑝𝑡−𝜏₀. (A.4)

Further manipulation of the loglikelihood yields to log𝑝𝐫_𝑛∣𝝊_𝑛∝^𝐾_𝑘=1𝑟_0,𝑛(𝑘)𝑝^∗𝑡−𝝊_𝑛²−𝑟_0,𝑛(𝑘)²∝𝑦_𝑛𝝊−𝝊_𝑛², (A.5) 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 signal-to-noise 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 𝝊̆=(𝜏̆₀,𝑓̆_𝑑0,𝜙̆₀)^𝑇. (ii)Few samples are sufficient to infer estimates of 𝝊. After correlation an integration over a certain interval is already done, 𝑇_int, and therefore the signal-to-noise 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 log-likelihood under the Gaussian assumption is log𝑝𝑦_𝑛∣𝝊_𝑛∝𝑦_𝑛(Δ𝝊)−ℎ_𝑛⁽²⁾𝝊_𝑛−𝝊̆²∝ℜ𝑦_𝑛(Δ𝝊)ℎ_𝑛⁽²⁾𝝊_𝑛−𝝊̆^∗, (A.6) with ℎ_𝑛⁽²⁾ being the postcorrelation signal model ℎ_𝑛⁽²⁾()=𝑎Δ𝝊⁰2𝐾sin𝜋Δ𝑓𝑇_int𝜋Δ𝑓𝑇_int𝑑[𝑛]_{𝑇𝑏/𝑇int}⋅𝑅_𝑝𝑞Δ𝜏₀cos𝜋Δ𝑓𝑇_int+Δ𝜙, (A.7) and 𝝊_𝑛 the unknown parameter we want to estimate at 𝑛.

If we set 𝝊_𝑛=𝝊̆, we can identify that log𝑝𝑦_𝑛∣𝝊_𝑛∝ℜ𝑦_𝑛𝝊−𝝊_𝑛ℎ_𝑛⁽²⁾(0)^∗∝𝑦_𝑛𝝊−𝝊_𝑛². (A.8)

From the latter mathematical derivations, we can conclude an important result: log𝑝𝐫_𝑛∣𝝊_𝑛∝log𝑝𝑦_𝑛∣𝝊_𝑛 (A.9) for a given integration interval 𝑇_int 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 post-correlation samples. Another important conclusions is that since in the pre-correlation approach we also need to integrate in order to increase the signal-to-noise ratio, effects happening faster than 𝑇_int will not be captured by the estimation algorithm. The same happens in the post-correlation 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., pre-or postcorrelated samples).

Acknowledgment

P. Closas and C. Fernández-Prades were supported by the European Commission under COST Action IC0803 (RFCSET).