Table of Contents Author Guidelines Submit a Manuscript
ISRN Signal Processing
Volume 2011 (2011), Article ID 148242, 11 pages
Research Article

Bayesian Change-Points Estimation Applied to GPS Signal Tracking

Laboratoire d'Informatique, Signaux et Images de la Côte d'Opale (LISIC), Université Lille Nord de France, Université du Littoral Côte d'Opale (ULCO) 50, rue Ferdinand Buisson BP 699, 62228 CALAIS Cedex, France

Received 24 February 2011; Accepted 20 April 2011

Academic Editors: G. Camps-Valls, C. S. Lin, and K. Sivakumar

Copyright © 2011 G. Stienne et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A hierarchical Bayesian model is applied to off-line segmentation of the GPS signal discriminator. The purpose of this work is to estimate the code delay of the receiving GPS CDMA code in order to retime the local receiver code and to estimate the pseudorange satellite receiver. The goal of our approach is to obtain a high-rate accurate positioning in the dynamic navigation case. We show that the behaviour of the coherent discriminator of a GPS pilot channel can be modelized by a piecewise stationary process. In our approach the discriminator behaviour in each stationary segment is approximated by a constant acceleration model, and the code delay at each end of the segments is known. The interest of this approach is that we use the coherent values of the discriminator in each segment to estimate the change instants of the process and to get in this case an accurate estimation of the code delays. In this context, a simultaneous estimation of the change instants is considered. We define the a posteriori distribution which integrates in its expression the signal change instants and the parameters of its statistical model. The proposed model leads after marginalization to a penalized contrast function that we minimize to estimate the discriminator change instants. The interest of the proposed model is that we can integrate in our estimate prior information on the roughly known values of the signal-to-noise ratio and relative speed satellite receiver. The potential of the proposed method is shown on experimentations realized on synthetic and real data for millisecond receiver localization.

1. Introduction

A GPS receiver can provide a position every millisecond at frequencies 𝑓𝐿1 and 𝑓𝐿5 or every 20 ms at frequency 𝑓𝐿2 [1, 2]. In the static case, the GPS receiver is at a fixed position and its localization processed by the receiver is averaged in time. In this case better positioning is obtained for a longer integration time. In the dynamic case, integration times are smaller in order to take into account the receiver trajectory. So we will have in the navigation case an accuracy that will depend on the receiver speed. The accuracy will decrease with the integration time as the speed of the receiver increases.

To improve the navigation accuracy, the GPS receiver is coupled with dead reckoning sensors. These sensors provide information used to integrate GPS signal most of the time in an extended Kalman filter. Unfortunately these sensors are in general inaccurate and drift with time, so they cannot be used for a long time integration. They are principally used to provide a position when the GPS receiver fails [3].

Most of the published works about GPS signal integration concern the problem of localization for low signal-to-noise ratio. An application in this case is for example the navigation in cover environment like indoor localization. Long coherent integration is then realized on several dozens of milliseconds for which the navigation message is known, and the GPS signal is supposed to be stationary. In these works an extended Kalman filter is used for the on-line carrier and code tracking [4]. We also found in the published works Bayesian methods using Markov chain Monte Carlo (MCMC) techniques. The goal of these methods is to estimate code delay and frequency offset of the GPS signal in the case of interference and jamming protection [5].

The GPS signal discriminator is a measurement that allows to compute the time delay between the CDMA code of the received signal and the local CDMA code generated by the receiver. In the classical match filter approach, the code delay is fixed, and the code is retimed for each new value of discriminator. In this context a coherent integration of the GPS signal on several milliseconds is used to get less noisy discriminator values. Unfortunately the coherent integration length is limited by the navigation message and by the phase, frequency, and code delay of the received signal that changes in time. The new GPS civil signals 𝐿2𝐶 and 𝐿5 have a dataless pilot channel (the 𝐶𝐿 code at 𝑓𝐿2 frequency, the 𝑄5 code at 𝑓𝐿5 frequency). These new channels allow to perform longer coherent integration of the signal, but they are limited with the evolution of the received signal parameters that change in time.

In our approach we estimate the code delay with several coherent discriminator values obtained by a code tracking open loop. We show that for a pilot channel we can suppose the discriminator piecewise linear instead of the data channel for which the discriminator sign is reversed when the bits of the navigation message change. Furthermore, we show that the change instants of the piecewise evolution of the discriminator are associated with known values of the code delay. In this context, we propose to estimate the change instants in order to compute the code delays. The proposed estimate is off line. We considered samples of the GPS signal grabbed during several seconds, and the code delays (instants of change in the piecewise evolution of the discriminator) are estimated simultaneously on the whole signal. When the change instants in the discriminator are estimated, we use a constant acceleration dynamic model in order to compute the pseudorange in each segment.

We can find in the published works a great number of contributions about off-line segmentation [6, 7]. These segmentation techniques are in general model selection methods, which are based on the minimisation of a penalized criterion. We differentiate in the published works two types of criteria: the efficient and the consistent criteria. Efficient criterion is based on the unbiased estimate of a quadratic risk, and the purpose of the method is to choose a model that minimizes this risk [8]. Consistent criterions suppose an existing true model of minimal size. Then, the criterion allows to define this model with a probability of 1 when the number of data tends towards infinity [9]. In the standard Bayesian approach, applied to segmentation of piecewise stationary process, we construct an a posteriori distribution of the process where we introduce a random vector of change instants 𝑅. This vector takes the value one at the change instant and the value zero between two changes [10]. In this case, to estimate the change instants, we have to find the vector 𝑟, which is a realization of 𝑅, that minimizes a penalized contrast function. This penalized contrast function is expressed from the a posteriori distribution of the process. In general, the Bayesian estimates are deduced from the maximum a posteriori (MAP) or the minimum mean-square error (MMSE) and need MCMC techniques to estimate the change times [10, 11].

