Research Article  Open Access
G. Stienne, S. Reboul, M. Azmani, S. Boutoille, J. B. Choquel, M. Benjelloun, "Bayesian ChangePoints Estimation Applied to GPS Signal Tracking", International Scholarly Research Notices, vol. 2011, Article ID 148242, 11 pages, 2011. https://doi.org/10.5402/2011/148242
Bayesian ChangePoints Estimation Applied to GPS Signal Tracking
Abstract
A hierarchical Bayesian model is applied to offline 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 highrate 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 signaltonoise 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 and or every 20 ms at frequency [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 signaltonoise 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 online 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 and have a dataless pilot channel (the code at frequency, the code at 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 offline 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 meansquare 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 satellitereceiver pseudorange evolution. In order to compute the navigation solution, the receiver multiplies the received signal with inphase and quadphase replica of the GPS signal:
In these expressions is the signal amplitude; it is normalized to drive the noise variance of , to 1. It is a function of the signaltonoise 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] In these expressions and are, respectively, the noise clock disturbances on the phase and frequency. 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 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 In this expression 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],
In these expressions is the normalized correlation function of the C/A CDMA code. 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 and . The values of are estimated in a filter from the inphase and quad phase observations , and , [4]. is the length of one bit of the CDMA code. Finally the value of is obtained when is maximum, so for , and .
2.2. Problem Positioning
The discriminator “earlyminuslate” is defined by the following expression for which we fix to the code delay between early and late code generated by the receiver: 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: 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 earlyminuslate strictly increases for . The variation strictly decreases for and stays equal to zero for . 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 , the discriminator behaviour is proportional to the code delay and is defined by When the code is retiming with a code delay of every pseudorange variation of , 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)
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 . When the code is retimed we use the following constant acceleration model to estimate the pseudorange satellite receiver : In this expression is the relative speed satellitereceiver. and 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 , where , as shown on Figure 1, is the instant such as . Then we have the following expression for in the first segment: 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 onedimensional process. Let be the observed discriminator values model as a series with conditional law . We suppose that the discriminator variations are piecewise stationary. Then we note the instants of change for which the code is retimed of . 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 . The observed series is composed of segments, then we have . We notice the number of samples in the segment . We notice the change instants in the observed series, associated with the change times configuration ( is the number of changes). We define the origin instant at and the last change at .
Let the following expression for the discriminator at instant in the segment : In the simplified (11), is the maximum discriminator value, and it depends on the SNR. We note 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], where is the prior law of the change configuration . Moreover we have where 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 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, 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 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 with . 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 We show in Appendix A that with 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 : with for each segment : The model is initialized by 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] 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 where is the probability to have a change at instant . 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 by the following normal distribution: where , the change instant for the segment , is processed with the Doppler frequency as follow: 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 : 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 . We show in Figure 2 the evolution of the penalty and its associate probability. In this example we have .
(a)
(b)
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.
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 (), called temperature (in our application ). . (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 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 and with probability elsewhere. Decrease the value of (in our application ). (6). Go to (2).
When has been estimated with the simulated annealing algorithm, we use the following equations to compute the pseudorange: with for each segment : 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 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 temperaturecompensated crystal oscillator clocks. The values of the phase and frequency randomwalk intensities are and .

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 signaltonoise ratio of 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 signaltonoise 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 signaltonoise ratio decreases.

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.
5. Conclusions and Future Works
In this paper, a Bayesian offline 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 This expression can be rewritten as follows if we introduce prior information: Then we have defined with (17) and (18) as follows: Let: Then we have then we can deduce the following distributions:
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: In this expression is the negative log of the likelihood distribution defined by (19) and is given by In the same way can be rewritten as follows: Finally the expression of the penalized contrast function is 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
References
 G. Gerten, “Protecting the global positioning system,” IEEE Aerospace and Electronic Systems Magazine, vol. 20, no. 11, pp. 3–8, 2005. View at: Publisher Site  Google Scholar
 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, LongBeach, Calif, USA, September 2004. View at: Google Scholar
 M. G. Petovello, C. O'Driscoll, and G. Lachapell, “Weak signal carrier tracking using extended coherent integration with an ultratight GNSS/IMU receiver,” in Proceedings of the European Navigation Conference (ENCGNSS '08), Toulouse, France, April 2008. View at: Google Scholar
 N. I. Ziedan, GNSS Receivers for Weak Signals, Artech House, London, UK, 2006.
 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
 S. Boutoille, S. Reboul, and M. Benjelloun, “A hybrid fusion system applied to offline detection and changepoints estimation,” Information Fusion, vol. 11, no. 4, pp. 325–337, 2010. View at: Publisher Site  Google Scholar
 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 Site  Google Scholar
 E. Lebarbier, “Detecting multiple changepoints in the mean of Gaussian process by model selection,” Signal Processing, vol. 85, no. 4, pp. 717–736, 2005. View at: Publisher Site  Google Scholar
 Y. C. Yao, “Estimating the number of changepoints via Schwarz criterion,” Statistics and Probability Letters, vol. 6, no. 3, pp. 181–189, 1988. View at: Google Scholar
 M. Lavielle, “Optimal segmentation of random processes,” IEEE Transactions on Signal Processing, vol. 46, no. 5, pp. 1365–1373, 1998. View at: Google Scholar
 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 Site  Google Scholar
 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. View at: Google Scholar
 E. Kaplan, Understanding GPS: Principles and Applications, Artech House, Norwood, Mass, USA, 1996.
 M. Lavielle and E. Lebarbier, “An application of MCMC methods for the multiple changepoints problem,” Signal Processing, vol. 81, no. 1, pp. 39–53, 2001. View at: Publisher Site  Google Scholar
 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 Site
Copyright
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.