# Bayesian Change-Points Estimation Applied to GPS Signal Tracking

Authors: G. Stienne; S. Reboul; M. Azmani; S. Boutoille; J. B. Choquel; M. Benjelloun
Journal: ISRN Signal Processing
 (2011)

---

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


---

## Body

## 1. Introduction

A GPS receiver can provide a position every millisecond at frequenciesfL1 and fL5 or every 20 ms at frequency fL2 [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 signalsL2C and L5 have a dataless pilot channel (the CL code at fL2 frequency, the Q5 code at fL5 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 RatioS/N0. 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)Cs(⋯) modulated at frequency fs. 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: (1)Is(t)=AsCs(t-τs)Cs(t-τ̂s)×cos⁡(2π(fs-f̂s)t+ϕs-ϕ̂s)+ηsI(t),Qs(t)=AsCs(t-τs)Cs(t-τ̂s)×sin(2π(fs-f̂s)t+ϕs-ϕ̂s)+ηsQ(t).In these expressionsAs is the signal amplitude; it is normalized to drive the noise variance of ηsI(t), ηsQ(t) to 1. It is a function of the signal-to-noise ratio (S/N0). 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, fe=fs-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)fei=fei-1+αi-1Tci-1+wfi,ϕei=ϕei-1+2πfeiTci+wϕi.
In these expressions wfi 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. Tci, the code length defined at the receiver, depends on the carrier frequency, so we have (3)Tci=Tcfsf̂s+fei.
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 (4)τei=τei-1+Tcifr-(f̂s+f̃ei)fr.
In this expression f̃ei=∑j=1iαj-1Tcj-1 is the error on the frequency associated with the movement of the receiver and fr 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.LetN be the number of samples integrated at the ith period Tc, we define for each satellite s [4],  (5)Ii(Δ)=1N∑k=1NIs(t)=2SN0dtR(τei-Δ)×sinc(fei+αiTci2)cos⁡(ϕei)+ηiI,Qi(Δ)=1N∑k=1NQs(t)=2SN0dtR(τei-Δ)×sinc(fei+αiTci2)sin(ϕei)+ηiQ.In these expressionsR(⋯) is the normalized correlation function of the C/A CDMA code. sinc(⋯) is a function that models the error on the estimate of the carrier frequency. In general the error values of ϕei and fei are estimated in a filter with Ii(0) and Qi(0). The values of τei are estimated in a filter from the in-phase and quad phase observations Iei=Ii(-Tb2),Ili=Ii(Tb2), and Qei=Qi(-Tb/2), Qli=Qi(Tb/2) [4]. Tb is the length of one bit of the CDMA code. Finally the value of S/N0 is obtained when Ii(0) is maximum, so for fei=0, τei=0 and ϕei=0.


### 2.2. Problem Positioning