In this paper, we propose a penalized contrast function deduced from the MAP estimate of the posterior distribution of the change instants, conditionally to the discriminator process parameters. The posterior distribution is defined as the product of a likelihood function with a prior law. For this application the parameter of the likelihood distribution is the Signal to Noise Ratio 𝑆/𝑁0. In our approach we consider that we have an approximate value of this parameter provided by the acquisition process and we model this information with a Gaussian distribution. The final expression of the likelihood is obtained after marginalization of this parameter. In order to define the prior law, we suppose the configuration of changes to be a sequence of Bernoulli variables [10]. For this model, the probabilities to have a change are different for each instant. We set a high probability to the instants of change supposed to be the real instants of change. These probabilities are derived from the relative speed satellite receiver computed with the Doppler frequency provided by the FLL tracking process. Finally the penalized contrast function is deduced from the negative logarithm of the posterior distribution. The change instants are estimated by minimizing the penalized contrast function with a simulated annealing algorithm.

The paper is organized as follows: the second part is dedicated to the GPS signal modelling and to the problem positioning. The estimate of the discriminator change instants is introduced in the third part. The fourth part describes experimentations on synthetic and real data. Finally a conclusion ends the paper.

2. Problem Positioning

2.1. GPS Model

The GPS satellite signal consists of a code division multiple access (CDMA) 𝐶𝑠() modulated at frequency 𝑓𝑠. The signal is binary phase shift keying (BPSK) modulated (the L5 signal is Quadra phase shift keying modulated). The code and carrier are assumed to be both time and frequency shifted. A GPS receiver realizes the signal acquisition in a first stage. It estimates the carrier frequency ̂𝑓𝑠, the code delay ̂𝜏𝑠, and the phase delay 𝜙𝑠 of the received signal. At the end of this first stage, the receiver tracks the variations of these parameters as a function of time in a tracking module. The frequency and phase indeed change with the Doppler effect (the relative speed satellite receiver), and the code delay with the satellite-receiver pseudorange evolution. In order to compute the navigation solution, the receiver multiplies the received signal with in-phase and quad-phase replica of the GPS signal: 𝐼𝑠(𝑡)=𝐴𝑠𝐶𝑠𝑡𝜏𝑠𝐶𝑠𝑡̂𝜏𝑠×cos2𝜋𝑓𝑠̂𝑓𝑠𝑡+𝜙𝑠𝜙𝑠+𝜂𝐼𝑠(𝑡),𝑄𝑠(𝑡)=𝐴𝑠𝐶𝑠𝑡𝜏𝑠𝐶𝑠𝑡̂𝜏𝑠×sin2𝜋𝑓𝑠̂𝑓𝑠𝑡+𝜙𝑠𝜙𝑠+𝜂𝑄𝑠(𝑡).(1)

In these expressions 𝐴𝑠 is the signal amplitude; it is normalized to drive the noise variance of 𝜂𝐼𝑠(𝑡), 𝜂𝑄𝑠(𝑡) to 1. It is a function of the signal-to-noise ratio (𝑆/𝑁0). The noise is supposed to be white, Gaussian, and centered. We can notice, in these expressions, the absence of the navigation message because in this paper we consider the dataless pilot channel. We then define 𝜏𝑒=𝜏𝑠̂𝜏𝑠, 𝑓𝑒=𝑓𝑠̂𝑓𝑠, 𝜙𝑒=𝜙𝑠𝜙𝑠 the respective code delay, frequency, and phase errors. The acquisition module provides initial values of the GPS signal parameters, so these errors are supposed to be close to zero at the beginning of the tracking stage. The goal of the tracking module is to maintain these errors as small as possible. Parameters are estimated every millisecond, the period 𝑇𝑐 of the C/A CDMA code. The GPS signal is down converted (in order to be sampled) at an intermediate frequency by means of a downconverter driven by a local clock oscillator. In this context the clock noise disturbances are modelled as normal random walks. Then we have at the 𝑖th millisecond [12]𝑓𝑒𝑖=𝑓𝑒𝑖1+𝛼𝑖1𝑇𝑐𝑖1+𝑤𝑓𝑖,𝜙𝑒𝑖=𝜙𝑒𝑖1+2𝜋𝑓𝑒𝑖𝑇𝑐𝑖+𝑤𝜙𝑖.(2) In these expressions 𝑤𝑓𝑖 and 𝑤𝜙𝑖 are, respectively, the noise clock disturbances on the phase and frequency. 𝛼𝑖1 models the speed variation of the frequency linked to the movements of the receiver and clock noise disturbance. 𝑇𝑐𝑖, the code length defined at the receiver, depends on the carrier frequency, so we have 𝑇𝑐𝑖=𝑇𝑐𝑓𝑠̂𝑓𝑠+𝑓𝑒𝑖.(3) Finally the error expression of the code delay, which is a function of the relative speed satellite receiver at the 𝑖th millisecond, is given by 𝜏𝑒𝑖=𝜏𝑒𝑖1+𝑇𝑐𝑖𝑓𝑟̂𝑓𝑠+̃𝑓𝑒𝑖𝑓𝑟.(4) In this expression ̃𝑓𝑒𝑖=𝑖𝑗=1𝛼𝑗1𝑇𝑐𝑗1 is the error on the frequency associated with the movement of the receiver and 𝑓𝑟 is the intermediate frequency of the emitted satellite signal. Parameters of the received signal are estimated with the values of 𝐼(𝑡) and 𝑄(𝑡) integrated on one or several periods 𝑇𝑐 of the signal.

