Bayesian Change-Points Estimation Applied to GPS Signal Tracking

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.


Introduction
A GPS receiver can provide a position every millisecond at frequencies f L1 and f L5 or every 20 ms at frequency f L2 [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-tonoise 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 L2C and L5 have a dataless pilot channel (the CL code at f L2 frequency, the Q5 code at f L5 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 R.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 r, which is a realization of R, 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 S/N 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.

Problem Positioning
2.1.GPS Model.The GPS satellite signal consists of a code division multiple access (CDMA) C s (• • • ) modulated at frequency f s .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 f s , the code delay τ s , and the phase delay φ s 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: In these expressions A s is the signal amplitude; it is normalized to drive the noise variance of η I s (t), η Q s (t) to 1.It is a function of the signal-to-noise ratio (S/N 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 τe = τ s − τ s , f e = f s − f s , φe = φ s − φ s 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 Tc 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 ith millisecond [12] (2) In these expressions w f i and wφ i are, respectively, the noise clock disturbances on the phase and frequency.α i−1 models the speed variation of the frequency linked to the movements of the receiver and clock noise disturbance.Tc i , the code length defined at the receiver, depends on the carrier frequency, so we have Finally the error expression of the code delay, which is a function of the relative speed satellite receiver at the ith millisecond, is given by In this expression f e i = i j=1 α j−1 Tc j−1 is the error on the frequency associated with the movement of the receiver and f r is the intermediate frequency of the emitted satellite signal.Parameters of the received signal are estimated with the values of I(t) and Q(t) integrated on one or several periods Tc of the signal.
Let N be the number of samples integrated at the ith period Tc, we define for each satellite s [4], ( In these expressions R(• • • ) is the normalized correlation function of the C/A CDMA code.sin c(• • • ) is a function that models the error on the estimate of the carrier frequency.In general the error values of φe i and f e i are estimated in a filter with I i (0) and Q i (0).The values of τe i are estimated in a filter from the in-phase and quad phase observations . Tb is the length of one bit of the CDMA code.Finally the value of S/N 0 is obtained when I i (0) is maximum, so for f e i = 0, τe i = 0 and φe i = 0.

Problem Positioning.
The discriminator "early-minuslate" is defined by the following expression for which we fix to Tb the code delay between early and late code generated by the receiver: In this expression η D i is a white Gaussian noise of variance 2. The function of correlation of the CDMA code is given by the following expression: 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 |τe i | < Tb/2.The variation strictly decreases for 3Tb/2 > |τe i | > Tb/2 and stays equal to zero for |τe i | > 3Tb/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 |τe i | ≤ Tb/2, the discriminator behaviour is proportional to the code delay τe i and is defined by When the code is retiming with a code delay of Tb/2 every pseudorange variation of CTb/2, the discriminator behaviour is approximately piecewise linear (C 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.We suppose here that the receiver trajectory does not modify the strict increasing (when the satellite moves away) i (ms) 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 Tb/2.When the code is retimed we use the following constant acceleration model to estimate the pseudorange satellite receiver d i : In this expression v i is the relative speed satellite-receiver.d t0 and v t0 are, respectively, the initial pseudo range and speed, processed in the acquisition step of the GPS signal.In this expression a is the acceleration supposed to be constant in each stationary piece.The value of a is defined in order to have d t1 = C Tb/2 + d t0 , where t 1 , as shown on Figure 1, is the instant such as τe t1 = Tb/2.Then we have the following expression for a in the first segment: The phase and frequency values that minimise f e i and φe i 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.

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 y = (y i , 1 ≤ i ≤ n) be the observed discriminator values model as a series with conditional law P(y).We suppose that the discriminator variations are piecewise stationary.Then we note r = (r i , 1 ≤ i ≤ n) the instants of change for which the code is retimed of Tb/2.r i 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 r n = 1.The observed series is composed of K r segments, then we have K r = n i=1 r i .We notice (n k , k ≥ 0) the number of samples in the segment k.We notice t = (t k , k ≥ 0) the change instants in the observed series, associated with the change times configuration r (K r is the number of changes).We define the origin instant at t 0 = 0 and the last change at t N = n.
Let the following expression for the discriminator at instant i in the segment k: In the simplified (11), s k is the maximum discriminator value, and it depends on the SNR.We note s = (s k , k ≥ 0) the set of maximum values of the discriminator.A i,k is a strictly increasing function (or decreasing function according to the satellite trajectory) which depends on the values of f e i , φe i , τe i .The distribution that describes the noisy values of the discriminator can be defined by [14], P r, y, s = P s y , r P y r π(r), (12) where π(r) is the prior law of the change configuration r.
Moreover we have where h(• • • ) is the likelihood function.Finally the maximum a posteriori estimate of the change configuration is defined by 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 s is defined as a random variable of mean μ and variance V .The mean value of the SNR μ is defined in the acquisition process and V 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).σ d represents the inaccuracy of this prior information.So, where μ, V and σ d 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 r (the number and position of changes) is estimated by minimizing the penalized contrast function defined by 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.

Contrast Function.
The likelihood distribution, of the discriminator defined in ( 8) and ( 9), is given by with σ 2 = 2.The stage of acquisition provides an inaccurate estimation s k of the SNR.Furthermore, the SNR varies with time.Then we suppose s k as normally distributed random variables with mean μ and variance V/n k [14].The value of μ is given by the acquisition module, and the value of V is fixed by the user.This value will be chosen as high as μ is inaccurate.Then we have We show in Appendix A that with To define the values of A i,k , we assume the constant acceleration dynamic model in the stationary pieces defined by (9).This model describes the variations of the discriminator y i in the stationary pieces.Then we have the following equations for A i,k in the segment k and for i ∈ {t k−1 + 1, t k }: with for each segment k: The model is initialized by where v 0 is the relative speed satellite receiver and μ the SNR provided by the acquisition process.

Penalty Definition.
We define r 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 r is given by [10] 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 π(r) to have a change configuration r as where λ i is the probability to have a change at instant i.Then we have the following penalty term: 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 t by the following normal distribution: where td k , the change instant for the segment k, is processed with the Doppler frequency f d as follow: σ 2 d 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 λ i at instant i: where I i (x) is the modified Bessel function, x is a parameter of concentration, and * the function of convolution.In this expression we define K r discrete Gaussian distributions centered on the instants td = (td k , k ≥ 0).We show in Figure 2 the evolution of the penalty and its associate probability.In this example we have td = (td 0 = 0, td 1 = 30, td 2 = 70, td 3 = 90).Finally we show in Appendix B that the penalized contrast function to minimize is given by 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.

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 r solution (that minimizes (30)) with probability one [15].For our application, the simulated annealing algorithm proceeds as follows.
(1) An initial random sequence r is drawn.Choose an initial value for the parameter T (T j , j = 0), called temperature (in our application T 0 = 500).j = 0.
(2) A new configuration r is drawn from r in one of these three cases: (i) add a change at a random position i in r, (ii) remove a change at a random position i in r, (iii) translate a change at a random position i in r.
(3) Process the penalized contrast function in three steps: (i) compute the A i,k with the change configuration r, (ii) compute the y i , φe i , and f e i for i ∈ {1, . . ., n} with the change configuration r.We use in this case the architecture tracking process described in Figure 3, (iii) compute U y ( r).
(5) Set r = r with probability one if ΔU y < 0 and with probability exp{−ΔU y /T j } elsewhere.Decrease the value of T j+1 = f (T j ) (in our application T j+1 = T j * 0.99).
When r has been estimated with the simulated annealing algorithm, we use the following equations to compute the pseudorange: with for each segment k: The initial relative speed satellite receiver v 0 is provided by the acquisition process.

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 L1 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 x and V 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 x will increase with the relative speed satellite receiver in the dynamic case.The value of V will increase when the signal to noise decreases because in this case the acquisition process provides inaccurate estimate of S/N 0 .
We show in Figure 4 positions computed every millisecond in the static noiseless case, for a signal-to-noise ratio of S/N 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.
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-tonoise ratios with the standard and the proposed method.
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.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.

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.

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.
Positions obtained with novatel GPS Positions obtained with our method 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.

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 r (the number and position of changes) is estimated by minimizing the following penalized contrast function: ( r) = Argmin In this expression V y (r; μ, V ) is the negative log of the likelihood distribution defined by ( 19) and is given by

Figure 5 :
Figure 5: Estimate dynamic trajectory for a S/N 0 of 50 dB.Hz.

Table 3 :Figure 6 :
Figure 6: Positions calculated in the static real case.

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

2 )( 1 − 4 )
In the same way − log π(r; σ d ) can be rewritten as follows:− log π(r; σ d ) λ i ).(B.3)Finally the expression of the penalized contrast function isV y r; μ, V − log π(r; σ d )In this expression the last three terms are independent of the change configuration r.In this context the penalized contrast function to minimize in order to estimate r is given by ( r) = Argmin (r)∈{0,

Table 2 :
Error on estimated static positions.