The discriminator “early-minus-late” is defined by the following expression for which we fix toTb the code delay between early and late code generated by the receiver:(6)Di=Iei-Ili=2SN0dtsinc(fei+αiTci2)×cos⁡(ϕei)(R(τei-Tb2)-R(τei+Tb2))+ηiD.
In this expression ηiD is a white Gaussian noise of variance 2. The function of correlation of the CDMA code is given by the following expression: (7)R(τ)=(1-|τTb|).
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|τei|<Tb/2. The variation strictly decreases for 3Tb/2>|τei|>Tb/2 and stays equal to zero for |τei|>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 |τei|≤Tb/2, the discriminator behaviour is proportional to the code delay τei and is defined by(8)Di=22SN0dtsinc(fei+αiTci2)×cos⁡(ϕei)(|τeiTb+12|-|τeiTb-12|)+ηiD,Di=22SN0dtsinc(fei+αiTci2)cos⁡(ϕei)τeiTb+ηiD.
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.Delay open loop discriminator behavior.
(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 ofTb/2.  When the code is retimed we use the following constant acceleration model to estimate the pseudorange satellite receiver di: (9)di=di-1+Tcvi-1+Tc22a,vi=vi-1+Tca.
In this expression vi is the relative speed satellite-receiver. dt0 and vt0 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 dt1=CTb/2+dt0, where t1, as shown on Figure 1, is the instant such as τet1=Tb/2. Then we have the following expression for a in the first segment: (10)a=CTb-2(t1-t0)Tcv0Tc2(t1-t0)2.
The phase and frequency values that minimise fei  and ϕei  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 Figure1, 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.


## 2.1. GPS Model

The GPS satellite signal consists of a code division multiple access (CDMA)Cs(⋯) modulated at frequency fs. 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: (1)Is(t)=AsCs(t-τs)Cs(t-τ̂s)×cos⁡(2π(fs-f̂s)t+ϕs-ϕ̂s)+ηsI(t),Qs(t)=AsCs(t-τs)Cs(t-τ̂s)×sin(2π(fs-f̂s)t+ϕs-ϕ̂s)+ηsQ(t).In these expressionsAs is the signal amplitude; it is normalized to drive the noise variance of ηsI(t), ηsQ(t) to 1. It is a function of the signal-to-noise ratio (S/N0). 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, fe=fs-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)fei=fei-1+αi-1Tci-1+wfi,ϕei=ϕei-1+2πfeiTci+wϕi.
In these expressions wfi 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. Tci, the code length defined at the receiver, depends on the carrier frequency, so we have (3)Tci=Tcfsf̂s+fei.
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 (4)τei=τei-1+Tcifr-(f̂s+f̃ei)fr.
In this expression f̃ei=∑j=1iαj-1Tcj-1 is the error on the frequency associated with the movement of the receiver and fr 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.LetN be the number of samples integrated at the ith period Tc, we define for each satellite s [4],  (5)Ii(Δ)=1N∑k=1NIs(t)=2SN0dtR(τei-Δ)×sinc(fei+αiTci2)cos⁡(ϕei)+ηiI,Qi(Δ)=1N∑k=1NQs(t)=2SN0dtR(τei-Δ)×sinc(fei+αiTci2)sin(ϕei)+ηiQ.In these expressionsR(⋯) is the normalized correlation function of the C/A CDMA code. sinc(⋯) is a function that models the error on the estimate of the carrier frequency. In general the error values of ϕei and fei are estimated in a filter with Ii(0) and Qi(0). The values of τei are estimated in a filter from the in-phase and quad phase observations Iei=Ii(-Tb2),Ili=Ii(Tb2), and Qei=Qi(-Tb/2), Qli=Qi(Tb/2) [4]. Tb is the length of one bit of the CDMA code. Finally the value of S/N0 is obtained when Ii(0) is maximum, so for fei=0, τei=0 and ϕei=0.


## 2.2. Problem Positioning