Let 𝑁 be the number of samples integrated at the 𝑖th period 𝑇𝑐, we define for each satellite 𝑠 [4], 𝐼𝑖(Δ)=1𝑁𝑁𝑘=1𝐼𝑠(𝑡)=2𝑆𝑁0𝑑𝑡𝑅𝜏𝑒𝑖Δ×sinc𝑓𝑒𝑖+𝛼𝑖𝑇𝑐𝑖2cos𝜙𝑒𝑖+𝜂𝐼𝑖,𝑄𝑖(Δ)=1𝑁𝑁𝑘=1𝑄𝑠(𝑡)=2𝑆𝑁0𝑑𝑡𝑅𝜏𝑒𝑖Δ×sinc𝑓𝑒𝑖+𝛼𝑖𝑇𝑐𝑖2sin𝜙𝑒𝑖+𝜂𝑄𝑖.(5)

In these expressions 𝑅() is the normalized correlation function of the C/A CDMA code. sin𝑐() is a function that models the error on the estimate of the carrier frequency. In general the error values of 𝜙𝑒𝑖 and 𝑓𝑒𝑖 are estimated in a filter with 𝐼𝑖(0) and 𝑄𝑖(0). The values of 𝜏𝑒𝑖 are estimated in a filter from the in-phase and quad phase observations 𝐼𝑒𝑖=𝐼𝑖𝑇𝑏2,𝐼𝑙𝑖=𝐼𝑖𝑇𝑏2, and 𝑄𝑒𝑖=𝑄𝑖(𝑇𝑏/2), 𝑄𝑙𝑖=𝑄𝑖(𝑇𝑏/2) [4]. 𝑇𝑏 is the length of one bit of the CDMA code. Finally the value of 𝑆/𝑁0 is obtained when 𝐼𝑖(0) is maximum, so for 𝑓𝑒𝑖=0, 𝜏𝑒𝑖=0 and 𝜙𝑒𝑖=0.

2.2. Problem Positioning

The discriminator “early-minus-late” is defined by the following expression for which we fix to 𝑇𝑏 the code delay between early and late code generated by the receiver:𝐷𝑖=𝐼𝑒𝑖𝐼𝑙𝑖=2𝑆𝑁0𝑑𝑡sinc𝑓𝑒𝑖+𝛼𝑖𝑇𝑐𝑖2×cos𝜙𝑒𝑖𝑅𝜏𝑒𝑖𝑇𝑏2𝑅𝜏𝑒𝑖+𝑇𝑏2+𝜂𝐷𝑖.(6) In this expression 𝜂𝐷𝑖 is a white Gaussian noise of variance 2. The function of correlation of the CDMA code is given by the following expression: 𝑅(𝜏)=1|||𝜏𝑇𝑏|||.(7) In the classical approach the code tracking is realized with a matched filter [13]. Every millisecond the discriminator value is processed, and the local code is delayed or advanced in order to minimize the shift with the received code. The goal of this process of retiming is to obtain discriminator values close to zero. In this context the local code is retiming every millisecond with a fix delay. The code is advanced or delayed as a function of the discriminator sign.

When the satellite moves away from the receiver, the theoretical variation of the discriminator early-minus-late strictly increases for |𝜏𝑒𝑖|<𝑇𝑏/2. The variation strictly decreases for 3𝑇𝑏/2>|𝜏𝑒𝑖|>𝑇𝑏/2 and stays equal to zero for |𝜏𝑒𝑖|>3𝑇𝑏/2. We show in Figure 1(a) the DOL (delay open loop) discriminator variations obtained when the code is not retiming. On Figure 1 the discriminator is shown for a frequency and phase error equal to zero. For |𝜏𝑒𝑖|𝑇𝑏/2, the discriminator behaviour is proportional to the code delay 𝜏𝑒𝑖 and is defined by𝐷𝑖=22𝑆𝑁0𝑑𝑡sinc𝑓𝑒𝑖+𝛼𝑖𝑇𝑐𝑖2×cos𝜙𝑒𝑖|||𝜏𝑒𝑖𝑇𝑏+12||||||𝜏𝑒𝑖𝑇𝑏12|||+𝜂𝐷𝑖,𝐷𝑖=22𝑆𝑁0𝑑𝑡sinc𝑓𝑒𝑖+𝛼𝑖𝑇𝑐𝑖2cos𝜙𝑒𝑖𝜏𝑒𝑖𝑇𝑏+𝜂𝐷𝑖.(8) When the code is retiming with a code delay of 𝑇𝑏/2 every pseudorange variation of 𝐶𝑇𝑏/2, the discriminator behaviour is approximately piecewise linear (𝐶 is the speed of light). We show in Figure 1(b) the piecewise linear variation of the discriminator in this case. In fact we model this evolution in each stationary piece with a constant acceleration dynamic model.

Figure 1: Delay open loop discriminator behavior.

