Bayesian Change-Points Estimation Applied to GPS Signal Tracking

Abstract

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 𝑆/𝑁₀. 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: 𝐼_𝑠(𝑡)=𝐴_𝑠𝐶_𝑠𝑡−𝜏_𝑠𝐶_𝑠𝑡−𝜏̂_𝑠×cos2𝜋𝑓_𝑠−𝑓_𝑠̂‌𝑡+𝜙_𝑠−𝜙_𝑠+𝜂^𝐼_𝑠(𝑡),𝑄_𝑠(𝑡)=𝐴_𝑠𝐶_𝑠𝑡−𝜏_𝑠𝐶_𝑠𝑡−𝜏̂_𝑠×sin2𝜋𝑓_𝑠−𝑓_𝑠̂‌𝑡+𝜙_𝑠−𝜙_𝑠+𝜂^𝑄_𝑠(𝑡).(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 (𝑆/𝑁₀). 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𝑆𝑁₀𝑑𝑡𝑅𝜏𝑒_𝑖−Δ×sinc𝑓𝑒_𝑖+𝛼_𝑖𝑇𝑐_𝑖2cos𝜙𝑒_𝑖+𝜂^𝐼_𝑖,𝑄_𝑖(Δ)=1√𝑁^𝑁_𝑘=1𝑄_𝑠(𝑡)=2𝑆𝑁₀𝑑𝑡𝑅𝜏𝑒_𝑖−Δ×sinc𝑓𝑒_𝑖+𝛼_𝑖𝑇𝑐_𝑖2sin𝜙𝑒_𝑖+𝜂^𝑄_𝑖.(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 𝑆/𝑁₀ 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𝑆𝑁₀𝑑𝑡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𝐷_𝑖=22𝑆𝑁₀𝑑𝑡sinc𝑓𝑒_𝑖+𝛼_𝑖𝑇𝑐_𝑖2×cos𝜙𝑒_𝑖|||𝜏𝑒_𝑖𝑇𝑏+12|||−|||𝜏𝑒_𝑖𝑇𝑏−12|||+𝜂^𝐷_𝑖,𝐷_𝑖=22𝑆𝑁₀𝑑𝑡sinc𝑓𝑒_𝑖+𝛼_𝑖𝑇𝑐_𝑖2cos𝜙𝑒_𝑖𝜏𝑒_𝑖𝑇𝑏+𝜂^𝐷_𝑖.(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.

(a)

(b)

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+𝑇𝑐²2𝑎,𝑣_𝑖=𝑣_𝑖−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 𝑡₁, as shown on Figure 1, is the instant such as 𝜏𝑒_𝑡1=𝑇𝑏/2. Then we have the following expression for 𝑎 in the first segment: 𝑎=𝐶𝑇𝑏−2𝑡₁−𝑡₀𝑇𝑐𝑣₀𝑇𝑐²𝑡₁−𝑡₀².(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 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 ℎ𝑦∣𝑟,𝑠;𝜎²=^𝐾𝑟_𝑘=1^𝑡𝑘_{𝑖=𝑡𝑘−1+1}1√2𝜋𝜎²exp⎛⎜⎝−𝑦_𝑖−𝑠_𝑘𝐴_𝑖,𝑘²2𝜎²⎞⎟⎠,(17) with 𝜎²=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 𝑃𝑠∣𝑟;𝜇,𝑉=^𝐾𝑟_𝑘=1√2𝜋𝑉𝑛_𝑘exp⎛⎜⎝−𝑛_𝑘𝑠_𝑘−𝜇²2𝑉⎞⎟⎠.(18) We show in Appendix A thatℎ𝑦∣𝑟;𝜇,𝑉=14𝜋^{𝑛/2𝐾𝑟}_𝑘=1𝑉_𝑘𝑛_𝑘𝑉×exp⎛⎜⎝−12⎛⎜⎝12^𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖+𝑛_𝑘𝜇²𝑉−𝜇²_𝑘𝑉_𝑘⎞⎟⎠⎞⎟⎠,(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𝑎_𝑘,𝑣_𝑖̃‌=𝑣̃‌_𝑖−1+𝑇𝑐𝑎_𝑘,(21) with for each segment 𝑘:𝑎_𝑘=2−2𝑇𝑐𝑛_𝑘𝑣̃‌_𝑡𝑘−1𝑇𝑐²𝑛²_𝑘,𝐴_{𝑡𝑘−1+1,𝑘}=𝑣̃‌_𝑡𝑘−1+𝑇𝑐𝑎_𝑘.(22) The model is initialized by 𝑣̃‌_𝑡0=𝑣₀𝜇,(23) where 𝑣₀ 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/𝑟𝑖=01−𝜆_𝑖,(25) where 𝜆_𝑖 is the probability to have a change at instant 𝑖. Then we have the following penalty term: ln𝜋𝑟=^𝑛_𝑖=1𝑟_𝑖ln𝜆_𝑖1−𝜆_𝑖+^𝑛_𝑖=1ln1−𝜆_𝑖.(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: 𝜆(𝑡)=12𝜋𝜎²_𝑑^𝐾𝑟_𝑘=1exp⎛⎜⎝−𝑡−𝑡𝑑_𝑘²2𝜎²_𝑑⎞⎟⎠,(27) where 𝑡𝑑_𝑘, the change instant for the segment 𝑘, is processed with the Doppler frequency 𝑓_𝑑 as follow: 𝑡𝑑_𝑘=𝑡𝑑_𝑘−1+𝑇𝑐𝑓_𝐿22𝑓_𝑑,𝑡𝑑₀=0.(28)𝜎²_𝑑 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,𝑡𝑑₁=30,𝑡𝑑₂=70,𝑡𝑑₃=90).

(a)

(b)

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^𝐾𝑟_𝑘=1⎧⎨⎩12^𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖−𝜇²_𝑘𝑉_𝑘−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 𝑇₀=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𝑎_𝑘,𝑣_𝑖=𝑣_𝑖−1+𝑇𝑐𝑎_𝑘,(31) with for each segment 𝑘: 𝑎_𝑘=𝐶𝑇𝑏−2𝑇𝑐𝑛_𝑘𝑣_𝑡𝑘−1𝑇𝑐²𝑛²_𝑘.(32) The initial relative speed satellite receiver 𝑣₀ 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×10⁻²¹𝑠 and 𝑆_𝑔=5.9×10⁻²⁰𝑠⁻¹.

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 𝑆/𝑁₀.

We show in Figure 4 positions computed every millisecond in the static noiseless case, for a signal-to-noise ratio of 𝑆/𝑁₀=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.

Appendices

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: 𝑃⎛⎜⎝𝑠𝑦,𝑟⎞⎟⎠ℎ𝑦𝑟=ℎ𝑦𝑟,𝑠;𝜎²𝑃𝑠𝑟;𝜇,𝑉.(A.2) Then we have Γ=ℎ(𝑦/𝑟,𝑠;𝜎²)𝑃(𝑠/𝑟;𝜇,𝑉) defined with (17) and (18) as follows: Γ=^𝑆𝑟_𝑘=112𝜋𝜎²^𝑛𝑘/2exp⎛⎜⎝−^𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦_𝑖−𝑠_𝑘𝐴_𝑖,𝑘²2𝜎²⎞⎟⎠×𝑛_𝑘2𝜋𝑉^1/2exp⎛⎜⎝−𝑛_𝑘𝑠_𝑘−𝜇²2𝑉⎞⎟⎠=^𝑆𝑟_𝑘=112𝜋𝜎²^𝑛𝑘/2𝑛_𝑘2𝜋𝑉^1/2×exp−12⎛⎜⎝1𝜎^2𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖+𝑛_𝑘𝑉𝜇²⎞⎟⎠×exp−12⎛⎜⎝𝑠²_𝑘⎛⎜⎝1𝜎^2𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝐴²_𝑖,𝑘+𝑛_𝑘𝑉⎞⎟⎠−2𝑠_𝑘⎛⎜⎝1𝜎^2𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦_𝑖𝐴_𝑖,𝑘+𝜇𝑛_𝑘𝑉⎞⎟⎠⎞⎟⎠.(A.3) Let:𝑇_𝑘=^𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝐴²_𝑖,𝑘,𝑉_𝑘=1𝑇_𝑘/𝜎²+𝑛_𝑘/𝑉,𝜇_𝑘=𝑉_𝑘⎛⎜⎝𝜇𝑛_𝑘𝑉+1𝜎^2𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦_𝑖𝐴_𝑖,𝑘⎞⎟⎠.(A.4) Then we haveΓ=^𝑆𝑟_𝑘=112𝜋𝑉_𝑘^1/2exp−12𝑉_𝑘𝑠_𝑘−𝜇_𝑘²×12𝜋𝜎²^𝑛𝑘/2𝑛_𝑘𝑉_𝑘𝑉^1/2×exp−12⎛⎜⎝1𝜎^2𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖+𝑛_𝑘𝑉𝜇²−𝜇²_𝑘𝑉_𝑘⎞⎟⎠,(A.5) then we can deduce the following distributions:𝑃⎛⎜⎝𝑠𝑦,𝑟;𝜇,𝑉,𝜎²⎞⎟⎠=^𝑆𝑟_𝑘=112𝜋𝑉_𝑘^1/2exp−12𝑉_𝑘𝑠_𝑘−𝜇_𝑘²,ℎ𝑦𝑟;𝜇,𝑉,𝜎²=12𝜋𝜎²^{𝑛/2𝑆𝑟}_𝑘=1𝑉𝑛_𝑘𝑉_𝑘^−1/2×exp−12⎛⎜⎝1𝜎^2𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖+𝑛_𝑘𝑉𝜇²−𝜇²_𝑘𝑉_𝑘⎞⎟⎠.(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𝜋)+^𝐾𝑟_𝑘=1−12log𝑉_𝑘+log𝑛_𝑘𝑉×^𝐾𝑟_𝑘=112⎧⎨⎩12^𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖−𝜇²_𝑘𝑉_𝑘+𝑛_𝑘𝜇²𝑉⎫⎬⎭.(B.2) In the same way −log𝜋(𝑟;𝜎_𝑑) can be rewritten as follows:−log𝜋𝑟;𝜎_𝑑=−^𝐾𝑟_𝑘=1log𝜆_𝑡𝑘−^𝐾𝑟_𝑘=1^𝑡𝑘−1_{𝑖=𝑡𝑘−1+1}1−𝜆_𝑖−^𝐾𝑟_𝑘=1log𝜆_𝑡𝑘+^𝐾𝑟_𝑘=1log1−𝜆_𝑡𝑘−^𝐾𝑟_𝑘=1^𝑡𝑘_{𝑖=𝑡𝑘−1+1}1−𝜆_𝑖−^𝐾𝑟_𝑘=1log𝜆_𝑡𝑘1−𝜆_𝑡𝑘−^𝐾𝑟_𝑘=1^𝑡𝑘_{𝑖=𝑡𝑘−1+1}1−𝜆_𝑖.(B.3) Finally the expression of the penalized contrast function is 𝑉_𝑦𝑟;𝜇,𝑉−log𝜋𝑟;𝜎_𝑑=^𝐾𝑟_𝑘=112⎧⎨⎩12^𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖−𝜇²_𝑘𝑉_𝑘−log𝑉_𝑘−2∗log𝜆_𝑡𝑘1−𝜆_𝑡𝑘⎫⎬⎭×𝑛2log(4𝜋)+^𝐾𝑟_𝑘=1log𝑛_𝑘𝑉−^𝐾𝑟_𝑘=1^𝑡𝑘_{𝑖=𝑡𝑘−1+1}1−𝜆_𝑖.(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 by𝑟̂=Argmin_{𝑟∈{0,1}}12^𝐾𝑟_𝑘=1⎧⎨⎩12^𝑡𝑘_{𝑖=𝑡𝑘−1+1}𝑦²_𝑖−𝜇²_𝑘𝑉_𝑘−log𝑉_𝑘⎫⎬⎭−^𝐾𝑟_𝑘=1log𝜆_𝑡𝑘1−𝜆_𝑡𝑘.(B.5)