The discriminator “early-minus-late” is defined by the following expression for which we fix toTb the code delay between early and late code generated by the receiver:(6)Di=Iei-Ili=2SN0dtsinc(fei+αiTci2)×cos⁡(ϕei)(R(τei-Tb2)-R(τei+Tb2))+ηiD.
In this expression ηiD is a white Gaussian noise of variance 2. The function of correlation of the CDMA code is given by the following expression: (7)R(τ)=(1-|τTb|).
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|τei|<Tb/2. The variation strictly decreases for 3Tb/2>|τei|>Tb/2 and stays equal to zero for |τei|>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 |τei|≤Tb/2, the discriminator behaviour is proportional to the code delay τei and is defined by(8)Di=22SN0dtsinc(fei+αiTci2)×cos⁡(ϕei)(|τeiTb+12|-|τeiTb-12|)+ηiD,Di=22SN0dtsinc(fei+αiTci2)cos⁡(ϕei)τeiTb+ηiD.
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.Delay open loop discriminator behavior.
(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 ofTb/2.  When the code is retimed we use the following constant acceleration model to estimate the pseudorange satellite receiver di: (9)di=di-1+Tcvi-1+Tc22a,vi=vi-1+Tca.
In this expression vi is the relative speed satellite-receiver. dt0 and vt0 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 dt1=CTb/2+dt0, where t1, as shown on Figure 1, is the instant such as τet1=Tb/2. Then we have the following expression for a in the first segment: (10)a=CTb-2(t1-t0)Tcv0Tc2(t1-t0)2.
The phase and frequency values that minimise fei  and ϕei  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 Figure1, 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 aone-dimensional process. Let y̲=(yi,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̲=(ri,1≤i≤n) the instants of change for which the code is retimed of Tb/2. ri 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 rn=1. The observed series is composed of Kr segments, then we have Kr=∑i=1nri. We notice (nk,k≥0) the number of samples in the segment k. We notice t̲=(tk,k≥0) the change instants in the observed series, associated with the change times configuration r̲ (Kr is the number of changes). We define the origin instant at t0=0 and the last change at tN=n.Let the following expression for the discriminator at instanti in the segment k: (11)yi=skAi,k.
In the simplified (11), sk is the maximum discriminator value, and it depends on the SNR. We note s̲=(sk,k≥0) the set of maximum values of the discriminator. Ai,k is a strictly increasing function (or decreasing function according to the satellite trajectory) which depends on the values of fei, ϕei, τei. The distribution that describes the noisy values of the discriminator can be defined by [14], (12)P(r̲,y̲,s̲)=P(s̲y̲,r̲)P(y̲r̲)π(r̲),
where π(r̲) is the prior law of the change configuration r̲. Moreover we have (13)P(s̲y̲,r̲)h(y̲r̲)=h(y̲r̲,s̲)P(s̲r̲),
where h(⋯) is the likelihood function. Finally the maximum a posteriori estimate of the change configuration is defined by (14)(r̲̂)=Argmax︸(r̲)P(r̲∣y̲).
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, (15)(r̲̂)=Argmin︸(r̲)-log⁡(h(y̲∣r̲;μ,V)π(r̲;σd)),
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 (16)Uy̲(r̲)=Vy̲(r̲;μ,V)-log⁡π(r̲;σd).
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  (17)h(y̲∣r̲,s̲;σ2)=∏k=1Kr∏i=tk-1+1tk12πσ2exp⁡(-(yi-skAi,k)22σ2),
with σ2=2. The stage of acquisition provides an inaccurate estimation sk of the SNR. Furthermore, the SNR varies with time. Then we suppose sk as normally distributed random variables with mean μ and variance V/nk [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 (18)P(s̲∣r̲;μ,V)=∏k=1Kr2πVnkexp⁡(-nk(sk-μ)22V).
We show in Appendix A that(19)h(y̲∣r̲;μ,V)=(14π)n/2∏k=1KrVknkV×exp⁡(-12(12∑i=tk-1+1tkyi2+nkμ2V-μk2Vk)),
with(20)Tk=∑i=tk-1+1tkAi,k,Vk=1(Tk/2+nk/V),μk=Vk(μnkV+12∑i=tk-1+1tkyiAi,k).
To define the values of Ai,k, we assume the constant acceleration dynamic model in the stationary pieces defined by (9). This model describes the variations of the discriminator yi in the stationary pieces. Then we have the following equations for Ai,k in the segment k and for i∈{tk-1+1,tk}: (21)Ai,k=Ai-1,k+Tcṽi-1+Tc22ak,ṽi=ṽi-1+Tcak,
with for each segment k:(22)ak=2-2Tcnkṽtk-1Tc2nk2,Atk-1+1,k=ṽtk-1+Tcak.
The model is initialized by (23)ṽt0=v0μ,
where v0 is the relative speed satellite receiver and μ the SNR provided by the acquisition process.


### 3.3. Penalty Definition

We definer̲ 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] (24)π(r̲)=λKr(1-λ)(n-Kr).
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 (25)π(r̲)=∏i=1/ri=1nλi∏i=1/ri=0n(1-λi),
where λi is the probability to have a change at instant i. Then we have the following penalty term: (26)ln⁡(π(r̲))=∑i=1nriln⁡(λi1-λi)+∑i=1nln⁡(1-λi).
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: (27)λ(t)=12πσd2∑k=1Krexp⁡(-(t-tdk)22σd2),
where tdk, the change instant for the segment k, is processed with the Doppler frequency fd as follow: (28)tdk=tdk-1+TcfL22fd,td0=0.σd2 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: (29)λi=∑k=1Kre-xIi(x)*δ(i-tdk),
where Ii(x) is the modified Bessel function, x is a parameter of concentration, and * the function of convolution. In this expression we define Kr discrete Gaussian distributions centered on the instants td̲=(tdk,k≥0). We show in Figure 2 the evolution of the penalty and its associate probability. In this example we have td̲=(td0=0,td1=30,td2=70,td3=90).(a) Probability evolution and (b) penalty evolution.
(a)(b)Finally we show in AppendixB that the penalized contrast function to minimize is given by (30)(r̲̂)=Argmin︸(r̲)∈{0,1}12∑k=1Kr{12∑i=tk-1+1tkyi2-μk2Vk-log⁡Vk}-∑k=1Krlog⁡λtk(1-λtk).
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 r̲̂ solution (that minimizes (30)) with probability one [15]. For our application, the simulated annealing algorithm proceeds as follows.  (1)
An initial random sequencer̲̂ is drawn. Choose an initial value for the parameter T (Tj,j=0), called temperature (in our application T0=500). j=0.(2)
A new configurationr̲̃ is drawn from r̲̂ in one of these three cases:
(i)
add a change at a random positioni in r̲̂,(ii)
remove a change at a random positioni in r̲̂,(iii)
translate a change at a random positioni in r̲̂.(3)
Process the penalized contrast function in three steps:
(i)
compute theAi,k with the change configuration r̲̃,(ii)
compute theyi, ϕei, and fei for i∈{1,…,n} with the change configuration r̲̃. We use in this case the architecture tracking process described in Figure 3,(iii)
computeUy̲(r̲̃).(4)
LetΔUy̲=Uy̲(r̲̃)-Uy̲(r̲̂).(5)
Setr̲̂=r̲̃ with probability one if ΔUy̲<0 and with probability exp{-ΔUy̲/Tj} elsewhere. Decrease the value of Tj+1=f(Tj) (in our application Tj+1=Tj*0.99).(6)
j=j+1. Go to (2).Figure 3
Architecture of the tracking process.Whenr̲̂ has been estimated with the simulated annealing algorithm, we use the following equations to compute the pseudorange:(31)di=di-1+Tcvi-1+Tc22ak,vi=vi-1+Tcak,
with for each segment k: (32)ak=CTb-2Tcnkvtk-1Tc2nk2.
The initial relative speed satellite receiver v0 is provided by the acquisition process.


## 3.1. Bayesian Approach

The problem we address in this paper can be defined as the detection of changes in the statistical distribution of aone-dimensional process. Let y̲=(yi,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̲=(ri,1≤i≤n) the instants of change for which the code is retimed of Tb/2. ri 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 rn=1. The observed series is composed of Kr segments, then we have Kr=∑i=1nri. We notice (nk,k≥0) the number of samples in the segment k. We notice t̲=(tk,k≥0) the change instants in the observed series, associated with the change times configuration r̲ (Kr is the number of changes). We define the origin instant at t0=0 and the last change at tN=n.Let the following expression for the discriminator at instanti in the segment k: (11)yi=skAi,k.
In the simplified (11), sk is the maximum discriminator value, and it depends on the SNR. We note s̲=(sk,k≥0) the set of maximum values of the discriminator. Ai,k is a strictly increasing function (or decreasing function according to the satellite trajectory) which depends on the values of fei, ϕei, τei. The distribution that describes the noisy values of the discriminator can be defined by [14], (12)P(r̲,y̲,s̲)=P(s̲y̲,r̲)P(y̲r̲)π(r̲),
where π(r̲) is the prior law of the change configuration r̲. Moreover we have (13)P(s̲y̲,r̲)h(y̲r̲)=h(y̲r̲,s̲)P(s̲r̲),
where h(⋯) is the likelihood function. Finally the maximum a posteriori estimate of the change configuration is defined by (14)(r̲̂)=Argmax︸(r̲)P(r̲∣y̲).
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, (15)(r̲̂)=Argmin︸(r̲)-log⁡(h(y̲∣r̲;μ,V)π(r̲;σd)),
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 (16)Uy̲(r̲)=Vy̲(r̲;μ,V)-log⁡π(r̲;σd).
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  (17)h(y̲∣r̲,s̲;σ2)=∏k=1Kr∏i=tk-1+1tk12πσ2exp⁡(-(yi-skAi,k)22σ2),
with σ2=2. The stage of acquisition provides an inaccurate estimation sk of the SNR. Furthermore, the SNR varies with time. Then we suppose sk as normally distributed random variables with mean μ and variance V/nk [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 (18)P(s̲∣r̲;μ,V)=∏k=1Kr2πVnkexp⁡(-nk(sk-μ)22V).
We show in Appendix A that(19)h(y̲∣r̲;μ,V)=(14π)n/2∏k=1KrVknkV×exp⁡(-12(12∑i=tk-1+1tkyi2+nkμ2V-μk2Vk)),
with(20)Tk=∑i=tk-1+1tkAi,k,Vk=1(Tk/2+nk/V),μk=Vk(μnkV+12∑i=tk-1+1tkyiAi,k).
To define the values of Ai,k, we assume the constant acceleration dynamic model in the stationary pieces defined by (9). This model describes the variations of the discriminator yi in the stationary pieces. Then we have the following equations for Ai,k in the segment k and for i∈{tk-1+1,tk}: (21)Ai,k=Ai-1,k+Tcṽi-1+Tc22ak,ṽi=ṽi-1+Tcak,
with for each segment k:(22)ak=2-2Tcnkṽtk-1Tc2nk2,Atk-1+1,k=ṽtk-1+Tcak.
The model is initialized by (23)ṽt0=v0μ,
where v0 is the relative speed satellite receiver and μ the SNR provided by the acquisition process.


## 3.3. Penalty Definition

We definer̲ 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] (24)π(r̲)=λKr(1-λ)(n-Kr).
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 (25)π(r̲)=∏i=1/ri=1nλi∏i=1/ri=0n(1-λi),
where λi is the probability to have a change at instant i. Then we have the following penalty term: (26)ln⁡(π(r̲))=∑i=1nriln⁡(λi1-λi)+∑i=1nln⁡(1-λi).
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: (27)λ(t)=12πσd2∑k=1Krexp⁡(-(t-tdk)22σd2),
where tdk, the change instant for the segment k, is processed with the Doppler frequency fd as follow: (28)tdk=tdk-1+TcfL22fd,td0=0.σd2 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: (29)λi=∑k=1Kre-xIi(x)*δ(i-tdk),
where Ii(x) is the modified Bessel function, x is a parameter of concentration, and * the function of convolution. In this expression we define Kr discrete Gaussian distributions centered on the instants td̲=(tdk,k≥0). We show in Figure 2 the evolution of the penalty and its associate probability. In this example we have td̲=(td0=0,td1=30,td2=70,td3=90).(a) Probability evolution and (b) penalty evolution.
(a)(b)Finally we show in AppendixB that the penalized contrast function to minimize is given by (30)(r̲̂)=Argmin︸(r̲)∈{0,1}12∑k=1Kr{12∑i=tk-1+1tkyi2-μk2Vk-log⁡Vk}-∑k=1Krlog⁡λtk(1-λtk).
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 r̲̂ solution (that minimizes (30)) with probability one [15]. For our application, the simulated annealing algorithm proceeds as follows.  (1)
An initial random sequencer̲̂ is drawn. Choose an initial value for the parameter T (Tj,j=0), called temperature (in our application T0=500). j=0.(2)
A new configurationr̲̃ is drawn from r̲̂ in one of these three cases:
(i)
add a change at a random positioni in r̲̂,(ii)
remove a change at a random positioni in r̲̂,(iii)
translate a change at a random positioni in r̲̂.(3)
Process the penalized contrast function in three steps:
(i)
compute theAi,k with the change configuration r̲̃,(ii)
compute theyi, ϕei, and fei for i∈{1,…,n} with the change configuration r̲̃. We use in this case the architecture tracking process described in Figure 3,(iii)
computeUy̲(r̲̃).(4)
LetΔUy̲=Uy̲(r̲̃)-Uy̲(r̲̂).(5)
Setr̲̂=r̲̃ with probability one if ΔUy̲<0 and with probability exp{-ΔUy̲/Tj} elsewhere. Decrease the value of Tj+1=f(Tj) (in our application Tj+1=Tj*0.99).(6)
j=j+1. Go to (2).Figure 3
Architecture of the tracking process.Whenr̲̂ has been estimated with the simulated annealing algorithm, we use the following equations to compute the pseudorange:(31)di=di-1+Tcvi-1+Tc22ak,vi=vi-1+Tcak,
with for each segment k: (32)ak=CTb-2Tcnkvtk-1Tc2nk2.
The initial relative speed satellite receiver v0 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 realL1 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 Table1 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 Sf=5×10-21s and Sg=5.9×10-20s-1.Table 1
Simulation parameters.
HDOPNumber of satellitesS/N0 (dB·Hz)x: Prior informationV: Prior information1.43750 / 40 / 301002x 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/N0.We show in Figure4 positions computed every millisecond in the static noiseless case, for a signal-to-noise ratio of S/N0=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 aS/N0 of 50 dB·Hz.We report in Table2 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.
Methodproposed/standardMean gap in east directionMean gap in north directionMean distanceS/N0=50 dB·Hz0.5/6.71.4/11.51.3/14.6S/N0=40 dB·Hz0.7/10.71/19.91.6/24.6S/N0=30 dB·Hz1.5/23.41.9/38.82.7/49.6For the dynamic case, we consider a circular trajectory of radius 10 m ended in 5 s. We show in Figure5, 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 aS/N0 of 50 dB.Hz.We report in Table3 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.
Methodproposed/standardMean gap in east directionMean gap in north directionMean distanceS/N0=50 dB.Hz0.5/8.41/14.80.8/18.5S/N0=40 dB.Hz0.6/11.51.1/21.81.2/26.7S/N0=30 dB.Hz1.6/2510.5/50.92.8/61.4


### 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 Figures6 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).


## 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 Table1 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 Sf=5×10-21s and Sg=5.9×10-20s-1.Table 1
Simulation parameters.
HDOPNumber of satellitesS/N0 (dB·Hz)x: Prior informationV: Prior information1.43750 / 40 / 301002x 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/N0.We show in Figure4 positions computed every millisecond in the static noiseless case, for a signal-to-noise ratio of S/N0=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 aS/N0 of 50 dB·Hz.We report in Table2 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.
Methodproposed/standardMean gap in east directionMean gap in north directionMean distanceS/N0=50 dB·Hz0.5/6.71.4/11.51.3/14.6S/N0=40 dB·Hz0.7/10.71/19.91.6/24.6S/N0=30 dB·Hz1.5/23.41.9/38.82.7/49.6For the dynamic case, we consider a circular trajectory of radius 10 m ended in 5 s. We show in Figure5, 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 aS/N0 of 50 dB.Hz.We report in Table3 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.
Methodproposed/standardMean gap in east directionMean gap in north directionMean distanceS/N0=50 dB.Hz0.5/8.41/14.80.8/18.5S/N0=40 dB.Hz0.6/11.51.1/21.81.2/26.7S/N0=30 dB.Hz1.6/2510.5/50.92.8/61.4


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



---

*Source: 148242-2011-06-23.xml*