We suppose here that the receiver trajectory does not modify the strict increasing (when the satellite moves away) or decreasing (when the satellite comes closer) of the discriminator. In our approach, the tracking process is off line, we then research simultaneously the set of instants associated with the maximum values of the discriminator. For each of these instants, the code is retimed, with a shift of 𝑇𝑏/2. When the code is retimed we use the following constant acceleration model to estimate the pseudorange satellite receiver 𝑑𝑖: 𝑑𝑖=𝑑𝑖1+𝑇𝑐𝑣𝑖1+𝑇𝑐22𝑎,𝑣𝑖=𝑣𝑖1+𝑇𝑐𝑎.(9) In this expression 𝑣𝑖 is the relative speed satellite-receiver. 𝑑𝑡0 and 𝑣𝑡0 are, respectively, the initial pseudo range and speed, processed in the acquisition step of the GPS signal. In this expression 𝑎 is the acceleration supposed to be constant in each stationary piece. The value of 𝑎 is defined in order to have 𝑑𝑡1=𝐶𝑇𝑏/2+𝑑𝑡0, where 𝑡1, as shown on Figure 1, is the instant such as 𝜏𝑒𝑡1=𝑇𝑏/2. Then we have the following expression for 𝑎 in the first segment: 𝑎=𝐶𝑇𝑏2𝑡1𝑡0𝑇𝑐𝑣0𝑇𝑐2𝑡1𝑡02.(10) The phase and frequency values that minimise 𝑓𝑒𝑖  and 𝜙𝑒𝑖  are estimated in each stationary piece with an extended Kalman filter [4].

With this modelization we need to know the change instants of the piecewise variation of the discriminator in order to estimate the pseudoranges. In practice we have discriminator values provided by the tracking open loop process. In the Bayesian approach we try to fit a model to the observations in order to estimate its parameters. In this context, as shown in Figure 1, the SNR and the change instants will be the parameters to estimate. We describe in the next section the Bayesian model proposed to estimate the change instants of the piecewise linear discriminator behaviour.

3. Bayesian Estimation

3.1. Bayesian Approach

The problem we address in this paper can be defined as the detection of changes in the statistical distribution of a one-dimensional process. Let 𝑦=(𝑦𝑖,1𝑖𝑛) be the observed discriminator values model as a series with conditional law 𝑃(𝑦). We suppose that the discriminator variations are piecewise stationary. Then we note 𝑟=(𝑟𝑖,1𝑖𝑛) the instants of change for which the code is retimed of 𝑇𝑏/2. 𝑟𝑖 takes the value 1 at the change instants and the value 0 between two changes. By convention there always is a change on the last sample, then we have 𝑟𝑛=1. The observed series is composed of 𝐾𝑟 segments, then we have 𝐾𝑟=𝑛𝑖=1𝑟𝑖. We notice (𝑛𝑘,𝑘0) the number of samples in the segment 𝑘. We notice 𝑡=(𝑡𝑘,𝑘0) the change instants in the observed series, associated with the change times configuration 𝑟 (𝐾𝑟 is the number of changes). We define the origin instant at 𝑡0=0 and the last change at 𝑡𝑁=𝑛.

Let the following expression for the discriminator at instant 𝑖 in the segment 𝑘: 𝑦𝑖=𝑠𝑘𝐴𝑖,𝑘.(11) In the simplified (11), 𝑠𝑘 is the maximum discriminator value, and it depends on the SNR. We note 𝑠=(𝑠𝑘,𝑘0) the set of maximum values of the discriminator. 𝐴𝑖,𝑘 is a strictly increasing function (or decreasing function according to the satellite trajectory) which depends on the values of 𝑓𝑒𝑖, 𝜙𝑒𝑖, 𝜏𝑒𝑖. The distribution that describes the noisy values of the discriminator can be defined by [14], 𝑃𝑟,𝑦,𝑠=𝑃𝑠𝑦,𝑟𝑃𝑦𝑟𝜋𝑟,(12) where 𝜋(𝑟) is the prior law of the change configuration 𝑟. Moreover we have 𝑃𝑠𝑦,𝑟𝑦𝑟=𝑦𝑟,𝑠𝑃𝑠𝑟,(13) where () is the likelihood function. Finally the maximum a posteriori estimate of the change configuration is defined by ̂𝑟=Argmax𝑟𝑃𝑟𝑦.(14) In our implementation we suppose that we have an approximate value of the SNR and a prior information on the instants of change. The approximate value of 𝑠 is defined as a random variable of mean 𝜇 and variance 𝑉. The mean value of the SNR 𝜇 is defined in the acquisition process and 𝑉 is the inaccuracy of this value. The prior information on the instants of change is a function of the relative speed satellite receiver (which is a function of the Doppler). 𝜎𝑑 represents the inaccuracy of this prior information. So, ̂𝑟=Argmin𝑟log𝑦𝑟;𝜇,𝑉𝜋𝑟;𝜎𝑑,(15) where 𝜇,𝑉 and 𝜎𝑑 are hyperparameters that, respectively, define the prior information on the SNR and on the change instants. When we take the negative logarithm of (14), we define a penalized empirical criterion, where 𝑟 (the number and position of changes) is estimated by minimizing the penalized contrast function defined by 𝑈𝑦𝑟=𝑉𝑦𝑟;𝜇,𝑉log𝜋𝑟;𝜎𝑑.(16) In this expression the first term measures the fidelity to the observation (contrast function), whereas the second term (penalty) is related to the number and the position of changes.

3.2. Contrast Function

The likelihood distribution, of the discriminator defined in (8) and (9), is given by 𝑦𝑟,𝑠;𝜎2=𝐾𝑟𝑘=1𝑡𝑘𝑖=𝑡𝑘1+112𝜋𝜎2exp𝑦𝑖𝑠𝑘𝐴𝑖,𝑘22𝜎2,(17) with 𝜎2=2. The stage of acquisition provides an inaccurate estimation 𝑠𝑘 of the SNR. Furthermore, the SNR varies with time. Then we suppose 𝑠𝑘 as normally distributed random variables with mean 𝜇 and variance 𝑉/𝑛𝑘 [14]. The value of 𝜇 is given by the acquisition module, and the value of 𝑉 is fixed by the user. This value will be chosen as high as 𝜇 is inaccurate. Then we have 𝑃𝑠𝑟;𝜇,𝑉=𝐾𝑟𝑘=12𝜋𝑉𝑛𝑘exp𝑛𝑘𝑠𝑘𝜇22𝑉.(18) We show in Appendix A that𝑦𝑟;𝜇,𝑉=14𝜋𝑛/2𝐾𝑟𝑘=1𝑉𝑘𝑛𝑘𝑉×exp1212𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖+𝑛𝑘𝜇2𝑉𝜇2𝑘𝑉𝑘,(19) with𝑇𝑘=𝑡𝑘𝑖=𝑡𝑘1+1𝐴𝑖,𝑘,𝑉𝑘=1𝑇𝑘/2+𝑛𝑘/𝑉,𝜇𝑘=𝑉𝑘𝜇𝑛𝑘𝑉+12𝑡𝑘𝑖=𝑡𝑘1+1𝑦𝑖𝐴𝑖,𝑘.(20) To define the values of 𝐴𝑖,𝑘, we assume the constant acceleration dynamic model in the stationary pieces defined by (9). This model describes the variations of the discriminator 𝑦𝑖 in the stationary pieces. Then we have the following equations for 𝐴𝑖,𝑘 in the segment 𝑘 and for 𝑖{𝑡𝑘1+1,𝑡𝑘}: 𝐴𝑖,𝑘=𝐴𝑖1,𝑘+𝑇𝑐̃𝑣𝑖1+𝑇2𝑐2𝑎𝑘,̃𝑣𝑖=̃𝑣𝑖1+𝑇𝑐𝑎𝑘,(21) with for each segment 𝑘:𝑎𝑘=22𝑇𝑐𝑛𝑘̃𝑣𝑡𝑘1𝑇𝑐2𝑛2𝑘,𝐴𝑡𝑘1+1,𝑘=̃𝑣𝑡𝑘1+𝑇𝑐𝑎𝑘.(22) The model is initialized by ̃𝑣𝑡0=𝑣0𝜇,(23) where 𝑣0 is the relative speed satellite receiver and 𝜇 the SNR provided by the acquisition process.

3.3. Penalty Definition

We define 𝑟 as a sequence of independent Bernoulli variables. If we suppose that we do not have any prior information on the positions of changes, the probability to have a sequence 𝑟 is given by [10] 𝜋𝑟=𝜆𝐾𝑟(1𝜆)(𝑛𝐾𝑟).(24) If we suppose that we have information on the possible location of changes, we associate to these locations a high probability to have a change. In this context we propose to define the prior probability 𝜋(𝑟) to have a change configuration 𝑟 as 𝜋𝑟=𝑛𝑖=1/𝑟𝑖=1𝜆𝑖𝑛𝑖=1/𝑟𝑖=01𝜆𝑖,(25) where 𝜆𝑖 is the probability to have a change at instant 𝑖. Then we have the following penalty term: ln𝜋𝑟=𝑛𝑖=1𝑟𝑖ln𝜆𝑖1𝜆𝑖+𝑛𝑖=1ln1𝜆𝑖.(26) In our application we approximately know the change instants thanks to the Doppler frequency given by the acquisition module and estimated in each stationary piece. Then we propose to define the probability to have a change at a given instant 𝑡 by the following normal distribution: 𝜆(𝑡)=12𝜋𝜎2𝑑𝐾𝑟𝑘=1exp𝑡𝑡𝑑𝑘22𝜎2𝑑,(27) where 𝑡𝑑𝑘, the change instant for the segment 𝑘, is processed with the Doppler frequency 𝑓𝑑 as follow: 𝑡𝑑𝑘=𝑡𝑑𝑘1+𝑇𝑐𝑓𝐿22𝑓𝑑,𝑡𝑑0=0.(28)𝜎2𝑑 defines the uncertainty on the change localization. This uncertainty will be as high as the dynamic trajectory of the receiver will change rapidly. However in practice the signal is sampled so the data are discrete. In this context we propose to use a discrete Gaussian kernel as an approximation of (27). Then we propose the following expression for the definition of 𝜆𝑖 at instant 𝑖: 𝜆𝑖=𝐾𝑟𝑘=1𝑒𝑥𝐼𝑖(𝑥)𝛿𝑖𝑡𝑑𝑘,(29) where 𝐼𝑖(𝑥) is the modified Bessel function, 𝑥 is a parameter of concentration, and the function of convolution. In this expression we define 𝐾𝑟 discrete Gaussian distributions centered on the instants 𝑡𝑑=(𝑡𝑑𝑘,𝑘0). We show in Figure 2 the evolution of the penalty and its associate probability. In this example we have 𝑡𝑑=(𝑡𝑑0=0,𝑡𝑑1=30,𝑡𝑑2=70,𝑡𝑑3=90).

Figure 2: (a) Probability evolution and (b) penalty evolution.

Finally we show in Appendix B that the penalized contrast function to minimize is given by ̂𝑟=Argmin𝑟{0,1}12𝐾𝑟𝑘=112𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖𝜇2𝑘𝑉𝑘log𝑉𝑘𝐾𝑟𝑘=1log𝜆𝑡𝑘1𝜆𝑡𝑘.(30) In the next subsection we describe the complete algorithm used to obtain the set of instants corresponding to the maximum value of the discriminator and how we proceed to obtain the pseudorange.

3.4. Minimization of the Contrast Function

To minimize this penalized contrast expression, we use a simulated annealing algorithm. It is known that searching for the real configuration of changes in all possible sequences is not computationally tractable. Nevertheless, by assuming statistical independence between the segments, the optimal segmenter can be solved with reduced computation via simulated annealing [10]. The simulated annealing algorithm is an iterative procedure that converges to the optimal ̂𝑟 solution (that minimizes (30)) with probability one [15]. For our application, the simulated annealing algorithm proceeds as follows. (1)An initial random sequence ̂𝑟 is drawn. Choose an initial value for the parameter 𝑇 (𝑇𝑗,𝑗=0), called temperature (in our application 𝑇0=500). 𝑗=0. (2)A new configuration ̃𝑟 is drawn from ̂𝑟 in one of these three cases: (i)add a change at a random position 𝑖 in ̂𝑟,(ii)remove a change at a random position 𝑖 in ̂𝑟,(iii)translate a change at a random position 𝑖 in ̂𝑟.(3)Process the penalized contrast function in three steps: (i)compute the 𝐴𝑖,𝑘 with the change configuration ̃𝑟, (ii)compute the 𝑦𝑖, 𝜙𝑒𝑖, and 𝑓𝑒𝑖 for 𝑖{1,,𝑛} with the change configuration ̃𝑟. We use in this case the architecture tracking process described in Figure 3, (iii)compute 𝑈𝑦(̃𝑟).(4)Let Δ𝑈𝑦=𝑈𝑦(̃𝑟)𝑈𝑦(̂𝑟). (5)Set ̂𝑟=̃𝑟 with probability one if Δ𝑈𝑦<0 and with probability exp{Δ𝑈𝑦/𝑇𝑗} elsewhere. Decrease the value of 𝑇𝑗+1=𝑓(𝑇𝑗) (in our application 𝑇𝑗+1=𝑇𝑗0.99). (6)𝑗=𝑗+1. Go to (2).

Figure 3: Architecture of the tracking process.

When ̂𝑟 has been estimated with the simulated annealing algorithm, we use the following equations to compute the pseudorange:𝑑𝑖=𝑑𝑖1+𝑇𝑐𝑣𝑖1+𝑇2𝑐2𝑎𝑘,𝑣𝑖=𝑣𝑖1+𝑇𝑐𝑎𝑘,(31) with for each segment 𝑘: 𝑎𝑘=𝐶𝑇𝑏2𝑇𝑐𝑛𝑘𝑣𝑡𝑘1𝑇𝑐2𝑛2𝑘.(32) The initial relative speed satellite receiver 𝑣0 is provided by the acquisition process.

4. Experimentation

In this experimentation we want to assess the proposed algorithm in two different cases. In the first case, the GPS signal is simulated in a real context. In the second case, the proposed method is assessed on a real 𝐿1 GPS signal. These experiments will be realized and described for a static receiver and for a circular dynamic trajectory. Positions are obtained, every millisecond, for each C/A code correlation measurement, and the signal is sampled with a 20 MHz frequency. For the real signal we extract in a first step the navigation message in order to compensate the discriminator changes associated with a variation in the sign of the navigation message. The proposed method is then applied in a second step to the real signal multiplied by the estimated navigation message. In this experimentation we compare the results obtained with the proposed method and the results obtained with a classical matched filter [13].

4.1. Experimentation on Synthetic Data

We simulate a real GPS constellation the third day of the GPS week 1291 at 12 hours, 00 minute, and 00 second. The receiver antenna localization is supposed to be on the roof of the laboratory. We report in Table 1 the different GPS simulation parameters. The developed algorithm is demonstrated using simulated GPS signals with C/A codes and using temperature-compensated crystal oscillator clocks. The values of the phase and frequency random-walk intensities are 𝑆𝑓=5×1021𝑠 and 𝑆𝑔=5.9×1020𝑠1.

Table 1: Simulation parameters.

𝑥 and 𝑉 are parameters that model the accuracy of the prior information. These parameters are tuned by the user in order to obtain the best results. The value of 𝑥 will increase with the relative speed satellite receiver in the dynamic case. The value of 𝑉 will increase when the signal to noise decreases because in this case the acquisition process provides inaccurate estimate of 𝑆/𝑁0.

We show in Figure 4 positions computed every millisecond in the static noiseless case, for a signal-to-noise ratio of 𝑆/𝑁0=50 dB.Hz. On these figures, we show positions processed with the proposed method in black and with a standard method in grey. The standard method is the classical match filter [13]. Positions calculated with the standard method are scattered around the real position. In the other hand, with the proposed method the calculated positions in black are closed to the real position.

Figure 4: Static positions computed for a 𝑆/𝑁0 of 50 dB·Hz.

We report in Table 2 the mean gap values and variance, according to North and East axis, between the real and estimate positions. We also describe the mean Euclidian distance between the real and estimated receiver positions. The values reported in this table are obtained for different signal-to-noise ratios with the standard and the proposed method.

Table 2: Error on estimated static positions.

For the dynamic case, we consider a circular trajectory of radius 10 m ended in 5 s. We show in Figure 5, the estimated positions obtained in the dynamic experimental case. We notice for the standard method that points are scattered around the real trajectory like in the static case. On the other hand, for the proposed method the trajectory is described by a succession of straight lines close to the real trajectory.

Figure 5: Estimate dynamic trajectory for a 𝑆/𝑁0 of 50 dB.Hz.

We report in Table 3 errors on computed positions in the dynamic case. Compared to the static case, the errors we obtain for the standard method are practically the same. However the errors obtained with the proposed method are smaller than those for the standard case but larger than those in the static case. Furthermore, we notice, like in the static case, a smaller variation of the precision, described in Table 3, for the proposed method when the signal-to-noise ratio decreases.

Table 3: Error on estimated dynamic positions.
4.2. Experimentation on Real Data

In this experimentation we want to show the feasibility of the proposed method. In the static case, the antenna is situated on the roof of the laboratory. For the dynamic case, the antenna is situated on the roof of a car. For this experimentation, we make a circular dynamic trajectory around a roundabout. We show in Figures 6 and 7 the positions, respectively, computed in the static and dynamic real cases during 7 s. For the static case, we show on Figure 6 the positions obtained with the proposed method in black and with the classical matched filter in grey. For the dynamic case, we show on Figure 7 the positions obtained with the proposed method in black, and we show in grey the vehicle trajectory obtained with an accurate inertial station coupled with a differential GPS from Novatel. The results obtained with real data show the feasibility of the proposed method. Finally, compared to the matched filter, the proposed method provides more accurate code delay estimation. Actually, the proposed method allows to integrate coherent discriminator values on longer time. This integration is performed with a constant acceleration model that is fit to the discriminator measurements. In this case the interpolation is also used to compute the pseudorange and so to smooth the measurements. The computed positions are in this case less noisy and close to the true trajectory in the dynamic case. Nevertheless, the proposed method is time consuming and needs high storage capacity for the signal samples.

Figure 6: Positions calculated in the static real case.
Figure 7: Positions calculated in the dynamic real case (shown on an aerial photo).

5. Conclusions and Future Works

In this paper, a Bayesian off-line estimation technique is applied to the implementation of an open loop code tracking of the GPS signal discriminator. The goal of this work is to estimate the localization of a dynamic receiver every millisecond. The estimation technique is based on a hierarchical statistical model that leads to the definition of a MAP estimate. This MAP estimate is a penalized contrast function. The contrast function is constructed from the distribution of the piecewise linear variation of the discriminator. The penalized function is deduced from the a priori probability to have a change. It is a function of time. Experimentation and test of this method on synthetic and real data are presented for a static case and for a dynamic case. We show that we obtain a better positioning accuracy with our method, compared to a standard matched filter technique.

The prospect of this work is about information fusion supplied by the multiband GPS signal that will be available in the modern GPS systems.


A. Likelihood of the Change Instants

In this appendix we want to derive the likelihood distribution of the change instants (𝑦𝑟;𝜇,𝑉). We have 𝑃𝑠𝑦,𝑟𝑦𝑟=𝑦𝑟,𝑠𝑃𝑠𝑟.(A.1) This expression can be rewritten as follows if we introduce prior information: 𝑃𝑠𝑦,𝑟𝑦𝑟=𝑦𝑟,𝑠;𝜎2𝑃𝑠𝑟;𝜇,𝑉.(A.2) Then we have Γ=(𝑦/𝑟,𝑠;𝜎2)𝑃(𝑠/𝑟;𝜇,𝑉) defined with (17) and (18) as follows: Γ=𝑆𝑟𝑘=112𝜋𝜎2𝑛𝑘/2exp𝑡𝑘𝑖=𝑡𝑘1+1𝑦𝑖𝑠𝑘𝐴𝑖,𝑘22𝜎2×𝑛𝑘2𝜋𝑉1/2exp𝑛𝑘𝑠𝑘𝜇22𝑉=𝑆𝑟𝑘=112𝜋𝜎2𝑛𝑘/2𝑛𝑘2𝜋𝑉1/2×exp121𝜎2𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖+𝑛𝑘𝑉𝜇2×exp12𝑠2𝑘1𝜎2𝑡𝑘𝑖=𝑡𝑘1+1𝐴2𝑖,𝑘+𝑛𝑘𝑉2𝑠𝑘1𝜎2𝑡𝑘𝑖=𝑡𝑘1+1𝑦𝑖𝐴𝑖,𝑘+𝜇𝑛𝑘𝑉.(A.3) Let:𝑇𝑘=𝑡𝑘𝑖=𝑡𝑘1+1𝐴2𝑖,𝑘,𝑉𝑘=1𝑇𝑘/𝜎2+𝑛𝑘/𝑉,𝜇𝑘=𝑉𝑘𝜇𝑛𝑘𝑉+1𝜎2𝑡𝑘𝑖=𝑡𝑘1+1𝑦𝑖𝐴𝑖,𝑘.(A.4) Then we haveΓ=𝑆𝑟𝑘=112𝜋𝑉𝑘1/2exp12𝑉𝑘𝑠𝑘𝜇𝑘2×12𝜋𝜎2𝑛𝑘/2𝑛𝑘𝑉𝑘𝑉1/2×exp121𝜎2𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖+𝑛𝑘𝑉𝜇2𝜇2𝑘𝑉𝑘,(A.5) then we can deduce the following distributions:𝑃𝑠𝑦,𝑟;𝜇,𝑉,𝜎2=𝑆𝑟𝑘=112𝜋𝑉𝑘1/2exp12𝑉𝑘𝑠𝑘𝜇𝑘2,𝑦𝑟;𝜇,𝑉,𝜎2=12𝜋𝜎2𝑛/2𝑆𝑟𝑘=1𝑉𝑛𝑘𝑉𝑘1/2×exp121𝜎2𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖+𝑛𝑘𝑉𝜇2𝜇2𝑘𝑉𝑘.(A.6)

B. Penalized Contrast Function

In this appendix we want to derive the penalized contrast function from the likelihood and prior distributions. We define a penalized empirical criterion, and 𝑟 (the number and position of changes) is estimated by minimizing the following penalized contrast function: ̂𝑟=Argmin𝑟{0,1}𝑉𝑦𝑟;𝜇,𝑉log𝜋𝑟;𝜎𝑑.(B.1) In this expression 𝑉𝑦(𝑟;𝜇,𝑉) is the negative log of the likelihood distribution defined by (19) and is given by𝑉𝑦𝑟;𝜇,𝑉=𝑛2log(4𝜋)+𝐾𝑟𝑘=112log𝑉𝑘+log𝑛𝑘𝑉×𝐾𝑟𝑘=11212𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖𝜇2𝑘𝑉𝑘+𝑛𝑘𝜇2𝑉.(B.2) In the same way log𝜋(𝑟;𝜎𝑑) can be rewritten as follows:log𝜋𝑟;𝜎𝑑=𝐾𝑟𝑘=1log𝜆𝑡𝑘𝐾𝑟𝑘=1𝑡𝑘1𝑖=𝑡𝑘1+11𝜆𝑖𝐾𝑟𝑘=1log𝜆𝑡𝑘+𝐾𝑟𝑘=1log1𝜆𝑡𝑘𝐾𝑟𝑘=1𝑡𝑘𝑖=𝑡𝑘1+11𝜆𝑖𝐾𝑟𝑘=1log𝜆𝑡𝑘1𝜆𝑡𝑘𝐾𝑟𝑘=1𝑡𝑘𝑖=𝑡𝑘1+11𝜆𝑖.(B.3) Finally the expression of the penalized contrast function is 𝑉𝑦𝑟;𝜇,𝑉log𝜋𝑟;𝜎𝑑=𝐾𝑟𝑘=11212𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖𝜇2𝑘𝑉𝑘log𝑉𝑘2log𝜆𝑡𝑘1𝜆𝑡𝑘×𝑛2log(4𝜋)+𝐾𝑟𝑘=1log𝑛𝑘𝑉𝐾𝑟𝑘=1𝑡𝑘𝑖=𝑡𝑘1+11𝜆𝑖.(B.4) In this expression the last three terms are independent of the change configuration 𝑟. In this context the penalized contrast function to minimize in order to estimate ̂𝑟 is given bŷ𝑟=Argmin𝑟{0,1}12𝐾𝑟𝑘=112𝑡𝑘𝑖=𝑡𝑘1+1𝑦2𝑖𝜇2𝑘𝑉𝑘log𝑉𝑘𝐾𝑟𝑘=1log𝜆𝑡𝑘1𝜆𝑡𝑘.(B.5)


  1. G. Gerten, “Protecting the global positioning system,” IEEE Aerospace and Electronic Systems Magazine, vol. 20, no. 11, pp. 3–8, 2005. View at Publisher · View at Google Scholar · View at Scopus
  2. J. Issler, L. Ries, J. Bourgeade, L. Lestarquit, and C. Macabiau, “Probabilistic approach of frequency diversity as interference mitigation means,” in Proceedings of the 17th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS '04), pp. 2136–2145, Long-Beach, Calif, USA, September 2004. View at Scopus
  3. M. G. Petovello, C. O'Driscoll, and G. Lachapell, “Weak signal carrier tracking using extended coherent integration with an ultra-tight GNSS/IMU receiver,” in Proceedings of the European Navigation Conference (ENC-GNSS '08), Toulouse, France, April 2008.
  4. N. I. Ziedan, GNSS Receivers for Weak Signals, Artech House, London, UK, 2006.
  5. I. Progri, M. Bromberg, and J. Wang, “Markov chain, Monte Carlo global search and integration for Bayesian, GPS, parameter estimation,” Journal of the Institute of Navigation, vol. 56, no. 3, pp. 195–204, 2009. View at Google Scholar
  6. S. Boutoille, S. Reboul, and M. Benjelloun, “A hybrid fusion system applied to off-line detection and change-points estimation,” Information Fusion, vol. 11, no. 4, pp. 325–337, 2010. View at Publisher · View at Google Scholar · View at Scopus
  7. S. Reboul and M. Benjelloun, “Joint segmentation of the wind speed and direction,” Signal Processing, vol. 86, no. 4, pp. 744–759, 2006. View at Publisher · View at Google Scholar · View at Scopus
  8. E. Lebarbier, “Detecting multiple change-points in the mean of Gaussian process by model selection,” Signal Processing, vol. 85, no. 4, pp. 717–736, 2005. View at Publisher · View at Google Scholar · View at Scopus
  9. Y. C. Yao, “Estimating the number of change-points via Schwarz criterion,” Statistics and Probability Letters, vol. 6, no. 3, pp. 181–189, 1988. View at Google Scholar · View at Scopus
  10. M. Lavielle, “Optimal segmentation of random processes,” IEEE Transactions on Signal Processing, vol. 46, no. 5, pp. 1365–1373, 1998. View at Google Scholar · View at Scopus
  11. E. Punskaya, C. Andrieu, A. Doucet, and W. J. Fitzgerald, “Bayesian curve fitting using MCMC with applications to signal segmentation,” IEEE Transactions on Signal Processing, vol. 50, no. 3, pp. 747–758, 2002. View at Publisher · View at Google Scholar · View at Scopus
  12. A. J. Van Dierendonck and J. B. McGraw, “Relationship between Allan variances and Kalman filter parameters,” in Proceedings of the 16th Annual Precise Time and Time Interval Application and Planning Meeting (PTTI '84), pp. 274–293, Greenbelt, Md, USA, 1984.
  13. E. Kaplan, Understanding GPS: Principles and Applications, Artech House, Norwood, Mass, USA, 1996.
  14. M. Lavielle and E. Lebarbier, “An application of MCMC methods for the multiple change-points problem,” Signal Processing, vol. 81, no. 1, pp. 39–53, 2001. View at Publisher · View at Google Scholar · View at Scopus
  15. D. Geman, Random Fields and Inverse Problems in Imaging, vol. 1427/1990 of Lecture Notes in Mathematics, New York, NY, USA, Springer, 1990. View at Publisher · View at Google Scholar